text:
% Parallel HYperslab Projection along Affine SubSpace (PHYPASS) algorithm
%
% M. Takizawa and M. Yukawa, "An Efficient Data-Reusing Kernel Adaptive
% Filtering Algorithm Based on Parallel Hyperslab Projection Along Affine
% Subspace," in Proc. ICASSP, May. 2013, pp.3557-3561.
%
% Remark: version01, August 2013
% Contributor for this code: Masa-aki Takizawa
%
% This file is part of the Kernel Adaptive Filtering Toolbox for Matlab.
% https://github.com/steven2358/kafbox/

classdef phypass < kernel_adaptive_filter

    properties (GetAccess = 'public', SetAccess = 'private')
        mu = 0.5; % step size
        s = 1; % number of update coefficients (dictionary size)
        sigma = 0.95; % threshold for dictionary
        p = 8; % number of hyper slabs
        kerneltype = 'gauss'; % kernel type
        kernelpar = 1; % kernel parameter
        omega = 1/8; % weight
    end

    properties (GetAccess = 'public', SetAccess = 'private')
        mem = []; % memory
        dict = []; % dictionary
        alpha = []; % expansion coefficients
        d = []; % output window
    end

    methods
        function kaf = phypass(parameters) % constructor
            if (nargin > 0) % copy valid parameters
                for fn = fieldnames(parameters)'
                    if ismember(fn,fieldnames(kaf))
                        kaf.(fn{1}) = parameters.(fn{1});
                    end
                end
            end
        end

        function y_est = evaluate(kaf,x) % evaluate the algorithm
            if size(kaf.dict,1)>0
                k = kernel(kaf.dict,x,kaf.kerneltype,kaf.kernelpar);
                y_est = k'*kaf.alpha;
            else
                y_est = zeros(size(x,1),1);
            end
        end

        function train(kaf,x,y) % train the algorithm
            if (length(kaf.d) < kaf.p) % grow the memory
                kaf.mem = [kaf.mem; x];
                kaf.d = [kaf.d; y];
            else % sliding window
                kaf.mem = [kaf.mem(2:end,:); x];
                kaf.d = [kaf.d(2:end); y];
            end
            if size(kaf.dict,1) == 0 % initialize
                kaf.dict = x;
                kaf.alpha = 0;
            else
                k = kernel(x,kaf.dict,kaf.kerneltype,kaf.kernelpar);
                if max(k) < kaf.sigma % coherence criterion
                    kaf.dict = [kaf.dict ; x]; % add the current input x
                    kaf.alpha = [kaf.alpha ; 0];
                end
            end
            num = length(kaf.d); % length of memory

            num_d = min(size(kaf.dict,1),kaf.s); % number of dictionary elements to update

            dict_id = zeros(num_d,num); % memory of dictionary index

            Kmemdict = kernel(kaf.mem,kaf.dict,kaf.kerneltype,kaf.kernelpar);
            Kdict = kernel(kaf.dict,kaf.dict,kaf.kerneltype,kaf.kernelpar);

            for k=1:num
                d_check = Kmemdict(k,:); % kernel between k'th memory element and full dictionary
                [mm,ii] = sort(d_check,'descend'); %#ok<ASGLU>
                dict_id(:,k) = ii(1:num_d); % memory indexes
            end

            k_n = zeros(num_d,num_d,num);
            y_n = zeros(num_d,num);
            alpha_new = zeros(num_d,num);
            for k=1:num
                k_n(:,:,k) = Kdict(dict_id(:,k),dict_id(:,k));
                y_n(:,k) = Kmemdict(k,dict_id(:,k))';
                alpha_new(:,k) = k_n(:,:,k)\y_n(:,k); % compute projections onto dictionary subspaces
            end
            ln_nume = 0;
            numerator = zeros(num,1);
            denominator = zeros(num,1);
            beta = zeros(num,1);

            for k=1:num
                numerator(k) = Kmemdict(k,:)*kaf.alpha;
                denominator(k) = Kmemdict(k,dict_id(:,k))*alpha_new(:,k);
                beta(k) = (kaf.d(k) - numerator(k))/(denominator(k)); %  progress of the projection onto hyperplanes
                ln_nume = ln_nume + kaf.omega * beta(k)^2 *...
                    alpha_new(:,k)' * k_n(:,:,k) * alpha_new(:,k); % numerator of the extrapolation coefficcient
            end

            G = zeros(num,num);
            for k = 1:num
                Grow =  alpha_new(:,k)'*Kdict(dict_id(:,k),dict_id);
                sti = 1;
                ndi = num_d;
                for l = 1:num
                    G(k,l) = Grow(1,sti:ndi)*alpha_new(:,l);
                    sti = ndi + 1;
                    ndi = ndi + num_d;
                end
            end

            ln_deno = kaf.omega^2 * beta' * G * beta; % denominator of the extrapolation coefficient
            ln = ln_nume/(ln_deno); % extrapolation coefficient
            alpha_new2 = zeros(length(kaf.alpha),num);
            for k=1:num
                alpha_new2(dict_id(:,k),k) = alpha_new(:,k);
            end
            for k=1:num
                kaf.alpha = kaf.alpha + kaf.mu * ln *...
                    kaf.omega * beta(k) * alpha_new2(:,k); % update expansion coefficients
            end

        end
    end
end

repo: kafbox


text:
function pass = test_isequal( pref )

% Grab some preferences
if ( nargin == 0 )
    pref = chebfunpref();
end

% Test with function 1
f = ballfun(ones(21,20,22));
g = f+f-f;

pass(1) = (isequal(f,g) && isequal(g,f));
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% EX_ARTICLE_SECTION_513: short example with macro-element integration rule.
%
% Example to solve the problem
%
%    - div ( grad (u)) = (8-9*sqrt(x.^2+y.^2)).*sin(2*atan(y./x))./(x.^2+y.^2)  in Omega
%                    u = 0                                                      on Gamma
%
% with                Omega = (1 < x^2+y^2 < 4) & (x > 0) & (y > 0)
% and exact solution      u = (x.^2+y.^2-3*sqrt(x.^2+y.^2)+2).*sin(2.*atan(y./x))
%
% using a quadrature rule on macro-elements as proposed in
%
%
% [1] Hughes, Reali, Sangalli. Comput. Methods Appl. Mech. Engrg.
% Volume 199, Issues 5-8, 1 January 2010, Pages 301-313
%
% This is the example of Section 5.1.3 in the article
%
% C. De Falco, A. Reali, R. Vazquez
% GeoPDEs: a research tool for IsoGeometric Analysis of PDEs
%
% Copyright (C) 2009, 2010 Carlo de Falco
% Copyright (C) 2011, 2015 Rafael Vazquez
%
%    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/>.

nodes   = [5.168367524056075e-02, 2.149829914261059e-01, ...
           3.547033685486441e-01, 5.000000000000000e-01, ...
           6.452966314513557e-01, 7.850170085738940e-01, ...
           9.483163247594394e-01];
weights = [1.254676875668223e-01, 1.708286087294738e-01, ...
           1.218323586744639e-01, 1.637426900584793e-01, ...
           1.218323586744638e-01, 1.708286087294741e-01, ...
           1.254676875668223e-01];
rule    = [nodes; weights];

geometry = geo_load ({@ring_polar_map, @ring_polar_map_der});

[knots, breaks] = kntuniform ([10 10], [2 2], [1 1]);
breaks{1} ([2:3:end-2,3:3:end-1]) = [];
breaks{2} ([2:3:end-2,3:3:end-1]) = [];
[qn, qw]  = msh_set_quad_nodes (breaks, {rule, rule}, [0 1]);
msh = msh_cartesian (breaks, qn, qw, geometry);

space = sp_bspline (knots, [2, 2], msh);

mat = op_gradu_gradv_tp (space, space, msh);
rhs = op_f_v_tp (space, msh, @(x, y) (8-9*sqrt(x.^2+y.^2)).*sin(2*atan(y./x))./(x.^2+y.^2));

drchlt_dofs = [];
for iside = 1:4
  drchlt_dofs = union (drchlt_dofs, space.boundary(iside).dofs);
end
int_dofs = setdiff (1:space.ndof, drchlt_dofs);

u = zeros (space.ndof, 1);
u(int_dofs) = mat(int_dofs, int_dofs) \ rhs(int_dofs);

sp_to_vtk (u, space, geometry, [20 20], 'laplace_solution.vts', 'u')
err = sp_l2_error (space, msh, u, @(x,y)(x.^2+y.^2-3*sqrt(x.^2+y.^2)+2).*sin(2.*atan(y./x)))

%!demo
%! ex_article_section_513;

%!test
%! nodes   = [5.168367524056075e-02, 2.149829914261059e-01, ...
%!            3.547033685486441e-01, 5.000000000000000e-01, ...
%!            6.452966314513557e-01, 7.850170085738940e-01, ...
%!            9.483163247594394e-01];
%! weights = [1.254676875668223e-01, 1.708286087294738e-01, ...
%!            1.218323586744639e-01, 1.637426900584793e-01, ...
%!            1.218323586744638e-01, 1.708286087294741e-01, ...
%!            1.254676875668223e-01];
%! rule    = [nodes; weights];
%! geometry = geo_load ({@ring_polar_map, @ring_polar_map_der});
%! [knots, breaks] = kntuniform ([10 10], [2 2], [1 1]);
%! breaks{1} ([2:3:end-2,3:3:end-1]) = [];
%! breaks{2} ([2:3:end-2,3:3:end-1]) = [];
%! [qn, qw]  = msh_set_quad_nodes (breaks, {rule, rule}, [0 1]);
%! msh = msh_cartesian (breaks, qn, qw, geometry);
%! space = sp_bspline (knots, [2, 2], msh);
%! mat = op_gradu_gradv_tp (space, space, msh, @(x, y) ones (size (x)));
%! rhs = op_f_v_tp (space, msh, @(x, y) (8-9*sqrt(x.^2+y.^2)).*sin(2*atan(y./x))./(x.^2+y.^2));
%! drchlt_dofs = [];
%! for iside = 1:4
%! drchlt_dofs = union (drchlt_dofs, space.boundary(iside).dofs);
%! end
%! int_dofs = setdiff (1:space.ndof, drchlt_dofs);
%! u = zeros (space.ndof, 1);
%! u(int_dofs) = mat(int_dofs, int_dofs) \ rhs(int_dofs);
%! err = sp_l2_error (space, msh, u, @(x,y)(x.^2+y.^2-3*sqrt(x.^2+y.^2)+2).*sin(2.*atan(y./x)));
%! assert (msh.nel, 9)
%! assert (space.ndof, 121)
%! assert (err, 4.71881595641687e-05, 1e-16)

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function plotData(X, y)
%PLOTDATA Plots the data points X and y into a new figure
%   PLOTDATA(x,y) plots the data points with + for the positive examples
%   and o for the negative examples. X is assumed to be a Mx2 matrix.

% Create New Figure
figure; hold on;

% ====================== YOUR CODE HERE ======================
% Instructions: Plot the positive and negative examples on a
%               2D plot, using the option 'k+' for the positive
%               examples and 'ko' for the negative examples.
%

% Find Indices of Positive and Negative Examples
pos = find(y==1);
neg = find(y == 0);

plot(X(pos, 1), X(pos, 2), 'k+','LineWidth', 2, 'MarkerSize', 7);
plot(X(neg, 1), X(neg, 2), 'ko', 'MarkerFaceColor', 'y', 'MarkerSize', 7);

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

hold off;

end

author: Borye


text:
%%   FDI Wind Turbines
% This script containing parameters for Bench Mark Model By Peter Fogh
% Odgaard, kk-electronic a/s, Viby J Denmark -Date 7.11.2008-
% with little modification by Nassim Laouti Date 13.11.2010
% Latest modification Date 16.02.2012

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c3=10;c4=1000;P_r=4.8e6;
% Simple time and time
Ts=1/100;Time=Ts:Ts:4400;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Wind data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
D_v_wind=[Time', 9+4*sin(0.01*Time)'];

load winddata
D_v_wind=[windtime windspeed];

%% Wind Model
seed1 = 256;
seed2 = 894;
turbulence_seeds=[seed1 seed2];
R=57.5;
H=87;
k =4.7;
a=2.2;
alpha = 0.1;
D_rotor=2*R;
r = D_rotor/2;
m = 1+alpha*(alpha-1)*R^2/(8*H^2);
Length_scale = 600; %[m]
L=Length_scale;
Turbulence_intensity = 12;
turb_int=Turbulence_intensity;
T_sample=0.05;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Pitch and Blade Model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
omega_n=11.11;xi=0.6;rho=1.225;
R=57.5;r0 = 1.5;

% cq table
load AeroDynamics

% Fault models
% Fault 6
xi2=0.45;omega_n2=5.73;
% Fault 7
xi3=0.9;omega_n3=3.42;

%transfers to ss models
[Apb,Bpb,Cpb,Dpb]=tf2ss([omega_n^2],[1 2*xi*omega_n omega_n^2]);
[Apb1,Bpb1,Cpb1,Dpb1]=tf2ss([omega_n2^2],[1 2*xi2*omega_n2 omega_n2^2]);
[Apb2,Bpb2,Cpb2,Dpb2]=tf2ss([omega_n3^2],[1 2*xi3*omega_n3 omega_n3^2]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Drive Train Model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B_dt=775.49;
B_r=7.11;
B_g=45.6;
N_g=95;
K_dt=2.7e9
eta_dt=0.97;
J_g=390
J_r=55e6

Addt=[-(B_dt+B_r)/J_r B_dt/N_g/J_r -K_dt/J_r; eta_dt*B_dt/N_g/J_g -(eta_dt*B_dt/N_g^2+B_g)/J_g eta_dt*K_dt/N_g/J_g; 1 -1/N_g 0];
Bddt=[1/J_r 0; 0 -1/J_g;0 0];
Cddt=[1 0 0;0 1 0];
Dddt=[0 0;0 0];

% Fault models
% Fault 9
eta_dt2=.91

Addt2=[-(B_dt+B_r)/J_r B_dt/N_g/J_r -K_dt/J_r; eta_dt2*B_dt/N_g/J_g -(eta_dt2*B_dt/N_g^2+B_g)/J_g eta_dt2*K_dt/N_g/J_g; 1 -1/N_g 0];
Bddt2=[1/J_r 0; 0 -1/J_g;0 0];
Cddt2=[1 0 0;0 1 0];
Dddt2=[0 0;0 0];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Generator & Converter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
alpha_gc=1/20e-3;
eta_gc=0.98;

%Fault models
%fault 8
Constant_tau_gc=100;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Controller
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
K_opt=.5* rho* pi*R^2*eta_dt*R^3*0.4554/(N_g*8)^3

K_i=1;
K_p=4;
Omega_nom=162;
Omega_delta=5;
P_r=4.8e6;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Sensors
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m_vw=1.5;
sigma_vm=0.5;
m_omega_r=0;
sigma_omega_r=0.004*2*pi;
m_omega_g=0;
sigma_omega_g=0.05;
m_tau_g=0;
sigma_tau_g=90;
m_P_g=0;
sigma_P_g=1e3;
m_Beta=0;
sigma_Beta=0.2;

%% Fault models (choose the scenario of fault parameters and magnitude)

% %1st fault scenario
% Constant_Beta_1_m1=5;Gain_Beta_2_m2=1.2;Constant_Beta_3_m1=10;Constant_Omega_r_m1=1.4;Gain_Omega_r_m2=1.1;Gain_Omega_g_m2=0.9;
% xi2=0.45;omega_n2=5.73;xi3=0.9;omega_n3=3.42;Constant_tau_gc=100;eta_dt2=.92;Constant_Omega_g_m1=140;

% %2nd fault scenario
% Constant_Beta_1_m1=6;Gain_Beta_2_m2=1.5;Constant_Beta_3_m1=8;Constant_Omega_r_m1=0.2;Gain_Omega_r_m2=1.2;Gain_Omega_g_m2=0.8;
% xi2=0.45;omega_n2=5.73;xi3=0.9;omega_n3=3.42;Constant_tau_gc=100;eta_dt2=.71;Constant_Omega_g_m1=80;

%3rd fault scenario
Constant_Beta_1_m1=4;Gain_Beta_2_m2=1.8;Constant_Beta_3_m1=12;Constant_Omega_r_m1=1.2;Gain_Omega_r_m2=0.7;Gain_Omega_g_m2=1.7;
xi2=0.45;omega_n2=5.73;xi3=0.9;omega_n3=3.42;Constant_tau_gc=700;eta_dt2=0.3;Constant_Omega_g_m1=100;

%4th fault scenario
% Constant_Beta_1_m1=2;Gain_Beta_2_m2=3;Constant_Beta_3_m1=15;Constant_Omega_r_m1=1.7;Gain_Omega_r_m2=1.7;Gain_Omega_g_m2=0.7;
% xi2=0.45;omega_n2=5.73;xi3=0.9;omega_n3=3.42;Constant_tau_gc=-500;eta_dt2=.6;Constant_Omega_g_m1=130;

% %5th fault scenario
% Constant_Beta_1_m1=-5;Gain_Beta_2_m2=10;Constant_Beta_3_m1=3;Constant_Omega_r_m1=2.5;Gain_Omega_r_m2=0.9;Gain_Omega_g_m2=0.2;
% xi2=0.45;omega_n2=5.73;xi3=0.9;omega_n3=3.42;Constant_tau_gc=1500;eta_dt2=.25;Constant_Omega_g_m1=120;

%6th fault scenario
% Constant_Beta_1_m1=1;Gain_Beta_2_m2=2;Constant_Beta_3_m1=6;Constant_Omega_r_m1=0.5;Gain_Omega_r_m2=0.2;Gain_Omega_g_m2=2.8;
% xi2=0.45;omega_n2=5.73;xi3=0.9;omega_n3=3.42;Constant_tau_gc=2500;eta_dt2=0.1;Constant_Omega_g_m1=150;

% % %7th fault scenario
% Constant_Beta_1_m1=24;Gain_Beta_2_m2=5;Constant_Beta_3_m1=20;Constant_Omega_r_m1=3.5;Gain_Omega_r_m2=3;Gain_Omega_g_m2=5;
% xi2=0.25;omega_n2=8.73;xi3=0.95;omega_n3=1.42;Constant_tau_gc=900;eta_dt2=0.4;Constant_Omega_g_m1=50;

%8th fault scenario
% Constant_Beta_1_m1=-3;Gain_Beta_2_m2=5;Constant_Beta_3_m1=7;Constant_Omega_r_m1=2;Gain_Omega_r_m2=0.5;Gain_Omega_g_m2=1.5;
% xi2=0.45;omega_n2=5.73;xi3=0.9;omega_n3=3.42;Constant_tau_gc=-1000;eta_dt2=0.22;Constant_Omega_g_m1=100;

% %9th fault scenario
% Constant_Beta_1_m1=2;Gain_Beta_2_m2=3;Constant_Beta_3_m1=7;Constant_Omega_r_m1=2;Gain_Omega_r_m2=0.5;Gain_Omega_g_m2=1.5;
% xi2=0.45;omega_n2=5.73;xi3=0.9;omega_n3=3.42;Constant_tau_gc=800;eta_dt2=0.42;Constant_Omega_g_m1=110;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[Apb1,Bpb1,Cpb1,Dpb1]=tf2ss([omega_n2^2],[1 2*xi2*omega_n2 omega_n2^2]);
[Apb2,Bpb2,Cpb2,Dpb2]=tf2ss([omega_n3^2],[1 2*xi3*omega_n3 omega_n3^2]);
%
Addt2=[-(B_dt+B_r)/J_r B_dt/N_g/J_r -K_dt/J_r; eta_dt2*B_dt/N_g/J_g -(eta_dt2*B_dt/N_g^2+B_g)/J_g eta_dt2*K_dt/N_g/J_g; 1 -1/N_g 0];
Bddt2=[1/J_r 0; 0 -1/J_g;0 0];
Cddt2=[1 0 0;0 1 0];
Dddt2=[0 0;0 0];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Fault signals (choose the scenario of fault appearance)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %1st fault scenario
% D_Fault1=[Time' 1-[zeros(1,2000/Ts) ones(1,100/Ts) zeros(1,2300/Ts)]'];
% D_Fault2=[Time' 1-[zeros(1,2300/Ts) ones(1,100/Ts) zeros(1,2000/Ts)]'];
% D_Fault3=[Time' 1-[zeros(1,2600/Ts) ones(1,100/Ts) zeros(1,1700/Ts)]'];
% D_Fault4=[Time' 1-[zeros(1,1500/Ts) ones(1,100/Ts) zeros(1,2800/Ts)]'];
% D_Fault5=[Time' 1-[zeros(1,1000/Ts) ones(1,100/Ts) zeros(1,3300/Ts)]'];
% D_Fault6=[Time' 1-[zeros(1,2900/Ts) ones(1,100/Ts) zeros(1,1400/Ts)]'];
% D_Fault7=[Time' 1-[zeros(1,3400/Ts) (0:1/3000:2999/3000) ones(1,40/Ts) (2999/3000:-1/3000:0) zeros(1,900/Ts)]'];
% D_Fault8=[Time' 1-[zeros(1,3800/Ts) ones(1,100/Ts) zeros(1,500/Ts)]'];
% D_Fault9=[Time' 1-[zeros(1,4000/Ts) ones(1,200/Ts) zeros(1,200/Ts)]'];
% D_Fault10=[Time' 1-[zeros(1,3600/Ts) ones(1,100/Ts) zeros(1,700/Ts) ]'];

%%2nd fault scenario
% D_Fault1=[Time' 1-[zeros(1,1000/Ts) ones(1,100/Ts) zeros(1,3300/Ts)]'];
% D_Fault2=[Time' 1-[zeros(1,3900/Ts) ones(1,100/Ts) zeros(1,400/Ts)]'];
% D_Fault3=[Time' 1-[zeros(1,1600/Ts) ones(1,100/Ts) zeros(1,2700/Ts)]'];
% D_Fault4=[Time' 1-[zeros(1,500/Ts) ones(1,100/Ts) zeros(1,3800/Ts)]'];
% D_Fault5=[Time' 1-[zeros(1,4000/Ts) ones(1,100/Ts) zeros(1,300/Ts)]'];
% D_Fault6=[Time' 1-[zeros(1,2900/Ts) ones(1,100/Ts) zeros(1,1400/Ts)]'];
% D_Fault7=[Time' 1-[zeros(1,3400/Ts) (0:1/3000:2999/3000) ones(1,40/Ts) (2999/3000:-1/3000:0) zeros(1,900/Ts)]'];
% D_Fault8=[Time' 1-[zeros(1,1800/Ts) ones(1,100/Ts) zeros(1,2500/Ts)]'];
% D_Fault9=[Time' 1-[zeros(1,3000/Ts) ones(1,200/Ts) zeros(1,1200/Ts)]'];
% D_Fault10=[Time' 1-[zeros(1,2600/Ts) ones(1,100/Ts) zeros(1,1700/Ts) ]'];

%3nd fault scenario
D_Fault1=[Time' 1-[zeros(1,100/Ts) ones(1,100/Ts) zeros(1,4200/Ts)]'];
D_Fault2=[Time' 1-[zeros(1,2500/Ts) ones(1,100/Ts) zeros(1,1800/Ts)]'];
D_Fault3=[Time' 1-[zeros(1,900/Ts) ones(1,100/Ts) zeros(1,3400/Ts)]'];
D_Fault4=[Time' 1-[zeros(1,1200/Ts) ones(1,100/Ts) zeros(1,3100/Ts)]'];
D_Fault5=[Time' 1-[zeros(1,1700/Ts) ones(1,100/Ts) zeros(1,2600/Ts)]'];
D_Fault6=[Time' 1-[zeros(1,3200/Ts) ones(1,100/Ts) zeros(1,1100/Ts)]'];
D_Fault7=[Time' 1-[zeros(1,3400/Ts) (0:1/3000:2999/3000) ones(1,40/Ts) (2999/3000:-1/3000:0) zeros(1,900/Ts)]'];
D_Fault8=[Time' 1-[zeros(1,4200/Ts) ones(1,100/Ts) zeros(1,100/Ts)]'];
D_Fault9=[Time' 1-[zeros(1,2000/Ts) ones(1,200/Ts) zeros(1,2200/Ts)]'];
D_Fault10=[Time' 1-[zeros(1,3900/Ts) ones(1,100/Ts) zeros(1,400/Ts)]'];

% %4th fault scenario
% D_Fault1=[Time' 1-[zeros(1,250/Ts) ones(1,100/Ts) zeros(1,4050/Ts)]'];
% D_Fault2=[Time' 1-[zeros(1,3200/Ts) ones(1,100/Ts) zeros(1,1100/Ts)]'];
% D_Fault3=[Time' 1-[zeros(1,1000/Ts) ones(1,100/Ts) zeros(1,3300/Ts)]'];
% D_Fault4=[Time' 1-[zeros(1,800/Ts) ones(1,100/Ts) zeros(1,3500/Ts)]'];
% D_Fault5=[Time' 1-[zeros(1,2200/Ts) ones(1,100/Ts) zeros(1,2100/Ts)]'];
% D_Fault6=[Time' 1-[zeros(1,2800/Ts) ones(1,100/Ts) zeros(1,1500/Ts)]'];
% D_Fault7=[Time' 1-[zeros(1,3400/Ts) (0:1/3000:2999/3000) ones(1,40/Ts) (2999/3000:-1/3000:0) zeros(1,900/Ts)]'];
% D_Fault8=[Time' 1-[zeros(1,1200/Ts) ones(1,100/Ts) zeros(1,3100/Ts)]'];
% D_Fault9=[Time' 1-[zeros(1,3900/Ts) ones(1,200/Ts) zeros(1,300/Ts)]'];
% D_Fault10=[Time' 1-[zeros(1,4200/Ts) ones(1,100/Ts) zeros(1,100/Ts)]'];

% %5th fault scenario
% D_Fault1=[Time' 1-[zeros(1,3000/Ts) ones(1,100/Ts) zeros(1,1300/Ts)]'];
% D_Fault2=[Time' 1-[zeros(1,1300/Ts) ones(1,100/Ts) zeros(1,3000/Ts)]'];
% D_Fault3=[Time' 1-[zeros(1,500/Ts) ones(1,100/Ts) zeros(1,3800/Ts)]'];
% D_Fault4=[Time' 1-[zeros(1,200/Ts) ones(1,100/Ts) zeros(1,4100/Ts)]'];
% D_Fault5=[Time' 1-[zeros(1,2000/Ts) ones(1,100/Ts) zeros(1,2300/Ts)]'];
% D_Fault6=[Time' 1-[zeros(1,3600/Ts) ones(1,100/Ts) zeros(1,700/Ts)]'];
% D_Fault7=[Time' 1-[zeros(1,3400/Ts) (0:1/3000:2999/3000) ones(1,40/Ts) (2999/3000:-1/3000:0) zeros(1,900/Ts)]'];
% D_Fault8=[Time' 1-[zeros(1,2400/Ts) ones(1,100/Ts) zeros(1,1900/Ts)]'];
% D_Fault9=[Time' 1-[zeros(1,2700/Ts) ones(1,200/Ts) zeros(1,1500/Ts)]'];
% D_Fault10=[Time' 1-[zeros(1,3300/Ts) ones(1,100/Ts) zeros(1,1000/Ts)]'];

% % %6th fault scenario
% D_Fault1=[Time' 1-[zeros(1,3800/Ts) ones(1,100/Ts) zeros(1,500/Ts)]'];
% D_Fault2=[Time' 1-[zeros(1,2600/Ts) ones(1,100/Ts) zeros(1,1700/Ts)]'];
% D_Fault3=[Time' 1-[zeros(1,4000/Ts) ones(1,100/Ts) zeros(1,300/Ts)]'];
% D_Fault4=[Time' 1-[zeros(1,2900/Ts) ones(1,100/Ts) zeros(1,1400/Ts)]'];
% D_Fault5=[Time' 1-[zeros(1,600/Ts) ones(1,100/Ts) zeros(1,3700/Ts)]'];
% D_Fault6=[Time' 1-[zeros(1,1500/Ts) ones(1,100/Ts) zeros(1,2800/Ts)]'];
% D_Fault7=[Time' 1-[zeros(1,1800/Ts) (0:1/3000:2999/3000) ones(1,40/Ts) (2999/3000:-1/3000:0) zeros(1,2500/Ts)]'];
% D_Fault8=[Time' 1-[zeros(1,100/Ts) ones(1,100/Ts) zeros(1,4200/Ts)]'];
% D_Fault9=[Time' 1-[zeros(1,2000/Ts) ones(1,200/Ts) zeros(1,2200/Ts)]'];
% D_Fault10=[Time' 1-[zeros(1,3500/Ts) ones(1,100/Ts) zeros(1,800/Ts)]'];
%
% %7th fault scenario
% D_Fault1=[Time' 1-[zeros(1,4200/Ts) ones(1,100/Ts) zeros(1,100/Ts)]'];
% D_Fault2=[Time' 1-[zeros(1,4000/Ts) ones(1,100/Ts) zeros(1,300/Ts)]'];
% D_Fault3=[Time' 1-[zeros(1,3700/Ts) ones(1,100/Ts) zeros(1,600/Ts)]'];
% D_Fault4=[Time' 1-[zeros(1,3500/Ts) ones(1,100/Ts) zeros(1,800/Ts)]'];
% D_Fault5=[Time' 1-[zeros(1,3000/Ts) ones(1,100/Ts) zeros(1,1300/Ts)]'];
% D_Fault6=[Time' 1-[zeros(1,2700/Ts) ones(1,100/Ts) zeros(1,1600/Ts)]'];
% D_Fault7=[Time' 1-[zeros(1,2500/Ts) (0:1/3000:2999/3000) ones(1,40/Ts) (2999/3000:-1/3000:0) zeros(1,1800/Ts)]'];
% D_Fault8=[Time' 1-[zeros(1,2100/Ts) ones(1,100/Ts) zeros(1,2200/Ts)]'];
% D_Fault9=[Time' 1-[zeros(1,1400/Ts) ones(1,200/Ts) zeros(1,2800/Ts)]'];
% D_Fault10=[Time' 1-[zeros(1,1000/Ts) ones(1,100/Ts) zeros(1,3300/Ts)]'];

% % 8th fault scenario
% D_Fault1=[Time' 1-[zeros(1,100/Ts) ones(1,100/Ts) zeros(1,4200/Ts)]'];
% D_Fault2=[Time' 1-[zeros(1,500/Ts) ones(1,100/Ts) zeros(1,3800/Ts)]'];
% D_Fault3=[Time' 1-[zeros(1,900/Ts) ones(1,100/Ts) zeros(1,3400/Ts)]'];
% D_Fault4=[Time' 1-[zeros(1,1200/Ts) ones(1,100/Ts) zeros(1,3100/Ts)]'];
% D_Fault5=[Time' 1-[zeros(1,1700/Ts) ones(1,100/Ts) zeros(1,2600/Ts)]'];
% D_Fault6=[Time' 1-[zeros(1,2900/Ts) ones(1,100/Ts) zeros(1,1400/Ts)]'];
% D_Fault7=[Time' 1-[zeros(1,3400/Ts) (0:1/3000:2999/3000) ones(1,40/Ts) (2999/3000:-1/3000:0) zeros(1,900/Ts)]'];
% D_Fault8=[Time' 1-[zeros(1,4200/Ts) ones(1,100/Ts) zeros(1,100/Ts)]'];
% D_Fault9=[Time' 1-[zeros(1,300/Ts) ones(1,200/Ts) zeros(1,3900/Ts)]'];
% D_Fault10=[Time' 1-[zeros(1,2200/Ts) ones(1,100/Ts) zeros(1,2100/Ts)]'];

% 9th fault scenario

% D_Fault2=[Time' 1-[zeros(1,100/Ts) ones(1,10/Ts) zeros(1,5/Ts) ones(1,10/Ts) zeros(1,5/Ts) ones(1,10/Ts) zeros(1,5/Ts) ones(1,10/Ts) zeros(1,5/Ts) ones(1,10/Ts) zeros(1,2330/Ts) ones(1,10/Ts) zeros(1,5/Ts) ones(1,10/Ts) zeros(1,5/Ts) ones(1,10/Ts) zeros(1,5/Ts) ones(1,10/Ts) zeros(1,5/Ts) ones(1,10/Ts) zeros(1,1830/Ts) ]'];
% Offset_Beta_2_m2=[Time' [zeros(1,100/Ts) 0.5*ones(1,10/Ts) zeros(1,5/Ts) 1*ones(1,10/Ts) zeros(1,5/Ts) 1.5*ones(1,10/Ts) zeros(1,5/Ts) 2*ones(1,10/Ts) zeros(1,5/Ts) 3*ones(1,10/Ts) zeros(1,2330/Ts) 0.5*ones(1,10/Ts) zeros(1,5/Ts) 1*ones(1,10/Ts) zeros(1,5/Ts) 1.5*ones(1,10/Ts) zeros(1,5/Ts) 2*ones(1,10/Ts) zeros(1,5/Ts) 3*ones(1,10/Ts) zeros(1,1830/Ts) ]'];
%
% D_Fault5=[Time' 1-[zeros(1,500/Ts) ones(1,10/Ts) zeros(1,5/Ts) ones(1,10/Ts) zeros(1,5/Ts) ones(1,10/Ts) zeros(1,5/Ts) ones(1,10/Ts) zeros(1,5/Ts) ones(1,10/Ts) zeros(1,2330/Ts) ones(1,10/Ts) zeros(1,5/Ts) ones(1,10/Ts) zeros(1,5/Ts) ones(1,10/Ts) zeros(1,5/Ts) ones(1,10/Ts) zeros(1,5/Ts) ones(1,10/Ts) zeros(1,1430/Ts) ]'];
% Offset_Omega_r_m2=[Time' [zeros(1,500/Ts) 0.1*ones(1,10/Ts) zeros(1,5/Ts) 0.2*ones(1,10/Ts) zeros(1,5/Ts) 0.5*ones(1,10/Ts) zeros(1,5/Ts) 1*ones(1,10/Ts) zeros(1,5/Ts) 2*ones(1,10/Ts) zeros(1,2330/Ts) 0.1*ones(1,10/Ts) zeros(1,5/Ts) 0.2*ones(1,10/Ts) zeros(1,5/Ts) 0.5*ones(1,10/Ts) zeros(1,5/Ts) 1*ones(1,10/Ts) zeros(1,5/Ts) 2*ones(1,10/Ts) zeros(1,1430/Ts) ]'];
% Offset_Omega_g_m2=[Time' [zeros(1,500/Ts) 2.5*ones(1,10/Ts) zeros(1,5/Ts) 4*ones(1,10/Ts) zeros(1,5/Ts) 5*ones(1,10/Ts) zeros(1,5/Ts) 6*ones(1,10/Ts) zeros(1,5/Ts) 8*ones(1,10/Ts) zeros(1,2330/Ts) 2.5*ones(1,10/Ts) zeros(1,5/Ts) 4*ones(1,10/Ts) zeros(1,5/Ts) 5*ones(1,10/Ts) zeros(1,5/Ts) 6*ones(1,10/Ts) zeros(1,5/Ts) 8*ones(1,10/Ts) zeros(1,1430/Ts) ]'];

% D_Fault1=[Time' 1-[zeros(1,300/Ts) ones(1,100/Ts) zeros(1,4000/Ts)]'];
% D_Fault3=[Time' 1-[zeros(1,700/Ts) ones(1,100/Ts) zeros(1,3600/Ts)]'];
% D_Fault4=[Time' 1-[zeros(1,1000/Ts) ones(1,100/Ts) zeros(1,3300/Ts)]'];
% D_Fault6=[Time' 1-[zeros(1,1900/Ts) ones(1,100/Ts) zeros(1,2400/Ts)]'];
% D_Fault7=[Time' 1-[zeros(1,1400/Ts) (0:1/3000:2999/3000) ones(1,40/Ts) (2999/3000:-1/3000:0) zeros(1,2900/Ts)]'];
% D_Fault8=[Time' 1-[zeros(1,4200/Ts) ones(1,100/Ts) zeros(1,100/Ts)]'];
% D_Fault9=[Time' 1-[zeros(1,3300/Ts) ones(1,200/Ts) zeros(1,900/Ts)]'];
% D_Fault10=[Time' 1-[zeros(1,2200/Ts) ones(1,100/Ts) zeros(1,2100/Ts)]'];

% %without Fault
% D_Fault1=[Time' 1-[zeros(1,4400/Ts)]'];
% D_Fault2=[Time' 1-[zeros(1,4400/Ts)]'];
% D_Fault3=[Time' 1-[zeros(1,4400/Ts)]'];
% D_Fault4=[Time' 1-[zeros(1,4400/Ts)]'];
% D_Fault5=[Time' 1-[zeros(1,4400/Ts)]'];
% D_Fault6=[Time' 1-[zeros(1,4400/Ts)]'];
% D_Fault7=[Time' 1-[zeros(1,4400/Ts)]'];
% D_Fault8=[Time' 1-[zeros(1,4400/Ts)]'];
% D_Fault9=[Time' 1-[zeros(1,4400/Ts)]'];
%  D_Fault10=[Time' 1-[zeros(1,4400/Ts)]'];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Noise Seeds
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

min_seed = 0;
max_seed = 999;
% Generate random seeds for 13 'Random Number' Generators blocks
seed = min_seed + (max_seed-min_seed).*rand(13, 1);
% Round up to integer
seed = ceil(seed);

%% Simulation
open_system('FDI_measures.mdl');

tic
sim('FDI_measures.mdl');
toc

author: Sable

repo: mcbench-benchmarks


text:
function ft_progress(varargin)

% FT_PROGRESS shows a graphical or non-graphical progress indication similar to the
% standard WAITBAR function, but with the extra option of printing it in the command
% window as a plain text string or as a rotating dial. Alternatively, you can also
% specify it not to give feedback on the progress.
%
% Prior to the for-loop, you should call either
%   ft_progress('init', 'none',    'Please wait...')
%   ft_progress('init', 'text',    'Please wait...')
%   ft_progress('init', 'textbar', 'Please wait...')      % ascii progress bar
%   ft_progress('init', 'dial',    'Please wait...')      % rotating dial
%   ft_progress('init', 'etf',     'Please wait...')      % estimated time to finish
%   ft_progress('init', 'gui',     'Please wait...')
%
% In each iteration of the for-loop, you should call either
% ft_progress(x)                                          % only show percentage
% ft_progress(x, 'Processing event %d from %d', i, N)     % show string, x=i/N
%
% After finishing the for-loop, you should call
%   ft_progress('close')
%
% Here is an example for the use of a progress indicator
%   ft_progress('init', 'etf', 'Please wait...');
%   for i=1:100
%     ft_progress(i/100, 'Processing event %d from %d', i, 100);
%     pause(0.03);
%   end
%   ft_progress('close')
%
% See also WAITBAR

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

persistent p         % the previous value of the progress
persistent c         % counter for the number of updates that is done
persistent h         % the handle of the dialog (in case of type=gui)
persistent t0        % initial time, required for ETF
persistent p0        % initial percentage, required for ETF
persistent type      % type of feedback, string with none, gui, text, textcr, textnl
persistent title     % the string containing the title
persistent angle     % the angle in degrees, for dial or textbar
persistent strlen    % the length of the previously printed string, used to remove it by \b
persistent stopwatch % the time of previous invocation, used to restrict number of updates
persistent closing
persistent lastargin

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin>1 && ischar(varargin{1}) && strcmp(varargin{1}, 'init')
  % reset these to the defaults
  p = 0;
  c = 0;
  h = 0;
  t0 = [];
  p0 = [];
  type = 'none'; % set below
  title = ''; % set below
  angle = 0;
  strlen = 0;
  stopwatch = tic();
  lastargin = [];
  closing = false;

  % determine the type of feedback
  type = varargin{2};
  if strcmp(type, 'textcr') || strcmp(type, 'textnl')
    type = 'text';
  end
  % determine the title of the dialog
  if nargin>2
    title = varargin{3};
  else
    title = '';
  end
  switch type
    case 'gui'
      % initialise the waitbar dialog
      if ~isempty(title)
        h = waitbar(0, title);
      else
        h = waitbar(0, 'Please wait');
      end
    case {'text', 'textnl', 'textcr'}
      if ~isempty(title)
        % print the title to the screen and go to the next line
        fprintf('%s\n', title)
      end
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif nargin==1 && ischar(varargin{1}) && strcmp(varargin{1}, 'close')

  if ~isempty(lastargin)
    % the last input argument is used when ft_progress('close') is called but the
    % previous invocation was not processed yet due to the restriction in the number of
    % updates to once every 100 ms
    closing = true;
    ft_progress(lastargin{:});
  end

  switch type
    case 'gui'
      % close the waitbar dialog
      close(h);
    case {'text', 'etf', 'dial', 'textbar'}
      % finish by going to the next line
      fprintf('\n');
  end

  % reset these to the defaults
  p = 0;
  c = 0;
  h = 0;
  t0 = [];
  p0 = [];
  type = 'none';
  title = '';
  angle = 0;
  strlen = 0;
  stopwatch = [];
  lastargin = [];
  closing = false;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
else

  % make sure we don't update more than once every 100 ms, otherwise certain
  % conditions result in a significant performance hit
  if ~isempty(stopwatch) && toc(stopwatch) < 0.1 && ~closing
    lastargin = varargin;
    return;
  end
  stopwatch = tic();
  lastargin = [];

  if strcmp(type, 'dial')
    % display should always be updated for the dial
  elseif (varargin{1}-p)<0.01 && strcmp(type, 'gui')
    % display should not be updated it the difference is less than one percent
    return;
  elseif (varargin{1}-p)<0.01 && strcmp(type, 'textbar')
    % display should not be updated it the difference is less than one percent
    return;
  elseif (varargin{1}-p)<0.01 && strcmp(type, 'etf')
    % display should not be updated it the difference is less than one percent
    return;
  end

  % count the number of updates, for debugging
  c = c+1;

  % remember the current value for the next function call
  p = varargin{1};

  switch type
    case 'gui'
      % update the the length of the bar in the waitbar dialog
      waitbar(varargin{1}, h);

    case 'etf'
      % compute the estimated time that the computation still needs to finish
      if isempty(t0) || isempty(p0)
        t0 = clock;
        p0 = p;
      end
      elapsed = etime(clock, t0);
      if nargin>1 && ~isempty(varargin{2})

        varargin{2} = [repmat(sprintf('\b'), [1 strlen]) varargin{2}];

        part1 = sprintf(varargin{2:end});
        part2 = sprintf(' - estimated time to finish is %d seconds', round(elapsed*(1-p)/(p-p0)));

        fprintf([part1 part2]);

        % record actual string length that was printed (subtracting all the \b's)
        strlen = length(part1) + length(part2) - strlen;
      else
        % only print the estimated time to finish
        part1 = sprintf([repmat(sprintf('\b'), [1 strlen]) ' - estimated time to finish is %d seconds'], round(elapsed*(1-p)/(p-p0)));
        fprintf(part1);
        strlen = length(part1) - strlen + 1;
      end

    case 'dial'
      dial = '|/-\|/-\';
      if ~isempty(title)
        % print the title and draw a new hand of the rotating dial
        part1 = sprintf([repmat(sprintf('\b'), [1 strlen]) '%s %s'], title, dial(1+angle/45));
        strlen = length(title) + 2;
      else
        % draw a new hand of the rotating dial
        part1 = sprintf([repmat(sprintf('\b'), [1 strlen]) '%s'], dial(1+angle/45));
        strlen = 1;
      end
      if part1(end) == '\'
        % a single backslash is not valid, add another one
        part1(end+1) = '\';
      end
      fprintf(part1);
      % increment the angle with 45 degrees
      angle = angle + 45;
      if angle==360
        % reset the angle to 0 degrees
        angle = 0;
      end

    case 'textbar'
      dial = '|/-\|/-\';
      % construct the line looking like this [------/ ]
      len = 75 - length(title) - 3;
      len1 = round(p*len); % number of '-' characters before the dial
      len2 = len - len1; % number of ' ' characters after the dial
      line = [title, ' [' repmat('-', 1, len1), dial(1+angle/45), repmat(' ', 1, len2), ']'];
      if c~=1
        backline = repmat('\b', [1 length(line)]);
      else
        backline = '';
      end
      % don't use carriage return, it sometimes leads to a new line
      fprintf([backline '%s'], line);
      % increment the angle with 45 degrees
      angle = angle + 45;
      if angle==360
        % reset the angle to 0 degrees
        angle = 0;
      end

    case 'text'
      if nargin>1
        % ensure the string does not end with a newline or carriage return, since
        % either would break compatibility with a -nodesktop matlab environment
        if length(varargin{2})>1 && (all(varargin{2}((end-1):end) == '\r')...
            || all(varargin{2}((end-1):end) == '\n'))
          varargin{2} = varargin{2}(1:end-2);
        end

        varargin{2} = [repmat(sprintf('\b'), [1 strlen]) varargin{2}];
        if usejava('desktop')
          % a newline is appropriate when using the desktop environment
          varargin{2} = [varargin{2} '\n'];
        end

        part1 = sprintf(varargin{2:end});
        fprintf(part1);
        strlen = length(part1) - strlen;

      else
        part1 = sprintf([repmat(sprintf('\b'), [1 strlen]) '%6.2f %%'], 100*varargin{1});
        fprintf(part1);
        strlen = length(part1) - strlen;
      end

  end % case gui, etf, dial, textbar, text
end % updating the displayed value

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% some test code follows
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ft_progress('init', 'gui'); for i=1:100, ft_progress(i/100, '%d/%d\r', i, 100); pause(0.03); end; ft_progress('close')
% ft_progress('init', 'textbar'); for i=1:100, ft_progress(i/100, '%d/%d\r', i, 100); pause(0.03); end; ft_progress('close')
% ft_progress('init', 'dial'); for i=1:100, ft_progress(i/100, '%d/%d\r', i, 100); pause(0.03); end; ft_progress('close')
% ft_progress('init', 'none'); for i=1:100, ft_progress(i/100, '%d/%d\r', i, 100); pause(0.03); end; ft_progress('close')

% ft_progress('init', 'text'); for i=1:100, ft_progress(i/100, '%d/%d' , i, 100); pause(0.03); end; ft_progress('close')
% ft_progress('init', 'text'); for i=1:100, ft_progress(i/100, '%d/%d\n', i, 100); pause(0.03); end; ft_progress('close')
% ft_progress('init', 'text'); for i=1:100, ft_progress(i/100, '%d/%d\r', i, 100); pause(0.03); end; ft_progress('close')

% ft_progress('init', 'textnl'); for i=1:100, ft_progress(i/100, '%d/%d' , i, 100); pause(0.03); end; ft_progress('close')
% ft_progress('init', 'textnl'); for i=1:100, ft_progress(i/100, '%d/%d\n', i, 100); pause(0.03); end; ft_progress('close')
% ft_progress('init', 'textnl'); for i=1:100, ft_progress(i/100, '%d/%d\r', i, 100); pause(0.03); end; ft_progress('close')

% ft_progress('init', 'textcr'); for i=1:100, ft_progress(i/100, '%d/%d' , i, 100); pause(0.03); end; ft_progress('close')
% ft_progress('init', 'textcr'); for i=1:100, ft_progress(i/100, '%d/%d\n', i, 100); pause(0.03); end; ft_progress('close')
% ft_progress('init', 'textcr'); for i=1:100, ft_progress(i/100, '%d/%d\r', i, 100); pause(0.03); end; ft_progress('close')

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function a = wishart_unit_sample ( m, df )

%*****************************************************************************80
%
%% WISHART_UNIT_SAMPLE samples the unit Wishart distribution.
%
%  Discussion:
%
%    This function requires functions from the PDFLIB and RNGLIB libraries.
%
%    The "initialize()" function from RNGLIB must be called before using
%    this function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 October 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Patrick Odell, Alan Feiveson,
%    A numerical procedure to generate a sample covariance matrix,
%    Journal of the American Statistical Association,
%    Volume 61, Number 313, March 1966, pages 199-203.
%
%    Stanley Sawyer,
%    Wishart Distributions and Inverse-Wishart Sampling,
%    Washington University,
%    30 April 2007, 12 pages.
%
%  Parameters:
%
%    Input, integer M, the order of the matrix.
%
%    Input, integer DF, the number of degrees of freedom.
%    M <= DF.
%
%    Output, real A(M,M), the sample matrix from the unit Wishart distribution.
%
  if ( df < m )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'WISHART_UNIT_SAMPLE - Fatal error!\n' );
    fprintf ( 1, '  DF = %d < M = %d.\n', df, m );
    error ( 'WISHART_UNIT_SAMPLE - Fatal error!\n' );
  end

  c = zeros ( m, m );

  for i = 1 : m
    df_chi = df - i + 1;
    c(i,i) = sqrt ( r8_chi_sample ( df_chi ) );
    for j = i + 1 : m
      c(i,j) = r8_normal_01_sample ( );
    end
  end

  a = c' * c;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%
% Evaluate multivariate polynomial f at x
%
%    Syntax:   >> y = mvPolynEvaluate(f,x)
%
%     Input:    f -- multivariate polynomial (as a coeff. matrix)
%               x -- vector value for the variables
%
%    Output:  the value of f(x)

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [p,w, fVal, exitFlag] = fitBezierCurve(t,x,nPoint,tSpan, xBnd)
% [p,w, fVal, exitFlag] = fitBezierCurve(t,x,nPoint,tSpan, xBnd)
%
% This function fits a rational bezier curve to the data x(t), using a
% bezier curve with n points.
%
% INPUTS:
%   t = [1 x nTime] monotonically increasing time vector
%   x = [nCurve x nTime] vector function at each point in t
%   nPoint = number of points to use for bezier curve
%   tSpan = [1 x 2] domain of t
%   xBnd = [nCurve x 2] = [lower, upper] bound on function
%
% OUTPUTS:
%   p = [nCurve x nPoint] = bezier control points
%   w = [1 x nPoint] = weights for each control point
%   fVal = mean-square-error in function fit
%   exitFlag = fmincon exit flag
%
% See also: bezierCurve, rationalBezierCurve
%

nCurve = size(x,1);

tNormalized = (t-tSpan(1))/diff(tSpan);
wGuess = 0.5*ones(1,nPoint);
tGuess = linspace(0,1,nPoint);
pGuess = interp1(tNormalized',x',tGuess,'linear','extrap')';
zGuess = packDecVar(pGuess,wGuess);

wLow = zeros(1,nPoint);
wUpp = ones(1,nPoint);
pLow = xBnd(:,1)*ones(1,nPoint);
pUpp = xBnd(:,2)*ones(1,nPoint);

zLow = packDecVar(pLow,wLow);
zUpp = packDecVar(pUpp,wUpp);

options = optimset(...
    'Display','off');

problem.Aineq = []; problem.Aeq = [];
problem.bineq = []; problem.beq = [];
problem.lb = zLow;
problem.ub = zUpp;
problem.x0 = zGuess;
problem.options = options;
problem.objective = @(z)evalFit(z,t,x,nCurve,nPoint, tSpan);
problem.solver = 'fmincon';

[zSoln,fVal,exitFlag] = fmincon(problem);

[p,w] = unpackDecVar(zSoln,nCurve,nPoint);

end

function [p,w] = unpackDecVar(z,nCurve,nPoint)

w = z(1:nPoint)';

p = reshape(z(nPoint+1:end),nCurve,nPoint);

end

function z = packDecVar(p,w)

[nCurve, nPoint] = size(p);

z = [w'; reshape(p,nCurve*nPoint,1)];

end

function mse = evalFit(z,t,x,nCurve,nPoint,tSpan)

[p,w] = unpackDecVar(z,nCurve,nPoint);
xFit = rationalBezierCurve(p,w,t,tSpan);

mse = mean(mean((xFit-x).^2));

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function [s] = vector2struct( v, template )

  if isnumeric(template)

    s = reshape(v, template);

  elseif isfield(template, 'sz') && isfield(template, 'shapes')

    s = vector2cell(v, template);

  else
    count = 0;

    s = [];
    for fi = 1:length(template.fields)
      f = template.fields{fi};
      sz = template.shapes{fi};

      if isfield(sz, 'fields')

        sub_struct = vector2struct(v(count+1:end), sz);
        c = length(struct2vector(sub_struct));
        s.(f) = sub_struct;

      elseif isstruct(sz)

        c = sum(cellfun(@(x) prod(x), sz.shapes));
        if isempty(v)
          next = zeros(c,1);
        else
          next = v(count + [1:c]);
        end
        s.(f) = vector2cell(next, sz);

      else

        c = prod(sz);
        if isempty(v)
          next = zeros(c,1);
        else
          next = v(count + [1:c]);
        end
        s.(f) = reshape(next, sz);

      end

      count = count + c;
    end
  end
end

repo: CategoryShapes


text:
%%****************************************************************
%% orthbasis: Compute an orthonomal basis from
%%            the power basis {I,A,A^2, .... ,A^m}
%%            via an Arnoldi-type iteration.
%%
%%   [Q,H,C] = orthbasis(A,m);
%%
%%   Output:  Q = a cell array containing the orthonormal basis
%%                obtained from {I,A,A^2, .... ,A^m}.
%%
%% use in chebymat.m, igmres.m
%%
%% SDPT3: version 3.0
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last modified: 2 Feb 01
%%****************************************************************

  function [Q,H,C] = orthbasis(A,m);

  N = length(A);
  C(1,1) = 1/sqrt(N);
  Q = cell(1,m+1);
  Q{1} = eye(N)/sqrt(N);
  for k = 1:m
      V = A*Q{k};
      Vold = V;
      for j = 1:k
          H(j,k) = sum(sum(conj(Q{j}).*V));
          V = V - H(j,k)*Q{j};
      end;
      if (norm(V,'fro') < norm(Vold,'fro'));
         for j = 1:k
            s(j,1) = sum(sum(conj(Q{j}).*V));
            V = V - s(j)*Q{j};
         end;
         H(1:k,k) = H(1:k,k) + s;
      end;
      nrm = norm(V,'fro');
      H(k+1,k) = nrm;
      Q{k+1} = V/nrm;
      C(1:k+1,k+1) = ([0; C(1:k,k)] - [C(1:k,1:k)*H(1:k,k); 0])/nrm;
  end;
%%================================================================


text:
function tests = test_dictionary_comparison
  tests = functiontests(localfunctions);
end

function test_1(testCase)
    d = spx.dict.simple.dirac_fourier_mtx(16);
    ratio = spx.dict.comparison.matching_atoms_ratio(d, d);
    verifyEqual(testCase, ratio, 1);
    % Let's mess up with one of the atoms
    d2 = d;
    d2(1:4, 1) = 1/2;
    ratio = spx.dict.comparison.matching_atoms_ratio(d, d2);
    verifyEqual(testCase, ratio, 31/32);
    d2(2:5, 2) = 1/2;
    ratio = spx.dict.comparison.matching_atoms_ratio(d, d2);
    verifyEqual(testCase, ratio, 30/32);
end

author: indigits

repo: sparse-plex


text:
function n = numberFrames(view,scan)
%
%    n = numberFrames(view,scan)
%
% Returns the number of frames of a selected scan
%
% gb 03/30/2005
if ieNotDefined('scan')
    scan = viewGet(view,'curScan');
end

global dataTYPES;
n = dataTYPES(view.curDataType).scanParams(scan).nFrames;

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function dicom_write_volume(Volume,filename,volscale,info)
% This function DICOM_WRITE_VOLUME will write a Matlab 3D volume as
% a stack of 2D slices in separate dicom files.
%
% dicom_write_volume(Volume,Filename,Scales,Info)
%
% inputs,
%   Volume: The 3D Matlab volume
%   Filename: The name of the dicom files
%   Scales: The dimensions of every voxel/pixel
%   Info: A struct with dicom tags and values
%
% Function is written by D.Kroon University of Twente (May 2009)

% Check inputs
if(exist('filename','var')==0), filename=[]; end
if(exist('info','var')==0), info=[]; end
if(exist('volscale','var')==0), volscale=[1 1 1]; end

% Show file dialog if no file name specified
if(isempty(filename))
    [filename, pathname] = uiputfile('*.dcm', 'Save to Dicom');
    filename= [pathname filename];
end

% Add dicom tags to info structure
if(~isstruct(info))
    info=struct;
    % Make random series number
    SN=round(rand(1)*1000);
    % Get date of today
    today=[datestr(now,'yyyy') datestr(now,'mm') datestr(now,'dd')];
    info.SeriesNumber=SN;
    info.AcquisitionNumber=SN;
    info.StudyDate=today;
    info.StudyID=num2str(SN);
    info.PatientID=num2str(SN);
    info.PatientPosition='HFS';
    info.AccessionNumber=num2str(SN);
    info.StudyDescription=['StudyMAT' num2str(SN)];
    info.SeriesDescription=['StudyMAT' num2str(SN)];
    info.Manufacturer='Matlab Convert';
    info.SliceThickness=volscale(3);
    info.PixelSpacing=volscale(1:2);
    info.SliceLocation=0;
end

% Remove filename extention
pl=find(filename=='.'); if(~isempty(pl)), filename=filename(1:pl-1); end

% Read Volume data
disp('Writing Dicom Files...');
for slicenum=1:size(Volume,3)
    filenamedicom=[filename number2string(slicenum) '.dcm'];
    % Add slice specific dicom info
    info.InstanceNumber = slicenum;
    info.SliceLocation = info.SliceLocation+volscale(3);
    % Write the dicom file
    disp(['Writing : ' filenamedicom]);
    dicomwrite(Volume(:,:,slicenum), filenamedicom, info)
end

function numstr=number2string(num)
    num=num2str(num);
    numzeros='000000';
    numstr=[numzeros(length(num):end) num];

author: tomdoel

repo: pulmonarytoolkit


text:
function [X,rho,eta,F] = nu(A,b,k,nu,s)
%NU Brakhage's nu-method.
%
% [X,rho,eta,F] = nu(A,b,k,nu,s)
%
% Performs k steps of Brakhage's nu-method for the problem
%    min || A x - b || .
% The routine returns all k solutions, stored as columns of
% the matrix X.  The solution norm and residual norm are returned
% in eta and rho, respectively.
%
% If nu is not specified, nu = .5 is the default value, which gives
% the Chebychev method of Nemirovskii and Polyak.
%
% If the singular values s are also provided, nu computes the
% filter factors associated with each step and stores them
% columnwise in the matrix F.

% Reference: H. Brakhage, "On ill-posed problems and the method of
% conjugate gradients"; in H. W. Engl & G. W. Groetsch, "Inverse and
% Ill-Posed Problems", Academic Press, 1987.

% Martin Hanke, Institut fuer Praktische Mathematik, Universitaet
% Karlsruhe and Per Christian Hansen, IMM, 03/21/92.

% Set parameter.
l_steps = 3;      % Number of Lanczos steps for est. of || A ||.
fudge   = 0.99;   % Scale A and b by fudge/|| A*L_p ||.
fudge_thr = 1e-4; % Used to prevent filter factors from exploding.

% Initialization.
if (k < 1), error('Number of steps k must be positive'), end
if (nargin==3), nu = .5; end
[m,n] = size(A); X = zeros(n,k);
if (nargout > 1)
  rho = zeros(k,1); eta = rho;
end;
if (nargin==5)
  F = zeros(n,k); Fd = zeros(n,1); s = s.^2;
end
V = zeros(n,l_steps); B = zeros(l_steps+1,l_steps);
v = zeros(n,1); eta = zeros(l_steps+1,1);

% Compute a rough estimate of the norm of A by means of a few
% steps of Lanczos bidiagonalization, and scale A and b such
% that || A || is slightly less than one.
beta = norm(b); u = b/beta;
for i=1:l_steps
  r = A'*u - beta*v;
  alpha = norm(r); v = r/alpha;
  B(i,i) = alpha; V(:,i) = v;
  p = A*v - alpha*u;
  beta = norm(p); u = p/beta;
  B(i+1,i) = beta;
end
scale = fudge/norm(B); A = scale*A; b = scale*b;
if (nargin==5), s = scale^2*s; end

% Prepare for iteration.
x = zeros(n,1);
d = A'*b;
r = d;
if (nargout>1), z = b; end

% Iterate.
for j=0:k-1

  % Updates.
  alpha = 4*(j+nu)*(j+nu+0.5)/(j+2*nu)/(j+2*nu+0.5);
  beta  = (j+nu)*(j+1)*(j+0.5)/(j+2*nu)/(j+2*nu+0.5)/(j+nu+1);
  Ad = A*d; AAd = A'*Ad;
  x  = x + alpha*d;
  r  = r - alpha*AAd;
  d  = r + beta*d;
  X(:,j+1) = x;
  if (nargout>1)
    z = z - alpha*Ad; rho(j+1) = norm(z)/scale;
  end;
  if (nargout>2), eta(j+1) = norm(x); end;

  % Filter factors.
  if (nargin==5)
    if (j==0)
      F(:,1) = alpha*s;
      Fd = s - s.*F(:,1) + beta*s;
    else
      F(:,j+1) = F(:,j) + alpha*Fd;
      Fd = s - s.*F(:,j+1) + beta*Fd;
    end
    if (j > 1)
      f = find(abs(F(:,j)-1) < fudge_thr & abs(F(:,j-1)-1) < fudge_thr);
      if ~isempty(f), F(f,j+1) = ones(length(f),1); end
    end
  end

end

repo: REST


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nick.battista@unc.edu
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting %%	lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us %to add a specific muscle model, please let Nick (nick.battista@unc.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the RUBBERBAND-EXAMPLE geometry and prints associated input files
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Dear_KC()

%
% Grid Parameters (MAKE SURE MATCHES IN input2d !!!)
%
Nx =  32;        % # of Eulerian Grid Pts. in x-Direction (MUST BE EVEN!!!)
Ny =  32;        % # of Eulerian Grid Pts. in y-Direction (MUST BE EVEN!!!)
Lx = 1.0;        % Length of Eulerian Grid in x-Direction
Ly = 1.0;        % Length of Eulerian Grid in y-Direction
ds = Lx/(2*Nx);  % Spatial step!

% Immersed Structure Geometric / Dynamic Parameters %
N = 2*Nx;        % Number of Lagrangian Pts. (2x resolution of Eulerian grid)
ds_Rest = 0;     % Resting length of springs
struct_name = 'dear_KC'; % Name for .vertex, .spring, etc files.

% Line on top and bottom of geometry
ds2 = 0.125*ds;
xLine = Lx/20:ds2:0.95*Lx;
yLineB = 0.025*Lx*ones(1,length(xLine));
yLineT= 0.975*Lx*ones(1,length(xLine));

xSq = [xLine xLine yLineB yLineT];
ySq = [yLineB yLineT xLine xLine];

% Call function to construct geometry
[x1,y1] = give_Me_Immsersed_Boundary_Geometry_1(Lx,Nx,ds);
x1 = [x1 xSq];
y1 = [y1 ySq];

% Call function to construct geometry
[x2,y2] = give_Me_Immsersed_Boundary_Geometry_2(Lx,Nx,ds,x1);
x2 = [x2 xSq];
y2 = [y2 ySq];

% Call function to construct geometry
[x3,y3] = give_Me_Immsersed_Boundary_Geometry_3(Lx,Nx,ds);
x3 = [x3 xSq];
y3 = [y3 ySq];

% Plot Geometry to test BEFORE taking out pts.
figure(1)
plot(x1,y1,'r*'); hold on;
axis([0 Lx 0 Ly]);

figure(2)
plot(x2,y2,'b*'); hold on;
axis([0 Lx 0 Ly]);

figure(3)
plot(x3,y3,'k*'); hold on;
axis([0 Lx 0 Ly]);

% Print files to .txt files
please_Print_Vertices_To_File(x1,y1,x2,y2,x3,y3)

% Prints .vertex file!
print_Lagrangian_Vertices(x1,y1,struct_name);

% Prints .spring file!
%k_Spring = 1e7;
%print_Lagrangian_Springs(xLag,yLag,k_Spring,ds_Rest,struct_name);

% Prints .beam file!
%k_Beam = 0.5; C = 0.0;
%print_Lagrangian_Beams(xLag,yLag,k_Beam,C,struct_name);

% Prints .target file!
k_Target = 1e6;
print_Lagrangian_Target_Pts(x1,k_Target,struct_name);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints VERTEX points to a file called struct.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Vertices(xLag,yLag,struct_name)

    N = length(xLag);

    vertex_fid = fopen([struct_name '.vertex'], 'w');

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

    %Loops over all Lagrangian Pts.
    for s = 1:N
        X_v = xLag(s);
        Y_v = yLag(s);
        fprintf(vertex_fid, '%1.16e %1.16e\n', X_v, Y_v);
    end

    fclose(vertex_fid);

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

function print_Lagrangian_Target_Pts(xLag,k_Target,struct_name)

    N = length(xLag);

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

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

    %Loops over all Lagrangian Pts.
    for s = 1:N
        if s > 284
            k_Target = 2.5e9;
        end
        fprintf(target_fid, '%d %1.16e\n', s, k_Target);
    end

    fclose(target_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints BEAM (Torsional Spring) points to a file called struct.beam
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Beams(xLag,yLag,k_Beam,C,struct_name)

    % k_Beam: beam stiffness
    % C: beam curvature

    N = length(xLag); % NOTE: Total number of beams = Number of Total Lag Pts. - 2

    beam_fid = fopen([struct_name '.beam'], 'w');

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

    %spring_force = kappa_spring*ds/(ds^2);

    %BEAMS BETWEEN VERTICES
    for s = 2:N-1
            if  s <= N-1
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, s+1, k_Beam, C);
            else
                %Case s=N
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, 1,   k_Beam, C);
            end
    end
    fclose(beam_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints SPRING points to a file called struct.spring
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Springs(xLag,yLag,k_Spring,ds_Rest,struct_name)

    N = length(xLag);

    spring_fid = fopen([struct_name '.spring'], 'w');

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

    %spring_force = kappa_spring*ds/(ds^2);

    %SPRINGS BETWEEN VERTICES
    for s = 1:N
            if s < N
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+1, k_Spring, ds_Rest);
            else
                %Case s=N
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, 1,   k_Spring, ds_Rest);
            end
    end
    fclose(spring_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the Lagrangian structure geometry for PHASE 1
%           msg: "Hi KC!!"
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag] = give_Me_Immsersed_Boundary_Geometry_1(Lx,Nx,ds)

% The immsersed structure message #1: Hi Nghieng!! %
len = Lx/6.5;

xC = 0;%-Lx/4;
yC = 0;%-Lx/2;
[xH,yH] = give_Me_The_Letter_Please_2(ds,len,'H',xC,yC); % (0,0)
[xi1,yi1] = give_Me_The_Letter_Please_2(ds,len,'i',xC,yC-0.05);
%
[xN,yN] = give_Me_The_Letter_Please_2(ds,len,'N',xC,yC);
[xg1,yg1] = give_Me_The_Letter_Please_2(ds,len,'g',xC,yC);
[xh,yh] = give_Me_The_Letter_Please_2(ds,len,'h',xC,yC);
[xi2,yi2] = give_Me_The_Letter_Please_2(ds,len,'i',xC,yC-0.05);
[xe,ye] = give_Me_The_Letter_Please_2(ds,len,'e',xC,yC);
[xn,yn] = give_Me_The_Letter_Please_2(ds,len,'n',xC,yC);
[xg2,yg2] = give_Me_The_Letter_Please_2(ds,len,'g',xC,yC);

% No room in domain for !'s
%[xEx1,yEx1] = give_Me_The_Letter_Please_2(ds,1.2*len,'!',xC-0.55,yC);
%[xEx2,yEx2] = give_Me_The_Letter_Please_2(ds,1.2*len,'!',xC-0.59,yC);

% HI
xH = xH + 0.45;   yH = yH + 0.7;
xi1 = xi1 + 0.55; yi1 = yi1 + 0.65;

% Nghieng
xN =  xN + 0.15;     yN = yN + 0.375;
xg1 = xg1 + 0.325;   yg1 = yg1 + 0.375;
xh =  xh + 0.435;      yh = yh + 0.375;
xi2 = xi2 + 0.515;     yi2 = yi2 + 0.325;
xe = xe + 0.59;       ye = ye + 0.375;
xn = xn + 0.7;       yn = yn + 0.375;
xg2 = xg2 + 0.8;     yg2 = yg2 + 0.375;

%
% COMBINE GEOMETRY
%
xLag = [xH xi1 xN xg1 xh xi2 xe xn xg2];
yLag = [yH yi1 yN yg1 yh yi2 ye yn yg2];

%
% Test Plot
%
% figure(1)
% plot(xH,yH,'.'); hold on;
% plot(xi1,yi1,'.'); hold on;
% plot(xN,yN,'.'); hold on;
% plot(xg1,yg1,'.'); hold on;
% plot(xh,yh,'.'); hold on;
% plot(xi2,yi2,'.'); hold on;
% plot(xe,ye,'.'); hold on;
% plot(xn,yn,'.'); hold on;
% plot(xg2,yg2,'.'); hold on;

%
% Construct Horizontal Lines for Added Pts. in Domain
%
xRow = (0:ds:Lx/3) + 1/3;
yRow = ones(size(xRow));

xRow2 = (0:ds:Lx/6+ds) + 0.41666;
yRow2 = ones(size(xRow2));

xRow_T1 = xRow; yRow_T1 = 0.925*yRow;
xRow_T2 = xRow2; yRow_T2 = 0.85*yRow2;

xRow_B1 = [xRow xRow(end)+ds]; yRow_B1 = 0.075*[yRow 1];
xRow_B2 = xRow2; yRow_B2 = 0.15*yRow2;

%
% Combine Letters + Horizontal Lines
%
xLag = [xLag xRow_T1 xRow_T2 xRow_B1 xRow_B2];
yLag = [yLag yRow_T1 yRow_T2 yRow_B1 yRow_B2];

%
% Test Plot 2
%
% figure(11)
% plot(xLag,yLag,'.'); hold on;
% axis([0 1 0 1]);

%
% CHECK TO MAKE SURE EQUAL NUMBER OF MOVING LAGS AS OTHER PHASES
% HERE: 284
%
%size(xLag)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the Lagrangian structure geometry for PHASE 2
%           msg: "Would you like to ..."
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag] = give_Me_Immsersed_Boundary_Geometry_2(Lx,Nx,ds,x1)

% The immsersed structure message #2: Woud you like to... %
len = Lx/8;

xC = -Lx/6;
yC = -2*Lx/3;
[xW,yW] = give_Me_The_Letter_Please(ds,len,'W',xC,yC);
[xo,yo] = give_Me_The_Letter_Please(ds,len,'o',xC-0.1,yC);
[xu,yu] = give_Me_The_Letter_Please(ds,len,'u',xC-0.18,yC);
[xl,yl] = give_Me_The_Letter_Please(ds,len,'l',xC-0.27,yC);
[xd,yd] = give_Me_The_Letter_Please(ds,1.2*len,'d',xC-0.29,yC);

[xy,yy] = give_Me_The_Letter_Please(ds,1.2*len,'y',xC-0.47,yC);
[xo2,yo2] = give_Me_The_Letter_Please(ds,len,'o',xC-0.55,yC);
[xu2,yu2] = give_Me_The_Letter_Please(ds,len,'u',xC-0.63,yC);

xC = -Lx/4;
yC = -1*Lx/4 - Lx/6;
[xl2,yl2] = give_Me_The_Letter_Please(ds,len,'l',xC,yC);
[xi,yi] = give_Me_The_Letter_Please(ds,len,'i',xC-0.01,yC);
[xk,yk] = give_Me_The_Letter_Please(ds,len,'k',xC-0.07,yC);
[xe,ye] = give_Me_The_Letter_Please(ds,len,'e',xC-0.13,yC);

[xt,yt] = give_Me_The_Letter_Please(ds,len,'t',xC-0.27,yC);
[xo3,yo3] = give_Me_The_Letter_Please(ds,len,'o',xC-0.32,yC);

yC = -1*Lx/5 - Lx/6;
[xo4,yo4] = give_Me_The_Letter_Please(ds/3,len/6,'o',xC-0.44,yC);
[xo5,yo5] = give_Me_The_Letter_Please(ds/3,len/6,'o',xC-0.48,yC);
[xo6,yo6] = give_Me_The_Letter_Please(ds/3,len/6,'o',xC-0.52,yC);

xLag = [xW xo xu xl xd xy xo2 xu2 xl2 xi xk xe xt xo3 xo4 xo5 xo6];
yLag = [yW yo yu yl yd yy yo2 yu2 yl2 yi yk ye yt yo3 yo4 yo5 yo6];

%x1Len = length(x1)
%xLagLen=length(xLag)
%PtsLeft = abs( length(xLag) - length(x1) )
%ds1 = (Lx/2.5)/( PtsLeft/2 - 1 );

xBL = 0:ds:Lx/2.5; xBL = xBL + 3/10*Lx;
yTL = 0.84*ones(1,length(xBL));
yBL = 0.18*ones(1,length(xBL));

xLag = [xLag xBL xBL];
yLag = [yLag yTL yBL];

% plot(xW,yW,'r*'); hold on;
% plot(xo,yo,'r*'); hold on;
% plot(xu,yu,'r*'); hold on;
% plot(xl,yl,'r*'); hold on;
% plot(xd,yd,'r*'); hold on;
%
% plot(xy,yy,'r*'); hold on;
% plot(xo2,yo2,'r*'); hold on;
% plot(xu2,yu2,'r*'); hold on;
%
% plot(xl2,yl2,'r*'); hold on;
% plot(xi,yi,'r*'); hold on;
% plot(xk,yk,'r*'); hold on;
% plot(xe,ye,'r*'); hold on;
%
% plot(xt,yt,'r*'); hold on;
% plot(xo3,yo3,'r*'); hold on;
%
% plot(xo4,yo4,'r*'); hold on;
% plot(xo5,yo5,'r*'); hold on;
% plot(xo6,yo6,'r*'); hold on;
% axis([0 Lx 0 Lx]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the Lagrangian structure geometry for PHASE 2
%           msg: "Would you like to ..."
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag] = give_Me_Immsersed_Boundary_Geometry_3(Lx,Nx,ds)

% The immsersed structure message #3: ...go on a date with me? %
len = Lx/9;

xC = -Lx/8;
yC = -3*Lx/4;
shift = 0.2;

[xo4,yo4] = give_Me_The_Letter_Please(ds/3,len/6,'o',xC,yC+0.09);
[xo5,yo5] = give_Me_The_Letter_Please(ds/3,len/6,'o',xC-0.04,yC+0.09);
[xo6,yo6] = give_Me_The_Letter_Please(ds/3,len/6,'o',xC-0.08,yC+0.09);

[xg,yg] = give_Me_The_Letter_Please(ds,len,'g',xC-shift,yC);
[xo,yo] = give_Me_The_Letter_Please(ds,len,'o',xC-0.1-shift,yC);

[xo2,yo2] = give_Me_The_Letter_Please(ds,1.2*len,'o',xC-0.25-shift,yC);
[xn,yn] = give_Me_The_Letter_Please(ds,1.2*len,'n',xC-0.34-shift,yC);

[xa,ya] = give_Me_The_Letter_Please(ds,len,'a',xC-0.48-shift,yC);

xC = -Lx/5;
yC = -Lx/2;
[xd,yd] = give_Me_The_Letter_Please(ds,len,'d',xC,yC);
[xa2,ya2] = give_Me_The_Letter_Please(ds,len,'a',xC-0.08,yC);
[xt,yt] = give_Me_The_Letter_Please(ds,len,'t',xC-0.16,yC);
[xe,ye] = give_Me_The_Letter_Please(ds,len,'e',xC-0.22,yC);

[xw,yw] = give_Me_The_Letter_Please(ds,1.1*len,'w',xC-0.37,yC);
[xi,yi] = give_Me_The_Letter_Please(ds,len,'i',xC-0.45,yC);
[xt2,yt2] = give_Me_The_Letter_Please(ds,len,'t',xC-0.505,yC);
[xh,yh] = give_Me_The_Letter_Please(ds,len,'h',xC-0.58,yC);

xC = -Lx/2.5;
yC = -Lx/4;
[xm,ym] = give_Me_The_Letter_Please(ds,1.1*len,'m',xC,yC);
[xe2,ye2] = give_Me_The_Letter_Please(ds,len,'e',xC-0.12,yC);
[xQU,yQU] = give_Me_The_Letter_Please(ds,1.2*len,'?',xC-0.22,yC);

xLag = [xo4 xo5 xo6 xg xo xo2 xn xa xd xa2 xt xe xw xi xt2 xh xm xe2 xQU];
yLag = [yo4 yo5 yo6 yg yo yo2 yn ya yd ya2 yt ye yw yi yt2 yh ym ye2 yQU];

% plot(xo4,yo4,'r*'); hold on;
% plot(xo5,yo5,'r*'); hold on;
% plot(xo6,yo6,'r*'); hold on;
%
% plot(xg,yg,'r*'); hold on;
% plot(xo,yo,'r*'); hold on;
%
% plot(xo2,yo2,'r*'); hold on;
% plot(xn,yn,'r*'); hold on;
%
% plot(xa,ya,'r*'); hold on;
%
% plot(xd,yd,'r*'); hold on;
% plot(xa2,ya2,'r*'); hold on;
% plot(xt,yt,'r*'); hold on;
% plot(xe,ye,'r*'); hold on;
%
% plot(xw,yw,'r*'); hold on;
% plot(xi,yi,'r*'); hold on;
% plot(xt2,yt2,'r*'); hold on;
% plot(xh,yh,'r*'); hold on;
%
% plot(xm,ym,'r*'); hold on;
% plot(xe2,ye2,'r*'); hold on;
% plot(xQU,yQU,'r*'); hold on;
%
% axis([0 Lx 0 Lx]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%s%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the Lagrangian structure geometry for the Heart
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag] = give_Me_Immsersed_Boundary_Geometry_Heart(Lx,Nx,frac,Nb)

% The immsersed structure is a heart %
ds = 2*pi/(Nb-1);
tVec = 0:ds:2*pi;

for i=1:length(tVec)
    t = tVec(i);
    xLag(i)	=	16*sin(t)^3;
    yLag(i)	=	13*cos(t)-5*cos(2*t)-2*cos(3*t)-cos(4*t);
end

xLag = frac*xLag + 0.5;
yLag = frac*yLag + 0.5;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the Lagrangian structure geometry for PHASE 2
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag] = give_Me_Immsersed_Boundary_Geometry_4(Lx,Nx,frac,Nb)

% The immsersed structure is a heart %
ds = 2*pi/(Nb-1);
tVec = 0:ds:2*pi;
for i=1:length(tVec)
    t = tVec(i);
    r(i) = 1 - sin(t);
    xLag(i)	=	r(i)*cos(t);
    yLag(i)	=	r(i)*sin(t);
end

xLag = frac*xLag;
yLag = frac*yLag;

minY = min(yLag);
yLag = yLag - 3*minY/8;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: take out Lagrangian Pts.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x1,y1,x2,y2] = please_Take_Out_Points(x1,y1,x2,y2)

n1 = 50;
n2 = 70;

x1 = [x1(1:n1) x1(n2:end)];
y1 = [y1(1:n1) y1(n2:end)];

x2 = [x2(1:n1) x2(n2:end)];
y2 = [y2(1:n1) y2(n2:end)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints all Vertices to File
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function please_Print_Vertices_To_File(X1,Y1,X2,Y2,X3,Y3)

fileID = fopen('All_Positions.txt','w');
for j=1:length(X1)
    fprintf(fileID,'%1.16e %1.16e %1.16e %1.16e %1.16e %1.16e %1.16e %1.16e\n', X1(j),Y1(j),X2(j),Y2(j),X3(j),Y3(j),X3(j),Y3(j) );
end
fclose(fileID);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: give me the letter!!!!!
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x,y] = give_Me_The_Letter_Please(ds,len,letter,xC,yC)

if strcmp(letter,'a')
    r = len/4.1;
    dt = ds/r;
    theta = 0:dt:2*pi;
    for i=1:length(theta)
       xA(i) = r*cos(theta(i)) - xC;
       yA(i) = r*sin(theta(i)) - yC - len/4;
    end

    yLine = [-r:ds:r r]; yLine = yLine - yC - len/4;
    xLine = r*ones(1,length(yLine)) - (xC-ds/6);

    x = [xA xLine];
    y = [yA yLine];

elseif strcmp(letter,'C')

    r = len/2;
        dt = ds/r;
    theta = pi/2:dt:pi-0.25*dt;
    theta = [theta pi -theta];
    for i=1:length(theta)
       x(i) = 3*r/4*cos(theta(i)) - xC;
       y(i) = r*sin(theta(i)) - yC;
    end

    xTB = 0:3*ds/4:len/6;
    xTB = xTB - xC + 0.01;
    yT = r*ones(1,length(xTB))-yC;
    yB = -r*ones(1,length(xTB))-yC;

    x = [xTB x xTB];
    y = [yT y yB];

    length(x)

elseif strcmp(letter,'d')

    r = len/4.1;
    dt = ds/r;
    theta = 0:dt:2*pi;
    for i=1:length(theta)
       xA(i) = r*cos(theta(i)) - xC;
       yA(i) = r*sin(theta(i)) - yC - len/4;
    end

    yLine = [-len/2:ds:len/2.25 r]; yLine = yLine - yC;
    xLine = r*ones(1,length(yLine)) - (xC-ds/6);

    x = [xA xLine];
    y = [yA yLine];

 elseif strcmp(letter,'e')

    r = len/4.1;
    yV = [-r:ds:r r];
    xH = yV - xC;
    yH = zeros(1,length(xH)) - yC;
    yH2 = yH - len/4;
    yH3 = yH - len/2;

    yV = yV - yC - len/4;
    xV = r*ones(1,length(yV));

    yV3 = 0:-ds:-len/4;
    yV3 = yV3 - yC;
    xV3 = r*ones(1,length(yV3));

    x = [-xV-(xC-ds/6) xH xV3-(xC-ds/6) xH xH];
    y = [yV yH yV3 yH2 yH3];

elseif strcmp(letter,'g')

    r = len/4.1;
    dt = ds/r;
    theta = 0:dt:2*pi;
    for i=1:length(theta)
       xA(i) = r*cos(theta(i)) - xC;
       yA(i) = r*sin(theta(i)) - yC - len/4;
    end

    yLine = [-len/2:ds:len/2.25 r]; yLine = yLine - yC - len/3;
    xLine = r*ones(1,length(yLine)) - (xC-ds/6);

    yV = [-r:ds:r r];
    xH = yV - xC;
    yH = zeros(1,length(xH)) - yC - 0.85*len;

    x = [xA xLine xH];
    y = [yA yLine yH];

elseif strcmp(letter,'h')

   yLine = -len/2:ds:0-ds/4;
   yLine = [yLine 0 -yLine];
   yLine1 = yLine - yC;
   xLine1 = 1/2*len/2*ones(1,length(yLine1));

   yLine2 = -len/2:ds:0-ds/4;
   yLine2 = [yLine2 0];
   yLine2 = yLine2 - yC;
   xLine2 = 1/2*len/2*ones(1,length(yLine2));

   xH = (-1/2*len/2+ds/2:ds:1/2*len/2-ds/4);
   xH = xH - xC;
   yH = zeros(1,length(xH)) - yC;

   x = [-xLine1-xC xH xLine2-xC];
   y = [yLine1 yH yLine2];

elseif strcmp(letter,'H')
   yLine = -len/2:ds:0-ds/4;
   yLine = [yLine 0 -yLine];
   yLine = yLine - yC;
   xLine = 1/2*len/2*ones(1,length(yLine));
   xH = (-1/2*len/2+ds/2:ds:1/2*len/2-ds/4);
   xH = xH - xC;
   yH = zeros(1,length(xH)) - yC;

   x = [-xLine-xC xH xLine-xC];
   y = [yLine yH yLine];

elseif strcmp(letter,'i')

   yLine = -len/2:ds:0-ds/4;
   yLine = yLine - yC;
   xLine = zeros(1,length(yLine)) - xC;

    r = len/21;
    dt = ds/(2*r);
    theta = 0:dt:2*pi;
    for i=1:length(theta)
       xDot(i) = r*cos(theta(i)) - xC;
       yDot(i) = r*sin(theta(i)) - yC + len/7;
    end

    x = [xLine xDot];
    y = [yLine yDot];

elseif strcmp(letter,'k')

   yLine = -len/2:ds:len/2;
   yLine = yLine - yC;
   xLine = zeros(1,length(yLine)) -xC-len/3.9;

   XD = 0:ds:len/sqrt(2);
   YD = zeros(1,length(XD));

   %rotate!
   for i=1:length(XD)
      xDt(i) = XD(i)*cos(pi/4)-YD(i)*sin(pi/4);
      yDt(i) = XD(i)*sin(pi/4)+YD(i)*cos(pi/4);

      xDb(i) = XD(i)*cos(pi/4)-YD(i)*sin(-pi/4);
      yDb(i) = XD(i)*sin(-pi/4)+YD(i)*cos(pi/4);
   end
   xDt = xDt - xC - len/4;
   yDt = yDt - yC;

   xDb = xDb - xC - len/4;
   yDb = yDb - yC;

   x = [xLine(1:ceil(0.75*end)) xDt(1:ceil(2*end/3)) xDb(1:ceil(2*end/3))];
   y = [yLine(1:ceil(0.75*end)) yDt(1:ceil(2*end/3))-len/8 yDb(1:ceil(2*end/3))-len/5.5];

elseif strcmp(letter,'K')

   yLine = -len/2:ds:len/2;
   yLine = yLine - yC;
   xLine = zeros(1,length(yLine)) -xC-len/3.9;

   XD = 0:ds:len/sqrt(2);
   YD = zeros(1,length(XD));

   %rotate!
   for i=1:length(XD)
      xDt(i) = XD(i)*cos(pi/4)-YD(i)*sin(pi/4);
      yDt(i) = XD(i)*sin(pi/4)+YD(i)*cos(pi/4);

      xDb(i) = XD(i)*cos(pi/4)-YD(i)*sin(-pi/4);
      yDb(i) = XD(i)*sin(-pi/4)+YD(i)*cos(pi/4);
   end
   xDt = xDt - xC - len/4.5;
   yDt = yDt - yC;

   xDb = xDb - xC - len/5.5;
   yDb = yDb - yC;

   x = [xLine xDt xDb];
   y = [yLine yDt yDb];

elseif strcmp(letter,'l')

   yLine = -len/2:ds:len/2.25;
   yLine = yLine - yC;
   xLine = 1/2*len/2*ones(1,length(yLine));

   x = -xLine-xC;
   y = yLine;

elseif strcmp(letter,'m')

    r = len/2;
    xH = -r:ds:r;  xH = xH - xC;
    yH = zeros(1,length(xH)) - yC;
    yV = 0:-ds:-len/2; yV = yV - yC;
    xV = -r*ones(1,length(yV)) - xC;
    xV2 = zeros(1,length(yV)) - xC;
    xV3 = r*ones(1,length(yV)) - xC;

    yVS = 0:ds:len/12;
    yVS = yVS - yC;
    xVS = -r*ones(1,length(yVS)) - xC;

    x = [xVS xH xV xV2 xV3];
    y = [yVS yH yV yV yV];

elseif strcmp(letter,'n')

    r = len/4.1;
    yV = [-r:ds:r r];
    xH = yV - xC;
    yH = zeros(1,length(xH)) - yC;

    yV = yV - yC - len/4;
    xV = r*ones(1,length(yV));

    yV2 = 0:ds:len/12;
    yV2 = yV2 - yC;
    xV2 = r*ones(1,length(yV2));

    x = [-xV-(xC-ds/6) -xV2-(xC-ds/6) xH xV-(xC-ds/6)];
    y = [yV yV2 yH yV];

elseif strcmp(letter,'o')

    r = len/4.1;
    dt = ds/r;
    theta = 0:dt:2*pi;
    for i=1:length(theta)
       x(i) = r*cos(theta(i)) - xC;
       y(i) = r*sin(theta(i)) - yC - len/4;
    end

 elseif strcmp(letter,'r')

    r = len/4.1;
    yV = [-r:ds:r r];
    xH = yV - xC;
    yH = zeros(1,length(xH)) - yC;

    yV = yV - yC - len/4;
    xV = r*ones(1,length(yV));

    yV2 = 0:ds:len/12;
    yV2 = yV2 - yC;
    xV2 = r*ones(1,length(yV2));

    yV3 = 0:-ds:-len/10;
    yV3 = yV3 - yC;
    xV3 = r*ones(1,length(yV3));

    x = [-xV-(xC-ds/6) -xV2-(xC-ds/6) xH xV3-(xC-ds/6)];
    y = [yV yV2 yH yV3];

elseif strcmp(letter,'t')

   yLine = -len/2:ds:len/2.25;
   yLine = yLine - yC;
   xLine = zeros(1,length(yLine));

   xH = -len/4:ds:len/4;
   xH = xH-xC;
   yH = zeros(1,length(xH)) - yC + len/9;

   x = [xLine-xC xH];
   y = [yLine yH];

elseif strcmp(letter,'u')

    r = len/4.1;
    yV = [-r:ds:r r];
    xH = yV - xC;
    yH = zeros(1,length(xH)) - yC - len/2;

    yV = yV - yC - len/4;
    xV = r*ones(1,length(yV));

    x = [-xV-(xC-ds/6) xH xV-(xC-ds/6)];
    y = [yV yH yV];

elseif strcmp(letter,'v')

   XD = 0:ds:sqrt(5)*len/4;
   YD = zeros(1,length(XD));

   %rotate!
   ang = atan(2);
   for i=1:length(XD)
      xDr(i) = XD(i)*cos(ang)-YD(i)*sin(ang);
      yDr(i) = XD(i)*sin(ang)+YD(i)*cos(ang);

      xDL(i) = XD(i)*cos(pi-ang)-YD(i)*sin(pi-ang);
      yDL(i) = XD(i)*sin(pi-ang)+YD(i)*cos(pi-ang);
   end
   xDr = xDr - xC;
   yDr = yDr - yC - len/2;

   xDL = xDL - xC;
   yDL = yDL - yC - len/2;

   x = [xDL xDr(2:end)];
   y = [yDL yDr(2:end)];

elseif strcmp(letter,'w')

   XD = 0:ds:sqrt(5)*len/4;
   YD = zeros(1,length(XD));

   %rotate!
   ang = atan(2);
   for i=1:length(XD)
      xDr(i) = XD(i)*cos(ang)-YD(i)*sin(ang);
      yDr(i) = XD(i)*sin(ang)+YD(i)*cos(ang);

      xDL(i) = XD(i)*cos(pi-ang)-YD(i)*sin(pi-ang);
      yDL(i) = XD(i)*sin(pi-ang)+YD(i)*cos(pi-ang);
   end
   xDr = xDr - xC - len/4;
   yDr = yDr - yC - len/2;

   xDL = xDL - xC - len/4;
   yDL = yDL - yC - len/2;

   x = [xDL xDr(2:end) xDL+len/2 xDr(2:end)+len/2];
   y = [yDL yDr(2:end) yDL yDr(2:end)];

elseif strcmp(letter,'W')

   XD = 0:ds:sqrt(5)*len/2;
   YD = zeros(1,length(XD));

   %rotate!
   ang = atan(2);
   for i=1:length(XD)
      xDr(i) = XD(i)*cos(ang)-YD(i)*sin(ang);
      yDr(i) = XD(i)*sin(ang)+YD(i)*cos(ang);

      xDL(i) = XD(i)*cos(pi-ang)-YD(i)*sin(pi-ang);
      yDL(i) = XD(i)*sin(pi-ang)+YD(i)*cos(pi-ang);
   end
   xDr = xDr - xC - len/4;
   yDr = yDr - yC - len/2;

   xDL = xDL - xC - len/4;
   yDL = yDL - yC - len/2;

   x = [xDL xDr(2:floor(end/2)) xDL(1:floor(end/2))+len/2 xDr(2:end)+len/2];
   y = [yDL yDr(2:floor(end/2)) yDL(1:floor(end/2))       yDr(2:end)];

elseif strcmp(letter,'y')

   XD = 0:ds:sqrt(5)*len/4;
   YD = zeros(1,length(XD));

   %rotate!
   ang = atan(2);
   for i=1:length(XD)
      xDr(i) = XD(i)*cos(ang)-YD(i)*sin(ang);
      yDr(i) = XD(i)*sin(ang)+YD(i)*cos(ang);

      xDL(i) = XD(i)*cos(pi-ang)-YD(i)*sin(pi-ang);
      yDL(i) = XD(i)*sin(pi-ang)+YD(i)*cos(pi-ang);

      xB(i) = XD(i)*cos(pi+ang)-YD(i)*sin(pi+ang);
      yB(i) = XD(i)*sin(pi+ang)+YD(i)*cos(pi+ang);
   end
   xDr = xDr - xC;
   yDr = yDr - yC - len/2;

   xDL = xDL - xC;
   yDL = yDL - yC - len/2;

   xB = xB - xC;
   yB = yB - yC - len/2;

   x = [xDL xDr(2:end) xB(2:end)];
   y = [yDL yDr(2:end) yB(2:end)];

elseif strcmp(letter,'!')

   yLine = -len/4:ds:len/2;
   yLine = yLine - yC;
   xLine = zeros(1,length(yLine)) - xC;

    r = len/21;
    dt = ds/(2*r);
    theta = 0:dt:2*pi;
    for i=1:length(theta)
       xDot(i) = r*cos(theta(i)) - xC;
       yDot(i) = r*sin(theta(i)) - yC - len/2.1;
    end

    x = [xLine xDot];
    y = [yLine yDot];

elseif strcmp(letter,'?')

    %xT = 0:ds:len/8-ds; xT = xT - xC - len/4;
    %yT = zeros(1,length(xT)) - yC + len/2;

    r = len/4;
    dt = ds/r;
    theta = 1.25*pi/2:-dt:-pi/2;
    for i=1:length(theta)
       xL(i) = 1.5*r*cos(theta(i)) - xC - len/8;
       yL(i) = r*sin(theta(i)) - yC + len/4;
    end

    yV = -ds:-ds:-len/3.25;
    xV = zeros(1,length(yV)) - xC - len/8;

    r = len/21;
    dt = ds/(2*r);
    theta = 0:dt:2*pi;
    for i=1:length(theta)
       xDot(i) = r*cos(theta(i)) - xC -len/8;
       yDot(i) = r*sin(theta(i)) - yC - len/2.1;
    end

    x = [xL xV xDot];
    y = [yL yV yDot];

end

author: nickabattista


text:
function FF = computeAllStatistics(fileName, win, step)

% This function computes the average and std values for the following audio
% features:
% - energy entropy
% - short time energy
% - spectral rolloff
% - spectral centroid
% - spectral flux
%
% ARGUMENTS:
% fileName: the name of the .wav file in which the signal is stored
% win: the processing window (in seconds)
% step: the processing step (in seconds)
%
% RETURN VALUE:
% F: a 12x1 array containing the 12 feature statistics
%

[x, fs] = wavread(fileName);

EE = Energy_Entropy_Block(x, win*fs, step*fs, 10);
E = ShortTimeEnergy(x, win*fs, step*fs);
Z = zcr(x, win*fs, step*fs, fs);
R = SpectralRollOff(x, win*fs, step*fs, 0.80, fs);
C = SpectralCentroid(x, win*fs, step*fs, fs);
F = SpectralFlux(x, win*fs, step*fs, fs);

FF(1) = statistic(EE, 1, length(EE), 'std');
FF(2) = statistic(Z, 1, length(Z), 'stdbymean');
FF(3) = statistic(R, 1, length(R), 'std');
FF(4) = statistic(C, 1, length(C), 'std');
FF(5) = statistic(F, 1, length(F), 'std');
FF(6) = statistic(E, 1, length(E), 'stdbymean');


text:
function [str, tag, cind, ccnt] = gencode_item(item, tag, tagctx, stoptag, tropts)

% function [str, tag, cind, ccnt] = gencode_item(item, tag, tagctx, stoptag, tropts)
% Generate code to recreate a generic item. This code does not deal with
% arrays of cfg_items, such a configuration should not exist with the
% current definition of a configuration tree.
%
% Traversal options
% struct with fields
% stopspec - match spec to stop forced setting of eflag
% dflag    - (not used here)
% clvl     - current level in tree
% mlvl     - maximum level to force settings - range 1 (top level only) to
%            Inf (all levels)
% cnt      - item count - used for unique tags
% mcnt     - (not evaluated here)
%
% 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: gencode_item.m 5678 2013-10-11 14:58:04Z volkmar $

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

%% Parent object
% Generate generic object
[str, tag, cind, ccnt] = gencode_item(item.cfg_item, tag, tagctx, stoptag, tropts);
% Check whether to generate code - ccnt == 0 means that generic object did
% not return code
if (tropts.clvl > tropts.mlvl || (~isempty(tropts.stopspec) && match(item, tropts.stopspec))) || ccnt == 0
    str = {};
    cind = [];
    ccnt = 0;
    return;
end;
% Reclassify generic object
str{cind} = sprintf('%s         = %s;', tag, class(item));
%% Strtype
% Generate strtype field
if ~isempty(item.strtype)
    str1 = gencode(item.strtype, sprintf('%s.strtype', tag), tagctx);
    str = [str(:)' str1(:)'];
end;
%% Num
% Generate num field
str{end+1} = sprintf('%s.num     = [%s];', tag, num2str(item.num));
%% Extras
% Generate extras field
if ~isempty(item.extras)
    str1 = gencode(item.extras, sprintf('%s.extras', tag), tagctx);
    str = [str(:)' str1(:)'];
end;


text:
function [b1, b2] = isManifoldMesh(varargin)
%ISMANIFOLDMESH Check whether the input mesh may be considered as manifold.
%
%   B = isManifoldMesh(V, F)
%   B = isManifoldMesh(V, E, F)
%   Checks if the specified mesh is a manifold. When mesh is a manifold,
%   all edges are connected to either 2 or 1 faces.
%
%   [B, HASBORDER] = isManifoldMesh(V, E, F)
%   Also checks whether the mesh contains border faces. Border faces
%   contains at least one edge which is ajacent to only one face.
%
%   Example
%     [V, F] = createOctahedron;
%     isManifoldMesh(V, F)
%     ans =
%       logical
%        1
%
%   See also
%     meshes3d, ensureManifoldMesh, trimMesh

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2019-01-31, using Matlab 9.5.0.944444 (R2018b)
% Copyright 2019-2022 INRA - Cepia Software Platform

vertices = varargin{1};
faces = varargin{2};

% compute edge to vertex array
if nargin == 3
    edges = faces;
    faces = varargin{3};
else
    % compute edge to vertex array
    edges = meshEdges(faces);
end

% compute face to edge indices array
% as a nFaces-by-3 array (each face connected to exactly three edges)
faceEdgeInds = meshFaceEdges(vertices, edges, faces);

% compute number of faces incident each edge
edgeFaces = trimeshEdgeFaces(faces);
edgeFaceDegrees = sum(edgeFaces > 0, 2);

% for each face, concatenate the face degree of each edge
faceEdgeDegrees = zeros(size(faces, 1), 3);
for iFace = 1:size(faces, 1)
    edgeInds = faceEdgeInds{iFace};
    faceEdgeDegrees(iFace, :) = edgeFaceDegrees(edgeInds);
end

regFaces = sum(ismember(faceEdgeDegrees, [1 2]), 2) == 3;
innerFaces = sum(faceEdgeDegrees == 2, 2) == 3;
borderFaces = regFaces & ~innerFaces;

% check if mesh is manifold: all faces are either regular or border
b1 = all(regFaces);

% check if some faces are border
b2 = any(borderFaces);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function net = net_move_to_device(net, device)

isDag = isa(net, 'dagnn.DagNN');

if isDag
   net.move(device)
else
   net = vl_simplenn_move(net, device) ;
   if strcmp(device, 'gpu')
       net.useGpu = true;
   else
       net.useGpu = false;
   end
end

repo: DBCNN


text:
function hydro = excitationIRF(hydro,tEnd,nDt,nDw,wMin,wMax)
% Calculates the normalized excitation impulse response function:
%
% 	:math:`\overline{K}_{e,i,\theta}(t) = {\frac{1}{2\pi}}\intop_{-\infty}^{\infty}{\frac{X_i(\omega,\theta)e^{i{\omega}t}}{{\rho}g}}d\omega`
%
% Default parameters can be used by inputting [].
% See ``WEC-Sim\examples\BEMIO`` for examples of usage.
%
% Parameters
% ----------
%     hydro : struct
%         Structure of hydro data
%
%     tEnd : float
%         Calculation range for the IRF, where the IRF is calculated from t
%         = 0 to tEnd, and the default is 100 s
%
%     nDt : float
%         Number of time steps in the IRF, the default is 1001
%
%     nDw : float
%         Number of frequency steps used in the IRF calculation
%         (hydrodynamic coefficients are interpolated to correspond), the
%         default is 1001
%
%     wMin : float
%         Minimum frequency to use in the IRF calculation, the default is
%         the minimum frequency from the BEM data
%
%     wMax : float
%         Maximum frequency to use in the IRF calculation, the default is
%         the maximum frequency from the BEM data
%
% Returns
% -------
%     hydro : struct
%         Structure of hydro data with excitation IRF
%

p = waitbar(0,'Calculating excitation IRFs...');  % Progress bar

% Set defaults if empty
if isempty(tEnd)==1;  tEnd = 100;           end
if isempty(nDt)==1;    nDt = 1001;            end
if isempty(nDw)==1;    nDw = 1001;            end
if isempty(wMin)==1;  wMin = min(hydro.w);  end
if isempty(wMax)==1;  wMax = max(hydro.w);  end

% Interpolate to the given t and w
t = linspace(-tEnd,tEnd,nDt);
w = linspace(wMin,wMax,nDw);
N = sum(hydro.dof)*hydro.Nh;

% Calculate the impulse response function for excitation
n = 0;
for i = 1:sum(hydro.dof)
    for j = 1:hydro.Nh
        ex_re = interp1(hydro.w,squeeze(hydro.ex_re(i,j,:)),w);
        ex_im = interp1(hydro.w,squeeze(hydro.ex_im(i,j,:)),w);
        hydro.ex_K(i,j,:) = (1/pi)*trapz(w,ex_re.*cos(w.*t(:))-ex_im.*sin(w.*t(:)),2);
        n = n+1;
    end
    waitbar(n/N)
end

hydro.ex_t = t;
hydro.ex_w = w;
close(p);

end

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
function file_name = file_name_inc ( file_name )

%*****************************************************************************80
%
% FILE_NAME_INC generates the next filename in a series.
%
%  Discussion:
%
%    It is assumed that the digits in the name, whether scattered or
%    connected, represent a number that is to be increased by 1 on
%    each call.  If this number is all 9's on input, the output number
%    is all 0's.  Non-numeric letters of the name are unaffected..
%
%    If the name is empty, then the routine stops.
%
%    If the name contains no digits, the empty string is returned.
%
%  Example:
%
%      Input            Output
%      -----            ------
%      'a7to11.txt'     'a7to12.txt'  (typical case.  Last digit incremented)
%      'a7to99.txt'     'a8to00.txt'  (last digit incremented, with carry.)
%      'a9to99.txt'     'a0to00.txt'  (wrap around)
%      'cat.txt'        ' '           (no digits in input name.)
%      ' '              STOP!         (error.)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 September 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character FILE_NAME(*), the string to be incremented.
%
%    Output, character FILE_NAME_NEW(*), the incremented string.
%
  lens = length ( file_name );

  if ( lens <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_NAME_INC - Fatal error!\n' );
    fprintf ( 1, '  The input filename is empty.\n' );
    error ( 'FILE_NAME_INC - Fatal error!' );
  end

  change = 0;

  for i = lens : -1 : 1

    c = file_name(i);

    if ( '0' <= c & c <= '8' )

      change = change + 1;

      c = c + 1;

      file_name(i) = c;

      return

    elseif ( c == '9' )

      change = change + 1;

      c = '0';

      file_name(i) = c;

    end

  end

  if ( change == 0 )
    file_name = ' ';
  end

  return
end


text:
function dist = dist_euclidean(metric, x)
% DIST_EUCLIDEAN   is a mex file that returns the euclidean distance matrix
% for metric_euclidean.  If we end up here the mex file has not been
% compiled. In this case the function returns NaN and the distance is
% evaluated in metric_euclidean.

% Copyright (c) 2008 Jarno Vanhatalo

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

dist = NaN;


text:
function arrivalObj = antelope(dbname, subset_expr)
% LOAD_ARRIVALS Load arrivals from a CSS3.0 database
% Modifed by Glenn Thompson 2017-11-28 to include some fields from assoc table
% if present, necessary to support hankelq code
ARRIVAL_TABLE_PRESENT = antelope.dbtable_present(dbname, 'arrival');
ASSOC_TABLE_PRESENT = antelope.dbtable_present(dbname, 'assoc');
ORIGIN_TABLE_PRESENT = antelope.dbtable_present(dbname, 'origin');
    if ~ARRIVAL_TABLE_PRESENT
        fprintf('No arrival table belonging to %s\n',dbname);
        return
    end

    fprintf('Loading arrivals from %s\n',dbname);

    % Open database
    db = dbopen(dbname,'r');
    disp('- database opened');

    % Apply subset expression
    db = dblookup_table(db,'arrival');
    disp('- arrival table opened');
    if exist('subset_expr','var')
        db = dbsubset(db,subset_expr);
        disp('- subsetted database')
    end

    nrows = dbnrecs(db);
    if nrows > 0

        % Sort by arrival time
        db = dbsort(db,'time');
        disp('- sorted arrival table')

        seaz =[]; deltim=[]; delta=[]; orid=[]; evid=[]; timeres=[]; otime=[]; depth=[];
        fprintf('- reading %d rows\n',nrows);
        [sta,chan,time,amp,signal2noise,iphase, arid] = dbgetv(db,'sta','chan','time','amp','snr','iphase','arid');

        % Join with assoc table if present
        if ASSOC_TABLE_PRESENT
            db2 = dblookup_table( db , 'assoc');
            db = dbjoin(db, db2);

            % Get the values
            [sta,chan,time,amp,signal2noise,iphase, arid, seaz, deltim, delta, orid, timeres] = ...
                dbgetv(db,'sta','chan','time','amp','snr','iphase','arid', ...
                'seaz', 'deltim', 'delta', 'orid', 'timeres');

            % Join with origin table if present
            if ORIGIN_TABLE_PRESENT
                db3 = dblookup_table( db , 'origin');
                db = dbjoin(db, db3);
                [evid,otime,depth] = dbgetv(db,'evid','time','depth');
            end
        end

        % Create arrival object
        disp('- creating arrival object')
        arrivalObj = Arrival(cellstr(sta), cellstr(chan), epoch2datenum(time), cellstr(iphase), ...
            'amp', amp, 'signal2noise', signal2noise, 'arid', arid, 'seaz', seaz, ...
            'deltim', deltim, 'delta', delta, 'otime', otime, 'orid', orid, ...
            'evid', evid, 'timeres', timeres, 'depth', depth);

        % Close database link
        dbclose(db);
        disp('- database closed')

        disp('- complete!')
    else
        fprintf('No arrivals found matching request\n');
    end

end

author: geoscience-community-codes

repo: GISMO


text:
function [ x, dp2, p1 ] = jacobi_ss_root ( x, n, alpha, beta, b, c )

%*****************************************************************************80
%
%% JACOBI_SS_ROOT improves an approximate root of a Jacobi polynomial.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 October 2005
%
%  Author:
%
%    Original FORTRAN77 version by Arthur Stroud, Don Secrest.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Arthur Stroud, Don Secrest,
%    Gaussian Quadrature Formulas,
%    Prentice Hall, 1966,
%    LC: QA299.4G3S7.
%
%  Parameters:
%
%    Input, real X, the approximate root.
%
%    Input, integer N, the order of the polynomial to be computed.
%
%    Input, real ALPHA, BETA, the exponents of (1-X) and
%    (1+X) in the quadrature rule.
%
%    Input, real B(N), C(N), the recursion coefficients.
%
%    Output, real X, the improved approximate root.
%
%    Output, real DP2, the value of J'(N)(X).
%
%    Output, real P1, the value of J(N-1)(X).
%
  maxstep = 10;

  for i = 1 : maxstep

    [ p2, dp2, p1 ] = jacobi_ss_recur ( x, n, alpha, beta, b, c );

    d = p2 / dp2;
    x = x - d;

    if ( abs ( d ) <= eps * ( abs ( x ) + 1.0 ) )
      return
    end

  end

  return
end


text:
classdef CEC2020_F9 < PROBLEM
% <single> <real>
% Composition function 2

%------------------------------- Reference --------------------------------
% C .T. Yue, K. V. Price, P. N. Suganthan, J. J. Liang, M. Z. Ali, B. Y.
% Qu, N. H. Awad, and P. P Biswas, Problem definitions and evaluation
% criteria for the CEC 2020 special session and competition on single
% objective bound constrained numerical optimization, Zhengzhou University,
% China and Nanyang Technological University, Singapore, 2019.
%------------------------------- 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".
%--------------------------------------------------------------------------

    properties
        O;      % Optimal decision vector
        Mat;	% Rotation matrix
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'CEC2020.mat'),'Data');
            obj.O = Data{9}.o;
            obj.M = 1;
            if isempty(obj.D) || obj.D < 10
                obj.D   = 5;
                obj.Mat = Data{9}.M_5;
            elseif obj.D < 15
                obj.D   = 10;
                obj.Mat = Data{9}.M_10;
            elseif obj.D < 20
                obj.D   = 15;
                obj.Mat = Data{9}.M_15;
            else
                obj.D   = 20;
                obj.Mat = Data{9}.M_20;
            end
            obj.lower    = zeros(1,obj.D) - 100;
            obj.upper    = zeros(1,obj.D) + 100;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            lambda = [10 1e-6 10 1];
            delta  = [10 20 30 40];
            bias   = [0 100 200 300];
            func   = {@Ackley,@Elliptic,@Griewank,@Rastrigin};
            W      = zeros(size(PopDec,1),4);
            F      = zeros(size(W));
            for i = 1 : size(W,2)
                tmp    = sum((PopDec-repmat(obj.O(i,1:size(PopDec,2)),size(PopDec,1),1)).^2,2);
                W(:,i) = 1./(sqrt(tmp)+1e-10).*exp(-tmp/2/obj.D/delta(i)^2);
                F(:,i) = func{i}((PopDec-repmat(obj.O(i,1:size(PopDec,2)),size(PopDec,1),1))*obj.Mat((i-1)*obj.D+1:i*obj.D,:)');
            end
            W = W./repmat(sum(W,2),1,size(W,2));
            PopObj = 2400 + sum(W.*(repmat(lambda,size(F,1),1).*F+repmat(bias,size(F,1),1)),2);
        end
    end
end

function F = Ackley(X)
    F = -20*exp(-0.2*sqrt(mean(X.^2,2))) - exp(mean(cos(2*pi*X),2)) + 20 + exp(1);
end

function F = Elliptic(X)
    F = sum((1e6).^(repmat(0:size(X,2)-1,size(X,1),1)/(size(X,2)-1+1e-6)).*X.^2,2);
end

function F = Griewank(X)
    X = 6*X;
    F = sum(X.^2,2)/4000 - prod(cos(X)./repmat(sqrt(1:size(X,2)),size(X,1),1),2) + 1;
end

function F = Rastrigin(X)
    X = 0.0512*X;
    F = sum(X.^2-10*cos(2*pi*X)+10,2);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO