text:
function simplenn_caffe_deploy(net, caffeModelBaseName, varargin)
% SIMPLENN_CAFFE_DEPLOY export a simplenn network to Caffe model
%   SIMPLENN_CAFFE_DEPLOY(NET, CAFFE_BASE_MODELNAME)
%   Export a simplenn network NET to a Caffe model.
%   The caffe model is stored in the following files:
%
%     [CAFFE_BASE_MODELNAME '.prototxt']   - Network definition file
%     [CAFFE_BASE_MODELNAME '.caffemodel'] - Binary Caffe model
%     [CAFFE_BASE_MODELNAME '_mean_image.binaryproto'] (optional) -
%         The average image (if set in net.normalization.averageImage)
%
%   Compiled MatCaffe (usually located in `<caffe_dir>/matlab`, built
%   with the `matcaffe` target) must be in path.
%
%   Only a limited subset of layers is currently supported and those are:
%
%     Conv, ReLU, Pool, LRN, SoftMax, SoftMaxLogLoss, DropOut
%
%   Please note that thanks to different implementations, the outputs of
%   simplenn and Caffe models are not neccessarily identical.
%
%   SIMPLENN_CAFFE_DEPLOY(NET, CAFFE_BASE_MODELNAME, 'OPT', VAL, ...)
%   takes the following options:
%
%   `removeDropout`:: `true`
%      When true, do not deploy dropout layers.
%
%   `replaceSoftMaxLoss`:: `true`
%      Replace SoftMax log loss with SoftMax.
%
%   `doTest`:: `true`
%      Compare the caffe model to the simplenn model using
%      `simplenn_caffe_compare`.
%
%   `testData`:: Random
%      Perform the test on the given data.
%
%   `inputBlobName`:: 'data'
%      Name of the input data blob in the final protobuf.
%
%   `labelBlobName`:: 'label'
%      Name of the input label blob in the final protobuf.
%
%   `outputBlobName`:: 'prob'
%      Name of the output blob in the resulting protobuf.
%
%   `silent`:: false
%      When true, suppresses all output to stdout.
%
%  See Also: simplenn_caffe_compare

% Copyright (C) 2015-16 Zohar Bar-Yehuda, 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).

opts.inputBlobName = 'data';
opts.outputBlobName = 'prob';
opts.labelBlobName = 'label';
opts.removeDropout = true;
opts.replaceSoftMaxLoss = true;
opts.doTest = true;
opts.testData = [];
opts.silent = false;
opts = vl_argparse(opts, varargin);
if ~exist('caffe.Net', 'class'), error('MatCaffe not in path.'); end

info = @(varargin) fprintf(1, varargin{:});
if opts.silent, info = @(varargin) []; end;

info('Exporting simplenn model to caffe model %s\n', caffeModelBaseName);
[modelDir, name] = fileparts(caffeModelBaseName);
[~,~,~] = mkdir(modelDir);

% -------------------------------------------------------------------------
%                                                          Tidy the network
% -------------------------------------------------------------------------
net = vl_simplenn_tidy(net);
net = vl_simplenn_move(net, 'cpu');
% Remove dropout layers
if opts.removeDropout
  net.layers(cellfun(@(l) strcmp(l.type, 'dropout'), net.layers)) = [];
end

if opts.replaceSoftMaxLoss
  % If last layer is softmax loss, replace it with softmax
  ll = net.layers{end};
  if strcmp(ll.type, 'softmaxloss') || ...
      (strcmp(ll.type, 'loss') && strcmp(ll.loss, 'softmaxlog'))
    net.layers{end}.type = 'softmax';
  elseif isequal(net.layers{end}.type, 'loss')
    error('Unsupported loss function: %s', net.layers{end}.loss);
  end
end

for idx = 1:numel(net.layers)
  % Add missing layer names
  if ~isfield(net.layers{idx}, 'name')
    net.layers{idx}.name = sprintf('layer%d', idx);
  end
end

avImage = [];
if isfield(net.meta, 'normalization') && ...
    isfield(net.meta.normalization, 'imageSize')
  imSize = net.meta.normalization.imageSize;
  if isfield(net.meta.normalization, 'averageImage')
    avImage = net.meta.normalization.averageImage;
    if numel(avImage) == imSize(3)
      avImage = reshape(avImage, 1, 1, imSize(3));
    end
  end
else
  error('Missing image size. Please set `net.normalization.imageSize`.');
end

% -------------------------------------------------------------------------
%                                                           Export prototxt
% -------------------------------------------------------------------------
prototxtFilename = [caffeModelBaseName '.prototxt'];
fid = fopen(prototxtFilename, 'w');

fprintf(fid, 'name: "%s"\n\n', name); % Network name

% Export input dimensions
fprintf(fid, 'input: "data"\n');
fprintf(fid, 'input_dim: 1\n');
fprintf(fid, 'input_dim: %d\n', imSize(3));
fprintf(fid, 'input_dim: %d\n', imSize(1));
fprintf(fid, 'input_dim: %d\n\n', imSize(2));

dummyData = zeros(imSize, 'single'); % Keep track of data size at each layer;

isFullyConnected = false(size(net.layers));
for idx = 1:numel(net.layers)
  % write layers
  fprintf(fid,'layer {\n');
  fprintf(fid,'  name: "%s"\n', net.layers{idx}.name); % Layer name
  layerInputSize = size(dummyData);
  if numel(layerInputSize) == 2
    layerInputSize(3) = 1;
  end
  switch net.layers{idx}.type
    case 'conv'
      filtH = size(net.layers{idx}.weights{1},1);
      filtW = size(net.layers{idx}.weights{1},2);
      if filtH < layerInputSize(1) || filtW < layerInputSize(2)
        % Convolution layer
        fprintf(fid, '  type: "Convolution"\n');
        write_connection(fid, net.layers, idx);
        fprintf(fid, '  convolution_param {\n');
        write_kernel(fid, [filtH, filtW]);
        fprintf(fid, '    num_output: %d\n', size(net.layers{idx}.weights{1},4));
        write_stride(fid, net.layers{idx}.stride);
        if isfield(net.layers{idx}, 'pad') && numel(net.layers{idx}.pad) == 4
          % Make sure pad is symmetrical
          if any(net.layers{idx}.pad([1, 3]) ~= net.layers{idx}.pad([2, 3]))
            error('Caffe only supports symmetrical padding');
          end
        end
        write_pad(fid, net.layers{idx}.pad);
        numGroups = layerInputSize(3) / size(net.layers{idx}.weights{1}, 3);
        assert(mod(numGroups, 1) == 0);
        if numGroups > 1
          fprintf(fid, '    group: %d\n', numGroups);
        end

        fprintf(fid, '  }\n');
      elseif filtH == layerInputSize(1) && filtW == layerInputSize(2)
        % Fully connected layer
        isFullyConnected(idx) = true;
        fprintf(fid, '  type: "InnerProduct"\n');
        write_connection(fid, net.layers, idx);
        fprintf(fid, '  inner_product_param {\n');
        fprintf(fid, '    num_output: %d\n', ...
          size(net.layers{idx}.weights{1}, 4));
        fprintf(fid, '  }\n');
      else
        error('Filter size (%d,%d) is larger than input size (%d,%d)', ...
          filtH, filtW, layerInputSize(1), layerInputSize(2))
      end

    case 'relu'
      fprintf(fid, '  type: "ReLU"\n');
      write_connection(fid, net.layers, idx);

    case 'pool'
      fprintf(fid, '  type: "Pooling"\n');
      % Check padding compatability with caffe. See:
      % http://www.vlfeat.org/matconvnet/matconvnet-manual.pdf
      % for more details.
      if numel(net.layers{idx}.pad) == 1
        net.layers{idx}.pad = repmat(net.layers{idx}.pad, 1, 4);
      end
      if numel(net.layers{idx}.stride) == 1
        net.layers{idx}.stride = repmat(net.layers{idx}.stride, 1, 2);
      end
      if numel(net.layers{idx}.pool) == 1
        net.layers{idx}.pool = repmat(net.layers{idx}.pool, 1, 2);
      end
      pad = net.layers{idx}.pad;
      if pad([2, 4]) == net.layers{idx}.pool - 1
        pad([2, 4]) = 0;
      else
        pad([2, 4]) = pad([2, 4]) - net.layers{idx}.stride + 1;
      end
      % Some older versions did not use these upper bounds
      pad = max(pad, 0);

      write_connection(fid, net.layers, idx);
      fprintf(fid, '  pooling_param {\n');
      switch (net.layers{idx}.method)
        case 'max'
          caffe_pool = 'MAX';
        case 'avg'
          caffe_pool = 'AVE';
        otherwise
          error('Unknown pooling type');
      end
      fprintf(fid, '    pool: %s\n', caffe_pool);
      write_kernel(fid, net.layers{idx}.pool);
      write_stride(fid, net.layers{idx}.stride);
      write_pad(fid, pad);
      fprintf(fid, '  }\n');

    case {'normalize', 'lrn'}
      % MATLAB param = [local_size, kappa, alpha/local_size, beta]
      fprintf(fid, '  type: "LRN"\n');
      write_connection(fid, net.layers, idx);
      fprintf(fid, '  lrn_param {\n');
      fprintf(fid, '    local_size: %d\n', net.layers{idx}.param(1));
      fprintf(fid, '    k: %f\n', net.layers{idx}.param(2));
      fprintf(fid, '    alpha: %f\n', net.layers{idx}.param(3)*net.layers{idx}.param(1));
      fprintf(fid, '    beta: %f\n', net.layers{idx}.param(4));
      fprintf(fid, '  }\n');

    case 'softmax'
      fprintf(fid, '  type: "Softmax"\n');
      write_connection(fid, net.layers, idx);

    case {'loss', 'softmaxloss'}
      fprintf(fid, '  type: "SoftmaxWithLoss"\n');
      write_connection(fid, net.layers, idx, true);

    case 'dropout'
      fprintf(fid, '  type: "Dropout"\n');
      write_connection(fid, net.layers, idx);
      fprintf(fid, '  dropout_param {\n');
      fprintf(fid, '    dropout_ratio: %d\n', net.layers{idx}.rate);
      fprintf(fid, '  }\n');

    otherwise
      error('Unknown layer type: %s', net.layers{idx}.type);
  end
  fprintf(fid,'}\n\n');
  layer = struct('layers', {net.layers(idx)});
  res = vl_simplenn(layer, dummyData);
  dummyData = res(end).x;
end
fclose(fid);
info('Network definition exported to: %s.\n', prototxtFilename);

% -------------------------------------------------------------------------
%                                                         Export caffemodel
% -------------------------------------------------------------------------
caffe.set_mode_cpu();
caffeNet = caffe.Net(prototxtFilename, 'test');
firstConv = true;
for idx = 1:numel(net.layers)
  layer_type = net.layers{idx}.type;
  layer_name = net.layers{idx}.name;
  switch layer_type
    case 'conv'
      filters = net.layers{idx}.weights{1};
      % Convert from HxWxCxN to WxHxCxN per Caffe's convention
      filters = permute(filters, [2 1 3 4]);
      if firstConv
        if size(filters, 3) == 3
          % We assume this is RGB Conv., need to convert RGB to BGR
          filters = filters(:,:, [3 2 1], :);
        end
        firstConv = false; % Do this only for first convolution;
      end
      if isFullyConnected(idx)
        % Fully connected layer, squeeze to 2 dims
        filters = reshape(filters, [], size(filters, 4));
      end
      biases = net.layers{idx}.weights{2}(:);
      caffeNet.layers(layer_name).params(1).set_data(filters); % set weights
      caffeNet.layers(layer_name).params(2).set_data(biases); % set bias
    case {'relu', 'normalize', 'lrn', 'pool', 'softmax'}
      % No weights - nothing to do
    otherwise
      error('Unknown layer type %s', layer_type)
  end
end
modelFilename = [caffeModelBaseName '.caffemodel'];
caffeNet.save(modelFilename);
delete(caffeNet);
info('Model file exported to: %s.\n', modelFilename);

% -------------------------------------------------------------------------
%                                                        Export mean image
% -------------------------------------------------------------------------
if ~isempty(avImage)
  if size(avImage, 1) == 1 && size(avImage, 2) == 1
    % Single value, we'll duplicate it to im_size
    avImage = repmat(avImage, imSize(1), imSize(2));
  end
  avImage = matlab_img_to_caffe(avImage);
  meanFilename = [caffeModelBaseName, '_mean_image.binaryproto'];
  caffe.io.write_mean(avImage, meanFilename)
  info('Mean image exported to: %s.\n', meanFilename);
end

if opts.doTest
  simplenn_caffe_compare(net, caffeModelBaseName, opts.testData);
end

  function write_stride(fid, stride)
      if numel(stride) == 1
        fprintf(fid, '    stride: %d\n', stride);
      elseif numel(stride) == 2
        fprintf(fid, '    stride_h: %d\n', stride(1));
        fprintf(fid, '    stride_w: %d\n', stride(2));
      end
  end

  function write_kernel(fid, kernelSize)
    if numel(kernelSize) == 1
      fprintf(fid, '    kernel_size: %d\n', kernelSize);
    elseif numel(kernelSize) == 2
      fprintf(fid, '    kernel_h: %d\n', kernelSize(1));
      fprintf(fid, '    kernel_w: %d\n', kernelSize(2));
    end
  end

  function write_pad(fid, pad)
    if numel(pad) == 1
      fprintf(fid, '    pad: %d\n', pad);
    elseif numel(pad) == 4
      fprintf(fid, '    pad_h: %d\n', pad(1));
      fprintf(fid, '    pad_w: %d\n', pad(2));
    else
      error('pad vector size must be 1 or 4')
    end
  end

  function write_connection(fid, layers, idx, isLoss)
    if idx == 1
      bottom_name = opts.inputBlobName;
    else
      bottom_name = layers{idx-1}.name;
    end
    top_name = layers{idx}.name;
    if idx == numel(layers) && ~isempty(opts.outputBlobName)
      top_name = opts.outputBlobName;
    end
    fprintf(fid, '  bottom: "%s"\n', bottom_name);
    if nargin > 3 && isLoss
      fprintf(fid, '  bottom: "%s"\n', opts.labelBlobName);
    end
    fprintf(fid, '  top: "%s"\n', top_name);
  end

  function img = matlab_img_to_caffe(img)
    img = single(img);
    % Convert from HxWxCxN to WxHxCxN per Caffe's convention
    img = permute(img, [2 1 3 4]);
    if size(img,3) == 3
      % Convert from RGB to BGR channel order per Caffe's convention
      img = img(:,:, [3 2 1], :);
    end
  end
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by: Giulio Tagliaferro
%  Contributors:     ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

function time = getFileStTime(filename)
% Return the start time of the file using the standard naming convention

[~,name, ext] = fileparts(filename);
if strcmpi(ext,'.${YY}p') || ((strcmpi(ext,'.${YY}[n|N]') || ~isempty(regexpi(ext,'\.\d\d[n|N]'))) && isempty(strfind(name, 'CGIM')))  || strcmpi(ext,'.${YY}l') || ~isempty(regexpi(ext,'\.\d\d[p|P]')) || ~isempty(regexp(ext,'\.\d\d[l|L]', 'once')) %#ok<STREMP>
    % name should be : cccwwwwd
    % note: nav file rinex does not have first epoch at the beginning
    % it is too expesivo to look all the file
    if length(name) >= 8
        week = str2double(name(4:7));
        dow = str2double(name(8));
        if ~isnan(week) && ~isnan(dow)
            time = GPS_Time.fromWeekDow(week, dow);
        else
            time = [];
        end
    else
        time = [];
    end

elseif isempty(strfind(lower(ext),lower('eph'))) || isempty(strfind(lower(ext),lower('sp3')))
    % read first 50 lines SP3 headers is 24 (allowing some space for more line comment)
    fid = fopen(filename,'rt');
    if fid > 0
        txt = fread(fid,61*50,'*char')';
        fclose(fid);
        if isempty(txt)
            Core.getLogger.addWarning(sprintf('"%s" is empty ', filename));
            time = [];
        else
            % get new line separators
            nl = regexp(txt, '\n')';
            if nl(end) <  numel(txt)
                nl = [nl; numel(txt)];
            end
            lim = [[1; nl(1 : end - 1) + 1] (nl - 1)];
            lim = [lim lim(:,2) - lim(:,1)];
            if lim(end,3) < 3
                lim(end,:) = [];
            end
            % find time line
            idx_epoch = find(txt(lim(:,1)) == '*');
            txt = txt(lim(idx_epoch):end);
            if ~isempty(idx_epoch)
                time = GPS_Time([str2num(txt(4:7)) str2num(txt(9:10)) str2num(txt(12:13)) str2num(txt(15:16)) str2num(txt(18:19)) str2num(txt(21:31))]);
            else
                time = [];
            end
        end
    else
        time = [];
    end
elseif isempty(strfind(lower(ext),lower('clk')))
    fid = fopen(filename,'rt');
    tline = fgetl(fid);
    i = 0;
    eoh_found = false;
    while ischar(tline) && i < 300 && eoh_found
        eoh_found = strfind(tline,'END OF HEADER');
        tline = fgetl(fid);
        i = i + 1;
    end
    fclose(fid);
    if eoh_found
        tline = fgetl(fid);
        if ischar(tline) && length(tline) >  36
            time = GPS_Time([str2num(tline(9:12)) str2num(tline(14:15)) str2num(tline(17:18)) str2num(tline(20:21)) str2num(tline(24:25)) str2num(tline(26:34))]);
        else
            time = [];
        end
    else
        time = [];
    end

else
    time = [];
end
end

author: goGPS-Project


text:
function [ht, rho, theta] = hough(IM, f, t, M, N)
% HOUGH Hough transform for detection lines in images.
% [HT,RHO,THETA] = HOUGH(IM,M,N).
% From an image IM, computes the integration of the values
% of the image over all the lines. The origin of the coordinates is
% in the middle of x axis (-250<=x<250, 0<=y<=1). The formula for the
% transform is rho = x*cos(theta) + y*sin(theta). Only the values of IM exceeding 5 % of the maximum
%	are taken into account (to speed up the algorithm).
%
%	IM    : image to be analyzed (size Xmax x Ymax).
%	M     : desired number of samples along the radial axis
%					(default : Xmax).
%	N     : desired number of samples along the azimutal (angle) axis
%					(default : Ymax).
%	HT    : output matrix (MxN matrix).
%	RHO   : sequence of samples along the radial axis.
%	THETA : sequence of samples along the azimutal axis.

if (nargin < 3 ),
 error('At least three parameters required');
end;
[Xmax,Ymax] = size(IM);

if (nargin == 3),
M = Xmax; N = Ymax;
elseif (nargin == 4),
N=Ymax;
end;

rhomax = sqrt(Xmax^2 + Ymax^2)/2;
deltar=rhomax/(M-1);
deltat=2*pi/N;

ht=zeros(M,N);
Max=max(max(IM));
for x = 0:Xmax-1,
    for y = 0:Ymax-1,
        if abs(IM(x+1,y+1))>Max/20,
            for theta = 0:deltat:(2*pi-deltat),
                 rho = x*cos(theta) + y*sin(theta);
                if ((rho>=0)&&(rho<=rhomax)),
                    ht(round(rho/deltar)+1,round(theta/deltat)+1) =...
                        ht(round(rho/deltar)+1,round(theta/deltat)+1)+...
                        IM(x+1,y+1);
                end
            end
        end
    end
end
rho=0:deltar:rhomax;
theta=0:deltat:(2*pi-deltat);

author: Sable

repo: mcbench-benchmarks


text:
function display(t)
%DISPLAY Command window display of a tensor.
%
%   DISPLAY(X) displays a tensor with its name.
%
%   See also TENSOR, TENSOR/DISP.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

disp(t,inputname(1));

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
filename='Cantileverbeam_Hexahedra_Linear_Structured_SYM';
ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'holes';
cost = {'compliance'};
weights = [1];
constraint = {'volume'};
optimizer = 'HAMILTON-JACOBI';
incrementFactor = 1;
designVariable = 'LevelSet';
filterType = 'P1';
constraint_case = 'INEQUALITY';

nsteps = 10;
Vfrac_final = 0.1;
Perimeter_target=3.5;
optimality_final =1e-3;
constr_final =1e-3;

BCscale_factor = 0.3;
rotation_per_it = 20;
HJiter0 = 1;
e2 = 100;
N_holes = [12 5 5];
R_holes = 0.2;
phase_holes = [0 0 0];

Vfrac_initial = 0.7;
optimality_initial = 5e-2;
constr_initial = 5e-2;

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

% maxiter = 1;
TOL.nu_plus = 0.3;
TOL.nu_minus = 0.3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% Recording Particle Velocity Example
%
% This example demonstrates how to record the particle velocity using a
% Cartesian or binary sensor mask. It builds on the Homogeneous Propagation
% Medium and Heterogeneous Propagation Medium examples.
%
% author: Bradley Treeby
% date: 1st November 2010
% last update: 20th September 2012
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave 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.
%
% k-Wave 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 Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

clear all;

% =========================================================================
% SIMULATION
% =========================================================================

% create the computational grid
Nx = 128;           % number of grid points in the x (row) direction
Ny = 128;           % number of grid points in the y (column) direction
dx = 0.1e-3;        % grid point spacing in the x direction [m]
dy = 0.1e-3;        % grid point spacing in the y direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy);

% define the properties of the propagation medium
medium.sound_speed = 1500;  % [m/s]
medium.density = 1000;      % [kg/m^3]
medium.alpha_coeff = 0.75;  % [dB/(MHz^y cm)]
medium.alpha_power = 1.5;

% create time array
t_end = 6e-6;
kgrid.t_array = makeTime(kgrid, medium.sound_speed, [], t_end);

% create initial pressure distribution using makeDisc
disc_magnitude = 5; % [au]
disc_x_pos = Nx/2;  % [grid points]
disc_y_pos = Ny/2;  % [grid points]
disc_radius = 5;    % [grid points]
source.p0 = disc_magnitude*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);

% define four sensor points centered about source.p0
sensor_radius = 40; % [grid points]
sensor.mask = zeros(Nx, Ny);
sensor.mask(Nx/2 + sensor_radius, Ny/2) = 1;
sensor.mask(Nx/2 - sensor_radius, Ny/2) = 1;
sensor.mask(Nx/2, Ny/2 + sensor_radius) = 1;
sensor.mask(Nx/2, Ny/2 - sensor_radius) = 1;

% set the acoustic variables that are recorded
sensor.record = {'p', 'u'};

% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor);

% =========================================================================
% VISUALISATION
% =========================================================================

% plot the initial pressure and sensor distribution
figure;
imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, source.p0 + sensor.mask, [-1 1]);
colormap(getColorMap);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;

% plot the simulated sensor data
figure;
[t, t_sc, t_prefix] = scaleSI(kgrid.t_array(end));
mx = 5e-7;
for sensor_num = 1:4
    % plot the pressure
    subplot(4, 3, 3*sensor_num - 2), plot(t_sc*kgrid.t_array, sensor_data.p(sensor_num, :), 'k-');
    set(gca, 'YLim', [-0.75, 0.75], 'XLim', [0 t_end*t_sc]);
    xlabel(['time [' t_prefix 's]']);
    ylabel('p');

    % plot the particle velocity ux
    subplot(4, 3, 3*sensor_num - 1), plot(t_sc*kgrid.t_array, sensor_data.ux(sensor_num, :), 'k-');
    set(gca, 'YLim', [-mx, mx], 'XLim', [0 t_end*t_sc]);
    xlabel(['time [' t_prefix 's]']);
    ylabel('ux');

    % plot the particle velocity uz
    subplot(4, 3, 3*sensor_num), plot(t_sc*kgrid.t_array, sensor_data.uy(sensor_num, :), 'k-');
    set(gca, 'YLim', [-mx, mx], 'XLim', [0 t_end*t_sc]);
    xlabel(['time [' t_prefix 's]']);
    ylabel('uy');
end

repo: Aero-matlab


function d2img = draw_2Dskeleton_ITOP(db,fileId,coord_pixel)

    db_path = strcat('/home/gyeongsikmoon/workspace/Data/Human_pose_estimation/ITOP/',db,'_view/','ITOP_',db,'_test_depth_map.h5');
    jointNum =  15;
    cubicSz = 2.0;
    imgHeight = 240;
    imgWidth = 320;

    annot_path = strcat('/home/gyeongsikmoon/workspace/Data/Human_pose_estimation/ITOP/',db,'_view/','ITOP_',db,'_test_labels.h5');

    coord_pixel = squeeze(coord_pixel);
    rgb_img = zeros(imgHeight,imgWidth,3);
    refDepths = zeros(1,jointNum);
    line_width = 6;

    gt = h5read(annot_path,'/real_world_coordinates',[1 1 fileId],[3 jointNum 1]);
    for jid = 1:jointNum
        refDepths(1,jid) = gt(3,jid,1);
    end
    refDepth = (min(refDepths(:)) + max(refDepths(:)))/2;

    img = h5read(db_path,'/data',[1 1 fileId],[imgWidth imgHeight 1]);
    img = permute(img,[2,1]);
    img(img==0) = refDepth+cubicSz/2;

    if strcmp(db,'top')
        mask = h5read(annot_path,'/segmentation',[1 1 fileId],[imgWidth imgHeight 1]);
        mask = permute(mask,[2,1]);
        img(mask==-1) = refDepth+cubicSz/2;
    end

    if strcmp(db,'side')
        if 1 <= min(coord_pixel(2,:)-15)
            img(1:round(min(coord_pixel(2,:)))-15,:) = refDepth+cubicSz/2;
        end
        if max(coord_pixel(2,:))+15 <= imgHeight
            img(round(max(coord_pixel(2,:)))+15:imgHeight,:) = refDepth+cubicSz/2;
        end
        if 1 <= min(coord_pixel(1,:)) - 15
            img(:,1:round(min(coord_pixel(1,:)))-15) = refDepth+cubicSz/2;
        end
        if max(coord_pixel(1,:)) + 15 <= imgWidth
            img(:,round(max(coord_pixel(1,:)))+15:imgWidth) = refDepth+cubicSz/2;
        end

        img(img < min(coord_pixel(3,:))-0.15) = refDepth+cubicSz/2;
        img(img > max(coord_pixel(3,:))+0.15) = refDepth+cubicSz/2;
    end

    img(img>refDepth+cubicSz/2) = refDepth + cubicSz/2;
    img(img<refDepth-cubicSz/2) = refDepth - cubicSz/2;
    img = img - refDepth;
    img = img/(cubicSz/2);

    img = img + 1;
    img = img/2;

    rgb_img(:,:,1) = img*255/255;
    rgb_img(:,:,2) = img*240/255;
    rgb_img(:,:,3) = img*204/255;

    f = figure;
    set(f, 'visible', 'off');
    imshow(rgb_img);
    hold on;

    plot([coord_pixel(1,1),coord_pixel(1,2)],[coord_pixel(2,1),coord_pixel(2,2)],'Color',[255/255,178/255,102/255],'LineWidth',line_width) %neck to head

    plot([coord_pixel(1,2),coord_pixel(1,3)],[coord_pixel(2,2),coord_pixel(2,3)],'Color',[255/255,153/255,153/255],'LineWidth',line_width) %neck to r_shoulder
    plot([coord_pixel(1,3),coord_pixel(1,5)],[coord_pixel(2,3),coord_pixel(2,5)],'Color',[255/255,102/255,102/255],'LineWidth',line_width) %r_shoulder to r_elbow
    plot([coord_pixel(1,5),coord_pixel(1,7)],[coord_pixel(2,5),coord_pixel(2,7)],'Color',[255/255,51/255,51/255],'LineWidth',line_width) %r_elbow to r_hand

    plot([coord_pixel(1,2),coord_pixel(1,4)],[coord_pixel(2,2),coord_pixel(2,4)],'Color',[153/255,255/255,153/255],'LineWidth',line_width) %neck to l_shoulder
    plot([coord_pixel(1,4),coord_pixel(1,6)],[coord_pixel(2,4),coord_pixel(2,6)],'Color',[102/255,255/255,102/255],'LineWidth',line_width) %l_shoulder to l_elbow
    plot([coord_pixel(1,6),coord_pixel(1,8)],[coord_pixel(2,6),coord_pixel(2,8)],'Color',[51/255,255/255,51/255],'LineWidth',line_width) %l_elbow to l_hand

    plot([coord_pixel(1,2),coord_pixel(1,9)],[coord_pixel(2,2),coord_pixel(2,9)],'Color',[230/255,230/255,0/255],'LineWidth',line_width) %neck to torso

    plot([coord_pixel(1,9),coord_pixel(1,10)],[coord_pixel(2,9),coord_pixel(2,10)],'Color',[255/255,153/255,255/255],'LineWidth',line_width) %torso to r_hip
    plot([coord_pixel(1,10),coord_pixel(1,12)],[coord_pixel(2,10),coord_pixel(2,12)],'Color',[255/255,102/255,255/255],'LineWidth',line_width) %r_hip to r_knee
    plot([coord_pixel(1,12),coord_pixel(1,14)],[coord_pixel(2,12),coord_pixel(2,14)],'Color',[255/255,51/255,255/255],'LineWidth',line_width) %r_knee to r_foot

    plot([coord_pixel(1,9),coord_pixel(1,11)],[coord_pixel(2,9),coord_pixel(2,11)],'Color',[153/255,204/255,255/255],'LineWidth',line_width) %torso to l_hip
    plot([coord_pixel(1,11),coord_pixel(1,13)],[coord_pixel(2,11),coord_pixel(2,13)],'Color',[102/255,178/255,255/255],'LineWidth',line_width) %l_hip to l_knee
    plot([coord_pixel(1,13),coord_pixel(1,15)],[coord_pixel(2,13),coord_pixel(2,15)],'Color',[51/255,153/255,255/255],'LineWidth',line_width) %l_knee to l_foot

    colorList = [
    255/255 178/255 102/255;

    230/255 230/255 0/255;

    255/255 153/255 153/255;
    153/255 255/255 153/255;

    255/255 102/255 102/255;
    102/255 255/255 102/255;

    255/255 51/255 51/255;
    51/255 255/255 51/255;

    230/255 230/255 0/255;

    255/255 153/255 255/255;
    153/255 204/255 255/255;

    255/255 102/255 255/255;
    102/255 178/255 255/255;

    255/255 51/255 255/255;
    51/255 153/255 255/255];
    scatter(coord_pixel(1,:),coord_pixel(2,:),150,colorList,'filled');

    set(gca,'Units','normalized','Position',[0 0 1 1]);  %# Modify axes size
    set(gcf,'Units','pixels','Position',[200 200 2*imgWidth 2*imgHeight]);  %# Modify figure size

    frame = getframe(gcf);
    framedata = frame.cdata;

    xmin = min(coord_pixel(1,:));
    xmax = max(coord_pixel(1,:));
    ymin = min(coord_pixel(2,:));
    ymax = max(coord_pixel(2,:));

    len = max(xmax-xmin+1,ymax-ymin+1) + 70;
    xcenter = (xmin + xmax)/2;
    ycenter = (ymin + ymax)/2;

    xmin = max(round(xcenter - len/2),1);
    xmax = min(round(xmin + len),imgWidth);
    ymin = max(round(ycenter - len/2),1);
    ymax = min(round(ymin + len),imgHeight);
    framedata = framedata(2*ymin:2*ymax,2*xmin:2*xmax,:);

    hold off;
    close(f);

    d2img = framedata;

end


text:
function [percentCorrect,correctLogical] = prtScorePercentCorrect(dataSet1,dataSet2)
% PERCENTCORRECT = prtScorePercentCorrect(GUESS, TRUTH)
%
%   PERCENTCORRECT = prtScorePercentCorrect(GUESS, TRUTH) returns the
%   percent of correct guesses in GUESS as compared to the truth in TRUTH.
%   GUESS and TRUTH should both be Nx1 vectors. The elements of both
%   TRUTH and GUESS should be binary or interger class labels.
%
%   Example:
%   TestDataSet = prtDataGenUnimodal;       % Create some test and
%   TrainingDataSet = prtDataGenUnimodal;   % training data
%   classifier = prtClassMap;               % Create a classifier
%   % Use minimum probablity of error rule
%   classifier.internalDecider = prtDecisionBinaryMinPe;
%   classifier = classifier.train(TrainingDataSet);    % Train
%   classified = run(classifier, TestDataSet);         % Test
%   percentCorr = prtScorePercentCorrect(classified,TestDataSet)
%
%   See also prtScoreConfusionMatrix, prtScoreRoc, prtScoreRmse

if nargin < 2
    dataSet2 = dataSet1;
end

[guesses,targets] = prtUtilScoreParseFirstTwoInputs(dataSet1,dataSet2);

if size(guesses,2) ~= 1
    error('prt:prtScorePercentCorrect','GUESS must be a N x 1 integer vector of class guesses');
else
    percentCorrect = mean(guesses == targets);
    if nargout > 1
        correctLogical = guesses == targets;
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
% gerschdisc.m
%
% This function plots the Gershgorin Discs for the matrix A passed as an argument.
% It will also plot the centers of such discs, and the actual eigenvalues
% of the matrix.

function gershdisc(A)

error(nargchk(nargin,1,1));
if size(A,1) ~= size(A,2)
    error('Matrix should be square');
    return;
end

% For each row, we say:
for i=1:size(A,1)
    % The circle has center in (h,k) where h is the real part of A(i,i) and
    % k is the imaginary part of A(i,i)   :
    h=real(A(i,i)); k=imag(A(i,i));

    % Now we try to compute the radius of the circle, which is nothing more
    % than the sum of norm of the elements in the row where i != j
    r=0;
    for j=1:size(A,1)
       if i ~= j
           r=r+(norm(A(i,j)));
       end
    end

    % We try to make a vector of points for the circle:
    N=256;
    t=(0:N)*2*pi/N;

    % Now we're able to map each of the elements of this vector into a
    % circle:
    plot( r*cos(t)+h, r*sin(t)+k ,'-');

    % We also plot the center of the circle for better undesrtanding:
    hold on;
    plot( h, k,'+');
end

% For the circles to be better graphed, we would like to have equal axis:
axis equal;

% Now we plot the actual eigenvalues of the matrix:
ev=eig(A);
for i=1:size(ev)
    rev=plot(real(ev(i)),imag(ev(i)),'ro');
end
hold off;
legend(rev,'Actual Eigenvalues');

end


text:
function mv = mv_ampsSubplots(mv, sz, runs, subtractMeans)
%
% mv = mv_ampsSubplots([mv], [sz], [runs], [subtractMeans=0]);
%
% Plot amplitude response vectors for each condition in a
% separate subplot for comparison.
%
% If sz = [rows cols] is provided, will arrange the subplots in the
% sepcified number of subplots. Otherwise, arranges the subplots
% in an approximately square array.
%
% If runs is provided, only computes mv_amps for those runs. Otherwise
% uses all runs.
%
% subtractMeans is a flag to remove the mean response across conditions
% for each voxel (as in Haxby et al., 2001). Default is 0.
%
% ras, 01/2007.
if notDefined('mv'),    mv = get(gcf, 'UserData');      end
if notDefined('sz'),    sz = [];                        end
if notDefined('runs'),  runs = unique(mv.trials.run);   end
if notDefined('subtractMeans'),	subtractMeans = 0;		end

amps = mv_amps(mv, runs);
nVoxels = size(amps, 1);
nConds = size(amps, 2);
sel = find(tc_selectedConds(mv));

if subtractMeans==1
	amps = amps - repmat(nanmean(amps, 2), [1 nConds]);
end

if checkfields(mv, 'ui', 'fig')
	delete(findobj('Parent', mv.ui.fig, 'Type', 'axes'));
end

% compute an order of axes to plot
if ~isempty(sz)
    if length(sz)==2
        nRows = sz(1);
        nCols = sz(2);
    elseif length(sz)==nConds
        % explicit specify axes order
    end
else
    nRows = ceil(sqrt(nConds));
    nCols = ceil(nConds/nRows);
end

ysz = [min(amps(:)) max(amps(:))];

% plot
for c = 1:nConds
    hax(c) = subplot(nRows, nCols, c);

    mv_sparkline(1:nVoxels, amps(:,c), ysz);

    if c == ((nRows-1)*nCols + 1) % lower left-hand corner
        axis on
        xlabel('Voxels', 'FontName', mv.params.font);
        ylabel('% Signal', 'FontName', mv.params.font);
        axis tight
    end
    title(mv.trials.condNames{sel(c)}, 'FontName', mv.params.font, ...
          'FontSize', mv.params.fontsz);
end
normAxes;

axes( hax( (nRows-1)*nCols + 1 ) );
% tuftify;

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ x, know ] = p00_sol ( problem, m, know )

%*****************************************************************************80
%
%% P00_SOL returns known solutions for a problem specified by index.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROBLEM, the problem index.
%
%    Input, integer M, the order of the problem.
%
%    Input/output, integer KNOW.
%    On input, KNOW is 0, or the index of the previously returned solution.
%    On output, KNOW is 0 if there are no more solutions, or it is the
%    index of the next solution.
%
%    Output, real X(M), the solution.
%
  if ( problem == 1 )
    [ x, know ] = p01_sol ( m, know );
  elseif ( problem == 2 )
    [ x, know ] = p02_sol ( m, know );
  elseif ( problem == 3 )
    [ x, know ] = p03_sol ( m, know );
  elseif ( problem == 4 )
    [ x, know ] = p04_sol ( m, know );
  elseif ( problem == 5 )
    [ x, know ] = p05_sol ( m, know );
  elseif ( problem == 6 )
    [ x, know ] = p06_sol ( m, know );
  elseif ( problem == 7 )
    [ x, know ] = p07_sol ( m, know );
  elseif ( problem == 8 )
    [ x, know ] = p08_sol ( m, know );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_SOL - Fatal error!\n' );
    fprintf ( 1, '  Problem index out of bounds.\n' );
    error ( 'P00_SOL - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function test_suite = test_conjugate_gradients
    clear all;
    initTestSuite;
end

function test_1
    A = [3 2; 2 6];
    X = [2 1 0 4; -2 2 0 4];
    B = A * X;
    solver = SPX_ConjugateDescent(A, B);
    x = solver.solve();
    %solver.printResults();
    verifyTrue(testCase, solver.hasConverged());
end

function test_2
    A  = [
    2 -1 0;
    -1 2 -1;
    0 -1 2];
    X = randi(10, 3,3);
    B = A * X;
    solver = SPX_ConjugateDescent(A, B);
    solver.MaxIterations = 50;
    x = solver.solve();
    verifyTrue(testCase, solver.hasConverged());
end

author: indigits

repo: sparse-plex


text:
function n=length(X)
%LENGTH (overloaded)

n = max(X.dim);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [all_bboxes_out] = localize_bboxes_of_all_imgs(...
    model, image_paths, all_bboxes_in, dst_directory, image_set_name, varargin)
% localize_bboxes_of_all_imgs: given a bounding box localization model and
% a set of images with their input candidate bounding boxes, for each
% image it predicts new bounding box coordinates such that the new boxes
% will be close (i.e. better localize) the objects of interest.
%
% INPUTS:
% 1) model: (type struct) the bounding box localization model
% 2) image_paths: NI x 1 cell array with the image paths of the set in the
% form of strings. NI is the number of images.
% 3) all_bboxes_in: a cell array with the candidate bounding
% boxes that are going to be processed. It can be of two forms:
%   a) (the flag is_per_class must being set to false) a NI x 1 cell array
%   where NI is the number of images. The all_bboxes_in{i} element
%   contains the candidate bounding boxes of the i-th image in the form of a
%   NB_i x 4 array where NB_i is the number of candidate bounding boxes of
%   the i-th image. The 4 columns of this array contain the bounding box
%   coordinates in the form of [x0,y0,x1,y1] (where the (x0,y0) and (x1,y1)
%   are the top-left and bottom-right corners).
%   b) (the flag is_per_class must being set to true) a C x 1 cell array
%   where C is the number of categories. The all_bboxes_in{j} element
%   in this case is NI x 1 cell array, where NI is the number of images,
%   with the candidate bounding boxes of the j-th category for each image.
%   The element all_bboxes_in{j}{i} is a NB_{i,j} x 4 array, where
%   NB_{i,j} is the number candidate bounding box of the i-th
%   image for the j-th category. The 4 columns of this array are the
%   bounding boxes in the form of [x0,y0,x1,y1] (where the (x0,y0) and (x1,y1)
%   are the top-left and bottom-right corners).
% 4) dst_directory: string with the path of the destination directory where
% the results are going to be cached.
% 5) image_set_name: string with the name of the image set that is being
% processed. It will be used on the name of the file where the results will
% be saved.
% 6) The rest input arguments are given in the form of Name,Value pair
% arguments and are:
% 'is_per_class': boolean value that if set to false, then the 3.a) form
% of the all_bboxes_in input parameter will be expected; otherwise, if set to
% true then the 3.b) form of the all_bboxes_in input parameter will be
% expected. Default value: false
% 'suffix': string with a suffix that it will be used on the name of the
% file where the results will be saved. Default value: ''
% 'all_bbox_gt': (optional) a NI cell array with the ground truth bounding
% boxes of each image; it will be used (if given) for printing the mAP on
% regular (see param 'checkpoint_step') step during processing the image set
% 'checkpoint_step': scalar value; during processing the NI images, after each
% checkpoint_step number of images the results are cached and the mAP
% results are printed. Default value: 500
%
% This file is part of the code that implements the following paper:
% Title      : "LocNet: Improving Localization Accuracy for Object Detection"
% Authors    : Spyros Gidaris, Nikos Komodakis
% Institution: Universite Paris Est, Ecole des Ponts ParisTech
% ArXiv link : http://arxiv.org/abs/1511.07763
% code       : https://github.com/gidariss/LocNet
%
% AUTORIGHTS
% --------------------------------------------------------
% Copyright (c) 2016 Spyros Gidaris
%
% Title     : "LocNet: Improving Localization Accuracy for Object Detection"
% ArXiv link: http://arxiv.org/abs/1511.07763
% Licensed under The MIT License [see LICENSE for details]
% ---------------------------------------------------------

ip = inputParser;
ip.addParamValue('is_per_class',    true,  @islogical);
ip.addParamValue('suffix',          '',    @ischar);
ip.addParamValue('all_bbox_gt',     {},    @iscell);
ip.addParamValue('checkpoint_step', 500,   @isnumeric);

ip.parse(varargin{:});
opts = ip.Results;

mkdir_if_missing(dst_directory);

filepath_results     = [dst_directory, filesep, 'loc_boxes_', opts.suffix, image_set_name, '.mat'];
in_progress_filepath = [dst_directory, filesep, 'loc_boxes_', opts.suffix, image_set_name, '_in_progress.mat'];

timestamp       = datestr(datevec(now()), 'yyyymmdd_HHMMSS');
log_file        = fullfile(dst_directory, 'output', ...
    ['log_file_localization',opts.suffix, image_set_name, '_', timestamp, '.txt']);
mkdir_if_missing(fileparts(log_file));

t_start = tic();
try
    all_bboxes_out = load_bboxes_scores(filepath_results);
catch
    diary(log_file);
    all_bboxes_out = process_images(...
        model, image_paths, all_bboxes_in, opts.all_bbox_gt, ...
        opts.is_per_class, in_progress_filepath, opts.checkpoint_step);

    save_bboxes_scores(filepath_results, all_bboxes_out);

    delete(in_progress_filepath);
    diary off;
end
fprintf('Localize bounding box proposals in %.4f minutes.\n', toc(t_start)/60);
end

function save_bboxes_scores(filename, bboxes_scores)
save(filename, 'bboxes_scores', '-v7.3');
end

function bboxes_scores = load_bboxes_scores(filename)
load(filename, 'bboxes_scores');
end

function [all_bboxes_out] = process_images(...
    model, image_paths, all_bboxes_in, all_bbox_gt, is_per_class, ...
    in_progress_filepath, checkpoint_step)

num_imgs      = length(image_paths);
nms_iou_thrs  = 0.3;
max_per_set   = 5 * num_imgs;
max_per_image = 100;
num_classes   = length(model.classes);
in_progress_filepath_prev = [in_progress_filepath, '.prev'];

try
    % load (if any) cached subset of results from images that have being
    % already processed from a previous run of the script that was stopped
    % before finish processing all the images.
    try
        [fist_img_idx, aboxes, all_bboxes_out, thresh] = load_progress(in_progress_filepath);
    catch
        [fist_img_idx, aboxes, all_bboxes_out, thresh] = load_progress(in_progress_filepath_prev);
    end
    mAP   = check_progress_on_mAP(aboxes, all_bbox_gt, fist_img_idx-1, model.classes);
    mAP_i = fist_img_idx-1;
catch exception
    fprintf('Exception message %s\n', getReport(exception));
    aboxes = cell(num_classes, 1);

    for i = 1:num_classes, aboxes{i} = cell(num_imgs, 1); end
    if is_per_class
        all_bboxes_out = cell(num_classes,1);
        for i = 1:num_classes, all_bboxes_out{i} = cell(num_imgs, 1); end
    else
        all_bboxes_out = cell(num_imgs,1);
    end
    thresh = -2.5 * ones(num_classes, 1);
    fist_img_idx = 1; mAP = 0; mAP_i = 0;
end

model.max_rois_num_in_gpu = 500; %find_max_rois_num_in_gpu(model);
num_chars = 0;
total_el_time = 0;
for i = fist_img_idx:num_imgs
    th = tic;
    % get the bounding box proposals of this image
    [bbox_proposals, bbox_scores] = get_this_img_bbox_proposals(all_bboxes_in, i, is_per_class);
    image = get_image(image_paths{i}); % read the image
    % refine the location of each bounding box with the localization model
    bbox_refined = localize_bboxes_of_image(model, image, bbox_proposals);
    localization_time = toc(th);

    th = tic;
    % prepare the candidate boundng box detections of this image
    bbox_cand_dets = prepare_bbox_cand_dets(bbox_refined, bbox_scores, num_classes, is_per_class);
    % prepare the output for this image
    all_bboxes_out = prepare_this_img_output(all_bboxes_out, i, is_per_class, bbox_cand_dets);
    % perform the post-processing step of non-max-suppresion on the
    % candidate boundng box detections of this image
    bbox_detections = post_process_candidate_detections(bbox_cand_dets, ...
        'thresholds',thresh, 'nms_iou_thrs',nms_iou_thrs,'use_gpu',true,...
        'max_per_image',max_per_image);

    for j = 1:num_classes, aboxes{j}{i} = bbox_detections{j}; end
    postprocessing_time = toc(th);

    if mod(i, checkpoint_step) == 0
        th = tic;
        for j = 1:num_classes, [aboxes{j}, thresh(j)] = keep_top_k(aboxes{j}, i, max_per_set, thresh(j)); end

        % save the till now progress
        save_progress(aboxes, all_bboxes_out, thresh, i, in_progress_filepath, in_progress_filepath_prev);

        % evaluate the mAP of the images 1 till i
        mAP = check_progress_on_mAP(aboxes, all_bbox_gt, i, model.classes);
        mAP_i = i;
        diary; diary; % flush diary

        %disp(thresh(:)');
        checkpoint_time = toc(th);
        fprintf('check_point time %.2fs\n',checkpoint_time)
        num_chars = 0;
    end

    elapsed_time = localization_time + postprocessing_time;
    [total_el_time, avg_time, est_rem_time] = timing_process(elapsed_time, total_el_time, fist_img_idx, i, num_imgs);
    fprintf(repmat('\b',[1, num_chars]));
    num_chars = fprintf('%s: bbox loc %d/%d:| ET %.3fs + %.3fs | AT: %.3fs | TT %.4fmin | ERT %.4fmin | mAP[%d/%d] = %.4f\n', ...
        procid(), i, num_imgs, localization_time, postprocessing_time, ...
        avg_time, total_el_time/60, est_rem_time/60, mAP_i, num_imgs, mAP);
end

delete(in_progress_filepath_prev);
end

function [total_el_time, ave_time, est_rem_time] = timing_process(...
    elapsed_time, total_el_time, fist_img_idx, i, num_imgs)
total_el_time   = total_el_time + elapsed_time;
ave_time        = total_el_time / (i-fist_img_idx+1);
est_rem_time    = ave_time * (num_imgs - i);
end

function [bbox_proposals, bbox_scores] = get_this_img_bbox_proposals(all_bboxes_in, img_id, is_per_class)
if is_per_class
    bbox_proposals_per_class = cellfun(@(x) x{img_id}(:,1:4), all_bboxes_in, 'UniformOutput', false);
    class_indices = [];
    for c = 1:length(all_bboxes_in)
        class_indices = [class_indices; ones(size(all_bboxes_in{c}{img_id},1),1,'single')*c];
    end
    bbox_scores_per_class = cellfun(@(x) x{img_id}(:,5:end), all_bboxes_in, 'UniformOutput', false);

    num_bbox_per_class  = cellfun(@(x) size(x,1), bbox_proposals_per_class,  'UniformOutput', true);
    bbox_proposals      = cell2mat(bbox_proposals_per_class(num_bbox_per_class>0));
    bbox_scores         = cell2mat(bbox_scores_per_class(num_bbox_per_class>0));
    bbox_proposals      = [bbox_proposals, class_indices];
else
    class_indices  = ones(size(all_bboxes_in{img_id},1),1,'single');
    bbox_proposals = [all_bboxes_in{img_id}(:,1:4),class_indices];
    bbox_scores    = all_bboxes_in{img_id}(:,5:end);
end
end

function bbox_cand_dets = prepare_bbox_cand_dets(bbox_refined, bbox_scores, num_classes, is_per_class)
if is_per_class
    class_indices      = bbox_refined(:,5);
    bbox_cand_dets_tmp = single([bbox_refined(:,1:4),bbox_scores]);
    bbox_cand_dets = cell(num_classes,1);
    for c = 1:num_classes
        bbox_cand_dets{c} = bbox_cand_dets_tmp(class_indices==c,:);
        if isempty(bbox_cand_dets{c}), bbox_cand_dets{c} = zeros(0,5,'single'); end
    end
else
    assert(size(bbox_refined,2)==5);
    assert(all(bbox_refined(:,5) == 1));
    bbox_cand_dets = single([bbox_refined(:,1:4), bbox_scores]);
end
end

function all_bboxes_out = prepare_this_img_output(all_bboxes_out, img_idx, is_per_class, bbox_this_img)
if is_per_class
    for j = 1:length(bbox_this_img), all_bboxes_out{j}{img_idx} = bbox_this_img{j}; end
else
    all_bboxes_out{img_idx} = bbox_this_img;
end
end

function [fist_img_idx, aboxes, abbox_scores, thresh] = load_progress(in_progress_filepath)
ld = load(in_progress_filepath);
fist_img_idx = ld.progress_state.img_idx + 1;
aboxes       = ld.progress_state.aboxes;
abbox_scores = ld.progress_state.abbox_scores;
thresh       = ld.progress_state.thresh;
end

function save_progress(aboxes, abbox_scores, thresh, img_idx, in_progress_filepath, in_progress_filepath_prev)
progress_state              = struct;
progress_state.img_idx      = img_idx;
progress_state.aboxes       = aboxes;
progress_state.abbox_scores = abbox_scores;
progress_state.thresh       = thresh;
if exist(in_progress_filepath, 'file')
    % in case it crash during updating in_progress_filepath
    copyfile(in_progress_filepath, in_progress_filepath_prev);
end
save(in_progress_filepath, 'progress_state', '-v7.3');
end

function mAP = check_progress_on_mAP(aboxes, all_bbox_gt, img_idx, classes)
mAP = 0;
if ~isempty(all_bbox_gt)
    aboxes      = cellfun(@(x) x(1:img_idx), aboxes, 'UniformOutput', false);
    mAP_result  = evaluate_average_precision_pascal( all_bbox_gt(1:img_idx), aboxes, classes );
    printAPResults(classes, mAP_result);
    mAP = mean([mAP_result(:).ap]');
end
end

function [boxes, thresh] = keep_top_k(boxes, end_at, top_k, thresh)
% ------------------------------------------------------------------------
% Keep top K
X = cat(1, boxes{1:end_at});
if isempty(X), return; end

scores = sort(X(:,end), 'descend');
thresh = scores(min(length(scores), top_k));
for image_index = 1:end_at
    bbox = boxes{image_index};
    keep = find(bbox(:,end) >= thresh);
    boxes{image_index} = bbox(keep,:);
end

end

author: gidariss

repo: LocNet

save_path: github-repos/MATLAB/gidariss-LocNet


text:
classdef Tonemap < handle
    %TONEMAP  Tonemapping algorithm used to map HDR image to 8-bit range
    %
    % See also: cv.TonemapDrago, cv.TonemapDurand, cv.TonemapReinhard,
    %  cv.TonemapMantiuk, tonemap
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    properties (Dependent)
        % positive value for gamma correction.
        %
        % Gamma value of 1.0 implies no correction, gamma equal to 2.2 is
        % suitable for most displays. Generally gamma > 1 brightens the image
        % and gamma < 1 darkens it.
        Gamma
    end

    %% Tonemap
    methods
        function this = Tonemap(varargin)
            %TONEMAP  Creates simple linear mapper with gamma correction
            %
            %     obj = cv.Tonemap()
            %     obj = cv.Tonemap('OptionName',optionValue, ...)
            %
            % ## Options
            % * __Gamma__ default 1.0
            %
            % See also: cv.Tonemap.process
            %
            this.id = Tonemap_(0, 'new', varargin{:});
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.Tonemap
            %
            if isempty(this.id), return; end
            Tonemap_(this.id, 'delete');
        end
    end

    %% Tonemap
    methods
        function dst = process(this, src)
            %PROCESS  Tonemaps image
            %
            %     dst = obj.process(src)
            %
            % ## Input
            % * __src__ source RGB image, 32-bit `single` 3-channel array.
            %
            % ## Output
            % * __dst__ destination image of same size as `src`, 32-bit
            %   `single` 3-channel array with values in [0,1] range.
            %
            % See also: cv.Tonemap.Tonemap
            %
            dst = Tonemap_(this.id, 'process', src);
        end
    end

    %% Algorithm
    methods (Hidden)
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.Tonemap.empty, cv.Tonemap.load
            %
            Tonemap_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Returns true if the algorithm is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the object is empty (e.g in the
            %   very beginning or after unsuccessful read).
            %
            % See also: cv.Tonemap.clear, cv.Tonemap.load
            %
            b = Tonemap_(this.id, 'empty');
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.Tonemap.save, cv.Tonemap.load
            %
            name = Tonemap_(this.id, 'getDefaultName');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm parameters to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in the specified
            % XML or YAML file.
            %
            % See also: cv.Tonemap.load
            %
            Tonemap_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from the specified XML or
            % YAML file (either from disk or serialized string). The previous
            % algorithm state is discarded.
            %
            % See also: cv.Tonemap.save
            %
            Tonemap_(this.id, 'load', fname_or_str, varargin{:});
        end
    end

    %% Getters/Setters
    methods
        function value = get.Gamma(this)
            value = Tonemap_(this.id, 'get', 'Gamma');
        end
        function set.Gamma(this, value)
            Tonemap_(this.id, 'set', 'Gamma', value);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function a = r8cb_indicator ( m, n, ml, mu )

%*****************************************************************************80
%
%% R8CB_INDICATOR sets up a R8CB indicator matrix.
%
%  Discussion:
%
%    The R8CB storage format is appropriate for a compact banded matrix.
%    It is assumed that the matrix has lower and upper bandwidths ML and MU,
%    respectively.  The matrix is stored in a way similar to that used
%    by LINPACK and LAPACK for a general banded matrix, except that in
%    this mode, no extra rows are set aside for possible fillin during pivoting.
%    Thus, this storage mode is suitable if you do not intend to factor
%    the matrix, or if you can guarantee that the matrix can be factored
%    without pivoting.
%
%    The original M by N matrix is "collapsed" downward, so that diagonals
%    become rows of the storage array, while columns are preserved.  The
%    collapsed array is logically ML+MU+1 by N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows of the matrix.
%    M must be positive.
%
%    Input, integer N, the number of columns of the matrix.
%    N must be positive.
%
%    Input, integer ML, MU, the lower and upper bandwidths.
%    ML and MU must be nonnegative, and no greater than min(M,N)-1.
%
%    Output, real A(ML+MU+1,N), the R8CB matrix.
%
  fac = 10^( i4_log_10 ( n ) + 1 );
  k = 0;

  for j = 1 : n
    for diag = 1 : ml + mu + 1

      i = diag + j - mu - 1;

      if ( 1 <= i & i <= m & i - ml <= j & j <= i + mu )
        value = fac * i + j;
      else
        k = k + 1;
        value = -k;
      end

      a(diag,j) = value;

    end
  end

  return
end


text:
function gridplot(xy,ev,bound,ebound)
%gridplot   quadrilateral grid verification
%   gridplot(xy,ev,bound,ebound);
%
%   IFISS function: DJS; 7 May 2006.
% Copyright (c) 2006 D.J. Silvester, H.C. Elman, A. Ramage
fprintf('\nGrid logistics ..\n')
nvtx=length(xy(:,1));
fprintf('  %g nodes \n',nvtx)
nelement=length(ev(:,1));
fprintf('  %g elements \n',nelement)
nboundvtx=length(bound);
fprintf('  %g nodes on Dirichlet boundary \n',nboundvtx)
nboundedge=length(ebound(:,1));
fprintf('  %g element edges on Dirichlet boundary \n\n',nboundedge)
%
%
adj=sparse(nvtx,nvtx); adx=sparse(nvtx,nvtx);
for i=1:nelement
   adj(ev(i,1),ev(i,2)) =1;
   adj(ev(i,2),ev(i,3)) =1;
   adj(ev(i,3),ev(i,4)) =1;
   adj(ev(i,4),ev(i,1)) =1;
end
%
%% define element edges
adjb=sparse(nvtx,nvtx);
% bottom boundary edges
k1=find(ebound(:,2)==1)';
for k=ebound(k1)
   adjb(ev(k,1),ev(k,2))=1;
end
% right boundary edges
k2=find(ebound(:,2)==2)';
for k=ebound(k2)
   adjb(ev(k,2),ev(k,3))=1;
end
% top boundary edges
k3=find(ebound(:,2)==3)';
for k=ebound(k3)
   adjb(ev(k,3),ev(k,4))=1;
end
% left boundary edges
k4=find(ebound(:,2)==4)';
for k=ebound(k4)
   adjb(ev(k,4),ev(k,1))=1;
end
%
figure(1)
gplot(adj,xy,'b')
hold on
stnode=int2str([1:nvtx]');
text(xy(:,1),xy(:,2),stnode)
axis('equal'),axis('off')
title('Indices of nodes of the element grid')
hold off
figure(2)
gplot(adj,xy,'b')
hold on
gplot(adjb,xy,'r')
xybd=xy(bound,:);
stbd=int2str([1:nboundvtx]');
text(xybd(:,1),xybd(:,2),stbd,'color','black')
title('Indices of nodes on the Dirichlet boundary')
axis('equal'),axis('off')
hold off
return


text:
function FBC = FbCrop( FB, delta )
% Crop a 2D filterbank (adjusting filter norms).
%
% Takes a filter bank and crops it by cropping off delta pixels from each
% side. Ensures that the mean response of each filter is 0 and that the L1
% norm is 1, i.e. sum(sum(abs(F)))==1.
%
% USAGE
%  FBC = FbCrop( FB, delta )
%
% INPUTS
%  FB       - original filterbank
%  delta    - amount to crop by
%
% OUTPUTS
%  FBC      - cropped filterbank
%
% EXAMPLE
%  load FbDoG.mat;  FBC=FbCrop(FB,4);
%  figure(1); montage2(FB,struct('extraInfo',1));
%  figure(2); montage2(FBC,struct('extraInfo',1));
%
% See also FBAPPLY2D
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.0
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

nd = ndims(FB);
if( nd~=2 && nd~=3 ); error('I must an MxNxK array'); end

cropsiz = size(FB);
cropsiz = [cropsiz(1:2)-2*delta, cropsiz(3)];
FBC = arrayToDims( FB, cropsiz );

for f=1:size(FB,3)
  FC = FBC(:,:,f);
  FC = FC - sum(sum(FC)) / numel(FC); % 0 mean
  FC = FC / sum(sum(abs(FC))); % L1 norm == 1
  FBC(:,:,f) = FC;
end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function K = gaussian(kern,d1,d2,ind1,ind2,kerParam),

% K = gaussian(d1,d2,ind1,ind2,param), compute the kernel
%     matrix between d1 and d2
% for a gaussian kernel exp(-||x-z||^2/(2*param^2))
%     where x is from d1 and z from d2

   K=get_x(d2,ind2)*get_x(d1,ind1)';
   kertmp=kernel;
   Kdn = get_norm(kertmp,d1,ind1).^2;
   Kn = get_norm(kertmp,d2,ind2).^2;
   K = ones(length(Kn),1)*Kdn' + Kn*ones(1,length(Kdn)) - 2*K;

   [numEx vDim oDim]=get_dim(d1);
   sigma=kerParam;

%   K = 1/(2*(2*pi)^(vDim/2)*sigma) * exp(-0.5*sigma^2*K);

K= 1/(sigma*(2*pi)^(vDim/2))*exp(-K/(2*sigma^2));

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [Bx,By,Bz] = magnetic_field_current_loop(x,y,z,x_p,y_p,z_p,a,I0)

%Kilian O'Donoghue

%30th July 2013

%This function calculates the magnetic field resulting from a single
%circular coil of radius a, carrying a currnet I0. The coil points in the z
%direction. The centre of the coil is located at {x_p,y_p,z_p} and the
%magnetic field is calculated at the point or points stored in the arrays
%x, y and z. This code accepts meshgrid inputs or single values along with
%vectors.

%These equations are based on those given in "General relation for the vector
%magnetic field of a circular current loop: a closer look" by Robert
%Schill, this is available from IEEE.

global u0 %permeability of free space is a global variable

rc=((x-x_p).^2+(y-y_p).^2).^.5; %Radial component is required for cylindrical coordinate system.
m=(4.*a.*rc).*(((rc+a).^2)+((z-z_p).^2)).^(-1); %This is a parameter for calculating the Elliptical integrals
kofkc=(pi/2)+(pi/8).*m+(9*pi/128).*m.^2; %K(k) elliptical function, this is a taylor expansion of the K elliptical integral.
eofkc=(pi/2)+(-pi/8).*m+(-3*pi/128).*m.^2;%E(k) elliptical function this is a taylor expansion of the E elliptical integral.

%Note for improved accuracy, Matlab has built in elliptical integral
%calculation but these expressions here are still very accurate when rc < a

Brc=(u0.*I0./(2.*pi.*rc)).*(z-z_p).*((((rc+a).^2)+((z-z_p).^2)).^(-.5)).*(-kofkc+eofkc.*((rc.^2+a.^2+(z-z_p).^2)./(((rc-a).^2)+((z-z_p).^2)))); %radial component of B%
Bz=(u0.*I0./(2.*pi)).*((((rc+a).^2)+((z-z_p).^2)).^(-.5)).*(kofkc-eofkc.*((rc.^2-a.^2+(z-z_p).^2)./(((rc-a).^2)+((z-z_p).^2)))); %axial component of B
Bx=Brc.*(x-x_p)./rc; %This converts the polar component into cartesian form.
By=Brc.*(y-y_p)./rc;

%The following sets any terms that result in Inf to zero, this occurs at
%the points near the coil itself.
Bx(isnan(Bx)) = 0 ;
By(isnan(By)) = 0 ;
Bz(isnan(Bz)) = 0 ;

path: github-repos/MATLAB/Sable-mcbench-benchmarks/mcbench-benchmarks-ba13b2f0296ef49491b95e3f984c7c41fccdb6d8/43394-magnetic-field-of-circular-coils-using-elliptical-integrals/Magnetic field of circular coil using elliptical integrals/magnetic_field_current_loop.m


text:
function [f,g] = hs71F(x)

f = x(1)*x(4)*sum(x(1:3)) + x(3);

if(nargout > 1)
    g = [ x(1)*x(4) + x(4)*sum(x(1:3))
                    x(1)*x(4)
                    x(1)*x(4) + 1
                    x(1)*sum(x(1:3)) ];
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function M = perform_lic(v, w, options)

% perform_lic - perform line integral convolution
%
%   M = perform_lic(v, w, options);
%
%   v is a vector field (should be approximately of unit norm).
%   w is the length of the convolution (in pixels)
%   M is an image of filtered noise that illustrate well the flow of v.
%
%   options.spot_size set the size of the features.
%
%   Set options.flow_correction=1 in order to fix problem
%   around singularity points of the flow.
%
%   The method is described in
%       Imaging vector fields using line integral convolution
%       Brian Cabral, Leith Casey Leedom
%       Siggraph 1993
%
%   If the vector field is not oriented (e.g. if it is an eigenvector field
%   of a tensor field) then set options.isoriented=1. It will perform two
%   LIC (horizontal and vertical) and then average them. This trick is
%   explained in
%       Interactive Tensor Field Design and Visualization on Surfaces
%       Eugene Zhang, James Hays, and Greg Turk
%       IEEE Transactions on Visualization and Computer Graphics
%       Volume 13 ,  Issue 1, Pages: 94-107, 2007.
%
%   Copyright (c) Gabriel Peyre 2007

n = size(v,1);

options.null = 0;
if isfield(options, 'M0')
    M0 = options.M0;
else
    M0 = randn(n);
    if isfield(options, 'spot_size')
        sigma = options.spot_size;
    else
        sigma = 2;
    end
    M0 = perform_blurring(M0, sigma, options);
end

if size(M0,3)>1
    % lic on each channel
    for i=1:size(M0,3)
        options.M0 = M0(:,:,i);
        M(:,:,i) = perform_lic(v, w, options);
    end
    return;
end

isoriented = getoptions(options, 'isoriented', 1);

if isoriented==0
    % run twice the lic in each direction
    options.isoriented = 1;
    options.method = 'xproj';
    v = perform_vf_reorientation(v, options);
    M1 = perform_lic(v, w, options);
    options.method = 'yproj';
    v = perform_vf_reorientation(v, options);
    M2 = perform_lic(v, w, options);
    % weight
    v = perform_vf_normalization(v);
    T = v(:,:,1).^2;
    M = T.*M1 + (1-T).*M2;
    return;
end

if isfield(options, 'niter_lic') && options.niter_lic>1
    M = options.M0;
    niter_lic = options.niter_lic;
    options.niter_lic = 1;
    for i=1:niter_lic
        options.M0 = M;
        M = perform_lic(v, w, options);
    end
    return;
end

histogram = getoptions(options, 'histogram', 'gaussian');
flow_correction = getoptions(options, 'flow_correction', 1);

if isstr(histogram)
    switch histogram
        case 'linear'
            hist = linspace(0,1, 100^2);
        case 'gaussian'
            hist = randn(100); hist = hist(:);
        otherwise
            error('Unkown kind of histograms');
    end
else
    hist = histogram;
end

if isfield(options, 'dt')
    dt = options.dt;
else
    dt = 0.5;
end

% perform integration of the vector field: Forward
T_list = 0:dt:w/2;
H = perform_vf_integration(v, dt, T_list, options );

% perform integration of the vector field: Backward
T_list(1) = [];
H1 = perform_vf_integration(-v, dt, T_list, options );
H = cat(4, H1(:,:,:,end:-1:1), H );
p = size(H,4);

% try to remove sampling problems
if flow_correction
    A = H(:,:,:,(end+1)/2);
    dX = H(:,:,1,:)-repmat(A(:,:,1,:), [1 1 1 p]);
    dY = H(:,:,2,:)-repmat(A(:,:,2,:), [1 1 1 p]);
    dX(dX>n/2) = dX(dX>n/2)-(n-1); dX(dX<-n/2) = dX(dX<-n/2)+(n-1);
    dY(dY>n/2) = dY(dY>n/2)-(n-1); dY(dY<-n/2) = dY(dY<-n/2)+(n-1);
    d = sqrt(dX.^2 + dY.^2);
    d = repmat( mean(mean(d)), [n n 1] ) - d;
    % threshold on the distance map deviation
    eta = 0.5;
end

% compute averaging
M = zeros(n);
[Y,X] = meshgrid(1:n,1:n);
W = zeros(n);
for i=1:p
    A = interp2(1:n,1:n, M0, H(:,:,2,i), H(:,:,1,i) );
    if flow_correction
        w = d(:,:,i)<eta;
    else
        w = ones(n);
    end
    M = M + A.*w;
    W = W+w;
end
W(W==0) = 1;
M = M./W;

if not(isempty(hist))
    M = perform_histogram_equalization(M, hist);
end

author: gpeyre

repo: matlab-toolboxes


text:
function gX = invcmpndKernGradX(kern, X, X2)

% INVCMPNDKERNGRADX
% See: cmpndKernGradX...
% SEEALSO cmpndKernParamInit, kernGradX, cmpndKernDiagGradX
%

% KERN

error('invcmpndKernGradX not yet implemented!')

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [TrainData, TestData, TrainLabel, TestLabel, numX] = createData()
    ncFeature = 36;  % features of color
    ntFeature = 51;  % features of text

    % load source domain 1
    A = textread(strcat('data/img/11_color36.txt'));
    B = textread(strcat('data/img/11_texture51.txt'));
    n1 = size(A,1)/ncFeature;
    for i = 1:n1
        Cc(:,i) = A((ncFeature*(i-1)+1):ncFeature*i,1);
        D(:,i) = B((ntFeature*(i-1)+1):ntFeature*i,1);
        TrainY1(1,i) = 1;
    end
    A = textread(strcat('data/img/15_color36.txt'));
    B = textread(strcat('data/img/15_texture51.txt'));
    n2 = size(A,1)/ncFeature;
    for i = 1:n2
        Cc(:,i+n1) = A((ncFeature*(i-1)+1):ncFeature*i,1);
        D(:,i+n1) = B((ntFeature*(i-1)+1):ntFeature*i,1);
        TrainY1(1,i+n1) = 2;
    end
    numX = [n1,n2];
    TrainX1 = [Cc;D];
    clear Cc;
    clear D;
    % load target domain
    A = textread(strcat('data/img/12_color36.txt'));
    B = textread(strcat('data/img/12_texture51.txt'));
    n1 = size(A,1)/ncFeature;
    for i = 1:n1
        Cc(:,i) = A((ncFeature*(i-1)+1):ncFeature*i,1);
        D(:,i) = B((ntFeature*(i-1)+1):ntFeature*i,1);
        TestY(1,i) = 1;
    end
    A = textread(strcat('data/img/16_color36.txt'));
    B = textread(strcat('data/img/16_texture51.txt'));
    n2 = size(A,1)/ncFeature;
    for i = 1:n2
        Cc(:,i+n1) = A((ncFeature*(i-1)+1):ncFeature*i,1);
        D(:,i+n1) = B((ntFeature*(i-1)+1):ntFeature*i,1);
        TestY(1,i+n1) = 2;
    end
    TestX = [Cc;D];

    clear A B Cc D;

    %% normalization
    TrainData = TrainX1;
    TrainLabel = TrainY1;
    TestData = TestX;
    TestLabel = TestY;
    column = size(TrainData,2);
    mode_TrainX = sqrt(sum(TrainData.*TrainData,1));
    for i = 1 : column
        TrainData(:,i) = TrainData(:,i)/mode_TrainX(1,i);
    end
    column = size(TestData,2);
    mode_TestX = sqrt(sum(TestData.*TestData,1));
    for i = 1 : column
        TestData(:,i) = TestData(:,i)/mode_TestX(1,i);
    end
    save TrainData.mat TrainData
    save TrainLabel.mat TrainLabel
    save TestData.mat TestData
    save TestLabel.mat TestLabel
    clear TrainX1 TrainY1 TestX TestY
end

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
%% Testing PSwarm

%% Rosenbrock [x = 1,1, fval = 0]
clc
fun = @(x) (1-x(1))^2 + 100 *(x(2) - x(1)^2)^2;
x0 = [0 0]';

opts = [];
opts.display = 2;

[x,fval,ef,iter,feval] = pswarm(fun,x0,-[10;10],[10;10],[],[],opts)

%% Constrained Rosenbrock [x = .9488,.9, fval = .0026]
clc
fun = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
x0 = [0 0]';
lb = [0;0];
ub = [1;0.9];

opts.display = 2;
opts.maxfeval = 1e6;

[x,fval,ef,iter] = pswarm(fun,x0,lb,ub,[],[],opts)


text:
classdef HBoxFlex < uix.HBox & uix.mixin.Flex
    %uix.HBoxFlex  Flexible horizontal box
    %
    %  b = uix.HBoxFlex(p1,v1,p2,v2,...) constructs a flexible horizontal
    %  box and sets parameter p1 to value v1, etc.
    %
    %  A horizontal box lays out contents from left to right.  Users can
    %  resize contents by dragging the dividers.
    %
    %  See also: uix.VBoxFlex, uix.GridFlex, uix.HBox, uix.HButtonBox

    %  Copyright 2009-2016 The MathWorks, Inc.
    %  $Revision: 1682 $ $Date: 2018-06-11 16:57:09 +0100 (Mon, 11 Jun 2018) $

    properties( Access = public, Dependent, AbortSet )
        DividerMarkings % divider markings [on|off]
    end

    properties( Access = private )
        ColumnDividers = uix.Divider.empty( [0 1] ) % column dividers
        FrontDivider % front divider
        DividerMarkings_ = 'on' % backing for DividerMarkings
        MousePressListener = event.listener.empty( [0 0] ) % mouse press listener
        MouseReleaseListener = event.listener.empty( [0 0] ) % mouse release listener
        MouseMotionListener = event.listener.empty( [0 0] ) % mouse motion listener
        ActiveDivider = 0 % active divider index
        ActiveDividerPosition = [NaN NaN NaN NaN] % active divider position
        MousePressLocation = [NaN NaN] % mouse press location
        BackgroundColorListener % background color listener
    end

    methods

        function obj = HBoxFlex( varargin )
            %uix.HBoxFlex  Flexible horizontal box constructor
            %
            %  b = uix.HBoxFlex() constructs a flexible horizontal box.
            %
            %  b = uix.HBoxFlex(p1,v1,p2,v2,...) sets parameter p1 to value
            %  v1, etc.

            % Create front divider
            frontDivider = uix.Divider( 'Parent', obj, ...
                'Orientation', 'vertical', ...
                'BackgroundColor', obj.BackgroundColor * 0.75, ...
                'Visible', 'off' );

            % Create listeners
            backgroundColorListener = event.proplistener( obj, ...
                findprop( obj, 'BackgroundColor' ), 'PostSet', ...
                @obj.onBackgroundColorChange );

            % Store properties
            obj.FrontDivider = frontDivider;
            obj.BackgroundColorListener = backgroundColorListener;

            % Set Spacing property (may be overwritten by uix.set)
            obj.Spacing = 5;

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

        end % constructor

    end % structors

    methods

        function value = get.DividerMarkings( obj )

            value = obj.DividerMarkings_;

        end % get.DividerMarkings

        function set.DividerMarkings( obj, value )

            % Check
            assert( ischar( value ) && any( strcmp( value, {'on','off'} ) ), ...
                'uix:InvalidArgument', ...
                'Property ''DividerMarkings'' must be ''on'' or ''off'.' )

            % Set
            obj.DividerMarkings_ = value;

            % Mark as dirty
            obj.Dirty = true;

        end % set.DividerMarkings

    end % accessors

    methods( Access = protected )

        function onMousePress( obj, source, eventData )
            %onMousePress  Handler for WindowMousePress events

            % Check whether mouse is over a divider
            loc = find( obj.ColumnDividers.isMouseOver( eventData ) );
            if isempty( loc ), return, end

            % Capture state at button down
            divider = obj.ColumnDividers(loc);
            obj.ActiveDivider = loc;
            obj.ActiveDividerPosition = divider.Position;
            root = groot();
            obj.MousePressLocation = root.PointerLocation;

            % Make sure the pointer is appropriate
            obj.updateMousePointer( source, eventData );

            % Activate divider
            frontDivider = obj.FrontDivider;
            frontDivider.Position = divider.Position;
            divider.Visible = 'off';
            frontDivider.Parent = [];
            frontDivider.Parent = obj;
            frontDivider.Visible = 'on';

        end % onMousePress

        function onMouseRelease( obj, ~, ~ )
            %onMousePress  Handler for WindowMouseRelease events

            % Compute new positions
            loc = obj.ActiveDivider;
            if loc > 0
                root = groot();
                delta = root.PointerLocation(1) - obj.MousePressLocation(1);
                iw = loc;
                jw = loc + 1;
                ic = loc;
                jc = loc + 1;
                divider = obj.ColumnDividers(loc);
                contents = obj.Contents_;
                ip = uix.getPosition( contents(ic), 'pixels' );
                jp = uix.getPosition( contents(jc), 'pixels' );
                oldPixelWidths = [ip(3); jp(3)];
                minimumWidths = obj.MinimumWidths_(iw:jw,:);
                if delta < 0 % limit to minimum distance from left neighbor
                    delta = max( delta, minimumWidths(1) - oldPixelWidths(1) );
                else % limit to minimum distance from right neighbor
                    delta = min( delta, oldPixelWidths(2) - minimumWidths(2) );
                end
                oldWidths = obj.Widths_(iw:jw);
                newPixelWidths = oldPixelWidths + delta * [1;-1];
                if oldWidths(1) < 0 && oldWidths(2) < 0 % weight, weight
                    newWidths = oldWidths .* newPixelWidths ./ oldPixelWidths;
                elseif oldWidths(1) < 0 && oldWidths(2) >= 0 % weight, pixels
                    newWidths = [oldWidths(1) * newPixelWidths(1) / ...
                        oldPixelWidths(1); newPixelWidths(2)];
                elseif oldWidths(1) >= 0 && oldWidths(2) < 0 % pixels, weight
                    newWidths = [newPixelWidths(1); oldWidths(2) * ...
                        newPixelWidths(2) / oldPixelWidths(2)];
                else % sizes(1) >= 0 && sizes(2) >= 0 % pixels, pixels
                    newWidths = newPixelWidths;
                end
                obj.Widths_(iw:jw) = newWidths;
            else
                return
            end

            % Deactivate divider
            obj.FrontDivider.Visible = 'off';
            divider.Visible = 'on';

            % Reset state at button down
            obj.ActiveDivider = 0;
            obj.ActiveDividerPosition = [NaN NaN NaN NaN];
            obj.MousePressLocation = [NaN NaN];

            % Mark as dirty
            obj.Dirty = true;

        end % onMouseRelease

        function onMouseMotion( obj, source, eventData )
            %onMouseMotion  Handler for WindowMouseMotion events

            loc = obj.ActiveDivider;
            if loc == 0 % hovering, update pointer
                obj.updateMousePointer( source, eventData );
            else % dragging column divider
                root = groot();
                delta = root.PointerLocation(1) - obj.MousePressLocation(1);
                iw = loc;
                jw = loc + 1;
                ic = loc;
                jc = loc + 1;
                contents = obj.Contents_;
                ip = uix.getPosition( contents(ic), 'pixels' );
                jp = uix.getPosition( contents(jc), 'pixels' );
                oldPixelWidths = [ip(3); jp(3)];
                minimumWidths = obj.MinimumWidths_(iw:jw,:);
                if delta < 0 % limit to minimum distance from left neighbor
                    delta = max( delta, minimumWidths(1) - oldPixelWidths(1) );
                else % limit to minimum distance from right neighbor
                    delta = min( delta, oldPixelWidths(2) - minimumWidths(2) );
                end
                obj.FrontDivider.Position = ...
                    obj.ActiveDividerPosition + [delta 0 0 0];
            end

        end % onMouseMotion

        function onBackgroundColorChange( obj, ~, ~ )
            %onBackgroundColorChange  Handler for BackgroundColor changes

            backgroundColor = obj.BackgroundColor;
            highlightColor = min( [backgroundColor / 0.75; 1 1 1] );
            shadowColor = max( [backgroundColor * 0.75; 0 0 0] );
            columnDividers = obj.ColumnDividers;
            for jj = 1:numel( columnDividers )
                columnDivider = columnDividers(jj);
                columnDivider.BackgroundColor = backgroundColor;
                columnDivider.HighlightColor = highlightColor;
                columnDivider.ShadowColor = shadowColor;
            end
            frontDivider = obj.FrontDivider;
            frontDivider.BackgroundColor = shadowColor;

        end % onBackgroundColorChange

    end % event handlers

    methods( Access = protected )

        function redraw( obj )
            %redraw  Redraw contents
            %
            %  c.redraw() redraws the container c.

            % Call superclass method
            redraw@uix.HBox( obj )

            % Create or destroy column dividers
            b = numel( obj.ColumnDividers ); % current number of dividers
            c = max( [numel( obj.Widths_ )-1 0] ); % required number of dividers
            if b < c % create
                for ii = b+1:c
                    divider = uix.Divider( 'Parent', obj, ...
                        'Orientation', 'vertical', ...
                        'BackgroundColor', obj.BackgroundColor );
                    obj.ColumnDividers(ii,:) = divider;
                end
            elseif b > c % destroy
                % Destroy dividers
                delete( obj.ColumnDividers(c+1:b,:) )
                obj.ColumnDividers(c+1:b,:) = [];
                % Update pointer
                if c == 0 && strcmp( obj.Pointer, 'left' )
                    obj.unsetPointer()
                end
            end

            % Compute container bounds
            bounds = hgconvertunits( ancestor( obj, 'figure' ), ...
                [0 0 1 1], 'normalized', 'pixels', obj );

            % Retrieve size properties
            widths = obj.Widths_;
            minimumWidths = obj.MinimumWidths_;
            padding = obj.Padding_;
            spacing = obj.Spacing_;

            % Compute column divider positions
            xColumnSizes = uix.calcPixelSizes( bounds(3), widths, ...
                minimumWidths, padding, spacing );
            xColumnPositions = [cumsum( xColumnSizes(1:c,:) ) + padding + ...
                spacing * transpose( 0:c-1 ) + 1, repmat( spacing, [c 1] )];
            yColumnPositions = [padding + 1, max( bounds(4) - 2 * padding, 1 )];
            yColumnPositions = repmat( yColumnPositions, [c 1] );
            columnPositions = [xColumnPositions(:,1), yColumnPositions(:,1), ...
                xColumnPositions(:,2), yColumnPositions(:,2)];

            % Position column dividers
            for ii = 1:c
                columnDivider = obj.ColumnDividers(ii);
                columnDivider.Position = columnPositions(ii,:);
                switch obj.DividerMarkings_
                    case 'on'
                        columnDivider.Markings = columnPositions(ii,4)/2;
                    case 'off'
                        columnDivider.Markings = zeros( [0 1] );
                end
            end

        end % redraw

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

            % Update listeners
            if isempty( newFigure )
                mousePressListener = event.listener.empty( [0 0] );
                mouseReleaseListener = event.listener.empty( [0 0] );
                mouseMotionListener = event.listener.empty( [0 0] );
            else
                mousePressListener = event.listener( newFigure, ...
                    'WindowMousePress', @obj.onMousePress );
                mouseReleaseListener = event.listener( newFigure, ...
                    'WindowMouseRelease', @obj.onMouseRelease );
                mouseMotionListener = event.listener( newFigure, ...
                    'WindowMouseMotion', @obj.onMouseMotion );
            end
            obj.MousePressListener = mousePressListener;
            obj.MouseReleaseListener = mouseReleaseListener;
            obj.MouseMotionListener = mouseMotionListener;

            % Call superclass method
            reparent@uix.HBox( obj, oldFigure, newFigure )

            % Update pointer
            if ~isempty( oldFigure ) && ~strcmp( obj.Pointer, 'unset' )
                obj.unsetPointer()
            end

        end % reparent

    end % template methods

    methods( Access = protected )

        function updateMousePointer ( obj, source, eventData  )

            oldPointer = obj.Pointer;
            if any( obj.ColumnDividers.isMouseOver( eventData ) )
                newPointer = 'left';
            else
                newPointer = 'unset';
            end
            switch newPointer
                case oldPointer % no change
                    % do nothing
                case 'unset' % change, unset
                    obj.unsetPointer()
                otherwise % change, set
                    obj.setPointer( source, newPointer )
            end

        end % updateMousePointer

    end % helpers methods

end % classdef

author: goGPS-Project


text:
function [status, message]=xlwrite(filename,A,sheet, range)
% XLWRITE Write to Microsoft Excel spreadsheet file using Java
%   XLWRITE(FILE,ARRAY) writes ARRAY to the first worksheet in the Excel
%   file named FILE, starting at cell A1. It aims to have exactly the same
%   behaviour as XLSWRITE. See also XLSWRITE.
%
%   XLWRITE(FILE,ARRAY,SHEET) writes to the specified worksheet.
%
%   XLWRITE(FILE,ARRAY,RANGE) writes to the rectangular region
%   specified by RANGE in the first worksheet of the file. Specify RANGE
%   using the syntax 'C1:C2', where C1 and C2 are opposing corners of the
%   region.
%
%   XLWRITE(FILE,ARRAY,SHEET,RANGE) writes to the specified SHEET and
%   RANGE.
%
%   STATUS = XLWRITE(FILE,ARRAY,SHEET,RANGE) returns the completion
%   status of the write operation: TRUE (logical 1) for success, FALSE
%   (logical 0) for failure.  Inputs SHEET and RANGE are optional.
%
%   Input Arguments:
%
%   FILE    String that specifies the file to write. If the file does not
%           exist, XLWRITE creates a file, determining the format based on
%           the specified extension. To create a file compatible with Excel
%           97-2003 software, specify an extension of '.xls'. If you do not
%           specify an extension, XLWRITE applies '.xls'.
%   ARRAY   Two-dimensional logical, numeric or character array or, if each
%           cell contains a single element, a cell array.
%   SHEET   Worksheet to write. One of the following:
%           * String that contains the worksheet name.
%           * Positive, integer-valued scalar indicating the worksheet
%             index.
%           If SHEET does not exist, XLWRITE adds a new sheet at the end
%           of the worksheet collection.
%   RANGE   String that specifies a rectangular portion of the worksheet to
%           read. Not case sensitive. Use Excel A1 reference style.
%           * If you specify a SHEET, RANGE can either fit the size of
%             ARRAY or specify only the first cell (such as 'D2').
%           * If you do not specify a SHEET, RANGE must include both
%             corners and a colon character (:), even for a single cell
%             (such as 'D2:D2').
%           * If RANGE is larger than the size of ARRAY, Excel fills the
%             remainder of the region with #N/A. If RANGE is smaller than
%             the size of ARRAY, XLWRITE writes only the subset that fits
%             into RANGE to the file.
%
%   Note
%   * This function requires the POI library to be in your javapath.
%     To add the Apache POI Library execute commands:
%     (This assumes the POI lib files are in folder 'poi_library')
%       javaaddpath('poi_library/poi-3.8-20120326.jar');
%       javaaddpath('poi_library/poi-ooxml-3.8-20120326.jar');
%       javaaddpath('poi_library/poi-ooxml-schemas-3.8-20120326.jar');
%       javaaddpath('poi_library/xmlbeans-2.3.0.jar');
%       javaaddpath('poi_library/dom4j-1.6.1.jar');
%   * Excel converts Inf values to 65535. XLWRITE converts NaN values to
%     empty cells.
%
%   EXAMPLES
%   % Write a 7-element vector to testdata.xls:
%   xlwrite('testdata.xls', [12.7, 5.02, -98, 63.9, 0, -.2, 56])
%
%   % Write mixed text and numeric data to testdata2.xls
%   % starting at cell E1 of Sheet1:
%   d = {'Time','Temperature'; 12,98; 13,99; 14,97};
%   xlwrite('testdata2.xls', d, 1, 'E1')
%
%
%   REVISIONS
%   20121004 - First version using JExcelApi
%   20121101 - Modified to use POI library instead of JExcelApi (allows to
%           generate XLSX)
%   20121127 - Fixed bug: use existing rows if present, instead of
%           overwrite rows by default. Thanks to Dan & Jason.
%   20121204 - Fixed bug: if a numeric sheet is given & didn't exist,
%           an error was returned instead of creating the sheet. Thanks to Marianna
%   20130106 - Fixed bug: use existing cell if present, instead of
%           overwriting. This way original XLS formatting is kept & not
%           overwritten.
%   20130125 - Fixed bug & documentation. Incorrect working of NaN. Thanks Klaus
%   20130227 - Fixed bug when no sheet number given & added Stax to java
%               load. Thanks to Thierry
%
%   Copyright 2012-2013, Alec de Zegher
%==============================================================================

% Check if POI lib is loaded
if exist('org.apache.poi.ss.usermodel.WorkbookFactory', 'class') ~= 8 ...
    || exist('org.apache.poi.hssf.usermodel.HSSFWorkbook', 'class') ~= 8 ...
    || exist('org.apache.poi.xssf.usermodel.XSSFWorkbook', 'class') ~= 8

    error('xlWrite:poiLibsNotLoaded',...
        'The POI library is not loaded in Matlab.\nCheck that POI jar files are in Matlab Java path!');
end

% Import required POI Java Classes
import org.apache.poi.ss.usermodel.*;
import org.apache.poi.hssf.usermodel.*;
import org.apache.poi.xssf.usermodel.*;
import org.apache.poi.ss.usermodel.*;

import org.apache.poi.ss.util.*;

status=0;

% If no sheet & xlrange is defined, attribute an empty value to it
if nargin < 3; sheet = []; end
if nargin < 4; range = []; end

% Check if sheetvariable contains range data
if nargin < 4 && ~isempty(strfind(sheet,':'))
    range = sheet;
    sheet = [];
end

% check if input data is given
if isempty(A)
    error('xlwrite:EmptyInput', 'Input array is empty!');
end
% Check that input data is not bigger than 2D
if ndims(A) > 2
	error('xlwrite:InputDimension', ...
        'Dimension of input array should not be higher than two.');
end

% Set java path to same path as Matlab path
java.lang.System.setProperty('user.dir', pwd);

% Open a file
xlsFile = java.io.File(filename);

% If file does not exist create a new workbook
if xlsFile.isFile()
    % create XSSF or HSSF workbook from existing workbook
    fileIn = java.io.FileInputStream(xlsFile);
    xlsWorkbook = WorkbookFactory.create(fileIn);
else
    % Create a new workbook based on the extension.
    [~,~,fileExt] = fileparts(filename);

    % Check based on extension which type to create. If no (valid)
    % extension is given, create XLSX file
    switch lower(fileExt)
        case '.xls'
            xlsWorkbook = HSSFWorkbook();
        case '.xlsx'
            xlsWorkbook = XSSFWorkbook();
        otherwise
            xlsWorkbook = XSSFWorkbook();

            % Also update filename with added extension
            filename = [filename '.xlsx'];
    end
end

% If sheetname given, enter data in this sheet
if ~isempty(sheet)
    if isnumeric(sheet)
        % Java uses 0-indexing, so take sheetnumer-1
        % Check if the sheet can exist
        if xlsWorkbook.getNumberOfSheets() >= sheet && sheet >= 1
            xlsSheet = xlsWorkbook.getSheetAt(sheet-1);
        else
            % There are less number of sheets, that the requested sheet, so
            % return an empty sheet
            xlsSheet = [];
        end
    else
        xlsSheet = xlsWorkbook.getSheet(sheet);
    end

    % Create a new sheet if it is empty
    if isempty(xlsSheet)
        warning('xlwrite:AddSheet', 'Added specified worksheet.');

        % Add the sheet
        if isnumeric(sheet)
            xlsSheet = xlsWorkbook.createSheet(['Sheet ' num2str(sheet)]);
        else
            % Create a safe sheet name
            sheet = WorkbookUtil.createSafeSheetName(sheet);
            xlsSheet = xlsWorkbook.createSheet(sheet);
        end
    end

else
    % check number of sheets
    nSheets = xlsWorkbook.getNumberOfSheets();

    % If no sheets, create one
    if nSheets < 1
        xlsSheet = xlsWorkbook.createSheet('Sheet 1');
    else
        % Select the first sheet
        xlsSheet = xlsWorkbook.getSheetAt(0);
    end
end

% if range is not specified take start row & col at A1
% locations are 0 indexed
if isempty(range)
    iRowStart = 0;
    iColStart = 0;
    iRowEnd = size(A, 1)-1;
    iColEnd = size(A, 2)-1;
else
    % Split range in start & end cell
    iSeperator = strfind(range, ':');
    if isempty(iSeperator)
        % Only start was defined as range
        % Create a helper to get the row and column
        cellStart = CellReference(range);
        iRowStart = cellStart.getRow();
        iColStart = cellStart.getCol();
        % End column calculated based on size of A
        iRowEnd = iRowStart + size(A, 1)-1;
        iColEnd = iColStart + size(A, 2)-1;
    else
        % Define start & end cell
        cellStart = range(1:iSeperator-1);
        cellEnd = range(iSeperator+1:end);

        % Create a helper to get the row and column
        cellStart = CellReference(cellStart);
        cellEnd = CellReference(cellEnd);

        % Get start & end locations
        iRowStart = cellStart.getRow();
        iColStart = cellStart.getCol();
        iRowEnd = cellEnd.getRow();
        iColEnd = cellEnd.getCol();
    end
end

% Get number of elements in A (0-indexed)
nRowA = size(A, 1)-1;
nColA = size(A, 2)-1;

% If data is a cell, convert it
if ~iscell(A)
    A = num2cell(A);
end

% Iterate over all data
for iRow = iRowStart:iRowEnd
    % Fetch the row (if it exists)
    currentRow = xlsSheet.getRow(iRow);
    if isempty(currentRow)
        % Create a new row, as it does not exist yet
        currentRow = xlsSheet.createRow(iRow);
    end

    % enter data for all cols
    for iCol = iColStart:iColEnd
        % Check if cell exists
        currentCell = currentRow.getCell(iCol);
        if isempty(currentCell)
            % Create a new cell, as it does not exist yet
            currentCell = currentRow.createCell(iCol);
        end

        % Check if we are still in array A
        if (iRow-iRowStart)<=nRowA && (iCol-iColStart)<=nColA
            % Fetch the data
            data = A{iRow-iRowStart+1, iCol-iColStart+1};

            if ~isempty(data)
                % if it is a NaN value, convert it to an empty string
                if isnumeric(data) && isnan(data)
                    data = '';
                end

                % Write data to cell
                currentCell.setCellValue(data);
            end

        else
            % Set field to NA
            currentCell.setCellErrorValue(FormulaError.NA.getCode());
        end
    end
end

% Write & close the workbook
fileOut = java.io.FileOutputStream(filename);
xlsWorkbook.write(fileOut);
fileOut.close();

status = 1;

end

author: Sable

repo: mcbench-benchmarks


text:
function ns3de_test02 ( )

%*****************************************************************************80
%
%% NS3DE_TEST02 samples the residual at the initial time.
%
%  Location:
%
%    http://people.sc.fsu.edu/~jburkardt/m_src/navier_stokes_3d_exact/ns3de_test02.m
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 January 2015
%
%  Author:
%
%    John Burkardt
%
  a = pi / 4.0;
  d = pi / 2.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'NS3DE_TEST02\n' );
  fprintf ( 1, '  Sample the Navier-Stokes residuals' );
  fprintf ( 1, '  at the initial time T = 0, using a region that is \n' );
  fprintf ( 1, '  the cube centered at (0,0,0) with "radius" 1.0,\n' );
  fprintf ( 1, '  Parameter A = %g\n', a );
  fprintf ( 1, '  Parameter D = %g\n', d );

  n = 1000;
  x = 2.0 * rand ( n, 1 ) - 1.0;
  y = 2.0 * rand ( n, 1 ) - 1.0;
  z = 2.0 * rand ( n, 1 ) - 1.0;
  t = 0.0;

  [ ur, vr, wr, pr ] = resid_ethier ( a, d, n, x, y, z, t );

  fprintf ( 1, '\n' );
  fprintf ( 1, '           Minimum       Maximum\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Ur:  %14.6g  %14.6g\n', min ( abs ( ur ) ), max ( abs ( ur ) ) );
  fprintf ( 1, '  Vr:  %14.6g  %14.6g\n', min ( abs ( vr ) ), max ( abs ( vr ) ) );
  fprintf ( 1, '  Wr:  %14.6g  %14.6g\n', min ( abs ( wr ) ), max ( abs ( wr ) ) );
  fprintf ( 1, '  Pr:  %14.6g  %14.6g\n', min ( abs ( pr ) ), max ( abs ( pr ) ) );

  return
end


text:
function cumsum = cfg_example_cumsum1
% Example script that creates an cfg_exbranch to sum two numbers. The
% inputs are entered as vector, the output is a vector containing the
% cumulative sums. This function differs from cfg_example_sum (except from
% names) only in the specification of the output subscript.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_example_cumsum1.m 1716 2008-05-23 08:18:45Z volkmar $

rev = '$Rev: 1716 $'; %#ok

%% Input Item
input1         = cfg_entry; % This is the generic data entry item
input1.name    = 'Input a'; % The displayed name
input1.tag     = 'a';       % The name appearing in the harvested job structure. This name must be unique among all items in the val field of the superior node
input1.strtype = 'e';       % No restriction on what type of data is entered. This could be used to restrict input to real numbers, integers ...
input1.num     = [1 Inf];   % Number of inputs required (2D-array with exactly one row and two column)
input1.help    = {'This is the input vector.','The elements will be added together.'}; % help text displayed

%% Executable Branch
cumsum      = cfg_exbranch;       % This is the branch that has information about how to run this module
cumsum.name = 'cumsum1';             % The display name
cumsum.tag  = 'cfg_example_cumsum1'; % The name appearing in the harvested job structure. This name must be unique among all items in the val field of the superior node
cumsum.val  = {input1};    % The items that belong to this branch. All items must be filled before this branch can run or produce virtual outputs
cumsum.prog = @cfg_example_run_cumsum1;  % A function handle that will be called with the harvested job to run the computation
cumsum.vout = @cfg_example_vout_cumsum1; % A function handle that will be called with the harvested job to determine virtual outputs
cumsum.help = {'Compute the cumulative sum of two numbers.'};

%% Local Functions
% The cfg_example_vout_cumsum1 function can go here, it is not useful outside
% the batch environment.
function vout = cfg_example_vout_cumsum1(job)
% Determine what outputs will be present if this job is run. In this case,
% the structure of the inputs is fixed, and the output is always a single
% number. Note that input items may not be numbers, they can also be
% dependencies.

vout = cfg_dep;                        % The dependency object
vout.sname      = 'cumsum(a)';       % Displayed dependency name
vout.src_output = substruct('.','cs'); % The output subscript reference. The length of the output vector depends on the unknown length of the input vector. Therefore, the vector output needs to be assigned to a struct field.

author: spm