text:
function results = traverseDirectoryStructure(rootPath, style, functionName, varargin)
%"traverseDirectoryStructure"
%   Given a path, <rootPath>, execute the function given by functionname
%   string <functionName> in rootPath and all of its subdirectories. Compile
%   the output of functionName's call in a datastructure. Must specify 'pre'
%   or 'post' in the style field for pre or post evaluation. varargin contains
%   any arguments for the function <functionname>.
%
%   JRA 11/7/03
%
% Usage: results = traverseDirectoryStructure(rootPath, style, functionName, varargin)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

oldpwd = pwd;
try
	cd(rootPath);
    rootPath = pwd;
	numArgs = nargout(functionName);

	%Establish result's struct fields but clear it of data.
	results.dir = '';
	results.data = '';
	results(1) = [];

	if strcmpi(style,'pre')
        [outputs{1:numArgs}] = feval(functionName, varargin{:});
        results(1).dir = rootPath;
        results(1).data = outputs;
	end

	dirContents = dir;
	for i = 1:length(dirContents)
        fileInfo = dirContents(i);
        if (strcmpi(fileInfo.name,'.') | strcmpi(fileInfo.name,'..'))
            continue
        end

        if fileInfo.isdir
            if nargin > 3
                output = traverseDirectoryStructure(fullfile(rootPath, fileInfo.name), style, functionName, varargin);
            else
                output = traverseDirectoryStructure(fullfile(rootPath, fileInfo.name), style, functionName);
            end
            results(end+1:end+length(output)) = output;
        end
	end

	if strcmpi(style,'post')
        [outputs{1:numArgs}] = feval(functionName,varargin{:});
        results(end+1).dir = rootPath;         %Add to end of array.
        results(end).data = outputs;           %Just end since array grew from last step.
    end
    cd(oldpwd);
catch
   cd(oldpwd)
   error(['Error in traverseDirectoryStructure: ' lasterr '.']);
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [net, info] = cnn_train_adagrad_ms(net, imdb, getBatch, varargin)
% CNN_TRAIN   Demonstrates training a CNN
%    CNN_TRAIN() is an example learner implementing stochastic gradient
%    descent with momentum to train a CNN for image classification.
%    It can be used with different datasets by providing a suitable
%    getBatch function.

opts.train = [] ;
opts.val = [] ;
opts.numEpochs = 300 ;
opts.batchSize = 256 ;
opts.useGpu = false ;
opts.learningRate = 0.001 ;
opts.continue = false ;
opts.expDir = fullfile('data','exp') ;
opts.conserveMemory = false ;
opts.sync = true ;
opts.prefetch = false ;
opts.weightDecay = 0.0005 ;
opts.errorType = 'multiclass' ;
opts.plotDiagnostics = false ;
opts.delta = 1e-8;
opts.display = 1;
opts.snapshot = 1;
opts.test_interval = 1;
opts = vl_argparse(opts, varargin) ;

if ~exist(opts.expDir, 'dir'), mkdir(opts.expDir); end
if isempty(opts.train), opts.train = find(imdb.images.set==1); end
if isempty(opts.val), opts.val = find(imdb.images.set==2); end
if isnan(opts.train), opts.train = []; end

% -------------------------------------------------------------------------
%                                                    Network initialization
% -------------------------------------------------------------------------

for i=1:numel(net.layers)
    if ~strcmp(net.layers{i}.type,'conv'), continue; end
    net.layers{i}.filtersMomentum = zeros(size(net.layers{i}.filters), ...
        class(net.layers{i}.filters)) ;
    net.layers{i}.biasesMomentum = zeros(size(net.layers{i}.biases), ...
        class(net.layers{i}.biases)) ; %#ok<*ZEROLIKE>
    if ~isfield(net.layers{i}, 'filtersLearningRate')
        net.layers{i}.filtersLearningRate = 1 ;
    end
    if ~isfield(net.layers{i}, 'biasesLearningRate')
        net.layers{i}.biasesLearningRate = 1 ;
    end
    if ~isfield(net.layers{i}, 'filtersWeightDecay')
        net.layers{i}.filtersWeightDecay = 1 ;
    end
    if ~isfield(net.layers{i}, 'biasesWeightDecay')
        net.layers{i}.biasesWeightDecay = 1 ;
    end
end

if opts.useGpu
    net = vl_simplenn_move(net, 'gpu') ;
    for i=1:numel(net.layers)
        if ~strcmp(net.layers{i}.type,'conv'), continue; end
        net.layers{i}.filtersMomentum = gpuArray(net.layers{i}.filtersMomentum) ;
        net.layers{i}.biasesMomentum = gpuArray(net.layers{i}.biasesMomentum) ;
    end
end

G_f = cell(numel(net.layers), 1);
G_b = cell(numel(net.layers), 1);

for l=1:numel(net.layers)

    if ~strcmp(net.layers{l}.type, 'conv'), continue ; end

    G_f{l} = zeros(size(net.layers{l}.filters), 'single');
    G_b{l} = zeros(size(net.layers{l}.biases), 'single');

end

% -------------------------------------------------------------------------
%                                                        Train and validate
% -------------------------------------------------------------------------

rng(0) ;

if opts.useGpu
    one = gpuArray(single(1)) ;
else
    one = single(1) ;
end

info.train.objective = [] ;
info.train.error = [] ;
info.train.topFiveError = [] ;
info.train.speed = [] ;
info.val.objective = [] ;
info.val.error = [] ;
info.val.topFiveError = [] ;
info.val.speed = [] ;

lr = opts.learningRate ;
res = [] ;
for epoch=1:opts.numEpochs

    % fast-forward to where we stopped
    modelPath = @(ep) fullfile(opts.expDir, sprintf('net-epoch-%d.mat', ep));
    modelFigPath = fullfile(opts.expDir, 'net-train.pdf') ;
    if opts.continue
        if exist(modelPath(epoch),'file')
            if epoch == opts.numEpochs
                load(modelPath(epoch), 'net', 'info') ;
            end
            continue ;
        end
        if epoch > 1
            fprintf('resuming by loading epoch %d\n', epoch-1) ;
            load(modelPath(epoch-1), 'net', 'info') ;
        end
    end

    train = opts.train(randperm(numel(opts.train))) ;
    val = opts.val ;

    info.train.objective(end+1) = 0 ;
    info.train.error(end+1) = 0 ;
    info.train.topFiveError(end+1) = 0 ;
    info.train.speed(end+1) = 0 ;
    info.val.objective(end+1) = 0 ;
    info.val.error(end+1) = 0 ;
    info.val.topFiveError(end+1) = 0 ;
    info.val.speed(end+1) = 0 ;

    for t=1:opts.batchSize:numel(train)
        % get next image batch and labels
        batch = train(t:min(t+opts.batchSize-1, numel(train))) ;
        batch_time = tic ;
        fprintf('training: epoch %02d: processing batch %3d of %3d ...', epoch, ...
            fix(t/opts.batchSize)+1, ceil(numel(train)/opts.batchSize)) ;
        [im_ori, labels_ori] = getBatch(imdb, batch) ;
        if opts.prefetch
            nextBatch = train(t+opts.batchSize:min(t+2*opts.batchSize-1, numel(train))) ;
            getBatch(imdb, nextBatch) ;
        end

        for ss = 1 : numel(imdb.scales)

            [im, labels] = rescale_im(im_ori, labels_ori, imdb.scales(ss),...
                imdb.mask, imdb.half_size, imdb.meanPixel);

            if opts.useGpu
                im = gpuArray(im) ;
            end

            % backprop
            net.layers{end}.class = labels ;
            res = vl_simplenn(net, im, one, res, ...
                'conserveMemory', opts.conserveMemory, ...
                'sync', opts.sync) ;

            % gradient step
            for l=1:numel(net.layers)
                if ~strcmp(net.layers{l}.type, 'conv'), continue ; end

                g_f = (net.layers{l}.filtersLearningRate) * ...
                    (opts.weightDecay * net.layers{l}.filtersWeightDecay) * net.layers{l}.filters + ...
                    (net.layers{l}.filtersLearningRate) / numel(batch) * res(l).dzdw{1};
                g_b = (net.layers{l}.biasesLearningRate) * ...
                    (opts.weightDecay * net.layers{l}.biasesWeightDecay) * net.layers{l}.biases + ...
                    (net.layers{l}.biasesLearningRate) / numel(batch) * res(l).dzdw{2};

                G_f{l} = G_f{l} + g_f .^ 2;
                G_b{l} = G_b{l} + g_b .^ 2;

                net.layers{l}.filters = net.layers{l}.filters - lr ./ (opts.delta + sqrt(G_f{l})) .* g_f;
                net.layers{l}.biases  = net.layers{l}.biases - lr ./ (opts.delta + sqrt(G_b{l})) .* g_b;
            end
        end

        % print information
        batch_time = toc(batch_time) ;
        speed = numel(batch)/batch_time ;
        info.train = updateError(opts, info.train, net, res, batch_time) ;

        fprintf(' %.2f s (%.1f images/s)', batch_time, speed) ;
        n = t + numel(batch) - 1 ;
        switch opts.errorType
            case 'multiclass'
                fprintf(' err %.1f err5 %.1f', ...
                    info.train.error(end)/n*100, info.train.topFiveError(end)/n*100) ;
                fprintf('\n') ;
            case 'binary'
                fprintf(' err %.1f', ...
                    info.train.error(end)/n*100) ;
                fprintf('\n') ;
            case 'euclideanloss'
                fprintf(' err %.1f', info.train.error(end) / n);
                fprintf('\n') ;
        end

        % debug info
        if opts.plotDiagnostics
            figure(2) ; vl_simplenn_diagnose(net,res) ; drawnow ;
        end
    end % next batch

    % evaluation on validation set
    if epoch == 1 || rem(epoch, opts.test_interval) == 0 || epoch == opts.numEpochs
        for t=1:opts.batchSize:numel(val)
            batch_time = tic ;
            batch = val(t:min(t+opts.batchSize-1, numel(val))) ;
            fprintf('validation: epoch %02d: processing batch %3d of %3d ...', epoch, ...
                fix(t/opts.batchSize)+1, ceil(numel(val)/opts.batchSize)) ;
            [im, labels] = getBatch(imdb, batch) ;
            if opts.prefetch
                nextBatch = val(t+opts.batchSize:min(t+2*opts.batchSize-1, numel(val))) ;
                getBatch(imdb, nextBatch) ;
            end
            if opts.useGpu
                im = gpuArray(im) ;
            end

            net.layers{end}.class = labels ;
            res = vl_simplenn(net, im, [], res, ...
                'disableDropout', true, ...
                'conserveMemory', opts.conserveMemory, ...
                'sync', opts.sync) ;

            % print information
            batch_time = toc(batch_time) ;
            speed = numel(batch)/batch_time ;
            info.val = updateError(opts, info.val, net, res, batch_time) ;

            fprintf(' %.2f s (%.1f images/s)', batch_time, speed) ;
            n = t + numel(batch) - 1 ;
            switch opts.errorType
                case 'multiclass'
                    fprintf(' err %.1f err5 %.1f', ...
                        info.val.error(end)/n*100, info.val.topFiveError(end)/n*100) ;
                    fprintf('\n') ;
                case 'binary'
                    fprintf(' err %.1f', ...
                        info.val.error(end)/n*100) ;
                    fprintf('\n') ;
                case 'euclideanloss'
                    fprintf(' err %.1f', info.val.error(end) / n);
                    fprintf('\n') ;
            end
        end
    end

    % save
    info.train.objective(end) = info.train.objective(end) / numel(train) ;
    info.train.error(end) = info.train.error(end) / numel(train)  ;
    info.train.topFiveError(end) = info.train.topFiveError(end) / numel(train) ;
    info.train.speed(end) = numel(train) / info.train.speed(end) ;
    info.val.objective(end) = info.val.objective(end) / numel(val) ;
    info.val.error(end) = info.val.error(end) / numel(val) ;
    info.val.topFiveError(end) = info.val.topFiveError(end) / numel(val) ;
    info.val.speed(end) = numel(val) / info.val.speed(end) ;
    if epoch == 1 || rem(epoch, opts.snapshot) == 0 || epoch == opts.numEpochs
        save(modelPath(epoch), 'net', 'info') ;
    end

    if epoch == 1 || rem(epoch, opts.display) == 0 || epoch == opts.numEpochs
        figure(1) ; clf ;
        subplot(1,2,1) ;
        semilogy(1:epoch, info.train.objective, 'k') ; hold on ;
        semilogy([1 opts.test_interval : opts.test_interval : epoch epoch], info.val.objective([1 opts.test_interval : opts.test_interval : epoch epoch]), 'b') ;
        xlabel('training epoch') ; ylabel('energy') ;
        grid on ;
        h=legend('train', 'val') ;
        set(h,'color','none');
        title('objective') ;
        subplot(1,2,2) ;
        switch opts.errorType
            case 'multiclass'
                plot(1:epoch, info.train.error, 'k') ; hold on ;
                plot(1:epoch, info.train.topFiveError, 'k--') ;
                plot([1 opts.test_interval : opts.test_interval : epoch epoch], info.val.error([1 opts.test_interval : opts.test_interval : epoch epoch]), 'b') ;
                plot([1 opts.test_interval : opts.test_interval : epoch epoch], info.val.topFiveError([1 opts.test_interval : opts.test_interval : epoch epoch]), 'b--') ;
                h=legend('train','train-5','val','val-5') ;
            case 'binary'
                plot(1:epoch, info.train.error, 'k') ; hold on ;
                plot([1 opts.test_interval : opts.test_interval : epoch epoch], info.val.error([1 opts.test_interval : opts.test_interval : epoch epoch]), 'b') ;
                h=legend('train','val') ;
            case 'euclideanloss'
                plot(1 : epoch, info.train.error, 'k'); hold on;
                plot([1 opts.test_interval : opts.test_interval : epoch epoch], info.val.error([1 opts.test_interval : opts.test_interval : epoch epoch]), 'b') ;
                h = legend('train', 'val') ;
        end
        grid on ;
        xlabel('training epoch') ; ylabel('error') ;
        set(h,'color','none') ;
        title('error') ;
        drawnow ;
        print(1, modelFigPath, '-dpdf') ;
    end

end
end

% -------------------------------------------------------------------------
function info = updateError(opts, info, net, res, speed)
% -------------------------------------------------------------------------
predictions = gather(res(end-1).x) ;
sz = size(predictions) ;
n = prod(sz(1:2)) ;

labels = net.layers{end}.class ;
info.objective(end) = info.objective(end) + sum(double(gather(res(end).x))) ;
info.speed(end) = info.speed(end) + speed ;
switch opts.errorType
    case 'multiclass'
        [~,predictions] = sort(predictions, 3, 'descend') ;
        error = ~bsxfun(@eq, predictions, reshape(labels, 1, 1, 1, [])) ;
        info.error(end) = info.error(end) +....
            sum(sum(sum(error(:,:,1,:))))/n ;
        info.topFiveError(end) = info.topFiveError(end) + ...
            sum(sum(sum(min(error(:,:,1:5,:),[],3))))/n ;
    case 'binary'

        labels = labels(:,:,1,:);
        [~,predictions] = sort(predictions, 3, 'descend') ;

        predictions = predictions(:,:,1,:);
        error = ~bsxfun(@eq, predictions, labels) ;
        info.error(end) = info.error(end) + sum(error(:))/n ;
    case 'euclideanloss'
        error = euclideanloss(sigmoid(predictions), labels);
        info.error(end) = info.error(end) + error;
end
end

function [im, labels] = rescale_im(im_ori, label_ori, scale, mask, half_size, meanPixel)

mask = imresize(mask, scale, 'nearest');

for i = 1:size(im_ori,4)
    im_ii = imresize(im_ori(:,:,:,i), scale, 'nearest');
    im_large = padarray(im_ii, [half_size, half_size], 'symmetric');
    im_ii = bsxfun(@minus, im_large, meanPixel);

    label_ii = imresize(label_ori(:,:,:,i), scale, 'nearest');

    im(:, :, :, i) = im_ii;
    labels(:, :, 1, i) = label_ii;
    labels(:, :, 2, i) = double(mask);

end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [imgOut,SH] = DiffuseConvolutionSH(img, falloff)
%
%
%        [imgOut,SH]=DiffuseConvolutionSH(img,falloff)
%
%
%        Input:
%           -img: an environment map in the latitude-longitude mapping
%           -falloff: a flag. If it is set 1, it means that fall-off will
%                     be taken into account
%
%        Output:
%           -imgOut: a diffuse convolved version of img
%           -SH: a [3,9] vector where spherical harmonics for img are
%           encoded
%
%     Copyright (C) 2011  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

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

%falloff compensation
if(falloff)
    img = FallOffEnvMap(img);
end

[r,c,col]=size(img);

SH = zeros(col, 9);

%projection constants
y00 = 0.282095;
y1x = 0.488603;
y2x = 1.092548;
y20 = 0.315392;
y22 = 0.546274;

%generation of directions

[X,Y] = meshgrid(1:c, 1:r);
phi   = pi * 2 * (X / c);
theta = pi * (Y / r);
sinTheta = sin(theta);

Dx = cos(phi) .* sinTheta;
Dy = cos(theta);
Dz = sin(phi) .* sinTheta;

for i=1:col
    img(:,:,i) = img(:,:,i) .* sinTheta;
end

%Environment projection on SH
for i=1:col
    %SH 0
    SH(i,1) = mean(mean(img(:,:,i) .* y00));
    %SH 1 -1 y
    SH(i,2) = mean(mean(img(:,:,i) .* Dy * y1x));
    %SH 1  0 z
    SH(i,3) = mean(mean(img(:,:,i) .* Dz * y1x));
    %SH 1  1 x
    SH(i,4) = mean(mean(img(:,:,i) .* Dx * y1x));
    %SH 2 -2 xy
    SH(i,5) = mean(mean(img(:,:,i) .* Dx .* Dy * y2x));
    %SH 2 -1 yz
    SH(i,6) = mean(mean(img(:,:,i) .* Dy .* Dz * y2x));
    %SH 2  1 xz
    SH(i,7) = mean(mean(img(:,:,i) .* Dx .* Dz * y2x));
    %SH 2  0 3z^2-1
    SH(i,8) = mean(mean(img(:,:,i) .* (3 * (Dz.^2) - 1) * y20));
    %SH 2  2 x^2-y^2
    SH(i,9) = mean(mean(img(:,:,i) .* (Dx.^2 - Dy.^2) * y22));
end

%scaling
SH = SH * pi * pi * 2;

%convolution
imgOut = EvaluationSH(SH, Dx, Dy, Dz);

end

author: banterle


text:
function b = chpsl ( ap, n, ipvt, b )

%*****************************************************************************80
%
%% CHPSL solves a complex hermitian system factored by CHPFA.
%
%  Discussion:
%
%    A division by zero may occur if CHPCO set RCOND to 0.0
%    or CHPFA set INFO nonzero.
%
%    To compute
%
%      inverse ( A ) * C
%
%    where C is a matrix with P columns
%
%      call chpfa(ap,n,ipvt,info)
%
%      if ( info == 0 )
%        do j = 1, p
%          call chpsl(ap,n,ipvt,c(1,j))
%        end do
%      end
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 May 2007
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1,
%    LC: QA214.L56.
%
%  Parameters:
%
%    Input, complex AP(N*(N+1)/2), the output from CHPFA.
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer IPVT(N), the pivot vector from CHPFA.
%
%    Input, complex B(N), the right hand side.
%
%    Output, complex B(N), the solution.
%

%
%  Loop backward applying the transformations and inverse ( D ) to B.
%
  k = n;
  ik = ( n * ( n - 1 ) ) / 2;

  while ( 0 < k )

    kk = ik + k;
%
%  1 x 1 pivot block.
%
    if ( 0 <= ipvt(k) )

      if ( k ~= 1 )

        kp = ipvt(k);

        if ( kp ~= k )
          t     = b(k);
          b(k)  = b(kp);
          b(kp) = t;
        end

        b(1:k-1) = b(1:k-1) + b(k) * ap(ik+1:ik+k-1);

      end
%
%  Apply D inverse.
%
      b(k) = b(k) / ap(kk);
      k = k - 1;
      ik = ik - k;

    else
%
%  2 x 2 pivot block.
%
      ikm1 = ik - ( k - 1 );

      if ( k ~= 2 )

        kp = abs ( ipvt(k) );

        if ( kp ~= k - 1 )
          t      = b(k-1);
          b(k-1) = b(kp);
          b(kp)  = t;
        end

        b(1:k-2) = b(1:k-2) + b(k)   * ap(ik+1:ik+k-2);
        b(1:k-2) = b(1:k-2) + b(k-1) * ap(ikm1+1:ikm1+k-2);

      end
%
%  Apply D inverse.
%
      km1k = ik + k - 1;
      kk = ik + k;
      ak = ap(kk) / conj ( ap(km1k) );
      km1km1 = ikm1 + k - 1;
      akm1 = ap(km1km1) / ap(km1k);
      bk = b(k) / conj ( ap(km1k) );
      bkm1 = b(k-1) / ap(km1k);
      denom = ak * akm1 - 1.0;
      b(k) = ( akm1 * bk - bkm1 ) / denom;
      b(k-1) = ( ak * bkm1 - bk ) / denom;
      k = k - 2;
      ik = ik - ( k + 1 ) - k;

    end

  end
%
%  Loop forward applying the transformations.
%
  k = 1;
  ik = 0;

  while ( k <= n )
%
%  1 x 1 pivot block.
%
    if ( 0 <= ipvt(k) )

      if ( k ~= 1 )

        b(k) = b(k) + conj ( ap(ik+1:ik+k-1) ) * transpose ( b(1:k-1) );
        kp = ipvt(k);

        if ( kp ~= k )
          t     = b(k);
          b(k)  = b(kp);
          b(kp) = t;
        end

      end

      ik = ik + k;
      k = k + 1;
%
%  2 x 2 pivot block.
%
    else

      if ( k ~= 1 )

        b(k) = b(k) + conj ( ap(ik+1:ik+k-1) ) * transpose ( b(1:k-1) );
        ikp1 = ik + k;
        b(k+1) = b(k+1) ...
          + conj ( ap(ikp1+1:ikp1+k-1) ) * transpose ( b(1:k-1) );
        kp = abs ( ipvt(k) );

        if ( kp ~= k )
          t     = b(k);
          b(k)  = b(kp);
          b(kp) = t;
        end

      end

      ik = ik + k + k + 1;
      k = k + 2;

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [y,ysup] = erf_rnd(x,rnd)
% input x real non-negative column vector
% rnd  -1  y = lower bound for erf(x)
%       1  y = upper bound for erf(x)
%      []  [y,ysup] inclusion of erf(x)
% rounding may be altered after leaving erf_rnd
%

% written  05/30/13     S.M. Rump
%

  x2 = hex2num('4017744f8f74e94b');     % ~5.86: erf(x)>pred(1) for x>=x2

  y = x;
  if isempty(rnd)
    ysup = x;
  end

  index = ( x<=0.5 );                   % first method
  if any(index)                         % x in [0,0.5]
    if isempty(rnd)
      [y(index),ysup(index)] = erf_rnd1(x(index),rnd,6);
    else
      y(index) = erf_rnd1(x(index),rnd,6);
    end
  end
  Index = index;                        % store finished indices

  index = ( ~Index ) & ( x<x2 );        % second method
  if any(index)                         % x in (0.5,x2)
    y_index = 1 - erfc_rnd2(x(index));
    if rnd==-1
      y(index) = y_index.inf;
    elseif rnd==1
      y(index) = y_index.sup;
    else
      y(index) = y_index.inf;
      ysup(index) = y_index.sup;
    end
  end

  index = ( x>=x2 );
  if any(index)                         % x in [x2,inf]
    if isempty(rnd)                     % inclusion [y,ysup] of erf(x)
      y(index) = hex2num('3fefffffffffffff');   % pred(1)
      ysup(index) = 1;
    elseif rnd==1                       % y upper bound for erf(x)
      y(index) = 1;
    else                                % y lower bound for erf(x)
      y(index) = hex2num('3fefffffffffffff');   % pred(1)
    end
  end

repo: OpenMAS


text:
function test08 ( )

%*****************************************************************************80
%
%% TEST08 tests UNICYCLE_INDEX, UNICYCLE_INDEX_TO_SEQUENCE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 June 2012
%
%  Author:
%
%    John Burkardt
%
  n = 6;
  test_num = 5;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST08\n' );
  fprintf ( 1, '  UNICYCLE_INDEX converts a unicycle to index form.\n' );
  fprintf ( 1, '  UNICYCLE_INDEX_TO_SEQUENCE converts an index to unicycle form.\n' );

  for test = 1 : test_num

    [ u, seed ] = unicycle_random ( n, seed );

    unicycle_print ( n, u, '  The unicycle:' );

    u_index = unicycle_index ( n, u );

    unicycle_index_print ( n, u_index, '  The index form:' );

    u2 = unicycle_index_to_sequence ( n, u_index );

    unicycle_print ( n, u2, '  The unicycle recovered:' );

  end

  return
end


text:
function F = tag(F,t)

if nargin > 1
    F = tag(lmi(F),t);
else
    F = tag(lmi(F));
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% op_phaseAlignAverages_fd.m
% Jamie Near, McGill University 2014.
%
% USAGE:
% [out,phs]=op_phaseAlignAverages_fd(in,minppm,maxppm,Npts,avg,weighting)
%
% DESCRIPTION:
% Perform time-domain spectral registration using a limited range of
% frequencies and using only phase adjustment (no frequency adjustment).
% This is rarely used.
%
% INPUTS:
% in         = Input data structure.
% minppm     = Minimum of frequency range (ppm).
% maxppm     = Maximum of frequnecy range (ppm).
% Npts       = Number of points in time domain to use for alignment.
% avg        = Align averages to the average of the averages ('y'), or the
%              first average in the series ('n');
% weighting	 = (Optional) Apply less weight to the later points of the fid?
%
% OUTPUTS:
% out        = Output following alignment of averages.
% phs        = Vector of phases (in degrees) used for alignment.

function [out,phs]=op_phaseAlignAverages_fd(in,minppm,maxppm,Npts,avg,weighting)

if ~in.flags.addedrcvrs
    error('ERROR:  I think it only makes sense to do this after you have combined the channels using op_addrcvrs.  ABORTING!!');
end

if nargin<6
    weighting='n';
end

if in.dims.subSpecs==0
    B=1;
else
    B=in.sz(in.dims.subSpecs);
end

phs=zeros(in.sz(in.dims.averages),B);
in_avg=op_averaging(in);
fids=zeros(size(in.fids));
if weighting=='y' || weighting=='Y'
    temp=op_freqrange(in_avg,minppm,maxppm);
    wgt_func=abs(temp.fids).^2;
else
    temp=op_freqrange(in_avg,minppm,maxppm);
    wgt_func=ones(size(temp.fids));
end

for m=1:B
    if avg=='y'||avg=='Y'
        base=in_avg;
        base=op_freqrange(base,minppm,maxppm);
        base=base.fids(1:Npts,m);
        begin=1;
    else
        base=op_freqrange(in,minppm,maxppm);
        base=base.fids(1:Npts,1,m);
        begin=2;
        fids(:,1,m)=in.fids(:,1,m);
    end
    for n=begin:in.sz(in.dims.averages)
        datarange=op_freqrange(in,minppm,maxppm);
        phsdiffs=(phase(base)-phase(datarange.fids(1:Npts,n,m)))*180/pi;
        phs(n,m)=mean(phsdiffs.*wgt_func(1:Npts,m))/mean(wgt_func(1:Npts,m));
        fids(:,n,m)=addphase(in.fids(:,n,m),phs(n,m));
    end
end

%re-calculate Specs using fft
specs=fftshift(ifft(fids,[],in.dims.t),in.dims.t);

%FILLING IN DATA STRUCTURE
out=in;
out.fids=fids;
out.specs=specs;

%FILLING IN THE FLAGS
out.flags=in.flags;
out.flags.writtentostruct=1;
out.flags.freqcorrected=1;

end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function Xrle=rle(XZv)
L=length(XZv);
j=1;
k=1;
i=1;
while i<2*L
    comp=1;
    for j=j:L
        if j==L
            break
        end;
         if XZv(j)==XZv(j+1)
            comp=comp+1;
        else
            break
        end;
    end;
        Xrle(k+1)=comp;
        Xrle(k)=XZv(j);
        if j==L & XZv(j-1)==XZv(j)
            break
        end;
        i=i+1;
        k=k+2;
        j=j+1;
        if j==L
            if mod(L,2)==0
            Xrle(k+1)=1;
            Xrle(k)=XZv(j);
            else
            Xrle(k+1)=1;
            Xrle(k)=XZv(j);
                       end;
             break
        end;
    end;

author: Sable

repo: mcbench-benchmarks


text:
function w = perform_vf_reorientation(v, options)

% perform_vf_reorientation - try to reorient the vf.
%
%   w = perform_vf_reorientation(v, options);
%
%   Flip the orientation of some vectors to ensure a coherent vector field.
%   v and w are (n,n,2) 2D vector field matrices.
%
%   options.method can be 'xproj', 'yproj', 'circproj', 'localproj',
%   'custproj', 'randomized', 'laplacian'.
%
%   'randomized' uses a slow simulated annealing method to optimize
%   an Ising energy.
%   'laplacian' solves a big linear system and is the best method (but slow).
%
%   For 'custproj' you have to provide an additional vector options.w.
%
%   See also: perform_vf_normalization.
%
%   Copyright (c) 2007 Gabriel Peyre

options.null = 0;
w = getoptions(options, 'w', [1 0]);
method = getoptions(options, 'method', 'xproj');

if size(v,3)~=2
    error('Works only for 2D vector fields.');
end

n = size(v,1);
p = size(v,2);

if strcmp(method, 'laplacian')

    % Try to solve the following linear system for some s(k) at each
    % location k:
    %   s(k) = 1/|V_k| * sum_{i \in V_k} s(k)*cos(theta(k)-theta(i))
    % where V_k are the neighbors of point k and theta(k) is the angle
    % theta(k)=atan2(v_y(k),v_x(k));

    % normalize
    d = sqrt(sum(v.^2,3)); d(d<eps) = 1;
    v = v ./ repmat(d,[1 1 2]);

    vx = v(:,:,1); vy = v(:,:,2);
    theta = atan2(vy,vx);

    i = []; j = []; s = [];
    sdiag = zeros(n^2,1);   % weights along the diagonal
    [Y,X] = meshgrid(1:n,1:n);
    I = X(:)+(Y(:)-1)*n;

    deplx = [1 -1 0 0];
    deply = [0 0 1 -1];
    for k=1:length(deplx)
        dx = deplx(k); dy = deply(k);
        J = find( X+dx<=n & X+dx>0 & Y+dy<=n & Y+dy>0);
        I1 = I(J);
        J1 = X(J)+dx+(Y(J)+dy-1)*n;
        i = [i; I1]; j = [j; J1];
        s = [s; cos( theta(I1)-theta(J1) )];
        sdiag(J) = sdiag(J)+1;
    end
    % add diagonal terms
    i = [i; I]; j = [j; I]; s = [s; -sdiag];
    % enforce constraint
    I = find(i>1);
    i = i(I); j = j(I); s = s(I);
    i = [i; 1]; j = [j; 1]; s = [s; 1];
    L = sparse(i,j,s);
    % solve
    y = zeros(n^2,1); y(1)=1;
    use_cg = getoptions(options, 'use_cg', 0);
    % use_cg = 0;
    if use_cg
        if not(isfield(options, 'niter_max'))
            options.niter_max = 500;
        end
        options.x = y*0+1;
        [s,err] = perform_conjugate_gradient(L'*L,L'*y,options);
    else
        s = L\y;
    end

    s = reshape(s,n,n);
    s = sign(s);
    w = v .* repmat(d, [1 1 2]) .* repmat(s, [1 1 2]);

    return;
end

if strcmp(method, 'laplacian1')
    %%%%%% OLD %%%%%%%%%
    % normalize
    d = sqrt(sum(v.^2,3)); d(d<eps) = 1;
    v = v ./ repmat(d,[1 1 2]);
    % padd with zeros
    v0 = zeros(n+2,n+2,2);
    v0(2:end-1,2:end-1,:) = v; v = v0;

    [dY,dX] = meshgrid(-1:1, -1:1); dX = dX(:); dY = dY(:);
    dX((end+1)/2) = []; dY((end+1)/2) = [];

    i = [n^2+1]; j = [n/2+(n/2-1)*n]; z = [1];
    u = zeros(n^2,1); % diagonal term
    for t=1:length(dX)
        dx = dX(t); dy = dY(t);
        % interaction
        a = sum( v(2:end-1,2:end-1,:) .* v(2+dx:end-1+dx,2+dy:end-1+dy,:), 3);
        [y,x] = meshgrid(1:n,1:n);
        I = find(x+dx>=1 & x+dx<=n & y+dy>=1 & y+dy<=n);
        x = x(I); y = y(I);
        w = a( max(1-dx,1):min(end-dx,end), max(1-dy,1):min(end-dy,end) );
        w(abs(w)<1e-3) = 1e-3;
        i = [i; x+(y-1)*n ];
        j = [j; x+dx+(y+dy-1)*n ];
        z = [ z; w(:) ];
        % add to diagonal
        u(x+(y-1)*n) = u(x+(y-1)*n) + 1;
    end
    % add diagonal term
    [y,x] = meshgrid(1:n,1:n); x = x(:); y = y(:);
    i = [i; x+(y-1)*n ];
    j = [j; x+(y-1)*n ];
    z = [ z; -u(:) ];

    A = sparse(i,j,z);
    b = zeros(n^2+1, 1); b(end)=1;
    use_cg = getoptions(options, 'use_cg', 1);
    if use_cg
        if not(isfield(options, 'niter_max'))
            options.niter_max = 500;
        end
        y = b;
        A1 = A'*A; y = A'*b;
        options.x = y*0+1;
        [s,err] = perform_conjugate_gradient(A1,y,options);
    else
        s = A\b;
    end
    s = reshape(s,n,n);
    s = sign(s);

    w = v(2:end-1,2:end-1,:);
    w = w .* repmat(d, [1 1 2]) .* repmat(s, [1 1 2]);
    return;

end

if strcmp(method, 'propagation')

    % normalize
    d = sqrt(sum(v.^2,3)); d(d<eps) = 1;
    v = v ./ repmat(d,[1 1 2]);
    % padd with zeros
    v0 = zeros(n+2,n+2,2);
    v0(2:end-1,2:end-1,:) = v; v = v0;

    neigh = {[1 0] [-1 0] [0 1] [0 -1]};
    % state
    n1 = n+2;
    S = zeros(n1);
    S(1,:) = -1; S(:,1) = -1; S(end,:) = -1; S(:,end) = -1;
    % initial point
    list = getoptions(options, 'start_points', floor(rand(2,1)*n+1));
    list = list+1;
    Prio = zeros(n1,n1)+Inf;
    for k=1:size(list,2)
        Prio(list(1,k),list(2,k))=0;
    end
    iter = 0;
    while not(isempty(list))
        iter = iter+1;
        progressbar(iter,n^2);
        % clf; imageplot(S); drawnow;
        I = list(1,:) + (list(2,:)-1)*n1;
        prio = Prio(I);
        % extract best match
        [tmp,I] = min(prio);
        i = list(1,I); j = list(2,I);
        prio(I) = []; list(:,I) = [];
        % set to dead
        S(i,j) = -1;
        % compute average vector
        tau = [S(i+1,j); S(i-1,j); S(i,j+1); S(i,j-1)]; tau = tau==-1;
        mtau = sum(tau);
        if mtau==0 && Prio(i,j)~=0 % iter>1
            error('Problem with computations.');
        end
        if mtau>0
            tau = repmat(tau, [1 1 2]);
            av = v(i+1,j,:).*tau(1,1,:) + v(i-1,j,:).*tau(2,1,:) + v(i,j+1,:).*tau(3,1,:) + v(i,j-1,:).*tau(4,1,:);
            av = av / mtau;
            if sum(av.*v(i,j,:))<0
                v(i,j,:) = -v(i,j,:);
            end
        end
        for s=1:length(neigh)
            is = i+neigh{s}(1); js = j+neigh{s}(2);
            w = squeeze(v(i,j,:));
            Prio(is,js) = min( Prio(is,js), Prio(i,j) + 1-abs( sum(w.*neigh{s}')) );
            if S(is,js)==0
                    % add to the front
                    list(:,end+1) = [is;js];
                    S(is,js) = 1; % the front
            end
        end
    end
    progressbar(iter,iter);

    w = v(2:end-1,2:end-1,:);
    w = w .* repmat(d, [1 1 2]);
    return;
end

if strcmp(method, 'localproj')
    % special case.
    for i=1:n
        for j=1:n
            m = zeros(1,1,2);
            if i>1
                m = m + v(i-1,j,:);
            end
            if j>1
                m = m + v(i,j-1,:);
            end
            if i>1 && j>1
                m = m + v(i-1,j-1,:);
            end
            s = dot( m, v(i,j,:) );
            if s<0
                v(i,j,:) = -v(i,j,:);
            end

        end
    end

    w = v;
    return;
end

if strcmp(method, 'randomized')

    % normalize
    d = sqrt(sum(v.^2,3)); d(d<eps) = 1;
    v = v ./ repmat(d,[1 1 2]);

    % padd with zeros
    v0 = zeros(n+2,n+2,2);
    v0(2:end-1,2:end-1,:) = v; v = v0;
    vx = v(:,:,1);
    vy = v(:,:,2);

    % number of iterations
    if isfield(options, 'niter_reorient')
        niter = options.niter_reorient;
    else
        niter = 2000;
    end
    % make 2 sub-grids
    [Y,X] = meshgrid(1:n,1:n);
    I = find( mod(X(:)+Y(:),2)==0 );
    [i1,j1] = ind2sub([n n], I);
    I = find( mod(X(:)+Y(:),2)==1 );
    [i2,j2] = ind2sub([n n], I);

    delta1 = 0.9; delta2 = 0;
    tlist = linspace(0,1,niter);
%    tlist = tlist.^3;

    err = [];
    for k=1:niter
        progressbar(k,niter);
        if mod(k,2)==1
            i = i1+1; j=j1+1;
        else
            i = i2+1; j=j2+1;
        end

        wx = vx( i + (j-1)*(n+2) );
        wy = vy( i + (j-1)*(n+2) );

        zx =    vx( i+1 + (j-1)*(n+2) ) + ...
            vx( i   + (j  )*(n+2) ) + ...
            vx( i-1 + (j-1)*(n+2) ) + ...
            vx( i   + (j-2)*(n+2) );
        zy =    vy( i+1 + (j-1)*(n+2) ) + ...
            vy( i   + (j  )*(n+2) ) + ...
            vy( i-1 + (j-1)*(n+2) ) + ...
            vy( i   + (j-2)*(n+2) );
        % normalize
        dd = sqrt(zx.^2+zy.^2); dd(dd<eps) = 1;
        zx = zx./dd; zy = zy./dd;

        delta = (1-tlist(k))*delta1 + tlist(k)*delta2;

        s = wx.*zx + wy.*zy;
        s = sign(s-delta);
        vx( i + (j-1)*(n+2) ) = wx.*s;
        vy( i + (j-1)*(n+2) ) = wy.*s;
        err(end+1) = sum(s<0);
    end
    v = cat(3,vx,vy);
    w = v(2:end-1,2:end-1,:);
    w = w .* repmat(d, [1 1 2]);
    return;
end

switch lower(method)
    case 'xproj'
        s = v(:,:,1);
    case 'yproj'
        s = v(:,:,2);
    case 'circproj'
        n = size(v,1);
        p = size(v,2);
        [Y,X] = meshgrid(0:p-1, 0:n-1);
        s = v(:,:,1).*X + v(:,:,2).*Y;
    case 'custproj'
        s = v(:,:,1)*w(1) + v(:,:,2)*w(2);
    otherwise
        error('Unknown method');
end

s = sign(s);
I = find(s==0); s(I) = 1;
w = v;
w(:,:,1) = v(:,:,1).*s;
w(:,:,2) = v(:,:,2).*s;

author: gpeyre

repo: matlab-toolboxes


text:
classdef rbm
% Restricted Boltzmann Machine Model object:
%----------------------------------------------------------------------------
% Initialize and train an RBM energy-based model. Supports binary, Gaussian,
% and multinomial inputs (still in development).
%
% Supports binary and noisy rectified linear units (NReLU) in the hidden layer.
%
% Supports joint modeling of binary and multinomial variables for classification.
%
% Model regularizers include L2 and L1 (via subgradients) weight decay, hidden
% unit sparsity (binary units only), and hidden unit dropout.
%----------------------------------------------------------------------------
% Dustin E. Stansbury
% stan_s_bury@berkeley.edu

properties
	class = 'rbm'; 		% GENERAL CLASS OF MODEL
	inputType = 'binary';% TYPE OF RBM ('BB','GB')
	classifier = false; % TRAIN MULTINOMIAL UNITS FOR CLASSIFICATION IN PARALLEL
	nClasses;			% # OF OUTPUT CLASSES (.classifier = true)
	nObs;				% # OF TRAINING OBSERVATIONS
	nVis;				% # OF VISIBLE UNITS (DIMENSIONS)
	nHid = 100;			% # OF HIDDEN UNITS
	rlu = 0;			% USE RECTIFIED LINEAR UNITS
	W;					% CONNECTION WEIGHTS
	dW;					% GRADIENT FOR CONN. WEIGHTS
	b;					% VISIBLE UNIT BIASES
	db;					% GRADIENT FOR VIS. BIAS
	c;					% HIDDEN UNIT BIASES
	dc;					% GRADIENT FOR HID. BIAS
	classW=0			% CLASSFIER OUTPUT WEIGHTS
	dClassW; 			% CLASS. WEIGHT GRADIENTS
	d=0;				% CLASSIFIER BIAS
	dd;					% GRADIENT FOR classifier BIASES
	log;				% ERROR AND LEARNING RATE LOGS
	aVis;				% VISIBLE LAYER ACTIVATIONS
	aHid;				% HIDDEN LAYER ACTIVATION
	pHid;				% HIDDEN LAYER PROBS
	pHid0;				% INITIAL HIDDEN LAYER PROBS
	aHid0;				% INITIAL HIDDEN LAYER ACTIVATION
	eHid;				% EXPECTATION OF HIDDEN STATE (POST TRAINING)
	pClass;				% PROBABILITIES OF CLASSIFIER UNITS (.classifier = 1)
	aClass;				% ACTIVATION OF CLASSIFIER UNITS (.classifier = 1)
	docLen;				% DOCUMENT LENGTH (MULTINOMIAL INPUTS)
	lRate = 0.1;		% DEFAULT LEARNING RATE
	batchIdx = [];		% BATCH INDICES INTO TRAINING DATA
	sampleVis = 0;		% SAMPLE THE VISIBLE UNITS
	sampleHid = 1;		% SAMPLE HIDDEN UNITS
	momentum = 0;		% DEFAULT MOMENTUM TERM
	nEpoch = 100;		% # OF FULL PASSES THROUGH TRIANING DATA
	wPenalty = 0;		% CURRENT WEIGHT PENALTY
	sparsity = 0;		% SPARSENESS FACTOR
	dropout = 0;		% HIDDEN UNIT DROPOUT
	doMask = 1;			% DROPOUT MASK
	topoMask = [];		% MASK FOR TOPOGRAPHY
	sparseGain=5;		% LEARNING RATE GAIN FOR SPARSITY
	batchSz = 100;		% # OF TRAINING POINTS PER BATCH
	nGibbs = 1;			% CONTRASTIVE DIVERGENCE (1)
	beginAnneal = Inf;	% # OF EPOCHS TO START SIMULATED ANNEALING
	beginWeightDecay=1% # OF EPOCHS TO START WEIGHT PENALTY
	verbose = 1;		% DISPLAY PROGRESS
	saveEvery = 0;		% # OF EPOCHS TO SAVE INTERMEDIATE MODELS
	displayEvery=500;	% # WEIGHT UPDATES TO VISUALIZE LEARNING
	visFun = [];		% USER-DEFINED FUNCTION ('@myFun')
	auxVars = []; 		% AUXILLARY VARIABLES, JUST IN CASE
	useGPU = 0; 		% USE CUDA, IF AVAILABLE
	gpuDevice;			% GPU DEVICE STRUCTURE
	interactive = 0;	% STORE A GLOBAL COPY FOR INTERACTIVE TRAINING
	saveFold
end % END PROPERTIES

methods

	function self = rbm(arch)
	% r = rbm(arch)
	%--------------------------------------------------------------------------
	% RBM constructor
	%--------------------------------------------------------------------------
	% INPUT:
	%  <arch>:  - a set of arguments defining the RBM architecture.
	%
	% OUTPUT:
	%     <r>:  - an RBM model object.
	%--------------------------------------------------------------------------
		if notDefined('arch')
		else
			self = self.init(arch);
		end
	end

	function print(self)
	% print()
	%--------------------------------------------------------------------------
	% Display the properties and methods for the RBM object class.
	%--------------------------------------------------------------------------
		properties(self)
		methods(self)
	end

	function self = train(self,X,targets)
	% self = train(X,[targets])
	%--------------------------------------------------------------------------
	% Train an RBM using Contrastive Divergence
	%--------------------------------------------------------------------------
	% INPUT:
	%        <X>:  - to-be modeled data |X| = [#Obs x #Vis]
	%  <targets>:  - category targets. can be either [#Obs x #Class] as 1 of
	%                K representation or [#Obs x 1], where each entry is a
	%                numerical category label. (optional)
	% OUTPUT:
	%     <self>:  - trained RBM object.
	%--------------------------------------------------------------------------
		lRate0 = self.lRate;
		if self.sparsity > 0, meanAct = 0; end
		if notDefined('targets')
			targets = 0;
		end

		self.nObs = size(X,1);

		% INITIALIZE CLASSIFIER UNITS
		if self.classifier || any(targets)
			self.classifier = true;
			[self,targets] = self.initClassifer(targets);
		end

		self.log.err = zeros(1,self.nEpoch);
		self.log.lRate = self.log.err;

		wPenalty = self.wPenalty;
		dCount = 1;
		mCount = 1;

		self.batchIdx = self.createBatches(X);

		% SEND DATA TO GPU
		% USER WILL NEED TO IMPLEMENT gpuDistribute.m
		if self.useGPU
			self = gpuDistribute(self);
			X = gpuArray(X);
			targets = gpuArray(targets);
		end

		iE = 1;
		while 1
			sumErr = 0;

			% SIMULATED ANNEALING
			if iE > self.beginAnneal
				self.lRate = max(lRate0*((iE-self.beginAnneal)^(-.25)),1e-8);
			end

			% WEIGHT DECAY
			if iE < self.beginWeightDecay
				self.wPenalty = 0;
			else
				self.wPenalty = wPenalty;
			end

			% LOOP OVER BATCHES
			for jB = 1:numel(self.batchIdx)
				batchX = X(self.batchIdx{jB},:);
				if self.classifier
					batchTargets = targets(self.batchIdx{jB},:);
				else
					batchTargets = 0;
				end

				if strcmp(self.inputType,'multinomial')
					self.docLen = sum(batchX,2);
				end

				if self.dropout > 0 % SAMPLE DROPOUT MASK ONECE PER MINIBATCH
					self.doMask = rand(size(batchX,1),self.nHid)>self.dropout;
				else
					self.doMask = 1;
				end

				% SAMPLE FROM MODEL, CALCULATE GRADIENTS, UPDATE PARAMS
				self = self.runGibbs(batchX,batchTargets);
				self = self.updateParams(batchX,batchTargets);
				sumErr = self.accumErr(batchX,sumErr);

				if ~isempty(self.visFun) & ~mod(dCount,self.displayEvery)&iE>1;
					self.auxVars.batchX = batchX;
					self.auxVars.batchTargets = targets;
					self.auxVars.error = self.log.err(1:iE-1);
					self.auxVars.lRate = self.log.lRate(1:iE-1);
					self.visLearning;
				end
				dCount = dCount+1;
				% MOVING AVERAGE OF HIDDEN ACTIVATIONS FOR SPARSITY
				if self.sparsity & ~self.rlu
					meanAct = mean(self.aHid)*0.1 + 0.9*meanAct;
				end
			end

			% SPARSITY (BINARY HIDDENS ONLY)
			if self.sparsity & ~self.rlu
				dcSparse = -self.lRate*self.sparseGain*(meanAct-self.sparsity);
				self.c = self.c + dcSparse;
			end

 			% ENDUCE SPARSITY
 			% if self.sparsity & ~self.rlu
 			% 	dcSparse = -self.lRate*self.sparseGain*(mean(self.pHid)-self.sparsity);
 			% 	self.c = self.c + dcSparse;
 			% end

			% LOG ERRORS, ETC.
			self.log.err(iE) = sumErr/self.nObs;
			self.log.lRate(iE) = self.lRate;

			if self.verbose
				self.printProgress(sumErr,iE,jB);
			end

			% SAVE?
			if iE > 1
				if ~mod(iE, self.saveEvery) & ~isempty(self.saveFold)
					r = self;
					if ~exist(r.saveFold,'dir')
						mkdir(r.saveFold);
					end
					save(fullfile(r.saveFold,sprintf('Epoch_%d.mat',iE)),'r'); clear r;
				end
			end

			% CLEAN UP IF USING INTERACTIVE TRAINING
			if iE >= self.nEpoch
				if self.interactive
					clear global r
				end
				break
			else % STORE GLOBAL COPY FOR INTERACTIVE TRAINING
				if self.interactive
					global r
					r = self;
				end
			end
			iE = iE + 1;
		end

		% PULL DATA FROM GPU, IF NEEDED
		% USER WILL NEED TO IMPLEMENT gpuGather.m
		if self.useGPU
			self = gpuGather(self);
			reset(self.gpuDevice);
			self.gpuDevice = [];
		end
		fprintf('\n');
	end

	function self = runGibbs(self,X,targets)
	% r = runGibbs(X,targets)
	%--------------------------------------------------------------------------
	% Draw MCMC samples from the current model via Gibbs sampling.
	%--------------------------------------------------------------------------
	% INPUT:
	%       <X>:  - minibatch data.
	% <targets>:  - possible categorical targets. if no categorization,
	%               <targets> should be empty ([]).
	%
	% OUTPUT:
	%    <r>:  - RBM object with updated states
	%--------------------------------------------------------------------------
		nObs = size(X,1);
		iC = 1;
		% GO UP
		self = self.hidGivVis(X,targets,1);
		% LOG INITIAL STATES FOR GRADIENT CALCULATION
		self.pHid0 = self.pHid;
		self.aHid0 = self.aHid;

		while 1
			% GO DOWN
			self = self.visGivHid(self.aHid,self.sampleVis);
			X = self.aVis;

			% GO BACK UP
			if iC >= self.nGibbs
				self = self.hidGivVis(self.aVis,targets,0);
				break
			else
				self = self.hidGivVis(self.aVis,targets,1);
			end
			iC = iC + 1;
		end
	end

	function self = hidGivVis(self,X,targets,sampleHid)
	% r = hidGivVis(X,targets,[sampleHid])
	%--------------------------------------------------------------------------
	% Update hidden unit probabilities and states conditioned on the current
	% states of the visible units.
	%--------------------------------------------------------------------------
	% INPUT:
	%         <X>:  - batch data.
	%   <targets>:  - possible target variables.
	% <sampleHid>:  - flag indicating to sample the states of the hidden units.
	%
	% OUTPUT:
	%      <r>:  - RBM object with updated hidden unit probabilities/states.
	%--------------------------------------------------------------------------
		hidBias = self.c;
		if strcmp(self.inputType,'multinomial')
			% WEIGHT BIASES BY DOCUMENT LENGTH
			hidBias = bsxfun(@times,self.docLen,hidBias);
		end

		if self.classifier
			% JOINTLY MODEL MULTINOMIAL OVER INPUT CLASSES
			pHid = self.sigmoid(bsxfun(@plus,X*self.W + targets*self.classW ,hidBias));
		else
			pHid = self.sigmoid(bsxfun(@plus,X*self.W, hidBias));
		end

		if self.rlu % RECTIFIED LINEAR UNITS
			self.aHid = max(0,pHid + randn(size(pHid)).*sqrt(self.sigmoid(pHid)));
		else
			if sampleHid
				self.aHid = single(pHid>rand(size(X,1),self.nHid));
			else
				self.aHid = pHid;
			end
		end
		self.pHid = pHid;

		% ENDUCE TOPOGRAPHY?
		if ~isempty(self.topoMask)
			for iH = 1:size(self.aHid,1)
				self.aHid(iH,:) = sum(bsxfun(@times,self.topoMask,self.aHid(iH,:)),2);
			end
		end
		% DEAL WITH DROPOUT
		self.aHid = self.aHid.*self.doMask;
	end

	function self = visGivHid(self,aHid,sampleVis)
	% r = hidGivVis(aHid,[sampleVis])
	%--------------------------------------------------------------------------
	% Update visible unit states conditioned on the current states of the hidden
	% units.
	%--------------------------------------------------------------------------
	% INPUT:
	%      <aHid>:  - current hidden unit states (activations)
	% <sampleVis>:  - flag indicating to sample the states of the visible units.
	%
	% OUTPUT:
	%      <r>:  - RBM object with updated visible unit probabilities/states.
	%--------------------------------------------------------------------------
		nObs = size(aHid,1);
		switch self.inputType
			case 'binary'
				pVis = self.sigmoid(bsxfun(@plus,aHid*self.W',self.b));
				if sampleVis
					self.aVis = pVis>rand(nObs,self.nVis);
				else
					self.aVis = pVis;
				end

			case 'gaussian'
				mu = bsxfun(@plus,aHid*self.W',self.b);

				if sampleVis
					self.aVis = self.drawNormal(mu);
				else
					self.aVis = mu;
				end

			case 'multinomial'
				pVis = self.softMax(bsxfun(@plus,aHid*self.W',self.b));
				self.aVis = zeros(size(pVis));
				for iO  = 1:nObs
					% DRAW D SEPARATE MULTINOMIALS FOR EACH INPUT
					self.aVis(iO,:) = mnrnd(self.docLen(iO),pVis(iO,:));
				end
		end

		if self.classifier
			self.pClass = self.softMax(bsxfun(@plus,self.aHid*self.classW',self.d));
			self.aClass = self.sampleClasses(self.pClass);
		end
	end

	function self = updateParams(self,X,targets);
	% self = updateParams(X,targets)
	%--------------------------------------------------------------------------
	% Update current model parameters based on the states of hidden and visible
	% units. Weight decay is applied here.
	%--------------------------------------------------------------------------
	% INPUT:
	%       <X>:  - minibatch data.
	% <targets>:  - possible categorical targets. if no categorization,
	%               <targets> should be empty ([]).
	%
	% OUTPUT:
	%    <self>:  - RBM object with updated states
	%--------------------------------------------------------------------------
		nObs = size(X,1);

		dW = (X'*self.pHid0 - self.aVis'*self.pHid)/nObs;
		self.dW = self.momentum*self.dW + ... % MOMENTUM
			          (1-self.momentum)*dW;         % NEW GRADIENT

		self.W = self.W + self.lRate*self.dW;

		% WEIGHT DECAY
		if self.wPenalty > 0 % L2
			self.W = self.W - self.lRate*self.wPenalty*self.W;
		elseif self.wPenalty < 0 % L1 SUBGRADIENT AT 0
			self.W = self.W + self.lRate*self.wPenalty*sign(self.W);
		end

		db = mean(X) - mean(self.aVis);
		dc = mean(self.pHid0) - mean(self.pHid);

		self.db = self.momentum*self.db + self.lRate*db;
		self.b = self.b + self.db;

		self.dc = self.momentum*self.dc + self.lRate*dc;
		self.c = self.c + self.dc;

		if self.classifier
			% CLASSIFIER WEIGHTS
			dClassW=(targets'*self.pHid0 - self.aClass'*self.pHid)/nObs;
			self.dClassW=self.momentum*self.dClassW+self.lRate*(dClassW-self.wPenalty*self.classW);
			self.classW = self.classW + self.dClassW;

			% CLASSIFIER BIASES
			dd = (sum(targets) - sum(self.aClass))/nObs;
			self.dd = self.momentum*self.dd + self.lRate*dd;
			self.d = self.d + self.dd;
		end
	end

	function err = accumErr(self,X,err0);
	% err = updateParams(X,err0)
	%--------------------------------------------------------------------------
	% Add reconstruction error (squared difference) for the current batch to
	% the total error.
	%--------------------------------------------------------------------------
	% INPUT:
	%       <X>:  - minibatch data.
	%    <err0>:  - current resevoir of error
	%
	% OUTPUT:
	%     <err>:  - updated error resevoir
	%--------------------------------------------------------------------------
		err = sum(sum((X-self.aVis).^2));
		err = err + err0;
	end

	function self = init(self,arch)
	% r = init(arch)
	%--------------------------------------------------------------------------
	% Initialize an rbm object based on provided architecture, <arch>.
	%--------------------------------------------------------------------------
	% INPUT:
	%    <arch>:  - a structure of architecture options. Possible fields are:
	%               .size       --> [#Visible x #Hidden] units
	%               .inputType  --> string indicating input type (i.e. 'binary',
	%                              'gaussian','multinomial')
	%               .classifyer --> flag indicating whether to train a class-
	%                               ifier model in parallel.
	%               .opts       --> additional cell array of options, defined
	%                               in argument-value pairs.
	%
	% OUTPUT:
	%     <r>:  - an initialized RBM object.
	%--------------------------------------------------------------------------
		arch = self.ensureArchitecture(arch);

		self.nVis = arch.size(1);
		self.nHid = arch.size(2);
		self.inputType = arch.inputType;
		self.classifier = arch.classifier;

		% PARSE ANY ADDITIONAL OPTIONS
		if isfield(arch,'opts');
			opts = arch.opts;
			fn = fieldnames(self);
			for iA = 1:2:numel(opts)
				if ~isstr(opts{iA})
					error('<opts> must be a cell array string-value pairs.');
				elseif sum(strcmp(fn,opts{iA}))
					self.(opts{iA})=opts{iA+1};
				end
			end
		end

		% INIT. WEIGHTS (A' LA BENGIO)
		range = sqrt(6/(2*self.nVis));
		self.W = single(2*range*(rand(self.nVis,self.nHid)-.5));

		self.dW = zeros(size(self.W),'single');
		self.b = zeros(1,self.nVis,'single');
		self.db = zeros(size(self.b),'single');
		self.c = zeros(1,self.nHid,'single');
		self.dc = zeros(size(self.c),'single');
	end

	function arch = ensureArchitecture(self,arch)
	% arch = ensureArchitecture(self,arch)
	%--------------------------------------------------------------------------
	% Preprocess the provided architecture structure <arch>.
	%--------------------------------------------------------------------------
	% INPUT:
	% <arch>:  - is either a [2 x 1] vector giving the [#vis x # hid], in which
	%            case we use the default model parameters, or it is a strucure
	%            with the fields:
	%             .size                  --> Network size; [#Vis x # Hid]
	%             .inputType (optional)  --> ['binary'] or 'gaussian'
	%             .classifier (optional) --> true or [false]
	%             .opt (optional)        --> a cell array of {'parameter',paramValue}
	%                                        of global options
	%--------------------------------------------------------------------------
		% IF arch IS GIVEN AS A VECTOR
		if ~isstruct(arch),arch.size = arch; end

		if ~isfield(arch,'size');
			error('need to supply the # of Visible and Hidden units')
		end

		if ~isfield(arch,'inputType');
			arch.inputType = self.inputType;
		end

		if ~isfield(arch,'classifier');
			arch.classifier = self.classifier;
		end

		if isfield(arch,'opts');
			opts = arch.opts;
			if ~iscell(opts) || mod(numel(opts),2)
				error('arch.opts must be a cell array of string-value pairs.')
			end
		end
	end

	function [self,targets] = initClassifer(self,targets);
	% [self,targets] = initClassifer(targets);
	%--------------------------------------------------------------------------
	% Initialize classifier units based on form of <targets>. Can also be used
	% to preprocess targets to be 1-of-K.
	%--------------------------------------------------------------------------

		% IF SUPPLIED TARGETS ARE A LIST OF LABELS
		if isvector(targets)
			targets = self.oneOfK(targets);
		end
		self.nClasses = size(targets,2);
		self.classW = single(0.1*randn(self.nClasses,self.nHid));
		self.dClassW = zeros(size(self.classW),'single');
		self.d = zeros(1,self.nClasses,'single');
		self.dd = self.d;
	end

	function batchIdx = createBatches(self,X)
	% batchIdx = createBatches(X)
	%--------------------------------------------------------------------------
	% Create minibatches based on dimensions of inputs <X>. Returns a cell
	% array with each entry giving the indices in to the rows of <X> for a
	% single batch.
	%--------------------------------------------------------------------------
		nObs = size(X,1);
		nBatches = ceil(nObs/self.batchSz);
		tmp = repmat(1:nBatches, 1, self.batchSz);
		tmp = tmp(1:nObs);
		randIdx=randperm(nObs);
		tmp = tmp(randIdx);
		for iB=1:nBatches
		    batchIdx{iB} = find(tmp==iB);
		end
	end

	function visLearning(self);
	% visLearning(self);
	%--------------------------------------------------------------------------
	% Deal with any visualiztions. Note, must set self.visFun to appropriate
	% visualization function handle.
	%--------------------------------------------------------------------------
		if ~isempty(self.visFun)
			try
				self.visFun(self);
			catch
				fprintf('\nVisualization failed...')
			end
		end
	end

	function printProgress(self,sumErr,iE,jB)
	% printProgress(sumErr,iE,jB)
	%--------------------------------------------------------------------------
	% Utility function to display std out.
	%--------------------------------------------------------------------------
			if iE > 1
			if self.log.err(iE) > self.log.err(iE-1) & iE > 1
				indStr = '(UP)    ';
			else
				indStr = '(DOWN)  ';
			end
		else
			indStr = '';
		end
		fprintf('Epoch %d/%d --> Recon. error: %f %s\r', ...
		iE,self.nEpoch,sumErr,indStr);
	end

	function E = hidExpect(self,X);
	% E = hidExpect(X);
	%--------------------------------------------------------------------------
	% Calculate hidden unit expectations for network input <X>
	%--------------------------------------------------------------------------
		switch self.inputType
		case 'multinomial'
			docLen = sum(X,2);
			E = self.sigmoid(bsxfun(@plus,X*self.W,bsxfun(@times,docLen,self.c)));
		otherwise
			if self.rlu
 				% E = max(0,bsxfun(@plus,X*self.W,self.c));
				E = max(0,X*self.W);
				E = bsxfun(@minus,E,mean(E));
				E = bsxfun(@rdivide,E,std(E));
			else
				E = self.sigmoid(bsxfun(@plus,X*self.W,self.c));
			end
		end
	end

	function samps = sample(self,X0,nSamples,nSteps)
	% samps = sample(X0,[nSamples],[nSteps])
	%--------------------------------------------------------------------------
	% Draw sample(s) from model using a Markov chain.
	%--------------------------------------------------------------------------
	% INPUT:
	%       <X0>:  - initial state from which to begin the Markov chain
	% <nSamples>:  - the number of individual samples to generate. default is
	%                to draw a single sample.
	%   <nSteps>:  - the number of transitions/steps to take in the Markov
	%                chain. default is to take 50 steps.
	%--------------------------------------------------------------------------

		self.sampleVis = 0;
		self.sampleHid = 0;
		if notDefined('nSamples');nSamples = 1;end
		if notDefined('nSteps'),nSteps = 50; end

		[nObs,nVis] = size(X0);
		samps = zeros(nObs,nVis,nSamples);

		if self.useGPU
			samps = gpuArray(single(samps));
			self = gpuDistribute(self);
		end

		for iS = 1:nSamples
			vis = X0;
			for iI = 1:nSteps
				hid = self.sigmoid(bsxfun(@plus,binornd(1,vis)* ...
									self.W,self.c));
				switch self.inputType
				case 'binary'
					vis = self.sigmoid(bsxfun(@plus,binornd(1,hid)* ...
					self.W',self.b));
				case 'gaussian'
					vis = bsxfun(@plus,hid*self.W',self.b);

				case 'multinomial'
					vis = bsxfun(@plus,hid*self.W',self.b)
					pVis = softMax(vis);
					% ENSURE PDF
					pVis = bsxfun(@rdivide,eVis,sum(eVis,2));
					vis = zeros(size(pVis));
					docLen = sum((find(pVis)),2);
					for iO  = 1:nObs
						% DRAW D SEPARATE MULTINOMIALS FOR EACH INPUT
						vis(iO,:) = mnrnd(docLen(iO),pVis(iO,:));
					end
				end
			end
			samps(:,:,iS) = vis;
		end

		if self.useGPU
			samps = gather(samps);
		end
	end

	function F = freeEnergy(self,X)
	% F = freeEnergy(X)
	%--------------------------------------------------------------------------
	% Calculate model free energy for an input <X>. Currently only available for
	% binary and gaussian inputs and binary hidden units.
	%--------------------------------------------------------------------------
	% INPUT:
	%  <X>:  - network input.
	%
	% OUTPUT:
	%  <F>:  - free energy, summing over all binary hidden units.
	%--------------------------------------------------------------------------

		nSamps = size(X,1);
		upBound = 50; loBound = -50; % FOR POSSIBLE OVERFLOW
		switch self.inputType
			case 'binary'
				H = ones(nSamps,1)*self.c + X*self.W;
				H = max(min(H,upBound),loBound);
				% SAMPLE ENERGIES
				sampE = -X*self.b' - sum(-log(1./exp(H)),2);
				F = mean(sampE);
			case 'gaussian'
				H = ones(nSamps,1)*self.c + X*self.W;
				H = max(min(H,upBound),loBound);
				% SAMPLE ENERGIES
				sampE = bsxfun(@minus,X,self.b);
				sampE = sum(sampE.^2,2)/2;
				sampE = sampE - sum(-log(1./(1+exp(H))),2);
				F = mean(sampE);
			case 'multinomial'
				error('need to implement multinomial free energy')
		end
	end

	function targets = oneOfK(self,labels)
	% targets = oneOfK(labels)
	%--------------------------------------------------------------------------
	% Create a 1-of-k representation of a set of labels, where <labels> is a
	% vector of integer class labels.
	%--------------------------------------------------------------------------
		classes = unique(labels);
		targets = zeros(numel(labels),max(classes));
		for iL = 1:numel(classes)
			targets(labels==classes(iL),iL)=1;
		end
	end

	function p = sigmoid(self,X)
	% p = sigmoid(X)
	%--------------------------------------------------------------------------
	% Sigmoid activation function
	%--------------------------------------------------------------------------
		if self.useGPU
			p = arrayfun(@(x)(1./(1 + exp(-x))),X);
		else
			p = 1./(1 + exp(-X));
		end
	end

	function p = drawNormal(self,mu);
	% p = drawNormal(mu);
	%--------------------------------------------------------------------------
	% Draw samples from a multivariate normal  with mean <mu> and identity
	% covariance.
	%--------------------------------------------------------------------------
		% ASSUMES UNIT VARIANCE OF ALL VISIBLES
		% (I.E. YOU SHOULD STANDARDIZE INPUTS)
		p = mvnrnd(mu,ones(1,self.nVis));
	end

	function c = softMax(self,X)
	% c = softMax(X)
	%--------------------------------------------------------------------------
	% Siftmax activation function for input X. Returns a vector of category
	% probabilities.
	%--------------------------------------------------------------------------
		c = bsxfun(@rdivide,exp(X),sum(exp(X),2));
	end

	function classes = sampleClasses(self,pClass)
	%--------------------------------------------------------------------------
	% classes = sampleClasses(pClass)
	%--------------------------------------------------------------------------
	% Sample class labels <classes> given a set of class probabilities, <pClass>
	%--------------------------------------------------------------------------
		[nObs,nClass]=size(pClass);
		classes = zeros(size(pClass));
		% ENSURE NORMALIZED
		pClass = bsxfun(@rdivide,pClass,sum(pClass,2));
		for iC = 1:nObs
			probs = pClass(iC,:);
			samp = cumsum(probs)>rand;
			idx = min(find(max(samp)==samp));
			classSamp = zeros(size(probs));
			classSamp(idx) = 1;
			classes(iC,:) = classSamp;
		end
	end

	function [pred,error,misClass] = classify(self,X,targets)
	% [pred,error,misClass] = classify(X,[targets])
	%--------------------------------------------------------------------------
	% Classify inputs <X> and calculate misclassification and error when comp-
	% ared to <targets>.
	%--------------------------------------------------------------------------
	% INPUT:
	%        <X>:  - a set of inputs to the network.
	%  <targets>:  - (optional). a set of target classes. Can be a 1-of-K vector
	%                or can a vector with equal length to the number of inputs
	%                where each entry is an interger class label.
	% OUTPUT:
	%     <pred>:  - predicted class label.
	%    <error>:  - the classification error, based on provided <targets>
	% <misClass>:  - the indices of the inputs that were misclassified.
	%--------------------------------------------------------------------------

		if notDefined('targets'),
			targets = [];
		else
			if isvector(targets);
				targets = self.oneOfK(targets);
			end
		end

		nObs = size(X,1);

		% ACTIVATE HIDDEN UNITS USING INPUT ONLY
		pHid = self.sigmoid(bsxfun(@plus,X*self.W,self.c));

		% CALCULATE CLASS PROBABILITY
		pClass = self.softMax(bsxfun(@plus,pHid*self.classW',self.d));
		% WINNER-TAKE-ALL CLASSIFICATION
		[~, pred] = max(pClass,[],2);
		if ~notDefined('targets') && nargout > 1
			[~,targets] = max(targets,[],2);

			% CALCULATE MISSCLASSIFICATION RATE
			misClass = find(pred~=targets);
			error = numel(misClass)/nObs;
		end
	end
end % END METHODS
end % END CLASSDEF

author: dustinstansbury

repo: medal


text:
classdef Panel < matlab.ui.container.Panel & uix.mixin.Panel
    %uix.Panel  Standard panel
    %
    %  b = uix.Panel(p1,v1,p2,v2,...) constructs a standard panel and sets
    %  parameter p1 to value v1, etc.
    %
    %  A card panel is a standard panel (uipanel) that shows one its
    %  contents and hides the others.
    %
    %  See also: uix.CardPanel, uix.BoxPanel, uipanel

    %  Copyright 2009-2016 The MathWorks, Inc.
    %  $Revision: 1480 $ $Date: 2017-02-15 16:56:13 +0100 (Wed, 15 Feb 2017) $

    methods

        function obj = Panel( varargin )
            %uix.Panel  Standard panel constructor
            %
            %  p = uix.Panel() constructs a standard panel.
            %
            %  p = uix.Panel(p1,v1,p2,v2,...) sets parameter p1 to value
            %  v1, etc.

            % Set properties
            try
                uix.set( obj, varargin{:} )
            catch e
                delete( obj )
                e.throwAsCaller()
            end

        end % constructor

    end % structors

    methods( Access = protected )

        function redraw( obj )

            % Compute positions
            bounds = hgconvertunits( ancestor( obj, 'figure' ), ...
                [0 0 1 1], 'normalized', 'pixels', obj );
            padding = obj.Padding_;
            xSizes = uix.calcPixelSizes( bounds(3), -1, 1, padding, 0 );
            ySizes = uix.calcPixelSizes( bounds(4), -1, 1, padding, 0 );
            position = [padding+1 padding+1 xSizes ySizes];

            % Redraw contents
            selection = obj.Selection_;
            if selection ~= 0
                uix.setPosition( obj.Contents_(selection), position, 'pixels' )
            end

        end % redraw

    end % template methods

end % classdef

author: goGPS-Project


text:
% FILTERSD  Solve a DENSE NLP using FilterSD
%
% THIS IS A LOW LEVEL FUNCTION - USE opti_filtersd() INSTEAD!
%
% FilterSD is a NLP solver written in FORTRAN and available from:
% https://projects.coin-or.org/filterSD
%
%   [x,fval,exitflag,stats,lambda] = filtersd(fun, grad, x0, lb, ub, nlcon, nljac, cl, cu, opts)
%
%   Input arguments:
%       fun - nonlinear function handle
%       grad - gradient of nonlinear function handle
%       x0 - initial solution guess
%       lb - decision variable lower bounds
%       ub - decision variable upper bounds
%       nlcon - nonlinear constraints handle
%       nljac - Jacobian of nonlinear constraints handle (dense matrix)
%       cl - nonlinear constraints lower bound
%       cu - nonlinear constraints upper bound
%       opts - solver options (see below)
%
%   Return arguments:
%       x - solution vector
%       fval - objective value at the solution
%       exitflag - exit status (see below)
%       stats - statistics structure
%       lambda - Lagrange multipliers at solution [bounds;nl constraints]
%
%   Option Fields (all optional):
%       maxiter - maximum major solver iterations
%       maxtime - maximum solver execution time [won't return last solution on bounded problems]
%       maxfeval - maximum function evaluations [won't return last solution on bounded problems]
%       display - solver display level [0-2]
%       rho - initial trust region radius
%       htol - tolerance allowed in sum h(x) of constraint infeasibilities
%       rgtol - tolerance allowed in reduced gradient L2 norm
%       iterfun - iteration callback of the form ifun(feval,fval,x)
%
%   Return Status:
%       0 - Successful run
%       1 - Unbounded NLP (f(x) <= fmin at htol feasible point)
%       2 - Bounds on x are inconsistent
%       3 - Local minimum of feasibility problem and h(x) > htol (nonlinear constraints are locally inconsistent)
%       4 - Initial point x has h(x) > ubd
%       5 - Iteration limit reached
%       6 - Termination with rho <= htol
%       7 - Not enough workspace memory in ws or lws
%       8 - Insufficient space for filter
%       9 and above - LCP solver problem
%       105 - User Exit
%
%   Copyright (C) 2013 Jonathan Currie (I2C2)

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function net = mlpunpak(net, w)
%MLPUNPAK Separates weights vector into weight and bias matrices.
%
%	Description
%	NET = MLPUNPAK(NET, W) takes an mlp network data structure NET and  a
%	weight vector W, and returns a network data structure identical to
%	the input network, except that the first-layer weight matrix W1, the
%	first-layer bias vector B1, the second-layer weight matrix W2 and the
%	second-layer bias vector B2 have all been set to the corresponding
%	elements of W.
%
%	See also
%	MLP, MLPPAK, MLPFWD, MLPERR, MLPBKP, MLPGRAD
%

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

% Check arguments for consistency
errstring = consist(net, 'mlp');
if ~isempty(errstring);
  error(errstring);
end

if net.nwts ~= length(w)
  error('Invalid weight vector length')
end

nin = net.nin;
nhidden = net.nhidden;
nout = net.nout;

mark1 = nin*nhidden;
net.w1 = reshape(w(1:mark1), nin, nhidden);
mark2 = mark1 + nhidden;
net.b1 = reshape(w(mark1 + 1: mark2), 1, nhidden);
mark3 = mark2 + nhidden*nout;
net.w2 = reshape(w(mark2 + 1: mark3), nhidden, nout);
mark4 = mark3 + nout;
net.b2 = reshape(w(mark3 + 1: mark4), 1, nout);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function GRIDobj2geotiff(A,file)

%GRIDobj2geotiff Exports an instance of GRIDobj to a geotiff file
%
% Syntax
%
%     GRIDobj2geotiff(DEM)
%     GRIDobj2geotiff(DEM,filename)
%
% Description
%
%     GeoTIFF is a common image file format that stores coordinates and
%     projection information to be read by most GIS software.
%     GRIDobj2geotiff writes an instance of GRIDobj to a GeoTIFF file.
%
%     GRIDobj2geotiff requires the function geotiffwrite available with
%     the Mapping Toolbox. If geotiffwrite does not exist on the search
%     path, the function will write a standard tif together with a
%     '.tfw'-file (worldfile, http://en.wikipedia.org/wiki/World_file ) to
%     the disk.
%
%     GRIDobj2geotiff(DEM) opens a dialogue box to save the GeoTIFF
%
%     GRIDobj2geotiff(DEM,filename) saves the DEM to the specified
%     filename
%
% Input arguments
%
%     DEM        instance of GRIDobj
%     filename   absolute or relative path and filename
%
% See also: GRIDobj
%
% Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
% Date: 17. August, 2017

narginchk(1,2)

% if only 1 argument, open file dialog box
if nargin == 1
    [FileName,PathName] = uiputfile({'*.tif'});
    if FileName == 0
        disp('     no output written to disk')
        return
    end
    file = [PathName FileName];
end

% try to use geotiffwrite, which comes with the Mapping Toolbox
try
    if isempty(A.georef);
        geotiffwrite(file,A.Z,A.refmat);
    else
        geotiffwrite(file,A.Z,A.georef.SpatialRef,...
            'GeoKeyDirectoryTag',A.georef.GeoKeyDirectoryTag);
    end
catch ME
    warning('TopoToolbox:GRIDobj',...
        ['GRIDobj2geotiff is unable to write a geotiff. Either you don''t \n'...
         'have the mapping toolbox, or there was another issue with geotiffwrite. \n'...
         'GRIDobj2geotiff instead writes a tif-image together with a world \n'...
         'file (*.tfw) which contains data on spatial referencing of the \n' ...
         'image, yet which lacks information on the type of projection used.']);

    % if geotiffwrite is not available or any other error occurs
    % a tif file will be written to the disk together with a worldfile
    % .tfw-file.
    [pathstr, name, ~] = fileparts(file);
    k = refmat2worldfile(A.refmat);
    dlmwrite(fullfile(pathstr,[name '.tfw']),k,'precision', '%.10f');
    A = A.Z;

    siz = size(A);
    cla = class(A);

    switch cla;
        case 'double'
            BpS = 64;
            TSF = Tiff.SampleFormat.IEEEFP;
        case 'single'
            BpS = 32;
            TSF = Tiff.SampleFormat.IEEEFP;
        otherwise
            if islogical(A);
                A = uint32(A);
                cla = 'uint32';
            end
            BpS = round(log2(double(intmax(cla))));
            TSF = Tiff.SampleFormat.UInt;

    end

    t = Tiff(file,'w');
    tagstruct.ImageLength = siz(1);
    tagstruct.ImageWidth = siz(2);
    tagstruct.BitsPerSample = BpS;
    tagstruct.SampleFormat = TSF;
    tagstruct.SamplesPerPixel = 1;
    tagstruct.RowsPerStrip = 16;
    tagstruct.PlanarConfiguration = Tiff.PlanarConfiguration.Chunky;
    tagstruct.Software = 'MATLAB';
    tagstruct.Photometric = 0;

    t.setTag(tagstruct);
    t.write(A);
    t.close;
end

end

function k = refmat2worldfile(r)
% does not support rotation

k(1,1) = r(2,1);
k(4,1) = r(1,2);
k(5,1) = r(3,1)+k(1);
k(6,1) = r(3,2)+k(4);
end

author: GERSL

repo: Fmask

save_path: github-repos/MATLAB/GERSL-Fmask


text:
function linplus_test265 ( )

%*****************************************************************************80
%
%% TEST265 tests R8GB_TO_R8GE, R8GE_TO_R8GB.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  m = 5;
  n = 8;
  ml = 2;
  mu = 1;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST265\n' );
  fprintf ( 1, '  R8GB_TO_R8GE copies a R8GB matrix to a R8GE matrix.\n' );
  fprintf ( 1, '  R8GE_TO_R8GB copies a R8GE matrix to a R8GB matrix.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix rows M =      %d\n', m );
  fprintf ( 1, '  Matrix columns N =   %d\n', n );
  fprintf ( 1, '  Lower bandwidth ML = %d\n', ml );
  fprintf ( 1, '  Upper bandwidth MU = %d\n', mu );

  a = r8gb_indicator ( m, n, ml, mu );

  r8gb_print ( m, n, ml, mu, a, '  The R8GB matrix:' );

  b = r8gb_to_r8ge ( m, n, ml, mu, a );

  r8ge_print ( m, n, b, '  The R8GE matrix:' );

  c = r8ge_to_r8gb ( m, n, ml, mu, b );

  r8gb_print ( m, n, ml, mu, c, '  The recovered R8GB matrix:' );

  return
end


text:
function [AerobicGrowth, AnaerobicGrowth] = testGrowth(model, biomassReaction)
% Tests growth, both on unlimited media and complex medium, both aerobic and
% anaerobic conditions. In anaerobic conditions, both oxygen uptake
% (exchange reaction) and cytosolic oxygen-utilizing reactions are blocked.
%
% USAGE
%           [AerobicGrowth, AnaerobicGrowth] = testGrowth(model, biomassReaction)
%
% INPUT
% model             COBRA model structure
% biomassReaction   String listing the biomass reaction
%
% OUTPUT
% AerobicGrowth     Numeric vector showing the flux through the biomass
%                   reaction (aerobic conditions):
%                   Column 1: "unlimited" media (all exchanges to -1000
%                   mmol/gDW/h)
%                   Column 2: complex medium as defined in
%                   InputFiles\ComplexMedium.txt
% AnaerobicGrowth   Numeric vector showing the flux through the biomass
%                   reaction (anaerobic conditions):
%                   Column 1: "unlimited" media (all exchanges to -1000
%                   mmol/gDW/h)
%                   Column 2: complex medium as defined in
%                   InputFiles\ComplexMedium.txt
%
% Stefania Magnusdottir, Nov 2017

tol=0.000001;

AerobicGrowth = zeros(1, 2);
AnaerobicGrowth = zeros(1, 2);

% Test if model can grow
% set "unlimited" constraints
model = changeRxnBounds(model, model.rxns(strncmp('EX_', model.rxns, 3)), -1000, 'l');
model = changeRxnBounds(model, model.rxns(strncmp('EX_', model.rxns, 3)), 1000, 'u');

% set objective
if nargin < 2
    error('Please provide biomass reaction')
end
if ~any(ismember(model.rxns, biomassReaction))
  %  error(['Biomass reaction "', biomassReaction, '" not found in model.'])
else
    model = changeObjective(model, biomassReaction);
end

% aerobic environment
modelO2 = changeRxnBounds(model, 'EX_o2(e)', -10, 'l');

% simulate
FBA = optimizeCbModel(modelO2, 'max');
if FBA.stat==1
    AerobicGrowth(1, 1) = FBA.f;
    if FBA.f > tol
        fprintf('Model grows on unlimited media (aerobic), flux through BOF: %d mmol/gDW/h\n', FBA.f)
    else
        warning('Model cannot grow on unlimited media (aerobic)')
    end
else
    warning('Model cannot grow on unlimited media (aerobic)')
end

% anaerobic environment
modelNoO2 = changeRxnBounds(model, 'EX_o2(e)', 0, 'l');

% simulate
FBA = optimizeCbModel(modelNoO2, 'max');
if FBA.stat==1
    AnaerobicGrowth(1, 1) = FBA.f;
    if FBA.f > tol
        fprintf('Model grows on unlimited media (anaerobic), flux through BOF: %d mmol/gDW/h\n', FBA.f)
    else
        warning('Model cannot grow on unlimited media (anaerobic)')
    end
else
    warning('Model cannot grow on unlimited media (anaerobic)')
end

% implement complex medium
constraints = readtable('ComplexMedium.txt', 'Delimiter', 'tab');
constraints=table2cell(constraints);
constraints=cellstr(string(constraints));

% apply complex medium
model = useDiet(model,constraints);

% aerobic environment
modelO2 = changeRxnBounds(model, 'EX_o2(e)', -10, 'l');

% simulate
FBA = optimizeCbModel(modelO2, 'max');
if FBA.stat==1
    AerobicGrowth(1, 2) = FBA.f;
    if FBA.f > tol
        fprintf('Model grows on complex medium (aerobic), flux through BOF: %d mmol/gDW/h\n', FBA.f)
    else
        warning('Model cannot grow on complex medium (aerobic)')
    end
else
    warning('Model cannot grow on complex medium (aerobic)')
end

% anaerobic environment
modelNoO2 = changeRxnBounds(model, 'EX_o2(e)', 0, 'l');

% simulate
FBA = optimizeCbModel(modelNoO2, 'max');
if FBA.stat==1
    AnaerobicGrowth(1, 2) = FBA.f;
    if FBA.f > tol
        fprintf('Model grows on complex medium (anaerobic), flux through BOF: %d mmol/gDW/h\n', FBA.f)
    else
        warning('Model cannot grow on complex medium (anaerobic)')
    end
else
    warning('Model cannot grow on complex medium (anaerobic)')
end

end

author: opencobra

repo: cobratoolbox


text:
function [VM, hrf, fit, e, param] = Det_Logit(V0,t,tc,Run)
%
% [VM, h, fit, e, param] = Det_Logit_allstim(V0,t,tc,Run)
%
% Estimate inverse logit (IL) HRF model
% Creates fitted curve - 3 logistic functions to be summed together - from parameter estimates
%
% INPUT: V0, t, tc, Run
% Run = stick function
% tc = time course
% t = vector of time points
% V0 = initial value for the parameter vector
%
% By Martin Lindquist, Christian Waugh and Tor Wager
% Created by Martin Lindquist on 10/02/09
% Last edited: 05/26/10 (ML)

numstim = length(Run);
len = length(Run{1});

% LB = [0.05, 1, 0, 0.05, 4, 0, 10];      % Lower bounds for parameters
% UB = [10, 15, 10, 10, 15, 5, 50];           % Upper bounds for parameters
% LB = repmat(LB, 1, numstim);
% UB = repmat(UB, 1, numstim);

% Remove intercept

b0 = pinv(ones(length(tc),1))*tc;
tc = tc - b0;

% Find optimal values

options = optimset('MaxFunEvals',10000000,'Maxiter',10000000,'TolX',1e-8,'TolFun',1e-8,'Display','off');

%VM = fminsearchbnd(@cost_allstim, V0, LB,UB,options,t,tc,Run);
VM = fminsearch(@msq_logit,V0,options,Run,t,tc);

% Use optimal values to fit hemodynamic response functions
hrf =zeros(length(t),numstim);
fitt = zeros(len,numstim);
param = zeros(3,numstim);

for g = 1:numstim
    hrf(:,g) = il_hdmf_tw2(t,VM(((g-1)*7+1):(g*7)));                   % Calculate HRF estimate (fit, given theta)
    param(:,g) = get_parameters2(hrf(:,g),t(end));
    fits(:,g) = conv(Run{g}, hrf(:,g));
    fitt(:,g) = fits(1:len,g);
end

fit = sum(fitt,2);
e = tc-fit;
fit = fit + b0;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   SUBFUNCTIONS
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function m=msq_logit(V,Run, t, tc)

numstim = length(Run);
len = length(Run{1});
h = zeros(length(t),numstim);
yhatt =zeros(len,numstim);

for k = 1:numstim
    h(:,k) = il_hdmf_tw2(t,V(((k-1)*7+1):(k*7)));           % Get IL model corresponding to parameters V
    yhat(:,k) = conv(Run{k}, h(:,k));                     % Convolve IL model with stick function
    yhatt(:,k) = yhat(1:len,k);
end

yhat2 = sum(yhatt,2); %Sum models together to get overall estimate

m = sum((tc-yhat2).^2);              % Calculate cost function

end

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

function [h,base] = il_hdmf_tw2(t,V)
% inverse logit -- creates fitted curve from parameter estimates
%
% t = vector of time points
% V = parameters

% 3 logistic functions to be summed together
base = zeros(length(t),3);
A1 = V(1);
T1 = V(2);
d1 = V(3);
A2 = V(4);
T2 = V(5);
A3 = V(6);
T3 = V(7);
d2 = -d1*(ilogit(A1*(1-T1)) - ilogit(A3*(1-T3)))/(ilogit(A2*(1-T2)) + ilogit(A3*(1-T3)));
d3 = abs(d2)-abs(d1);

base(:,1)= d1*ilogit(A1*(t-T1))';
base(:,2)= d2*ilogit(A2*(t-T2))';
base(:,3)= d3*ilogit(A3*(t-T3))';
h = sum(base,2)';

end

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

function [L] = ilogit(t)
L = exp(t)./(1+exp(t));
end

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

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%% DEMO_febio_0047_cylinder_embedded_probe_01
% Below is a demonstration for:
%
% * Building geometry for a tissue segment with an embedded probe
% * Defining the boundary conditions
% * Coding the febio structure
% * Running the model
% * Importing and visualizing the displacement results

%% Keywords
%
% * febio_spec version 3.0
% * febio, FEBio
% * probe
% * rigid body constraints
% * tetrahedral elements, tet4
% * triangular elements, tri3
% * slab, block, rectangular
% * sphere
% * static, solid
% * hyperelastic, Ogden
% * displacement logfile
% * stress logfile

%%
clear; close all; clc;

%%
% Plot settings
fontSize=15;
faceAlpha=1;
lineWidth1=1.5;
lineWidth2=3;
markerSize1=15;
markerSize2=30;
edgeWidth=2;
edgeColor='k';
faceAlpha1=1;

%% Control parameters

% Path names
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
savePath=fullfile(defaultFolder,'data','temp');

% Defining file names
febioFebFileNamePart='tempModel';
febioFebFileName=fullfile(savePath,[febioFebFileNamePart,'.feb']); %FEB file name
febioLogFileName=[febioFebFileNamePart,'.txt']; %FEBio log file name
febioLogFileName_disp=[febioFebFileNamePart,'_disp_out.txt']; %Log file name for exporting displacement
febioLogFileName_force=[febioFebFileNamePart,'_force_out.txt']; %Log file name for exporting force
febioLogFileName_strainEnergy=[febioFebFileNamePart,'_energy_out.txt']; %Log file name for exporting strain energy density
febioLogFileName_strain=[febioFebFileNamePart,'_strain_out.txt']; %Log file name for exporting strain

probeHeight=75;

probeRadius=2; % The radius of the hemi-spher portion
nRefine=0; % Number of |subtri| refinements for icosahedron

pointSpacing=3;
dAdd=7;
tissueRadius=probeRadius+dAdd;
tissueHeight=probeHeight+dAdd;
volumeFactor=5;

displacementMagnitude=-1;

%Material parameter set
c1=1e-3; %Shear-modulus-like parameter
m1=2; %Material parameter setting degree of non-linearity
k_factor=1e2; %Bulk modulus factor
k=c1*k_factor; %Bulk modulus

% FEA control settings
numTimeSteps=10; %Number of time steps desired
max_refs=25; %Max reforms
max_ups=0; %Set to zero to use full-Newton iterations
opt_iter=6; %Optimum number of iterations
max_retries=5; %Maximum number of retires
dtmin=(1/numTimeSteps)/100; %Minimum time step size
dtmax=1/numTimeSteps; %Maximum time step size

%% Build probe

probeMeshInputStruct.sphereRadius=probeRadius;% => The radius of the hemi-spher portion
probeMeshInputStruct.nRefine=nRefine;% => Number of |subtri| refinements for icosahedron
probeMeshInputStruct.cylinderHeight=probeHeight-probeRadius;% => height of the cylinder part
probeMeshInputStruct.cylinderStepSize=[];% => Aproximate node spacing for cylinder portion
probeMeshInputStruct.patchType='tri';

[Fp,Vp,Cp]=hemiSphereCylMesh(probeMeshInputStruct);
Fp=fliplr(Fp); %Invert face orientation

Vp(:,3)=Vp(:,3)-max(Vp(:,3));

%Get top curve
Eb=patchBoundary(Fp);
indProbeTop=edgeListToCurve(Eb);
indProbeTop=indProbeTop(1:end-1);
Vst=Vp(indProbeTop,:);

%%
cFigure; hold on;
gpatch(Fp,Vp,'gw','k');
plotV(Vst,'b.-','lineWidth',lineWidth1,'MarkerSize',markerSize1);
patchNormPlot(Fp,Vp);
axisGeom(gca,fontSize);
drawnow;

%% Build tissue

%Sketching profile
ns=150;
t=linspace(0,2*pi,ns);
t=t(1:end-1);

x=tissueRadius*cos(t);
y=tissueRadius*sin(t);
z=zeros(size(x));
Vc=[x(:) y(:) z(:)];
np=ceil(max(pathLength(Vc))./pointSpacing);
[Vc]=evenlySampleCurve(Vc,np,'pchip',1);

% Extruding model
cPar.numSteps=round(tissueHeight/pointSpacing);
cPar.depth=tissueHeight;
cPar.patchType='tri';
cPar.dir=-1;
cPar.closeLoopOpt=1;
[Fg,Vg]=polyExtrude(Vc,cPar);
Fg=fliplr(Fg);

Vgb=Vg(cPar.numSteps:cPar.numSteps:end,:);

Vgt=Vg(1:cPar.numSteps:end,:);

%% Cap ends

regionCell={Vgt(:,[1 2]),Vst(:,[1 2])};
[Ft,Vt]=regionTriMesh2D(regionCell,pointSpacing,0,0);
Vt(:,3)=mean(Vgt(:,3));

regionCell={Vgb(:,[1 2])};
[Fb,Vb]=regionTriMesh2D(regionCell,pointSpacing,0,0);
Fb=fliplr(Fb); %flip face orientation
Vb(:,3)=mean(Vgb(:,3));

%%
% Visualize

cFigure; hold on;

gpatch(Fp,Vp,'rw','k',0.5);
gpatch(Fg,Vg,'gw','k',0.5);
gpatch(Fb,Vb,'bw','k',0.5);
gpatch(Ft,Vt,'bw','k',0.5);

plotV(Vgb,'b.-','lineWidth',lineWidth1,'MarkerSize',markerSize1);
plotV(Vgt,'b.-','lineWidth',lineWidth1,'MarkerSize',markerSize1);
plotV(Vst,'b.-','lineWidth',lineWidth1,'MarkerSize',markerSize1);

axisGeom(gca,fontSize);
drawnow;

%% Merge model components

[F,V,C]=joinElementSets({Fg,Ft,Fb,Fp},{Vg,Vt,Vb,Vp});
[F,V]=mergeVertices(F,V);

%%

cFigure;
subplot(1,2,1); hold on;
gpatch(F,V,C,'none',0.5);
axisGeom(gca,fontSize);
colormap gjet; icolorbar;

subplot(1,2,2); hold on;
gpatch(F,V,C);
patchNormPlot(F,V,2);
plotV(Vst,'b.-','lineWidth',lineWidth1,'MarkerSize',markerSize1);
axisGeom(gca,fontSize);
colormap gjet; icolorbar;

drawnow;

%% Mesh solid using tetgen

%%
% Create tetgen meshing input structure

[regionA]=tetVolMeanEst(F,V); %Volume for a regular tet based on edge lengths
V_inner=getInnerPoint(F,V); %Interior point for region

inputStruct.stringOpt='-pq1.2AaY';
inputStruct.Faces=F;
inputStruct.Nodes=V;
inputStruct.holePoints=[];
inputStruct.faceBoundaryMarker=C; %Face boundary markers
inputStruct.regionPoints=V_inner; %region points
inputStruct.regionA=regionA*volumeFactor; %Desired volume for tets
inputStruct.minRegionMarker=2; %Minimum region marker

%%
% Mesh model using tetrahedral elements using tetGen

[meshOutput]=runTetGen(inputStruct); %Run tetGen

%%
% Visualize mesh

meshView(meshOutput);

%%
% Access model element and patch data
F=meshOutput.faces;
V=meshOutput.nodes;
C=meshOutput.faceMaterialID;
E=meshOutput.elements;
elementMaterialID=meshOutput.elementMaterialID;

Fb=meshOutput.facesBoundary;
Cb=meshOutput.boundaryMarker;

%% Define boundary condition node sets

logicRigid= Cb==1 | Cb==3;
bcSupportList=Fb(logicRigid,:);
bcSupportList=unique(bcSupportList(:));

logicIndentor= Cb==4;
bcPrescribeList=Fb(logicIndentor,:);
bcPrescribeList=unique(bcPrescribeList(:));

%%
% Visualize boundary conditions

cFigure; hold on;
gpatch(Fb,V,'bw','none',0.5);
hp(1)=plotV(V(bcSupportList,:),'k.','lineWidth',lineWidth1,'MarkerSize',markerSize1);
hp(2)=plotV(V(bcPrescribeList,:),'r.','lineWidth',lineWidth1,'MarkerSize',markerSize1);
legend(hp,{'BC Full support','BC Prescribed displacement'});
axisGeom(gca,fontSize);
drawnow;

%% Defining the FEBio input structure
% See also |febioStructTemplate| and |febioStruct2xml| and the FEBio user
% manual.

%Get a template with default settings
[febio_spec]=febioStructTemplate;

%febio_spec version
febio_spec.ATTR.version='3.0';

%Module section
febio_spec.Module.ATTR.type='solid';

%Create control structure for use by all steps
stepStruct.Control.time_steps=numTimeSteps;
stepStruct.Control.step_size=1/numTimeSteps;
stepStruct.Control.solver.max_refs=max_refs;
stepStruct.Control.solver.max_ups=max_ups;
stepStruct.Control.time_stepper.dtmin=dtmin;
stepStruct.Control.time_stepper.dtmax=dtmax;
stepStruct.Control.time_stepper.max_retries=max_retries;
stepStruct.Control.time_stepper.opt_iter=opt_iter;

%Add template based default settings to proposed control section
[stepStruct.Control]=structComplete(stepStruct.Control,febio_spec.Control,1); %Complement provided with default if missing

%Remove control field (part of template) since step specific control sections are used
febio_spec=rmfield(febio_spec,'Control');

febio_spec.Step.step{1}.Control=stepStruct.Control;
febio_spec.Step.step{1}.ATTR.id=1;
febio_spec.Step.step{2}.Control=stepStruct.Control;
febio_spec.Step.step{2}.ATTR.id=2;

%Material section
materialName1='Material1';
febio_spec.Material.material{1}.ATTR.name=materialName1;
febio_spec.Material.material{1}.ATTR.type='Ogden';
febio_spec.Material.material{1}.ATTR.id=1;
febio_spec.Material.material{1}.c1=c1;
febio_spec.Material.material{1}.m1=m1;
febio_spec.Material.material{1}.c2=c1;
febio_spec.Material.material{1}.m2=-m1;
febio_spec.Material.material{1}.k=k;
% Mesh section
% -> Nodes
febio_spec.Mesh.Nodes{1}.ATTR.name='Object1'; %The node set name
febio_spec.Mesh.Nodes{1}.node.ATTR.id=(1:size(V,1))'; %The node id's
febio_spec.Mesh.Nodes{1}.node.VAL=V; %The nodel coordinates

% -> Elements
partName1='Part1';
febio_spec.Mesh.Elements{1}.ATTR.name=partName1; %Name of this part
febio_spec.Mesh.Elements{1}.ATTR.type='tet4'; %Element type
febio_spec.Mesh.Elements{1}.elem.ATTR.id=(1:1:size(E,1))'; %Element id's
febio_spec.Mesh.Elements{1}.elem.VAL=E; %The element matrix

% -> NodeSets
nodeSetName1='bcSupportList';
nodeSetName2='bcPrescribeList';

febio_spec.Mesh.NodeSet{1}.ATTR.name=nodeSetName1;
febio_spec.Mesh.NodeSet{1}.node.ATTR.id=bcSupportList(:);

febio_spec.Mesh.NodeSet{2}.ATTR.name=nodeSetName2;
febio_spec.Mesh.NodeSet{2}.node.ATTR.id=bcPrescribeList(:);

%MeshDomains section
febio_spec.MeshDomains.SolidDomain.ATTR.name=partName1;
febio_spec.MeshDomains.SolidDomain.ATTR.mat=materialName1;

%Boundary condition section
%-> Fix boundary conditions
febio_spec.Boundary.bc{1}.ATTR.type='fix';
febio_spec.Boundary.bc{1}.ATTR.node_set=nodeSetName1;
febio_spec.Boundary.bc{1}.dofs='x,y,z';

% -> Prescribe boundary conditions
%STEP 1 Up/down
febio_spec.Step.step{1}.Boundary.bc{1}.ATTR.type='prescribe';
febio_spec.Step.step{1}.Boundary.bc{1}.ATTR.node_set=nodeSetName2;
febio_spec.Step.step{1}.Boundary.bc{1}.dof='z';
febio_spec.Step.step{1}.Boundary.bc{1}.scale.ATTR.lc=1;
febio_spec.Step.step{1}.Boundary.bc{1}.scale.VAL=displacementMagnitude;
febio_spec.Step.step{1}.Boundary.bc{1}.relative=1;

febio_spec.Step.step{1}.Boundary.bc{2}.ATTR.type='fix';
febio_spec.Step.step{1}.Boundary.bc{2}.ATTR.node_set=nodeSetName2;
febio_spec.Step.step{1}.Boundary.bc{2}.dofs='x,y';

%STEP 2 Sideways
febio_spec.Step.step{2}.Boundary.bc{1}.ATTR.type='prescribe';
febio_spec.Step.step{2}.Boundary.bc{1}.ATTR.node_set=nodeSetName2;
febio_spec.Step.step{2}.Boundary.bc{1}.dof='x';
febio_spec.Step.step{2}.Boundary.bc{1}.scale.ATTR.lc=2;
febio_spec.Step.step{2}.Boundary.bc{1}.scale.VAL=displacementMagnitude;
febio_spec.Step.step{2}.Boundary.bc{1}.relative=1;

febio_spec.Step.step{2}.Boundary.bc{2}.ATTR.type='fix';
febio_spec.Step.step{2}.Boundary.bc{2}.ATTR.node_set=nodeSetName2;
febio_spec.Step.step{2}.Boundary.bc{2}.dofs='y,z';

%LoadData section
% -> load_controller
febio_spec.LoadData.load_controller{1}.ATTR.id=1;
febio_spec.LoadData.load_controller{1}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{1}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{1}.points.point.VAL=[0 0; 0.25 1; 0.5 0; 0.75 -1; 1 0];

febio_spec.LoadData.load_controller{2}.ATTR.id=2;
febio_spec.LoadData.load_controller{2}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{2}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{2}.points.point.VAL=[1 0; 1.25 1; 1.5 0; 1.75 -1; 2 0];

%Output section
% -> log file
febio_spec.Output.logfile.ATTR.file=febioLogFileName;
febio_spec.Output.logfile.node_data{1}.ATTR.file=febioLogFileName_disp;
febio_spec.Output.logfile.node_data{1}.ATTR.data='ux;uy;uz';
febio_spec.Output.logfile.node_data{1}.ATTR.delim=',';

febio_spec.Output.logfile.node_data{2}.ATTR.file=febioLogFileName_force;
febio_spec.Output.logfile.node_data{2}.ATTR.data='Rx;Ry;Rz';
febio_spec.Output.logfile.node_data{2}.ATTR.delim=',';

febio_spec.Output.logfile.element_data{1}.ATTR.file=febioLogFileName_strain;
febio_spec.Output.logfile.element_data{1}.ATTR.data='E1;E2;E3';
febio_spec.Output.logfile.element_data{1}.ATTR.delim=',';

%% Quick viewing of the FEBio input file structure
% The |febView| function can be used to view the xml structure in a MATLAB
% figure window.

%%
% |febView(febio_spec); %Viewing the febio file|

%% Exporting the FEBio input file
% Exporting the febio_spec structure to an FEBio input file is done using
% the |febioStruct2xml| function.

febioStruct2xml(febio_spec,febioFebFileName); %Exporting to file and domNode

%% Running the FEBio analysis
% To run the analysis defined by the created FEBio input file the
% |runMonitorFEBio| function is used. The input for this function is a
% structure defining job settings e.g. the FEBio input file name. The
% optional output runFlag informs the user if the analysis was run
% succesfully.

febioAnalysis.run_filename=febioFebFileName; %The input file name
febioAnalysis.run_logname=febioLogFileName; %The name for the log file
febioAnalysis.disp_on=1; %Display information on the command window
febioAnalysis.runMode='external';%'internal' or 'external';

[runFlag]=runMonitorFEBio(febioAnalysis);%START FEBio NOW!!!!!!!!

%% Import FEBio results

if runFlag==1 %i.e. a succesful run

    % Importing nodal displacements from a log file
    [time_mat, N_disp_mat,~]=importFEBio_logfile(fullfile(savePath,febioLogFileName_disp)); %Nodal displacements
    time_mat=[0; time_mat(:)]; %Time

    N_disp_mat=N_disp_mat(:,2:end,:);
    sizImport=size(N_disp_mat);
    sizImport(3)=sizImport(3)+1;
    N_disp_mat_n=zeros(sizImport);
    N_disp_mat_n(:,:,2:end)=N_disp_mat;
    N_disp_mat=N_disp_mat_n;

    DN_MAG=sqrt(sum(N_disp_mat.^2,2));
    DN=N_disp_mat(:,:,end);
    DN_magnitude=sqrt(sum(DN(:,3).^2,2));
    V_def=V+DN;
    V_DEF=N_disp_mat+repmat(V,[1 1 size(N_disp_mat,3)]);
    X_DEF=V_DEF(:,1,:);
    Y_DEF=V_DEF(:,2,:);
    Z_DEF=V_DEF(:,3,:);
    [CF_def]=vertexToFaceMeasure(Fb,DN_magnitude);

     %%
    % Importing element data from a log file
    [~,E_data,~]=importFEBio_logfile(fullfile(savePath,febioLogFileName_strain)); %Element data

    %Remove nodal index column
    E_data=E_data(:,2:end,:);

    E_data=sqrt(0.5*( (E_data(:,1,:)-E_data(:,2,:)).^2 + (E_data(:,2,:)-E_data(:,3,:)).^2 + (E_data(:,3,:)-E_data(:,1,:)).^2 ));

    %Add initial state i.e. zero
    sizImport=size(E_data);
    sizImport(3)=sizImport(3)+1;
    E_data_mat_n=zeros(sizImport);
    E_data_mat_n(:,:,2:end)=E_data;
    E_data=E_data_mat_n;

    VE=patchCentre(E,V);
    logicCutElements=VE(:,2)>=0;

    [F_cut,CF_cut_data]=element2patch(E(logicCutElements,:),E_data(logicCutElements,:,1));
    [indBoundary]=tesBoundary(F_cut);

    CV=faceToVertexMeasure(F_cut(indBoundary,:),V,CF_cut_data(indBoundary,:));

    %%
    % Plotting the simulated results using |anim8| to visualize and animate
    % deformations

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; %Open figure
    gtitle([febioFebFileNamePart,': Press play to animate']);
%     subplot(1,2,1);
    hp1=gpatch(Fb,V_def,'kw','none',0.25); %Add graphics object to animate
    hp2=gpatch(F_cut(indBoundary,:),V_def,CV,'k',1); %Add graphics object to animate
    hp2.FaceColor='interp';
%     gpatch(Fb,V,0.5*ones(1,3),'k',0.25); %A static graphics object
    colormap(gjet(250)); colorbar;
    caxis([0 max(E_data(:))/3]);
    axisGeom(gca,fontSize);
    axis([min(X_DEF(:)) max(X_DEF(:)) min(Y_DEF(:)) max(Y_DEF(:)) min(Z_DEF(:)) max(Z_DEF(:))]);
    axis manual;
    camlight headlight;

    drawnow;

    % Set up animation features
    animStruct.Time=time_mat; %The time vector
    for qt=1:1:size(N_disp_mat,3) %Loop over time increments
        DN=N_disp_mat(:,:,qt); %Current displacement
        DN_magnitude=sqrt(sum(DN.^2,2)); %Current displacement magnitude
        V_def=V+DN; %Current nodal coordinates
%         [CF_def]=vertexToFaceMeasure(Fb,DN_magnitude); %Current color data to use

        [~,CF_cut_data]=element2patch(E(logicCutElements,:),E_data(logicCutElements,:,qt));
        CV=faceToVertexMeasure(F_cut(indBoundary,:),V,CF_cut_data(indBoundary,:));

        %Set entries in animation structure
        animStruct.Handles{qt}=[hp1 hp2 hp2]; %Handles of objects to animate
        animStruct.Props{qt}={'Vertices','Vertices','CData'}; %Properties of objects to animate
        animStruct.Set{qt}={V_def,V_def,CV}; %Property values for to set in order to animate
    end
    anim8(hf,animStruct); %Initiate animation feature
    drawnow;

end

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

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON