text:
% PROBABILITY DISTRIBUTION FUNCTIONS (in the dist-folder):
%
% Priors
%   PRIOR_FIXED       Fix parameter to its current value
%   PRIOR_GAMMA       Gamma prior structure
%   PRIOR_GAUSSIAN    Gaussian prior structure
%   PRIOR_INVGAMMA    Inverse-gamma prior structure
%   PRIOR_INVT        Inverse Student-t prior structure
%   PRIOR_INVUNIF     Inverse uniform prior structure
%   PRIOR_LAPLACE     Laplace (double exponential) prior structure
%   PRIOR_LOGGAUSSIAN Log-Gaussian prior structure
%   PRIOR_LOGLOGUNIF  Uniform prior structure for the log(log(parameter))
%   PRIOR_LOGT        Student-t prior structure for the logarithm of the parameter
%   PRIOR_LOGUNIF     Uniform prior structure for the logarithm of the parameter
%   PRIOR_SINVCHI2    Scaled inverse-chi-square prior structure
%   PRIOR_SQINVGAMMA  Gamma prior structure for square inverse of the parameter
%   PRIOR_SQINVLOGUNIF Uniform prior structure for the log of the square
%                      inverse of parameter
%   PRIOR_SQINVSINVCHI2 Scaled-Inv-Chi^2 prior structure for square inverse
%                       of the parameter
%   PRIOR_SQINVUNIF   Uniform prior structure for the square inverse of the
%                     parameter
%   PRIOR_SQRTINVT    Student-t prior structure for the square root of
%                     inverse of the parameter
%   PRIOR_INVSQRTUNIF Uniform prior structure for the square root of
%                     inverse of the parameter
%   PRIOR_SQRTT       Student-t prior structure for the square root of the
%                     parameter
%   PRIOR_SQRTUNIF    Uniform prior structure for the square root of the
%                     parameter
%   PRIOR_T           Student-t prior structure
%   PRIOR_UNIF        Uniform prior structure
%
% Probability (log/cumulative) density functions
%   BETA_CDF      - Beta cumulative distribution function
%   BETA_INV      - Inverse of the beta cumulative distribution function (cdf)
%   BETA_LPDF     - Beta log-probability density function (lpdf)
%   BETA_PDF      - Beta probability density function (pdf)
%   DIR_LPDF      - Log probability density function of uniform Dirichlet
%                   distribution
%   DIR_PDF       - Probability density function of uniform Dirichlet
%                   distribution
%   GAM_CDF       - Cumulative of Gamma probability density function (cdf)
%   GAM_LPDF      - Log of Gamma probability density function (lpdf)
%   GAM_PDF       - Gamma probability density function (pdf)
%   GEO_LPDF      - Geometric log probability density function (lpdf)
%   INVGAM_LPDF   - Inverse-Gamma log probability density function
%   INVGAM_PDF    - Inverse-Gamma probability density function
%   LAPLACE_LPDF  - Laplace log-probability density function (lpdf)
%   LAPLACE_PDF   - Laplace probability density function (pdf)
%   LOGN_LPDF     - Log normal log-probability density function (lpdf)
%   LOGT_LPDF     - Log probability density function (lpdf) for log Student's T
%   MNORM_LPDF    - Multivariate-Normal log-probability density function (lpdf)
%   MNORM_PDF     - Multivariate-Normal log-probability density function (lpdf)
%   NBIN_CDF      - Negative binomial cumulative distribution function (cdf)
%   NBIN_INV      - Inverse of Negative binomial cumulative distribution function (inv)
%   NBIN_PDF      - Negative binomial probability density function (pdf)
%   NEGBIN_LPDF   - Negative binomial log probability density function
%   NEGBIN_PDF    - Negative binomial probability density function
%   NEGBINZTR_LPDF - Zero trunc. negative binomial log probability density function
%   NEGBINZTR_PDF - Zero trunc. negative binomial log probability density function
%   NORM_CDF      - Normal cumulative probability density function (cdf)
%   NORM_INV      - Inverse of the normal cumulative distribution function (cdf)
%   NORM_LPDF     - Normal log-probability density function (lpdf)
%   NORM_PDF      - Normal probability density function (pdf)
%   POISS_LPDF    - Poisson log-probability density function
%   POISS_PDF     - Poisson probability density function
%   SINVCHI2_LPDF - Scaled inverse-chi log-probability density function
%   SINVCHI2_PDF  - Scaled inverse-chi probability density function
%   T_CDF         - Student's t cumulative distribution function (cdf)
%   T_INV         - Inverse of Student's T cumulative distribution function (cdf)
%   T_LPDF        - Student's T log-probability density function (lpdf)
%   T_PDF         - Student's T probability density function (pdf)
%
% Random numbers
%   BETARAND      - Random matrices from beta distribution
%   CATRAND       - Random matrices from categorical distribution
%   DIRRAND       - Uniform dirichlet random vectors
%   EXPRAND       - Random matrices from exponential distribution
%   GAMRAND       - Random matrices from gamma distribution
%   GAMRAND1      - Random matrices from gamma distribution (mex)
%   INTRAND       - Random matrices from uniform integer distribution
%   INVGAMRAND    - Random matrices from inverse gamma distribution
%   INVGAMRAND1   - Random matrices from inverse gamma distribution (mex)
%   INVWISHRND    - Random matrices from inverse Wishart distribution
%   NORMLTRAND    - Random draws from a left-truncated normal
%                   distribution, with mean = mu, variance = sigma2
%   NORMRTRAND    - Random draws from a right-truncated normal
%                   distribution, with mean = mu, variance = sigma2
%   NORMTRAND     - Random draws from a normal truncated to interval
%   NORMTZRAND    - Random draws from a normal distribution truncated by zero
%   SINVCHI2RAND  - Random matrices from scaled inverse-chi distribution
%   TRAND         - Random numbers from Student's t-distribution
%   UNIFRAND      - Generate unifrom random numberm from interval [A,B]
%   WISHRND       - Random matrices from Wishart distribution
%
% Others
%   HAMMERSLEY    - Hammersley quasi-random sequence
%

author: gpstuff-dev

repo: gpstuff

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


text:
function create_tsne_Callback(hObject, eventdata, handles)
% Creates a 2-dimensional t-sne image of calls.

padding = 1000; % Pad the temp image by this amount, so that calls near the border still fit.
blackLevel = 10; % Subtract this value from each call image to make a nicer picture.

% Select embedding type
embeddingType = questdlg('Embed with UMAP or t-SNE?', 'Embedding Method', 't-SNE' , 'UMAP', 't-SNE');
if isempty(embeddingType); return; end
if strcmp(embeddingType, 'UMAP')
    if ~exist('run_umap.m', 'file')
        msgbox('Please download UMAP and add it to MATLAB''s path and try again')
        web('   https://www.mathworks.com/matlabcentral/fileexchange/71902-uniform-manifold-approximation-and-projection-umap');
        return
    end
end

%% Choose clustering dimensions
inputParameters= questdlg('Select the input dimensions', 'Input Dimensions', 'Variational autoencoder embeddings', 'Contour shape, frequency, and duration', 'Contour shape, frequency, and duration');
if isempty(inputParameters); return; end

% Choose to assign colors by the call classification or by pitch.
colorType = questdlg({'Color the calls by frequecy (pitch), or by cluster identity?','If coloring by cluster, you may not use pre-extracted contours'}, 'Choose Color', 'Frequency' , 'Cluster', 'Frequency');
if isempty(colorType); return; end

switch inputParameters
    case 'Contour shape, frequency, and duration'
        % Get the clustering parameters, prepare data as if performing k-means
        clusterParameters = inputdlg({'Images height (pixels)','Image width (pixels)','Perplexity (if using t-SNE, ignore for UMAP)','Max number of calls to plot (set to 0 to plot everything)', 'Shape weight','Frequency weight','Duration weight'},...
            'Choose cluster parameters:',1,...
            {'6000','6000','30','2000','1','1','1'});
        if isempty(clusterParameters); return; end

        slope_weight = str2double(clusterParameters{5});
        freq_weight = str2double(clusterParameters{6});
        duration_weight = str2double(clusterParameters{7});

        % Get the data
        [ClusteringData, clustAssign] = CreateClusteringData(handles, 'forClustering', true);

        %% Extract features
        ReshapedX   = cell2mat(cellfun(@(x) imresize(x',[1 9]) ,ClusteringData.xFreq,'UniformOutput',0));
        slope       = diff(ReshapedX,1,2);
        slope       = zscore(slope);
        freq        = cell2mat(cellfun(@(x) imresize(x',[1 8]) ,ClusteringData.xFreq,'UniformOutput',0));
        freq        = zscore(freq);
        duration    = repmat(ClusteringData.Duration,[1 8]);
        duration    = zscore(duration);

        data = [
            freq     .*  freq_weight,...
            slope    .*  slope_weight,...
            duration .*  duration_weight,...
            ];

    case 'Variational autoencoder embeddings'
        clusterParameters = inputdlg({'Images height (pixels)','Image width (pixels)','Perplexity (if using t-SNE, ignore for UMAP)','Max number of calls to plot (set to 0 to plot everything)'},...
            'Choose cluster parameters:',1,...
            {'6000','6000','30','2000'});
        if isempty(clusterParameters); return; end

        FromExisting = questdlg('From existing model?','Cluster','Yes','No','No');
        if isempty(FromExisting); return; end
        switch FromExisting % Load Model
            case 'No'
                [encoderNet, decoderNet, options, ClusteringData] = create_VAE_model(handles);
                [FileName, PathName] = uiputfile(fullfile(handles.data.squeakfolder, 'Clustering Models', 'Variational Autoencoder Model.mat'), 'Save clustering model (optional)');
                if ~isnumeric(FileName) % Save the new model
                    save(fullfile(PathName, FileName), 'encoderNet', 'decoderNet', 'options');
                end
            case 'Yes' % Load the VAE model
                [FileName, PathName] = uigetfile(fullfile(handles.data.squeakfolder, 'Clustering Models', '*.mat'), 'Select VAE model');
                if isnumeric(FileName);return;end
                load(fullfile(PathName,FileName),'encoderNet','options');
                [ClusteringData, clustAssign] = CreateClusteringData(handles, 'spectrogramOptions', options.spectrogram, 'scale_duration', options.maxDuration, 'freqRange', options.freqRange, 'forClustering', true);
                if isempty(ClusteringData); return; end
        end
        data = extract_VAE_embeddings(encoderNet, options, ClusteringData);

end

imsize = str2double(clusterParameters(1:2))';
perplexity = str2double(clusterParameters{3});
NumberOfCalls = str2double(clusterParameters{4});

% Plot all calls if number of calls is set to 0
if NumberOfCalls == 0
    NumberOfCalls = size(ClusteringData,1);
end
NumberOfCalls = min(size(ClusteringData,1), NumberOfCalls);

%% Get parameters

% Calculate embeddings
rng default;

% Run embedding
switch embeddingType
    case 't-SNE'
        embed = tsne(data,'Verbose',1,'Perplexity',perplexity);
    case 'UMAP'
        embed = run_umap(data);
end
% Rescale values between 0 and 1
embed = (embed - min(embed)) ./ (max(embed)-min(embed));

ClusteringData.embedY = 1-embed(:,2); % flip Y coordinates so the images looks like the UMAP figure
ClusteringData.embedX = embed(:,1);

switch colorType
    case 'Frequency'
        minfreq = prctile(ClusteringData.MinFreq, 1);
        maxfreq = prctile(ClusteringData.MinFreq + ClusteringData.Bandwidth, 99);
        ColorData = HSLuv_to_RGB(256, 'H',  [0 270], 'S', 100, 'L', 75); % Make a color map for each category
        ColorData = reshape(ColorData,size(ColorData,1),1,size(ColorData,2));
    case 'Cluster'
        [clustAssignID, cName] = findgroups(clustAssign); % Convert categories into numbers
        ClusteringData.Cluster = clustAssignID; % Append the category number to clustering data

        % make it so that adjacent clusters are generally different colors.
        % it turns out that this isn't trivial, so try 200 different color
        % orders and use the best one.
        embedings=[ClusteringData.embedY,  ClusteringData.embedX];
        z=unique(clustAssignID);
        for c=1:height(z)
            idx=clustAssignID==z(c);
            clusterCentroids(c,:)=mean(embedings(idx,:));
        end
        % = splitapply(@mean,[ClusteringData.embedY,  ClusteringData.embedX], clustAssignID);
        clusterCentroids = pdist2(clusterCentroids, clusterCentroids);
        hueAngle = [
            sin(linspace(0,2*pi, length(cName)))
            cos(linspace(0,2*pi, length(cName)))]';
        hueDistance = pdist2(hueAngle, hueAngle, 'cosine');
        colorSeperation = zeros(200,1);
        colorOrders = zeros(200,length(cName));
        for i = 1:200
            colorOrders(i,:) = randperm(length(cName)); % Randomize the order of the colors
            colorSeperation(i) = sum((1-clusterCentroids).* hueDistance(colorOrders(i,:),colorOrders(i,:)), 'all');
        end
        [~, idx] = max(colorSeperation)
        colorOrder = colorOrders(idx,:)

        hueRange = [0 360-360/length(cName)];
        cMap = HSLuv_to_RGB(length(cName), 'H', hueRange, 'S', 100, 'L', 75); % Make a color map for each category
        cMap = cMap(colorOrder,:);
        figure('Color','w') % Display the colors
        h = image(reshape(cMap,[],1,3));
        yticklabels(h.Parent, cellstr(cName));
        yticks(h.Parent,1:length(cName));
end

%% Create the image
im = zeros([imsize+padding*2,3],'uint8');

% Only plot the X number of calls
calls2plot = randsample(size(ClusteringData, 1), NumberOfCalls, false);

for i = calls2plot'
    call = ClusteringData(i, :);
    % Get x and y coordinates to place with image
    iy = imsize(1) * call.embedY;
    iy = iy:iy + size(call.Spectrogram{:},1) - 1;
    iy = round(iy - mean(size(call.Spectrogram{:},1))) + padding;
    ix = imsize(2) * call.embedX;
    ix = ix:ix + size(call.Spectrogram{:},2) - 1;
    ix = round(ix - mean(size(call.Spectrogram{:},2))) + padding;

    % Either use the call pitch or the cluster id to apply a color mask
    switch colorType
        case 'Frequency'
            % Interpolate the color data
            freqdata = linspace(call.MinFreq, call.MinFreq + call.Bandwidth, size(call.Spectrogram{:},1));
            colorMask = interp1(linspace(minfreq, maxfreq, size(ColorData,1)), ColorData, freqdata, 'nearest', 'extrap');
        case 'Cluster'
            colorMask = reshape(cMap(call.Cluster,:),1,1,3);
    end
    im(iy,ix,:) = max(im(iy,ix,:),uint8(single(call.Spectrogram{:}) .* colorMask - blackLevel));
end

% Crop the image at the first and last non-empty pixel
[y1,~] = find(max(max(im,[],3),[],2),1,'first');
[y2,~] = find(max(max(im,[],3),[],2),1,'last');

[~,x1] = find(max(max(im,[],3),[],1),1,'first');
[~,x2] = find(max(max(im,[],3),[],1),1,'last');

[fname,fpath] = uiputfile({'*.jpg'; '*.png'},'Save image', 'embeddings');
imwrite(im2uint8(im(y1:y2,x1:x2,:)),fullfile(fpath,fname))

% Open the image in a file manager
if ispc % Open the file in windows explorer
    system(['explorer.exe /select,"' fullfile(fpath,fname) '"']);
    % winopen(fullfile(fpath,fname));
elseif ismac % Open the file in finder - UNTESTED (I don't have a mac)
    system(['open -a Finder "' fullfile(fpath,fname) '"']);
elseif isunix % Open the file in linux file manager - UNTESTED
    system(['xdg-open "' fullfile(fpath) '"']); % open folder in file manager
    % system(['xdg-open "' fullfile(fpath, fname) '"']); % open image in default viewer
end

author: DrCoffey

repo: DeepSqueak

save_path: github-repos/MATLAB/DrCoffey-DeepSqueak


text:
classdef MatNatModality
    % MatNatModality Enumeration representing a modality for imaging data
    %
    %     Licence
    %     -------
    %     Part of MatNat. https://github.com/tomdoel/matnat
    %     Author: Tom Doel, 2015.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        Name
        XnatType
    end

    methods
        function obj = MatNatModality(name, xnatType)
            obj.Name = name;
            obj.XnatType = xnatType;
        end
    end

    enumeration
        MR ('MR', 'xnat:mrScanData')
    end

    properties
    end

    methods (Static)
        function modality = getModalityFromXnatString(xnatType)
            % Get the modality type from the Xnat type string

            allEnums = enumeration('MatNatModality');
            for enum = allEnums
                if strcmp(xnatType, enum.XnatType)
                    modality = enum;
                    return;
                end
            end
            modality = [];
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function Y = au_qmr_ilu_s(tr,X,i)
%
%  Solves shifted linear systems with the matrix A or its transpose A' by
%  the QMR iteration with ILU preconditioning:
%
%  for tr = 'N':
%
%    Y = inv(A+p(i)*I)*X;
%
%  for tr = 'T':
%
%    Y = inv(A.'+p(i)*I)*X.
%
%  The incomplete LU factors of the ILU preconditioner for A+p(i)*I are
%  provided as global data if mc has been set to 'C' when calling
%  'au_qmr_ilu_m_i'. In this case, this data must be generated by calling
%  'au_qmr_ilu_s_i' before calling this routine! If mc has been set to 'M',
%  the incomplete factorizations are computed within this routine to save
%  memory. However, the computational cost is increased in this case.
%
%  Calling sequence:
%
%    Y = au_qmr_ilu_s(tr,X,i)
%
%  Input:
%
%    tr        (= 'N' or 'T') determines whether shifted systems with
%              A or A' should be solved;
%    X         a matrix of proper size;
%    i         the index of the shift parameter.
%
%  Output:
%
%    Y         the resulting solution matrix.
%
%
%   LYAPACK 1.0 (Thilo Penzl, August 1999)

if nargin~=3
  error('Wrong number of input arguments.');
end

global LP_A LP_P LP_MC LP_MAX_IT_QMR LP_TOL_QMR LP_TOL_ILU LP_INFO_QMR

if ~length(LP_A) | ~length(LP_P) | ~length(LP_MC) | ~length(LP_MAX_IT_QMR) ...
  | ~length(LP_TOL_QMR) | ~length(LP_TOL_ILU) | ~length(LP_INFO_QMR)
  error('This routine needs global data which must be generated by calling ''au_qmr_ilu_m_i'' first.');
end

[n,m] = size(X);

if LP_MC=='C'

  eval(lp_e( 'global LP_L',i,' LP_U',i ));

  is_init1 = eval(lp_e( 'length(LP_L',i,')' ));
  is_init2 = eval(lp_e( 'length(LP_U',i,')' ));
  if ~is_init1 | ~is_init2
    error('This routine needs global data which must be generated by calling ''au_qmr_ilu_s_i'' first.');
  end

else

  eval(lp_e( '[LP_L',i,',LP_U',i,'] = luinc(LP_A+LP_P(i)*speye(n),LP_TOL_ILU);' ));

%
%  This shows the amount of fill-in created by the ILU preconditioner:
%
  if LP_INFO_QMR >= 5
    nnA = nnz(LP_A);
    eval(lp_e( 'nnLU = nnz(LP_L',i,')+nnz(LP_U',i,');'  ));
    disp(sprintf('ILU-Ratio (nnz(L)+nnz(U))/nnz(A) = %6g',nnLU/nnA));
  end

end

Y = zeros(n,m);

if tr=='N'
  for j = 1:m
    eval(lp_e(  '[Y(:,j),fl,rr] = qmr(LP_A+LP_P(i)*speye(n),X(:,j),LP_TOL_QMR,LP_MAX_IT_QMR,LP_L',i,',LP_U',i,');'  ));
    if LP_INFO_QMR & fl
      disp('WARNING in ''au_qmr_ilu_s'': QMR terminated not properly!');
      disp(sprintf('  ( normalized precond. QMR residual = %6g )',rr));
      rr0 = norm(LP_A*Y(:,j)+LP_P(i)*Y(:,j)-X(:,j))/norm(X(:,j));
      disp(sprintf('  ( normalized UNprecond. QMR residual = %6g )',rr0));
    end
  end
elseif tr=='T'
  for j = 1:m
    eval(lp_e(  '[Y(:,j),fl,rr] = qmr((LP_A.''+LP_P(i)*speye(n)),X(:,j),LP_TOL_QMR,LP_MAX_IT_QMR,LP_U',i,'.'',LP_L',i,'.'');'  ));
    if LP_INFO_QMR & fl
      disp('WARNING in ''au_qmr_ilu_s'': QMR terminated not properly!');
      disp(sprintf('  ( normalized QMR residual = %6g )',rr));
      rr0 = norm(LP_A.'*Y(:,j)+LP_P(i)*Y(:,j)-X(:,j))/norm(X(:,j));
      disp(sprintf('  ( normalized UNprecond. QMR residual = %6g )',rr0));
    end
  end
else
  error('tp must be either ''N'' or ''T''.');
end

author: Sable

repo: mcbench-benchmarks


text:
function r = sin(a)
%SIN          Hessian (elementwise) sine
%

% written  04/04/04     S.M. Rump
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  K = prod(size(a.x));
  if K==1                   % scalar hessian

    r.x = sin(a.x);
    cosax = cos(a.x);
    r.dx = cosax * a.dx;
    r.hx = cosax * a.hx - reshape( ((0.5*r.x)*a.dx) * a.dx.' , size(a.hx) );

  else                      % matrix hessian

    N = getappdata(0,'INTLAB_HESSIAN_NUMVAR');
    N2 = N^2;

    r.x = sin(a.x);
    if issparse(a.hx)               % input sparse

      ax = cos(full(a.x(:)));
      sizeax = length(ax);
      [ia,ja,sa] = find(a.dx);
      % check for emptyness: cures Matlab bug
      % a=sparse([],[],[],2,1), [i,j,s]=find(a), s(i).*s(:)  yields error
      if isempty(ia)
        r.dx = sparse([],[],[],N,sizeax);
        r.hx = sparse([],[],[],N2,sizeax);
      else
        adx1 = (-0.5)*r.x(:);
        if isa(a.x,'intval')          % sparse intval
          rdx = times(ax(ja),sa(:),0);
          adx1 = times(adx1(ja),sa(:),0);
          if rdx.complex
            r.dx = intval( sparse(ia,ja,rdx.mid,N,sizeax) , sparse(ia,ja,rdx.rad,N,sizeax) , 'midrad' );
          else
            r.dx = intval( sparse(ia,ja,rdx.inf,N,sizeax) , sparse(ia,ja,rdx.sup,N,sizeax) , 'infsup' );
          end
          if adx1.complex
            adx1 = intval( sparse(ia,ja,adx1.mid,N,sizeax) , sparse(ia,ja,adx1.rad,N,sizeax) , 'midrad' );
          else
            adx1 = intval( sparse(ia,ja,adx1.inf,N,sizeax) , sparse(ia,ja,adx1.sup,N,sizeax) , 'infsup' );
          end
        else                          % sparse point
          r.dx = sparse(ia,ja,ax(ja).*sa(:),N,sizeax);
          adx1 = sparse(ia,ja,adx1(ja).*sa(:),N,sizeax);
        end
        r.hx = adx2rhx(N,sizeax,adx1,a.dx);
      end
      [ia,ja,sa] = find(a.hx);        % sparse point or intval
      % check for emptyness: cures Matlab bug
      % a=sparse([],[],[],2,1), [i,j,s]=find(a), s(i).*s(:)  yields error
      if ~isempty(ia)
        if isa(a.x,'intval')
          rhx = times(ax(ja),sa(:),0);
          if rhx.complex
            r.hx = r.hx + intval( sparse(ia,ja,rhx.mid,N2,sizeax) , sparse(ia,ja,rhx.rad,N2,sizeax) , 'midrad' );
          else
            r.hx = r.hx + intval( sparse(ia,ja,rhx.inf,N2,sizeax) , sparse(ia,ja,rhx.sup,N2,sizeax) , 'infsup' );
          end
        else
          r.hx = r.hx + sparse(ia,ja,ax(ja).*sa(:),N2,sizeax);
        end
      end

    else                            % input full

      r.x = sin(a.x);
      rx = r.x(:).';
      cosax = cos(a.x(:).');
      cosax = cosax(ones(N*N,1),:);
      r.dx = a.dx .* cosax(1:N,:);
      adx = repmat(0.5*rx,N,1) .* a.dx;
      r.hx = a.hx .* cosax - adx(repmat(1:N,N,1),:) .* a.dx(repmat(1:N,1,N),:);

    end

  end

  r = class(r,'hessian');

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
%% boxconstrain
% Below is a demonstration of the features of the |boxconstrain| function

%% Syntax
% |[xb]=boxconstrain(x,lb,ub,m,s);|

%% Description
% The |boxconstrain| function can be used to constrain parameters from
% [-inf,inf] to the range [lb ub]. The sigmoidal tanh function is used to
% do this mapping. The middle of the sigmoid, denoted by the parameter m
% need not be (lb+ub)/2. If x=m xb=m.
% The optional parameter s sets the slope of the sigmoid at
% x=m, the default slope is 1.

%% Examples

clear; close all; clc;

%% Example: Constraining parameters

m=5; %The middle
lb=0; %The lower bound
ub=10; %The upper bound
x=linspace(m-20,m+20,1000); % x paremeters for plotting
[y]=boxconstrain(x,lb,ub,m);

%%

cFigure; hold on;
xlabel('x'); ylabel('Constrained x');
h1=plot([min(x) max(x)],[ub ub],'b--','LineWidth',3);
h2=plot([min(x) max(x)],[lb lb],'r--','LineWidth',3);
h3=plot(m,m,'k.','MarkerSize',50);
h4=plot(x,y,'k-','MarkerSize',25,'LineWidth',3);

legend([h1;h2;h3;h4],{'Upper limit','Lower limit','Middle','Constrained'},'Location','Northwest');

axis tight;
set(gca,'FontSize',25);
grid on; box on;
drawnow;

%% Example: Assymetric mappings

m=1; %The middle
lb=0; %The lower bound
ub=10; %The upper bound

x=linspace(m-20,m+20,1000); % x paremeters for plotting

[y]=boxconstrain(x,lb,ub,m);

%%

cFigure; hold on;
xlabel('x'); ylabel('Constrained x');
h1=plot([min(x) max(x)],[ub ub],'b--','LineWidth',3);
h2=plot([min(x) max(x)],[lb lb],'r--','LineWidth',3);
h3=plot(m,m,'k.','MarkerSize',50);
h4=plot(x,y,'k-','MarkerSize',25,'LineWidth',3);

legend([h1;h2;h3;h4],{'Upper limit','Lower limit','Middle','Constrained'},'Location','Northwest');

axis tight;
set(gca,'FontSize',25);
grid on; box on;
drawnow;

%% Example: Custom slop

m=5; %The middle
lb=0; %The lower bound
ub=10; %The upper bound
s=0.5; %Slope at the middle

x=linspace(m-20,m+20,1000); % x paremeters for plotting

[y]=boxconstrain(x,lb,ub,m,s);

%%

cFigure; hold on;
xlabel('x'); ylabel('Constrained x');
h1=plot([min(x) max(x)],[ub ub],'b--','LineWidth',3);
h2=plot([min(x) max(x)],[lb lb],'r--','LineWidth',3);
h3=plot(m,m,'k.','MarkerSize',50);
h4=plot(x,y,'k-','MarkerSize',25,'LineWidth',3);

legend([h1;h2;h3;h4],{'Upper limit','Lower limit','Middle','Constrained'},'Location','Northwest');

axis tight;
set(gca,'FontSize',25);
grid on; box on;
drawnow;

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

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%  INTERNAL FUCNTION
%
%

author: jmaih


text:
function planC = CERRImport(optName, archiveDir, saveFile)
%function planC = CERRImport(optName,archiveDir,saveDir)
%
%Description: Wrapper file to import AAPM/RTOG directories.  The main work is done
%by importRTOGDir.  If called without arguments, the user selects an archive file and
%an options file via GUI dialogs, and the name of the file to save the
%CERR plan to is input via the file save GUI.  Can be called with 1, 2, or 3
%arguments.
%
%Inputs:
%optName -- is the name of an options file (an example file is initCERRImport).
%archiveDir -- is the path of the archive directory, with trailing slash.
%saveFile -- the name, including path, of a file to which to save the CERR plan.
%
%Output: None, except the saved file.
%
%Globals:  None.
%
%Author: J. O. Deasy, deasy@radonc.wustl.edu
%
%Reference:
%W. Harms, Specifications for Tape/Network Format for Exchange of Treatment Planning Information,
%version 3.22., RTOG 3D QA Center, (http://rtog3dqa.wustl.edu), 1997.
%
%Last Modified:  14 Feb 02, by JOD.
%
%Version:  1.x. (Incremented by 0.001 for every modification of this file or a called file.)
%
%Copyright:
%This software is copyright J. O. Deasy and Washington University in St Louis.
%A free license is granted to use or modify but only for non-commercial non-clinical use.
%Any user-modified software must retain the original copyright and notice of changes if redistributed.
%Any user-contributed software retains the copyright terms of the contributor.

% copyright (c) 2001-2006, Washington University in St. Louis.
% Permission is granted to use or modify only for non-commercial,
% non-treatment-decision applications, and further only if this header is
% not removed from any file. No warranty is expressed or implied for any
% use whatever: use at your own risk.  Users can request use of CERR for
% institutional review board-approved protocols.  Commercial users can
% request a license.  Contact Joe Deasy for more information
% (radonc.wustl.edu@jdeasy, reversed).

fname1 = '';
%Interactive or not:
if nargin == 0
    [fname1, pathname] = uigetfile('*.*','Select 0000 file to begin importing.');
    archiveDir = pathname;

    [fname, pathname] = uigetfile('*.m','Select CERR options .m file');
    optName = [pathname fname];

    %   % Would the user like to store the plan as a compressed file?
    %   % The following loop investigates this.
    %   % Compression (gz) is made through GNUZIP (freeware)
    %   % Default compression options are 'verbose' and 'best'
    %   gzip 1.2.4 Win32 (02 Dec 97)
    %   usage: gzip [-acdfhlLnNrtvV19] [-S suffix] [file ...]
    %       -a --ascii       ascii text; convert end-of-lines using local conventions
    %       -c --stdout      write on standard output, keep original files unchanged
    %       -d --decompress  decompress
    %       -f --force       force overwrite of output file and compress links
    %       -h --help        give this help
    %       -l --list        list compressed file contents
    %       -L --license     display software license
    %       -n --no-name     do not save or restore the original name and time stamp
    %       -N --name        save or restore the original name and time stamp
    %       -q --quiet       suppress all warnings
    %       -r --recursive   operate recursively on directories
    %       -S .suf  --suffix .suf     use suffix .suf on compressed files
    %       -t --test        test compressed file integrity
    %       -v --verbose     verbose mode
    %       -V --version     display version number
    %       -1 --fast        compress faster
    %       -9 --best        compress better
    %       file...          files to (de)compress. If none given, use standard input.
    % Logic for compressing files:
    %   Matlab binaries will be created first in all cases; subsequently,
    %       they will be replaced by the corresponding compressed binary IF
    %       the user has selected this option in 'CERROptions.m'

    %fid = fopen(optName,'r');
    %while(1)
    %    tline = fgetl(fid);
    %    if findstr((tline), 'optS.saveZip'), break, end
    %end
    %fclose(fid);
    %index = findstr(tline, '=');
    %toZipOrNot = lower(midstring(tline, (index + 1), length(tline) - 1));
    %if findstr(toZipOrNot, 'yes')
    %    toZipOrNot = 'yes';
    %elseif findstr(toZipOrNot, 'no')
    %    toZipOrNot = 'no';
    %else
    warndlg('Output file compression not set: defaulting to full.', ...
        'CERR Configuration');
end
optS = opts4Exe(optName);

if strcmp(optS.saveZip,'yes')

    %end        % APA commented this end (may be placed incorrectly, gave
    %an error during crosslisting)

    saveFlag = 'ui';

elseif nargin == 1 %Assume options given
    [fname1, pathname] = uigetfile('*.*','Select aapm0000 file to begin importing.');
    archiveDir = pathname;
    fid = fopen(optName,'r');
    while(1)
        tline = fgetl(fid);
        if strfind((tline), 'optS.saveZip'), break, end
    end
    fclose(fid);
    index = strfind(tline, '=');
    toZipOrNot = lower(midstring(tline, (index + 1), length(tline) - 1));
    if strfind(toZipOrNot, 'yes')
        toZipOrNot = 'yes';
    elseif strfind(toZipOrNot, 'no')
        toZipOrNot = 'no';
    else
        warndlg('Output file compression not set: defaulting to full.', ...
            'CERR Configuration');
    end
    optS = opts4Exe(optName);

    saveFlag = 'ui';

elseif nargin == 2
    fid = fopen(optName,'r');
    while(1)
        tline = fgetl(fid);
        if strfind((tline), 'optS.saveZip'), break, end
    end
    fclose(fid);
    index = strfind(tline, '=');
    toZipOrNot = lower(midstring(tline, (index + 1), length(tline) - 1));
    if strfind(toZipOrNot, 'yes')
        toZipOrNot = 'yes';
    elseif strfind(toZipOrNot, 'no')
        toZipOrNot = 'no';
    else
        warndlg('Output file compression not set: defaulting to full.', ...
            'CERR Configuration');
    end
    optS = opts4Exe(optName);

    saveFlag = 'ui';

elseif nargin == 3
    fid = fopen(optName,'r');
    while(1)
        tline = fgetl(fid);
        if strfind((tline), 'optS.saveZip'), break, end
    end
    fclose(fid);
    index = strfind(tline, '=');
    toZipOrNot = lower(midstring(tline, (index + 1), length(tline) - 1));
    if strfind(toZipOrNot, 'yes')
        toZipOrNot = 'yes';
    elseif strfind(toZipOrNot, 'no')
        toZipOrNot = 'no';
    else
        warndlg('Output file compression not set: defaulting to full.', ...
            'CERR Configuration');
    end
    optS = opts4Exe(optName);
    if ~strcmp(saveFile,'no')
        saveFlag = 'spec';
    else
        saveFlag = 'no';
    end
    fname1 = 'aapm0000';

end

%Call the import tool:
planC = importRTOGDir(optS, archiveDir, fname1)

switch saveFlag

    case 'ui'

        [fname pname] = uiputfile( '*.mat', ...
            'Save the plan data as:');
        saveFile = [pname fname];
        if ~strcmp(midstring(saveFile, ...
                length(saveFile)-3, length(saveFile)), '.mat')
            % '*.mat' file type not typed by the user
            % this situation can create a file error when compressing
            %       will add the '*.mat' suffix prior to archive being written
            saveFile = strcat(saveFile, '.mat');
        end
        save(saveFile,'planC')
        if strcmp(toZipOrNot, 'yes') % User wants to compress output file
            outstr = gnuCERRCompression(saveFile,'compress');
            h = dos(outstr);
            if (h) % Function returned as zero when compression complete
                errordlg('File not found by compression routine.', 'CERR compression');
            end
        end

    case 'spec'

        if isempty(saveFile) %save within the subdirectory
            saveFile = [archiveDir 'planC.mat'];
            if strcmp(toZipOrNot, 'yes') % User wants to compress output file
                outstr = gnuCERRCompression(saveFile,'compress');
                h = dos(outstr);
                if (h) % Function returned as zero when compression complete
                    errordlg('File not found by compression routine.', 'CERR compression');
                end
            end
        end

        try
            save(saveFile,'planC')
            if strcmp(toZipOrNot, 'yes') % User wants to compress output file
                outstr = gnuCERRCompression(saveFile,'compress');
                h = dos(outstr);
                if (h) % Function returned as zero when compression complete
                    errordlg('File not found by compression routine.', 'CERR compression');
                end
            end

        catch
            warndlg('That file name was not valid.', 'CERR file save')
            [fname pname] = uiputfile('*.mat','Save the plan data as:');
            saveFile = [pname fname];
            save(saveFile,'planC')
            if strcmp(toZipOrNot, 'yes') % User wants to compress output file
                outstr = gnuCERRCompression(saveFile,'compress');
                h = dos(outstr);
                if (h) % Function returned as zero when compression complete
                    errordlg('File not found by compression routine.', 'CERR compression');
                end
            end

        end
    case 'no'

end

% Unzip routines
% outstr = ['gzip.exe -vd ', fmat]
% dos(outstr)

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function angle_axis = RotationMatrix2AngleAxis(R)

% The conversion of a rotation matrix to the angle-axis form is
% numerically problematic when then rotation angle is close to zero
% or to Pi. The following implementation detects when these two cases
% occurs and deals with them by taking code paths that are guaranteed
% to not perform division by a small number.

% x = k * 2 * sin(theta), where k is the axis of rotation.
angle_axis(1) = R(1+2, 1+1) - R(1+1, 1+2);
angle_axis(2) = R(1+0, 1+2) - R(1+2, 1+0);
angle_axis(3) = R(1+1, 1+0) - R(1+0, 1+1);

% Since the right hand side may give numbers just above 1.0 or
% below -1.0 leading to atan misbehaving, we threshold.
costheta = min(max((R(1+0, 1+0) + R(1+1, 1+1) + R(1+2, 1+2) - 1.0) / 2.0,  -1.0), 1.0);

% sqrt is guaranteed to give non-negative results, so we only
% threshold above.
sintheta = min(sqrt(angle_axis(1) * angle_axis(1) + angle_axis(2) * angle_axis(2) + angle_axis(3) * angle_axis(3)) / 2.0, 1.0);

% Use the arctan2 to get the right sign on theta
theta = atan2(sintheta, costheta);

% Case 1: sin(theta) is large enough, so dividing by it is not a
% problem. We do not use abs here, because while jets.h imports
% std::abs into the namespace, here in this file, abs resolves to
% the int version of the function, which returns zero always.
%
% We use a threshold much larger then the machine epsilon, because
% if sin(theta) is small, not only do we risk overflow but even if
% that does not occur, just dividing by a small number will result
% in numerical garbage. So we play it safe.
kThreshold = 1e-12;
if ((sintheta > kThreshold) || (sintheta < -kThreshold))
    r = theta / (2.0 * sintheta);
    angle_axis = angle_axis * r;
    return;
end

% Case 2: theta ~ 0, means sin(theta) ~ theta to a good
% approximation.
if (costheta > 0.0)
    angle_axis = angle_axis * 0.5;
    return;
end

% Case 3: theta ~ pi, this is the hard case. Since theta is large,
% and sin(theta) is small. Dividing by theta by sin(theta) will
% either give an overflow or worse still numerically meaningless
% results. Thus we use an alternate more complicated formula
% here.

% Since cos(theta) is negative, division by (1-cos(theta)) cannot
% overflow.
inv_one_minus_costheta = 1.0 / (1.0 - costheta);

% We now compute the absolute value of coordinates of the axis
% vector using the diagonal entries of R. To resolve the sign of
% these entries, we compare the sign of angle_axis[i]*sin(theta)
% with the sign of sin(theta). If they are the same, then
% angle_axis[i] should be positive, otherwise negative.

for i=1:3
    angle_axis(i) = theta * sqrt((R(i, i) - costheta) * inv_one_minus_costheta);
    if (((sintheta < 0.0) && (angle_axis(i) > 0.0)) || ((sintheta > 0.0) && (angle_axis(i) < 0.0)))
        angle_axis(i) = -angle_axis(i);
    end
end

author: jianxiongxiao

repo: ProfXkit


text:
function s = rdivide(f, g)
%./   Divide DELTAFUNS with DELTAFUNS
%   RDIVIDE(F, G) computes the pointwise division F./G. The operation is only
%   defined if G does not have any delta functions and G has no roots.
%
% See also LDIVIDE, TIMES.
%
% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Case of empty arguments:
if ( isempty(f) || isempty(g) )
    % Return an empty DELTAFUN:
    s = deltafun();
    return
end

% Check if inputs are other than DELTAFUNs, CLASSICFUNs or doubles.
if ( (~isa(f, 'deltafun') && ~isa(f, 'classicfun') && ~isa(f, 'double')) || ...
     (~isa(g, 'deltafun') && ~isa(g, 'classicfun') && ~isa(g, 'double')) )
    error('CHEBFUN:DELTAFUN:rdivide:rdivide', ...
        'Input can only be a DELTAFUN, a CLASSICFUN or a double')
end

%% Reciprocal: ( CLASSICFUN or DOUBLE ). / DELTAFUN
if ( isa(f, 'double') || isa( f, 'classicfun') && isa(g, 'deltafun') )
    if ( anyDelta(g) )
        error('CHEBFUN:DELTAFUN:rdivide:rdivide', ...
            'Division by delta functions is not defined.');
    end
    % A smooth function is returned in this case:
    s = f ./ g.funPart;
    return
end

%% DELTAFUN./DOUBLE
if ( isa(f, 'deltafun') && isa(g,'double') )
    % Copy the DELTAFUN in the output:
    s = f;
    % Divide the funPart with the double g:
    s.funPart = 1/g * f.funPart;
    % Update the delta functions:
    s.deltaMag = 1/g * s.deltaMag;
    if ( ~anyDelta(s) )
        s = s.funPart;
    end
    return
end

%% DELTAFUN./CLASSICFUN
if ( isa(f, 'deltafun') && isa(g, 'classicfun') )
    % Take reciprocal and make a DELTAFUN:
    g = deltafun( 1./g, [], [] );
    % Now, multiply:
    s = f .* g;
    return
end

%% DELTAFUN./DELTAFUN
if ( isa(f, 'deltafun') && isa(g, 'deltafun') )
    if ( anyDelta(g) )
        error('CHEBFUN:DELTAFUN:rdivide:rdivide', ...
            'Division by delta functions is not defined.');
    end
    % Take reciprocal of the funPart only and make a DELTAFUN:
    g = deltafun( 1 ./ g.funPart, [], []);
    % Now, multiply:
    s = f .* g;
    return
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function sp=findConmatFromEdges(nodeList,edgeList)
% function connectionMatrix=findConmatFromEdges(edgeList)
% Generates a connection matrix from a list of edges
% The returned connection matrix is max(edgeList(:)) square
%
% nodeList is an 8xn array
% nodes:  8xN array of
%    Nx(x,y,z,num_edges,edge_offset,layer,dist,pqindex).
% edges:  1xM array of node indices.  The edge_offset of
%    each node points into the starting location of its set
%    of edges.

nNodes=size(nodeList,2);

mxNode=max(edgeList(:));
sp=spalloc(nNodes,nNodes,fix(nNodes/40));

for t=1:nNodes

	nEdges=nodeList(4,t);
	offset=nodeList(5,t);

	theseNodes=edgeList(offset:(offset+nEdges-1));
	sp(t,theseNodes)=1;
end

sp=((sp+sp')~=0); % Matrix is symmetrical

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function triplot
    % triplot plots a triangular diagram
    % points on it are plotted by tripts
    % written by Gerry Middleton, 1995
    xscale = 2/sqrt(3);
    x = [0 xscale/2 xscale 0];		 % define coordinates of apices
    y = [-1 0 -1 -1];
    plot(x,y,'k');
    axis('equal');
    axis off;						 % turn off rectangular axes
    hold on;						 % hold, to permit adding to plot later
    for i = 1:9						 % plot dotted lines parallel to sides,
        x1 = i*xscale/20;			 % at 10 percent intervals
        x2 = (20-i)*xscale/20;
        x3 = i*xscale/10;
        y1 = -1 + i/10;
        y2 = -i/10;
        xx = [x1 x2];
        yy = [y1 y1];
        xx2 = [x1 x3];
        yy2 = [y1 -1];
        xx3 = [x3 (xscale/2 + x1)];
        yy3 = [-1 y2];
        plot(xx,yy,'k:',xx2,yy2,'k:',xx3,yy3,'k:');
    end
    set(gcf,'color','w')

author: CelsoReyes


text:
function varargout = erfcinv(varargin)

switch class(varargin{1})

    case 'sdpvar'
        varargout{1} = InstantiateElementWise(mfilename,varargin{:});

    case 'char'

        operator = CreateBasicOperator('decreasing','callback');
        operator.inverse = @(x)(erfc(x));
        operator.derivative = @(x)-1./(exp(-((erfcinv(x))).^2)*2/sqrt(pi));
        operator.inflection = [1 -1];
        operator.domain = [0 2];

        varargout{1} = [];
        varargout{2} = operator;
        varargout{3} = varargin{3};

    otherwise
        error(['SDPVAR/' upper(mfilename) ' called with weird argument']);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


function [dy edy dx edx y ey x ex tmpw tmph] = pad(total_boxes,w,h)
	%compute the padding coordinates (pad the bounding boxes to square)
	tmpw=total_boxes(:,3)-total_boxes(:,1)+1;
	tmph=total_boxes(:,4)-total_boxes(:,2)+1;
	numbox=size(total_boxes,1);

    dx=ones(numbox,1);dy=ones(numbox,1);
	edx=tmpw;edy=tmph;

	x=total_boxes(:,1);y=total_boxes(:,2);
	ex=total_boxes(:,3);ey=total_boxes(:,4);

	tmp=find(ex>w);
	edx(tmp)=-ex(tmp)+w+tmpw(tmp);ex(tmp)=w;

	tmp=find(ey>h);
	edy(tmp)=-ey(tmp)+h+tmph(tmp);ey(tmp)=h;

	tmp=find(x<1);
	dx(tmp)=2-x(tmp);x(tmp)=1;

	tmp=find(y<1);
	dy(tmp)=2-y(tmp);y(tmp)=1;
end


text:
function C = glm_determine_noise(Y,E,tr)
% Estimate the noise covariance for each voxel using its time series
%
%    C = glm_determine_noise(Y,E,tr):
%
% This routine is called when we try to whiten the noise in the
% event-related GLM tools.
%
% Given a matrix of time series, Y, and an estimated error for each voxel,
% E, and a frame period, tr, estimate the temporal noise covariance matrix
% for each voxel, C.
%
% Size of each argument:
% Input: Y: # Time Points x # Voxels x # Runs
%        E:
%        tr: single integer
% Output: C: # Predictors x # Predictors x # Runs
%         The # of predictors is the total # of
%         different predictors in a GLM, including
%         DC baseline estimate predictors and, if
%         a selective-averaging GLM is being applied,
%         all nh hemodynamic predictors for each condition.
%         (The # predictors is the # of columns in the
%         design matrix produced by er_createDesMtx.)
%
%
% original code by gb, 11/04:
% updated by ras, 05/05

if nargin<3, help(mfilename); error('Insufficient args.'); end

% params
delayRange = 20; % # of seconds over which to estimate noise
kmax = floor(delayRange / tr);
nFrames = size(E,1);
nSlices = size(E,2);

% use a heuristic to select voxels within the skull:
withinSkullVoxels = glm_choose_voxels(Y);

% estimate the autocorrelation within these
% voxels over the selected delay range (e.g., 20 secs):
B = E(:,withinSkullVoxels); % voxels in brain
for k = 1:kmax
    corr(k,:) = sum(B(1:nFrames-k,:).*B(k+1:nFrames,:))/sum(B.^2);
end
corr = mean(corr,2);

% Determine the parameters alpha and rho using in estimating the
% autocorrelation of residual errors (Re). These parameters will be used in
% the estimation function:
% Re(k) = 1,  k=0;
%         (1-alpha) * rho^k, 0 < |k| <= kmax
%         0, |k| > kmax
% This is equation (5) in the Burock and Dale 2000 paper, and eq. (13) in
% the Greve FS-FAST theory paper.
kmax = length(corr);

[x eval] = fminsearch(@glm_fit_model, [0, 1/2], optimset('TolX',1e-8),kmax,corr);
alpha = x(1);
rho   = x(2);

% finally, compute the noise covariance matrix:
tmp = [1 ((1 - alpha)*(rho.^(1:kmax))) zeros(1,nFrames - kmax - 1)];
C = toeplitz(tmp);

return
% /-----------------------------------------------------/ %

% sum(E(1:nFrames-k,:).*e(k+1:nFrames,:))/sum(E.^2);


text:
function [connectivity] = channelconnectivity(cfg, data)

% CHANNELCONNECTIVIY creates a NxN matrix that describes whether channels
% are connected as neighbours
%
% See also FT_PREPARE_NEIGHBOURS, TRIANGLE2CONNECTIVITY

if (isfield(cfg, 'avgoverchan') && strcmp(cfg.avgoverchan, 'yes')) || isempty(cfg.neighbours)
  if nargin < 2
    nchan = numel(cfg.channel);
  else
    nchan = numel(data.label);
  end
  connectivity = false(nchan,nchan);

else
  % check whether to use channels in the cfg or data
  if nargin < 2 && isfield(cfg, 'channel')
    chans = cfg.channel;
  elseif nargin == 2
    chans = data.label;
  else
    ft_error('either the cfg needs to have both cfg.channel and cfg.neighbours, or a second (data) input argument needs to be specified');
  end

  nchan        = length(chans);
  connectivity = false(nchan,nchan);
  neighbours   = struct(cfg.neighbours);

  try
    % the original loop was the following, which is very slow:
    %   for chan=1:length(cfg.neighbours)
    %     [seld] = match_str(chans, cfg.neighbours(chan).label);
    %     [seln] = match_str(chans, cfg.neighbours(chan).neighblabel);
    %     if isempty(seld)
    %       % this channel was not present in the data
    %       continue;
    %     else
    %       % add the neighbours of this channel to the matrix
    %       connectivity(seld, seln) = true;
    %     end
    %   end
    % the above is still used in case the fast version fails (catch below)
    % the following code should be equivalent:

    if numel(neighbours) < nchan/4
      ft_error('this only works when at least 1/4th of all channels has neighbours defined, or with neighbours = []');
      % FIXME the above error is true (or maybe with <1/5th) but I don't
      % understand why, ES 25-nov-2013
      % the above error will cause the code to fall through to the catch
      % below, so no big worries
    end

    % avoid using match_str inside the loop for performance reasons
    [sel1,sel2] = match_str(chans, {neighbours.label}');

    % make sure we only have the neighbours present in the data
    neighbours = neighbours(sel2);

    % make big list of all neighbour labels...
    allneighb = {neighbours.neighblabel};

    % determine whether neighbours are stored as row or column vector
    siz1 = cellfun(@size, allneighb, repmat({1},size(allneighb)));
    siz2 = cellfun(@size, allneighb, repmat({2},size(allneighb)));
    if sum(siz2) > sum(siz1)
      catdim = 2;
    else
      catdim = 1;
    end

    % ...store dimensionality...
    numAllNeighb = cellfun(@numel, allneighb);
    allneighb = cat(catdim, allneighb{:});
    % ...so that we can do with one call of match_str (with fullout=1)...
    selneighb = match_str(chans, allneighb, 1);
    % ...and then we can reshape the indices again to get the neighbours
    selneighb = mat2cell(selneighb, numAllNeighb);
    for k = 1:numel(sel1)
      % remove zeroes from selneighb (these correspond to channels absent)
      seln = selneighb{k};
      seln = seln(seln > 0);
      % add the neighbours of this channel to the matrix
      connectivity(sel1(k), seln) = true;
    end

  catch
    % the fast version failed, use slow instead

    for chan=1:length(neighbours)
      [seld] = match_str(chans, neighbours(chan).label);
      [seln] = match_str(chans, neighbours(chan).neighblabel);
      if isempty(seld)
        % this channel was not present in the data
        continue;
      else
        % add the neighbours of this channel to the matrix
        connectivity(seld, seln) = true;
      end
    end

  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function s = sem(x,varargin)

%sem - Compute standard error of the mean (SEM).
%
%  USAGE
%
%    s = sem(x)
%
%    x              vector or matrix over which the SEM should be computed
%    dim            dimension to run SEM across (default = 1)

% Copyright (C) 2008-2013 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

if nargin < 1,
  error('Incorrect number of parameters (type ''help sem'' for details).');
end
if nargin == 1
    dim = 1;
else
    dim = varargin{1};
end
if ~isnumeric(x) & ~isdvector(x),
  error('Incorrect input - use vector or matrix (type ''help sem'' for details).');
end

% if any(size(x)==1), x = x(:); end

% s = std(x,[],dim)./sqrt(size(x,1));
s = nanstd(x,[],dim) ./ sqrt(sum(~isnan(x),dim));

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [y] = spm_gx_state_fmri(x,u,P,M)
% Simulated BOLD response and copied state vector
% FORMAT [y] = spm_gx_state_fmri(x,u,P,M)
% y          - BOLD response and copied state vector
%
% x          - state vector     (see spm_fx_fmri)
% P          - Parameter vector (see spm_fx_fmri)
% M          - model specification structure (see spm_nlsi)
%
% The `copied state vector' passes the first hidden variable in each region
% to the output variable y, so that 'neural activities' can be plotted
% by spm_dcm_generate.m
%
% See spm_fx_fmri.m and spm_dcm_generate.m
%__________________________________________________________________________
% Copyright (C) 2011 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_gx_state_fmri.m 6262 2014-11-17 13:47:56Z karl $

y = spm_gx_fmri(x,u,P,M);

% Copy first hidden state (neural activity) from each region
nregions=size(x,1);
for i=1:nregions,
    y=[y;x(i,1)];
end
y=full(y);

author: spm


text:
function pstruct = tapas_hgf_ar1_binary_mab_namep(pvec)
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2017 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

pstruct = struct;

l = (length(pvec)+1)/6;

if l ~= floor(l)
    error('tapas:hgf:UndetNumLevels', 'Cannot determine number of levels');
end

pstruct.mu_0      = pvec(1:l);
pstruct.sa_0      = pvec(l+1:2*l);
pstruct.phi       = pvec(2*l+1:3*l);
pstruct.m         = pvec(3*l+1:4*l);
pstruct.ka        = pvec(4*l+1:5*l-1);
pstruct.om        = pvec(5*l:6*l-1);

return;

author: translationalneuromodeling

repo: tapas


text:
function polygon_triangulate_test ( )

%*****************************************************************************80
%
%% POLYGON_TRIANGULATE_TEST tests the POLYGON_TRIANGULATE library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 May 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'POLYGON_TRIANGULATE_TEST\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Test the POLYGON_TRIANGULATE library.\n' );

  polygon_triangulate_test01 ( );

  polygon_triangulate_test02 ( 'comb' );
  polygon_triangulate_test02 ( 'hand' );
  polygon_triangulate_test02 ( 'i18' );

  polygon_triangulate_test03 ( 'comb' );
  polygon_triangulate_test03 ( 'hand' );
  polygon_triangulate_test03 ( 'i18' );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'POLYGON_TRIANGULATE_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function polygon_triangulate_test01 ( )

%*****************************************************************************80
%
%% POLYGON_TRIANGULATE_TEST01 tests the "comb_10" polygon.
%
%  Discussion:
%
%    There are N-3 triangles in the triangulation.
%
%    For the first N-2 triangles, the first edge listed is always an
%    internal diagonal.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 May 2014
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'POLYGON_TRIANGULATE_TEST01\n' );
  fprintf ( 1, '  Triangulate the comb_10 polygon.\n' );

  x = [ ...
    8; ...
    7; ...
    6; ...
    5; ...
    4; ...
    3; ...
    2; ...
    1; ...
    0; ...
    4];

  y = [ ...
    0; ...
   10; ...
    0; ...
   10; ...
    0; ...
   10; ...
    0; ...
   10; ...
    0; ...
   -2 ];

  n = 10;

  triangles = polygon_triangulate ( n, x, y );

  i4mat_print ( 3, n - 2, triangles, '  Triangles' );

  return
end
function polygon_triangulate_test02 ( prefix )

%*****************************************************************************80
%
%% POLYGON_TRIANGULATE_TEST02 triangulates a polygon described in a file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 May 2014
%
%  Author:
%
%    John Burkardt.
%

%
%  Create filenames.
%
  node_filename = strcat ( prefix, '_nodes.txt' );
  element_filename = strcat ( prefix, '_elements.txt' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'POLYGON_TRIANGULATE_TEST02\n' );
  fprintf ( 1, '  Read polygon coordinates in "%s"\n', node_filename );
%
%  Read the node coordinates.
%
  [ ~, n ] = r8mat_header_read ( node_filename );

  xy = r8mat_data_read ( node_filename, 2, n );
%
%  Get the triangulation.
%
  triangle_num = n - 2;
  triangles = polygon_triangulate ( n, xy(1,:), xy(2,:) );
%
%  Write the triangulation to a file.
%
  i4mat_write ( element_filename, 3, triangle_num, triangles );

  fprintf ( 1, '  Write triangulation to "%s"\n', element_filename );

  return
end
function polygon_triangulate_test03 ( prefix )

%*****************************************************************************80
%
%% POLYGON_TRIANGULATE_TEST03 plots a triangulation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 May 2014
%
%  Author:
%
%    John Burkardt.
%

%
%  Create filenames.
%
  node_filename = strcat ( prefix, '_nodes.txt' );
  plot_filename = strcat ( prefix, '.png' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'POLYGON_TRIANGULATE_TEST03\n' );
  fprintf ( 1, '  Read polygon coordinates in "%s"\n', node_filename );
%
%  Read the node coordinates.
%
  [ ~, n ] = r8mat_header_read ( node_filename );

  xy = r8mat_data_read ( node_filename, 2, n );
%
%  Get the triangulation.
%
  triangle_num = n - 2;
  triangles = polygon_triangulate ( n, xy(1,:), xy(2,:) );
%
%  Plot the triangulation.
%
  figure ( 1 )
  clf
  grid on
  hold on

  plot ( xy(1,1:n), xy(2,1:n), 'k.', 'MarkerSize', 25 )

  for i = 1 : n

    ip1 = mod ( i, n ) + 1;

    line ( [ xy(1,i), xy(1,ip1) ], ...
           [ xy(2,i), xy(2,ip1) ], ...
           'LineWidth', 3, 'Color', 'b' )

  end

  for i = 1 : n - 3
    i1 = triangles(1,i);
    i2 = triangles(2,i);
    line ( [ xy(1,i1), xy(1,i2) ], ...
           [ xy(2,i1), xy(2,i2) ], ...
           'LineWidth', 3, 'Color', 'r' )
  end

  xlabel ( '<-- X -->', 'Fontsize', 16 )
  ylabel ( '<-- Y -->', 'Fontsize', 16 )
  title ( prefix, 'Fontsize', 24 );
  axis equal

  hold off

  print ( '-dpng', plot_filename );
  fprintf ( 1, '  Plot triangulation, save to graphics file "%s"\n', plot_filename );

  return
end
function column_num = file_column_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_COLUMN_COUNT counts the columns in the first line of a file.
%
%  Discussion:
%
%    The file is assumed to be a simple text file.
%
%    Most lines of the file are presumed to consist of COLUMN_NUM words,
%    separated by spaces.  There may also be some blank lines, and some
%    comment lines, which have a "#" in column 1.
%
%    The routine tries to find the first non-comment non-blank line and
%    counts the number of words in that line.
%
%    If all lines are blanks or comments, it goes back and tries to analyze
%    a comment line.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the file.
%
%    Output, integer COLUMN_NUM, the number of columns in the file.
%
  FALSE = 0;
  TRUE = 1;
%
%  Open the file.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_COLUMN_COUNT - Error!' );
  end
%
%  Read one line, but skip blank lines and comment lines.
%  Use FGETL so we drop the newline character!
%
  got_one = FALSE;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( s_len_trim ( line ) == 0 )

    elseif ( line(1) == '#' )

    else
      got_one = TRUE;
      break;
    end

  end

  fclose ( input_unit );

  if ( got_one == FALSE )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Warning!\n' );
    fprintf ( 1, '  The file does not seem to contain any data.\n' );
    column_num = -1;
    return
  end

  column_num = s_word_count ( line );

  return
end
function row_num = file_row_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_ROW_COUNT counts the number of row records in a file.
%
%  Discussion:
%
%    Each input line is a "RECORD".
%
%    The records are divided into three groups:
%
%    * BLANK LINES (nothing but blanks)
%    * COMMENT LINES (begin with a '#')
%    * DATA RECORDS (anything else)
%
%    The value returned by the function is the number of data records.
%
%    By the way, if the MATLAB routine FGETS is used, instead of
%    FGETL, then the variable LINE will include line termination
%    characters, which means that a blank line would not actually
%    have zero characters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the input file.
%
%    Output, integer ROW_NUM, the number of rows found.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_ROW_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_ROW_COUNT - Error!' );
  end

  blank_num = 0;
  comment_num = 0;
  row_num = 0;

  record_num = 0;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    record_num = record_num + 1;
    record_length = s_len_trim ( line );

    if ( record_length <= 0 )
      blank_num = blank_num + 1;
    elseif ( line(1) == '#' )
      comment_num = comment_num + 1;
    else
      row_num = row_num + 1;
    end

  end

  fclose ( input_unit );

  return
end
function i4mat_print ( m, n, a, title )

%*****************************************************************************80
%
%% I4MAT_PRINT prints an I4MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 September 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of the matrix.
%
%    Input, integer A(M,N), an M by N matrix to be printed.
%
%    Input, string TITLE, a title.
%
  i4mat_print_some ( m, n, a, 1, 1, m, n, title );

  return
end
function i4mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% I4MAT_PRINT_SOME prints out a portion of an I4MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of the matrix.
%
%    Input, integer A(M,N), an M by N matrix to be printed.
%
%    Input, integer ILO, JLO, the first row and column to print.
%
%    Input, integer IHI, JHI, the last row and column to print.
%
%    Input, string TITLE, a title.
%
  incx = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );

  if ( m <= 0 || n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  (None)\n' );
    return
  end

  for j2lo = max ( jlo, 1 ): incx : min ( jhi, n )

    j2hi = j2lo + incx - 1;
    j2hi = min ( j2hi, n );
    j2hi = min ( j2hi, jhi );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col: ' );

    for j = j2lo : j2hi
      fprintf ( 1, '%7d  ', j );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row\n' );
    fprintf ( 1, '\n' );

    i2lo = max ( ilo, 1 );
    i2hi = min ( ihi, m );

    for i = i2lo : i2hi

      fprintf ( 1, '%5d: ', i );
      for j = j2lo : j2hi
        fprintf ( 1, '%7d  ', a(i,j) );
      end

      fprintf ( 1, '\n' );

    end

  end

  return
end
function i4mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% I4MAT_WRITE writes an I4MAT file.
%
%  Discussion:
%
%    An I4MAT is an array of I4's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 June 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, integer TABLE(M,N), the points.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'I4MAT_WRITE - Error!' );
  end
%
%  Write the data.
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %12d', round ( table(i,j) ) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function table = r8mat_data_read ( input_filename, m, n )

%*****************************************************************************80
%
%% R8MAT_DATA_READ reads data from an R8MAT file.
%
%  Discussion:
%
%    An R8MAT is an array of R8's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer M, N, the number of rows and columns of data.
%
%    Output, real TABLE(M,N), the point coordinates.
%

%
%  Build up the format string for reading M real numbers.
%
  string = ' ';

  for i = 0 : m
    string = strcat ( string, ' %f' );
  end

  input_unit = fopen ( input_filename );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_DATA_READ - Error!\n' );
    fprintf ( 1, '  Could not open the file.\n' );
    error ( 'R8MAT_DATA_READ - Error!' );
  end

  table = zeros(m,n);

  i = 0;

  while ( i < n )

    line = fgets ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( line(1) == '#' )

    elseif ( s_len_trim ( line ) == 0 )

    else

      [ x, count ] = sscanf ( line, string );

      if ( count == m )
        i = i + 1;
        table(1:m,i) = x(1:m);
      end

    end

  end

  fclose ( input_unit );

  return
end
function [ m, n ] = r8mat_header_read ( input_filename )

%*****************************************************************************80
%
%% R8MAT_HEADER_READ reads the header from an R8MAT file.
%
%  Discussion:
%
%    An R8MAT is an array of R8's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Output, integer M, the spatial dimension.
%
%    Output, integer N, the number of points.
%
  m = file_column_count ( input_filename );

  if ( m <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data columns in\n' );
    fprintf ( 1, '  the file %s.\n', input_filename );
  end

  n = file_row_count ( input_filename );

  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data rows in\n' );
    fprintf ( 1, '  the file %s\n', input_filename );
  end

  return
end
function len = s_len_trim ( s )

%*****************************************************************************80
%
% S_LEN_TRIM returns the length of a character string to the last nonblank.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be measured.
%
%    Output, integer LENGTH, the length of the string up to the last nonblank.
%
  len = length ( s );

  while ( 0 < len )
    if ( s(len) ~= ' ' )
      return
    end
    len = len - 1;
  end

  return
end
function word_num = s_word_count ( s )

%*****************************************************************************80
%
%% S_WORD_COUNT counts the number of "words" in a string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be examined.
%
%    Output, integer WORD_NUM, the number of "words" in the string.
%    Words are presumed to be separated by one or more blanks.
%
  FALSE = 0;
  TRUE = 1;

  word_num = 0;
  s_length = length ( s );

  if ( s_length <= 0 )
    return;
  end

  blank = TRUE;

  for i = 1 : s_length

    if ( s(i) == ' ' )
      blank = TRUE;
    elseif ( blank == TRUE )
      word_num = word_num + 1;
      blank = FALSE;
    end

  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
lambda = 5;
limits = [0 10; 0 10];

clutter_model = PoissonRateUniformPositionX('ClutterRate',lambda,'Limits',limits);

lambda_samples = clutter_model.random(1000,'cardinality');
spatial_samples = clutter_model.random();

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function varargout = transformPoint3d(varargin)
%TRANSFORMPOINT3D Transform a point with a 3D affine transform
%
%   PT2 = transformPoint3d(PT1, TRANS);
%   PT2 = transformPoint3d(X1, Y1, Z1, TRANS);
%   where PT1 has the form [xp yp zp], and TRANS is a [3x3], [3x4], [4x4]
%   matrix, return the point transformed according to the affine transform
%   specified by TRANS.
%
%   Format of TRANS is a 4-by-4 matrix.
%
%   The function accepts transforms given using the following formats:
%   [a b c]   ,   [a b c j] , or [a b c j]
%   [d e f]       [d e f k]      [d e f k]
%   [g h i]       [g h i l]      [g h i l]
%                                [0 0 0 1]
%
%   PT2 = transformPoint3d(PT1, TRANS)
%   also work when PT1 is a [Nx3xMxPxETC] array of double. In this case,
%   PT2 has the same size as PT1.
%
%   PT2 = transformPoint3d(X1, Y1, Z1, TRANS);
%   also work when X1, Y1 and Z1 are 3 arrays with the same size. In this
%   case, PT2 will be a 1-by-3 cell containing {X Y Z} outputs of size(X1).
%
%   [X2 Y2 Z2] = transformPoint3d(...);
%   returns the result in 3 different arrays the same size as the input.
%   This form can be useful when used with functions like meshgrid or warp.
%
%   MESH2 = transformPoint3d(MESH, TRANS) transforms the field 'vertices'
%   of the struct MESH and returns the same struct with the transformed
%   vertices
%
%   See also:
%   points3d, transforms3d, translation3d, meshgrid
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 10/02/2005.
%

%   23/03/2006 add support for non vector point data
%   26/10/2006 better support for large data handling: iterate on points
%       in the case of a memory lack.
%   20/04/2007 add link to rotationXX functions
%   29/09/2010 fix bug in catch case
%   12/03/2011 slightly reduce memory usage

% parse input arguments
meshStructSwitch = false;
if length(varargin) == 2
    if isstruct(varargin{1}) && isfield(varargin{1}, 'vertices')
        % If first argument is a struct with the field 'vertices' the
        % output will be the same struct, but with the transformed vertices
        meshStructSwitch = true;
        meshStruct = varargin{1};
        varargin{1}=varargin{1}.vertices;
    end
    % Point coordinates are given in a single N-by-3-by-M-by-etc argument.
    % Preallocate x, y, and z to size N-by-1-by-M-by-etc, then fill them in
    dim = size(varargin{1});
    dim(2) = 1;
    [x, y, z] = deal(zeros(dim, class(varargin{1})));
    x(:) = varargin{1}(:,1,:);
    y(:) = varargin{1}(:,2,:);
    z(:) = varargin{1}(:,3,:);
    trans  = varargin{2};
elseif length(varargin) == 4
    % Point coordinates are given in 3 different arrays
    x = varargin{1};
    y = varargin{2};
    z = varargin{3};
    dim = size(x);
    trans = varargin{4};
end

% eventually add null translation
if size(trans, 2) == 3
    trans = [trans zeros(size(trans, 1), 1)];
end

% eventually add normalization
if size(trans, 1) == 3
    trans = [trans;0 0 0 1];
end

% convert coordinates
NP  = numel(x);
try
    % vectorial processing, if there is enough memory
    %res = (trans*[x(:) y(:) z(:) ones(NP, 1)]')';
    %res = [x(:) y(:) z(:) ones(NP, 1)]*trans';
    res = [x(:) y(:) z(:) ones(NP,1,class(x))] * trans';

    % Back-fill x,y,z with new result (saves calling costly reshape())
    x(:) = res(:,1);
    y(:) = res(:,2);
    z(:) = res(:,3);
catch ME
    disp(ME.message)
    % process each point one by one, writing in existing array
    for i = 1:NP
        res = [x(i) y(i) z(i) 1] * trans';
        x(i) = res(1);
        y(i) = res(2);
        z(i) = res(3);
    end
end

% process output arguments
if nargout <= 1
    % results are stored in a unique array
    if length(dim) > 2 && dim(2) > 1
        warning('geom3d:shapeMismatch',...
            'Shape mismatch: Non-vector xyz input should have multiple x,y,z output arguments. Cell {x,y,z} returned instead.')
        varargout{1} = {x,y,z};
    else
        if meshStructSwitch
            meshStruct.vertices = [x y z];
            varargout{1}=meshStruct;
        else
            varargout{1} = [x y z];
        end
    end

elseif nargout == 3
    varargout{1} = x;
    varargout{2} = y;
    varargout{3} = z;
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [neglogli, dL, H] = neglogli_poissGLM(prs,XX,YY,dtbin)
% [neglogli, dL, H] = Loss_GLM_logli_exp(prs,XX);
%
% Compute negative log-likelihood of data undr Poisson GLM model with
% exponential nonlinearity
%
% Inputs:
%   prs [d x 1] - parameter vector
%    XX [T x d] - design matrix
%    YY [T x 1] - response (spike count per time bin)
% dtbin [1 x 1] - time bin size used
%
% Outputs:
%   neglogli   = negative log likelihood of spike train
%   dL [d x 1] = gradient
%   H  [d x d] = Hessian (second deriv matrix)

% Compute GLM filter output and condititional intensity
vv = XX*prs; % filter output
rr = exp(vv)*dtbin; % conditional intensity (per bin)

% ---------  Compute log-likelihood -----------
Trm1 = -vv'*YY; % spike term from Poisson log-likelihood
Trm0 = sum(rr);  % non-spike term
neglogli = Trm1 + Trm0;

% ---------  Compute Gradient -----------------
if (nargout > 1)
    dL1 = -XX'*YY; % spiking term (the spike-triggered average)
    dL0 = XX'*rr; % non-spiking term
    dL = dL1+dL0;
end

% ---------  Compute Hessian -------------------
if nargout > 2
    H = XX'*bsxfun(@times,XX,rr); % non-spiking term
end

author: pillowlab

repo: GLMspiketraintutorial


text:
function [f_constr,label,auxdata] = getFrictionCone(obj, f, fric_coef)
    % returns the symbolic expression of the friction cone
    % constraints of the contact
    %
    % @note The friction cone constraints does not includes the
    % ZMP/CWS constraints.
    %
    % Parameters:
    %  f: the SymVariable of the constraint wrenches
    %  @type SymVariable
    %  fric_coef: the coefficient of friction. @type struct
    %
    % Return values:
    %  f_constr: symbolic function of the friction cone
    %  @type SymFunction
    %  label: the label of the constraints @type cellstr
    %  auxdata: constant data used in the constraint function
    %  @type cell
    %
    % Optional fields of fric_coef:
    %  mu: the (static) coefficient of friction. @type double
    %  gamma: the coefficient of torsional friction @type double

    assert(size(obj.WrenchBase,2) == length(f),...
        ['The dimension of the constraint wrenchs is incorrect.\n',...
        'Expected %d, instead %d'], size(obj.WrenchBase,2), length(f));

    mu = SymVariable('mu');
    gamma = SymVariable('gamma');

    fun_name = ['u_friction_cone_', obj.Name];
    switch obj.Type
        case 'PlanarLineContactWithFriction'
            % x, y, z
            constr = [f(2); % fz >= 0
                f(1) + (mu/sqrt(2))*f(2);
                -f(1) + (mu/sqrt(2))*f(2)];

            % create a symbolic function object
            f_constr = SymFunction(fun_name,...
                constr,{f},{mu});

            % create the label text
            label = {'normal_force';
                'friction_x_pos';
                'friction_x_neg';
                };

            % validate the provided static friction coefficient
            validateattributes(fric_coef.mu,{'double'},...
                {'scalar','real','>=',0},...
                'ContactFrame.getFrictionCone','mu');
            auxdata = {fric_coef.mu};
        case 'PlanarPointContactWithFriction'
            % x, y, z
            constr = [f(2); % fz >= 0
                f(1) + (mu/sqrt(2))*f(2);
                -f(1) + (mu/sqrt(2))*f(2)];

            % create a symbolic function object
            f_constr = SymFunction(fun_name,...
                constr,{f},{mu});

            % create the label text
            label = {'normal_force';
                'friction_x_pos';
                'friction_x_neg';
                };

            % validate the provided static friction coefficient
            validateattributes(fric_coef.mu,{'double'},...
                {'scalar','real','>=',0},...
                'ContactFrame.getFrictionCone','mu');
            auxdata = {fric_coef.mu};
        case 'PointContactWithFriction'
            % x, y, z
            constr = [f(3); % fz >= 0
                f(1) + (mu/sqrt(2))*f(3);
                -f(1) + (mu/sqrt(2))*f(3);
                f(2) + (mu/sqrt(2))*f(3);
                -f(2) + (mu/sqrt(2))*f(3)];
%                 f(4) + gamma * f(3);       % -gamma * fz < wz
%                 -f(4) + gamma * f(3)];     % wz < gamma * fz

            % create a symbolic function object
            f_constr = SymFunction(fun_name,...
                constr,{f},{[mu]});

            % create the label text
            label = {'normal_force';
                'friction_x_pos';
                'friction_x_neg';
                'friction_y_pos';
                'friction_y_neg'};
%                 'tor_firction_neg';
%                 'tor_firction_pos';
%                 };

            % validate the provided static friction coefficient
            validateattributes(fric_coef.mu,{'double'},...
                {'scalar','real','>=',0},...
                'ContactFrame.getFrictionCone','mu');
            % validate the provided torsional friction coefficient
            %             validateattributes(fric_coef.gamma,{'double'},...
            %                 {'scalar','real','>=',0},...
            %                 'ContactFrame.getFrictionCone','gamma');
            auxdata = [fric_coef.mu];
        case 'PointContactWithoutFriction'
            % z
            constr = f;

            % create a symbolic function object
            f_constr = SymFunction(fun_name,...
                constr,{f});

            % create the label text
            label = {'normal_force';
                };
            auxdata = [];
        case 'LineContactWithFriction'
            % x, y, z, roll, yaw
            constr = [f(3); % fz >= 0
                f(1) + (mu/sqrt(2))*f(3);  % -mu/sqrt(2) * fz < fx
                -f(1) + (mu/sqrt(2))*f(3); % fx < mu/sqrt(2) * fz
                f(2) + (mu/sqrt(2))*f(3);  % -mu/sqrt(2) * fz < fu
                -f(2) + (mu/sqrt(2))*f(3); % fy < mu/sqrt(2) * fz
                f(5) + gamma * f(3);       % -gamma * fz < wy
                -f(5) + gamma * f(3)];     % wy < gamma * fz

            % create a symbolic function object
            f_constr = SymFunction(fun_name,...
                constr,{f},{[mu;gamma]});

            % create the label text
            label = {'normal_force';
                'friction_x_pos';
                'friction_x_neg';
                'friction_y_pos';
                'friction_y_neg';
                'tor_firction_neg';
                'tor_firction_pos';
                };

            % validate the provided static friction coefficient
            validateattributes(fric_coef.mu,{'double'},...
                {'scalar','real','>=',0},...
                'ContactFrame.getFrictionCone','mu');

            % validate the provided torsional friction coefficient
            validateattributes(fric_coef.gamma,{'double'},...
                {'scalar','real','>=',0},...
                'ContactFrame.getFrictionCone','gamma');
            auxdata = [fric_coef.mu; fric_coef.gamma];
        case 'LineContactWithoutFriction'
            % z, roll
            constr = f;
            % create a symbolic function object
            f_constr = SymFunction(fun_name,...
                constr,{f});

            % create the label text
            label = {'normal_force';
                };
            auxdata = [];
        case 'PlanarContactWithFriction'
            % x, y, z, roll, pitch, yaw
            constr = [f(3); % fz >= 0
                f(1) + (mu/sqrt(2))*f(3);  % -mu/sqrt(2) * fz < fx
                -f(1) + (mu/sqrt(2))*f(3); % fx < mu/sqrt(2) * fz
                f(2) + (mu/sqrt(2))*f(3);  % -mu/sqrt(2) * fz < fu
                -f(2) + (mu/sqrt(2))*f(3); % fy < mu/sqrt(2) * fz
                f(6) + gamma * f(3);       % -gamma * fz < wz
                -f(6) + gamma * f(3)];     % wz < gamma * fz

            % create a symbolic function object
            f_constr = SymFunction(fun_name,...
                constr,{f},{[mu;gamma]});

            % create the label text
            label = {'normal_force';
                'friction_x_pos';
                'friction_x_neg';
                'friction_y_pos';
                'friction_y_neg';
                'tor_firction_neg';
                'tor_firction_pos';
                };

            % validate the provided static friction coefficient
            validateattributes(fric_coef.mu,{'double'},...
                {'scalar','real','>=',0},...
                'ContactFrame.getFrictionCone','mu');

            % validate the provided torsional friction coefficient
            validateattributes(fric_coef.gamma,{'double'},...
                {'scalar','real','>=',0},...
                'ContactFrame.getFrictionCone','gamma');
            auxdata = [fric_coef.mu; fric_coef.gamma];
        case 'PlanarContactWithoutFriction'
            % z, roll, pitch,
            constr = f;
            % create a symbolic function object
            f_constr = SymFunction(fun_name,...
                constr,{f});

            % create the label text
            label = {'normal_force';
                };
            auxdata = [];
    end

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function oev = getoe(ioev)

% interactive request of classical orbital elements

% NOTE: all angular elements are returned in radians

% input

%  ioev = request array (1 = yes, 0 = no)

% output

%  oev(1) = semimajor
%  oev(2) = orbital eccentricity
%  oev(3) = orbital inclination
%  oev(4) = argument of perigee
%  oev(5) = right ascension of the ascending node (ioev(5) = 1)
%  oev(5) = east longitude of the ascending node (ioev(5) = 2)
%  oev(6) = true anomaly

% Orbital Mechanics with MATLAB

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dtr = pi / 180;

oev = zeros(6,1);

if (ioev(1) == 1)
    % semimajor axis (kilometers)
    % (semimajor axis > req)

    for itry = 1:1:5
        fprintf('\nplease input the semimajor axis (kilometers)');
        fprintf('\n(semimajor axis > 0)\n');

        sma = input('? ');

        if (sma > 0)
            break;
        end
    end
    oev(1) = sma;
end

if (ioev(2) == 1)
    % orbital eccentricity (non-dimensional)
    % (0 <= eccentricity < 1)

    for itry = 1:1:5
        fprintf('\nplease input the orbital eccentricity (non-dimensional)');
        fprintf('\n(0 <= eccentricity < 1)\n');

        ecc = input('? ');

        if (ecc >= 0 && ecc < 1)
            break;
        end
    end
    oev(2) = ecc;
end

if (ioev(3) == 1)
    % orbital inclination (degrees)
    % (0 <= inclination <= 180)

    for itry = 1:1:5
        fprintf('\nplease input the orbital inclination (degrees)');
        fprintf('\n(0 <= inclination <= 180)\n');

        orbinc = input('? ');

        if (orbinc >= 0 && orbinc <= 180)
            break;
        end
    end
    oev(3) = dtr * orbinc;
end

if (ioev(4) == 1)
    % argument of perigee (degrees)
    % (0 <= argument of perigee <= 360)

    if (ecc > 0)

        for itry = 1:1:5

            fprintf('\nplease input the argument of perigee (degrees)');
            fprintf('\n(0 <= argument of perigee <= 360)\n');

            argper = input('? ');

            if (argper >= 0 && argper <= 360)
                break;
            end
        end
    else
        argper = 0;
    end
    oev(4) = dtr * argper;
end

if (ioev(5) >= 1)

    if (ioev(5) == 1)
        % raan (degrees)
        % (0 <= raan <= 360)

        if (oev(3) > 0)
            for itry = 1:1:5
                fprintf('\nplease input the right ascension of the ascending node (degrees)');
                fprintf('\n(0 <= raan <= 360)\n');

                raan = input('? ');

                if (raan >= 0 && raan <= 360)
                    break;
                end
            end
            oev(5) = dtr * raan;
        else
            oev(5) = 0;
        end
    end

    if (ioev(5) == 2)
        % elan (degrees)
        % (0 <= elan <= 360)

        if (oev(3) > 0)
            for itry = 1:1:5
                fprintf('\nplease input the east longitude of the ascending node (degrees)');
                fprintf('\n(0 <= elan <= 360)\n');

                elan = input('? ');

                if (elan >= 0 && elan <= 360)
                    break;
                end
            end
            oev(5) = dtr * elan;
        else
            oev(5) = 0;
        end
    end
end

if (ioev(6) == 1)
    % true anomaly (degrees)
    % (0 <= true anomaly <= 360)

    for itry = 1:1:5
        fprintf('\nplease input the true anomaly (degrees)');
        fprintf('\n(0 <= true anomaly <= 360)\n');

        tanom = input('? ');

        if (tanom >= 0 && tanom <= 360)
            break;
        end
    end
    oev(6) = dtr * tanom;
end

author: Sable


text:
% simple demo to convert a polygon mesh into a voxel representation

% code is 99.99% based on
% [1] http://www.mathworks.com/matlabcentral/fileexchange/24086-polygon2voxel
% [2] http://www.mathworks.com/matlabcentral/fileexchange/21044-3d-voxelizer

%   % Compile the c-coded function
%   mex polygon2voxel_double.c -v

load model;

vertices = vertices - repmat(mean(vertices,1),size(vertices,1),1);

FV.faces = faces;
FV.vertices = vertices;

Volume=polygon2voxel(FV,[100 100 100],'auto');

%% visualization 1
figure
[X,Y,Z]=ind2sub(size(Volume),find(Volume(:)));
plot3(X,Y,Z,'.');
axis equal;
xlabel('x');
ylabel('y');
zlabel('z');

%% visualization 2
% 3d pirnter style visualization to add layer by layer
plot3D(Volume,'timed', 0.1)

figure
%% visualization 3
for i=1:size(Volume,1)
    imagesc(squeeze(Volume(i,:,:)));
    axis equal;
    axis tight;
    axis off
    title(i);
    pause(0.1);
end

for i=1:size(Volume,2)
    imagesc(squeeze(Volume(:,i,:)));
    axis equal;
    axis tight;
    axis off
    title(i);
    pause(0.1);
end

for i=1:size(Volume,3)
    imagesc(squeeze(Volume(:,:,i)));
    axis equal;
    axis tight;
    axis off
    title(i);
    pause(0.1);
end


text:
function words = ExtractWordsFromString(inputString)

inWord = 0;
words = {};

for k=1:length(inputString)
    if inputString(k) == ' '
        isSpace = 1;
    elseif ~isempty( regexp(inputString(k), '\s', 'once') )
        isSpace = 1;
    else
        isSpace = 0;
    end

    if inWord == 1
        if isSpace == 1
            words{end+1} = inputString(word_start:k-1);
            inWord = 0;
        end

    else
        if isSpace == 0
            inWord = 1;
            word_start = k;
        end
    end
end

% if the last symbol is not space and we are still in a word.
if inWord==1
    words{end+1} = inputString(word_start:end);
end

author: singaxiong

repo: SignalGraph


text:
function [ s ] = kpath_length( k, A, mode )
%KPATH_LENGTH Find the distance along a path in k-space.
%   [ s ] = kpath_length( k, A, mode ) calculates the legnth along a series
%   of N points in k-space whose fractional coordinates are given by the
%   3xN array k. This is useful for creating band structure plots. The 3x3
%   matrix A gives the lattice vectors (mode = 'direct') or the reciprocal
%   lattice vectors (mode = 'reciprocal'). If mode is not specified, the
%   default is reciprocal.

    if nargin ==3
        switch mode
            case 'direct'
                B = reciprocal_lattice(A);
            case 'reciprocal'
                B = A;
        end
    end

    k = k*B; % convert to Cartesian coordinates

    s = zeros(1,size(k,1));

    for i = 1:size(k,1)-1
        d = sqrt(sum((k(i,:)-k(i+1,:)).^2));
        s(i+1) = s(i) + d;
    end

end

repo: mcbench-benchmarks


text:
function f = times(f, g, varargin)
%.*   TRIGTECH multiplication.
%   F.*G multiplies TRIGTECH objects F and G or a TRIGTECH by a scalar if either
%   F or G is a scalar.
%
%   If F is an array-valued TRIGTECH, then F.*C is supported if C is a row
%   vector of doubles with the same number of columns as F.
%
% See also MTIMES, RDIVIDE.

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

% TRIGTECH * [] = []:
if ( isempty(f) || isempty(g) )
    f = [];
    return
end

if ( ~isa(f, 'trigtech') )      % Ensure F is a TRIGTECH.

    f = times(g, f, varargin{:});
    return

elseif ( isa(g, 'double') )     % TRIGTECH .* double.

    % Do the multiplication:
    if ( size(g, 2) > 1 )
        f.values = bsxfun(@times, f.values, g);
        f.coeffs = f.vals2coeffs(f.values);
    else
        f.values = f.values*g;
        f.coeffs = f.coeffs*g;
    end
    f.isReal = f.isReal & isreal(g);
    return

elseif ( ~isa(f, 'trigtech') || ~isa(g, 'trigtech') )
    % Don't know how to do the operation.

    error('CHEBFUN:TRIGTECH:times:typeMismatch', ...
        ['Incompatible operation between objects.\n', ...
         'Make sure functions are of the same type.']);

elseif ( size(f.values, 1) == 1 )
    % If we have (constant TRIGTECH).*TRIGTECH, reverse the order and call TIMES
    % again:
    f = times(g, f.values);
    return

elseif ( size(g.values, 1) == 1)
    % If we have TRIGTECH.*(constant TRIGTECH), convert the (constant TRIGTECH)
    % to a scalar and call TIMES again:
    f = times(f, g.values);
    return
end

% Get the size of each TRIGTECH:
[fn, fm] = size(f.values);
[gn, gm] = size(g.values);

% The length of the product is known:
fNew = prolong(f, fn + gn - 1);

% Check dimensions:
if ( fm ~= gm )
    if ( fm == 1 )
        % Allow [Inf x 1] .* [Inf x m].
        fNew.values = repmat(fNew.values, 1, gm);
        fNew.coeffs = repmat(fNew.coeffs, 1, gm);
    elseif ( gm == 1 )
        % Allow [Inf x m] .* [Inf x 1].
        g.values = repmat(g.values, 1, fm);
        g.coeffs = repmat(g.coeffs, 1, fm);
    else
        error('CHEBFUN:TRIGTECH:times:dim2', ...
            'Inner matrix dimensions must agree.');
    end
end

% Check for two cases where the output is known in advance to be positive,
% namely F == conj(G) or F == G and isreal(F).
pos = false;

% Multiply values:
if ( isequal(f, g) )
   values = fNew.values.^2;
   if ( isreal(f) )
       pos = true;
   end
elseif ( isequal(conj(f), g) )
   values = conj(fNew.values).*fNew.values;
   pos = true;
else
   gNew = prolong(g, fn + gn - 1);
   values = fNew.values.*gNew.values;
end

% Assign values and coefficients back to f:
f.values = values;
f.coeffs = f.vals2coeffs(values);

% Update ishappy:
f.ishappy = f.ishappy && g.ishappy;

% Simplify.
f = simplify(f);

if ( pos )
    % Here we know that the product of F and G should be positive. However,
    % SIMPLIFY may have destroyed this property, so we enforce it.
    f.values = abs(f.values);
    f.coeffs = f.vals2coeffs(f.values);
    f.isReal = true(1, size(f.coeffs, 2));
else
    f.isReal = f.isReal & g.isReal;
end

% If f and g are real then make the result real.
f.values(:,f.isReal) = real(f.values(:,f.isReal));

end


text:
function [vhat,iteration]=decode_ldpc_matlab(rx_waveform,No,h,rows,cols,ind,r,c,max_iter)

% bit to check
% check to bit
% update APP LLR
% hard decision
% if valid codework, then exit

vhat(1,1:cols)=0;   %length of decoded word = number of colums of H

s=struct('a',sparse(rows,cols,0),'b',sparse(rows,cols,0),'g',sparse(rows,cols,0));
% a:check to bit messages of previous iteration, b: current check to bit messages,
% g:updated codebit APP LLR

gamma_n=(4/No)*rx_waveform;

s.g(ind)=gamma_n(c);

bitmessage=sparse(rows,cols,0);
bitmessage_2=sparse(rows,cols,0);

for iteration=1:max_iter       %iteration loop

    bitmessage(ind)=tanh( ( -(s.g(ind))+ s.a(ind))  /2);    %bit to check messages

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %loop for computing check-to-bit messages
    for i=1:rows
        cc=c(find(r==i));   %
        rowtemp=bitmessage(i,:);    %ith row of bitmessage

        rowtemp_int=zeros(1,length(cc));
        rowtemp_int=rowtemp(cc);    %interesting elements corresponding to 1's in H
        ze=find(rowtemp_int==0);    %number of zeros in interesting elements

        if length(ze)>=2    %if #of zeros >=2, all interesting elements are zero
            bitmessage_2(i,cc)=0;

        elseif isempty(ze)==1
            rowtemp_int1=zeros(1,length(cc));
            rowtemp_int1(1:length(cc))=prod(rowtemp_int);
            bitmessage_2(i,cc)=rowtemp_int1./rowtemp_int;

        else  %one zero in interesting elements
            zuss=find(rowtemp_int==0);
            rowtemp_int2=rowtemp_int;
            rowtemp_int2(zuss)=1;
            rowtemp_int3=zeros(1,length(cc));
            rowtemp_int3(zuss)=prod(rowtemp_int2);
            bitmessage_2(i,cc)=rowtemp_int3;
        end

    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    s.b(ind)=-2*atanh(bitmessage_2(ind));   %check to bit messages

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    s.b(find(s.b==inf))=1000;    % limiting the checknode to bitnode LLR values to +-50
    s.b(find(s.b==-inf))=-1000;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    s.a(ind)=s.b(ind);                      %current check to bit messages for use in next iteration

    sum_of_b=sum(s.b,1);

    tempsum=sum_of_b+gamma_n;               % updating APP LLR
    s.g(ind)=tempsum(c);

    vhat(find(tempsum>=0))=1;
    vhat(find(tempsum<0))=0;

    if mod(vhat*h',2)==0
        break
    end

end      %iteration loop -end

%   Copyright (c) 2005 by Shaikh Faisal Zaheer, faisal.zaheer@gmail.com
%   $Revision: 1.0 $  $Date: 10/11/2005 $


text:
function trans = createScaling(varargin)
%CREATESCALING Create the 3*3 matrix of a scaling in 2 dimensions.
%
%   TRANS = createScaling(SX, SY);
%   return the matrix corresponding to scaling by SX and SY in the 2
%   main directions.
%   The returned matrix has the form:
%   [SX  0  0]
%   [0  SY  0]
%   [0   0  1]
%
%   TRANS = createScaling(SX);
%   Assume SX and SY are equals.
%
%   TRANS = createScaling(CENTER, ...);
%   Specifies the center of the scaling transform. The argument CENTER
%   should be a 1-by-2 array representing coordinates of center.
%
%   See also
%   transforms2d, transformPoint, createTranslation, createRotation
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2004-04-07
% Copyright 2004-2022 INRA - TPV URPOI - BIA IMASTE

% defined default arguments
sx = 1;
sy = 1;
cx = 0;
cy = 0;

% process input arguments
if nargin == 1
    % the argument is either the scaling factor in both direction,
    % or a 1x2 array containing scaling factor in each direction
    var = varargin{1};
    sx = var(1);
    sy = var(1);
    if length(var)>1
        sy = var(2);
    end
elseif nargin == 2
    % the 2 arguments are the scaling factors in each dimension
    sx = varargin{1};
    sy = varargin{2};
elseif nargin == 3
    % first argument is center, 2nd and 3d are scaling factors
    center = varargin{1};
    cx = center(1);
    cy = center(2);
    sx = varargin{2};
    sy = varargin{3};
end

% concatenate results in a 3-by-3 matrix
trans = [sx 0 cx*(1-sx); 0 sy cy*(1-sy); 0 0 1];

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom