text:
function polygon_properties_test02 ( )

%*****************************************************************************80
%
%% POLYGON_PROPERTIES_TEST02 tests POLYGON_AREA.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2014
%
%  Author:
%
%    John Burkardt
%
  test_num = 2;
  area_exact_test = [ 2.0, 6.0 ];
  n_test = [ 4, 8 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'POLYGON_PROPERTIES_TEST02\n' );
  fprintf ( 1, '  For a polygon:\n' );
  fprintf ( 1, '  POLYGON_AREA computes the area.\n' );

  for test = 1 : test_num

    n = n_test(test);
    area_exact = area_exact_test(test);

    if ( test == 1 )

      v = [ ...
        1.0, 0.0; ...
        2.0, 1.0; ...
        1.0, 2.0; ...
        0.0, 1.0 ]';

    elseif ( test == 2 )

      v = [ ...
        0.0, 0.0; ...
        3.0, 0.0; ...
        3.0, 3.0; ...
        2.0, 3.0; ...
        2.0, 1.0; ...
        1.0, 1.0; ...
        1.0, 2.0; ...
        0.0, 2.0 ]';

    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Number of polygonal vertices = %d\n', n );

    r8mat_transpose_print ( 2, n, v, '  The polygon vertices:' );

    area = polygon_area ( n, v );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Exact area is        %g\n', area_exact );
    fprintf ( 1, '  The computed area is %g\n', area );

  end

  return
end


text:
function res = badsamples(this, chanind, sampind, trialind)
% Returns an array of 0/1 marking bad data based on artefact events and bad flags
% FORMAT res = badsamples(this, chanind, sampind, trialind)
% _______________________________________________________________________
% Copyright (C) 2013 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: badsamples.m 7199 2017-11-01 16:42:12Z vladimir $

if ischar(chanind) && isequal(chanind, ':')
    chanind = 1:nchannels(this);
end

if ischar(sampind) && isequal(sampind, ':')
    sampind = 1:nsamples(this);
end

if ischar(trialind) && isequal(trialind, ':')
    trialind = 1:ntrials(this);
end

if ~isequal(type(this), 'continuous') && ~any(trialonset(this))
    error('Trial onset information is not available. Cannot map artefact events to samples.');
end

res = false(length(chanind), nsamples(this), length(trialind));
for i = 1:length(trialind)

    ev = events(this, trialind(i));
    if iscell(ev)
        ev = ev{1};
    end

    if ~isempty(ev)
        ev = ev(intersect(strmatch('artefact', {ev.type}),...
            find(cellfun(@ischar, {ev.value}))));
        for k = 1:numel(ev)
            [dum, chan] = intersect(chanind, selectchannels(this, ev(k).value));
            samples     = find((trialonset(this, trialind(i))+time(this))>=ev(k).time & ...
                (trialonset(this, trialind(i))+time(this))<=(ev(k).time+ev(k).duration));
            res(chan, samples, i) = true;
        end
    end
end

res = res(:, sampind, :);
res(badchannels(this, chanind), :, :) = true;
res(:, :, badtrials(this, trialind))  = true;

author: spm


text:
function [CnoHelp] = SkipMatlabHelp (C)
%
%   [CnoHelp] = SkipMatlabHelp (C)
%
%Purpose:
%   removes the first set of comments from a character array
%   that would normally consiture matlab's help text (see example)
%
%
%Input Parameters:
%   C : a character array
%
%
%Output Parameters:
%   CnoHelp : C without the help comments
%
%
%
%More Info :
%   C = sprintf ('%%Hello Baby\n%%Couci Couci Cooo\nAbc = 100\n%%Waldo\njon=6')
%   [CnoHelp] = SkipMatlabHelp (C)
%
%  PurgeComments
%	NextString
%	FindChar
%
%     Author : Ziad Saad
%     Date : Sat Mar 27 14:27:36 CST 1999

%Define the function name for easy referencing
FuncName = 'SkipMatlabHelp';

%Debug Flag
DBG = 1;

%initailize return variables
err = 1;
CnoHelp = [];

if (C(1) == '%'),
	cpos = 1;
	Loc = 1;
	while (~isempty(Loc) & Loc(1) == 1 & ~isempty(C)),
		%get the 1st line
			[err, Cnext, cend] = NextString (C, 'NewLine', 1);
		%get ridd of 1st line
			nC = length(C);
			C = C(cend+1:nC);
		%search for first %
			[err, Loc] = FindChar (C, '%');
	end
	CnoHelp = C;
else
	CnoHelp = C;
end

err = 0;
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [Jul1,Jul2]=Cal2TAI(year,month,day,hour,minute,second)
%%CAL2TT  Convert dates in terms of the Gregorian calendar in years,
%          months, days, hours, minutes and seconds with the time in
%          universal coordinated time (UTC) to a two-part Julian
%          date in international atomic time (TAI).
%
%INPUTS:   year    A matrix of integer  years in the Gregorian
%                  calendar under UTC time.
%          month   A matrix of integer months in the Gregorian calendar
%                  under UTC time. 1<=month<=12
%          day     A matrix of integer days in the Gregorian calendar
%                  under UTC time. Days count from 1.
%          hour    A matrix of integer hour under the Gregorian calendar.
%                  UTC time 0<=hour<=23
%          minute  A matrix of integer minutes in the Gregorian calendar
%                  under UTC time. 0<=minute<=59.
%          second  A matrix of double floating point second in the
%                  Gregorian calendar under UTC time. This is normally
%                  less than 60, but can be a value less than 61 or 59 at
%                  the right hour on a day with a leap second.
%
%OUTPUTS:   Jul1, Jul2  Matrices of the time as pseudo-Julian dates in TAI
%                       where each row/column corresponds to the values in
%                       the same row/column of the input matrices.
%
%This function just calls the function Cal2UTC and then UTC2TAI.
%
%March 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

[Jul1,Jul2]=Cal2UTC(year,month,day,hour,minute,second);
[Jul1,Jul2]=UTC2TAI(Jul1,Jul2);
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
clear all
        % Desemnarea unor variabile simbolice
syms x u v t
        % Generarea unei functii de o variabila
f=sin(x)^3;
        % Deschiderea unei ferestre grafice noi
figure(1)
        % Reprezentarea functiei f intre limitele implicite
ezplot(f)
        % Reprezentarea functiei f intre limitele [0, 8*pi]
figure(2)
ezplot(f,[0,8*pi])
        % Generarea unei functii de doua variabile
g=u^2-v^2+1;
        % Reprezentarea functiei g de doua variabile intre limitele date
figure(3)
ezplot(g,[-2,2,-5,5])
        % Generarea unei functii parametrice
x=sin(t);
y=cos(t);
        % Reprezentarea functiei parametrice intre limitele impuse
figure(4)
ezplot(x,y)

author: Sable

repo: mcbench-benchmarks


text:
%% Example of using KST class for interfacing with KUKA iiwa robots

% An example script, it is used to show how to use the different
% functions of the KUKA Sunrise Toolbox

% First start the server on the KUKA iiwa controller
% Then run the following script in Matlab

% This example works with Sunrise application version KST_1.7  and higher.
% Copyright Mohammad SAFEEA, 26th-June-2018

close all;clear;clc;
warning('off')
%% Create the robot object
ip='172.31.1.147'; % The IP of the controller
arg1=KST.LBR7R800; % choose the robot iiwa7R800 or iiwa14R820
arg2=KST.Medien_Flansch_elektrisch; % choose the type of flange
Tef_flange=eye(4); % transofrm matrix of EEF with respect to flange
iiwa=KST(ip,arg1,arg2,Tef_flange); % create the object

%% Start a connection with the server
flag=iiwa.net_establishConnection();
if flag==0
    return;
end
pause(1);
disp('Doing some stuff')

%% move to initial position
jPos={0,pi*20/180,0,-pi*70/180,0,pi*60/180,0}; % initial confuguration
relVel=0.15; % relative velocity
iiwa.movePTPJointSpace(jPos, relVel); % point to point motion in joint space

%% Linear relative motion of end effector, relative to base frame
deltaX=0;deltaY=0;deltaZ=100.;
Pos{1}=deltaX;
Pos{2}=deltaY;
Pos{3}=deltaZ;

vel=250; % linear velocity of end effector, mm/sec
iiwa.movePTPLineEefRelBase(Pos, vel);
Pos{3}=-deltaZ;
iiwa.movePTPLineEefRelBase(Pos, vel);

%% Linear relative motion of end effector, relative to EEF frame
Pos{1}=deltaX;
Pos{2}=deltaY;
Pos{3}=deltaZ;

vel=150; % linear velocity of end effector, mm/sec

iiwa.movePTPLineEefRelEef(Pos, vel);
Pos{3}=-deltaZ;
iiwa.movePTPLineEefRelEef(Pos, vel);

%% turn off the server
iiwa.net_turnOffServer();
warning('on');

repo: KST-Kuka-Sunrise-Toolbox


text:
function [Xgp Ygp]=projectToGroundPlane(Xi, Yi, sceneInfo)
%
% (C) Anton Andriyenko, 2012
%
% The code may be used free of charge for non-commercial and
% educational purposes, the only requirement is that this text is
% preserved within the derivative work. For any other purpose you
% must contact the authors for permission. This code may not be
% redistributed without written permission from the authors.

[F N]=size(Xi);
Xgp=zeros(size(Xi));
Ygp=zeros(size(Xi));

for t=1:F
    extar=find(Xi(t,:));
    for id=extar
        [Xgp(t,id) Ygp(t,id) zw]=imageToWorld(Xi(t,id), Yi(t,id), sceneInfo.camPar);
    end
end

end

author: VisDrone

repo: DroneCrowd

save_path: github-repos/MATLAB/VisDrone-DroneCrowd


text:
function G = slpwmetricgraph(X, varargin)
%SLPWMETRICGRAPH Constructs a graph based on pairwise metrics
%
% $ Syntax $
%   - G = slpwmetricgraph(X, ...)
%   - G = slpwmetricgraph(X, Xt, ...)
%
% $ Arguments $
%   - X:        The sample matrix with each column as a (source) node
%   - Xt:       The sample matrix with each column as a (target) node
%   - G:        The constructed graph
%
% $ Description $
%   - G = slpwmetricgraph(X, ...) constructs a graph based on computation
%     of pairwise metric between vector samples. You can specify the
%     following properties:
%     \*
%     \t   The Properties of Graph Matrix construction           \\
%     \h      name       &      description
%            'sparse'    & whether the target graph G is sparse
%                          (default = true)
%            'valtype'   & The type of values in G: 'logical'|'numeric'
%                          (default = 'numeric')
%                          The value output by evalfunctor should conform
%                          to the specified valtype
%            'maxblk'    & The maximum number of elements that can be
%                          computed in each batch. (default = 1e7)
%            'mfunctor'  & The functor to compute the metrics. In the form:
%                          V = f(X1, X2, ...)
%                          default = {@slmetric_pw, 'eucdist')
%            'tfunctor'  & The functor to transform the metric values,
%                          like tv = f(sv, ...)
%                          The transform is taken before the threshold
%                          filtering. default = []
%            'thres'     & The threshold values. The values not in the
%                          valid range set by the threshold is regarded as
%                          zeros.
%                          (default = [], means not to use threshold)
%            'threstype' & The type of threshold value, (default = 'ub')
%                           - 'lb': a lower bound scalar
%                                   valid_range: x >= thres
%                           - 'ub': a upper bound scalar
%                                   valid_range: x <= thres
%                           - 'lub': a vector of lower and upper bounds
%                                    valid_range: thres(1) <= x <= thres(2)
%     \*
%
%   - G = slpwmetricgraph(X, Xt, ...) constructs a bigraph with different
%     source and target node set. The properties mentioned above are
%     all applicable here.
%
% $ Remarks $
%   - The implementation is based on slpwgraph.
%
% $ History $
%   - Created by Dahua Lin, on Sep 8th, 2006
%   - Modified by Dahua Lin, on Sep 10th, 2006
%       - Base on upgraded slpwgraph function
%       - Add support for bigraph
%

%% parse and verify input

if nargin >= 2
    if isnumeric(varargin{1})
        Xt = varargin{1};
        if nargin == 2
            params = {};
        else
            params = varargin(2:end);
        end
    else
        Xt = X;
        params = varargin;
    end
else
    Xt = X;
    params = {};
end

opts.sparse = true;
opts.valtype = 'numeric';
opts.maxblk = 1e7;
opts.mfunctor = {@slmetric_pw, 'eucdist'};
opts.tfunctor = [];
opts.thres = [];
opts.threstype = 'ub';
opts = slparseprops(opts, params{:});

if ~ismember(opts.threstype, {'lb', 'ub', 'lub'})
    error('sltoolbox:invalidarg', ...
        'Invalid type of threshold: %s', opts.threstype);
end

%% Main skeleton

n = size(X, 2);
nt = size(Xt, 2);
mfunctor = opts.mfunctor;
tfunctor = opts.tfunctor;
filter = getvaluefilter(opts);
evalfunctor = {@metricevalfunc, mfunctor, tfunctor, filter};

G = slpwgraph(X, Xt, n, nt, evalfunctor, ...
    'sparse', opts.sparse, ...
    'valtype', opts.valtype, ...
    'maxblk', opts.maxblk);

%% Core functions

function V = metricevalfunc(X, Xt, inds1, inds2, mfunctor, tfunctor, filter)

X1 = X(:, inds1);
X2 = Xt(:, inds2);
V = slevalfunctor(mfunctor, X1, X2);

if ~isempty(tfunctor)
    V = slevalfunctor(tfunctor, V);
end

if ~isempty(filter)
    zero_range = ~filter(V);
    V(zero_range) = 0;
end

function filter = getvaluefilter(opts)

if isempty(opts.thres)
    filter = [];
else
    thr = opts.thres;
    switch opts.threstype
        case 'lb'
            if ~isscalar(thr)
                error('sltoolbox:invalidarg', 'For lb type, the thres should be a scalar');
            end
            filter = @(x) x >= thr;
        case 'ub'
            if ~isscalar(thr)
                error('sltoolbox:invalidarg', 'For ub type, the thres should be a scalar');
            end
            filter = @(x) x <= thr;
        case 'lub'
            if ~isvector(thr) || length(thr) ~= 2
                error('sltoolbox:invalidarg', 'For lub type, the thres should be a length-2 vector');
            end
            filter = @(x) x >= thr(1) & x <= thr(2);
    end
end

author: lmthang


text:
%% Copyright (C) 2016 Abhinav Tripathi
%% Copyright (C) 2016, 2018-2019, 2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym chebyshevU (@var{n}, @var{x})
%% Find the nth symbolic Chebyshev polynomial of the second kind.
%%
%% If @var{n} is a vector then it returns a vector with Chebyshev polynomials
%% of the second kind for each element of @var{n}.
%%
%% Examples:
%% @example
%% @group
%% syms x
%% chebyshevU(1, x)
%%   @result{} (sym) 2⋅x
%% chebyshevU(2, x)
%%   @result{} (sym)
%%          2
%%       4⋅x  - 1
%% syms n
%% chebyshevU(n, x)
%%   @result{} (sym) chebyshevu(n, x)
%% @end group
%% @end example
%%
%% The inputs can be vectors, for example:
%% @example
%% @group
%% syms x
%% chebyshevU([0 1 2], x)
%%   @result{} (sym 1×3 matrix)
%%       ⎡           2    ⎤
%%       ⎣1  2⋅x  4⋅x  - 1⎦
%% @end group
%% @end example
%%
%% @seealso{@@sym/chebyshevT, @@double/chebyshevU}
%% @end defmethod

function y = chebyshevU(n, x)
  if (nargin ~= 2)
    print_usage ();
  end
  y = elementwise_op ('chebyshevu', sym(n), sym(x));
end

%!error chebyshevU (sym(1))
%!error chebyshevU (sym(1), 2, 3)

%!assert (isequaln (chebyshevU (2, sym(nan)), sym(nan)))

%!shared x
%! syms x

%!assert(isequal(chebyshevU(0, x), sym(1)))
%!assert(isequal(chebyshevU(1, x), 2*x))
%!assert(isequal(chebyshevU(2, x), 4*x*x - 1))
%!assert(isequal(chebyshevU([0 1 2], x), [sym(1) 2*x (4*x*x-1)]))

%!test
%! % round trip
%! syms n z
%! f = chebyshevU (n, z);
%! h = function_handle (f, 'vars', [n z]);
%! A = h (1.1, 2.2);
%! B = chebyshevU (1.1, 2.2);
%! assert (A, B)

repo: octsympy


text:
% Script demonstrating usage of the cbpdndlcns function.
%
% Author: Brendt Wohlberg <brendt@lanl.gov>
%
% This file is part of the SPORCO library. Details of the copyright
% and user license can be found in the 'Copyright' and 'License' files
% distributed with the library.

% Training images
S0 = zeros(512, 512, 5, 'single');
S0(:,:,1) = single(stdimage('lena.grey')) / 255;
S0(:,:,2) = single(stdimage('barbara.grey')) / 255;
S0(:,:,3) = single(stdimage('kiel.grey')) / 255;
S0(:,:,4) = single(rgb2gray(stdimage('mandrill'))) / 255;
tmp = single(stdimage('man.grey')) / 255;
S0(:,:,5) = tmp(101:612, 101:612);

%Reduce images size to speed up demo script
tmp = zeros(256, 256, 5, 'single');
for k = 1:size(S0,3),
  tmp(:,:,k) = imresize(S0(:,:,k), 0.5);
end
S0 = tmp;

% Filter input images and compute highpass images
npd = 16;
fltlmbd = 5;
[Sl, Sh] = lowpass(S0, fltlmbd, npd);

% Construct initial dictionary
D0 = zeros(8,8,32, 'single');
D0(3:6,3:6,:) = single(randn(4,4,32));

% Set up cbpdndl parameters
lambda = 0.2;
opt = [];
opt.Verbose = 1;
opt.MaxMainIter = 250;
opt.rho = 50*lambda + 0.5;
opt.sigma = size(Sh,3);
opt.AutoRho = 1;
opt.AutoRhoPeriod = 10;
opt.AutoSigma = 1;
opt.AutoSigmaPeriod = 10;
opt.XRelaxParam = 1.8;
opt.DRelaxParam = 1.8;

% Do dictionary learning
[D, X, optinf] = cbpdndlcns(D0, Sh, lambda, opt);

% Display learned dictionary
figure;
imdisp(tiledict(D));

% Plot functional value evolution
figure;
plot(optinf.itstat(:,2));
xlabel('Iterations');
ylabel('Functional value');

author: thfylsty


%clear;
clear;
pp = '/home/yutian/re-id/dataset/market1501/bounding_box_test/*.jpg';
%% change ID-level attribute to instance-level attribute
load('market_attribute.mat');
file = dir(pp);
counter_data=1;
counter_last = 1;
class = 0;
c_last = '';
cc = [];
for i=1:length(file)
    c = strsplit(file(i).name,'_');
    if(~isequal(c{1},c_last))
        a = str2num(c{1});
        if a>0
            class = class + 1;
        end
        fprintf('%d::%d\n',class,counter_data-counter_last);
        cc=[cc;counter_data-counter_last];
        c_last = c{1};
        counter_last = counter_data;
    end
    if a>0
        dtupblack(counter_data) = market_attribute.test.upblack(class);
        dtupwhite(counter_data) = market_attribute.test.upwhite(class);
        dtupred(counter_data) = market_attribute.test.upred(class);
        dtuppurple(counter_data) = market_attribute.test.uppurple(class);
        dtupyellow(counter_data) = market_attribute.test.upyellow(class);
        dtupblue(counter_data) = market_attribute.test.upblue(class);
        dtupgreen(counter_data) = market_attribute.test.upgreen(class);
        dtupgray(counter_data) = market_attribute.test.upgray(class);

        dtdownblack(counter_data) = market_attribute.test.downblack(class);
        dtdownwhite(counter_data) = market_attribute.test.downwhite(class);
        dtdownpink(counter_data) = market_attribute.test.downpink(class);
        dtdownpurple(counter_data) = market_attribute.test.downpurple(class);
        dtdownyellow(counter_data) = market_attribute.test.downyellow(class);
        dtdownblue(counter_data) = market_attribute.test.downblue(class);
        dtdowngreen(counter_data) = market_attribute.test.downgreen(class);
        dtdowngray(counter_data) = market_attribute.test.downgray(class);
        dtdownbrown(counter_data) = market_attribute.test.downbrown(class);

        dtup(counter_data) = market_attribute.test.up(class);
        dtgender(counter_data) = market_attribute.test.gender(class);
        dtdown(counter_data) = market_attribute.test.down(class);
        dtclothes(counter_data) = market_attribute.test.clothes(class);
        dthat(counter_data) = market_attribute.test.hat(class);
        dtbackpack(counter_data) = market_attribute.test.backpack(class);
        dthandbag(counter_data) = market_attribute.test.handbag(class);
        dtbag(counter_data) = market_attribute.test.bag(class);
        dtage(counter_data) = market_attribute.test.age(class);
        dthair(counter_data) = market_attribute.test.hair(class);
        counter_data = counter_data + 1;
    end
end
%% evaluate
load('gallery_market.mat')
temp1 = 0;temp2 = 0;temp3 = 0;temp4 = 0;temp5 = 0;temp6 = 0;temp7 = 0;temp8 = 0;
temp9 = 0; temp10 = 0;temp11 = 0; temp12 = 0;
for j = 1:13115
    if (gallery(j,1)==dtgender(j))
        temp1 = temp1+1;
    end
    if (gallery(j,2)==dthair(j))
        temp3 = temp3+1;
    end
    if (gallery(j,3)==dtage(j))
        temp2 = temp2+1;
    end
    if (gallery(j,4)==dtup(j))
        temp4 = temp4+1;
    end
    if (gallery(j,5)==dtdown(j))
        temp5 = temp5+1;
    end
    if (gallery(j,6)==dtclothes(j))
        temp6 = temp6+1;
    end
    if (gallery(j,7)==dthat(j))
        temp10 = temp10+1;
    end
    if (gallery(j,8)==dtbackpack(j))
        temp7 = temp7+1;
    end
    if (gallery(j,9)==dtbag(j))
        temp9 = temp9+1;
    end
    if (gallery(j,10)==dthandbag(j))
        temp8 = temp8+1;
    end
    if(gallery(j,11)==1 && dtupblack(j)==2)
        temp11 = temp11+1;
    end
    if(gallery(j,11)==2 && dtupwhite(j)==2)
        temp11 = temp11+1;
    end
    if(gallery(j,11)==3 && dtupred(j)==2)
        temp11 = temp11+1;
    end
    if(gallery(j,11)==4 && dtuppurple(j)==2)
        temp11 = temp11+1;
    end
    if(gallery(j,11)==5 && dtupgray(j)==2)
        temp11 = temp11+1;
    end
    if(gallery(j,11)==6 && dtupblue(j)==2)
        temp11 = temp11+1;
    end
    if(gallery(j,11)==7 && dtupgreen(j)==2)
        temp11 = temp11+1;
    end
    if(gallery(j,11)==8 && dtupyellow(j)==2)
        temp11 = temp11+1;
    end
    if(gallery(j,12)==1 && dtdownblack(j)==2)
        temp12 = temp12+1;
    end
    if(gallery(j,12)==2 && dtdownwhite(j)==2)
        temp12 = temp12+1;
    end
    if(gallery(j,12)==3 && dtdownpink(j)==2)
        temp12 = temp12+1;
    end
    if(gallery(j,12)==4 && dtdowngray(j)==2)
        temp12 = temp12+1;
    end
    if(gallery(j,12)==5 && dtdownblue(j)==2)
        temp12 = temp12+1;
    end
    if(gallery(j,12)==6 && dtdowngreen(j)==2)
        temp12 = temp12+1;
    end
    if(gallery(j,12)==7 && dtdownbrown(j)==2)
        temp12 = temp12+1;
    end
    if(gallery(j,12)==8 && dtdownyellow(j)==2)
        temp12 = temp12+1;
    end
    if(gallery(j,12)==9 && dtdownpurple(j)==2)
        temp12 = temp12+1;
    end
end

precision1 = temp1 /j
precision2 = temp2 /j
precision3 = temp3 /j
precision4 = temp4 /j
precision5 = temp5 /j
precision6 = temp6 /j
precision7 = temp7 /j
precision8 = temp8 /j
precision9 = temp9 /j
precision10 = temp10 /j
precision11 = temp11 /j
precision12 = temp12 /j
average = (precision1+precision2+precision3+precision4+precision5+precision6+precision7+precision8+precision9+precision10+precision11+precision12)/12


text:
% close all; clc;

processStatus = [processStatus, newline, ...
                     datestr(now), ' >> ', '[Start text outputing].'];

targetFile = 'in.txt';
targetFile = strcat(inputFilesFolder, '\', targetFile);

processStatus = [processStatus, newline, ...
                     datestr(now), ' >> ', 'Reading text file ...'];

% Read the contents an array
fid = fopen(targetFile);
linesFromFile = fread(fid, '*char')';
fclose(fid);

% Read the contents an array in ascii dec values
fid = fopen(targetFile);
linesFromFileInAscii = fread(fid, '*uint8')';
fclose(fid);

% '\r\t' == [13   10]
newlineInAscii1 = [13 10]; % carriage return + newline (used in most text files)
% '\r\t' unaddressed yet
newlineInAscii2 = [13 9]; % newlineChar = '\r\t';

% scan for both to be sure
locationOfCarriageReturnChar = strfind(linesFromFileInAscii, newlineInAscii1(1));
% locationOfHorizontalTabChar = strfind(linesFromFileInAscii, newlineInAscii1(2));

fileContents = [];
startPos = 1;

for i=1:length(locationOfCarriageReturnChar)
    if(linesFromFileInAscii(locationOfCarriageReturnChar(i)+1) == newlineInAscii1(2))
        fileContents = [fileContents {linesFromFile(startPos:(locationOfCarriageReturnChar(i)-1))}];
        startPos = locationOfCarriageReturnChar(i)+2; % skip over 2 character
    end
end

fileContents = [fileContents {linesFromFile(startPos:end)}];

% locationOfCarriageReturnChar = strfind(linesFromFileInAscii, newlineInAscii2(1));
% locationOfHorizontalTabChar = strfind(linesFromFileInAscii, newlineInAscii2(2));

lineCount = 0;

for i = 1:length(fileContents)
   if (length(char(fileContents(i))) ~= 0)
        lineCount = lineCount + 1;
        nonEmptyLines(lineCount) = fileContents(i);
   end
end

nonEmptyLinesInAscii = [];

nonEmptyLinesInAscii = uint8(char(nonEmptyLines(1)));

for i = 2:lineCount
    lineToBeAddressed = uint8(char(nonEmptyLines(i)));
    nonEmptyLinesInAscii = [nonEmptyLinesInAscii newlineInAscii1 lineToBeAddressed]; % fit in the newlineInAscii1
end

targetFile = 'out.txt';
targetFile = strcat(outputFilesFolder, '\', targetFile);

processStatus = [processStatus, newline, ...
                     datestr(now), ' >> ', 'Writing text file ...'];

% write to file
fid = fopen(targetFile, 'w');
fwrite(fid, char(nonEmptyLinesInAscii), '*char');
fclose(fid);

processStatus = [processStatus, newline, ...
                     datestr(now), ' >> ', '[END text outputing].'];


text:
function ySliderCallL()
%
%APA, 04/12/2010
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

global stateS
hs = stateS.legendSlider;
%hs = gcbo;
ud = get(hs,'userdata');
if isempty(ud)
    return;
end
ha = ud{1};
posYall = ud{2};
%posOldY = axis(ha);
posOldY = [get(ha,'xLim'), get(ha,'yLim')];
val = get(hs,'value');
% posYallO = posYall;
% posYallO(1) = posYall(1) -(posOldY(4)-posOldY(3));
valT = posYall(2) - val + posYall(1);
if valT<=posYall(1)
    posNewY = [posYall(1)-(posOldY(4)-posOldY(3)) posYall(1)];
elseif valT>=posYall(2)
    posNewY = [posYall(2)-(posOldY(4)-posOldY(3)) posYall(2)];
else
    posNewY = [posYall(2)-val posYall(2)-val+posYall(1)];
end

% posNewXY = posOldY;
% posNewXY(3:4) = posNewY;
% axis(ha,posNewXY)
%ha.YLim = posNewY;
set(ha,'YLim',posNewY);

% %%% Make lines/text outside axis not visible
% Iall = posYallO(1)+1:posYallO(2)-1;
% InotVisivle1 = find(Iall<posNewY(1));
% InotVisivle2 = find(Iall>posNewY(2));
% InotVisivle = [InotVisivle1 InotVisivle2];
% set(stateS.handle.legend.text(InotVisivle),'string','');
% set(stateS.handle.legend.lines(InotVisivle),'visible','off');

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
classdef MW8 < PROBLEM
% <multi/many> <real> <large/none> <constrained>
% Constrained benchmark MOP proposed by Ma and Wang

%------------------------------- Reference --------------------------------
% Z. Ma and Y. Wang, Evolutionary constrained multiobjective optimization:
% Test suite construction and performance comparisons. IEEE Transactions on
% Evolutionary Computation, 2019, 23(6): 972-986.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            if isempty(obj.M); obj.M = 3; end
            if isempty(obj.D); obj.D = 15; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values and constraint violations
        function Population = Evaluation(obj,varargin)
            X = varargin{1};
            X = max(min(X,repmat(obj.upper,size(X,1),1)),repmat(obj.lower,size(X,1),1));
            z = 1 - exp(-10*(X(:,obj.M:end) - ((obj.M:obj.D) - 1)/obj.D).^2);
            g = sum((1.5 + (0.1/obj.D)*z.^2 - 1.5*cos(2*pi*z)),2);
            PopObj = repmat(1+g,1,obj.M).*flip(cumprod([ones(size(X,1),1),cos(X(:,1:obj.M-1)*pi/2)],2),2).*[ones(size(X,1),1),sin(X(:,obj.M-1:-1:1)*pi/2)];
            l      = asin(PopObj(:,end)./sqrt(sum(PopObj.^2,2)));
            PopCon = sum(PopObj.^2,2) - (1.25 - 0.5*sin(6*l).^2).^2;
            Population  = SOLUTION(X,PopObj,PopCon,varargin{2:end});
            obj.FE      = obj.FE + length(Population);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
             R = UniformPoint(N,obj.M);
             R = R./repmat(sqrt(sum(R.^2,2)),1,obj.M);
             R(1-(1.25 - 0.5*sin(6*asin(R(:,end))).^2).^2>0,:) = [];
        end
        %% Generate the feasible region
        function R = GetPF(obj)
            if obj.M == 2
                [x,y] = meshgrid(linspace(0,1.2,400));
                z     = nan(size(x));
                fes   = x.^2 + y.^2 - (1.25-0.5*sin(6*asin(y./sqrt(x.^2+y.^2))).^2).^2 <= 0;
                z(fes & x.^2+y.^2>=1) = 0;
                R = {x,y,z};
            elseif obj.M == 3
                a = linspace(0,pi/2,40)';
                x = sin(a)*cos(a');
                y = sin(a)*sin(a');
                z = cos(a)*ones(size(a'));
                fes     = 1 - (1.25-0.5*sin(6*asin(z)).^2).^2 <= 0;
                z(~fes) = nan;
                R = {x,y,z};
            else
                R = [];
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function plotc(x,y,v,marker)
%FUNCTION PLOTC(X,Y,V,'MARKER') plots the values of v colour coded
% at the positions specified by x and y, and v (z-axis) in a 3-D axis
% system. A colourbar is added on the right side of the figure.
%
% The colorbar strectches from the minimum value of v to its
% maximum in 9 steps (10 values).
%
% The last argument is optional to define the marker being used. The
% default is a point. To use a different marker (such as circles, ...) send
% its symbol to the function (which must be enclosed in '; see example).
%
% The plot is actually a 3D plot but the orientation of the axis is set
% such that it appears to be a plane 2D plot. However, you can toggle
% between 2D and 3D view either by using the command 'view(3)' (for 3D
% view) or 'view(2)' (for 2D), or by interactively rotating the axis
% system.
%
% Example:
% Define three vectors
%    x=1:10;y=1:10;p=randn(10,1);
%    plotc(x,y,p)
%
%    x=randn(100,1);
%    y=2*x+randn(100,1);
%    p=randn(100,1);
%    plotc(x,y,p,'d')
%    view(3)
%
% Uli Theune, University of Alberta, 2004
% modified by Stephanie Contardo, British OCeanographic Data Centre, 2006
%

delete(gca)
if nargin <4
    marker='.';
end

map=colormap;
miv=min(v);
mav=max(v);
clrstep = (mav-miv)/size(map,1) ;
% Plot the points
hold on
for nc=1:size(map,1)
    iv = find(v>miv+(nc-1)*clrstep & v<=miv+nc*clrstep) ;
    plot3(x(iv),y(iv),v(iv),marker,'color',map(nc,:),'markerfacecolor',map(nc,:))
end
hold off

% Re-format the colorbar
h=colorbar;

set(h,'ylim',[1 length(map)]);
yal=linspace(1,length(map),10);
set(h,'ytick',yal);
% Create the yticklabels
ytl=linspace(miv,mav,10);
s=char(10,4);
for i=1:10
    if min(abs(ytl)) >= 0.001
        B=sprintf('%-4.3f',ytl(i));
    else
        B=sprintf('%-3.1E',ytl(i));
    end
    s(i,1:length(B))=B;
end
set(h,'yticklabel',s);
grid on
view(2)


text:
function [u,Du,eqn,info] = PoissonP3(node,elem,pde,bdFlag,option)
%% POISSONP2 Poisson equation: P3 quadratic element.
%
% u = PoissonP3(node,elem,pde,bdFlag,option) produces the cubic
%   finite element approximation of the Poisson equation
%
%       -div(d*grad(u))=f  in \Omega, with
%       Dirichlet boundary condition u=g_D on \Gamma_D,
%       Neumann boundary condition   d*grad(u)*n=g_N on \Gamma_N,
%       Robin boundary condition     g_R*u + d*grad(u)*n=g_N on \Gamma _R
%
% [u,Du,eqn,info] = PoissonP3(node,elem,pde,bdFlag,option)
%
% The usage is the same as Poisson. For quadratic elements, middle points
% of each edge are degree of freedom. See dofP2doc for detail.
%
% Example
%   femrateP2
%
% See also Poisson, Poisson3, Poisson3P2,PoissonP2
% Created by Jie Zhou
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

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

tic;
%% Construct Data Structure
[elem2dof,elem2edge,edge,bdDof]  = dofP3(elem);
N = size(node,1);  NT = size(elem,1); NE = size(edge,1);
Ndof = N + 2*NE + NT;

%% Compute geometric quantities and gradient of local basis
[Dlambda,area] = gradbasis(node,elem);

%% Assemble stiffness matrix
% Since Dphi_i*Dphi_j is four degree, so four order numerical quadrature rule is used here
if ~isfield(pde,'d'), pde.d = []; end
if ~isfield(option,'quadorder')
    option.quadorder = 5;   % default order 2(p-1)+1
end
[lambda, w] = quadpts(option.quadorder);
nQuad = size(lambda,1);
ii = zeros(55*NT,1); jj = zeros(55*NT,1); sA = zeros(55*NT,nQuad);
% generate sparse pattern
index = 0;
for i = 1:10
    for j = i:10
        ii(index+1:index+NT) = double(elem2dof(:,i));
        jj(index+1:index+NT) = double(elem2dof(:,j));
        index = index + NT;
    end
end
% compute non-zeros
for p = 1:nQuad
    % Dphi at quadrature points
    Dphip(:,:,1) = (27/2*lambda(p,1)*lambda(p,1)-9*lambda(p,1)+1).*Dlambda(:,:,1);
    Dphip(:,:,2) = (27/2*lambda(p,2)*lambda(p,2)-9*lambda(p,2)+1).*Dlambda(:,:,2);
    Dphip(:,:,3) = (27/2*lambda(p,3)*lambda(p,3)-9*lambda(p,3)+1).*Dlambda(:,:,3);
    Dphip(:,:,4) = 9/2*((3*lambda(p,2)*lambda(p,2)-lambda(p,2)).*Dlambda(:,:,3)+...
                   lambda(p,3)*(6*lambda(p,2)-1).*Dlambda(:,:,2));
    Dphip(:,:,5) = 9/2*((3*lambda(p,3)*lambda(p,3)-lambda(p,3)).*Dlambda(:,:,2)+...
                   lambda(p,2)*(6*lambda(p,3)-1).*Dlambda(:,:,3));
    Dphip(:,:,6) = 9/2*((3*lambda(p,3)*lambda(p,3)-lambda(p,3)).*Dlambda(:,:,1)+...
                   lambda(p,1)*(6*lambda(p,3)-1).*Dlambda(:,:,3));
    Dphip(:,:,7) = 9/2*((3*lambda(p,1)*lambda(p,1)-lambda(p,1)).*Dlambda(:,:,3)+...
                   lambda(p,3)*(6*lambda(p,1)-1).*Dlambda(:,:,1));
    Dphip(:,:,8) = 9/2*((3*lambda(p,1)*lambda(p,1)-lambda(p,1)).*Dlambda(:,:,2)+...
                   lambda(p,2)*(6*lambda(p,1)-1).*Dlambda(:,:,1));
    Dphip(:,:,9)  = 9/2*((3*lambda(p,2)*lambda(p,2)-lambda(p,2)).*Dlambda(:,:,1)+...
                   lambda(p,1)*(6*lambda(p,2)-1).*Dlambda(:,:,2));
    Dphip(:,:,10) = 27*(lambda(p,1)*lambda(p,2)*Dlambda(:,:,3)+lambda(p,1)*lambda(p,3)*Dlambda(:,:,2)+...
                   lambda(p,3)*lambda(p,2)*Dlambda(:,:,1));
    index = 0;
    for i = 1:10
        for j = i:10
            Aij = 0;
            if isempty(pde.d) || isnumeric(pde.d)
                Aij = Aij + w(p)*dot(Dphip(:,:,i),Dphip(:,:,j),2);
            else
                pxy = lambda(p,1)*node(elem(:,1),:) ...
                    + lambda(p,2)*node(elem(:,2),:) ...
                    + lambda(p,3)*node(elem(:,3),:);
                Aij = Aij + w(p)*dot(Dphip(:,:,i),Dphip(:,:,j),2).*pde.d(pxy);
            end
             if ~isempty(pde.d) && isnumeric(pde.d) % d is piecewise constant
                 Aij = pde.d.*Aij;
             end
            Aij = Aij.*area;
            sA(index+1:index+NT,p) = Aij;
            index = index + NT;
        end
    end
end
sA = sum(sA,2);
diagIdx = (ii == jj);   upperIdx = ~diagIdx;
A = sparse(ii(diagIdx),jj(diagIdx),sA(diagIdx),Ndof,Ndof);
AU = sparse(ii(upperIdx),jj(upperIdx),sA(upperIdx),Ndof,Ndof);
A = A + AU + AU';
clear Aij ii jj sA

%% Assemble right hand side by high order quadrature rule
b = zeros(Ndof,1);
if ~isfield(option,'fquadorder')
    option.fquadorder = 6;   % default order
end
if ~isfield(pde,'f') || (isreal(pde.f) && (pde.f==0))
    pde.f = [];
end
if ~isempty(pde.f)
    % quadrature points in the barycentric coordinate
    [lambda,w] = quadpts(option.fquadorder);
    nQuad = size(lambda,1);
    phi(:,1) = 0.5*(3*lambda(:,1)-1).*(3*lambda(:,1)-2).*lambda(:,1);
    phi(:,2) = 0.5*(3*lambda(:,2)-1).*(3*lambda(:,2)-2).*lambda(:,2);
    phi(:,3) = 0.5*(3*lambda(:,3)-1).*(3*lambda(:,3)-2).*lambda(:,3);
    phi(:,4) = 9/2*lambda(:,3).*lambda(:,2).*(3*lambda(:,2)-1);
    phi(:,5) = 9/2*lambda(:,3).*lambda(:,2).*(3*lambda(:,3)-1);
    phi(:,6) = 9/2*lambda(:,1).*lambda(:,3).*(3*lambda(:,3)-1);
    phi(:,7) = 9/2*lambda(:,1).*lambda(:,3).*(3*lambda(:,1)-1);
    phi(:,8) = 9/2*lambda(:,1).*lambda(:,2).*(3*lambda(:,1)-1);
    phi(:,9) = 9/2*lambda(:,1).*lambda(:,2).*(3*lambda(:,2)-1);
    phi(:,10) = 27*lambda(:,1).*lambda(:,2).*lambda(:,3);

    bt = zeros(NT,10);
    for p = 1:nQuad
        % quadrature points in the x-y coordinate
        pxy = lambda(p,1)*node(elem(:,1),:) ...
            + lambda(p,2)*node(elem(:,2),:) ...
            + lambda(p,3)*node(elem(:,3),:);
        if isfield(pde,'f') && isnumeric(pde.f)
            fp = pde.f;        % piecewise constant
        else
            fp = pde.f(pxy);   % function handle
        end
        for j = 1:10
            bt(:,j) = bt(:,j) + w(p)*phi(p,j)*fp;
        end
    end
    bt = bt.*repmat(area,1,10);
    b = accumarray(elem2dof(:),bt(:),[Ndof 1]);
end

%% Boundary Conditions
if nargin<=3, bdFlag = []; end
[AD,b,u,freeDof,isPureNeumann] = getbdP3(b);

%% Record assembeling time
assembleTime = toc;
if ~isfield(option,'printlevel'), option.printlevel = 1; end
if option.printlevel >= 2
    fprintf('Time to assemble matrix equation %4.2g s\n',assembleTime);
end

%% Solve the system of linear equations
if isempty(freeDof), return; end
% Set up solver
if isempty(option) || ~isfield(option,'solver')    % no option.solver
    if Ndof <= 2e3  % Direct solver for small size systems
        option.solver = 'direct';
    else            % Multigrid-type  solver for large size systems
        option.solver = 'mg';
    end
end

solver = option.solver;
% solve
switch solver
    case 'direct'
        tic;
        u(freeDof) = AD(freeDof,freeDof)\b(freeDof);
        residual = norm(b - AD*u);
        info = struct('solverTime',toc,'itStep',0,'err',residual,'flag',2,'stopErr',residual);
    case 'none'
        info = struct('solverTime',[],'itStep',0,'err',[],'flag',3,'stopErr',[]);
    case 'mg'
        option.x0 = u;
        option.solver = 'CG';
        option.tol = Ndof^(-2);
        [u,info] = mg(AD,b,elem,option,edge);
    case 'amg'
        option.solver = 'CG';
        option.tol = Ndof^(-2);
        [u(freeDof),info] = amg(AD(freeDof,freeDof),b(freeDof),option);
end
clear phi lambda
% post-process for pure Neumann problem
if isPureNeumann
    intguh = double(sparse(NT,1));
    [lambda,w] = quadpts(3);  %This is P3 element.
    nQuad = size(lambda,1);
    phi(:,1) = 0.5*(3*lambda(:,1)-1).*(3*lambda(:,1)-2).*lambda(:,1);
    phi(:,2) = 0.5*(3*lambda(:,2)-1).*(3*lambda(:,2)-2).*lambda(:,2);
    phi(:,3) = 0.5*(3*lambda(:,3)-1).*(3*lambda(:,3)-2).*lambda(:,3);
    phi(:,4) = 9/2*lambda(:,3).*lambda(:,2).*(3*lambda(:,2)-1);
    phi(:,5) = 9/2*lambda(:,3).*lambda(:,2).*(3*lambda(:,3)-1);
    phi(:,6) = 9/2*lambda(:,1).*lambda(:,3).*(3*lambda(:,3)-1);
    phi(:,7) = 9/2*lambda(:,1).*lambda(:,3).*(3*lambda(:,1)-1);
    phi(:,8) = 9/2*lambda(:,1).*lambda(:,2).*(3*lambda(:,1)-1);
    phi(:,9) = 9/2*lambda(:,1).*lambda(:,2).*(3*lambda(:,2)-1);
    phi(:,10) = 27*lambda(:,1).*lambda(:,2).*lambda(:,3);
     for     p = 1:nQuad
     intguh(:) = intguh(:) + w(p)*(phi(p,1)*u(elem2dof(:,1))+phi(p,2)*u(elem2dof(:,2))+phi(p,3)*u(elem2dof(:,3))...
                           + phi(p,4)*u(elem2dof(:,4))+phi(p,5)*u(elem2dof(:,5))+phi(p,6)*u(elem2dof(:,6))...
                           + phi(p,7)*u(elem2dof(:,7))+phi(p,8)*u(elem2dof(:,8))+phi(p,9)*u(elem2dof(:,9))...
                           + phi(p,10)*u(elem2dof(:,10)));
     end                 %compute the intgrable of uh.
     intguh = intguh.*area;
         uc = sum(intguh)/sum(area);
          u = u - uc;    % normalization for pure Neumann problem
end

% Output information
eqn = struct('A',AD,'b',b,'edge',edge,'freeDof',freeDof);
info.assembleTime = assembleTime;

%% Compute Du
Du = [];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunctions getbdP3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function [AD,b,u,freeDof,isPureNeumann] = getbdP3(b)
    %% Boundary conditions for Poisson equation: P3 quadratic FEM.
    %
    % The set up of boundary condition consists of two parts:
    %
    % 1) Modify the matrix for Dirichlet boundary nodes, which are not degree
    % of freedom. Values at these nodes are evaluatation of pde.g_D. The
    % original stiffness matrix A is turn into the matrix AD by enforcing
    % AD(fixedDof,fixedDof)=I, AD(fixedDof,freeDof)=0, AD(freeDof,fixedDof)=0.
    %
    % 2) Modify the right hand side b. The Neumann boundary integral is added
    % to b. For Dirichlet boundary ndoes, b(fixedDof) is the evaluation of
    % pde.g_D.
    %
    % Special attentation should be given for the pure Neumann boundary
    % condition. To enforce the compatible condition, the vector b should have
    % mean value zero. To avoid a singular matrix, the 1st node is chosen as
    % fixedDof.
    %
    % The order of assigning Neumann and Dirichlet boundary condition is
    % important to get the right setting at the intersection nodes of Dirichlet
    % and Neumann boundary edges.

    u = zeros(Ndof,1);

    %% Set up boundary and basic parameter
    if ~isfield(pde,'g_D'), pde.g_D = []; end
    if ~isfield(pde,'g_N'), pde.g_N = []; end
    if ~isfield(pde,'g_R'), pde.g_R = []; end

    %% Part 1: Modify the matrix for Dirichlet and Robin condition
    % Robin boundary condition
    Robin = [];
    idxR = (bdFlag(:) == 3);      %index of Robin edges in bdFlag
    if any(idxR)
        isRobin = false(NE,1);
        isRobin(elem2edge(idxR)) = true;
        Robin = edge(isRobin,:);  % Robin edges
    end
    if ~isempty(Robin) && ~isempty(pde.g_R) && ~(isnumeric(pde.g_R) && (pde.g_R == 0))
        if ~isfield(option,'gRquadorder')
            option.gRquadorder = 8;   % we should use six order  rule.
        end
        [lambdagR,weightgR] = quadpts1(option.gRquadorder);
        nQuadgR = size(lambdagR,1);
        % cubic bases (1--3--4--2)
        bdphi = zeros(nQuadgR,4);
        bdphi(:,1) = 0.5*(3*lambdagR(:,1)-1).*(3*lambdagR(:,1)-2).*lambdagR(:,1);
        bdphi(:,2) = 0.5*(3*lambdagR(:,2)-1).*(3*lambdagR(:,2)-2).*lambdagR(:,2);
        bdphi(:,3) = 9/2*lambdagR(:,1).*lambdagR(:,2).*(3*lambdagR(:,1)-1);
        bdphi(:,4) = 9/2*lambdagR(:,1).*lambdagR(:,2).*(3*lambdagR(:,2)-1);
        % length of edge
        el = sqrt(sum((node(Robin(:,1),:) - node(Robin(:,2),:)).^2,2));
        NR = size(Robin,1);
        ss = zeros(NR,4,4);
        for pp = 1:nQuadgR
            ppxy = lambdagR(pp,1)*node(Robin(:,1),:) ...
                 + lambdagR(pp,2)*node(Robin(:,2),:);
            gRp = pde.g_R(ppxy);
            for iR = 1:4
                for jR = iR:4   % only compute half of the off-diagonal part
                    ss(:,iR,jR) = ss(:,iR,jR) + ...
                    weightgR(pp)*gRp*bdphi(pp,iR).*bdphi(pp,jR);
                end
            end
        end
        ss(:) = ss(:).*repmat(el,16,1);
        Robin(:,3) = 2*find(isRobin)+N-1;  % the third one maps to corresponding dof
        Robin(:,4) = 2*find(isRobin)+N;    % the fourth one maps to corresponding dof
        index = 0;
        for iR = 1:4
            for jR = 1:4
                iiR(index+1:index+NR) = double(Robin(:,iR));
                jjR(index+1:index+NR) = double(Robin(:,jR));
                if jR>=iR
                    ssR(index+1:index+NR) = ss(:,iR,jR);
                else
                    ssR(index+1:index+NR) = ss(:,jR,iR);
                end
                index = index + NR;
            end
        end
        A = A + sparse(iiR,jjR,ssR,Ndof,Ndof);
    end

    % Find Dirichlet boundary dof: fixedDof
    fixedDof = []; freeDof = [];
    isFixedDof = false(Ndof,1);
    if ~isempty(bdFlag)
        isDirichlet(elem2edge(bdFlag(:)==1)) = true;
        isFixedDof(edge(isDirichlet,:)) = true;
        isFixedDof(N + 2*find(isDirichlet')-1) = true;
        isFixedDof(N + 2*find(isDirichlet')) = true;
        fixedDof = find(isFixedDof);
        freeDof = find(~isFixedDof);
    end
    if isempty(bdFlag) && ~isempty(pde.g_D) && isempty(pde.g_N) && isempty(pde.g_R)
        fixedDof = bdDof;
        isFixedDof(fixedDof) = true;
        freeDof = find(~isFixedDof);
    end
    isPureNeumann = false;
    if isempty(fixedDof) && isempty(Robin)  % pure Neumann boundary condition
        % pde.g_N could be empty which is homogenous Neumann boundary condition
        isPureNeumann = true;
        fixedDof = 1;
        freeDof = 2:Ndof;    % eliminate the kernel by enforcing u(1) = 0;
    end
    % Modify the matrix
    % Build Dirichlet boundary condition into the matrix AD by enforcing
    % |AD(fixedDof,fixedDof)=I, AD(fixedDof,freeDof)=0,
    % AD(freeDof,fixedDof)=0|.
    if ~isempty(fixedDof)
        bdidx = zeros(Ndof,1);
        bdidx(fixedDof) = 1;
        Tbd = sparse(1:Ndof,1:Ndof,bdidx,Ndof,Ndof);
        T = sparse(1:Ndof,1:Ndof,1-bdidx,Ndof,Ndof);
        AD = T*A*T + Tbd;
    else
        AD = A;
    end

    %% Part 2: Find boundary edges and modify the load b
    % Find boundary edges: Neumann and Robin
    Neumann = [];
    if ~isempty(bdFlag)
        idxN = (bdFlag(:) == 2);      % all Neumann edges in bdFlag
        Neumannidx = elem2edge(idxN | idxR); % index of Neumann and Robin edges
        % since boundary integral is also needed for Robin edges
        Neumann   = edge(Neumannidx,:);
    end
    if isempty(bdFlag) && (~isempty(pde.g_N) || ~isempty(pde.g_R))
        % no bdFlag, only pde.g_N or pde.g_R is given in the input
        Neumannidx = find(bdDof>N);
        Neumann = edge(Neumannidx,:);
    end

    % Neumann boundary condition
    if ~isempty(pde.g_N) && ~isempty(Neumann) && ~(isnumeric(pde.g_N) && (pde.g_N == 0))
        if ~isfield(option,'gNquadorder')
            option.gNquadorder = 6;
        end
        [lambdagN,weightgN] = quadpts1(option.gNquadorder);
        nQuadgN = size(lambdagN,1);
        % quadratic bases (1---3---4--2)
        bdphi = zeros(nQuadgN,4);
        bdphi(:,1) = 0.5*(3*lambdagN(:,1)-1).*(3*lambdagN(:,1)-2).*lambdagN(:,1);
        bdphi(:,2) = 0.5*(3*lambdagN(:,2)-1).*(3*lambdagN(:,2)-2).*lambdagN(:,2);
        bdphi(:,3) = 9/2*lambdagN(:,1).*lambdagN(:,2).*(3*lambdagN(:,1)-1);
        bdphi(:,4) = 9/2*lambdagN(:,1).*lambdagN(:,2).*(3*lambdagN(:,2)-1);
        % length of edge
        el = sqrt(sum((node(Neumann(:,1),:) - node(Neumann(:,2),:)).^2,2));
        ge = zeros(size(Neumann,1),4);
        for pp = 1:nQuadgN
            ppxy = lambdagN(pp,1)*node(Neumann(:,1),:) ...
                 + lambdagN(pp,2)*node(Neumann(:,2),:);
            gNp = pde.g_N(ppxy);
            ge(:,1) = ge(:,1) + weightgN(pp)*gNp*bdphi(pp,1);
            ge(:,2) = ge(:,2) + weightgN(pp)*gNp*bdphi(pp,2);
            ge(:,3) = ge(:,3) + weightgN(pp)*gNp*bdphi(pp,3);
            ge(:,4) = ge(:,4) + weightgN(pp)*gNp*bdphi(pp,4);
        end
        ge = ge.*repmat(el,1,4);
        b(1:N) = b(1:N) + accumarray(Neumann(:), [ge(:,1); ge(:,2)],[N,1]);
        b(N+2*Neumannidx-1) = b(N+2*Neumannidx-1) + ge(:,3);
        b(N+2*Neumannidx)   = b(N+2*Neumannidx) + ge(:,4);
    end

    % Dirichlet boundary conditions
    if ~isPureNeumann && ~isempty(fixedDof) && ...
       ~isempty(pde.g_D) && ~(isnumeric(pde.g_D) && (pde.g_D == 0))
        % interpolation
        idx = (fixedDof > N);  % index of edge nodes
        u(fixedDof(~idx)) = pde.g_D(node(fixedDof(~idx),:)); % bd value at vertex dofs
        % for P3,  we should divide the points of edge into two parts.
        bdEdgeIdx = fixedDof(idx) - N;
        %  First parts, the points  * is in  1---*------2
        bdEdgeMid = node(edge(isDirichlet,1),:)+(node(edge(isDirichlet,2),:) ...
                  - node(edge(isDirichlet,1),:))/3;
        u(N + bdEdgeIdx(1:2:end)) = pde.g_D(bdEdgeMid);
      %  Second parts, the points * is in  1------*---2
        bdEdgeMid = node(edge(isDirichlet,1),:)+2*(node(edge(isDirichlet,2),:)...
                  - node(edge(isDirichlet,1),:))/3;
        u(N + bdEdgeIdx(2:2:end)) = pde.g_D(bdEdgeMid);
        % modify the right hand side
        b = b - A*u;
    end
    if ~isPureNeumann % non-empty Dirichlet boundary condition
        b(fixedDof) = u(fixedDof);
    end

    % Pure Neumann boundary condition
    if isPureNeumann
        b = b - mean(b);   % compatilbe condition: sum(b) = 0
        b(1) = 0;
    end
    end % end of getbdP3
end % end of function PoissonP3

repo: Aero-matlab


text:
% 样条插值举例
clear; clc

X = [27.7, 28, 29, 30];
Y = [4.1, 4.3, 4.1, 3.0];
df0 = 3.0; dfn = -4.0;  % 第一类边界条件

n = length(X) - 1;
H = diff(X); % 每个插值小区间的长度
mu = H(1:n-1) ./ (H(1:n-1) + H(2:n));
lambda = 1 - mu;
% 计算二阶差商
Y1 = diff(Y) ./ diff(X);  % 一阶差商
Y2 = diff(Y1) ./ ( X(3:end) - X(1:end-2) );  % 二阶差商
% 计算右端项
d = 6*Y2;
d0 = 6/H(1) * (Y1(1) - df0);
dn = 6/H(end) * (dfn - Y1(end));
% 给出系数矩阵
A = 2*eye(n+1) + diag([mu(:);1],-1) +diag([1;lambda(:)],1);
b = [d0; d(:); dn];
M = A\b;

% 计算 s_k(x) 的系数, 两种形式
p = zeros(n,4); p0 = zeros(n,4);
for k = 1 : n
    p(k,1) = (M(k+1) - M(k))/(6*H(k));
    p(k,2) = M(k)/2;
    p(k,3) = (Y(k+1) - Y(k))/ H(k) - H(k)/6 * (2*M(k) + M(k+1));
    p(k,4) = Y(k);
    p0(k,1) = M(k)/(6*H(k));
    p0(k,2) = M(k+1)/(6*H(k));
    p0(k,3) = (Y(k) - M(k)*H(k)*H(k)/6) / H(k);
    p0(k,4) = (Y(k+1) - M(k+1)*H(k)*H(k)/6) / H(k);
end

% 调用 Matlab 的三次样条插值函数
pp = spline(X, [df0;Y(:);dfn]);

%输出结果
fprintf('按形式一输出: \n'); disp(p0);
fprintf('按形式二输出: \n'); disp(p);
fprintf('spline的结果(形式二): \n'); disp(pp.coefs);

repo: NumericalAnalysis


text:
function W = gsp_vec2adj(A,x,type)
%GSP_VEC2ADJ Create the matrix W with sparsity pattern A and entries x
%
%   Usage: W = gsp_vec2adj(A,x)
%              gsp_vec2adj(A,x,type)
%
%   Input parameters:
%       A       : Matrix of sparsity pattern (e.g. binary adjacency matrix)
%       x       : Vector of entries
%       type    : Type of matrix (string) (default 'sym' if A is symmetric)
%                    *  sym: W will be symmetric (x must be of size nnz(A)/2)
%                    * asym: W is asymmetric (x must be of size nnz(A))
%   Output parameters
%       W       : Weighted matrix
%
%   Create the matrix W with sparsity pattern A and entries x%
%

% Author: Francesco Grassi
% Date   : July 2016

if nargin<3
    if issymmetric(full(double(A)))
        type='sym';
    else
        type='asym';
    end
end

if ~nnz(A)==length(x) && ~nnz(A)/2==length(x)
    error('x must have the same size of nnz(A) or nnz(A)/2 if A is symmetric')
end

[N1,N2] = size(A);

switch type
    case 'sym'
        [i,j] = find(triu(A));
        W = sparse(i,j,x,N1,N2)+sparse(i,j,x,N1,N2)';
    case 'asym'
        [i,j] = find(A);
        W = sparse(i,j,x,N1,N2);
    otherwise
        error('Unknow type');
end

repo: gspbox


text:
%% inequality restrictions
% use the parameter restrictions block
% a >= b
% cos(tan(b))<=log(c)

%% equality restrictions
% One can use the parameter restrictions block as the interface
% to equality restrictions. Those restrictions can be nonlinear
% a = 3*b+4*c
% a = tan(cos(log(sin(b+c))))
% in the example above, if a is not estimated while b and c can be
%
% An alternative and more general solution is to use a steady state file or
% the steady_state_model block.

%% inequality restrictions as equality restrictions
% suppose we have the following inequality restriction
% a>=b
% this implies that there is a c>=0 such that a=b+c. One can exploit this
% insight to transform inequality restrictions into equality restrictions.
% We would need to define b and c as estimated while a will not appear in
% the parameterization block. The implied equality restriction will then be
% written either in the parameter_restrictions block (when this is formally
% implemented) or in a steady state file

%% fmincon and estimation with inequality restrictions
% lately, my experience with fmincon is that it does not always handle
% inequality restrictions well. Estimation may break down completely while
% fmincon tries to compute the derivatives of the inequality restrictions
% at a particular point. Unfortunately the function doing that step is
% p-coded and it is difficult to implement a workaround within fmincon
% itself (The Mathworks should be able to help on that...)
% in the models we estimate, this is sometimes the case for transition
% probabilities. Fmincon is known to work well when problems are smooth and
% differentiable. So in order to avoid the breakdown of fmincon, one has to
% be smart and RISE is:
% 1- use priors to smooth the shape of the objective function
% 2- use block_wise optimization
% 3- use a different optimization algorithm like bee_gate and its brothers,
% which I still need to push into the current version of RISE
%{
myblocks={
    {'beta_trans','kappa','rho_d','rho_s','rho_o','Lambda','Gamma',...
    'sig_d(vol,1)','sig_d(vol,2)','sig_s(vol,1)','sig_s(vol,2)',...
    'sig_r(vol,1)','sig_r(vol,2)','sig_o(oil,1)','sig_o(oil,2)','phi_pi',...
    'phi_y','rho_r'}
    {'vol_tp_1_2','vol_tp_2_1','oil_tp_1_2','oil_tp_2_1'}
};
m = estimate(m, 'optimizer', 'fmincon',...
    'estim_blocks',myblocks);
%}

author: jmaih


text:
function trueFalse = is2DData( Data )
%
% Purpose : Tests if the input argument represents valid 2D homogeneous
% coordinates.
%
% Uses (syntax) :
%   is2DData( Data )
%
% Input Parameters :
%   Data := the variable to be tested (should be 3xn, n greater than 0)
%
% Return Parameters :
%   trueFalse := 0 or 1 (false or true)
%
% Description and algorithms:
%   Tests the size of the input argument
%
% References :
%
% Cite this as :
%
% Author : Matthew Harker
% Date : July 13, 2005
% Version : 1.0
%--------------------------------------------------------------------------
% (c) 2005, O'Leary, Harker, University of Leoben, Leoben, Austria
% email: automation@unileoben.ac.at, url: automation.unileoben.ac.at
%--------------------------------------------------------------------------
% History:
%   Date:           Comment:
%   July 13, 2005   Original Version 1.0
%--------------------------------------------------------------------------
%
[m,n] = size( Data ) ;
%
if (m == 3) & (n > 0)
    %
    trueFalse = 1 ;
    %
else
    %
    trueFalse = 0 ;
    %
end

author: tobycollins

repo: IPPE

save_path: github-repos/MATLAB/tobycollins-IPPE


%==============================================================================
% This code is part of the Matlab-based toolbox
% LagLDDDM - A Lagrangian Gauss--Newton--Krylov Solver for Mass- and
%                        Intensity-Preserving Diffeomorphic Image Registration
%
% For details and license info see
% - https://github.com/C4IR/FAIR.m/tree/master/add-ons/LagLDDMM
%
% function [T,dT] = getPICMatrixAnalyticIntegral(omega,m,hp,xp,eps,varargin)
%
% Building push-forward matrix for a Particle-In-Cell method using anlytic
% integration. See also description in Sec. 3 of the paper:
%
%  ##2
%
%
% Input:
%  omega    - representation of computational domain
%  mc       - discretization size of sampling grid of data
%  mf       - discretization size of high-resolution image
%  hp       - cell size on particle grid used for computing particle's mass
%             which is rho_i \approx prod(hp)* \rho(x_i)
%  xp       - particle positions, size(xp) = [dim*np,1];
%  eps      - width of particles
%  rho      - particles mass, size(rho)=[np,1]
%  varargin - optional additional input
%
% Output:
%  T        - push-forward matrix, i.e. rho(xp) = C*rho
%  dT       - derivative of (C(xp)*rho) with respect to xp.
%
% =========================================================================
function [T,dT] = getPICMatrixAnalyticIntegral(omega,mc,mf,xp,varargin)

if nargin==0, help(mfilename); runMinimalExample; return; end

hp = (omega(2:2:end)-omega(1:2:end))./mf;   % cell-size in particle mesh
epsP = hp;                                  % width of particles
doDerivative = (nargout==2);
for k=1:2:length(varargin) % overwrites defaults
    eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;

dT = [];
dim = numel(omega)/2;
xp  = reshape(xp,[],dim);
np  = size(xp,1);         % number of particles
n   = prod(mc);            % number of voxels in sampling grid

h   = (omega(2:2:end)-omega(1:2:end))./mc;
pwidth = (ceil(epsP./h));    % upper bound for support of basis functions

% map particle positions to the domain [0,m(1)] x [0,m(dim)]
for i=1:dim
    xp(:,i) = (xp(:,i)-omega(2*i-1))/h(i);
end
% omega = zeros(1,2*dim);omega(2:2:end) = mc;

% get cell index of particles center of mass
P = ceil(xp);
w = xp-(P-1);

switch dim
    case 1
        B = reshape(int1D(w,pwidth,epsP,h,hp),[],1);
        J = repmat(  (1:np)',[2*pwidth+1,1]);
        I = reshape(bsxfun(@plus,P(:,1),-pwidth:pwidth),[],1);

        valid = (1<=I) & (I<=mc(1));
        I = I(valid); J = J(valid); B = B(valid);

        T = sparse(I,J,B,n,np);
        if doDerivative
            bx = reshape(diff1D(w,pwidth,epsP,h,hp),[],1);
            ids = (1:numel(bx)); ids = ids(valid);
            dT = @(rho) sparse(I,J,rho(J).*bx(ids),n,np)/h;
        end

    case 2
        s2i  = @(i) i(:,1)+ mc(1)   *(i(:,2)-1);   % sub2ind for cell-centered grid

        B1 = int1D(w(:,1),pwidth(1),epsP(1),h(1),hp(1));
        B2 = int1D(w(:,2),pwidth(2),epsP(2),h(2),hp(2));
        if doDerivative
            b1 = diff1D(w(:,1),pwidth(1),epsP(1),h(1),hp(1));
            b2 = diff1D(w(:,2),pwidth(2),epsP(2),h(2),hp(2));
        end

        nVoxel = prod([size(B1,2),size(B2,2)]);
        J = repmat((1:np)',[nVoxel 1]);
        I = zeros([np*nVoxel,2]); B = zeros(np*nVoxel,1); pp = 1;
        bx = B; by = B;
        for px = -pwidth(1):pwidth(1)
            for py = -pwidth(2):pwidth(2)
                idx = (pp-1)*np+(1:np); pp = pp+1;
                I(idx,:) = [P(:,1)+px, P(:,2)+py];

                B(idx) =  B1(:,px+pwidth(1)+1).*B2(:,py+pwidth(2)+1);

                if doDerivative
                    bx(idx) = B2(:,py+pwidth(2)+1).*b1(:,px+pwidth(1)+1);
                    by(idx) = B1(:,px+pwidth(1)+1).*b2(:,py+pwidth(2)+1);
               end
            end
        end
        valid = (1<=I(:,1)) & (I(:,1)<=mc(1)) & (1<=I(:,2)) & (I(:,2)<=mc(2));
        I = s2i(I(valid,:));
        J = J(valid);
        B = B(valid);
        T = sparse(I,J,B,n,np);

        if doDerivative
            ids = (1:numel(bx))'; ids = ids(valid);
            dT = @(rho) [ sparse(I,J,rho(J).*bx(ids)/h(1),n,np), sparse(I,J,rho(J).*by(ids)/h(2),n,np)];
        end
    case 3
        % sub2ind for cell-centered grid
        s2i  = @(i) i(:,1)+ mc(1)   *(i(:,2)-1)+(mc(1)  )*(mc(2)  )*(i(:,3)-1);

        B1 = int1D(w(:,1),pwidth(1),epsP(1),h(1),hp(1));
        B2 = int1D(w(:,2),pwidth(2),epsP(2),h(2),hp(2));
        B3 = int1D(w(:,3),pwidth(3),epsP(3),h(3),hp(3));

        nVoxel = prod([size(B1,2),size(B2,2),size(B3,2)]);
        I = zeros([np*nVoxel,3]); J = repmat( (1:np)',[nVoxel,1]); B = zeros(np*nVoxel,1);
        pp = 1;
        if doDerivative
            b1 = diff1D(w(:,1),pwidth(1),epsP(1),h(1),hp(1));
            b2 = diff1D(w(:,2),pwidth(2),epsP(2),h(2),hp(2));
            b3 = diff1D(w(:,3),pwidth(3),epsP(3),h(3),hp(3));

            bx = B; by = B; bz = B;
        end
        for pz = -pwidth(3):pwidth(3)
            i3 = P(:,3)+pz;
            B3t = B3(:,pz+pwidth(3)+1);
            for py = -pwidth(2):pwidth(2)
                i2 = P(:,2)+py;
                B2t = B2(:,py+pwidth(2)+1);
                B23 = B3t.*B2t;
                for px = -pwidth(1):pwidth(1)
                    idx = (pp-1)*np+(1:np); pp = pp+1;
                    I(idx,1) = P(:,1)+px;
                    I(idx,2)=i2;
                    I(idx,3)=i3;

                    % remove cells that lie out of the domain
                    B1t = B1(:,px+pwidth(1)+1);

                    B(idx)   =  B1t.*B23;

                    if doDerivative
                       % compute derivatives of weights
                        bx(idx)     = B23.*b1(:,px+pwidth(1)+1);
                        by(idx)     = B1t.*B3t.*b2(:,py+pwidth(2)+1);
                        bz(idx)     = B1t.*B2t.*b3(:,pz+pwidth(3)+1);
                    end
                end
            end
        end
        valid =    (1<=I(:,1)) & (I(:,1)<=mc(1)) ...
            & (1<=I(:,2)) & (I(:,2)<=mc(2)) ...
            & (1<=I(:,3)) & (I(:,3)<=mc(3));

        I = s2i(I(valid,:));
        J = J(valid);
        B = B(valid);
        T = sparse(I,J,B,n,np);

        if doDerivative
            ids = (1:numel(bx))'; ids = ids(valid);
            D1 = @(rho) sparse(I,J,rho(J).*bx(ids)/h(1),n,np);
            D2 = @(rho) sparse(I,J,rho(J).*by(ids)/h(2),n,np);
            D3 = @(rho) sparse(I,J,rho(J).*bz(ids)/h(3),n,np);
            dT = @(rho) [D1(rho), D2(rho), D3(rho)];
        end

    otherwise
        error('dimension must be 1,2,3')
end

function Bij = int1D(w,pwidth,eps,h,hp)
Bij = zeros(numel(w),2*pwidth+1);
Bleft = B(-pwidth-w,eps,h);
for p = -pwidth:pwidth
    Bright = B(1+p-w,eps,h);
    Bij(:,p+pwidth+1)  = hp*(Bright - Bleft);
    Bleft = Bright;
end
function Bij = diff1D(w,pwidth,eps,h,hp)
Bij = zeros(numel(w),2*pwidth+1);
Bleft = b(-pwidth-w,eps,h);
for p = -pwidth:pwidth
    Bright = b(1+p-w,eps,h);
    Bij(:,p+pwidth+1)  = hp*(Bright - Bleft);
    Bleft = Bright;
end
Bij = - Bij;

function bij = b(x,eps,h)
bij = zeros(numel(x),1);

ind1 = (-eps/h<=x)&(x<=0);
ind2 = (0<x)&(x<=eps/h);

bij(ind1)  = 1 + h*x(ind1)./eps;
bij(ind2)  = 1 - h*x(ind2)./eps;
bij = bij /eps;

function Bij = B(x,eps,h)
Bij = zeros(numel(x),1);

ind1 = (-eps/h<=x)&(x<=0);
ind2 = (0<x)&(x<=eps/h);
ind3 = (eps/h<x);

Bij(ind1) = x(ind1) + 1./(2*eps/h).*x(ind1).^2+eps/(h*2);
Bij(ind2) = x(ind2) - 1./(2*eps/h).*x(ind2).^2+eps/(h*2);
Bij(ind3) = eps/h;
Bij = Bij/eps;

function [T,dT] = derivativeTestFctn(omega,m,hc,xp,rho)
[T,dT] = feval(mfilename,omega,m,hc,xp);
T = T*rho;
dT = dT(rho);

function runMinimalExample

%  ========== 1 D ==================
omega = [0 1]; mc = 64; mf = 64;
hc     = (omega(2:2:end)-omega(1:2:end))./mc;
hf     = (omega(2:2:end)-omega(1:2:end))./mf;
rho   = ones(mc,1); rho([1:8,end-8:end]) = 0;

xp = getCellCenteredGrid(omega,mc);
xp = xp +.02;

% transport rho
C = feval(mfilename,omega,mf,mc,xp);
rhonew = C*rho;

% visualize result
figure(1);clf;
subplot(1,2,1)
plot(getCellCenteredGrid(omega,mc),rho);
title(sprintf('rho, mass:%e',prod(hc)*sum(rho)));

subplot(1,2,2)
plot(getCellCenteredGrid(omega,mf),rhonew)
title(sprintf('rhonew, mass:%e',prod(hf)*sum(rhonew)));
% check derivative
fctn = @(xp) derivativeTestFctn(omega,mf,mc,xp,rho);
checkDerivative(fctn,xp(:));
%  ========== 2 D ==================
omega = [-2 2 -2 2]; mc = [32 32]; mf = mc*2;
hc     = (omega(2:2:end)-omega(1:2:end))./mc;
hf     = (omega(2:2:end)-omega(1:2:end))./mf;
f      = @(x) 0 + (sqrt(x(:,1).^2 + .4*x(:,2).^2)<.7);
rho    = f(reshape(getCellCenteredGrid(omega,mc),[],2));

% choose parameters
xp = getCellCenteredGrid(omega,mc);
xp = xp+.1*hc(1) ;

% transport rho
C = feval(mfilename,omega,mf,mc,xp);
rhonew = C*rho;

% visualize result
figure(2); clf;
subplot(1,2,1);
viewImage2Dsc(rho,omega,mc);
title(sprintf('rho, mass:%e',prod(hc)*sum(rho)));

subplot(1,2,2);
viewImage2Dsc(rhonew,omega,mf);
title(sprintf('rho, mass:%e',prod(hf)*sum(rhonew)));

% check derivative
fctn = @(xp) derivativeTestFctn(omega,mf,mc,xp,rho);
checkDerivative(fctn,xp(:));

%  ========== 3 D ==================
omega = [-2 2 -2 2 -3 3]; mc = [16 32 24]; mf = mc;
hc     = (omega(2:2:end)-omega(1:2:end))./mc;
hf     = (omega(2:2:end)-omega(1:2:end))./mf;
f      = @(x) 0 + (sqrt(x(:,1).^2 + .3*x(:,2).^2 + .2*x(:,3).^2)<.7);
rho    = f(reshape(getCellCenteredGrid(omega,mc),[],3));

% choose parameters
xp = getCellCenteredGrid(omega,mc);
xp = xp+0.1*hc(1) ;

% transport rho
C = feval(mfilename,omega,mf,mc,xp);
rhonew = C*rho;

% visualize result
figure(2); clf;
subplot(1,2,1);
imgmontage(rho,omega,mc);
title(sprintf('rho, mass:%e',prod(hc)*sum(rho)));

subplot(1,2,2);
imgmontage(rhonew,omega,mf);
title(sprintf('rho, mass:%e',prod(hf)*sum(rhonew)));

% check derivative
fctn = @(xp) derivativeTestFctn(omega,mf,mc,xp,rho);
checkDerivative(fctn,xp(:));


text:
%
% testing the gaver-stehfest module
function y=testgs(L)
% L is the number of coefficients
% (examples: L=8, 10, 12, 14, 16, so on..)
sum=0.0;
for l=1:20
    t(l)=l * 0.1;
    calcv(l)=gavsteh('fun1',t(l),L);
    exactv(l)=t(l);
    sum=sum+(1- exactv/calcv)^2;
end
result1=[exactv' calcv' calcv'-exactv']
relerr=sqrt(sum)

%another example
sum=0.0;
for l=1:20
    t(l)=l * 0.1;
    calcv(l)=gavsteh('fun2',t(l),L);
    exactv(l)=sin(t(l));
    sum=sum+(1- exactv/calcv)^2;
end
result2=[exactv' calcv' calcv'-exactv']
relerr=sqrt(sum)


text:
function [grp_phase, cep, ts] = modified_group_delay_feature(file_name, rho, gamma, num_coeff, frame_shift)

%input:
%     file_name: path for the waveform. The waveform should have a header
%     rho: a parameter to control the shape of modified group delay spectra
%     gamma: a parameter to control the shape of the modified group delay spectra
%     num_coeff: the desired feature dimension
%     [frame_shift]:
%
%output:
%     grp_phase: modifed gropu delay spectrogram
%     cep: modified group delay cepstral feature.
%     ts: time instants at the center of each analysis frame.
%
%Example:
%     [grp_phase, cep, ts] = modified_group_delay_feature('./100001.wav', 0.4, 0.9, 12);
% Please tune rho and gamma for better performance
%     See also: howtos/HOWTO_features.m
%
% by Zhizheng Wu (zhizheng.wu@ed.ac.uk)
% http://www.zhizheng.org
%
% The code has been used in the following three papers:
% Zhizheng Wu, Xiong Xiao, Eng Siong Chng, Haizhou Li, "Synthetic speech detection using temporal modulation feature", IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP) 2013.
% Zhizheng Wu, Tomi Kinnunen, Eng Siong Chng, Haizhou Li, Eliathamby Ambikairajah, "A study on spoofing attack in state-of-the-art speaker verification: the telephone speech case", Asia-Pacific Signal and Information Processing Association Annual Summit and Conference (APSIPA ASC) 2012.
% Zhizheng Wu, Eng Siong Chng, Haizhou Li, "Detecting Converted Speech and Natural Speech for anti-Spoofing Attack in Speaker Recognition", Interspeech 2012.
%
% feel free to modify the code and welcome to cite above papers :)

[speech,fs]  = wavread(file_name);

if nargin<2;
    rho = 0.4;
end
if nargin<3;
    gamma = 0.9;
end
if nargin<4;
    num_coeff = 12;
end
frame_length = 0.025; %msec
if nargin<5;
    frame_shift  = 0.010; %msec
end
NFFT         = 512;
pre_emph     = true;

%%% Pre-emphasis + framing
if (pre_emph)
    speech = filter([1 -0.97], 1, speech);
end;
frame_length = round((frame_length)*fs);
frame_shift = round((frame_shift)*fs);
[frames, ts] = enframe(speech, hamming(frame_length), frame_shift);
ts = (ts-1)/fs;

frame_num    = size(frames, 1);
frame_length = size(frames, 2);
delay_vector = [1:1:frame_length];
delay_matrix = repmat(delay_vector, frame_num, 1);

delay_frames = frames .* delay_matrix;

x_spec = fft(frames', NFFT);
y_spec = fft(delay_frames', NFFT);
x_spec = x_spec(1:NFFT/2+1, :);
y_spec = y_spec(1:NFFT/2+1, :);

temp_x_spec = abs(x_spec);

dct_spec = dct(medfilt1(log(temp_x_spec), 5));
smooth_spec = idct(dct_spec(1:30,:), NFFT/2+1);

grp_phase1 = (real(x_spec).*real(y_spec) + imag(y_spec) .* imag(x_spec)) ./(exp(smooth_spec).^ (2*rho));
grp_phase = (grp_phase1 ./ abs(grp_phase1)) .* (abs(grp_phase1).^ gamma);
grp_phase = grp_phase ./ (max(max(abs(grp_phase))));

grp_phase(isnan(grp_phase)) = 0.0;

cep = dct(grp_phase);
cep = cep(2:num_coeff+1, :)';

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep