text:
function [clusters,CLU,subclusters] = mask2clusters(P,varargin)
% Extracts clusters and con img data from mask
%
% Use with *mask_intersection.m*
%
% To get clusters but not extract data, enter only one argument.
%
% To get clusters and choose extraction imgs with the GUI, enter an empty [] 2nd argument.
%
% :Usage:
% ::
%
%    [clusters,CLU,subclusters] = mask2clusters(img mask file with voxels,[imgs to extract data from],[df])
%
%
% DOES *NOT* CONVERT BETWEEN DIFFERENT VOXEL SIZES AND POSITIONS BETWEEN IMNAMES AND SPM/VOL STRUCTS
%
% :See also: roi_probe
%
% If no imgs are entered, Z-scores are values from mask
%
% If df is entered, values in mask img are converted to Z-scores with spm_t2z.m
%
% If extract img names are empty and df is entered, assume we're using
% values from mask as t-values and convert to Z-scores
%
% WARNING: for spm2 compatibility, ABSOLUTE VALUES of voxel sizes are
% returned; e.g., ignores analyze flipping in SPM2.
%
% % Matlab 6.5/OSX bug gives seg fault or something if mask is too big.
%
% :Example:
% ::
%
%    cl = mask2clusters('myimage.img',[img string mtx],[]); % no z-score
%    conversion, extracts data from [img string mtx]
%
%    cl = mask2clusters('rob_tmap_0002_filt_t_3-05_k10_neg.img')
%
%    % This one works with already-loaded image data and a mat matrix:
%    V = spm_vol('rob_tmap_0002_filt_t_3-05_k10_neg.img'); dat = spm_read_vols(V);
%    cl = mask2clusters(dat,V.mat);
%
%
% ..
%    tor wager
%
%    modification 2/27/03
% ..

% ..
%    set up inputs
% ..
clusters = [];
df = []; imP = [];  resl = 0; vmat = [];

for i = 1:length(varargin)
    tmp = varargin{i};
    if isstr(tmp), imP = tmp;       % images to extract data from
    elseif isstruct(tmp), vmat = tmp.mat;  % it's an spm_vol V struct
    elseif isscalar(tmp), df = tmp; % df to convert from t to z-scores
    elseif all(size(tmp) > 1), vmat = tmp;      % it's a SPM .mat matrix to be used with raw data instead of string img name
    end
end

if isempty(P)
    P = spm_get(1,'*img','Select mask image with voxels to extract');
end

if ischar(P)
    P = deblank(P);

    if isempty(which(P)) && ~exist(P,'file')
        disp(['Cannot find file: ' P]);
        return
    end

    if ~isempty(which(P)), P = which(P); end
end

% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% get coordinates and height values from the mask
% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

% img2voxel can take raw data, if vmat is entered; if string, vmat not
% used.
[CLU.XYZ,CLU.XYZmm,CLU.Z,CLU.V] = img2voxel(P,vmat);
CLU.XYZ = CLU.XYZ(1:3,:);

% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% change pseudo-T  or T values to Z scores
% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
if ~isempty(imP)
    if isempty(df) & length(varargin) > 1, df = size(imP,1) - 1;, end     % to get df
    if df == 0, df = [];, end
end

if ~isempty(df)
    disp(['Converting to Z scores based on ' num2str(df) ' df.'])
    [CLU.Z] = spm_t2z(CLU.Z,df);
else
    df = 0;
    %disp('Saving values in mask file in clusters.Z (no z-score conversion)')
    %CLU.Z = ones(1,size(CLU.XYZ,2));
end

if size(CLU.Z,1) > size(CLU.Z,2), CLU.Z = CLU.Z';, end

% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% fill in other fields
% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

% adjust crit_t from 0 to 1 for compatibility with fixed_TSU (Talairach)
% needs to be a non-zero value to display on Talairach
% ------------------------------------------------------
CLU.crit_t = 1;

CLU.cl_size = 0;
CLU.df = df;

% for compatibility with SPM struct and cluster analysis
% ------------------------------------------------------
CLU.voxSize = diag(CLU.V.mat)';
CLU.voxSize = CLU.voxSize(1:3);
CLU.VOX = CLU.voxSize;		% compatible with VOL structure
CLU.M = CLU.V(1).mat;

CLU.u = CLU.crit_t;
CLU.k = CLU.cl_size;

if isfield(CLU.V,'fname')
    [a,b,c] = fileparts(CLU.V.fname);
    CLU.title = [b];
else
    CLU.title = 'Analyze image file.';
end

if resl, imP = spm_get(Inf,'*img','Select imgs for data extraction.');
else, % input imP              % to get P
end

%if ~exist(P(1,:))
%    warning('Image files cannot be found in ind. subject directories: Path has changed?')
    % tries to re-find P - only works for my specific directory
%    for i = 1:size(P,1)
%        [d,f,e] = fileparts(P(1,:));
%        [dummy,d] = fileparts(d);
%        newP{i,1}=fullfile('..',d,[f e]);
%    end
%    P = cell2mat(newP);
%    disp(['Found: ' P(1,:) ' etc.'])
%end

if isempty(CLU),
    warning('EMPTY ... mask2clusters found no eligible voxels.')
    clusters = [];
    return
end

if isempty(CLU.XYZ) | isempty(CLU.Z)
    warning('EMPTY ... mask2clusters found no eligible voxels.')
    clusters = [];
    return
end

[clusters] = tor_extract_rois(imP,CLU,CLU,1);
for i = 1:length(clusters),
    clusters(i).P = P;,
    clusters(i).imP = imP;,
    if size(imP,1) == 1 & df == 0,
        %disp('Saving values in mask file in clusters.Z')
        clusters(i).Z = clusters(i).all_data;
    end

    % for SPM2 compatibility
    clusters(i).voxSize = abs(clusters(i).voxSize);
end

if nargout > 2, [subclusters] = cluster_princomp(clusters); end

return

author: canlab

repo: CanlabCore


text:
function regression_table(y,X,names,varargin)
%regression_table(y,X,names,[do figures],[do robust])
%
% performs regression using the regress command, and prints an output table
% X should not have an intercept
%
% Tor Wager

dofigs = 1; dorobust = 0;
if length(varargin) > 0, dofigs = varargin{1}; end
if length(varargin) > 1, dorobust = varargin{2}; end

[n, k] = size(X);

if dorobust
    [b,stat]=robustfit(X,y);
    BINT = NaN * zeros(n + 1,2);
    STATS = [NaN NaN NaN];
    partialr = NaN * zeros(1, k + 1);

else
    [b,dev,stat]=glmfit(X,y);
    [B,BINT,R,RINT,STATS] = regress(y,[ones(n ,1) X]);

    % partial correlations
    partialr(1) = NaN;  % intercept
    for i = 1:k, [x,y,partialr(i+1),partialp] = partialcor(X,y,i); end
end

names = [{'Intercept'}, names];

fprintf(1,'Parameter\tb-hat\tt\tp\tConf. interval\t\tPartial Corr\tRob. Partial Corr.\n')

if dofigs
    create_figure('Regression Table output', k, 2);

end

for i = 1:length(b)

    if i == 1, %intercept
            fprintf(1,'%s\t%3.2f\t%3.2f\t%3.4f\t%3.2f\t%3.2f\t\n', ...
        names{i},b(i),stat.t(i),stat.p(i),BINT(i,1),BINT(i,2));

    else
        if dofigs
            subplot(k, 2, 2 * (i - 2) + 1); [r,str,sig] = prplot(y,X,i-1,0); xlabel([names{i}])
            subplot(k, 2, 2 * (i - 2) + 2); [r2,str,sig] = prplot(y,X,i-1,1); xlabel([names{i} ': Robust IRLS'])
        else
            r = partialr(i);
            r2 = NaN;
        end

        fprintf(1,'%s\t%3.2f\t%3.2f\t%3.4f\t%3.2f\t%3.2f\t%3.2f\t%3.2f\t\n', ...
        names{i},b(i),stat.t(i),stat.p(i),BINT(i,1),BINT(i,2),r,r2);
    end
end

fprintf(1,'\nOverall R2 = %3.2f, F = %3.2f, p = %3.4f\t\n',STATS(1),STATS(2),STATS(3));

return

save_path: github-repos/MATLAB/canlab-CanlabCore


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

% Stationary partially periodic covariance function for a stationary
% covariance function k0 such as covMaternard, covPPard, covRQard and covSEard.
% Partially periodic means that the covariance function is periodic only
% along some dimensions (as opposed to all dimensions).
% Stationary means that the covariance function k0(x,z) depends on the
% data points x,z only through the squared distance
% dxz = (x-z)'*inv(P)*(x-z) where the P matrix is diagonal with ARD parameters
% ell_1^2,...,ell_D^2, where D is the dimension of the input space.
% The covariance function is parameterized as:
%
% k(x,z) = k0(u(x),u(z)),
% [u(x)]([i,i+D]) = [sin(pi*x/p); cos(pi*x/p)]      if p_i is finite
% [u(x)]([i,i+D]) = [x; 0]                          if p_i is not finite
%
% where the period p belongs to covPPERard and hyp0 belong to k0:
%
% hyp = [ log(p_1)
%         log(p_2)
%          .
%         log(p_D)
%         hyp0 ]
%
% The first D hyperparameters of k0 are the log lengthscales such that
% hyp0(i) = log(ell(i)) for i=1..D.
% Note that for k0 = covSEard and D = 1, a faster alternative is covPeriodic.
%
% Copyright (c) by Luigi Acerbi, 2015-12-20.
% Original code by Hannes Nickisch.
%
% See also COVFUNCTIONS, COVPERARD.

nocov = false;                   % default case when no cov argument is provided
if nargin==0, cov = {@covSEard}; nocov = true; end                % default case
if isnumeric(cov)       % detect old version where the cov parameter was missing
  % i <- z, z <- x, x <- hyp, hyp <- cov
  if nargin>3, i = z; end
  if nargin>2, z = x; end
  if nargin>1, x = hyp; end
  hyp = cov; cov = {@covSEard}; nocov = true;
end

if nocov && nargin<2 || ~nocov && nargin<3         % report number of parameters
  K = ['(D+',feval(cov{:}),')']; return
end
if nocov && nargin<3 || ~nocov && nargin<4, z = []; end    % make sure, z exists
xeqz = isempty(z); dg = strcmp(z,'diag');                       % determine mode

[n,D] = size(x);
p = exp(hyp(1:D));
lell = hyp(D+(1:D));
hyp0 =  [lell; lell; hyp(2*D+1:end)];

if nocov && nargin<4 || ~nocov && nargin<5
  [x,z] = u(x,z,p,dg);                      % apply the embedding u:IR^D->IR^2*D
  K = feval(cov{:},hyp0,x,z);
else
  if i<=D
    if isfinite(p(i))                             % periodic dimension
        if dg                                                   % compute distance d
          di = zeros([n,1]);
        else
          if xeqz                                             % symmetric matrix Kxx
            di = repmat(reshape(x(:,i),n, 1),[1, n])...
                -repmat(reshape(x(:,i),1, n),[n, 1]);
          else                                               % cross covariances Kxz
            nz = size(z,1);
            di = repmat(reshape(x(:,i),n, 1),[1,nz])...
                -repmat(reshape(z(:,i),1,nz),[n, 1]);
          end
        end
        di = 2*pi*di/p(i); dD2_dlpi = -2*sin(di).*di;   % derivative dD2i/dlog(p(i))
        [x,z] = u(x,z,p,dg);                    % apply the embedding u:IR^D->IR^2*D
        if dg                                            % compute squared distances
          D2 = zeros(n,1);
        else
          if xeqz
            D2 = sq_dist(x(:,[i,i+D])');
          else
            D2 = sq_dist(x(:,[i,i+D])',z(:,[i,i+D])');
          end
        end
        % reconstruct derivative w.r.t. D2i from derivative w.r.t. log(ell(i))
        dK_dD2 = feval(cov{:},hyp0,x,z,i) + feval(cov{:},hyp0,x,z,i+D);
        dK_dD2 = dK_dD2./(-2*D2); dK_dD2(D2<1e-12) = 0;
        K = dK_dD2.*dD2_dlpi;                                     % apply chain rule
    else                                            % non-periodic dimension
      K = zeros(n);
    end
  else
    [x,z] = u(x,z,p,dg);                    % apply the embedding u:IR^D->IR^2*D
    if i<=2*D
        if isfinite(p(i-D))
          K = feval(cov{:},hyp0,x,z,i-D) + feval(cov{:},hyp0,x,z,i);
        else
          K = feval(cov{:},hyp0,x,z,i-D);
        end
    else
      K = feval(cov{:},hyp0,x,z,i);
    end
  end
end

function [x,z] = u(x,z,p,dg)                % apply the embedding u:IR^D->IR^2*D
  [~,D] = size(x);
  pd = find(isfinite(p));     % periodic dimensions
  x(:,pd) = 2*pi*x(:,pd)*diag(1./p(pd));
  x = [x, zeros(size(x))];
  x(:,pd+D) = cos(x(:,pd)); x(:,pd) = sin(x(:,pd));
  if numel(z)>0 && ~dg
    z(:,pd) = 2*pi*z(:,pd)*diag(1./p(pd));
    z = [z, zeros(size(z))];
    z(:,pd+D) = cos(z(:,pd)); z(:,pd) = sin(z(:,pd));
  end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
function opts = getReweightedOptimizationParams

    % Maximum number of majorization-minimization (MM) iterations used for
    % iteratively re-weighted minimization (should be between 5 and 15).
    opts.maxMMIter = 10;

    % Maximum number of scaled conjugate gradients (SCG) iterations in the
    % inner optimization loop used for image reconstruction (should be
    % between 3 and 10).
    opts.maxSCGIter = 5;

    % Maximum number of cross validation (CV) iterations for hyperparameter
    % selection (should be between 10 and 30).
    opts.maxCVIter = 20;

    % Fraction of the observations used for the training stage in
    % hyperparameter selection (should be between 0.90 and 0.95).
    opts.fractionCVTrainingObservations = 0.95;

    % Initial search range (lower and upper bound) for the regularization
    % hyperparameter on a logarithmic scale (should be between 10-12 and
    % 10^0).
    opts.hyperparameterCVSearchRange = [-12 0];

    % Termination tolerance for max(x(t) - x(t-1)) for the estimate of the
    % high-resolution image over the MM and SCG iterations (should be
    % approx. 1e-3).
    opts.terminationTol = 1e-3;

    % Number of levels used for coarse-to-fine optimization. If empty, the
    % number of levels is set automatically (should be between 2 and 5).
    opts.numCoarseToFineLevels = [];

    % The weighting function that is used to determine the observation
    % confidence weights. This function needs to be defined as a struct P
    % with the following fields:
    %   - P.function is a function handle to weighting function to compute
    %     the confidence weights based on the residual error. The residual
    %     errors for the different frames need to be stored as a cell
    %     array.
    %   - P.parameters are additional parameters that are passed to the
    %     function P.function.
    opts.observationWeightingFunction = estimateHuberObservationWeights;

    % The weighting function that is used to determine the adpative prior
    % weights. This function needs to be defined as a struct P with the
    % following fields:
    %   - P.function is a function handle to weighting function to compute
    %     the prior weights based on a given image in a certain transform
    %     domain.
    %   - P.parameters are additional parameters that are passed to the
    %     function P.function.
    opts.priorWeightingFunction = estimateWBTVPriorWeights;

author: thomas-koehler

repo: SupER


text:
function setup_convnet_burden()
%SETUP_CONVNET_BURDEN Sets up convnet-burden, by adding its folders
% to the Matlab path
%
% Copyright (C) 2017 Samuel Albanie
% Licensed under The MIT License [see LICENSE.md for details]

  check_dependency('mcnExtraLayers') ;
  root = fileparts(mfilename('fullpath')) ;
  addpath(root, [root '/matlab'], [root '/core']) ;

% -----------------------------------
function check_dependency(moduleName)
% -----------------------------------

  name2path = @(name) strrep(name, '-', '_') ;
  setupFunc = ['setup_', name2path(moduleName)] ;
  if exist(setupFunc, 'file')
    vl_contrib('setup', moduleName) ;
  else
    % try adding the module to the path, supressing the warning
    warning('off', 'MATLAB:dispatcher:pathWarning') ;
    addpath(fullfile(vl_rootnn, 'contrib', moduleName)) ;
    warning('on', 'MATLAB:dispatcher:pathWarning') ;

    if exist(setupFunc, 'file')
      vl_contrib('setup', moduleName) ;
    else
      waiting = true ;
      msg = ['module %s was not found on the MATLAB path. Would you like ' ...
             'to install it now? (y/n)\n'] ;
      prompt = sprintf(msg, moduleName) ;
      while waiting
        str = input(prompt,'s') ;
        switch str
          case 'y'
            vl_contrib('install', moduleName) ;
            vl_contrib('compile', moduleName) ;
            vl_contrib('setup', moduleName) ;
            return ;
          case 'n'
            throw(exception) ;
          otherwise
            fprintf('input %s not recognised, please use `y` or `n`\n', str) ;
        end
      end
    end
  end

author: albanie

repo: convnet-burden


text:
function [C_pos,C_neg,Ctot_pos,Ctot_neg] = clustering_coef_wu_sign(W,coef_type)
%CLUSTERING_COEF_WU_SIGN     Multiple generalizations of the clustering coefficient
%
%   [C_pos,C_neg,Ctot_pos,Ctot_neg] = clustering_coef_wu_sign(W,coef_type);
%
%   The weighted clustering coefficient is the average weight or intensity
%   of all triangles associated with each node.
%
%   Inputs:
%       W,
%           Weighted undirected connection matrix
%
%       corr_type,
%           Desired type of clustering coefficient.
%           Options:
%           1,  (default) Onnela et al. formula, used in original
%               clustering_coef_wu.m. Computed separately for positive &
%               negative weights.
%           2,  Zhang & Horvath formula, similar to Onnela formula except
%               denominator of Onnela formula relies on binarizing the
%               network whereas this denominator is based on weight value,
%               which reduces the sensitivity of this measure to the
%               weights directly connected to the node of interest.
%               Computed separately for positive & negative weights.
%           3,  Constantini & Perugini's generalization of the Zhang &
%               Horvath formula. This formula takes both positive &
%               negative weights into account simultaneously, & is
%               particularly sensitive to non-redundancy in path
%               information based on sign (i.e., when two weights are
%               positive & one negative, or all three are negative, both of
%               which indicate that the weight of the third path is not
%               redundant information). Produces only one value.
%
%
%   Outputs:
%       C_pos/C_neg,
%           Clustering coefficient vector for positive/negative weights.
%           For the third option, only one vector is outputted (as C_pos).
%       Ctot_pos/Ctot_neg,
%           Mean clustering coefficient for positive and negative weights.
%
%   References:
%       Onnela et al. (2005) Phys Rev E 71:065103
%       Zhang & Horvath (2005) Stat Appl Genet Mol Biol 41:1544-6115
%       Costantini & Perugini (2014) PLOS ONE 9:e88669
%
%
%   Contributor: Jeff Spielberg, Boston University, 2014-2015
%                (script based on clustering_coef_wu.m)

%
%   Modification History:
%   May 2014: Added computation of pos & neg weights separately &
%             computation of mean coefficient (Jeff Spielberg)
%   May 2015: Added computation of Zhang & Horvath and Constantini &
%             Perugini formulas (Jeff Spielberg)
%   May 2016: Bugfix in computation of the denominator of the Costantini &
%             Perugini (flag 3) version (Chiara Pintossi)

if ~exist('coef_type','var')
    coef_type = 1;
end

n            = length(W);                   %number of nodes
W(1:n+1:end) = 0;

switch coef_type
    case 1
        W_pos                = W.*(W>0);
        K_pos                = sum(W_pos~=0,2);
        cyc3_pos             = diag((W_pos.^(1/3))^3);
        K_pos(cyc3_pos == 0) = inf;             %if no 3-cycles exist, make C=0 (via K=inf)
        C_pos                = cyc3_pos./(K_pos.*(K_pos-1));         %clustering coefficient
        Ctot_pos             = mean(C_pos);

        W_neg                = -W.*(W<0);
        K_neg                = sum(W_neg~=0,2);
        cyc3_neg             = diag((W_neg.^(1/3))^3);
        K_neg(cyc3_neg == 0) = inf;             %if no 3-cycles exist, make C=0 (via K=inf)
        C_neg                = cyc3_neg./(K_neg.*(K_neg-1));         %clustering coefficient
        Ctot_neg             = mean(C_neg);
    case 2
        W_pos    = W.*(W>0);
        cyc3_pos = zeros(n,1);
        cyc2_pos = zeros(n,1);
        for i = 1:n
            for j = 1:n
                for q = 1:n
                    cyc3_pos(i) = cyc3_pos(i)+(W_pos(j,i)*W_pos(i,q)*W_pos(j,q));
                    if j~=q
                        cyc2_pos(i) = cyc2_pos(i)+(W_pos(j,i)*W_pos(i,q));
                    end
                end
            end
        end
        cyc2_pos(cyc3_pos == 0) = inf;             %if no 3-cycles exist, make C=0 (via K=inf)
        C_pos                   = cyc3_pos./cyc2_pos;         %clustering coefficient
        Ctot_pos                = mean(C_pos);

        W_neg    = -W.*(W<0);
        cyc3_neg = zeros(n,1);
        cyc2_neg = zeros(n,1);
        for i = 1:n
            for j = 1:n
                for q = 1:n
                    cyc3_neg(i) = cyc3_neg(i)+(W_neg(j,i)*W_neg(i,q)*W_neg(j,q));
                    if j~=q
                        cyc2_neg(i) = cyc2_neg(i)+(W_neg(j,i)*W_neg(i,q));
                    end
                end
            end
        end
        cyc2_neg(cyc3_neg == 0) = inf;             %if no 3-cycles exist, make C=0 (via K=inf)
        C_neg                   = cyc3_neg./cyc2_neg;         %clustering coefficient
        Ctot_neg                = mean(C_neg);
    case 3
        cyc3         = zeros(n,1);
        cyc2         = zeros(n,1);

        for i = 1:n
            for j = 1:n
                for q = 1:n
                    cyc3(i) = cyc3(i)+(W(j,i)*W(i,q)*W(j,q));
                    if j~=q
                        cyc2(i) = cyc2(i)+abs(W(j,i)*W(i,q));
                    end
                end
            end
        end

        cyc2(cyc3 == 0) = inf;             %if no 3-cycles exist, make C=0 (via K=inf)
        C_pos           = cyc3./cyc2;         %clustering coefficient
        Ctot_pos        = mean(C_pos);
        C_neg           = nan(size(C_pos));
        Ctot_neg        = nan(size(Ctot_pos));
end


text:
classdef (Abstract) mme_buslink_pf_ac < mp.mme_buslink

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

%     properties
%         name = 'buslink';
%     end

    methods
        function obj = add_vars(obj, mm, nm, dm, mpopt)
            nme = obj.network_model_element(nm);

            mm.init_indexed_name('var', 'Plink', {3});
            mm.init_indexed_name('var', 'Qlink', {3});
            mmx_i1 = mm.var.N + 1;
            for p = 1:3
                mm.add_var('Plink', {p}, nme.nk, 0, -Inf, Inf);
            end
            mmx_iN = mm.var.N;
            mm.aux_data.var_map{end+1} = ...
                {'zr', nm.zr.idx.i1.Plink(1), nm.zr.idx.iN.Plink(end), [], mmx_i1, mmx_iN, []};

            mmx_i1 = mm.var.N + 1;
            for p = 1:3
                mm.add_var('Qlink', {p}, nme.nk, 0, -Inf, Inf);
            end
            mmx_iN = mm.var.N;
            mm.aux_data.var_map{end+1} = ...
                {'zi', nm.zi.idx.i1.Qlink(1), nm.zi.idx.iN.Qlink(end), [], mmx_i1, mmx_iN, []};
        end

        function obj = add_constraints(obj, mm, nm, dm, mpopt)
            nme = obj.network_model_element(nm);

            %% add Qlink constraints for PV node on 3-phase side

            %% indices of buslinks connected to PV nodes via port 2 (3&4)
            ad = mm.aux_data;
            [~, k, ~] = find(nme.C(ad.pv, nme.nk+1:2*nme.nk));
            n = length(k);
            if n
                I = sparse(1:n, k, 1, n, nme.nk);
                zz = sparse(n, nme.nk);
                A = [I -I zz; I zz -I];
                b = zeros(2*n, 1);
                vs = struct('name', {'Qlink', 'Qlink', 'Qlink'}, ...
                            'idx', {{1}, {2}, {3}} );
                mm.add_lin_constraint('buslink_qlink', A, b, b, vs);
            end
        end

        function [A_va_pq, A_va_pv, b_va, A_vm, b_vm] = voltage_constraints(obj, nme, ad)
            %% form constraint matrices for matching
            %%  voltage angles for pv and pq nodes
            %%  voltage magnitudes for pq nodes
            %% columns of A_va_pv correspond to ...
            %%  'Va_pv', 'Va3_pv', 'Va3_pv', 'Va3_pv',
            %% columns of A_va_pq correspond to ...
            %%  'Va_pq', 'Va3_pq', 'Va3_pq', 'Va3_pq'
            %% columns of A_vm correspond to ...
            %%  'Vm_pq', 'Vm3_pq', 'Vm3_pq', 'Vm3_pq'
            nk = nme.nk;

            %% basic constraint matrix for voltage equality (all nodes)
            [A, b_va, b_vm, Istack] = nme.voltage_constraints();

            %% sub-matrices by node type
            A_ref = A(:, ad.ref);
            A_pv = A(:, ad.pv);
            A_pq = A(:, ad.pq);

            %% voltage angle constraints (all buslinks)
            A_va_pq = A_pq;
            A_va_pv = A_pv;

            %% voltage magnitude constraints (all buslinks)
            A_vm = A_pq;

            %% indices of buslinks connected to REF nodes (fixed va)
            [~, k_va1, ~] = find(nme.C(ad.ref, 1:nk));      %% via port 1
            [~, k_va2, ~] = find(nme.C(ad.ref, nk+1:2*nk)); %% via port 2 (3&4)

            %% adjust RHS of va constraints involving fixed voltage angles
            if ~isempty(k_va1) || ~isempty(k_va2)
                [va, vm] = nme.aux_data_va_vm(ad);
                va1 = nme.C(:, 1:nk)' * va;         %% port 1 voltage angles
                va2 = nme.C(:, nk+1:2*nk)' * va;    %% port 2 voltage angles
                va_adj = zeros(nk, 1);
                va_adj(k_va1) = -va1(k_va1);
                va_adj(k_va2) =  va2(k_va2);
                b_va = b_va + Istack * va_adj;

                %% delete redundate (and currently incorrect) constraints
                %% corresponding to k_va2 and ports 3 and 4
                if ~isempty(k_va2)
                    d_va = [nk+k_va2; 2*nk+k_va2];
                    A_va_pq(d_va, :) = [];
                    A_va_pv(d_va, :) = [];
                    b_va(d_va) = [];
                end
            end

            %% indices of buslinks connected to REF/PV nodes (fixed vm)
            rpv = [ad.ref;ad.pv];   %% indices of REF/PV nodes
            [~, k_vm1, ~] = find(nme.C(rpv, 1:nk));         %% via port 1
            [~, k_vm2, ~] = find(nme.C(rpv, nk+1:2*nk));    %% via port 2 (3&4)

            %% adjust RHS of vm constraints involving fixed voltage magnitudes
            if ~isempty(k_vm1) || ~isempty(k_vm2)
                [va, vm] = nme.aux_data_va_vm(ad);
                vm1 = nme.C(:, 1:nk)' * vm;         %% port 1 voltage magnitudes
                vm2 = nme.C(:, nk+1:2*nk)' * vm;    %% port 2 voltage magnitudes
                vm_adj = zeros(nk, 1);
                vm_adj(k_vm1) = -vm1(k_vm1);
                vm_adj(k_vm2) =  vm2(k_vm2);
                b_vm = b_vm + Istack * vm_adj;

                %% delete redundate (and currently incorrect) constraints
                %% corresponding to k_vm2 and ports 3 and 4
                if ~isempty(k_vm2)
                    d_vm = [nk+k_vm2; 2*nk+k_vm2];
                    A_vm(d_vm, :) = [];
                    b_vm(d_vm) = [];
                end
            end
        end

        function obj = data_model_update(obj, mm, nm, dm, mpopt)
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
% This is the kernel of FAIR
%
% the kernel contains
%   contents.m	this file
%
%   data                data and examples included in the toolbox
%   distances           various distance measures (SSD etc.) and their administration
%   examples            numerous examples (organized by chapters of the FAIR book)
%   imgModel            various image models (interpolation etc.) and their administration
%   landmarks           tools for landmark registration
%   matrixfree	        tools for matrixfree operations
%   numerics            numerical tools including optimization
%   regularizers        various regularizer (elastic etc.) and their administration
%   tools               numerous tools
%   transformations     various transformation models (affine etc.) and their administration
%==============================================================================
function debit = contents
if nargout == 0, help(mfilename); return; end;

debit = {
  'contents.m'
};
%==============================================================================


text:
function []=coreg_pos(firstL,lastL,firstP,lastP,transL,transP,slc_osf)
% COREG_POS create a position file for fine coregistration
%    []=coreg_pos(firstL,lastL,firstP,lastP,transL,transP)
%
%   Andy Hooper, Aug 2005
%
% ======================================================================
% 07/2006 AH:  changed spacing of windows to give same number (3600)
%              independent of image size
% 11/2009 AH:  code to use translations put back in
% 03/2010 AH:  revert to just 3600 windows
% 09/2009 MA:  minor update for oversampled data + translation_lp
% 10/2010 JCM: drop negative values at the edges of the scene
% ======================================================================
%

if nargin<7
 slc_osf=1 ;              % default oversampling factor is 1
end
slc_osf
winSize=64*slc_osf ;      % correlation window size default: 64 should match with coreg.dorisin
                          %                                : 128 for oversampled data

nWinx=30    % # of windows in range
nWiny=120   % # of windows in azimuth

x=round(linspace(firstP+winSize+transP,lastP-winSize+transP,nWinx));
y=round(linspace(firstL+winSize+transL,lastL-winSize+transL,nWiny));

% drop negative values at the edges of the scene JCM
keep_ix=x>0;
keep_iy=y>0;
x=x(keep_ix);
y=y(keep_iy);

[Y,X]=meshgrid(y,x);

fid=fopen('fc_pos.in','w');
fprintf(fid,'%i  %i\n',[Y(:),X(:)]');
fclose(fid);

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
function pi = CPD_to_pi(CPD, msg_type, n, ps, msg, evidence)
% CPD_TO_PI Compute the pi vector (gaussian)
% function pi = CPD_to_pi(CPD, msg_type, n, ps, msg, evidence)

switch msg_type
 case 'd',
  error('gaussian_CPD can''t create discrete msgs')
 case 'g',
  dps = ps(CPD.dps);
  k = evidence{dps};
  if isempty(k)
    error('gmux node must have observed discrete parent')
  end
  m = msg{n}.pi_from_parent{k};
  B = CPD.weights(:,:,k);
  pi.mu = CPD.mean(:,k) + B * m.mu;
  pi.Sigma = CPD.cov(:,:,k) + B * m.Sigma * B';
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function val = get(c,prop_name)

   % GET - Get correlation properties
   %
   % VAL = get(C,PROP_NAME)
   % See HELP CORRELATION for description of primary property names:
   %   Waveforms, Trig, Corr, Lag, Stat, Link, Clust
   %
   % Additional scalar properties:
   %   Traces:         number of traces
   %   Data_Length:    number of samples in each trace
   %   Fs:             Frequency
   %   Period:         Period
   %   Nyq:            Nyquist frequency
   %   **NOTE - The scalar properties listed above exist for each waveform
   %   (except Traces). However, the correlation object requires them to be
   %   the same. GET only reads the value associated with the first waveform
   %   Normally this will not be an issue. If the correlation object
   %   was created manually it is possible for these fields to be different.
   %
   % Additional vector properties (nx1):
   %   Start, Trig, End:             same as Start_Matlab, End_Matlab
   %   Start_Str, Trig_Str, End_Str:           String times
   %   Start_Matlab, Trig_Matlab, End_Matlab:  Matlab-format times
   %   Start_Epoch, Trig_Epoch, End_Epoch:     Epoch-format times
   %
   % Additional vector properties (nx1 cell vector)
   %   Sta (or Station):       Cell vector of station names
   %   Chan (or Component):    Cell vector of channel names
   %
   % Additional matrix properties (nxm matrix)
   %   Data:                    Matrix of raw trace data
   %                           (n traces) x (m samples)

   % Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   warning('using NewCorrelation/get, alternate ways exist of accessing data.')

   narginchk(2,inf)

   % got rid of test for c

   switch upper(prop_name)
      case {'WAVEFORMS', 'WAVEFORM', 'WAVES'}
         val = c.W;
      case {'TRIG'}
         val = c.trig;
      case {'TRIG_STR'}
         val = cellstr(datestr(c.trig,'dd-mmm-yyyy HH:MM:SS.FFF'));
      case {'TRIG_MATLAB'}
         val = c.trig;
      case {'CORR'}
         val = c.corrmatrix;
      case {'LAG'}
         val = c.lags;
      case {'STAT'}
         val = c.stat;
      case {'LINK'}
         val = c.link;
      case {'CLUST'}
         val = c.clust;
      case {'TRACES'}
         val = c.ntraces;

         % FROM WAVEFORM/GET (SCALAR OUTPUT)
      case {'DATA_LENGTH'}
         val = c.traces(1).nsamples();
      case {'FS'}
         val = c.traces(1).samplerate;

         %     % FROM WAVEFORM/GET (VECTOR)
      case {'STATION', 'STA'}
         val = c.stations;
      case {'COMPONENT', 'CHAN'}
         val = c.channels;

         % OTHER ROUTINES
      case {'DATA'}
         val = double(c.traces);

      otherwise
         try
            % see if it is a valid waveform property
            warning('accessing a waveform property %s. ', prop_name);
            val = get(c.W,prop_name);
         catch er
            error('''%s'' is not a valid argument for correlation/get\n%s\n', prop_name, er.message);
         end
   end;
end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CHECKVALS
% This function checks to see if each waveform has the same parameters

function checkvals(vals)

   same = all(vals(:) == vals(1));

   if ~same
      % Warning is disabled because it gets called repeatedly (annoyingly)
      % from within other scripts
      % warning('Waveforms have different frequencies or different numbers of samples. Consider VERIFY function');
   end;
end

author: geoscience-community-codes

repo: GISMO


text:
function cost = leaveoneout(model, estfct,combinefct)

% Estimate the performance of a trained model with leave-one-out crossvalidation
%
% CAUTION!! Use this function only to obtain the value of the leave-one-out score
% function given the tuning parameters. Do not use this function together with
% 'tunelssvm', but use 'leaveoneoutlssvm' instead. The latter is a faster
% implementation which uses previously computed results.
%
% >> leaveoneout({X,Y,type,gam,sig2})
% >> leaveoneout(model)
%
% In each iteration, one leaves one point, and fits a model on the
% other data points. The performance of the model is estimated
% based on the point left out. This procedure is repeated for each
% data point. Finally, all the different estimates of the
% performance are combined (default by computing the mean). The
% assumption is made that the input data is distributed independent
% and identically over the input space.
%
%
% Full syntax
%
%     1. Using the functional interface for the LS-SVMs:
%
% >> cost = leaveoneout({X,Y,type,gam,sig2,kernel,preprocess})
% >> cost = leaveoneout({X,Y,type,gam,sig2,kernel,preprocess}, estfct)
% >> cost = leaveoneout({X,Y,type,gam,sig2,kernel,preprocess}, estfct, combinefct)
%
%       Outputs
%         cost          : Cost estimated by leave-one-out crossvalidation
%
%       Inputs
%         X             : Training input data used for defining the LS-SVM and the preprocessing
%         Y             : Training output data used for defining the LS-SVM and the preprocessing
%         type          : 'function estimation' ('f') or 'classifier' ('c')
%         gam           : Regularization parameter
%         sig2          : Kernel parameter (bandwidth in the case of the 'RBF_kernel')
%         kernel(*)     : Kernel type (by default 'RBF_kernel')
%         preprocess(*) : 'preprocess'(*) or 'original'
%         estfct(*)     : Function estimating the cost based on the residuals (by default mse)
%         combinefct(*) : Function combining the estimated costs on the different folds (by default mean)
%
%
%     2. Using the object oriented interface for the LS-SVMs:
%
% >> cost = leaveoneout(model)
% >> cost = leaveoneout(model, estfct)
% >> cost = leaveoneout(model, estfct, combinefct)
%
%       Outputs
%         cost          : Cost estimated by leave-one-out crossvalidation
%
%       Inputs
%         model         : Object oriented representation of the model
%         estfct(*)     : Function estimating the cost based on the residuals (by default mse)
%         combinefct(*) : Function combining the estimated costs on the different folds (by default mean)
%
%
%
% See also:
%    crossvalidate, trainlssvm, simlssvm

% Copyright (c) 2011,  KULeuven-ESAT-SCD, License & help @ http://www.esat.kuleuven.be/sista/lssvmlab

% LS-SVMlab
eval('model = initlssvm(model{:});',' ');
eval('estfct;','estfct=''mse'';');
eval('combinefct;','combinefct=''mean'';');

%
%initialize: no incremental  memory allocation
%
p = randperm(model.nb_data);
px = model.xtrain(p,:);
py = model.ytrain(p,:);
[~,Y] = postlssvm(model,[],py); % Y is raw data, non preprocessed

% kernel matrix computation
K = kernel_matrix(px,model.kernel_type,model.kernel_pars);

Ka = pinv([K+eye(model.nb_data)./model.gam ones(model.nb_data,1);ones(1,model.nb_data) 0]);
sol = Ka*[py;0]; model.alpha = sol(1:end-1); model.b = sol(end);
yh = py - model.alpha./diag(Ka(1:model.nb_data,1:model.nb_data));

[~,yh] = postlssvm(model,[],yh);
if ~(model.type(1)=='c')
    cost = feval(estfct,yh-Y);
else
    cost = feval(estfct,Y,sign(yh));
end

author: peterhcharlton

repo: RRest


text:
function varargout = tanh(varargin)
%TANH   Hyperbolic tangent of a DISKFUN.
%   TANH(F) returns the hyperbolic tangent of a DISKFUN object F.
%
% See also DISKFUN/TAN.

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

[varargout{1:nargout}] = tanh@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function fx2 = p12_fx2 ( x )

%*****************************************************************************80
%
%% P12_FX2 evaluates the second derivative of the function for problem 12.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the abscissa.
%
%    Output, real FX2, the second derivative of the function at X.
%
  factor = 1000.0;

  if ( x == 1.0 )
    fx2 = 0.0;
  else
    y = x - 1.0;
    fx2 = - 2.0 * factor * exp ( - 1.0 / y.^2 ) * ( y.^2 - 2.0 ) / y.^5;
  end

author: johannesgerer

repo: jburkardt-m


text:
%% This function is used to turn off the server on the robot
% function [ output_args ] = net_turnOffServer( t )
% t: is the TCP/IP connection object
% Copy right, Mohammad SAFEEA, 3rd of May 2017

repo: KST-Kuka-Sunrise-Toolbox


text:
function [h, compUpJV] =  lfmComputeH3JV(gamma1_p, gamma1_m, sigma2, t1, ...
    t2, preFactor, mode)

% LFMCOMPUTEH3JV Helper function for computing part of the LFMJV kernel.
% FORMAT
% DESC computes a portion of the LFMJV kernel.
% ARG gamma1 : Gamma value for first system.
% ARG gamma2 : Gamma value for second system.
% ARG sigma2 : length scale of latent process.
% ARG t1 : first time input (number of time points x 1).
% ARG t2 : second time input (number of time points x 1).
% ARG preFactor : precomputed constants.
% ARG mode: indicates the correct derivative.
% RETURN h : result of this subcomponent of the kernel for the given values.
%
% COPYRIGHT : Mauricio Alvarez, 2010

% KERN

% Evaluation of h

if nargout>1
    compUpJV{1} = lfmjvComputeUpsilonMatrix(gamma1_p,sigma2, t1,t2, mode);
    compUpJV{2} = lfmjvComputeUpsilonMatrix(gamma1_m,sigma2, t1,t2, mode);
    h = preFactor(1)*compUpJV{1} + preFactor(2)*compUpJV{2};
else
    h = preFactor(1)*lfmjvComputeUpsilonMatrix(gamma1_p,sigma2, t1,t2, mode) ...
        + preFactor(2)*lfmjvComputeUpsilonMatrix(gamma1_m,sigma2, t1,t2, mode);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [posterior,out] = VBA_wrapup(posterior,options,dim,suffStat,u,y,it,fromPause)
% wraps up the VB inversion
% function [posterior,out] = VBA_wrapup(posterior,options,dim,suffStat,u,y,it)
% IN:
%   - [VBA internal variables; see VBA_NLStateSpaceModel.m]
% OUT:
%   - posterior: the final 'posterior' structure
%   - out: the final 'out' structure

try; fromPause; catch; fromPause = 0; end

if options.DisplayWin
    display = options.display;
end

%------- Summary output structure -----%
% Restore posterior from ODE limit if required
[posterior,options,dim,suffStat] = VBA_odeLim2NLSS(posterior,options,dim,suffStat,u,1);
% Store options, sufficient statistics, etc...
out.F = suffStat.F(end);
out.options = options;
out.u = VBA_getU(u,options,dim,'back2micro');
out.y = y;
out.dim = dim;
out.it = it;
out.suffStat = suffStat;
out.date = clock;
out.dt = toc(options.tStart);
out.fit = VBA_getFit(posterior,out);
if fromPause
    return
end

% Display inversion status
if ~options.OnLine
    [st,i] = dbstack;
    try,ifInit=isequal(st(i+2).name,'VBA_Initialize');catch,ifInit=0;end
    if ~ifInit % main VB inversion has converged
        status = 'inversion';
        try
            if options.DisplayWin
                % display diagnostics
                figure(display.hfp)
                set(display.hm(1),'string','Deriving diagnostics... ');
                set(display.hm(2),'string','0%');
                out.options.display = display;
                out.diagnostics = VBA_getDiagnostics(posterior,out);
                set(display.hm(2),'string','OK');
                VBA_ReDisplay(posterior,out);
            end
            out.options = rmfield(out.options,'display');
            out.options.hf = display.hfp;
        end
    else
        status = 'initialization';
        if options.DisplayWin
            out.options.display = display;
        end
    end

    % display overall inversion time
    if floor(out.dt./60) == 0
        timeString = [num2str(floor(out.dt)),' sec'];
    else
        timeString = [num2str(floor(out.dt./60)),' min'];
    end
    VBA_disp(['VB ',status,' complete (took ~',timeString,').'],options)
    if ~ifInit
        str = VBA_summary(out,1);
        VBA_disp(str(setdiff(1:6,1:2)),options)
    end
%     VBA_disp(' ',options)
    drawnow
end

author: MBB-team

repo: VBA-toolbox


text:
function [mmps2] = cmps22mmps2(cmps2)
% Convert acceleration from centimeters per square centimeter to
% millimeters per second-squared
% Chad A. Greene 2012
mmps2 = cmps2*1e+1;
end

author: Sable

repo: mcbench-benchmarks


text:
function filter = spm_cfg_eeg_filter
% configuration file for EEG Filtering
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Stefan Kiebel
% $Id: spm_cfg_eeg_filter.m 5377 2013-04-02 17:07:57Z vladimir $

rev = '$Rev: 5377 $';

D = cfg_files;
D.tag = 'D';
D.name = 'File Name';
D.filter = 'mat';
D.num = [1 1];
D.help = {'Select the EEG mat file.'};

type = cfg_menu;
type.tag = 'type';
type.name = 'Type';
type.labels = {'Butterworth', 'FIR'};
type.values = {'butterworth', 'fir'};
type.val = {'butterworth'};
type.help = {'Select the filter typee.'};

band = cfg_menu;
band.tag = 'band';
band.name = 'Band';
band.labels = {'Lowpass', 'Highpass', 'Bandpass', 'Stopband'};
band.values = {'low' 'high' 'bandpass' 'stop'};
band.val = {'low'};
band.help = {'Select the filter band.'};

freq = cfg_entry;
freq.tag = 'freq';
freq.name = 'Cutoff(s)';
freq.strtype = 'r';
freq.num = [1 inf];
freq.help = {'Enter the filter cutoff'};

dir = cfg_menu;
dir.tag = 'dir';
dir.name = 'Direction';
dir.labels = {'Zero phase', 'Forward', 'Backward'};
dir.values = {'twopass', 'onepass', 'onepass-reverse'};
dir.val = {'twopass'};
dir.help = {'Select the filter direction.'};

order = cfg_entry;
order.tag = 'order';
order.name = 'Order';
order.val = {5};
order.strtype = 'n';
order.num = [1 1];
order.help = {'Enter the filter order'};

prefix         = cfg_entry;
prefix.tag     = 'prefix';
prefix.name    = 'Filename Prefix';
prefix.help    = {'Specify the string to be prepended to the filenames of the filtered dataset. Default prefix is ''f''.'};
prefix.strtype = 's';
prefix.num     = [1 Inf];
prefix.val     = {'f'};

filter = cfg_exbranch;
filter.tag = 'filter';
filter.name = 'Filter';
filter.val = {D type band freq dir order, prefix};
filter.help = {'Filters EEG/MEG data.'};
filter.prog = @eeg_filter;
filter.vout = @vout_eeg_filter;
filter.modality = {'EEG'};

function out = eeg_filter(job)
% construct the S struct
S = job;
S.D = S.D{1};

out.D = spm_eeg_filter(S);
out.Dfname = {fullfile(out.D.path, out.D.fname)};

function dep = vout_eeg_filter(job)
% Output is always in field "D", no matter how job is structured
dep = cfg_dep;
dep.sname = 'Filtered Data';
% reference field "D" from output
dep.src_output = substruct('.','D');
% this can be entered into any evaluated input
dep.tgt_spec   = cfg_findspec({{'strtype','e'}});

dep(2) = cfg_dep;
dep(2).sname = 'Filtered Datafile';
% reference field "Dfname" from output
dep(2).src_output = substruct('.','Dfname');
% this can be entered into any file selector
dep(2).tgt_spec   = cfg_findspec({{'filter','mat'}});

author: spm


text:
%% this is an implementation of kalman filter using belief propagation.
% author: Shuang Wang
% email: shw070@ucsd.edu
% Division of Biomedical Informatics, University of California, San Diego.

clear;
addpath(genpath('.'));
%%
% x_k = A x_t-1 + q; N(0, V_Q);
% y_k = B x_k + w; N(0, V_W);
% x_1 ~ N(mu_0, V_0);

%% generating fake data.
s = RandStream('mt19937ar','Seed',4);
RandStream.setGlobalStream(s);
T = 40;
x = zeros(2, T);
y = zeros(2, T);
d_y = size(y, 1);
d_x = size(x, 1);
mu_0 = [4; 4];
V_0 = [1, 0
       0, 1];
V_Q = [0.1,   0
         0, 0.3];
V_W = [0.05,   0
          0, 0.1];
A = [   1, 0.04
     0.02,    1];
B = [1, 0
     0, 1];

 x(:, 1) = mvnrnd(mu_0, V_0, 1)';
 y(:, 1) = B * x(:, 1) + mvnrnd([0, 0], V_W, 1)';
 for t = 2:T
     x(:, t) = A * x(:, t-1) + mvnrnd([0, 0], V_Q, 1)';
     y(:, t) = B * x(:, t) + mvnrnd([0, 0], V_W, 1)';
 end

[msg, varNode_x, factorNode_x, factorNode_y, belief] = initializeMessage(T, mu_0, V_0, d_x, y, B, V_W);

for bp_iter = 1:2*T
    %% variable node update
    for i = 1:T
        if (i < T)
            msg{varNode_x{i}.backwardNeighborMsgID}.toFactorNode ...
                = GaussianMultiply(msg{varNode_x{i}.forwardNeighborMsgID}.toVarNode,...
                msg{varNode_x{i}.lowerNeighborMsgID}.toVarNode);

            msg{varNode_x{i}.forwardNeighborMsgID}.toFactorNode ...
                = GaussianMultiply(msg{varNode_x{i}.backwardNeighborMsgID}.toVarNode,...
                msg{varNode_x{i}.lowerNeighborMsgID}.toVarNode);
        else
            msg{varNode_x{i}.backwardNeighborMsgID}.toFactorNode = msg{varNode_x{i}.lowerNeighborMsgID}.toVarNode;
        end
    end

    %% factor node update
    for i = 2:T
        msg{factorNode_x{i}.forwardNeighborMsgID}.toVarNode.mu ...
            = A * msg{factorNode_x{i}.backwardNeighborMsgID}.toFactorNode.mu;
        if(~isempty(msg{factorNode_x{i}.backwardNeighborMsgID}.toFactorNode.V))
            msg{factorNode_x{i}.forwardNeighborMsgID}.toVarNode.V ...
                = A * msg{factorNode_x{i}.backwardNeighborMsgID}.toFactorNode.V * (A') + V_Q;
            msg{factorNode_x{i}.forwardNeighborMsgID}.toVarNode.iV = [];
        else
            msg{factorNode_x{i}.forwardNeighborMsgID}.toVarNode.V = [];
            tmp_iAtiviA = (A')\msg{factorNode_x{i}.backwardNeighborMsgID}.toFactorNode.iV/A;
            msg{factorNode_x{i}.forwardNeighborMsgID}.toVarNode.iV ...
                = tmp_iAtiviA -(eye(d_x) + tmp_iAtiviA*V_Q)\tmp_iAtiviA*V_Q*tmp_iAtiviA;
        end
        msg{factorNode_x{i}.backwardNeighborMsgID}.toVarNode.mu ...
            = A\msg{factorNode_x{i}.forwardNeighborMsgID}.toFactorNode.mu;
        if(~isempty(msg{factorNode_x{i}.forwardNeighborMsgID}.toFactorNode.V))
            msg{factorNode_x{i}.backwardNeighborMsgID}.toVarNode.iV = [];
            msg{factorNode_x{i}.backwardNeighborMsgID}.toVarNode.V ...
                = A\(msg{factorNode_x{i}.forwardNeighborMsgID}.toFactorNode.V + V_Q)/(A');
        else
            tmp_iv = msg{factorNode_x{i}.forwardNeighborMsgID}.toFactorNode.iV;
            msg{factorNode_x{i}.backwardNeighborMsgID}.toVarNode.V = [];
            msg{factorNode_x{i}.backwardNeighborMsgID}.toVarNode.iV ...
                = (A')*(tmp_iv -(eye(d_x) + tmp_iv*V_Q)\tmp_iv*V_Q*tmp_iv)*A;
        end
    end
end % bp_iter

%% update belief
for i = 1:T
    belief{i} = GaussianMultiply(msg{varNode_x{i}.backwardNeighborMsgID}.toFactorNode,...
                msg{varNode_x{i}.backwardNeighborMsgID}.toVarNode);
    mu_tmp_sm(:,i) = belief{i}.mu;
end
%%
figure(1); clf;
plot(x(1, :), x(2, :), '.'); hold on;
plot(y(1, :), y(2, :), 'ro'); hold on;
plot(mu_tmp_sm(1, :), mu_tmp_sm(2, :), 'g*'); hold on;
h = legend('Hidden state', 'observation', 'kalman smoother');
set(h, 'box', 'off', 'location', 'best');
fprintf('SAD:x_z = %f\n', sum(abs((y(:) - x(:)))));
fprintf('SAD:kf_sm_z = %f\n', sum(abs((mu_tmp_sm(:) - x(:)))));


text:
function cosh_values_test ( )

%*****************************************************************************80
%
%% COSH_VALUES_TEST demonstrates the use of COSH_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 June 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'COSH_VALUES_TEST:\n' );
  fprintf ( 1, '  COSH_VALUES stores values of \n' );
  fprintf ( 1, '  the hyperbolic cosine function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X           FX\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = cosh_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, fx );

  end

  return
end


text:
function msm_to_mm_test07 ( )

%*****************************************************************************80
%
%% MSM_TO_MM_TEST07 tests MSM_TO_MM_ARRAY_INTEGER_SYMMETRIC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 November 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MSM_TO_MM_TEST07\n' );
  fprintf ( 1, '  Convert an MSM to MM array integer symmetric format.\n' );

  output_filename = 'msm_to_mm_test07.mm';

  a = i4mat_indicator ( 4, 4 );
  a = a + a';
%
%  Have MSM_TO_MM write the matrix to a file.
%
  msm_to_mm ( output_filename, a, 'array', 'integer', 'symmetric' );

  return
end


text:
function res = spm_eeg_reduce_cva(S)
% Plugin for data reduction using PCA
% FORMAT res = spm_eeg_reduce_pca(S)
%
% S                     - input structure
% fields of S:
%    S.ncomp            - number of PCA components
%
% Output:
%  res -
%   If no input is provided the plugin returns a cfg branch for itself
%
%   If input is provided:
%      montage struct implementing projection to PCA subspace
%______________________________________________________________________________________
% Copyright (C) 2012 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: spm_eeg_reduce_cva.m 5675 2013-10-09 14:27:17Z vladimir $

if nargin == 0

    cvachan = cfg_branch;
    cvachan.tag = 'cvachan';
    cvachan.name = 'Channels to reduce';
    cvachan.val = {spm_cfg_eeg_channel_selector};

    refchan = cfg_branch;
    refchan.tag = 'refchan';
    refchan.name = 'Reference channels';
    refchan.val = {spm_cfg_eeg_channel_selector};

    ncomp = cfg_entry;
    ncomp.tag = 'ncomp';
    ncomp.name = 'Number of components';
    ncomp.strtype = 'n';
    ncomp.num = [1 1];
    ncomp.val = {1};
    ncomp.help = {'Number of components to retain'};

    outlabel = cfg_entry;
    outlabel.tag = 'outlabel';
    outlabel.name = 'Output channel label';
    outlabel.strtype = 's';
    outlabel.num = [1 Inf];
    outlabel.help = {'Label for the output channel(s).',...
        'Numbers are added for multiple channels.'};

    foi = cfg_entry;
    foi.tag = 'foi';
    foi.name = 'Frequency band of interest';
    foi.strtype = 'r';
    foi.num = [1 2];
    foi.val = {[0 Inf]};
    foi.help = {'Frequency window to optimize for'};

    tshiftwin = cfg_entry;
    tshiftwin.tag = 'tshiftwin';
    tshiftwin.name = 'Time shift window';
    tshiftwin.strtype = 'r';
    tshiftwin.num = [1 2];
    tshiftwin.val = {[0 0]};
    tshiftwin.help = {'Time shift window (ms). [0 0] - instantaneous'};

    tshiftres = cfg_entry;
    tshiftres.tag = 'tshiftres';
    tshiftres.name = 'Time shift resolution';
    tshiftres.strtype = 'r';
    tshiftres.num = [1 1];
    tshiftres.val = {5};
    tshiftres.help = {'Time shift resolution (ms)'};

    chanset = cfg_branch;
    chanset.tag = 'chanset';
    chanset.name = 'Set';
    chanset.val = {cvachan, refchan, ncomp, outlabel, foi, tshiftwin, tshiftres};

    chansets = cfg_repeat;
    chansets.tag = 'chansets';
    chansets.name = 'Channel sets';
    chansets.values = {chanset};
    chansets.num = [1 Inf];
    chansets.val = {chanset};

    cva = cfg_branch;
    cva.tag = 'cva';
    cva.name = 'CVA';
    cva.val = {chansets};

    res = cva;

    return
end

D = S.D;

nsets = numel(S.chanset);
badind = D.badchannels;

% Assuming projecting to columns
montage = [];
montage.labelorg    = D.chanlabels;
montage.labelnew    = {};
montage.chantypenew = {};
montage.tra         = zeros(0, D.nchannels);

spm('Pointer', 'Watch');drawnow;
spm_progress_bar('Init', nsets, 'Channel sets processed'); drawnow;
if nsets > 100, Ibar = floor(linspace(1, nsets,100));
else Ibar = 1:nsets; end

for i = 1:nsets

    spm_progress_bar('Set','ylabel','preparing data...');

    cvaind = setdiff(D.selectchannels(spm_cfg_eeg_channel_selector(S.chanset(i).cvachan.channels)), badind);
    refind = setdiff(D.selectchannels(spm_cfg_eeg_channel_selector(S.chanset(i).refchan.channels)), badind);

    if any(S.chanset(i).tshiftwin)
        tshiftind = S.chanset(i).tshiftwin(1):S.chanset(i).tshiftres:S.chanset(i).tshiftwin(2);
        tshiftind = [0 round(1e-3*D.fsample*tshiftind)];
        tshiftind = repmat(1:D.nsamples, length(tshiftind), 1)+repmat(tshiftind(:), 1, D.nsamples);
        tshiftind = tshiftind(:, all(tshiftind>0 & tshiftind<=D.nsamples));
    else
        tshiftind = repmat(1:D.nsamples, 2, 1);
    end

    nshift = size(tshiftind, 1)-1;
    nt     = size(tshiftind, 2);

    Y  = zeros(length(cvaind), nt, D.ntrials);
    Yr = zeros(length(refind)*nshift, nt, D.ntrials);

    for j = 1:D.ntrials
        cY  = D(cvaind, :, j);
        cYr = D(refind, :, j);

        if S.chanset(i).foi(1) > 0
            cY  = ft_preproc_highpassfilter(cY, D.fsample, S.chanset(i).foi(1),...
                5, 'but', 'twopass', 'reduce');
            cYr = ft_preproc_highpassfilter(cYr, D.fsample, S.chanset(i).foi(1),...
                5, 'but', 'twopass', 'reduce');
        end

        if isfinite(S.chanset(i).foi(2))
            cY  = ft_preproc_lowpassfilter(cY, D.fsample, S.chanset(i).foi(2),...
                5, 'but', 'twopass', 'reduce');
            cYr = ft_preproc_lowpassfilter(cYr, D.fsample, S.chanset(i).foi(2),...
                5, 'but', 'twopass', 'reduce');
        end

        Y(:, :, j) = cY(:, tshiftind(1, :));

        for k = 1:size(cYr, 1)
            Yr(((k-1)*nshift+1):k*nshift, :, j) = reshape(cYr(k, tshiftind(2:end, :)'),[], nshift)';
        end
    end

    Y  = reshape(Y,  size(Y, 1), []);
    Yr = reshape(Yr, size(Yr, 1), []);

    if D.fsample>2.5*S.chanset(i).foi(2)
        dec = floor(D.fsample/(2.5*S.chanset(i).foi(2)));
        Y   =  Y(:, 1:dec:end);
        Yr  = Yr(:, 1:dec:end);
    end

    spm_progress_bar('Set','ylabel','running CVA...');

    CVA = spm_cva(Y', Yr');

    ncomp = min(S.chanset(i).ncomp, size(CVA.V, 2));
    for j = 1:ncomp
        if ncomp == 1
            montage.labelnew{end+1, 1} = S.chanset(i).outlabel;
        else
            montage.labelnew{end+1, 1} = [S.chanset(i).outlabel num2str(j)];
        end

        montage.tra(end+1, end)  = 0;
        montage.tra(end, cvaind) = CVA.V(:, j)';

        montage.chantypenew{end+1}='LFP';
    end

    if ismember(i, Ibar)
        spm_progress_bar('Set', i); drawnow;
    end
end

spm_progress_bar('Clear');

if ~isempty(S.chanind)
    montage.labelnew = [montage.labelnew; D.chanlabels(S.chanind)'];
    I = eye(D.nchannels);
    montage.tra = [montage.tra; I(S.chanind, :)];
    montage.chantypenew = [montage.chantypenew, D.chantype(S.chanind)];
end

res = montage;


text:
function [pos, tri] = mesh_icosahedron

% MESH_ICOSAHEDRON returns the vertices and triangle of a 12-vertex icosahedral
% mesh.
%
% Use as
%   [pos, tri] = mesh_icosahedron
%
% See also MESH_TETRAHEDRON, MESH_OCTAHEDRON, MESH_SPHERE

% Copyright (C) 2002, Robert Oostenveld
% Copyright (C) 2019, Robert Oostenveld and Jan-Mathijs Schoffelen
%
% 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/>.
%

tri = [
  1   2   3
  1   3   4
  1   4   5
  1   5   6
  1   6   2
  2   8   3
  3   9   4
  4  10   5
  5  11   6
  6   7   2
  7   8   2
  8   9   3
  9  10   4
  10  11   5
  11   7   6
  12   8   7
  12   9   8
  12  10   9
  12  11  10
  12   7  11
  ];

pos = zeros(12, 3);

rho = 0.4*sqrt(5);
phi = 2*pi*(0:4)/5;

pos( 1, :) = [0 0  1];          % top point

pos(2:6, 1) = rho*cos(phi)';
pos(2:6, 2) = rho*sin(phi)';
pos(2:6, 3) = rho/2;

pos(7:11, 1) = rho*cos(phi - pi/5)';
pos(7:11, 2) = rho*sin(phi - pi/5)';
pos(7:11, 3) = -rho/2;

pos(12, :) = [0 0 -1];          % bottom point

% scale all vertices to the unit sphere
pos = pos ./ repmat(sqrt(sum(pos.^2,2)), 1,3);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ r, seed ] = r8_uniform_01 ( seed )

%*****************************************************************************80
%
%% R8_UNIFORM_01 returns a unit pseudorandom R8.
%
%  Discussion:
%
%    This routine implements the recursion
%
%      seed = 16807 * seed mod ( 2**31 - 1 )
%      r8_uniform_01 = seed / ( 2**31 - 1 )
%
%    The integer arithmetic never requires more than 32 bits,
%    including a sign bit.
%
%    If the initial seed is 12345, then the first three computations are
%
%      Input     Output      R8_UNIFORM_01
%      SEED      SEED
%
%         12345   207482415  0.096616
%     207482415  1790989824  0.833995
%    1790989824  2035175616  0.947702
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Springer Verlag, pages 201-202, 1983.
%
%    Pierre L'Ecuyer,
%    Random Number Generation,
%    in Handbook of Simulation,
%    edited by Jerry Banks,
%    Wiley Interscience, page 95, 1998.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, pages 362-376, 1986.
%
%    Peter Lewis, Allen Goodman, James Miller,
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, pages 136-143, 1969.
%
%  Parameters:
%
%    Input, integer SEED, the integer "seed" used to generate
%    the output random number.  SEED should not be 0.
%
%    Output, real R, a random value between 0 and 1.
%
%    Output, integer SEED, the updated seed.  This would
%    normally be used as the input seed on the next call.
%
  if ( seed == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_UNIFORM_01 - Fatal error!\n' );
    fprintf ( 1, '  Input SEED = 0!\n' );
    error ( 'R8_UNIFORM_01 - Fatal error!' );
  end

  seed = floor ( seed );

  seed = mod ( seed, 2147483647 );

  if ( seed < 0 )
    seed = seed + 2147483647;
  end

  k = floor ( seed / 127773 );

  seed = 16807 * ( seed - k * 127773 ) - k * 2836;

  if ( seed < 0 )
    seed = seed + 2147483647;
  end

  r = seed * 4.656612875E-10;

  return
end


text:
function b = includes(s1, s2)

% Sep 4, 2020: Removed verLessThan because it's slowing down Homer3 startup

% if verLessThan('matlab','9.1')
%     b = ~isempty(strfind(s1,s2));
% else
%     b = contains(s1,s2);
% end

if iscell(s1)
    b = ~isempty(find(strcmp(s1, s2)));
else
    b = ~isempty(strfind(s1,s2)); %#ok<STREMP>
    % b = contains(s1,s2);
end

author: BUNPC


text:
function wc = FWT2_PO(x,L,qmf)
% FWT2_PO -- 2-d MRA wavelet transform (periodized, orthogonal)
%  Usage
%    wc = FWT2_PO(x,L,qmf)
%  Inputs
%    x     2-d image (n by n array, n dyadic)
%    L     coarse level
%    qmf   quadrature mirror filter
%  Outputs
%    wc    2-d wavelet transform
%
%  Description
%    A two-dimensional Wavelet Transform is computed for the
%    array x.  To reconstruct, use IWT2_PO.
%
%  See Also
%    IWT2_PO, MakeONFilter
%
	[n,J] = quadlength(x);
	wc = x;
	nc = n;
	for jscal=J-1:-1:L, % from fine (J) to coarse (L)
		top = (nc/2+1):nc; bot = 1:(nc/2);
		for ix=1:nc,
			row = wc(ix,1:nc);
			wc(ix,bot) = DownDyadLo(row,qmf);
			wc(ix,top) = DownDyadHi(row,qmf);
		end
		for iy=1:nc,
			row = wc(1:nc,iy)';
			wc(top,iy) = DownDyadHi(row,qmf)';
			wc(bot,iy) = DownDyadLo(row,qmf)';
		 end
		nc = nc/2; % downsampling
	end

%
% Copyright (c) 1993. David L. Donoho
%

%
%  Part of Wavelab Version 850
%  Built Tue Jan  3 13:20:40 EST 2006
%  This is Copyrighted Material
%  For Copying permissions see COPYING.m
%  Comments? e-mail wavelab@stat.stanford.edu

author: thomaskuestner


text:
function J=calcCartRRJacob(components,xState,useHalfRange,lTx,lRx)
%%CALCCARTRRJACOB Compute the Jacobian matrix for a Cartesian measurement
%           possibly with a bistatic range rate component in 1D, 2D or 3D
%           for a Cartesian target state. This function can be useful when
%           using an extended Kalman filter (EKF) with Cartesian-converted
%           measurements and range rate. Without the range rate component,
%           the result is just the measurement matrix, which, when
%           multiplied by the target state, extracts the position
%           components and is often designated by H in the discrete-time
%           Kalman filter.
%
%INPUTS: components This specified whether the output should contain a
%               range-rate component. Possible values are:
%               0 The output contains position and range rate rows. This is
%                 the default if an empty matrix is passed.
%               1 The output just has position rows.
%        xState The xDimX1 target state vector in the global coordinate
%               system with [x;y;z;xDot;yDot;zDot] components in 3D or
%               [x;y;xDot;yDot] components in 2D.
%  useHalfRange A boolean value specifying whether the bistatic range
%               (and thus the range rate) value has been divided by two.
%               This normally comes up when operating in monostatic mode,
%           lTx The transmitter state vector in the global coordinate
%               system with [x;y;z;xDot;yDot;zDot] components in 3D or
%               [x;y;xDot;yDot] components in 2D. If omitted, then a vector
%               of zeros is used.
%           lRx The receiver state vector in the global coordinate system
%               with [x;y;z;xDot;yDot;zDot] components in 3D or
%               [x;y;xDot;yDot] components in 2D. If omitted, then a vector
%               of zeros is used.
%
%OUTPUTS: J The (xDim/2+1)XxDim Jacobian matrix. Each row is a component
%           of x-y-z-range rate and each column is the derivative with
%           respect to x,y,z,xDot,yDot,zDot. If range rate is not desired,
%           then the matrix is (xDim/2)XxDim in size.
%
%The derivatives of the position components with themselves are clear (1
%for common components and zero for different components). The derivatives
%of the range rate components are obtained using the rangeRateGradient
%function.
%
%February 2017 David F.Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

xDim=size(xState,1);
posDim=xDim/2;
measDim=posDim+1;

if(isempty(components))
    components=0;
end

if(components==0)
    if(nargin<3||isempty(useHalfRange))
        useHalfRange=false;
    end

    if(nargin<4||isempty(lTx))
        lTx=zeros(xDim,1);
    end

    if(nargin<5||isempty(lRx))
        lRx=zeros(xDim,1);
    end

    J=zeros(measDim,xDim);
    %Position components are just passed through.
    J(1:posDim,1:posDim)=eye(posDim,posDim);

    J(posDim+1,:)=rangeRateGradient(xState,useHalfRange,lTx,lRx);
else
    J=zeros(posDim,xDim);
    J(1:posDim,1:posDim)=eye(posDim,posDim);
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 dx = simpleSysWrapper(~,x,w)
% System dynamics part
u = w'*Psi_fun(x);
dx = susp_sys(x,u);       % dx as the first 1-4 states of the wrapper
end

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function fx = annulus ( n, x )

%*****************************************************************************80
%
%% ANNULUS evaluates the annulus function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 December 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points.
%
%    Input, real X(N,2), the point coordinates.
%
%    Output, real FX(N), the function values.
%
  fx = zeros ( n, 1 );

  c1 = 0.25 * ones ( n, 2 );
  r1 = 0.5 * ones ( n, 1 );

  c2 = 0.0 * ones ( n, 2 );
  r2 = 1.0 * ones ( n, 1 );

  a = ( sqrt ( sum ( ( x - c1 ) .* ( x - c1 ), 2 ) ) ) - r1;
  b = r2 - ( sqrt ( sum ( ( x - c2 ) .* ( x - c2 ), 2 ) ) );

  fx = a .* b;

  return
end


text:
function segNum = get_segNum
%Get the previously highest value of segNum for all displayed active contours
%LM:  18 April 02, JOD.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

%Get handles to all the existing plotted segments
hSegsV = findobj('tag','newSegment');

segNum = 0;
for i = 1 : length(hSegsV)
  segS = get(hSegsV(i),'userdata');
  if segS.segNum > segNum
    segNum = segS.segNum;
  end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [Ioutput]=zerocrossing(Izerosmooth)

[rr,cc]=size(Izerosmooth);
Ioutput=zeros([rr,cc]);

for i=2:rr-1
    for j=2:cc-1
        if (Izerosmooth(i,j)>0)
             if (Izerosmooth(i,j+1)>=0 && Izerosmooth(i,j-1)<0) || (Izerosmooth(i,j+1)<0 && Izerosmooth(i,j-1)>=0)

                Ioutput(i,j)= Izerosmooth(i,j+1);

            elseif (Izerosmooth(i+1,j)>=0 && Izerosmooth(i-1,j)<0) || (Izerosmooth(i+1,j)<0 && Izerosmooth(i-1,j)>=0)
                    Ioutput(i,j)= Izerosmooth(i,j+1);
            elseif (Izerosmooth(i+1,j+1)>=0 && Izerosmooth(i-1,j-1)<0) || (Izerosmooth(i+1,j+1)<0 && Izerosmooth(i-1,j-1)>=0)
                  Ioutput(i,j)= Izerosmooth(i,j+1);
            elseif (Izerosmooth(i-1,j+1)>=0 && Izerosmooth(i+1,j-1)<0) || (Izerosmooth(i-1,j+1)<0 && Izerosmooth(i+1,j-1)>=0)
                  Ioutput(i,j)=Izerosmooth(i,j+1);
            end

        end

    end
end


text:
%% Optical flow evaluation demo
%
% Computes flow field between two images using various methods and display it
% (deepflow, simpleflow, sparsetodenseflow, Farneback, TV-L1).
%
% Sources:
%
% * <https://github.com/opencv/opencv_contrib/blob/3.1.0/modules/optflow/samples/optical_flow_evaluation.cpp>
%

%% Input images
% a pair of 8-bit color images
im1 = imread(fullfile(mexopencv.root(),'test','RubberWhale1.png'));
im2 = imread(fullfile(mexopencv.root(),'test','RubberWhale2.png'));
assert(isequal(size(im1), size(im2)), ...
    'Dimension mismatch between input images');
if ~mexopencv.isOctave() && mexopencv.require('images')
    %HACK: IMSHOWPAIR not implemented in Octave
    imshowpair(im1, im2);
end

%% Compare the different methods
algorithms = {'farneback', 'simpleflow', 'tvl1', 'deepflow', ...
    'sparsetodenseflow', 'pcaflow', ...
    'DISflow_ultrafast', 'DISflow_fast', 'DISflow_medium', 'variational'};
for i=1:numel(algorithms)
    % prepare images
    if any(strcmp(algorithms{i}, {'farneback', 'tvl1', 'deepflow', ...
            'DISflow_ultrafast', 'DISflow_fast', 'DISflow_medium', ...
            'variational'})) && size(im1,3)==3
        % 1-channel images are expected
        img1 = cv.cvtColor(im1, 'RGB2GRAY');
        img2 = cv.cvtColor(im2, 'RGB2GRAY');
    elseif strcmp(algorithms{i}, 'simpleflow') && size(im1,3)==1
        % 3-channel images expected
        img1 = cv.cvtColor(im1, 'GRAY2RGB');
        img2 = cv.cvtColor(im2, 'GRAY2RGB');
    else
        % sparsetodenseflow/pcaflow handle both 1- or 3-channels
        img1 = im1;
        img2 = im2;
    end

    % compute flow field between img1 and img2 using current method
    tic
    switch lower(algorithms{i})
        case 'farneback'
            %{
            obj = cv.FarnebackOpticalFlow();
            flow = obj.calc(img1, img2);
            %}
            flow = cv.calcOpticalFlowFarneback(img1, img2);
        case 'simpleflow'
            flow = cv.calcOpticalFlowSF(img1, img2);
        case 'deepflow'
            flow = cv.calcOpticalFlowDF(img1, img2);
        case 'sparsetodenseflow'
            flow = cv.calcOpticalFlowSparseToDense(img1, img2);
        case 'tvl1'
            obj = cv.DualTVL1OpticalFlow();
            flow = obj.calc(img1, img2);
        case 'pcaflow'
            %obj = cv.OpticalFlowPCAFlow(prior);  % path to a prior file for PCAFlow
            obj = cv.OpticalFlowPCAFlow();
            flow = obj.calc(img1, img2);
        case 'disflow_ultrafast'
            obj = cv.DISOpticalFlow('Preset','UltraFast');
            flow = obj.calc(img1, img2);
        case 'disflow_fast'
            obj = cv.DISOpticalFlow('Preset','Fast');
            flow = obj.calc(img1, img2);
        case 'disflow_medium'
            obj = cv.DISOpticalFlow('Preset','Medium');
            flow = obj.calc(img1, img2);
        case 'variational'
            obj = cv.VariationalRefinement();
            flow = obj.calc(img1, img2);
    end
    fprintf('%18s: ', algorithms{i});
    toc

    % display the flow
    if true
        [mag, ang] = cv.cartToPolar(flow(:,:,1), flow(:,:,2), 'Degrees',true);
    else
        [ang, mag] = cart2pol(flow(:,:,1), flow(:,:,2));
        if mexopencv.isOctave()
            %HACK: RAD2DEG not implemented in Octave
            ang = (ang + pi) * (180 / pi);
        else
            ang = rad2deg(ang + pi);
        end
    end
    mag = cv.normalize(mag, 'Alpha',0, 'Beta',1, 'NormType','MinMax');
    hsv = cat(3, ang, ones(size(ang),class(ang)), mag); % H=[0,360], S,V=[0,1]
    rgb = cv.cvtColor(hsv, 'HSV2RGB');                  % R,G,B=[0,1]
    figure, imshow(rgb)
    title(sprintf('Computed flow: %s',algorithms{i}), 'Interpreter','none')
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [resRates] = ma_getResRates(massLoss)
%getResRates Summary of this function goes here
%   Detailed explanation goes here
    maxNumRes = max(horzcat(massLoss.lossConvert.resLost, massLoss.lossConvert.resConvert));
    resRates = zeros(1,maxNumRes);

    for(i=1:length(massLoss.lossConvert)) %#ok<*NO4LP>
        lossConvert = massLoss.lossConvert(i);
        resRates(lossConvert.resLost) = resRates(lossConvert.resLost) - lossConvert.resLostRates(lossConvert.resLost);

        sumRates = sum(lossConvert.resLostRates);
        resRates(lossConvert.resConvert) = resRates(lossConvert.resConvert) + sumRates*lossConvert.resConvertPercent(lossConvert.resConvert);
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [x,fval,exitflag,info] = opti_filtersd(fun,grad,x0,lb,ub,nlcon,nljac,nljacstr,cl,cu,opts)
%OPTI_FILTERSD Solve a NLP using FILTERSD
%
%   min f(x)       subject to:    cl <= nlcon(x) <= cu
%    x                            lb <= x <= ub
%
%   x = opti_filtersd(fun,grad,x0) solves a NLP where fun is the objective
%   function, grad is the gradient of the objective and x0 is a starting
%   guess.
%
%   x = opti_filtersd(fun,grad,x0,lb,ub) solves subject to decision
%   variable bounds lb <= x <= ub. Infinite bounds are OK.
%
%   x = opti_filtersd(fun,x0,lb,ub,nlcon,nljac,[],cl,cu) solves subject to
%   the nonlinear inequality constraints cl <= nlcon(x) <= cu, and nljac
%   returns a DENSE Jacobian of the constraints.
%
%   x = opti_filtersd(fun,x0,lb,ub,nlcon,nljac,nljacstr,cl,cu) solves as
%   above, except nljac returns a SPARSE Jacobian of the constraints, and
%   nljacstr supplies a sparse matrix of ALL possible non-zero locations
%   within the Jacobian.
%
%   x = opti_filtersd(fun,...,cu,opts) uses opts to pass optiset options to
%   the solver.
%
%   [x,fval,exitflag,info] = opti_filtersd(...) returns the objective value
%   at the solution, together with the solver exitflag, and an information
%   structure.
%
%   THIS IS A WRAPPER FOR FILTERSD
%   See supplied Eclipse Public License

%   Copyright (C) 2013 Jonathan Currie (IPL)

if(nargin < 11), opts = optiset; else opts = optiset(opts); end
if(nargin < 10), cu = []; end
if(nargin < 9),  cl = []; end
if(nargin < 8),  nljacstr = []; end
if(nargin < 7),  nljac = []; end
if(nargin < 6),  nlcon = []; end
if(nargin < 5),  ub = []; end
if(nargin < 4),  lb = []; end
if(nargin < 3),  error('FILTERSD requires at least 3 arguments'); end

%Default is not sparse
sp = false;

%Determine we have a sparse problem
if(~isempty(nljacstr))
    sp = true; %must be sparse if we passed an arg here
    if(~isa(nljacstr,'function_handle'))
        error('Jacobian Structure must be a function handle');
    end
end

%Check we have a valid x0
if(isempty(x0) || any(isnan(x0)))
    error('FILTERSD requires an initial guess, x0!');
end

%Setup display level
opts.display = dispLevel(opts.display);
opts.optiver = optiver;

t = tic;
if(sp)
    % Run FILTERSD [Sparse]
    jacstr = nljacstr(); %filtersd just requires the matrix
    [x, fval, exitflag, stats, lambda] = filtersdsp(fun,grad,x0,lb,ub,nlcon,nljac,jacstr,cl,cu,opts);
else
    % Run FILTERSD [Dense]
    [x, fval, exitflag, stats, lambda] = filtersd(fun,grad,x0,lb,ub,nlcon,nljac,cl,cu,opts);
end

%Collect Results
info.Iterations = stats.niter;
info.FuncEvals = stats.nfval;
info.GradEvals = stats.ngval;
info.Time = toc(t);
if(sp)
    info.Algorithm = 'FILTERSD: Nonlinear Optimization using the Filter Search & Trust Region [Sparse]';
else
    info.Algorithm = 'FILTERSD: Nonlinear Optimization using the Filter Search & Trust Region [Dense]';
end

switch(exitflag)
    case 0
        info.Status = 'Locally Optimal';
        exitflag = 1;
    case {5,101,102}
        info.Status = 'Exceeded Iterations / Function Evaluations / Time';
        exitflag = 0;
    case {1,2,4,7,8}
        info.Status = 'Initialization Error';
        exitflag = -2;
    case {3,6}
        info.Status = 'Infeasible / Inconsistent Constraints';
        exitflag = -1;
    case {9,10}
        info.Status = 'LCP Solver Error';
        exitflag = -3;
    case 105
        info.Status = 'User Exited';
        exitflag = -5;
    otherwise
        info.Status = 'FILTERSD Error';
end

info.Lambda = lambda;


text:
function failed_bug2359

% MEM 2gb
% WALLTIME 00:30:00
% DEPENDENCY ft_prepare_mesh ft_prepare_sourcemodel

cd(dccnpath('/home/common/matlab/fieldtrip/data/test/bug2359'));

cortex = ft_read_headshape('cortex_20484.surf.gii');
iskull = ft_read_headshape('iskull_2562.surf.gii');
oskull = ft_read_headshape('oskull_2562.surf.gii');
scalp  = ft_read_headshape('scalp_2562.surf.gii');

figure
ft_plot_mesh(cortex, 'facecolor', 'b');
ft_plot_mesh(iskull, 'facecolor', [0.5 0.5 0.5], 'edgecolor', 'none', 'facealpha', 0.5);

%% PART 1, implement the correction of the cortical sheet by moving points that are too close to the surface inward

% usually the "moveinward" cortex modification would be done for EEG-BEM, but it is
% easier to test it for a MEG singleshell model (for which it is actually not
% needed).

cfg = [];
cfg.method = 'singleshell';
vol = ft_prepare_headmodel(cfg, iskull);

cfg = [];
cfg.headmodel = vol;
cfg.sourcemodel = cortex;    % this is in mm
cfg.inwardshift = 0;  % this should be expressed in the units consistent with cfg.unit
cfg.moveinward = 0;
gridorig = ft_prepare_sourcemodel(cfg);

cfg = [];
cfg.headmodel = vol;
cfg.sourcemodel = cortex;      % this is in mm
cfg.inwardshift = -5;  % outward shifted
cfg.moveinward = 0;
gridoutward = ft_prepare_sourcemodel(cfg);

cfg = [];
cfg.headmodel = vol;
cfg.sourcemodel = cortex;    % this is in mm
cfg.inwardshift = 5; % inward shifted
cfg.moveinward = 0;
gridinward = ft_prepare_sourcemodel(cfg);

cfg = [];
cfg.headmodel = vol;
cfg.sourcemodel = cortex;    % this is in mm
cfg.inwardshift = 0;  % keep this at the original place
cfg.moveinward = 5;  % dipoles moved inwards
gridcorrect = ft_prepare_sourcemodel(cfg);

assert(isempty(gridorig.outside));
assert(isempty(gridoutward.outside));
assert(~isempty(gridinward.outside)); % this should have a few vertices outside the inward shifted surface
assert(isempty(gridcorrect.outside))

figure
ft_plot_mesh(iskull, 'facecolor', [0.5 0.5 0.5], 'edgecolor', 'none', 'facealpha', 0.5);
ft_plot_mesh(gridorig.pos(gridinward.outside,:), 'vertexcolor', 'b');
ft_plot_mesh(gridcorrect.pos(gridinward.outside,:), 'vertexcolor', 'r');

%% PART 2, implement the spherify of the cortical sheet

mesh = cat(1, iskull, oskull, scalp);

cfg = [];
cfg.method = 'concentricspheres';
vol = ft_prepare_headmodel(cfg, mesh);

cfg = [];
cfg.headmodel = vol;
cfg.sourcemodel = cortex;    % this is in mm
cfg.spherify = 'yes';
gridsphere = ft_prepare_sourcemodel(cfg);

assert(isempty(gridsphere.outside));

figure
ft_plot_headmodel(vol, 'edgecolor', 'none', 'facecolor', 'skin', 'facealpha', 0.5);
ft_plot_mesh(gridsphere)

%% this is a weird modification, I am just curious to see how it works

cfg = [];
cfg.method = 'singlesphere';
vol = ft_prepare_headmodel(cfg, iskull);

cfg = [];
cfg.headmodel = vol;
cfg.sourcemodel.xgrid = -200:10:200;    % this is in mm
cfg.sourcemodel.ygrid = -200:10:200;    % this is in mm
cfg.sourcemodel.zgrid =  -50:10:150;    % this is in mm
cfg.spherify = 'yes';
gridsphere = ft_prepare_sourcemodel(cfg);

figure
ft_plot_headmodel(vol, 'edgecolor', 'none', 'facecolor', 'skin', 'facealpha', 0.5);
ft_plot_mesh(gridsphere)


text: any(cellfun(@exist,{'lsl_pull_sample_buf', 'lsl_push_sample'})==3)

repo: REST


text:
function [feat,idxs_bbox_pair,cb1,cb2,o] = get_spatial_features_same_img(locs,idxs_bbox_pair,scores,class_id,bVis)

if (nargin < 3)
    scores = [];
end

if (nargin < 5)
    bVis = false;
end

if (nargin < 6)
    lab = ones(size(locs,1),1);
end

half_size = 132;

xs = locs(:, 1);
ys = locs(:, 2);
boxes = [xs-half_size, ys-half_size, xs+half_size, ys+half_size];

cb1 = [mean(boxes(idxs_bbox_pair(:,1),[1 3]),2) mean(boxes(idxs_bbox_pair(:,1),[2 4]),2)];
cb2 = [mean(boxes(idxs_bbox_pair(:,2),[1 3]),2) mean(boxes(idxs_bbox_pair(:,2),[2 4]),2)];

h1 = boxes(idxs_bbox_pair(:,1),4) - boxes(idxs_bbox_pair(:,1),2);
h2 = boxes(idxs_bbox_pair(:,2),4) - boxes(idxs_bbox_pair(:,2),2);
assert(any(h1>0));
assert(any(h2>0));

h_m = (h1+h2)/2;
offset = abs(cb1 - cb2)./[h_m h_m];
sc = abs(h1 - h2)./h_m;
offsetM = offset(:,1).*offset(:,2);
[o,o2,o3] = boxoverlapMx(boxes(idxs_bbox_pair(:,1),:), boxes(idxs_bbox_pair(:,2),:));

feat = cat(2, offset, offsetM, sc, o, o2, o3, ...
            offset.^2, offsetM.^2, sc.^2, o.^2, o2.^2, o3.^2, ...
            exp(-offset), exp(-offsetM), exp(-sc), exp(-o), exp(-o2), exp(-o3));
if (~isempty(scores))
    scores1 = scores(idxs_bbox_pair(:,1),:);
    scores2 = scores(idxs_bbox_pair(:,2),:);
    feat = cat(2, feat, scores1, scores2);
end

if (bVis)
    idxs = 1:min(size(cb1,1),100);

    plot(cb1(idxs,1),cb1(idxs,2),'b+','MarkerSize',10);
    plot(cb2(idxs,1),cb2(idxs,2),'g+','MarkerSize',10);

    for i = 1:length(idxs)
        if (lab(idxs(i)) == 1)
            plot([cb1(idxs(i),1); cb2(idxs(i),1)],[cb1(idxs(i),2); cb2(idxs(i),2)],'r-','lineWidth',1);
        else
            plot([cb1(idxs(i),1); cb2(idxs(i),1)],[cb1(idxs(i),2); cb2(idxs(i),2)],'b-','lineWidth',1);
        end
    end

    legendName = {['pidx' num2str(class_id)]};
    legend(legendName);
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function e = end(f, k, n)
%END   Rightmost point of a CHEBFUN domain (or last row/col of quasimatrix).
%   END(F, K, N) This function is called for indexing expressions involving a
%   CHEBFUN F when END is part of the K-th index out of N indices.
%
%   If F is a column CHEBFUN and K is 1, this function returns the last point
%   in F's domain. If K is 2, it returns the index of F's last column.
%
%   If F is a row CHEBFUN and K is 1, this function returns the index of F's
%   last row. If K is 2, it returns the last point in F's domain.

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

if ( n > 2 )
    error('CHEBFUN:CHEBFUN:end:ngt2', 'Index exceeds CHEBFUN dimensions.');
end

if ( ((k == 2) && ~f(1).isTransposed) || ...
     ((k == 1) && f(1).isTransposed && (n > 1)) )
    % 'end' row/column of the array-valued CHEBFUN.
    if ( isempty(f) )
        e = 0;
    else
        e = numColumns(f);
    end

else
    % 'end' of the domain.
    e = f(1).domain(end);

end

end


text:
function [B,A] = octdsgn(Fc,Fs,N);
% OCTDSGN  Design of an octave filter.
%    [B,A] = OCTDSGN(Fc,Fs,N) designs a digital octave filter with
%    center frequency Fc for sampling frequency Fs.
%    The filter are designed according to the Order-N specification
%    of the ANSI S1.1-1986 standard. Default value for N is 3.
%    Warning: for meaningful design results, center values used
%    should preferably be in range Fs/200 < Fc < Fs/5.
%    Usage of the filter: Y = FILTER(B,A,X).
%
%    Requires the Signal Processing Toolbox.
%
%    See also OCTSPEC, OCT3DSGN, OCT3SPEC.

% Author: Christophe Couvreur, Faculte Polytechnique de Mons (Belgium)
%         couvreur@thor.fpms.ac.be
% Last modification: Aug. 22, 1997, 9:00pm.

% References:
%    [1] ANSI S1.1-1986 (ASA 65-1986): Specifications for
%        Octave-Band and Fractional-Octave-Band Analog and
%        Digital Filters, 1993.

if (nargin > 3) | (nargin < 2)
  error('Invalide number of arguments.');
end
if (nargin == 2)
  N = 3;
end
if (Fc > 0.70*(Fs/2))
  error('Design not possible. Check frequencies.');
end

% Design Butterworth 2Nth-order octave filter
% Note: BUTTER is based on a bilinear transformation, as suggested in [1].
%W1 = Fc/(Fs/2)*sqrt(1/2);
%W2 = Fc/(Fs/2)*sqrt(2);
pi = 3.14159265358979;
beta = pi/2/N/sin(pi/2/N);
alpha = (1+sqrt(1+8*beta^2))/4/beta;
W1 = Fc/(Fs/2)*sqrt(1/2)/alpha;
W2 = Fc/(Fs/2)*sqrt(2)*alpha;
[B,A] = butter(N,[W1,W2]);


text:
function Cs = WENStoClockwise(D);
% WENStoClockwise: Converts a data structure where wedge orientations
%                  are organized by cardinal point (in the order West
%                  East North South) into another structure with
%                  orientations starting from top left wedge and
%                  increasing in a clockwise fashion.
%
% See Also
%   SeparateAngles, DetailCurveCoeff, ClockwisetoWENS

   [nq,nw,nl,ns] = size(D);
   Cs = cell(1,nq*nw);
   cnt = 1;

    for w = 1:nw
      Cs{cnt} = squeeze(D(3,w,:,:)); cnt = cnt+1;
    end
    for w=nw:-1:1
      Cs{cnt} = squeeze(D(2,w,:,:)).'; cnt = cnt + 1;
    end
    for w=1:nw
      Cs{cnt} = squeeze(D(4,w,:,:)); cnt = cnt+1;
    end
    for w=nw:-1:1
      Cs{cnt} = squeeze(D(1,w,:,:)).'; cnt = cnt + 1;
    end


text:
%Simple script to quickly demo use of CalCamArm.m

%path to images of checkerboards
imageFolder = './Example Data/Images/';
%loading arm transformations
load('./Example Data/armMat.mat');
%checkerboard square widths in mm
squareSize = 13.6;

%run calibration
[TBase, TEnd, cameraParams, TBaseStd, TEndStd, pixelErr] = CalCamArm(imageFolder, armMat, squareSize,'maxBaseOffset',0.5);

%print results
fprintf('\nFinal camera to arm base transform is\n')
disp(TBase);

fprintf('Final end effector to checkerboard transform is\n')
disp(TEnd);

fprintf('Final camera matrix is\n')
disp(cameraParams.IntrinsicMatrix');

fprintf('Final camera radial distortion parameters are\n')
disp(cameraParams.RadialDistortion);

fprintf('Final camera tangential distortion parameters are\n')
disp(cameraParams.TangentialDistortion);

author: ZacharyTaylor

repo: Camera-to-Arm-Calibration


text:
function partition_problem_test02 ( n, w )

%*****************************************************************************80
%
%% PARTITION_PROBLEM_TEST02 tests PARTITION_COUNT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 May 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of weights.
%
%    Input, integer W(N), a set of weights.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PARTITION_PROBLEM_TEST02:\n' );
  fprintf ( 1, '  PARTITION_COUNT counts the number of exact solutions\n' );
  fprintf ( 1, '  of the partition problem.\n' );

  count = partition_count ( n, w );

  fprintf ( 1, '\n' );
  fprintf ( 1, '     I        W\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %4d  %8d\n', i, w(i) );
  end
  fprintf ( 1, '\n' );
  fprintf ( 1, 'Number of solutions = %d\n', count );

  return
end


text:
function table = hn_exponential_product ( p, b )

%*****************************************************************************80
%
%% HN_EXPONENTIAL_PRODUCT: exponential products exp(b*x)*Hn(i,x)*Hn(j,x).
%
%  Discussion:
%
%    Hn(i,x) is the normalized physicist's Hermite polynomial of degree I.
%
%    For polynomial chaos applications, it is of interest to know the
%    value of the integrals of products of exp(B*X) with every possible pair
%    of basis functions.  That is, we'd like to form
%
%      Tij = Integral ( -oo < X < +oo ) exp(B*X) * Hn(I,X) * Hn(J,X) exp(-X*X) dx
%
%    We will estimate these integrals using Gauss-Hermite quadrature.
%    Because of the exponential factor exp(B*X), the quadrature will not
%    be exact.
%
%    However, when B = 0, the quadrature is exact, and moreoever, the
%    table will be the identity matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer P, the maximum degree of the polyonomial factors.
%    0 <= P.
%
%    Input, real B, the coefficient of X in the exponential factor.
%
%    Output, real TABLE(P+1,P+1), the table of integrals.  TABLE(I,J)
%    represents the weighted integral of exp(B*X) * Hn(I+1,X) * Hn(J+1,X).
%
  table(1:p+1,1:p+1) = 0.0;

  order = floor ( ( 3 * p + 4 ) / 2 );

  [ x_table, w_table ] = h_quadrature_rule ( order );
  r8vec2_print ( order, x_table, w_table, '  Quadrature rule:' );

  for k = 1 : order

    x = x_table(k);
    h_table = hn_polynomial_value ( 1, p, x );
%
%  The following formula is an outer product in H_TABLE.
%
    table(1:p+1,1:p+1) = table(1:p+1,1:p+1) ...
      + w_table(k) * exp ( b * x ) * h_table(1:p+1)' * h_table(1:p+1);

  end

  return
end