text:
% EEG_MULTIEEGPLOT - Produce an EEGPLOT of a the average of an epoched dataset
%                  (with optional pre-labelling of specific trials).
% Usage:
%   >> eeg_multieegplot( data,trialrej, elecrej, ...
%                                'key1', value, 'key2', value ... );
% Inputs:
%   data        - input data (channels x points or channels x points x trials).
%   trialrej    - array of 0s and 1s (depicting rejected trials) (size sweeps)
%   elecrej     - array of 0s and 1s (depicting electrodes rejected in
%                 all trials) (size nbelectrodes x sweeps )
%   oldtrialrej - array of 0s and 1s (depicting rejected trials) (size sweeps)
%   oldelecrej  - array of 0s and 1s (depicting electrodes rejected in
%                 all trials) (size nbelectrodes x sweeps )
%
% Note: 1) {'Key', value } Arguments are passed on to EEGPLOT
%       2) To ignore previous rejections simply set 'oldtrialrej' and
%          'oldelecrej' arguments to empty ([]).
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: EEGPLOT, EEGPLOT2EVENT, EEGPLOT2TRIAL, EEGLAB

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% 01-25-02 reformated help & license -ad
% 03-07-02 corrected help -ad

function eeg_multieegplot( data, rej, rejE, oldrej, oldrejE, varargin);

if nargin < 1
   help eeg_multieegplot;
   return;
end
if nargin < 4
   oldrej = [];
end
if nargin < 5
   oldrejE = [];
end
if ~exist('command')
	command = '';
end
chans = size(data,1);
pnts  = size(data,2);

colnew = [0.8 0.8 1];
colold = [0.8 1 0.8];

if ndims(data) > 2 % --------------- considering epoched EEG data
   rejeegplot = [];
   % concatenate old rejection if not empty
   % --------------------------------------
   if ~isempty( oldrej )
   		oldrej 	 = find( oldrej > 0);
   		oldrejE  = oldrejE(:, oldrej)';
   		rejeegplot = trial2eegplot( oldrej, oldrejE, pnts, colold);
   end;
   % convert for eegplot
   % -------------------
   if ~isempty(	rej )
   		rej    = find( rej > 0);
   		rejE   = rejE(:,rej)';
   		secondrejeegplot = trial2eegplot( rej, rejE, pnts, colnew); % see bottom of code for this function
   		rejeegplot = [ rejeegplot' secondrejeegplot' ]';

   		% remove duplicates
   		% -----------------
		%[tmp I] = unique_bc( rejeegplot(:,1) );
		%rejeegplot = rejeegplot(I,:);
   end
else % ---------------------------------------- considering continuous EEG
   % for continuous EEG, electrodes (rejE and oldrejE) are not considered yet
   % because there would be a format problem (these rejection are stored in
   % the event array).

   rejeegplot = [];
   if ~isempty(rej) %assuming nrejection x 2 (2 = begin end)
      s = size(rej, 1);
      rejeegplot = [ rej(3:4) colnew*ones(s,1)  zeros(s, chans) ];
   end;

   % pooling previous rejections
   if ~isempty(oldrej)
       tmp = [ oldrej(:, 3:4) colold*ones(s,1) zeros(size(oldrej,1), chans) ];
       rejeegplot = [rejeegplot; tmp];
   end
end

if isempty(varargin)
	eegplot(data, 'winlength', 5, 'position', [100 300 800 500], 'winrej', rejeegplot, 'xgrid', 'off');
else
	eegplot(data, 'winlength', 5, 'position', [100 300 800 500], 'winrej', rejeegplot, 'xgrid', 'off', varargin{:} );
end;
return;

% convert eeglab format to eeplot format of rejection window
% ----------------------------------------------------------
function rejeegplot = trial2eegplot( rej, rejE, pnts, color)
   	rejeegplot = zeros(length(rej), size(rejE,2)+5);
   	rejeegplot(:, 6:end) = rejE;
   	rejeegplot(:, 1) = (rej(:)-1)*pnts;
   	rejeegplot(:, 2) = rej(:)*pnts-1;
   	rejeegplot(:, 3:5) = ones(size(rejeegplot,1),1)*color;
return

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function status = GE_convertVolume(inBaseName,outName,numImages)
%
% status = GE_convertVolume(inDir,outName,[numImages])
%
%GE_convertVolume
%
% Version 3.1
%
% Converts a series of GE slices into Analyze format for SPM.
% inDir is the name of the directory containing the first file of
% the series, e.g. 003
% outName is used for naming the Analyze files
% status = 1 if there is an error.
%
% Modified to use spm functions
%
% Souheil J. Inati
% Dartmouth College
% September 2001
% souheil.inati@dartmouth.edu
%

status = 1;

if (nargin < 2)
        error('Not enough input arguments.')
        return
end

% Create the name of the first file in inDir
firstfile = fullfile('',[inBaseName,'.001']);

% Read the Header from the first file
disp(['Reading first file: ',firstfile]);
[su_hdr,ex_hdr,se_hdr,im_hdr,pix_hdr,im_offset] = GE_readHeader(firstfile);

% The image Dimensions
Nx =  im_hdr.imatrix_X; % X Voxels
Ny =  im_hdr.imatrix_Y; % Y Voxels
if(exist('numImages','var') & ~isempty(numImages))
    Nz = numImages;
else
    Nz =  im_hdr.slquant;   % Z Voxels
end
volSize = [Nx Ny Nz];

disp(['Reading volume (volsize: ',num2str(volSize),')...']);

% Read in the volume
[imageVol, lastfile] = GE_readVolume(inBaseName, volSize, pix_hdr.img_depth, im_offset);

disp(['Converting...']);
% Is the first image the first or the last?
if (se_hdr.end_loc - se_hdr.start_loc) > 0
  scandir = 1;
else
  scandir = -1;
end

% Compute the M matrix
M = GE_createSPMmat(im_hdr,scandir);

% Create the SPM volume
V.fname = outName;
V.dim = [volSize spm_type('int16')]; % short ints
V.mat = M;
V.pinfo = [1 0 0]';
V = spm_create_image(V);

% Write out the SPM Volume
% Don't use spm_write_vol because that sets a scale factor.
for i = 1:Nz
  V = spm_write_plane(V,squeeze(imageVol(:,:,i)),i);
end

% Done with vol
fprintf('Wrote %s.img.\n',outName);

% Done
fprintf('Conversion Finished. \n\n');

status = 0;

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function W = lssvc(A, TYPE, PAR, C)
%LSSVC Least-Squares Support Vector Classifier
%
%       W = lssvc(A,TYPE,PAR,C);
%
% INPUT
%   A       dataset
%   TYPE    Type of the kernel (optional; default: 'p')
%   PAR     Kernel parameter (optional; default: 1)
%   C       Regularization parameter (optional; default: 1)
%
% OUTPUT
%   W       Mapping: Least-Squares Support Vector Classifier
%
% DESCRIPTION
% Optimizes a least-squares support vector classifier for the dataset A by
% quadratic programming. The classifier can be of one of the types
% as defined by PROXM. Default is linear (TYPE = 'p', PAR = 1). The
% regularization parameter C allows for seeting the level of overfitting.
% A smaller value for C allows for more class overlap. Default C = 10.
%
% NOTE
% This implementation uses the LS-SVMLab toolbox. This toolbox is available
% from http://www.esat.kuleuven.ac.be/sista/lssvmlab/. Make sure the path
% to LS-SVMLab is set in Matlab.
%
% See also MAPPINGS, DATASETS, PROXM

% Copyright: L.J.P. van der Maaten, l.vandermaaten@micc.unimaas.nl
% MICC-IKAT, Maastricht University, Maastricht, The Netherlands

    name = 'LSSVC';

    % Check whether LS-SVMLab is installed and set in the path
    if ~exist('trainlssvm.m', 'file')
        error('Could not find LS-SVMLab. Make sure LS-SVMLab is installed and set in the Matlab path.');
    end

    % Perform LS-SVM training on dataset A
    if nargin <= 1 || (nargin >= 2 && ~isa(TYPE, 'prmapping'))

        % Set LS-SVM parameters
        if ~exist('TYPE', 'var') || ~isstruct(TYPE)
            if nargin < 2 || isempty(TYPE), data.kernel = 'p';    else data.kernel = TYPE; end       % kernel function
            if nargin < 3 || isempty(PAR),  data.par = 1;         else data.par = PAR; end           % kernel parameters
            if nargin < 4 || isempty(C),    data.gamma = 10;      else data.gamma = C; end           % regularization parameter
            data.par2 = -1;

            % LS-SVMLab uses different naming for kernels; correct this
            switch lower(data.kernel)
                case {'p', 'polynomial'},       data.kernel = 'poly_kernel'; data.par2 = 1;
                case {'h', 'homogeneous'},      data.kernel = 'poly_kernel'; data.par2 = 0;
                case {'e', 'exponential'},      error('This kernel is not supported by LS-SVMLab.');
                case {'r', 'radial_basis'},     data.kernel = 'RBF_kernel'; data.par = data.par^2;          % LS-SVMLab uses sigma^2 as input
                case {'s', 'sigmoid'},          data.kernel = 'MLP_kernel';
                case {'d', 'distance'},         error('This kernel is not supported by LS-SVMLab.');
                case {'m', 'minkowski'},        error('This kernel is not supported by LS-SVMLab.');
                case {'c', 'city-block'},       error('This kernel is not supported by LS-SVMLab.');
                case {'o', 'cosine'},           error('This kernel is not supported by LS-SVMLab.');
                otherwise, error('Unknown kernel function.');
            end
        else
            data = TYPE;
        end

        % Handle the case in which the dataset is empty
        if ~exist('A', 'var') || isempty(A)
            W = prmapping('lssvc', 'untrained', data);
            W = setname(W, name);
            return;
        end

        % Perform some checks on the data
        islabtype(A, 'crisp');          % allow crisp labels only
        isvaldset(A, 1, 2);             % at least one object per class, two objects
        if ischar(A.lablist), error('Only numerical labels are allowed in LSSVC.'); end

        % Train binary LS-SVM
        [m, k, c] = getsize(A);         % size of the training set
		if length(A.lablist) <= 2
            if strcmp(data.kernel, 'poly_kernel')
                [alpha, b] = trainlssvm({A.data, A.labels, 'classification', data.gamma, [data.par2; data.par], data.kernel});
            else
                [alpha, b] = trainlssvm({A.data, A.labels, 'classification', data.gamma, data.par, data.kernel});
            end

		% Train multi-class LS-SVM
		else
			[labels_code, codebook, old_codebook] = code(A.labels, 'code_MOC');
            if strcmp(data.kernel, 'poly_kernel')
                [alpha, b] = trainlssvm({A.data, labels_code, 'classification', data.gamma, [data.par2; data.par], data.kernel});
            else
                [alpha, b] = trainlssvm({A.data, labels_code, 'classification', data.gamma, data.par, data.kernel});
            end
		end

        % Store results of training in data-struct
        data.A = A;
        data.alpha = alpha;
        data.b = b;
		if length(A.lablist) > 2
			data.labels_code = labels_code;
			data.codebook = codebook;
			data.old_codebook = old_codebook;
		end
        W = prmapping('lssvc', 'trained', data, getlablist(A), k, c);
        W = setname(W, name);

    % Apply LS-SVM mapping B on dataset A
    elseif nargin == 2 && isa(TYPE, 'prmapping')

        % Initialize some variables
        B = TYPE;
        [m, k] = getsize(A);        % size of the test set
        [k, c] = size(B);           % K features with C classes
        W = zeros(m, c);            % output: C class densities for M objects

        % Perform binary classification
		if length(B.data.A.lablist) <= 2
            if strcmp(B.data.kernel, 'poly_kernel')
                C = simlssvm({B.data.A.data, B.data.A.labels, 'classification', B.data.gamma, [B.data.par2; B.data.par], B.data.kernel, 'preprocess'}, ...
                             {B.data.alpha, B.data.b}, A.data);
            else
				C = simlssvm({B.data.A.data, B.data.A.labels, 'classification', B.data.gamma, B.data.par, B.data.kernel, 'preprocess'}, ...
                             {B.data.alpha, B.data.b}, A.data);
            end

		% Perform multi-class classification (only numerical labels!)
        else
            if strcmp(B.data.kernel, 'poly_kernel')
                C = simlssvm({B.data.A.data, B.data.labels_code, 'classification', B.data.gamma, [B.data.par2; B.data.par], B.data.kernel, 'preprocess'}, ...
                             {B.data.alpha, B.data.b}, A.data);
            else
                C = simlssvm({B.data.A.data, B.data.labels_code, 'classification', B.data.gamma, B.data.par, B.data.kernel, 'preprocess'}, ...
                             {B.data.alpha, B.data.b}, A.data);
            end
			C = code(C, B.data.old_codebook, [], B.data.codebook, 'codedist_hamming');
        end

        % Make sure all labels are numbered in the same way as in A
        for i=1:length(C)
            [tf, loc] = ismember(B.data.A.lablist, C(i));
            W(i, find(loc)) = 1;
        end

        % Return new labeled dataset
        W = setdat(A, W, B);

    % Should not happen
    else
        error('Illegal call.');
    end

author: marianux

repo: ecg-kit

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


text:
clear all; close all; clc

addpath('apm')

s = 'http://byu.apmonitor.com';
b = 'mhe';

% Connect to Arduino
tclab;

% Run time in minutes
run_time = 10.0;

% Number of cycles (1 cycle per 3 seconds)
loops = round(20*run_time);

% milli-volts input
Q1 = zeros(1,loops);
Q2 = zeros(1,loops);
Q1(3:end) = 100.0;
Q1(50:end) = 0.0;
Q1(100:end) = 80.0;

Q2(25:end) = 60.0;
Q2(75:end) = 100.0;
Q2(125:end) = 25.0;

for i = 130:180
    if mod(i,10)==0
        Q1(i:i+10) = rand(1) * 100;
    end
    if mod(i+5,10)==0
        Q2(i:i+10) = rand(1) * 100;
    end
end

% Temperature (degC)
T1 = ones(1,loops) * T1C(); % measured T
T2 = ones(1,loops) * T2C(); % measured T
T1mhe = ones(1,loops) * T1C(); % measured T
T2mhe = ones(1,loops) * T2C(); % measured T
Umhe = ones(1,loops) * 10.0;
taumhe = ones(1,loops) * 5.0;
amhe1 = ones(1,loops) * 0.01;
amhe2 = ones(1,loops) * 0.0075;
time = zeros(1,loops);

% time
tm = zeros(1,loops);

% moving horizon estimation
mhe_init();

start_time = clock;
prev_time = start_time;

% dynamic plot (note: subplots needs to be declared here first)
figure(1)
subplot(3,1,1)
hold on, grid on
anexp1 = animatedline('LineStyle','-', 'Color', 'k', 'LineWidth', 2);
anpred1 = animatedline('LineStyle','--','Color', 'r','LineWidth', 2);
anexp2 = animatedline('LineStyle','-', 'Color', 'b', 'LineWidth', 2);
anpred2 = animatedline('LineStyle','--','Color', 'g','LineWidth', 2);
ylabel('Temperature \circC')
legend('T_1 Measured', 'T_1 Predicted', ...
    'T_2 Measured', 'T_2 Predicted', ...
    'Location', 'northwest')
title('Temperature Estimation')
subplot(3,1,2)
hold on, grid on
anQ1 = animatedline('LineStyle','-', 'Color', 'k', 'LineWidth', 2);
anQ2 = animatedline('LineStyle','--', 'Color', 'b', 'LineWidth', 2);
ylabel('Power Level Q (%)')
legend('Q_1', 'Q_2', 'Location', 'northwest')
subplot(3,1,3)
hold on, grid on
anU = animatedline('LineStyle','-', 'Color', 'r', 'LineWidth', 2);
antau = animatedline('LineStyle','--', 'Color', 'b', 'LineWidth', 2);
ana1 = animatedline('LineStyle','-', 'Color', 'g', 'LineWidth', 2);
ana2 = animatedline('LineStyle','-', 'Color', 'k', 'LineWidth', 2);
ylabel('Parameters')
legend('U', 'tau', 'a1 x 1000', 'a2 x 1000', 'Location', 'northwest')
xlabel('Time (sec)')

for ii = 1:loops
    % adjust power level
    h1(Q1(ii));
    h2(Q2(ii));

    % Pause Sleep time
    pause_max = 3.0;
    pause_time = pause_max - etime(clock,prev_time);
    if pause_time >= 0.0
        pause(pause_time - 0.01)
    else
        pause(0.01)
    end

    % Record time and change in time
    t = clock;
    dt = etime(t,prev_time);
    if ii>=2
        time(ii) = time(ii-1) + dt;
    end
    prev_time = t;

    % read and record from temperature controller
    T1(ii) = T1C();
    T2(ii) = T2C();

    % Start estimating parameters after 10 cycles (30 sec)
    if ii==10
       apm_option(s,b,'U.STATUS',1);
       apm_option(s,b,'tau.STATUS',1);
       apm_option(s,b,'a1.STATUS',1);
       apm_option(s,b,'a2.STATUS',1);
    end

    % non-linear energy balance
    jj = ii+1;
    params = mhe(T1(ii),T2(ii),Q1(ii),Q2(ii));
    Umhe(jj) = params(1);
    taumhe(jj) = params(2);
    a1(jj) = params(3);
    a2(jj) = params(4);
    T1mhe(jj) = params(5);
    T2mhe(jj) = params(6);

    % plot
    addpoints(anexp1,time(ii),T1(ii))
    addpoints(anpred1,time(ii),T1mhe(ii))
    addpoints(anexp2,time(ii),T2(ii))
    addpoints(anpred2,time(ii),T2mhe(ii))
    addpoints(anQ1,time(ii),Q1(ii))
    addpoints(anQ2,time(ii),Q2(ii))
    addpoints(anU,time(ii),Umhe(ii))
    addpoints(antau,time(ii),taumhe(ii))
    addpoints(ana1,time(ii),1000*a1(ii))
    addpoints(ana2,time(ii),1000*a2(ii))
    drawnow

    % open web-interface
    if ii==20
        apm_web(s,b);
    end
end

h1(0);
h2(0);
disp('Heaters off')
% turn off heater but keep LED on if T > 50
if (T1C() || T2C()) > 50
    led(1)
    disp(['Warning, heater temperature 1 =', num2str(T1C())])
    disp(['Warning, heater temperature 2 =', num2str(T2C())])
else
    led(0)
end

% save txt file with data
data = [time',Q1',Q2',T1',T2'];
csvwrite('data.txt',data);

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%          Inverse Short-Time Fourier Transform        %
%               with MATLAB Implementation             %
%                                                      %
% Author: M.Sc. Eng. Hristo Zhivomirov       12/26/13  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x, t] = istft(stft, h, nfft, fs)

% function: [x, t] = istft(stft, h, nfft, fs)
% stft - STFT matrix (only unique points, time across columns, freq across rows)
% h - hop size
% nfft - number of FFT points
% fs - sampling frequency, Hz
% x - signal in the time domain
% t - time vector, s

% estimate the length of the signal
coln = size(stft, 2);
xlen = nfft + (coln-1)*h;
x = zeros(1, xlen);

% form a periodic hamming window
win = hamming(nfft, 'periodic');

% perform IFFT and weighted-OLA
if rem(nfft, 2)                     % odd nfft excludes Nyquist point
    for b = 0:h:(h*(coln-1))
        % extract FFT points
        X = stft(:, 1 + b/h);
        X = [X; conj(X(end:-1:2))];

        % IFFT
        xprim = real(ifft(X));

        % weighted-OLA
        x((b+1):(b+nfft)) = x((b+1):(b+nfft)) + (xprim.*win)';
    end
else                                % even nfft includes Nyquist point
    for b = 0:h:(h*(coln-1))
        % extract FFT points
        X = stft(:, 1+b/h);
        X = [X; conj(X(end-1:-1:2))];

        % IFFT
        xprim = real(ifft(X));

        % weighted-OLA
        x((b+1):(b+nfft)) = x((b+1):(b+nfft)) + (xprim.*win)';
    end
end

W0 = sum(win.^2);                   % find W0
x = x.*h/W0;                        % scale the weighted-OLA

% calculate the time vector
actxlen = length(x);                % find actual length of the signal
t = (0:actxlen-1)/fs;               % generate time vector

end

author: dynamicslab


text:
function [results_table_pos, results_table_neg, r, excluded_region_table, r_excluded, region_list] = table_of_atlas_regions_covered(r, varargin)
% Make a table of which atlas parcels are covered by a set of regions or map identified in a study
%
% [results_table_pos, results_table_neg, r, excluded_region_table, r_excluded, region_list] = table_of_atlas_regions_covered(r, [atlas_obj])
%
% fMRI activation maps often include large suprathreshold areas that span
% multiple brain regions.  Traditional tables divide areas by contiguous
% regions ("blobs"), but this type of division is not useful if blobs cannot
% be easily summarized by a single anatomical label. A complementary approach
% is to identify which areas defined in a standard brain atlas ("parcels")
% are covered by the activation map. Parcels can be defined based on anatomy
% (e.g., cytoarchitecture) or function (e.g., resting-state fMRI).
% For example, a large blob may span the insula, claustrum, and putamen.
% This can be described in a table that lists all of these regions, along with
% the percentage of each parcel covered by the activation map.
%
% This function uses the object method subdivide_by_atlas() to create such
% a table, which is dividided into two tables with positive and negative average
% statistic values (results_table_pos, results_table_neg). These are Matlab
% table-class objects.  It also returns a region-class object for reference
% and rendering the subdivided blobs on brain slices or surfaces.
%
% The default threshold for coverage is 25% or more of the atlas parcel
% covered by the activation map.
%
% You can enter any atlas-class object as a 2nd argument. If you don't,
% a default atlas object will be used.
%
% Examples:
% -----------------------------------------------------------------
% % Load a thresholded map:
% % (Note: You must have Neuroimaging_Pattern_Masks repo with subfolders on your path)
% % This is a predictive map for drug craving (Koban et al. 2022, Nat Neurosci)
%
% ncsthr = fmri_data(which('NCS_multithr_001k5_005_05_pruned.nii'), 'noverbose');
%
% % Convert to region object:
% r = region(ncsthr);
%
% [results_table_pos, results_table_neg, r, excluded_region_table, table_legend_text] = table_of_atlas_regions_covered(r);
% montage(r, 'regioncenters', 'colormap');
%
% % Also see large blobs that didn't sufficiently cover any one atlas parcel:
% montage(r_excluded(cat(1, r_excluded.numVox) > 20), 'regioncenters', 'colormap');
%
% % You can use the image_vector method to apply this to a single-image
% % thresholded map:
% [results_table_pos, results_table_neg, r, excluded_region_table, r_excluded, region_labels] = table_of_atlas_regions_covered(ncsthr);
%
% Note: For a key to labels when using the Glasser 2016 Nature atlas, or the
% CANlab combined atlas (which uses Glasser), see GlasserTableS1.pdf in the
% original paper or CANlab github, or http://braininfo.rprc.washington.edu/
%
% see also region.table, for autolabeling of regions with an atlas and more
% detailed table output.
%
% Tor Wager, Feb 2023

n_cols = 140;                       % 140 good for HTML reports
sep_str = repmat('_', 1, n_cols);   % see textwrap

dosep = true;
dolegend = true;
table_legend_text = '';

[results_table_pos, results_table_neg] = deal(table());

if isempty(r) || length(r) == 0
    disp('No regions to display');
    return
end

if ~isempty(varargin)

    atlas_obj = varargin{1};
    if ~isa(atlas_obj, 'atlas'), error('2nd input must be an atlas-class object'), end
else

    atlas_obj = load_atlas('canlab2018_2mm');

end

[r, percent_coverage] = subdivide_by_atlas(r, atlas_obj);

% attach this so we can deal with reordering of regions
for i = 1:length(r)
    r(i).custom_info2 = percent_coverage(i);
    r(i).custom_info2_descrip = 'atlas parcel coverage';
end

% Separate subregions with positive and negative values if requested
% -------------------------------------------------------------------------
if dosep
    % separate pos and neg
    [poscl, negcl] = posneg_separate(r);

    r = [poscl negcl];
    ispos = [true(1, length(poscl)) false(1, length(negcl))]; % logical for splitting combined cl later

    clear poscl negcl

    fprintf('\n%s\nPositive Effects\n', sep_str)

else

    fprintf('\n%s\nTable of all regions\n', sep_str)
end

% Volume in mm^3
% count per cubic mm - voxel count * voxel volume
get_vox_volume_in_mm = @(r) size(r.XYZ, 2) .* prod(abs(diag(r.M(1:3, 1:3))));

Volume = zeros(length(r), 1);

for i = 1:length(r)

    Volume(i, 1) = get_vox_volume_in_mm(r(i));

end

results_table = table({r.shorttitle}', 'VariableNames', {'Region'});

results_table.Volume = Volume;

results_table.XYZ = round(cat(1, r.mm_center));

results_table = [results_table get_signed_max(r, 'Z', 'maxZ')];  % use function because may be empty, handle if so

results_table.Properties.VariableNames{4} = 'maxVal';

atlas_region_coverage = cat(1, r.custom_info2);

results_table = addvars(results_table, round(atlas_region_coverage), 'NewVariableNames', 'Coverage');

% Exclude
wh_exclude = cat(1, r.numVox) < 2 | atlas_region_coverage < 25;

excluded_region_table = results_table(wh_exclude, :);

results_table(wh_exclude, :) = [];

ispos(wh_exclude) = [];
r_excluded = r(wh_exclude);

r(wh_exclude) = [];

if dosep

    results_table_pos = results_table(ispos, :);
    results_table_neg = results_table(~ispos, :);

end

% Legend info

voxvol = prod(abs(diag(r(1).M(1:3, 1:3))));

table_legend_text = {'Note: XYZ values are Montreal Neurologic Institute coordinates in mm.'};
table_legend_text{2} = sprintf('Voxel volume is %3.2fmm^3', voxvol);
table_legend_text{3} = 'maxVal is maximum positive or negatively-valued value in the .Z field. Its meaning varies across analyses. For a t-test, it usually reflects the maximum t-value.';
table_legend_text{4} = 'Coverage: Percentage of voxels in the labeled atlas parcel covered by the identified activation region.';
table_legend_text{5} = 'Regions shown are those with >1 voxel and >=25% coverage of the atlas parcel. See excluded_region_table for additional regions excluded in this way.';

% Now split into positive and neg sub-tables and display
if dosep
    % Tables for pos and neg results separately

    if any(ispos)
        disp(results_table_pos)
    else
        disp('No regions to display');
    end

    fprintf('\nNegative Effects\n')
    if any(~ispos)
        disp(results_table_neg)
    else
        disp('No regions to display');
    end

else
    % No separation

    if size(results_table, 1) > 0
        disp(results_table)
    else
        disp('No regions to display');
    end

end % dosep

% Get region list
% -------------------------------------------------------------
region_list = cell(4, 1);
for i = 1:4, region_list{i} = ''; end
region_list{1} = 'Positive effects:';
region_list{3} = 'Negative effects:';

for i = 1:size(results_table_pos, 1)
    region_list{2} = [region_list{2} sprintf('%s (%2.0f%%), ', results_table_pos.Region{i}, results_table_pos.Coverage(i))];
end

for i = 1:size(results_table_neg, 1)
    region_list{4} = [region_list{4} sprintf('%s (%2.0f%%), ', results_table_neg.Region{i}, results_table_neg.Coverage(i))];
end

if dolegend == false || isempty(table_legend_text)

    return

else
    % Print legend text

    canlab_print_legend_text(table_legend_text{:});
end

disp('')
disp('List of regions covered:')
for i = 1:4
    disp(region_list{i});
end

% clean up text for legend
% if length(table_legend_text) > 2
%     table_legend_text(2:3) = [];
% end

end % main function

function val_table = get_signed_max(cl, myfield, tablevarname)
% Returns a table with var "MaxZ" nregions x 1, or empty if cl.Z is empty

%maxZ = @(i) max(double(cl(i).Z));

smax = @(i) signedmax(cl(i).(myfield));

for i = 1:length(cl)

    if ~isempty(cl(i).(myfield))
        myZ(i, 1) = smax(i);
    else
        myZ(i, 1) = NaN;
    end

    %if isinf(smax(i)), keyboard, end

end

% Fix infinite vals - only for .Z . so this is not generalizable beyond
% this function without modifications:
maxZ = norminv(1 - 1E-12);
myZ(myZ > maxZ) = maxZ;

if all(isnan(myZ))
    val_table = [];
else
    val_table = table(myZ,  'VariableNames', {tablevarname});
end

end % function

function val = signedmax(vals)

vals = double(vals);
[maxabs, wh] = max(abs(vals));

val = sign(vals(wh)) .* maxabs;

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function dep = cfg_vout_mkdir(job)

% Define virtual outputs for cfg_run_mkdir. The directory name can either
% be assigned to a cfg_files directory input or to a evaluated cfg_entry.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_vout_mkdir.m 1716 2008-05-23 08:18:45Z volkmar $

rev = '$Rev: 1716 $'; %#ok
dep            = cfg_dep;
if isa(job.name, 'cfg_dep')
    dep.sname  = sprintf('Make Directory');
else
    dep.sname  = sprintf('Make Directory ''%s''', job.name);
end;
dep.src_output = substruct('.','dir');
dep.tgt_spec   = cfg_findspec({{'filter','dir', 'strtype','e'}});

author: spm


text:
function C=circulantMatrix(v,reverseDir)
%%CIRCULANTMATRIX Create a circulant matrix where the first row is given by
%            v. A circulatn matrix is a matrix where each row going down
%            from the top rotates the elements of the previous row by 1.
%
%INPUTS: v A 1XN or NX1 vector containing the elemtns that go into the
%          first row.
% reverseDir If true, this boolean parameter indicates that each subsequent
%          row should be left-shifted from the previous one. If this
%          parameter is omitted or an empty matrix is passed, then the
%          default is false.
%
%OUTPUTS: C The NXN circulant matrix whose first row is v.
%
%EXAMPLE:
% CF=circulantMatrix(1:5,false)
% CR=circulantMatrix(1:5,true)
%The results are:
%CF=[1,2,3,4,5;
%    5,1,2,3,4;
%    4,5,1,2,3;
%    3,4,5,1,2;
%    2,3,4,5,1];
%CR=[1,2,3,4,5;
%    2,3,4,5,1;
%    3,4,5,1,2;
%    4,5,1,2,3;
%    5,1,2,3,4];
%
%October 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    if(nargin<2||isempty(reverseDir))
        reverseDir=false;
    end

    N=length(v);
    v=v(:);
    x=[v(N:-1:1,1);flipud(v(2:end))];
    if(reverseDir)
        idx=bsxfun(@plus,mod(0:-1:-((N-1)),N)',(N:-1:1));
    else
        idx=bsxfun(@plus,(0:(N-1))',(N:-1:1));
    end

    C=x(idx);
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 rendering = spm_cfg_render
% SPM Configuration file for Render
%__________________________________________________________________________
% Copyright (C) 2013-2016 Wellcome Trust Centre for Neuroimaging

% $Id: spm_cfg_render.m 6925 2016-11-09 17:23:40Z guillaume $

%==========================================================================
% Extract
%==========================================================================

%--------------------------------------------------------------------------
% data Data
%--------------------------------------------------------------------------
data         = cfg_files;
data.tag     = 'data';
data.name    = 'Data';
data.help    = {'Images to create rendering/surface from (usually grey and white matter segmentations).'};
data.filter  = 'image';
data.ufilter = '.*';
data.num     = [1 Inf];

%--------------------------------------------------------------------------
% mode Output
%--------------------------------------------------------------------------
mode        = cfg_menu;
mode.tag    = 'mode';
mode.name   = 'Output';
mode.help   = {'Operation mode.'};
mode.labels = {'Save Rendering'
               'Save Extracted Surface'
               'Save Rendering and Surface'}';
mode.values = {1 2 3};
mode.val    = {3};

%--------------------------------------------------------------------------
% thresh Surface isovalue(s)
%--------------------------------------------------------------------------
thresh         = cfg_entry;
thresh.tag     = 'thresh';
thresh.name    = 'Surface isovalue(s)';
thresh.help    = {...
    'Enter one or more values at which isosurfaces through the input images will be computed.'
    'This is only relevant for extracting surfaces, not rendering.'};
thresh.strtype = 'r';
thresh.val     = {0.5};
thresh.num     = [1 Inf];

%--------------------------------------------------------------------------
% extract Extract Surface
%--------------------------------------------------------------------------
extract      = cfg_exbranch;
extract.tag  = 'extract';
extract.name = 'Extract Surface';
extract.val  = {data mode thresh};
extract.help = {'Surface extraction.'};
extract.prog = @spm_surf;
extract.vout = @vout_extract;

%==========================================================================
% Render
%==========================================================================

%--------------------------------------------------------------------------
% spmmat Select SPM.mat
%--------------------------------------------------------------------------
spmmat         = cfg_files;
spmmat.tag     = 'spmmat';
spmmat.name    = 'Select SPM.mat';
spmmat.help    = {'Select the SPM.mat file that contains the design specification.'};
spmmat.filter  = 'mat';
spmmat.ufilter = '^SPM\.mat$';
spmmat.num     = [1 1];

%--------------------------------------------------------------------------
% contrasts Contrast(s)
%--------------------------------------------------------------------------
contrasts         = cfg_entry;
contrasts.tag     = 'contrasts';
contrasts.name    = 'Contrast(s)';
contrasts.help    = {'Index of contrast(s). If more than one number is entered, analyse a conjunction hypothesis.'};
contrasts.strtype = 'n';
contrasts.num     = [1 Inf];

%--------------------------------------------------------------------------
% threshdesc Threshold type
%--------------------------------------------------------------------------
threshdesc        = cfg_menu;
threshdesc.tag    = 'threshdesc';
threshdesc.name   = 'Threshold type';
threshdesc.help   = {''};
threshdesc.labels = {'FWE' 'none' 'FDR'};
threshdesc.values = {'FWE' 'none' 'FDR'};
threshdesc.val    = {'FWE'};

%--------------------------------------------------------------------------
% thresh Threshold
%--------------------------------------------------------------------------
thresh         = cfg_entry;
thresh.tag     = 'thresh';
thresh.name    = 'Threshold';
thresh.help    = {''};
thresh.strtype = 'r';
thresh.num     = [1 1];
thresh.val     = {0.05};

%--------------------------------------------------------------------------
% extent Extent (voxels)
%--------------------------------------------------------------------------
extent         = cfg_entry;
extent.tag     = 'extent';
extent.name    = 'Extent (voxels)';
extent.help    = {''};
extent.strtype = 'w';
extent.num     = [1 1];
extent.val     = {0};

%--------------------------------------------------------------------------
% contrasts Contrast(s)
%--------------------------------------------------------------------------
contrasts1         = cfg_entry;
contrasts1.tag     = 'contrasts';
contrasts1.name    = 'Contrast(s)';
contrasts1.help    = {'Index of contrast(s) for masking - leave empty for no masking.'};
contrasts1.strtype = 'n';
contrasts1.num     = [1 Inf];

%--------------------------------------------------------------------------
% thresh Mask threshold
%--------------------------------------------------------------------------
thresh1         = cfg_entry;
thresh1.tag     = 'thresh';
thresh1.name    = 'Mask threshold';
thresh1.help    = {''};
thresh1.strtype = 'r';
thresh1.num     = [1 1];
thresh1.val     = {0.05};

%--------------------------------------------------------------------------
% mtype Nature of mask
%--------------------------------------------------------------------------
mtype        = cfg_menu;
mtype.tag    = 'mtype';
mtype.name   = 'Nature of mask';
mtype.help   = {''};
mtype.labels = {'Inclusive' 'Exclusive'};
mtype.values = {0 1};

%--------------------------------------------------------------------------
% mask Mask definition
%--------------------------------------------------------------------------
mask      = cfg_branch;
mask.tag  = 'mask';
mask.name = 'Mask definition';
mask.val  = {contrasts1 thresh1 mtype};
mask.help = {''};

%--------------------------------------------------------------------------
% generic Masking
%--------------------------------------------------------------------------
generic1        = cfg_repeat;
generic1.tag    = 'generic';
generic1.name   = 'Masking';
generic1.help   = {''};
generic1.values = {mask};
generic1.num    = [0 1];

%--------------------------------------------------------------------------
% conspec Contrast query
%--------------------------------------------------------------------------
conspec      = cfg_branch;
conspec.tag  = 'conspec';
conspec.name = 'Contrast query';
conspec.val  = {spmmat contrasts threshdesc thresh extent generic1};
conspec.help = {''};

%--------------------------------------------------------------------------
% generic Contrasts
%--------------------------------------------------------------------------
generic        = cfg_repeat;
generic.tag    = 'generic';
generic.name   = 'Contrasts';
generic.help   = {''};
generic.values = {conspec};
generic.num    = [1 3];

%--------------------------------------------------------------------------
% rendfile Render File
%--------------------------------------------------------------------------
rendfile         = cfg_files;
rendfile.tag     = 'rendfile';
rendfile.name    = 'Render File';
rendfile.help    = {'File containing the images to render on to.'};
rendfile.filter  = {'mat','mesh'};
rendfile.ufilter = '.*';
rendfile.num     = [1 1];

%--------------------------------------------------------------------------
% render Display Surface
%--------------------------------------------------------------------------
render      = cfg_exbranch;
render.tag  = 'display';
render.name = 'Display Surface';
render.val  = {rendfile generic};
render.help = {'Surface rendering.'};
render.prog = @run_render;

%==========================================================================
% rendering Rendering
%==========================================================================
rendering        = cfg_choice;
rendering.tag    = 'render';
rendering.name   = 'Rendering';
rendering.help   = {'Rendering utilities.'};
rendering.values = {extract render};

%==========================================================================
function run_render(job)
for i=1:numel(job.conspec)
    xSPM.swd       = spm_file(job.conspec(i).spmmat{1},'fpath');
    xSPM.Ic        = job.conspec(i).contrasts;
    xSPM.u         = job.conspec(i).thresh;
    xSPM.Im        = [];
    if ~isempty(job.conspec(i).mask)
        xSPM.Im    = job.conspec(i).mask.contrasts;
        xSPM.pm    = job.conspec(i).mask.thresh;
        xSPM.Ex    = job.conspec(i).mask.mtype;
    end
    xSPM.thresDesc = job.conspec(i).threshdesc;
    xSPM.k         = job.conspec(i).extent;
    %xSPM.n        = 1; % conjunction
    xSPM.units     = {'mm' 'mm' 'mm'};
    [SPM, xSPM]    = spm_getSPM(xSPM);
    dat(i) = struct('XYZ', xSPM.XYZ,...
                    't',   xSPM.Z',...
                    'mat', xSPM.M,...
                    'dim', xSPM.DIM);
end
% Force non-interactive mode...
global prevrend
prevrend = struct('rendfile',job.rendfile{1}, 'brt',1, 'col',eye(3));
spm_render(dat,1,job.rendfile{1});

%==========================================================================
function dep = vout_extract(job)

cdep = 1;
if any(job.mode==[1 3])
    dep(cdep)            = cfg_dep;
    dep(cdep).sname      = 'Render .mat File';
    dep(cdep).src_output = substruct('.','rendfile');
    dep(cdep).tgt_spec   = cfg_findspec({{'filter','mat','strtype','e'}});
    cdep = cdep + 1;
end

if any(job.mode==[2 3])
    for k=1:numel(job.thresh)
        dep(cdep)            = cfg_dep;
        dep(cdep).sname      = sprintf('Surface .gii File (thr=%.02f)', ...
            job.thresh(k));
        dep(cdep).src_output = substruct('.','surffile', '()',{k});
        dep(cdep).tgt_spec   = cfg_findspec({{'filter','mesh','strtype','e'}});
        cdep = cdep + 1;
    end
end


text:
%kvclorec 'Sup-Reconstruction'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros vclorec.pane file
%
% Parameters:
% InputFile: i1 'Input Image #1', required: 'first input image'
% InputFile: i2 'Input Image #2', required: 'second input image'
% InputFile: str 'Struct Element', required: 'Structuring Element'
% OutputFile: o 'Output Image ', required: 'resulting output image'
%
% Example: o = kvclorec({i1, i2, str}, {'i1','';'i2','';'str','';'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% vclorec - Sup-Reconstruction
%
%  DESCRIPTION
% .I vclorec
% will do the closing by reconstruction of an input image by a set structuring element and a conditional image.
%
% The two input images must be of the same size and have the same
% number of data bands.  The input
% images may be of the same or different data types.  If the input
% images are of different data types, then all input images will be
% upcast to the highest input image data type.  The output
% image data type will be the same as that of the highest data type
% of the input images.
%
%
%
%  EXAMPLES
% vclorec -i1 ball.xv -i2 feath.xv -str b272.str -o outimage.xv
%
% Will do the closing by reconstruction of the image 1 "ball.xv", using "feath.xv" as the conditional image, by the structuring set b272.str, with the resulting image written to "outimage.xv".
%
%  "SEE ALSO"
% vcast1(1)
%
%  RESTRICTIONS
% .I vclorec
% can be defined for all data types supported by Khoros, but at the moment it has been implemented just for the bit and unsigned char types.
% The structuring elements are subsets of the 3x3 matrix and the origin is always at the center of this matrix.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993-1997 Junior Barrera, Roberto Lotufo.  All rights reserved.
%

function varargout = kvclorec(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kvclorec(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i1', '__input';'i2', '__input';'str', '__input';'o', '__output'};
maxval={0,0,0,0};
minval={0,0,0,0};
istoggle=[0,0,0,0];
was_set=istoggle * 0;
paramtype={'InputFile','InputFile','InputFile','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'vclorec"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function varargout = mvncdf(varargin)

switch class(varargin{1})

    case {'sdpvar','ndsdpvar'}

        if nargin > 1
            error('sdpvar/mvncdf currently only supports 1 argument, i.e. assumed zero mean and unit variance');
        end
        varargin{1} = reshape(varargin{1},[],1);
        varargout{1} = yalmip('define',mfilename,varargin{1});

    case 'char'

        X = varargin{3};
        F = [];

        operator = CreateBasicOperator('positive','increasing','callback');
        if length(X) == 1
            operator.convexity = @convexity;
            operator.derivative = @(x)(1/2)*(1/sqrt(2))*exp(-(x/sqrt(2)).^2)*2/sqrt(pi);
            operator.inverse = @(x)(sqrt(2)*erfinv(2*x-1));
        else
            operator.bounds = @bounds;
        end

        operator.range = [0 1];
        varargout{1} = F;
        varargout{2} = operator;
        varargout{3} = X;

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

function [L,U] = bounds(xL,xU)
L = mvncdf(xL);
U = mvncdf(xU);

function vexity = convexity(xL,xU)
if xL >= 0
    vexity = 'concave';
elseif xU <= 0
    vexity = 'convex';
else
    vexity = 'none';
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function gamma = omp(varargin)
%OMP Sparsity-constrained Orthogonal Matching Pursuit.
%  GAMMA = OMP(D,X,G,T) solves the optimization problem
%
%       min  |X - D*GAMMA|_2     s.t.  |GAMMA|_0 <= T
%      gamma
%
%  for each of the signals in X, using Batch Orthogonal Matching Pursuit.
%  Here, D is a dictionary with normalized columns, X is a matrix
%  containing column signals, T is the # of non-zeros in each signal
%  representation, and G is the Gramm matrix D'*D. The output GAMMA is a
%  matrix containing the sparse representations as its columns.
%
%  GAMMA = OMP(D,X,[],T) performs the same operation, but without the
%  matrix G, using OMP-Cholesky. This call produces the same output as
%  Batch-OMP, but is significantly slower. Using this syntax is only
%  recommended when available memory is too small to store G.
%
%  GAMMA = OMP(DtX,G,T) is the fastest implementation of OMP, but also
%  requires the most memory. Here, DtX stores the projections D'*X. In this
%  case Batch-OMP is used, but without having to compute D'*X in advance,
%  which slightly improves runtime. Note that in general, the call
%
%    GAMMA = OMP(D'*X,G,T);
%
%  will be faster than the call
%
%    GAMMA = OMP(D,X,G,T);
%
%  due to optimized matrix multiplications in Matlab. However, when the
%  entire matrix D'*X cannot be stored in memory, one of the other two
%  versions can be used. Both compute D'*X for just one signal at a time,
%  and thus require much less memory.
%
%  GAMMA = OMP(...,PARAM1,VAL1,PARAM2,VAL2,...) specifies additional
%  parameters for OMP. Available parameters are:
%
%    'gammamode' - Specifies the representation mode for GAMMA. Can be
%                  either 'full' or 'sparse', corresponding to a full or
%                  sparse matrix, respectively. By default, GAMMA is
%                  returned as a sparse matrix.
%    'messages'  - Specifies whether progress messages should be displayed.
%                  When positive, this is the number of seconds between
%                  status prints. When negative, indicates that no messages
%                  should be displayed (this is the default).
%    'checkdict' - Specifies whether dictionary normalization should be
%                  verified. When set to 'on' (default) the dictionary
%                  atoms are verified to be of unit L2-norm. Setting this
%                  parameter to 'off' disables verification and accelerates
%                  function performance. Note that an unnormalized
%                  dictionary will produce invalid results.
%    'profile'   - Can be either 'on' or 'off'. When 'on', profiling
%                  information is displayed at the end of the funciton
%                  execution.
%
%
%  Summary of OMP versions:
%
%    version      |   speed     |   memory
%  --------------------------------------------------
%   OMP(DtX,G,T)  |  very fast  |  very large
%   OMP(D,X,G,T)  |  fast       |  moderate
%   OMP(D,X,[],T) |  very slow  |  small
%  --------------------------------------------------
%
%
%  References:
%  [1] M. Elad, R. Rubinstein, and M. Zibulevsky, "Efficient Implementation
%      of the K-SVD Algorithm using Batch Orthogonal Matching Pursuit",
%      Technical Report - CS, Technion, April 2008.
%
%  See also OMP2.

%  Ron Rubinstein
%  Computer Science Department
%  Technion, Haifa 32000 Israel
%  ronrubin@cs
%
%  April 2009

% default options

sparse_gamma = 1;
msgdelta = -1;
checkdict = 1;
profile = 0;

% determine number of parameters

paramnum = 1;
while (paramnum<=nargin && ~ischar(varargin{paramnum}))
  paramnum = paramnum+1;
end
paramnum = paramnum-1;

% parse options

for i = paramnum+1:2:length(varargin)
  paramname = varargin{i};
  paramval = varargin{i+1};

  switch lower(paramname)

    case 'gammamode'
      if (strcmpi(paramval,'sparse'))
        sparse_gamma = 1;
      elseif (strcmpi(paramval,'full'))
        sparse_gamma = 0;
      else
        error('Invalid GAMMA mode');
      end

    case 'messages'
      msgdelta = paramval;

    case 'checkdict'
      if (strcmpi(paramval,'on'))
        checkdict = 1;
      elseif (strcmpi(paramval,'off'))
        checkdict = 0;
      else
        error('Invalid checkdict option');
      end

    case 'profile'
      if (strcmpi(paramval,'on'))
        profile = 1;
      elseif (strcmpi(paramval,'off'))
        profile = 0;
      else
        error('Invalid profile mode');
      end

    otherwise
      error(['Unknown option: ' paramname]);
  end

end

% determine call type

if (paramnum==3)
  DtX = varargin{1};
  G = varargin{2};
  T = varargin{3};
  D = [];
  X = [];
elseif (paramnum==4)
  D = varargin{1};
  X = varargin{2};
  G = varargin{3};
  T = varargin{4};
  DtX = [];
else
  error('Invalid number of parameters');
end

% verify dictionary normalization

if (checkdict)
  if (isempty(G))
    atomnorms = sum(D.*D);
  else
    atomnorms = diag(G);
  end
  if (any(abs(atomnorms-1) > 1e-2))
    error('Dictionary columns must be normalized to unit length');
  end
end

% omp

gamma = ompmex(D,X,DtX,G,T,sparse_gamma,msgdelta,profile);

author: thfylsty


text:
function [ w, x ] = line_rule ( a, b, order )

%*****************************************************************************80
%
%% LINE_RULE returns a quadrature rule for a line segment in 1D.
%
%  Discussion:
%
%    The integration region is:
%      A <= X <= B
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 September 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Carlos Felippa,
%    A compendium of FEM integration formulas for symbolic work,
%    Engineering Computation,
%    Volume 21, Number 8, 2004, pages 867-890.
%
%  Parameters:
%
%    Input, real A, B, the lower and upper limits.
%
%    Input, integer ORDER, the order of the rule.
%
%    Output, real W(ORDER), the weights.
%
%    Output, real X(ORDER), the abscissas.
%
  if ( order == 1 )
    [ w, x ] = line_unit_o01 ( );
  elseif ( order == 2 )
    [ w, x ] = line_unit_o02 ( );
  elseif ( order == 3 )
    [ w, x ] = line_unit_o03 ( );
  elseif ( order == 4 )
    [ w, x ] = line_unit_o04 ( );
  elseif ( order == 5 )
    [ w, x ] = line_unit_o05 ( );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LINE_RULE - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of ORDER.\n' );
    error ( 'LINE_RULE - Fatal error!' );
  end
%
%  Transform from [-1,+1] to [A,B]
%
  for j = 1 : order
    w(j) = w(j) * ( b - a ) / 2.0;
    x(j) = ( ( 1.0 - x(j) ) * a   ...
           + ( 1.0 + x(j) ) * b ) ...
           /   2.0;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function lagrange_nd_test09 ( )

%*****************************************************************************80
%
%% LAGRANGE_ND_TEST09 tests LAGRANGE_PARTIAL in 3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 January 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LAGRANGE_ND_TEST09\n' );
  fprintf ( 1, '  LAGRANGE_PARTIAL determines\n' );
  fprintf ( 1, '  the Lagrange interpolating polynomials L(x)\n' );
  fprintf ( 1, '  for ND points in D dimensions, assuming that\n' );
  fprintf ( 1, '  the number of points is less than or equal to\n' );
  fprintf ( 1, '  R = Pi(D,N), the number of monomials of degree N or less\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  For this example, the data points are the same as those\n' );
  fprintf ( 1, '  used by the level 2 Clenshaw Curtis sparse grid in 3D.\n' );

  sq2h = sqrt ( 2.0 ) / 2.0;

  d = 3;
  n = 4;
  r = mono_upto_enum ( d, n );
  nd = 25;
  xd = [ ...
      0.0,   0.0,   0.0; ...
     -1.0,   0.0,   0.0; ...
      1.0,   0.0,   0.0; ...
      0.0,  -1.0,   0.0; ...
      0.0,   1.0,   0.0; ...
      0.0,   0.0,  -1.0; ...
      0.0,   0.0,   1.0; ...
     -sq2h,  0.0,   0.0; ...
      sq2h,  0.0,   0.0; ...
     -1.0,  -1.0,   0.0; ...
      1.0,  -1.0,   0.0; ...
     -1.0,   1.0,   0.0; ...
      1.0,   1.0,   0.0; ...
      0.0,  -sq2h,  0.0; ...
      0.0,   sq2h,  0.0; ...
     -1.0,   0.0,  -1.0; ...
      1.0,   0.0,  -1.0; ...
     -1.0,   0.0,   1.0; ...
      1.0,   0.0,   1.0; ...
      0.0,  -1.0,  -1.0; ...
      0.0,   1.0,  -1.0; ...
      0.0,  -1.0,   1.0; ...
      0.0,   1.0,   1.0; ...
      0.0,   0.0,  -sq2h; ...
      0.0,   0.0,   sq2h ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension D = %d\n', d );
  fprintf ( 1, '  Maximum degree N = %d\n', n );
  fprintf ( 1, '  Number of monomials R = %d\n', r );
  fprintf ( 1, '  Number of data points ND = %d\n', nd );

  r8mat_transpose_print ( d, nd, xd, '  Data points XD:' );

  [ po, pc, pe ] = lagrange_partial ( d, n, r, nd, xd );
%
%  Print the polynomials.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Lagrange polynomials for XD data points:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : nd
    o = po(i);
    label = sprintf ( '  P(%d)(x) =', i );
    polynomial_print ( d, o, pc(i,1:o), pe(i,1:o), label );
  end
%
%  Evaluate the polynomials at XD.
%
  value = zeros ( nd, nd );

  for j = 1 : nd
    o = po(j);
    label = sprintf ( '  P(%d)(x) =', j );
    value(1:nd,j) = polynomial_value ( d, o, pc(j,1:o), pe(j,1:o), nd, xd );
  end

  err = r8mat_is_identity ( nd, value );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Frobenius norm of Lagrange matrix error = %g\n', err );

  return
end


text:
function gf = gradddivHexp(Hkt_r, varargin)
% function gf = gradddivH(Hkt_r, varargin)
% gf is 1 x k*t row vector
% Hkt_r = reshape(Hkt,1,k*t) -> row vector
% Vxt = varargin{1};  %data
% Wxk = varargin{2};  %W matrix
% Wxk_fix = varargin{3}; %fixed part of the Wxk matrix (e.g. background) ->rows
% Hkt_fix = varargin{4}; %fixed part (lines) of the H matrix (e.g. background)

alphaH=1; %for now....
Vxt = varargin{1};  %data
Wxk_tmp = varargin{2};  %W matrix
Wxk_fix = varargin{3}; %fixed part of the Wxk matrix (e.g. background) ->rows
Hkt_fix = varargin{4}; %fixed part (lines) of the H matrix (e.g. background)
peval = varargin{5}; %parameters

if ~isfield(peval, 'w_lambda') peval.w_lambda=0; end

t=size(Vxt,2);
k=length(peval.h_dovec);

Hkt_tmp = exp(reshape(Hkt_r,k,t));

Wxk = zeros(peval.numpix, peval.ncomp+1);
Hkt = zeros(peval.ncomp+1, peval.nt);

Wxk(:,peval.w_dovec)=Wxk_tmp;
Hkt(peval.h_dovec,:)=Hkt_tmp;

Wxk(:,peval.w_fixvec)=Wxk_fix;
Hkt(peval.h_fixvec,:)=Hkt_fix;

deltasum=sum(sum(Wxk_tmp,1))-k;
if and(~isempty(Wxk_tmp), abs(deltasum)>10^-6)
    error('Wxk is not correctly normalized! (sum(Wxk_tmp,1)<>1)\n sum(Wxk_tmp,1)=%f',deltasum)
end

ap = peval.alphapen;
sumH_t = sum(Hkt,2);
sumH_t_sq = sum(sumH_t.^ap);
sumH = sum(sumH_t);

penaltyterm = (1-ap)*(1/(sumH)^ap)*sum(sumH_t_sq) + ap*(1/sumH^(ap-1))*(sumH_t).^(ap-1);
penaltyterm_kt = repmat(penaltyterm, 1, t);

gfkt = (1-Wxk'*(Vxt./(Wxk*Hkt)))*alphaH.*Hkt + penaltyterm_kt; %d/dh(d-divergence)
% one is tehre because Wxt is normalized: sum(Wxt,1)=1
gf=reshape(gfkt(1:k,:),1,k*t); %making row vector


text:
classdef chebtech2 < chebtech
%CHEBTECH2   Approximate smooth functions on [-1,1] with Chebyshev interpolants.
%
%   Class for approximating smooth functions on the interval [-1,1] using
%   function values at 2nd-kind Chebyshev points and coefficients of the
%   corresponding 1st-kind Chebyshev series expansion.
%
% Constructor inputs:
%   CHEBTECH2(OP) constructs a CHEBTECH2 object from the function handle OP. OP
%   should be vectorised (i.e., accept a vector input) and ouput a vector of
%   the same length. CHEBTECH2 objects allow for array-valued construction
%   (i.e., of an array-valued function), in which case OP should accept a column
%   vector of length N and return a matrix of size NxM.
%
%   CHEBTECH2(OP, DATA) constructs a CHEBTECH2 using the additional data
%   supplied in the DATA structure.  Fields currently recognized are:
%     DATA.VSCALE    (Default:  0)
%     DATA.HSCALE    (Default:  1)
%         The constructor builds a CHEBTECH2 with 'happiness' (see
%         HAPPINESSCHECK.m) relative to the maximum of the given vertical scale
%         DATA.VSCALE and the (column-wise) infinity norm of the sampled
%         function values of OP, and the fixed horizontal scale DATA.HSCALE. If
%         not given (or given as empty), the VSCALE defaults to 0 initially,
%         and HSCALE defaults to 1.
%   If any fields in DATA are empty or not supplied, or if DATA itself is empty
%   or not supplied, appropriate default values are set.
%
%   CHEBTECH2(OP, DATA, PREF) overrides the default behavior with that given by
%   the preference structure PREF. See CHEBTECH.TECHPREF for details.
%
%   CHEBTECH2(VALUES, ...) returns a CHEBTECH2 object which interpolates the
%   values in the columns of VALUES at 2nd-kind Chebyshev points and
%   CHEBTECH2({VALUES, COEFFS}, ... ) uses the Chebyshev coefficients passed in
%   COEFFS rather than computing them from VALUES. If VALUES, is empty then the
%   CHEBTECH2 is constructed directly from the COEFFS. If COEFFS are passed,
%   the resulting CHEBTECH2 is always deemed 'happy'.
%
% Examples: % Basic construction: f = chebtech2(@(x) sin(x))
%
%   % Construction with preferences:
%   p.sampleTest = 0; % See CHEBTECH.TECHPREF for details
%   f = chebtech2(@(x) sin(x), [], [], p)
%
%   % Array-valued construction:
%   f = chebtech2(@(x) [sin(x), cos(x), exp(x)])
%
% See also CHEBTECH, CHEBTECH.TECHPREF, CHEBPTS, HAPPINESSCHECK, REFINE.

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CHEBTECH2 Class Description:
%
% The CHEBTECH2 class represents smooth functions on the interval [-1,1] using
% function values at 2nd-kind Chebyshev points and coefficients of the
% corresponding 1st-kind Chebyshev series expansion.
%
% The constructor is supplied with a handle that evaluates a given function on
% an increasingly fine Chebyshev 2nd-kind grid (see REFINE.m) until the
% representation is deemed 'happy' (see HAPPINESSCHECK.m). The resulting object
% can be used to evaluate and operate on the input function.
%
% More information can be found in the CHEBTECH class definition.
%
% Class diagram: [<<CHEBTECH>>] <-- [CHEBTECH2]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS CONSTRUCTOR:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = false )

        function obj = chebtech2(op, data, pref)
            % Parse inputs.
            if ( (nargin == 0) || isempty(op) )
                % Return an empty CHEBTECH2 on null input:
                return
            end

            if ( (nargin < 2) || isempty(data) )
                    data = struct();
            end

            if ( (nargin < 3) || isempty(pref) )
                pref = chebtech.techPref();
            else
                pref = chebtech.techPref(pref);
            end

            data = chebtech.parseDataInputs(data, pref);

            % Force nonadaptive construction if PREF.FIXEDLENGTH is numeric and
            % we're not using contour integrals.
            if ( ~(isnumeric(op) || iscell(op)) && ...
                    ~isnan(pref.fixedLength) && ~pref.useTurbo )
                % Evaluate op on the Chebyshev grid of given size:
                op = feval(op, chebtech2.chebpts(pref.fixedLength));
            end

            % Actual construction takes place here:
            [obj, values] = populate(obj, op, data, pref);

            if ( isnumeric(op) || iscell(op) )
                % Set length of obj to PREF.FIXEDLENGTH (if it is non-trivial).
                if ( ~isnan(pref.fixedLength ) )
                    obj = prolong(obj, pref.fixedLength);
                end

                % No need to error check when constructing from discrete data.
                return
            elseif ( obj.ishappy )
                % Use contour integrals if requested.
                if ( pref.useTurbo )
                    obj = constructorTurbo(obj, op, pref);
                end

                % No need to error check if we are happy:
                return
            end

            % Check for NaNs (if not happy):
            if ( pref.extrapolate )
                % Check for NaNs in interior only (because extrapolate was on):
                if ( any(any(isnan(obj.coeffs(2:end-1,:)))) )
                    error('CHEBFUN:CHEBTECH2:chebtech2:nanEval', ...
                        'Function returned NaN when evaluated.')
                end
                % We make sure not to return NaNs at +1 and -1.
                valuesTemp = extrapolate(obj, values);
                valuesTemp([1 end], :) = valuesTemp([1,end],:);
                obj.coeffs = obj.vals2coeffs(valuesTemp);
            elseif ( any(isnan(obj.coeffs(:))) )
                % Here we throw an error if NaNs were encountered anywhere.
                error('CHEBFUN:CHEBTECH2:chebtech2:nanEval2', ...
                    'Function returned NaN when evaluated.')
            end
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% STATIC METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = true )

        % Aliasing:
        coeffs = alias(coeffs, m)

        % Angles of Chebyshev points. (i.e., acos(chebpts(n))
        t = angles(n)

        % Evaluate a Chebyshev interpolant using the barycentric formula:
        out = bary(x, values)

        % Compute Chebyshev barycentric weights:
        w = barywts(n)

        % Compute Chebyshev points (x) and optionally quadrature (w)
        % and barycentric (v) weights:
        [x, w, v, t] = chebpts(n);

        % Tensor product grid of Chebyshev points in 1D, 2D or 3D:
        out = tensorGrid(N, dom)

        % Convert coefficients to values:
        values = coeffs2vals(coeffs);

        % Make a CHEBTECH2 (constructor shortcut):
        f = make(varargin);

        % Compute Chebyshev quadrature weights:
        w = quadwts(n)

        % Refinement function for CHEBTECH2 construction (evaluates OP on grid):
        [values, points, giveUp] = refine(op, values, pref)

        % Return the value-based discretization class which uses CHEBTECH2:
        function disc = returnValsDisc()
            disc = @chebcolloc2;
        end

        % Convert values to coefficients:
        coeffs = vals2coeffs(values)

    end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function padua_test04 ( )

%*****************************************************************************80
%
%% PADUA_TEST04 tests PADUA_POINTS and PADUA_POINTS_SET.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 June 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PADUA_TEST04\n' );
  fprintf ( 1, '  PADUA_POINTS computes the points of a Padua rule.\n' );
  fprintf ( 1, '  PADUA_POINTS_SET looks them up in a table.\n' );

  for l = 3 : 4
    n = padua_order ( l );
    xy1 = padua_points ( l );
    [ x2, y2 ] = padua_points_set ( l );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Level %d  Padua points\n', l );
    fprintf ( 1, '\n' );
    for j = 1 : n
      fprintf ( 1, '  %4d  %14.6g  %14.6g\n', j, xy1(1,j), xy1(2,j) );
      fprintf ( 1, '        %14.6g  %14.6g\n',    x2(j), y2(j) );
    end
  end

  return
end


text:
function [dLdp,iCpY,L] = mci_linsqr_deriv (P,M,U,Y)
% Gradient of likelihood for linear regression
% FORMAT [dLdp,iCpY,L] = mci_linsqr_deriv (P,M,U,Y)
%
% P         parameters
% M         model
% U         inputs
% Y         data
%
% dLdp      gradient of log joint
% iCpY      curvature (Fisher Information)
% L         log joint
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: mci_linsqr_deriv.m 6548 2015-09-11 12:39:47Z will $

G = mci_linsqr_gen (P,M,U);
if isstruct(Y)
    e = Y.y-G;
else
    e = Y-G;
end
X = U.X;

N=size(X,1);
dydp=2*(ones(N,1)*P').*X;

%dydp = spm_diff(M.IS,P,M,U,1);

dLdp=dydp'*M.iCe*e;
dLdp=dLdp';

iCpY=dydp'*M.iCe*dydp;

if nargout > 2
    L=mci_linsqr_like (P,M,U,Y);
end


text:
function ninv = norm_inv(x,mu,sigma)
%NORM_INV Inverse of the normal cumulative distribution function (cdf).
%   Y = NORM_INV(X,MU,SIGMA) Returns inverse of the normal cdf with
%   mean, MU, and standard deviation, SIGMA, at the values in X.
%
%   The size of X is the common size of the input arguments. A scalar input
%   functions as a constant matrix of the same size as the other inputs.
%
%   Default values for MU and SIGMA are 0 and 1 respectively.
%
% Copyright (c) 1995-1997, 2005-2007 Kurt Hornik

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

if nargin < 3
  sigma = 1;
end
if nargin < 2
  mu = 0;
end

if (~isscalar(mu) || ~isscalar(sigma))
  if ~((size(x,1) == size(mu,1)) && (size(mu,1) == size(sigma,1)))
    error ('norm_inv: x, mu and sigma must be of common size or scalars');
  end
end

[n,nin] = size (x);
ninv = zeros (n,nin);

if (isscalar (mu) && isscalar (sigma))
  if (find (isinf (mu) || isnan (mu) || (sigma < 0) || isinf(sigma)))
    ninv = NaN*ones(n,nin);
  else
    ninv =  mu+sigma.*sqrt(2)*erfinv(2*x-1);
  end
else
  k = find(isinf(mu) || isnan(mu) || (sigma < 0) || isinf(sigma));
  if (any(k))
    ninv(k) = NaN;
  end

  k = find(~isinf(mu) && ~isnan(mu) && (sigma > 0) && (sigma < Inf));
  if (any (k))
    ninv(k) = mu(k)+sigma(k).*sqrt(2)*erfinv(2*x(k)-1);
  end
end

k = find ((sigma == 0) && (x > 0) && (x < 1));
if (any(k))
  ninv(k) = mu(k);
end

ninv((sigma == 0) & (x == 0)) = -Inf;
ninv((sigma == 0) & (x == 1)) = Inf;

end

author: gpstuff-dev

repo: gpstuff

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


text:
% pop_loadascinstep() - import an INStep ASC (ASCII) file
%
% Usage:
%   >> OUTEEG = pop_loadascinstep( filename );
%
% Inputs:
%   filename       - file name
%
% Outputs:
%   OUTEEG         - EEGLAB data structure
%
% Author: Arnaud Delorme, SCCN, INC, UCSD, July 2006

% Copyright (C) 2006 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [EEG, com] = pop_loadascinstep( filename );

if nargin < 1
    [filename, filepath] = uigetfile('*.*', 'Choose an INStep ASC file -- pop_loadascinstep');
    drawnow;
    if filename == 0 return; end;
    filename = fullfile(filepath, filename);
end;

EEG = eeg_emptyset;
fid = fopen(filename, 'r');
if fid == -1, error('Cannot open file'); end;
num1 = fscanf(fid, '%f', 1);tmp   = fgetl(fid);
num2 = fscanf(fid, '%f', 1);tmp   = fgetl(fid);
num3 = fscanf(fid, '%f', 1);tmp   = fgetl(fid);

if ~isempty(num3)
    EEG.nbchan = num1;
    EEG.pnts   = num2;
    EEG.srate  = 1000/num3;
    tline = fgetl(fid);
    tline = fgetl(fid);
else
    EEG.nbchan = num1;
    EEG.srate  = 1000/num2;
    tline = tmp;
end;

allf  = parsetxt(tline);

EEG.data = fscanf(fid, '%f', [EEG.nbchan+5 Inf]);
fclose(fid);

EEG.xmin = EEG.data(3,1)/1000;
EEG.data(4,:) = EEG.data(4,:).*EEG.data(2,:); % stimulus times category
EEG.data([1:3],:) = []; % trial column
EEG.chanlocs = struct('labels', allf(4:end));
EEG = eeg_checkset(EEG);

EEG = pop_chanevent(EEG, 1, 'edge', 'leading', 'delchan', 'on');
EEG = pop_chanevent(EEG, 1, 'edge', 'leading', 'delchan', 'on', 'delevent', 'off', 'nbtype', 1, 'typename', 'resp' );
EEG = eeg_checkset(EEG, 'eventconsistency');

com = sprintf('EEG = pop_loadascinstep(''%s'');',filename);

author: PatternRecognition

repo: OpenBMI


text:
function [EEG, postOut] = prepPostProcess(EEG, postIn)
% PREP post-processing: HP filter, remove interpolated, cleanup metadata
%
%
%
%% Setup the options for EEGLab to use double precision
    [backupOptionsFile, currentOptionsFile, warningsState] = setupForEEGLAB();
    finishup = onCleanup(@() cleanup(backupOptionsFile, currentOptionsFile, ...
        warningsState));

%% Check the parameters
postOut = struct();

%% Initial setup
    try
        %% Check the input
        if nargin < 1 || ~isstruct(EEG)
            error('postProcess:NotEnoughArguments', 'first argument must be a structure');
        elseif nargin < 2 || ~exist('postIn', 'var') || isempty(postIn)
            postIn = struct();
        end
        if ~isstruct(postIn)
            error('postProcess:NoData', 'second argument must be a structure')
        end
        if ~isfield(EEG, 'etc') || ~isfield(EEG.etc, 'noiseDetection') || ...
                hasPrepErrors(EEG.etc.noiseDetection)
            error('postProcess:NoPrep', 'Prep has not executed or has errors');
        end

        %% Set up the output structure
        postOut = struct('keepFiltered', [], 'removeInterpolatedChannels', [], ...
            'cleanupReference', []);
        defaults = getPrepDefaults(EEG, 'postprocess');

        [postOut, errors] = checkPrepDefaults(postIn, postOut, defaults);
        if ~isempty(errors)
            error('postProcess:BadParameters', ['|' sprintf('%s|', errors{:})]);
        end
        defaults = getPrepDefaults(EEG, 'general');
        [postOut, errors] = checkPrepDefaults(postOut, postOut, defaults);
        if ~isempty(errors)
            error('postProcess:BadGeneralParameters', ['|' sprintf('%s|', errors{:})]);
        end
        EEG.etc.noiseDetection.postProcess = postOut;

    catch mex
        errorMessages.setup = ...
            ['postProcessing failed due to previous issues: ' ...
            getReport(mex, 'basic', 'hyperlinks', 'off')];
        EEG.etc.noiseDetection.errors.postProcess = errorMessages;
        return;
    end

%% Perform filtering if requested
     try
        if postOut.keepFiltered
          fprintf('Refiltering the data in post processing\n');
          [EEG, detrend] = removeTrend(EEG, postOut);
          EEG.etc.noiseDetection.postProcess.detrend = detrend;
        end
    catch mex
        errorMessages.detrend = ...
            ['postProcessing failed removeTrend: ' ...
            getReport(mex, 'basic', 'hyperlinks', 'off')];
        EEG.etc.noiseDetection.errors.status = 'post-processing error';
        EEG.etc.noiseDetection.errors.postProcess = errorMessages;
        if strcmpi(postOut.errorMsgs, 'verbose')
            warning('[%s]\n%s', mex.identifier, ...
                getReport(mex, 'extended', 'hyperlinks', 'on'));
        end
        return;
     end

%% Remove interpolated channels if requested
     try
        interpolatedChannels = EEG.etc.noiseDetection.reference.interpolatedChannels.all;
        if postOut.removeInterpolatedChannels && ~isempty(interpolatedChannels)
          EEG.etc.noiseDetection.postProcess.removedChannelNumbers = ...
               interpolatedChannels;
          EEG.etc.noiseDetection.postProcess.removedChannelData = ...
               EEG.data(interpolatedChannels, :);
          EEG.etc.noiseDetection.postProcess.removedChanlocs = ...
               EEG.chanlocs(interpolatedChannels);
          EEG.chanlocs(interpolatedChannels) = [];
          EEG.data(interpolatedChannels, :) = [];
          EEG.nbchan = length(EEG.chanlocs);
        else
           EEG.etc.noiseDetection.postProcess.removedChannelNumbers = [];
           EEG.etc.noiseDetection.postProcess.removedChannelLocs = [];
        end
    catch mex
        errorMessages.removeInterpolated = ...
            ['postProcessing failed to remove interpolated channels: ' ...
              getReport(mex, 'basic', 'hyperlinks', 'off')];
        EEG.etc.noiseDetection.errors.status = 'post-processing error';
        EEG.etc.noiseDetection.errors.postProcess = errorMessages;
        if strcmpi(postOut.errorMsgs, 'verbose')
            warning('[%s]\n%s', mex.identifier, ...
                getReport(mex, 'extended', 'hyperlinks', 'on'));
        end
        return;
     end
     postOut = EEG.etc.noiseDetection.postProcess;

 %% Cleanup if requested
     try
        if postOut.cleanupReference
           EEG = cleanupReference(EEG);
        end
    catch mex
        errorMessages.cleanupReference = ...
            ['postProcessing failed to cleanup reference information: ' ...
            getReport(mex, 'basic', 'hyperlinks', 'off')];
        EEG.etc.noiseDetection.errors.status = 'post-processing error';
        EEG.etc.noiseDetection.errors.postProcess = errorMessages;
        if strcmpi(postOut.errorMsgs, 'verbose')
            warning('[%s]\n%s', mex.identifier, ...
                getReport(mex, 'extended', 'hyperlinks', 'on'));
        end
        return;
     end

     %% If we got this far everything was good
     EEG.etc.noiseDetection.errors.postProcess = 0;
end

%% Cleanup callback
function cleanup(backupFile, currentFile, warningsState)
% Restore EEGLAB options file and warning settings
   restoreEEGOptions(backupFile, currentFile);
   warning(warningsState);
end % cleanup

author: VisLab

repo: EEG-Clean-Tools


text:
function ConvHDRtoStack(fmtIn, fmtOut, bSampling, ldr_gamma)
%
%        ConvHDRtoStack(fmtIn, fmtOut, bSampling, ldr_gamma)
%
%
%        For example:
%           ConvLDRtoLDR('hdr', 'jpg', 2.2);
%
%        This lines tonemaps all the .hdr files in the folder using the
%        Reinhard et al.'s operator and it saves them as .jpg files using
%        gamma 2.2
%
%        Input:
%           -fmtIn: an input string represeting the LDR format of the images
%           to be converted. This can be: 'hdr', 'pfm'
%           -fmtOut: an input string represeting the LDR format of
%           converted images. This can be: 'jpeg', 'jpg', 'png', etc.
%           -bSampling: using or not histogram sampling
%           -ldr_gamma: the encoding gamma for the LDR images. The default
%           value is 2.2
%
%        Output:
%           -ret: a boolean value, true or 1 if the method succeeds
%
%     Copyright (C) 2012-15  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

if(~exist('bSampling', 'var'))
    bSampling = 0;
end

if(~exist('ldr_gamma', 'var'))
    ldr_gamma = 2.2;
end

lst = dir(['*.', fmtIn]);

for i=1:length(lst)
    disp(lst(i).name);

    tmp_name = lst(i).name;
    img = hdrimread(tmp_name);
    L = lum(img);
    L = imresize(L, 0.5, 'bilinear');

    if(bSampling)
        fstops = ExposureHistogramSampling(L, 8, 2);
    else
        minL = round(log2(min(L(:))));
        maxL = round(log2(max(L(:))));
        fstops = -maxL:-minL;
    end

    for j=1:length(fstops)
        disp(fstops(j));
        img_exp_j = GammaTMO(img, ldr_gamma, fstops(j));

        bSkip = 0;
        if(~bSampling)
            tImg1 = ClampImg(round(255 * img_exp_j) / 255, 0, 1);
            val = mean(tImg1(:)); %mean value

            if((val < 0.1) || (val > 0.9))
                bSkip = 1;
            end

        end

        if(~bSkip)
            tmp_name_we = RemoveExt(tmp_name);
            tmp_name_out = [tmp_name_we, '_fstop_', num2str(j), '.', fmtOut];
            imwrite(img_exp_j, tmp_name_out);
        end
    end
end

end

author: banterle


text:
function next = next_boundary_node_q9 ( node )

%*****************************************************************************80
%
%% NEXT_BOUNDARY_NODE_Q9 returns the next boundary node in a Q9 element.
%
%  Reference Element Q9:
%
%    |
%    1  4--7--3
%    |  |     |
%    |  |     |
%    S  8  9  6
%    |  |     |
%    |  |     |
%    0  1--5--2
%    |
%    +--0--R--1-->
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE, the index of the current node.  An input
%    value of 0 (or any "unusual" value") indicates that the
%    first edge node is desired.
%
%    Output, integer NEXT, the index of the next edge node.
%
  if ( node == 1 )
    next = 5;
  elseif ( node == 5 )
    next = 2;
  elseif ( node == 2 )
    next = 6;
  elseif ( node == 6 )
    next = 3;
  elseif ( node == 3 )
    next = 7;
  elseif ( node == 7 )
    next = 4;
  elseif ( node == 4 )
    next = 8;
  else
    next = 1;
  end

  return
end


text:
%compute dytailcentral_mm

function [data,units]=compute_dytailcentral_mm(trx,n)

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
dytailcentral_mm=cell(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
    dytailcentral_mm{1,i}=(trx(larva).ytailcentral_mm(2:end)-trx(larva).ytailcentral_mm(1:end-1))./trx(larva).dt;
end

units=parseunits('mm/s');
data=dytailcentral_mm;

author: kristinbranson

repo: JAABA


text:
function r8_round2_test ( )

%*****************************************************************************80
%
%% R8_ROUND2_TEST tests R8_ROUND2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  x = pi;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8_ROUND2_TEST\n' );
  fprintf ( 1, '  R8_ROUND2 rounds a number to a\n' );
  fprintf ( 1, '  specified number of base 2 digits.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Test effect on PI:\n' );
  fprintf ( 1, '  X = %f\n', x );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  NPLACE  XROUND\n' );
  fprintf ( 1, '\n' );

  for i = 0 : 20
    nplace = i;
    xround = r8_round2 ( nplace, x );
    fprintf ( 1, '  %8d  %f\n', i, xround );
  end

  return
end


text:
function [y,t] = ramp(A,T,t0,c)
%
%RAMP unit ramp function
% function [y,t] = ramp(A,T,t0,c)
% This function generates the unit ramp function of
%amplitude A
%width 2T.
%shift of t0
%Scaling factor of c
%If no output is mentioned the function plots the signal.
%
%Example:
%ramp(1,2,3,4) plots the rampl signal of amplitude 1, stimewidth of 4
%time-shift of 3 and scaling factor of 4.
%[y,t]=ramp(1,2,3,4) is same as above with no plot.
%
% Author: B.H.Sri Hari, IITMadras
% Copyright B.H.Sri Hari.
%$Date:03/28/2006 13:04:15$
if nargin < 3,
    t0 = 0;
    c=1;
elseif nargin==3,
    c=1;
end
ra = abs(t0/c);
if ra < 1 | ra==0
    t = -2*T:.01:2*T;
elseif ra>=1 | ra==0
    t = -2*T*(ra):.01:2*T*(ra);
end
y = A*(c*t-t0).*(c*t>=t0);
y1 = zeros(1,length(t));
if nargout == 0,
    plot(t,y,'k','LineWidth',2);
    hold on;
    plot(t,y1,':k');
    hold on;
    plot(y1,t,':k')
    hold off
    title('Ramp signal')
   if A<0 & ra ~= 0
        axis([-2*T*(ra) 2*T*(ra) 4*A/3 -A/4]);
    elseif A<0  & ra == 0
        axis([-2*T 2*T 4*A/3 -A/4])
    elseif A>0 & ra~=0
        axis([-2*T*(ra) 2*T*(ra) -A/4 4*A/3]);
    else
            axis([-2*T 2*T -A/4 4*A/3]);
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function [ft2] = ha2ft2(ha)
% Convert area from hectares to square feet.
% Chad A. Greene 2012
ft2 = ha*107639.1041671;


text:
function [data, count] = freadbkj(infile, lines, bkformat, machinefmt ,r0, rN, c0, cN)
% FREADBK  --  Read binary data file.
%
%   DATA = FREADBK pops up a menu to select a filename, asks for the number
%   of lines and format of the input file, and reads in the data (double).
%
%   DATA = FREADBK(FILENAME) fopens the specified file with read permission and
%   returns a row vector DATA.  Format is assumed to be float32 format.
%
%   [DATA, COUNT] = FREADBK(FILENAME) optionally returns the number of elements
%   successfully read.
%
%   [DATA, COUNT] = FREADBK(FILENAME, NUMLINES) includes an optional NUMLINES
%   argument for the number of lines in the file.  Matrix DATA is returned of
%   size NUMLINES rows and the appropriate number of columns.
%
%   [DATA, COUNT] = FREADBK(FILENAME, NUMLINES, BKFORMAT) reads in a file in
%   the format BKFORMAT.  Which is one of the formats specified in FREAD (see
%   help fread), or one of these with a prepended 'cpx' for complex data types.
%   (DATA must be stored major row order, pixel interleaved):
%	'cpxfloat32'	complex floating point, 32 bits.
%	'cpxint16'	complex signed 16 bit short integers, ...
%	'cpx...'	...
%      ('mph'           format: use 'cpxfloat32')
%      ('hgt'           format: use 'freadhgt' function)
%
%   [DATA, COUNT] = FREADBK(FILENAME, NUMLINES, BKFORMAT, MACHINEFMT) additionally
%   specifies the order for reading or writing bytes or bits in the file using the machinefmt argument
%   MACHINEFMT:
%   'n' or 'native'         Your system byte ordering (default)
%   'b' or 'ieee-be'        Big-endian ordering (GAMMA software, for example)
%   'l' or 'ieee-le'        Little-endian ordering
%   's' or 'ieee-be.l64'    Big-endian ordering, 64-bit long data type
%   'a' or 'ieee-le.l64'    Little-endian ordering, 64-bit long data type
%
%   [DATA, COUNT] = FREADBK(FILENAME, NUMLINES, BKFORMAT, r0, rN, c0, cN) only
%   read in the part between: row r0 to row rN, and between column c0 to cN.
%   Indexes start at 1, stop at numlines/width. r0=0 all rows, c0=0 means all
%   columns.
%
%   Examples:
%     To read in a file with 2 channels (complex), stored pixel interleaved, ie,
%     (row1: RE1 - IM1 - RE2 - IM2 - RE3 - ...)
%     (row2: RE1 - IM1 - RE2 - IM2 - RE3 - ...)
%     (      ...                              )
%     in float32 format, which has a height of 100 lines, and a width of 250
%     complex pixels, use:
%       D = freadbk('filename',100,'cpxfloat32');
%
%     To crop this file while reading, between rows 1:10, columns 101:110, use:
%       D = freadbk('filename',100,'cpxfloat32','n',1, 10, 101, 110);
%
%   See also FWRITEBK, FOPEN, FREAD, FWRITE, LOAD, SAVE, FPRINTF,
%   FSEEK, FTELL, FREADHGT, FWRITEHGT
%

% $Revision: 1.12 $  $Date: 2001/09/28 14:24:32 $
% Bert Kampes, 04-Mar-2000

%   ======================================================================
%   04/2009 MJ add machinefmt
%   ======================================================================

%%% Handle input. (varargin?)
%if (nargin >  3) error('working on this... not yet ok, sorry.'); end;
false=0; true=1;
complextype=false; readwholefile=false;%		defaults

if (nargin <  8) cN = 0; end;%				fall through (not...)
if (nargin <  7) c0 = 0; end;%				fall through
if (nargin <  6) rN = 0; end;%				fall through
if (nargin <  5) r0 = 0; readwholefile=true; end;%	fall through
if (nargin <  4) machinefmt = 'n'; end;%	fall through
if (nargin <  3) bkformat='float32'; end;%		fall through
if (nargin <  2) lines=1; end;%				fall through
if (nargin <  1)
    [infile, inpath] = uigetfile('*', 'Select binary input float file', 0,0);
    infile   = [inpath,infile];
    lines    = input('Enter number of lines in file: ');
    bkformat = input('Format (enter between single quotes): ');% make a gui...
end;
%%% Check if whole file should be read regardless of r0,rNc0,cN.
if (c0==1 & r0==1 & rN==lines)
  bytesperelem = datatypesize(bkformat);%	in bytes
  filesize     = fsize(infile);%		in bytes
  filewidth    = filesize/(bytesperelem*lines);%ok for complex
  if (cN==filewidth)
    readwholefile=true;
  end;
end
% if (readwholefile==true)
%   disp(['% Reading whole file: ', infile]);
% end;

%%% Check bkformat for complex type: 'cpx*'
if (~ischar(bkformat)) error('FREADBK: bkformat must be string.'); end;
if (~ischar(infile))   error('FREADBK: infile must be string.'); end;
if (length(bkformat)==3)
  if (bkformat=='mph')
    disp('changing mph format to cpxfloat32');
    bkformat='cpxfloat32';
  end
  if (bkformat=='hgt')
    error('please use freadhgt for hgt format.');
  end
end
%%% complex types defined as prepended 'cpx'
if (length(bkformat)>6)
  if (bkformat(1:3)=='cpx')
    complextype = true;
    bkformat=bkformat(4:length(bkformat));
  end;
end;

%%% Read from file in column vector.
fid = fopen(infile,'r',machinefmt);
if (fid<0)%						try one more time
  [infile, inpath] = uigetfile('*', 'Select binary input file', 0,0);
  infile = [inpath,infile];
  fid    = fopen(infile,'r',machinefmt);
  if (fid<0) error(ferror(fid)); end;
end;

if (readwholefile==true)
  [data,count]=fread(fid,bkformat);%		read data in column vector
				%		count is number of elements, not bytes...
else
  bytesperelem = datatypesize(bkformat);%	in bytes
  filesize     = fsize(fid);%			in bytes
  filewidth    = filesize/(bytesperelem*lines);%ok for complex
  if (r0 > rN) error('r0 > rN'); end;
  if (c0 > cN) error('c0 > cN'); end;
  if (c0==0) c0=1; cN=filewidth; end;
  if (r0==0) r0=1; rN=lines; end;
  if (r0 < 1)            error('r0 < 1'); end;
  if (rN > lines)        error('rN > lines'); end;
  if (c0 < 1)            error('c0 < 1'); end;
  if (cN > filewidth)    error('cN > width'); end;
  if (~isint(filewidth)) error('numlines file seems to be wrong(?)'); end;

  data    = [];
  count   = 0;
  offset  = bytesperelem*(c0-1);%		in bytes
  mywidth = cN-c0+1;%				number of elems 2b read
  lines   = rN-r0+1;%				number of lines in matrix
  if (complextype==true)
    mywidth=mywidth*2;%				correction for pix interleaved complex
    offset =offset*2;%				correction for pix interleaved complex
  end;
  for ii=r0-1:rN-1
    start   = ii*filewidth*bytesperelem+offset;
    status  = fseek(fid,start,'bof');%			beginning of file
    if (status ~= 0)
      disp(['r0 rN c0 cN: ',num2str(r0),' ',num2str(rN),...
			' ',num2str(c0),' ',num2str(cN)]);
      disp(['ii offset start: ',...
	     num2str(ii),' ',num2str(offset),' ',num2str(start)]);
      error(ferror(fid));
    end;
    [d cnt] = fread(fid,mywidth,bkformat);%		column vector
    count   = count+cnt;
    data    = [data;d];%				augment column vector
  end;
end;
fclose(fid);

if (complextype==true)
  % should be possible by something like, but for speed makes no big diff.
  %data = reshape(data,width,2*lines);%		transpose band interleaved if complex
  %data = complex(data(:,1:lines),data(:,lines+1:2*lines)).';
  data = complex(data(1:2:count),data(2:2:count));
  count=count/2;%					correction for complex types
end;

width = count/lines;%
data  = reshape(data,width,lines).';

%%% Export to type if requested
% data = int8(data) ...

%%% EOF

author: DinhHoTongMinh

repo: TomoSAR


text:
% The following vectorized Matlab code implements Algorithm 1 from the
% paper
% Projection onto the probability simplex: An efficient algorithm with a
% simple proof, and an application, by W. Wang and M.A. Carreira-Perpinan,
% see https://arxiv.org/abs/1309.1541
%
% It projects each column vector in the D N matrix Y onto the probability
% simplex in D dimensions.

function X = SimplexColProj(Y)

Y = Y';
[N,D] = size(Y);
X = sort(Y,2,'descend');
Xtmp = (cumsum(X,2)-1)*diag(sparse(1./(1:D)));
X = max(bsxfun(@minus,Y,Xtmp(sub2ind([N,D],(1:N)',sum(X>Xtmp,2)))),0);
X = X';

author: hiroyuki-kasai

repo: NMFLibrary


text:
%DEMSOM1 Demonstrate SOM for visualisation.
%
%	Description
%	 This script demonstrates the use of a SOM with  a two-dimensional
%	grid to map onto data in  two-dimensional space.  Both on-line and
%	batch training algorithms are shown.
%
%	See also
%	SOM, SOMPAK, SOMTRAIN
%

%	Copyright (c) Ian T Nabney (1996-2001)

randn('state', 42);
rand('state', 42);
nin = 2;
ndata = 300;
% Give data an offset so that network has something to learn.
x = rand(ndata, nin) + ones(ndata, 1)*[1.5 1.5];

clc;
disp('This demonstration of the SOM, or Kohonen network, shows how the')
disp('network units after training lie in regions of high data density.')
disp('First we show the data, which is generated uniformly from a square.')
disp('Red crosses denote the data and black dots are the initial locations')
disp('of the SOM units.')
disp(' ')
disp('Press any key to continue.')
pause
net = som(nin, [8, 7]);
c1 = sompak(net);
h1 = figure;
plot(x(:, 1), x(:, 2), 'r+');
hold on
plot(c1(:,1), c1(:, 2), 'k.');
drawnow;  % Force figure to be drawn before training starts
options = foptions;

% Ordering phase
options(1) = 1;
options(14) = 50;
%options(14) = 5; % Just for testing
options(18) = 0.9;  % Initial learning rate
options(16) = 0.05; % Final learning rate
options(17) = 8;    % Initial neighbourhood size
options(15) = 1;    % Final neighbourhood size

disp('The SOM network is trained in two phases using an on-line algorithm.')
disp('Initially the neighbourhood is set to 8 and is then reduced')
disp('linearly to 1 over the first 50 iterations.')
disp('Each iteration consists of a pass through the complete')
disp('dataset, while the weights are adjusted after each pattern.')
disp('The learning rate is reduced linearly from 0.9 to 0.05.')
disp('This ordering phase puts the units in a rough grid shape.')
disp('Blue circles denote the units at the end of this phase.')
disp(' ')
disp('Press any key to continue.')
pause
net2 = somtrain(net, options, x);
c2 = sompak(net2);
plot(c2(:, 1), c2(:, 2), 'bo');
drawnow;

% Convergence phase
options(1) = 1;
options(14) = 400;
options(18) = 0.05;
options(16) = 0.01;
options(17) = 0;
options(15) = 0;

disp('The second, convergence, phase of learning just updates the winning node.')
disp('The learning rate is reduced from 0.05 to 0.01 over 400 iterations.')
disp('Note how the error value does not decrease monotonically; it is')
disp('difficult to decide when training is complete in a principled way.')
disp('The units are plotted as green stars.')
disp(' ')
disp('Press any key to continue.')
pause
net3 = somtrain(net2, options, x);
c3 = sompak(net3);
plot(c3(:, 1), c3(:, 2), 'g*');
drawnow;

% Now try batch training
options(1) = 1;
options(6) = 1;
options(14) = 50;
options(17) = 3;
options(15) = 0;
disp('An alternative approach to the on-line algorithm is a batch update')
disp('rule.  Each unit is updated to be the average weights')
disp('in a neighbourhood (which reduces from 3 to 0) over 50 iterations.');
disp('Note how the error is even more unstable at first, though eventually')
disp('it does converge.')
disp('The final units are shown as black triangles.')
disp(' ')
disp('Press any key to continue.')
pause
net4 = somtrain(net, options, x);
c4 = sompak(net4);
plot(c4(:, 1), c4(:, 2), 'k^')
legend('Data', 'Initial weights', 'Weights after ordering', ...
    'Weights after convergence', 'Batch weights', 2);
drawnow;

disp(' ')
disp('Press any key to end.')
disp(' ')
pause

close(h1);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt