text:
function line=plotLineROI(view, smoothStep, flip, titlename, plotFlag)
% plotLineROI - plot current map along a line ROI
%
% line=plotLineROI(view, smoothStep, flip, titlename, plotFlag);
%
% Arguments:
%  view : mrVista view struct
%  smoothStep : averaging steps across neighboring vertices [default=0]
%  flip       : flip line [default=false]
%  titlename  : title for plot
%
% Outputs:
%  line.x     : distance across cortical surface
%  line.y     : values in current view
%
% 2007/09 SOD: wrote it.
if ~exist('view','var') || isempty(view),
    view = getSelectedGray;
end
if ~exist('smoothStep','var') || isempty(smoothStep),
    smoothStep = 0;
else
    if smoothStep==-1, % aks subject
        ttltxt = sprintf('Enter number of smoothing steps [0 = no smoothing]');
        smoothStep = inputdlg('Number larger than 0, recommended range [0-5]:',ttltxt,1,{'0'});
        smoothStep = str2double(smoothStep{1});
        if isempty(smoothStep), smoothStep = 0; end
   end
end
if ~strcmp(view.viewType,'Gray'),
    error('This function is only defined for the Gray view');
end
if ~exist('flip','var') || isempty(flip),
    flip = false;
end
if ~exist('titlename','var') || isempty(titlename),
   titlename = [];
end
if ~exist('plotFlag','var') || isempty(plotFlag),
    plotFlag = 1;
end

mrGlobals;

% get ROI coords
roi = view.ROIs(view.selectedROI);

% get dimensions
mmPerPix = readVolAnatHeader(vANATOMYPATH);

% get all gray coords, nodes and edges
%coords  = double(viewGet(view,'coords'));
nodes   = double(view.nodes);
edges   = double(view.edges);
numNeighbors = double(view.nodes(4,:));
edgeOffsets  = double(view.nodes(5,:));

% Get nearest gray node
allNodeIndices = zeros(1,size(roi.coords,2));
for ii=1:size(roi.coords,2)
    grayNode = find(nodes(2,:) == roi.coords(1,ii) & ...
                    nodes(1,:) == roi.coords(2,ii) & ...
                    nodes(3,:) == roi.coords(3,ii));

    % Catch errors.
    if(isempty(grayNode))
        error('No gray nodes were found!');
    end
    if(length(grayNode)>1)
        fprintf('[%s]: WARNING- coord %d - more than one grayNode found!\n',mfilename,ii);
        grayNode = grayNode(1);
    end
    allNodeIndices(ii) = grayNode;
end

% HACK: find starting point:
% Define as the point that has the largest distance to
% one other point in the line with the fewest direct neighbors.
% This will in most cases be the end or the beginning.
allDist = zeros(numel(allNodeIndices));
for n=1:numel(allNodeIndices),
    tmp = mrManDist(nodes, edges, allNodeIndices(n), mmPerPix, -1, 0);
    allDist(:,n)    = tmp(allNodeIndices);
end

%figure;imagesc(allDist)

% nodes distances
distNodes   = zeros(numel(allNodeIndices),1);
for n=1:numel(distNodes)-1
    distNodes(n+1,1) = allDist(n,n+1);
end
distNodes =cumsum(abs(distNodes));

% get current map
curMap =  double(view.(view.ui.displayMode){view.curScan});

% smooth
if smoothStep,
    numNeighbors = numNeighbors(:);
    if ~strcmpi('ph',view.ui.displayMode),
        for n=1:smoothStep,
            curMap = sumOfNeighbors(curMap,edges,edgeOffsets,numNeighbors)./(numNeighbors+1);
        end
    else
        % phase data, so go complex, smooth and go back
        curMap = -exp(i*curMap).*double(view.co{view.curScan});
        cm_r   = real(curMap);
        cm_i   = imag(curMap);
        for n=1:smoothStep,
            cm_r = sumOfNeighbors(cm_r,edges,edgeOffsets,numNeighbors)./(numNeighbors+1);
            cm_i = sumOfNeighbors(cm_i,edges,edgeOffsets,numNeighbors)./(numNeighbors+1);
        end
        curMap = angle(cm_r+cm_i*i)+pi;
    end;
end

%%%% PLOT CODE BELOW

% line data
line.y   = curMap(allNodeIndices);
line.x   = distNodes;
line.y = line.y(:);
line.x = line.x(:);
if flip,
    line.x   = abs(line.x-max(line.x));
end

% ras 05/2008: if the user is plotting phase data, and if there are
% retinotopy params set (e.g. for mapping polar angle or eccentricity),
% auto-convert the phase values to the appropriate parameter
if strcmp(view.ui.displayMode, 'ph')
	rparams = retinoGetParams(view);
	if ~isempty(rparams)  % these parameters have been set
		if isequal(lower(rparams.type), 'polar_angle')
			% map from phase to polar angle
			line.y = polarAngle(line.y, rparams);
		else
			% assume eccentricity: map to eccentricity
			line.y = eccentricity(line.y, rparams);
		end

		% temporarily set the display mode to map, and the map name and units
		% to the relevant values (this should not propagate to the global view
		% variable)
		view.ui.displayMode = 'map';
		view.mapName = rparams.type;
		view.mapUnits = 'degrees';
	end
end

% Plot if requested
if plotFlag==1
    if strcmp(view.ui.displayMode,'ph')
        if (max(line.y)-min(line.y))>6
            % line data
            line.y   = pi-line.y;
            %line.coords = coords(:,sampleNodes);
            %line.index  = sampleNodes;
            % plot
            figure;
            hold on;
            plot(line.x,mod(line.y,2*pi)-pi,'ko-');
            plot([0 ceil(max(line.x)/10)*10],[0 0],'k');

            ylim([-pi pi])
        else
            line.y   = line.y-pi;
            % plot
            newGraphWin;
            hold on;
            plot(line.x,line.y,'ko-');
            plot([0 ceil(max(line.x)/10)*10],[0 0],'k');
            ylim([-pi pi])
        end
    else
        newGraphWin;
        hold on;
        plot(line.x,line.y,'ko-');
    end
    if ~isempty(titlename),
        title(titlename);
    end

    if isequal(view.ui.displayMode, 'map')
        ylabelText = sprintf('%s (%s)', view.mapName, view.mapUnits);
    else
        ylabelText = view.ui.displayMode;
    end
    ylabel(ylabelText, 'FontSize', 14, 'Interpreter', 'none');
    xlabel('distance (mm)', 'FontSize', 14, 'Interpreter', 'none');

    % kalanit wanted to add this -- don't blame me!
    set(gcf, 'Color', 'w');

    % add a button to flip the line along the X axis, in case the hack
    % didn't correctly identify the start and end points
    cb = ['tmp = get(gca, ''Children''); ' ...
          'set(tmp(end), ''XData'', fliplr(get(tmp(end), ''XData''))); ' ...
          'clear tmp '];
    uicontrol('Style', 'pushbutton', 'Units', 'norm', ...
              'Position', [.7 .02 .1 .05], 'String', 'Flip L/R', ...
              'Callback', cb);
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ a, info ] = dpofa ( a, lda, n )

%*****************************************************************************80
%
%% DPOFA factors a real symmetric positive definite matrix.
%
%  Discussion:
%
%    DPOFA is usually called by DPOCO, but it can be called
%    directly with a saving in time if RCOND is not needed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 June 2005
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Moler, Bunch and Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, real A(LDA,N), the symmetric matrix to be  factored.  Only the
%    diagonal and upper triangle are used.
%
%    Input, integer LDA, the leading dimension of the array A.
%
%    Input, integer N, the order of the matrix.
%
%    Output, real A(LDA,N), an upper triangular matrix R so that A = R'*R
%    where R' is the transpose.  The strict lower triangle is unaltered.
%    If INFO /= 0, the factorization is not complete.
%
%    Output, integer INFO, error flag.
%    0, for normal return.
%    K, signals an error condition.  The leading minor of order K is not
%    positive definite.
%
  for j = 1 : n

    s = 0.0;

    for k = 1 : j-1
      t = a(k,j) - ddot ( k-1, a(1:k-1,k), 1, a(1:k-1,j), 1 );
      t = t / a(k,k);
      a(k,j) = t;
      s = s + t * t;
    end

    s = a(j,j) - s;

    if ( s <= 0.0 )
      info = j;
      return
    end

    a(j,j) = sqrt ( s );

  end

  info = 0;

  return
end

author: johannesgerer

repo: jburkardt-m


function [env, versionString] = getEnvironment()
% Determine environment (Octave, MATLAB) and version string
% TODO: Unify private `getEnvironment` functions
    persistent cache

    if isempty(cache)
        isOctave = exist('OCTAVE_VERSION', 'builtin') ~= 0;
        if isOctave
            env = 'Octave';
            versionString = OCTAVE_VERSION;
        else
            env = 'MATLAB';
            vData = ver(env);
            versionString = vData.Version;
        end

        % store in cache
        cache.env = env;
        cache.versionString = versionString;

    else
        env = cache.env;
        versionString = cache.versionString;
    end
end


text:
function id=satno2id(sat)

global glc

[sys,prn]=satsys(sat);

switch sys
    case glc.SYS_GPS,id=sprintf('G%02d',prn-glc.MINPRNGPS+1);return;
    case glc.SYS_GLO,id=sprintf('R%02d',prn-glc.MINPRNGLO+1);return;
    case glc.SYS_GAL,id=sprintf('E%02d',prn-glc.MINPRNGAL+1);return;
    case glc.SYS_BDS,id=sprintf('C%02d',prn-glc.MINPRNBDS+1);return;
    case glc.SYS_QZS,id=sprintf('J%02d',prn-glc.MINPRNQZS+1);return;
end

id='';

return

repo: GINav


text:
function add_rx_tx_io(hRD)

% add AXI4 and AXI4-Lite slave interfaces
hRD.addAXI4SlaveInterface( ...
    'InterfaceConnection', 'axi_cpu_interconnect/M11_AXI', ...
    'BaseAddress',         '0x43C00000', ...
    'MasterAddressSpace',  'sys_ps7/Data');

% % AGC control input for transceiver
% hRD.addInternalIOInterface( ...
%     'InterfaceID',    'Enable AGC', ...
%     'InterfaceType',  'OUT', ...
%     'PortName',       'en_agc', ...
%     'PortWidth',      1, ...
%     'InterfaceConnection', 'gpio_en_agc', ...
%     'IsRequired',     false);

% GPIO status output for transceiver
hRD.addInternalIOInterface( ...
    'InterfaceID',    'CTRL_STATUS', ...
    'InterfaceType',  'IN', ...
    'PortName',       'gpio_status', ...
    'PortWidth',      8, ...
    'InterfaceConnection', 'gpio_status', ...
    'IsRequired',     false);

% GPIO Control input for transceiver
hRD.addInternalIOInterface( ...
    'InterfaceID',    'AD9361 CTRL IN', ...
    'InterfaceType',  'OUT', ...
    'PortName',       'gpio_ctl', ...
    'PortWidth',      4, ...
    'InterfaceConnection', 'gpio_ctl', ...
    'IsRequired',     false);

% DMA Ready signal
hRD.addInternalIOInterface( ...
    'InterfaceID',    'DMA Ready', ...
    'InterfaceType',  'IN', ...
    'PortName',       'dma_rdy', ...
    'PortWidth',      1, ...
    'InterfaceConnection', 'axi_ad9361_adc_dma/s_axis_ready', ...
    'IsRequired',     false);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Rx Reference design interfaces
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hRD.addInternalIOInterface( ...
    'InterfaceID',    'IP Data Valid OUT', ...
    'InterfaceType',  'OUT', ...
    'PortName',       'dut_data_valid', ...
    'PortWidth',      1, ...
    'InterfaceConnection', 'util_ad9361_adc_pack/adc_valid_0', ...
    'IsRequired',     false);

hRD.addInternalIOInterface( ...
    'InterfaceID',    'IP Data 0 OUT', ...
    'InterfaceType',  'OUT', ...
    'PortName',       'dut_data_0', ...
    'PortWidth',      16, ...
    'InterfaceConnection', 'util_ad9361_adc_pack/adc_data_0', ...
    'IsRequired',     false);

hRD.addInternalIOInterface( ...
    'InterfaceID',    'IP Data 1 OUT', ...
    'InterfaceType',  'OUT', ...
    'PortName',       'dut_data_1', ...
    'PortWidth',      16, ...
    'InterfaceConnection', 'util_ad9361_adc_pack/adc_data_1', ...
    'IsRequired',     false);

hRD.addInternalIOInterface( ...
    'InterfaceID',    'IP Data 2 OUT', ...
    'InterfaceType',  'OUT', ...
    'PortName',       'dut_data_2', ...
    'PortWidth',      16, ...
    'InterfaceConnection', 'util_ad9361_adc_pack/adc_data_2', ...
    'IsRequired',     false);

hRD.addInternalIOInterface( ...
    'InterfaceID',    'IP Data 3 OUT', ...
    'InterfaceType',  'OUT', ...
    'PortName',       'dut_data_3', ...
    'PortWidth',      16, ...
    'InterfaceConnection', 'util_ad9361_adc_pack/adc_data_3', ...
    'IsRequired',     false);

hRD.addInternalIOInterface( ...
    'InterfaceID',    'AD9361 ADC Data I0', ...
    'InterfaceType',  'IN', ...
    'PortName',       'sys_wfifo_0_dma_wdata', ...
    'PortWidth',      16, ...
    'InterfaceConnection', 'util_ad9361_adc_fifo/dout_data_0', ...
    'IsRequired',     false);

hRD.addInternalIOInterface( ...
    'InterfaceID',    'AD9361 ADC Data Q0', ...
    'InterfaceType',  'IN', ...
    'PortName',       'sys_wfifo_1_dma_wdata', ...
    'PortWidth',      16, ...
    'InterfaceConnection', 'util_ad9361_adc_fifo/dout_data_1', ...
    'IsRequired',     false);

hRD.addInternalIOInterface( ...
    'InterfaceID',    'AD9361 ADC Data I1', ...
    'InterfaceType',  'IN', ...
    'PortName',       'sys_wfifo_2_dma_wdata', ...
    'PortWidth',      16, ...
    'InterfaceConnection', 'util_ad9361_adc_fifo/dout_data_2', ...
    'IsRequired',     false);

hRD.addInternalIOInterface( ...
    'InterfaceID',    'AD9361 ADC Data Q1', ...
    'InterfaceType',  'IN', ...
    'PortName',       'sys_wfifo_3_dma_wdata', ...
    'PortWidth',      16, ...
    'InterfaceConnection', 'util_ad9361_adc_fifo/dout_data_3', ...
    'IsRequired',     false);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Tx Reference design interfaces
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hRD.addInternalIOInterface( ...
    'InterfaceID',    'AD9361 DAC Data I0', ...
    'InterfaceType',  'OUT', ...
    'PortName',       'axi_ad9361_dac_data_i0', ...
    'PortWidth',      16, ...
    'InterfaceConnection', 'axi_ad9361_dac_fifo/din_data_0', ...
    'IsRequired',     false);

hRD.addInternalIOInterface( ...
    'InterfaceID',    'AD9361 DAC Data Q0', ...
    'InterfaceType',  'OUT', ...
    'PortName',       'axi_ad9361_dac_data_q0', ...
    'PortWidth',      16, ...
    'InterfaceConnection', 'axi_ad9361_dac_fifo/din_data_1', ...
    'IsRequired',     false);

hRD.addInternalIOInterface( ...
    'InterfaceID',    'AD9361 DAC Data I1', ...
    'InterfaceType',  'OUT', ...
    'PortName',       'axi_ad9361_dac_data_i1', ...
    'PortWidth',      16, ...
    'InterfaceConnection', 'axi_ad9361_dac_fifo/din_data_2', ...
    'IsRequired',     false);

hRD.addInternalIOInterface( ...
    'InterfaceID',    'AD9361 DAC Data Q1', ...
    'InterfaceType',  'OUT', ...
    'PortName',       'axi_ad9361_dac_data_q1', ...
    'PortWidth',      16, ...
    'InterfaceConnection', 'axi_ad9361_dac_fifo/din_data_3', ...
    'IsRequired',     false);

hRD.addInternalIOInterface( ...
    'InterfaceID',    'IP Data 0 IN', ...
    'InterfaceType',  'IN', ...
    'PortName',       'util_dac_unpack_dac_data_00', ...
    'PortWidth',      16, ...
    'InterfaceConnection', 'util_ad9361_dac_upack/dac_data_0', ...
    'IsRequired',     false);

hRD.addInternalIOInterface( ...
    'InterfaceID',    'IP Data 1 IN', ...
    'InterfaceType',  'IN', ...
    'PortName',       'util_dac_unpack_dac_data_01', ...
    'PortWidth',      16, ...
    'InterfaceConnection', 'util_ad9361_dac_upack/dac_data_1', ...
    'IsRequired',     false);

hRD.addInternalIOInterface( ...
    'InterfaceID',    'IP Data 2 IN', ...
    'InterfaceType',  'IN', ...
    'PortName',       'util_dac_unpack_dac_data_02', ...
    'PortWidth',      16, ...
    'InterfaceConnection', 'util_ad9361_dac_upack/dac_data_2', ...
    'IsRequired',     false);

hRD.addInternalIOInterface( ...
    'InterfaceID',    'IP Data 3 IN', ...
    'InterfaceType',  'IN', ...
    'PortName',       'util_dac_unpack_dac_data_03', ...
    'PortWidth',      16, ...
    'InterfaceConnection', 'util_ad9361_dac_upack/dac_data_3', ...
    'IsRequired',     false);

hRD.addInternalIOInterface( ...
    'InterfaceID',    'IP Load Tx Data OUT', ...
    'InterfaceType',  'OUT', ...
    'PortName',       'util_dac_unpack_dac_valid_00', ...
    'PortWidth',      1, ...
    'InterfaceConnection', 'util_ad9361_dac_upack/dac_valid_0', ...
    'IsRequired',     false);

hRD.addInternalIOInterface( ...
    'InterfaceID',    'IP Valid Tx Data IN', ...
    'InterfaceType',  'IN', ...
    'PortName',       'util_dac_unpack_upack_valid_00', ...
    'PortWidth',      1, ...
    'InterfaceConnection', 'util_ad9361_dac_upack/upack_valid_0', ...
    'IsRequired',     false);

author: analogdevicesinc


text:
function cp = cornerprocess(c, T, q)
%CORNERPROCESS Processes the output of function cornermetric.
%	CP = CORNERPROCESS(C, T, Q) postprocesses C, the output of
%	function CORNERMETRIC, with the objective of reducing the
%	number of irrelevant corner points (with respect to threshold T)
%	and the number of multiple corners in a neighborhood of size
%	Q-by-Q. If there are multiple corner points contained within
%	that neighborhood, they are eroded morphologically to one corner
%	point.
%
%   A corner point is said to have been found at coordinates (I, J)
%   if C(I,J) > T.
%
%	A good practice is to normalize the values of C to the range [0
%	1], in im2double format before inputting C into this function.
%	This facilitates interpretation of the results and makes
%	thresholding more intuitive.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% Peform thresholding.
cp = c > T;

% Dilate CP to incorporate close neighbors.
B = ones(q);
cp = imdilate(cp, B);

% Shrink connnected components to single points.
cp = bwmorph(cp, 'shrink','Inf');

author: dipum

repo: dipum-toolbox

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


text:
function [s] = find_first_path(guesses,quiet)
  % [s] = find_first_path(guesses,quiet)
  si = find(cellfun(@(guess) exist(guess,'file'),guesses),1,'first');
  if isempty(si)
    if nargin>1 && quiet
      s = [];
      return;
    else
      error('Could not find path');
    end
  end
  s = guesses{si};
end

author: alecjacobson

repo: gptoolbox


text:
function model = update_integer_bounds(model);

if ~isempty(model.integer_variables)
    % Clean up things like 0.00001 <= x <= 0.999999 to 0 <= x <= 1
    % however, don't kill equalities 1 <= x <= 1
    lbfixed = fix(model.lb(model.integer_variables));
    ubfixed = fix(model.ub(model.integer_variables));
    fixlb = find(model.lb(model.integer_variables) == lbfixed);
    fixub = find(model.ub(model.integer_variables) == ubfixed);
    model.lb(model.integer_variables) = ceil(model.lb(model.integer_variables)-1e-4);
    model.ub(model.integer_variables) = floor(model.ub(model.integer_variables)+1e-4);
    if ~isempty(fixlb)
        model.lb(model.integer_variables(fixlb)) = lbfixed(fixlb);
    end
    if ~isempty(fixub)
        model.ub(model.integer_variables(fixub)) = ubfixed(fixub);
    end
end
if ~isempty(model.binary_variables)
    lbfixed = fix(model.lb(model.binary_variables));
    ubfixed = fix(model.ub(model.binary_variables));
    fixlb = find(model.lb(model.binary_variables) == lbfixed);
    fixub = find(model.ub(model.binary_variables) == ubfixed);
    model.lb(model.binary_variables) = ceil(model.lb(model.binary_variables)-1e-4);
    model.ub(model.binary_variables) = floor(model.ub(model.binary_variables)+1e-4);
    if ~isempty(fixlb)
        model.lb(model.binary_variables(fixlb)) = lbfixed(fixlb);
    end
    if ~isempty(fixub)
        model.ub(model.binary_variables(fixub)) = ubfixed(fixub);
    end
end
if any(model.lb(model.binary_variables) > model.ub(model.binary_variables))
    model.feasible = 0;
end
if any(model.lb(model.integer_variables) > model.ub(model.integer_variables))
    model.feasible = 0;
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function canlab_glm_publish(varargin)
% :Usage:
% ::
%
%     canlab_glm_publish(directory_specifications [options])
%
% :Directory Specification:
%
%   **'s', dirs:**
%        Generates HTML reports of the output from scn_spm_design_check for
%        directories in cell array dires (string allowable for single dir).
%        If a directory is a subject level analysis, the HTML will be generated
%        for that subject in the analysis directory.
%        Ex:
%        ::
%
%        canlab_glm_publish('s',{'model1/1011' 'model1/1014'})
%
%        If a directory contains subject level analyses, an HTML will be
%        generated with output for each subject level analysis.
%        Ex:
%        ::
%
%        canlab_glm_publish('s','model1')
%
%        ASSUMPTION: lower level analyses contain an SPM.mat file.
%
%
%   **'g', dirs:**
%        For each "robfit directory" in cell array dirs, will run robust_results_batch
%        on all contrast directories (e.g., robust0001/) (string allowable for single dir).
%        ("robfit directories" contain robfit contrast directories (like robust0001))
%        EITHER directories contain EXPT.mat files (see help robfit)
%        OR an EXPT struct is loaded in the workspace and a single directory is specified
%        OR will do best to figure out info normally contained in EXPT.mat
%        Ex:
%        ::
%
%        canlab_glm_publish('g', {'group_n35' 'group_anxiety_n35' 'group_sadness_n35'})
%
% :Note: directory paths may be absolute or relative (to working directory)
%
%
% :Options:
%
%   **'t', {[pthresh clustersize] ...}:**
%        Use the paired voxelwise_pthresh and minimum_cluster_size thresholds
%        with which to produce robfit results maps.
%        This option must follow immediately after a 'g' option (see above) and
%        will only apply to the analyses specified in that option.
%        ONLY applies to robfit directories (no bearing on lower level design checks)
%
%        DEFAULT: {[.001 5] [.005 1] [.05 1]}
%        Ex:
%        ::
%
%        canlab_glm_publish('g', pwd, 't', {[.001 1] [.005 10] [.05 10] [.01 25]})
%
%
%   **'email', address:**
%        send notification email to address when done running
%        Ex:
%        ::
%
%        canlab_glm_publish('g', pwd, 'email', 'ruzic@colorado.edu')
%

thresh = [.001 .005 .05];
size = [5 1 1];

%% setup
STARTINGDIR = pwd;

%% parse arguments

i=1;
while i<=numel(varargin)
    if ischar(varargin{i})
        switch varargin{i}
            case {'s' 'l'}
                i=i+1;
                if isstr(varargin{i})
                    varargin{i} = {varargin{i}};
                elseif ~iscellstr(varargin{i})
                    error('[Argument to ''s'' must a be string or cell array.]')
                end
                spmlower_report(STARTINGDIR,varargin{i});

            case {'g' 'r'}
                i=i+1;
                if isstr(varargin{i})
                    varargin{i} = {varargin{i}};
                elseif ~iscellstr(varargin{i})
                    error('[Argument to ''g'' must a be string or cell array.]')
                end
                robfitdirs = varargin{i};

                if i < nargin && strcmp(varargin{i+1},'t')
                    i=i+2;
                    t = cell2mat(varargin{i});
                    thresh = t(1:2:end);
                    size = t(2:2:end);
                end

                robfit_report(STARTINGDIR,robfitdirs,thresh,size);

            case {'email'}
                i=i+1;
                address = varargin{i};

            otherwise
                error(['Unrecognized argument: ' varargin{i}])

        end
    else
        error(['Unrecognized argument: ' varargin{i}])
    end
    i=i+1;
end

%% clean up
cd(STARTINGDIR)

%% email notification
if exist('address','var')
    try
        [ignore output] = system(sprintf('printf "canlab_glm_publish has finished running.\n" | mail -v -s "canlab_glm_publish done" %s',address)); %#ok
    catch  %#ok
        fprintf('The notification email failed to send.\n');
    end
end

end

function spmlower_report(STARTINGDIR,analysisdirs)

for i = 1:numel(analysisdirs)
    if regexp(analysisdirs{i},'^/')
        targdir = analysisdirs{i};
    else
        targdir = fullfile(STARTINGDIR, analysisdirs{i});
    end

    if ~exist(targdir,'dir')
        warning('No such directory: %s', targdir)
        continue
    end

    cd(targdir);

    outputdir = fullfile(targdir, 'design_review_html');
    mkdir(outputdir)

    p = struct('useNewFigure', false, 'maxHeight', 1500, 'maxWidth', 1200, ...
        'outputDir', outputdir, 'showCode', false);

    fout = publish('canlab_glm_publish_subject_levels.m', p);
    fprintf('Created subject level design review:\n\t%s\n', fout);
end

end

%% robfit group analyses (robust_results_batch)
function robfit_report(STARTINGDIR,robfitdirs,thresh,size)

assignin('base','thresh',thresh)
assignin('base','size',size)

for i = 1:numel(robfitdirs)
    if regexp(robfitdirs{i},'^/')
        targdir = robfitdirs{i};
    else
        targdir = fullfile(STARTINGDIR, robfitdirs{i});
    end

    if ~exist(targdir,'dir')
        fprintf('*** WARNING: skipping robfit directory (not found): %s\n', targdir)
        continue
    end

    cd(targdir);
    canlab_glm_publish_group_levels
end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function cnn_toy_data_generator(dataDir)
%CNN_TOY_DATA_GENERATOR
%   Generates toy data in the given path: random image of triangles,
%   squares and circles.
%
%   The directory format is: '<dataDir>/<set>/<label>/<sample>.png', where
%   <set> is 'train' or 'val', <label> is an integer between 1 and 3, and
%   <sample> is the sample index.

% Copyright (C) 2017 Joao F. Henriques.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

  % Set random seed
  rng('default') ;
  rng(0) ;

  % The sets, and number of samples per label in each set
  sets = {'train', 'val'} ;
  numSamples = [1500, 150] ;

  % Number of polygon points in each class. The circle is created with 50
  % points.
  numPoints = [3, 4, 50] ;

  for s = 1:2  % Iterate sets
    for label = 1:3  % Iterate labels
      fprintf('Generating images for set %s, label %i...\n', sets{s}, label) ;

      mkdir(sprintf('%s/%s/%i', dataDir, sets{s}, label)) ;

      for i = 1:numSamples(s)  % Iterate samples
        % Points of a regular polygon, with random rotation and scale
        radius = randi([11, 14]) ;
        angles = rand(1) * 2 * pi + (0 : 2 * pi / numPoints(label) : 2 * pi) ;
        xs = 16.5 + cos(angles) * radius ;
        ys = 16.5 + sin(angles) * radius ;

        % Generate image
        image = poly2mask(xs, ys, 32, 32) ;

        % Save it
        imwrite(image, sprintf('%s/%s/%i/%04i.png', dataDir, sets{s}, label, i)) ;
      end
    end
  end

end

repo: LapSRN


text:
% msTrain

LOAD = 1;

nsegments = [5 7 10 15 20 25 35 50 60 70 80 90 100 Inf]; % number of segments per segmentation
ncv = 1; % number of cross-validation sets
labeltol = 0.9; % required percentage of single-label pixels for segment to be good
nclasses = 23;

datadir = '~/data/eccv08/msrc21/';
imdir = '~/data/msrc/MSRC_ObjCategImageDatabase_v2/Images/';
outdir = '~/data/eccv08/msrc21/';

if LOAD
    load(fullfile(datadir, 'msrc_imsegs.mat'));
    load(fullfile(datadir, 'trainLabels.mat'));
    load(fullfile(datadir, 'trainTestFn_tu.mat'));
    load(fullfile(outdir, 'msSuperpixelLabels.mat'));
    load(fullfile(outdir, 'msSuperpixelData.mat'));
    load(fullfile(outdir, 'msEdgeData.mat'));
    load(fullfile(outdir, 'msEdgeClassifier_tu.mat'));
end

allind = union(train, test);
trainind1 = train(1:50);
trainind = train;
trainind2 = setdiff(trainind, trainind1);
testind = test;

nimages = numel(imsegs);

disp('Converting labels to superpixels')
if ~exist('splabels', 'var')
    splabels = msLabelMap2Sp(labels, {imsegs.segimage});
    save(fullfile(outdir, 'msSuperpixelLabels.mat'), 'splabels');
end

disp('Getting superpixel features')
if ~exist('spfeatures', 'var')
    spfeatures = mcmcGetAllSuperpixelData(imdir, imsegs);
    save(fullfile(outdir, 'msSuperpixelData.mat'), 'spfeatures');
end

disp('Getting edge features')
if ~exist('efeatures', 'var')
    [efeatures, adjlist] = mcmcGetAllEdgeData(spfeatures, imsegs);
    save(fullfile(outdir, 'msEdgeData.mat'), 'efeatures', 'adjlist');
end

disp('Training edge classifier')
if ~exist('eclassifier', 'var')
    eclassifier = msTrainEdgeClassifier(efeatures(trainind1), ...
        adjlist(trainind1), imsegs(trainind1), splabels(trainind1));
    ecal = msCalibrateEdgeClassifier(efeatures(trainind2), adjlist(trainind2), ...
        imsegs(trainind2), eclassifier, splabels(trainind2), ncv);
    ecal = ecal{1};
    save(fullfile(outdir, 'msEdgeClassifier_tu.mat'), 'eclassifier', 'ecal');
end

disp('Getting multiple segmentations')
if ~exist('smaps', 'var')
    for f = allind
        pE{f} = test_boosted_dt_mc(eclassifier, efeatures{f});
        pE{f} = 1 ./ (1+exp(ecal(1)*pE{f}+ecal(2)));
        smaps{f} = msCreateMultipleSegmentations(pE{f}, adjlist{f}, ...
            imsegs(f).nseg, nsegments);
    end
    save(fullfile(outdir, 'msMultipleSegmentations_tu_tmp.mat'), 'pE', 'smaps');
end

disp('Getting segment features')
if ~exist('segfeatures', 'var')
    for f = allind
        if mod(f, 50)==0, disp(num2str(f)), end
        im = im2double(imread([imdir '/' imsegs(f).imname]));
        imdata = mcmcComputeImageData(im, imsegs(f));

        for k = 1:numel(nsegments)
            segfeatures{f, k} = mcmcGetSegmentFeatures(imsegs(f), ...
                spfeatures{f}, imdata, smaps{f}(:, k), (1:max(smaps{f}(:, k))));
            [mclab{f, k}, mcprc{f, k}, allprc{f, k}, trainw{f, k}] = ...
                msSegs2labels(imsegs(f), smaps{f}(:, k), splabels{f}, nclasses);
            seglabel{f, k} = mclab{f, k}.*(mcprc{f, k}>labeltol);
            seggood{f,k} =  1*(mcprc{f, k}>labeltol) + (-1)*(mcprc{f, k}<labeltol);
        end
    end
    save(fullfile(outdir, 'msSegmentFeatures_tu_tmp.mat'), 'segfeatures', 'seglabel', 'seggood', 'trainw');
end

disp('Training segmentation classifier')
if ~exist('segclassifier', 'var')
    traindata = segfeatures(trainind, :);
    trainlabels = seggood(trainind, :);
    trainweights = trainw(trainind, :);
    segclassifier = mcmcTrainSegmentationClassifier2(traindata, trainlabels, trainweights);
    save(fullfile(outdir, 'msSegmentationClassifier_tu_tmp.mat'), 'segclassifier');
end

disp('Training label classifier')
if ~exist('labelclassifier', 'var')
    traindata = segfeatures(trainind, :);
    trainlabels = seglabel(trainind, :);
    trainweights = trainw(trainind, :);
    labelclassifier = msTrainLabelClassifier(traindata, trainlabels, trainweights, classnames, Inf);
    %labelclassifier = mcmcTrainSegmentClassifier2(traindata, trainlabels, trainweights);
    save(fullfile(outdir, 'msLabelClassifier_tu_tmp.mat'), 'labelclassifier');
end

save(fullfile(outdir, 'msClassifiers_tu_tmp.mat'), 'eclassifier', 'segclassifier', 'labelclassifier', 'ecal');

author: zouchuhang

repo: LayoutNet


text:
% STD_SELECTDESIGN - select an existing STUDY design.
%                      Use STD_MAKEDESIGN to add a new STUDY.design.
%
% Usage:
%   >> [STUDY] = std_selectdesign(STUDY, ALLEEG, designind);
%
% Inputs:
%   STUDY     - EEGLAB STUDY structure
%   STUDY     - EEGLAB ALLEEG structure
%   designind - desired (existing) design index
%
% Outputs:
%   STUDY     - EEGLAB STUDY structure with currentdesign set to the input design index.
%
% Author: Arnaud Delorme, Institute for Neural Computation, UCSD, 2010-

% Copyright (C) Arnaud Delorme, 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.

function STUDY = std_selectdesign(STUDY, ALLEEG, designind);

if nargin < 3
    help std_selectdesign;
    return;
end

if designind < 1 || designind > length(STUDY.design) || isempty(STUDY.design(designind).name)
    disp('Cannot select an empty STUDY.design');
    return;
end
STUDY.currentdesign = designind;
STUDY = std_rmalldatafields( STUDY );

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
% StackExchange Signal Processing Q60916
% https://dsp.stackexchange.com/questions/60916
% What Is the Bilateral Filter Category: LPF, HPF, BPF or BSF?
% References:
%   1.  A
% Remarks:
%   1.  B
% TODO:
% 	1.  C
% Release Notes
% - 1.0.001     02/03/2022
%   *   Using 'sgtitle()' instead of 'suptitle()' (Requires R2018b and above).
% - 1.0.000     19/08/2019
%   *   First release.

%% General Parameters

subStreamNumberDefault = 79;

run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = ON;

%% Simulation Parameters

inputImageFileName  = 'Lenna.png';

kernelRadius    = 7;
spatialStd      = 2;
vRangeStd       = [0.001, 0.01, 0.05, 0.5, 1e10];

%% Generate Data

mI = im2double(imread(inputImageFileName));
mI = mean(mI, 3);
mO = repmat(mI, [1, 1, 3]);

mRefPixels  = [84, 432; 82, 442; 267, 159; 420, 320];
mColors     = [1, 0, 0; 0, 1, 0; 0, 0, 1; 1, 0, 1];

%% Analysis

numPatches      = size(mRefPixels, 1);
kernelLength    = (2 * kernelRadius) + 1;

for ii = 1:numPatches
    mO(mRefPixels(ii, 1), mRefPixels(ii, 2), :) = reshape(mColors(ii, :), [1, 1, 3]);
    mO(mRefPixels(ii, 1) - kernelRadius - 1, mRefPixels(ii, 2) - kernelRadius - 1:mRefPixels(ii, 2) + kernelRadius + 1, :) = repmat(reshape(mColors(ii, :), [1, 1, 3]), [1, kernelLength + 2, 1]);
    mO(mRefPixels(ii, 1) + kernelRadius + 1, mRefPixels(ii, 2) - kernelRadius - 1:mRefPixels(ii, 2) + kernelRadius + 1, :) = repmat(reshape(mColors(ii, :), [1, 1, 3]), [1, kernelLength + 2, 1]);
    mO(mRefPixels(ii, 1) - kernelRadius - 1:mRefPixels(ii, 1) + kernelRadius + 1, mRefPixels(ii, 2) - kernelRadius - 1, :) = repmat(reshape(mColors(ii, :), [1, 1, 3]), [kernelLength + 2, 1, 1]);
    mO(mRefPixels(ii, 1) - kernelRadius - 1:mRefPixels(ii, 1) + kernelRadius + 1, mRefPixels(ii, 2) + kernelRadius + 1, :) = repmat(reshape(mColors(ii, :), [1, 1, 3]), [kernelLength + 2, 1, 1]);
end

figureIdx = figureIdx + 1;

hFigure     = figure('Position', figPosLarge);
hAxes       = axes;
hImageObj   = imshow(mO);
set(get(hAxes, 'Title'), 'String', {['Lenna Image with Patches']}, ...
    'FontSize', fontSizeTitle);

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

for ii = 1:length(vRangeStd)
    rangeStd = vRangeStd(ii);

    figureIdx = figureIdx + 1;
    hFigure = figure('Position', figPosLarge);

    plotIdx = 0;
    patchIdx = 0;

    for jj = 1:4
        patchIdx = patchIdx + 1;
        plotIdx = plotIdx + 1;

        hAxes       = subplot(4, 5, plotIdx);
        hImageObj   = imshow(mI(mRefPixels(patchIdx, 1) - kernelRadius:mRefPixels(patchIdx, 1) + kernelRadius, mRefPixels(patchIdx, 2) - kernelRadius:mRefPixels(patchIdx, 2) + kernelRadius));
        set(get(hAxes, 'Title'), 'String', {['Input Patch']}, ...
            'FontSize', fontSizeTitle);

        plotIdx = plotIdx + 1;

        % Spatial Weight
        mWs = CalcSpatialWeights(kernelRadius, spatialStd);

        hAxes       = subplot(4, 5, plotIdx);
        hImageObj   = imshow(mWs, []);
        set(get(hAxes, 'Title'), 'String', {['Spatial Weight']}, ...
            'FontSize', fontSizeTitle);

        plotIdx = plotIdx + 1;

        % Spatial Weight
        mWr = CalcRangeWeights(mI, mRefPixels(patchIdx, :), kernelRadius, rangeStd);

        hAxes       = subplot(4, 5, plotIdx);
        hImageObj   = imshow(mWr, []);
        set(get(hAxes, 'Title'), 'String', {['Range Weight']}, ...
            'FontSize', fontSizeTitle);

        plotIdx = plotIdx + 1;

        % Kernel Weight
        mW = mWr .* mWs;
        mW = mW / sum(mW(:));

        hAxes       = subplot(4, 5, plotIdx);
        hImageObj   = imshow(mW, []);
        set(get(hAxes, 'Title'), 'String', {['Kernel Weight']}, ...
            'FontSize', fontSizeTitle);

        plotIdx = plotIdx + 1;

        hAxes       = subplot(4, 5, plotIdx);
        hImageObj   = imshow(log(1 + abs(fftshift(fft2(mW)))), []);
        set(get(hAxes, 'Title'), 'String', {['Kernel - Frequency']}, ...
            'FontSize', fontSizeTitle);

    end

    sgtitle(hFigure, ['Kernel Analysis for rangeStd - ', num2str(rangeStd)]);

    if(generateFigures == ON)
        saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
    end

end

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
function [OBrientest] = OBrientest(X,alpha)
%O'Brien's Test for Homogeneity of Variances.
%[In the Obrien's test the data are transforming to
%yij = ((nj-1.5)*nj*((xij-mean(xj))**2)-((0.5)*(var(xj))*(nj-1)))/((nj-1)*(nj-2))
%and uses the F distribution performing an one-way ANOVA using y as the
%dependent variable (O'Brien, 1979)].
%
%   Syntax: function [OBrientest] = OBrientest(X,alpha)
%
%     Inputs:
%          X - data matrix (Size of matrix must be n-by-2; data=column 1, sample=column 2).
%       alpha - significance level (default = 0.05).
%     Outputs:
%          - Sample variances vector.
%          - Whether or not the homoscedasticity was met.
%
%    Example: From the example 10.1 of Zar (1999, p.180), to test the O'Brien's
%             homoscedasticity of data with a significance level = 0.05.
%
%                                 Diet
%                   ---------------------------------
%                       1       2       3       4
%                   ---------------------------------
%                     60.8    68.7   102.6    87.9
%                     57.0    67.7   102.1    84.2
%                     65.0    74.0   100.2    83.1
%                     58.6    66.3    96.5    85.7
%                     61.7    69.8            90.3
%                   ---------------------------------
%
%           Data matrix must be:
%            X=[60.8 1;57.0 1;65.0 1;58.6 1;61.7 1;68.7 2;67.7 2;74.0 2;66.3 2;69.8 2;
%            102.6 3;102.1 3;100.2 3;96.5 3;87.9 4;84.2 4;83.1 4;85.7 4;90.3 4];
%
%     Calling on Matlab the function:
%             OBrientest(X)
%
%       Answer is:
%
% The number of samples are: 4
%
% ----------------------------
% Sample    Size      Variance
% ----------------------------
%   1        5         9.3920
%   2        5         8.5650
%   3        4         7.6567
%   4        5         8.3880
% ----------------------------
%
% O'Brien's Test for Equality of Variances F=0.0171, df1= 3, df2=15
% Probability associated to the F statistic = 0.9968
% The associated probability for the F test is larger than 0.05
% So, the assumption of homoscedasticity was met.
%

%  Created by A. Trujillo-Ortiz and R. Hernandez-Walls
%             Facultad de Ciencias Marinas
%             Universidad Autonoma de Baja California
%             Apdo. Postal 453
%             Ensenada, Baja California
%             Mexico.
%             atrujo@uabc.mx
%
%  April 18, 2003.
%
%  To cite this file, this would be an appropriate format:
%  Trujillo-Ortiz, A. and R. Hernandez-Walls. (2003). Obrientest: O'Brien's test for
%    homogeneity of variances. A MATLAB file. [WWW document]. URL http://www.mathworks.com/
%    matlabcentral/fileexchange/loadFile.do?objectId=3335&objectType=FILE
%
%  References:
%
%  O'Brien, R. G. (1979), A General ANOVA Method for Robust Tests of
%           Additive Models for Variances. Journal of the American
%           Statistical Association, 74:877-880.
%  Zar, J. H. (1999), Biostatistical Analysis (2nd ed.).
%           NJ: Prentice-Hall, Englewood Cliffs. p. 180.
%

if nargin < 2,
   alpha = 0.05;
end

Y=X;
k=max(Y(:,2));
fprintf('The number of samples are:%2i\n\n', k);

%O'Brien Procedure.
n=[];s2=[];Z=[];
indice=Y(:,2);
for i=1:k
   Ye=find(indice==i);
   eval(['Y' num2str(i) '=Y(Ye,1);']);
   eval(['mY' num2str(i) '=mean(Y(Ye,1));']);
   eval(['n' num2str(i) '=length(Y' num2str(i) ') ;'])
   eval(['s2' num2str(i) '=(std(Y' num2str(i) ').^2) ;'])
   eval(['xn= n' num2str(i) ';'])
   eval(['xs2= s2' num2str(i) ';'])
   eval(['Z' num2str(i) '= ((n' num2str(i) ' - 1.5)*(n' num2str(i) ')*((Y' num2str(i) ' - mY' num2str(i) ').^2)-((0.5)*(s2' num2str(i) ')*(n' num2str(i) ' - 1)))/((n' num2str(i) ' - 1)*(n' num2str(i) ' - 2));']);
   eval(['x= Z' num2str(i) ';']);
   n=[n;xn];s2=[s2;xs2];Z=[Z;x];
end

for i=1:k
   if n(i)==2
      error('Requires sample sizes greater than two. Please, redefine the data matrix.');
   end
end

Y=[Z Y(:,2)];

fprintf('-----------------------------\n');
disp(' Sample    Size      Variance')
fprintf('-----------------------------\n');
for i=1:k
   fprintf('   %d       %2i         %.4f\n',i,n(i),s2(i))
end
fprintf('-----------------------------\n');
disp(' ')

C=(sum(Y(:,1)))^2/length(Y(:,1)); %correction term.
SST=sum(Y(:,1).^2)-C; %total sum of squares.
dfT=length(Y(:,1))-1; %total degrees of freedom.

indice=Y(:,2);
for i=1:k
   Ye=find(indice==i);
   eval(['A' num2str(i) '=Y(Ye,1);']);
end

A=[];
for i=1:k
   eval(['x =((sum(A' num2str(i) ').^2)/length(A' num2str(i) '));']);
   A=[A,x];
end

SSA=sum(A)-C; %sample sum of squares.
dfA=k-1; %sample degrees of freedom.
SSE=SST-SSA; %error sum of squares.
dfE=dfT-dfA; %error degrees of freedom.
MSA=SSA/dfA; %sample mean squares.
MSE=SSE/dfE; %error mean squares.
F=MSA/MSE; %F-statistic.
v1=dfA;df1=v1;
v2=dfE;df2=v2;

P = 1 - fcdf(F,v1,v2);  %probability associated to the F-statistic.

fprintf('O''Brien''s Test for Equality of Variances F=%3.4f, df1=%2i, df2=%2i\n', F,df1,df2);
fprintf('Probability associated to the F statistic = %3.4f\n', P);

if P >= alpha;
  fprintf('The associated probability for the F test is equal or larger than% 3.2f\n', alpha);
  fprintf('So, the assumption of homoscedasticity was met.\n');
else
  fprintf('The associated probability for the F test is smaller than% 3.2f\n', alpha);
  fprintf('So, the assumption of homoscedasticity was not met.\n');
end

author: Sable

repo: mcbench-benchmarks


text:
function can = mpot_to_cpot(mom)
% MPOT_TO_CPOT Convert a moment potential to canonical form.
% mom = mpot_to_cpot(can)

[g, h, K] = moment_to_canonical(mom.logp, mom.mu, mom.Sigma);
can = cpot(mom.domain, mom.sizes, g, h, K);

%%%%%%%%%%%

function [g, h, K] = moment_to_canonical(logp, mu, Sigma)
% MOMENT_TO_CANONICAL Convert moment characteristics to canonical form.
% [g, h, K] = moment_to_canonical(logp, mu, Sigma)

K = inv(Sigma);
h = K*mu;
n = length(K);
if isempty(mu)
  g = logp + 0.5*(log(det(K)) - n*log(2*pi));
else
  g = logp + 0.5*(log(det(K)) - n*log(2*pi) - mu'*K*mu);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [l] = read_label(sname, lname)
% l = read_label(<sname>, lname)
%
% reads the label file 'lname' from the subject 'sname'
% in the subject's label directory into the vector l
% l will be nvertices-by-5, where each column means:
% (1) vertex number, (2-4) xyz at each vertex, (5) stat
%
% IMPORTANT: the vertex number is 0-based.
%

%
% read_label.m
%
% Original Author: Bruce Fischl
% CVS Revision Info:
%    $Author: nicks $
%    $Date: 2011/03/02 00:04:12 $
%    $Revision: 1.7 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

l = [];

if(nargin ~= 2)
  fprintf('l = read_label(<sname>, lname)\n');
  return;
end

if(~isempty(sname))
  sdir = getenv('SUBJECTS_DIR') ;
  fname = sprintf('%s/%s/label/%s.label', sdir, sname, lname) ;
else
	ind = findstr(lname, '.label') ;
	if (length(ind) > 0)
        fname = lname ;
    else
		fname = sprintf('%s.label', lname);
	end
end

% open it as an ascii file
fid = fopen(fname, 'r') ;
if(fid == -1)
  fprintf('ERROR: could not open %s\n',fname);
  return;
end

fgets(fid) ;
if(fid == -1)
  fprintf('ERROR: could not open %s\n',fname);
  return;
end

line = fgets(fid) ;
nv = sscanf(line, '%d') ;
l = fscanf(fid, '%d %f %f %f %f\n') ;
l = reshape(l, 5, nv) ;
l = l' ;

fclose(fid) ;


text:
function out = isDiag(~)
%ISDIAG   Returns 1 for a SPINOP.
%
% See also SPINOP.

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

out = 1;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [in] = mm2in(mm)
% Convert length from millimeters to inches.
% Chad A. Greene 2012
in = mm/25.4;


text:
function [u, Ps, ue] = spm_uc_clusterFDR(q,df,STAT,R,n,Z,XYZ,V2R,ui,G)
% Cluster False Discovery critical extent threshold
% FORMAT [u, Ps, ue] = spm_uc_clusterFDR(q,df,STAT,R,n,Z,XYZ,ui[,G])
%
% q     - prespecified upper bound on False Discovery Rate
% df    - [df{interest} df{residuals}]
% STAT  - statistical field
%         'Z' - Gaussian field
%         'T' - T - field
%         'X' - Chi squared field
%         'F' - F - field
% R     - RESEL Count {defining search volume}
% n     - conjunction number
% Z     - height {minimum over n values}
%         or mapped statistic image(s)
% XYZ   - locations [x y x]' {in voxels}
%         or vector of indices of elements within mask
%         or mapped mask image
% V2R   - voxel to resel
% ui    - feature-inducing threshold
% G     - patch structure (for surface-based inference)
%
% u     - critical extent threshold
% Ps    - sorted p-values
% ue    - critical extent threshold for FWE
%__________________________________________________________________________
%
% References
%
% J.R. Chumbley and K.J. Friston, "False discovery rate revisited: FDR and
% topological inference using Gaussian random fields". NeuroImage,
% 44(1):62-70, 2009.
%
% J.R. Chumbley, K.J. Worsley, G. Flandin and K.J. Friston, "Topological
% FDR for NeuroImaging". NeuroImage, 49(4):3057-3064, 2010.
%__________________________________________________________________________
% Copyright (C) 2009-2013 Wellcome Trust Centre for Neuroimaging

% Justin Chumbley & Guillaume Flandin
% $Id: spm_uc_clusterFDR.m 5809 2013-12-20 14:30:22Z guillaume $

% Read statistical value from disk if needed
%--------------------------------------------------------------------------
if isstruct(Z)
    Vs         = Z;
    Vm         = XYZ;
    [Z, XYZmm] = spm_read_vols(Vs(1),true);
    for i=2:numel(Vs)
        Z      = min(Z, spm_read_vols(Vs(i),true));
    end
    Z          = Z(:)';
    XYZ        = round(Vs(1).mat \ [XYZmm; ones(1, size(XYZmm, 2))]);
    try
        M      = spm_vol(spm_file(Vs(1).fname,'basename','mask'));
        I      = spm_get_data(M,XYZ,false) > 0;
    catch
        I      = ~isnan(Z) & Z~=0;
    end
    XYZ        = XYZ(1:3,I);
    Z          = Z(I);
    if isstruct(Vm)
        Vm     = spm_get_data(Vm,XYZ,false) > 0;
    end
    XYZ        = XYZ(:,Vm);
    Z          = Z(:,Vm);
end

% Threshold the statistical field
%--------------------------------------------------------------------------
XYZ      = XYZ(:,Z >= ui);

% Extract size of excursion sets
%--------------------------------------------------------------------------
if nargin == 9
    N    = spm_clusters(XYZ);
else
    tmp  = false(size(G.vertices,1),1);
    tmp(XYZ(1,:)) = true;
    N    = spm_mesh_clusters(G,tmp);
    N    = N(XYZ(1,:));
end
if ~isempty(N)
    N    = histc(N,(0:max(N))+0.5);
end
N        = N(1:end-1);
N        = N .* V2R;

% Compute uncorrected p-values based on N using Random Field Theory
%--------------------------------------------------------------------------
Ps       = zeros(1,numel(N));
Pk       = zeros(1,numel(N));
ws       = warning('off','SPM:outOfRangePoisson');
for i = 1:length(N)
    [Pk(i), Ps(i)] = spm_P_RF(1,N(i),ui,df,STAT,R,n);
end
warning(ws);
[Ps, J]  = sort(Ps, 'ascend');

S        = length(Ps);

% Calculate FDR inequality RHS
%--------------------------------------------------------------------------
cV       = 1;    % Benjamini & Yeuketeli cV for independence/PosRegDep case
Fi       = (1:S)/S*q/cV;

% Find threshold
%--------------------------------------------------------------------------
I        = find(Ps <= Fi, 1, 'last');
if isempty(I)
    u    = Inf;
else
    u    = N(J(I)) / V2R;
end

% As we're there, also determine the FWE critical extent threshold
%--------------------------------------------------------------------------
if nargout == 3
    [Pk, J] = sort(Pk, 'ascend');
    I       = find(Pk <= q);
    if isempty(I)
        ue  = Inf;
    else
        ue  = min(N(J(I))) / V2R;
    end
end

author: spm


text:
function [] = display_stop_reason(epoch, infos, options, disp_name, reason, max_reached_flag)
% Function to display stop reason
%
% This file is part of NMFLibrary.
%
% Created by H.Kasai on July 14, 2022
%

    if options.verbose > 0

        if options.verbose == 1 && ~max_reached_flag
            fprintf(' (incomplete)\n');
        end

        fprintf('# %s: %s', disp_name, reason);
    end

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function r = save_mgh2(vol, fname, M, mr_parms);
%
% save_mgh2(vol,fname, M, <mr_parms>);
%
% M is the 4x4 vox2ras transform such that
% y(i1,i2,i3), xyz = M*[i1 i2 i3 1] where the
% indicies are 0-based
%
% mr_parms = [tr flipangle te ti]
%
% See also: load_mgh, vox2ras_0to1
%
%

%
% save_mgh2.m
%
% Original Author: Doug Greve
% CVS Revision Info:
%    $Author: nicks $
%    $Date: 2011/03/02 00:04:13 $
%    $Revision: 1.4 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

r = 1;

if(nargin < 2 | nargin > 4)
  msg = 'USAGE: save_mgh2(vol,fname,M)';
  return;
end

if(exist('mr_parms')~=1) mr_parms = []; end
if(isempty(mr_parms))   mr_parms = [0 0 0 0]; end
if(length(mr_parms) ~= 4)
  fprintf('ERROR: mr_parms length = %d, must be 4\n', ...
	  lenght(mr_parms));
  return;
end

% These dont appear to be used %
MRI_UCHAR =  0 ;
MRI_INT =    1 ;
MRI_LONG =   2 ;
MRI_FLOAT =  3 ;
MRI_SHORT =  4 ;
MRI_BITMAP = 5 ;
MRI_TENSOR = 6 ;

fid = fopen(fname, 'wb', 'b') ;
if(fid == -1)
  fprintf('ERROR: could not open %s for writing\n',fname);
  return;
end

[ndim1,ndim2,ndim3,rows,cols] = size(vol) ;
fwrite(fid, 1, 'int') ;		% magic #
fwrite(fid, ndim1, 'int') ;
fwrite(fid, ndim2, 'int') ;
fwrite(fid, ndim3, 'int') ;
fwrite(fid, 1, 'int') ;		% # of frames
if (ndims(vol) == 5)
  is_tensor = 1 ;
  fwrite(fid, MRI_TENSOR, 'int') ; % type = MRI_TENSOR
else
  is_tensor = 0 ;
  fwrite(fid, MRI_FLOAT, 'int') ;  % type = MRI_FLOAT
end

%%?????????????%%%
fwrite(fid, 1, 'int') ;          % dof (not used)
dof = fread(fid, 1, 'int') ;

UNUSED_SPACE_SIZE= 256;
USED_SPACE_SIZE = (3*4+4*3*4);  % space for ras transform

MdcD = M(1:3,1:3);
delta = sqrt(sum(MdcD.^2));

Mdc = MdcD./repmat(delta,[3 1]);
Pcrs_c = [ndim1/2 ndim2/2 ndim3/2 1]'; %'
Pxyz_c = M*Pcrs_c;
Pxyz_c = Pxyz_c(1:3);

fwrite(fid, 1,      'short') ;       % ras_good_flag = 1
fwrite(fid, delta,  'float32') ;
fwrite(fid, Mdc,    'float32') ;
fwrite(fid, Pxyz_c, 'float32') ;

unused_space_size = UNUSED_SPACE_SIZE-2 ;
unused_space_size = unused_space_size - USED_SPACE_SIZE ;
fwrite(fid, zeros(unused_space_size,1), 'char') ;

fwrite(fid,vol,'float32');

fwrite(fid, mr_parms, 'float32') ;
fclose(fid) ;

r = 0;

return;


text:
function skelVargplvmModify(handle, yVal, chan, skel, padding, zeroIndices)

% SKELMODIFY Update visualisation of skeleton data.
% FORMAT
% DESC updates a skeleton representation in a 3-D plot.
% ARG handle : a vector of handles to the structure to be updated.
% ARG channels : the channels to update the skeleton with.
% ARG skel : the skeleton structure.
%
% SEEALSO : skelVisualise
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006
% MODIFICATIONS: Andreas C. Damianou, 2014

% SHEFFIELDML

if nargin<5
  padding = 0;
end
channels = demCmu35VargplvmLoadChannels(yVal,skel);
channels = [channels zeros(1, padding)];
vals = skel2xyz(skel, channels);
connect = skelConnectionMatrix(skel);

indices = find(connect);
[I, J] = ind2sub(size(connect), indices);

set(handle(1), 'Xdata', vals(:, 1), 'Ydata', vals(:, 3), 'Zdata', ...
                 vals(:, 2));
%/~
%set(handle(1), 'visible', 'on')
%~/

for i = 1:length(indices)
  set(handle(i+1), 'Xdata', [vals(I(i), 1) vals(J(i), 1)], ...
            'Ydata', [vals(I(i), 3) vals(J(i), 3)], ...
            'Zdata', [vals(I(i), 2) vals(J(i), 2)]);
end

function [vals, connect] = wrapAround(vals, lim, connect);

quot = lim(2) - lim(1);
vals = rem(vals, quot)+lim(1);
nVals = floor(vals/quot);
for i = 1:size(connect, 1)
  for j = find(connect(i, :))
    if nVals(i) ~= nVals(j)
      connect(i, j) = 0;
    end
  end
end

author: SheffieldML

repo: vargplvm


text:
function [mu sigma2] = estimateGaussian(X)
%ESTIMATEGAUSSIAN This function estimates the parameters of a
%Gaussian distribution using the data in X
%   [mu sigma2] = estimateGaussian(X),
%   The input X is the dataset with each n-dimensional data point in one row
%   The output is an n-dimensional vector mu, the mean of the data set
%   and the variances sigma^2, an n x 1 vector
%

% Useful variables
[m, n] = size(X);

% You should return these values correctly
mu = zeros(n, 1);
sigma2 = zeros(n, 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the mean of the data and the variances
%               In particular, mu(i) should contain the mean of
%               the data for the i-th feature and sigma2(i)
%               should contain variance of the i-th feature.
%

mu = mean(X);
sigma2 = var(X, 1);

% =============================================================

end

author: merwan

repo: ml-class

save_path: github-repos/MATLAB/merwan-ml-class


text:
function [change,error_handle]=dicomrt_sortct(filelist,xlocation,ylocation,zlocation,filename)
% dicomrt_sortct(filelist,xlocation,ylocation,zlocation,filename)
%
% Sort CT data set specified in filename.
%
% See also dicomrt_loaddose, dicomrt_loadct, dicomrt_loadctlist
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

% Set parameters
counter=0;
error_handle=0;
match=0;

% sort ct slices with respect to their z location
[zlocation_sorted,sorted_index]=sort(zlocation);

% Progress bar
h = waitbar(0,['Sorting progress:']);
set(h,'Name','dicomrt_sortct: sorting CT images');

diff_location=dicomrt_mmdigit(zlocation_sorted-zlocation,6);

if find(diff_location) % ct slices are not sorted
    disp('Message: slices are not sorted!');
    change=1;
    % sort filelist accordingly
    filelist_sorted=' ';
    for i=1:length(sorted_index)
        counter=counter+1;
        filelist_sorted=char(filelist_sorted,filelist(sorted_index(i),:));
        xlocation_sorted(i)=xlocation(sorted_index(i));
        ylocation_sorted(i)=ylocation(sorted_index(i));
        waitbar(i/length(sorted_index),h);
    end

    xlocation_sorted=xlocation_sorted';
    ylocation_sorted=ylocation_sorted';
    filelist_sorted(1,:)=[];
    disp('Message: ct slices are now sorted!');
    % finds duplicates
    duplicate_index=0;
    [duplicate]=diff(zlocation_sorted);
    [duplicate_index]=find(duplicate==0);
    if length(duplicate_index)~=0
        warning('dicomrt_sortct: Images with the same zlocation were found. You can exit now or proceed to delete duplicates.');
        leaveoption=input('Exit now ? (Y/N) [N]','s');
        if leaveoption == 'Y' | leaveoption == 'y';
            disp('OK. Check your file list and try again. ');
            error_handle=1;return;
        else
            duplicates_name=' ';
            filelist_sorted_scrub=' ';
            for i=1:length(duplicate_index)
                duplicates_name=char(duplicates_name,filelist_sorted(duplicate_index(i),:));
                filelist_sorted(duplicate_index(i),1:10)='duplicate-'; % mark duplicates
                xlocation_sorted(duplicate_index(i))=pi;
                ylocation_sorted(duplicate_index(i))=pi;
            end
            for i=1:length(filelist_sorted)
                if isequal(filelist_sorted(i,1:10),'duplicate-')~=1
                    filelist_sorted_scrub=char(filelist_sorted_scrub,filelist_sorted(i,:));
                else
                end
            end % at this point duplicates are deleted in new filelist
                        filelist_sorted_scrub(1,:)=[];
            disp('The following duplicates has been deleted from file list:');
            duplicates_name(1,:)=[];
            duplicates_name
            %
            % now sorting xlocation and ylocation accordingly to what done before
            % (prepare for scout images identification)
            %
            for i=1:length(xlocation_sorted)
                if xlocation_sorted(i)~=pi
                    xlocation_sorted_scrub=xlocation_sorted(i);
                end
            end
            for i=1:length(ylocation_sorted)
                if ylocation_sorted(i)~=pi
                    ylocation_sorted_scrub=ylocation_sorted(i);
                end
            end
        end
    else
        filelist_sorted_scrub=filelist_sorted;
    end

    % Export sorted/scrubed filelist
    newfilename=[filename,'.sort.txt'];
    newfile=fopen(newfilename,'w');

    for i=1:size(filelist_sorted_scrub,1)
        fprintf(newfile,'%s',filelist_sorted_scrub(i,:)); fprintf(newfile,'\n');
    end

    fclose(newfile);
    disp(['A new file list has been written by dicomrt_sortct with name: ',newfilename]);
    disp('This file will be used to import ct data instead');

else % ct slices are already sorted
    change=0;
end

% Close progress bar
close(h);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function varargout = scl_inter(varargin)
% Format
% For getting the value
% dat = scl_inter(obj)
%
% For setting the value
% obj = scl_inter(obj,dat)
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id$

if nargin==2,
    varargout{1} = asgn(varargin{:});
elseif nargin==1,
    varargout{1} = ref(varargin{:});
else
    error('Wrong number of arguments.');
end;
return;

function dat = ref(obj)
dat = obj.scl_inter;
return;

function obj = asgn(obj,dat)
if isnumeric(dat), % && numel(dat)<=1,
    obj.scl_inter = double(dat);
else
    error('"scl_inter" must be numeric.');
end;
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef L2Norm < dagnn.ElementWise

  properties
    % keep the name "param" for  compatibility with previous implementations
    param = 1e-10 ;
  end

  methods
    function outputs = forward(obj, inputs, params)
      outputs{1} = vl_nnl2norm(inputs{1}, 'epsilon', obj.param) ;
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
      derInputs{1} = vl_nnl2norm(inputs{1}, derOutputs{1}, 'epsilon', obj.param) ;
      derParams = {} ;
    end

    function rfs = getReceptiveFields(obj)
      rfs.size = [1 1] ;
      rfs.stride = [1 1] ;
      rfs.offset = [1 1] ;
    end

    function obj = L2Norm(varargin)
      obj.load(varargin) ;
    end
  end
end

repo: MFT


text:
function diaphony ( input_filename )

%*****************************************************************************80
%
%% MAIN is the main program for DIAPHONY.
%
%  Discussion:
%
%    DIAPHONY reads a table dataset and applies the diaphony test.
%
%  Usage:
%
%    diaphony ( 'table_file' )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2012
%
%  Author:
%
%    John Burkardt
%
  verbose = 0;

  if ( verbose )
    timestamp ( );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'DIAPHONY:\n' );
    fprintf ( 1, '  MATLAB version\n' );
    fprintf ( 1, '  Compute the diaphony of a point set.\n' );
  end
%
%  Get the filename.
%
  if ( nargin < 1 )
    fprintf ( 1, '\n' );
    input_filename = input ( 'Enter the name of the input file.' );
  end
%
%  Get the data size.
%
  [ dim_num, point_num ] = r8mat_header_read ( input_filename );

  if ( verbose )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  The spatial dimension is %d\n', dim_num );
    fprintf ( 1, '  The number of points is %d\n', point_num );
  end
%
%  Read the data.
%
  points = r8mat_data_read ( input_filename, dim_num, point_num );

  if ( min ( min ( points ) ) < 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'DIAPHONY - Fatal error!\n' );
    fprintf ( 1, '  At least one coordinate of a point is less than 0!\n' );
    error ( 'DIAPHONY - Fatal error!' )
  elseif ( 1.0 < max ( max ( points ) ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'DIAPHONY - Fatal error!\n' );
    fprintf ( 1, '  At least one coordinate of a point is greater than 1!\n' );
    error ( 'DIAPHONY - Fatal error!' )
  end
%
%  Analyze the data.
%
  d = diaphony_compute ( dim_num, point_num, points );

  e = 1.0 / sqrt ( point_num );
  de = d / e;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  File  M  N  Diaphony  1/sqrt(N)  D/sqrt(N)\n' );
  fprintf ( 1, '  %s  %2d  %5d  %14.6g  %14.6g  %14.6g\n', ...
    input_filename, dim_num, point_num, d, e, de );
%
%  Terminate.
%
  if ( verbose )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'DIAPHONY:\n' );
    fprintf ( 1, '  Normal end of execution.\n' );
    fprintf ( 1, '\n' );
    timestamp ( )
  end

  return
end
function d = diaphony_compute ( dim_num, point_num, x )

%*****************************************************************************80
%
%% DIAPHONY_COMPUTE evaluates the diaphony of a N-dimensional point set.
%
%  Discussion:
%
%    The diaphony is analogous to, and related to, the discrepancy,
%    and is a measure of how well spread a set of point is.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 May 2003
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Peter Heelekalek, Harald Niederreiter,
%    The Weighted Spectral Test: Diaphony,
%    ACM Transactions on Modeling and Computer Simulation,
%    Volume 8, Number 1, January 1998, pages 43-60.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, real X(DIM_NUM,POINT_NUM), the point set, which is
%    presumed to lie in the DIM_NUM dimensional unit hypercube.
%
%    Output, real D, the value of the diaphony.
%
  d = 0.0;

  for i = 1 : point_num
    for j = 1 : point_num
      z(1:dim_num) = x(1:dim_num,i) - x(1:dim_num,j);
      for k = 1 : dim_num
        z(k) = r8_modp ( z(k), 1.0 );
      end
      d = d - 1.0 + prod ( ...
        ( 1.0 + 2.0 * pi^2 * ( z(1:dim_num).^2 - z(1:dim_num) + 1.0 / 6.0 ) ) );
    end
  end

  bot = point_num^2 * ( ( 1.0 + pi^2 / 3.0 )^dim_num - 1.0 );

  d = d / bot;

  d = sqrt ( d );

  return
end
function column_num = file_column_count ( input_file_name )

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

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

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( s_len_trim ( line ) == 0 )

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

    else
      got_one = TRUE;
      break;
    end

  end

  fclose ( input_unit );

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

  column_num = s_word_count ( line );

  return
end
function row_num = file_row_count ( input_file_name )

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

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

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

  record_num = 0;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

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

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

  end

  fclose ( input_unit );

  return
end
function value = r8_modp ( x, y )

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

  value = mod ( x, y );

  if ( value < 0.0 )
    value = value + abs ( y );
  end

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

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

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

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

  input_unit = fopen ( input_filename );

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

  table = zeros(m,n);

  i = 0;

  while ( i < n )

    line = fgets ( input_unit );

    if ( line == -1 )
      break;
    end

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

    elseif ( s_len_trim ( line ) == 0 )

    else

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

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

    end

  end

  fclose ( input_unit );

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

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

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

  n = file_row_count ( input_filename );

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

  return
end
function len = s_len_trim ( s )

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

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

  return
end
function word_num = s_word_count ( s )

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

  word_num = 0;
  s_length = length ( s );

  if ( s_length <= 0 )
    return;
  end

  blank = TRUE;

  for i = 1 : s_length

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

  end

  return
end
function timestamp ( )

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

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef MimModelController < CoreBaseClass

    properties (Access = private)
        Mim
        ModelMap
    end

    methods
        function obj = MimModelController(mim)
            obj.Mim = mim;
            obj.ModelMap = MimModelMap(mim);
        end

        function [value, hash] = getValue(obj, modelId)
            [value, hash] = obj.ModelMap.getValue(modelId);
            obj.ModelMap.autoUpdate();
        end

        function setValue(obj, modelId, value)
            obj.ModelMap.setValue(modelId, value);
            obj.ModelMap.autoUpdate();
        end

        function modelId = buildModelId(obj, modelClassName, parameters)
            modelId = obj.ModelMap.buildModelId(modelClassName, parameters);
        end

        function clear(obj)
            obj.ModelMap.clear();
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [minconpos, mincon]=API(fun1,LB,UB,NumAnts,Nmoves,LocalMoves,nonlinfun,rpmax,RES,varargin)
%
% APA, 2003

warning off
%  API Finds a unconstrained or constrained minimum of a function of several variables.
%     API solves problems of the form:
%         min F(X)  subject to:  g(x) <= 0,  h(x) = Beq (inequality & equality constraints)
%          X                     LB <= X <= UB (boundary constraints)
%
%     X=API(FUN,LB,UB,NumAnts,Nmoves,LocalMoves,NONLNCON,rpmax) minimizes FUN
%     subject to inequality and equality constraints in NONLINCON. (Set NONLINCON=[]
%     and rpmax=[] if no inequality and equality constraints exist).
%     NumAnts = Number of API Ants to be used.
%     Nmoves = Number of Nest Moves.
%     LocalMoves = Number of Local Moves for API Ants.
%     rpmax = Augmented Lagrangian parameter.
%
%     X=API(FUN,LB,UB,NumAnts,Nmoves,LocalMoves,NONLNCON,rpmax,P1,P2,...) passes the
%     problem-dependent parameters P1,P2,... directly to the functions FUN
%     and NONLCON: feval(FUN,X,P1,P2,...) and feval(NONLCON,X,P1,P2,...).
%
% Example:
% FUN='testbanana';LB=[-5 -5];UB=[3 3];NumAnts=20;Nmoves=15;LocalMoves=25;
% NONLNCON='nonlinbanana';rpmax=5;
% X=API(FUN,LB,UB,NumAnts,Nmoves,LocalMoves,NONLNCON,rpmax)
%
% Listing of Objective and Constraint for the above exxample
%
% function f=testbanana(x)
% f=100*(x(1)^2-x(2))^2 + (1-x(1))^2;
%
% function [g,h]=nonlinbanana(x)
% x1=x(1); x2=x(2);
% %%% Equality Constraints %%%%%%%%
% h=0;
% %%% Inequality Constraints %%%%%%
% g(1,:)=-(x1.^2+x2.^2-1.2);
% g(2,:)=-(1.5-(x1-.2)^2-x2.^2);

if isempty(nonlinfun)

    xb=[LB(:) UB(:)];

    % Specify the starting point
    Ninit = (LB+UB)/2;

    n=length(Ninit);

    %%% **************--- Unconstrained Loop ---****************

    N=Ninit;% Initial Nest Location
    min1=inf;
    disp('Nest Move    Objective')
    for ii=1:Nmoves   % Number of Nest Moves

        Smem=[];
        for num=1:NumAnts
            if(num==1)
                Asite=0.01;
            else
                Asite=0.01*(1/0.01)^(num/NumAnts);
            end
            for i=1:2
                for j=1:n
                    r=rand;
                    S(j)=N(j)+((-0.5+r)*Asite*(xb(j,2)-xb(j,1)));
                    while(~((xb(j,1)<=S(j))&(xb(j,2)>=S(j))))
                        r=rand;
                        S(j)=N(j)+((-0.5+r)*Asite*(xb(j,2)-xb(j,1)));
                    end
                end
                Smem=[Smem;S];
            end
        end
        % Fun1='fun';
        %%%%%%%%%%%% Tandem Running %%%%%%%%%%%%%%%
        a1=1+floor(NumAnts*rand);
        a2=1+floor(NumAnts*rand);
        f1=feval(fun1,Smem(2*a1-1,:),varargin{:});
        f2=feval(fun1,Smem(2*a1,:),varargin{:});
        f3=feval(fun1,Smem(2*a2-1,:),varargin{:});
        f4=feval(fun1,Smem(2*a2,:),varargin{:});

        if(f1<f2)
            fpos1=2*a1-1;
            fmin1=f1;
        else
            fpos1=2*a1;
            fmin1=f2;
        end
        if(f3<f4)
            fpos2=2*a2-1;
            fmin2=f3;
        else
            fpos2=2*a2;
            fmin2=f4;
        end
        if(fmin1<fmin2)
            Smem(fpos2,:)=Smem(fpos1,:);
        else
            Smem(fpos1,:)=Smem(fpos2,:);
        end
        %%%%%%%%%%%%%%tandem running loop ends%%%%%%%%%%%%%%%%%%

        for num=1:2:(2*NumAnts-1)

            if(num==1)
                Alocal=0.1*0.01;
            else
                Alocal=(1/10)*0.01*(1/0.01)^(((num+1)/2)/NumAnts);
            end
            pos1=num;
            pos2=num+1;

            for i=1:LocalMoves   % Number of Local Sites % Local Search Loop Starts here

                Spos=Smem(num,:);
                for j=1:n
                    r=rand;
                    Sloc(j)=Spos(j)+(-0.5+r)*Alocal*(xb(j,2)-xb(j,1));
                    while(~((xb(j,1)<=Sloc(j))&(xb(j,2)>=Sloc(j))))
                        r=rand;
                        Sloc(j)=Spos(j)+(-0.5+r)*Alocal*(xb(j,2)-xb(j,1));
                    end
                end

                V1=feval(fun1,Spos,varargin{:});
                V2=feval(fun1,Sloc,varargin{:});
                if(V1>V2)
                    Smem(pos1,:)=Sloc;
                    if(min1>V2)
                        minpos=Sloc;
                        min1=V2;
                    end
                else
                    temppos=Smem(pos1,:);
                    Smem(pos1,:)=Smem(pos2,:);
                    Smem(pos2,:)=temppos;
                    if(min1>V1)
                        minpos=Spos;
                        min1=V1;
                    end
                end

            end

        end
        %! Put-in the residue
        if min1<RES
            return
        end
        N=minpos;
        % disp('[Nest Move    Objective]')
        disp(['    ',num2str(ii),'        ',num2str(min1)])

    end

    %%%***********--- Unconstrained Loop Ends ---***************
    minconpos=minpos;
    mincon=min1;

else

    %API1Gen.m
    ftemp=fun1;nonlinfuntemp=nonlinfun;
    global rp mu1 si1 fun1 nonlinfun
    fun1=ftemp;nonlinfun=nonlinfuntemp;
    xb=[LB(:) UB(:)];

    % Specify the starting point
    Ninit = (LB+UB)/2;

    n=length(Ninit);
    rp=1;
    [g1, h1]=feval('nonlinAPI',Ninit,varargin{:});
    mu1=zeros(length(g1),1);
    si1=zeros(length(h1),1);

    mincon=inf;
    NITER=0;
    while(rp<rpmax)

        NITER=NITER+1;

        %%% **************--- Unconstrained Loop ---****************

        N=Ninit;% Initial Nest Location
        min1=inf;

        for ii=1:Nmoves   % Number of Nest Moves

            Smem=[];
            for num=1:NumAnts
                if(num==1)
                    Asite=0.01;
                else
                    Asite=0.01*(1/0.01)^(num/NumAnts);
                end
                for i=1:2
                    for j=1:n
                        r=rand;
                        S(j)=N(j)+((-0.5+r)*Asite*(xb(j,2)-xb(j,1)));
                        while(~((xb(j,1)<=S(j))&(xb(j,2)>=S(j))))
                            r=rand;
                            S(j)=N(j)+((-0.5+r)*Asite*(xb(j,2)-xb(j,1)));
                        end
                    end
                    Smem=[Smem;S];
                end
            end
            %%%%%%%%%%%% Tandem Running %%%%%%%%%%%%%%%
            a1=1+floor(NumAnts*rand);
            a2=1+floor(NumAnts*rand);
            f1=feval(@funAPI,Smem(2*a1-1,:),varargin{:});
            f2=feval(@funAPI,Smem(2*a1,:),varargin{:});
            f3=feval(@funAPI,Smem(2*a2-1,:),varargin{:});
            f4=feval(@funAPI,Smem(2*a2,:),varargin{:});

            if(f1<f2)
                fpos1=2*a1-1;
                fmin1=f1;
            else
                fpos1=2*a1;
                fmin1=f2;
            end
            if(f3<f4)
                fpos2=2*a2-1;
                fmin2=f3;
            else
                fpos2=2*a2;
                fmin2=f4;
            end
            if(fmin1<fmin2)
                Smem(fpos2,:)=Smem(fpos1,:);
            else
                Smem(fpos1,:)=Smem(fpos2,:);
            end
            %%%%%%%%%%%%%%tandem running loop ends%%%%%%%%%%%%%%%%%%

            for num=1:2:(2*NumAnts-1)

                if(num==1)
                    Alocal=0.1*0.01;
                else
                    Alocal=(1/10)*0.01*(1/0.01)^(((num+1)/2)/NumAnts);
                end
                pos1=num;
                pos2=num+1;

                for i=1:LocalMoves   % Number of Local Sites % Local Search Loop Starts here

                    Spos=Smem(num,:);
                    for j=1:n
                        r=rand;
                        Sloc(j)=Spos(j)+(-0.5+r)*Alocal*(xb(j,2)-xb(j,1));
                        while(~((xb(j,1)<=Sloc(j))&(xb(j,2)>=Sloc(j))))
                            r=rand;
                            Sloc(j)=Spos(j)+(-0.5+r)*Alocal*(xb(j,2)-xb(j,1));
                        end
                    end

                    V1=feval(@funAPI,Spos,varargin{:});
                    V2=feval(@funAPI,Sloc,varargin{:});
                    if(V1>V2)
                        Smem(pos1,:)=Sloc;
                        if(min1>V2)
                            minpos=Sloc;
                            min1=V2;
                        end
                    else
                        temppos=Smem(pos1,:);
                        Smem(pos1,:)=Smem(pos2,:);
                        Smem(pos2,:)=temppos;
                        if(min1>V1)
                            minpos=Spos;
                            min1=V1;
                        end
                    end

                end

            end

            N=minpos;

        end

        %%%***********--- Unconstrained Loop Ends ---***************

        %%% --- Augmented Lagrangian Method ---
        if(rp>1)
            mincon=feval(@funAPI,minconpos,varargin{:});
        end
        if(min1<mincon)
            minconpos=minpos;
            mincon=min1;
        end

        [temp1, temp2]=feval(nonlinfun, minconpos,varargin{:});
        mu1 = mu1 + rp.*temp1(:);
        si1 = si1 + rp.*temp2(:);
        %%% ------ End of Augmented Lagrangian -----

        rp=rp*1.15;
        if NITER==1
            disp('Iteration     Objective    Maximum Constraint')
            disp(['     ',num2str(NITER),'         ',num2str(mincon),'         ',num2str(max(temp1))])
        else
            disp(['     ',num2str(NITER),'         ',num2str(mincon),'         ',num2str(max(temp1))])
        end

    end %%% End of Constrained Loop (i.e. rp)

    %%%% function API1gen(fun1,nonlinAPI) ends

end

function Faug=funAPI(x,varargin)

global mu1 si1 rp nonlinfun fun1
if ~isempty(varargin)
    f=feval(fun1,x,varargin{:});
else
    f=feval(fun1,x);
end
if ~isempty(varargin)
    [zeta1, h1]=feval(nonlinfun,x,varargin{:});
else
    [zeta1, h1]=feval(nonlinfun,x);
end
Fp = f + rp/2*(sum(zeta1.^2) + sum(h1.^2));
Faug = Fp + sum(mu1.*zeta1) + sum(si1.*h1);

%%%%% fun1(x) ends %%%%

function [zeta1, h1]=nonlinAPI(x,varargin)
global mu1 si1 rp h1 nonlinfun

if ~isempty(varargin)
    [g1, h1]=feval(nonlinfun,x,varargin{:});
else
    [g1, h1]=feval(nonlinfun,x);
end
if(rp==1)
    zeta1=g1;
else
    zeta1=max(g1,-mu1/rp);
end

%%%%% function nonlinAPI(x) ends