text:
function y = vl_nnloss(x,c,varargin)
%VL_NNLOSS CNN categorical or attribute loss.
%   Y = VL_NNLOSS(X, C) computes the loss incurred by the prediction
%   scores X given the categorical labels C.
%
%   The prediction scores X are organised as a field of prediction
%   vectors, represented by a H x W x D x N array. The first two
%   dimensions, H and W, are spatial and correspond to the height and
%   width of the field; the third dimension D is the number of
%   categories or classes; finally, the dimension N is the number of
%   data items (images) packed in the array.
%
%   While often one has H = W = 1, the case W, H > 1 is useful in
%   dense labelling problems such as image segmentation. In the latter
%   case, the loss is summed across pixels (contributions can be
%   weighed using the `InstanceWeights` option described below).
%
%   The array C contains the categorical labels. In the simplest case,
%   C is an array of integers in the range [1, D] with N elements
%   specifying one label for each of the N images. If H, W > 1, the
%   same label is implicitly applied to all spatial locations.
%
%   In the second form, C has dimension H x W x 1 x N and specifies a
%   categorical label for each spatial location.
%
%   In the third form, C has dimension H x W x D x N and specifies
%   attributes rather than categories. Here elements in C are either
%   +1 or -1 and C, where +1 denotes that an attribute is present and
%   -1 that it is not. The key difference is that multiple attributes
%   can be active at the same time, while categories are mutually
%   exclusive. By default, the loss is *summed* across attributes
%   (unless otherwise specified using the `InstanceWeights` option
%   described below).
%
%   DZDX = VL_NNLOSS(X, C, DZDY) computes the derivative of the block
%   projected onto the output derivative DZDY. DZDX and DZDY have the
%   same dimensions as X and Y respectively.
%
%   VL_NNLOSS() supports several loss functions, which can be selected
%   by using the option `type` described below. When each scalar c in
%   C is interpreted as a categorical label (first two forms above),
%   the following losses can be used:
%
%   Classification error:: `classerror`
%     L(X,c) = (argmax_q X(q) ~= c). Note that the classification
%     error derivative is flat; therefore this loss is useful for
%     assessment, but not for training a model.
%
%   Top-K classification error:: `topkerror`
%     L(X,c) = (rank X(c) in X <= K). The top rank is the one with
%     highest score. For K=1, this is the same as the
%     classification error. K is controlled by the `topK` option.
%
%   Log loss:: `log`
%     L(X,c) = - log(X(c)). This function assumes that X(c) is the
%     predicted probability of class c (hence the vector X must be non
%     negative and sum to one).
%
%   Softmax log loss (multinomial logistic loss):: `softmaxlog`
%     L(X,c) = - log(P(c)) where P(c) = exp(X(c)) / sum_q exp(X(q)).
%     This is the same as the `log` loss, but renormalizes the
%     predictions using the softmax function.
%
%   Multiclass hinge loss:: `mhinge`
%     L(X,c) = max{0, 1 - X(c)}. This function assumes that X(c) is
%     the score margin for class c against the other classes.  See
%     also the `mmhinge` loss below.
%
%   Multiclass structured hinge loss:: `mshinge`
%     L(X,c) = max{0, 1 - M(c)} where M(c) = X(c) - max_{q ~= c}
%     X(q). This is the same as the `mhinge` loss, but computes the
%     margin between the prediction scores first. This is also known
%     the Crammer-Singer loss, an example of a structured prediction
%     loss.
%
%   When C is a vector of binary attribures c in (+1,-1), each scalar
%   prediction score x is interpreted as voting for the presence or
%   absence of a particular attribute. The following losses can be
%   used:
%
%   Binary classification error:: `binaryerror`
%     L(x,c) = (sign(x - t) ~= c). t is a threshold that can be
%     specified using the `threshold` option and defaults to zero. If
%     x is a probability, it should be set to 0.5.
%
%   Binary log loss:: `binarylog`
%     L(x,c) = - log(c(x-0.5) + 0.5). x is assumed to be the
%     probability that the attribute is active (c=+1). Hence x must be
%     a number in the range [0,1]. This is the binary version of the
%     `log` loss.
%
%   Logistic log loss:: `logistic`
%     L(x,c) = log(1 + exp(- cx)). This is the same as the `binarylog`
%     loss, but implicitly normalizes the score x into a probability
%     using the logistic (sigmoid) function: p = sigmoid(x) = 1 / (1 +
%     exp(-x)). This is also equivalent to `softmaxlog` loss where
%     class c=+1 is assigned score x and class c=-1 is assigned score
%     0.
%
%   Hinge loss:: `hinge`
%     L(x,c) = max{0, 1 - cx}. This is the standard hinge loss for
%     binary classification. This is equivalent to the `mshinge` loss
%     if class c=+1 is assigned score x and class c=-1 is assigned
%     score 0.
%
%   VL_NNLOSS(...,'OPT', VALUE, ...) supports these additionals
%   options:
%
%   InstanceWeights:: []
%     Allows to weight the loss as L'(x,c) = WGT L(x,c), where WGT is
%     a per-instance weight extracted from the array
%     `InstanceWeights`. For categorical losses, this is either a H x
%     W x 1 or a H x W x 1 x N array. For attribute losses, this is
%     either a H x W x D or a H x W x D x N array.
%
%   TopK:: 5
%     Top-K value for the top-K error. Note that K should not
%     exceed the number of labels.
%
%   See also: VL_NNSOFTMAX().

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

if ~isempty(varargin) && ~ischar(varargin{1})  % passed in dzdy
  dzdy = varargin{1} ;
  varargin(1) = [] ;
else
  dzdy = [] ;
end

opts.instanceWeights = [] ;
opts.classWeights = [] ;
opts.threshold = 0 ;
opts.loss = 'softmaxlog' ;
opts.topK = 5 ;
opts = vl_argparse(opts, varargin, 'nonrecursive') ;

inputSize = [size(x,1) size(x,2) size(x,3) size(x,4)] ;

% Form 1: C has one label per image. In this case, get C in form 2 or
% form 3.
c = gather(c) ;
if numel(c) == inputSize(4)
  c = reshape(c, [1 1 1 inputSize(4)]) ;
  c = repmat(c, inputSize(1:2)) ;
end

hasIgnoreLabel = any(c(:) == 0);

% --------------------------------------------------------------------
% Spatial weighting
% --------------------------------------------------------------------

% work around a bug in MATLAB, where native cast() would slow
% progressively
if isa(x, 'gpuArray')
  switch classUnderlying(x) ;
    case 'single', cast = @(z) single(z) ;
    case 'double', cast = @(z) double(z) ;
  end
else
  switch class(x)
    case 'single', cast = @(z) single(z) ;
    case 'double', cast = @(z) double(z) ;
  end
end

labelSize = [size(c,1) size(c,2) size(c,3) size(c,4)] ;
assert(isequal(labelSize(1:2), inputSize(1:2))) ;
assert(labelSize(4) == inputSize(4)) ;
instanceWeights = [] ;
switch lower(opts.loss)
  case {'classerror', 'topkerror', 'log', 'softmaxlog', 'mhinge', 'mshinge'}
    % there must be one categorical label per prediction vector
    assert(labelSize(3) == 1) ;

    if hasIgnoreLabel
      % null labels denote instances that should be skipped
      instanceWeights = cast(c(:,:,1,:) ~= 0) ;
    end

  case {'binaryerror', 'binarylog', 'logistic', 'hinge'}

    % there must be one categorical label per prediction scalar
    assert(labelSize(3) == inputSize(3)) ;

    if hasIgnoreLabel
      % null labels denote instances that should be skipped
      instanceWeights = cast(c ~= 0) ;
    end

  otherwise
    error('Unknown loss ''%s''.', opts.loss) ;
end

if ~isempty(opts.instanceWeights)
  % important: this code needs to broadcast opts.instanceWeights to
  % an array of the same size as c
  if isempty(instanceWeights)
    instanceWeights = bsxfun(@times, onesLike(c), opts.instanceWeights) ;
  else
    instanceWeights = bsxfun(@times, instanceWeights, opts.instanceWeights);
  end
end

% --------------------------------------------------------------------
% Do the work
% --------------------------------------------------------------------

switch lower(opts.loss)
  case {'log', 'softmaxlog', 'mhinge', 'mshinge'}
    % from category labels to indexes
    numPixelsPerImage = prod(inputSize(1:2)) ;
    numPixels = numPixelsPerImage * inputSize(4) ;
    imageVolume = numPixelsPerImage * inputSize(3) ;

    n = reshape(0:numPixels-1,labelSize) ;
    offset = 1 + mod(n, numPixelsPerImage) + ...
             imageVolume * fix(n / numPixelsPerImage) ;
    ci = offset + numPixelsPerImage * max(c - 1,0) ;
end

if nargin <= 2 || isempty(dzdy)
  switch lower(opts.loss)
    case 'classerror'
      [~,chat] = max(x,[],3) ;
      t = cast(c ~= chat) ;
    case 'topkerror'
      [~,predictions] = sort(x,3,'descend') ;
      t = 1 - sum(bsxfun(@eq, c, predictions(:,:,1:opts.topK,:)), 3) ;
    case 'log'
      t = - log(x(ci)) ;
    case 'softmaxlog'
      Xmax = max(x,[],3) ;
      ex = exp(bsxfun(@minus, x, Xmax)) ;
      t = Xmax + log(sum(ex,3)) - x(ci) ;
    case 'mhinge'
      t = max(0, 1 - x(ci)) ;
    case 'mshinge'
      Q = x ;
      Q(ci) = -inf ;
      t = max(0, 1 - x(ci) + max(Q,[],3)) ;
    case 'binaryerror'
      t = cast(sign(x - opts.threshold) ~= c) ;
    case 'binarylog'
      t = -log(c.*(x-0.5) + 0.5) ;
    case 'logistic'
      %t = log(1 + exp(-c.*X)) ;
      a = -c.*x ;
      b = max(0, a) ;
      t = b + log(exp(-b) + exp(a-b)) ;
    case 'hinge'
      t = max(0, 1 - c.*x) ;
  end
  if ~isempty(instanceWeights)
    y = instanceWeights(:)' * t(:) ;
  else
    y = sum(t(:));
  end
else
  if ~isempty(instanceWeights)
    dzdy = dzdy * instanceWeights ;
  end
  switch lower(opts.loss)
    case {'classerror', 'topkerror'}
      y = zerosLike(x) ;
    case 'log'
      y = zerosLike(x) ;
      y(ci) = - dzdy ./ max(x(ci), 1e-8) ;
    case 'softmaxlog'
      Xmax = max(x,[],3) ;
      ex = exp(bsxfun(@minus, x, Xmax)) ;
      y = bsxfun(@rdivide, ex, sum(ex,3)) ;
      y(ci) = y(ci) - 1 ;
      y = bsxfun(@times, dzdy, y) ;
    case 'mhinge'
      y = zerosLike(x) ;
      y(ci) = - dzdy .* (x(ci) < 1) ;
    case 'mshinge'
      Q = x ;
      Q(ci) = -inf ;
      [~, q] = max(Q,[],3) ;
      qi = offset + numPixelsPerImage * (q - 1) ;
      W = dzdy .* (x(ci) - x(qi) < 1) ;
      y = zerosLike(x) ;
      y(ci) = - W ;
      y(qi) = + W ;
    case 'binaryerror'
      y = zerosLike(x) ;
    case 'binarylog'
      y = - dzdy ./ (x + (c-1)*0.5) ;
    case 'logistic'
      % t = exp(-Y.*X) / (1 + exp(-Y.*X)) .* (-Y)
      % t = 1 / (1 + exp(Y.*X)) .* (-Y)
      y = - dzdy .* c ./ (1 + exp(c.*x)) ;
    case 'hinge'
      y = - dzdy .* c .* (c.*x < 1) ;
  end
end

% --------------------------------------------------------------------
function y = zerosLike(x)
% --------------------------------------------------------------------
if isa(x,'gpuArray')
  y = gpuArray.zeros(size(x),classUnderlying(x)) ;
else
  y = zeros(size(x),'like',x) ;
end

% --------------------------------------------------------------------
function y = onesLike(x)
% --------------------------------------------------------------------
if isa(x,'gpuArray')
  y = gpuArray.ones(size(x),classUnderlying(x)) ;
else
  y = ones(size(x),'like',x) ;
end

repo: LapSRN


text:
% reload paths etc.
% NOTE!
% The mex lock has been removed from caffe_.cpp!
%munlock('caffe_') ;
clear all ; %#ok
caffePath = '/users/albanie/coding/libs/caffes/se-caffe/matlab' ;
%munlock('caffe_') ;
%clear caffe_.mexa64 ;
clear caffe_ ;
rehash path ;

% hack to avoid cuda errors
a = gpuArray(1); clear a ; %#ok

% refresh/set up ssd-caffe
addpath(caffePath) ;

caffeRoot = '/users/albanie/coding/libs/caffes/se-caffe' ;

% set paths and load caffe ssd-model
modelDir = '/users/albanie/data/models/caffe/SE-nets/proto' ;
dataDir = '/users/albanie/data/models/caffe/SE-nets/weights' ;
model = fullfile(modelDir, 'SE-ResNet-50-debug.prototxt') ;
weights = fullfile(dataDir, 'SE-ResNet-50.caffemodel') ;
caffe.set_mode_cpu() ;

% to use the relative paths defined in the prototxt, we change into the
% caffe root to load the network
cd(caffeRoot) ;

% load model
caffeNet = caffe.Net(model, weights, 'test') ;

% load mcn model for comparison

% set path to trunk model (VGG-VD-16)
opts.trunkModelPath = fullfile(vl_rootnn, 'data', 'models-import', ...
                                'SE-ResNet-50-mcn.mat') ;
opts.imageSize = [224 224] ;
net = load(opts.trunkModelPath) ; net = dagnn.DagNN.loadobj(net) ;

repo: MFT


text:
function [data] = blc(data, interval, optional)

% BLC does a baseline correction using the prestimulus interval of the data
%
%   [data] = baseline(data, interval);
%   [data] = baseline(data, begin, end);
%
% If no begin and end are specified, the whole timeinterval is
% used for baseline correction.

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

% determine the dimension of the data
if length(size(data))==3
  % multiple epochs
  dim=3;
  nepoch = size(data,1);
  nchan  = size(data,2);
  ntime  = size(data,3);
else
  % single epoch with multiple channels
  dim=2;
  nchan  = size(data,1);
  ntime  = size(data,2);
end

% determine the interval to use for baseline correction
if nargin==1
  % default is to use the whole timeinterval for baseline correction
  interval = 1:ntime;
elseif nargin==2
  % use the interval specified by the user
elseif nargin==3
  % create the interval from the specified begin and end point
  interval = interval:optional;
end

if dim==3
  % the data contains multiple epochs
  for epoch=1:nepoch
    baseline = mean(data(epoch,:,interval), 3);
    for chan=1:nchan
      data(epoch,chan,:) = data(epoch,chan,:) - baseline(chan);
    end
  end
else
  % the data contains a single epoch
  baseline = mean(data(:,interval), 2);
  for chan=1:nchan
    data(chan,:) = data(chan,:) - baseline(chan);
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef testEnergyEquivalenceVoigtAndTensorNotationForIsoTensor  ...
         < testEnergyEquivalenceVoigtAndTensorNotation

    methods (Access = public)

        function obj = testEnergyEquivalenceVoigtAndTensorNotationForIsoTensor()
            obj@testEnergyEquivalenceVoigtAndTensorNotation();
        end

    end

    methods
        function generateFourthOrderTensor(obj)
             obj.Ch = IsotropicConstitutiveTensor(1,1/3);
        end
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [Xc3,tri3,xc3,xp3,dc3,xc_texture,nc3,conf_nc3,Nn3] = Meshing(Xc2,xc2,xp2,Thresh_connect,N_smoothing,om,T,N_x,N_y,fc,cc,kc,alpha_c,fp,cp,kp,alpha_p),

% scaled connection threshold

T_connect = Thresh_connect; 	% scaled threshold

fprintf(1,'Organizing the data...\n');

xp_frac = rem(xp2,1);
xc_frac = rem(xc2(1,:),1);

if std(xp_frac) > std(xc_frac),
    disp('Dense depth map in the image coordinates');
    temporal = 1;
    spatial = 0;
else
    disp('Dense depth map in the cross image and projector frame');
    temporal = 0;
    spatial = 1;
end;

if spatial,

    % Something to fix the organization:
    xpmin = min(xp2);
    xp_ind = round(unique(xp2 - xpmin));
    step_xp = min(diff(xp_ind)); %xp_ind(2) - xp_ind(1);
    xp4 = (xp2 - xpmin)/step_xp + xpmin;

    xpmin = min(xp4);
    xpmax = max(xp4);

    xcmin = min(xc2(2,:));
    xcmax = max(xc2(2,:));

else

    % Something to fix the organization:

    xp4 = xc2(1,:);

    xpmin = min(xp4);
    xpmax = max(xp4);

    xcmin = min(xc2(2,:));
    xcmax = max(xc2(2,:));

end;

Nrow = xcmax - xcmin + 1;
Ncol = xpmax - xpmin + 1;

Xmesh = zeros(Nrow,Ncol);
Ymesh = zeros(Nrow,Ncol);
Zmesh =  zeros(Nrow,Ncol);
Cmesh = zeros(Nrow,Ncol);
xcmesh = zeros(Nrow,Ncol);
ycmesh = zeros(Nrow,Ncol);

ind_col = round(xp4 - xpmin + 1);
ind_row = xc2(2,:) - xcmin + 1;
ind = ind_row + (ind_col-1)*Nrow;
ni = length(ind);

% Good format for ivview:
Xmesh(ind) = Xc2(1,:);
Ymesh(ind) = Xc2(2,:); 		% tries to have it in the right position
Zmesh(ind) = Xc2(3,:); 		% tries to have it in the right position
Cmesh(ind) = ones(1,ni);
xcmesh(ind) = xc2(1,:);
ycmesh(ind) = xc2(2,:);

% Hypothesis on the triangles:

D1 = Cmesh(2:Nrow,1:(Ncol-1)) & Cmesh(1:(Nrow-1),2:Ncol);
F1 = Cmesh(1:(Nrow-1),1:(Ncol-1)) & D1;
F2 = Cmesh(2:Nrow,2:Ncol) & D1;
C1 =  (F1 | F2);

D2 = Cmesh(1:(Nrow-1),1:(Ncol-1)) & Cmesh(2:Nrow,2:Ncol);
F3 = Cmesh(1:(Nrow-1),2:Ncol) & D2;
F4 = Cmesh(2:Nrow,1:(Ncol-1)) & D2;
C2 = (F3 | F4);

Ambi = C1 & C2; 			% ambiguous zones
Div1 = C1 & ~C2; 			% needs to check relative distance of points
Div2 = ~C1 & C2; 			% needs to check relative distance of points

% diagonal measure:

%Dm1 = abs(Zmesh(2:Nrow,1:(Ncol-1)) -  Zmesh(1:(Nrow-1),2:Ncol));

Dm1 =( Xmesh(2:Nrow,1:(Ncol-1)) -  Xmesh(1:(Nrow-1),2:Ncol)).^2 + (Ymesh(2:Nrow,1:(Ncol-1)) -  Ymesh(1:(Nrow-1),2:Ncol)).^2 + (Zmesh(2:Nrow,1:(Ncol-1)) -  Zmesh(1:(Nrow-1),2:Ncol)).^2;

%Dm2 = abs(Zmesh(1:(Nrow-1),1:(Ncol-1)) - Zmesh(2:Nrow,2:Ncol));

Dm2 = (Xmesh(1:(Nrow-1),1:(Ncol-1)) - Xmesh(2:Nrow,2:Ncol)).^2 + (Ymesh(1:(Nrow-1),1:(Ncol-1)) - Ymesh(2:Nrow,2:Ncol)).^2 + (Zmesh(1:(Nrow-1),1:(Ncol-1)) - Zmesh(2:Nrow,2:Ncol)).^2 ;

Div1n = Div1 | ((Dm1 <= Dm2)&Ambi);
Div2n = Div2 | ((Dm2 < Dm1)&Ambi);

Div11 = Div1n & F1;
Div12 = Div1n & F2;
Div21 = Div2n & F3;
Div22 = Div2n & F4;

% look at local difference:

dZ_r = abs(Zmesh(:,2:Ncol)-Zmesh(:,1:(Ncol-1)));
dZ_c = abs(Zmesh(2:Nrow,:)-Zmesh(1:(Nrow-1),:));

Div11 = Div11 & (dZ_r(1:(Nrow-1),:)<T_connect) & (dZ_c(:,1:(Ncol-1))<T_connect); % & (Dm1 < T_connect);
Div12 = Div12 & (dZ_r(2:Nrow,:)<T_connect) & (dZ_c(:,2:Ncol)<T_connect); %& (Dm1 < T_connect);

Div21 = Div21 & (dZ_r(1:(Nrow-1),:)<T_connect) & (dZ_c(:,2:Ncol)<T_connect); % & (Dm2 < T_connect);
Div22 = Div22 & (dZ_r(2:Nrow,:)<T_connect) & (dZ_c(:,1:(Ncol-1))<T_connect); % & (Dm2 < T_connect);

%%% Smoothing:

for i = 1:N_smoothing,

    fprintf(1,'Surface smoothing %d\n',i);

    % first find the neighbor points of every point:

    t = zeros(Nrow,Ncol);
    b = zeros(Nrow,Ncol);
    l = zeros(Nrow,Ncol);
    r = zeros(Nrow,Ncol);
    tr = zeros(Nrow,Ncol);
    br = zeros(Nrow,Ncol);
    tl = zeros(Nrow,Ncol);
    bl = zeros(Nrow,Ncol);

    t(2:Nrow,2:Ncol) = (Div12 | Div21);
    t(2:Nrow,1:(Ncol-1)) = t(2:Nrow,1:(Ncol-1)) | (Div11 | Div22);
    b(1:(Nrow-1),2:Ncol) = (Div12 | Div21);
    b(1:(Nrow-1),1:(Ncol-1)) = b(1:(Nrow-1),1:(Ncol-1)) | (Div11 | Div22);
    r(2:Nrow,1:(Ncol-1)) = (Div12 | Div22);
    r(1:(Nrow-1),1:(Ncol-1)) = r(1:(Nrow-1),1:(Ncol-1)) | (Div11 | Div21);
    l(2:Nrow,2:Ncol) = (Div12 | Div22);
    l(1:(Nrow-1),2:Ncol) = l(1:(Nrow-1),2:Ncol) | (Div11 | Div21);
    tr(2:Nrow,1:(Ncol-1)) = Div11 | Div12;
    br(1:(Nrow-1),1:(Ncol-1)) = Div21 | Div22;
    tl(2:Nrow,2:Ncol) = Div21 | Div22;
    bl(1:(Nrow-1),2:Ncol) = Div11 | Div12;

    Nn = t + b + l + r + tr + br + tl + bl;

    XX = Xmesh; 		     	% zeros(Nrow,2); zeros(2,Ncol+2)];
    YY = Ymesh; 				% zeros(Nrow,2); zeros(2,Ncol+2)];
    ZZ = Zmesh; 				% zeros(Nrow,2); zeros(2,Ncol+2)];

    [is,js] = find(Nn);

    indd = find((is > 1) & (is < Nrow) & (js > 1) & (js < Ncol));

    is = is(indd);
    js = js(indd);

    sm =  is + (js-1)*(Nrow);
    sm_t = is + (js-1)*(Nrow) - 1;
    sm_b = is + (js-1)*(Nrow) + 1;
    sm_r = is + (js)*(Nrow);
    sm_l = is + (js-2)*(Nrow);

    sm_tr =  is + (js)*(Nrow) - 1;
    sm_br =  is + (js)*(Nrow) + 1;
    sm_tl =  is + (js-2)*(Nrow) -1;
    sm_bl =  is + (js-2)*(Nrow) + 1;

    XX(sm) = 0.5 * XX(sm) + 0.5 * ((XX(sm_t).*t(sm)+XX(sm_b).*b(sm)+XX(sm_r).*r(sm)+XX(sm_l).*l(sm)+XX(sm_tr).*tr(sm)+XX(sm_br).*br(sm)+XX(sm_tl).*tl(sm)+XX(sm_bl).*bl(sm))./Nn(sm));

    YY(sm) = 0.5 * YY(sm) + 0.5 * ((YY(sm_t).*t(sm)+YY(sm_b).*b(sm)+YY(sm_r).*r(sm)+YY(sm_l).*l(sm)+YY(sm_tr).*tr(sm)+YY(sm_br).*br(sm)+YY(sm_tl).*tl(sm)+YY(sm_bl).*bl(sm))./Nn(sm));

    ZZ(sm) = 0.5 * ZZ(sm) + 0.5 * ((ZZ(sm_t).*t(sm)+ZZ(sm_b).*b(sm)+ZZ(sm_r).*r(sm)+ZZ(sm_l).*l(sm)+ZZ(sm_tr).*tr(sm)+ZZ(sm_br).*br(sm)+ZZ(sm_tl).*tl(sm)+ZZ(sm_bl).*bl(sm))./Nn(sm));

    Xmesh = XX(1:Nrow,1:Ncol);
    Ymesh = YY(1:Nrow,1:Ncol);
    Zmesh = ZZ(1:Nrow,1:Ncol);

    %%% reconnect after smoothing:

    % diagonal measure:

    Dm1 =( Xmesh(2:Nrow,1:(Ncol-1)) -  Xmesh(1:(Nrow-1),2:Ncol)).^2 + (Ymesh(2:Nrow,1:(Ncol-1)) -  Ymesh(1:(Nrow-1),2:Ncol)).^2 + (Zmesh(2:Nrow,1:(Ncol-1)) -  Zmesh(1:(Nrow-1),2:Ncol)).^2;

    Dm2 = (Xmesh(1:(Nrow-1),1:(Ncol-1)) - Xmesh(2:Nrow,2:Ncol)).^2 + (Ymesh(1:(Nrow-1),1:(Ncol-1)) - Ymesh(2:Nrow,2:Ncol)).^2 + (Zmesh(1:(Nrow-1),1:(Ncol-1)) - Zmesh(2:Nrow,2:Ncol)).^2 ;

    Div1n = Div1 | ((Dm1 <= Dm2)&Ambi);
    Div2n = Div2 | ((Dm2 < Dm1)&Ambi);

    Div11 = Div1n & F1;
    Div12 = Div1n & F2;
    Div21 = Div2n & F3;
    Div22 = Div2n & F4;

    % look at local difference:

    dZ_r = abs(Zmesh(:,2:Ncol)-Zmesh(:,1:(Ncol-1)));
    dZ_c = abs(Zmesh(2:Nrow,:)-Zmesh(1:(Nrow-1),:));

    Div11 = Div11 & (dZ_r(1:(Nrow-1),:)<T_connect) & (dZ_c(:,1:(Ncol-1))<T_connect); % & (Dm1 < T_connect);
    Div12 = Div12 & (dZ_r(2:Nrow,:)<T_connect) & (dZ_c(:,2:Ncol)<T_connect); %& (Dm1 < T_connect);

    Div21 = Div21 & (dZ_r(1:(Nrow-1),:)<T_connect) & (dZ_c(:,2:Ncol)<T_connect); % & (Dm2 < T_connect);
    Div22 = Div22 & (dZ_r(2:Nrow,:)<T_connect) & (dZ_c(:,1:(Ncol-1))<T_connect); % & (Dm2 < T_connect);

end; 				% of smoothing

% At that point the Divij are the final connections

% Find the number of neighbors per points:

t = zeros(Nrow,Ncol);
b = zeros(Nrow,Ncol);
l = zeros(Nrow,Ncol);
r = zeros(Nrow,Ncol);
tr = zeros(Nrow,Ncol);
br = zeros(Nrow,Ncol);
tl = zeros(Nrow,Ncol);
bl = zeros(Nrow,Ncol);

t(2:Nrow,2:Ncol) = (Div12 | Div21);
t(2:Nrow,1:(Ncol-1)) = t(2:Nrow,1:(Ncol-1)) | (Div11 | Div22);
b(1:(Nrow-1),2:Ncol) = (Div12 | Div21);
b(1:(Nrow-1),1:(Ncol-1)) = b(1:(Nrow-1),1:(Ncol-1)) | (Div11 | Div22);
r(2:Nrow,1:(Ncol-1)) = (Div12 | Div22);
r(1:(Nrow-1),1:(Ncol-1)) = r(1:(Nrow-1),1:(Ncol-1)) | (Div11 | Div21);
l(2:Nrow,2:Ncol) = (Div12 | Div22);
l(1:(Nrow-1),2:Ncol) = l(1:(Nrow-1),2:Ncol) | (Div11 | Div21);
tr(2:Nrow,1:(Ncol-1)) = Div11 | Div12;
br(1:(Nrow-1),1:(Ncol-1)) = Div21 | Div22;
tl(2:Nrow,2:Ncol) = Div21 | Div22;
bl(1:(Nrow-1),2:Ncol) = Div11 | Div12;

Nn = t + b + l + r + tr + br + tl + bl;

% build up the matrix of used points: (for renumbering)
% Number of neighbor triangles:

top_left = Div12 + Div21 + Div22;
top_right = Div11 + Div12 + Div22;
bot_left = Div11 + Div12 + Div21;
bot_right = Div11 + Div21 + Div22;

Used_points = zeros(Nrow,Ncol);
Used_points(2:Nrow,2:Ncol) = Used_points(2:Nrow,2:Ncol)+top_left;
Used_points(2:Nrow,1:(Ncol-1)) = Used_points(2:Nrow,1:(Ncol-1))+top_right;
Used_points(1:(Nrow-1),2:Ncol) = Used_points(1:(Nrow-1),2:Ncol)+bot_left;
Used_points(1:(Nrow-1),1:(Ncol-1)) = Used_points(1:(Nrow-1),1:(Ncol-1))+bot_right;

[xc3_2, xp3] = find(Used_points);

ind_points = (xp3-1)*Nrow + xc3_2;

N_vertices = length(ind_points);

Ind_Mat = -ones(Nrow,Ncol);
Ind_Mat(ind_points) = (1:N_vertices)-1;

% Regenere the 3D coordinates, the image location, and the projector coordinates:

Xc3 =  [Xmesh(ind_points) Ymesh(ind_points) Zmesh(ind_points)]';

% number of neighbors:

Nn3 = Nn(ind_points)';

xc3 =  project_points2(Xc3,zeros(3,1),zeros(3,1),fc,cc,kc,alpha_c);% project2(Xc3,eye(3),[0;0;0],fc,cc,kc);

%xc3(2,:) = xc3_2' + xcmin - 1;

%xp3_2D = project_points2(Xc3,om,T,fp,cp,kp,alpha_p);

if spatial,
    xp3 = step_xp * (xp3'-1) + xpmin;
else
    xp3_2D = project_points2(Xc3,om,T,fp,cp,kp,alpha_p);
    xp3 = xp3_2D(1,:);
end;

% Texture coordinates:
xc_texture = [(xc3(1,:)+.5)/(N_x);1-((xc3(2,:)+.5)/(N_y))];

% The boundaries faces (not fully connected):

Div11_u = Div11;
Div12_u = Div12;
Div21_u = Div21;
Div22_u = Div22;

[r11,c11] = find(Div11_u);
[r12,c12] = find(Div12_u);
[r21,c21] = find(Div21_u);
[r22,c22] = find(Div22_u);

% Work with Div11:
ind11_1 =  Nrow*c11+r11;
ind11_2 =  Nrow*(c11-1)+r11;
ind11_3 =  Nrow*(c11-1)+r11+1;

% Work with Div12:
ind12_1 =  Nrow*c12+r12;
ind12_2 =  Nrow*(c12-1)+r12+1;
ind12_3 =  Nrow*c12+r12+1;

% Work with Div21:
ind21_1 =  Nrow*c21+r21;
ind21_2 =  Nrow*(c21-1)+r21;
ind21_3 =  Nrow*c21+r21+1;

% Work with Div22:
ind22_1 =  Nrow*(c22-1)+r22;
ind22_2 =  Nrow*(c22-1)+r22+1;
ind22_3 =  Nrow*c22+r22+1;

% FaceSet
Faces11 = [Ind_Mat(ind11_1) Ind_Mat(ind11_2) Ind_Mat(ind11_3)];
Faces12 = [ Ind_Mat(ind12_1) Ind_Mat(ind12_2) Ind_Mat(ind12_3)];
Faces21 = [Ind_Mat(ind21_1) Ind_Mat(ind21_2) Ind_Mat(ind21_3)];
Faces22 = [ Ind_Mat(ind22_1) Ind_Mat(ind22_2) Ind_Mat(ind22_3)];

Faces = [Faces11;Faces12;Faces21;Faces22]';

N_triangles = size(Faces,2);

% matlab formulation:
tri3 = Faces' + 1;

% Direction of observation:
dc3 = -Xc3;
dc3 = dc3 ./ (ones(3,1) * (sqrt(sum(dc3.^2))));

% use Used_points to keep the number of neighbor triangles:

% compute the surface normals for Div11, Div12, Div21 and Div22, and then
% average them using Used_points as number of them.

nx11 = zeros(Nrow-1,Ncol-1);
ny11 = zeros(Nrow-1,Ncol-1);
nz11 = zeros(Nrow-1,Ncol-1);

nx12 = zeros(Nrow-1,Ncol-1);
ny12 = zeros(Nrow-1,Ncol-1);
nz12 = zeros(Nrow-1,Ncol-1);

nx21 = zeros(Nrow-1,Ncol-1);
ny21 = zeros(Nrow-1,Ncol-1);
nz21 = zeros(Nrow-1,Ncol-1);

nx22 = zeros(Nrow-1,Ncol-1);
ny22 = zeros(Nrow-1,Ncol-1);
nz22 = zeros(Nrow-1,Ncol-1);

u11 = [(Xmesh(ind11_2)-Xmesh(ind11_1))';(Ymesh(ind11_2)-Ymesh(ind11_1))';(Zmesh(ind11_2)-Zmesh(ind11_1))'];
v11 = [(Xmesh(ind11_3)-Xmesh(ind11_1))';(Ymesh(ind11_3)-Ymesh(ind11_1))';(Zmesh(ind11_3)-Zmesh(ind11_1))'];

nx11(r11+(c11-1)*(Nrow-1)) = u11(2,:).*v11(3,:) - u11(3,:).*v11(2,:);
ny11(r11+(c11-1)*(Nrow-1)) = u11(3,:).*v11(1,:) - u11(1,:).*v11(3,:);
nz11(r11+(c11-1)*(Nrow-1)) = u11(1,:).*v11(2,:) - u11(2,:).*v11(1,:);

u12 = [(Xmesh(ind12_2)-Xmesh(ind12_1))';(Ymesh(ind12_2)-Ymesh(ind12_1))';(Zmesh(ind12_2)-Zmesh(ind12_1))'];
v12 = [(Xmesh(ind12_3)-Xmesh(ind12_1))';(Ymesh(ind12_3)-Ymesh(ind12_1))';(Zmesh(ind12_3)-Zmesh(ind12_1))'];

nx12(r12+(c12-1)*(Nrow-1)) = u12(2,:).*v12(3,:) - u12(3,:).*v12(2,:);
ny12(r12+(c12-1)*(Nrow-1)) = u12(3,:).*v12(1,:) - u12(1,:).*v12(3,:);
nz12(r12+(c12-1)*(Nrow-1)) = u12(1,:).*v12(2,:) - u12(2,:).*v12(1,:);

u21 = [(Xmesh(ind21_2)-Xmesh(ind21_1))';(Ymesh(ind21_2)-Ymesh(ind21_1))';(Zmesh(ind21_2)-Zmesh(ind21_1))'];
v21 = [(Xmesh(ind21_3)-Xmesh(ind21_1))';(Ymesh(ind21_3)-Ymesh(ind21_1))';(Zmesh(ind21_3)-Zmesh(ind21_1))'];

nx21(r21+(c21-1)*(Nrow-1)) = u21(2,:).*v21(3,:) - u21(3,:).*v21(2,:);
ny21(r21+(c21-1)*(Nrow-1)) = u21(3,:).*v21(1,:) - u21(1,:).*v21(3,:);
nz21(r21+(c21-1)*(Nrow-1)) = u21(1,:).*v21(2,:) - u21(2,:).*v21(1,:);

u22 = [(Xmesh(ind22_2)-Xmesh(ind22_1))';(Ymesh(ind22_2)-Ymesh(ind22_1))';(Zmesh(ind22_2)-Zmesh(ind22_1))'];
v22 = [(Xmesh(ind22_3)-Xmesh(ind22_1))';(Ymesh(ind22_3)-Ymesh(ind22_1))';(Zmesh(ind22_3)-Zmesh(ind22_1))'];

nx22(r22+(c22-1)*(Nrow-1)) = u22(2,:).*v22(3,:) - u22(3,:).*v22(2,:);
ny22(r22+(c22-1)*(Nrow-1)) = u22(3,:).*v22(1,:) - u22(1,:).*v22(3,:);
nz22(r22+(c22-1)*(Nrow-1)) = u22(1,:).*v22(2,:) - u22(2,:).*v22(1,:);

% Sum all the relevant normal components for each vertice:

nx = zeros(Nrow,Ncol);
ny = zeros(Nrow,Ncol);
nz = zeros(Nrow,Ncol);

nx(1:(Nrow-1),1:(Ncol-1)) = nx11 + nx21 + nx22;
nx(2:Nrow,1:(Ncol-1)) = nx(2:Nrow,1:(Ncol-1)) + nx11 + nx12 + nx22;
nx(1:(Nrow-1),2:Ncol) = nx(1:(Nrow-1),2:Ncol) + nx11 + nx12 + nx21;
nx(2:Nrow,2:Ncol) = nx(2:Nrow,2:Ncol) + nx12 + nx21 + nx22;

ny(1:(Nrow-1),1:(Ncol-1)) = ny11 + ny21 + ny22;
ny(2:Nrow,1:(Ncol-1)) = ny(2:Nrow,1:(Ncol-1)) + ny11 + ny12 + ny22;
ny(1:(Nrow-1),2:Ncol) = ny(1:(Nrow-1),2:Ncol) + ny11 + ny12 + ny21;
ny(2:Nrow,2:Ncol) = ny(2:Nrow,2:Ncol) + ny12 + ny21 + ny22;

nz(1:(Nrow-1),1:(Ncol-1)) = nz11 + nz21 + nz22;
nz(2:Nrow,1:(Ncol-1)) = nz(2:Nrow,1:(Ncol-1)) + nz11 + nz12 + nz22;
nz(1:(Nrow-1),2:Ncol) = nz(1:(Nrow-1),2:Ncol) + nz11 + nz12 + nz21;
nz(2:Nrow,2:Ncol) = nz(2:Nrow,2:Ncol) + nz12 + nz21 + nz22;

% Vertice normals:
nc3 =  [nx(ind_points)';ny(ind_points)';nz(ind_points)'];

% normalization of the normals:

conf_nc3 = sum(nc3.^2);
Norms = sqrt(conf_nc3);

nc3 = nc3 ./ (ones(3,1)*Norms);

% update of the normal components:

nx(ind_points) = nc3(1,:);
ny(ind_points) = nc3(2,:);
nz(ind_points) = nc3(3,:);

%clear Ambi C1 C2 Cmesh D1 D2 Div Div1 Div11 Div11_u Div12 Div12_u Div2 Div21 Div21_u Div22 Div22_u Div1n Div2n Dm1 Dm2 F1 F2 F3 F4 Faces Faces11 Faces12 Faces21 Faces22 Header Ind_Mat Light Ncol Nn Norms Nrow Stripe T_connect Used_points Vertice_norm XX Xc3 Xmesh YY Ymesh ZZ Zmesh b bl bot_left bot_right br c11 c12 c21 c22 coeff dZ_c dZ_r day dot file i ind ind11_1 ind11_2 ind11_3 ind12_1 ind12_2 ind12_3 ind21_3 ind21_1 ind21_2 ind21_3 ind22_1 ind22_2 ind22_3 ind_col ind_row indd is js l ni nx nx11 nx12 nx21 nx22 ny ny11 ny12 ny21 ny22 nz nz11 nz12 nz21 nz22 r r11 r12 r21 r22 rasterimage sm sm_b sm_bl sm_br sm_l sm_r sm_t sm_tl sm_tr t tl top_left top_right tr u11 u12 u21 u22 unshading v v11 v12 v21 v22 Xc3 xc3_2 xc_texture xcmax xcmesh xcmin xp3 xpmax xpmin ycmesh xc3 ind_points

repo: ekfmonoslam


text:
function hermite_cubic_test09 ( )

%*****************************************************************************80
%
%% HERMITE_CUBIC_TEST09 tests HERMITE_CUBIC_SPLINE_INTEGRATE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 February 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HERMITE_CUBIC_TEST09:\n' );
  fprintf ( 1, '  HERMITE_CUBIC_SPLINE_INTEGRATE integrates a Hermite\n' );
  fprintf ( 1, '  cubic spline from A to B.\n' );
%
%  Define the cubic spline.
%
  x1 = 0.0;
  x2 = 10.0;

  nn = 11;
  xn = linspace ( x1, x2, nn );
  [ fn, dn, sn, tn ] = cubic_value ( xn );

  n = 25;
  a = linspace ( 2.5, 2.5, n );
  b = linspace ( x1 - 1.0, x2 + 1.0, n );

  q = hermite_cubic_spline_integrate ( nn, xn, fn, dn, n, a, b );

  fprintf ( 1, '\n' );
  fprintf ( 1, '                                 Exact       Computed\n' );
  fprintf ( 1, '    I      A           B         Integral    Integral\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n

    q_exact = cubic_integrate ( a(i), b(i) );

    fprintf ( 1, '  %3d  %10f  %10f  %10.6g  %10.6g\n', ...
      i, a(i), b(i), q_exact, q(i) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [x_best, psi_best, out] = LevMar(mapp, lin_sym_solver, x0, options)
% LevMar is a Levenberg-Marquardt trust-region algorithm for solving
% systems of nonlinear equations :math:`h(x) = 0`, `x` in :math:`R^m`
% using the nonlinear unconstrained minimization :math:`\textrm{min}\ \psi(x) = 1/2 ||h(x)||^2`
% s.t. `x` in :math:`R^m`.
%
% USAGE:
%
%    [x_best, psi_best, out] = LevMar(mapp, lin_sym_solver, x0, options)
%
% INPUTS:
%    mapp:              function handle provides `h(x)` and gradient `h(x)`
%    lin_sym_solver:    function handle for solving the linear system
%    x0:                initial point
%    options:           structure including the parameteres of scheme
%
%                         * .eta - parameter of the scheme
%                         * .MaxNumIter - maximum number of iterations
%                         * .MaxNumMapEval - maximum number of function evaluations
%                         * .MaxNumGmapEval - maximum number of subgradient evaluations
%                         * .TimeLimit - maximum running time
%                         * .epsilon - accuracy parameter
%                         * .x_opt - optimizer
%                         * .psi_opt - optimum
%                         * .adaptive - update lambda adaptively
%                         * .flag_x_error - 1: saves :math:`x_{error}`, 0: do not saves :math:`x_{error}` (default)
%                         * .flag_psi_error - 1: saves :math:`\psi_{error}`, 0: do not saves :math:`\psi_{error}` (default)
%                         * .flag_time - 1: saves :math:`\psi_{error}`, 0: do not saves :math:`\psi_{error}` (default)
%                         * .Stopping_Crit - stopping criterion
%
%                           1. stop if :math:`||grad|| \leq \epsilon`
%                           2. stop if :math:`||nhxk|| \leq \epsilon`
%                           3. stop if `MaxNumIter` is reached
%                           4. stop if `MaxNumMapEval` is reached
%                           5. stop if `MaxNumGmapEval` is reached
%                           6. stop if `TimeLimit` is reached
%                           7. stop if :math:`||grad|| \leq \textrm{max}(\epsilon, \epsilon^2 * ngradx0)`
%                           8. stop if :math:`||nhxk|| \leq \textrm{max}(\epsilon, \epsilon^2 * nhx0)`
%                           9. stop if (default) :math:`||hxk|| \leq \epsilon` or `MaxNumIter` is reached
%
% OUTPUT:
%    x_best:            the best approximation of the optimizer
%    psi_best:          the best approximation of the optimum
%    out:               structure including more output information
%
%                         * .T - running time
%                         * .Niter - total number of iterations
%                         * .Nmap - total number of mapping evaluations
%                         * .Ngmap - total number of mapping gradient evaluations
%                         * .merit_func - array including all merit function values
%                         * .x_error - relative error :math:`\textrm{norm}(x_k(:)-x_{opt}(:))/\textrm{norm}(x_{opt})`
%                         * .psi_error - relative error :math:`(\psi_k-\psi_{opt})/(\psi_0-\psi_{opt}))`
%                         * .Status - reason of termination
%
% .. REFERENCE:
% .. 1. I.C.F. Ipsen, C.T. Kelley, S.R. Pope, Rank-deficient nonlinear least squares problems and subset selection, SIAM Journal on Numerical Analysis, 49(3), 1244-1266 (2011
% .. 2. C.T. Kelley, Iterative Methods for Optimization. SIAM Press, Philadelphia, 1999.
% .. Author: - Masoud Ahookhosh, System Biochemistry Group, Luxembourg Center for System Biomedicine, University of Luxembourg, Luxembourg
%            - Update: July 2017, M. Ahookhosh

format longG ;

% ================ Error messages for input and output =================
if nargin > 4
    error('The number of input arguments is more than what is needed');
elseif nargin < 4
    error('The number of input arguments is not enough');
end;

if isempty(mapp)
    error('the function handle mapp has to be defined');
elseif ~isa(mapp,'function_handle')
    error('mapp should be a function handle');
end

if isempty(lin_sym_solver)
    error('the function handle lin_sym_solver has to be defined');
elseif ~isa(lin_sym_solver,'function_handle')
    error('lin_sym_solver should be a function handle');
end

if isempty(x0)
    error('The starting point x0 has to be defined');
elseif ~isa(x0,'numeric')
    error('x0 should be a numeric vector');
end

% =================== initializing the parameters ======================
% ===== user has requested viewing the default values of "options" =====
[eta,epsilon,MaxNumIter,MaxNumMapEval,MaxNumGmapEval,adaptive, ...
    TimeLimit,flag_x_error,flag_psi_error,flag_time,Stopping_Crit] ...
    = Initialization(options);

if ~isa(eta,'numeric') || (eta <= 0)
    error('eta should be numeric and eta in (0,4*delta)');
end

if isfield(options,'x_opt')
    x_opt=options.x_opt;
elseif flag_x_error==1
    error('x_error requires to x_opt be specified');
end

if flag_x_error == 1
    Nxopt      = sqrt(sum(x_opt(:).^2));
    x_error(1) = sqrt(sum((x0(:)-x_opt(:)).^2))/Nxopt;
end

if flag_psi_error == 1
    psi_error(1) = 1;
end

if flag_time == 1
    Time(1) = 0;
end

mu0        = 1e-4;
mulow      = 0.25;
muhigh     = 0.75;
wdown      = 0.5;
wup        = 2;
nu0        = 1e-3;

xk         = x0;
Niter      = 1;
[hxk,ghxk] = mapp(x0);
Nmap       = 1;
Ngmap      = 1;
grad       = ghxk*hxk;
ngradx0    = norm(grad);
nhx0       = norm(hxk);
nhxk       = nhx0;
I          = eye(length(xk));
psik       = 0.5*nhxk^2;
merit_func = psik;
nuk        = ngradx0;
StopFlag   = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%% Main body of LevMar.m %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T0 = tic;

% ======================= start of the main loop =======================
while ~StopFlag

    Gk    = ghxk*ghxk';
    Hk    = Gk+nuk*I;
    dk    = lin_sym_solver(Hk,grad);
    xk1   = xk+dk;
    hxk1  = mapp(xk1);

    Nmap  = Nmap+1;

    flag   = 0;
    initer = 0;
    while ~flag
        nhxk1  = norm(hxk1);
        psik1 = 0.5*nhxk1^2;
        ared  = psik-psik1;
        pred  = -0.5*(grad'*dk);
        rk    = ared/pred;
        if rk < mu0
            nuk    = max(nuk*wup,nu0);
		    Hk     = Gk+nuk*I;
            dk     = lin_sym_solver(Hk,grad);
            xk1    = xk+dk;
            hxk1   = mapp(xk1);
            Nmap   = Nmap+1;
            initer = initer+1;
            if initer>30
                flag = 1;
            end
	    elseif rk < mulow
	    	xk    = xk1;
            Niter = Niter+1;
            nuk   = max(nuk*wup,nu0);
            flag  = 1;
        else
		    xk    = xk1;
            Niter = Niter+1;
            if rk > muhigh
		        nuk = wdown*nuk;
                %if nuk < 10^(-3)
                if nuk < 10^(-3)
                    nuk =1e-8;
                end
            end
            flag = 1;
        end
    end

    [hxk,ghxk] = mapp(xk);
    Nmap       = Nmap+1;
    grad       = ghxk*hxk;
    nhxk       = norm(hxk);

    % ================= Gathering output information ===================
    psik              = 0.5*nhxk^2;
    merit_func(Niter) = psik;
    if flag_time == 1
        Time(Niter+1) = toc(T0);
    end

    if flag_x_error == 1
        Nx_opt = norm(x_opt);
        x_error(Niter+1) = sqrt(sum((xk(:)-x_opt(:)).^2))/Nx_opt;
    end

    if flag_psi_error == 1
        psi_error(Niter+1) = (psik-psi_opt)/(psi0-psi_opt);
    end

    % ================== checking stopping criteria ====================
    T = toc(T0);

    [StopFlag,Status] = StopCriterion(grad,nhxk,Niter,Nmap, ...
    Ngmap,MaxNumIter,MaxNumMapEval,MaxNumGmapEval,T,TimeLimit, ...
    epsilon,nhx0,ngradx0,Stopping_Crit);

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Outputs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Status
x_best         = xk;
psi_best       = psik;
out.T          = T;
out.nhx        = nhxk;
out.merit_func = merit_func';
out.Niter      = Niter;
out.Nmap       = Nmap;
out.Ngmap      = Ngmap;
out.Status     = Status;

if flag_x_error == 1
    out.x_error = x_error;
end
if flag_psi_error == 1
    out.psi_error = psi_error;
end
if flag_time == 1
    out.Time = Time;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%% End of LevMar.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: opencobra

repo: cobratoolbox


text:
function blas1_z_test10 ( )

%*****************************************************************************80
%
%% TEST10 tests ZDROT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 May 2006
%
%  Author:
%
%    John Burkardt
%
  n = 6;

  for j = 1 : n
    x(j) = 10 * j + j * i;
  end

  for j = 1 : n
    y(j) = 20 * j + 2 * j * i;
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST10\n' );
  fprintf ( 1, '  ZDROT carries out a Givens rotation\n' );
  fprintf ( 1, '  on a complex vector.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  X and Y\n' );
  fprintf ( 1, '\n' );
  for j = 1 : n
    fprintf ( 1, '  %6d  %10f  %10f  %10f  %10f\n', ...
      j, real ( x(j) ), imag ( x(j) ), real ( y(j) ), imag ( y(j) ) );
  end

  c = 0.5;
  s = sqrt ( 1.0 - c * c );
  [ x, y ] = zdrot ( n, x, 1, y, 1, c, s );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  ZDROT ( N, X, 1, Y, 1, %f, %f )\n', c, s );
  fprintf ( 1, '\n' );
  for j = 1 : n
    fprintf ( 1, '  %6d  %10f  %10f  %10f  %10f\n', ...
      j, real ( x(j) ), imag ( x(j) ), real ( y(j) ), imag ( y(j) ) );
  end

  return
end


text:
function xyzl_data_print ( point_num, line_num, line_data_num, ...
  line_pointer, line_data )

%*****************************************************************************80
%
%% XYZL_DATA_PRINT prints the data of an XYZL file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, integer LINE_NUM, the number of lines.
%
%    Input, integer LINE_DATA_NUM, the number of line items.
%
%    Input, integer LINE_POINTER(LINE_NUM+1), pointers to the
%    first line item for each line.
%
%    Input, integer LINE_DATA(LINE_DATA_NUM), indices
%    of points that form lines.
%
  fprintf ( 1, '\n' );
  for line = 1 : line_num
    fprintf ( 1, '  %4d  %8d  %8d\n', ...
      line, line_pointer(line) : line_pointer(line+1) - 1 );
  end

  fprintf ( 1, '\n' );
  for line = 1 : line_num
    for i = line_pointer(line) : line_pointer(line+1) - 1
      fprintf ( 1, '  %d', line_data(i) );
    end
    fprintf ( 1, '\n' );
  end

  return
end


text:
function ssimLabelFeatures(ftrDir,vocabPath,labelDir,imgList,parms)

  load(vocabPath); % Gives the variable V

  for i=1:length(imgList),
    IMGname=imgList{i};
    onlyImgName=extractImgName(IMGname);
    fprintf('Labeling ssim descriptors of %s\n',onlyImgName);
    load([ftrDir onlyImgName '.mat']);
    textonlabels=uint16(MEXfindnearestl2(features,V));
    %save([RESdirs(j,:),IMGname],'textonlabels');

    origImg=ssimReadImg(IMGname,parms);
    [labelImage,cmap]=drawLabels(textonlabels,size(V,2),ftrCoords,salientCoords,uniformCoords,origImg);

    %height=size(labelImage,1);width=size(labelImage,2);
    %labelImage=labelImage( (1+parms.margin):(height-parms.margin),(1+parms.margin):(width-parms.margin));
    save([labelDir onlyImgName '.mat'],'labelImage');

  end;
end

function [labelImage,cmap]=drawLabels(labels,numClusters,dCoords,sCoords,uCoords,origImg)
  numLabels=numClusters;

  if size(labels,2)~=size(dCoords,2)
    fprintf('Something wrong with dimensions of arrays of labels and drawing coordinates\n');keyboard;exit;
  end;

  labelImage=zeros(size(origImg,1),size(origImg,2)); % Intialize with the unlabeled tag

  %for i=1:size(dCoords,2)
    %labelImage(dCoords(2,i),dCoords(1,i))=labels(i);
  %end
  labelImage(sub2ind(size(labelImage),dCoords(2,:),dCoords(1,:)))=labels;

  cmap=colorcube(numLabels);

end

%function [labelImage,cmap]=drawLabels(labels,numClusters,dCoords,sCoords,uCoords,origImg)
  %numLabels=numClusters+2; % Other than texton labels, we have salient,uniform
  %labels=labels+2; % The first 2 labels are for uniform, and salient (in order)
%
  %if size(labels,2)~=size(dCoords,2)
    %fprintf('Something wrong with dimensions of arrays of labels and drawing coordinates\n');keyboard;exit;
  %end;
%
  %labelImage=zeros(size(origImg,1),size(origImg,2)); % Intialize with the unlabeled tag
%
  %%for i=1:size(dCoords,2)
    %%labelImage(dCoords(2,i),dCoords(1,i))=labels(i);
  %%end
  %labelImage(sub2ind(size(labelImage),dCoords(2,:),dCoords(1,:)))=labels;
%
  %for i=1:size(uCoords,2)
    %labelImage(uCoords(2,i),uCoords(1,i))=1;
  %end
%
  %for i=1:size(sCoords,2)
    %labelImage(sCoords(2,i),sCoords(1,i))=2;
  %end
  %
  %cmap=colorcube(numLabels);
  %
%end

author: adikhosla

repo: feature-extraction


text:
function newimg = GPReduce2(img,blocksize,displayflag)

if ~exist('displayflag')
	displayflag = 1;
end

dim = size(img);

newimg = [];
for b = 1:blocksize:dim(3)
    if b<dim(3)
        img2 = img(:,:,b:b+blocksize-1);
    else
        img2 = img(:,:,b);
    end
	newimg2 = GPReduce(img2,displayflag);
	newimg = cat(3,newimg,newimg2);
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function testbed_stationplot(Q, varargin)

options = struct(...
    'stations', 1:size(Q.W,1), ...
    'linespec', 'x');

% Check arguments
[options, errmsg, remopts] = argparse(options, varargin{:});
error(errmsg);

mapplot(Q.options.user_data.coordinates(options.stations,1:2)', options.linespec, ...
        remopts{:});

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function imageWrapper = DMLoadImagesFromMetadataGrouping(metadataGrouping, dicomLibrary, reporting)
    % DMLoadImagesFromMetadataGrouping. Loads metadata from a series of DICOM files
    %
    %     Syntax
    %     ------
    %
    %         image_volume = DMLoadImagesFromMetadataGrouping(metadata_list, reporting)
    %
    %             imageWrapper     a CoreWrapper object containing the image volume
    %
    %             metadataGrouping a DMFileGrouping containing metadata
    %
    %             dicomLibrary     (Optional) An object implementing DMDicomLibraryInterface
    %
    %             reporting       (Optional) A CoreReporting or other implementor of CoreReportingInterface,
    %                             for error and progress reporting. Create a CoreReporting
    %                             with no arguments to hide all reporting. If no
    %                             reporting object is specified then a default
    %                             reporting object with progress dialog is
    %                             created
    %
    %
    %     Licence
    %     -------
    %     Part of DicoMat. https://github.com/tomdoel/dicomat
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the BSD 3-Clause license. Please see the file LICENSE for details.
    %

    if nargin < 3
        reporting = CoreReportingDefault;
    end

    if nargin < 2
        dicomLibrary = DMDicomLibrary.getLibrary;
    end

    reporting.ShowProgress('Reading pixel data');
    reporting.UpdateProgressValue(0);

    imageWrapper = CoreWrapper;

    num_slices = length(metadataGrouping.Metadata);

    % Load image slice
    first_image_slice = dicomLibrary.dicomread(metadataGrouping.Metadata{1});
    if isempty(first_image_slice)
        return;
    end

    % Pre-allocate image matrix
    size_i = metadataGrouping.Metadata{1}.Rows;
    size_j = metadataGrouping.Metadata{1}.Columns;
    size_k = num_slices;
    samples_per_pixel = metadataGrouping.Metadata{1}.SamplesPerPixel;

    % Pre-allocate image matrix
    data_type = whos('first_image_slice');
    data_type_class = data_type.class;
    if (strcmp(data_type_class, 'char'))
        reporting.ShowMessage('DMLoadImagesFromMetadataGrouping:SettingDatatypeToInt8', 'Char datatype detected. Setting to int8');
        data_type_class = 'int8';
    end
    imageWrapper.RawImage = zeros([size_i, size_j, size_k, samples_per_pixel], data_type_class);
    imageWrapper.RawImage(:, :, 1, :) = first_image_slice;

    for file_index = 2 : num_slices
        imageWrapper.RawImage(:, :, file_index, :) = dicomLibrary.dicomread(metadataGrouping.Metadata{file_index});
        reporting.UpdateProgressValue(round(100*(file_index)/num_slices));
    end

    reporting.CompleteProgress;
end

author: tomdoel

repo: pulmonarytoolkit


text:
function captureFigure( figh, filename, trim )
%captureFigure: use print to capture a figure with over-sampling
%
%   gpubench.captureFigure(figh,filename,trim)

%   Copyright 2011 The MathWorks, Inc.

if nargin<3
    trim = true;
end

% Have to hard-code the DPI since Windows doesn't report it correctly
% dpi = get(0,'ScreenPixelsPerInch');
%dpi = 192;
dpi = 160;

% Copy the position to the paper position to get the size right
pos = get( figh, 'Position' );
set( figh, 'PaperUnits', 'Points', 'PaperPosition', pos );

% Print the image to file at screen resolution
print( figh, '-zbuffer', sprintf('-r%d',round(0.5*dpi)), '-dpng', filename );

% Load and trim the top whitespace
gpubench.trimImage( filename, trim );

repo: Aero-matlab


text:
% VL_DEMO_SIFT_VS_UBC Compare VLFeat SIFT with Lowe's original

% --------------------------------------------------------------------
%                            Load a figure and original SIFT keypoints
% --------------------------------------------------------------------

im = imread(fullfile(vl_root, 'data', 'box.pgm')) ;
[f1,d1] = vl_ubcread(fullfile(vl_root, 'data', 'box.sift')) ;

% --------------------------------------------------------------------
%                                    Compare with VLFeat SIFT detector
% --------------------------------------------------------------------

f2 = vl_sift(im2single(im), ...
             'firstoctave', -1, ...
             'peakthresh', .01, ...
             'edgethresh', 10, ...
             'windowsize', 2, ...
             'verbose') ;

D = sqrt(vl_alldist2(f1(1:2,:), f2(1:2,:))) ;
[d12,m12] = min(D,[],2) ;
[d21,m21] = min(D,[],1) ;

matches = [1:size(f1,2), m21 ; m12',  1:size(f2,2)] ;
proxim  = [d12', d21] ;

h = histc(proxim, [0 .01 .05 +inf]) ;
h = h / sum(h) * 100 ;
h = h(1:end-1) ;

figure(1) ; clf ;
imagesc(im) ; colormap(gray(256)) ; hold on ;
vl_plotframe(f1, 'linewidth', 3, 'color', 'r') ;
vl_plotframe(f2, 'linewidth',  1, 'color', 'b') ;
axis image off ;
vl_demo_print('sift_vs_ubc_1', 0.7) ;

figure(2) ; clf ;
pie(h) ;
colormap(hot(3)) ;
legend({'0.01 pixels', '0.05 pixels', 'others'}, ...
       'location', 'northeastoutside') ;
set(findobj(2, '-property', 'fontsize'), 'fontsize', 11) ;
vl_demo_print('sift_vs_ubc_2') ;

% --------------------------------------------------------------------
%                                  Compare with VLFeat SIFT descriptor
% --------------------------------------------------------------------

[drop,d2]=vl_sift(im2single(im), 'frames', f1, 'verbose', 'firstoctave', -1) ;

D = sqrt(double(vl_alldist2(d1, d2, 'l2'))) ;
[d12,m12] = min(D,[],2) ;
[d21,m21] = min(D,[],1) ;

matches = [1:size(d1,2), m21 ; m12',  1:size(d2,2)] ;
proxim  = [d12', d21] ;

meanDist = mean(D(:)) ;

h = histc(proxim, meanDist * [0 .05 .10 .20 +inf]) ;
h = h / sum(h) * 100 ;
h = h(1:end-1) ;

figure(3) ; clf ;
pie(h) ;
colormap(hot(4)) ;
legend({'5% differnce', '10% difference', '20% difference', 'others'}, ...
       'location', 'northeastoutside') ;
set(findobj(3, '-property', 'fontsize'), 'fontsize', 11) ;
vl_demo_print('sift_vs_ubc_3') ;

author: jianxiongxiao

repo: ProfXkit


text:
function tetris(action)

% TETRIS for MATLAB
%
%	To begin, simply type: tetris
%
%	Only one game can be run at a time.
%	UN-PAUSE (I know that's not a word) the game by hitting RETURN
%	in the Matlab command window.
%
%	Please feel free to modify, improve, rewrite, cleanup the code etc.
%	Written and tested on a SUN running Matlab 4.2c, but should work with
%	Matlab 4.0, if not you can debug/modify it.
%	This is FREEWARE, enjoy!   -  Peace & Love

%
%	Author: Maurice Smith, 1/19/95 (msmith@bme.jhu.edu)

global CENTER_X CENTER_Y POS STUCK USED_SPACE H;
global WIDTH1 HEIGHT1 LEVEL TIME_DELAY SQ_X SQ_Y;
global H_SLIDER_LEVEL H_SLIDER_WIDTH H_SLIDER_HEIGHT;
global H_TEXT_LEVEL H_TEXT_WIDTH H_TEXT_HEIGHT H_TEXT_SCORE;
global START_HANDLE STOP_HANDLE STOP_GAME EXIT_HANDLE;
global H_A H_B BOX H_LINE MIN_DELAY MAX_DELAY;

if nargin < 1, action = 'initialize'; end;

if strcmp(action,'initialize'),

default_width = 10;		% default width of tetris play-field
default_height = 25;		% default hieght of tetris play-field
default_level = 5;		% default skill level
default_box = 1;		% default play field outline (on=1, off=0)

WIDTH1 = default_width;
HEIGHT1 = default_height;
LEVEL = default_level;
BOX = default_box;

MIN_DELAY = 0.1;		% delay corresponding to skill level 9
MAX_DELAY = 1;			% delay corresponding to skill level 0
TIME_DELAY = MAX_DELAY*exp((LEVEL/9)*log(MIN_DELAY/MAX_DELAY));
fig=figure( ...
        'Name','Tetris for Matlab', 'NumberTitle','off', ...
        'Visible','off', 'BackingStore','off');

SQ_X = [0 -1 -1 0]';	% vertex coordinates of unit square (for "patch" command)
SQ_Y = [0 0 -1 -1]';

on = 'set(go_handle,''back'',''cyan''), ';
off = 'set(go_handle,''back'',''default''), ';

figure(fig);
H_A = axes( ...
        'Units','normalized', 'Position', [.1 .18 .4 .8], ...
        'Visible','off', 'DrawMode','fast', ...
        'NextPlot','replace', 'XLim',[1 WIDTH1+1],'YLim',[0 HEIGHT1]);
H_B = axes( ...
        'Units','normalized', 'Position', [.1 .18 .4 .8],...
        'Visible','off', 'DrawMode','fast', ...
        'NextPlot','replace', 'XLim',[1 WIDTH1+1],'YLim',[0 HEIGHT1]);

if BOX == 1,			% draw outline box initially
	BOX = ~BOX;
	tetris('box');
end

pause_handle = uicontrol('units','normalized',...
	'position', [.7 .04 .1 .05],'string','pause', ...
	'callback','disp(''press RETURN to continue''); pause;', ...
	'interruptible','no');

box_handle = uicontrol('units','normalized',...
	'position', [.7 .15 .1 .05],'string','box', ...
	'callback','tetris(''box'');', ...
	'interruptible','yes');

STOP_HANDLE = uicontrol('units','normalized',...
	'position', [.55 .04 .1 .05],'string','stop', ...
	'callback','tetris(''stop'')', ...
	'interruptible','yes','visible','off');

EXIT_HANDLE = uicontrol('units','normalized',...
	'position', [.85 .04 .1 .05],'string','exit', ...
	'callback','delete(gcf)', ...
	'interruptible','no');

START_HANDLE = uicontrol('units','normalized',...
	'position', [.55 .04 .1 .05],'string','start', ...
	'callback','tetris(''start'')', ...
	'interruptible','yes','visible','on');

left_handle = uicontrol('units','normalized',...
	'position', [.18 .07 .1 .05 ],'string','<---', ...
	'callback','tetris(''left'')', ...
	'interruptible','yes');

right_handle = uicontrol('units','normalized',...
	'position', [.32 .07 .1  .05 ],'string','--->', ...
	'callback','tetris(''right'')', ...
	'interruptible','yes');

drop_handle = uicontrol('units','normalized',...
	'position',[.25 .01 .1 .05 ],'string','drop', ...
	'callback','tetris(''drop'')');

rotate_handle = uicontrol('units','normalized',...
	'position',[.25 .13 .1 .05],'string','rotate', ...
	'callback','tetris(''rotate'')');

H_TEXT_SCORE = uicontrol(...
		'ForegroundColor',[ 1 1 1 ],...
		'Position',[ 0.65 0.85 0.2 0.07 ],...
		'String','Score: xxxxx',...
		'Units','normalized');

v = version;
if str2num(v(1:3)) >= 4.2,	% if Matlab 4.2 or later, create slider controls

H_SLIDER_LEVEL = uicontrol(...
		'CallBack', 'tetris(''level''); ', 'Max',[ 9 ], 'Min',[ 0 ],...
		'val', LEVEL, 'Position', [ 0.6 0.3 0.04 0.4 ],...
		'Style', 'slider', 'Units', 'normalized');

H_SLIDER_HEIGHT = uicontrol(...
		'CallBack', 'tetris(''height''); ', 'Max', [ 50 ], 'Min', [ 10 ],...
		'val', HEIGHT1, 'Position', [ 0.78 0.3 0.04 0.4 ],...
		'Style', 'slider', 'Units', 'normalized');

H_SLIDER_WIDTH = uicontrol(...
		'CallBack','tetris(''width''); ','Max',[ 20 ],'Min',[ 6 ],...
		'val',WIDTH1,'Position',[ 0.9 0.3 0.04 0.4 ],...
		'Style','slider','Units','normalized');

H_label_width = uicontrol(...
		'ForegroundColor',[ 0 1 1 ], 'Position',[0.87 .72 .1 .06],...
		'String','Width', 'Style','text', 'Units','normalized');

H_label_Height = uicontrol(...
		'ForegroundColor',[ 1 1 0 ], 'Position',[ 0.75 0.72 0.1 0.06 ],...
		'String','Height', 'Style','text', 'Units','normalized');

H_label_level = uicontrol(...
		'ForegroundColor',[ 1 0 0 ], 'Position',[ 0.57 0.72 0.1 0.06 ],...
		'String','Level', 'Style','text', 'Units','normalized');

H_TEXT_WIDTH = uicontrol(...
		'ForegroundColor',[ 0 1 1 ], 'Position',[0.87 .22 .1 .06],...
		'String',num2str(WIDTH1), 'Style','text', 'Units','normalized');

H_TEXT_HEIGHT = uicontrol(...
		'ForegroundColor',[ 1 1 0 ], 'Position',[ 0.75 0.22 0.1 0.06 ],...
		'String',num2str(HEIGHT1), 'Style','text', 'Units','normalized');

H_TEXT_LEVEL = uicontrol(...
		'ForegroundColor',[ 1 0 0 ], 'Position',[ 0.57 0.22 0.1 0.06 ],...
		'String',num2str(LEVEL), 'Style','text', 'Units','normalized');

end;    % end version 4.2 extras (slider controls & labels)

end; 	% end initialization routine

if strcmp(action,'level'),				% set skill level
	LEVEL = round(get(H_SLIDER_LEVEL,'val'));
	set(H_TEXT_LEVEL,'string',LEVEL);
	TIME_DELAY = MAX_DELAY*exp((LEVEL/9)*log(MIN_DELAY/MAX_DELAY));
end;

if strcmp(action,'height'),			% set play-field height
	HEIGHT1 = round(get(H_SLIDER_HEIGHT,'val'));
	set(H_TEXT_HEIGHT,'string',HEIGHT1);
end;

if strcmp(action,'width'),			% set play-field width
	WIDTH1 = round(get(H_SLIDER_WIDTH,'val'));
	set(H_TEXT_WIDTH,'string',WIDTH1);
end;

if strcmp(action,'box'),		% toggle display of play-field outline
	BOX = ~BOX;
	if BOX == 1,
		H_LINE = line ([1 WIDTH1+1 WIDTH1+1 1 1], [1 1 HEIGHT1 HEIGHT1 1]);
	else, delete(H_LINE);
	end
end;

if strcmp(action,'left'),		% Move Left
	stopped = 0;
	for i=1:4, if USED_SPACE(POS(i,1)-1,POS(i,2)) == 1, stopped=1; end; end;
	if stopped == 0,
		POS(:,1) = POS(:,1) - 1;
		CENTER_X = CENTER_X - 1;
		for i=1:4, set(H(i),'XData',SQ_X+POS(i,1),'YData',SQ_Y+POS(i,2)); end;
		drawnow;
	end;
end;

if strcmp(action,'right'),		% Move Right
	stopped = 0;
	for i=1:4, if USED_SPACE(POS(i,1)+1,POS(i,2)) == 1, stopped=1; end; end;

	if stopped == 0,
		POS(:,1) = POS(:,1) + 1;
		CENTER_X = CENTER_X + 1;
		for i=1:4, set(H(i),'XData',SQ_X+POS(i,1),'YData',SQ_Y+POS(i,2)); end;
		drawnow;
	end;
end;

if strcmp(action,'drop'),		% Drop Piece Down
	while STUCK==0,
		for i=1:4, if USED_SPACE(POS(i,1),POS(i,2)-1) == 1, STUCK=1; end; end;
		if STUCK == 0, POS(:,2) = POS(:,2) -1; end;
	end;
	for i=1:4, set(H(i),'XData',SQ_X+POS(i,1),'YData',SQ_Y+POS(i,2)); end;
	drawnow;

end;

if strcmp(action,'rotate'),			% Rotate
	new_x = -(POS(:,2) - CENTER_Y) + CENTER_X;
	new_y = (POS(:,1) - CENTER_X) + CENTER_Y;
	stopped = 0;
	for i=1:4,
		if ((new_x>1) & (new_x<(WIDTH1+2)) & (new_y>1)),
			if USED_SPACE(new_x(i), new_y(i)) == 1, stopped=1; end;
		else, stopped=1;
		end;
	end;
	if stopped == 0,			% space not occupied
		POS = [new_x new_y];
		for i=1:4, set(H(i),'XData',SQ_X+POS(i,1),'YData',SQ_Y+POS(i,2)); end;
		drawnow;
	end;
end;

if strcmp(action,'stop'),             % stop game
	STOP_GAME = 1;
end;

% start tetris game

if strcmp(action,'start'),             % start game

STOP_GAME = 0;
score = 0;
set(START_HANDLE, 'visible', 'off');
set(STOP_HANDLE, 'visible', 'on');
set(EXIT_HANDLE, 'visible', 'off');
set(H_TEXT_SCORE, 'string', ['Score: ' num2str(score)]);

h_matb = zeros(WIDTH1+2,HEIGHT1);
USED_SPACE = zeros(WIDTH1+2,HEIGHT1);
USED_SPACE(:,1) = ones(WIDTH1+2,1);
USED_SPACE ([1 WIDTH1+2],:) = ones(2,HEIGHT1);
%USED_SPACE
drawnow;
axes(H_A);

% shape definitions

a_shape = [-1 0 1 2; 0 0 0 0]';		% rectangle 1x4
b_shape = [-1 0 0 1; 0 0 1 0]';		% Z mirror
c_shape = [-1 0 0 1; 0 0 1 1]';		% Z normal
d_shape = [-1 0 0 1; 1 1 0 0]';		% truncated T
e_shape = [0 0 1 1; 0 1 0 1]';		% square 2x2
f_shape = [-1 0 1 1; 0 0 0 1]';		% L normal
g_shape = [-1 -1 0 1; 1 0 0 0]';	% L mirror

shape_matrix = [a_shape b_shape c_shape d_shape e_shape f_shape g_shape];
s = size(shape_matrix);
num_shapes = s(2) / 2;

color_matrix = [ 0 0 1; 0 1 0; 0 1 1; 1 0 0; 1 0 1; 1 1 0; 1 1 1];
home_pos = ones(4,1) * [(WIDTH1/2) (HEIGHT1-2)];
POS = a_shape + home_pos;
for i=1:4, H(i) = patch(SQ_X+POS(i,1),SQ_Y+POS(i,2),'r'); end;

while ((sum(USED_SPACE(2:WIDTH1+1,HEIGHT1-4)) == 0) & (STOP_GAME == 0)),

	shape_index = ceil(num_shapes*rand(1)); 	% choose random shape
	current_color = color_matrix(shape_index,:);	% associated color
	POS = shape_matrix(:,(2*shape_index - 1):(2*shape_index)) + home_pos;	% set initial postion of shape
	CENTER_X = home_pos(1,1);			% set origin of rotation
	CENTER_Y = home_pos(1,2);
	for i = 1:4,
		set(H(i),'XData',SQ_X+POS(i,1), 'YData',SQ_Y+POS(i,2), ...
			'FaceColor', current_color, 'Visible', 'on');
	end;
	drawnow;
	STUCK = 0;
	t0 = clock;
	while ((STUCK == 0) & (STOP_GAME == 0)),
		drawnow;
		if etime(clock,t0) > TIME_DELAY,	% move piece down one unit
			%etime(clock,t0)
			for i=1:4,
				if USED_SPACE(POS(i,1),POS(i,2)-1) == 1, STUCK=1; end;
			end;

			if STUCK == 0,
				t0 = clock;
				POS(:,2) = POS(:,2) - 1;
				CENTER_Y = CENTER_Y - 1;
				for i = 1:4,
					set(H(i),'YData',SQ_Y+POS(i,2));
				end;
				drawnow;
				%etime(clock,t0)

			end
		end;
		drawnow;
	end
	for i=1:4, USED_SPACE(POS(i,1),POS(i,2)) = 1; end;	% retire piece
	for i = 1:4, set(H(i),'Visible', 'off'); end;		% hide active piece
	axes(H_B);
	for i=1:4, h_matb(POS(i,1),POS(i,2)) = patch(SQ_X+POS(i,1),SQ_Y+POS(i,2),current_color); end;

	i=2;
	while i < HEIGHT1,  				% check if any lines are complete
		if sum(USED_SPACE(2:WIDTH1+1,i)) == WIDTH1,
			score = score + 10*(LEVEL+1);	% add score for completing line
			set(H_TEXT_SCORE, 'string', ['Score: ' num2str(score)]);
			for k1=2:WIDTH1+1, delete(h_matb(k1,i)); end;
			for k2=(i+1):HEIGHT1,
				for k1=2:WIDTH1+1,
					if h_matb(k1,k2) ~= 0,
						%disp([k1 k2 h_matb(k1,k2)]);
						set(h_matb(k1,k2),'YData',k2-1+SQ_Y);
					end;
				end;
			end;
			h_matb(2:WIDTH1+1,i:(HEIGHT1-1)) = h_matb(2:WIDTH1+1,(i+1):HEIGHT1);
			h_matb(2:WIDTH1+1,HEIGHT1) = zeros(WIDTH1,1);
			USED_SPACE(2:WIDTH1+1,i:(HEIGHT1-1)) = USED_SPACE(2:WIDTH1+1,(i+1):HEIGHT1);
			USED_SPACE(2:WIDTH1+1,HEIGHT1) = zeros(WIDTH1,1);
		else i = i+1;
		end

	end
	axes(H_A);

end  % end while loop

for k2=2:HEIGHT1,		% delete objects and clear play-field
	for k1=2:WIDTH1+1,
		if h_matb(k1,k2) ~= 0,
			delete(h_matb(k1,k2));
		end;
	end;
end;
set(START_HANDLE, 'visible', 'on');
set(STOP_HANDLE, 'visible', 'off');
set(EXIT_HANDLE, 'visible', 'on');

end  	% end 'start' command (end game)

author: Sable

repo: mcbench-benchmarks


text:
% NOTE: DO NOT REMOVE THIS LINE XMAKEFILE_TOOL_CHAIN_CONFIGURATION
function toolChainConfiguration = arduinomega2560()
%ARDUINO Defines a tool chain configuration

% Copyright 2012 The MathWorks, Inc.

try
    hSetup = realtime.setup.PackageInfo;
    spPkg = hSetup.getSpPkgInfo('Arduino Mega 2560');
    arduinoRootDir = hSetup.getTpPkgRootDir('ARDUINO', spPkg);
catch me
    toolChainConfiguration = {};
    return;
end

% General
toolChainConfiguration.Configuration = 'Arduinomega2560';
toolChainConfiguration.Version = '2.0';
toolChainConfiguration.Description = 'MS Visual Studio';
toolChainConfiguration.Operational = true;
toolChainConfiguration.Decorator = 'linkfoundation.xmakefile.decorator.eclipseDecorator';
% Make
toolChainConfiguration.MakePath = @(src) fullfile(arduinoRootDir, 'hardware', 'tools', 'avr', 'utils', 'bin', 'make');
toolChainConfiguration.MakeFlags = '-f "[|||MW_XMK_GENERATED_FILE_NAME[R]|||]" [|||MW_XMK_ACTIVE_BUILD_ACTION_REF|||]';
toolChainConfiguration.MakeInclude = '';
% Compiler
toolChainConfiguration.CompilerPath = @(src) fullfile(arduinoRootDir, 'hardware', 'tools', 'avr', 'bin', 'avr-gcc');
toolChainConfiguration.CompilerFlags = '-c -x none';
toolChainConfiguration.SourceExtensions = '.c,.cpp';
toolChainConfiguration.HeaderExtensions = '.h';
toolChainConfiguration.ObjectExtension = '.o';
% Linker
toolChainConfiguration.LinkerPath = @(src) fullfile(arduinoRootDir, 'hardware', 'tools', 'avr', 'bin', 'avr-gcc');
toolChainConfiguration.LinkerFlags = '-o [|||MW_XMK_GENERATED_TARGET_REF|||]';
toolChainConfiguration.LibraryExtensions = '.lib,.a';
toolChainConfiguration.TargetExtension = '.elf';
toolChainConfiguration.TargetNamePrefix = '';
toolChainConfiguration.TargetNamePostfix = '';
% Archiver
toolChainConfiguration.ArchiverPath = @(src) fullfile(arduinoRootDir, 'hardware', 'tools', 'avr', 'bin', 'avr-ar');
toolChainConfiguration.ArchiverFlags = '-crs $(TARGET_FILE)';
toolChainConfiguration.ArchiveExtension = '.a';
toolChainConfiguration.ArchiveNamePrefix = 'lib';
toolChainConfiguration.ArchiveNamePostfix = '';
% Pre-build
toolChainConfiguration.PrebuildEnable = false;
toolChainConfiguration.PrebuildToolPath = '';
toolChainConfiguration.PrebuildFlags = '';
% Post-build
toolChainConfiguration.PostbuildEnable = true;
toolChainConfiguration.PostbuildToolPath = @(src) fullfile(arduinoRootDir, 'hardware', 'tools', 'avr', 'bin', 'avr-objcopy');
toolChainConfiguration.PostbuildFlags = '-O ihex -R .eeprom $(TARGET)  [|||MW_XMK_MODEL_NAME|||].hex';
% Execute
toolChainConfiguration.ExecuteDefault = false;
toolChainConfiguration.ExecuteToolPath = '';
toolChainConfiguration.ExecuteFlags = '';
% Directories
toolChainConfiguration.DerivedPath = '[|||MW_XMK_SOURCE_PATH_REF|||]';
toolChainConfiguration.OutputPath = '';
% Custom
toolChainConfiguration.Custom1 = '';
toolChainConfiguration.Custom2 = '';
toolChainConfiguration.Custom3 = '';
toolChainConfiguration.Custom4 = '';
toolChainConfiguration.Custom5 = '';
end


text:
function a = elastix_affine_struct2matrix(tf)
% ELASTIX_AFFINE_STRUCT2MATRIX  Convert any type of affine transform from
% elastix struct format to homogeneous coordinates matrix format.
%
% A = ELASTIX_AFFINE_STRUCT2MATRIX(TF)
%
%   TF is a struct in elastix format of a 2D transform. The transform may
%   be of any affine type:
%
%     'AffineTransform'
%     'SimilarityTransform'
%     'EulerTransform'
%     'TranslationTransform'
%
%   It may also have an arbitrary TF.CenterOfRotationPoint, it doesn't need
%   to be 0. Subsequent transforms pointed at from
%   TF.InitialTransformParametersFileName are ignored.
%
%   A is an affine transform matrix in homegeneous coordinates, referred to
%   the centre of coordinates = 0. The transform is computed as
%
%     [y 1] = [x 1] * A
%
%   where A = [a11 a12 0]
%             [a21 a22 0]
%             [t1  t2  1]
%
%   [t1 t2] is the translation vector, whereas
%
%   for a general affine transform: [a11 a12]
%                                   [a21 a22]
%
%   for a similarity transform:     [ cos(theta) sin(theta)] * s
%                                   [-sin(theta) cos(theta)]
%
%   for a rigid transform:          [ cos(theta) sin(theta)]
%                                   [-sin(theta) cos(theta)]
%
%   for a translation transform:    [1 0]
%                                   [0 1]
%
%   If TF is a vector of structures, A is an array where each A(:, :, I)
%   corresponds to TP(I).
%
% See also: elastix_affine_matrix2struct.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2015 University of Oxford
% Version: 0.2.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
narginchk(1, 1);
nargoutchk(0, 1);

% number of input transforms
N = length(tf);

% init output
a = zeros(3, 3, N);

% loop each transform
for I = 1:N

    a(:, :, I) = one_struct2matrix(tf(I));

end

end

%% auxiliary function
% one_struct2matrix: convert a single transform to matrix form
function a = one_struct2matrix(tf)

switch (tf.Transform)

    case 'AffineTransform'

        % transformation nomenclature
        a11 =   tf.TransformParameters(1);
        a21 =   tf.TransformParameters(2);
        a12 =   tf.TransformParameters(3);
        a22 =   tf.TransformParameters(4);
        t =     tf.TransformParameters(5:6);

        % matrix block that doesn't contain the translation
        A = [a11 a12; a21 a22];

    case 'SimilarityTransform'

        % transformation nomenclature
        s =     tf.TransformParameters(1);
        theta = tf.TransformParameters(2);
        t =     tf.TransformParameters(3:4);

        % matrix block that doesn't contain the translation
        A = [cos(theta) sin(theta);...
            -sin(theta) cos(theta)] * s;

    case 'EulerTransform'

        % transformation nomenclature
        theta = tf.TransformParameters(1);
        t =     tf.TransformParameters(2:3);

        % matrix block that doesn't contain the translation
        A = [cos(theta) sin(theta);...
            -sin(theta) cos(theta)];

    case 'TranslationTransform'

        % transformation nomenclature
        t =     tf.TransformParameters(1:2);

        % matrix block that doesn't contain the translation
        A = [1 0;...
            0 1];

    otherwise

        error('Transform not implemented')

end

% centre of rotation
c = tf.CenterOfRotationPoint;

% center transform on origin
t = t + c * (eye(2) - A);

% create affine transform matrix
a = [A [0;0]; t 1];

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function s=phaseran(x,c)
%Syntax: s=phaseran(x,c)
%_______________________
%
% Makes c phase randomized surrogates of a time series x.
%
% s is the phase randomized time series.
% x is the original time series.
% c is the number of surrogates.
%
%
% References:
%
% Theiler J, Galdrikian B, Longtin A, Eubank S, Farmer D J (1992): Using
% Surrogate Data to Detect Nonlinearity in Time Series. In Nonlinear Modeling
% and Forecasting, eds. Casdagli M & Eubank S. 163-188. Addison-Wesley
%
% Theiler J, Eubank S,Galdrikian B, Longtin A,  Farmer D J (1992): Testing
% for nonlinearity in time series: the method of surrogate data. Physica D
% 58: 77-94
%
%
% Alexandros Leontitsis
% Department of Education
% University of Ioannina
% 45110 - Dourouti
% Ioannina
% Greece
%
% University e-mail: me00743@cc.uoi.gr
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% 12 Apr 2002

if nargin<1 | isempty(x)==1
   error('You should provide a time series.');
else
   % x must be a vector
   if min(size(x))>1
      error('Invalid time series.');
   end
   x=x(:);
   % N is the time series length
   N=length(x);
end

if nargin<2 | isempty(c)==1
   c=1;
else
   % c must be scalar
   if sum(size(c))>2
      error('c must be scalar.');
   end
   % c must be greater or equal than 1
   if c<1
      error('c must be greater or equal than 1.');
   end
end

% FFT on x
y=fft(x);
% Magnitudes
m=abs(y);
% Angles
p=angle(y);
% The imaginary unit
i=sqrt(-1);
% Half of the data points
h=floor(N/2);

for j=1:c
   % Randomized phases
   if rem(N,2)==0
      p1=rand(h-1,1)*2*pi;
      p(2:N)=[p1' p(h+1) -flipud(p1)'];
      % Adjust the magnitudes
      m=[m(1:h+1);flipud(m(2:h))];
   else
      p1=rand(h,1)*2*pi;
      p(2:N)=[p1 -flipud(p1)];
   end
   % Back to the complex numbers
   s(:,j)=m.*exp(i*p);
   % Back to the time series (phase randomized surrogates)
   s(:,j)=real(ifft(s(:,j)));
end


text:
classdef LossSmoothL1 < dagnn.Loss
%LossSmoothL1  Smooth L1 loss
%  `LossSmoothL1.forward({x, x0, w})` computes the smooth L1 distance
%  between `x` and `x0`, weighting the elements by `w`.
%
%  Here the smooth L1 loss between two vectors is defined as:
%
%     Loss = sum_i f(x_i - x0_i) w_i.
%
%  where f is the function (following the Faster R-CNN definition):
%
%              { 0.5 * sigma^2 * delta^2,         if |delta| < 1 / sigma^2,
%   f(delta) = {
%              { |delta| - 0.5 / sigma^2,         otherwise.
%
%  In practice, `x` and `x0` can pack multiple instances as 1 x 1 x C
%  x N arrays (or simply C x N arrays).

  properties
    sigma = 1.
  end

  methods
    function outputs = forward(obj, inputs, params)
      sigma2 = obj.sigma^2 ;
      delta = inputs{1} - inputs{2} ;
      absDelta = abs(delta) ;

      linearRegion = (absDelta > 1. / sigma2) ;
      absDelta(linearRegion) = absDelta(linearRegion) - 0.5/sigma2 ;
      absDelta(~linearRegion) = 0.5 * sigma2 * absDelta(~linearRegion).^2 ;

      % Mutliply by instance weights and sum.
      outputs{1} = inputs{3}(:)' * absDelta(:) ;

      % Accumulate loss statistics.
      if obj.ignoreAverage, return; end;
      n = obj.numAveraged ;
      m = n + gather(sum(inputs{3}(:))) + 1e-9 ;
      obj.average = (n * obj.average + gather(outputs{1})) / m ;
      obj.numAveraged = m ;
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
    % Function derivative:
    %
    %          { sigma^2 * x,             if |x| < 1 / sigma^2,
    %  f'(x) = {
    %          { sign(x),                 otherwise.

      sigma2 = obj.sigma^2 ;
      delta = inputs{1} - inputs{2} ;
      absDelta = abs(delta) ;

      linearRegion = (absDelta > 1. / sigma2) ;
      delta(linearRegion) = sign(delta(linearRegion));
      delta(~linearRegion) = sigma2 * delta(~linearRegion) ;

      derInputs = {inputs{3} .* delta .* derOutputs{1}, [], []} ;
      derParams = {} ;
    end

    function obj = LossSmoothL1(varargin)
      obj.load(varargin) ;
      obj.loss = 'smoothl1';
    end
  end
end


text:
function save_model_path = mask_train(conf, imdb_train, roidb_train, varargin)
% save_model_path = fast_rcnn_train(conf, imdb_train, roidb_train, varargin)
% --------------------------------------------------------
% Fast R-CNN
% Reimplementation based on Python Fast R-CNN (https://github.com/rbgirshick/fast-rcnn)
% Copyright (c) 2015, Shaoqing Ren
% Licensed under The MIT License [see LICENSE for details]
% --------------------------------------------------------

%% inputs
    ip = inputParser;
    ip.addRequired('conf',                              @isstruct);
    ip.addRequired('imdb_train',                        @iscell);
    ip.addRequired('roidb_train',                       @iscell);
    ip.addParamValue('do_val',          false,          @isscalar);
    ip.addParamValue('imdb_val',        struct(),       @isstruct);
    ip.addParamValue('roidb_val',       struct(),       @isstruct);
    ip.addParamValue('val_iters',       100,            @isscalar);
    ip.addParamValue('val_interval',    400,           @isscalar);
    ip.addParamValue('snapshot_interval',...
                                        2000,          @isscalar);
    ip.addParamValue('solver_def_file', fullfile(pwd, 'models', 'mask_rcnn_prototxts/resnet_50layers_conv', 'solver_30k40k.prototxt'), ...
                                                        @isstr);
    ip.addParamValue('net_file',        fullfile(pwd, 'models', 'mask_rcnn_prototxts/resnet_50layers_conv', 'solver_30k40k.prototxt'), ...
                                                        @isstr);
    ip.addParamValue('cache_name',      'mask_RES50_VOC2007', ...
                                                        @isstr);

    ip.parse(conf, imdb_train, roidb_train, varargin{:});
    opts = ip.Results;

%% try to find trained model
    imdbs_name = cell2mat(cellfun(@(x) x.name, imdb_train, 'UniformOutput', false));
    cache_dir = fullfile(pwd, 'output', 'mask_rcnn_cachedir', opts.cache_name, imdbs_name);
    save_model_path = fullfile(cache_dir, 'final')
    if exist(save_model_path, 'file')
        fprintf('Found existing model\n');
        return;
    end
%% init
    % init caffe solver
    mkdir_if_missing(cache_dir);
    caffe_log_file_base = fullfile(cache_dir, 'caffe_log');
    caffe.init_log(caffe_log_file_base);
    caffe_solver = caffe.Solver(opts.solver_def_file);
    caffe_solver.net.copy_from(opts.net_file);
    % init log
    timestamp = datestr(datevec(now()), 'yyyymmdd_HHMMSS');
    mkdir_if_missing(fullfile(cache_dir, 'log'));
    log_file = fullfile(cache_dir, 'log', ['train_', timestamp, '.txt']);
    diary(log_file);

    % set random seed
    prev_rng = seed_rand(conf.rng_seed);
    caffe.set_random_seed(conf.rng_seed);

    % set gpu/cpu
    if conf.use_gpu
        caffe.set_mode_gpu();
    else
        caffe.set_mode_cpu();
    end
    %disp(conf);
    %disp(opts);
   fprintf('Preparing data\n');
%% making tran/val data
    [image_roidb_train, bbox_means, bbox_stds]...
                            = mask_prepare_image_roidb_voc(conf, opts.imdb_train, opts.roidb_train);

  if opts.do_val
        [image_roidb_val]...
                                = mask_prepare_image_roidb_voc(conf, opts.imdb_val, opts.roidb_val, bbox_means, bbox_stds);
        % fix validation data
        shuffled_inds_val = generate_random_minibatch([], image_roidb_val, conf.ims_per_batch);
        shuffled_inds_val = shuffled_inds_val(randperm(length(shuffled_inds_val), opts.val_iters));
    end
   fprintf('Done\n');
%%  try to train/val with images which have maximum size potentially, to validate whether the gpu memory is enough
    num_classes = size(image_roidb_train(1).overlap, 2);
    fprintf('Checking GPU Memory\n');
    check_gpu_memory(conf, caffe_solver, num_classes, opts.do_val);
    fprintf('Done\n');
%% training
    shuffled_inds = [];
    train_results = [];
    val_results = [];
    iter_ = caffe_solver.iter();
    max_iter = caffe_solver.max_iter();

    while (iter_ < max_iter)
        caffe_solver.net.set_phase('train');
        % generate minibatch training data
        [shuffled_inds, sub_db_inds] = generate_random_minibatch(shuffled_inds, image_roidb_train, conf.ims_per_batch);
        [im_blob, rois_blob, labels_blob, bbox_targets_blob, bbox_loss_weights_blob, mask_targets_blob] = ...
            mask_get_minibatch(conf, image_roidb_train(sub_db_inds));

        net_inputs = {im_blob, rois_blob, labels_blob, bbox_targets_blob, bbox_loss_weights_blob, mask_targets_blob};

        caffe_solver.net.reshape_as_input(net_inputs);

        % one iter SGD update
        caffe_solver.net.set_input_data(net_inputs);
        caffe_solver.step(1);
        rst = caffe_solver.net.get_output();
        train_results = parse_rst(train_results, rst);
        % do valdiation per val_interval iterations
        if ~mod(iter_, opts.val_interval)
            if opts.do_val
                caffe_solver.net.set_phase('test');
                for i = 1:length(shuffled_inds_val)
                    sub_db_inds = shuffled_inds_val{i};
                    [im_blob, rois_blob, labels_blob, bbox_targets_blob, bbox_loss_weights_blob, mask_targets_blob] = ...
                        mask_get_minibatch(conf, image_roidb_val(sub_db_inds));

                    % Reshape net's input blobs
                    net_inputs = {im_blob, rois_blob, labels_blob, bbox_targets_blob, bbox_loss_weights_blob, mask_targets_blob};
                    caffe_solver.net.reshape_as_input(net_inputs);

                    caffe_solver.net.forward(net_inputs);
                    rst = caffe_solver.net.get_output();
                    val_results = parse_rst(val_results, rst);
                end
            end
            show_state(iter_, train_results, val_results);
            train_results = [];
            val_results = [];
            diary; diary; % flush diary
        end

        % snapshot
        if ~mod(iter_, opts.snapshot_interval)
            snapshot(caffe_solver, bbox_means, bbox_stds, cache_dir, sprintf('iter_%d', iter_));
        end
        iter_ = caffe_solver.iter();
    end

    % final snapshot
    snapshot(caffe_solver, bbox_means, bbox_stds, cache_dir, sprintf('iter_%d', iter_));
    save_model_path = snapshot(caffe_solver, bbox_means, bbox_stds, cache_dir, 'final');

    diary off;
    caffe.reset_all();
    rng(prev_rng);
end

function [shuffled_inds, sub_inds] = generate_random_minibatch(shuffled_inds, image_roidb_train, ims_per_batch)

    % shuffle training data per batch
    if isempty(shuffled_inds)
        % make sure each minibatch, only has horizontal images or vertical
        % images, to save gpu memory

        hori_image_inds = arrayfun(@(x) x.im_size(2) >= x.im_size(1), image_roidb_train, 'UniformOutput', true);
        vert_image_inds = ~hori_image_inds;
        hori_image_inds = find(hori_image_inds);
        vert_image_inds = find(vert_image_inds);

        % random perm
        lim = floor(length(hori_image_inds) / ims_per_batch) * ims_per_batch;
        hori_image_inds = hori_image_inds(randperm(length(hori_image_inds), lim));
        lim = floor(length(vert_image_inds) / ims_per_batch) * ims_per_batch;
        vert_image_inds = vert_image_inds(randperm(length(vert_image_inds), lim));

        % combine sample for each ims_per_batch
        hori_image_inds = reshape(hori_image_inds, ims_per_batch, []);
        vert_image_inds = reshape(vert_image_inds, ims_per_batch, []);
        shuffled_inds = [hori_image_inds, vert_image_inds];
        shuffled_inds = shuffled_inds(:, randperm(size(shuffled_inds, 2)));

        shuffled_inds = num2cell(shuffled_inds, 1);
    end

    if nargout > 1
        % generate minibatch training data
        sub_inds = shuffled_inds{1};
        assert(length(sub_inds) == ims_per_batch);
        shuffled_inds(1) = [];
    end
end

function check_gpu_memory(conf, caffe_solver, num_classes, do_val)
%%  try to train/val with images which have maximum size potentially, to validate whether the gpu memory is enough
    % generate pseudo training data with max size
    im_blob = single(zeros(max(conf.scales), conf.max_size, 3, conf.ims_per_batch));
    rois_blob = single(repmat([0; 0; 0; max(conf.scales)-1; conf.max_size-1], 1, conf.batch_size));
    rois_blob = permute(rois_blob, [3, 4, 1, 2]);
    labels_blob = single(ones(conf.batch_size, 1));
    labels_blob = permute(labels_blob, [3, 4, 2, 1]);
    bbox_targets_blob = zeros(4 * (num_classes+1), conf.batch_size, 'single');
    bbox_targets_blob = single(permute(bbox_targets_blob, [3, 4, 1, 2]));
    bbox_loss_weights_blob = bbox_targets_blob;

    mask_targets_blob = zeros(conf.batch_size, num_classes+1, 28, 28, 'single');
    mask_targets_blob(:) = -1;
    mask_targets_blob = single(permute(mask_targets_blob, [4, 3, 2, 1]));

    net_inputs = {im_blob, rois_blob, labels_blob, bbox_targets_blob, bbox_loss_weights_blob, mask_targets_blob};
    % Reshape net's input blobs
    caffe_solver.net.reshape_as_input(net_inputs);
    % one iter SGD update
    caffe_solver.net.set_input_data(net_inputs);
    caffe_solver.step(1);
    if do_val
        % use the same net with train to save memory
        caffe_solver.net.set_phase('test');
        caffe_solver.net.forward(net_inputs);
        caffe_solver.net.set_phase('train');
    end
end

function model_path = snapshot(caffe_solver, bbox_means, bbox_stds, cache_dir, file_name)
    bbox_stds_flatten = reshape(bbox_stds', [], 1);
    bbox_means_flatten = reshape(bbox_means', [], 1);

    % merge bbox_means, bbox_stds into the model
    bbox_pred_layer_name = 'bbox_pred';
    weights = caffe_solver.net.params(bbox_pred_layer_name, 1).get_data();
    biase = caffe_solver.net.params(bbox_pred_layer_name, 2).get_data();
    weights_back = weights;
    biase_back = biase;

    weights = ...
        bsxfun(@times, weights, bbox_stds_flatten'); % weights = weights * stds;
    biase = ...
        biase .* bbox_stds_flatten + bbox_means_flatten; % bias = bias * stds + means;

    caffe_solver.net.set_params_data(bbox_pred_layer_name, 1, weights);
    caffe_solver.net.set_params_data(bbox_pred_layer_name, 2, biase);

    model_path = fullfile(cache_dir, file_name);
    caffe_solver.net.save(model_path);
    fprintf('Saved as %s\n', model_path);

    % restore net to original state
    caffe_solver.net.set_params_data(bbox_pred_layer_name, 1, weights_back);
    caffe_solver.net.set_params_data(bbox_pred_layer_name, 2, biase_back);
end

function show_state(iter, train_results, val_results)
    fprintf('\n------------------------- Iteration %d -------------------------\n', iter);
    fprintf('Training : error %.3g, loss (cls %.3g, reg %.3g, mask %.3g)\n', ...
        1 - mean(train_results.accuarcy.data), ...
        mean(train_results.loss_cls.data), ...
        mean(train_results.loss_bbox.data), ...
        mean(train_results.loss_mask.data));
    if exist('val_results', 'var') && ~isempty(val_results)
        fprintf('Testing  : error %.3g, loss (cls %.3g, reg %.3g, mask %.3g)\n', ...
            1 - mean(val_results.accuarcy.data), ...
            mean(val_results.loss_cls.data), ...
            mean(val_results.loss_bbox.data), ...
            mean(val_results.loss_mask.data));
    end
end

author: jasjeetIM

repo: Mask-RCNN

save_path: github-repos/MATLAB/jasjeetIM-Mask-RCNN


text:
function y = tanh(x)
exp_term = exp(x);
exp_term2 = exp(-x);
y = (exp_term - exp_term2) ./ (exp_term + exp_term2);

end

author: singaxiong

repo: SignalGraph


text:
function [fusion,w0] = qfuser_v5(w,scores,wfuse)

if nargin==0
    test_this();
    return;
end

% block 1
f1 = linear_fuser([],scores.scores);
w1 = wfuse;
[whead,wtail] = splitvec_fh(length(w1));
f1 = f1(whead);

% block 2
modelQ = scores.modelQ;
[q,n1] = size(modelQ);
modelQ = [modelQ;ones(1,n1)];
segQ = scores.segQ;
[q2,n2] = size(segQ);
segQ = [segQ;ones(1,n2)];
assert(q==q2);
q = q + 1;

wq = q*(q+1)/2;
f2 = AWB_fh(modelQ',segQ,tril_to_symm_fh(q,wtail));
w2 = zeros(wq,1);

% assemble
fusion = sum_of_functions(w,[1,1],f1,f2);
w0 = [w1;w2];

end

function test_this()

m = 5;
k = 2;
n1 = 4;
n2 = 5;

scores.sindx = [1,2,3];
scores.qindx = [4,5];

scores.scores = randn(m,n1*n2);
scores.modelQ = randn(k,n1);
scores.segQ = randn(k,n2);

wfuse = [1,2,3,4]';

[fusion,w0] = qfuser_v4([],scores,wfuse);

%test_MV2DF(fusion,w0);

[fusion(w0),linear_fuser(wfuse,scores.scores(scores.sindx,:))]

%fusion(w0)

end

author: nesl


text:
% Copyright (C) 1993-2014, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

%%begin
% In the field of robotics there are many possible ways of representing
% positions and orientations, but the homogeneous transformation is well
% matched to MATLABs powerful tools for matrix manipulation.
%
% Homogeneous transformations describe the relationships between Cartesian
% coordinate frames in terms of translation and orientation.

%  A pure translation of 0.5m in the X direction is represented by

transl(0.5, 0.0, 0.0)

% a rotation of 90degrees about the Y axis by

troty(pi/2)

% and a rotation of -90degrees about the Z axis by

trotz(-pi/2)

%  these may be concatenated by multiplication

t = transl(0.5, 0.0, 0.0) * troty(pi/2) * trotz(-pi/2)

%
% If this transformation represented the origin of a new coordinate frame with respect
% to the world frame origin (0, 0, 0), that new origin would be given by

t * [0 0 0 1]'

% the orientation of the new coordinate frame may be expressed in terms of
% Euler angles

tr2eul(t)

% or roll/pitch/yaw angles

tr2rpy(t)

% It is important to note that tranform multiplication is in general not
% commutative as shown by the following example

trotx(pi/2) * trotz(-pi/8)
trotz(-pi/8) * trotx(pi/2)

author: Allopart

repo: rbpf-gmapping


text:
%RELMTEST
%script to test two earthqauake rate hypotheses using earthquake data
%
region='Northern California Aftershock model';
nsim=100;
mt=2;
 park1=vRatesH;
 park2=vRatesN;
%load H1;
%load H2;
%park1=H1; %Hypothesis with variable b-values
%park2=H2; %Null hypothesis, uniform b-value
clear test null;
%     xmin(i)=park1(j,1);
%     xmax(i)=park1(j,2);
%     ymin(i)=park1(j,3);
%     ymax(i)=park1(j,4);
%     zmin(i)=park1(j,5);
%     zmax(i)=park1(j,6);
%     magmin(i)=park1(j,7);
%     magmax(i)=park1(j,8);
%     lamda1(i)=park1(j,9);
%     weight(i)=park1(j,10);
nquake=park1(:,11);
[m,n] = size(park1);
magmin=park1(:,7);
lamda1=park1(:,9);
lamda2=park2(:,9);
weight1=park1(:,10);
weight2=park2(:,10);
weight=weight1.*weight2.*(magmin>mt);
%
% Remove rows of matrix for which weight is zero
%
j=0;
for i = 1:m
    if weight(i)>0
        j=j+1;
        w(j)=weight(i);
        nq(j)=nquake(i);
        lam1(j)=lamda1(i);
        lam2(j)=lamda2(i);
        mmin(j)=magmin(i);
    end
end
nq=w.*nq;
Nquake=sum(nq);
lam1=w.*lam1;
lam2=w.*lam2;
clear park1 park2 lamda1 lamda2 nquake magmin weight weight1 weight2;
%
%make a weighted magnitude-frequency plot
%
mf=[mmin;nq;lam1;lam2]';
mfsort=sortrows(mf);
mag=mfsort(:,1);

Fobs=flip(cumsum(flip(mfsort(:,2))));
Fth1=flip(cumsum(flip(mfsort(:,3))));
Fth2=flip(cumsum(flip(mfsort(:,4))));
figure%(1)
semilogy(mag,Fobs,'r',mag,Fth1,'g',mag,Fth2,'b');
grid;
axis([3,8,.0001,100]);

%
%    Evaluate whether total number of quakes is consistent with H1
%
%
Nhat=sum(lam1)
peq=poisspdf(Nquake, Nhat); % probability of exactly Nquake
Ple=poisscdf(Nquake, Nhat); % probability of less than or equal to Nquake
Pless=Ple-peq;              % probability of less than Nquake
Pmore=1-Ple                 % probability of more than Nquake
P1_equal=peq;
P1_less=Pless;
P1_more=Pmore;
Nhat1=Nhat;
lamcum1=cumsum(lam1)/Nhat1;

%   Evaluate whether total number of quakes is consistent with H2

Nhat=sum(lam2)
peq=poisspdf(Nquake, Nhat); % probability of exactly Nquake
Ple=poisscdf(Nquake, Nhat); % probability of less than or equal to Nquake
Pless=Ple-peq;              % probability of less than Nquake
Pmore=1-Ple                 % probability of more than Nquake
P2_equal=peq;
P2_less=Pless;
P2_more=Pmore;
Nhat2=Nhat;
lamcum2=cumsum(lam2)/Nhat2;
%
%   simulate catalogs according to H1,
%   and evaluate likelihood scores of nsquake1 and real catalog using lamda1 and lamda2
%
nsquake=simulate(Nquake,lam1, nsim);
[LLR1, rank11,rank12] = Rtest(lam1, lam2, nq, nsquake, w);
%
%   simulate catalogs according to H2,
%   and evaluate likelihood scores of nsquake1 and real catalog using lamda1 and lamda2
%
nsquake=simulate(Nquake,lam2, nsim);
[LLR2, rank21,rank22] = Rtest(lam1, lam2, nq, nsquake, w);
%
%Plot cumulative likelihood scores for two hypotheses
%
alpha = sum(LLR2>0)/nsim
beta = sum(LLR1<0)/nsim
index=[1:nsim]/nsim;
x=[0,0];y=[0,1];
figure_w_normalized_uicontrolunits(2);
plot(LLR1,index,'g',LLR2,index,'r',x,y,'b')';
xlabel('Likelihood ratio (Variable b/Constant b)')';
ylabel('Fraction of cases');
title('Green assumes variable-b hypothesis; Red assumes constant=b hypothesis');

region, mt,Nquake, Nhat1,Nhat2,P1_less,P1_more,P2_less,P2_more, alpha, beta, rank11,rank12, rank21, rank22
%[rank1,rank2]
%[Nhat1,Nhat2]
%[P1_less, P2_less]
%[P1_more, P2_more]

author: CelsoReyes


text:
%  演示调用 LU 分解的紧凑方式
clear;
A=[1 2 3; 2 5 2; 3 1 5];
b=[14;18;20];

% LU三角分解法
[L, U, x] = my_lu(A, b)

% 列主元LU三角分解法
%[L, U, x] = my_lu_with_column_pivoting(A, b)

% 直接调用matlab左除运算符
x_ = A\b

repo: NumericalAnalysis


text:
% Small technical example.
% Shows how to retrieve the details on a quincunx grid, after a multi-
% resolution decomposition has been performed.
%
disp('Small technical example.');
disp('Shows how to retrieve the details on a quincunx grid, after a multi-');
disp('resolution decomposition has been performed.');
disp('FOR MORE INFORMATION:  help retrieveQ1001');
disp(' ');
disp('See also the report http://repository.cwi.nl:8888/cwi_repository/docs/IV/04/04178D.pdf');
disp('Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>');
disp(' (C) 1998-2006 Stichting CWI, Amsterdam, The Netherlands');
disp(' ');
%---PARAMETERS-----------------------------------------------------------------
% How to execute, set parameters
N = 6;                   %  maximum level (even number) in lifting scheme
filtername = 'Neville4';
%
%---INSERT YOUR IMAGE HERE-----------------------------------------------------
if exist('imread','file') == 2
  Orig = double(imread('zenithgray.TIF','tiff'));
else
  load zenithgray; Orig = zenithgray; clear zenithgray;
end
%
disp([' Dimensions of original      ' int2str( size(Orig) )]);
%---DECOMPOSITION--------------------------------------------------------------
disp([' Filter type is ' filtername]);
[C,S] = QLiftDec2(Orig,N,filtername);
%
%---RETRIEVE DETAIL AT ODD LEVEL-----------------------------------------------
level = N-1;
disp([' Level ' int2str(level)]);
[F10, F01] = retrieveQ1001(level, 'd', C, S);
%
disp('The dimensions of the detail function read as follows');
disp([int2str(size(F10)) ' and ' int2str(size(F01))]);
%


text:
function len = edgeLength(varargin)
%EDGELENGTH Return length of an edge
%
%   L = edgeLength(EDGE);
%   Returns the length of an edge, with parametric representation:
%   [x1 y1 x2 y2].
%
%   The function also works for several edges, in this case input is a
%   N-by-4 array, containing parametric representation of each edge, and
%   output is a N-by-1 array containing length of each edge.
%
%   See also:
%   edges2d, edgeAngle
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 19/02/2004
%

%   HISTORY
%   15/04/2005 changes definition for edge, uses [x1 y1 x2 y2] instead of
%       [x0 y0 dx dy].

%   TODO : specify norm (euclidian, taxi, ...).

if nargin == 1
    % input is an edge [X1 Y1 X2 Y2]
    edge = varargin{1};
    len = hypot(edge(:,3)-edge(:,1), edge(:,4)-edge(:,2));

elseif nargin == 2
    % input are two points [X1 Y1] and [X2 Y2]
    p1 = varargin{1};
    p2 = varargin{2};
    len = hypot(p2(:,1)-p1(:,1), p2(:,2)-p1(:,2));

end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips