text:
function geometry_test0126 ( )

%*****************************************************************************80
%
%% TEST0126 tests SPHERE_CAP_VOLUME_3D, SPHERE_CAP_VOLUME_ND.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  dim_num = 3;
  ntest = 12;

  center(1:3) = [ 0.0, 0.0, 0.0 ];
  r = 1.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0126\n' );
  fprintf ( 1, '  SPHERE_CAP_VOLUME_3D computes the volume of a\n' );
  fprintf ( 1, '    spherical cap, defined by a plane that cuts the\n' );
  fprintf ( 1, '    sphere to a thickness of H units.\n' );
  fprintf ( 1, '  SPHERE_CAP_VOLUME_ND does the same operation,\n' );
  fprintf ( 1, '    but in N dimensions.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Using a radius R = %f\n', r );

  fprintf ( 1, '\n' );
  fprintf ( 1, '        H           Cap        Cap\n' );
  fprintf ( 1, '                    volume_3d  volume_nd\n' );
  fprintf ( 1, '\n' );

  for i = 0 : ntest

    h = 2.0 * r * i / ntest;

    volume1 = sphere_cap_volume_3d ( r, h );

    volume2 = sphere_cap_volume_nd ( dim_num, r, h );

    fprintf ( 1, '  %10f  %10f  %10f\n', h, volume1, volume2 );

  end

  return
end


text:
function [ringWedgeScans ringWedgeFields] = readRingWedgeScans(dlgTitle);
%
%   [ringWedgeScans ringWedgeFields]  = readRingWedgeScans(dlgTitle);
%
% Author: Wandell, Brewer
% Purpose:
%    Utility to read ring and wedge scan numbers
%    ringWedgeScans(1) is the scan num of the ring
%    ringWedgeScans(2) is the scan num of the wege
%
% ras, 04/2009: updated with slots for the data field (co, amp, ph, map)
% for each scan as well. (Important for pRF data.)
if notDefined('dlgTitle'), dlgTitle='Read ring and wedge scan numbers'; end

dlg(1).fieldName = 'ringScan';
dlg(end).style = 'number';
dlg(end).value = 1;
dlg(end).string = 'Eccentricity (ring) data scan?';

dlg(end+1).fieldName = 'ringField';
dlg(end).style = 'popup';
dlg(end).list = {'co' 'amp' 'ph' 'map'};
dlg(end).value = 'map';
dlg(end).string = 'Polar Angle (wedge) data field?';

dlg(end+1).fieldName = 'wedgeScan';
dlg(end).style = 'number';
dlg(end).value = 1;
dlg(end).string = 'Eccentricity (ring) data scan?';

dlg(end+1).fieldName = 'wedgeField';
dlg(end).style = 'popup';
dlg(end).list = {'co' 'amp' 'ph' 'map'};
dlg(end).value = 'ph';
dlg(end).string = 'Polar Angle (wedge) data field?';

[resp ok] = generalDialog(dlg, dlgTitle);

if ~ok, error('User Aborted.');		end

ringWedgeScans = [resp.ringScan resp.wedgeScan];
ringWedgeFields = {resp.ringField resp.wedgeField};

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% BOX_INTERSECT Given sets of axis-aligne boxes determine which pairs overlap.
%
% I = box_intersect(A1,A2)
%
% Inputs:
%   A1  #A by dim list of minimum corners
%   A2  #A by dim list of maximum corners
% Outputs:
%   I  #I by 2 list of indices into A such that box I(i,1) of set A intersects
%     with box I(i,2) of set A.
%
% I = box_intersect(A1,A2,B1,B2)
%
% Inputs:
%   A1  #A by dim list of minimum corners
%   A2  #A by dim list of maximum corners
%   B1  #B by dim list of minimum corners
%   B2  #B by dim list of maximum corners
% Outputs:
%   I  #I by 2 list of indices into A such that box I(i,1) of set A intersects
%     with box I(i,2) of set B.

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function out=factorial(x)

precision=x(1).precision;
out_rval=cell(size(x));
out_ival=cell(size(x));

for ii=1:numel(x)
 imag=false;
 [xrval,xival]=getVals(x,ii);
 if hasimag(xival), imag=true; end
 if imag
  error('Factorial must be on nonegative integers');
 else
  out_rval{ii}=mpfr_factorial(precision,xrval);
 end
end % for ii=1:max(ex,
out=class(struct('rval',out_rval,...
                  'ival',out_ival,...
                  'precision',precision),'mp');

author: opencobra

repo: cobratoolbox


text:
function test15
%TEST15 test cs_amd
%
% Example:
%   test15
% See also: testall

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

rand ('state', 0) ;
randn ('state', 0) ;
clf

for trials = 1:100
    n = fix (200 * rand (1)) ;
    d = 0.05 * rand (1) ;
    A = sprandn (n, n, d) ;

    % add a randomly placed dense column
    k = fix (n * rand (1)) ;
    k = max (1, k) ;
    k = min (n, k) ;
    A (:,k) = 1 ;

    try
        p0 = amd (A) ;
    catch
        p0 = symamd (A) ;
    end
    p1 = cs_amd (A) ;

    if (any (sort (p1) ~= 1:n))
        error ('not perm!') ;
    end

    C = A+A' + speye (n) ;
    lnz0 = sum (symbfact (C (p0,p0))) ;
    lnz1 = sum (symbfact (C (p1,p1))) ;
    subplot (1,3,1) ; spy (C)
    subplot (1,3,2) ; spy (C (p0,p0))
    subplot (1,3,3) ; spy (C (p1,p1))
    fprintf ('n %4d nz %6d lnz %6d %6d\n', n, nnz(A), lnz0, lnz1) ;
    drawnow

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef ImdbFCN < ImdbCalvin
    % ImdbFCN
    %
    % Standard Imdb for FCN experiments.
    %
    % Copyright by Holger Caesar, 2016

    properties
        % Set in constructor
        dataset

        batchOpts = struct();
    end
    methods
        function obj = ImdbFCN(dataset, dataDir, nnOpts)
            % Call default constructor
            obj = obj@ImdbCalvin();
            obj.dataset = dataset;

            % FCN-specific
            obj.batchOpts.imageSize = [512, 512] - 128;
            obj.batchOpts.labelStride = 1;

            obj.batchOpts.imageFlipping = true;
            obj.batchOpts.rgbMean = single([128; 128; 128]);
            obj.batchOpts.classWeights = [];
            obj.batchOpts.imageNameToLabelMap = @(imageName) obj.dataset.getImLabelMap(imageName);
            obj.batchOpts.translateLabels = false;
            obj.batchOpts.maskThings = false;

            obj.batchOpts.useInvFreqWeights = false;

            % Dataset-specific
            obj.batchOpts.vocAdditionalSegmentations = true;
            obj.batchOpts.vocEdition = '11';
            obj.batchOpts.dataDir = dataDir;

            % Load VOC-style IMDB
            obj.loadImdb(nnOpts);
        end

        function loadImdb(obj, nnOpts)
            % loadImdb(obj, nnOpts)
            %
            % Creates or load the current VOC-style IMDB.

            %%% VOC specific
            if strStartsWith(obj.dataset.name, 'VOC')
                % Get PASCAL VOC segmentation dataset plus Berkeley's additional segmentations
                imdbPath = fullfile(nnOpts.expDir, 'imdbVoc.mat');
                if exist(imdbPath, 'file')
                    % Load imdbVoc
                    load(imdbPath, 'imdbVoc', 'rgbStats');
                else
                    imdbVoc = vocSetup('dataDir', obj.batchOpts.dataDir, ...
                        'edition', obj.batchOpts.vocEdition, ...
                        'includeTest', false, ...
                        'includeSegmentation', true, ...
                        'includeDetection', false);
                    if obj.batchOpts.vocAdditionalSegmentations
                        imdbVoc = vocSetupAdditionalSegmentations(imdbVoc, 'dataDir', obj.batchOpts.dataDir);
                    end
                    rgbStats = getDatasetStatistics(imdbVoc);

                    save(imdbPath, 'imdbVoc', 'rgbStats');
                end
                obj.batchOpts = structOverwriteFields(obj.batchOpts, imdbVoc);

                obj.batchOpts.rgbMean = rgbStats.rgbMean;
                obj.batchOpts.rgbMean = reshape(obj.batchOpts.rgbMean, [1 1 3]);
                obj.batchOpts.translateLabels = true;
                obj.batchOpts.imageNameToLabelMap = @(imageName) imread(sprintf(obj.batchOpts.paths.classSegmentation, imageName));
            else
                %%% Other datasets
                % Get labels and image path
                obj.batchOpts.classes.name = obj.dataset.getLabelNames();
                obj.batchOpts.paths.image = fullfile(obj.dataset.getImagePath(), sprintf('%%s%s', obj.dataset.imageExt));

                % Get trn + tst/val images
                imageListTrn = obj.dataset.getImageListTrn();
                imageListTst = obj.dataset.getImageListTst();

                % Remove images without labels
                missingImageIndicesTrn = obj.dataset.getMissingImageIndices('train');
                imageListTrn(missingImageIndicesTrn) = [];
                % TODO: is it a good idea to remove test images?
                % (only doing it on non-competitive EdiStuff
                if isa(obj.dataset, 'EdiStuffDataset') || isa(obj.dataset, 'EdiStuffSubsetDataset')
                    missingImageIndicesTst = obj.dataset.getMissingImageIndices('test');
                    imageListTst(missingImageIndicesTst) = [];
                end
                imageCountTrn = numel(imageListTrn);
                imageCountTst = numel(imageListTst);

                obj.batchOpts.images.name = [imageListTrn; imageListTst];
                obj.batchOpts.images.segmentation = true(imageCountTrn+imageCountTst, 1);
                obj.batchOpts.images.set = nan(imageCountTrn+imageCountTst, 1);
                obj.batchOpts.images.set(1:imageCountTrn) = 1;
                obj.batchOpts.images.set(imageCountTrn+1:end) = 2;

                obj.batchOpts.rgbMean = obj.dataset.getMeanColor();
                obj.batchOpts.rgbMean = reshape(obj.batchOpts.rgbMean, [1 1 3]);
                obj.batchOpts.translateLabels = false;
                obj.batchOpts.imageNameToLabelMap = @(imageName) obj.dataset.getImLabelMap(imageName);
            end

            % Specify level of supervision for each train image
            if ~nnOpts.misc.weaklySupervised
                % Full supervision
                obj.batchOpts.images.isFullySupervised = true(numel(obj.batchOpts.images.name), 1);
            elseif ~semiSupervised
                % Weak supervision
                obj.batchOpts.images.isFullySupervised = false(numel(obj.batchOpts.images.name), 1);
            else
                % Semi supervision: Set x% of train and all val to true
                obj.batchOpts.images.isFullySupervised = true(numel(obj.batchOpts.images.name), 1);
                if isa(obj.dataset, 'EdiStuffDataset')
                    selWS = find(~ismember(obj.batchOpts.images.name, obj.dataset.datasetFS.getImageListTrn()));
                    assert(numel(selWS) == 18431);
                else
                    selTrain = find(obj.batchOpts.images.set == 1);
                    selTrain = selTrain(randperm(numel(selTrain)));
                    selWS = selTrain((selTrain / numel(selTrain)) >= semiSupervisedRate);
                end
                obj.batchOpts.images.isFullySupervised(selWS) = false;

                if nnOpts.misc.semiSupervisedOnlyFS
                    % Keep x% of train and all val
                    selFS = obj.batchOpts.images.isFullySupervised(:) | obj.batchOpts.images.set(:) == 2;
                    obj.batchOpts.images.name = obj.batchOpts.images.name(selFS);
                    obj.batchOpts.images.set = obj.batchOpts.images.set(selFS);
                    obj.batchOpts.images.segmentation = obj.batchOpts.images.segmentation(selFS);
                    obj.batchOpts.images.isFullySupervised = obj.batchOpts.images.isFullySupervised(selFS);

                    if strStartsWith(obj.dataset.name, 'VOC')
                        obj.batchOpts.images.id = obj.batchOpts.images.id(selFS);
                        obj.batchOpts.images.classification = obj.batchOpts.images.classification(selFS);
                        obj.batchOpts.images.size = obj.batchOpts.images.size(:, selFS);
                    end
                end
            end

            % Make sure val images are always fully supervised
            obj.batchOpts.images.isFullySupervised(obj.batchOpts.images.set == 2) = true;

%             % Print overview of the fully and weakly supervised number of training
%             % images
%             fsCount = sum( obj.batchOpts.images.isFullySupervised(:) & obj.batchOpts.images.set(:) == 1);
%             wsCount = sum(~obj.batchOpts.images.isFullySupervised(:) & obj.batchOpts.images.set(:) == 1);
%             fsRatio = fsCount / (fsCount+wsCount);
%             wsRatio = 1 - fsRatio;
%             fprintf('Images in train: %d FS (%.1f%%), %d WS (%.1f%%)...\n', fsCount, fsRatio * 100, wsCount, wsRatio * 100);

            % Get training and test/validation subsets
            % We always validate and test on val
            trainInds = find(obj.batchOpts.images.set == 1 & obj.batchOpts.images.segmentation);
            valInds   = find(obj.batchOpts.images.set == 2 & obj.batchOpts.images.segmentation);
            obj.data.train = obj.batchOpts.images.name(trainInds); %#ok<FNDSB>
            obj.data.val   = obj.batchOpts.images.name(valInds); %#ok<FNDSB>
            obj.data.train = obj.data.train(:);
            obj.data.val   = obj.data.val(:);

            % Dataset-independent imdb fields
            obj.numClasses = obj.dataset.labelCount;
        end

        function[inputs, numElements] = getBatch(obj, batchIdx, net, nnOpts)
            % [inputs, numElements] = getBatch(obj, batchIdx, net, nnOpts)
            %
            % Returns a single image, its labels and boxes to be used as a sub-batch in training, validation or testing.
            % Note: Currently train and val batches are treated the same.
            %
            % Copyright by Holger Caesar, 2015

            % Check inputs
            assert(~isempty(obj.datasetMode));
            assert(numel(batchIdx) == 1);
            imageCount = numel(batchIdx);
            imageIdx = batchIdx;

            % Check settings
            assert(~isempty(obj.batchOpts.rgbMean));

            % Determine whether testing
            testMode = strcmp(obj.datasetMode, 'test');

            % Init labels
            if ~testMode
                lx = 1 : obj.batchOpts.labelStride : obj.batchOpts.imageSize(2);
                ly = 1 : obj.batchOpts.labelStride : obj.batchOpts.imageSize(1);
                labels = zeros(numel(ly), numel(lx), 1, imageCount, 'double'); % must be double for to avoid numerical precision errors in vl_nnloss, when using many classes
                if nnOpts.misc.weaklySupervised
                    labelsImageCell = cell(imageCount, 1);
                end
                if nnOpts.misc.maskThings
                    assert(isa(obj.dataset, 'EdiStuffDataset'));
                    datasetIN = ImageNetDataset();
                end
                masksThingsCell = cell(imageCount, 1); % by default this is empty
            end

            % Get image
            if true
                % Get image
                imageName = obj.data.(obj.datasetMode){imageIdx};
                imageOri = single(obj.dataset.getImage(imageName)) * 255;

                if ~testMode
                    image = zeros(obj.batchOpts.imageSize(1), obj.batchOpts.imageSize(2), 3, imageCount, 'single');

                    % Crop and rescale image
                    h = size(imageOri, 1);
                    w = size(imageOri, 2);
                    sz = obj.batchOpts.imageSize(1 : 2);
                    scale = max(h / sz(1), w / sz(2));
                    scale = scale .* (1 + (rand(1) - .5) / 5);
                    sy = round(scale * ((1:sz(1)) - sz(1)/2) + h/2);
                    sx = round(scale * ((1:sz(2)) - sz(2)/2) + w/2);

                    % Flip image
                    if obj.batchOpts.imageFlipping && rand > 0.5
                        sx = fliplr(sx);
                    end

                    % Get image indices in valid area
                    okx = find(1 <= sx & sx <= w);
                    oky = find(1 <= sy & sy <= h);
                    image(oky, okx, :, 1) = imageOri(sy(oky), sx(okx), :);
                else
                    image = imageOri;

                    % Workaround: Limit image size to avoid running out of RAM
                    maxImSize = 700;
                    if ~isempty(maxImSize),
                        maxSize = max(size(image, 1), size(image, 2));
                        if maxSize > maxImSize,
                            factor = maxImSize / maxSize;
                            image = imresize(image, factor);
                        end;
                    end;

                    % Some networks requires the image to be a multiple of 32 pixels
                    imageNeedsToBeMultiple = true;
                    if imageNeedsToBeMultiple
                        sz = [size(image, 1), size(image, 2)];
                        sz = round(sz / 32) * 32;
                        image = imresize(image, sz);
                    end
                end

                % Subtract mean image
                image = bsxfun(@minus, image, obj.batchOpts.rgbMean);
            end

            % Get labels
            if true
                % Fully supervised: Get pixel level labels
                if ~testMode
                    % Get pixel-level GT
                    if obj.dataset.annotation.hasPixelLabels || obj.batchOpts.images.isFullySupervised(imageIdx)
                        anno = uint16(obj.batchOpts.imageNameToLabelMap(imageName));

                        % Translate labels s.t. 255 is mapped to 0
                        if obj.batchOpts.translateLabels,
                            % Before: 255 = ignore, 0 = bkg, 1:n = classes
                            % After : 0 = ignore, 1 = bkg, 2:n+1 = classes
                            anno = mod(anno + 1, 256);
                        end
                        % 0 = ignore, 1:n = classes
                    else
                        anno = [];
                    end

                    if ~isempty(anno)
                        tlabels = zeros(sz(1), sz(2), 'double');
                        tlabels(oky,okx) = anno(sy(oky), sx(okx));
                        tlabels = single(tlabels(ly, lx));
                        labels(:, :, 1, 1) = tlabels; % 0: ignore
                    end
                end

                if ~testMode
                    % Weakly supervised: extract image-level labels
                    if nnOpts.misc.weaklySupervised
                        if ~isempty(anno) && ~all(anno(:) == 0)
                            % Get image labels from pixel labels
                            % These are already translated (if necessary)
                            curLabelsImage = unique(anno);
                        else
                            curLabelsImage = obj.dataset.getImLabelInds(imageName);

                            % Translate labels s.t. 255 is mapped to 0
                            if obj.batchOpts.translateLabels
                                curLabelsImage = mod(curLabelsImage + 1, 256);
                            end

                            if obj.dataset.annotation.labelOneIsBg
                                % Add background label
                                curLabelsImage = unique([0; curLabelsImage(:)]);
                            end
                        end

                        % Remove invalid pixels
                        curLabelsImage(curLabelsImage == 0) = [];

                        % Store image-level labels
                        labelsImageCell{1} = single(curLabelsImage(:));
                    end
                end

                % Optional: Mask out thing pixels
                if ~testMode
                    if nnOpts.misc.maskThings
                        % Get mask
                        longName = datasetIN.shortNameToLongName(imageName);
                        mask = datasetIN.getImLabelBoxesMask(longName);

                        % Resize it if necessary
                        if size(mask, 1) ~= size(image, 1) || ...
                                size(mask, 2) ~= size(image, 2)
                            mask = imresize(mask, [size(image, 1), size(image, 2)]);
                        end
                        masksThingsCell{1} = mask;
                    end
                end
            end

            % Extract inverse class frequencies from dataset
            if ~testMode
                if obj.batchOpts.useInvFreqWeights,
                    if nnOpts.misc.weaklySupervised,
                        classWeights = obj.dataset.getLabelImFreqs('train');
                    else
                        classWeights = obj.dataset.getLabelPixelFreqs('train');
                    end

                    % Inv freq and normalize classWeights
                    classWeights = classWeights ./ sum(classWeights);
                    nonEmpty = classWeights ~= 0;
                    classWeights(nonEmpty) = 1 ./ classWeights(nonEmpty);
                    classWeights = classWeights ./ sum(classWeights);
                    assert(~any(isnan(classWeights)));
                else
                    classWeights = [];
                end
                obj.batchOpts.classWeights = classWeights;
            end

            % Move image to GPU
            if strcmp(net.device, 'gpu')
                image = gpuArray(image);
            end

            % Store in output struct
            inputs = {'input', image};
            if ~testMode
                %TODO: fix imageIdx relative to subset
                if obj.dataset.annotation.hasPixelLabels || obj.batchOpts.images.isFullySupervised(imageIdx)
                    inputs = [inputs, {'label', double(labels)}];  % has to be double to avoid problems in loss
                end
                if nnOpts.misc.weaklySupervised
                    assert(~any(cellfun(@(x) isempty(x), labelsImageCell)));
                    inputs = [inputs, {'labelsImage', labelsImageCell}];
                end

                % Instance/pixel weights, can be left empty
                inputs = [inputs, {'classWeights', obj.batchOpts.classWeights}];

                % Decide which level of supervision to pick
                if nnOpts.misc.semiSupervised
                    % SS
                    isWeaklySupervised = ~obj.batchOpts.images.isFullySupervised(imageIdx);
                else
                    % FS or WS
                    isWeaklySupervised = nnOpts.misc.weaklySupervised;
                end
                inputs = [inputs, {'isWeaklySupervised', isWeaklySupervised}];
                inputs = [inputs, {'masksThingsCell', masksThingsCell}];
            end
            numElements = 1; % One image
        end

        function[allBatchInds] = getAllBatchInds(obj)
            % Obtain the indices and ordering of all batches (for this epoch)

            batchCount = numel(obj.data.(obj.datasetMode));
            if strcmp(obj.datasetMode, 'train')
                allBatchInds = randperm(batchCount)';
            elseif strcmp(obj.datasetMode, 'val'),
                allBatchInds = (1:batchCount)';
            elseif strcmp(obj.datasetMode, 'test'),
                allBatchInds = (1:batchCount)';
            else
                error('Error: Unknown datasetMode!');
            end
        end

        function initEpoch(obj, epoch)
            % Call default method
            initEpoch@ImdbCalvin(obj, epoch);
        end
    end
end

author: nightrome

repo: matconvnet-calvin


text:
function Population = EnvironmentalSelection(Population,T)
% The environmental selection of PeEA

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Li Li

    [Fitness,extreme] = CalFitness(Population.objs);
    extreme = unique(extreme); %eliminate the identical solutions
    [~,mm] = size(extreme);

	%% Calculate the angle between each two solutions
    Angle = acos(1-pdist2(Population.objs,Population.objs,'cosine'));
    Angle(logical(eye(length(Population)))) = inf;

    %% Angle-based tournament selection
    Remain = 1 : length(Population);
    Remain(extreme) = []; %eliminate extremes first
    while length(Remain) > (T-mm) %mins mm first
        % Identify the two solutions A and B with the minimum angle
        [sortA,rank1] = sort(Angle(Remain,Remain),2); %sortA is the sorted angle, rank1 NXN index of matric angle
        [~,rank2]     = sortrows(sortA);              %rank2 Nx1 index of sorted angle
        A = rank2(1);                                 %A, with minimum angle
        B = rank1(A,1);                               %B corrponding to A with minimum angle
        % Eliminate one of A and B
        if  Fitness(Remain(A)) > Fitness(Remain(B))
            Remain(A) = [];
        else
            Remain(B) = [];
        end
    end
    % Population for next generation
    Population = Population([Remain,extreme]); %add extremes into pop
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [model] = tapas_sem_multiv_model(data, ptheta, pars)
%% Set up the model.
%
% Input
%       hgf         -- Hgf model.
%       pars        -- pars structure.
% Output
%       model       -- Model structre.
%

% aponteeduardo@gmail.com
% copyright (C) 2016
%
%% Define the model.
%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

model = struct('graph', []);
model.graph = cell(4, 1);

for i = 1:4
    model.graph{i} = struct('llh', [], 'htheta', []);
end

model.graph{1}.llh = @tapas_sem_multiv_llh;
model.graph{2}.llh = @tapas_mdlinear_hier_llh;
model.graph{3}.llh = @tapas_mdlinear_llh;
model.graph{4}.llh = [];

% Computes the likelihood for a single node.
model.graph{2}.llh_sn = @tapas_mdlinear_hier_llh_sn;

model.graph{1}.htheta = struct('pe', 0.5, 'T', pars.T, 'model', ptheta);

model.graph{2}.htheta = struct('T', ones(size(pars.T, 2)));
model.graph{3}.htheta = struct('T', ones(size(pars.T, 2)));

% The last level is a dummy used to store the hyperpriors.

model.graph{4}.htheta = struct('y', [], 'u', []);

end

author: translationalneuromodeling

repo: tapas


text:
% reader struct
%    array: specify whether input is array or single channel
%    multiArrayFiles: if input is array, specify whether each channel is
%    one file.
%    useChannel: if input is array, specify which channel(s) to use.
%    if not defined, use all channels.
%
function feat = Reader_wavfile2spectrum(files, reader, useGPU, precision)
if nargin<4
    precision = 'single';
end
for i=1:length(files)
    [wav, fs] = Reader_waveform(files(i), reader);
    if strcmpi(precision, 'single')
        wav = single(wav{1});
    else
        wav = double(wav{1});
    end

    reader = SetDefaultValue(reader, 'frame_len', fs*0.025);
    reader = SetDefaultValue(reader, 'frame_shift', fs*0.01);
    reader = SetDefaultValue(reader, 'window_type', 'hamming');
    reader = SetDefaultValue(reader, 'removeDC', 0);
    reader = SetDefaultValue(reader, 'useGPU', 0);
    reader = SetDefaultValue(reader, 'doDithering', 0);

    FFT_length = 2^nextpow2(reader.frame_len);
    nCh = size(wav',2);
    tmp = sfft_multi(wav', reader.frame_len, reader.frame_shift, FFT_length, reader.window_type, reader.removeDC, reader.useGPU);
    % take the first half of the Fourier coefficients
    tmp = tmp(1:FFT_length/2+1,:,:);
    % reshape multi-channel Fourier coefficients into a vector
    tmp_feat = reshape(tmp, nCh*(FFT_length/2+1), size(tmp,3));

    if strcmpi(precision, 'single')
        feat{i} = single(gather(tmp_feat));
    else
        feat{i} = double(gather(tmp_feat));
    end
end
end

author: singaxiong

repo: SignalGraph


text:
function LabelAxis(ha, xlbls, ylbls)
%Puts text labels on axis of graph
%example:
%figure(1); plot ([1 2], [4 5]);
%LabelAxis(gca,{'Low' 'Hi'}, {'four' 'five'})

if (~isempty(xlbls) & ~iscellstr(xlbls)),
   fprintf(2,'Need array of cellstrings for x input');
end
if (~isempty(ylbls) & ~iscellstr(ylbls)),
   fprintf(2,'Need array of cellstrings for y input');
end

lx = length(xlbls);
ly = length(ylbls);

if (lx > 0),
   set(ha,'XTick',[1:1:lx]);
   set(ha,'XTickLabel',char(xlbls));
end
if (ly > 0),
   set(ha,'YTick',[1:1:ly]);
   set(ha,'YTickLabel',char(ylbls));
end

return;


text:
function pre_install (desc)
  dir =  (fileparts (mfilename ("fullpath")));
  putenv ("GEOPDES_INCLUDE_DIR", [ "-I" canonicalize_file_name(dir) filesep() "inst" filesep() "include"]);
endfunction

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
filename='Bridge_quad_fine';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance';'perimeter'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'PROJECTED GRADIENT'; incrementFactor = 1;
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.2;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-2;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function VisualizeHogFeatures(bdir,fly,fnum,varargin)
% function VisualizeHogFeatures(expdir,fly,fnum,'jabfile',jabfile)

[moviename,trxfilename,params,jabfile] = myparse(varargin,...
  'moviename','movie.ufmf','trxfilename','trx.mat','params',getSTParams,...
  'jabfile','');

if ~isempty(jabfile)
   J = load(jabfile,'-mat');
   moviename = J.x.file.moviefilename;
   trxfilename = J.x.file.trxfilename;
   params = J.x.stInfo;
end

moviename = fullfile(bdir,moviename);
trackfilename = fullfile(bdir,trxfilename);

fname = 'hf';
%% params
psize = params.psize;
nbins = params.nbins;
npatches_x = params.npatches_x;
npatches_y = params.npatches_y;
patchsz_x = psize*npatches_x;
patchsz_y = psize*npatches_y;
wd = params.wd;
scale = params.scale;

%% compute the bins

% this seems to be what the centers correspond to
bincenters = linspace(0,pi,nbins+1);
bincenters = bincenters(1:nbins);
dt = mean(diff(bincenters));
binedges = [bincenters(1)-dt/2,(bincenters(1:end-1)+bincenters(2:end))/2,bincenters(end)+dt/2];

%%

tracks = load(trackfilename);
tracks = tracks.trx;

[readfcn,nframes,fid,headerinfo] = get_readframe_fcn(moviename);
im1 = readfcn(fnum);

trackndx = fnum - tracks(fly).firstframe + 1;
% locy = round(tracks(fly).y(trackndx));
% locx = round(tracks(fly).x(trackndx));
% im1 = extractPatch(im1,...
%   locy,locx,tracks(fly).theta(trackndx),patchsz);
locy = tracks(fly).y(trackndx);
locx = tracks(fly).x(trackndx);
im1 = CropImAroundTrx(im1,...
  locx,locy,tracks(fly).theta(trackndx)-pi/2,(patchsz_x-1)/2,(patchsz_y-1)/2);

H = zeros(npatches_y,npatches_x,nbins);
firstframe = tracks(fly).firstframe;
parfor yy = 1:npatches_y
  for xx = 1:npatches_x
    for oo = 1:nbins
      pfname = fullfile(bdir,'perframe',sprintf('st_%s_%02d_%02d_%d.mat',fname,yy,xx,oo));
      q = load(pfname);
      trackndx = fnum - firstframe + 1;
      H(yy,xx,oo) = q.data{fly}(trackndx);

    end
  end
end

% plot
hfig = figure();%'Visible','off');
clf;
hax = axes('Position',[0,0,1,1]);
set(hfig,'Units','pixels','Position',get(0,'ScreenSize'));

im1curr = im1;

him = imshow(imresize(im1curr,scale));
axis image;
truesize;
colormap gray;
hold on;
axis off;

colors = hsv(nbins);
colors = colors([ (end/2+1):end 1:end/2],:);

[nr,nc,~] = size(im1);
% maxv2 = max(H(:));
maxv2 = 0.03;

hogpatch = [wd wd -wd -wd wd;-psize psize psize -psize -psize]/2;
h = [];
for xi = 1:ceil(nc/psize),
  cx = (psize/2 + (xi-1)*psize)*scale+ 1 ;
  if cx+psize/2 > nc*scale,
    break;
  end
  for yi = 1:ceil(nr/psize),
    cy = (psize/2 + (yi-1)*psize)*scale+ 1 ;
    if cy+psize/2 > nr*scale,
      break;
    end

    for bini = 1:nbins,
      tmp = bincenters(bini);
      curpatch = [cos(tmp) -sin(tmp); sin(tmp) cos(tmp)]*hogpatch;
      xcurr = cx + curpatch(1,:)*scale;
      ycurr = cy + curpatch(2,:)*scale;
      h(yi,xi,bini) = patch(xcurr,ycurr,colors(bini,:),'LineStyle','none','FaceAlpha',min(1,H(yi,xi,bini)/maxv2));
    end

  end
end
truesize(hfig);
% im = getframe(hax);

author: kristinbranson

repo: JAABA


text:
function temp_obj_eval = my_compute_precision_recall_fmeasure(u, u0_GT, u0_SKL_GT)
% 121227: Reza FARRAHI MOGHADDAM and Hossein ZIAEI NAFCHI
% 090519: Reza FARRAHI MOGHADDAM
% used in objective_evaluation_core.m

%
[xm ym] = size(u);
% figure, imshow(u)

%
if (numel(u0_GT) == 0)
    u0_GT = NaN * ones([xm ym]);
end

%
if (nargin == 2)
    % u0_SKL_GT = NaN * ones([xm ym]);
	u0_SKL_GT = ~bwmorph(~u0_GT, 'thin', 'inf');
end

% TP pixels
temp_tp = [u == 0] & [u0_GT == 0];

% FP pixels
temp_fp = [u == 0] & [u0_GT ~= 0];

% FN pixels
temp_fn = [u ~= 0] & [u0_GT == 0];

% TN pixels
temp_tn = [u ~= 0] & [u0_GT ~= 0];

% SKL TP / FN pixels
temp_skl_tp = [u == 0] & [u0_SKL_GT == 0];
temp_skl_fp = [u == 0] & [u0_SKL_GT ~= 0];
temp_skl_fn = [u ~= 0] & [u0_SKL_GT == 0];
temp_skl_tn = [u ~= 0] & [u0_SKL_GT ~= 0];

% counts
count_tp = sum(sum(temp_tp));
count_fp = sum(sum(temp_fp));
count_fn = sum(sum(temp_fn));
count_tn = sum(sum(temp_tn));
count_skl_tp = sum(sum(temp_skl_tp));
count_skl_fp = sum(sum(temp_skl_fp));
count_skl_fn = sum(sum(temp_skl_fn));
count_skl_tn = sum(sum(temp_skl_tn));

% precision
temp_p = count_tp / (count_fp + count_tp);

% recall
temp_r = count_tp / (count_fn + count_tp);
if (temp_r == 0)
    temp_r = NaN;
end

% p-recall
temp_pseudo_p = count_skl_tp / (count_skl_fp + count_skl_tp);
temp_pseudo_r = count_skl_tp / (count_skl_fn + count_skl_tp);
if (temp_pseudo_r == 0)
    temp_pseudo_r = NaN;
end

% f-measure
temp_f = 100 * 2 * (temp_p * temp_r) / (temp_p + temp_r);

% p-f-measure
temp_pseudo_f = 100 * 2 * (temp_p * temp_pseudo_r) / (temp_p + temp_pseudo_r);

% sensetivity
temp_sens = count_tp / (count_tp + count_fn);
if (temp_sens == 0)
    temp_sens = NaN;
end

% specificity
temp_spec = count_tn / (count_tn + count_fp);
if (temp_spec == 0)
    temp_spec = NaN;
end

% BCR: Balanced Classification Rate
temp_BCR = 0.5 * (temp_sens + temp_spec);

% AUC: Area Under the Curve
temp_AUC = 0.5 * (temp_sens + temp_spec);

% BER: Balanced Error Rate
temp_BER = 100 * (1 - temp_BCR);

% S-F-measure: harmonic mean of sensetivity and specificity
temp_s_f_measure = 100 * 2 * (temp_sens * temp_spec) / (temp_sens + temp_spec);

% Accuracy: mean of sensetivity and specificity
temp_accu = (count_tp + count_tn) / (count_tp + count_tn + count_fp + count_fn);

% gAccuracy: Geometric mean of sensetivity and specificity
temp_g_accu = sqrt(temp_sens * temp_spec);

% NRM (Negative Rate Metric) (*10^-2)
NR_FN = count_fn / (count_fn + count_tp);
NR_FP = count_fp / (count_fp + count_tn);
temp_NRM = (NR_FN + NR_FP) / 2;

% PSNR
err=sum(sum(temp_fp | temp_fn)) / (xm * ym);
temp_PSNR = 10 * log10( 1 / err);

% DRD: Distance Reciprocal Distortion Metric
blkSize=8; % even number
MaskSize=5; % odd number
u0_GT1 = false(xm + 2, ym + 2);
u0_GT1(2 : xm + 1, 2 : ym + 1) = u0_GT;
intim = cumsum(cumsum(u0_GT1, 1), 2);
NUBN = 0; blkSizeSQR = blkSize ^ 2;
for i= 2 : blkSize : (xm - blkSize + 1)
    for j = 2 : blkSize : (ym - blkSize + 1)
        blkSum=intim(i + blkSize - 1, j + blkSize - 1) - intim(i - 1, j + blkSize - 1) - intim(i + blkSize - 1, j - 1) + intim(i - 1, j - 1);
        if blkSum == 0 || blkSum == blkSizeSQR
        else
            NUBN = NUBN + 1;
        end
    end
end

wm = zeros(MaskSize, MaskSize);
ic = (MaskSize + 1) / 2; jc = ic; % center coordinate
for i = 1 : MaskSize
    for j = 1 : MaskSize
        wm(i, j) = 1 / (sqrt((i - ic) .^ 2 + (j - jc) .^ 2));
    end
end
wm(ic, jc) = 0;
wnm = wm ./ (sum(wm(:))); % Normalized weight matrix

u0_GT_Resized = zeros(xm + ic + 1, ym + jc + 1);
u0_GT_Resized(ic : xm + ic - 1, jc : ym + jc - 1) = u0_GT;
u_Resized = zeros(xm + ic + 1, ym + jc + 1);
u_Resized(ic : xm + ic - 1, jc : ym + jc - 1) = u;
temp_fp_Resized = [u_Resized == 0] & [u0_GT_Resized ~= 0];
temp_fn_Resized = [u_Resized ~= 0] & [u0_GT_Resized == 0];
Diff = temp_fp_Resized | temp_fn_Resized;
[xm2 ym2] = size(Diff);
SumDRDk = 0;
for i = ic : xm2 - ic + 1
    for j = jc : ym2 - jc + 1
        if Diff(i, j) == 1
            Local_Diff = my_xor_infile(u0_GT_Resized(i - ic + 1 : i + ic -1 , j - ic + 1 : j + ic - 1), u_Resized(i, j));
            DRDk = sum(sum(Local_Diff.* wnm));
            SumDRDk = SumDRDk + DRDk;
        end
    end
end
temp_DRD = SumDRDk / NUBN;

% MPM: Misclassification penalty metric
Contour = bwmorph(~u0_GT, 'remove');
Dist = bwdist(Contour, 'Chessboard');
D = sum(Dist(:));
MP_FN = sum(Dist(temp_fn)) / D;
MP_FP = sum(Dist(temp_fp)) / D;
temp_MPM = (MP_FN + MP_FP) / 2;

% output
%
temp_obj_eval.Precision = temp_p;
temp_obj_eval.Recall = temp_r;
temp_obj_eval.Fmeasure = temp_f;
temp_obj_eval.P_Precision = temp_pseudo_p;
temp_obj_eval.P_Recall = temp_pseudo_r;
temp_obj_eval.P_Fmeasure = temp_pseudo_f;
temp_obj_eval.Sensitivity = temp_sens;
temp_obj_eval.Specificity = temp_spec;
temp_obj_eval.BCR = temp_BCR;
temp_obj_eval.AUC = temp_AUC;
temp_obj_eval.BER = temp_BER;
temp_obj_eval.SFmeasure = temp_s_f_measure;
temp_obj_eval.Accuracy = temp_accu;
temp_obj_eval.GAccuracy = temp_g_accu;
temp_obj_eval.NRM = temp_NRM;
temp_obj_eval.PSNR = temp_PSNR;
temp_obj_eval.DRD = temp_DRD;
temp_obj_eval.MPM = temp_MPM;

end

function temp_xor_infile = my_xor_infile(u_infile, u0_GT_infile)
% Reza

temp_fp_infile = [u_infile == 0] & [u0_GT_infile ~= 0];
temp_fn_infile = [u_infile ~= 0] & [u0_GT_infile == 0];
temp_xor_infile = temp_fp_infile | temp_fn_infile;

end

author: Sable

repo: mcbench-benchmarks


text:
function param = gpReversibleDynamicsExtractParam(model)

% GPREVERSIBLEDYNAMICSEXTRACTPARAM Extract parameters from the GP reversible dynamics model.

% FGPLVM

param = gpDynamicsExtractParam(model);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function view = installTalairachCoordinates(view);
%loadTalairachCoordinates() Load Tal. coordinates from file and convert the
%to Volume ROIs
%The Gray (Volume) window must be open.
%
% Authors RFD and JR 09/18/01
% 2002.01.17 RFD modified initialization to allow GRAY view as well as
% volume.

global mrSESSION

if ~strcmp(view.viewType,'Volume') & ~strcmp(view.viewType,'Gray');

    error('view must be VOLUME or GRAY')

end

transform = loadTalairachXform(mrSESSION.subject);

talDat = loadTalairachCoordinates;

xyz=[talDat{1}',talDat{2}',talDat{3}'];

roiNames = cellstr(talDat{4});

roiColors = cellstr(talDat{5});

% Let the user select from the available coords

[idx,ok] = listdlg('PromptString','Select Talairach coordinates',...

       'SelectionMode','multiple',...

       'ListString',roiNames);

if ~ok %user has aborted selection

    return

end

xyz=xyz(idx,:);

roiNames=roiNames(idx);

roiColors=roiColors(idx);

vCoords = talairachToVol(xyz, transform.vol2Tal);

for(ii=[1:size(vCoords,1)])

    view = newROI(view, roiNames{ii}, 1, roiColors{ii});

    view.ROIs(view.selectedROI).coords = round(vCoords(ii,:))';

%     roi.name = roiNames{ii};%['Talairach_',num2str(tc(ii,1)),'_',num2str(tc(ii,2)),'_',num2str(tc(ii,3)),''];

%     roi.coords = round(vCoords(ii,:))';

%     roi.color = roiColors{ii};

%     roi.viewType = view.viewType;

%     view = addROI(view,roi);

end


text:
function [keypointsAll, imgidxs_missing]= det2keypointsImgMulticutBBoxMultiPoses(expidx,parts,nImgs)

fprintf('det2keypointsImgMulticutBBoxMultiPoses()\n');

p = rcnn_exp_params(expidx);
load([p.detDir '/annolist_poses'],'annolist_poses');

load(p.testGT,'annolist');

assert(length(annolist) == length(annolist_poses));

nmissing = 0;
bVis = false;
imgidxs_missing = [];

if (isfield(p,'nPeopleMax'))
    nPeopleMax = p.nPeopleMax;
else
    nPeopleMax = inf;
end

if (isfield(p,'nPeopleMin'))
    nPeopleMin = p.nPeopleMin;
else
    nPeopleMin = 0;
end

nTotalBBox = 0;
nTotalBBoxDet = 0;
nTotalImgs = 0;
nTotalImgsNotDet = 0;
idxmap = [8 7 6 12 13 14 -1 -1 2 1 5 4 3 9 10 11];
for imgidx = 1:nImgs

    fprintf('.');
    keypointsAll(imgidx).imgname = annolist_poses(imgidx).image.name;

    ndet = 0; clusidx = 0;
    bNeedNumberOfPeople = false;
    if (length(annolist(imgidx).annorect) <= nPeopleMax && ...
        length(annolist(imgidx).annorect) >= nPeopleMin)
        nTotalImgs = nTotalImgs + 1;
        bNeedNumberOfPeople = true;
        for ridx = 1:length(annolist_poses(imgidx).annorect)
            nTotalBBox = nTotalBBox + 1;
            rect = annolist_poses(imgidx).annorect(ridx);
            if (isfield(rect,'joints') && ~isempty(rect.joints))
                nTotalBBoxDet = nTotalBBoxDet + 1;
                clusidx = clusidx + 1;
                joints = rect.joints;
                %             img = imread(annolist_poses(imgidx).image.name);
                %             figure(100); clf; imagesc(img); axis equal; hold on;
                %             labels = cell(0);
                %             for i = 1:14
                %                 labels{i} = i;
                %             end
                %             joints = double(joints);
                %             plot(joints(:,1),joints(:,2),'r+','MarkerSize',5);
                %             text(joints(:,1)+10,joints(:,2),labels,'FontSize',6,'BackgroundColor','w');
                for i = 1:length(p.pidxs)
                    pidx = p.pidxs(i);
                    assert(parts(pidx+1).pos(1) == parts(pidx+1).pos(2));
                    jidx = parts(pidx+1).pos(1);
                    ndet = ndet + 1;
                    keypointsAll(imgidx).det{clusidx}{jidx+1} = joints(idxmap(jidx+1),:);
                end
            end
        end
        if (ndet == 0)
            nTotalImgsNotDet = nTotalImgsNotDet + 1;
        end
    end
    if (ndet == 0)
        for i = 1:length(p.pidxs)
            pidx = p.pidxs(i);
            jidx = parts(pidx+1).pos(1);
            keypointsAll(imgidx).det{1}{jidx+1} = [];
        end
        if (~bNeedNumberOfPeople)
            imgidxs_missing = [imgidxs_missing; imgidx];
            nmissing = nmissing + 1;
        end
        if (~mod(imgidx, 100))
            fprintf(' %d/%d\n',imgidx,length(annolist_poses));
        end
        continue;
    end

    if (bVis)
        colors = {'r','g','b','c','m','y'};
        markers = {'+','o','s'}; %,'x','.','-','s','d'
        img = imread(annolist_poses(imgidx).image.name);
        figure(100); clf; imagesc(img); axis equal; hold on;
        for j = 1:length(keypointsAll(imgidx).det)
            labels = cell(0);
            n = 0;
            points = zeros(1,2);
            for i = 1:length(p.pidxs)
                pidx = p.pidxs(i);
                jidx = parts(pidx+1).pos(1);
                pp = keypointsAll(imgidx).det{j}{jidx+1};
                if (~isempty(pp))
                    n = n  + 1;
                    labels{n} = num2str(i);
%                     labels{n} = sprintf('%1.2f',pp(3));
                    points(n,:) = pp(1:2);
                end
            end
            lp = j;
            if (lp <= 5)
                m = markers{1};
            elseif (lp > 5 && lp <= 11)
                m = markers{2};
            else
                m = markers{3};
            end
            plot(points(:,1),points(:,2),[colors{mod(lp,6)+1} m],'MarkerSize',5);
            text(points(:,1)+10,points(:,2),labels,'FontSize',6,'BackgroundColor','w');
            axis off;
        end
    end

    if (~mod(imgidx, 100))
        fprintf(' %d/%d\n',imgidx,length(annolist_poses));
    end
end
fprintf(' done\n');
fprintf('nmissing: %d\n',nmissing);
fprintf('nTotal: %d\n',nTotalBBox);
fprintf('nTotalDet: %d\n',nTotalBBoxDet);
fprintf('nTotalImgs: %d\n',nTotalImgs);
fprintf('nTotalImgsNotDet: %d\n',nTotalImgsNotDet);
nTotalImgsNotDet/nTotalImgs
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
filename = 'bridge_10_2';
ptype = 'MACRO';
method = 'SIMP_P3';
materialType = 'ISOTROPIC';
initial_case = 'holes';
cost = {'volume'};
weights = [1];
constraint = {'complianceConstraintC2','complianceConstraintC2','complianceConstraintC2'};
constraint_case = {'INEQUALITY','INEQUALITY','INEQUALITY'};
optimizerUnconstrained = 'SLERP';%'PROJECTED GRADIENT';
optimizer = 'AlternatingPrimalDual';%'DualNestedInPrimal';'DualNestedInPrimal';%'AlternatingPrimalDual';%'AlternatingPrimalDual';
incrementFactor = 1;
designVariable = 'LevelSet';%'Density';
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.4;
optimality_final =1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
Perimeter_target = 5;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

% For all tests
plotting = true;
printing = true;
printing_physics = false;
monitoring = true;
monitoring_interval = 1;
maxiter = 3e4;


text:
function masks = txt2img( strings, h, txtPrp )
% Convert text string to a binary image.
%
% For each input string, creates a binary matrix of dimensions [hxw] that
% when displayed as an image shows the string. See example for usage.
% First, uses text() to display the text using Matlab's display
% functionality. Next uses getframe() to capture the screen and crops the
% resulting image appropriately.
%
% USAGE
%  masks = txt2img( strings, h, [txtPrp] )
%
% INPUTS
%  strings  - {1xn} text string or cell array of text strings to convert
%  h        - font height in pixels
%  txtPrp   - [] additional properites for call to "text"
%
% OUTPUTS
%  masks    - {1xn} binary image masks of height h for each string
%
% EXAMPLE
%  masks=txt2img('hello world',100); im(masks{1})
%
% See also CHAR2IMG, GETFRAME, TEXT
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.35
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

% set up blank figure to which to write text
if(nargin<3 || isempty(txtPrp)), txtPrp={}; end
txtPrp = [txtPrp 'units','pixels', 'position',[2 2], 'FontUnits', ...
  'pixels', 'fontsize',h-2, 'VerticalAlignment','Bottom', 'string'];
ss=get(0,'ScreenSize'); w=ss(3)-250;
hf=figure; clf; imshow(ones([h-1 w]));
% write each string and use screen capture to save image
if(~iscell(strings)), strings={strings}; end
n=length(strings); masks=cell(1,n);
for s=1:n
  ht=text(txtPrp{:},strings{s}); pos=get(ht,'Extent');
  if(pos(3)>w), error('string does not fit on screen'); end
  M=getframe(gca); M=M.cdata; M=M(:,1:pos(3),1); M=uint8(M>100);
  k=find(sum(M==0)>0); if(k), M=M(:,max(1,k(1)-1):min(end,k(end)+1)); end
  masks{s}=M; delete(ht);
end; close(hf);
end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function IQR=interquartileRange(x,dim,definition)
%%INTERQUARTILERANGE Determine the interquartile range of the values in x.
%             This is notionally the difference between the value in the
%             75th percentile (the third quartile) and the 25th percentile
%             (the first quartile) of the data. However, the strict
%             definition varies, which is why there are multiple options.
%
%INPUTS: x A vector, matrix or hypermatrix of real values. The
%          interquartile range is only taken over a particular dimension.
%      dim An optional parameter specifying the dimensions across which
%          the interquartile range is found. If this value is omitted or an
%          empty matrix is passed, then the interquartile range is computed
%          across the first non-singleton dimension of x.
% definition This selects the definition fo the interquartile range that is
%          to be used. Let n be the number of points in dimension dim of x.
%          The following definitions are considering each independent
%          vector in dimension dim being processed. Possible values are:
%          0 (The default if omitted or an empty matrix is passed) If an
%            odd number of points is passed, then the first quartile is the
%            median of ordered (ascending) points up to floor(n/2) and the
%            third quartile is the median of ordered points from ceil(n/2)
%            to n. Thus, the median is assigned to the upper half. If an
%            even number of points is passed, then the first quartile is
%            the median of the ordered points up to n/2 and the third
%            quartile is the median of ordered points n/2+1 to n. Thus, the
%            halves are split.
%          1 This is the same as 0, except when given an odd number of
%            points, the middle point is omitted, not assigned to the upper
%            quartile. If only 1 point is given, 0 is returned.
%          2 This is the same as 0 and 1 except when given an odd number of
%            points, the median point is included in both the lower and the
%            upper quartile. This type of interquartile range is also known
%            as the midhinge.
%          3 The first and third quartiles are given by percentileVal with
%            its definition=0.
%          4 The first and third quartiles are given by percentileVal with
%            its definition=1.
%          5 This uses the kthOrderStat to get the 1+floor(0.25*(n-1)) and
%            1+ceil(0.75*(n-1)) order statistics.
%          6 In the set of ordered points, let the split point be
%            k50=ceil(n/2). The lower quartile value is the one at index
%            ceil(k50/2) and the upper quartile value is the one at index
%            k50+ceil((n-k50)/2). This is identical to considering the
%            points part of an empirical distribution and choosing the
%            percentiles using EmpiricalD.invCDF. If none aling directly
%            with 1/4 and 3/4, choose the next highest point.
%
%OUTPUTS: IQR The interquartile range values. This is a matrix with the
%             same dimensions as x except dimension dim has become unitary.
%             If dimension dim is a singleton, then IRQ will be zero for
%             all definitions.
%
%EXAMPLE:
%Here, one sees some of the differences:
% for def=0:6
%     IQR=interquartileRange([1,2,3,4,5],[],def)
% end
% one will get 2.5, 3, 2, 2.5, 2, 2, and 2.
%
%September 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(isempty(x))
    IQR=[];
    return;
end

%Select the first non-singleton dimension if dim is not provided.
if(nargin<2||isempty(dim))
    dim=find(size(x)>1,1);
    if(isempty(dim))
        dim=1;
    end
end

if(nargin<3||isempty(definition))
    definition=0;
end

n=size(x,dim);

if(n==1)%If the interquartile range is a singleton, then zero is returned
        %for all methods.
    IQR=zeros(size(x));
    return;
end

switch(definition)
    case 0%The definition used in Mathematica and in Matlab's IQR function.
          %The middle point for an odd number of points is assigned to the
          %third quartile.
        idxVec=repmat({':'},1,ndims(x));

        x=sort(x,dim,'ascend');
        if(mod(n,2)==0)%Split the halves
            idxVec{dim}=1:(n/2);
            Q1=median(x(idxVec{:}),dim);

            idxVec{dim}=(n/2+1):n;
            Q3=median(x(idxVec{:}),dim);
        else%Assign the median point to the third quartile.
            idxVec{dim}=1:floor(n/2);
            Q1=median(x(idxVec{:}),dim);

            idxVec{dim}=ceil(n/2):n;
            Q3=median(x(idxVec{:}),dim);
        end
    case 1%Alternative successive median definition 1.
        idxVec=repmat({':'},1,ndims(x));

        x=sort(x,dim,'ascend');
        if(mod(n,2)==0)%Split the halves
            idxVec{dim}=1:(n/2);
            Q1=median(x(idxVec{:}),dim);

            idxVec{dim}=(n/2+1):n;
            Q3=median(x(idxVec{:}),dim);
        else%Omit the median point.
            idxVec{dim}=1:floor(n/2);
            Q1=median(x(idxVec{:}),dim);

            idxVec{dim}=(ceil(n/2)+1):n;
            Q3=median(x(idxVec{:}),dim);
        end
    case 2%Alternative successive median definition 2.
        idxVec=repmat({':'},1,ndims(x));

        if(mod(n,2)==0)%Split the halves
            idxVec{dim}=1:(n/2);
            Q1=median(x(idxVec{:}),dim);

            idxVec{dim}=(n/2+1):n;
            Q3=median(x(idxVec{:}),dim);
        else%Include the median point in both the upper and the lower
            %quartiles.
            idxVec{dim}=1:ceil(n/2);
            Q1=median(x(idxVec{:}),dim);

            idxVec{dim}=(ceil(n/2)):n;
            Q3=median(x(idxVec{:}),dim);
        end

    case 3%The percentile definition, option 0.
        Q1=percentileVal(x,25,dim,0);
        Q3=percentileVal(x,75,dim,0);
    case 4%The percentile definition, option 1.
        Q1=percentileVal(x,25,dim,1);
        Q3=percentileVal(x,75,dim,1);
    case 5%The order statistic definition.
        [Q1,x]=kthOrderStat(x,1+floor(0.25*(n-1)),dim);
        Q3=kthOrderStat(x,1+ceil(0.75*(n-1)),dim);
    case 6
        idxVec=repmat({':'},1,ndims(x));

        x=sort(x,dim,'ascend');

        k50=ceil(n/2);
        k25=ceil(k50/2);
        k75=k50+ceil((n-k50)/2);

        idxVec{dim}=k25;
        Q1=x(idxVec{:});

        idxVec{dim}=k75;
        Q3=x(idxVec{:});
    otherwise
        error('Unknown definition of the interquartile range specified.')
end

IQR=Q3-Q1;

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 kern = simKernParamInit(kern)

% SIMKERNPARAMINIT SIM kernel parameter initialisation.
% The single input motif (SIM) kernel is specifically designed for
% working with gene networks where there is assumed to be a single
% transcription factor controlling several genes. If each gene is
% related to the transcription factor through the following
% differential equation,
%
% dx(t)/dt = B + S f(t-delta) - D x(t),
%
% where D is a decay term, S is a response term, delta is a time delay
% and B is an initial level. Then if f(t) is assumed to come from a
% Gaussian process with an RBF covariance function x(t) is a Gaussian
% process with a covariance function provided by the single input
% motif kernel.
%
% The kernel is designed to interoperate with the multiple output
% block kernel so that f(t) can be inferred given several different
% instantiations of x(t) (associated with different genes).
%
% By default the parameters (B, S, delta and D) are constrained positive. If
% kern.options.isNegativeS is set true then the parameter S is allowed to go
% negative.
%
% FORMAT
% DESC initialises the single input motif kernel structure with some
% default parameters.
% ARG kern : the kernel structure which requires initialisation.
% RETURN kern : the kernel structure with the default parameters placed in.
%
% SEEALSO : kernCreate, kernParamInit, simKernCompute
%
% COPYRIGHT : Neil D. Lawrence, 2006, 2009

% KERN

if kern.inputDimension > 1
  error('SIM kernel only valid for one-D input.')
end

if isfield(kern, 'options') && isfield(kern.options, 'gaussianInitial') && ...
      kern.options.gaussianInitial,
  kern.gaussianInitial = 1;
  kern.initialVariance = 1;
else
  kern.gaussianInitial = 0;
end

kern.delay = 0;
kern.decay = 1;
kern.initVal = 1;
kern.variance = 1;
kern.inverseWidth = 1;

if kern.gaussianInitial,
  kern.nParams = 4;
else
  kern.nParams = 3;
end

if isfield(kern, 'options') ...
      && isfield(kern.options, 'isNegativeS') ...
      && kern.options.isNegativeS,
  kern.isNegativeS = true;
  kern.transforms.index = setdiff(1:kern.nParams, 3);
  kern.sensitivity = 1;
else
  kern.isNegativeS = false;
  kern.transforms.index = 1:kern.nParams;
end
kern.transforms.type = optimiDefaultConstraint('positive');

kern.isStationary = false;
kern.isNormalised = false;
kern.positiveTime = true;


text:
%% Schmid Factor Analysis
% This script describes how to analyze Schmid factors.
%
% Let us assume a Nickel crystal

CS = crystalSymmetry('cubic',[3.523,3.523,3.523],'mineral','Nickel')

%%
% Since Nickel is fcc a dominant slip system is given by the slip plane
% normal

n = Miller(1,1,1,CS,'hkl')

%%
% and the slip direction (which needs to be orthogonal)

d = Miller(0,-1,1,CS,'uvw')

%%
% For tension in direction 123
r = normalize(vector3d(1,2,3))

%%
% the Schmid factor for the slip system [0-11](111) is defined by
%
tau = dot(d,r,'noSymmetry') * dot(n,r,'noSymmetry')

%%
% The same computation can be performed by defining the slip system as an
% MTEX variable

sS = slipSystem(d,n)

%%
% and using the command <slipSystem.SchmidFactor.html SchmidFactor>

sS.SchmidFactor(r)

%%
% Ommiting the tension direction r the command
% <slipSystem.SchmidFactor.html SchmidFactor> returns the Schmid factor as
% a <S2FunHarmonic.S2FunHarmonic.html spherical function>

SF = sS.SchmidFactor

% plot the Schmid factor in dependency of the tension direction
plot(SF)

% find the tension directions with the maximum Schmid factor
[SFMax,pos] = max(SF)

% and annotate them
annotate(pos)

%% Stress Tensor
% Instead by the tension direction the stress might be specified by a
% stress tensor

sigma = stressTensor.uniaxial(vector3d.Z)

%%
% Then the Schmid factor for the slip system |sS| and the stress tensor
% |sigma| is computed by

sS.SchmidFactor(sigma)

%% Active Slip System
% In general a crystal contains not only one slip system but at least all
% symmetrically equivalent ones. Those can be computed with

sSAll = sS.symmetrise('antipodal')

%%
% The option |antipodal| indicates that Burgers vectors in oposite
% direction should not be distinguished.
% Now

tau = sSAll.SchmidFactor(r)

%%
% returns a list of Schmid factors and we can find the slip system with the
% largest Schmid factor using

[tauMax,id] = max(abs(tau))

sSAll(id)

%%
% The above computation can be easily extended to a list of tension
% directions

% define a grid of tension directions
r = plotS2Grid('resolution',0.5*degree,'upper');

% compute the Schmid factors for all slip systems and all tension
% directions
tau = sSAll.SchmidFactor(r);

% tau is a matrix with columns representing the Schmid factors for the
% different slip systems. Lets take the maximum rhowise
[tauMax,id] = max(abs(tau),[],2);

% vizualize the maximum Schmid factor
contourf(r,tauMax)
mtexColorbar

%%
% We may also plot the index of the active slip system
pcolor(r,id)

mtexColorMap black2white

%%
% and observe that within the fundamental sectors the active slip system
% remains the same. We can even visualize the the plane normal and the slip
% direction

% if we ommit the option antipodal we can distinguish
% between the oposite burger vectors
sSAll = sS.symmetrise

% take as directions the centers of the fundamental regions
r = symmetrise(CS.fundamentalSector.center,CS);

% compute the Schmid factor
tau = sSAll.SchmidFactor(r);

% here we do not need to take the absolut value since we consider both
% burger vectors +/- b
[~,id] = max(tau,[],2);

% plot active slip plane in red
hold on
quiver(r,sSAll(id).n,'LineWidth',2,'Color','r');

% plot active slip direction in green
hold on
quiver(r,sSAll(id).b.normalize,'LineWidth',2,'Color','g');
hold off

%%
% If we perform this computation in terms of spherical functions we obtain

% ommiting |r| gives us a list of 12 spherical functions
tau = sSAll.SchmidFactor

% now we take the max of the absolute value over all these functions
contourf(max(abs(tau),[],1),'upper')
mtexColorbar

%% The Schmid factor for EBSD data
% So far we have always assumed that the stress tensor is already given
% relatively to the crystal coordinate system. Next, we want to examine the
% case where the stress is given in specimen coordinates and we know the
% orientation of the crystal. Lets import some EBSD data and computet the
% grains

mtexdata csl

% take some subset
ebsd = ebsd(ebsd.inpolygon([0,0,200,50]))

grains = calcGrains(ebsd);
grains = smooth(grains,5);

plot(ebsd,ebsd.orientations,'micronbar','off')
hold on
plot(grains.boundary,'linewidth',2)
hold off

%%
% We want to consider the following slip systems

sS = slipSystem.fcc(ebsd.CS)
sS = sS.symmetrise;

%%
% Since, those slip systems are in crystal coordinates but the stress
% tensor is in specimen coordinates we either have to rotate the slip
% systems into specimen coordinates or the stress tensor into crystal
% coordinates. In the following sections we will demonstrate both ways.
% Lets start with the first one

% rotate slip systems into specimen coordinates
sSLocal = grains.meanOrientation * sS

%%
% These slip systems are now arranged in matrix form
% where the rows corrspond to the crystal reference frames of the different
% grains and the rows are the symmetrically equivalent slip systems.
% Computing the Schmid faktor we end up with a matrix of the same size

% compute Schmid factor
sigma = stressTensor.uniaxial(vector3d.X)
SF = sSLocal.SchmidFactor(sigma);

% take the maxium allong the rows
[SFMax,active] = max(SF,[],2);

% plot the maximum Schmid factor
plot(grains,SFMax,'micronbar','off','linewidth',2)
mtexColorbar location southoutside

%%
% Next we want to visualize the active slip systems.

% take the active slip system and rotate it in specimen coordinates
sSactive = grains.meanOrientation .* sS(active);

hold on
% visualize the trace of the slip plane
quiver(grains,sSactive.trace,'color','b')

% and the slip direction
quiver(grains,sSactive.b,'color','r')
hold off

%%
% We observe that the Burgers vector is in most case aligned with the
% trace. In those cases where trace and Burgers vector are not aligned the
% slip plane is not perpendicular to the surface and the Burgers vector
% sticks out of the surface.

%%
% Next we want to demonstrate the alternative route

% rotate the stress tensor into crystal coordinates
sigmaLocal = inv(grains.meanOrientation) * sigma

%%
% This becomes a list of stress tensors with respect to crystal coordinates
% - one for each grain. Now we have both the slip systems as well as the
% stress tensor in crystal coordiantes and can compute the Schmid factor

% the resulting matrix is the same as above
SF = sS.SchmidFactor(sigmaLocal);

% and hence we may proceed analogously
% take the maxium allong the rows
[SFMax,active] = max(SF,[],2);

% plot the maximum Schmid factor
plot(grains,SFMax)
mtexColorbar

% take the active slip system and rotate it in specimen coordinates
sSactive = grains.meanOrientation .* sS(active);

hold on
% visualize the trace of the slip plane
quiver(grains,sSactive.trace,'color','b')

% and the slip direction
quiver(grains,sSactive.b,'color','r')

hold off

%% Strain based analysis on the same data set

eps = strainTensor(diag([1,0,-1]))

epsCrystal = inv(grains.meanOrientation) * eps

[M, b, W] = calcTaylor(epsCrystal, sS);

plot(grains,M,'micronbar','off')
mtexColorbar southoutside

%%

[ bMax , bMaxId ] = max( b , [ ] , 2 ) ;
sSGrains = grains.meanOrientation .* sS(bMaxId) ;
hold on
quiver ( grains , sSGrains.b)
quiver ( grains , sSGrains.trace)
hold off

author: mtex-toolbox

repo: mtex

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


text:
filename='Bridge_hexahedra_fine';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'PROJECTED GRADIENT'; incrementFactor = 1;
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.05;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-2;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;


text:
%kvsupgen 'Sup-Generating Operator (Hit-Miss)'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros vsupgen.pane file
%
% Parameters:
% InputFile: i 'Input Image', required: 'input image'
% InputFile: str1 'Struct Element 1', required: 'Structuring Element'
% InputFile: str2 'Struct Element 2', required: 'Structuring Element'
% OutputFile: o 'Output Image ', required: 'resulting output image'
%
% Example: o = kvsupgen({i, str1, str2}, {'i','';'str1','';'str2','';'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% vsupgen - Sup-Generating Operator (Hit-Miss)
%
%  DESCRIPTION
% .I vsupgen
% will do the sup-generating of an input image by two set structuring elements.
%
% The output image data type will be the same as that of the highest data type
% of the input image.
%
%
%
%  EXAMPLES
% vsupgen -i ball.xv -str1 b272.str -str2 b273.str -o outimage.xv
%
% Will do the sup-generating of image "ball.xv" by the structuring set b272.str and b273.str, with the resulting image written to "outimage.xv".
%
%  "SEE ALSO"
% vcast1(1)
%
%  RESTRICTIONS
% .I vsupgen
% 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 = kvsupgen(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,..] = kvsupgen(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'str1', '__input';'str2', '__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 'vsupgen"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function calpak_test33 ( )

%*****************************************************************************80
%
%% CALPAK_TEST33 tests MONTH_LENGTH_COMMON.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 March 2013
%
%  Author:
%
%    John Burkardt
%
  n_test = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST33\n' );
  fprintf ( 1, '  For the Common calendar:\n' );
  fprintf ( 1, '  MONTH_LENGTH_COMMON returns month lengths.\n' );

  y_test(1) = 1582;
  y_test(2) = 1752;
  y_test(3) = 1900;
  y_test(4) = 2000;

  for i_test = 1 : n_test

    y = y_test(i_test);
    sy = y_to_s_common ( y );
    months = year_length_months_common ( y );
    days = year_length_common ( y );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %d\n',y );
    fprintf ( 1, '  %s\n', sy );
    fprintf ( 1, '    Year length in months = %d\n', months );
    fprintf ( 1, '    Year length in days = %d\n', days );
    fprintf ( 1, '\n' );

    for m = 1 : months
      month_name = month_to_month_name_common ( m );
      fprintf ( 1, '  %10s  %2d\n', month_name, month_length_common ( y, m ) );
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [B_hat,mapping] = train_linear_classifier(X,labels,w)
%TRAIN_LINEAR_CLASSIFIER

 [B_hat,mapping] = linear_classifier('train',X,labels,w,'classification',[]);

end

author: epfl-lasa


text:
% modulator
function [data_out]=mod_d(data_in,rate_id)
switch (rate_id)
    case 0
       M=2;  %Size of signal constellation for BFSK
    case {1,2}
       M=4;  %Size of signal constellation QPSK
    case {3,4}
       M=16;  %Size of signal constellation 16QAM
    case {5,6}
       M=64;  %Size of signal constellation 64QAM
    otherwise
       display('error in constellation modulator give proper rate_id')
end

k=log2(M); % no. of bits per symbol
x=data_in;
% bits to symbol mapping
xnew1=reshape(x,k,length(x)/k);
xsym = bi2de(xnew1.','left-msb');
% modulating % the symbole signal xsym should be
%  a column vector containing integers between 0 to M-1
switch (rate_id)
    case {0,1,2}
       y = pskmod(xsym,M);
    case {3,4,5,6}
       y = qammod(xsym,M);
    otherwise
       display('error in modulation give proper rate_id')
end
% scatterplot(y)
%% Transmitted Signal
data_out = y;


text:
function c = subsasgn(c,s,b)
%SUBSASGN     Implements subscripted assignments for long numbers
%
%   example   c(2:3) = b
%

% written  12/30/98     S.M. Rump
% modified 11/11/99     S.M. Rump  precb<=precc
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 08/26/12     S.M. Rump  global variables removed
%

  INTLAB_LONG_ERROR = getappdata(0,'INTLAB_LONG_ERROR');

  if length(s)>1
    error('multiple indexing for long assignment not allowed')
  end

  if strcmp(s.type,'()')     % assignment c(i) = b
    if isempty(b)
      c.sign(s(1).subs{:}) = [];
      c.exponent(s(1).subs{:}) = [];
      c.mantissa(s(1).subs{:},:) = [];
      if INTLAB_LONG_ERROR
        c.error.mant(s(1).subs{:}) = [];
        c.error.exp(s(1).subs{:}) = [];
      end
      return
    else
      if isa(b,'double')
        b = long(b);
      end
    end
    if isempty(c)
      c.sign = [];
      c.exponent = [];
      c.mantissa = [];
      c.error.mant = [];
      c.error.exp = [];
      c = class(c,'long');
    end
    c.sign(s(1).subs{:},1) = b.sign;
    c.exponent(s(1).subs{:},1) = b.exponent;
    precb = size(b.mantissa,2);
    precc = size(c.mantissa,2);
    if precb<=precc
      c.mantissa(s(1).subs{:},precb+1:end) = 0;
      c.mantissa(s(1).subs{:},1:precb) = b.mantissa;
    else
      c.mantissa(end,precb) = 0;
      c.mantissa(s(1).subs{:},:) = b.mantissa;
    end
    if INTLAB_LONG_ERROR
      c.error.mant(s(1).subs{:},1) = b.error.mant;
      c.error.exp(s(1).subs{:},1) = b.error.exp;
    end
  else
    error('invalid call of subsasgn')
  end

repo: OpenMAS


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
%
%
% problem 1- Unilateral Laplace Transform computation

syms s t
f=-1.25+3.5*t*exp(-2*t) +1.25*exp(-2*t);

F=laplace(f,s);
simplify (F);
pretty(ans)

%verification
f=ilaplace(F,t);
simplify(f);
pretty(ans)


text:
function images = loadPRFsize(vw, images, curSlice)
%
% load 'ecc' into images.M1 and 'pRF size' into images.M3 field (ecc is
% reloaded since it was scaled to 0-2pi at the initial loading)
%
% images = loadPRFsize(vw, images, curSlice)

global dataTYPES

if notDefined('vw'), vw=getCurView;  end
if notDefined('images'), error('images should be defined');  end
if notDefined('curSlice'), curSlice = viewGet(vw, 'Current Slice'); end

for i=1:size(dataTYPES,2)
    if strcmp(dataTYPES(i).name, 'Averages')
        vw.curDataType = i;
    end
end
vw = selectDataType(vw,vw.curDataType);
vw = loadCorAnal(vw);

% reload ecc (initially loaded ecc was scaled into the range
% between 0 and 2*pi for atlas fitting)
images.M1=vw.map{1}(:,:,curSlice);
% load pRF size
images.M3=vw.amp{1}(:,:,curSlice);

author: vistalab

repo: vistasoft


text:
function [Mean, Covar] = ecmninit(Data, InitMethod)
%ECMNINIT Calculate initial mean and covariance for ECMNMLE.
%	Initial estimates for the Mean and Covariance of Data, where Data has
%	NUMSAMPLES samples of NUMSERIES random variables with missing data.
%
%	[Mean, Covar] = ecmninit(Data);
%	[Mean, Covar] = ecmninit(Data, InitMethod);
%
% Inputs:
%	Data - NUMSAMPLES x NUMSERIES matrix with NUMSAMPLES samples of a
%		NUMSERIES-dimensional random vector. Missing values are indicated
%		by NaNs. Data is the only required argument.
%
% Optional Inputs:
%	InitMethod - String to identify one of three initialization methods to
%		compute initial estimates for the mean and covariance of the data.
%		The default method is 'nanskip'. The initialization methods are:
%		'nanskip' - (default) Skip all records with NaNs.
%		'twostage' - Estimate mean, fill NaNs with mean, then estimate covar.
%		'diagonal' - Form a diagonal covar.
%
% Outputs:
%	Mean - NUMSERIES x 1 column vector initial estimate for mean of Data.
%	Covar - NUMSERIES x NUMSERIES matrix initial estimate for covariance of
%		Data.
%
% See also ECMNMLE.

%	Author(s): R.Taylor, 4-11-2005
%	Copyright 2005 The MathWorks, Inc.
%	$Revision: $ $Date: $

% Step 1 - check arguments

if nargin < 2
	InitMethod = 'NANSKIP';
end
if nargin < 1
	error('Finance:ecmninit:MissingInputArg', ...
		'The required input argument Data is missing.');
end

if isempty(Data)
    error('Finance:ecmninit:EmptyInputData', ...
		'The required input argument Data is empty.');
end

% Step 2 - initialization

[NumSamples, NumSeries] = size(Data);

if any(sum(isnan(Data),1) == NumSamples)
	error('Finance:ecmninit:TooManyNaNs', ...
		'One or more data series has all NaN values.');
end

if sum(sum(isinf(Data)))
	error('Finance:ecmninit:InfiniteValue', ...
		'One or more infinite values found in data.');
end

InitMethod = upper(InitMethod);
if ~any(strcmp(InitMethod,{'NANSKIP','TWOSTAGE','DIAGONAL'}))
	warning('Finance:ecmninit:UnknownInitMethodString', ...
		'Unknown InitMethod string. Will use default NANSKIP.');
	InitMethod = 'NANSKIP';
end

if strcmp(InitMethod,'NANSKIP') && (NumSamples - sum(any(isnan(Data),2)) <= NumSeries)
	warning('Finance:ecmninit:ChangeInitMethod', ...
		'Unable to do NANSKIP initialization. Switching to TWOSTAGE method ...');
	InitMethod = 'TWOSTAGE';
end

% Step 3 - estimate initial mean and covariance according to selected method

if strcmp(InitMethod,'TWOSTAGE')
	Mean = nanmean(Data)';
	Covar = zeros(NumSeries,NumSeries);
	for k = 1:NumSamples
		Z = zeros(NumSeries,1);
		for i = 1:NumSeries
			if ~isnan(Data(k,i))
				Z(i) = Data(k,i) - Mean(i);
			end
		end
		Covar = Covar + Z * Z';
	end
	Covar = (1.0/NumSamples) .* Covar;
elseif strcmp(InitMethod,'DIAGONAL')
	Mean = nanmean(Data)';
	Covar = diag(nanvar(Data,1));
else									% default method is NANSKIP
	P = ~any(isnan(Data),2);
	Mean = mean(Data(P,:))';
	Covar = cov(Data(P,:),1);
end


text:
%IM_BDILATION Fixed mapping for binary dilation (DIP_Image)
%
%	B = IM_BDILATION(A,N,CONNECTIVITY,EDGE_CONDITION)
%	B = A*IM_BDILATION([],N,CONNECTIVITY,EDGE_CONDITION)
%	B = A*IM_BDILATION(N,CONNECTIVITY,EDGE_CONDITION)
%
% INPUT
%   A        Dataset with binary object images dataset (possibly multi-band)
%   N        Number of iterations (default 1)
%   CONNECTIVITY    See BDILATION
%   EDGE_CONDITION  Value of edge, default 0
%
% OUTPUT
%   B        Dataset with dilated images
%
% SEE ALSO (PRTools Guide)
% DATASETS, DATAFILES, DIP_IMAGE, BDILATION

% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Physics, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

function b = im_bdilation(varargin)

	argin = shiftargin(varargin,'scalar');
  argin = setdefaults(argin,[],1,-2,0);
  if mapping_task(argin,'definition')
    b = define_mapping(argin,'fixed');
    b = setname(b,'Image dilation');
  else
    [a,n,connect,edgecon] = deal(argin{:});
    if isa(a,'prdataset') % allows datafiles too
      isobjim(a);
      b = filtim(a,mfilename,{n,connect,edgecon});
    elseif isa(a,'double') || isa(a,'dip_image') % here we have a single image
      if checktoolbox('dipimage')
        a = dip_image(a,'bin');
        b = bdilation(a,n,connect,edgecon);
      else
        diplibwarn
        b = bwmorph(a,'dilate',n);
      end
    else
      error('Illegal input')
    end
	end

return

author: marianux

repo: ecg-kit

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


text:
function [x,lambda_opt,svd_state,Y_hat] = ridge_gcv(varargin)
%[x_hat,lambda_opt] = ridgeGCV(Y,A,P)
%
% Estimates a ridge regression model, also know as Tikhonov regularization,
% or minimum norm with L2 prior.
%
% x_hat = argmin(x) ||Y-A*x||^2 + lambda*||P*x||^2
% with lambda > 0
%
%
% [..., svd_state] = ridge_gcv(...)  returns the svd decomposition. This can
% be reused if the target vector does not change.
%
% [..., Y_hat] = ridge_gcv(...)  returns the predicted target vector.
% Residuals can be computed via E = Y-Y_hat;
%
% This code is based on a previous implementation used in Valdes-Hernandez
% et al. (2009), written by Alejandro Ojeda and Pedro Valdez-Hernandez at
% the Cuban Neuroscience Center in 2009.
%
% Author: Alejandro Ojeda, SCCN/INC/UCSD, Jul-2012
%         Tim Mullen, SCCN/INC/UCSD, Jan-2013, Apr-2013
%
% References:
%   Pedro A. Valdes-Hernandez, Alejandro Ojeda, Eduardo Martinez-Montes, Agustin
%       Lage-Castellanos, Trinidad Virues-Alba, Lourdes Valdes-Urrutia, Pedro A.
%       Valdes-Sosa, 2009. White matter architecture rather than
%       cortical surface area correlates with the EEG alpha rhythm. NeuroImage 49
%       (2010) 2328–2339

arg_define([0 Inf],varargin, ...
    arg_norep({'Y','TargetVector','y'},mandatory,[],'The target vector'), ...
    arg_norep({'A','DesignMatrix'},mandatory,[],'The design matrix. This is the data matrix (ie X).'), ...
    arg({'P','PriorInvCov'},[],[],'Prior inverse covariance (precision) matrix for params. Can be an [nc x nc] matrix, where nc is the number of columns of A. Can also be a scalar, P, specifying the prior inverse variance (precision) of each parameter (diagonal covariance matrix). If empty, identity covariance matrix assumed. A sparse matrix is advised if precision matrix is not dense.'), ...
    arg_nogui({'blksz','DesignMatrixBlockSize','designMatrixBlockSize'},[],[],'Design matrix structure. Can be a tuple [numrows numcols], in which case A consists of identical blocks of this size, along the main diagonal'), ...
    arg_subswitch({'lambdaMode','LambdaSelectionMode'},'grid_gcv', { ...
        'manual' { ...
            arg({'lambda','RegularizationParam'},1,[0 Inf],'Regularization parameter (lambda)','type','denserealdouble') ...
            }, ...
         'grid_gcv' { ...
            arg({'gridSize','GridSize'},100,[0 Inf],'Grid size for regularization param search. This is used to automatically select the regularization parameter which minimizes the Generalized Cross-Validation (GCV) criteria.') ...
            arg({'plotGCV','PlotGCV'},false,[],'Plot GCV curve'), ...
            } ...
     },'Selection mode for lambda. Automatic (GCV grid search) or Manual (must provide lambda)'), ...
    arg_nogui('svd_state',[],[],'SVD decomposition structure'), ...
    arg({'verb','Verbosity'},false,[],'Verbose output') ...
    );

[nr,nc] = size(A);

% if A is not sufficiently sparse, convert to full
if issparse(A) && nnz(A)/numel(A) > 0.9
    A = full(A);
end

if isempty(svd_state)
    if verb, fprintf('Computing SVD of design matrix.\nr'); end

    % init prior covmat
    if isscalar(P)
        P=P*speye(nc);
    end
    if ~isempty(P)
        Pinv = inverse(P);
    end
    % compute SVD
    if ~isempty(P)
        APinv = A*Pinv;
    else
        APinv = A;
    end
    if isempty(blksz)
        [U,S,V] = svd_wrapper(APinv);
    else
        %warning('block optimization not yet implemented');
        [U,S,V] = svd_wrapper(APinv);
%         [U,S,V] = svd_wrapper(APinv(1:blksz(1),1:blksz(2)));
    end
    if isempty(P)
        iPV = V;
    else
        iPV = Pinv*V;
    end
    s   = diag(S);
    s2  = s.^2;
    Ut  = U';
else
    iPV = svd_state.iPV;
    s   = svd_state.s;
    s2  = svd_state.s2;
    Ut  = svd_state.Ut;
end

UtY = Ut*Y;

switch lambdaMode.arg_selection
    case 'grid_gcv'
        % search over a grid of lambda values for the value that minimizes
        % the Generalized Cross Validation (GCV) criteria
        % lambdaOpt = argmin(lambda) { GCV(lambda) }

        % automatically determine lambda range based on singular values
        tol     = max([nr nc])*eps(max(s));
        lgrid   = logspace(log10(tol),log10(max(s)),lambdaMode.gridSize);
        gcv     = zeros(lambdaMode.gridSize,1);
        for it=1:lambdaMode.gridSize
            % compute GCV criteria
            d       = lgrid(it)./(s2+lgrid(it));
            f       = diag(d)*UtY;
            gcv(it) = dot(f,f,1)/sum(d)^2;
        end
        loc = getMinima(gcv);
        if isempty(loc),
            % no minimum found, search for elbow instead
            if verb
                fprintf('no GCV minimum, finding elbow...\nr');
            end
            [val loc] = hlp_findElbow(gcv);  % min(gcv)
        end
        loc         = loc(end);
        lambda_opt  = lgrid(loc);
        if verb
            fprintf('lambda: %0.5g, GCV: %0.5g\nr',lambda_opt,gcv(loc));
        end
        if lambdaMode.plotGCV
            plotLambdaGCV(lgrid,gcv,lambda_opt,loc);
        end
    case 'manual'
        lambda_opt = lambdaMode.lambda;
    otherwise
        error('SIFT:ridge_gcv:BadLambdaRule', ...
              'Unknown lambda learning rule %s',lambdaMode.arg_selection);
end

% solve system for parameters x
x = iPV*bsxfun(@times,(s./(s2+lambda_opt^2)),UtY);

if nargout > 2
    Y_hat = A*x;
end
if nargout > 3
    svd_state.iPV = iPV;
    svd_state.s   = s;
    svd_state.s2  = s2;
    svd_state.UtY = UtY;
end

% Helper functions
% -------------------------------------------------------------------------
function [U,S,V] = svd_wrapper(A)
% compute SVD
if issparse(A)
    [U,S,V] = svds(A,min(size(A)));
else
    [U,S,V] = svd(A,'econ');
end

function plotLambdaGCV(lgrid,gcv,lambda_opt,loc)
% plot lambda versus GCV
figure;
semilogx(lgrid,gcv)
xlabel('log-lambda');
ylabel('GCV');
hold on;
plot(lambda_opt,gcv(loc),'rx','linewidth',2);
hold off; grid on;

function indmin = getMinima(x)
% get minimum of a function
fminor = diff(x)>=0;
fminor = ~fminor(1:end-1, :) & fminor(2:end, :);
fminor = [0; fminor; 0];
indmin = find(fminor);

repo: REST


text:
function [predictedCluster,centroids,eigens] = kernelkmeans(Kn, k)

    [H, ~] = eigs(Kn, k);
    H_normalized = H ./ repmat(sqrt(sum(H.^2, 2)), 1, k);
    predictedCluster = kmeans(H_normalized, k);

    eigens = H;
    centroids = zeros(k,k);

    for i = 1:k
        centroids(i,:) = sum(H_normalized(predictedCluster==i,:),1)/size(H_normalized(predictedCluster==i,:),1);
    end
end


text:
function save(qual,outfilename)
% Saves a Seg_Quality object to a file.  The file type is
% determined by the extension.
% Inputs:
%   qual: The Seg_Quality object to be saved.
%   outfilename: The name for the output file.

assert(nargin==2)
assert(isa(qual,'Seg_Quality'))
assert(isa(outfilename,'char'))
assert(qual.validate())

dotpos = strfind(outfilename,'.');
assert(~isempty(dotpos))
extension = outfilename(dotpos(end)+1:end);
assert(~isempty(extension))

if strcmp(extension,'hdf5') || strcmp(extension,'h5')
    qual.save_hdf5(outfilename);
elseif strcmp(extension,'mat')
    qual.save_mat(outfilename);
elseif strcmp(extension,'txt')
    qual.save_txt(outfilename);
else
    error('Unknown extension "%s"\n',extension)
end

author: nesl


text:
function res = smoothPolyline(poly, M)
%SMOOTHPOLYLINE Smooth a polyline using local averaging.
%
%   RES = smoothPolygon(POLY, M)
%   POLY contains the polyline vertices, and M is the size of smoothing
%   (given as the length of the convolution window).
%   Extremities of the polyline are smoothed with reduced window (last and
%   first vertices are kept identical, second and penultimate vertices are
%   smoothed with 3 values, etc.).
%
%   Example
%     img = imread('circles.png');
%     img = imfill(img, 'holes');
%     contours = bwboundaries(img');
%     poly = contours{1}(201:500,:);
%     figure; drawPolyline(poly, 'b'); hold on;
%     poly2 = smoothPolyline(poly, 21);
%     drawPolygon(poly2, 'm');
%
%   See also
%     polygons2d, smoothPolygon, simplifyPolyline, resamplePolyline

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2015-02-17, using Matlab 8.4.0.150421 (R2014b)
% Copyright 2015-2022 INRA - Cepia Software Platform

% compute the number of elements before and after
M1 = floor((M - 1) / 2);
M2 = ceil((M - 1) / 2);

% create convolution vector
v2 = ones(M, 1) / M;

% apply filtering on central part of the polyline
res(:,1) = conv(poly(:,1), v2, 'same');
res(:,2) = conv(poly(:,2), v2, 'same');

% need to recompute the extremities
for i = 1:M1
    i2 = 2 * i - 1;
    res(i, 1) = mean(poly(1:i2, 1));
    res(i, 2) = mean(poly(1:i2, 2));
end
for i = 1:M2
    i2 = 2 * i - 1;
    res(end - i + 1, 1) = mean(poly(end-i2+1:end, 1));
    res(end - i + 1, 2) = mean(poly(end-i2+1:end, 2));
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [f] = psi(z)
%Psi     Psi (or Digamma) function valid in the entire complex plane.
%
%                 d
%        Psi(z) = --log(Gamma(z))
%                 dz
%
%usage: [f] = psi(z)
%
%tested under versions 6.0 and 5.3.1
%
%        Z may be complex and of any size.
%
%        This program uses the analytical derivative of the
%        Log of an excellent Lanczos series approximation
%        for the Gamma function.
%
%References: C. Lanczos, SIAM JNA  1, 1964. pp. 86-96
%            Y. Luke, "The Special ... approximations", 1969 pp. 29-31
%            Y. Luke, "Algorithms ... functions", 1977
%            J. Spouge,  SIAM JNA 31, 1994. pp. 931
%            W. Press,  "Numerical Recipes"
%            S. Chang, "Computation of special functions", 1996
%
%
%see also:   GAMMA GAMMALN GAMMAINC PSIN
%see also:   mhelp psi
%see also:   mhelp GAMMA

%Paul Godfrey
%pgodfrey@conexant.com
%July 13, 2001
%see gamma for calculation details...

siz = size(z);
z=z(:);
zz=z;

f = 0.*z; % reserve space in advance

%reflection point
p=find(real(z)<0.5);
if ~isempty(p)
   z(p)=1-z(p);
end

%Lanczos approximation for the complex plane

g=607/128; % best results when 4<=g<=5

c = [  0.99999999999999709182;
      57.156235665862923517;
     -59.597960355475491248;
      14.136097974741747174;
      -0.49191381609762019978;
        .33994649984811888699e-4;
        .46523628927048575665e-4;
       -.98374475304879564677e-4;
        .15808870322491248884e-3;
       -.21026444172410488319e-3;
        .21743961811521264320e-3;
       -.16431810653676389022e-3;
        .84418223983852743293e-4;
       -.26190838401581408670e-4;
        .36899182659531622704e-5];

n=0;
d=0;
for k=size(c,1):-1:2
    dz=1./(z+k-2);
    dd=c(k).*dz;
    d=d+dd;
    n=n-dd.*dz;
end
d=d+c(1);
gg=z+g-0.5;
%log is accurate to about 13 digits...

f = log(gg) + (n./d - g./gg) ;

if ~isempty(p)
   f(p) = f(p)-pi*cot(pi*zz(p));
end

p=find(round(zz)==zz & real(zz)<=0 & imag(zz)==0);
if ~isempty(p)
   f(p) = Inf;
end

f=reshape(f,siz);

return

%A demo of this routine is:
clc
clear all
close all
x=-4:1/16:4.5;
y=-4:1/16:4;
[X,Y]=meshgrid(x,y);
z=X+i*Y;
f=psi(z);
p=find(abs(f)>10);
f(p)=10;

mesh(x,y,abs(f),phase(f));
view([45 10]);
rotate3d;

figure(2);
ezplot psi;
grid on;

One=psi(2)-psi(1)
EulerGamma=-psi(1)

return


text:
% pop_importegimat() - import EGI Matlab segmented file
%
% Usage:
%   >> EEG = pop_importegimat(filename, srate, latpoint0);
%
% Inputs:
%   filename    - Matlab file name
%   srate       - sampling rate
%   latpoint0   - latency in sample ms of stimulus presentation.
%                 When data files are exported using Netstation, the user specify
%                 a time range (-100 ms to 500 ms for instance). In this
%                 case, the latency of the stimulus is 100 (ms). Default is 0 (ms)
%
% Output:
%   EEG        - EEGLAB dataset structure
%
% Authors: Arnaud Delorme, CERCO/UCSD, Jan 2010

% Copyright (C) Arnaud Delorme, 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_importegimat(filename, srate, latpoint0);

    if nargin < 3, latpoint0 = 0; end;
    if nargin < 4, lookupchanfile = 0; end;
    if nargin < 1
        % ask user
        [filename, filepath] = uigetfile('*.mat', 'Choose a Matlab file from Netstation -- pop_importegimat()');
        if filename == 0 return; end;
        filename = fullfile(filepath, filename);

        promptstr    = { { 'style' 'text'       'string' 'Sampling rate (Hz)' } ...
                         { 'style' 'edit'       'string' '250' } ...
                         { 'style' 'text'       'string' 'Sample latency for stimulus (ms)' } ...
                         { 'style' 'edit'       'string' '0' } };
        geometry = { [2 1] [2 1] };
        result       = inputgui( 'geometry', geometry, 'uilist', promptstr, ...
                                 'helpcom', 'pophelp(''pop_importegimat'')', ...
                                 'title', 'Import a Matlab file from Netstation -- pop_runica()');

        if length(result) == 0 return; end;
        srate = str2num(result{1});
        latpoint0 = str2num(result{2});
        if isempty(latpoint0), latpoint0 = 0; end;
    end;

    EEG = eeg_emptyset;
    fprintf('Reading EGI Matlab file %s\n', filename);
    tmpdata = load('-mat', filename);

    % get data types
    % --------------
    allfields = fieldnames(tmpdata);
    for index = 1:length(allfields)
        allfields{index} = allfields{index}(1:findstr(allfields{index}, 'Segment')-2);
    end;
    datatypes = unique(allfields);
    datatypes(cellfun(@isempty, datatypes)) = [];

    % read all data
    % -------------
    counttrial = 1;
    EEG.srate  = srate;
    latency = (latpoint0/1000)*EEG.srate+1;
    for index = 1:length(datatypes)
        tindex = 1;
        for tindex = 1:length(allfields)
            if isfield(tmpdata, sprintf('%s_Segment%d', datatypes{index}, tindex))
                datatrial = getfield(tmpdata, sprintf('%s_Segment%d', datatypes{index}, tindex));
                if counttrial == 1
                    EEG.pnts = size(datatrial,2);
                    EEG.data = repmat(single(0), [size(datatrial,1), size(datatrial,2), 1000]);
                end;
                EEG.data(:,:,counttrial) = datatrial;

                EEG.event(counttrial).type    = datatypes{index};
                EEG.event(counttrial).latency = latency;
                EEG.event(counttrial).epoch   = counttrial;

                counttrial = counttrial+1;
                latency = latency + EEG.pnts;
            end;
        end;
    end;
    fprintf('%d trials read\n', counttrial-1);
    EEG.data(:,:,counttrial:end) = [];

    EEG.setname  = filename(1:end-4);
    EEG.nbchan   = size(EEG.data,1);
    EEG.trials   = counttrial-1;
    if latpoint0 ~= 1
        EEG.xmin     = -latpoint0/1000;
    end;
    EEG = eeg_checkset(EEG);

    % channel location
    % ----------------
    if all(EEG.data(end,1:10) == 0)
        disp('Deleting empty data reference channel (reference channel location is retained)');
        EEG.data(end,:)   = [];
        EEG.nbchan        = size(EEG.data,1);
        EEG = eeg_checkset(EEG);
    end;
    EEG = readegilocs(EEG);

    com = sprintf('EEG = pop_importegimat(''%s'', %3.2f, %3.2f, %d);', filename, srate, latpoint0, lookupchanfile);

author: PatternRecognition

repo: OpenBMI