text:
%% Microstrip low-pass filter analysis using 3D FDTD code with UPML
%% absorbing borders (ABC)
%
% Here we use FDTD 3D with UPML to calculate scattering coefficients S_{11}
% and S_{21} for planar microstrip low-pass filter following by the original
% paper by D. Sheen, S. Ali, M. Abouzahra, J. Kong "Application of the
% Three-Dimensional Finite-Difference Time-Domain Method to the Analysis of
% planar Microstrip Circuits", IEEE Trans. on Microwave Theory and Techniques
% (http://dx.doi.org/10.1109/22.55775).
% Also |S_{21}| dependence can be taken from the paper "Computational
% electromagnetic method for interconnects and small structures" by C.
% Balanis, A. Policarpou and S. Georgakopoulos
% (http://dx.doi.org/10.1006/spmi.2000.0865)
%% Current code includes some improvements in comparison with the original
%% calculations:
% 1) imposing UPML instead of Mur ABCs;
% 2) using real metal (copper) as a patch conductor material instead of PEC;
% 3) applying matched load at the ends of filter's transmission microstrip
%    lines to prevent physical reflections;
% 4) no "magnetic wall" or "electric wall" conditions at the Ez source plane.
%
function FDTD_3D_Lowpass
close all; clear; clc;
%% Physical constants
   epsilon0 = 8.85418782e-12; mu0 = 1.25663706e-6;
   c = 1.0/sqrt(mu0*epsilon0);

%% Gaussian half-width
   t_half = 15.0e-12;

%% Microstrip transmission lines parameters
   lineW = 2.413e-3;
   lineH = 1.0e-3;
   lineEr = 2.2;
   Z0 = 49.2526;

%% End time
   t_end = 1.5e-9;

%% Total mesh dimensions and grid cells sizes (without PML)
   nx = 80; ny = 100; nz = 16;
   dx = 0.4064e-3; dy = 0.4233e-3; dz = 0.2650e-3;

%% Number of PML layers
   PML = 5;

%% Matrix of material's constants
   number_of_materials = 4;
   % For material of number x = 1,2,3... :
   % Material(x,1) - relative permittivity, Material(x,2) - relative permeability,
   % Material(x,3) - specific conductivity
   % Vacuum
   Material(1,1) = 1.0;   Material(1,2) = 1.0;   Material(1,3) = 0.0;
   % Metal (Copper)
   Material(2,1) = 1.0;   Material(2,2) = 1.0;   Material(2,3) = 5.88e+7;
   % Substrate material (RT/Duroid 5880)
   Material(3,1) = lineEr;   Material(3,2) = 1.0;   Material(3,3) = 0.0;
   % Matched load material is calculated from transmission line parameters
   Material(4,1) = 1.0;   Material(4,2) = 1.0;   Material(4,3) = lineH/(Z0*lineW*dy);

% Add PML layers
   nx = nx + 2*PML; ny = ny + 2*PML; nz = nz + 2*PML;
% Calculate dt
   dt = (1.0/c/sqrt( 1.0/(dx^2) + 1.0/(dy^2) + 1.0/(dz^2)))*0.9999;
   number_of_iterations = ceil(t_end/dt);

%% 3D array for geometry
   Index = ones(nx, ny, nz);

%% Define of low-pass filter geometry
   % Ground plane
   Index((1+PML):(nx-PML), (1+PML):(ny-PML), PML+1) = 2;
   % Rectangular patch (one cell thickness)
   Index((nx/2-25):(nx/2+25), (ny/2-3):(ny/2+3), PML+5) = 2;
   % Transmission line from port 1 to rectangular patch
   Index((nx/2-10):(nx/2-5), (PML+1):ny/2, PML+5) = 2;
   % Transmission line from rectangular patch to port 2
   Index((nx/2+5):(nx/2+10), ny/2:(ny-PML), PML+5) = 2;
   % Dielectric substrate between ground plane and filter patch
   Index((1+PML):(nx-PML), (1+PML):(ny-PML), (PML+2):(PML+4)) = 3;
   % Matched load before port 1
   Index((nx/2-10):(nx/2-5), PML+1, (PML+2):(PML+4)) = 4;
   % Matched load after port 2
   Index((nx/2+5):(nx/2+10), ny-PML, (PML+2):(PML+4)) = 4;

%% 3D FDTD physical (fields) and additional arrays are defined as 'single'
%% to increase performance
   Ex = zeros(nx, ny+1, nz+1, 'single');
   Gx = zeros(nx, ny+1, nz+1, 'single');
   Fx = zeros(nx, ny+1, nz+1, 'single');
   Ey = zeros(nx+1, ny, nz+1, 'single');
   Gy = zeros(nx+1, ny, nz+1, 'single');
   Fy = zeros(nx+1, ny, nz+1, 'single');
   Ez = zeros(nx+1, ny+1, nz, 'single');
   Gz = zeros(nx+1, ny+1, nz, 'single');
   Fz = zeros(nx+1, ny+1, nz, 'single');
   Hx = zeros(nx+1, ny, nz, 'single');
   Bx = zeros(nx+1, ny, nz, 'single');
   Hy = zeros(nx, ny+1, nz, 'single');
   By = zeros(nx, ny+1, nz, 'single');
   Hz = zeros(nx, ny, nz+1, 'single');
   Bz = zeros(nx, ny, nz+1, 'single');

%% FDTD PML coefficients arrays. Here they are already filled with values
%% corresponding to free space
   m = 4; ka_max = 1.0; R_err = 1.0e-16;
   eta = sqrt(mu0/epsilon0*Material(1,1)/Material(1,2));
   k_Ex_c = ones(nx, ny, nz, 'single')*2.0*epsilon0;
   k_Ex_d = ones(nx, ny, nz, 'single')*(-2.0*epsilon0);
   k_Ey_a = ones(nx+1, ny, nz, 'single');
   k_Ey_b = ones(nx+1, ny, nz, 'single')/(2.0*epsilon0);
   k_Gz_a = ones(nx+1, ny, nz, 'single');
   k_Gz_b = ones(nx+1, ny, nz, 'single');
   k_Hy_a = ones(nx, ny, nz, 'single');
   k_Hy_b = ones(nx, ny, nz, 'single')/(2.0*epsilon0);
   k_Hx_c = ones(nx+1, ny, nz, 'single')*2.0*epsilon0/mu0;
   k_Hx_d = ones(nx+1, ny, nz, 'single')*(-2.0*epsilon0/mu0);
   k_Bz_a = ones(nx, ny, nz, 'single');
   k_Bz_b = ones(nx, ny, nz, 'single')*dt;
   k_Gx_a = ones(nx, ny+1, nz, 'single');
   k_Gx_b = ones(nx, ny+1, nz, 'single');
   k_Ey_c = ones(nx, ny, nz, 'single')*2.0*epsilon0;
   k_Ey_d = ones(nx, ny, nz, 'single')*(-2.0*epsilon0);
   k_Ez_a = ones(nx, ny+1, nz, 'single');
   k_Ez_b = ones(nx, ny+1, nz, 'single')/(2.0*epsilon0);
   k_Bx_a = ones(nx, ny, nz, 'single');
   k_Bx_b = ones(nx, ny, nz, 'single')*dt;
   k_Hy_c = ones(nx, ny+1, nz, 'single')*2.0*epsilon0/mu0;
   k_Hy_d = ones(nx, ny+1, nz, 'single')*(-2.0*epsilon0/mu0);
   k_Hz_a = ones(nx, ny, nz, 'single');
   k_Hz_b = ones(nx, ny, nz, 'single')/(2.0*epsilon0);
   k_Ex_a = ones(nx, ny, nz+1, 'single');
   k_Ex_b = ones(nx, ny, nz+1, 'single')/(2.0*epsilon0);
   k_Gy_a = ones(nx, ny, nz+1, 'single');
   k_Gy_b = ones(nx, ny, nz+1, 'single');
   k_Ez_c = ones(nx, ny, nz, 'single')*2.0*epsilon0;
   k_Ez_d = ones(nx, ny, nz, 'single')*(-2.0*epsilon0);
   k_Hx_a = ones(nx, ny, nz, 'single');
   k_Hx_b = ones(nx, ny, nz, 'single')/(2.0*epsilon0);
   k_By_a = ones(nx, ny, nz, 'single');
   k_By_b = ones(nx, ny, nz, 'single')*dt;
   k_Hz_c = ones(nx, ny, nz+1, 'single')*2.0*epsilon0/mu0;
   k_Hz_d = ones(nx, ny, nz+1, 'single')*(-2.0*epsilon0/mu0);

%% General FDTD coefficients
   I = 1:number_of_materials;
   K_a(I) = (2.0*epsilon0*Material(I,1) - Material(I,3)*dt)./...
            (2.0*epsilon0*Material(I,1) + Material(I,3)*dt);
   K_b(I) = 2.0*dt./(2.0*epsilon0*Material(I,1) + Material(I,3)*dt);
   K_c(I) = Material(I,2);
   Ka = single(K_a(Index)); Kb = single(K_b(Index)); Kc = single(K_c(Index));

%% PML coefficients along x-axis
   sigma_max = -(m + 1.0)*log(R_err)/(2.0*eta*PML*dx);
   for I=0:(PML-1)
        sigma_x = sigma_max*((PML - I)/PML)^m;
        ka_x = 1.0 + (ka_max - 1.0)*((PML - I)/PML)^m;
        k_Ey_a(I+1,:,:) = (2.0*epsilon0*ka_x - sigma_x*dt)/...
                          (2.0*epsilon0*ka_x + sigma_x*dt);
        k_Ey_a(nx-I,:,:) = k_Ey_a(I+1,:,:);
        k_Ey_b(I+1,:,:) = 1.0/(2.0*epsilon0*ka_x + sigma_x*dt);
        k_Ey_b(nx-I,:,:) = k_Ey_b(I+1,:,:);
        k_Gz_a(I+1,:,:) = (2.0*epsilon0*ka_x - sigma_x*dt)/...
                          (2.0*epsilon0*ka_x + sigma_x*dt);
        k_Gz_a(nx-I,:,:) = k_Gz_a(I+1,:,:);
        k_Gz_b(I+1,:,:) = 2.0*epsilon0/(2.0*epsilon0*ka_x + sigma_x*dt);
        k_Gz_b(nx-I,:,:) = k_Gz_b(I+1,:,:);
        k_Hx_c(I+1,:,:) = (2.0*epsilon0*ka_x + sigma_x*dt)/mu0;
        k_Hx_c(nx-I,:,:) = k_Hx_c(I+1,:,:);
        k_Hx_d(I+1,:,:) = -(2.0*epsilon0*ka_x - sigma_x*dt)/mu0;
        k_Hx_d(nx-I,:,:) = k_Hx_d(I+1,:,:);

        sigma_x = sigma_max*((PML - I - 0.5)/PML)^m;
        ka_x = 1.0 + (ka_max - 1.0)*((PML - I - 0.5)/PML)^m;
        k_Ex_c(I+1,:,:) = 2.0*epsilon0*ka_x + sigma_x*dt;
        k_Ex_c(nx-I-1,:,:) = k_Ex_c(I+1,:,:);
        k_Ex_d(I+1,:,:) = -(2.0*epsilon0*ka_x - sigma_x*dt);
        k_Ex_d(nx-I-1,:,:) = k_Ex_d(I+1,:,:);
        k_Hy_a(I+1,:,:) = (2.0*epsilon0*ka_x - sigma_x*dt)/...
                          (2.0*epsilon0*ka_x + sigma_x*dt);
        k_Hy_a(nx-I-1,:,:) = k_Hy_a(I+1,:,:);
        k_Hy_b(I+1,:,:) = 1.0/(2.0*epsilon0*ka_x + sigma_x*dt);
        k_Hy_b(nx-I-1,:,:) = k_Hy_b(I+1,:,:);
        k_Bz_a(I+1,:,:) = (2.0*epsilon0*ka_x - sigma_x*dt)/...
                          (2.0*epsilon0*ka_x + sigma_x*dt);
        k_Bz_a(nx-I-1,:,:) = k_Bz_a(I+1,:,:);
        k_Bz_b(I+1,:,:) = 2.0*epsilon0*dt/(2.0*epsilon0*ka_x + sigma_x*dt);
        k_Bz_b(nx-I-1,:,:) = k_Bz_b(I+1,:,:);
   end

%% PML coefficients along y-axis
   sigma_max = -(m + 1.0)*log(R_err)/(2.0*eta*PML*dy);
   for J=0:(PML-1)
        sigma_y = sigma_max*((PML - J)/PML)^m;
        ka_y = 1.0 + (ka_max - 1.0)*((PML - J)/PML)^m;
        k_Gx_a(:,J+1,:) = (2.0*epsilon0*ka_y - sigma_y*dt)/...
                          (2.0*epsilon0*ka_y + sigma_y*dt);
        k_Gx_a(:,ny-J,:) = k_Gx_a(:,J+1,:);
        k_Gx_b(:,J+1,:) = 2.0*epsilon0/(2.0*epsilon0*ka_y + sigma_y*dt);
        k_Gx_b(:,ny-J,:) = k_Gx_b(:,J+1,:);
        k_Ez_a(:,J+1,:) = (2.0*epsilon0*ka_y - sigma_y*dt)/...
                          (2.0*epsilon0*ka_y + sigma_y*dt);
        k_Ez_a(:,ny-J,:) = k_Ez_a(:,J+1,:);
        k_Ez_b(:,J+1,:) = 1.0/(2.0*epsilon0*ka_y + sigma_y*dt);
        k_Ez_b(:,ny-J,:) = k_Ez_b(:,J+1,:);
        k_Hy_c(:,J+1,:) = (2.0*epsilon0*ka_y + sigma_y*dt)/mu0;
        k_Hy_c(:,ny-J,:) = k_Hy_c(:,J+1,:);
        k_Hy_d(:,J+1,:) = -(2.0*epsilon0*ka_y - sigma_y*dt)/mu0;
        k_Hy_d(:,ny-J,:) = k_Hy_d(:,J+1,:);

        sigma_y = sigma_max*((PML - J - 0.5)/PML)^m;
        ka_y = 1.0 + (ka_max - 1.0)*((PML - J - 0.5)/PML)^m;
        k_Ey_c(:,J+1,:) = 2.0*epsilon0*ka_y+sigma_y*dt;
        k_Ey_c(:,ny-J-1,:) = k_Ey_c(:,J+1,:);
        k_Ey_d(:,J+1,:) = -(2.0*epsilon0*ka_y-sigma_y*dt);
        k_Ey_d(:,ny-J-1,:) = k_Ey_d(:,J+1,:);
        k_Bx_a(:,J+1,:) = (2.0*epsilon0*ka_y-sigma_y*dt)/...
                          (2.0*epsilon0*ka_y+sigma_y*dt);
        k_Bx_a(:,ny-J-1,:) = k_Bx_a(:,J+1,:);
        k_Bx_b(:,J+1,:) = 2.0*epsilon0*dt/(2.0*epsilon0*ka_y+sigma_y*dt);
        k_Bx_b(:,ny-J-1,:) = k_Bx_b(:,J+1,:);
        k_Hz_a(:,J+1,:) = (2.0*epsilon0*ka_y-sigma_y*dt)/...
                          (2.0*epsilon0*ka_y+sigma_y*dt);
        k_Hz_a(:,ny-J-1,:) = k_Hz_a(:,J+1,:);
        k_Hz_b(:,J+1,:) = 1.0/(2.0*epsilon0*ka_y+sigma_y*dt);
        k_Hz_b(:,ny-J-1,:) = k_Hz_b(:,J+1,:);
   end

%% PML coefficients along z-axis
   sigma_max = -(m + 1.0)*log(R_err)/(2.0*eta*PML*dz);
   for K=0:(PML-1)
        sigma_z = sigma_max*((PML - K)/PML)^m;
        ka_z = 1.0 + (ka_max - 1.0)*((PML-K)/PML)^m;
        k_Ex_a(:,:,K+1) = (2.0*epsilon0*ka_z - sigma_z*dt)/...
                          (2.0*epsilon0*ka_z + sigma_z*dt);
        k_Ex_a(:,:,nz-K) = k_Ex_a(:,:,K+1);
        k_Ex_b(:,:,K+1) = 1.0/(2.0*epsilon0*ka_z + sigma_z*dt);
        k_Ex_b(:,:,nz-K) = k_Ex_b(:,:,K+1);
        k_Gy_a(:,:,K+1) = (2.0*epsilon0*ka_z - sigma_z*dt)/...
                          (2.0*epsilon0*ka_z + sigma_z*dt);
        k_Gy_a(:,:,nz-K) = k_Gy_a(:,:,K+1);
        k_Gy_b(:,:,K+1) = 2.0*epsilon0/(2.0*epsilon0*ka_z + sigma_z*dt);
        k_Gy_b(:,:,nz-K) = k_Gy_b(:,:,K+1);
        k_Hz_c(:,:,K+1) = (2.0*epsilon0*ka_z + sigma_z*dt)/mu0;
        k_Hz_c(:,:,nz-K) = k_Hz_c(:,:,K+1);
        k_Hz_d(:,:,K+1) = -(2.0*epsilon0*ka_z - sigma_z*dt)/mu0;
        k_Hz_d(:,:,nz-K) = k_Hz_d(:,:,K+1);

        sigma_z = sigma_max*((PML - K - 0.5)/PML)^m;
        ka_z = 1.0 + (ka_max - 1.0)*((PML - K - 0.5)/PML)^m;
        k_Ez_c(:,:,K+1) = 2.0*epsilon0*ka_z + sigma_z*dt;
        k_Ez_c(:,:,nz-K-1) = k_Ez_c(:,:,K+1);
        k_Ez_d(:,:,K+1) = -(2.0*epsilon0*ka_z - sigma_z*dt);
        k_Ez_d(:,:,nz-K-1) = k_Ez_d(:,:,K+1);
        k_Hx_a(:,:,K+1) = (2.0*epsilon0*ka_z - sigma_z*dt)/...
                          (2.0*epsilon0*ka_z + sigma_z*dt);
        k_Hx_a(:,:,nz-K-1) = k_Hx_a(:,:,K+1);
        k_Hx_b(:,:,K+1) = 1.0/(2.0*epsilon0*ka_z + sigma_z*dt);
        k_Hx_b(:,:,nz-K-1) = k_Hx_b(:,:,K+1);
        k_By_a(:,:,K+1) = (2.0*epsilon0*ka_z - sigma_z*dt)/...
                          (2.0*epsilon0*ka_z + sigma_z*dt);
        k_By_a(:,:,nz-K-1) = k_By_a(:,:,K+1);
        k_By_b(:,:,K+1) = 2.0*epsilon0*dt/(2.0*epsilon0*ka_z + sigma_z*dt);
        k_By_b(:,:,nz-K-1) = k_By_b(:,:,K+1);
   end

%% Main 3D FDTD+UPML routine (operates with 'singles' to increase speed)
   hhh = waitbar(0, 'Calculations in progress...');
   tic;
   for T=0:(number_of_iterations-1)
        %% Calculate Fx -> Gx -> Ex
        I = 1:nx; J = 2:ny; K = 2:nz;
        Fx_r = Fx(I,J,K);
        Fx(I,J,K) = Ka(I,J,K).*Fx(I,J,K) + Kb(I,J,K).*...
                    ((Hz(I,J,K) - Hz(I,J-1,K))/dy - (Hy(I,J,K) - Hy(I,J,K-1))/dz);
        Gx_r = Gx(I,J,K);
        Gx(I,J,K) = k_Gx_a(I,J,K).*Gx(I,J,K) + k_Gx_b(I,J,K).*(Fx(I,J,K) - Fx_r);
        Ex(I,J,K) = k_Ex_a(I,J,K).*Ex(I,J,K) + k_Ex_b(I,J,K).*...
                    (k_Ex_c(I,J,K).*Gx(I,J,K) + k_Ex_d(I,J,K).*Gx_r);

        %% Calculate Fy -> Gy -> Ey
        I = 2:nx; J = 1:ny; K = 2:nz;
        Fy_r = Fy(I,J,K);
        Fy(I,J,K) = Ka(I,J,K).*Fy(I,J,K) + Kb(I,J,K).*...
                    ((Hx(I,J,K) - Hx(I,J,K-1))/dz - (Hz(I,J,K) - Hz(I-1,J,K))/dx);
        Gy_r = Gy(I,J,K);
        Gy(I,J,K) = k_Gy_a(I,J,K).*Gy(I,J,K) + k_Gy_b(I,J,K).*(Fy(I,J,K) - Fy_r);
        Ey(I,J,K) = k_Ey_a(I,J,K).*Ey(I,J,K) + k_Ey_b(I,J,K).*...
                    (k_Ey_c(I,J,K).*Gy(I,J,K) + k_Ey_d(I,J,K).*Gy_r);

        %% Calculate Fz -> Gz -> Ez
        I = 2:nx; J = 2:ny; K = 1:nz;
        Fz_r = Fz(I,J,K);
        Fz(I,J,K) = Ka(I,J,K).*Fz(I,J,K) + Kb(I,J,K).*...
                    ((Hy(I,J,K) - Hy(I-1,J,K))/dx - (Hx(I,J,K) - Hx(I,J-1,K))/dy);
        Gz_r = Gz(I,J,K);
        Gz(I,J,K) = k_Gz_a(I,J,K).*Gz(I,J,K) + k_Gz_b(I,J,K).*(Fz(I,J,K) - Fz_r);
        Ez(I,J,K) = k_Ez_a(I,J,K).*Ez(I,J,K) + k_Ez_b(I,J,K).*...
                    (k_Ez_c(I,J,K).*Gz(I,J,K) + k_Ez_d(I,J,K).*Gz_r);

        %% Source of vertical electric field Ez applied to the whole face
        %% plane at the port 1
        if (T*dt<=10.0*t_half)
            Ez((nx/2-10):(nx/2-5), PML+2, (PML+2):(PML+4)) = Source(t_half, 1.0, T*dt);
        end

        %% Save reflected Ez at port 1, passed through Ez at port 2 and Ez
        %% distribution between electrodes
        Incident(T+1) = Source(t_half, 1.0, T*dt);
        Reflected(T+1) = Ez(nx/2-7, PML+2, PML+3) - Incident(T+1);
        Passed(T+1) = Ez(nx/2+7, ny-PML-1, PML+3);
        Ez_saved(1:(ny-2*PML-1), 1:(nx-2*PML-1), T+1) = ...
                Ez((PML+1):(nx-PML-1), (PML+1):(ny-PML-1), PML+3)';

        %% Calculate Bx -> Hx
        I = 2:nx; J = 1:ny; K = 1:nz;
        Bx_r = Bx(I,J,K);
        Bx(I,J,K) = k_Bx_a(I,J,K).*Bx(I,J,K) + k_Bx_b(I,J,K).*...
                    ((Ey(I,J,K+1) - Ey(I,J,K))/dz - (Ez(I,J+1,K) - Ez(I,J,K))/dy);
        Hx(I,J,K) = k_Hx_a(I,J,K).*Hx(I,J,K) + k_Hx_b(I,J,K).*...
                    (k_Hx_c(I,J,K).*Bx(I,J,K) + k_Hx_d(I,J,K).*Bx_r)./Kc(I,J,K);

        %% Calculate By -> Hy
        I = 1:nx; J = 2:ny; K = 1:nz;
        By_r = By(I,J,K);
        By(I,J,K) = k_By_a(I,J,K).*By(I,J,K) + k_By_b(I,J,K).*...
                    ((Ez(I+1,J,K) - Ez(I,J,K))/dx - (Ex(I,J,K+1) - Ex(I,J,K))/dz);
        Hy(I,J,K) = k_Hy_a(I,J,K).*Hy(I,J,K) + k_Hy_b(I,J,K).*...
                    (k_Hy_c(I,J,K).*By(I,J,K) + k_Hy_d(I,J,K).*By_r)./Kc(I,J,K);

        %% Calculate Bz -> Hz
        I = 1:nx; J = 1:ny; K = 2:nz;
        Bz_r = Bz(I,J,K);
        Bz(I,J,K) = k_Bz_a(I,J,K).*Bz(I,J,K) + k_Bz_b(I,J,K).*...
                    ((Ex(I,J+1,K) - Ex(I,J,K))/dy - (Ey(I+1,J,K) - Ey(I,J,K))/dx);
        Hz(I,J,K) = k_Hz_a(I,J,K).*Hz(I,J,K) + k_Hz_b(I,J,K).*...
                    (k_Hz_c(I,J,K).*Bz(I,J,K) + k_Hz_d(I,J,K).*Bz_r)./Kc(I,J,K);

        %% Progress bar updates at every time percent
        if ( mod(T, ceil(number_of_iterations/100)) == 0 )
            waitbar((T+1)/number_of_iterations, hhh);
        end

    end
    toc
    close(hhh);

%% Plots: signals, scattering parameters and Ez field between plates
   figure('units','normalized','outerposition',[0 0 1 1]);
   set(gcf, 'doublebuffer', 'on');
   % Input signal, reflected (at port 1) and passed through (port 2)
     subplot(2,2,1);
     plot(1e9*dt*(0:T), Incident, 'k-', 1e9*dt*(0:T), Reflected, 'r-', 1e9*dt*(0:T), Passed, 'b-', 'LineWidth', 2);
     title('Input, reflected (Port 1) and passed (Port 2) signals', 'FontSize', 18);
     xlabel('Time, [ns]', 'FontSize', 16);
     ylabel('Strength, [V/cm]', 'FontSize', 16);
     legend('Incident (Port 1)' ,'Reflected (Port 1)', 'Passed (Port 2)');
     grid on;
   % |S_11| and |S_21|
     NFFT = 2^nextpow2(T);
     S_11 = fft(Reflected, NFFT)./fft(Incident, NFFT);
     S_21 = fft(Passed, NFFT)./fft(Incident, NFFT);
     subplot(2,2,3);
     f = (0.5/dt)*linspace(0, 1, NFFT/2);
     plot(f*1e-9, 20*log10(abs(S_11(1:NFFT/2))), 'r', ...
          f*1e-9, 20*log10(abs(S_21(1:NFFT/2))), 'b', 'LineWidth', 2);
     title('Scattering parameters - |S_{11}| and |S_{21}|', 'FontSize', 18);
     xlabel('Frequency, [GHz]', 'FontSize', 16);
     ylabel('Magnitude, [dB]', 'FontSize', 16);
     legend('Reflected (Port 1)', 'Passed (Port 2)', 'Location', 'SouthEast');
     xlim([0 20]);
     grid on;
   % Electric field strength Ez between electrodes
     for T=1:20:number_of_iterations
         subplot(2,2,[2 4]);
         contour(1e3*dx*(1:(nx-2*PML-1)), 1e3*dy*(1:(ny-2*PML-1)), ...
         Index((PML+1):(nx-PML-1), (PML+1):(ny-PML-1), PML+5)', 1);
         hold on;
         subplot(2,2,[2 4]);
         pcolor(1e3*dx*(1:(nx-2*PML-1)), 1e3*dy*(1:(ny-2*PML-1)), double(Ez_saved(:, :, T)));
         title(['Time t = ', num2str(round(T*dt*1e12)),' ps'], 'FontSize', 18);
         xlabel('[mm]', 'FontSize', 16);
         ylabel('[mm]', 'FontSize', 16);
         shading interp;
         axis image;
         caxis([-1 1]);
         colorbar;
         drawnow;
     end

%% Gauss function for voltage source
function [res] = Source(t_half, amplitude, t)
% Pulse delay
t0 = 3.0*t_half;
res = amplitude*exp(-0.5*((t-t0)/(t_half/2.35482))^2.0);


text:
%
% Implementation of Exposure Fusion
%
% written by Tom Mertens, Hasselt University, August 2007
% e-mail: tom.mertens@gmail.com
%
% This work is described in
%   "Exposure Fusion"
%   Tom Mertens, Jan Kautz and Frank Van Reeth
%   In Proceedings of Pacific Graphics 2007
%
%
% Usage:
%   result = exposure_fusion(I,m);
%   Arguments:
%     'I': represents a stack of N color images (at double
%       precision). Dimensions are (height x width x 3 x N).
%     'm': 3-tuple that controls the per-pixel measures. The elements
%     control contrast, saturation and well-exposedness, respectively.
%
% Example:
%   'figure; imshow(exposure_fusion(I, [0 0 1]);'
%   This displays the fusion of the images in 'I' using only the well-exposedness
%   measure
%

function [R, W]= exposure_fusion(I,m)

r = size(I,1);
c = size(I,2);
N = size(I,4);

W = ones(r,c,N);

%compute the measures and combines them into a weight map
contrast_parm = m(1);
sat_parm = m(2);
wexp_parm = m(3);

if (contrast_parm > 0)
    W = W.*contrast(I).^contrast_parm;
end
if (sat_parm > 0)
    W = W.*saturation(I).^sat_parm;
end
if (wexp_parm > 0)
    W = W.*well_exposedness(I).^wexp_parm;
end

%normalize weights: make sure that weights sum to one for each pixel
W = W + 1e-12; %avoids division by zero
W = W./repmat(sum(W,3),[1 1 N]);

% create empty pyramid
pyr = gaussian_pyramid_(zeros(r,c,3));
nlev = length(pyr);

% multiresolution blending
for i = 1:N
    % construct pyramid from each input image
	pyrW = gaussian_pyramid_(W(:,:,i));
	pyrI = laplacian_pyramid_(I(:,:,:,i));

    % blend
    for l = 1:nlev
        w = repmat(pyrW{l},[1 1 3]);
        pyr{l} = pyr{l} + w.*pyrI{l};
    end
end

% reconstruct
R = reconstruct_laplacian_pyramid_(pyr);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% contrast measure
function C = contrast(I)
h = [0 1 0; 1 -4 1; 0 1 0]; % laplacian filter
N = size(I,4);
C = zeros(size(I,1),size(I,2),N);
for i = 1:N
    mono = rgb2gray(I(:,:,:,i));
    C(:,:,i) = abs(imfilter(mono,h,'replicate'));
end

% saturation measure
function C = saturation(I)
N = size(I,4);
C = zeros(size(I,1),size(I,2),N);
for i = 1:N
    % saturation is computed as the standard deviation of the color channels
    R = I(:,:,1,i);
    G = I(:,:,2,i);
    B = I(:,:,3,i);
    mu = (R + G + B)/3;
    C(:,:,i) = sqrt(((R - mu).^2 + (G - mu).^2 + (B - mu).^2)/3);
end

% well-exposedness measure
function C = well_exposedness(I)
sig = .2;
N = size(I,4);
C = zeros(size(I,1),size(I,2),N);
for i = 1:N
    R = exp(-.5*(I(:,:,1,i) - .5).^2/sig.^2);
    G = exp(-.5*(I(:,:,2,i) - .5).^2/sig.^2);
    B = exp(-.5*(I(:,:,3,i) - .5).^2/sig.^2);
    C(:,:,i) = R.*G.*B;
end

author: mahmoudnafifi


text:
function varargout = process_test_permutation2( varargin )
% PROCESS_TEST_PERMUTATION2: Permutation two-sample tests (independent).

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, Dimitrios Pantazis, 2016

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Permutation test: Independent';
    sProcess.Category    = 'Stat2';
    sProcess.SubGroup    = 'Test';
    sProcess.Index       = 104;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Statistics';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data',  'results',  'timefreq',  'matrix'};
    sProcess.OutputTypes = {'pdata', 'presults', 'ptimefreq', 'pmatrix'};
    sProcess.nInputs     = 2;
    sProcess.nMinFiles   = 2;

    % === GENERIC EXTRACT OPTIONS
    % Label
    sProcess.options.extract_title.Comment    = 'Select data to test:';
    sProcess.options.extract_title.Type       = 'label';
    % Options
    sProcess = process_extract_values('DefineExtractOptions', sProcess);
    % DISABLE ABSOLUTE VALUE
    sProcess.options.isabs.Value = 0;
    sProcess.options.isnorm.Value = 0;
    sProcess.options.isabs.Hidden = 1;
    sProcess.options.isnorm.Hidden = 1;

    % === EXCLUDE ZERO VALUES
    sProcess.options.iszerobad.Comment = 'Exclude the zero values from the computation';
    sProcess.options.iszerobad.Type    = 'checkbox';
    sProcess.options.iszerobad.Value   = 1;
    % === OUTPUT COMMENT
    sProcess.options.Comment.Comment = 'Comment (empty=default): ';
    sProcess.options.Comment.Type    = 'text';
    sProcess.options.Comment.Value   = '';

    % === TEST: title
    sProcess.options.test_title.Comment    = '
Test statistic:';
    sProcess.options.test_title.Type       = 'label';
    % === TEST: type
    sProcess.options.test_type.Comment = {['Student''s t-test   (equal variance)
t = (mean(A)-mean(B)) / (Sx * sqrt(1/nA + 1/nB))
' ...
                                           'Sx = sqrt(((nA-1)*var(A) + (nB-1)*var(B)) / (nA+nB-2))'], ...
                                          ['Student''s t-test   (unequal variance)
', ...
                                           't = (mean(A)-mean(B)) / sqrt(var(A)/nA + var(B)/nB)'], ...
                                          ['Absolute mean test:    (works with unconstrained sources)
' ...
                                           'T = (|mean(A)|-|mean(B)|) / sqrt(|var(A)|/nA + |var(B)|/nB)']; ...
                                          ...  ['Wilcoxon rank-sum test
', ...
                                          ...  'R = tiedrank([A,B]),    W = sum(R(1:nA))']; ...
                                          'ttest_equal', 'ttest_unequal', 'absmean'}; % , 'wilcoxon'};
    sProcess.options.test_type.Type    = 'radio_label';
    sProcess.options.test_type.Value   = 'ttest_equal';

    % ===== STATISTICAL TESTING OPTIONS =====
    sProcess.options.label2.Comment  = '
Statistical testing (Monte-Carlo):';
    sProcess.options.label2.Type     = 'label';
    % === NUMBER OF RANDOMIZATIONS
    sProcess.options.randomizations.Comment = 'Number of randomizations:';
    sProcess.options.randomizations.Type    = 'value';
    sProcess.options.randomizations.Value   = {1000, '', 0};
    % === TAIL FOR THE TEST STATISTIC
    sProcess.options.tail.Comment  = {'One-tailed (-)', 'Two-tailed', 'One-tailed (+)', ''; ...
                                      'one-', 'two', 'one+', ''};
    sProcess.options.tail.Type     = 'radio_linelabel';
    sProcess.options.tail.Value    = 'two';
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = process_test_parametric2('FormatComment', sProcess);
    Comment = ['Perm ' Comment];
end

%% ===== RUN =====
function sOutput = Run(sProcess, sInputsA, sInputsB) %#ok<DEFNU>
    % Initialize returned variables
    sOutput = [];

    % ===== GET OPTIONS =====
    % Get generic extract options
    OPTIONS = process_extract_values('GetExtractOptions', sProcess, sInputsA(1));
    % Exclude zero values
    if isfield(sProcess.options, 'iszerobad') && isfield(sProcess.options.iszerobad, 'Value') && ~isempty(sProcess.options.iszerobad.Value)
        OPTIONS.isZeroBad = sProcess.options.iszerobad.Value;
    else
        OPTIONS.isZeroBad = 1;
    end
    % Get test type
    OPTIONS.TestType = sProcess.options.test_type.Value;
    OPTIONS.TestTail = sProcess.options.tail.Value;
    Randomizations = sProcess.options.randomizations.Value{1};
    % Invalid test/tail combinations
    if ismember(OPTIONS.TestType, {'ttest_onesample'}) && OPTIONS.isAbsolute && ismember(OPTIONS.TestTail, {'two', 'one-'})
        bst_report('Warning', sProcess, [], 'Testing |X|>0: Using a positive one-tailed test (one+) instead.');
        OPTIONS.TestTail = 'one+';
    elseif strcmpi(OPTIONS.TestType, 'chi2_onesample') && ismember(OPTIONS.TestTail, {'two', 'one-'})
        bst_report('Warning', sProcess, [], 'Testing |X|>0: Using a positive one-tailed test (one+) instead.');
        OPTIONS.TestTail = 'one+';
    elseif strcmpi(OPTIONS.TestType, 'signtest') && ismember(OPTIONS.TestTail, {'two', 'one-'})
        bst_report('Warning', sProcess, [], 'The sign test statistic produces positive values only: Using a positive one-tailed test (one+) instead.');
        OPTIONS.TestTail = 'one+';
    end

    % ===== CHECK INPUT FILES =====
    % Make sure that file type is indentical for both sets
    if ~isempty(sInputsA) && ~isempty(sInputsB) && ~strcmpi(sInputsA(1).FileType, sInputsB(1).FileType)
        bst_report('Error', sProcess, [], 'Cannot process inputs from different types.');
        return;
    end
    % Check the number of files in input
    if (length(sInputsA) < 2)
        bst_report('Error', sProcess, [], 'Not enough files in input.');
        return;
    end
    % Load time vector from the first file: if same as input, discard input
    TimeVector = in_bst(sInputsA(1).FileName, 'Time');
    if ~isempty(OPTIONS.TimeWindow) && (abs(TimeVector(1) - OPTIONS.TimeWindow(1)) < 1e-4) && (abs(TimeVector(end) - OPTIONS.TimeWindow(2)) < 1e-4)
        OPTIONS.TimeWindow = [];
    end
    % Load freq range from the first file: if same as input, discard input
    if ~isempty(OPTIONS.FreqRange)
        % Load Freqs field from the input file
        TfMat = in_bst_timefreq(sInputsA(1).FileName, 0, 'Freqs');
        if iscell(TfMat.Freqs)
            BandBounds = process_tf_bands('GetBounds', TfMat.Freqs);
            FreqList = unique(BandBounds(:));
        else
            FreqList = TfMat.Freqs;
        end
        if (abs(OPTIONS.FreqRange(1) - FreqList(1)) < 1e-4) && (abs(OPTIONS.FreqRange(2) - FreqList(end)) < 1e-4)
            OPTIONS.FreqRange = [];
        end
    end

    % ===== INPUT DATA =====
    % Concatenante values in the 4th dimension
    OPTIONS.Dim = 4;
    % Set all the bad values to 0
    OPTIONS.isBadZero = 1;
    % Call extraction process: FilesA
    [sDataA, OutputType, matName] = process_extract_values('Extract', sProcess, sInputsA, OPTIONS);
    if isempty(sDataA)
        bst_report('Error', sProcess, [], 'No data read from FilesA.');
        return;
    end
    % Read FilesB
    sDataB = process_extract_values('Extract', sProcess, sInputsB, OPTIONS);
    if isempty(sDataB)
        bst_report('Error', sProcess, [], 'No data read from FilesB.');
        return;
    elseif (size(sDataA.(matName),1) ~= size(sDataB.(matName),1)) || (size(sDataA.(matName),2) ~= size(sDataB.(matName),2)) || (size(sDataA.(matName),3) ~= size(sDataB.(matName),3))
        bst_report('Error', sProcess, [], 'Files A and B do not have the same number of signals or time samples.');
        return;
    end
% WARNING ONLY APPLIES TO PARAMETRIC TESTS
%     % Time-frequency: Warning if processing power
%     if strcmpi(sInputsA(1).FileType, 'timefreq') && isfield(sDataA, 'Measure') && isequal(sDataA.Measure, 'power')
%         bst_report('Warning', sProcess, [], ['You are testing power values, while a more standard analysis is to test the magnitude (ie. sqrt(power)).' 10 ...
%             'Option #1: Recompute the time-frequency maps using the option "Measure: Magnitude".' 10 ...
%             'Option #2: Run the process "Extract > Measure from complex values", with option "Magntiude".']);
%     end

    % ===== UNCONSTRAINED SOURCES =====
    % Detect if the source model is unconstrained
    isUnconstrained = panel_scout('isUnconstrained', sDataA);
    % Do not allow unconstrained sources without a norm
    if isUnconstrained && ~OPTIONS.isAbsolute
        % Unconstrained models: Ok if using
        if strcmpi(OPTIONS.TestType, 'absmean') % && (sDataA.nComponents == 3)
            OPTIONS.TestType = 'absmean_unconstr';
        else
            bst_report('Error', sProcess, [], ['Cannot run this test on unconstrained sources:' 10 'you must compute the norm of the three orientations first.']);
            return;
        end
    end

    % === COMPUTE TEST ===
    % Run the permutation test
    [pmap, tmap, nA, nB] = bst_permtest(sDataA.(matName), sDataB.(matName), OPTIONS.TestType, OPTIONS.Dim, Randomizations, OPTIONS.TestTail, OPTIONS.isZeroBad);
    % Finished processing
    bst_progress('text', 'Saving the results...');

    % Bad channels: For recordings, keep only the channels that are good in BOTH A and B sets
    switch lower(sInputsA(1).FileType)
        case 'data'
            ChannelFlag = sDataA.ChannelFlag;
            ChannelFlag(sDataB.ChannelFlag == -1) = -1;
        case {'results', 'timefreq', 'matrix'}
            ChannelFlag = [];
    end

    % === DISPLAYED UNITS ===
    switch (OPTIONS.TestType)
        case 'ttest_equal',      DisplayUnits = 't';
        case 'ttest_unequal',    DisplayUnits = 't';
        case 'ttest_paired',     DisplayUnits = 't';
        case 'signtest',         DisplayUnits = 'N';
        case 'wilcoxon_paired',  DisplayUnits = 'W';
        case 'absmean',          DisplayUnits = 'T';
        case 'absmean_unconstr', DisplayUnits = 'T';
        % case 'wilcoxon',        DisplayUnits = 'W';
        otherwise,              error('Invalid statistic.');
    end

    % === CONVERT BACK MATRIX => DATA ===
    % If processing recordings with only some sensor types selected
    if strcmpi(sInputsA(1).FileType, 'data') && strcmpi(OutputType, 'matrix') && ~isempty(OPTIONS.SensorTypes) && ~OPTIONS.isAvgTime && ~OPTIONS.isAvgRow && ~OPTIONS.isAvgFreq
        % Get the list of selected sensors
        dataTypes = strtrim(str_split(OPTIONS.SensorTypes, ',;'));
        % If only major data types were selected: save results in "data" format
        if ~isempty(dataTypes) && all(ismember(dataTypes, {'MEG','EEG','MEG MAG''MEG GRAD','MEG GRAD2','MEG GRAD3','SEEG','ECOG','NIRS'}))
            % Load channel file
            ChannelMat = in_bst_channel(sInputsA(1).ChannelFile);
            % Find channel names in the output row names
            iChannels = channel_find(ChannelMat.Channel, OPTIONS.SensorTypes);
            % Convert output data matrices
            % tmap
            tmap_tmp = zeros(length(ChannelMat.Channel), size(tmap,2), size(tmap,3));
            tmap_tmp(iChannels,:,:) = tmap;
            tmap = tmap_tmp;
            % pmap
            pmap_tmp = zeros(size(tmap_tmp));
            pmap_tmp(iChannels,:,:) = pmap;
            pmap = pmap_tmp;
            % nA
            if ~isempty(nA)
                nA_tmp = zeros(size(tmap_tmp));
                nA_tmp(iChannels,:,:) = nA;
                nA = nA_tmp;
            end
            % nB
            if ~isempty(nB)
                nB_tmp = zeros(size(tmap_tmp));
                nB_tmp(iChannels,:,:) = nB;
                nB = nB_tmp;
            end
            % New channel flag
            tmpChannelFlag = -1 .* ones(length(ChannelMat.Channel), 1);
            if ~isempty(ChannelFlag) && (length(ChannelFlag) == length(iChannels))
                tmpChannelFlag(iChannels) = ChannelFlag;
            else
                tmpChannelFlag(iChannels) = 1;
            end
            ChannelFlag = tmpChannelFlag;
            % Convert Stat structure
            OutputType = 'data';
            sDataA.RowNames = [];
        end
    end

    % === OUTPUT STRUCTURE ===
    % Initialize output structure
    sOutput = db_template('statmat');
    sOutput.pmap         = pmap;
    sOutput.tmap         = tmap;
    sOutput.df           = [];
    sOutput.Correction   = 'no';
    sOutput.Type         = OutputType;
    sOutput.ChannelFlag  = ChannelFlag;
    sOutput.Time         = sDataA.Time;
    sOutput.ColormapType = 'stat2';
    sOutput.DisplayUnits = DisplayUnits;
    if strcmpi(OPTIONS.TestType, 'absmean_unconstr')
        sOutput.nComponents = 1;
    elseif isfield(sDataA, 'nComponents')
        sOutput.nComponents = sDataA.nComponents;
    end
    if isfield(sDataA, 'GridAtlas')
        sOutput.GridAtlas = sDataA.GridAtlas;
    end
    if isfield(sDataA, 'Freqs')
        sOutput.Freqs = sDataA.Freqs;
    end
    if isfield(sDataA, 'TFmask')
        sOutput.TFmask = sDataA.TFmask;
    end
    % Row names
    if isfield(sDataA, 'Description') && ~isempty(sDataA.Description)
        RowNames = sDataA.Description;
    elseif isfield(sDataA, 'RowNames') && ~isempty(sDataA.RowNames)
        RowNames = sDataA.RowNames;
    else
        RowNames = [];
    end
    if ~isempty(RowNames)
        if strcmpi(OutputType, 'matrix')
            sOutput.Description = RowNames;
        elseif strcmpi(OutputType, 'timefreq')
            sOutput.RowNames = RowNames;
        end
    end
    % Save options
    sOutput.Options = OPTIONS;
    % Save the number of good samples used for both sets:
    % In compressed format (keeping only one value per row, if all the other dimensions are the same)
    if isequal(nA(1:size(nA,1))', mean(mean(nA,2),3))
        sOutput.Options.nGoodSamplesA = nA(1:size(nA,1))';
        sOutput.Options.nGoodSamplesB = nB(1:size(nB,1))';
    % Or saving the full list of good samples
    else
        sOutput.Options.nGoodSamplesA = nA;
        sOutput.Options.nGoodSamplesB = nB;
    end
end

author: brainstorm-tools


text:
function [phi,lambda] = LC_froca(x,y,maxlat,minlat,maxlon,minlon)

    %LC_FROM_CARTESIAN
    %
    %	[phi,lambda] = LC_from_cartesian(x,y,maxlat,minlat,maxlon,minlon)
    %
    %	Function to compute the angular coordinates PHI (latitude) and
    %	LAMBDA (longitude) from the cartesian coordinates X and Y
    %	given the MAXLAT, MINLAT, MAXLON & MINLON of the Lambert
    %	conformal map on which these points have to be mapped.
    %
    %	where * phi: current location latitude
    %	      * lambda: current location longitude
    %	      * x & y : current cartesian coordinates
    %	      * maxlat: maximum latitude limit of the map
    %	      * minlat: minimum latitude limit of the map
    %	      * maxlon: maximum longitude limit of the map
    %	      * minlon: minimum longitude limit of the map
    %	               (remember: West longitude is < 0!)
    %
    %	If the LC_MAP function has been called before this function
    %	AND the same map limits are used, then it is not neccessary to
    %	enter the last 4 arguments:
    %
    %	[phi,lambda] = LC_to_cartesian(x,y)
    %
    %	Source: Equations taken from "Map Projections Used by the
    %	        U.S. Geological Survey" by John P. Snyder
    %	        Geological Survey Bulletin 1532, pg: 101-109.

    report_this_filefun();

    todeg = 180 / pi;

    % set the global variables
    global scale
    global phi0 lambda0 phi1 phi2
    global maxlatg minlatg maxlong minlong
    ZG = ZmapGlobal.Data;
    torad = ZG.torad;
    Re = ZG.Re;

    if nargin == 2
        %get data from global variables
        maxlat = maxlatg; minlat = minlatg;
        maxlon = maxlong; minlon = minlong;
    elseif nargin == 6
        % set the global variable for later use
        maxlatg = maxlat; minlatg = minlat;
        maxlong = maxlon; minlong = minlon;

        % set some constants
        scale = 1;

        % get the Standard Parallels and Center Coordinates
        phi2 = (minlat + ((maxlat + minlat) / 4)) * torad;
        phi1 = (maxlat - ((maxlat + minlat) / 4)) * torad;
        phi0 = (phi1 + phi2) / 2;
        lambda0 = ((minlon + maxlon) / 2) * torad;
    else
        disp('This function requires 2 or 6 input arguments!')
        help LC_to_cartesian
        return
    end

    % compute the constant of the cone: sine_phi0
    tan1 = tan((pi/4) + (phi1/2));
    tan2 = tan((pi/4) + (phi2/2));
    sine_phi0 = log(cos(phi1)/cos(phi2)) / log(tan2/tan1);

    % compute the auxiliary function: psi
    psi = (cos(phi1) * (tan1.^sine_phi0)) / sine_phi0;

    % compute the polar radius to the origin: rho0
    tan0 = tan((pi/4) + (phi0/2));
    rho0 = (Re * psi) / (tan0.^sine_phi0);

    % compute the polar angles: theta
    theta = atan(x ./ (rho0 - y));

    % compute rho (inverse)
    rho = sign(sine_phi0) * sqrt((x.^2) + (rho0 - y).^2);

    % store the latitudes and longitudes in output variables
    arctan = atan((Re * psi ./ rho).^(1/sine_phi0));
    phi = ((2 * arctan) - (pi/2)) * todeg;
    lambda = ((theta / sine_phi0) + lambda0) * todeg;
end

author: CelsoReyes


text:
function I = movieToImages( M )
% Creates a stack of images from a matlab movie M.
%
% Repeatedly calls frame2im. Useful for playback with playMovie.
%
% USAGE
%  I = movieToImages( M )
%
% INPUTS
%  M   - a matlab movie
%
% OUTPUTS
%  I   - MxNxT array (of images)
%
% EXAMPLE
%  load( 'images.mat' );  [X,map]=gray2ind(video(:,:,1));
%  M = fevalArrays( video, @(x) im2frame(gray2ind(x),map) );
%  I = movieToImages(M); playMovie(I);
%
% See also PLAYMOVIE
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.0
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

I = fevalArrays( M, @frame2Ii );

function I = frame2Ii( F )

[I,map] = frame2im( F );
if( isempty(map) )
  if( size(I,3)==3 )
    classname = class( I );
    I = sum(I,3)/3;
    I = feval( classname, I );
  end
else
  I = ind2gray( I, map );
end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function a_sorted = r8vec_sort_heap_a ( n, a )

%*****************************************************************************80
%
%% R8VEC_SORT_HEAP_A ascending sorts an R8VEC using heap sort.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 April 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    A Nijenhuis and H Wilf,
%    Combinatorial Algorithms,
%    Academic Press, 1978, second edition,
%    ISBN 0-12-519260-6.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the array.
%
%    Input, real A(N), the array to be sorted;
%
%    Output, real A_SORTED(N), the sorted array.
%
  if ( n < 1 )
    a_sorted = [];
    return
  end

  if ( n == 1 )
    a_sorted(1) = a(1);
    return
  end
%
%  1: Put A into descending heap form.
%
  a_sorted = r8vec_heap_d ( n, a );
%
%  2: Sort A.
%
%  The largest object in the heap is in A(1).
%  Move it to position A(N).
%
  temp = a_sorted(1);
  a_sorted(1) = a_sorted(n);
  a_sorted(n) = temp;
%
%  Consider the diminished heap of size N1.
%
  for n1 = n-1 : -1 : 2
%
%  Restore the heap structure of A(1) through A(N1).
%
    a_sorted(1:n1) = r8vec_heap_d ( n1, a_sorted );
%
%  Take the largest object from A(1) and move it to A(N1).
%
    temp = a_sorted(1);
    a_sorted(1) = a_sorted(n1);
    a_sorted(n1) = temp;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function rgb = cortex

% returns a predefined color as [red green blue] values
%
% skin_surface         = [255 213 119]/255;
% outer_skull_surface  = [140  85  85]/255;
% inner_skull_surface  = [202 100 100]/255;
% cortex = [255 213 119]/255;
% black  = [0   0   0  ]/255;
% white  = [255 255 255]/255;
% red    = [255 0   0  ]/255;
% green  = [0   192 0  ]/255;
% blue   = [0   0   255]/255;
% yellow = [255 255 0  ]/255;
% cortex_light = [199 194 169]/255;
% cortex_dark  = [100 97 85]/255;

rgb = [255 213 119]/255;

author: spm


text:
%
%  T = build_toep( c, d, n );
%
%  Given:
%    c - the nonzero part of a central column of a banded Toeplitz
%        matrix
%    d - index of c containing the diagonal element of T
%    n - dimension of the banded Toeplitz matrix
%
%  The banded Toeplitz matrix is constructed explicitly.
%
function T = build_toep( c, d, n )

m = length( c );

col = zeros(n,1);
row = col';
col(1:m-d+1,1) = c(d:m);
row(1,1:d) = c(d:-1:1)';
T = toeplitz( col, row );

repo: IRtools


text:
% POLARCONT Polar contour plot
%
% Richard Rieber
% rrieber@gmail.com
% April 4, 2007
% Updated June 15, 2007
%
% function [C,h] = polarcont(r,theta,z,N,s)
%
% Purpose: This function creates polar contour plots on the current active
%          figure
%
% Inputs:  o r     - Radius vector of length m
%          o theta - Angle vector in radians of length n
%          o z     - Magnitude at the points specified in r and theta of
%                    size m x n
%          o N     - The number of contours to plot [OPTIONAL]
%          o s     - Linespec as described in PLOT [OPTIONAL]
%
% Outputs: o C     - returns contour matrix C as described in CONTOURC
%          o h     - Column vector H of handles to LINE or PATCH objects,
%                    one handle per line.
%
% OTHER NOTES:
% - Both C and h can be used as inputs to CLABEL
% - Colors are defined in colormap
% - Treat this function as a standard contour plot

function [C,h] = polarcont(r,theta,z,N,s)

[a,b] = size(z);

if a ~= length(r)
    error('r is not the same length as the first dimension of z')
end

if b ~= length(theta)
    error('theta is not the same length as the second dimension of z')
end

x = zeros(a,b);
y = zeros(a,b);

for j = 1:a
    for k = 1:b
        x(j,k) = r(j)*cos(theta(k));
        y(j,k) = r(j)*sin(theta(k));
    end
end

if nargin == 3
    [C,h] = contour(x,y,z);
elseif nargin == 4
    [C,h] = contour(x,y,z,N);
elseif nargin == 5
    [C,h] = contour(x,y,z,N,s);
else
    error('Incorrect number of inputs')
end

author: Sable

repo: mcbench-benchmarks


text:
% THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION IS RELEASED "AS IS."  THE U.S. GOVERNMENT MAKES NO WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, CONCERNING THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION, INCLUDING, WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.  IN NO EVENT WILL THE U.S. GOVERNMENT BE LIABLE FOR ANY DAMAGES, INCLUDING ANY LOST PROFITS, LOST SAVINGS OR OTHER INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE, OR INABILITY TO USE, THIS SOFTWARE OR ANY ACCOMPANYING DOCUMENTATION, EVEN IF INFORMED IN ADVANCE OF THE POSSIBILITY OF SUCH DAMAGES.
%
% file: prop_diff_ci.m

% 030110 tdr created
% 031223 tdr added test for NaN in inputs
% 040203 tdr corrected error found by Steve Walls in check_inputs (was using n in place of n1 when ensuring that we had integers)

function metrics = prop_diff_ci(x1,n1,x2,n2,alpha,method,verbose)

if (nargin ~= 7 & nargin ~= 6 & nargin ~= 5),
    error('Requires 5, 6, or 7 input arguments (x1,n1,x2,n2,alpha), (x1,n1,x2,n2,alpha,method), or (x1,n1,x2,n2,alpha,method,verbose); default method = 3');
end
if nargin == 5, method = 3; verbose = 0; end;
if nargin == 6, verbose = 0; end;
inputsOkay = check_inputs(x1,n1,x2,n2,alpha,method,verbose);
if ~inputsOkay
    metrics = [NaN, NaN, NaN];
    return
end;

tic;

% call selected interval estimator method
switch method
case 1, % one-sided CI - note provides both upper and lower one-sided CIs
    ci = get_prop_diff_ci1(x1,n1,x2,n2,alpha,method,verbose);
case 2, % minimum length CI
    error('Method 2 Not Available - choose Method 1, 3, 4, 5 or 6');
case 3, % symmetrical about estimate, expect when outside range
    ci = get_prop_diff_ci3(x1,n1,x2,n2,alpha,method,verbose);
case 4, % equal tail masses, note - may not enclose estimate
    ci = get_prop_diff_ci1(x1,n1,x2,n2,alpha/2,method,verbose);
case 5, % Using Fagan's application of "exact" method
    ci = get_ci_fagan(x1,n1,x2,n2,alpha,method,verbose);
case 6, % Normal Approximation
    ci = get_ci_rss(x1,n1,x2,n2,alpha,method,verbose);
case 7, % Law of Large Numbers Approximation
    error('Method 7 Not Available - choose Method 1, 3, 4, 5 or 6');
otherwise
   error('Not a valid method: 1 (1-sided), 2 (min length), 3 (symm.), 4 (equal tail), 5 (MatLab-binofit), 6 (Normal Approx)');
end;
rt = toc;

if verbose
    p1_hat = x1/n1; p2_hat = x2/n2; delta_p_hat = p1_hat - p2_hat;
    length = ci(3)-ci(2);
    lower_tail = 1 - prop_diff(x1,n1,x2,n2,ci(2));
    upper_tail = prop_diff(x1,n1,x2,n2,ci(3));
    if (method == 1)
        actual_alpha = (lower_tail + upper_tail)/2;
    else
        actual_alpha = lower_tail + upper_tail;
    end;
    if (abs(actual_alpha - alpha) > 0.00005) % close_enough is 0.00001 throughout, but spec is 0.0005.
        warning('Interval not to spec: abs(desired_alpha - actual_alpha) < 0.00005)');
        warning_stats = [alpha actual_alpha alpha-actual_alpha]
    end;
    disp('p_diff_hat, Lower CI Bound, Upper CI Bound, x1, n1, x2, n2, Desired alpha, Method, Length, Lower Tail, Upper Tail, Actual alpha, Delta alpha, Run Time')
    metrics = [ci x1 n1 x2 n2 alpha method length lower_tail upper_tail actual_alpha (alpha - actual_alpha) rt];
    disp(sprintf('%8.6g, %8.6g, %8.6g, %d, %d, %d, %d, %8.6g, %d, %8.6g, %8.6g, %8.6g, %8.6g, %8.6g, %8.2g',metrics));
else
    metrics = ci;
end;

% ------------------------------------------------------------
% Difference CI based on root sum of squares of individual CI's, see
% Fagan'99 in Computers in Biology and Medicine.

function  ci = get_ci_fagan(x1,n1,x2,n2,alpha,method,verbose);
    ci1 = prop_ci(x1,n1,alpha,method);
    ci2 = prop_ci(x2,n2,alpha,method);
    diff = ci1(1) - ci2(1);
    correction_factor = 0.18 / max(n1,n2); % see Fagan'99, p.85
    rss_lower = diff - sqrt((ci1(2)-ci1(1))^2 + (ci2(3)-ci2(1))^2);
    rss_upper = diff + sqrt((ci1(3)-ci1(1))^2 + (ci2(2)-ci2(1))^2);
    diff_lower = max((ci1(2)-ci2(3)), (rss_lower - correction_factor));
    diff_upper = min((ci1(3)-ci2(2)), (rss_upper + correction_factor));
    ci = [diff diff_lower diff_upper];

% ------------------------------------------------------------
% Difference CI based on root sum of squares of individual CI's

function  ci = get_ci_rss(x1,n1,x2,n2,alpha,method,verbose);
    ci1 = prop_ci(x1,n1,alpha,method);
    ci2 = prop_ci(x2,n2,alpha,method);
    diff = ci1(1) - ci2(1);
    rss_lower = diff - sqrt((ci1(2)-ci1(1))^2 + (ci2(3)-ci2(1))^2);
    rss_upper = diff + sqrt((ci1(3)-ci1(1))^2 + (ci2(2)-ci2(1))^2);
    ci = [diff rss_lower rss_upper];

% ------------------------------------------------------------
% start checking inputs
function inputsOkay = check_inputs(x1,n1,x2,n2,alpha,method,verbose)

inputsOkay = 1;

% inputs must all be scalars
if (max(max([size(x1); size(n1); size(x2); size(n2); size(alpha)])) ~= 1)
    error('Non-scalar input');
end;

if isnan(x1) | isnan(n1) | isnan(x2) | isnan(n2) | isnan(alpha)
    warning('NaN input')
    inputsOkay = 0;
    return
end;

% x and n must be integers
if (round(x1) ~= x1)
    x1 = round(x1);
    warning('Non-integer input x1');
end;
if (round(x2) ~= x2)
    x2 = round(x2);
    warning('Non-integer input x2');
end;
if (round(n1) ~= n1)
    n1 = round(n1);
    warning('Non-integer input n1');
end;
if (round(n2) ~= n2)
    n2 = round(n2);
    warning('Non-integer input n2');
end;

% n must be > 0
if (n1 <= 0 | n2 <= 0),
   inputsOkay = 0;
   warning('n <= 0');
   return;
end;

% n should be < 10^6
if (n1 > 10^6 | n2 > 10^6),
   warning('n > 10^6, Results may not be accurate.');
end;

% x must be >= 0 and <=n
if ( x1 < 0 | x1 > n1 | x2 < 0 | x2 > n2),
   inputsOkay = 0;
   warning('x < 0 or x > n');
   return;
end;

% alpha must be > 0.0 and < 1.0
if ( alpha < 0 | alpha > 1),
   inputsOkay = 0;
   warning('alpha < 0 or alpha > 1');
   return;
end;

% verbose must be 0 or 1
if ( verbose ~= 0 & verbose ~= 1),
   inputsOkay = 0;
   warning('verbose not zero or one');
   return;
end;
% end checking inputs
% ------------------------------------------------------------


text:
function info = get_performance_bqp(Xopt,yopt,Sopt,SDP,POP,path)

if iscell(path)
    for i = 1:length(path)
        addpath(genpath(path{i}))
    end
else
    addpath(genpath(path))
end

blk         = SDP.blk;
At          = SDP.At;
b           = SDP.b;
C           = svec(blk,SDP.C);
Xoptmat     = Xopt;
Xopt        = svec(blk,Xopt);
Sopt        = svec(blk,Sopt);
Amap  = @(X) AXfun(blk,At,X);
ATmap = @(y) Atyfun(blk,At,y);

%% standard SDP KKT residuals
Rp          = Fnorm(b - Amap(Xopt))/(1+Fnorm(b));
Rd          = Fnorm(ops(ops(C,'-',ATmap(yopt)),'-',Sopt))/(1+Fnorm(C));
pobj        = blktrace(blk,C,Xopt);
dobj        = b'*yopt;
gap         = abs(pobj - dobj)/(1+abs(pobj) + abs(dobj));

%% round a rank-one solution
Xmom        = Xoptmat{1};
[V,~]       = sorteig(Xmom);
v           = V(:,1)/V(1,1);
x_est       = v(2:1+POP.d);
f_est       = POP.f.coefficient * ( prod(x_est.^(POP.f.degmat),1) )';

S_mineig    = mineig( smat(blk, ops(C,'-',ATmap(yopt)) ) );
S_mineig_1  = min(S_mineig,0);
M           = SDP.M;
f_lb        = SDP.b'*yopt + M'*S_mineig_1;
eta         = abs(f_est - f_lb)/(1+abs(f_est)+abs(f_lb));

info.Rp     = Rp;
info.Rd     = Rd;
info.Rg     = gap;
info.Rs     = eta;
info.pobj   = pobj;
info.dobj   = dobj;
info.f_est      = f_est;
info.S_mineig   = S_mineig;
info.f_lb       = f_lb;
info.x_est      = x_est;

fprintf('\n========== Performance of Binary Quadratic Program ==========\n')
fprintf('Rp: %3.2e, Rd: %3.2e, Rg: %3.2e, Rs: %3.2e.\n',Rp,Rd,gap,eta);
fprintf('f_est: %3.4e, f_lb: %3.4e, pobj: %3.4e, dobj: %3.4e.\n',f_est,f_lb,pobj,dobj);
fprintf('==============================================================\n')

if iscell(path)
    for i = 1:length(path)
        rmpath(genpath(path{i}))
    end
else
    rmpath(genpath(path))
end
end

author: MIT-SPARK

repo: CertifiablyRobustPerception


text:
function LU = extract_bounds_from_milpsubsref_operator(LU,extstruct,extvariables,i)
arg = extstruct(i).arg{1};
epi = getvariables(extstruct(i).var);
[M,m] = derivebounds(reshape(arg,[],1),LU);
LU(epi,1) = min(m);
LU(epi,2) = max(M);
for j = 1:length(extstruct(i).arg{2}.subs)
    index = extstruct(i).arg{2}.subs{j};
    if isa(index,'sdpvar')
        if isequal(getbase(index),[0 1])
            LU(getvariables(index),1) = max(LU(getvariables(index),1),1);
            LU(getvariables(index),2) = min(LU(getvariables(index),2),numel(arg));
        end
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [model, varargout] = grlvq_core(trainSet, trainLab, varargin)
%grlvq_core.m - trains the Generalized Relevance LVQ algorithm
%NOTE: minimal requirement version 7.4.0.336 (R2007a)
%
%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  %% Use the wrapper GRLVQ.M to access the functionality in style of %%
%  %% the SOM Toolbox (i.e. with data structs).                       %%
%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  example for usage:
%  trainSet = [1,2,3;4,5,6;7,8,9];
%  trainLab = [1;1;2];
%  GRLVQ_model=GRLVQ_train(trainSet,trainLab); % minimal parameters required
%  estimatedTrainLabels = GRLVQ_classify(trainSet, GRLVQ_model);
%  trainError = mean( trainLab ~= estimatedTrainLabels );
%
% input:
%  trainSet : matrix with training samples in its rows
%  trainLab : vector with the labels of the training set
% optional parameters:
%  PrototypesPerClass: (default=1) the number of prototypes per class used. This could
%  be a number or a vector with the number for each class
%  initialPrototypes : (default=[]) a set of prototypes to start with. If not given initialization near the class means
%  initialRelevances : the relevances to start with. If not given all
%  relevances are equal at start
%  regularization    : (default=0) values usually between 0 and 1 treat with care.
%  Regularizes the relevances to be more homogeneous
%  testSet           : (default=[]) an optional test set used to compute
%  the test error. The last column is expected to be a label vector
%  comparable        : (default=0) a flag which resets the random generator
%  to produce comparable results if set to 1
%  optimization      : (default=fminlbfgs) indicates which optimization is used: sgd or fminlbfgs
% parameter for the stochastic gradient descent sgd
%  nb_epochs             : (default=100) the number of epochs for sgd
%  learningRatePrototypes: (default=[]) the learning rate for the prototypes.
%  Could be the start and end value used for a sigmoidal spectrum or a vector of length nb_epochs
%  learningRateRelevances: (default=[]) the learning rate for the matrix.
%  Could be the start and end value used for a sigmoidal spectrum or a vector of length nb_epochs
%  relevanceStart        : (default=1) the epoch to start the matrix training
% parameter for the build-in function fminlbfgs
%  threshstop       : (default=0) the training error for early stopping
%  useEarlyStopping : (default=1) use early stopping based on threshstop
%  Display          : (default='off') the optimization output 'iter' or 'off'
%  GradObj          : (default=on) use the gradient information or not
%  HessUpdate       : (default=lbfgs) the update can be 'lbfgs', 'bfgs' or 'steepdesc'
%  TolFun           : (default=1e-6) the tolerance
%  MaxIter          : (default=2500) the maximal number of iterations
%  MaxFunEvals      : (default=1000000) the maximal number of function evaluations
%  TolX             : (default=1e-10) tolerance
%  DiffMinChange    : (default=1e-10) minimal change
%
% output: the GRLVQ model with prototypes w their labels c_w and the relevances lamda
%  optional output:
%  initialization : a struct containing the settings
%  trainError     : error in the training set
%  testError      : error in the training set (only computed if 'testSet' is given)
%  costs          : the output of the cost function
%
% Citation information:
% B. Hammer and T. Villmann, Generalized relevance learning vector quantization, Neural Networks 15, 1059-1068, 2002.
%
% P. Schneider, K. Bunte, B. Hammer and M. Biehl, Regularization in Matrix Relevance Learning,
% IEEE Transactions on Neural Networks, vol. 21, nb. 5, pp. 831-840, 2010.
%
% Kerstin Bunte
% uses the Fast Limited Memory Optimizer fminlbfgs.m written by Dirk-Jan Kroon available at the MATLAB central
% kerstin.bunte@googlemail.com
% Fri Nov 09 14:13:52 CEST 2012
%
% Conditions of GNU General Public License, version 2 and BSD License apply.
% See file 'license-gpl2.txt' and 'BSD_license.txt' enclosed in this package.
% Programs are not for use in critical applications!
%

% Contributed to SOM Toolbox vs2, December 3rd, 2012 by Alexander Schulz
% Copyright (c) Kerstin Bunte
% http://www.cis.hut.fi/projects/somtoolbox/

nout = max(nargout,1)-1;
p = inputParser;   % Create an instance of the class.
p.addRequired('trainSet', @isfloat);
p.addRequired('trainLab', @(x) length(x)==size(trainSet,1) & isnumeric(x));

p.addOptional('testSet', [], @(x)(size(x,2)-1)==size(trainSet,2) & isfloat(x));
p.addParamValue('PrototypesPerClass', ones(1,length(unique(trainLab))), @(x)(sum(~(x-floor(x)))/length(x)==1 && (length(x)==length(unique(trainLab)) || length(x)==1)));
p.addParamValue('initialPrototypes',[], @(x)(size(x,2)-1==size(trainSet,2) && isfloat(x)));
p.addParamValue('initialRelevances',[], @(x)(numel(x)==size(trainSet,2) && isfloat(x)));
p.addOptional('comparable', 0, @(x)(~(x-floor(x))));
p.addParamValue('regularization',0, @(x)(isfloat(x) && x>=0));
p.addOptional('optimization', 'fminlbfgs', @(x)any(strcmpi(x,{'sgd','fminlbfgs'})));
% parameter for the stochastic gradient descent
p.addOptional('nb_epochs', 100, @(x)(~(x-floor(x))));
p.addParamValue('learningRatePrototypes', [], @(x)(isfloat(x) || isa(x,'function_handle')));
p.addParamValue('learningRateRelevances', [], @(x)(isfloat(x) || isa(x,'function_handle')));
p.addOptional('relevanceStart', 1, @(x)(~(x-floor(x))));
% parameter for the build-in function
p.addOptional('threshstop',0,@(x) isfloat(x));
p.addOptional('useEarlyStopping',1,@(x)(~(x-floor(x))));
p.addOptional('Display', 'off', @(x)any(strcmpi(x,{'iter','off'})));
p.addOptional('GradObj', 'on', @(x)any(strcmpi(x,{'on','off'})));
p.addOptional('HessUpdate', 'lbfgs', @(x)any(strcmpi(x,{'lbfgs', 'bfgs', 'steepdesc'})));
p.addOptional('TolFun',1e-6,@(x) isfloat(x));
p.addOptional('MaxIter', 2500, @(x)(~(x-floor(x))));
p.addOptional('MaxFunEvals', 1000000, @(x)(~(x-floor(x))));
p.addOptional('TolX',1e-10,@(x) isfloat(x));
p.addOptional('DiffMinChange',1e-10,@(x)isfloat(x));
p.CaseSensitive = true;
p.FunctionName = 'GRLVQ';
% Parse and validate all input arguments.
p.parse(trainSet, trainLab, varargin{:});

%%% check if results should be comparable
if p.Results.comparable,
    rng('default');
end
%%% set useful variables
nb_samples = size(trainSet,1);
nb_features = size(trainSet,2);
if size(trainLab,1)~=nb_samples, trainLab = trainLab';end

classes = unique(trainLab);
nb_classes = length(classes);
testSet = p.Results.testSet;

% global regularization;
regularization = p.Results.regularization;
if regularization, disp(['Regularize the relevances with ',num2str(regularization)]);end

initialization = rmfield(p.Results, 'trainSet');
initialization.trainSet = [num2str(nb_samples),'x',num2str(nb_features),' matrix'];
initialization = rmfield(initialization, 'trainLab');
initialization.trainLab = ['vector of length ',num2str(length(trainLab))];
if ~isempty(testSet)
    initialization = rmfield(initialization, 'testSet');
    initialization.testSet = [num2str(size(testSet,1)),'x',num2str(size(testSet,2)),' matrix'];
end
switch(p.Results.optimization)
case{'sgd'}
    initialization = rmfield(initialization, 'useEarlyStopping');
    initialization = rmfield(initialization, 'Display');
    initialization = rmfield(initialization, 'GradObj');
    initialization = rmfield(initialization, 'HessUpdate');
    initialization = rmfield(initialization, 'TolFun');
    initialization = rmfield(initialization, 'MaxIter');
    initialization = rmfield(initialization, 'MaxFunEvals');
    initialization = rmfield(initialization, 'TolX');
    initialization = rmfield(initialization, 'DiffMinChange');
case{'fminlbfgs'}
    disp('The fminlbfgs optimization uses some global variables:');
    disp('threshstop earlystopped useEarlyStopping');
    initialization = rmfield(initialization, 'nb_epochs');
    initialization = rmfield(initialization, 'learningRatePrototypes');
    initialization = rmfield(initialization, 'learningRateRelevances');
    initialization = rmfield(initialization, 'relevanceStart');
end
% Display all arguments.
%disp 'Settings for GRLVQ:'
%disp(initialization);

%%% check the number of prototypes per class if one integer is given and turn
%%% it into a vector
nb_ppc = p.Results.PrototypesPerClass;
if length(nb_ppc)~=nb_classes,
    nb_ppc = ones(1,nb_classes)*nb_ppc;
end

%%% initialize the prototypes
if isempty(p.Results.initialPrototypes)
    % initialize near the class centers
    w = zeros(sum(nb_ppc),nb_features);
    c_w = zeros(sum(nb_ppc),1);
    actPos = 1;
    for actClass=1:nb_classes
        nb_prot_c = nb_ppc(actClass);
        classMean = mean(trainSet(trainLab==classes(actClass),:));
        % set the prototypes to the class mean and add a random variation
        % between -0.1 and 0.1
        w(actPos:actPos+nb_prot_c-1,:) = classMean(ones(nb_prot_c,1),:)+(rand(nb_prot_c,nb_features)*2-ones(nb_prot_c,nb_features))/10;
        c_w(actPos:actPos+nb_prot_c-1) = classes(actClass);
        actPos = actPos+nb_prot_c;
    end
else
    % initialize with given w
    w = p.Results.initialPrototypes(:,1:end-1);
    c_w = p.Results.initialPrototypes(:,end);
end
%%% initialize the matrix
if isempty(p.Results.initialRelevances)
    lambda = ones(1,nb_features);
else
%     if isvector(p.Results.initialRelevances)
%         omega = diag(p.Results.initialRelevances);
%     else
    lambda = p.Results.initialRelevances;
%     end
end
% normalize the relevances
lambda = lambda / sum(lambda);
model = struct('w',w,'c_w',c_w,'lambda',lambda);
clear w c_w lambda;

dim = length(model.lambda);

switch(p.Results.optimization)
case{'sgd'}
    %%% gradient descent variables
    nb_epochs = p.Results.nb_epochs;
    relevanceStart = p.Results.relevanceStart;
    % compute the vector of nb_epochs learning rates alpha for the prototype learning
    if isa(p.Results.learningRatePrototypes,'function_handle')
        % with a given function specified from the user
        alphas = arrayfun(p.Results.learningRatePrototypes, 1:nb_epochs);
    elseif length(p.Results.learningRatePrototypes)>2
        if length(p.Results.learningRatePrototypes)==nb_epochs
            alphas = p.Results.learningRatePrototypes;
        else
            disp('The learning rate vector for the prototypes does not fit the nb of epochs');
            return;
        end
    else
        % or use an decay with a start and a decay value
        if isempty(p.Results.learningRatePrototypes)
            initialization.learningRatePrototypes = [nb_features/100, nb_features/10000];
        end
        alpha_start = initialization.learningRatePrototypes(1);
        alpha_end = initialization.learningRatePrototypes(2);
        alphas = arrayfun(@(x) alpha_start * (alpha_end/alpha_start)^(x/nb_epochs), 1:nb_epochs);
    %     alphas = arrayfun(@(x) alpha_start / (1+(x-1)*alpha_end), 1:nb_epochs);
    end
    % compute the vector of nb_epochs learning rates epsilon for the Matrix learning
    epsilons = zeros(1,nb_epochs);
    if isa(p.Results.learningRateRelevances,'function_handle')
        % with a given function specified from the user
    % 	epsilons = arrayfun(p.Results.learningRateRelevances, 1:nb_epochs);
        epsilons(relevanceStart:nb_epochs) = arrayfun(p.Results.learningRateRelevances, relevanceStart:nb_epochs);
    elseif length(p.Results.learningRateRelevances)>2
        if length(p.Results.learningRateRelevances)==nb_epochs
            epsilons = p.Results.learningRateRelevances;
        else
            disp('The learning rate vector for the Matrix does not fit the nb of epochs');
            return;
        end
    else
        % or use an decay with a start and a decay value
        if isempty(p.Results.learningRateRelevances)
            initialization.learningRateRelevances = [nb_features/1000, nb_features/100000];
        end
        eps_start = initialization.learningRateRelevances(1);
        eps_end = initialization.learningRateRelevances(2);
    %     epsilons = arrayfun(@(x) eps_start * (eps_end/eps_start)^(x/nb_epochs), 1:nb_epochs);
        epsilons(relevanceStart:nb_epochs) = arrayfun(@(x) eps_start * (eps_end/eps_start)^((x-relevanceStart)/(nb_epochs-relevanceStart)), relevanceStart:nb_epochs);
    end

    %%% initialize requested outputs
    trainError = [];
    costs = [];
    testError = [];
    if nout>=2,
        % train error requested
        trainError = ones(1,nb_epochs+1);
        estimatedLabels = GRLVQ_classify(trainSet, model); % error after initialization
        trainError(1) = sum( trainLab ~= estimatedLabels )/nb_samples;
        if nout>=3,
            % test error requested
            if isempty(testSet)
                testError = [];
                disp('The test error is requested, but no labeled test set given. Omitting the computation.');
            else
                testError = ones(1,nb_epochs+1);
                estimatedLabels = GRLVQ_classify(testSet(:,1:end-1), model); % error after initialization
                testError(1) = sum( testSet(:,end) ~= estimatedLabels )/length(estimatedLabels);
            end
            if nout>=4,
                % costs requested
                disp('The computation of the costs is an expensive operation, do it only if you really need it!');
                costs(1) = GRLVQ_costfun(trainSet, trainLab, model, regularization);
            end
        end
    end

    %%% optimize with stochastic gradient descent
    for epoch=1:nb_epochs
        if mod(epoch,100)==0, disp(epoch); end
        % generate order to sweep through the trainingset
        order = randperm(nb_samples);

        % perform one sweep through trainingset
        for i=1:nb_samples
            % select one training sample randomly
            xi = trainSet(order(i),:);
            c_xi = trainLab(order(i));
%             dist = ((xi(ones(size(model.w,1),1),:))-model.w)*model.omega'*(model.omega*((xi(ones(size(model.w,1),1),:))-model.w)');
%             dist = diag(dist);
            dist = sum(bsxfun(@times,bsxfun(@minus,xi,model.w).^2,model.lambda), 2);
            % determine the two winning prototypes
            % nearest prototype with the same class
            [sortDist,sortIdx] = sort(dist);
            count = 1;
            J = sortIdx(count);
            while model.c_w(sortIdx(count)) ~= c_xi,
                count = count+1;
                J = sortIdx(count);
            end
            dJ = sortDist(count);
            count = 1;
            K = sortIdx(count);
            while model.c_w(sortIdx(count)) == c_xi,
                count = count+1;
                K = sortIdx(count);
            end
            dK = sortDist(count);

            wJ = model.w(J,:);
            wK = model.w(K,:);
            % prototype update
            norm_factor = (dJ + dK)^2;
            DJ = (xi-wJ);
            DK = (xi-wK);

            dwJ = (2*dK/norm_factor)*2*model.lambda.*DJ;
            dwK = (2*dJ/norm_factor)*2*model.lambda.*DK;

            model.w(J,:) = wJ + alphas(epoch) * dwJ;
            model.w(K,:) = wK - alphas(epoch) * dwK;

            % update relevances
            if epsilons(epoch)>0, % epoch >= RelevanceStart
%                 f1 = (2*dK/norm_factor)*2*(model.omega*DJ')*DJ;
%                 f2 = (2*dJ/norm_factor)*2*(model.omega*DK')*DK;
                f1 = (2*dK/norm_factor)*DJ.^2;
                f2 = (2*dJ/norm_factor)*DK.^2;
                % update lambda
                if regularization,
                    f3 = diag(pinv(sqrt(diag(model.lambda))))';
%                     test = bsxfun(@times,model.lambda,ones(1,length(model.lambda))');
%                     test(1:nb_features+1:nb_features*nb_features) = 1;
%                     f3 = 1/(prod(model.lambda)+eps)*prod(test,2)';
                else
                    f3 = 0;
                end
                model.lambda = model.lambda-epsilons(epoch) * (f1-f2  - regularization * f3);
                % normalization
                model.lambda (model.lambda < 0) = 0;
                model.lambda = model.lambda / sum(model.lambda);
            end
        end
        if nout>=2,
            % train error requested
            estimatedLabels = GRLVQ_classify(trainSet, model); % error after epoch
            trainError(epoch+1) = sum( trainLab ~= estimatedLabels )/nb_samples;
            if nout>=3,
                % test error requested
                if ~isempty(testSet)
                    estimatedLabels = GRLVQ_classify(testSet(:,1:end-1), model); % error after initialization
                    testError(epoch+1) = sum( testSet(:,end) ~= estimatedLabels )/length(estimatedLabels);
                end
                if nout>=4,
                    % costs requested
%                     dist = computeDistance(trainSet, model.w, model);
%                     if regularization,
%                         regTerm = regularization * log(det(model.omega*model.omega'));
%                     else
%                         regTerm = 0;
%                     end
%                     costs(epoch+1) = sum(arrayfun(@(idx) GLVQ_costfun(min(dist(idx,model.c_w == trainLab(idx))),...
%                                                                       min(dist(idx,model.c_w ~= trainLab(idx)))), 1:size(dist,1)))-regTerm;
                    costs(epoch+1) = GRLVQ_costfun(trainSet, trainLab, model, regularization);
                end
            end
        end
    end
case{'fminlbfgs'}
    %%% optimization options
    options = struct( ...
      'Display',p.Results.Display, ...
      'GradObj',p.Results.GradObj, ...
      'GradConstr',false, ...
      'GoalsExactAchieve',0, ...
      'TolFun',p.Results.TolFun, ...
      'MaxIter',p.Results.MaxIter, ...
      'MaxFunEvals', p.Results.MaxFunEvals, ...
      'TolX',p.Results.TolX, ...
      'DiffMinChange',p.Results.DiffMinChange, ...
      'OutputFcn','LVQ_progresser', ...
      'HessUpdate',p.Results.HessUpdate ...
    );
    clear('progresser'); % memory therein might need reset
    global threshstop earlystopped useEarlyStopping % for LVQ_progresser.m datval labval n_vec
    useEarlyStopping = p.Results.useEarlyStopping; % use early stopping
%     useEarlyStopping = []; % if empty, no validation set given/used for early stopping
    earlystopped = false;
    threshstop = p.Results.threshstop; % stop if classification below this threshold for early stopping
%     global training_data training_label LabelEqualsPrototype LRrelevances LRprototypes prototypeLabel;
%     training_data = trainSet;
%     training_label = trainLab;
%     prototypeLabel = model.c_w;
    nb_prototypes = numel(model.c_w);
%     nb_samples = size(training_data,1);
%     LabelEqualsPrototype = trainLab*ones(1,nb_prototypes) == (model.c_w*ones(1,nb_samples))';
    LabelEqualsPrototype = bsxfun(@eq,trainLab,model.c_w');
    earlystopped = false; % don't change, assigned in progresser.m
    clear('progresser'); % memory therein might need reset
    newfval = realmax('single');
    % fminlbfgs optimizer courtesy of Dirk-Jan Kroon:
    % http://www.mathworks.de/matlabcentral/fileexchange/23245
    % early stopping to be implemented in progresser function
    variables = zeros(nb_prototypes+nb_features,size(trainSet,2));
    variables(1:nb_prototypes,:) = model.w;
    variables(nb_prototypes+1:end,:) = diag(sqrt(model.lambda));
    LRprototypes = 1; % learn prototype locations
    LRrelevances = 0; % don't learn metric
%     [variables,fval] = fminlbfgs(@GRLVQ_optfun,variables,options);
    [variables,fval] = fminlbfgs(@(variables) GMLVQ_optfun(variables,trainSet,LabelEqualsPrototype,LRrelevances,LRprototypes,model.c_w,regularization),variables,options);
    if not(isempty(fval))
        newfval = fval;
    end
    LRprototypes = 0; % don't learn prototype locations
    LRrelevances = 1; % learn metric
%     [variables,fval] = fminlbfgs(@GRLVQ_optfun,variables,options);
    [variables,fval] = fminlbfgs(@(variables) GMLVQ_optfun(variables,trainSet,LabelEqualsPrototype,LRrelevances,LRprototypes,model.c_w,regularization),variables,options);
    if not(isempty(fval))
        newfval = fval;
    end
    if not(isempty(fval))
      clear('progresser'); % memory therein might need reset
      LRprototypes = 1;
      for i = 1:100  % depending on data, re-iterations might further improve
%         [variables,fval] = fminlbfgs(@GRLVQ_optfun,variables,options);
        [variables,fval] = fminlbfgs(@(variables) GMLVQ_optfun(variables,trainSet,LabelEqualsPrototype,LRrelevances,LRprototypes,model.c_w,regularization),variables,options);
        if not(isempty(fval))
          if abs(fval - newfval) < 1e-3
            newfval = fval;
            break
          end
          newfval = fval;
        end
        if isempty(fval) || earlystopped
          break
        end
      end
    end
    model.w = variables(1:nb_prototypes,:);
    model.lambda = diag(variables(nb_prototypes+1:end,:)'*variables(nb_prototypes+1:end,:))';
%     model.lambda = model.lambda/sum(model.lambda);
    if nout>=2,
        % train error requested
        estimatedLabels = GRLVQ_classify(trainSet, model); % error after initialization
        trainError = mean( trainLab ~= estimatedLabels );
        if nout>=3,
            % test error requested
            if isempty(testSet)
                testError = [];
                disp('The test error is requested, but no labeled test set given. Omitting the computation.');
            else
                estimatedLabels = GRLVQ_classify(testSet(:,1:end-1), model); % error after initialization
                testError = mean( testSet(:,end) ~= estimatedLabels );
            end
            if nout>=4,
                % costs requested
%                 dist = computeDistance(trainSet, model.w, model);
%                 if regularization,
%                     regTerm = regularization * log(det(model.omega*model.omega'));
%                 else
%                     regTerm = 0;
%                 end
%                 costs = sum(arrayfun(@(idx) GLVQ_costfun(min(dist(idx,model.c_w == trainLab(idx))),...
%                                                          min(dist(idx,model.c_w ~= trainLab(idx)))), 1:size(dist,1)))-regTerm;
                costs = GRLVQ_costfun(trainSet, trainLab, model, regularization);
            end
        end
    end
end

%%% output of the training
varargout = cell(nout);
for k=1:nout
	switch(k)
		case(1)
			varargout(k) = {initialization};
		case(2)
			varargout(k) = {trainError};
		case(3)
			varargout(k) = {testError};
		case(4)
            varargout(k) = {costs};
	end
end

function cost = GRLVQ_costfun(trainSet, trainLab, model, regularization)
%GRLVQ_costfun.m - computes the costs for a given training set and GRLVQ
%model with or without regularization
%  example for usage:
%  trainSet = [1,2,3;4,5,6;7,8,9];
%  trainLab = [1;1;2];
%  GMLVQ_model=GMLVQ_train(trainSet,trainLab); % minimal parameters required
%  costs = GMLVQ_costfun(trainSet, trainLab, GMLVQ_model, 0);
%
% input:
%  trainSet : matrix with training samples in its rows
%  trainLab : a vector of training labels
%  model    : GMLVQ model with prototypes w their labels c_w and the matrix omega
%  regularization: the factor>=0 for the regularization
%
% output    : cost function value
%
% Kerstin Bunte (based on the code from Marc Strickert)
% kerstin.bunte@googlemail.com
% Mon Nov 05 09:05:52 CEST 2012
%
% Conditions of GNU General Public License, version 2 apply.
% See file 'license-gpl2.txt' enclosed in this package.
% Programs are not for use in critical applications!
%
nb_samples = length(trainLab);
% labels should be a row vector
if size(trainLab,1)~=nb_samples, trainLab = trainLab';end

% LabelEqPrototype = trainLab*ones(1,numel(model.c_w)) == (model.c_w*ones(1,nb_samples))';
LabelEqPrototype = bsxfun(@eq,trainLab,model.c_w');
dists = computeDistance(trainSet, model.w, model);
if regularization,
    regTerm = regularization * log(prod(model.lambda));
% if strcmp(p.Results.optimization,'sgd')
else
    regTerm = 0;
end
Dwrong = dists;
Dwrong(LabelEqPrototype) = realmax(class(Dwrong));   % set correct labels impossible
distwrong = min(Dwrong.'); % closest wrong
clear Dwrong;

Dcorrect = dists;
Dcorrect(~LabelEqPrototype) = realmax(class(Dcorrect)); % set wrong labels impossible
distcorrect = min(Dcorrect.'); % closest correct
clear Dcorrect;
clear dists;
distcorrectpluswrong = distcorrect + distwrong;
distcorrectminuswrong = distcorrect - distwrong;
mu = distcorrectminuswrong ./ distcorrectpluswrong;
if regularization,
    regTerm = regularization * log(prod(model.lambda));
else
    regTerm = 0;
end
cost = sum(mu)-regTerm;

function [f G]  = GMLVQ_optfun(variables,training_data,LabelEqualsPrototype,LRrelevances,LRprototypes,prototypeLabel,regularization)
% [f G] = GMLVQ_optfun(variables)
% function to be optimzed by matrix relevance learning vector quantization
% variables = [prototype matrix;omega matrix]
% global variables are
%   training_data        : data vectors as row vectors, i.e. attributes in columns
%   LabelEqualsPrototype : binary matrix indicating coocurrences of
%   training labels and prototype labels
%   prototypeLabel       : label vector for the prototypes
%   regularization       : the regularization parameter
%   LRrelevances         : learning rate for the relevance matrix
%   LRprototypes         : learning rate for the prototypes
%
% Kerstin Bunte (modified based on the code of Marc Strickert http://www.mloss.org/software/view/323/)
% kerstin.bunte@googlemail.com
% Fri Nov 09 14:13:52 CEST 2012
%
% Conditions of GNU General Public License, version 2 apply.
% See file 'license-gpl2.txt' enclosed in this package.
% Programs are not for use in critical applications!
%
if isempty(LRprototypes) % values between 1e-2,1e-3,... 1e-8 seem pragmatic
    LRprototypes = 1; % no relevance learning by default
end
if isempty(LRrelevances) % values between 1e-2,1e-3,... 1e-8 seem pragmatic
    LRrelevances = 0; % no relevance learning by default
end
[n_data, n_dim] = size(training_data);
nb_prototypes =  numel(prototypeLabel);
omegaT = variables(nb_prototypes+1:end,:)';
n_vec = size(variables,1) - nb_prototypes;

dists = squaredEuclidean(training_data*omegaT, variables(1:nb_prototypes,:)*omegaT);

Dwrong = dists;
Dwrong(LabelEqualsPrototype) = realmax(class(Dwrong));   % set correct labels impossible
[distwrong pidxwrong] = min(Dwrong.'); % closest wrong
clear Dwrong;

Dcorrect = dists;
Dcorrect(~LabelEqualsPrototype) = realmax(class(Dcorrect)); % set wrong labels impossible
[distcorrect pidxcorrect] = min(Dcorrect.'); % closest correct
clear Dcorrect;

distcorrectpluswrong = distcorrect + distwrong;
distcorrectminuswrong = distcorrect - distwrong;
mu = distcorrectminuswrong ./ distcorrectpluswrong;
% callitq = 1./(1 + exp(-squashsigmoid * mu)); % apply sigmoidal

if regularization,
    regTerm = regularization * log(det(omegaT'*omegaT));
else
    regTerm = 0;
end
f = sum(mu)-regTerm;
% f = mean(callitq);

if nargout > 1  % gradient needed not just function eval
    G = zeros(size(variables)); % initially no gradient
    %       callitq = squashsigmoid * callitq .* (1-callitq); % derivative of sigmoid
    %       distcorrectpluswrong = 2 * callitq ./ distcorrectpluswrong.^2; % degeneration?
    distcorrectpluswrong = 4 ./ distcorrectpluswrong.^2; % norm_factor for derivative for every data sample
    if LRrelevances > 0
        Gw = zeros(n_vec,n_dim);
    end
    for k=1:nb_prototypes%(n_vec+1):size(lambda,1) % update all prototypes
        idxc = (k == pidxcorrect);  % Js: idxs where actual prototype is nearest correct
        idxw = (k == pidxwrong);    % Ks: idxs where actual prototype is nearest wrong

        dcd =  distcorrect(idxw) .* distcorrectpluswrong(idxw);
        dwd =    distwrong(idxc) .* distcorrectpluswrong(idxc);
        if LRrelevances > 0
            % part of derivative of distance
            difc = bsxfun(@minus,training_data(idxc,:),variables(k,:)); % DJs
            difw = bsxfun(@minus,training_data(idxw,:),variables(k,:)); % DKs
            % update omega
            Gw = Gw - (bsxfun(@times,difw,dcd.') * omegaT).' * difw + ...
                      (bsxfun(@times,difc,dwd.') * omegaT).' * difc;
            if LRprototypes > 0
                G(k,:) = dcd * difw - dwd * difc;
            end
        else
            if LRprototypes > 0
                G(k,:) = dcd * training_data(idxw,:) - dwd * training_data(idxc,:) + (sum(dwd)-sum(dcd)) * variables(k,:);
            end
        end
    end
if regularization,
    f3 = (pinv(omegaT'))';
else
    f3 = 0;
end
    % some rescalings needed
    if LRrelevances > 0
        G(nb_prototypes+1:nb_prototypes+n_vec,:) = 2/n_data * LRrelevances * Gw - regularization*f3;
    end
    if LRprototypes > 0
        G(1:nb_prototypes,:) = 1./n_data * LRprototypes * G(1:nb_prototypes,:) * omegaT * omegaT.';
    end
    G = G .* (1 + .0001 * (rand(size(G))-.5)); % help break symmetries
end
% if 0,
% w = variables(1:nb_prototypes,:);
% dJs = zeros(1,nb_samples);
% dKs = zeros(1,nb_samples);
% Js = zeros(1,nb_samples);
% Ks = zeros(1,nb_samples);
% norm_factors = zeros(1,nb_samples);
% DJs = zeros(nb_samples,size(training_data,2));
% DKs = zeros(nb_samples,size(training_data,2));
% for i=1:nb_samples
%     % select one training sample randomly
%     xi = training_data(i,:);
%     c_xi = trainLab(i);
%
%     dist = ((xi(ones(size(w,1),1),:))-w)*omegaT*(omegaT'*((xi(ones(size(w,1),1),:))-w)');
%     dist = diag(dist);
%     % determine the two winning prototypes
%     % nearest prototype with the same class
%     [sortDist,sortIdx] = sort(dist);
%     count = 1;
%     J = sortIdx(count);
%     while prototypeLabel(sortIdx(count)) ~= c_xi,
%         count = count+1;
%         J = sortIdx(count);
%     end
%     dJ = sortDist(count);
%     dJs(i) = dJ;
%     Js(i) = J;
%     count = 1;
%     K = sortIdx(count);
%     while prototypeLabel(sortIdx(count)) == c_xi,
%         count = count+1;
%         K = sortIdx(count);
%     end
%     dK = sortDist(count);
%     dKs(i) = dK;
%     Ks(i) = K;
%
%     wJ = w(J,:);
%     wK = w(K,:);
%     % prototype update
%     norm_factors(i) = 4/((dJ + dK)^2);
%     DJ = (xi-wJ);
%     DK = (xi-wK);
%     DJs(i,:) = DJ;
%     DKs(i,:) = DK;
% end
% end

function D = squaredEuclidean(A, B)
% computes the sqared Euclidean distance
%
%   D = squaredEuclidean(X) returns the squared Euclidean distance matrix of data in rows of X
%   D = squaredEuclidean(X, Y) returns the distance matrix with all distances between the points in X and Y.
%
if nargin == 1 % means that one matrix
    D = bsxfun(@plus, sumsquared(A,2), bsxfun(@minus, sumsquared(A,2).', 2*A*A.')); % 2*(Y*Y.')
else
    D = bsxfun(@plus, sumsquared(A,2), bsxfun(@minus, sumsquared(B,2).', 2*A*B.'));
end
D = max(D,0);

function sq  = sumsquared(x,dim)
% sq  = sumsquared(x,dim)
% sum of all squared element of matrix x along dimension dim

persistent isoctave

if isempty(isoctave)
  isoctave = exist('OCTAVE_VERSION','builtin');
end

if nargin == 1
  dim = 1;
end

if isoctave
  sq = sumsq(x,dim);
else
  sq = sum(x.^2, dim);
end

author: ilarinieminen

repo: SOM-Toolbox


text:
function [x, fx, exitFlag] = bisection_s(f,lb,ub,target,options)%#codegen
% BISECTION Fast and robust root-finding method that handles n-dim arrays.
%   bisection_s => "s" = single or serial, for turning into a mex file
%   [x,fVal,ExitFlag] = BISECTION(f,LB,UB,target,options) finds x +/- TolX
%   (LB < x < UB) such that f(x) = target +/- TolFun.
%
%   x = BISECTION(f,LB,UB) finds the root(s) of function f on the interval
%   [LB, UB], i.e. finds x such that f(x) = 0 where LB <= x <= UB. f will
%   never be evaluated outside of the interval specified by LB and UB. f
%   should have only one root and f(UB) and f(LB) must bound it. Elements
%   of x are NaN for instances where a solution could not be found.
%
%   x = BISECTION(f,LB,UB,target) finds x such that f(x) = target.
%
%   x = BISECTION(f,LB,UB,target,TolX) will terminate the search when the
%   search interval is smaller than TolX (TolX must be positive).
%
%   x = BISECTION(f,LB,UB,target,options) solves with the default
%   parameters replaced by values in the structure OPTIONS, an argument
%   created with the OPTIMSET function. Used options are TolX and TolFun.
%   Note that OPTIMSET will not allow arrays for tolerances, so set the
%   fields of the options structure manually for non-scalar TolX or TolFun.
%
%   [x,fVal] = BISECTION(f,...) returns the value of f evaluated at x.
%
%   [x,fVal,ExitFlag] = BISECTION(...) returns an ExitFlag that describes
%   the exit condition of BISECTION. Possible values of elements of
%   ExitFlag and the corresponding exit conditions are
%
%       1   Search interval smaller than TolX.
%       2   Function value within TolFun of target.
%       3   Search interval smaller than TolX AND function value within
%           TolFun of target.
%      -1   No solution found.
%
%   Any or all of f(scalar), f(array), LB, UB, target, TolX, or TolFun may
%   be scalar or n-dim arrays. All non-scalar arrays must be the same size.
%   All outputs will be this size.
%
%   Default values are target = 0, TolX = 1e-6, and TolFun = 0.
%
%   There is no iteration limit. This is because BISECTION (with a TolX
%   that won't introduce numerical issues) is guaranteed to converge if f
%   is a continuous function on the interval [UB, LB] and f(x)-target
%   changes sign on the interval.
%
%   The bisection method is very robust root-finding method. The absolute
%   error is halved at each step so the method converges linearly. However,
%   Brent's method (such as implemented in FZERO) can converge
%   superlinearly and is as robust. FZERO also has more features and input
%   checking, so use BISECTION in cases where either the optimization
%   toolbox is unavailable or if FZERO would have to be implemented in a
%   loop to solve multiple cases, in which case BISECTION will be much
%   faster because of vectorization.
%
%   Define LB, UB, target, TolX, and TolFun for each specific application
%   using great care for the following reasons:
%     - There is no iteration limit, so given an unsolvable task, such as
%       TolX = TolFun = 0, BISECTION remains in an unending loop.
%     - Spacing between very large floating point numbers is likely to be
%       greater than TolX.
%     - There is no initial check to make sure that f(x) - target changes
%       sign between LB and UB.
%     - Very large or very small numbers can introduce numerical issues.
%
%   Example 1: Find cube root of array 'target' without using NTHROOT and
%   compare speed to using FZERO.
%       options = optimset('TolX', 1e-9);
%       target = [(-100:.1:100)' (-1000:1:1000)'];
%
%       tic;
%       xfz = zeros(size(target));
%       for ii = 1:numel(target)
%           xfz(ii) = fzero(@(x) x.^3-target(ii), [-20 20], options);
%       end
%       fzero_time = toc
%
%       tic;
%       xbis = bisection(@(x) x.^3, -20, 20, target, options);
%       bisection_time = toc
%
%       fprintf('FZERO took %0.0f times longer than BISECTION.\n',...
%                   fzero_time/bisection_time)
%
%   Example 2: Find roots by varying the function coefficients.
%       [A, B] = meshgrid(linspace(1,2,6), linspace(4,12,10));
%       f = @(x) A.*x.^0.2 + B.*x.^0.87 - 15;
%       xstar = bisection(f,0,5);
%
%   See also FZERO, FMINBND, OPTIMSET, FUNCTION_HANDLE.
%
%   [x,fVal,ExitFlag] = BISECTION(f,LB,UB,target,options)

%   Copyright 2010-2013 Sky Sartorius
%   Author  - Sky Sartorius
%   Contact - www.mathworks.com/matlabcentral/fileexchange/authors/101715

% --- Process inputs. ---
% Set default values
tolX    = 1e-6;
tolFun  = 0;
if nargin == 5
    if isstruct(options)
        if isfield(options,'TolX') && ~isempty(options.TolX)
            tolX = options.TolX;
        end
        if isfield(options,'TolFun') && ~isempty(options.TolFun)
            tolFun = options.TolFun;
        end
    else
        tolX = options;
    end
end
if nargin<4 || isempty(target); target=0; end

ub_in = ub; lb_in = lb;
f = @(x) f(x) - target;

% --- Flip UB and LB if necessary. ---
isFlipped = lb>ub;
if any(isFlipped(:))
    ub(isFlipped) = lb_in(isFlipped);
    lb(isFlipped) = ub_in(isFlipped);
    ub_in = ub; lb_in = lb;
end

% --- Make sure everything is the same size for a non-scalar problem. ---
if isscalar(lb) && isscalar(ub)
    % Test if f returns multiple outputs for scalar input.
    if ~isscalar(target)
        ub = ub + zeros(size(target));
    else
        jnk = f(ub);
        if ~isscalar(jnk)
            ub = ub + zeros(size(jnk));
        end
    end
end

% Check if lb and/or ub need to be made into arrays.
if isscalar(lb) && ~isscalar(ub)
    lb = lb + zeros(size(ub));
elseif ~isscalar(lb) && isscalar(ub)
    ub = ub + zeros(size(lb));
end

x = zeros(size(lb));
testconvergence
% --- Iterate ---
iterations = 0;
while any(stillNotDone(:))
    bigger  = fx.*f(ub) > 0;
    ub(bigger)= x(bigger);
    lb(~bigger)= x(~bigger);

    testconvergence;
    iterations = iterations + 1;

    if(iterations > 10000)
%         error('Too many iterations in bisection!');
        break;
    end
end

    function testconvergence
        x=(ub+lb)/2;
        fx=f(x);
        outsideTolFun =  abs(fx)  > tolFun;
        outsideTolX =   (ub - lb) > tolX;
        stillNotDone = outsideTolX & outsideTolFun;
    end

% --- Check that f(x+tolX) and f(x-tolX) have opposite sign. ---
fu = f(min(x+tolX,ub_in));
fl = f(max(x-tolX,lb_in));
unboundedRoot = (fu.*fl) > 0;

% Throw out unbounded results if not meeting TolFun convergence criteria.
x(unboundedRoot & outsideTolFun) = NaN;

% --- Catch NaN elements of UB, LB, target, or other funky stuff. ---
x(isnan(fx)) = NaN;

% --- Characterize results. ---
fx = fx + target;
if nargout > 2
    exitFlag                                    = +~outsideTolX;
    exitFlag(~outsideTolFun)                    =  2;
    exitFlag(~outsideTolFun & ~outsideTolX)     =  3;
    exitFlag(isnan(x))                          = -1;
end

end

% V2: July     2010
% V3: December 2012
% don't remember when
%   typo line 39; added fn handle to see also; made array in example 2
%   smaller; changed wording in example 1
% 2013-08-23
%  -changed scalar*ones(...) calls to scalar+zeros(...) calls based on
%   http://undocumentedmatlab.com/blog/allocation-performance-take-2/
%  -rearranged help block and formatted a tiny bit

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function transducer = makeTransducer(kgrid, transducer_properties)
%MAKETRANSDUCER Create k-Wave ultrasound transducer.
%
% DESCRIPTION:
%       makeTransducer creates an object of the kWaveTransducer class which
%       can be substituted for the source or sensor inputs when using
%       kspaceFirstOrder3D.
%
%       Note: This function will not work with older versions of MATLAB in
%       which custom class definitions are not supported.
%
% USAGE:
%       transducer = makeTransducer(kgrid, settings)
%
% INPUTS:
%
%       kgrid           - k-Wave grid structure returned by makeGrid
%       settings        - input structure used to define the properties of
%                         the transducer (see below)
%
%       The following parameters can be appended as fields to the settings
%       input structure. These parameters are fixed when the transducer is
%       initialised (they cannot be changed without re-creating the
%       transducer). All the settings are optional and are given their
%       default values if not defined.
%
%       number_elements - total number of transducer elements (default =
%                         128)
%       element_width   - width of each element in grid points
%                         (default = 1)
%       element_length  - length of each element in grid points
%                         (default = 10)
%       element_spacing - spacing (kerf width) between the transducer
%                         elements in grid points (default = 0)
%       position        - position of the corner of the transducer within
%                         the grid in grid points (default = [1, 1, 1])
%       radius          - radius of curvature of the transducer [m]
%                         (currently only inf is supported) (default = inf)
%       input_signal    - signal used to drive the ultrasound transducer
%                         (where the sampling rate is defined by kgrid.dt
%                         or kgrid.t_array) (default = [])
%
%       The following parameters can also be appended as fields to the
%       settings input structure, and can additionally be modified after
%       the transducer has been initialised.
%
%       active_elements - transducer elements that are currently active
%                         elements (default = all elements)
%       beamforming_delay_offset
%                       - beamforming delay offset (used to force
%                         beamforming delays to be positive) (default =
%                         'auto')
%       elevation_focus_distance
%                       - focus depth in the elevation direction [m]
%                         (default = inf)
%       focus_distance  - focus distance used to calculate beamforming
%                         delays [m] (default = inf)
%       receive_apodization
%                       - receive apodization; can be set to any of the
%                         window shapes supported by getWin, or given as a
%                         vector the same length as number_elements
%                         (default = 'Rectangular')
%       sound_speed     - sound speed used to calculate beamforming delays
%                         [m/s] (default = 1540)
%       steering_angle  - steering angle used to calculate beamforming
%                         delays [deg] (default = 0)
%       transmit_apodization
%                       - transmit apodization; can be set to any of the
%                         window shapes supported by getWin, or given as a
%                         vector the same length as number_elements
%                         (default = 'Rectangular')
%
% OUTPUTS:
%
%       transducer      - kWaveTransducer object which can be used to
%                         replace the source or sensor inputs of
%                         kspaceFirstOrder3D
%
%       In addition to the input parameters given above (which are also
%       accessible after the transducer has been created), the
%       kWaveTransducer object has a number dependent properties and
%       methods.
%
%       input_signal    - user defined input signal appended and prepended
%                         with additional zeros depending on the values of
%                         focus_distance, elevation_focus_distance, and
%                         steering_angle
%       mask            - binary mask of the active transducer elements
%       number_active_elements
%                       - current number of active transducer elements
%       transducer_width
%                       - total width of the transducer in grid points
%
%       active_elements_mask
%                       - return a binary mask of the active transducer
%                         elements (identical to mask)
%       all_elements_mask
%                       - return a binary mask of all the transducer
%                         elements (both active and inactive)
%       beamforming_delays
%                       - return a vector of the beam forming delays (in
%                         units of time samples) for each active element
%                         based on the focus and steering angle settings
%       combine_sensor_data(sensor_data)
%                       - combine the sensor data returned by
%                         kspaceFirstOrder3D-OMP and
%                         kspaceFirstOrder3D-CUDA to give a single time
%                         series per active transducer element (anlogous to
%                         the output from the MATLAB code), rather than a
%                         time series per grid point
%       delay_mask      - return a mask of the active transducer elements,
%                         where the mask values contain the beamforming
%                         delays (an integer input can also be given to
%                         control the beamforming delays used, where 1:
%                         both delays, 2: elevation only, 3: azimuth only)
%       elevation_beamforming_delays - return a vector of the elevation
%                         beam forming delays (in units of time samples)
%                         for each active element based on the elevation
%                         focus setting
%       get_receive_apodization
%                       - return the receive apodization
%       get_transmit_apodization
%                       - return the transmit apodization
%       indexed_active_elements_mask
%                       - return a mask of the active transducer elements,
%                         where the mask values indicate which transducer
%                         element each grid point corresponds to
%       indexed_elements_mask
%                       - return a mask of all the transducer elements
%                         (both active and inactive), where the mask values
%                         indicate which transducer element each grid point
%                         corresponds to
%       plot            - plot the transducer using voxelPlot
%       properties      - print a list of the transducer properties to the
%                         command line
%       scan_line(sensor_data)
%                       - combine the input sensor data using the current
%                         apodization and beamforming setting to generate a
%                         single scan line
%       transmit_apodization_mask
%                       - return a mask of the active transducer elements,
%                         where the mask values contain the apodization
%                         weights
%
% ABOUT:
%       author          - Bradley Treeby
%       date            - 28th July 2011
%       last update     - 4th November 2013
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox
%
% See also makeGrid, kspaceFirstOrder3D, kWaveTransducer

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

try
    % create a new instance of the kWaveTransducer class
    transducer = kWaveTransducer(kgrid, transducer_properties);
catch ME
    % if user defined classes aren't supported, throw an error
    if strcmp(ME.identifier, 'MATLAB:UndefinedFunction')
        error('The transducer cannot be created because user defined classes are not supported in your version of MATLAB. To use this functionality, please try using a newer MATLAB version.');
    end
    rethrow(ME);
end

repo: Aero-matlab