text:
function [net,stats, ders_iter] = cnn_train_dag_ridge(net, imdb, input, getBatch, varargin)
%CNN_TRAIN_DAG Demonstrates training a CNN using the DagNN wrapper
%    CNN_TRAIN_DAG() is similar to CNN_TRAIN(), but works with
%    the DagNN wrapper instead of the SimpleNN wrapper.

% Copyright (C) 2014-16 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

% if show test1, please assign true to is_show_test1
is_show_test1 = true;%false;%

addpath(fullfile(vl_rootnn, 'examples'));

opts.expDir = fullfile('data','exp') ;
opts.continue = true ;
opts.batchSize = 256 ;
opts.numSubBatches = 1 ;
opts.train = [] ;
opts.val = [] ;
opts.gpus = [] ;
opts.prefetch = false ;
opts.numEpochs = 300 ;
opts.learningRate = 0.001 ;
opts.weightDecay = 0.0005 ;

opts.solver = @solver.adam; % Empty array - optimised SGD solver
[opts, varargin] = vl_argparse(opts, varargin);
if isempty(opts.solver)
  opts.solverOpts.momentum = 0.9;
else
  assert(isa(opts.solver, 'function_handle') && nargout(opts.solver) == 2,...
    'Invalid solver - a function handle with two outputs expected.');
  % A call without any input arg - def opts
  opts.solverOpts = opts.solver();
end

opts.saveSolverState = true ;
opts.randomSeed = 0 ;
opts.profile = false ;
opts.parameterServer.method = 'mmap' ;
opts.parameterServer.prefix = 'mcn' ;

opts.derOutputs = {'objective_r', 1} ;
opts.extractStatsFn = @extractStats ;
opts.plotStatistics = false;
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
% if isnan(opts.val), opts.val = [] ; end

% -------------------------------------------------------------------------
%                                                            Initialization
% -------------------------------------------------------------------------

evaluateMode = isempty(opts.train) ;
if ~evaluateMode
  if isempty(opts.derOutputs)
    error('DEROUTPUTS must be specified when training.\n') ;
  end
end

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

modelPath = @(ep) fullfile(opts.expDir, sprintf('net-epoch-%d.mat', ep));
modelFigPath = fullfile(opts.expDir, 'net-train.pdf') ;

start = opts.continue * findLastCheckpoint(opts.expDir) ;
if start >= 1
  fprintf('%s: resuming by loading epoch %d\n', mfilename, start) ;
  [net, state, stats] = loadState(modelPath(start)) ;
else
  state = [] ;
end

global out;
out=false;

% if is_show_test1
%     act_num = zeros(512,1);%-test1
%     epoch=0;
    ders_iter = zeros(size(input{1},3),opts.numEpochs);
% end
for epoch=start+1:opts.numEpochs

    if out
        break;
    end

  % Set the random seed based on the epoch and opts.randomSeed.
  % This is important for reproducibility, including when training
  % is restarted from a checkpoint.

  rng(epoch + opts.randomSeed) ;
  %prepareGPUs(opts, epoch == start+1) ;

  % Train for one epoch.
  params = opts ;
  params.epoch = epoch ;
  params.learningRate = opts.learningRate(min(epoch, numel(opts.learningRate))) ;
  params.train = opts.train(randperm(numel(opts.train))) ; % shuffle
  params.val = opts.val(randperm(numel(opts.val))) ;
  %params.imdb = imdb ;
  params.getBatch = getBatch ;

  if numel(opts.gpus) <= 1
    [net, state, ders_iter(:,epoch)] = processEpoch(net, input, state, params, 'train') ;
    [net, state] = processEpoch(net, input, state, params, 'val') ;
    if ~evaluateMode
      %saveState(modelPath(epoch), net, state) ;
    end
    lastStats = state.stats ;
  else
    spmd
      [net, state, ders_iter(:,epoch)] = processEpoch(net, input, state, params, 'train') ;
      [net, state] = processEpoch(net, input, state, params, 'val') ;
      if labindex == 1 && ~evaluateMode
        %saveState(modelPath(epoch), net, state) ;
      end
      lastStats = state.stats ;
    end
    lastStats = accumulateStats(lastStats) ;
  end

  stats.train(epoch) = lastStats.train ;
  stats.val(epoch) = lastStats.val ;
  clear lastStats ;
  %saveStats(modelPath(epoch), stats) ;

  if opts.plotStatistics
    switchFigure(1) ; clf ;
    plots = setdiff(...
      cat(2,...
      fieldnames(stats.train)', ...
      fieldnames(stats.val)'), {'num', 'time'}) ;
    for p = plots
      p = char(p) ;
      values = zeros(0, epoch) ;
      leg = {} ;
      for f = {'train', 'val'}
        f = char(f) ;
        if isfield(stats.(f), p)
          tmp = [stats.(f).(p)] ;
          values(end+1,:) = tmp(1,:)' ;
          leg{end+1} = f ;
        end
      end
      subplot(1,numel(plots),find(strcmp(p,plots))) ;
      plot(1:epoch, values','o-') ;
      xlabel('epoch') ;
      title(p) ;
      legend(leg{:}) ;
      grid on ;
    end
    drawnow ;
    print(1, modelFigPath, '-dpdf') ;
  end
end

% if is_show_test1
%     %-test1
%     aa=zeros(epoch,1);
%      for i=1:epoch
%         aa(i) = sum(act_num>i-1 );
%      end
%      figure,plot(aa);
%      xlabel('adjust times');
%      ylabel('channel numbers');
%      title('adjust times vs channel numbers');
%     %-test1
% end
% With multiple GPUs, return one copy
if isa(net, 'Composite'), net = net{1} ; end

% -------------------------------------------------------------------------
function [net, state, ders] = processEpoch(net, input, state, params, mode)
% -------------------------------------------------------------------------
% Note that net is not strictly needed as an output argument as net
% is a handle class. However, this fixes some aliasing issue in the
% spmd caller.

% initialize with momentum 0
if isempty(state) || isempty(state.solverState)
  state.solverState = cell(1, numel(net.params)) ;
end

% move CNN  to GPU as needed
numGpus = numel(params.gpus) ;
if numGpus >= 1
  net.move('gpu') ;
  %state.momentum = cellfun(@gpuArray, state.momentum, 'uniformoutput', false) ;
end
if numGpus > 1
  parserv = ParameterServer(params.parameterServer) ;
  net.setParameterServer(parserv) ;
else
  parserv = [] ;
end

% profile
if params.profile
  if numGpus <= 1
    profile clear ;
    profile on ;
  else
    mpiprofile reset ;
    mpiprofile on ;
  end
end

num = 0 ;
epoch = params.epoch ;
subset = params.(mode) ;
adjustTime = 0 ;

stats.num = 0 ; % return something even if subset = []
stats.time = 0 ;

start = tic ;
for t=1:params.batchSize:numel(subset)
%   fprintf('%s: epoch %02d: %3d/%3d:', mode, epoch, ...
%           fix((t-1)/params.batchSize)+1, ceil(numel(subset)/params.batchSize)) ;
  batchSize = min(params.batchSize, numel(subset) - t + 1) ;

  for s=1:params.numSubBatches
    % get this image batch and prefetch the next
    batchStart = t + (labindex-1) + (s-1) * numlabs ;
    batchEnd = min(t+params.batchSize-1, numel(subset)) ;
    batch = subset(batchStart : params.numSubBatches * numlabs : batchEnd) ;
    num = num + numel(batch) ;
    if numel(batch) == 0, continue ; end

    feats=input{1};
    labels=input{2};

%lx     inputs ={'input1', gpuArray(ims),'input2', ...
%          gpuArray(ims), 'label', gpuArray(labels)};
    inputs ={'input', feats, 'label_gaussian', labels};

    if params.prefetch
      if s == params.numSubBatches
        batchStart = t + (labindex-1) + params.batchSize ;
        batchEnd = min(t+2*params.batchSize-1, numel(subset)) ;
      else
        batchStart = batchStart + numlabs ;
      end
      nextBatch = subset(batchStart : params.numSubBatches * numlabs : batchEnd) ;
      %params.getBatch(params.imdb, nextBatch) ;
    end

    if strcmp(mode, 'train')
      net.mode = 'normal' ;
      net.accumulateParamDers = (s ~= 1) ;
      net.eval(inputs, params.derOutputs, 'holdOn', s < params.numSubBatches) ;
%         if is_show_test1
%               showing gradient value in each depth channel %-test1
              ders = squeeze(gather(sum(sum(net.vars(1).der))));
%                 ders = squeeze(gather(sum(sum(net.params(1).der))));
%               figure(3);
%               plot(a);
%               xlabel('channel number'); ylabel('value');title('gradient weights');
%               [val, ind] = sort(a,'descend');
%               act_num(ind(1:64))= act_num(ind(1:64))+1;  %-test1
%         end
    else
      net.mode = 'test' ;
      net.eval(inputs) ;
    end
  end

  % Accumulate gradient.
  if strcmp(mode, 'train')
    if ~isempty(parserv), parserv.sync() ; end
    state = accumulateGradients(net, state, params, batchSize, parserv) ;
  end

  % Get statistics.
  time = toc(start) + adjustTime ;
  batchTime = time - stats.time ;
  stats.num = num ;
  stats.time = time ;
  stats = params.extractStatsFn(stats,net) ;
  currentSpeed = batchSize / batchTime ;
  averageSpeed = (t + batchSize - 1) / time ;
  if t == 3*params.batchSize + 1
    % compensate for the first three iterations, which are outliers
    adjustTime = 4*batchTime - time ;
    stats.time = time + adjustTime ;
  end

%   fprintf(' %.1f (%.1f) Hz', averageSpeed, currentSpeed) ;
  for f = setdiff(fieldnames(stats)', {'num', 'time'})
    f = char(f) ;
    fprintf(' %s: %.8f', f, stats.(f)) ;

    if strcmp(f,'objective_r')&& stats.(f)<0.02
        global out;
        out=true;
    end
    break;
  end
  fprintf('\n') ;
end

% Save back to state.
state.stats.(mode) = stats ;
if params.profile
  if numGpus <= 1
    state.prof.(mode) = profile('info') ;
    profile off ;
  else
    state.prof.(mode) = mpiprofile('info');
    mpiprofile off ;
  end
end
if ~params.saveSolverState
  state.solverState = [] ;
else
  state.solverState = cellfun(@gather, state.solverState, 'uniformoutput', false) ;
end

net.reset() ;
%net.move('cpu') ;

% -------------------------------------------------------------------------
function state = accumulateGradients(net, state, params, batchSize, parserv)
% -------------------------------------------------------------------------
numGpus = numel(params.gpus) ;
otherGpus = setdiff(1:numGpus, labindex) ;

for p=1:numel(net.params)

  if ~isempty(parserv)
    parDer = parserv.pullWithIndex(p) ;
  else
    parDer = net.params(p).der ;
  end

  switch net.params(p).trainMethod

    case 'average' % mainly for batch normalization
      thisLR = net.params(p).learningRate ;
      net.params(p).value = vl_taccum(...
          1 - thisLR, net.params(p).value, ...
          (thisLR/batchSize/net.params(p).fanout),  parDer) ;

    case 'gradient'
      thisDecay = params.weightDecay * net.params(p).weightDecay ;
      thisLR = params.learningRate * net.params(p).learningRate ;

      if isempty(params.solver)
        if isempty(state.solverState{p})
          state.solverState{p} = zeros(size(parDer), 'like', parDer);
        end

        state.solverState{p} = vl_taccum(...
          params.solverOpts.momentum,  state.solverState{p}, ...
          - (1 / batchSize), parDer) ;
        net.params(p).value = vl_taccum(...
          (1 - thisLR * thisDecay / (1 - params.solverOpts.momentum)),  ...
          net.params(p).value, ...
          thisLR, state.solverState{p}) ;
      else
        grad = (1 / batchSize) * parDer + thisDecay * net.params(p).value;
        % call solver function to update weights
        [net.params(p).value, state.solverState{p}] = ...
          params.solver(net.params(p).value, state.solverState{p}, ...
          grad, params.solverOpts, thisLR) ;
      end

    otherwise
      error('Unknown training method ''%s'' for parameter ''%s''.', ...
        net.params(p).trainMethod, ...
        net.params(p).name) ;
  end
end

% -------------------------------------------------------------------------
function stats = accumulateStats(stats_)
% -------------------------------------------------------------------------

for s = {'train', 'val'}
  s = char(s) ;
  total = 0 ;

  % initialize stats stucture with same fields and same order as
  % stats_{1}
  stats__ = stats_{1} ;
  names = fieldnames(stats__.(s))' ;
  values = zeros(1, numel(names)) ;
  fields = cat(1, names, num2cell(values)) ;
  stats.(s) = struct(fields{:}) ;

  for g = 1:numel(stats_)
    stats__ = stats_{g} ;
    num__ = stats__.(s).num ;
    total = total + num__ ;

    for f = setdiff(fieldnames(stats__.(s))', 'num')
      f = char(f) ;
      stats.(s).(f) = stats.(s).(f) + stats__.(s).(f) * num__ ;

      if g == numel(stats_)
        stats.(s).(f) = stats.(s).(f) / total ;
      end
    end
  end
  stats.(s).num = total ;
end

% -------------------------------------------------------------------------
function stats = extractStats(stats, net)
% -------------------------------------------------------------------------
sel = find(cellfun(@(x) isa(x,'dagnn.Loss'), {net.layers.block})) ;
for i = 1:numel(sel)
  stats.(net.layers(sel(i)).outputs{1}) = net.layers(sel(i)).block.average ;
end

% -------------------------------------------------------------------------
function saveState(fileName, net_, state)
% -------------------------------------------------------------------------
net = net_.saveobj() ;
save(fileName, 'net', 'state') ;

% -------------------------------------------------------------------------
function saveStats(fileName, stats)
% -------------------------------------------------------------------------
if exist(fileName)
  save(fileName, 'stats', '-append') ;
else
  save(fileName, 'stats') ;
end

% -------------------------------------------------------------------------
function [net, state, stats] = loadState(fileName)
% -------------------------------------------------------------------------
load(fileName, 'net', 'state', 'stats') ;
net = dagnn.DagNN.loadobj(net) ;
if isempty(whos('stats'))
  error('Epoch ''%s'' was only partially saved. Delete this file and try again.', ...
        fileName) ;
end

% -------------------------------------------------------------------------
function epoch = findLastCheckpoint(modelDir)
% -------------------------------------------------------------------------
list = dir(fullfile(modelDir, 'net-epoch-*.mat')) ;
tokens = regexp({list.name}, 'net-epoch-([\d]+).mat', 'tokens') ;
epoch = cellfun(@(x) sscanf(x{1}{1}, '%d'), tokens) ;
epoch = max([epoch 0]) ;

% -------------------------------------------------------------------------
function switchFigure(n)
% -------------------------------------------------------------------------
if get(0,'CurrentFigure') ~= n
  try
    set(0,'CurrentFigure',n) ;
  catch
    figure(n) ;
  end
end

% -------------------------------------------------------------------------
function clearMex()
% -------------------------------------------------------------------------
clear vl_tflow vl_imreadjpeg ;

% -------------------------------------------------------------------------
function prepareGPUs(opts, cold)
% -------------------------------------------------------------------------
numGpus = numel(opts.gpus) ;
if numGpus > 1
  % check parallel pool integrity as it could have timed out
  pool = gcp('nocreate') ;
  if ~isempty(pool) && pool.NumWorkers ~= numGpus
    delete(pool) ;
  end
  pool = gcp('nocreate') ;
  if isempty(pool)
    parpool('local', numGpus) ;
    cold = true ;
  end

end
if numGpus >= 1 && cold
  fprintf('%s: resetting GPU\n', mfilename)
  clearMex() ;
  if numGpus == 1
    gpuDevice(opts.gpus)
  else
    spmd
      clearMex() ;
      gpuDevice(opts.gpus(labindex))
    end
  end
end

author: XinLi-zn

repo: TADT

save_path: github-repos/MATLAB/XinLi-zn-TADT


text:
% Thanks for Carlos Lopez for this file and many other enhancements to the mp toolbox.
%quick and dirty tests
precision=1000;%notice that tol~=(.5)^precision;
a42=mp(rand(4,2),precision);%It will be nice to use mprand, but it doesn't work
b42=mp(rand(4,2),precision);

%2) direct/inverse pairs
c_mp=norm(asin(sin(a42))-a42);disp(['2.1 sin/asin' setstr(9) num2str(c_mp)])
c_mp=norm(acos(cos(a42))-a42);disp(['2.2 cos/acos' setstr(9) num2str(c_mp)])
c_mp=norm(atan(tan(a42))-a42);disp(['2.3 tan/atan' setstr(9) num2str(c_mp)])
c_mp=norm(asinh(sinh(a42))-a42);disp(['2.4 sinh/asinh' setstr(9) num2str(c_mp)])
c_mp=norm(acosh(cosh(a42))-a42);disp(['2.5 cosh/acosh' setstr(9) num2str(c_mp)])
c_mp=norm(atanh(tanh(a42))-a42);disp(['2.6 tanh/atanh' setstr(9) num2str(c_mp)])
c_mp=norm(exp(log(a42))-a42);disp(   ['2.7 exp/log   ' setstr(9) num2str(c_mp)])
c_mp=norm(sqrt(a42.^2)-a42);disp(    ['2.8 sqrt/^2   ' setstr(9) num2str(c_mp)])

%4) Catastrophic cancellation
z=sqrt(117)-sqrt(116);zExact=1/(sqrt(117)+sqrt(116));
disp('4.1 Improvement in the catastrophic cancellation of sqrt(117)-sqrt(116)')
disp([setstr(9) 'Standard double precision:' num2str(abs(z/zExact-1)) ])
for precision=300:100:1000
 zmp=sqrt(mp(117,precision))-sqrt(mp(116,precision));zmpExact=1/(sqrt(mp(117,precision))+sqrt(mp(116,precision)));
 disp([setstr(9) 'Multiple precision = ' num2str(precision) ':'  num2str(abs((zmp/zmpExact-1)))])
end

disp('4.2 Accuracy of tan(pi/4)-1')
for precision=300:100:1000
 zmp=mp('pi',precision)/4;
 disp([setstr(9) 'Multiple precision = ' num2str(precision) ':'  num2str(((tan(zmp)-1)))])
end
disp('4.3 Accuracy of 4*atan2(-1,1)+pi')
for precision=300:100:1000
 zmp=mp('pi',2*precision);%Request a higher precision pi in order to compare
 one=mp(1,precision);
 disp([setstr(9) 'Multiple precision = ' num2str(precision) ':'  ...
       num2str(((atan2(-one,one)*4+zmp)))])
end

%5) Test matrix; we will only test that the accuracy
%of an exact inverse times the matrix is closer to the identity in mp rather than in double
n=14;
disp(['5.1 Hilbert matrix and its exact inverse.'])
disp(['    norm in double: ' num2str([norm(hilb(n)*invhilb(n)-eye(n))])])
for precision=300:100:1000
 mp_set_defaults(precision);
 %5.1 Hilbert
 one=mp(1);

 J = 1:n;
 J = J(ones(n,1),:);
 I = J';
 E = mp(ones(n,n));
 mpHilb = E./(I+J-one);

 %Now the inverse; this is a slightly modified version of invhilb.m to produce mp values
 p = n;
 H = mp(zeros(n,n));
 for k = 1:n
  i=mp(k);
  if k > 1, p = ((n-i+one)*p*(n+i-one))/(k-one)^2; end
  r = p*p;
  H(k,k) = r/(2*i-one);
  for j = i+1:n
   %             r = -((n-j+one)*r*(n+j-one))/(j-one)^2;
   num = -((n-j+one)*r*(n+j-one));
   den= (j-one)^2;
   r = num/den;
   H(k,j) = r/(i+j-one);
   H(j,k) = r/(i+j-one);
  end
 end
 disp(['precision=' num2str(precision) ':' setstr(9) num2str(norm(H*mpHilb-eye(n)))]);
end
%check also some auxiliary functions (unrelated with precision!)
X = mp([2 8 4;7 3 9]);
%see "help min" for full documentation
if all(min(X,[],1) == [2 3 4]), disp('Test 6.1 passed'), else, disp('Test 6.1 failed'),end
if all(min(X,[],2) == [2;3]), disp('Test 6.2 passed'), else, disp('Test 6.2 failed'),end
if all(min(X,5)==[2 5 4;5 3 5]), disp('Test 6.3 passed'), else, disp('Test 6.3 failed'),end

if all(max(X,[],1)==[7 8 9]), disp('Test 6.4 passed'), else, disp('Test 6.4 failed'),end
if all(max(X,[],2)==[8;9]), disp('Test 6.5 passed'), else, disp('Test 6.5 failed'),end
if all(max(X,5)==[5 8 5;7 5 9]),disp('Test 6.6 passed'), else, disp('Test 6.6 failed'),end

X = mp([0 1 2;3 4 5]);
if all(sum(X,1) == [3 5 7]), disp('Test 6.7 passed'), else, disp('Test 6.7 failed'),end
if all(sum(X,2) == [ 3;12]), disp('Test 6.8 passed'), else, disp('Test 6.8 failed'),end

X = mp([3 7 5;0 4 2]);
if all(sort(X,1) == [0 4 2;3 7 5] ), disp('Test 6.9 passed'), else, disp('Test 6.9 failed'),end
if all(sort(X,2) == [3 5 7;0 2 4] ), disp('Test 6.10 passed'), else, disp('Test 6.10 failed'),end

author: opencobra

repo: cobratoolbox


text:
function [ STATEints ] = IDStateEpisode(STATEints,maxintdur,minSTATEdur)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
%
%INPUTS
%   STATEints
%   maxintdur   maximum interruption
%   minSTATEdur minimum state duration

%%
%SWSints = stateintervals{2};

[STATEints] = MergeSeparatedInts_ss(STATEints,maxintdur);
STATElengths = STATEints(:,2)-STATEints(:,1);
STATEints = STATEints(STATElengths>=minSTATEdur,:);

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function parts = strsplit_cerr(splitstr, str, option)
%STRSPLIT_CERR Split string into pieces.
%
%   STRSPLIT_CERR(SPLITSTR, STR, OPTION) splits the string STR at every occurrence
%   of SPLITSTR and returns the result as a cell array of strings.  By default,
%   SPLITSTR is not included in the output.
%
%   STRSPLIT_CERR(SPLITSTR, STR, OPTION) can be used to control how SPLITSTR is
%   included in the output.  If OPTION is 'include', SPLITSTR will be included
%   as a separate string.  If OPTION is 'append', SPLITSTR will be appended to
%   each output string, as if the input string was split at the position right
%   after the occurrence SPLITSTR.  If OPTION is 'omit', SPLITSTR will not be
%   included in the output.

%   Author:      Peter J. Acklam
%   Time-stamp:  2004-09-22 08:48:01 +0200
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam
%
% 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/>.

nargsin = nargin;
error(nargchk(2, 3, nargsin));
if nargsin < 3
	option = 'omit';
else
	option = lower(option);
end

splitlen = length(splitstr);
parts = {};

while 1

	k = strfind(str, splitstr);
	if isempty(k)
		parts{end+1} = str;
		break
	end

	switch option
		case 'include'
			parts(end+1:end+2) = {str(1:k(1)-1), splitstr};
		case 'append'
			parts{end+1} = str(1 : k(1)+splitlen-1);
		case 'omit'
			parts{end+1} = str(1 : k(1)-1);
		otherwise
			error(['Invalid option string -- ', option]);
	end

	str = str(k(1)+splitlen : end);

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%DEMO_NEURALNETCOV  Demonstration of Gaussian process with a neural
%                   network covariance function
%
%  Description
%    Infinite neural network solutions in 2D and 1D regression
%    problems with a comparison to Gaussian process solution given
%    by squared exponential covariance function. The noisy
%    observations y are assumed to satisfy
%
%         y = f + e,    where e ~ N(0, s^2)
%
%    where f is an unknown underlying function. A zero mean
%    Gaussian process prior is assumed for f
%
%         f ~ N(0, K),
%
%    where K is the covariance matrix whose elements are given by
%    neural network (or squared exponential) covariance function. A
%    prior is assumed for parameters of the covariance functions,
%    and the inference is done with a MAP estimate for parameter
%    values.
%
%    For more detailed discussion of infinite neural networks, see
%    e.g.
%
%      Neal, R. M. (1996). Bayesian Learning for Neural Networks.
%      Springer-Verlag.
%
%      Williams, C. K. I. (1996). Computing with infinite networks.
%      In Advances in Neural Information Processing Systems 9. MIT
%      Press, Cambridge, MA.
%
%
%  See also DEMO_REGRESSION1

% Copyright (c) 2010 Jaakko Riihim�ki, Aki Vehtari

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

% 2D REGRESSION DATA

% create 2D example data
x=rand(300,2)*2-1;
y=zeros(size(x,1),1); y(x(:,1)>0&x(:,2)>0)=1;
y=y+0.1*randn(size(y));

[n, nin] = size(x);

% --- Construct the model ---

% squared exponential covariance function
gpcf1 = gpcf_sexp('lengthScale', ones(1,nin), 'magnSigma2', 1);
% neural network covariance function
gpcf2 = gpcf_neuralnetwork('weightSigma2', ones(1,nin), 'biasSigma2', 1);
% Gaussian noise structures
lik = lik_gaussian('sigma2', 0.2^2);

% a prior structure for GP parameters
pt = prior_t('s2', 4);
gpcf1 = gpcf_sexp(gpcf1, 'lengthScale_prior', pt, 'magnSigma2_prior', pt);
gpcf2 = gpcf_neuralnetwork(gpcf2, 'weightSigma2_prior', pt, 'biasSigma2_prior', pt);

gp = gp_set('lik', lik, 'cf', gpcf1);
gp2 = gp_set('lik', lik, 'cf', gpcf2);

% --- MAP estimate using scaled conjugate gradient algorithm ---
%     (see scg for more details)

% Set the options for the optimization
opt=optimset('TolFun',1e-3,'TolX',1e-3);
% Optimize with the scaled conjugate gradient method
gp=gp_optim(gp,x,y,'opt',opt);
gp2=gp_optim(gp2,x,y,'opt',opt);

% create points where predictions are made
[xt1,xt2]=meshgrid(-1.5:0.05:1.5,-1.5:0.05:1.5);
xt=[xt1(:) xt2(:)];
% compute the predictions
[Eft_map, Varft_map] = gp_pred(gp, x, y, xt);
[Eft_map2, Varft_map2] = gp_pred(gp2, x, y, xt);

% Plot the predictions and data
figure, set(gcf, 'color', 'w')
mesh(xt1, xt2, reshape(Eft_map,size(xt1,1),size(xt1,2)));
hold on
plot3(x(:,1), x(:,2), y, '*')
axis on;
title('GP (squared exponential) predictions and the data points');

figure, set(gcf, 'color', 'w')
mesh(xt1, xt2, reshape(Eft_map2,size(xt1,1),size(xt1,2)));
hold on
plot3(x(:,1), x(:,2), y, '*')
axis on;
title('GP (neural network) predictions and the data points');

% 1D REGRESSION DATA

% create a 1D toy data
x=rand(100,1)*4-2;
y=norm_pdf(4*x)+0.05*randn(size(x));
[n, nin] = size(x);

gpcf1 = gpcf_sexp('lengthScale', ones(1,nin), 'magnSigma2', 1);
gpcf2 = gpcf_neuralnetwork('weightSigma2', ones(1,nin), 'biasSigma2', 1);

gpcf1 = gpcf_sexp(gpcf1, 'lengthScale_prior', pt, 'magnSigma2_prior', pt);
gpcf2 = gpcf_neuralnetwork(gpcf2, 'weightSigma2_prior', pt, 'biasSigma2_prior', pt);
lik = lik_gaussian();

gp = gp_set('lik', lik, 'cf', gpcf1);
gp2 = gp_set('lik', lik, 'cf', gpcf2);

% --- MAP estimate using scaled conjugate gradient algorithm ---
%     (see fminscg for more details)

% Set the options for the optimization
opt=optimset('TolFun',1e-3,'TolX',1e-3);
% Optimize with the scaled conjugate gradient method
gp=gp_optim(gp,x,y,'opt',opt);
gp2=gp_optim(gp2,x,y,'opt',opt);

% create points where predictions are made
xgrid=linspace(min(x)-1.5,max(x)+1.5,200)';
[Eft_map, Varft_map, lpyt_map, Eyt_map, Varyt_map] = gp_pred(gp, x, y, xgrid, 'yt', ones(200,1));
[Eft_map2, Varft_map2, lpyt_map2, Eyt_map2, Varyt_map2] = gp_pred(gp2, x, y, xgrid, 'yt', ones(200,1));

% Plot the predictions and data
color1=ones(1,3)*0.8; color2=ones(1,3)*0.5;
figure, set(gcf, 'color', 'w'), hold on
h1=fill([xgrid' fliplr(xgrid')], [(Eyt_map+1.96*sqrt(Varyt_map))' fliplr((Eyt_map-1.96*sqrt(Varyt_map))')], color1, 'edgecolor', color1);
% GP mean
h2=plot(xgrid, Eyt_map, 'color', color2, 'linewidth', 3);
% observations
h3=plot(x, y, 'xk', 'markersize', 10, 'linewidth', 2);
% true function
h4=plot(xgrid, norm_pdf(4*xgrid), 'color', 'r', 'linewidth', 2);
legend([h1 h2 h3 h4], 'GP 95% CI', 'GP mean', 'observations', 'true latent function')
title('GP (squared exponential) predictions and the data points');

figure, set(gcf, 'color', 'w'), hold on
h1=fill([xgrid' fliplr(xgrid')], [(Eyt_map2+1.96*sqrt(Varyt_map2))' fliplr((Eyt_map2-1.96*sqrt(Varyt_map2))')], color1, 'edgecolor', color1);
% GP mean
h2=plot(xgrid, Eyt_map2, 'color', color2, 'linewidth', 3);
% observations
h3=plot(x, y, 'xk', 'markersize', 10, 'linewidth', 2);
% true function
h4=plot(xgrid, norm_pdf(4*xgrid), 'color', 'r', 'linewidth', 2);
legend([h1 h2 h3 h4], 'GP 95% CI', 'GP mean', 'observations', 'true latent function')
title('GP (neural network) predictions and the data points');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [dat, ref] = ft_preproc_rereference(dat, refchan, method, handlenan)

% FT_PREPROC_REREFERENCE computes the average reference over all EEG channels
% or rereferences the data to the selected channels
%
% Use as
%   [dat] = ft_preproc_rereference(dat, refchan, method, handlenan)
% where
%   dat        data matrix (Nchans X Ntime)
%   refchan    vector with indices of the new reference channels, or 'all'
%   method     string, can be 'avg' or 'median'
%   handlenan  boolean, can be true or false
%
% If the new reference channel is not specified, the data will be
% rereferenced to the average of all channels.
%
% See also PREPROC

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

% determine the size of the data
[Nchans, Nsamples] = size(dat);

% determine the new reference channels
if nargin<2 || isempty(refchan) || (ischar(refchan) && strcmp(refchan, 'all'))
  refchan = 1:Nchans;
end

if nargin<3 || isempty(method)
  method = 'avg';
end

if nargin<4 || isempty(handlenan)
  handlenan = false;
end

hasnan = any(any(isnan(dat(refchan,:))));

if hasnan && handlenan
  % preprocessing works differently if channels contain NaN
  switch method
    case 'avg'
      ref = nanmean(dat(refchan,:), 1);
    case 'median'
      ref = nanmedian(dat(refchan,:), 1);
    otherwise
      ft_error('unsupported method')
  end % switch
else
  % preprocessing fails on channels that contain NaN
  if any(isnan(dat(:)))
    ft_warning('FieldTrip:dataContainsNaN', 'data contains NaN values');
  end
  % compute the average value over the reference channels
  switch method
    case 'avg'
      ref = mean(dat(refchan,:), 1);
    case 'median'
      ref = median(dat(refchan,:), 1);
    otherwise
      ft_error('unsupported method')
  end % switch
end

% subtract the new reference from the data
for chan=1:Nchans
  dat(chan,:) = dat(chan,:) - ref;
end

author: spm


text:
% function [IDX_SC, Uk, Dk, time_SC] = SC(G,lap_type)
%
% This is an implementation of the classical Spectral CLustering algorithm
% using either:
% - the combinatorial Laplacian if lap_type='combinatorial' (in this case,
% the features are not normalized);
% - or the normalized Laplacian if lap_type='normalized' (in this case,
% the features are normalized).
% See the paper cited below for many references on this algorithm.
%
% Copyright (C) 2016 Nicolas Tremblay, Gilles Puy.
% This file is part of the CSCbox (Compressive Spectral Clustering toolbox)
%
% The CSCbox 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.
%
% The CSCbox 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 you use this toolbox please kindly cite
%     N. Tremblay, G. Puy, R. Gribonval and P. Vandergheynst.
%     Compressive Spectral Clustering.
%     ArXiv e-prints:1602.02018 Feb. 2016.

function [IDX_SC, Uk, Dk, time_SC] = SC(G,lap_type)

normBOOL=strcmp(lap_type,'normalized');
G = gsp_create_laplacian(G,lap_type);

% eigs
if isfield(G, 'G.Dk'), % if it was already calculated, do not do it again
    Dk=G.Dk;
    Uk=G.Uk;
    time_SC.eigs=G.time_SC_eigs;
else % run eigs
    tic;
    opts.isreal=1;opts.issym=1;opts.maxit=10000;
    [Uk,Dk]=eigs(G.L,G.k,'SA',opts);Dk=diag(Dk);
    time_SC.eigs=toc;
end

if normBOOL % kmeans on Yk
    tic;
    Yk=Uk./repmat(sqrt(sum(Uk.^2,2)),1,G.k);
    IDX_SC = kmeans(Yk, G.k,'Replicates',20);
    time_SC.kmeans=toc;
else % kmeans on Uk
    tic;
    IDX_SC = kmeans(Uk, G.k,'Replicates',20);
    time_SC.kmeans=toc;
end

time_SC.total=time_SC.eigs+time_SC.kmeans;

repo: gspbox


text:
root_10 = 'ACS_10_5YR_DP02';
ACS_path = [root_10, '/', root_10, '_metadata.csv'];
[~, text_10] = xlsread(ACS_path);

root_11 = 'ACS_11_5YR_DP02';
ACS_path = [root_11, '/', root_11, '_metadata.csv'];
[~, text_11] = xlsread(ACS_path);

root_12 = 'ACS_12_5YR_DP02';
ACS_path = [root_12, '/', root_12, '_metadata.csv'];
[~, text_12] = xlsread(ACS_path);

root_13 = 'ACS_13_5YR_DP02';
ACS_path = [root_13, '/', root_13, '_metadata.csv'];
[~, text_13] = xlsread(ACS_path);

root_14 = 'ACS_14_5YR_DP02';
ACS_path = [root_14, '/', root_14, '_metadata.csv'];
[~, text_14] = xlsread(ACS_path);

root_15 = 'ACS_15_5YR_DP02';
ACS_path = [root_15, '/', root_15, '_metadata.csv'];
[~, text_15] = xlsread(ACS_path);

root_16 = 'ACS_11_5YR_DP02';
ACS_path = [root_16, '/', root_16, '_metadata.csv'];
[~, text_16] = xlsread(ACS_path);

repo: MCM-ICM


text:
% Test file for bndfun/diff.m

function pass = test_diff(pref)

% Get preferences.
if ( nargin < 1 )
    pref = chebfunpref();
end

% Set the domain.
dom = [-2 7];

% Generate a few random points to use as test values.
seedRNG(6178);
x = diff(dom) * rand(100, 1) + dom(1);

%%
% Spot-check derivatives for a couple of functions.

f = bndfun(@(x) exp(x/10) - x, struct('domain', dom), pref);
df = diff(f);
df_exact = @(x) exp(x/10)./10 - 1;
err = df_exact(x) - feval(df, x);
pass(1) = (norm(err, inf) < 1e3*get(f, 'vscale')*eps);

f = bndfun(@(x) atan(x), struct('domain', dom), pref);
df = diff(f);
df_exact = @(x) 1./(1 + x.^2);
err = df_exact(x) - feval(df, x);
pass(2) = (norm(err, inf) < 1e3*get(f, 'vscale')*eps);

f = bndfun(@(x) sin(x), struct('domain', dom), pref);
df = diff(f);
df_exact = @(x) cos(x);
err = df_exact(x) - feval(df, x);
pass(3) = (norm(err, inf) < 1e3*get(f, 'vscale')*eps);

z = exp(2*pi*1i/3);
f = bndfun(@(t) airy(z*t), struct('domain', dom), pref);
df = diff(f);
df_exact = @(t) z*airy(1, z*t);
err = df_exact(x) - feval(df, x);
pass(4) = (norm(err, inf) < 1e3*get(f, 'vscale')*eps);

%%
% Verify that calling diff() gives the same answer as direct construction.

f = bndfun(@(x) 0.5*x - 0.0625*sin(8*x), struct('domain', dom), pref);
df = bndfun(@(x) sin(4*x).^2, struct('domain', dom), pref);
err = diff(f) - df;
pass(5) = (get(err, 'vscale') < 1e4*get(f, 'vscale')*eps);

%%
% Verify basic differentiation rules.

f = bndfun(@(x) x.*sin(x.^2) - 1, struct('domain', dom), pref);
df = diff(f);
g = bndfun(@(x) exp(-x.^2), struct('domain', dom), pref);
dg = diff(g);

tol_f= 10*get(f, 'vscale')*eps;
tol_df = 10*get(df, 'vscale')*eps;
tol_g= 10*get(g, 'vscale')*eps;
tol_dg = 10*get(dg, 'vscale')*eps;

errfn = diff(f + g) - (df + dg);
err = feval(errfn, x);
pass(6) = (norm(err, inf) < max([tol_f ; tol_g ; tol_df ; tol_dg]));

errfn = diff(f.*g) - (f.*dg + g.*df);
err = feval(errfn, x);
pass(7) = (norm(err, inf) < 1e1*max([tol_f ; tol_g ; tol_df ; tol_dg]));

const = bndfun(@(x) ones(size(x)), struct('domain', dom), pref);
dconst = diff(const);
err = feval(dconst, x);
pass(8) = (norm(err, inf) <= 10*get(dconst, 'vscale')*eps);

%%
% Check higher-order derivatives.

f = bndfun(@(x) x.*atan(x) - x - 0.5*log(1 + x.^2), struct('domain', dom), ...
    pref);
df2 = diff(f, 2);
df2_exact = @(x) 1./(1 + x.^2);
err = df2_exact(x) - feval(df2, x);
pass(9) = (norm(err, inf) < 1e7*get(df2, 'vscale')^2*eps);

f = bndfun(@(x) sin(x), struct('domain', dom), pref);
df4 = diff(f, 4);
df4_exact = @(x) sin(x);
err = norm(df4_exact(x) - feval(df4, x), inf);
tol = 10*get(df4, 'vscale')*eps;
pass(10) = err < 1e6*tol;

f = bndfun(@(x) x.^5 + 3*x.^3 - 2*x.^2 + 4, struct('domain', dom), pref);
df6 = diff(f, 6);
df6_exact = @(x) zeros(size(x));
err = df6_exact(x) - feval(df6, x);
pass(11) = (norm(err, inf) <= get(df6, 'vscale')^6*eps);

%%
% Check operation for array-valued bndfun objects.

f = bndfun(@(x) [sin(x) x.^2 exp(1i*x)], struct('domain', dom), pref);
df_exact = @(x) [cos(x) 2*x 1i*exp(1i*x)];
err = feval(diff(f), x) - df_exact(x);
pass(12) = (norm(err(:), inf) < 1e3*max(get(f, 'vscale')*eps));

% DIM option.
dim2df = diff(f, 1, 2);
g = @(x) [(x.^2 - sin(x)) (exp(1i*x) - x.^2)];
err = feval(dim2df, x) - g(x);
pass(13) = (norm(err(:), inf) < 10*max(get(f, 'vscale')*eps));

dim2df2 = diff(f, 2, 2);
g = @(x) exp(1i*x) - 2*x.^2 + sin(x);
err = feval(dim2df2, x) - g(x);
pass(14) = (norm(err(:), inf) < 10*max(get(f, 'vscale')*eps));

% DIM option should return an empty bndfun for non-array-valued input.
f = bndfun(@(x) x.^3, struct('domain', dom), pref);
dim2df = diff(f, 1, 2);
pass(15) = (isempty(dim2df));

%% Test on singular function:

pow = -0.5;
op = @(x) (x - dom(1)).^pow.*sin(x);
pref.blowup = true;
data.domain = dom;
data.exponents = [pow 0];
f = bndfun(op, data, pref);
df = diff(f);
vals_df = feval(df, x);
df_exact = @(x) (x - dom(1)).^(pow-1).*(pow*sin(x)+(x - dom(1)).*cos(x));
vals_exact = feval(df_exact, x);
err = vals_df - vals_exact;
pass(16) = ( norm(err, inf) < 1e5*eps*norm(vals_exact, inf) );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function PlotFieldonMesh(coordinates,nodes,component)
%--------------------------------------------------------------------------
% Code written by : Siva Srinivas Kolukula                                |
%                   Senior Research Fellow                                |
%                   Structural Mechanics Laboratory                       |
%                   Indira Gandhi Center for Atomic Research              |
%                   India                                                 |
% E-mail : allwayzitzme@gmail.com                                         |
%          http://sites.google.com/site/kolukulasivasrinivas/             |
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% Purpose:
%         To plot the profile of a component on mesh
% Synopsis :
%           ProfileonMesh(coordinates,nodes,component)
% Variable Description:
%           coordinates - The nodal coordinates of the mesh
%           -----> coordinates = [node X Y ]
%           nodes - The nodal connectivity of the elements
%           -----> nodes = [node1 node2......]
%           component - The components whose profile to be plotted
%           -----> components = a column vector in the order of node
%                               numbers
%--------------------------------------------------------------------------

nel = length(nodes) ;                  % number of elements
nnode = length(coordinates) ;          % total number of nodes in system
nnel = size(nodes,2);                % number of nodes per element
%
% Initialization of the required matrices
X = zeros(nnel,nel) ;
Y = zeros(nnel,nel) ;
Z = zeros(nnel,nel) ;
profile = zeros(nnel,nel) ;
%
for iel=1:nel
     for i=1:nnel
     nd(i)=nodes(iel,i);         % extract connected node for (iel)-th element
     X(i,iel)=coordinates(nd(i),1);    % extract x value of the node
     Y(i,iel)=coordinates(nd(i),2);    % extract y value of the node
     end
     profile(:,iel) = component(nd') ;         % extract component value of the node
end

% Plotting the FEM mesh and profile of the given component
     f3 = figure ;
     set(f3,'name','Postprocessing','numbertitle','off') ;
     plot(X,Y,'k')
     fill(X,Y,profile)
     axis off ;
     % Colorbar Setting
     SetColorbar
 end

author: Sable

repo: mcbench-benchmarks


text:
function [w,mu,P]=EMAlgGaussClust(z,w,mu,P,numIter)
%%EMALGGAUSSCLUST Use the expectation maximization (EM) algorithm to
%                 refine estimates of the components of a Gaussian mixture
%                 with a known number of terms given a set of samples.
%
%INPUTS: z A zDim X numPoints set of samples of the Gaussian mixture.
%        w A KX1 set of initial weight estimates of the K Gaussians in the
%          mixture.
%       mu A zDim X K set of initial mean estimates of the component
%          Gaussians in the mixture.
%        P A zDim X zDim X K hypermatrix of initial covariance matrix
%          estimates of the components of the Gaussians in the mixture.
%  numIter The number of iterations of the EM algorithm to perform.
%
%OUTPUTS: w The refined weights.
%        mu The refined means.
%         P The refined covariance matrix estimates.
%
%The EM algorithm for Gaussian mixtures is an implementation of the
%algorithm described in Chapter 9.2.2 of [1].
%
%REFERENCES:
%[1] C. M. Bishop, Pattern Recognition and Machine Learning. Cambridge,
%    United Kingdom: Springer, 2007.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    %zDim=size(z,1);
    numPoints=size(z,2);
    K=size(mu,2);

    gamma=zeros(numPoints,K);
    for curIter=1:numIter
        %Calculate the posterior weights.
        for k=1:K
            gamma(:,k)=GaussianPDF(z,mu(:,k),P(:,:,k))*w(k);
        end
        %Normalize the weights for each measurement.
        gamma=bsxfun(@rdivide,gamma,sum(gamma,2));

        %Update the means, covariances and weights using the posterior
        %weights.
        for k=1:K
            Nk=sum(gamma(:,k));
            w(k)=Nk/numPoints;

            [mu(:,k), P(:,:,k)]=calcMixtureMoments(z,gamma(:,k)/Nk);
        end
    end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function quad_mesh_test05 ( )

%*****************************************************************************80
%
%% QUAD_MESH_TEST05 tests BOUNDARY_EDGE_COUNT_EULER_Q4_MESH.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'QUAD_MESH_TEST05\n' );
  fprintf ( 1, '  BOUNDARY_EDGE_COUNT_EULER_Q4_MESH counts the\n' );
  fprintf ( 1, '    boundary edges using Euler''s formula.\n' );

  [ node_num, element_num, hole_num ] = example1_q4_mesh_size ( );

  boundary_edge_num = boundary_edge_count_euler_q4_mesh ( node_num, ...
    element_num, hole_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of boundary edges = %d\n', boundary_edge_num );
  fprintf ( 1, '  Correct number =           %d\n', 22 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function blas1_d_test05 ( )

%*****************************************************************************80
%
%% TEST05 demonstrates DMACH.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 May 2006
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST05\n' );
  fprintf ( 1, '  DMACH returns some approximate machine numbers.\n' );
  fprintf ( 1, '\n' );
  job = 1;
  fprintf ( 1, '  DMACH(1) = EPS =  %e\n', dmach ( job ) );
  job = 2;
  fprintf ( 1, '  DMACH(2) = TINY = %e\n', dmach ( job ) );
  job = 3;
  fprintf ( 1, '  DMACH(3) = HUGE = %e\n', dmach ( job ) );

  return
end


text:
function n = dim(A)
%DIM          Dimension of a square matrix
%
%    n = dim(A)
%

% written  10/16/98     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  if A.complex
    [m n] = size(A.mid);
  else
    [m n] = size(A.inf);
  end

  if m ~= n
    error('function dim called with non-square matrix')
  end;

repo: OpenMAS


text:
%%vis3d(img, mycmap): a GUI-enabled visualization of the 3D volume img,
%%with an optional colormap argument.

%{
Joshua Stough
JHU, iacl
June 2012
April 2013

BSD Licence:
Copyright (c) 2013, Joshua Stough
All rights reserved.

Please comment on MATLAB fileexchange if you find this useful. Thank you.

This program offers a visualization of 3D volumes.

Obviously, this is not a new idea. This program might offer over other
programs:
    -change slice via click drag or scroll
    -orthogonal planes are visualized on each projection in a color-coded
    manner, allowing a better sense of where you are in the 3D space.
    -windowing capability
    -changing the colormap
    -neurological coordinate frame

However, after seeing the slice function, I think I must be missing a lot
of built-in functionality, making this program pointless. But I share in any
event as I obviously cannot see all ends (nerd).

References to really good and thoughtful 3D viewer implementations in
matlab (including buttons and all axis flipping that I can't do):
http://www.mathworks.com/matlabcentral/fileexchange/2256-orthoview
http://www.mathworks.com/matlabcentral/fileexchange/764-sliceomatic

For some 3d matrix of double or uint8, call
>> vis3d(img);
or, if you have some multilabel segmentation image, try:
>> vis3d(seg, 'multi');
or, if you want define your own colormap, you can send that instead:
>> mymap = colormap('lines');
>> vis3d(seg, mymap);
lastly, you can customize the figure title:
>> vis3d(seg, '', 'My Figure Title');

Updates/Modifications:
4/13:
-fixed slice number issue: the out-of-plane slice number was incorrect
in the initial window.
-changed interaction: now, scrolling changes slice number
-fixed the 3d volume visualization errors arising from the treatment of
images versus matrices in matlab.
-fixed tick mark misrepresentation.
-added windowing on the colormap
-added choice of colormap in addition to typical maps.
%}

function [] = vis3d(seg, mycmap, figName)

    %We check for the colormap after we have a figure open.
    if nargin < 3
        figName = 'vis3d: click and scroll/drag vertically (up=closer to eye).';
        userSpecFigName = false; %did the user specify the figure name.
    else
        userSpecFigName = true;
    end
    %arg 2 taken care of below.

    if ~(isa(seg, 'double') || isa(seg, 'uint8'))
        fprintf('The image is not double or uint8.  You may receive many warnings.\n');
        fprintf('Cast your image to either before sending here.');
    end
   %seg = seg/max(seg(:)); %Make 0-1.
    segRange = [min(seg(:)) max(seg(:))];
    origSegRange = segRange;
    %segRange is required to ensure that if you view a slice with more or
    %less range, you get colors consistent with other slices.

    sizeSeg = size(seg);
    sn = round(sizeSeg/2);
    %slice numbers for the three projections.

    titleLines = {'YX view: ', 'YZ view: ', 'XZ view: '};
    %titleColors = {'b', 'r', 'g'};
    titleColors = {[.3 .3 .8], [.8 .3 .3], [.3 .8 .3]};

    %The figure handle, with its callbacks.
    mainHandle = figure('Position', [50 50 800 800], ...
        'WindowButtonDownFcn', @buttonDownCallback, ...
        'WindowButtonUpFcn', @buttonUpCallback, ...
        'Name', figName); hold on;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Adding colormap popupmenu
    colormapText = uicontrol('Style', 'text', 'Position', [300 0 100, 20], ...
        'FontSize', 12, 'FontWeight', 'bold');
    set(colormapText, 'String', 'Colormap:', 'ForegroundColor', [1 1 1]);

    colormapPopup = uicontrol('Style', 'popupmenu', ...
        'Position', [400 0 200 20], ...
        'Callback', @specifyColormap);
    set(colormapPopup, 'String', {'gray'; 'multi-label'; 'jet'}, ...
        'Value', 1, 'FontSize', 12, 'FontWeight', 'bold');

    %define the colormaps.
    %A colormap I find useful for multi-label segmentation images, which is
    %why I made this thing to begin with (and why the image is called seg).
    multi = [0 0 0;
            1 0 1;
            0 .7 1;
            1 0 0;
            .3 1 0;
            0 0 1;
            1 1 1;
            1 .7 0];
    graymap = colormap('gray');
    jetmap = colormap('jet');
    the_cmaps = {graymap; multi; jetmap};

    %Is there a user-defined map?
    %Now that we're sure to have a figure open, get a colormap.
    if nargin < 2 || strcmp(mycmap, '')
        mycmap = gray;
    elseif strcmp(mycmap, 'multi')
        mycmap = multi;
        set(colormapPopup, 'Value', 2);
    else
        %User defined a colormap to use.  Set it as current.
        the_cmaps{end+1} =  mycmap;
        set(colormapPopup, 'String', {'gray'; 'multi-label'; 'jet'; 'user-spec'}, ...
        'Value', 4);
    end
    %End of colormap popup menu.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %flip the matrix to better organize views. This is memory-stupid!!!
    %Please help!
    xyview = seg;
    xzview = seg;
    yzview = seg;

    xyview = flipdim(xyview, 1);

    xzview = flipdim(xzview, 1);
    xzview = permute(xzview, [3 2 1]);
    %xzview = flipdim(xzview, 2);

    yzview = permute(yzview, [1 3 2]);
    yzview = flipdim(yzview, 1);

    %Now, we plot the three orthogonal slices and try to get the two lines
    %on each slice, which represent the other planes.
    %Notice we are using axis xy, so that the tick marks are correct, that
    %also means we don't need to flip z in yz and xz views...
    handles{1} = subplot(2,2,3);
    imHandles{1} = imagesc(squeeze(xyview(:,:,sn(3))), segRange); colormap(mycmap); axis image; axis xy;
    %set(handles{1}, 'LooseInset', get(gca,'TightInset')) %made smaller...
    titleHandles{1} = title(handles{1}, 'XY view', 'Color', titleColors{1}, 'FontSize', 14, 'FontWeight', 'bold');

    handles{2} = subplot(2,2,1);
    imHandles{2} = imagesc(squeeze(yzview(sn(1),:,:)), segRange); colormap(mycmap); axis image; axis xy;
    titleHandles{2} = title(handles{2}, 'YZ view', 'Color', titleColors{2}, 'FontSize', 14, 'FontWeight', 'bold');
    %set(handles{2}, 'XDir', 'rev');
    %I used reverse here to have appropriate neurological coord, but then I
    %thought, people like to see 'through the eyes of the patient' in
    %coronal (i.e., patient l-r equals l-r in the image).

    handles{3} = subplot(2,2,2);
    imHandles{3} = imagesc(squeeze(xzview(:,sn(2),:)), segRange); colormap(mycmap); axis image; axis xy;
    titleHandles{3} = title(handles{3}, 'XZ view', 'Color', titleColors{3}, 'FontSize', 14, 'FontWeight', 'demi');

    %Making sure the titles all have the initial slice number as well.
    offsliceCoord = [3 1 2];
    for i = 1:3
        set(titleHandles{i}, 'String', sprintf('%s%03d', titleLines{i}, sn(offsliceCoord(i))));
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %lower and upper bound on the intensity window and text box.

    %Adding uicontrol elements would remove the toolbar from the figure
    %without this statement.
    set(mainHandle, 'Toolbar', 'figure');

    IntWinSliderHandle{1} = uicontrol('Style', 'slider', ...
        'Position', [0 0 200 20], ...
        'Callback', @IntWinBound);
    set(IntWinSliderHandle{1}, 'Value', segRange(1));

    IntWinSliderHandle{2} = uicontrol('Style', 'slider', ...
        'Position', [0 25 200 20], ...
        'Callback', @IntWinBound);
    set(IntWinSliderHandle{2}, 'Value', segRange(2));

    for i = 1:2
        set(IntWinSliderHandle{i}, 'Min', segRange(1));
        set(IntWinSliderHandle{i}, 'Max', segRange(2));
    end

    IntWinTextBound{1} = uicontrol('Style', 'edit', ...
        'Position', [210 0 60 20], ...
        'Callback', @defineBound);
    set(IntWinTextBound{1}, 'String', segRange(1));

    IntWinTextBound{2} = uicontrol('Style', 'edit', ...
        'Position', [210 25 60 20], ...
        'Callback', @defineBound);
    set(IntWinTextBound{2}, 'String', segRange(2));

    IntWinText = uicontrol('Style', 'text', 'Position', [0 50 210 15], ...
        'FontSize', 12, 'FontWeight', 'bold');

    set(IntWinText, 'String', 'Windowing', ...
            'ForegroundColor', [1 1 1]);
    %End of intensity windowing elements.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %And the 3D view.
    handles{4} = subplot(2,2,4);

    %We're using the slice visualization function here, but it's screwy
    %because the x and y are transposed and axes must be flipped for the
    %appropriate perspective,

    %I wanted the callback to avoid duplicating code, but the camera position
    %is hard to get before slice is called...buttonUpCallback([]); so here
    %is init.
    %due to confusing coordinates for the slice function, we need to
    %permute and flip dimensions...
    visSeg = seg;
    visSeg = permute(visSeg, [2 1 3]);
    %visSeg = flipdim(visSeg, 1);
    %visSeg = flipdim(visSeg, 2);
    flippedx = inline(sprintf('%d - x + 1', size(visSeg, 2)));
    %The 2nd dimension of the matrix ends up being 'called' x by matlab.
    %wtf?

    hslc=slice(visSeg, flippedx(sn(1)) ,sn(2),sn(3));
    axis equal; axis vis3d; set(hslc(1:3),'LineStyle','none');
    xlabel 'p-a' ;ylabel 'l-r' ;zlabel 'i-s';
    %set(gca, 'XDir', 'reverse');
    %The above XDir didn't work, without screwing up the data.  So, I'll
    %just change the labels...WAIT!: that's wrong because the slice numbers
    %become wrong then.  Have to flip...
    myx = get(gca, 'XTickLabel');
    set(gca, 'XTickLabel', flippedx(str2num(myx)));

    title 'Manually enable/disable Rotate3D.'
    %To maintain the 3D perspective throughout.
    camPos3D = get(handles{4}, 'CameraPosition');

    %Too many problems trying to regain control from the rotate3d.
    %rotHandle = rotate3d(handles{4}); so it's up to the user.
    rotHandle = rotate3d;
    setAllowAxesRotate(rotHandle, handles{1}, false);
    setAllowAxesRotate(rotHandle, handles{2}, false);
    setAllowAxesRotate(rotHandle, handles{3}, false);
    %End of the 3D view elements
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Generate the lines representing the orthogonal planes.  This part is
    %so confusing, in part because line will use the bottom
    %left, whereas we're thinking of the top left (image) as being 0,0.
    %But it's mostly my fault, because coordinate systems are so HARD...

    lines{1}.x = [1 sizeSeg(2); sn(2) sn(2)]';
    lines{1}.y = [sn(1) sn(1); 1 sizeSeg(1)]';
    %the x coords of the horiz and vertical line respectively.
    %the y coords of the horiz and vertical line respectively.

    lines{2}.x = [1 sizeSeg(2); sn(2) sn(2)]';
    lines{2}.y = [sn(3) sn(3); 1 sizeSeg(3)]';

    lines{3}.x = [1 sizeSeg(1); sn(1) sn(1)]';
    lines{3}.y = [sn(3) sn(3); 1 sizeSeg(3)]';

    %We now instantiate the lines.  We also want the line colors to reflect
    %the plane they represent, which is given by the titleColors. The trick
    %is weird though: when we're visualizing the XY plane, the two lines on
    %it are the XZ plane (vertical), and the YZ plane, horizontal.  whew.
    %for each, first is horizontal, next is vertical.
    lineHandles = {};
    subplot(handles{1});
    lineHandles{1}(1) = line(lines{1}.x(:,1), lines{1}.y(:,1), 'Color', titleColors{2}, 'LineWidth', 1);
    lineHandles{1}(2) = line(lines{1}.x(:,2), lines{1}.y(:,2), 'Color', titleColors{3}, 'LineWidth', 1);

    subplot(handles{2});
    lineHandles{2}(1) = line(lines{2}.x(:,1), lines{2}.y(:,1), 'Color', titleColors{1}, 'LineWidth', 1);
    lineHandles{2}(2) = line(lines{2}.x(:,2), lines{2}.y(:,2), 'Color', titleColors{3}, 'LineWidth', 1);

    subplot(handles{3});
    lineHandles{3}(1) = line(lines{3}.x(:,1), lines{3}.y(:,1), 'Color', titleColors{1}, 'LineWidth', 1);
    lineHandles{3}(2) = line(lines{3}.x(:,2), lines{3}.y(:,2), 'Color', titleColors{2}, 'LineWidth', 1);
    %End of the line crap.  This is really confusing.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %Variables to keep track of mouse motion and which slice is being
    %maniupulated.
    lastPoint = [];
    whichView = 1;

    %At this point, time to define our callbacks.

    %This is the callback when the user clicks.  We want to determine which
    %if any projection they clicked on so that further vertical motion can
    %be associated with changing the slice for that projection.
    function buttonDownCallback(varargin)

        p = get(gca, 'CurrentPoint');
        lastPoint = [p(1); p(3)];

        %We know where they clicked with respect to the last active axes,
        %determine which axes that was.
        whichView = find(cell2mat(handles) == gca);

        if whichView == 4
            %Can't do anthing with the 3D view, unless it's rotate
            %selection.  The user has to do it themselves.
            %set(rotHandle, 'Enable', 'on', 'ActionPostCallback', @buttonUpOnRotateCallback);
            return
        end

        %By getting the camera position on a click, we can maintain the 3d
        %perspective in the case of changing slices.
        camPos3D = get(handles{4}, 'CameraPosition');

        %Now determine if they clicked on a valid spot, within the coords
        %of the chosen image/axes.
        xlim = get(gca, 'XLim');
        ylim = get(gca, 'YLim');
        if lastPoint(1) >= xlim(1) && lastPoint(1) <= xlim(2) && ...
                lastPoint(2) >= ylim(1) && lastPoint(2) <= ylim(2)
            %They clicked on a valid point, time to activate the motion
            %callback. Also let's tell them the intensity at the point they
            %click. But we also don't want to remove what they decided to
            %call the main figure.
            t_img = get(imHandles{whichView}, 'CData');
            p = round(p);
            if userSpecFigName
                set(mainHandle, 'Name', sprintf('%s: Pos: %d %d, Int: %f', ...
                    figName, p(3), p(1), t_img(p(3), p(1))));
            else
                set(mainHandle, 'Name', sprintf('Pos: %d %d, Int: %f', ...
                    p(3), p(1), t_img(p(3), p(1))));
            end
            set(mainHandle, 'WindowButtonMotionFcn', @dragCallback);
            set(mainHandle, 'WindowScrollWheelFcn', @scrollCallback);
            %fprintf('You clicked on %f, %f on gca %f (%d).\n', p(1), p(3), gca, whichView);
        %else
            %fprintf('Oops, you did''t click flippedx(sn(1))a valid point.\n');
        end
    end

    %At this point, it's established that the user clicked on one of the
    %images and is currently dragging.  We will change the image on the
    %relevant slice and update slice info (in the title).
    %And then, if this works, I'll try the red lines on the other two
    %projections.
    function dragCallback(varargin)
        p = get(gca, 'CurrentPoint');
        motionV = round([p(1); p(3)] - lastPoint);
        if abs(motionV(2)) < 2
            return
        end %Wait until the motion is noticeable.

        lastPoint = [p(1); p(3)];

        newslice = sn(offsliceCoord(whichView)) + round(motionV(2)/2);

        updateSliceViz(newslice);

    end

    function buttonUpCallback(varargin)
        %Nullify the motion and correct the 3D view.
        %Slice is too easy though, thanks Matlab.  I should be ashamed
        %http://www.mathworks.com/matlabcentral/fileexchange/2256-orthoview
        set(mainHandle, 'WindowButtonMotionFcn', '');

        axes(handles{4}); hslc=slice(visSeg, flippedx(sn(1)), sn(2), sn(3));%rotate3d on;
        axis equal; axis vis3d; set(hslc(1:3),'LineStyle','none');
        xlabel 'p-a' ;ylabel 'l-r' ;zlabel 'i-s';

        %Same stupid and incorrect label issue.
        myx = get(gca, 'XTickLabel');
        set(gca, 'XTickLabel', flippedx(str2num(myx)));

        title 'Manually enable/disable Rotate3D.'
        set(handles{4}, 'CameraPosition', camPos3D);
    end

    %I had problems with 'ActionPostCallback' and so I'm leaving rotation
    %of the 3d plot up to the user clicking on and off the button.
    %function buttonUpOnRotateCallback(varargin)
    %    set(rotHandle, 'ActionPostCallback', '');
    %    set(rotHandle, 'Enable', 'off');
    %end

    function scrollCallback(src,evnt)
        %Change whichView's slice by the vertical motion.
        offsliceCoord = [3 1 2];
        %if they're moving xy slice, that's a change in z, etc.

        if whichView == 4
            return
        end

        if evnt.VerticalScrollCount > 0
            newslice = sn(offsliceCoord(whichView)) - 1;
        else
            newslice = sn(offsliceCoord(whichView)) + 1;
        end

        updateSliceViz(newslice);

    end

    function updateSliceViz(newslice)
        %Change whichView's slice by the vertical motion.

        offsliceCoord = [3 1 2];
        %if they're moving xy slice, that's a change in z, etc.

        if newslice > 0 && newslice <= sizeSeg(offsliceCoord(whichView))
            sn(offsliceCoord(whichView)) = newslice;
        end
        subplot(handles{whichView});

        %This is a tough to comprehend part.  I want to move the lines
        %representing the current plane in the other orthogonal views.
        %I'll be changing line properties accordingly.
        if whichView == 1
            %XY view (axial)
            %imagesc(squeeze(seg(:,:,sn(3))), segRange); axis image; hold on
            %While you don't have to respecify the colormap, you do need to
            %remind imagesc of the appropriate range and ensure the axes
            %are good. In fact, I forgot that with matlab, everything is an
            %object.  Just change the data in the object.
            set(imHandles{1}, 'CData', squeeze(xyview(:,:,sn(3))));

            %Moving in z, so change the horizontal lines in the YZ and XZ
            %plots.
            set(lineHandles{2}(1), 'YData', [sn(3) sn(3)]');
            set(lineHandles{3}(1), 'YData', [sn(3) sn(3)]');

        elseif whichView == 2
            %YZ view (coronal)
            set(imHandles{2}, 'CData', squeeze(yzview(sn(1),:,:)));
            %imagesc(squeeze(seg(sn(1),:,:)), segRange); axis image

            %Moving in x, so change the horizontal line in the XY
            %(representing front-back) and the vertical in and XZ.
            %plots
            set(lineHandles{1}(1), 'YData', [sn(1) sn(1)]');
            set(lineHandles{3}(2), 'XData', [sn(1) sn(1)]');
        else
            %XZ view (sagittal)
            set(imHandles{3}, 'CData', squeeze(xzview(:,sn(2),:)));
            %imagesc(squeeze(seg(:,sn(2),:)), segRange); axis image

            %Moving in y, so change the vertical line in plot 1 and
            %horizontal line in plot 2.
            set(lineHandles{1}(2), 'XData', [sn(2) sn(2)]');
            set(lineHandles{2}(2), 'XData', [sn(2) sn(2)]');
        end
        %title(handles{whichView}, strcat(titleLines{whichView}, num2str(sn(offsliceCoord(whichView)))),...
        %    'Color', titleColors{whichView});
        %Why recreate the title if all we really want to do is change the
        %string.
        set(titleHandles{whichView}, 'String', ...
            sprintf('%s%03d', titleLines{whichView}, sn(offsliceCoord(whichView))));
    end

    function IntWinBound(varargin)
        %Here we perform the linear intensity windowing on the images.
        %First, determine the window selected by the user and update the
        %text appropriately.

        for whichSlider = 1:2
            slider_value = get(IntWinSliderHandle{whichSlider}, 'Value');
            segRange(whichSlider) = slider_value;
            set(IntWinTextBound{whichSlider}, 'String', segRange(whichSlider));
        end

        %Now, change the CLim of the four axes.  I couldn't find how to get
        %at just the axes (which have CLim defined), because subplot
        %apparently never returned a handle for my handles cell array.
        %Anyway, looping through the main figure's children did the trick.
        %(For now anyway, until I start using uipanels and such).
        for c = get(mainHandle, 'Children')'
            try
                set(c, 'CLim', segRange);
            catch exception
                continue
            end
        end

    end

    function defineBound(varargin)
        for whichBound = 1:2
            bound_str = get(IntWinTextBound{whichBound},'String');
            if strcmp(bound_str, 'max')
                bound_value = max(origSegRange);
            elseif strcmp(bound_str, 'min')
                bound_value = min(origSegRange);
            else
                bound_value = str2double(bound_str);
            end
            set(IntWinTextBound{whichBound},'String', num2str(bound_value));
            segRange(whichBound) = bound_value;
        end

        for c = get(mainHandle, 'Children')'
            try
                set(c, 'CLim', segRange);
            catch exception
                continue
            end
        end
    end

    function specifyColormap(hObject, ~)
        choice = get(hObject, 'Value');

        colormap(the_cmaps{choice});

    end

end


text:
function varargout = solver_sBPDN_WW( A, alpha, W1, beta, W2, b, epsilon, mu, x0, z0, opts, varargin )
% SOLVER_SBPDN_WW BPDN with two separate (weighted) l1-norm terms. Uses smoothing.
% [ x, out, opts ] = solver_sBPDN_WW( A, alpha, W_1, beta, W_2, b, epsilon, mu, x0, z0, opts )
%    Solves the smoothed basis pursuit denoising problem
%        minimize alpha*norm(W_1 x,1) + beta*norm(W_2 x, 1) + 0.5*mu*(x-x0).^2
%        s.t.     norm(A*x-b,2) <= epsilon
%    by constructing and solving the composite dual.
%    A, W_1 and W_2 must be a linear operator or matrix, and b must be a vector. The
%    initial points x0, z0 and the options structure opts are optional.
%    See also solver_sBPDN and solver_sBPDN_W

% Supply default values
error(nargchk(8,12,nargin));
if nargin < 9, x0 = []; end
if nargin < 10, z0 = []; end
if nargin < 11, opts = []; end
if ~isfield( opts, 'restart' ), opts.restart = 5000; end

if epsilon < 0
    error('TFOCS error: epsilon is negative');
end
if ~epsilon
    error('TFOCS error: cannot handle epsilon = 0.  Please call solver_sBP instead');
elseif epsilon < 100*builtin('eps')
    warning('TFOCS:badConstraint',...
        'TFOCS warning: epsilon is near zero; consider calling solver_sBP instead');
end

% Need to estimate the norms of A*A' and W*W' in order to be most efficient
if isfield( opts, 'noscale' ) && opts.noscale,
    normA2 = 1; normW12 = 1; normW22 = 1;
else
    normA2 = []; normW12 = []; normW22 = [];
    if isfield( opts, 'normA2'  )
        normA2 = opts.normA2;
        opts = rmfield( opts, 'normA2' );
    end
    if isfield( opts, 'normW12' )
        normW12 = opts.normW12;
        opts = rmfield( opts, 'normW12' );
    end
    if isfield( opts, 'normW22' )
        normW22 = opts.normW22;
        opts = rmfield( opts, 'normW22' );
    end
end
if isempty( normA2 ),
    normA2 = linop_normest( A ).^2;
end
if isempty( normW12 ),
    normW12 = linop_normest( W1 ).^2;
end
if isempty( normW22 ),
    normW22 = linop_normest( W2 ).^2;
end
if isempty(alpha),
    alpha = 1;
end
if isempty(beta),
    beta = 1;
end

proxScale1 = sqrt( normW12 / normA2 );
proxScale2 = sqrt( normW22 / normA2 );
prox       = { prox_l2( epsilon ), ...
               proj_linf( proxScale1 * alpha ),...
               proj_linf( proxScale2 * beta ) };
W1         = linop_compose( W1, 1 / proxScale1 );
W2         = linop_compose( W2, 1 / proxScale2 );
[varargout{1:max(nargout,1)}] = ...
    tfocs_SCD( [], { A, -b; W1, 0; W2, 0 }, prox, mu, x0, z0, opts, varargin{:} );

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function [kHz] = THz2kHz(THz)
% Convert frequency from terahertz to kilohertz.
% Chad A. Greene 2012
kHz = THz*1e+9;


text:
function [c] = cspace3(obs, bot, theta_steps, post_rotation)
if nargin < 4
  post_rotation = eye(2);
end
import iris.cspace.minkowski_sum;
if iscell(obs)
  padded = iris.pad_obstacle_points(obs);
  obs = cell2mat(reshape(padded, size(padded, 1), [], length(obs)));
end
bot = -bot;

if isscalar(theta_steps)
  th = linspace(-pi, pi, theta_steps);
else
  th = theta_steps;
end

c = zeros(3, 2*size(bot, 2), size(obs,3) * (length(th)-1));

if size(obs, 2) == 1
  idx = 0;
  for j = 1:(length(th)-1)
    rbot0 = post_rotation * iris.util.rotmat(th(j)) * bot;
    c_obs0 = bsxfun(@plus, rbot0, obs);

    rbot1 = post_rotation * iris.util.rotmat(th(j+1)) * bot;
    c_obs1 = bsxfun(@plus, rbot1, obs);

    c(:,:,idx+(1:size(obs,3))) = [c_obs0, c_obs1;
                                  th(j)-1e-3+zeros([1,size(c_obs0,2),size(obs,3)]), th(j+1)+1e-3+zeros([1,size(c_obs1,2),size(obs,3)])];
    idx = idx + size(obs,3);
  end
else
  idx = 1;
  for j = 1:(length(th)-1)
    for k = 1:size(obs,3)
      rbot0 = post_rotation * iris.util.rotmat(th(j)) * bot;
      c_obs0 = minkowski_sum(rbot0, obs(:,:,k));

      rbot1 = post_rotation * iris.util.rotmat(th(j+1)) * bot;
      c_obs1 = minkowski_sum(rbot1, obs(:,:,k));
      c_pts = [c_obs0, c_obs1;
              th(j)*ones(1,size(c_obs0,2)) - 1e-3, th(j+1)*ones(1,size(c_obs1,2)) + 1e-3];
      if size(c_pts, 2) > size(c, 2)
        c = [c, zeros([size(c, 1), size(c_pts,2)-size(c,2), size(c,3)])];
      end
      c(:,:,idx) = c_pts;
      idx = idx + 1;
    end
  end
end

author: rdeits

repo: iris-distro

save_path: github-repos/MATLAB/rdeits-iris-distro


text:
function [proj,geo,angles]=NikonDataLoader(filepath,varargin)
%NikonDataLoader(filepath) Loads Nikon uCT datasets into TIGRE standard
%
%   NikonDataLoader(filepath OPT,VAL,...) uses options and values.
%      These are options in case you don't want to load the entire
%      dataset, but only particular sets of projections.
%      The possible options in OPT are:
%           'sampling': type of sampling. default 'equidistant' Can be:
%                  'equidistant': equidistantly sample the entire set
%                       of angles. 'num_angles' should be set for partial
%                       loading of the data.
%                  'step': sample the entire set of projections every
%                         'sampling_step' angles.
%                  'continous': Load the first 'num_angles' amount of
%                             angles only.
%
%           'num_angles': Number of total angles to load. Default all of
%                    them. Useful for 'equidistant' and 'continous' loading
%
%           'sampling_step': step to load when loading projections.
%                   Default=1. Useful for 'step' loading.
%%

[geo,angles]=readXtekctGeometry(filepath);
[proj,geo,angles]=loadNikonProjections(filepath,geo,angles,varargin{:});

end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function X = model_2_data(N, K, S)
    % Data as per model 2 in eldar2010average paper
    gen  = spx.data.synthetic.SparseSignalGenerator(N, K, S);
    X = gen.complex_gaussian;
end

author: indigits

repo: sparse-plex


text:
function best_acc = LCKSVD_top(dataset, N_train, k, ...
        sparsitythres, valpha, vbeta)
    % Syntax: LCKSVD_top(dataset, N_train, k, sparsitythres, valpha, vbeta)
    addpath(genpath('utils'));  % add K-SVD box
    addpath(genpath('LCKSVD'));  % add K-SVD box
    addpath('build_spams');
    addpath('utils');
    best_acc = zeros(1, 2);
    if nargin == 0
        dataset = 'myYaleB';
        N_train = 10;
        k = 8;
        sparsitythres = 10;
        valpha = 0.01;
        vbeta = 0.01;
    end
    %% get data
    [dataset, Y_train, Y_test, label_train, label_test] = train_test_split(...
        dataset, N_train);
    %% main
    [acc, rt] = LCKSVD_wrapper(Y_train, label_train, Y_test, label_test,...
                    k, sparsitythres, valpha, vbeta);
    %% save data
    if ~exist('results', 'dir')
        mkdir('results');
    end
    if ~exist(fullfile('results', 'LCKSVD'), 'dir')
            mkdir('results', 'LCKSVD');
        end
    fn1 = fullfile('results', 'LCKSVD', strcat(dataset, '_N_', ...
        num2str(N_train), '_k_', num2str(k), '_a_', num2str(valpha), '_b_', ...
        num2str(vbeta), '_', getTimeStr(), '_1.mat'));
    disp(fn1);
    fn2 = fullfile('results', 'LCKSVD', strcat(dataset, '_N_', ...
        num2str(N_train), '_k_', num2str(k), '_a_', num2str(valpha), ...
        '_b_', num2str(vbeta), '_', getTimeStr(), '_2.mat'));
    best_acc = acc;
    rrt = rt;
    acc = best_acc(1);
    rt = rrt(1);
    save(fn1, 'acc', 'rt');
    acc = best_acc(2);
    rt = rrt(2);
    save(fn2, 'acc', 'rt');
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function test_suite = test_a_bit
initTestSuite

function test_now

function test_later

author: ricedsp


text:
function scanFileNameC = scan2imageOut(planC,scanNumV,tmpDirPath,reorientFlag,extn,dataType)

if ~exist('tmpDirPath','var') || isempty(tmpDirPath) || ~exist(tmpDirPath,'dir')
    tmpDirPath = fullfile(getCERRPath, 'ImageRegistration', 'tmpFiles');
end

if ~exist('reorientFlag','var') || isempty(reorientFlag)
    reorientFlag = 1;
end

if ~exist('extn','var') || isempty(extn)
    extn = 'nii';
end

if ~exist('dataType','var')
    dataType = [];
end

for i = 1:numel(scanNumV)
    scanNum = scanNumV(i);
    [scanUniqName, ~] = genScanUniqName(planC,scanNum);
    [affineMat,scan3M_RAS,voxel_size] = getPlanCAffineMat(planC, scanNum, reorientFlag);
    [~,orientationStr,~] = returnViewerAxisLabels(planC,scanNum);
    if ~isempty(dataType)
        eval(['scan3M_RAS = ' dataType '(scan3M_RAS);']);
    end
    qOffset = affineMat(1:3,end)';
    if strcmpi(extn,'nii')
        scanFileName = fullfile(tmpDirPath, ['scan_' num2str(scanNum) '_' scanUniqName '.nii']);
        niiC{i} = vol2nii(scan3M_RAS,affineMat,qOffset,voxel_size,[],scanFileName);
    elseif strcmpi(extn,'nrrd')
        scanFileName = fullfile(tmpDirPath, ['scan_' num2str(scanNum) '_' scanUniqName '.nrrd']);
        vol2nrrd(scan3M_RAS,affineMat,qOffset,voxel_size,orientationStr,scanFileName);
    end
    scanFileNameC{i} = scanFileName;
end


text:
%DISP Display dataset
%
%	DISP(A)
%
% Display dataset A and its content

author: marianux

repo: ecg-kit

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


text:
%------------------------------- getting ready ----------------------------
clc
clear all
close all
%--------------------------------------------------------------------------

global A E
tol  = 0.001;
max  = 20;
init = 0;
xn   = init*ones(2,78);          % estimation for initial trejctory for states
pn   = init*ones(2,78);          % estimation for initial trejctory for co-states
check = [0.01:0.01:0.78];

for loop = 1:max
    xh1 = [0 0 1 0];             % initial conditions for homogeneous solution1
    xh2 = [0 0 0 1];             % initial conditions for homogeneous solution2
    xp0 = [0.05 0 0 0];          % initial conditions for particular  solution0
    a = 0;
    b = 0.01;
    for I=1:78
        t       = [xn(:,I)' pn(:,I)'];
        pointer = 100*loop+I
        clc
        A  = calculate_matrix(t);
        x1 = [-2*(t(1)+0.25)+(t(2)+0.5)*exp(25*t(1)/(t(1)+2))-(t(1)+0.25)^2*t(3)/0.2    0.5-t(2)-(t(2)+0.5)*exp(25*t(1)/(t(1)+2))];
        x2 = [-2*t(1)+2*t(3)-t(3)*(t(2)+0.5)*(50/(t(1)+2)^2)*exp(25*t(1)/(t(1)+2))+t(3)^2*(t(1)+0.25)/0.2+t(4)*(t(2)+0.5)*(50/(t(1)+2)^2)*exp(25*t(1)/(t(1)+2)) -2*t(2)-t(3)*exp(25*t(1)/(t(1)+2))+t(4)*(1+exp(25*t(1)/(t(1)+2)))];
        e1 = [x1 x2]';
        e2 = A*t';
        E  = e1 - e2;

        clear t1 c
        x0  = [xp0 xh1 xh2]';
        [t1,c]  = ode45('stirred_model', [a b], x0);

        L       = length(t1);
        X1p(I)  = c(L,1);
        X2p(I)  = c(L,2);
        P1p(I)  = c(L,3);
        P2p(I)  = c(L,4);
        xp0     = [X1p(I) X2p(I) P1p(I) P2p(I)];

        X11h(I) = c(L,5);
        X12h(I) = c(L,6);
        P11h(I) = c(L,7);
        P12h(I) = c(L,8);
        xh1     = [X11h(I) X12h(I) P11h(I) P12h(I)];

        X21h(I) = c(L,9 );
        X22h(I) = c(L,10);
        P21h(I) = c(L,11);
        P22h(I) = c(L,12);
        xh2     = [X21h(I) X22h(I) P21h(I) P22h(I)];

        % simulation time for next interval of integrtion
        a = b;
        b = b + 0.01;

    end

    % calculting constant to satisfy boudary condition
    B1  = [0 0]'-[P1p(78) P2p(78)]';
    ph1 = [P11h(78) P12h(78)]';
    ph2 = [P21h(78) P22h(78)]';
    A1  = [ph1 ph2];
    const = A1\B1;

    % calculation of states and costates
    state   = [X1p;X2p] + [const(1)*X11h;const(1)*X12h]+[const(2)*X21h;const(2)*X22h];
    costate = [P1p;P2p] + [const(1)*P11h;const(1)*P12h]+[const(2)*P21h;const(2)*P22h];
    test = terminate(state,xn,costate,pn);
    fprintf('Terminate Check = %f\n',test);
    const
    pause
    clc
    if test < tol
        break
    end

    xn = state;
    pn = costate;

end

plot(check,state,'-');
title('Optimal Trejctory');
grid on
%--------------------------------------------------------------------------
for j = 1:78
   u(j)    = 0.5*costate(1,j)*(state(1,j)+0.25);
   intg(j) = state(1,j)^2+state(2,j)^2+0.1*u(j)^2;
end
figure
plot(check,u);
title('U^*');
grid on
J = integrate(check,intg);
fprintf('Optimal Cost = %f\n',J);


text:
function sum_wf(handles,einzelmat)
    imagefile = get(handles.im_file_edit, 'String');
    outpath = get(handles.out_dir_edit, 'String');
    if outpath(length(outpath))~='\'
        outpath = [outpath,'\'];
    end

    %Obtain info on image type, etc
    file_info = get_file_info2(imagefile, handles);

    load(einzelmat,'image_sum','x_size','y_size','x_offset','y_offset');

    exf       = str2double(get(handles.exf_edit,'String')); %expansion factor
    ishift     = str2double(get(handles.ishift,'String'));
    jshift     = str2double(get(handles.jshift,'String'));
    xw_render = str2double(get(handles.xw_render_edit,'String'));
    yw_render = str2double(get(handles.yw_render_edit,'String'));

    if xw_render ~= 0
        xw=xw_render;
    else
        xw=x_size;
    end
    if yw_render ~= 0
        yw=yw_render;
    else
        yw=y_size;
    end

    x_off=x_offset+ishift;
    y_off=y_offset+jshift;

    wf_sum=zeros(yw,xw);
    wf_sum=image_sum(y_off:y_off+yw-1,x_off:x_off+xw-1);
    wf_sum=wf_sum/max(max(wf_sum));
    wf_sum=wf_sum/max(max(wf_sum));

    imex=zeros(yw*exf,xw*exf);
    for j=1:yw
        for i=1:xw
            x0=(i-1)*exf+1;
            x1=x0+exf-1;
            y0=(j-1)*exf+1;
            y1=y0+exf-1;
            imex(y0:y1,x0:x1)=wf_sum(j,i);
        end
    end

    figure
    imagesc(imex)
    title(['Expanded Widefield Sum: Frames: ',num2str(file_info.start),'-',num2str(file_info.stop)])
    colormap gray
    axis image

    if get(handles.write_image_checkbox,'Value') %write to file if box is checked
        if (file_info.part_name(length(file_info.part_name)) == '_') %Avoid saving with repeat seperators
            wf_sum_file=[outpath,'\',file_info.part_name,num2str(file_info.start),'-',num2str(file_info.stop),'exp',num2str(exf),'_wf.tif'];
        else
            wf_sum_file=[outpath,'\',file_info.part_name,num2str(file_info.start),'-',num2str(file_info.stop),'_exp',num2str(exf),'_wf.tif'];
        end
        imwrite(uint8(imex*255),wf_sum_file,'Compression','none');
    end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function c=ref_dsti_1(f)
%REF_DSTI_1  Reference Discrete Sine Transform type I
%   Usage:  c=ref_dsti_1(f);
%
%

L=size(f,1);
W=size(f,2);

if L==1
  c=f;
  return;
end;

R=1/sqrt(2)*[zeros(1,L,assert_classname(f));...
	     eye(L);
	     zeros(1,L,assert_classname(f));...
	     -flipud(cast(eye(L),assert_classname(f)))];

c=i*R'*dft(R*f);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function g = modelLatentGradients(model)

% MODELLATENTGRADIENTS Gradients of the latent variables for dynamics models in the GPLVM.

% FGPLVM

fhandle = str2func([model.type 'LatentGradients']);
g = fhandle(model);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function differ_test01 ( )

%*****************************************************************************80
%
%% DIFFER_TEST01 tests DIFFER_MATRIX.
%
%  Discussion:
%
%    DIFFER_MATRIX computes a modified Vandermonde matrix A1.
%
%    The solution of a system A1 * X1 = B is related to the solution
%    of the system A2 * X2 = B, where A2 is the standard Vandermonde
%    matrix, simply by X2(I) = X1(I) * A(I,1).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2013
%
%  Author:
%
%    John Burkardt
%
  n = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'DIFFER_TEST01\n' );
  fprintf ( 1, '  Demonstrate that the DIFFER matrix is "really"\n' );
  fprintf ( 1, '  a Vandermonde matrix.\n' );

  stencil(1:n,1) = [ 2.5; 3.3; -1.3; 0.5 ];
  x1(1:n,1) = [ 1.0; 2.0; 3.0; 4.0 ];
  a = differ_matrix ( n, stencil );
  r8mat_print ( n, n, a, '  Stencil matrix:' );
  b(1:n,1) = a(1:n,1:n) * x1(1:n,1);
%
%  Set up and solve the DIFFER system.
%
  x1 = a \ b;

  r8vec_print ( n, x1, '  Solution of DIFFER system:' );
%
%  R8VM_SL solves the related Vandermonde system.
%  A simple transformation gives us the solution to the DIFFER system.
%
  job = 0;
  [ x2, info ] = r8vm_sl ( n, stencil, b, job );

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST01 - Warning!\n' );
    fprintf ( 1, '  VANDERMONDE system is singular.\n' );
    error ( 'TEST01 - Vandermonde system is singular.' );
    return
  end

  r8vec_print ( n, x2, '  Solution of VANDERMONDE system:' );

  x2(1:n,1) = x2(1:n,1) ./ stencil(1:n,1);
  r8vec_print ( n, x2, '  Transformed solution of VANDERMONDE system:' );

  return
end


text:
function [vdeg] = trideg2(pp,tt)
%TRIDEG2 calc. topological degree for vertices in a 2-simpl-
%ex triangulation.
%   [VDEG] = TRIDEG2(VERT,TRIA) returns the no. of triangles
%   incident to each vertex. VDEG is a V-by-1 array of vert-
%   ex degrees, VERT is a V-by-D array of XY coordinates,
%   and TRIA is a T-by-3 array of vertex indexing, where
%   each row defines a triangle, such that
%   VERT(TRIA(II,1),:), VERT(TRIA(II,2),:) and
%   VERT(TRIA(II,3),:) are the coordinates of the II-TH tri-
%   angle.
%
%   See also TRISCR2, TRIVOL2, TRIANG2, TRIBAL2

%   Darren Engwirda : 2017 --
%   Email           : de2363@columbia.edu
%   Last updated    : 10/07/2018

%---------------------------------------------- basic checks
    if (~isnumeric(pp) || ~isnumeric(tt) )
        error('trideg2:incorrectInputClass' , ...
            'Incorrect input class.') ;
    end

%---------------------------------------------- basic checks
    if (ndims(pp) ~= +2 || ndims(tt) ~= +2 )
        error('trideg2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end
    if (size(pp,2) < +2 || size(tt,2) < +3 )
        error('trideg2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end

    nvrt = size(pp,1) ;
    ntri = size(tt,1) ;

%---------------------------------------------- basic checks
    if (min(min(tt(:,1:3))) < +1 || ...
            max(max(tt(:,1:3))) > nvrt )
        error('trideg2:invalidInputs', ...
            'Invalid TRIA input array.') ;
    end

%------------------------------------- compute vertex degree
    vdeg = sum(sparse( ...
        tt(:,1:3),repmat( ...
            (1:ntri)',1,3),+1,nvrt,ntri),2) ;

end

author: CHLNDDEV


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Jellyfish Example Courtesy of Alexander P. Hoover, PhD
%
% Converted from IBAMR: 1/16/2018 by NAB.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Make_Jelly_Geometry()

close all;
clear all;
L = 8;                              % height of computational domain (m) for keeping desired resolution
Lh = 10;                            % actual height of computational domain (m) (MATCHES INPUT2D)
Lw = 3;                             % width of computational domain (m) (MATCHES INPUT2D)
N = 768;                            % number of Cartesian grid meshwidths at the finest level of the AMR grid
dx = L/N;                           % Cartesian mesh width (m)
ds = dx/2;

a=.5;                               % bell radius (semi-minor axis, horizontal axis, note width=2a)
b=.75;                              % bell semi-major axis
d=-0.25;
factor_a=.8;

F=1e5; %5e0

theta=zeros(1000,1);
theta_lim=asin(d/b);
theta_test=pi/2;

x_points=zeros(1000,1);
z_points=zeros(1000,1);
id_points=zeros(1000,1);
offset = 0;

kappa_spring = 1e7; %1e5               % spring constant (Newton)
kappa_beam = 2.5e5; %1e5    %5e3              % beam stiffness constant (Newton m^2)
%kappa_beam_flexible = kappa_beam/5;   % beam stiffness constant (Newton m^2)
kappa_target = kappa_spring;           % target point penalty spring constant (Newton)

c=0;
while(theta_test<(pi-theta_lim))
    c=c+1;
    theta(c)=theta_test;

    x_points(c)=a*cos(theta(c));
    z_points(c)=b*sin(theta(c));
    id_points(c)=c-1;

    theta_test=ds/((a*sin(theta(c)))^(2)+(b*cos(theta(c)))^(2))^(.5)+theta(c);

end

c_stiff=c;

npts=2*c-1;
npts_wing=floor(npts/2);
npts_musc=floor(npts_wing/4);

for j=(c+1):(npts)
    x_points(j)=-1*x_points(j-c+1);
    z_points(j)=z_points(j-c+1);
    id_points(j)=j-1;
end

mesh_name = 'jelly';
xShift = 1.5;
yShift = 2;

x_points=x_points(1:npts)+xShift;
z_points=z_points(1:npts)+yShift;
it_points=id_points(1:npts);

plot(x_points(:),z_points(:),'*'); hold on;
axis([0 8 0 8])

% Lag Pts to Mess up Flow At Edge
xBlock = ds:4*ds:Lw-ds;
yBlock = (Lh-5*ds)*ones(1,length(xBlock))+ds;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Print .vertex information
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

vertex_fid = fopen([mesh_name num2str(N) '.vertex'], 'w');

    fprintf(vertex_fid, '%d\n', npts + npts_musc*2 + length(xBlock));
    lag_ct = 0;

    %
    % bell
    %
    for j=1:npts
        fprintf(vertex_fid, '%1.16e %1.16e\n', x_points(j), z_points(j));
        lag_ct = lag_ct + 1;
    end

    %
    % muscles
    %
    for s = 1:npts_musc
        fprintf(vertex_fid, '%1.16e %1.16e\n', x_points(npts_wing+1-npts_musc+s), z_points(npts_wing+1-npts_musc+s));
        plot(x_points(npts_wing+1-npts_musc+s),z_points(npts_wing+1-npts_musc+s),'r*'); hold on;
        lag_ct = lag_ct + 1;
    end
    for s = 1:npts_musc
        fprintf(vertex_fid, '%1.16e %1.16e\n', x_points(npts-npts_musc+s), z_points(npts-npts_musc+s));
        plot(x_points(npts-npts_musc+s),z_points(npts-npts_musc+s),'r*'); hold on;
        lag_ct = lag_ct + 1;
    end

    for ii=1:length(xBlock)
        fprintf(vertex_fid, '%1.16e %1.16e\n', xBlock(ii), yBlock(ii));
    end

fclose(vertex_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Print .spring information
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

spring_fid = fopen([mesh_name num2str(N) '.spring'], 'w');

    npts_spring_type1=npts-1;

    fprintf(spring_fid, '%d\n', npts-1 + npts_musc);

    fprintf('\nNumber of springs before muscles: %d \n\n',npts-1)

    factor = 1;%ds^2/ds;

    %
    % bell
    %
    for s = 1:c-1
        resting=sqrt((x_points(s)-x_points(s+1))^(2)+(z_points(s)-z_points(s+1))^(2));
        fprintf(spring_fid, '%d %d %1.16e %1.16e %d\n', id_points(s)+1, id_points(s+1)+1, kappa_spring*ds/(ds^2)*factor, resting, 1);
    end
    for s = c+1:npts-1
        resting=sqrt((x_points(s)-x_points(s+1))^(2)+(z_points(s)-z_points(s+1))^(2));
        fprintf(spring_fid, '%d %d %1.16e %1.16e %d\n', id_points(s)+1, id_points(s+1)+1, kappa_spring*ds/(ds^2)*factor, resting, 1);
    end
    resting=sqrt((x_points(1)-x_points(c+1))^(2)+(z_points(1)-z_points(c+1))^(2));
    fprintf(spring_fid, '%d %d %1.16e %1.16e %d\n', id_points(1)+1, id_points(c+1)+1, kappa_spring*ds/(ds^2)*factor, resting, 1);

    %
    % muscles
    %
    for s = 1:npts_musc
        fprintf(spring_fid, '%d %d %1.16e %1.16e %d\n',npts+s-1+1, npts+s+npts_musc-1+1, F, 0, 1);
    end

    fclose(spring_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Print .nonInv_beam information
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

beam_fid = fopen([mesh_name num2str(N) '.nonInv_beam'], 'w');

    fprintf(beam_fid, '%d\n', npts-2);

    factor=1;% = (ds^4)/ds;

    for s = 2:c-1
        C1 = x_points(s-1)+x_points(s+1)-2*x_points(s);
        C2 = z_points(s-1)+z_points(s+1)-2*z_points(s);
        fprintf(beam_fid, '%d %d %d %1.16e %1.16e %1.16e\n', id_points(s-1)+1, id_points(s)+1, id_points(s+1)+1, kappa_beam*ds/(ds^4)*factor, C1, C2);
    end
    for s = c+2:npts-1
        C1 = x_points(s-1)+x_points(s+1)-2*x_points(s);
        C2 = z_points(s-1)+z_points(s+1)-2*z_points(s);
        fprintf(beam_fid, '%d %d %d %1.16e %1.16e %1.16e\n', id_points(s-1)+1, id_points(s)+1, id_points(s+1)+1, kappa_beam*ds/(ds^4)*factor, C1, C2);
    end

    C1 = x_points(c+2)+x_points(1)-2*x_points(c+1);
    C2 = z_points(c+2)+z_points(1)-2*z_points(c+1);
    fprintf(beam_fid, '%d %d %d %1.16e %1.16e %1.16e\n', id_points(c+2)+1, id_points(c+1)+1, id_points(1)+1, kappa_beam*ds/(ds^4)*factor, C1, C2);

    C1 = x_points(c+1)+x_points(2)-2*x_points(1);
    C2 = z_points(c+1)+z_points(2)-2*z_points(1);
    fprintf(beam_fid, '%d %d %d %1.16e %1.16e %1.16e\n', id_points(c+1)+1, id_points(1)+1, id_points(2)+1, kappa_beam*ds/(ds^4)*factor, C1, C2);

    fclose(beam_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% PRINT TARGET POINTS!!!
%
% print target points (flow blocker along edge)
k_Target = 2.5e6;
nBefore = lag_ct; % Counts pts in jellyfish for bookkeeping for .target file
struct_name = ['jelly' num2str(N)];
print_Lagrangian_Target_Pts(xBlock,k_Target,struct_name,nBefore)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints TARGET points to a file called 'struct_name'.target
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Target_Pts(xLag,k_Target,struct_name,nBefore)

    N = length(xLag);
    Nstart = nBefore+1;
    Nend = nBefore+N;

    target_fid = fopen([struct_name '.target'], 'w');

    fprintf(target_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.
    for s = Nstart:Nend
        fprintf(target_fid, '%d %1.16e\n', s, k_Target);
    end

    fclose(target_fid);

author: nickabattista


text:
function varargout = twister(varargin)
% VL_TWISTER Random number generator
%    VL_TWISTER() is essentially equivalent to MATLAB native RAND()
%    when using the Twister random number generator. VL_TWISTER(),
%    VL_TWISTER(M,N,P,...)  and VL_TWISTER([M N P ...]) are equivalent
%    to RAND(), RAND(M,N,P,...)  and RAND([M N P ...]) respectively.
%
%    The state of the random generator can be seeded by
%    VL_TWISTER('STATE', X), where X is a DOUBLE scalar (this is
%    equivalent to RAND('TWISTER', X)). The state can be read by
%    VL_TWISTER('STATE') (equivalent to RAND('TWISTER')) and set by
%    VL_TWISTER('STATE', STATE) (equivalent to RAND('TWISTER',
%    STATE)). Here STATE is a vector of 625 elements of class
%    UINT32. Finally VL_TWISTER('STATE',KEY) seeds the generator by a
%    vector of DOUBLE of length not greater than 624.
%
%    VL_TWISTER() is slightly faster than RAND(). Moreover it can be
%    used to control the state of the random number generator used by
%    all VLFEAT functions.
%
%    See also: VL_HELP().
[varargout{1:nargout}] = vl_twister(varargin{:});

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function adj_num = ns_adj_count ( node_num, triangle_num, ...
  variable_num, triangle_node, triangle_neighbor, node_u_variable, ...
  node_v_variable, node_p_variable )

%*****************************************************************************80
%
%% NS_ADJ_COUNT counts adjacencies in a Navier Stokes triangulation.
%
%  Discussion:
%
%    This routine is called to count the adjacencies, so that the
%    appropriate amount of memory can be set aside for storage when
%    the adjacency structure is created.
%
%    The value of ADJ_NUM computed and returned by this routine should
%    be identical to the value computed by NS_ADJ_COL_SET.
%
%    The triangulation is assumed to involve 6-node triangles.
%
%    Variables for the horizontal and vertical velocities are associated
%    with every node.  Variables for the pressure are associated only with
%    the vertex nodes.
%
%    We are interested in determining the number of nonzero entries in the
%    stiffness matrix of the Stokes equations, or the jacobian matrix of
%    the Navier Stokes equations.  To this end, we will say, somewhat
%    too broadly, that two variables are "adjacent" if their associated
%    nodes both occur in some common element.  This adjacency of variables
%    I and J is taken to be equivalent to the possible nonzeroness of
%    matrix entries A(I,J) and A(J,I).
%
%    A sparse compressed column format is used to store the counts for
%    the nonzeroes.  In other words, while the value ADJ_NUM reports the
%    number of adjacencies, the vector ADJ_COL is sufficient to allow us
%    to properly set up a sparse compressed matrix for the actual storage
%    of the sparse matrix, if we desire to proceed.
%
%  Local Node Numbering:
%
%       3
%    s  |\
%    i  | \
%    d  |  \
%    e  6   5  side 2
%       |    \
%    3  |     \
%       |      \
%       1---4---2
%
%         side 1
%
%  Variable Diagram:
%
%      UVP
%       |\
%       | \
%       |  \
%      UV   UV
%       |    \
%       |     \
%       |      \
%      UVP--UV--UVP
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer TRIANGLE_NUM, the number of triangles.
%
%    Input, integer VARIABLE_NUM, the number of variables.
%
%    Input, integer TRIANGLE_NODE(6,TRIANGLE_NUM), lists the nodes that
%    make up each triangle.  The first three nodes are the vertices,
%    in counterclockwise order.  The fourth value is the midside
%    node between nodes 1 and 2; the fifth and sixth values are
%    the other midside nodes in the logical order.
%
%    Input, integer TRIANGLE_NEIGHBOR(3,TRIANGLE_NUM), for each side of
%    a triangle, lists the neighboring triangle, or -1 if there is
%    no neighbor.
%
%    Input, integer NODE_U_VARIABLE(NODE_NUM), NODE_V_VARIABLE(NODE_NUM),
%    NODE_P_VARIABLE(NODE_NUM), the index of the horizontal velocity,
%    vertical velocity and pressure variables associated with a node,
%    or -1 if no such variable is associated with the node.
%
%    Output, integer ADJ_NUM, the number of Navier Stokes variable adjacencies.
%
  triangle_order = 6;

  adj_num = 0;
%
%  Set every variable to be adjacent to itself.
%
  adj_num = variable_num;
%
%  Set every variable to be adjacent to the other variables associated with
%  that node.
%
%  U <=> V
%  U <=> P (if there is a P variable)
%  V <=> P (if there is a P variable)
%
  for node = 1 : node_num

    adj_num = adj_num + 2;

    p1 = node_p_variable(node);

    if ( 0 < p1 )
      adj_num = adj_num + 4;
    end

  end
%
%  Examine each triangle.
%
  for triangle = 1 : triangle_num
%
%  For sure, we add the new adjacencies:
%
%    U5 V5 <=> U1 V1 P1
%    U6 V6 <=> U2 V2 P2
%    U4 V4 <=> U3 V3 P3
%    U5 V5 <=> U4 V4
%    U6 V6 <=> U4 V4
%    U6 V6 <=> U5 V5
%
    adj_num = adj_num + 60;
%
%  Add edges (1,2), (1,4), (2,4) if this is the first occurrence,
%  that is, if the edge (1,4,2) is on a boundary (TRIANGLE2 <= 0)
%  or if this triangle is the first of the pair in which the edge
%  occurs (TRIANGLE < TRIANGLE2).
%
%  Maybe add
%
%    U1 V1 P1 <=> U2 V2 P2
%    U1 V1 P1 <=> U4 V4
%    U2 V2 P2 <=> U4 V4
%
    triangle2 = triangle_neighbor(1,triangle);

    if ( triangle2 < 0 | triangle < triangle2 )
      adj_num = adj_num + 42;
    end
%
%  Maybe add
%
%    U2 V2 P2 <=> U3 V3 P3
%    U2 V2 P2 <=> U5 V5
%    U3 V3 P3 <=> U5 V5
%
    triangle2 = triangle_neighbor(2,triangle);

    if ( triangle2 < 0 | triangle < triangle2 )
      adj_num = adj_num + 42;
    end
%
%  Maybe add
%
%    U1 V1 P1 <=> U3 V3 P3
%    U1 V1 P1 <=> U6 V6
%    U3 V3 P3 <=> U6 V6
%
    triangle2 = triangle_neighbor(3,triangle);

    if ( triangle2 < 0 | triangle < triangle2 )
      adj_num = adj_num + 42;
    end

  end

  return
end


text:
function [val grad]=test_func(x)

val=(x-2).^2;
grad=2*(x-2);


text:
function [y,dy] = softbndloss(x,slb,sub,TolCon)
%SOFTBNDLOSS Loss function for soft bounds for function minimization.

% Penalization relative scale
if nargin < 4 || isempty(TolCon); TolCon = 1e-3; end

compute_grad = nargout > 1;     % Compute gradient only if requested

ell = (sub - slb).*TolCon;

y = 0;
dy = zeros(size(x));

idx = x < slb;
if any(idx)
    y = y + 0.5*sum(((slb(idx) - x(idx))./ell(idx)).^2);
    if compute_grad
        dy(idx) = (x(idx) - slb(idx))./ell(idx).^2;
    end
end

idx = x > sub;
if any(idx)
    y = y + 0.5*sum(((x(idx) - sub(idx))./ell(idx)).^2);
    if compute_grad
        dy(idx) = (x(idx) - sub(idx))./ell(idx).^2;
    end
end

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
% LOG    Natural logarithm.
%    LOG(X) is the natural logarithm of the elements of X.
%    Complex results are produced if X is not positive.
%
%    See also LOG1P, LOG2, LOG10, EXP, LOGM, REALLOG.
%
%    Reference page in Doc Center
%       doc log
%
%    Other functions named log
%
%       codistributed/log    gpuArray/log    sym/log    ts/log
%       fints/log
%

author: jmaih