text:
function fem1d_pack_test01 ( )

%*****************************************************************************80
%
%% FEM1D_PACK_TEST01 tests LOCAL_BASIS_1D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    None
%
  node_num = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM1D_PACK_TEST01:\n' );
  fprintf ( 1, '  LOCAL_BASIS_1D evaluates the local basis functions\n' );
  fprintf ( 1, '  for a 1D element.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Test that the basis functions, evaluated at the nodes,\n' );
  fprintf ( 1, '  form the identity matrix.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of nodes = %d\n', node_num );

  node_x(1:node_num) = [ 1.0, 2.0, 4.0, 4.5 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Node coordinates:\n' );
  fprintf ( 1, '\n' );
  for j = 1 : node_num
    fprintf ( 1, '  %8d  %7.3f\n', j, node_x(j) );
  end

  for j = 1 : node_num
    x = node_x(j);
    phi_matrix(1:node_num,j) = local_basis_1d ( node_num, node_x, x );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  A(I,J) = PHI(I) at node (J):\n' );
  fprintf ( 1, '\n' );
  for i = 1 : node_num
    for j = 1 : node_num
      fprintf ( 1, '  %7.3f', phi_matrix(i,j) );
    end
    fprintf ( 1, '\n' );
  end

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The PHI functions should sum to 1 at random X values:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '       X        Sum ( PHI(:)(X) )\n' );
  fprintf ( 1, '\n' );

  for j = 1 : 5
    [ x, seed ] = r8_uniform_ab ( 1.0, 4.5, seed );
    phi = local_basis_1d ( node_num, node_x, x );
    fprintf ( 1, '  %14.6g  %14.6g\n', x, sum ( phi(1:node_num) ) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function w = cplxdual2D(x, J, Faf, af)

% Dual-Tree Complex 2D Discrete Wavelet Transform
%
% USAGE:
%   w = cplxdual2D(x, J, Faf, af)
% INPUT:
%   x - 2-D array
%   J - number of stages
%   Faf{i}: first stage filters for tree i
%   af{i}:  filters for remaining stages on tree i
% OUTPUT:
%   w{j}{i}{d1}{d2} - wavelet coefficients
%       j = 1..J (scale)
%       i = 1 (real part); i = 2 (imag part)
%       d1 = 1,2; d2 = 1,2,3 (orientations)
%   w{J+1}{m}{n} - lowpass coefficients
%       d1 = 1,2; d2 = 1,2
% EXAMPLE:
%   x = rand(256);
%   J = 5;
%   [Faf, Fsf] = FSfarras;
%   [af, sf] = dualfilt1;
%   w = cplxdual2D(x, J, Faf, af);
%   y = icplxdual2D(w, J, Fsf, sf);
%   err = x - y;
%   max(max(abs(err)))
%
% WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
% http://taco.poly.edu/WaveletSoftware/

% normalization
x = x/2;

for m = 1:2
    for n = 1:2
        [lo w{1}{m}{n}] = afb2D(x, Faf{m}, Faf{n});
        for j = 2:J
            [lo w{j}{m}{n}] = afb2D(lo, af{m}, af{n});
        end
        w{J+1}{m}{n} = lo;
    end
end

for j = 1:J
    for m = 1:3
        [w{j}{1}{1}{m} w{j}{2}{2}{m}] = pm(w{j}{1}{1}{m},w{j}{2}{2}{m});
        [w{j}{1}{2}{m} w{j}{2}{1}{m}] = pm(w{j}{1}{2}{m},w{j}{2}{1}{m});
    end
end

author: thomaskuestner


text:
function kbinfo2(cmd)
% KBINFO keyboard information
%   This simple example assists you to understand how to get the keyboard information
%
%   Author: Fu-Sung Wang, E-mail: fusung@webmail.pme.nthu.edu.tw , Update: 23-Dec-2004 11:29:53
%   $Revision

%% Initialization
global Hfig Htxt;        % Htxt: text handle, Hfig: figure handle
global MAXWORDLEN;       % MATWORDLEN: maximum word length

%% Create figure and text, and MAXWORDLEN is also initialized
if nargin == 0           % no input argument needed
Hfig = figure('Name','Keyboard Infomation','Numbertitle','off','Menubar','none',...         % Hfig: figure handle
              'Color',[0.831373 0.815686 0.784314],'Resize','off','DoubleBuffer','on',...
              'Position',[150,150,220,400],'KeyPressFcn',[mfilename,'(''RUN'')']);
Htxt = uicontrol('Tag','typelist','Style', 'ListBox', ...                                      % Htxt: text handle
                 'Position',[0 0 220 400],'FontSize',10, ...
                 'HorizontalAlignment','Left','KeyPressFcn',[mfilename,'(''RUN'')']);
MAXWORDLEN = 200;       % maximum word length
end

%% Retrieve keyboard information
% keyIn = get(Hfig, 'CurrentCharacter');
keyIn = num2str(keyinfo');
rec = getappdata(Hfig,'rec');           % get value from figure object (Hfig)
rec = strvcat(rec,keyIn);                % string concatenation

%% Reset rec
if size(str2num(rec),1) > MAXWORDLEN
    rec = [];
    rec = keyIn;
end

%% Update data
set(Htxt,'String',rec);
len_str = size(str2num(rec),1);
set(Htxt,'Value',len_str);
set(Hfig,'Name',[num2str(len_str) ' / ' num2str(MAXWORDLEN)] );
setappdata(Hfig,'rec',rec);

author: Sable

repo: mcbench-benchmarks


text:
function planC = storeStructuresMatrices(planC, indicesM, structBitsM, scanNum)
%"storeStructuresMatrices"
%   Either 1) modifies planC by adding another cell just before the index cell, or
%          2) uses the preallocated structureArray cell in planC.
%
%this new cell includes the index matrix and the bit matrix for the structures.
%
% 16 Aug 02, V H Clark
% 18 Feb 05, JRA, Added support for multiple scanSets.
% 14 Feb 08, APA, Added UID association
%
%Usage:
%   function planC = storeStructuresMatrices(planC, indicesM, structBitsM)

fieldExists = isfield(planC{end}, 'structureArray'); %this will be true for newer versions, false for old versions
indexS = planC{end};
assScanUID = planC{indexS.scan}(scanNum).scanUID;
structSetUID = createUID('STRUCTURESET');
if fieldExists %most common
    if isempty(planC{planC{end}.structureArray})
        planC{planC{end}.structureArray} = struct('indicesArray', indicesM, 'bitsArray', structBitsM, 'assocScanUID', assScanUID, 'structureSetUID', structSetUID);
    else
        planC{planC{end}.structureArray}(scanNum) = struct('indicesArray', indicesM, 'bitsArray', structBitsM, 'assocScanUID', assScanUID, 'structureSetUID', structSetUID);
    end
else
    structIndex = length(planC);
    planC{end+1} = planC{end}; %moves the index to the end

    planC{structIndex}(scanNum) = struct('indicesArray', indicesM, 'bitsArray', structBitsM, 'assocScanUID', assScanUID, 'structureSetUID', structSetUID);
    %change index
    planC{end}.structureArray = structIndex;
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% MAKEQLD
%
% Script to compile the mex file qld.c and qldorg.c.

% Store the original path
oldpath = pwd;
% Move to the private dir of prtools
np = which('ldc');
[newpath,fname] = fileparts(np);
cd(fullfile(newpath,'private'));
% Compile
mex -c qldorg.c
if strcmp(computer,'PCWIN')
	mex qld.c qldorg.obj
else
	mex qld.c qldorg.o
end
% Return home
cd(oldpath);

author: marianux

repo: ecg-kit

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


text:
function SD = opto2homer(opto)

% OPTO2HOMER constructs a Homer-compatible sensor definition (SD) from a FieldTrip
% opto structure.
%
% See https://www.nitrc.org/plugins/mwiki/index.php/homer2:Homer_Input_Files#NIRS_data_file_format
%
% The Homer SD structure contains the source/detector geometry and has the following fields:
%
% nSrcs    - Number of lasers; scalar variable
% nDets    - Number of detectors; scalar variable
% SrcPos   - Array of probe coordinates of the lasers; dimensions <number of lasers> by 3
% DetPos   - Array of probe coordinates of the detectors; dimensions <number of detectors> by 3
% Lambda   - Wavelengths used for data acquisition; dimensions <number of wavelengths> by 1
% MeasList - List of source/detector/wavelength measurement channels. It’s an array with dimensions, <number of channels> by 4.The meaning of the 4 columns are as follows:
%   Column 1 index of the source from the SD.SrcPos list.
%   Column 2 index of the detector from the SD.DetPos list.
%   Column 3 is unused right now and contains all ones.
%   Column 4 index of the wavelength from SD.Lambda.
%
% The FieldTrip optode structure is defined in FT_DATATYPE_SENS
%
% See also HOMER2OPTO, FT_DATATYPE_SENS

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

[nchan, noptode] = size(opto.tra);
transmitters = find(strcmp(opto.optotype, 'transmitter'));
receivers    = find(strcmp(opto.optotype, 'receiver'));

SD = [];
SD.nSrcs = length(transmitters);
SD.nDets = length(receivers);
SD.SrcPos = opto.optopos(transmitters,:);
SD.DetPos = opto.optopos(receivers,:);
SD.Lambda = opto.wavelength;

SD.MeasList = zeros(nchan, 4);
for i=1:nchan
  optodes = opto.tra(i,:);
  SD.MeasList(i,1) = find(transmitters==find(optodes>0));     % the number of the transmitter
  SD.MeasList(i,2) = find(receivers==find(optodes<0));        % the number of the receiver
  SD.MeasList(i,4) = optodes(optodes>0);                      % the wavelength
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function plotGraph(vertices, edges, figureRef)

if nargin < 3
    figureRef = figure;
end

cla(figureRef);

figure(figureRef);
plot(vertices(:,1),vertices(:,2),'rx');
hold on;

nv = size(vertices,1);
ne = size(edges,1);

v1 = vertices(edges(:,1),:);
v2 = vertices(edges(:,2),:);

vList = ones(ne*3,2);

vList(1:3:ne*3,:) = v1;
vList(2:3:ne*3,:) = v2;
vList(3:3:ne*3,:) = NaN;

plot(vList(:,1),vList(:,2),'b-');

axis equal;
xlim([0 1]);
ylim([0 1]);

author: Doraemonzzz


text:
function obj = SIPEVD_wbmethod(obj)

global H W_mopt Nk Nt Nr  Ns Vn ifVFD  n;
t1 = clock;
j = 0;
w = zeros(Nk,1);
v = zeros(Nk,1);
%init
if (ifVFD)
    W_equal = W_mopt;
else
    W_equal = exp( 1i*unifrnd(0,2*pi,Nr,Ns,Nk) );
end
V_equal = zeros(Nt,Ns,Nk);
H_equal = zeros(Ns,Ns,Nk);
m_mse = zeros(Nk,1);
for i = 1:Nk
    w(i) = trace(W_equal(:,:,i)'*W_equal(:,:,i));
end

H1 = zeros(Nt,Ns,Nk);
H2 = zeros(Nr,Ns,Nk);
trigger = 1;
m_MSE_new = 100;

%limit the iterations number by i<10
while (trigger > 1e-5 && j<10)

    Vn1 = Vn * w;
    for i = 1: Nk
        H1(:,:,i) = H(:,:,i)'*W_equal(:,:,i);
    end
    [V_RF,V_U] = SIPEVD_method(H1,Vn1);

    for i = 1:Nk
        V_equal(:,:,i) = V_RF * V_U(:,:,i);
        v(i) = trace(V_equal(:,:,i)'*V_equal(:,:,i));
        H2(:,:,i) = H(:,:,i)*V_equal(:,:,i);
    end
    Vn2 = Vn * v;
    [W_RF,W_B] = SIPEVD_method(H2,Vn2);

    m_MSE_old = m_MSE_new;

    for k = 1:Nk
        W_equal(:,:,k) = W_RF * W_B(:,:,k);
        w(k) = trace(W_equal(:,:,k)'*W_equal(:,:,k));
        H_equal(:,:,k) = W_equal(:,:,k)'*H2(:,:,k);
        m_mse(k) = trace(H_equal(:,:,k) * H_equal(:,:,k)' - H_equal(:,:,k) - H_equal(:,:,k)')...
            + Vn * v(k) * w(k);
    end
    m_MSE_new = sum(m_mse)/Nk;
    trigger = m_MSE_old - m_MSE_new;
    j = j + 1;
    obj.modmse(j,n) = m_MSE_new + Ns;
end

for i = 1:Nk
    V_B(:,:,i)= V_U(:,:,i) /sqrt(v(i));
end

t2 = clock;
runtime  = etime(t2,t1);
obj.V_B = V_B;
obj.W_B = W_B;
obj.V_RF = V_RF;
obj.W_RF = W_RF;
obj.runtime = obj.runtime + runtime;
obj.iter = obj.iter + j;
obj = get_wbmetric(obj);

author: Zzhaoxingyu

repo: hybrid-beamforming-for-three-scenes


text:
function engine = var_elim_inf_engine(bnet, varargin)
% VAR_ELIM_INF_ENGINE Variable elimination inference engine
% engine = var_elim_inf_engine(bnet)
%
% For details on variable elimination, see
% - R. Dechter, "Bucket Elimination: A Unifying Framework for Probabilistic Inference", UA1 96, pp. 211-219.
% - Z. Li and B. D'Ambrosio, "Efficient inference in Bayes networks as a combinatorial
%     optimization problem", Intl. J. Approximate Reasoning, 11(1):55-81, 1994
% - R. McEliece and S. M. Aji, "The Generalized Distributive Law", IEEE Trans. Inform. Theory, 46(2), 2000

% This is where we will store the results between enter_evidence and marginal_nodes
engine.evidence = [];

engine = class(engine, 'var_elim_inf_engine', inf_engine(bnet));

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [dat, volInfo, cl] = iimg_threshold(image_names, varargin)
% Thresholds images to be at least thresh(1) and at most thresh(2)
%
% :Usage:
% ::
%
%     function [dat, volInfo, cl] = iimg_threshold(image_names, varargin)
%
% :Input types for image names:
%   1) String matrix of filenames
%   2) 4-D array of 3-D data volumes
%   3) voxels x images 2-D index array
%   4) image_vector object
%
% :Outputs:
%
%   **imdat:**
%        index vector of thresholded data
%
%   **volInfo:**
%        structure of info about volume
%
%   **cl:**
%        optional (slower) clusters structure from dat
%
% :Command strings:
%
%   **'imgtype':**
%        followed by 't', 'p', or 'data' (default)
%       specify type of values in image
%
%   **'threshtype':**
%        followed by 't' or 'p'
%
%   **'df':**
%        followed by degrees of freedom, for p- and t-image threshold calc
%
%   **'k':**
%        followed by extent threshold in voxels (slower)
%
%   **'abs':**
%        absolute value must be > threshold.  Use to get both pos. and neg.
%        results in two-tailed test.
%
%   **'intersect:**
%        Create intersection of images; uses abs, so + and - values
%        count as yesses for intersection
%
%   **'contrast':**
%        followed by contrasts across images
%
%   **'volInfo:**
%        followed by volInfo structure.  Necessary for extended output
%
% :Examples:
% ::
%
%    % Threshold an image set (p) to be at least zero
%    image_names =
%    /Users/tor/Documents/Tor_Documents/CurrentExperiments/Lab/Pre-appraisal/Results/rscalped_avg152T1_graymatter_smoothed.img
%    /Users/tor/Documents/Tor_Documents/CurrentExperiments/Lab/Pre-appraisal/Results/t_intercept.img
%
%    [dat, volInfo] = iimg_threshold(image_names, 'thr', 0);
%
%    % Do the same, but take the intersection and write an output image
%    [dat, volInfo] = iimg_threshold(image_names, 'thr', 3, 'outnames', 'intersct', 'masked_t.img');
%
%
%    % Threshold a t-image based on a p-value and a df
%    image_names = '/Users/tor/Documents/Tor_Documents/CurrentExperiments/Lab/Pre-appraisal/Results/t_intercept.img'
%    [dat, volInfo] = iimg_threshold(image_names, 'thr', .005, 'imgtype', 't', 'threshtype', 'p', 'df', 28, 'outnames', 'thresh_t.img');
%
%    cl = mask2clusters('masked_t.img');
%    cluster_orthviews(cl);
%
%    % The same, but threshold based on absolute value (+ and - values)
%    [dat, volInfo] = iimg_threshold(image_names, 'thr', .005, 'abs', 'imgtype', 't', 'threshtype', 'p', 'df', 28, 'outnames', 'thresh_t.img');
%
%    % Threshold a p-value image directly
%    [dat, volInfo] = iimg_threshold('X-M-Y_pvals.img', 'thr', [0 .05], 'outnames', 'X-M-Y_pvals_thresh.img');
%    [dat, volInfo, cl] = iimg_threshold(inname, 'thr', [0 .05], 'outnames', outname);
%
%    % Threshold a p-value image, with cluster sizes
%    [dat, volInfo, cl] = iimg_threshold(inname, 'thr', [0 .05], 'k', 10);
%
%    % Threshold and display a t-image using FDR, getting both positive and negative results
%    [dat, volInfo, cl] = iimg_threshold('contrast_t.img', 'imgtype', 't', 'df', 37, 'thr', .2, 'threshtype', 'fdr', 'k', 3, 'abs');
%    cluster_orthviews(cl);
%    spm_orthviews_hotcool_colormap(cat(2,cl.Z), 1.52);
%
% ..
%    Authorship and Updates:
%    Created by Tor Wager, edited by Matthew Davidson
%    Update April 2007 by TW : correct FDR-thresholding bug
% ..

    % ..
    %    Set up arguments
    % ..

    df = []; % defaults
    % maskname = deblank(image_names(1, :));
    thr = [0 Inf];                  % must be > 0, < Inf
    imgtype = 'data';               % could be data, t, p, or F
    threshtype = 'none';
    mask = [];
    k = 1;
    dointersect = 0;
    con = [];
    outnames = [];
    doabs = 0;                      % absolute value thresholding
    extended_output_flag = 2;    %*(nargout > 2); % yes if we've requested clusters to be created needed for extent also

    % inputs
    for i = 1:length(varargin)
        arg = varargin{i};
        if ischar(arg)
            switch lower(arg)
                case 'imgtype', imgtype = varargin{i+1};
                case 'threshtype', threshtype = varargin{i+1};
                case 'df', df = varargin{i+1};
                case {'thr', 'threshold'}, thr = varargin{i+1};
                case 'k', k = varargin{i+1};
                case 'intersect', dointersect = 1;
                case 'mask', mask = varargin{i+1};
                case 'contrast', con = varargin{i+1};
                case 'outnames', outnames = varargin{i+1};
                case 'abs', doabs = 1;
                case 'volinfo', volInfo = varargin{i+1};
            end
        end
    end

    % Make upper threshold infinity if not entered
    if(length(thr) < 2)
        switch threshtype
            case 'p'
                thr = [0 thr];
            case 'fdr'
                % threshold should be a single p-value, which is converted
                % to [0 fdr_thresh] by convert_threshold, below
                % check and make sure p-value
                if thr < eps || thr > 1
                    error('Invalid threshold (thr) for FDR: must be p-value threshold.');
                end

            case 't'
                thr = [thr Inf];

            case 'none'
                % we should have a data image
                if ~strcmp(imgtype, 'data'), error('You must use imgtype ''data'' (default) if you do not enter a threshold type (''t'' or ''p'')'); end

                thr = [thr Inf];

            otherwise
                error('Unknown threshold type (''threshtype'').');
        end
    end

    %% --------------------------------------
    % * Read image data
    % --------------------------------------
    if isa(image_names, 'statistic_image')
        tmpVolInfo = image_names.volInfo;
        image_names = replace_empty(image_names);
        dat = image_names.dat(:, 1);
    else
        [tmpVolInfo, dat] = iimg_read_img(image_names, extended_output_flag);
    end

    if issparse(dat), dat = full(dat); end
    clear image_names

    % use input volInfo if entered; this is so you can input an indexed image
    % instead of a filename and get extended output (cluster sizes)
    if ~exist('volInfo', 'var'), volInfo = tmpVolInfo; end
    clear tmpVolInfo

    if ~isempty(mask)
        dat = iimg_mask(mask,dat,volInfo);
    end

    if dointersect
        % If any of the input images has a zero (the designated "no value") for a voxel, then clear it out for all images
        wh = any(dat==0,2);
        dat(wh,:) = 0;
    end

    % --------------------------------------
    % * Convert and apply threshold
    % --------------------------------------
    thr = convert_threshold(imgtype, threshtype, thr, df, dat);

    if doabs
        whzero = abs(dat) < thr(1) | abs(dat) > thr(2);  % remove these values
    else
        whzero = dat < thr(1) | dat > thr(2);
    end
    dat(whzero) = 0;

    %% --------------------------------------
    % * Apply size threshold
    % --------------------------------------
    [dat, nvox] = iimg_cluster_extent(dat, volInfo, k);

    %% --------------------------------------
    % special operations
    % --------------------------------------
    if ~isempty(con)
        % Apply contrast(s) across images
        dat = apply_contrast(dat, con);
    end

    %% --------------------------------------
    % * Write out image names, if asked
    % --------------------------------------
    if ~isempty(outnames)
        iimg_write_images(dat, volInfo, outnames);
    end

    %% --------------------------------------
    % * clusters output
    % --------------------------------------
    if extended_output_flag

        % could enter u and k here, but don't need to because it's already
        % thresholded; also, neg. elements will be removed in some cases if
        % threshold is re-applied here.

        %cl = iimg_indx2clusters(dat, volInfo, thr, k);
        cl = iimg_indx2clusters(dat, volInfo);
    end
end

% Sub-functions

function thr = convert_threshold(imgtype, threshtype, thr, df, dat)
    % handle FDR
    switch imgtype
        case 'data'
            if(strcmp(threshtype, 'fdr'))
                error('FDR thresholding has no meaning for data images. Try using a t- or a p-image');
            end
        case 't'
            switch threshtype
                case 't'
                case 'p'
                    % convert p-value threshold to t-value threshold
                    if ~exist('df', 'var') || isempty(df)
                        error('You must enter df to use p-value based thresholding with a t img.');
                    end
                    thr(~isinf(thr)) = tinv(1-thr(~isinf(thr)), df);
                case 'fdr'
                    % convert t-values in image to p-values for FDR
                    % thresholding
                    if ~exist('df', 'var') || isempty(df)
                        error('You must enter df to use p-value based thresholding with a t img.');
                    end
                    pdat = 1 - tcdf(dat(dat ~= 0 & ~isnan(dat)), df); % zero is a special invalid value

                    fdr_threshold = FDR(pdat, thr); % get p-value threshold
                    if isempty(fdr_threshold), fdr_threshold = 0; end

                    fprintf('Computed FDR threshold. For p < %.3f corrected, p-threshold is %.7f, ', thr, fdr_threshold);
                    thr = [tinv(1-fdr_threshold, df) Inf];

                    fprintf('t-threshold is %.2f\n', thr(1));

                otherwise
                    error('Threshold type must be t or p values for t-images')
            end

        case 'p'
            switch threshtype
                case 't'
                    if ~exist('df', 'var') || isempty(df)
                        error('You must enter df to use t-value based thresholding with a p img.');
                    end
                    thr = 1 - tcdf(thr, df);
                case 'p' % do nothing
                case 'fdr'

                    dat = dat(dat ~= 0 & ~isnan(dat)); % zero is a special invalid value

                    fdr_threshold = FDR(dat, thr);
                    if isempty(fdr_threshold), fdr_threshold = 0; end

                    %thr(isinf(thr)) = fdr_threshold;
                    fprintf('Computed FDR threshold. For p < %.3f corrected, p-threshold is %.7f, ', thr, fdr_threshold);

                    thr = [0 fdr_threshold];  % from zero to new thresh
                otherwise
                    error('Threshold type must be t or p values for p-images')
            end
    end
end

function condat = apply_contrast(dat, con)
    if size(con, 2) ~= size(dat, 2)
        error('Contrast matrix must have as many columns as images.');
    end

    condat(:, i) = dat * con';
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function rmStimInfo(view, fid);
% rmStimInfo - get information of current loaded rm stimulus
% Useful if you want to predict the reponses to a scan using a previously
% solved rm and a new stimulus. Loading the model generally also loads the
% stimulus associated with the model, which may be different from the
% stimulus associated with the current scan.
%
% rmStimInfo(view, id);
%
% 11/2006 SOD: wrote it.
%  9/2008  JW: added 'fid' as input argument to permit printing to file

% to do: need more specific information?

if notDefined('view'), view = getCurView; end;
if notDefined('fid'), fid = 1; end

% load params
try
    params = viewGet(view,'rmparams');
catch,
    params = [];
end;

% if no parameters are loaded exit
if isempty(params),
    disp(sprintf('[%s]:No stimulus parameters loaded.',mfilename));
    return;
end;

nscans = numel(params.stim);

% give info:
fprintf(fid,'\n---------------------------------------\n');
fprintf(fid,'[%s]:There are stimulus descriptions for %d scan(s).\n',mfilename, nscans);

for n=1:nscans,
    fprintf(fid,'\tScan %d: %s\n',    n,params.stim(n).stimType);

    fprintf(fid,'\t\tStim size: %d deg (radius)\n',  params.stim(n).stimSize);

    fprintf(fid,'\t\tNumber of frames: %d\n', params.stim(n).nFrames);

    if strcmpi(params.stim(n).stimType, 'stimfromscan')
        fprintf(fid,'\t\tFilter Type: %s\n',     params.stim(n).imFilter);
        fprintf(fid,'\t\tImage file:\n\t\t\t%s\n',     params.stim(n).imFile);
        fprintf(fid,'\t\tParameter file:\n\t\t\t%s\n', params.stim(n).paramsFile);
    end

end
fprintf(fid,'---------------------------------------\n\n');

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% sub module for 1. order kinetics in simpltrans
c1 = exp(-lambda*dtdiff)*c1;


text:
function output = apply_net_filter(input_v, input_h)
    global config mem;
    input_size = size(input_v);
    %output = config.NEW_MEM(zeros(size(input_v, 1), size(input_v, 2), config.chs, config.batch_size));
    output = zeros(size(input_v, 1), size(input_v, 2), config.chs, config.batch_size);
    p_size = config.MEM.p_size;
    size_differ = config.size_differ;
    %count = config.NEW_MEM(zeros(size(input_v, 1), size(input_v, 2)));
    count = zeros(size(input_v, 1), size(input_v, 2));
    for v = 1 : length(config.MEM.start_rows)
        for h = 1 : length(config.MEM.start_cols)
            v_start = config.MEM.start_rows(v);
            v_end = v_start + p_size - 1;
            h_start = config.MEM.start_cols(h);
            h_end = h_start + p_size - 1;
            if(v_end > input_size(1))
                v_end = input_size(1);
                v_start = v_end - p_size + 1;
            end
            if(h_end > input_size(2))
                h_end = input_size(2);
                h_start = h_end - p_size + 1;
            end
            input_piece = cat(4, input_v(v_start:v_end, h_start:h_end,:), permute(input_h(v_start:v_end, h_start:h_end,:), [2 1 3]));

            op_test_pipe(input_piece, mem.fake_output_for_test);
            %output_piece = mem.output;
            output_piece = gather(mem.output);
            if(size_differ(1) ~= size_differ(2))
                output_piece = output_piece(size_differ(2)/2+1:size(output_piece,1)-size_differ(2)/2, size_differ(1)/2+1:size(output_piece,2)-size_differ(1)/2,:,:);
            end
            output_piece(:,:,:,2) = permute(output_piece(:,:,:,2), [2 1 3]);

            output(v_start+(max(size_differ)/2):v_end-(max(size_differ)/2), h_start+(max(size_differ)/2):h_end-(max(size_differ)/2),:,:) = ...
                    output(v_start+(max(size_differ)/2):v_end-(max(size_differ)/2), h_start+(max(size_differ)/2):h_end-(max(size_differ)/2),:,:) + output_piece;
            count(v_start+(max(size_differ)/2):v_end-(max(size_differ)/2), h_start+(max(size_differ)/2):h_end-(max(size_differ)/2)) = ...
                    count(v_start+(max(size_differ)/2):v_end-(max(size_differ)/2), h_start+(max(size_differ)/2):h_end-(max(size_differ)/2)) + 1;
        end
    end
    count = max(count, 1);
    output = bsxfun(@rdivide, output, count);
    %output = gather(output);

    %output = output_piece;
    %output = padarray(output, [1 1]);
end

author: jimmy-ren


text:
function imResult = blendMode_PinLight(A, B, offsetW, offsetH)
%% Pin Light blending mode: replaces colors depending on the brightness of
%   the blend color. If the blend color is more than 50% brightness and the
%   base color is darker than the blend color, then the base color is
%   replaced with the blend color. If the blend color is less than 50%
%   brightness and the base color is lighter than the blend color, then the
%   base color is replaced with the blend color.
%
% Input:
%       A       -       Base Image
%       B       -       Top Image
%   offsetW     -   move picture B horizontally in respect to the top-left
%                   corner of picture A. Default value = 1.
%   offsetH     -   move picture B vertically in respect to the top-left
%                   corner of picture A. Default value = 1.
%
% Output:
%       imResult    -   Result of the blending, having the same size of the
%                       Base Image A.
%

%% Check Input
a = size(A);
b = size(B);
blendMode_checkInput(nargin, a, b, func2str(@blendMode_PinLight));

if nargin < 3
    offsetW = 1;
    offsetH = 1;
end

if nargin < 4
    offsetH = 1;
end

%% Implementation
imResult = A;

if (((offsetW ~= 1) || (offsetH ~= 1)) || (sum(a == b) ~= length(a)))
    [A, B] = blendMode_ResizeImages(A, B, a, b, offsetW, offsetH);
end

ind = B < 0.5;
indCompl = abs(ind - 1);

tmp1 = blendMode_Darken(A, 2 * B);
tmp2 = blendMode_Lighten(A, 2 * (B - 0.5));
C = ind .* tmp1 + indCompl .* tmp2;

if (((offsetW ~= 1) || (offsetH ~= 1)) || (sum(a == b) ~= length(a)))
    imResult = blendMode_CreateResult(imResult, C, offsetW, offsetH);
else
    imResult = C;
end

end


text:
function printHyperlink( URL, OptionalText )
% Prints a hyperlink to the command window
%
% Syntax:	PRINTHYPERLINK( URL, OPTIONALTEXT)
%
% Inputs:
% 	URL - Address for the hyperlink
% 	OptionalText - The text to be printed which is hyperlinked
%
% See also: fprintf, cprintf

% Author: Jacob Donley
% University of Wollongong
% Email: jrd089@uowmail.edu.au
% Copyright: Jacob Donley 2017
% Date: 16 February 2017
% Revision: 0.2 (30 March 2017)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin < 2
    OptionalText = URL;
end
fprintf(['' OptionalText '\n']);
if exist('cprintf','file')
    cprintf('hyper', '');
end

end

author: jtkim-kaist

repo: Speech-enhancement


text:
function [Q, vertexBool, fileNameOut] = lrsReadRay(modelName,param)
% Read in a vertex representation (*.ext) of a polytope derived from lrs
% See http://cgm.cs.mcgill.ca/~avis/C/lrslib/USERGUIDE.html#file
%
% USAGE:
%
%    [Q, vertexBool, fileNameOut] = lrsReadRay(modelName,param)
%
% INPUT:
% modelName     string giving the prefix of the *.ext file that will contain the vertex representation
%               It is assumed the file is pwd/*.ine, otherwise provide the full path.
%
% OUTPUT:
% Q             m x n integer matrix where each row is a variable and each column is a vertex or ray
% vertexBool         n x 1 Boolean vector indicating which columns of Q are vertices
%                    By default, all columns of Q are assumed to be rays.
%
% Ronan Fleming 2021

if ~exist('param','var')
    param = struct();
end
if ~isfield(param,'positivity')
    param.positivity = 0;
end
if ~isfield(param,'inequality')
    param.positivity = 0;
end
if ~isfield(param,'sh')
    param.sh = 0;
end
if ~exist('modelName','var')
    modelName = 'test';
end

if param.inequality == 0
%     if param.positivity == 0
%         modelName = [modelName '_pos_eq'];
%     else
%         modelName = [modelName '_neg_eq'];
%     end
else
    if param.positivity == 1
        modelName = [modelName '_pos_ineq'];
    else
        modelName = [modelName '_neg_ineq'];
    end
end
fileNameOut = [modelName '.ext'];

fid = fopen(fileNameOut);
if fid<0
    disp(fileNameOut)
    error('Could not open lrs output file.');
end

% pause(eps)
while 1
    tline = fgetl(fid);
    if strcmp(tline, 'begin')
        break;
    elseif ~ischar(tline)
        error('Could not read lrs output file.');
    end
end

% find the number of columns
C = textscan(fid, '%s %f %s', 1);
nCols = C{2};

% move on pointer one line
fgetl(fid);

% count the number of rows in the file
nRows = 0;
while 1
    if ~strcmp(fgetl(fid), 'end')
        nRows = nRows + 1;
    else
        break;
    end
end

fclose(fid);

% pwd
fid = fopen(fileNameOut);
if fid<0
    disp(fileNameOut)
    error('Could not open lrs output file.');
end

while 1
    if strcmp(fgetl(fid), 'begin')
        break;
    end
end

% find the number of columns
C = textscan(fid, '%s %f %s', 1);
nCols = C{2};

% move on pointer one line
fgetl(fid);

% read rows into a matrix
P = sparse(nRows, nCols);

for r = 1:nRows
    line = fgetl(fid);
    if ~contains(line,'/')
        scannedLine = sscanf(line, '%d')';  % added transpose here for reading in LP solutions
        P(r, :) = scannedLine;
    else
        line = strrep(line, '/', '.');
        scannedLine = sscanf(line, '%f')';
        for c = 1:nCols
            M = mod(scannedLine(c), 1);
            if M ~= 0
                F = fix(scannedLine(c));
                scannedLine(c) = F / M;
            else
                scannedLine(c) = int16(scannedLine(c));
            end
        end
        % pause(eps);
    end
end
fclose(fid);

% Each vertex is given in the form
% 1   v0   v 1 ...   vn-1
V = P(P(:, 1) ~= 0, 2:end)';
%remove zero vertices
V = V(:,sum(V,1)~=0);

% Each ray is given in the form
% 0   r0   r 1 ...   rn-1
R = P(P(:, 1) == 0, 2:end)';  % not the transpose

% order the rays by the number of nnz
[mlt, nlt] = size(R);
nNonZero = zeros(nlt, 1);
for n = 1:nlt
    nNonZero(n) = nnz(R(:, n));
end
%remove zero rays
R = R(:,nNonZero~=0);
[B, IX] = sort(nNonZero);
R = R(:, IX);

%first vertices then rays
Q = [V, R];
vertexBool = false(size(Q,2),1);
vertexBool(1:size(V,2))=1;

author: opencobra

repo: cobratoolbox


text:
function jed = transition_to_jed_english ( )

%*****************************************************************************80
%
%% TRANSITION_TO_JED_ENGLISH returns the English calendar transition as a JED.
%
%  Discussion:
%
%    In the English calendar, the last moment of the Julian calendar was
%      11:59 pm, 2 September 1752 Julian/English,
%      11:59 pm, 13 September 1752 Gregorian/CE.
%    The first minute of the Gregorian calendar ended at
%      12:01 am, 3 September 1752 Julian,
%      12:01 am, 15 September 1752 Gregorian/CE/English.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 December 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real JED, the Julian Ephemeris Date of the date.
%
  jed = 2361221.5;

  return
end


text:
function C = complex( A, B )
%COMPLEX Construct complex SEPARABLEAPPROX from real and imaginary parts.
%   C = COMPLEX(A, B) returns the complex SEPARABLEAPPROX A + Bi, where A and B are
%   real valued SEPARABLEAPPROX objects with the same domain.
%
%   C = COMPLEX(A) for real SEPARABLEAPPROX A returns the complex result C with real
%   part A and all zero imaginary part. isreal(C) returns false.
%
% See also IMAG, CONJ, ABS, REAL.

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

if ( nargin == 2 )
    if ( ~isa(B, 'separableApprox') )
        error('CHEBFUN:SEPARABLEAPPROX:complex:inputs', ...
            'Second input must be a CHEBFUN2.');
    elseif ( ~isreal( A ) || ~isreal( B ) )
        error('CHEBFUN:SEPARABLEAPPROX:complex:notReal1', ...
            'Inputs must be real valued.');
    end
    C = A + 1i*B;
else
    if ( ~isreal( A ) )
        error('CHEBFUN:SEPARABLEAPPROX:complex:notReal2', ...
            'Input must be real valued.');
    end
    % Make complex.
    C = A + 0*1i;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% this script tests NGLDM features between CERR and pyradiomics on wavelet filtered image.
%
% RKP, 03/22/2018

%% Load image
gldmParamFileName = fullfile(fileparts(fileparts(getCERRPath)),...
    'Unit_Testing','tests_for_cerr','test_ngldm_radiomics_extraction_settings.json');
cerrFileName = fullfile(fileparts(fileparts(getCERRPath)),...
    'Unit_Testing','data_for_cerr_tests','CERR_plans','head_neck_ex1_20may03.mat.bz2');

planC = loadPlanC(cerrFileName,tempdir);
indexS = planC{end};

paramS = getRadiomicsParamTemplate(gldmParamFileName);
strNum = getMatchingIndex(paramS.structuresC{1},{planC{indexS.structures}.structureName});
scanNum = getStructureAssociatedScan(strNum,planC);

scanType = 'wavelet';
dirString = 'HHH';

%% NGLDM features CERR

ngldmS = calcGlobalRadiomicsFeatures...
            (scanNum, strNum, paramS, planC);
ngldmS = ngldmS.Wavelets_Coif1__HHH.ngldmFeatS;

cerrNgldmV = [ngldmS.lde, ngldmS.hde, ngldmS.lgce, ngldmS.hgce, ...
    ngldmS.ldlge, ngldmS.ldhge, ngldmS.hdlge, ngldmS.hdhge, ...
    ngldmS.gln, ngldmS.glnNorm, ngldmS.dcn, ngldmS.dcnNorm,...
    ngldmS.dcp, ngldmS.glv, ngldmS.dcv, ngldmS.entropy, ngldmS.energy];

% %% Calculate features using pyradiomics
%
% testM = single(planC{indexS.scan}(scanNum).scanArray) - ...
%     single(planC{indexS.scan}(scanNum).scanInfo(1).CTOffset);
% mask3M = zeros(size(testM),'logical');
% [rasterSegments, planC, isError] = getRasterSegments(strNum,planC);
% [maskBoundBox3M, uniqueSlices] = rasterToMask(rasterSegments, scanNum, planC);
% mask3M(:,:,uniqueSlices) = maskBoundBox3M;
%
% dx = planC{indexS.scan}(scanNum).scanInfo(1).grid1Units;
% dy = planC{indexS.scan}(scanNum).scanInfo(1).grid1Units;
% dz = mode(diff([planC{indexS.scan}(scanNum).scanInfo(:).zValue]));
% pixelSize = [dx dy dz]*10;
%
% teststruct = PyradWrapper(testM, mask3M, pixelSize, scanType, dirString);
% %teststruct = PyradWrapper(testM, mask3M, scanType, dirString);
%
% pyradNgldmNamC = {'SmallDependenceEmphasis', 'LargeDependenceEmphasis',...
%     'LowGrayLevelCountEmphasis', 'HighGrayLevelCountEmphasis',  'SmallDependenceLowGrayLevelEmphasis', ...
%     'SmallDependenceHighGrayLevelEmphasis', 'LargeDependenceLowGrayLevelEmphasis', ...
%     'LargeDependenceHighGrayLevelEmphasis', 'GrayLevelNonUniformity', 'GrayLevelNonUniformityNorm', ...
%     'DependenceNonUniformity', 'DependenceNonUniformityNormalized', ...
%     'DependencePercentage', 'GrayLevelVariance', 'DependenceVariance', ...
%     'DependenceEntropy', 'DependenceEnergy'};
%
%
% pyradNgldmNamC = strcat(['wavelet','_', dirString, '_gldm_'],pyradNgldmNamC);
%
% pyRadNgldmV = [];
% for i = 1:length(pyradNgldmNamC)
%     if isfield(teststruct,pyradNgldmNamC{i})
%         pyRadNgldmV(i) = teststruct.(pyradNgldmNamC{i});
%     else
%         pyRadNgldmV(i) = NaN;
%     end
% end
%
% %% Compare
% ngldmDiffV = (cerrNgldmV - pyRadNgldmV) ./ cerrNgldmV * 100

%% Compare with previously calculated pyradiomics NGLDM features
saved_pyRadNgldmV = [0.204701151831240,55.0338803599788,NaN,NaN,0.000365710992103713,196.983572569714,0.0594517626789526,51247.9161461091,1583.25410269984,NaN,817.939756484913,0.0866002918459410,NaN,11.2964611875247,18.6018880476280,6.32953752322967,NaN];
ngldmDiffV = (cerrNgldmV - saved_pyRadNgldmV) ./ cerrNgldmV * 100


text:
% Normalized Leaky Kernel Affine Projection Algorithm
%
% W. Liu and J.C. Principe, "Kernel Affine Projection Algorithms", EURASIP
% Journal on Advances in Signal Processing, Volume 2008, Article ID 784292,
% 12 pages. http://dx.doi.org/10.1155/2008/784292
%
% Remark: This implementation includes a maximum dictionary size M. With
% M=Inf this algorithm is equivalent to KAPA-4 from the publication. With
% M=Inf and lambda=0 it is equivalent to KAPA-2.
%
% This file is part of the Kernel Adaptive Filtering Toolbox for Matlab.
% https://github.com/steven2358/kafbox/

classdef nlkapa < kernel_adaptive_filter

    properties (GetAccess = 'public', SetAccess = 'private') % parameters
        eta = .05; % learning rate
        eps = 1E-4; % Newton regularization
        lambda = 1E-2; % Tikhonov regularization
        M = 1000; % maximum dictionary size
        P = 20; % number of regressors
        kerneltype = 'gauss'; % kernel type
        kernelpar = 1; % kernel parameter
    end

    properties (GetAccess = 'protected', SetAccess = 'private') % variables
        xmem = []; % input memory
        ymem = []; % output memory
        dict = []; % dictionary
        alpha = []; % expansion coefficients
    end

    methods
        function kaf = nlkapa(parameters) % constructor
            if (nargin > 0) % copy valid parameters
                for fn = fieldnames(parameters)'
                    if ismember(fn,fieldnames(kaf))
                        kaf.(fn{1}) = parameters.(fn{1});
                    end
                end
            end
        end

        function y_est = evaluate(kaf,x) % evaluate the algorithm
            if size(kaf.dict,1)>0
                k = kernel(kaf.dict,x,kaf.kerneltype,kaf.kernelpar);
                y_est = k'*kaf.alpha;
            else
                y_est = zeros(size(x,1),1);
            end
        end

        function train(kaf,x,y) % train the algorithm
            if size(kaf.dict,2)==0 % initialize
                kaf.dict = x;
                kaf.alpha = kaf.eta*y;
                kaf.xmem = x;
                kaf.ymem = y;
            else
                if size(kaf.dict,1) < kaf.M
                    if size(kaf.xmem,1)<kaf.P
                        % grow memory
                        kaf.xmem = [kaf.xmem; x];
                        kaf.ymem = [kaf.ymem; y];
                    else
                        % slide memory
                        kaf.xmem = [kaf.xmem(2:kaf.P,:); x];
                        kaf.ymem = [kaf.ymem(2:kaf.P); y];
                    end

                    ymem_est = kaf.evaluate(kaf.xmem);
                    e = kaf.ymem - ymem_est;
                    G = kernel(kaf.xmem,kaf.xmem,...
                        kaf.kerneltype,kaf.kernelpar);

                    kaf.dict = [kaf.dict; x]; % grow dictionary
                    kaf.alpha = [(1-kaf.lambda*kaf.eta)*kaf.alpha; 0];
                    % leak and grow coefficients

                    m = size(kaf.alpha,1);
                    p = size(kaf.xmem,1);

                    % update p last coefficients
                    kaf.alpha(m-p+1:m) = kaf.alpha(m-p+1:m) + ...
                        kaf.eta*inv(G+kaf.eps*eye(p))*e; %#ok<MINV>
                    % prefer inv to \ to avoid instability
                end
            end
        end
    end
end

repo: kafbox


text:
function [MCI] = spm_mci_mfx_dynamic (MCI)
% Mixed Effects Inference for Dynamical Systems
% FORMAT [MCI] = spm_mci_mfx_dynamic (MCI)
%
% MCI               Data structure containing fields:
%
% .M{n}             Model for nth of N replications (e.g. subjects)
% .U{n}             Inputs for nth replication
% .Y{n}             Data for nth replication
% .fixed            Gaussian prior (.pE and .pC) over FFX
% .S                Second level model describing population mean, m, and
%                   precision, Lambda. The parameters in S.prior
%                   define the sufficient statistics of p(Lambda) (.a and .B)
%                   and p(m|Lambda) (.beta and.m)
% .total_its        Total number of samples per subject
% .rinit            Proportion of samples to collect prior to use of
%                   Empirical (group) prior
% .verbose          Show progress of optimisation
%
% KNOWN, FIXED or RANDOM EFFECTS:
% The initial states, flow and output
% parameters can be fixed or random effects or take on known values:
%
% .assign.init_par  'fixed', 'random' or 'known'
% .assign.flow_par  'fixed', 'random' or 'known'
% .assign.out_par   'fixed', 'random' or 'known'
%
% .pinit0           Initial values of initial state parameters
%                   [Ninit x 1] for fixed, [Ninit x N] for random
% .pflow0           Initial values of flow parameters
%                   [Nflow x 1] for fixed, [Nflow x N] for random
% .pout0            Initial values of output parameters
%                   [Nout x 1] for fixed, [Nout x N] for random
%
% .update_obs_noise Update observation noise, Gamma ? [yes/no] (1/0), default=1
% .verbose          Show progress of optimisation
%
% The output fields are:
%
% POSTERIOR SAMPLES:
% .sv               [Nv x Nsamples] group fixed effects samples, v
% .sm               [Nw x Nsamples] group random effect means, m
% .sw               [Nw x N x Nsamples] subject random effects, w
% .Ce               [Ny x Ny x Nsamples] Obs noise covariance samples
% .postind          Indices for posterior (ie. excluding burn-in)
%
% POSTERIOR MEANS:
% .sv_mean          [Nv x 1] posterior mean over v
% .sm_mean          [Nw x 1] posterior mean over m
% .sw_mean          [Nw x N] posterior mean over w
%
% For Dynamical Systems models:
% .pinit            Estimated initial states
% .pflow            Estimated flow
% .pout             Estimated output params
%
% SUFFICIENT STATISTICS:
% .noise            Parameters of p(Gamma|Y,w,v): .c0,.D0,.cN,.DN
% .S.post           Parameters of p(Lambda|w) (.a and.B)
%                   and p(m|Lambda,w) (.beta and .m)
%
% W.Penny, M Klein-Flugge and B Sengupta (2015) Mixed Effects Langevin
% Monte Carlo, Submitted, 2015.
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_mci_mfx_dynamic.m 6697 2016-01-27 14:57:28Z spm $

try, verbose=MCI.verbose; catch, verbose=0; end

% Update observation noise ?
try, update_obs_noise=MCI.update_obs_noise; catch, update_obs_noise=1; end

% Numbers of iterations, with and without group model
try, total_its=MCI.total_its; catch, total_its=1024; end
try, rinit=MCI.rinit; catch, rinit=0.25; end
rfx_its=2;
mfx_its=(1-rinit)*total_its/rfx_its;
rfx1_its=rinit*total_its;

try, ffx_its=MCI.ffx_its; catch, ffx_its=4; end
try, ffx1_its=MCI.ffx1_its; catch, ffx1_its=64; end

M=MCI.M;U=MCI.U;Y=MCI.Y;
Np=length(spm_vec(M{1}.pE));
if isfield(M{1},'x0')
    % For dynamical systems
    d=size(M{1}.x0,1);
end

% Observation noise
Ny=size(Y{1}.y,2);
noise.c0=Ny;
noise.D0=eye(Ny);

% Prior over fixed effects
try
    fixed=MCI.fixed;
    fixed.vpE=spm_vec(fixed.pE); % m_v (vector)
catch
    fixed=[];
end

% Prior over random effects (second level model)
try
    S=MCI.S;
    Nrand=S.prior.a;
catch
    if strcmp(MCI.assign.init_par,'random')
        Nrand=d;
    else
        Nrand=0;
    end
    if strcmp(MCI.assign.flow_par,'random')
        Nrand=Nrand+M{1}.Npflow;
    end
    if strcmp(MCI.assign.out_par,'random')
        Nrand=Nrand+M{1}.Npout;
    end
    S.prior.P=Nrand;
    S.prior.a=Nrand/2;
    S.prior.B=eye(Nrand);
    S.prior.beta=1;
    S.prior.m=zeros(Nrand,1);
    S.N=length(M);
end
% Prior mean and cov of random effects from second level model
m=S.prior.m;
C=S.prior.B/S.prior.a;
rfx_prior.pE=m;
rfx_prior.pC=C;
MCI.S=S;

% Initialise fixed and random effects
[w_init,v_init,assign,update_ffx,update_rfx] = spm_mci_vw_init (MCI);

% Sampling params
fixed_mcmc.assign=assign;
subj_mcmc.assign=assign;
subj_mcmc.verbose=0;

% Set up sample matrix
Nfixed = length(v_init);
if Nfixed > 0
    %v = zeros(mfx_its,Nfixed);
    for j=1:rfx1_its,
        v(j,:) = v_init';
    end
end
w = w_init;

% Initialisation of each model
subj_mcmc.verbose=verbose;
subj_mcmc.maxits=rfx1_its;
for n=1:S.N,
    subj_mcmc.init=w(:,n);
    if rfx1_its > 0
        Psamp = spm_mci_random (subj_mcmc,rfx_prior,v_init,M{n},U{n},Y{n});
        sw(n,:,:) = Psamp;
        w(:,n) = Psamp(:,end);
    else
        sw(n,:,1) = w(:,n);
    end
end
sw=permute(sw,[2 1 3]);
sm=squeeze(mean(sw,2));
for j=1:rfx1_its,
    Ce(:,:,j)=M{1}.Ce;
end
subj_mcmc.maxits=rfx_its;
subj_mcmc.verbose=0;

% Main Loop
for it=1:mfx_its

    if verbose
        disp(sprintf('MCI iteration %d',it));
    end

    if it>1 && Nrand > 0
        % Update second level params
        S = spm_nwpost (S,w);
        [m,Lambda,C] = spm_nwrnd (S.post,1);
    end

    if update_ffx
        if verbose, disp('Updating fixed effects'); end
        % 2. Update estimate of group fixed effects

        % Start sampling where left off
        fixed_mcmc.init=v(it,:)';
        fixed_mcmc.update_obs_noise=update_obs_noise;

        if it==1
            fixed_mcmc.verbose=1;
            fixed_mcmc.maxits=ffx1_its;
        else
            fixed_mcmc.maxits=ffx_its;
            fixed_mcmc.verbose=0;
        end
        [Psamp,noise,M] = spm_mci_fixed (fixed_mcmc,w,fixed,noise,M,U,Y);
        P = Psamp(end,:);
        v(rfx1_its+it,:) = P;
        vit=P';
    else
        vit=[];
    end

    if update_rfx
        % Updated prior on random effects
        rfx_prior.pE=m;
        rfx_prior.pC=C;

        % 1. Update estimates of subject random effects
        if verbose, disp('Updating random effects'); end

        for n=1:S.N,
            if verbose
                disp(sprintf('Subject %d out of %d',n,S.N));
            end
            % Start sampling where left off
            subj_mcmc.init=w(:,n);
            Psamp = spm_mci_random (subj_mcmc,rfx_prior,vit,M{n},U{n},Y{n});
            w(:,n) = Psamp(:,end);
        end

    end

    % Update observation noise precision
    if update_obs_noise
        if verbose, disp('Updating observation noise precision'); end
        [noise,M] = spm_mci_obsnoise (w,vit,assign,noise,M,U,Y);
    end
    Ce(:,:,rfx1_its+it)=M{1}.Ce;

    % Store samples
    sm(:,rfx1_its+it)=m;
    sw(:,:,rfx1_its+it)=w;

end

% Define post burn-in
total_its=size(sm,2);
burn_in=round(0.3*total_its);
post_ind=[burn_in+1:total_its];

% Return samples
MCI.sm=sm;
MCI.sw=sw;
try, MCI.sv=v'; catch, MCI.sv=[]; end

% Return posterior means
if update_rfx
    MCI.sm_mean=mean(sm(:,post_ind),2);
    MCI.sw_mean=mean(sw(:,:,post_ind),3);
end
try, MCI.sv_mean=mean(v(post_ind,:)); catch, MCI.sv_mean=[]; end
MCI.post_ind=post_ind;

% Extract posterior means for pinit, pflow, pout
switch assign.init_par,
    case 'fixed',
        MCI.pinit=MCI.sv_mean(:,assign.v_init)';
    case 'random',
        MCI.pinit_sub=MCI.sw_mean(assign.w_init,:);
        MCI.pinit=MCI.sm_mean(assign.w_init,:);
    otherwise
        try, MCI.pinit=MCI.pinit0; catch, MCI.pinit=[]; end
end
switch assign.flow_par,
    case 'fixed',
        MCI.pflow=MCI.sv_mean(:,assign.v_flow)';
    case 'random',
        MCI.pflow_sub=MCI.sw_mean(assign.w_flow,:);
        MCI.pflow=MCI.sm_mean(assign.w_flow,:);
    otherwise
        try, MCI.pflow=MCI.flow0; catch, MCI.pflow=[]; end
end
switch assign.out_par,
    case 'fixed',
        MCI.pout=MCI.sv_mean(:,assign.v_out)';
    case 'random',
        MCI.pout_sub=MCI.sw_mean(assign.w_out,:);
        MCI.pout=MCI.sm_mean(assign.w_out,:);
    otherwise
        try, MCI.pout=MCI.out0; catch, MCI.pout=[]; end
end

MCI.M=M;
MCI.S=S;
MCI.noise=noise;
MCI.assign=assign;
MCI.Ce=Ce;

author: spm


text:
% cell2file   writes a cell array of strings, or number into a file
% ==This function has been obtained modifyng the file str2file.m I found in the FEX.==
% Last Modified: 2008/01/09 (yyyy/mm/dd).
%
% Syntax: cell2file(fid,cell,....) where ... is PARAMETER, VALUE (in pair).
%
% Possible Parameter are: 'EndOfLine' and 'Delimiter'
%
% cell2str can write to a file both sting and numeric format.
%
% You just have to pass the handle to the file
%(so that you can choice your opening attribute) and the cellarray.
%
% Ex:
% A={'prova',[1],'prova';
% [3],'txt',[4]};
%
% fid=fopen('prova.log','a+');
% cell2file(fid,A,'EndOfLine','\r\n');
% fclose(fid);

function cell2file(fid,str,varargin)

id=find(strcmpi('Delimiter',varargin)==1);
if isempty(id)
    delimiter=';';
else
    delimiter=varargin{id+1};
end

id=find(strcmpi('EndOfLine',varargin)==1);
if isempty(id)
    EndOfLine='\n';
else
    EndOfLine=varargin{id+1};
end

for k=1:size(str,1)
    for j=1:size(str,2)
        app=[str{k,j}];
        if isnumeric(app)
            app=num2str(app);
        end
        if not(j==size(str,2))
        fprintf(fid,['%s' delimiter],app);
        else
            fprintf(fid,'%s',app);
        end
    end

    if not(k==size(str,1))
    fprintf(fid,EndOfLine);
    end
end


text:
function mC = SpectralRollOff(signal,windowLength, step, c, fs)
signal = signal / max(abs(signal));
curPos = 1;
L = length(signal);
numOfFrames = (L-windowLength)/step + 1;
H = hamming(windowLength);
m = [0:windowLength-1]';
for (i=1:numOfFrames)
    window = (signal(curPos:curPos+windowLength-1));
    FFT = (abs(fft(window,512)));
    FFT = FFT(1:255);
    totalEnergy = sum(FFT);
    curEnergy = 0.0;
    countFFT = 1;
    while ((curEnergy<=c*totalEnergy) && (countFFT<=255))
        curEnergy = curEnergy + FFT(countFFT);
        countFFT = countFFT + 1;
    end
    mC(i) = ((countFFT-1))/(fs/2);
    curPos = curPos + step;
end


text:
function [rankVal,V,U,S]=matrixRank(A,algorithm)
%%MATRIXRANK Determine the rank of a matrix using a chosen criterion.
%
%INPUTS: A The MXN matrix whose rank is desired.
% algorithm An optional parameter specifying the criterion to use. All
%          methods compare the singular values to a bound. Possible values
%          are:
%          0 (The default if omitted or an empty matrix is passed) Use
%            eps(norm(A,1)) as the bound.
%          1 Use max(size(A))*eps(max(s)), where s is the vector of
%            singular values.
%          2 Use eps()*norm(A,1).
%          3 Use max(size(A))*eps()*max(s).
%
%OUTPUTS: rankVal The rank of matrix A.
%               V The right singular vectors from the SVD used here (The
%                 full SVD is computed, not the 'econ' or 0 options). Given
%                 the rank, one can extract a basis or a nullspace.
%               U The left singular vectors from the SVD.
%               S A diagonal matrix of singular values from the SVD.
%
%The notion of the matrix rank relating to a thresholding of the singular
%values is discussed in Chapter 5.4.1 of [1].
%
%REFERENCES:
%[1] G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed.
%    Baltimore: Johns Hopkins University Press, 2013.
%
%September 2019 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(algorithm))
    algorithm=0;
end

if(nargout>1)
    [U,S,V]=svd(A);
else
    [~,S,~]=svd(A);
end
s=diag(S);

switch(algorithm)
    case 0
        rankVal=sum(s>eps(norm(A,1)));
    case 1
        rankVal=sum(s>max(size(A))*eps(max(s)));
    case 2
        rankVal=sum(s>eps()*norm(A,1));
    case 3
        rankVal=sum(s>max(size(A))*eps()*max(s));
    otherwise
        error('Unknown rank algorithm specified.')
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function triangulation_display ( prefix, node_show, element_show, ...
  neighbor_show )

%*****************************************************************************80
%
%% TRIANGULATION_DISPLAY displays a triangulation.
%
%  Discussion:
%
%    TRIANGULATION_DISPLAY displays a triangulated set of nodes on
%    the MATLAB graphics screen.
%
%  Usage:
%
%    triangulation_display ( 'prefix', node_show, element_show, (neighbor_show) )
%
%    where:
%
%    'prefix' is the common prefix for the node and element files:
%
%    * prefix_nodes.txt,             the node coordinates.
%    * prefix_elements.txt,          the nodes that make up each element.
%    * prefix_element_neighbors.txt, (optional), the neighbors of each element.
%
%    'node_show' indicates the node visibility:
%
%    0: do not show the nodes;
%    1:        show the nodes;
%    2:        show the nodes, and label them.
%
%    'element_show' indicates the element visibility:
%
%    0: do not show the element;
%    1:        show the element;
%    2:        show the element, and color them.
%    3:        show the element, and color them, and label them.
%
%    'neighbor_show' indicates the element neighbor visibility,
%    if the element neighbor file exists:
%
%    0: do not show the element neighbor labels;
%    1:        show the element neighbor labels.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 April 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGULATION_DISPLAY\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read a node dataset of NODE_NUM points in 2 dimensions.\n' );
  fprintf ( 1, '  Read an associated triangulation dataset of ELEMENT_NUM\n' );
  fprintf ( 1, '  elements using 3, 4 or 6 nodes.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Display an image of the data on the MATLAB graphics screen.\n' );
%
%  First argument is the file prefix.
%
  if ( nargin < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGULATION_DISPLAY:\n' );
    prefix = input ( '  Enter the file prefix:  ' );
  end
%
%  Create the file names.
%
  node_filename = strcat ( prefix, '_nodes.txt' );
  element_filename = strcat ( prefix, '_elements.txt' );
  element_neighbor_filename = strcat ( prefix, '_element_neighbors.txt' );
%
%  Second argument is node visibility.
%
  if ( nargin < 2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Options for node visibility:\n' );
    fprintf ( 1, '  0: do not show the nodes;\n' );
    fprintf ( 1, '  1:        show the nodes;\n' );
    fprintf ( 1, '  2:        show the nodes, and label them.\n' );
    node_show = input ( '  Enter the node visibility option:  ' );
  end
%
%  Third argument is element visibility.
%
  if ( nargin < 3 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Options for element visibility:\n' );
    fprintf ( 1, '  0: do not show the elements;\n' );
    fprintf ( 1, '  1:        show the elements;\n' );
    fprintf ( 1, '  2:        show the elements, and color them.\n' );
    fprintf ( 1, '  3:        show the elements, and color them, and label them.\n' );
    element_show = input ( '  Enter the element visibility option:  ' );
  end
%
%  Fourth argument is neighbor visibility.
%
  if ( nargin < 4 && file_exist ( element_neighbor_filename ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Options for element neighbor visibility:\n' );
    fprintf ( 1, '  0: do not show the element neighbor labels;\n' );
    fprintf ( 1, '  1:        show the element neighbor labels;\n' );
    neighbor_show = input ( '  Enter the element neighbor visibility option:  ' );
  end
%
%  Read the node data.
%
  [ dim_num, node_num ] = r8mat_header_read ( node_filename );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the header of "%s".\n', node_filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension DIM_NUM = %d\n', dim_num );
  fprintf ( 1, '  Number of nodes NODE_NUM  = %d\n', node_num );

  node_xy = r8mat_data_read ( node_filename, dim_num, node_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the data in "%s".\n', node_filename );

  r8mat_transpose_print_some ( dim_num, node_num, node_xy, 1, 1, 2, 5, ...
    '  Initial portion of NODE_XY:' );
%
%  Read the element data.
%
  [ element_order, element_num ] = i4mat_header_read ( element_filename );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the header of "%s".\n', element_filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Element order      = %d\n', element_order );
  fprintf ( 1, '  Number of elements = %d\n', element_num );

  element_node = i4mat_data_read ( element_filename, ...
    element_order, element_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the data in "%s".\n', element_filename );

  i4mat_transpose_print_some ( element_order, element_num, element_node, ...
    1, 1, element_order, 5, '  Initial portion of ELEMENT_NODE:' );
%
%  If the neighbor file exists, read the neighbor data.
%
  if ( file_exist ( element_neighbor_filename ) && 0 < neighbor_show )

    [ three, element_num ] = i4mat_header_read ( element_neighbor_filename );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Read the header of "%s".\n', element_neighbor_filename );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Element order      = %d\n', element_order );
    fprintf ( 1, '  Number of elements = %d\n', element_num );

    element_neighbor = i4mat_data_read ( element_neighbor_filename, ...
      3, element_num );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Read the data in "%s".\n', element_neighbor_filename );

    i4mat_transpose_print_some ( 3, element_num, element_neighbor, ...
      1, 1, element_order, 5, '  Initial portion of ELEMENT_NEIGHBOR:' );

  else

    neighbor_show = 0;
    element_neighbor = [];

  end
%
%  Detect and correct zero-based node indexing.
%
  element_node = mesh_base_one ( node_num, element_order, element_num, ...
    element_node );
%
%  Display a plot.
%
  title_string = s_escape_tex ( prefix );

  triangulation_display_matlab ( title_string, node_num, node_xy, element_order, ...
    element_num, element_node, element_neighbor, node_show, element_show, ...
    neighbor_show );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGULATION_DISPLAY\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

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

  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:
%
%    21 February 2004
%
%  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 value = file_exist ( file_name )

%*****************************************************************************80
%
%% FILE_EXIST reports whether a file exists.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character FILE_NAME, the name of the file.
%
%    Output, logical FILE_EXIST, is TRUE if the file exists.
%
  fid = fopen ( file_name );

  if ( fid == -1 )
    value = 0;
  else
    fclose ( fid );
    value = 1;
  end

  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:
%
%    31 December 2006
%
%  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 value = i4_modp ( i, j )

%*****************************************************************************80
%
%% I4_MODP returns the nonnegative remainder of I4 division.
%
%  Discussion:
%
%    If
%      NREM = I4_MODP ( I, J )
%      NMULT = ( I - NREM ) / J
%    then
%      I = J * NMULT + NREM
%    where NREM is always nonnegative.
%
%    The MOD function computes a result with the same sign as the
%    quantity being divided.  Thus, suppose you had an angle A,
%    and you wanted to ensure that it was between 0 and 360.
%    Then mod(A,360) would do, if A was positive, but if A
%    was negative, your result would be between -360 and 0.
%
%    On the other hand, I4_MODP(A,360) is between 0 and 360, always.
%
%  Example:
%
%        I     J     MOD  I4_MODP    Factorization
%
%      107    50       7       7    107 =  2 *  50 + 7
%      107   -50       7       7    107 = -2 * -50 + 7
%     -107    50      -7      43   -107 = -3 *  50 + 43
%     -107   -50      -7      43   -107 =  3 * -50 + 43
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 1999
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the number to be divided.
%
%    Input, integer J, the number that divides I.
%
%    Output, integer VALUE, the nonnegative remainder when I is
%    divided by J.
%
  if ( j == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4_MODP - Fatal error!\n' );
    fprintf ( 1, '  Illegal divisor J = %d\n', j );
    error ( 'I4_MODP - Fatal error!' );
  end

  value = mod ( i, j );

  if ( value < 0 )
    value = value + abs ( j );
  end

  return
end
function value = i4_wrap ( ival, ilo, ihi )

%*****************************************************************************80
%
%% I4_WRAP forces an integer to lie between given limits by wrapping.
%
%  Example:
%
%    ILO = 4, IHI = 8
%
%    I   Value
%
%    -2     8
%    -1     4
%     0     5
%     1     6
%     2     7
%     3     8
%     4     4
%     5     5
%     6     6
%     7     7
%     8     8
%     9     4
%    10     5
%    11     6
%    12     7
%    13     8
%    14     4
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 October 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer IVAL, an integer value.
%
%    Input, integer ILO, IHI, the desired bounds for the integer value.
%
%    Output, integer I4_WRAP, a "wrapped" version of IVAL.
%
  jlo = min ( ilo, ihi );
  jhi = max ( ilo, ihi );

  wide = jhi - jlo + 1;

  if ( wide == 1 )
    value = jlo;
  else
    value = jlo + i4_modp ( ival - jlo, wide );
  end

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

%*****************************************************************************80
%
%% I4MAT_DATA_READ reads data from an I4MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer M, N, the number of rows and columns in the data.
%
%    Output, integer TABLE(M,N), the point coordinates.
%
  table = zeros ( m, n );
%
%  Build up the format string for reading M real numbers.
%
  string = ' ';

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

  input_unit = fopen ( input_filename );

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

  i = 0;

  while ( i < n )

    line = fgets ( input_unit );

    if ( line == -1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'I4MAT_DATA_READ - Error!\n' );
      fprintf ( 1, '  End of input while reading data.\n' );
      error ( 'I4MAT_DATA_READ - Error!' );
    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 ] = i4mat_header_read ( input_filename )

%*****************************************************************************80
%
%% I4MAT_HEADER_READ reads the header from an I4MAT file.
%
%  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, 'I4MAT_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, 'I4MAT_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 i4mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% I4MAT_TRANSPOSE_PRINT_SOME prints some of an I4MAT, transposed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 June 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    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 );

  for i2lo = max ( ilo, 1 ) : incx : min ( ihi, m )

    i2hi = i2lo + incx - 1;
    i2hi = min ( i2hi, m );
    i2hi = min ( i2hi, ihi );

    inc = i2hi + 1 - i2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row: ' );
    for i = i2lo : i2hi
      fprintf ( 1, '%7d  ', i );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col\n' );
    fprintf ( 1, '\n' );

    j2lo = max ( jlo, 1 );
    j2hi = min ( jhi, n );

    for j = j2lo : j2hi

      fprintf ( 1, '%5d  ', j );
      for i2 = 1 : inc
        i = i2lo - 1 + i2;
        fprintf ( 1, '%7d  ', a(i,j) );
      end
      fprintf ( 1, '\n' );

    end

  end

  return
end
function element_node = mesh_base_one ( node_num, element_order, ...
  element_num, element_node )

%*****************************************************************************80
%
%% MESH_BASE_ONE ensures that the element definition is one-based.
%
%  Discussion:
%
%    The ELEMENT_NODE array contains nodes indices that form elements.
%    The convention for node indexing might start at 0 or at 1.
%    Since a MATLAB program will naturally assume a 1-based indexing, it is
%    necessary to check a given element definition and, if it is actually
%    0-based, to convert it.
%
%    This function attempts to detect 0-based node indexing and correct it.
%
%    Thanks to Feifei Xu for pointing out that I was subtracting 1 when I
%    should have been adding 1!  29 November 2012.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 November 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer ELEMENT_ORDER, the order of the elements.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input/output, integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM), the element
%    definitions.
%
  node_min = min ( min ( element_node(1:element_order,1:element_num) ) );
  node_max = max ( max ( element_node(1:element_order,1:element_num) ) );

  if ( node_min == 0 && node_max == node_num - 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MESH_BASE_ONE:\n' );
    fprintf ( 1, '  The element indexing appears to be 0-based!\n' );
    fprintf ( 1, '  This will be converted to 1-based.\n' );
    element_node(1:element_order,1:element_num) = ...
      element_node(1:element_order,1:element_num) + 1;
  elseif ( node_min == 1 && node_max == node_num )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MESH_BASE_ONE:\n' );
    fprintf ( 1, '  The element indexing appears to be 1-based!\n' );
    fprintf ( 1, '  No conversion is necessary.\n' );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MESH_BASE_ONE - Warning!\n' );
    fprintf ( 1, '  The element indexing is not of a recognized type.\n' );
  end

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

%*****************************************************************************80
%
%% R8MAT_DATA_READ reads data from an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2006
%
%  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.
%
  table = zeros ( m, n );
%
%  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

  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.
%
%  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 r8mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, real 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 = 5;

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

  for i2lo = max ( ilo, 1 ) : incx : min ( ihi, m )

    i2hi = i2lo + incx - 1;
    i2hi = min ( i2hi, m );
    i2hi = min ( i2hi, ihi );

    inc = i2hi + 1 - i2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row: ' );
    for i = i2lo : i2hi
      fprintf ( 1, '%7d       ', i );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col\n' );

    j2lo = max ( jlo, 1 );
    j2hi = min ( jhi, n );

    for j = j2lo : j2hi

      fprintf ( 1, '%5d ', j );
      for i2 = 1 : inc
        i = i2lo - 1 + i2;
        fprintf ( 1, '%12f', a(i,j) );
      end
      fprintf ( 1, '\n' );

    end

  end

  return
end
function s2 = s_escape_tex ( s1 )

%*****************************************************************************80
%
%% S_ESCAPE_TEX de-escapes TeX escape sequences.
%
%  Discussion:
%
%    In particular, every occurrence of the characters '\', '_',
%    '^', '{' and '}' will be replaced by '\\', '\_', '\^',
%    '\{' and '\}'.  A TeX interpreter, on seeing these character
%    strings, is then likely to return the original characters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S1, the string to be de-escaped.
%
%    Output, string S2, a copy of the string, modified to avoid TeX escapes.
%
  s1_length = length ( s1 );

  s1_pos = 0;
  s2_pos = 0;
  s2 = [];

  while ( s1_pos < s1_length )

    s1_pos = s1_pos + 1;

    if ( s1(s1_pos) == '\' || ...
         s1(s1_pos) == '_' || ...
         s1(s1_pos) == '^' || ...
         s1(s1_pos) == '{' || ...
         s1(s1_pos) == '}' )
      s2_pos = s2_pos + 1;
      s2 = strcat ( s2, '\' );
    end

    s2_pos = s2_pos + 1;
    s2 = strcat ( s2, s1(s1_pos) );

  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 LEN, 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
function triangulation_display_matlab ( title_string, node_num, node_xy, ...
  element_order, element_num, element_node, element_neighbor, node_show, ...
  element_show, neighbor_show )

%*****************************************************************************80
%
%% TRIANGULATION_DISPLAY_MATLAB displays a triangulation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 April 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string TITLE_STRING, a title for the plot.
%
%    Input, integer NODE_NUM, the number of points.
%
%    Input, real NODE_XY(2,NODE_NUM), the nodes.
%
%    Input, integer ELEMENT_ORDER, the order of the elements.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM), lists,
%    for each element, the indices of the points that form the vertices
%    of the element.
%
%    Input, integer ELEMENT_NEIGHBOR(3,ELEMENT_NUM).  ELEMENT_NEIGHBOR(I,J)
%    is the neighbor adjacent to the size opposite node ELEMENT_NODE(I,J),
%    or -1 if no such neighbor.
%
%    Input, logical NODE_SHOW:
%    0, do not show the nodes.
%    1, show the nodes.
%    2, show the nodes, and label them.
%
%    Input, logical ELEMENT_SHOW,
%    0, do not show the elements.
%    1, show the elements.
%    2, show the elements, and color them.
%    3, show the elements, and color them, and label them.
%
%    Input, logical NEIGHBOR_SHOW,
%    0, do not show the element neighbor labels.
%    1, show the element neighbor labels.
%

%
%  Clear the graphics screen.
%
  clf;
%
%  Get the scale.
%
  xy_min(1) = min ( node_xy(1,:) );
  xy_max(1) = max ( node_xy(1,:) );

  xy_min(2) = min ( node_xy(2,:) );
  xy_max(2) = max ( node_xy(2,:) );

  xy_range(1:2) = xy_max(1:2) - xy_min(1:2);

  margin = 0.025 * max ( xy_range(1), xy_range(2) );

  x_min = xy_min(1) - margin;
  x_max = xy_max(1) + margin;
  y_min = xy_min(2) - margin;
  y_max = xy_max(2) + margin;
%
%  Draw the element faces.
%
  if ( 2 <= element_show )

    face_color = 'g';

    for element = 1 : element_num
      p = element_node(1:3,element);
      patch ( node_xy(1,p), node_xy(2,p), face_color );
    end

  end
%
%  Draw the element lines.
%
  if ( 1 <= element_show )

    line_color = 'r';

    for element = 1 : element_num
      p = [ element_node(1:3,element)', element_node(1,element) ];
      line ( 'XData', node_xy(1,p), 'YData', node_xy(2,p), ...
        'Color', line_color, 'LineWidth', 2 );
    end

  end
%
%  Label the elements
%
  if ( 3 <= element_show )

    for element = 1 : element_num

      ave_x = 0.0;
      ave_y = 0.0;

      for i = 1 : element_order
        node = element_node(i,element);
        ave_x = ave_x + node_xy(1,node);
        ave_y = ave_y + node_xy(2,node);
      end

      ave_x = ave_x / element_order;
      ave_y = ave_y / element_order;

      text ( ave_x, ave_y, num2str ( element ) );

    end

  end
%
%  Draw the nodes.
%
  hold on

  if ( 1 <= node_show )
    plot ( node_xy(1,:), node_xy(2,:), 'o', 'markersize', 4, ...
      'markerfacecolor', 'black' );
  end
%
%  Label the nodes.
%
  if ( 2 <= node_show )
    if ( node_num <= 100 )
      abit = margin;
    else
      abit = 0.0;
    end
    for j = 1 : node_num
      text ( node_xy(1,j) + abit, node_xy(2,j) + abit, num2str(j) );
    end
  end
%
%  Label the element neighbors.
%
  if ( 1 <= neighbor_show )
    for j = 1 : element_num
      for i = 1 : 3
        np0 = element_node ( i, j );
        ip1 = i4_wrap ( i + 1, 1, 3 );
        np1 = element_node ( ip1, j );
        ip2 = i4_wrap ( i + 2, 1, 3 );
        np2 = element_node ( ip2, j );
        x = 0.40 * node_xy(1,np1) + 0.40 * node_xy(1,np2) + 0.20 * node_xy(1,np0);
        y = 0.40 * node_xy(2,np1) + 0.40 * node_xy(2,np2) + 0.20 * node_xy(2,np0);
        text ( x, y, num2str ( element_neighbor(i,j) ) );
      end
    end
  end
%
%  Label the plot.
%
  xlabel ( '--X axis--' )
  ylabel ( '--Y axis--' )
  title ( title_string )

  axis ( [ x_min, x_max, y_min, y_max ] );
  axis equal

  hold off

  return
end

author: johannesgerer


text:
% VL_DEMO_GMM_2D_RAND  Demonstrate clustering points with a GMM

%% Create a random set of points
numPoints = 5000 ;
dimension = 2 ;
numClusters = 20 ;
data = rand(dimension, numPoints) ;

%% Learn a GMM: fit the points at maximum likelihood
vl_twister('state',0) ;
[means, covariances, priors] = ...
    vl_gmm(data, numClusters, ...
           'MaxNumIterations', 1000, ...
           'Verbose') ;

figure(1) ; clf ; hold on
plot(data(1,:),data(2,:),'r.');
for i=1:numClusters
  vl_plotframe([means(:,i)' covariances(1,i) 0 covariances(2,i)], ...
               'Color','blue','LineWidth',2);
end

title('GMM: Gaussian mixture initialized choosing random points') ;
axis equal ; axis off ;
vl_demo_print('gmm_2d_rand',0.6);

author: jianxiongxiao

repo: ProfXkit


text:
function d =  generate(a)

x=randn(a.l,length(a.mean));
n=length(a.mean);

if size(a.cov,2)==1
  cov=eye(n)*a.cov;   %% make cov matrix with only diagonal elems (all same)
end;
if size(a.cov,1)==1
  cov=diag(a.cov);   %% make cov matrix with only diagonal elems (different)
end

for i=1:a.l
  if length(a.cov)==1
    x(i,:)= (x(i,:)*sqrt(a.cov)) + a.mean; continue;
  else
    if size(a.cov,1)==1
      x(i,:)= (x(i,:).* sqrt(a.cov)) + a.mean; continue;
    else
      x(i,:)= x(i,:) * a.cov^0.5 + a.mean;
    end
  end
end

d=data(['gauss' ' l=' num2str(a.l)],x,[]);


text:
function [xs,phid,wlsden,timeiter,niter] = wls_gcd(A, W, yy, xhat, ng, niter, tol)
%	weighted least squares minimized by
%	grouped coordinate descent algorithm
%	Input
%		A	[nn,np]		system matrix
%		W	[nn,nn]		inverse data covariance
%		yy	[nn,1]		noisy data
%		xhat	[np,1] or [nx,ny]	initial estimate
%		ng			# groups, or, if array, groups
%		tol			stopping criterion for convergence
%		niter			max
%
%	Output
%		xs	[np,niter+1]	estimates each iteration
%		phid	[niter+1,1]	objective decrease
%		wlsden	[np,1]		curvature of surrogate
%		timeiter [niter+1,1]	CPU time each iteration
%		niter			# total iterations
%
% Saowapak Sotthivirat
% 8/15/00

[nx,ny] = size(xhat);

np = ncol(A);
yy = yy(:);
xhat = xhat(:);
xs(:,1) = xhat;

%
%	build groups
%
if numel(ng) == 1
	groups = zeros(np, ng);
	for ig = 1:ng
		gg(:,ig) = [ig:ng:np]';
		groups(gg(:,ig),ig) = ones(size(gg(:,ig)));
	end
elseif numel(ng) == 2
	groups	= group2d(nx, ny,ng,ones(ng(1),ng(2)),0);
	ng	= ncol(groups);
	for ig = 1:ng
		gg(:,ig) = find(groups(:,ig).*[1:np]');
	end
else
	groups	= ng;
	ng	= ncol(groups);
end

%	precompute quadratic part of denominator
%	change '1' to 'ig' for overlapping groups
wlsden = zeros(np, 1);
for ig = 1:ng
	gr = groups(:,ig);
	ggi = gg(:,ig);

	t = A * gr;	% faster than: sum(A(:,gg)',1)';
	t = W * t;	% only slightly slower than: diag(W) .* t;
			% faster than either: A(:,gg)'*t or A' * t !
	t = t' * A;
	wlsden(ggi) = t(ggi)';
end

wii = diag(W);
res = A * xhat - yy;
phi0 = res'*W*res/2;
ii = 1;
phid(ii) = 0;
timeiter(ii) = 0;
t0 = cputime;

while phid(ii) < tol & ii <= niter

	for ig = 1:ng
		ggi = gg(:,ig);
		xold = xhat(ggi);
		xnew = xold;
		ader = ((W*res)'*A(:,ggi))';
		diff = ader./wlsden(ggi);
		xnew = xnew - diff;
		xhat(ggi) = max(0,xnew);
		res = res - A(:,ggi) * diff;

	end
	ii = ii+1;
	xs(:,ii) = xhat;
	phid(ii) = phi0 - res'*W*res/2;
	timeiter(ii) = cputime - t0;
end

niter = ii;

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function legendre_symbol_test ( )

%*****************************************************************************80
%
%% LEGENDRE_SYMBOL_TEST tests LEGENDRE_SYMBOL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2008
%
%  Author:
%
%    John Burkardt
%
  ntest = 4;
  ptest = [ 7, 11, 13, 17 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'LEGENDRE_SYMBOL_TEST\n' );
  fprintf ( 1, '  LEGENDRE_SYMBOL computes the Legendre\n' );
  fprintf ( 1, '  symbol (Q/P) which records whether Q is \n' );
  fprintf ( 1, '  a quadratic residue modulo the prime P.\n' );

  for i = 1 : ntest
    p = ptest(i);
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Legendre Symbols for P = %d\n', p );
    fprintf ( 1, '\n' );
    for q = 0 : p
      l = legendre_symbol ( q, p );
      fprintf ( 1, '  %6d  %6d  %6d\n', p, q, l );
    end
  end

  return
end


text:
function evaluation(norm_pred_list,gt_dir,setting_name,setting_class,legend_name)
load(gt_dir);
if ~exist(sprintf('./plot/baselines/Val/%s/%s',setting_class,legend_name),'dir')
    mkdir(sprintf('./plot/baselines/Val/%s/%s',setting_class,legend_name));
end
IoU_thresh = 0.5;
event_num = 61;
thresh_num = 1000;
org_pr_cruve = zeros(thresh_num,2);
count_face = 0;

for i = 1:event_num
    img_list = file_list{i};
    gt_bbx_list = face_bbx_list{i};
    pred_list = norm_pred_list{i};
    sub_gt_list = gt_list{i};
    img_pr_info_list = cell(length(img_list),1);

    fprintf('%s, current event %d\n',setting_name,i);
    for j = 1:length(img_list)
        gt_bbx = gt_bbx_list{j};
        pred_info = pred_list{j};
        keep_index = sub_gt_list{j};
        count_face = count_face + length(keep_index);

        if isempty(gt_bbx) || isempty(pred_info)
            continue;
        end
        ignore = zeros(size(gt_bbx,1),1);
        if ~isempty(keep_index)
            ignore(keep_index) = 1;
        end

        [pred_recall, proposal_list] = image_evaluation(pred_info, gt_bbx, ignore, IoU_thresh);

        img_pr_info = image_pr_info(thresh_num, pred_info, proposal_list, pred_recall);
        img_pr_info_list{j} = img_pr_info;
    end
    for j = 1:length(img_list)
        img_pr_info = img_pr_info_list{j};
        if ~isempty(img_pr_info)
            org_pr_cruve(:,1) = org_pr_cruve(:,1) + img_pr_info(:,1);
            org_pr_cruve(:,2) = org_pr_cruve(:,2) + img_pr_info(:,2);
        end
    end
end
pr_cruve = dataset_pr_info(thresh_num, org_pr_cruve, count_face);
save(sprintf('./plot/baselines/Val/%s/%s/wider_pr_info_%s_%s.mat',setting_class,legend_name,legend_name,setting_name),'pr_cruve','legend_name','-v7.3');
end

function [pred_recall,proposal_list] = image_evaluation(pred_info, gt_bbx, ignore, IoU_thresh)
    pred_recall = zeros(size(pred_info,1),1);
    recall_list = zeros(size(gt_bbx,1),1);
    proposal_list = zeros(size(pred_info,1),1);
    proposal_list = proposal_list + 1;
    pred_info(:,3) = pred_info(:,1) + pred_info(:,3);
    pred_info(:,4) = pred_info(:,2) + pred_info(:,4);
    gt_bbx(:,3) = gt_bbx(:,1) + gt_bbx(:,3);
    gt_bbx(:,4) = gt_bbx(:,2) + gt_bbx(:,4);
    for h = 1:size(pred_info,1)
        overlap_list = boxoverlap(gt_bbx, pred_info(h,1:4));
        [max_overlap, idx] = max(overlap_list);
        if max_overlap >= IoU_thresh
            if (ignore(idx) == 0)
                recall_list(idx) = -1;
                proposal_list(h) = -1;
            elseif (recall_list(idx)==0)
                recall_list(idx) = 1;
            end
        end
        r_keep_index = find(recall_list == 1);
        pred_recall(h) = length(r_keep_index);
    end
end

function img_pr_info = image_pr_info(thresh_num, pred_info, proposal_list, pred_recall)
    img_pr_info = zeros(thresh_num,2);
    for t = 1:thresh_num
        thresh = 1-t/thresh_num;
        r_index = find(pred_info(:,5)>=thresh,1,'last');
        if (isempty(r_index))
            img_pr_info(t,2) = 0;
            img_pr_info(t,1) = 0;
        else
            p_index = find(proposal_list(1:r_index) == 1);
            img_pr_info(t,1) = length(p_index);
            img_pr_info(t,2) = pred_recall(r_index);
        end
    end
end

function pr_cruve = dataset_pr_info(thresh_num, org_pr_cruve, count_face)
    pr_cruve = zeros(thresh_num,2);
    for i = 1:thresh_num
        pr_cruve(i,1) = org_pr_cruve(i,2)/org_pr_cruve(i,1);
        pr_cruve(i,2) = org_pr_cruve(i,2)/count_face;
    end
end

author: peiyunh

repo: tiny

save_path: github-repos/MATLAB/peiyunh-tiny


text:
%%Given the ned position (L,lon, h) computes the local radious and gravity
function [Rn, Re, gn, sL, cL, WIE_E]=geoparam(pos_n)
%Constants
WIE_E=7292115e-11;    %earth's rotaion rate
SM_AXIS=6378137;
E_SQR=0.00669437999014;
NORMAL_GRV=9.7803253359;
GRV_CONS=0.00193185265241;
FLATTENING=0.00335281066475;
M_FAKTOR=0.00344978650684;

sL=sin(pos_n(1));
cL=cos(pos_n(1));
h=pos_n(3);
Rn=6335439.327292829/(sqrt(1.0-E_SQR*sL*sL)*(1.0-E_SQR*sL*sL));
Re=SM_AXIS/(sqrt(1.0-E_SQR*sL*sL));
g1=NORMAL_GRV*(1+GRV_CONS*sL*sL)/(sqrt(1.0-E_SQR*sL*sL));
gn=g1*(1.0-(2.0/SM_AXIS)*(1.0+FLATTENING+M_FAKTOR-2.0*FLATTENING*sL*sL)*h+3.0*h*h/SM_AXIS/SM_AXIS);

repo: ekfmonoslam


text:
%SETNAN Set nan values.
%
%   [ X ] = SETNAN( X, VAL ) Sets all nan values in the numeric array
%   X to VAL (default 0).

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function prob_thin = probThin(BandCirrus)
%PROBCIRRUS
    prob_thin = BandCirrus./400;
    clear BandCirrus;
%     prob_thin(prob_thin>1)=1;
    prob_thin(prob_thin<0)=0;
end

author: GERSL

repo: Fmask

save_path: github-repos/MATLAB/GERSL-Fmask


text:
function [sstate, si] = tapas_sem_hier_get_stored_state(data, model, ...
    inference, state)
%% Get the state that should be stored
%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

sstate = struct('graph', {cell(4, 1)}, 'llh', {cell(4, 1)}, 'v', []);

order = state.graph{2}.u.temperature_ordering;
sstate.graph{2} = state.graph{2}.y(:, order);

sstate.llh{1} = state.llh{1}(:, order);

sstate.v = state.v(:, order);

if state.nsample < inference.nburnin
    si = mod(state.nsample, inference.ndiag) + 1;
else
    si = state.nsample - inference.nburnin + 1;
end

end

author: translationalneuromodeling

repo: tapas