text:
function varargout = ft_version(command)

% FT_VERSION returns the version of FieldTrip and the path where it is installed
%
% FieldTrip is not released with version numbers as "2.0", "2.1", etc. Instead, we
% share our development version on http://github.com/fieldtrip/fieldtrip. You can use
% git to make a local clone of the development version. Furthermore, we make
% more-or-less daily releases of the code available on
% https://github.com/fieldtrip/fieldtrip/releases and as zip file on our FTP server.
%
% If you use git with the development version, the version is labeled with the hash
% of the latest commit like "128c693". You can access the specific version "XXXXXX"
% at https://github.com/fieldtrip/fieldtrip/commit/XXXXXX.
%
% If you download the daily released version from our FTP server, the version is part
% of the file name "fieldtrip-YYYYMMDD.zip", where YYY, MM and DD correspond to year,
% month and day.
%
% Use as
%   ft_version
% to display the latest revision number on screen, or
%   [ftver, ftpath] = ft_version
% to get the version and the installation root directory.
%
% When using git with the development version, you can also get additional information with
%   ft_version revision
%   ft_version branch
%   ft_version clean
%
% See also FT_PLATFORM_SUPPORTS, VERSION, VER, VERLESSTHAN

% Copyright (C) 2012-2019, Eelke Spaak, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/donders/fieldtrip
% 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$

persistent issvn
persistent isgit
persistent ftver
persistent ftpath
persistent is_pc

if isempty(is_pc)
  is_pc = ispc;
end

if isempty(ftpath)
  ftpath = fileparts(mfilename('fullpath'));
  ftpath = ftpath(1:end-10); % strip away '/utilities' where this function is located
end

% set the defaults
clean = 'no';
branch = 'unknown';

if is_pc
  % this requires a file extension
  ext = '.exe';
else
  ext = '';
end

if isempty(issvn)
  % are we dealing with an SVN working copy of FieldTrip?
  issvn = isfolder(fullfile(ftpath, '.svn'));
end

if isempty(isgit)
  % are we dealing with an GIT working copy of FieldTrip?
  isgit = exist(fullfile(ftpath, '.git'), 'file');
  % is the git command line executable available?
  if isgit
    [status, output] = system(sprintf('git%s --version', ext));
    if status>0
      if ~is_pc
        % the git command line executable will probably not be available on windows
        ft_warning('you seem to have an GIT development copy of FieldTrip, yet ''git'' does not work as expected');
      end
      isgit = false;
    end
  end
end

if nargin<1
  % this is only supported for git
  command = 'revision';
end

if ~isempty(ftver) && ~isempty(ftpath) && nargin<1
  % use the previously determined values

elseif issvn
  % use svn system call to determine latest revision
  [status, output] = system(sprintf('cd %s && svn%s info', ftpath, ext));
  if status > 0
    if ~is_pc
      % the command line tools will probably not be available on windows
      ft_warning('you seem to have an SVN development copy of FieldTrip, yet ''svn info'' does not work as expected');
    end
    ftver = 'unknown';
  else
    rev = regexp(output, 'Revision: (.*)', 'tokens', 'dotexceptnewline');
    rev = rev{1}{1};
    ftver = ['r' rev];
  end

elseif isgit
  % use git system call to determine latest revision
  switch command
    case 'revision'
      [status, output] = system(sprintf('cd %s && git%s rev-parse --short HEAD', ftpath, ext));
      ftver = strtrim(output); % remove trailing newline character
    case 'branch'
      [status, output] = system(sprintf('cd %s && git%s rev-parse --abbrev-ref HEAD', ftpath, ext));
      branch = strtrim(output); % remove trailing newline character
    case 'clean'
      [status, output] = system(sprintf('cd %s && git%s diff --quiet --exit-code', ftpath, ext));
      if status
        clean = 'no';
      else
        clean = 'yes';
      end
    otherwise
      ft_error('unsupported command "%s"', command);
  end % switch command

elseif isequal(regexp(ftpath, ['.*\' filesep '[fF]ieldtrip-fieldtrip-[[0-9][a-z]]{7}']), 1)
  % this corresponds with being downloaded from the Mathworks file exchange link to github
  % which results in a ftpath like /Users/robert/matlab/fieldtrip-fieldtrip-851478d
  ftver = ftpath(end-6:end);

elseif isequal(regexp(ftpath, ['.*\' filesep '[fF]ieldtrip-20[0-9]{6}']), 1)
  % this corresponds with the daily version from the ftp server
  % which results in a ftpath like /Users/robert/matlab/fieldtrip-20160317
  ftver = ftpath(end-7:end);

elseif isequal(regexp(ftpath, ['.*\' filesep '[fF]ieldtrip-lite20[0-9]{6}']), 1)
  % this corresponds with the daily version from the ftp server
  % which results in a ftpath like /Users/robert/matlab/fieldtrip-20160317
  ftver = ftpath(end-7:end);

else
  % get it from the Contents.m file in the FieldTrip directory
  if ~isdeployed
    tmp = ver(ftpath);
    ftver = tmp.Version;
  else
    ftver = 'deployed';
  end
end % if issvn, isgit or otherwise

if ~nargout
  switch command
    case 'revision'
      fprintf('\nThis is FieldTrip, revision %s.\n\n', ftver);
    case 'branch'
      fprintf('\nThis is FieldTrip, branch %s.\n\n', branch);
    case 'clean'
      fprintf('\nThis is FieldTrip, clean %s.\n\n', clean);
  end
else
  switch command
    case 'revision'
      varargout = {ftver, ftpath};
    case 'branch'
      varargout = {branch};
    case 'clean'
      varargout = {clean};
  end
end


text:
function [CV,CF] = spline_cage(P,C,varargin)
  % SPLINE_CAGE Build a triangle mesh "cage" surrounding an input spline curve
  % trying to keep control points inside the triangles.
  %
  % [CV,CF] = spline_cage(P,C)
  %
  % Inputs:
  %   P  #P by 2 list of control point locations
  %   C  #C by 4 list of cubic Bezier curves as indices into rows of P
  %   Optional:
  %     'GridSize' followed by number of cells along the x-axis in the grid used
  %     for unsigned distance calculation and bwmesh {1000}
  %     'Pad' followed by number of cells to pad grid by {4}
  %     'PolyTol' followed by tolerance parameter for `spline_to_poly` used for
  %     approximating spline distances {0.1}
  %     'Tol'  followed by tolerance parameter for `bwmesh` {1.0}
  %     'SmoothingIters'  followed by smoothing iters parameter for `bwmesh`
  %       {10}
  % Outputs:
  %   CV  #CV by 2 list of cage triangulation vertex positions
  %   CF  #CF by 3 list of cage triangulation triangle indices into rows of CV
  %

  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'GridSize','Pad','PolyTol','Tol','SmoothingIters'}, ...
    {'gs','pad','poly_tol','bwmesh_tol','bwmesh_smoothing_iters'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  gs = 1000;
  pad = 4;
  poly_tol = 0.1;
  bwmesh_tol = 1.0;
  bwmesh_smoothing_iters = 10;

  % Just use distance to a well sampled poly-line
  [V,E] = spline_to_poly(P,C,poly_tol);
  [BC,side,r] = voxel_grid(P,gs,'Pad',pad);
  X = reshape(BC(:,1),side([2 1]));
  Y = reshape(BC(:,2),side([2 1]));
  sqrD = point_mesh_squared_distance(BC,V,E);
  D = sqrt(reshape(sqrD,size(X)));
  % Find maximum distance sampled at control points to try to be sure they're
  % inside the extracted mesh
  Dmax = max(interp2(X,Y,reshape(D,size(X)),P(:,1),P(:,2)));
  % convert into fuzzy logistic looking function
  h = 2/sqrt(2)*r(1);
  Dmax = Dmax+h;
  Z = 1-min(max(D-Dmax,0),h)/h;
  % Ignore holes
  Z = imfill(Z);
  [CV,CF] = bwmesh(Z,'Tol',bwmesh_tol,'SmoothingIters',bwmesh_smoothing_iters);
  % bwmesh flips things around and doesn't know the scale
  CV(:,2) = size(Z,1)-CV(:,2);
  CV = (CV-0.0)/size(Z,2)*(max(X(:))-min(X(:)))+[min(X(:)) min(Y(:))];
end

author: alecjacobson

repo: gptoolbox


text:
function i4_to_pascal_test ( )

%*****************************************************************************80
%
%% I4_TO_PASCAL_TEST tests I4_TO_PASCAL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 April 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4_TO_PASCAL_TEST\n' );
  fprintf ( 1, '  I4_TO_PASCAL converts a linear index to\n' );
  fprintf ( 1, '  Pascal triangle indices.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     K  =>   I     J\n' );
  fprintf ( 1, '\n' );

  for k = 1 : 20

    [ i, j ] = i4_to_pascal ( k );

    fprintf ( 1, '  %4d    %4d  %4d\n', k, i, j );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function B=autocorr2d(H)

% Compute the 2D  autocorrelation of matrix ( Image )
% Algorithm based on Wiener - Khintchine Theorem*.
%
% * :http://mathworld.wolfram.com/Wiener-KhinchinTheorem.html
%
% Special Case :
% 2D zero mean Gaussian process with variance=1:
%           ==========================
%           | P=randn(400,600);      |
%           | CorrFunc=autocorr2d(H);|
%           ==========================
% CorrFunc is symmetric with Dirac impulsion in the center.
% with max(CorrFunc(:))=1
% July, 24, 2012
% KHMOU Youssef

[n m]=size(H);
% Divide by the size for normalization

B=abs(fftshift(ifft2(fft2(H).*conj(fft2(H)))))./(n*m);
figure, surf(B);
shading interp;

author: Sable

repo: mcbench-benchmarks


text:
function [model] = rvr_train(X, y, options)
% RVM_TRAIN Trains an RVM Model using SB1 Tipping Toolbox
%
%   input ----------------------------------------------------------------
%
%       o X        : (N x D), N  input data points of D dimensionality.
%
%       o y        : (N x 1), N  output data points
%
%       o options     : struct
%
%
%   output ----------------------------------------------------------------
%
%       o model       : struct.
%
%
%% %    RVM OPTIONS
%       ALPHA   Scalar initial value for hyperparameters
%       BETA    Initial value for inverse noise variance (in regression)
%               Set this negative to fix the value, rather than estimate
%       KERNEL  Kernel type: see SB1_KERNELFUNCTION for options
%       LEN     Kernel length scale
%       USEBIAS Set to non-zero to utilise a "bias" offset
%       MAXITS  Maximum iterations to run for.

%% Parse RVM Options

% Transform Data to Columns
X = X(:);
y = y(:);

% Parsing Parameter for RVR
N	= length(X);
useBias = options.useBias;
kernel_	= options.kernel_;
width   = options.width;
maxIts	= options.maxIts;
monIts	= round(maxIts/20);

% Set verbosity of output (0 to 4)
setEnvironment('Diagnostic','verbosity',3);
% Set file ID to write to (1 = stdout)
setEnvironment('Diagnostic','fid',1);

%% Set up initial hyperparameters - precise settings should not be critical
%
initAlpha	= (1/N)^2;
epsilon		= std(y) * 10/100; % Initial guess of 10% noise-to-signal
initBeta	= 1/epsilon^2;

%% Train RVR Model

% "Train" a sparse Bayes kernel-based model (relevance vector machine)
[weights, used, bias, marginal, alpha, beta, gamma] = ...
    SB1_RVM(X,y,initAlpha,initBeta,kernel_,width,useBias,maxIts,monIts);

%% Model Parameters
%       WEIGHTS Parameter values of estimated model (sparse)
%       USED    Index vector of "relevant" kernels (data points)
%       BIAS    Value of bias or offset parameter
%       ML      Log marginal likelihood of model
%       ALPHA   Estimated hyperparameter values (sparse)
%       BETA    Estimated inverse noise variance for regression
%       GAMMA   "Well-determinedness" factors for relevant kernels

model.weights  = weights;
model.kernel_  = kernel_;
model.width    = width;
model.RVs_idx  = used;
model.bias     = bias;
model.marginal = marginal;
model.alpha    = alpha;
model.beta     = beta;
model.gamma    = gamma;
model.RVs      = X(used,:);

end

author: epfl-lasa


text:
% HYBRJ  Solve a SNLE using HYBRJ
%
% THIS IS A LOW LEVEL FUNCTION - USE opti_hybrj() INSTEAD!
%
% hybrj uses the Minpack Powell Hybrid library.
%
%   [x,fval,exitflag,fevals] = hybrj(fun,grad,x0,opts)
%
%   Input arguments:
%       fun - vector of nonlinear functions (equations to be solved)
%       grad - gradient of nonlinear equations (optional)
%       x0 - initial solution guess
%       opts - solver options (see below)
%
%   Return arguments:
%       x - solution vector
%       fval - objective value at the solution
%       exitflag - exit status (see below)
%       fevals - number of function evaluations taken
%
%   Option Fields (all optional):
%       display - solver display level [0,1,2]
%       maxfeval - maximum function evaluations
%       maxtime - maximum solver execution time
%       iterfun - Iteration Callback Function, stop = iterfun(iter,fval,x)
%
%   Return Status:
%       1 - optimal
%       0 - maximum iterations exceeded
%      -1 - infeasible / could not converge
%      -2 - hybrj error / other
%      -5 - user exit
%
%
%   Copyright (C) 2011 Jonathan Currie (I2C2)

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function outn=le(x,y)

precAndSize

for ii=1:max(ex,ey)
 imag=false;
 if ex==1
  [xrval,xival]=getVals(x,1);
  [yrval,yival]=getVals(y,ii);
 elseif ey==1
  [xrval,xival]=getVals(x,ii);
  [yrval,yival]=getVals(y,1);
 else
  [xrval,xival]=getVals(x,ii);
  [yrval,yival]=getVals(y,ii);
 end
  outn(ii)=mpfr_le(precision,xrval,yrval);
  if outn(ii)~=0, outn(ii)=1; end
end % for ii=1:max(ex,

author: opencobra

repo: cobratoolbox


text:
function [y,nz] = ompdenoise1(params,msgdelta)
%OMPDENOISE1 OMP denoising of 1-D signals.
%  OMPDENOISE1 denoises a 1-dimensional signal using OMP denoising. The
%  function syntax is identical to OMPDENOISE, but it runs significantly
%  faster on 1-D signals. OMPDENOISE1 requires somewhat more memory than
%  OMPDENOISE (approximately the size of the input signal), so if memory is
%  limited, OMPDENOISE can be used instead.
%
%  See also OMPDENOISE.

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

% parse input arguments %

x = params.x(:);
D = params.dict;
blocksize = params.blocksize;

% maxval %
if (isfield(params,'maxval'))
  maxval = params.maxval;
else
  maxval = 1;
end

% gain %
if (isfield(params,'gain'))
  gain = params.gain;
else
  gain = 1.15;
end

% maxatoms %
if (isfield(params,'maxatoms'))
  maxatoms = params.maxatoms;
else
  maxatoms = floor(prod(blocksize)/2);
end

% stepsize %
if (isfield(params,'stepsize'))
  stepsize = params.stepsize;
else
  stepsize = 1;
end
if (any(stepsize<1))
  error('Invalid step size.');
end

% noise mode %
if (isfield(params,'noisemode'))
  switch lower(params.noisemode)
    case 'psnr'
      sigma = maxval / 10^(params.psnr/20);
    case 'sigma'
      sigma = params.sigma;
    otherwise
      error('Invalid noise mode specified');
  end
elseif (isfield(params,'sigma'))
  sigma = params.sigma;
elseif (isfield(params,'psnr'))
  sigma = maxval / 10^(params.psnr/20);
else
  error('Noise strength not specified');
end

% lambda %
if (isfield(params,'lambda'))
  lambda = params.lambda;
else
  lambda = maxval/(10*sigma);
end

% msgdelta %
if (nargin <2)
  msgdelta = 5;
end
if (msgdelta<=0)
  msgdelta = -1;
end

epsilon = sqrt(prod(blocksize)) * sigma * gain;   % target error for omp

MEM_LOW = 1;
MEM_NORMAL = 2;
MEM_HIGH = 3;

if (isfield(params,'memusage'))
  switch lower(params.memusage)
    case 'low'
      memusage = MEM_LOW;
    case 'normal'
      memusage = MEM_NORMAL;
    case 'high'
      memusage = MEM_HIGH;
    otherwise
      error('Invalid memory usage mode');
  end
else
  memusage = MEM_NORMAL;
end

% denoise the signal %

if (memusage >= MEM_NORMAL)
  G = D'*D;
end

% process the signal in batches to conserve memory
% choose batchsize so im2col returns a matrix of approximately the same
% size as the signal
batchsize = ceil(length(x)*stepsize/blocksize + blocksize);

y = zeros(size(x));
ids = 1:min(batchsize,length(x));
nz = 0;

tid = timerinit('ompdenoise', length(x));
while (length(ids)>=blocksize)

  % extract the signal blocks
  blocks = im2col(x(ids),[blocksize 1],'sliding');
  blocks = blocks(:,1:stepsize:end);

  % remove DC
  [blocks, dc] = remove_dc(blocks,'columns');

  % denoise the blocks
  if (memusage == MEM_LOW)
    gamma = omp2(D,blocks,[],epsilon,'maxatoms',maxatoms);
  else
    gamma = omp2(D,blocks,G,epsilon,'maxatoms',maxatoms);
  end
  nz = nz + nnz(gamma);
  cleanblocks = add_dc(D*gamma, dc, 'columns');

  y(ids) = y(ids) + col2imsum(cleanblocks, [blocksize 1], [length(ids) 1], [stepsize 1]);
  ids = ids + floor((batchsize-blocksize)/stepsize)*stepsize + stepsize;
  if (ids(end)>length(x))
    ids = ids(ids<=length(x));
  end

  % display status
  if (msgdelta>0 && ~isempty(ids)>0)
    timereta(tid, ids(1), msgdelta);
  end

end

if (msgdelta>0)
  timereta(tid, length(x));
end

cnt = countcover(size(x),[blocksize 1],[stepsize 1]);
y = (y+lambda*x)./(cnt + lambda);
y = reshape(y,size(params.x));

author: thomas-koehler

repo: SupER


text:
%% Machine Learning Online Class
%  Exercise 1: Linear regression with multiple variables
%
%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  linear regression exercise.
%
%  You will need to complete the following functions in this
%  exericse:
%
%     warmUpExercise.m
%     plotData.m
%     gradientDescent.m
%     computeCost.m
%     gradientDescentMulti.m
%     computeCostMulti.m
%     featureNormalize.m
%     normalEqn.m
%
%  For this part of the exercise, you will need to change some
%  parts of the code below for various experiments (e.g., changing
%  learning rates).
%

%% Initialization

%% ================ Part 1: Feature Normalization ================

%% Clear and Close Figures
clear all; close all; clc

fprintf('Loading data ...\n');

%% Load Data
data = csvread('ex1data2.txt');
X = data(:, 1:2);
y = data(:, 3);
m = length(y);

% Print out some data points
fprintf('First 10 examples from the dataset: \n');
fprintf(' x = [%.0f %.0f], y = %.0f \n', [X(1:10,:) y(1:10,:)]');

fprintf('Program paused. Press enter to continue.\n');
pause;

% Scale features and set them to zero mean
fprintf('Normalizing Features ...\n');

[X mu sigma] = featureNormalize(X);

% Add intercept term to X
X = [ones(m, 1) X];

%% ================ Part 2: Gradient Descent ================

% ====================== YOUR CODE HERE ======================
% Instructions: We have provided you with the following starter
%               code that runs gradient descent with a particular
%               learning rate (alpha).
%
%               Your task is to first make sure that your functions -
%               computeCost and gradientDescent already work with
%               this starter code and support multiple variables.
%
%               After that, try running gradient descent with
%               different values of alpha and see which one gives
%               you the best result.
%
%               Finally, you should complete the code at the end
%               to predict the price of a 1650 sq-ft, 3 br house.
%
% Hint: By using the 'hold on' command, you can plot multiple
%       graphs on the same figure.
%
% Hint: At prediction, make sure you do the same feature normalization.
%

fprintf('Running gradient descent ...\n');

% Choose some alpha value
alpha = 0.3;
num_iters = 100;

% Init Theta and Run Gradient Descent
theta = zeros(3, 1);
[theta, J_history] = gradientDescentMulti(X, y, theta, alpha, num_iters);

% Plot the convergence graph
figure;
plot(1:numel(J_history), J_history, '-b', 'LineWidth', 2);
xlabel('Number of iterations');
ylabel('Cost J');

% Display gradient descent's result
fprintf('Theta computed from gradient descent: \n');
fprintf(' %f \n', theta);
fprintf('\n');

% Estimate the price of a 1650 sq-ft, 3 br house
% ====================== YOUR CODE HERE ======================
% Recall that the first column of X is all-ones. Thus, it does
% not need to be normalized.
price = 0; % You should change this

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

fprintf(['Predicted price of a 1650 sq-ft, 3 br house ' ...
         '(using gradient descent):\n $%f\n'], price);

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ================ Part 3: Normal Equations ================

fprintf('Solving with normal equations...\n');

% ====================== YOUR CODE HERE ======================
% Instructions: The following code computes the closed form
%               solution for linear regression using the normal
%               equations. You should complete the code in
%               normalEqn.m
%
%               After doing so, you should complete this code
%               to predict the price of a 1650 sq-ft, 3 br house.
%

%% Load Data
data = csvread('ex1data2.txt');
X = data(:, 1:2);
y = data(:, 3);
m = length(y);

% Add intercept term to X
X = [ones(m, 1) X];

% Calculate the parameters from the normal equation
theta = normalEqn(X, y);

% Display normal equation's result
fprintf('Theta computed from the normal equations: \n');
fprintf(' %f \n', theta);
fprintf('\n');

% Estimate the price of a 1650 sq-ft, 3 br house
% ====================== YOUR CODE HERE ======================
price = 0; % You should change this

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

fprintf(['Predicted price of a 1650 sq-ft, 3 br house ' ...
         '(using normal equations):\n $%f\n'], price);

author: zhouxc

repo: Stanford-Machine-Learning-Course


text:
function mprob = convBonmin(prob,opts,warn)
%CONVBONMIN Convert OPTI problem to BONMIN Mixed Integer Nonlinear problem
%
%   mprob = convBonmin(prob,opts)

%   Copyright (C) 2011 Jonathan Currie (I2C2)

%COMPATIBLE CPLEX VERSION
CPLEXVER = '12.5.0.0';

%Ensure all args passed
if(nargin < 2)
    error('You must supply both the problem + options');
end
if(~isstruct(prob) || ~isstruct(opts))
    error('Both prob + opts must be structures!');
end

%Ensure we have defaults for bonmin
if(isempty(opts.solverOpts))
    opts.solverOpts = bonminset;
end

%First convert to IPOPT problem (used as the internal relaxed solver)
mprob = convIpopt(prob,opts);
%Remove options from IPOPT structure we don't need
mprob.options.ipopt = rmfield(mprob.options.ipopt,'max_cpu_time'); %don't want this set in the local solver
mprob.options.ipopt = rmfield(mprob.options.ipopt,'print_level'); %don't want this set in the local solver
%Remove options already at OPTI default (appears BONMIN sets some internally different from docs)
mprob.options.ipopt = optiRemoveDefaults(mprob.options.ipopt,ipoptset('bonmin'));

%Set BONMIN options
mprob.options.bonmin = bonminset(opts.solverOpts,'noIpopt'); %this will remove ipopt options from the bonminset structure
%Remove options set elsewhere
mprob.options.bonmin = rmfield(mprob.options.bonmin,'var_lin');
mprob.options.bonmin = rmfield(mprob.options.bonmin,'cons_lin');
%Remove options already at OPTI default (appears BONMIN sets some internally different from docs)
mprob.options.bonmin = optiRemoveDefaults(mprob.options.bonmin,bonminset());

%If user has specified Cplex, ensure we can load it
if(isfield(mprob.options.bonmin,'milp_solver') && strcmpi(mprob.options.bonmin.milp_solver,'cplex'))
    try
        %Forces to load Cplex on current platform (x86 or x64) if MATLAB interface present
        c = Cplex; %#ok<NASGU>
    catch %#ok<CTCH>
       %nothing to do
    end
    try
        %The following will fail with a "cannot load module" if it cannot find correct version of Cplex DLL on user's PC
        v = bonminCplex(); %#ok<NASGU>
    catch %#ok<CTCH>
        if(warn)
            optiwarn('opti:bonmin_Cplex','Could not load BONMIN with CPLEX - ensure CPLEX (v%s) is installed and licensed on your PC. Using CBC for now.',CPLEXVER);
        end
        mprob.options.bonmin = rmfield(mprob.options.bonmin,'milp_solver');
    end
end

%Convert OPTISET options to equivalent BONMIN ones
mprob.options.bonmin.node_limit = opts.maxnodes;
mprob.options.bonmin.time_limit = opts.maxtime;
mprob.options.bonmin.integer_tolerance = opts.tolint;

%Set Display Level
if(strcmpi(opts.display,'iter'))
    mprob.options.display = 1;
else
    mprob.options.display = 0;
end

%Setup Integer Vars
mprob.options.var_type = prob.int.ind; %uses -1, 0, 1 format
%Setup Nonlinear Vars
if(isempty(opts.solverOpts.var_lin) || any(isnan(opts.solverOpts.var_lin)))
    mprob.options.var_lin = zeros(prob.sizes.ndec,1); %for now assuming all nonlinear
else
    varlin = opts.solverOpts.var_lin;
    if(length(varlin) ~= prob.sizes.ndec)
        error('The decision variable linearity vector is the wrong size! Expected %d x 1',prob.sizes.ndec);
    else
        mprob.options.var_lin = varlin;
    end
end
%Setup Nonlinear Constraints Linearity
if(prob.sizes.nrow) %bug fix when bonmin solves with double-sided linear row constraints
    nlincon = prob.sizes.nrow;
else
    nlincon = prob.sizes.nineq + prob.sizes.neq;
end
if(prob.sizes.nnlrow) %bug fix when bonmin solves with double-sided nonlinear row constraints
    nnlcon = prob.sizes.nnlrow + prob.sizes.nqc;
else
    nnlcon = prob.sizes.nnlineq + prob.sizes.nnleq + prob.sizes.nqc;
end
if(isempty(opts.solverOpts.cons_lin) || any(isnan(opts.solverOpts.cons_lin)))
   mprob.options.cons_lin = [zeros(nnlcon,1); ones(nlincon,1)]; %augment both kinds
else
    conslin = opts.solverOpts.cons_lin;
    if(length(conslin) ~= (nlincon+nnlcon))
        error('The constraint linearity vector is the wrong size! Expected %d x 1',nlincon+nnlcon);
    else
        mprob.options.cons_lin = conslin;
    end
end


text:
%*****************THE MAIN PROGRAMME*****************
%====================================================
%Practise problem for the course of CAO
%Prof P.Beckers
%====================================================
%Porpose :Drawing the orthotomic surface of a Bezier surface  %
% Student: BUI QUOC TINH
% European Master Mechanices of Contructions	(EMMC)
% University the Liege. Belgium										    		     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%====================================================
clear all;

%Definition of the control points directly
P00=[-10 0 10];P01=[-10 5 5];P02=[-10 5 -5];P03=[-10 0 -10];
P10=[-5 5 10];P11=[-5 5 5];P12=[-5 5 -5];P13=[-5 5 -10];
P20=[5 5 10];P21=[5 5 5];P22=[5 5 -5];P23=[5 5 -10];
P30=[10 0 10];P31=[10 5 5];P32=[10 5 -5];P33=[10 5 -10];

%Matrix control points of Bezier surface
PP=[P00 P01 P02 P03;...
    P10 P11 P12 P13;...
    P20 P21 P22 P23;...
    P30 P31 P32 P33];

 %S source point
%disp('INPUT S SOURCE POINT ')
%xs=input('input xs = ');
%ys=input('input ys = ');
%zs=input('input zs = ');
xs=-5;
ys=3;
zs=5;
S=[xs;ys;zs];

%Calculation of the points in the surface
k=20;
for i=1:k
	if i>9 & mod(i,10)==0
		disp(i) % Using show is the progress of the calculations.
	end
for j=1:k
u=0.01*i;
v=0.01*j;

%=============================================
%the function of bezier surafce polynomial
q=P00*B(0,3,u)*B(0,3,v)+P01*B(0,3,u)*B(1,3,v)+P02*B(0,3,u)*B(2,3,v)+P03*B(0,3,u)*B(3,3,v)+...
+P10*B(1,3,u)*B(0,3,v)+P11*B(1,3,u)*B(1,3,v)+P12*B(1,3,u)*B(2,3,v)+P13*B(1,3,u)*B(3,3,v)+...
+P20*B(2,3,u)*B(0,3,v)+P21*B(2,3,u)*B(1,3,v)+P22*B(2,3,u)*B(2,3,v)+P23*B(2,3,u)*B(3,3,v)+...
+P30*B(3,3,u)*B(0,3,v)+P31*B(3,3,u)*B(1,3,v)+P32*B(3,3,u)*B(2,3,v)+P33*B(3,3,u)*B(3,3,v);

%the first derivatives u parametic of bezier surface function
q_u=P00*B_daoham(0,3,u)*B(0,3,v)+P01*B_daoham(0,3,u)*B(1,3,v)+P02*B_daoham(0,3,u)*B(2,3,v)+P03*B_daoham(0,3,u)*B(3,3,v)+...
+P10*B_daoham(1,3,u)*B(0,3,v)+P11*B_daoham(1,3,u)*B(1,3,v)+P12*B_daoham(1,3,u)*B(2,3,v)+P13*B_daoham(1,3,u)*B(3,3,v)+...
+P20*B_daoham(2,3,u)*B(0,3,v)+P21*B_daoham(2,3,u)*B(1,3,v)+P22*B_daoham(2,3,u)*B(2,3,v)+P23*B_daoham(2,3,u)*B(3,3,v)+...
+P30*B_daoham(3,3,u)*B(0,3,v)+P31*B_daoham(3,3,u)*B(1,3,v)+P32*B_daoham(3,3,u)*B(2,3,v)+P33*B_daoham(3,3,u)*B(3,3,v);

%%the first derivatives v parametic of bezier surface function
q_v=P00*B(0,3,u)*B_daoham(0,3,v)+P01*B(0,3,u)*B_daoham(1,3,v)+P02*B(0,3,u)*B_daoham(2,3,v)+P03*B(0,3,u)*B_daoham(3,3,v)+...
+P10*B(1,3,u)*B_daoham(0,3,v)+P11*B(1,3,u)*B_daoham(1,3,v)+P12*B(1,3,u)*B_daoham(2,3,v)+P13*B(1,3,u)*B_daoham(3,3,v)+...
+P20*B(2,3,u)*B_daoham(0,3,v)+P21*B(2,3,u)*B_daoham(1,3,v)+P22*B(2,3,u)*B_daoham(2,3,v)+P23*B(2,3,u)*B_daoham(3,3,v)+...
+P30*B(3,3,u)*B_daoham(0,3,v)+P31*B(3,3,u)*B_daoham(1,3,v)+P32*B(3,3,u)*B_daoham(2,3,v)+P33*B(3,3,u)*B_daoham(3,3,v);

%===============
qx(i,j)=q(1);%Saving points matrix into the remember computer
qy(i,j)=q(2);
qz(i,j)=q(3);

%===============
qx_u=q_u(1);
qy_u=q_u(2);
qz_u=q_u(3);

%===============
qx_v=q_v(1);
qy_v=q_v(2);
qz_v=q_v(3);

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

%The geomatric information matrix of tangent plane
R=[qx(i,j) qx_u qx_v;...
  qy(i,j) qy_u qy_v;...
  qz(i,j) qz_u qz_v];

%The normal vector of tangent plane as same as the orient vector of a traigh lines
%perpendicular with it and through source point
nn=[(qy_u.*qz_v)-(qz_u.*qy_v);...
    (qz_u.*qx_v)-(qz_v.*qx_v);...
    (qx_u.*qy_v)-(qx_v.*qy_u)];

%The geomatric information matrix of straigh lines d perpendicular with it and through source point
d=[xs nn(1);...
   ys nn(2);...
   zs nn(3)];

%Solving the equations
AA=[nn(1) -qx_u -qx_v;...
   nn(2) -qy_u -qy_v;...
   nn(3) -qz_u -qz_v];
syms t1 t2 t3 % parametric variables
T=[t1; t2; t3];
BB=[(qx(i,j)-xs);...
    (qy(i,j)-ys);...
    (qz(i,j)-zs)];
T=inv(AA)*BB;

%The intersection points between the tangent plane and the straigh d
H=S+T(1)*nn;

%The reflective points of the source point about the tangent plane
S_dx=2*H-S;
%==================
S_dxx(i,j)=S_dx(1);
S_dxy(i,j)=S_dx(2);
S_dxz(i,j)=S_dx(3);
end
end

% Drawing the bezier surface
view(3);
grid
surface(qx,qy,qz); xlabel('x');ylabel('y');zlabel('z');%
title('The orthotomic surface of a bezier surface');

%Drawing the source point
hold on
plot3(xs,ys,zs,'Pr','ButtonDownFcn','animator start','EraseMode','xor','LineWidth',2, ...
'Marker','>','MarkerSize',5,'Tag','pointer');
text(xs+0.2,ys+0.2,zs+0.2,'S')
hold on

%Drawing the orthotomic surface
surface(S_dxx,S_dxy,S_dxz); xlabel('x');ylabel('y');zlabel('z');

% Thank you very much for your course !


text:
function fire_serial ( forest_size, prob_spread )

%*****************************************************************************80
%
%% MAIN is the main program for FIRE_SERIAL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 July 2013.
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, integer FOREST_SIZE, specifies the forest as a square of
%    FOREST_SIZE by FOREST_SIZE trees.  Default is 20.
%
%    Input, real PROB_SPREAD, the probability that fire will spread from
%    a burning tree to its neighbor.  Default is 0.5.
%
  if ( nargin < 1 )
    forest_size = 20;
  elseif ( ischar ( forest_size ) )
    forest_size = str2num ( forest_size );
  end

  if ( nargin < 2 )
    prob_spread = 0.5;
  elseif ( ischar ( prob_spread ) )
    prob_spread = str2num ( prob_spread );
  end

  SMOLDERING = 1;

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FIRE_SERIAL\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  A probabilistic simulation of a forest fire.\n' );
  fprintf ( 1, '  The forest size is %d\n', forest_size );
  fprintf ( 1, '  The probability of tree-to-tree spread is %g\n', prob_spread );
%
%  Initialize the random number generator.
%
  seed = get_seed ( );
  fprintf ( 1, '  The random number generator is seeded by %d\n', seed );
%
%  Initialize the values in the forest.
%
  forest = forest_initialize ( forest_size );
%
%  Choose a tree at random where the fire will start.
%
  [ i_ignite, seed ] = i4_uniform_ab ( 1, forest_size, seed );
  [ j_ignite, seed ] = i4_uniform_ab ( 1, forest_size, seed );
  forest(i_ignite,j_ignite) = SMOLDERING;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Fire starts at tree(%d,%d)\n', i_ignite, j_ignite );
%
%  Let time run until nothing is burning any more.
%
  while ( forest_is_burning ( forest_size, forest ) )
    [ forest, seed ] = forest_burns ( forest_size, forest, seed, prob_spread );
  end
%
%  Display the final forest state.
%
  forest_print ( forest_size, forest, i_ignite, j_ignite );
%
%  Report the percentage of forest burned.
%
  percent = get_percent_burned ( forest_size, forest );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Percentage of forest burned = %g\n', percent );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FIRE_SERIAL:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function [ value, seed ] = fire_spreads ( seed, prob_spread )

%*****************************************************************************80
%
%% FIRE_SPREADS determines whether the fire spreads.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 July 2013
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input/output, integer SEED, a seed for the random
%    number generator.
%
%    Input, real PROB_SPREAD, the probability of spreading.
%
%    Output, logical VALUE, is TRUE if the fire spreads.
%
  [ u, seed ] = r8_uniform_01 ( seed );

  if ( u < prob_spread )
    value = 1;
  else
    value = 0;
  end

  return
end
function [ forest, seed ] = forest_burns ( forest_size, forest, seed, ...
  prob_spread )

%*****************************************************************************80
%
%% FOREST_BURNS models a single time step of the burning forest.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 July 2013
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, integer FOREST_SIZE, the linear dimension of the forest.
%
%    Input/output, integer FOREST(FOREST_SIZE,FOREST_SIZE), an
%    array with an entry for each tree in the forest.
%
%    Input/output, integer SEED, a seed for the random
%    number generator.
%
%    Input, real PROB_SPREAD, the probability that the fire will
%    spread from a burning tree to an unburnt one.
%
  BURNING = 2;
  BURNT = 3;
  SMOLDERING = 1;
  UNBURNT = 0;
%
%  Burning trees burn down;
%  Smoldering trees ignite;
%
  for j = 1 : forest_size
    for i = 1 : forest_size
      if ( forest(i,j) == BURNING )
        forest(i,j) = BURNT;
      elseif ( forest(i,j) == SMOLDERING )
        forest(i,j) = BURNING;
      end
    end
  end
%
%  Unburnt trees might catch fire.
%
  for j = 1 : forest_size
    for i = 1 : forest_size

      if ( forest(i,j) == BURNING )
%
%  North.
%
        if ( 1 < i )
          if ( forest(i-1,j) == UNBURNT )
            [ value, seed ] = fire_spreads ( seed, prob_spread );
            if ( value )
              forest(i-1,j) = SMOLDERING;
            end
          end
        end
%
%  South.
%
        if ( i < forest_size )
          if ( forest(i+1,j) == UNBURNT )
            [ value, seed ] = fire_spreads ( seed, prob_spread );
            if ( value )
              forest(i+1,j) = SMOLDERING;
            end
          end
        end
%
%  West.
%
        if ( 1 < j )
          if ( forest(i,j-1) == UNBURNT )
            [ value, seed ] = fire_spreads ( seed, prob_spread );
            if ( value )
              forest(i,j-1) = SMOLDERING;
            end
          end
        end
%
%  East.
%
        if ( j < forest_size )
          if ( forest(i,j+1) == UNBURNT )
            [ value, seed ] = fire_spreads ( seed, prob_spread );
            if ( value )
              forest(i,j+1) = SMOLDERING;
            end
          end
        end

      end

    end
  end

  return
end
function forest = forest_initialize ( forest_size )

%*****************************************************************************80
%
%% FOREST_INITIALIZE initializes the forest values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 July 2013
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, integer FOREST_SIZE, the linear dimension of the forest.
%
%    Output, integer FOREST(FOREST_SIZE,FOREST_SIZE), an array
%    with an entry for each tree in the forest.
%
  UNBURNT = 0;

  forest(1:forest_size,1:forest_size) = UNBURNT;

  return
end
function value = forest_is_burning ( forest_size, forest )

%*****************************************************************************80
%
%% FOREST_IS_BURNING reports whether any trees in the forest are burning.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 July 2013
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, integer FOREST_SIZE, the linear dimension of the forest.
%
%    Input, integer FOREST(FOREST_SIZE,FOREST_SIZE), an array
%    with an entry for each tree in the forest.
%
%    Output, logical FOREST_IS_BURNING, is TRUE if any tree in the forest
%    is in the SMOLDERING or BURNING state.
%
  BURNING = 2;
  SMOLDERING = 1;

  value = 0;

  for j = 1 : forest_size
    for i = 1 : forest_size
      if ( forest(i,j) == SMOLDERING || forest(i,j) == BURNING )
        value = 1;
      end
    end
  end

  return
end
function forest_print ( forest_size, forest, i_ignite, j_ignite )

%*****************************************************************************80
%
%% FOREST_PRINT prints the state of the trees in the forest.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 July 2013
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, integer FOREST_SIZE, the linear dimension of the forest.
%
%    Input, integer FOREST(FOREST_SIZE,FOREST_SIZE), an array
%    with an entry for each tree in the forest.
%
%    Input, integer I_IGNITE, J_IGNITE, the location of the start
%    of the fire.
%
  BURNT = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Map of fire damage.\n' );
  fprintf ( 1, '  Fire started at "*".\n' );
  fprintf ( 1, '  Burned trees are indicated by ".".\n' );
  fprintf ( 1, '  Unburned trees are indicated by "X".\n' );
  fprintf ( 1, '\n' );

  for i = 1 : forest_size
    fprintf ( 1, '  ' );
    for j = 1 : forest_size
      if ( i == i_ignite && j == j_ignite )
        fprintf ( 1, '*' );
      elseif ( forest(i,j) == BURNT )
        fprintf ( 1, '.' );
      else
        fprintf ( 1, 'X' );
      end
    end
    fprintf ( 1, '\n' );
  end

  return
end
function percent = get_percent_burned ( forest_size, forest )

%*****************************************************************************80
%
%% GET_PERCENT_BURNED computes the percentage of the forest that burned.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 July 2013
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, integer FOREST_SIZE, the linear dimension of the forest.
%
%   Input, integer FOREST(FOREST_SIZE,FOREST_SIZE), an array
%    with an entry for each tree in the forest.
%
%    Output, real PERCENT, the percentage of the forest
%    that burned.
%
  BURNT = 3;

  total = 0;
  for j = 1 : forest_size
    for i = 1 : forest_size
      if ( forest(i,j) == BURNT )
        total = total + 1;
      end
    end
  end

  percent = total / forest_size / forest_size;

  return
end
function seed = get_seed ( )

%*****************************************************************************80
%
%% GET_SEED returns a seed for the random number generator.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 November 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer SEED, a random seed value.
%
  I_MAX = 2147483647;

  time_array = clock;

  hour = time_array(4);
  minute = time_array(5);
  second = time_array(6);

  temp = ( second + 60 * ( minute + 60 * hour ) ) / ( 60.0 * 60.0 * 24.0 );

  if ( temp <= 0.0 )
    temp = temp + 1.0;
  end

  if ( 1.0 < temp )
    temp = temp - 1.0;
  end

  seed = 1 + floor ( I_MAX * temp );

  return
end
function [ c, seed ] = i4_uniform_ab ( a, b, seed )

%*****************************************************************************80
%
%% I4_UNIFORM_AB returns a scaled pseudorandom I4.
%
%  Discussion:
%
%    The pseudorandom number will be scaled to be uniformly distributed
%    between A and B.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Second Edition,
%    Springer, 1987,
%    ISBN: 0387964673,
%    LC: QA76.9.C65.B73.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, December 1986, pages 362-376.
%
%    Pierre L'Ecuyer,
%    Random Number Generation,
%    in Handbook of Simulation,
%    edited by Jerry Banks,
%    Wiley, 1998,
%    ISBN: 0471134031,
%    LC: T57.62.H37.
%
%    Peter Lewis, Allen Goodman, James Miller,
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, Number 2, 1969, pages 136-143.
%
%  Parameters:
%
%    Input, integer A, B, the minimum and maximum acceptable values.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, integer C, the randomly chosen integer.
%
%    Output, integer SEED, the updated seed.
%
  i4_huge = 2147483647;

  if ( seed == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4_UNIFORM_AB - Fatal error!\n' );
    fprintf ( 1, '  Input SEED = 0!\n' );
    error ( 'I4_UNIFORM_AB - Fatal error!' );
  end

  seed = floor ( seed );
  a = round ( a );
  b = round ( b );

  seed = mod ( seed, i4_huge );

  if ( seed < 0 )
    seed = seed + i4_huge;
  end

  k = floor ( seed / 127773 );

  seed = 16807 * ( seed - k * 127773 ) - k * 2836;

  if ( seed < 0 )
    seed = seed + i4_huge;
  end

  r = seed * 4.656612875E-10;
%
%  Scale R to lie between A-0.5 and B+0.5.
%
  r = ( 1.0 - r ) * ( min ( a, b ) - 0.5 ) ...
    +         r   * ( max ( a, b ) + 0.5 );
%
%  Use rounding to convert R to an integer between A and B.
%
  value = round ( r );

  value = max ( value, min ( a, b ) );
  value = min ( value, max ( a, b ) );

  c = value;

  return
end
function [ r, seed ] = r8_uniform_01 ( seed )

%*****************************************************************************80
%
%% R8_UNIFORM_01 returns a unit pseudorandom R8.
%
%  Discussion:
%
%    This routine implements the recursion
%
%      seed = 16807 * seed mod ( 2^31 - 1 )
%      r8_uniform_01 = seed / ( 2^31 - 1 )
%
%    The integer arithmetic never requires more than 32 bits,
%    including a sign bit.
%
%    If the initial seed is 12345, then the first three computations are
%
%      Input     Output      R8_UNIFORM_01
%      SEED      SEED
%
%         12345   207482415  0.096616
%     207482415  1790989824  0.833995
%    1790989824  2035175616  0.947702
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Springer Verlag, pages 201-202, 1983.
%
%    Pierre L'Ecuyer,
%    Random Number Generation,
%    in Handbook of Simulation,
%    edited by Jerry Banks,
%    Wiley Interscience, page 95, 1998.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, pages 362-376, 1986.
%
%    Peter Lewis, Allen Goodman, James Miller,
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, pages 136-143, 1969.
%
%  Parameters:
%
%    Input, integer SEED, the integer "seed" used to generate
%    the output random number.  SEED should not be 0.
%
%    Output, real R, a random value between 0 and 1.
%
%    Output, integer SEED, the updated seed.  This would
%    normally be used as the input seed on the next call.
%
  if ( seed == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_UNIFORM_01 - Fatal error!\n' );
    fprintf ( 1, '  Input SEED = 0!\n' );
    error ( 'R8_UNIFORM_01 - Fatal error!' );
  end

  seed = floor ( seed );

  seed = mod ( seed, 2147483647 );

  if ( seed < 0 )
    seed = seed + 2147483647;
  end

  k = floor ( seed / 127773 );

  seed = 16807 * ( seed - k * 127773 ) - k * 2836;

  if ( seed < 0 )
    seed = seed + 2147483647;
  end

  r = seed * 4.656612875E-10;

  return
end
function timestamp ( )

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

  return
end


text:
rng(1); % For reproducibility
r = sqrt(rand(100,1)); % Radius
t = 2*pi*rand(100,1);  % Angle
data1 = [r.*cos(t), r.*sin(t)]; % Points

r2 = sqrt(3*rand(100,1)+1); % Radius
t2 = 2*pi*rand(100,1);      % Angle
data2 = [r2.*cos(t2), r2.*sin(t2)]; % points

figure;
plot(data1(:,1),data1(:,2),'r.','MarkerSize',15)
hold on
plot(data2(:,1),data2(:,2),'b.','MarkerSize',15)
ezpolar(@(x)1);ezpolar(@(x)2);
axis equal
hold off

data3 = [data1;data2];
theclass = ones(200,1);
theclass(1:100) = -1;

%Train the SVM Classifier
cl = fitcsvm(data3,theclass,'KernelFunction','rbf',...
    'BoxConstraint',Inf,'ClassNames',[-1,1]);

% Predict scores over the grid
d = 0.02;
[x1Grid,x2Grid] = meshgrid(min(data3(:,1)):d:max(data3(:,1)),...
    min(data3(:,2)):d:max(data3(:,2)));
xGrid = [x1Grid(:),x2Grid(:)];
[~,scores] = predict(cl,xGrid);

% Plot the data and the decision boundary
figure;
h(1:2) = gscatter(data3(:,1),data3(:,2),theclass,'rb','.');
hold on
ezpolar(@(x)1);
h(3) = plot(data3(cl.IsSupportVector,1),data3(cl.IsSupportVector,2),'ko');
contour(x1Grid,x2Grid,reshape(scores(:,2),size(x1Grid)),[0 0],'k');
legend(h,{'-1','+1','Support Vectors'});
axis equal
hold off

cl2 = fitcsvm(data3,theclass,'KernelFunction','rbf');
[~,scores2] = predict(cl2,xGrid);

figure;
h(1:2) = gscatter(data3(:,1),data3(:,2),theclass,'rb','.');
hold on
ezpolar(@(x)1);
h(3) = plot(data3(cl2.IsSupportVector,1),data3(cl2.IsSupportVector,2),'ko');
contour(x1Grid,x2Grid,reshape(scores2(:,2),size(x1Grid)),[0 0],'k');
legend(h,{'-1','+1','Support Vectors'});
axis equal
hold off

author: Grootzz

repo: GLCM-SVM

save_path: github-repos/MATLAB/Grootzz-GLCM-SVM


text:
function Dist = sc_pdist2( X, Y, metric )
% Calculates the pairwise distance between sets of vectors.
%
% Let X be an D-by-M matrix representing m points in D-dimensional space
% and Y be an D-by-N matrix representing another set of points in the same
% space. This function computes the M-by-N distance matrix Dist where Dist(i,j)
% is the distance between X(:,j) and Y(:,j).  This function has been
% optimized where possible, with most of the distance computations
% requiring few or no loops.
%
% The metric can be one of the following:
%
% 'euclidean' / 'sqeuclidean':
%   Euclidean / SQUARED Euclidean distance.  Note that 'sqeuclidean'
%   is significantly faster.
%
% 'chisq' 'cs'
%   The chi-squared distance between two vectors is defined as:
%    d(x,y) = sum( (xi-yi)^2 / (xi+yi) ) / 2;
%   The chi-squared distance is useful when comparing histograms.
%
%  'hell' 'hik'
%   The hellinger's distance is dedined as:
%    d(x,y) = sum (x^.5 - y^.5) .^ 2
%
% 'cos'
%   Distance is defined as the cosine of the angle between two vectors.
%
% 'emd'
%   Earth Mover's Distance (EMD) between positive vectors (histograms).
%   Note for 1D, with all histograms having equal weight, there is a simple
%   closed form for the calculation of the EMD.  The EMD between histograms
%   x and y is given by the sum(abs(cdf(x)-cdf(y))), where cdf is the
%   cumulative distribution function (computed simply by cumsum).
%
% 'L1','LINF'
%   The L1 distance between two vectors is defined as:  sum(abs(x-y));
%   The L-Inf distance is defined as: max(abs(x-y));
%
% USAGE
%  Dist = sc_pdist2( X, Y, [metric] )
%
% INPUTS
%  X        - [D x M] matrix of M D-dimensional vectors
%  Y        - [D x N] matrix of N D-dimensional vectors
%  metric   - ['sqeuclidean'], 'chisq', 'cosine', 'emd', 'euclidean', 'L0'
%              'L1', 'Linf', 'hik'
%
% OUTPUTS
%  Dist        - [M x N] distance matrix
% % %
% (C) Shicai Yang, 2012
% Institute of Systems Engineering, Southeast University, Nanjing

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

switch metric
    case {0,'sqeuclidean','sqe','seu'}
        Dist = distEucSq(X,Y);
    case {'euclidean','euc','l2','L2'}
        Dist = sqrt(distEucSq(X,Y));
    case {'L0','hamming','ham','l0'}
        Dist = distL0(X,Y);
    case {'L1','cityblock','cit','l1'}
        Dist = distL1(X,Y);
    case {'Linf','linf','chebychev','che'}
        Dist = distLinf(X,Y);
    case {'cosine','cos'}
        Dist = distCosine(X,Y);
    case {'emd','EMD'}
        Dist = distEmd(X,Y);
    case {'chisq','cs','chi2'}
        Dist = distChiSq(X,Y);
    case {'hellingers','hell','hik'}
        Dist = distHell(X,Y);
    otherwise
        error(['pdist2 - unknown metric: ' metric]);
end
% Dist = max(0,Dist);
end

function D = distL0( X, Y )
    D=pdist2(X',Y','hamming');
end

function D = distL1( X, Y )
%     m = size(X,2);  n = size(Y,2);
%     mOnes = ones(1,m); D = zeros(m,n);
%     for i=1:n
%         yi = Y(:,i);  yi = yi(:,mOnes);
%         D(:,i) = sum( abs( X-yi));
%     end
    D=pdist2(X',Y','cityblock');
end

function D = distLinf( X, Y )
%     m = size(X,2);  n = size(Y,2);
%     mOnes = ones(1,m); D = zeros(m,n);
%     for i=1:n
%         yi = Y(:,i);  yi = yi(:,mOnes);
%         D(:,i) = max( abs( X-yi));
%     end
    D=pdist2(X',Y','chebychev');
end

function D = distCosine( X, Y )
    p=size(X,1);
    XX = sqrt(sum(X.*X))+eps; X = X ./ XX(ones(1,p),:);
    YY = sqrt(sum(Y.*Y))+eps; Y = Y ./ YY(ones(1,p),:);
    D = 1 - X'*Y;
end

function D = distEmd( X, Y )
    Xcdf = cumsum(X,1);
    Ycdf = cumsum(Y,1);
%     m = size(X,2);  n = size(Y,2);
%     mOnes = ones(1,m); D = zeros(m,n);
%     for i=1:n
%       ycdf = Ycdf(:,i);
%       ycdfRep = ycdf(:,mOnes);
%       D(:,i) = sum(abs(Xcdf - ycdfRep));
%     end
    D = zeros(size(X,2),size(Y,2),'double');
    for i=1:size(Y,2)
      D(:,i) = sum(abs(bsxfun(@minus,Xcdf,Ycdf(:,i))));
    end
end

function D = distChiSq( X, Y )
% note: supposedly it's possible to implement this without a loop!
%     m = size(X,2);  n = size(Y,2);
%     mOnes = ones(1,m); D = zeros(m,n);
%     for i=1:n
%       yi = Y(:,i);  yiRep = yi(:,mOnes);
%       s = yiRep + X;    d = yiRep - X;
%       D(:,i) = sum( d.^2 ./ (s+eps));
%     end
%     D = D/2;
    D = zeros(size(X,2),size(Y,2),'double');
    for i=1:size(Y,2)
      yi = Y(:,i);
      s=bsxfun(@plus,X,yi);
      d=bsxfun(@times,X,yi);
      D(:,i) = sum(s-4*d./(s+eps),1);
    end
    D = D/2;
end

function D = distEucSq( X, Y )
    XX = sum(X.*X);
    YY = sum(Y.*Y);
    D = bsxfun(@plus,XX',YY)-2*X'*Y;
end

function D = distHell( X, Y )
    m = size(X,2);  n = size(Y,2);
    D = zeros(m,n,'double');
    X=X.^.5;Y=Y.^.5;
    for i=1:n
      yi = Y(:,i);
      s=bsxfun(@minus,X,yi);
      D(:,i) = sum(s.^2,1);
    end
end

author: jindongwang

repo: activityrecognition


text:
function C = curvature(DEM,ctype,varargin)

%CURVATURE 8-connected neighborhood curvature of a digital elevation model
%
% Syntax
%
%     C = curvature(DEM)
%     C = curvature(DEM,type)
%     C = curvature(DEM,type,pn,pv,...)
%
% Description
%
%     curvature returns the second numerical derivative (curvature) of a
%     digital elevation model. By default, curvature returns the profile
%     curvature (profc).
%
% Input arguments
%
%     DEM    digital elevation model (GRIDobj)
%     type   'profc' (default) : profile curvature [m^(-1)]
%            'planc' : planform curvature or contour curvature [m^(-1)]
%            'tangc' : tangential curvature [m^(-1)]
%            'meanc' : mean curvature [m^(-1)]
%            'total' : total curvature [m^(-2)]
%
%     Parameter name value/pairs
%
%     'useblockproc'    true or {false}: use block processing
%                       (see function blockproc)
%     'useparallel'     true or {false}: use parallel computing toolbox
%     'blocksize'       blocksize for blockproc (default: 5000)
%
% Output arguments
%
%     C      curvature (GRIDobj)
%
% Remarks
%
%     Please note that curvature is not defined for cells with zero
%     gradient. Here, curvature is set to zero.
%
%     All formulas are according to Schmidt et al. (2003) on page 800.
%
% Example
%
%     DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
%     DEM = filter(DEM);
%     C = curvature(DEM,'planc');
%     imageschs(DEM,C,'percentclip',0.1)
%
% Reference
%
%     Schmidt, J., Evans, I.S., Brinkmann, J., 2003. Comparison of
%     polynomial models for land surface curvature calculation.
%     International Journal of Geographical Information Science 17,
%     797-814. doi:10.1080/13658810310001596058
%
% See also: GRIDobj/gradient8
%
% Author:  Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
% Date: 17. August, 2017

% check input arguments
narginchk(1,inf);
if nargin == 1
    ctype = 'profc';
else
    ctype = validatestring(ctype,{'profc','planc','tangc','meanc','total'});
end

p = inputParser;
p.FunctionName = 'GRIDobj/curvature';
addParamValue(p,'useblockproc',false,@(x) isscalar(x));
addParamValue(p,'blocksize',5000,@(x) isscalar(x));
addParamValue(p,'useparallel',false,@(x) isscalar(x));
parse(p,varargin{:});

% create a copy of the DEM instance
C = DEM;
c = class(DEM.Z);
switch c
    case 'double'
        C.Z = double.empty(0,0);
    otherwise
        C.Z = single.empty(0,0);
end

% Large matrix support. Break calculations in chunks using blockproc
% Parallisation for large grids using blockproc does in my experience with
% four cores hardly increase the speed.
if p.Results.useblockproc
    blksiz = bestblk(size(DEM.Z),p.Results.blocksize);
    cs  = C.cellsize;
    fun = @(x) curvaturesub(x,cs,ctype);
    C.Z = blockproc(DEM.Z,blksiz,fun,...
           'BorderSize',[1 1],...
           'Padmethod','symmetric',...
           'UseParallel',p.Results.useparallel);
else
    C.Z = curvaturesub(DEM.Z,C.cellsize,ctype);
end

C.name = ctype;

end
% subfunction

function curv = curvaturesub(dem,cs,ctype)

if isstruct(dem);
    dem = dem.data;
    % DEM has already been padded
    correctedges = false;
    shape = 'same';
else
    correctedges = true;
    shape = 'valid';
end

% First-order partial derivatives:
[fx,fy] = gradient(dem,cs);

if correctedges
    dem = padarray(dem,[1 1],'symmetric');
end
% Second order derivatives according to Evans method (see Olaya 2009)
%
% z1 z2 z3
% z4 z5 z6
% z7 z8 z9

% kernel for d2z/dx2
kernel = [1 -2 1; 1 -2 1; 1 -2 1]./(3*cs.^2);
fxx = conv2(dem,kernel,shape);
% kernel for d2z/dy2
kernel = kernel';
fyy = conv2(dem,kernel,shape);
% kernel for d2z/dxy
kernel = [-1 0 1; 0 0 0; 1 0 -1]./(4*cs.^2);
fxy = conv2(dem,kernel,shape);

%% Other options to calculate Second-order partial derivatives:
% r = gradient(p,cs);
% t = gradient(q',cs)';
% % Second-order mixed partial derivative:
% s = gradient(p',cs)';

switch ctype
    case 'profc'
        curv = - (fx.^2 .* fxx + 2*fx.*fy.*fxy + fy.^2.*fyy)./((fx.^2 + fy.^2).*(1 + fx.^2 + fy.^2).^(3/2));
    case 'tangc'
        curv = - (fy.^2 .* fxx + 2*fx.*fy.*fxy + fx.^2.*fyy)./((fx.^2 + fy.^2).*(1 + fx.^2 + fy.^2).^(1/2));
    case 'planc'
        curv = - (fy.^2 .* fxx + 2*fx.*fy.*fxy + fx.^2.*fyy)./((fx.^2 + fy.^2).^(3/2));
    case 'meanc'
        curv = - ((1+fy.^2).*fxx - 2.*fxy.*fx.*fy + (1+fx.^2).*fyy)./ ...
            (2.* (fx.^2+fy.^2+1).^(3/2));
%         curv = (fx.^2 .* fxx + 2*fx.*fy.*fxy + fy.^2.*fyy)./((fx.^2 + fy.^2).*(1 + fx.^2 + fy.^2)) ...
%             - ((1+fy).^2 .* fxx + 2*fx.*fy.*fxy + (1+fx).^2.*fyy)./(2.*(1 + fx.^2 + fy.^2).^(3/2));
    case 'total'
        curv = fxx.^2 + 2*fxy.^2+fyy.^2;
end

if correctedges
    dem = dem(2:end-1,2:end-1);
end
curv(isinf(curv) | isnan(curv)) = 0;
curv(isnan(dem)) = nan;
curv = reshape(curv,size(dem));

end

author: GERSL

repo: CCDC

save_path: github-repos/MATLAB/GERSL-CCDC


text:
function face_normal = stla_face_normal_compute ( node_num, face_num, ...
  node_xyz, face_node )

%*****************************************************************************80
%
%% STLA_FACE_NORMAL_COMPUTE computes normal vectors for an ASCII StereoLithography file.
%
%  Discussion:
%
%    This routine computes the normal vector to each triangular face
%    in the STLA solid.  If the nodes of each triangular face are
%    listed in counterclockwise order (as seen from outside the solid),
%    then the normal vectors will be properly outward facing.
%
%    The normal vectors will have unit Euclidean norm.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 September 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    3D Systems, Inc,
%    Stereolithography Interface Specification,
%    October 1989.
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer FACE_NUM, the number of faces.
%
%    Input, real NODE_XYZ(3,NODE_NUM), the node coordinates.
%
%    Input, integer FACE_NODE(3,FACE_NUM), the nodes making faces.
%
%    Input, integer FACE_MAX, the maximum number of faces.
%
%    Output, real FACE_NORMAL(3,FACE_NUM), the normal vector at each face.
%
  for face = 1 : face_num

    n1 = face_node(1,face);
    n2 = face_node(2,face);
    n3 = face_node(3,face);

    v1(1:3) = node_xyz(1:3,n2) - node_xyz(1:3,n1);
    v2(1:3) = node_xyz(1:3,n3) - node_xyz(1:3,n1);

    face_normal(1:3,face) = r8vec_cross_3d ( v1, v2 );

    norm = sqrt ( sum ( ( face_normal(1:3,face) ).^2 ) );

    if ( norm ~= 0.0 )
      face_normal(1:3,face) = face_normal(1:3,face) / norm;
    end

  end

  return
end


text:
%     Steffen Urban email: steffen.urban@kit.edu
%     Copyright (C) 2014  Steffen Urban
%
%     This program is free software; you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation; either version 2 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License along
%     with this program; if not, write to the Free Software Foundation, Inc.,
%     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

% 04.03.2014 by Steffen Urban
% error function for the center of distortion search

function [error] = errCenterUrban(x, calib_data)

error = 0;
xc = x(1);
yc = x(2);
% call calibration function
[RRfin, ss] = calibrate(calib_data.Xt, calib_data.Yt, calib_data.Xp_abs, ...
                                                        calib_data.Yp_abs, xc,yc, ...
                                                        calib_data.taylor_order, ...
                                                        calib_data.ima_proc);
lauf = 1;
M = [calib_data.Xt,calib_data.Yt,ones(size(calib_data.Xt))];

for i = 1:size(RRfin,3)
    % if calibration was not possible add a high value
    % to penalize the minimization away from that point
    if calib_data.RRfin(:,:,i)==0
      error= error+sum(ones(length(calib_data.Xp_abs),1)*sqrt( (calib_data.ocam_model.width/2)^2 + (calib_data.ocam_model.height/2)^2));
    else
        Mc = RRfin(:,:,i)*M';
        Xpp=calib_data.Xp_abs(:,:,i);
        Ypp=calib_data.Yp_abs(:,:,i);
        [xp1,yp1] = omni3d2pixel(ss, Mc, calib_data.ocam_model.width, calib_data.ocam_model.height);
        if (isinf(xp1) | isinf(yp1))
             error = error+sum(ones(length(calib_data.Xp_abs),1)*sqrt( (calib_data.ocam_model.width/2)^2 + (calib_data.ocam_model.height/2)^2));
        else
            xp = xp1 + xc;
            yp = yp1 + yc;
            lauf = lauf+length(Xpp);
            error = error + sum((Xpp-xp').^2) + sum((Ypp-yp').^2);
        end

    end
end
error = sqrt(error / lauf);

end

author: urbste

repo: ImprovedOcamCalib


text:
function asa066_test ( )

%*****************************************************************************80
%
%% ASA066_TEST tests ASA066.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ASA066_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the ASA066 library.\n' );

  asa066_test01 ( );
  asa066_test02 ( );
  asa066_test03 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ASA066_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function [C_e, C_d] = elec_distance_consec(dist_e, dist_d)

% elec_distance_consec - Identify distance between nearest neighbour electrodes
%
% Requires input data generated by 'elec_distance.m'
% Assumes using FUSA CNL electrode positions (Scan ECI-Cap)
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:54 $

% Licence:  GNU GPL, no implied or express warranties
% History:  05/2000, Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  fprintf('\n%s\n\n', 'Calculating consecutive electrode distances.');

  End_e = strvcat('3', '5', '7', '9', '13', '22', '30', '41', '49', '54', '59', '67', '78', '86', '97', '101', '108', '110', '115', '117', '120', '124');

  C_e = cell(0);
  C_d = [];

  for row = 1:length(dist_d)

     % Check that row in dist_d matrix is non-zero

     NonZeroColumns = find(dist_d(row,:));

     if ~isempty(NonZeroColumns)

	  % Calculate consecutive electrode distance,
	  % excluding end electrodes defined in End_e

		if( (row+1) <= length(dist_e) )

			e = str2num(char(dist_e(row,row+1)));

			if(strmatch(num2str(e(1)), End_e, 'exact'))

				% found end electrode, do nothing
			else

				C_d(end + 1) = dist_d(row,row+1);
                                C_e(end + 1) = {num2str(e(1))};
			end
		end
	end
  end
  clear r x;

author: PatternRecognition

repo: OpenBMI


text:
function ufmf2avi(infile,outfile,frms)

[readframe,nframes,movie_fid,headerinfo]= ...
  get_readframe_fcn(infile);

if nargin>2
  nframes = frms;
end

writerObj = VideoWriter(outfile);

writerObj.FrameRate = 30;
writerObj.Quality = 60;

open(writerObj);

for ndx = 1:nframes
  ii = readframe(ndx);
  if mod(ndx,1000) == 0,
      fprintf('Frame %d / %d, video %s\n',ndx,nframes,infile);
  end
   writeVideo(writerObj,ii);
end

close(writerObj);

if movie_fid>0
  fclose(movie_fid);

end

author: kristinbranson

repo: JAABA


text:
classdef Container < handle
    %uix.mixin.Container  Container mixin
    %
    %  uix.mixin.Container is a mixin class used by containers to provide
    %  various properties and template methods.

    %  Copyright 2009-2016 The MathWorks, Inc.
    %  $Revision: 1778 $ $Date: 2019-01-23 23:08:35 +0000 (Wed, 23 Jan 2019) $

    properties( Dependent, Access = public )
        Contents % contents in layout order
    end

    properties( Access = public, Dependent, AbortSet )
        Padding % space around contents, in pixels
    end

    properties( Access = protected )
        Contents_ = gobjects( [0 1] ) % backing for Contents
        Padding_ = 0 % backing for Padding
    end

    properties( Dependent, Access = protected )
        Dirty % needs redraw
    end

    properties( Access = private )
        Dirty_ = false % backing for Dirty
        FigureObserver % observer
        FigureListener % listener
        ChildObserver % observer
        ChildAddedListener % listener
        ChildRemovedListener % listener
        SizeChangedListener % listener
        ActivePositionPropertyListeners = cell( [0 1] ) % listeners
    end

    methods

        function obj = Container()
            %uix.mixin.Container  Initialize
            %
            %  c@uix.mixin.Container() initializes the container c.

            % Create observers and listeners
            figureObserver = uix.FigureObserver( obj );
            figureListener = event.listener( figureObserver, ...
                'FigureChanged', @obj.onFigureChanged );
            childObserver = uix.ChildObserver( obj );
            childAddedListener = event.listener( ...
                childObserver, 'ChildAdded', @obj.onChildAdded );
            childRemovedListener = event.listener( ...
                childObserver, 'ChildRemoved', @obj.onChildRemoved );
            sizeChangedListener = event.listener( ...
                obj, 'SizeChanged', @obj.onSizeChanged );

            % Store observers and listeners
            obj.FigureObserver = figureObserver;
            obj.FigureListener = figureListener;
            obj.ChildObserver = childObserver;
            obj.ChildAddedListener = childAddedListener;
            obj.ChildRemovedListener = childRemovedListener;
            obj.SizeChangedListener = sizeChangedListener;

            % Track usage
            obj.track()

        end % constructor

    end % structors

    methods

        function value = get.Contents( obj )

            value = obj.Contents_;

        end % get.Contents

        function set.Contents( obj, value )

            % For those who can't tell a column from a row...
            if isrow( value )
                value = transpose( value );
            end

            % Check
            [tf, indices] = ismember( value, obj.Contents_ );
            assert( isequal( size( obj.Contents_ ), size( value ) ) && ...
                numel( value ) == numel( unique( value ) ) && all( tf ), ...
                'uix:InvalidOperation', ...
                'Property ''Contents'' may only be set to a permutation of itself.' )

            % Call reorder
            obj.reorder( indices )

        end % set.Contents

        function value = get.Padding( obj )

            value = obj.Padding_;

        end % get.Padding

        function set.Padding( obj, value )

            % Check
            assert( isa( value, 'double' ) && isscalar( value ) && ...
                isreal( value ) && ~isinf( value ) && ...
                ~isnan( value ) && value >= 0, ...
                'uix:InvalidPropertyValue', ...
                'Property ''Padding'' must be a non-negative scalar.' )

            % Set
            obj.Padding_ = value;

            % Mark as dirty
            obj.Dirty = true;

        end % set.Padding

        function value = get.Dirty( obj )

            value = obj.Dirty_;

        end % get.Dirty

        function set.Dirty( obj, value )

            if value
                if obj.isDrawable() % drawable
                    obj.redraw() % redraw now
                else % not drawable
                    obj.Dirty_ = true; % flag for future redraw
                end
            end

        end % set.Dirty

    end % accessors

    methods( Access = private, Sealed )

        function onFigureChanged( obj, ~, eventData )
            %onFigureChanged  Event handler

            % Call template method
            obj.reparent( eventData.OldFigure, eventData.NewFigure )

            % Redraw if possible and if dirty
            if obj.Dirty_ && obj.isDrawable()
                obj.redraw()
                obj.Dirty_ = false;
            end

        end % onFigureChanged

        function onChildAdded( obj, ~, eventData )
            %onChildAdded  Event handler

            % Do nothing if add is internal tree surgery
            if isTreeSurgery( eventData.Child ), return, end

            % Call template method
            obj.addChild( eventData.Child )

        end % onChildAdded

        function onChildRemoved( obj, ~, eventData )
            %onChildRemoved  Event handler

            % Do nothing if container is being deleted
            if strcmp( obj.BeingDeleted, 'on' ), return, end

            % Do nothing if remove is internal tree surgery
            if isTreeSurgery( eventData.Child ), return, end

            % Call template method
            obj.removeChild( eventData.Child )

        end % onChildRemoved

        function onSizeChanged( obj, ~, ~ )
            %onSizeChanged  Event handler

            % Mark as dirty
            obj.Dirty = true;

        end % onSizeChanged

        function onActivePositionPropertyChanged( obj, ~, ~ )
            %onActivePositionPropertyChanged  Event handler

            % Mark as dirty
            obj.Dirty = true;

        end % onActivePositionPropertyChanged

    end % event handlers

    methods( Abstract, Access = protected )

        redraw( obj )

    end % abstract template methods

    methods( Access = protected )

        function addChild( obj, child )
            %addChild  Add child
            %
            %  c.addChild(d) adds the child d to the container c.

            % Add to contents
            obj.Contents_(end+1,:) = child;

            % Add listeners
            if isa( child, 'matlab.graphics.axis.Axes' )
                obj.ActivePositionPropertyListeners{end+1,:} = ...
                    event.proplistener( child, ...
                    findprop( child, 'ActivePositionProperty' ), ...
                    'PostSet', @obj.onActivePositionPropertyChanged );
            else
                obj.ActivePositionPropertyListeners{end+1,:} = [];
            end

            % Mark as dirty
            obj.Dirty = true;

        end % addChild

        function removeChild( obj, child )
            %removeChild  Remove child
            %
            %  c.removeChild(d) removes the child d from the container c.

            % Remove from contents
            contents = obj.Contents_;
            tf = contents == child;
            obj.Contents_(tf,:) = [];

            % Remove listeners
            obj.ActivePositionPropertyListeners(tf,:) = [];

            % Mark as dirty
            obj.Dirty = true;

        end % removeChild

        function reparent( obj, oldFigure, newFigure ) %#ok<INUSD>
            %reparent  Reparent container
            %
            %  c.reparent(a,b) reparents the container c from the figure a
            %  to the figure b.

        end % reparent

        function reorder( obj, indices )
            %reorder  Reorder contents
            %
            %  c.reorder(i) reorders the container contents using indices
            %  i, c.Contents = c.Contents(i).

            % Reorder contents
            obj.Contents_ = obj.Contents_(indices,:);

            % Reorder listeners
            obj.ActivePositionPropertyListeners = ...
                obj.ActivePositionPropertyListeners(indices,:);

            % Mark as dirty
            obj.Dirty = true;

        end % reorder

        function tf = isDrawable( obj )
            %isDrawable  Test for drawability
            %
            %  c.isDrawable() is true if the container c is drawable, and
            %  false otherwise.  To be drawable, a container must be
            %  rooted.

            tf = ~isempty( obj.FigureObserver.Figure );

        end % isDrawable

        function track( obj )
            %track  Track usage

            persistent TRACKED % single shot
            if isempty( TRACKED )
                v = ver( 'layout' );
                try %#ok<TRYNC>
                    % Disable tracking it's too hidden in the code
                    % it does not allow users to know it is happening
                    % (even if it is anonymous)
                    % uix.tracking( 'UA-82270656-2', v(1).Version, class( obj ) )
                end
                TRACKED = true;
            end

        end % track

    end % template methods

end % classdef

function tf = isTreeSurgery( child )
%isTreeSurgery  Test for internal tree surgery
%
%  Certain internal operations perform tree surgery, removing a child
%  temporarily and adding it back under an additional node.

if ~isa( child, 'matlab.graphics.axis.Axes' ) % only certain types
    tf = false;
elseif verLessThan( 'MATLAB', '9.5' ) % only certain versions
    tf = false;
else % check stack
    s = dbstack();
    tf = any( strcmp( 'ScribeStackManager.createLayer', {s.name} ) );
end

end % isTreeSurgery

author: goGPS-Project


text:
function radius = computeRadiusFromTrueAEcc(tru, sma, ecc)
%computeRadiusFromTrueAEcc Summary of this function goes here
%   Detailed explanation goes here
    p = sma.*(1-ecc.^2);
    radius = p./(1+ecc.*cos(tru));
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function r = checkTimingReport(varargin)
% Find timing report
if isempty(varargin)
    hdl_prj = 'hdl_prj';
else
    hdl_prj = varargin{1};
end
path = [hdl_prj,'/vivado_ip_prj/vivado_prj.runs/impl_1'];
filename = 'system_wrapper_timing_summary_routed.rpt';
full = [path,'/',filename];

if ~isfile(full)
    full1 = full;
    % Try ZYNQMP naming
    path = [hdl_prj,'/vivado_ip_prj/vivado_prj.runs/impl_1'];
    filename = 'system_top_timing_summary_routed.rpt';
    full = [path,'/',filename];
    if ~isfile(full)
        error(['No timing report found under: ',full,' or ',full1]);
    end
end

fid = fopen(full);
str = textscan(fid, '%s');
str = strcat(str{1}{:});
%%
if ~contains(str,'(VIOLATED)')
    disp('Passed Timing');
    r = true;
else
    disp('Failed Timing');
    r = false;
end

author: analogdevicesinc


text:
function a = le( x, y )

%   Disciplined convex programming information for LE (<=):
%      The left-hand side of a less-than constraint must be convex. The
%      right-hand side must be concave. Of course, real constant and
%      affine expressions are both convex and concave and can be used on
%      either side as well.
%
%   Disciplined geometric programming information for LE (<=):
%      The left-hand side of a less-than constraint must be log-convex,
%      including positive constants, monomials, posynomials, generalized
%      posynomials, and products thereof. The right-hand side must be
%      log-concave---including positive constants, monomials,
%      reciprocals of log-convex expressions, and products thereof.
%
%   Note that CVX does not distinguish between strict less-than (<) and
%   less-than-or-equal (<=) constraints; they are treated identically.
%   Feasible interior-point solvers tend to return points which satisfy
%   strict inequality, but not all solvers do.

b = newcnstr( evalin( 'caller', 'cvx_problem', '[]' ), x, y, '<=' );
if nargout, a = b; end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function r = length(a)
%LENGTH       Taylor length length(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
%

  r = max(a.size);

repo: OpenMAS


text:
function bw = bwfunc(re)
mm = 1;
m =mean2(re);
s = std2(re);
maxv = max(re(:));
switch mm
    case 1
        T = m + 0.5*(maxv - m);
    case 2
        ratio = 0.6;
        T = ratio * maxv;
    case 3
        ratio = 3;
        T = m+ ratio*s;
end
bw = re> T;
end

author: daxjuanxiong

repo: infrared-small-target-detection


text:
% Algorithm for finding connected components in a graph
% Note: Valid for undirected graphs only.
%
% INPUTS: adj - adjacency matrix, nxn
% OUTPUTS: a list of the components comp{i}=[j1,j2,...jk]
%
% Other routines used: findConnCompI.m, degrees.m
% GB: last updated, September 22, 2012

function comp_mat = findConnComp(adj)

[deg,~,~]=degrees(adj);            % degrees
comp_mat={};                       % initialize components matrix

for i=1:length(deg)
    if deg(i)>0
        done=0;
        for x=1:length(comp_mat)
            if length(find(comp_mat{x}==i))>0   % i in comp_mat(x).mat
                done=1;
                break
            end
        end
        if not(done)
            comp=findConnCompI(adj,i);
            comp_mat{length(comp_mat)+1}=comp;
        end

    elseif deg(i)==0
        comp_mat{length(comp_mat)+1}=[i];
    end

end

author: aeolianine

repo: octave-networks-toolbox


text:
function s=mergestruct(s1,s2)
%
% s=mergestruct(s1,s2)
%
% merge two struct objects into one
%
% authors:Qianqian Fang (q.fang <at> neu.edu)
% date: 2012/12/22
%
% input:
%      s1,s2: a struct object, s1 and s2 can not be arrays
%
% output:
%      s: the merged struct object. fields in s1 and s2 will be combined in s.
%
% license:
%     BSD License, see LICENSE_BSD.txt files for details
%
% -- this function is part of jsonlab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
%

if(~isstruct(s1) || ~isstruct(s2))
    error('input parameters contain non-struct');
end
if(length(s1)>1 || length(s2)>1)
    error('can not merge struct arrays');
end
fn=fieldnames(s2);
s=s1;
for i=1:length(fn)
    s=setfield(s,fn{i},getfield(s2,fn{i}));
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [pgx,pgy,area,perim] = minboundparallelogram(x,y,metric)
% minboundparallelogram: Compute the minimal bounding parallelogram of points in the plane
% usage: [pgx,pgy,area,perimeter] = minboundparallelogram(x,y,metric)
%
% While minboundrect allows you to choose between a minimal area or
% perimeter, it turns out that the minimal area paralellogram is not
% always unique. However, there will exist a unique parallelogram of
% minimal perimeter.
%
% arguments: (input)
%  x,y - vectors of points, describing points in the plane as
%        (x,y) pairs. x and y must be the same lengths.
%
%  metric - (OPTIONAL) - single letter character flag which
%        denotes the use of minimal area or perimeter as the
%        metric to be minimized. metric may be either 'a' or 'p',
%        capitalization is ignored. Any other contraction of 'area'
%        or 'perimeter' is also accepted.
%
%        DEFAULT: 'a'    ('area')
%
% arguments: (output)
%  pgx,pgy - 5x1 vectors of points that define the minimal
%        bounding parallelogram.
%
%  area - (scalar) area of the minimal parallelogram itself.
%
%  perimeter - (scalar) perimeter of the minimal parallelogram as found
%
%
% Example usage:
%  x = rand(50000,1);
%  y = rand(50000,1);
%  tic,[px,py,area,perim] = minboundparallelogram(x,y);toc
%
%  Elapsed time is 0.083340 seconds.
%
% [px,py]
% ans =
%       8.73336658824275e-05      3.51448852413649e-05
%           1.00002616513798       3.7428827532671e-05
%          0.999923888887252         0.999977663164396
%      -1.49425848464391e-05         0.999975379222104
%       8.73336658824275e-05      3.51448852413649e-05
%
% area
% area =
%          0.999879069698332
%
% perim
% perim =
%            3.9997581420842
%
%
% See also: minboundcircle, minboundtri, minboundrect
%
%
% Author: John D'Errico
% E-mail: woodchips@rochester.rr.com
% Release: 1.0
% Release date: 8/19/2010

% default for metric
if (nargin<3) || isempty(metric)
  metric = 'a';
elseif ~ischar(metric)
  error('metric must be a character flag if it is supplied.')
else
  % check for 'a' or 'p'
  metric = lower(metric(:)');
  ind = strmatch(metric,{'area','perimeter'});
  if isempty(ind)
    error('metric does not match either ''area'' or ''perimeter''')
  end

  % just keep the first letter.
  metric = metric(1);
end

% preprocess data
x=x(:);
y=y(:);

% not many error checks to worry about
n = length(x);
if n~=length(y)
  error('x and y must be the same sizes')
end

% start out with the convex hull of the points to
% reduce the problem dramatically. Note that any
% points in the interior of the convex hull are
% never needed, so we drop them.
if n>3
  edges = convhull(x,y);  % 'Pp' will silence the warnings

  % exclude those points inside the hull as not relevant
  % also sorts the points into their convex hull as a
  % closed polygon

  x = x(edges(1:(end-1)));
  y = y(edges(1:(end-1)));

  % probably fewer points now, unless the points are fully convex
  nedges = length(x) - 1;

elseif n > 1
  % n must be 2 or 3
  nedges = n;
  x(end+1) = x(1);
  y(end+1) = y(1);
else
  % n must be 0 or 1
  nedges = n;
end

% now we must find the bounding parallelogram of those
% that remain.

% special case small numbers of points. If we trip any
% of these cases, then we are done, so return.
switch nedges
  case 0
    % empty begets empty
    pgx = [];
    pgy = [];
    area = [];
    perimeter = [];
    return
  case 1
    % with one point, the rect is simple.
    pgx = repmat(x,1,5);
    pgy = repmat(y,1,5);
    area = 0;
    perimeter = 0;
    return
  case 2
    % only two points. also simple.
    pgx = x([1 2 2 1 1]);
    pgy = y([1 2 2 1 1]);
    area = 0;
    perimeter = 2*sqrt(diff(x).^2 + diff(y).^2);
    return
end
% 3 or more points.

% will need a 2x2 rotation matrix through an angle theta
Rmat = @(theta) [cos(theta) sin(theta);-sin(theta) cos(theta)];

% get the angle of each edge of the hull polygon.
nx = length(x);
ind1 = 1:nx;
ind2 = circshift(ind1,[0 -1]);
edgeangles = atan2(y(ind2) - y(ind1),x(ind2) - x(ind1));
% move the angle to be in [0,pi).
edgeangles = mod(edgeangles,pi);

% now just check each edge of the convex hull
nang = length(edgeangles);
area = inf;
perim = inf;
met = inf;
xy = [x,y];
for i = 1:nang
  % line that defines the base of the parallelogram
  p1 = xy(i,:);

  % we will rotate and translate the points so this
  % edge lies on the x axis. The rotation matrix is...
  rot = Rmat(-edgeangles(i));
  xyr = (xy - repmat(p1,nx,1))*rot;

  % the height of the parallelogram is
  if max(xyr(:,2)) >= max(-xyr(:,2))
    pgheight = max(xyr(:,2));
  else
    pgheight = min(xyr(:,2));
  end

  % rotate everything, but keep the angles in the interval [0,pi)
  anglesr = mod(edgeangles - edgeangles(i),pi);

  % what are the smallest and largest possible angles of the
  % remaining edges? This will define the limits of where the
  % secondary sides of the bounding parallelogram may lie.
  anglesr(i) = [];
  anglesr(anglesr == 0) = [];
  anglebounds = [min(anglesr), max(anglesr)];

  % use fminbnd to search over the family of parallelograms with
  % the given base.
  sideang = fminbnd(@(ang) pgramobj(ang,xyr,pgheight),anglebounds(1),anglebounds(2));

  % get the final parallelogram
  [pmin,amin,pgxy] = pgramobj(sideang,xyr,pgheight);

  if metric=='a'
    M_i = amin;
  else
    M_i = pmin;
  end

  % A new metric value for the current interval.
  % Is it better than the previous best?
  if M_i < met
    % keep this one
    met = M_i;
    area = amin;
    perim = pmin;

    % recover the parallelogram in the original coordinate system
    pgxy = pgxy*rot' + repmat(p1,5,1);

    pgx = pgxy(:,1);
    pgy = pgxy(:,2);

%    plot(x,y,'kx',pgx,pgy,'-ro')
%    title(num2str([perim,area]))
%    axis equal
%    pause

  end
end

% all done
% plot(x,y,'kx',pgx,pgy,'-ro')
% title(num2str([perim,area]))
% axis equal

% =================================================================
function [perim,area,pgxy] = pgramobj(ang,xyr,pgheight)
%  For a given angle (in radians) computes the bounding
%  parallelogram where the base is assumed to lie on the
%  x axis. Thus the upper face must be parallel to the x axis.

% ang will always be in the open interval (0,pi)

% The vector that points along the sides will be
sidevec = [cos(ang),sin(ang)];

% for this value of ang, what is the normal vector to
% that side of the parallelogram? It always points in
% the positive x direction, by the way it is constructed.
nvec = [sin(ang);-cos(ang)];

% What is the left most point that a side will hit against?
% The rightmost? get that information from a dot product of
% the points with the normal vector of the sides.
dp = xyr*nvec;
leftmost = min(dp);
rightmost = max(dp);

% a point on the line that contains the left side of the
% parallelogram is
leftpoint = leftmost*nvec.';

% and the right hand side line passes through this point.
rightpoint = rightmost*nvec.';

% we can define each line by the parametric equations
%  Pleft(t) = leftpoint + t*sidevec
%  Pright(t) = rightpoint + s*sidevec
%
% now find the values of s and t, such that those lines
% intersect the x axis (thus y == 0)
t0 = -leftpoint(2)./sidevec(2);
s0 = -rightpoint(2)./sidevec(2);

% the actual x axis intercepts are at
x1 = leftpoint(1) + t0*sidevec(1);
x2 = rightpoint(1) + s0*sidevec(1);

% the length of a side of the parallelogram is given by the
% definition of the trig function from triangle sides.
sidelength = pgheight./cos(pi/2 - ang);
baselength = abs(x2 - x1);

% so the perimeter of the parallellogram is given by
perim = 2*(abs(sidelength) + baselength);

% don't bother to do these unless necessary
if nargout > 1
  % the area is also trivial to compute, as
  area = abs(pgheight)*baselength;

  % and the vertices of the parallelogram are
  pgxy = [[x1,0]; [x2,0]; ...
    [x2,0] + sidelength*sidevec ; ...
    [x1,0] + sidelength*sidevec ; ...
    [x1,0]];

end

author: Sable

repo: mcbench-benchmarks


text:
gammaln(0.1) - 2.2527126517342059598697
gammaln(0.6) - .39823385806923489961685
gammaln(0.7) - .26086724653166651438573
gammaln(1.0)
gammaln(2.0)
gammaln(3.4) - 1.0923280598027415674947
gammaln(4.0) - 1.791759469228055000812477
gammaln(8.0) - 8.525161361065414300165531
gammaln(64.0) - 201.00931639928152667928
gammaln(256.0) - 1161.71210111840065079
if gammaln(0) ~= Inf
  error('gammaln(0) should be Inf');
end
if ~isnan(gammaln(-1))
  error('gammaln(-1) should be NaN');
end
if gammaln(Inf) ~= Inf
  error('gammaln(Inf) should be Inf');
end
% should be NaN?
gammaln(-Inf)
if ~isnan(gammaln(NaN))
  error('gammaln(NaN) should be NaN');
end

author: Cloud-CV

repo: object-proposals


text:
function fem = genFEM3D(mesh,femtype,bcind,option)

%% Usage: Form FEM global degrees of freedom on 3D tetrahedron mesh
%
% INPUTS:
% mesh --- a struct data contains rich mesh information.
% bcind --- [bc1,bc2,bc3,bc4,bc5,bc6] denotes the boundary condition of
%            bc1: x = xmin  bc2: x = xmax
%            bc3: y = ymin  bc4: y = ymax
%            bc5: z = zmin  bc6: z = zmax
%            if bc = 1: Dirichlet boundary condition
%               bc = 2: Neumann boundary condition
%               bc = 3: Robin boundary condition
% femtype --- the type of finite element methods
%             Conforming: P1, P2
%
% OUTPUTS:
% fem --- a struct data contains the following fields:
%         fem.p: (x,y,z) coordinate of each vertex w.r.t a global DoF
%         fem.t: indices of global DoF in each element
%         fem.bc: indices of global DoF on the boundary
%         fem.mapper: a vector to extract unknowns from global DoF.
%         fem.type: type of finite element methods
%         fem.ldof: number of local DoF on each element
%         fem.ng: number of Gaussian quadrature points on each element

% Last Modified: 08/07/2020 by Xu Zhang

%% 0. Inputs
if nargin == 2
    bcind = [1,1,1,1,1,1];  option.ng = 4;
end
if nargin == 3
    option.ng = 4;
end
if ~isfield(option,'ng')
    option.ng = 4;
end
ng = option.ng;

%% 1. Form fem p, t, and basis functions
if strcmp(femtype,'P1')
    fem = genP1FEM3D(mesh); option.ngauss = 4;
elseif strcmp(femtype,'P2')
    fem = genP2FEM3D(mesh); option.ngauss = 11;
end

%% 2. Form Gaussian Quadrature
p = mesh.p; t = mesh.t;
X1 = p(t(:,1),:); X2 = p(t(:,2),:); X3 = p(t(:,3),:); X4 = p(t(:,4),:);
gw = gaussWtetra(ng);
[gx,gy,gz] = gaussPtetra(X1,X2,X3,X4,ng);
A = tetraArea(X1,X2,X3,X4);
fem.gw = gw; fem.gx = gx; fem.gy = gy; fem.gz = gz; fem.area = A;

%% 3. Form bc mapper
[bc,mapper] = boundaryNode3D(mesh,fem.p,bcind);
fem.bc = bc; fem.mapper = mapper; fem.ng = option.ng; fem.bcind = bcind;

repo: ifem


text:
%% Machine Learning Online Class - Exercise 3 | Part 2: Neural Networks

%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  linear exercise. You will need to complete the following functions
%  in this exericse:
%
%     lrCostFunction.m (logistic regression cost function)
%     oneVsAll.m
%     predictOneVsAll.m
%     predict.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% Initialization
clear ; close all; clc

%% Setup the parameters you will use for this exercise
input_layer_size  = 400;  % 20x20 Input Images of Digits
hidden_layer_size = 25;   % 25 hidden units
num_labels = 10;          % 10 labels, from 1 to 10
                          % (note that we have mapped "0" to label 10)

%% =========== Part 1: Loading and Visualizing Data =============
%  We start the exercise by first loading and visualizing the dataset.
%  You will be working with a dataset that contains handwritten digits.
%

% Load Training Data
fprintf('Loading and Visualizing Data ...\n')

load('ex3data1.mat');
m = size(X, 1);

% Randomly select 100 data points to display
sel = randperm(size(X, 1));
sel = sel(1:100);

displayData(X(sel, :));

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ================ Part 2: Loading Pameters ================
% In this part of the exercise, we load some pre-initialized
% neural network parameters.

fprintf('\nLoading Saved Neural Network Parameters ...\n')

% Load the weights into variables Theta1 and Theta2
load('ex3weights.mat');

%% ================= Part 3: Implement Predict =================
%  After training the neural network, we would like to use it to predict
%  the labels. You will now implement the "predict" function to use the
%  neural network to predict the labels of the training set. This lets
%  you compute the training set accuracy.

pred = predict(Theta1, Theta2, X);

fprintf('\nTraining Set Accuracy: %f\n', mean(double(pred == y)) * 100);

fprintf('Program paused. Press enter to continue.\n');
pause;

%  To give you an idea of the network's output, you can also run
%  through the examples one at the a time to see what it is predicting.

%  Randomly permute examples
rp = randperm(m);

for i = 1:m
    % Display
    fprintf('\nDisplaying Example Image\n');
    displayData(X(rp(i), :));

    pred = predict(Theta1, Theta2, X(rp(i),:));
    fprintf('\nNeural Network Prediction: %d (digit %d)\n', pred, mod(pred, 10));

    % Pause with quit option
    s = input('Paused - press enter to continue, q to exit:','s');
    if s == 'q'
      break
    end
end

author: Ayatans

repo: Machine-Learning-homework


text:
name='bimodal.png'
raw=double(imread(name));
rawgray=raw(:,:,1);
smoothgray=imgaussfilt(rawgray,70);
figure()
surf(rawgray,'edgecolor','none')
figure()
surf(smoothgray,'edgecolor','none')
imwrite(uint16(smoothgray*255),['green' name '.tiff']);
testIN=imread(['green' name]);
figure()
surf(testIN,'edgecolor','none')

author: BrianHaidet

repo: AlphaPhoenix


text:
function outputImage = imdilate(this, structureElement, nD)
% Dilates image clusters slice-wise; mimicks imdilate in matlab functionality
%
%   Y = MrImage()
%   dilatedY = Y.imdilate(structureElement, nD)
%
% This is a method of class MrImage.
%
% IN
%   structureElement
%           morphological structuring element for dilation, e.g.
%           strel('disk', 2) for a disk of radius 2
%           default: strel('disk', 2)
%   nD      dimensionality to perform operation
%           '2d' = slicewise application, separate 2d images
%           '3d' = as volume
%
% OUT
%   outputImage
%           MrImage where data matrix is inflated
%
% EXAMPLE
%   Y = MrImage();
%   dilatedY = Y.imdilate()
%   dilatedY = Y.imdilate(strel('disk', 5))
%
%
%   See also MrImage imdilate MrImage.imerode perform_unary_operation

% Author:   Saskia Klein & Lars Kasper
% Created:  2014-08-04
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public Licence (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

% update geometry-header with flipping left-right
if nargin < 2
    structureElement = strel('disk', 2);
end

if nargin < 3
    nD = '2d';
end

if isreal(this)
    outputImage = this.perform_unary_operation(...
        @(x) imdilate(x, structureElement), nD);
else
    outputImage = this.abs.perform_unary_operation(...
        @(x) imdilate(x, structureElement), nD);
end

author: translationalneuromodeling

repo: tapas


text:
function [S] = readSEL(filename)
  % READSEL read selection from .sel file (as used by Mario Botsch and Olga
  % Sorkine in the accompanying data to their survey "On linear variational
  % surface deformation methods"
  % http://igl.ethz.ch/projects/deformation-survey/
  %
  % Input:
  %   filename  path to .sel file
  % Output:
  %   S #V list of ids:
  %     0  fixed boundary (silver, completely fixed handle-region)
  %     1  interior (blue, "solved-for", "free" region)
  %     2  handle boundary (gold, re-arrangable handle-region)
  %
  % Example:
  %  S = readSEL(filename);
  %  % Swap 0 (fixed boundary) with 1 (interior) so that selection mask
  %  % functions as handle ids (with 0 meaning interior)
  %  R = (S<1).*1 + (S==1).*0 + (S>1).*S;
  %

  fp = fopen(filename,'r');
  % read initial comments
  C = textscan(fp, '#%[^\n]');
  % remaining lines contain one integer per line
  S = fscanf(fp,'%d');
  % clean up file
  fclose(fp);
end

author: alecjacobson

repo: gptoolbox