text:
function sensor_data = kspaceSecondOrder(kgrid, medium, source, sensor, varargin)
%KSPACESECONDORDER  Fast time-domain simulation of wave propagation for homogeneous media.
%
% DESCRIPTION:
%       kspaceSecondOrder simulates the time-domain propagation of linear
%       compressional waves through a one, two, or three dimensional
%       homogeneous acoustic medium given four input structures: kgrid,
%       medium, source, and sensor. The computation is based on an exact
%       second-order k-space model for media with power law absorption. At
%       each time-step (defined by kgrid.t_array), the pressure at the
%       positions defined by sensor.mask are recorded and stored. If
%       kgrid.t_array is set to 'auto', this array is automatically
%       generated using makeTime. To prevent wave wrapping, the
%       computational domain can be automatically expanded by a factor of
%       two by setting the optional input 'ExpandGrid' to true.
%
%       An initial pressure distribution can be specified by assigning a
%       matrix (the same size as the computational grid) of arbitrary
%       numeric values to source.p0. An initial pressure gradient can
%       similarly be specified using source.dp0dt. The pressure is returned
%       as an array of time series at the sensor locations defined by
%       sensor.mask. This is specified as a binary matrix (i.e., a matrix
%       of 1's and 0's the same size as the computational grid)
%       representing the grid points within the computational grid that
%       will collect the data. The sensor_data is returned using MATLAB's
%       standard column-wise linear matrix index ordering with the recorded
%       data indexed as sensor_data(sensor_position, time). The final
%       pressure field over the complete computational grid can also be
%       obtained by setting sensor.record to {'p', 'p_final'}. In this
%       case, the output sensor_data is returned as a structure with the
%       outputs appended as the structure fields sensor_data.p and
%       sensor_data.p_final.
%
%       Compared to the first-order simulation functions
%       kspaceFirstOrder1D, kspaceFirstOrder2D, and kspaceFirstOrder3D,
%       kspaceSecondOrder is restricted to homogeneous media and has less
%       functionality. However, it is also more computationally efficient
%       and allows an initial pressure gradient to be specified.
%
% USAGE:
%       sensor_data = kspaceSecondOrder(kgrid, medium, source, sensor)
%       sensor_data = kspaceSecondOrder(kgrid, medium, source, sensor, ...)
%
% INPUTS:
%       kgrid               - k-Wave grid structure returned by makeGrid
%                             containing Cartesian and k-space grid fields
%       kgrid.t_array       - evenly spaced array of time values [s] (set
%                             to 'auto' by makeGrid)
%
%       medium.sound_speed  - homogeneous sound speed within the acoustic
%                             medium [m/s]
%       medium.alpha_power  - power law absorption exponent
%       medium.alpha_coeff  - power law absorption coefficient
%                             [dB/(MHz^y cm)]
%       medium.alpha_mode   - optional input to force either the absorption
%                             or dispersion terms in the equation of state
%                             to be excluded; valid inputs are
%                             'no_absorption' or 'no_dispersion'
%
%       source.p0           - initial pressure within the acoustic medium
%       source.dp0dt        - initial pressure gradient within the acoustic
%                             medium
%
%       sensor.mask         - binary grid specifying where the pressure is
%                             recorded at each time-step
%       sensor.record       - cell array of the acoustic parameters to
%                             record in the form sensor.record = {'p'};
%                             valid inputs are:
%                               'p'       (acoustic pressure)
%                               'p_final' (final pressure field)
%
% OPTIONAL INPUTS:
%       Optional 'string', value pairs that may be used to modify the
%       default computational settings.
%
%       'ExpandGrid'    - Boolean controlling whether the grid size is
%                         expanded on two sides to delay the time before
%                         wave wrapping occurs (default = false).
%       'MeshPlot'      - Boolean controlling whether mesh is used in place
%                         of imagesc to plot the pressure field (default =
%                         false).
%       'PlotFrames'    - Boolean controlling whether the pressure field
%                         for each time step is plotted in a new window
%                         (default = false).
%       'PlotFreq'      - The number of iterations which must pass before
%                         the simulation plot is updated (default = 10).
%       'PlotScale'     - [min, max] values used to control the scaling for
%                         imagesc (visualisation) (default = [-1, 1].
%       'PlotSim'       - Boolean controlling whether the simulation
%                         iterations are progressively plotted (default =
%                         true).
%       'Smooth'        - Boolean controlling whether source.p0 is smoothed
%                         using smooth before computation (default = true).
%
% OUTPUTS:
% If sensor.record is not defined by the user:
%       sensor_data         - time varying pressure recorded at the sensor
%                             positions given by sensor.mask
%
% If sensor.record is defined by the user:
%       sensor_data.p       - time varying pressure recorded at the sensor
%                             positions given by sensor.mask (returned if
%                             'p' is set)
%       sensor_data.p_final - final pressure field over the complete domain
%                             (returned if 'p_final' is set)
%
% ABOUT:
%       author      - Bradley Treeby and Ben Cox
%       date        - 21st August 2008
%       last update - 25th August 2014
%
% 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 kspaceFirstOrder1D, kspaceFirstOrder2D, kspaceFirstOrder3D,
% makeGrid, makeTime, smooth

% 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/>.

% start the timer
tic;

% =========================================================================
% DEFINE LITERALS
% =========================================================================

% minimum number of input variables
NUM_REQ_INPUT_VARIABLES = 4;

% optional input defaults
EXPAND_GRID_DEF = false;
MESH_PLOT_DEF = false;
PLOT_FRAMES_DEF = false;
PLOT_FREQ_DEF = 10;
PLOT_SCALE_DEF = [-1 1];
PLOT_SIM_DEF = true;
SMOOTH_P0_DEF = true;

% create colormap for visualisations
COLOR_MAP = getColorMap;

% =========================================================================
% CHECK INPUT STRUCTURES
% =========================================================================

% calculate t_array using makeTime if it is not given
if strcmp(kgrid.t_array, 'auto')
    kgrid.t_array = makeTime(kgrid, medium.sound_speed);
end

% check medium fields
checkFieldNames(medium, {'sound_speed', 'density', 'alpha_coeff', 'alpha_power', 'alpha_mode'});
enforceFields(medium, {'sound_speed'});
if isfield(medium, 'density')
    disp('WARNING: medium.density is not a valid input parameter');
end

% check source fields
checkFieldNames(source, {'p0', 'dp0dt'});

% check the absorption mode input is valid
if isfield(medium, 'alpha_mode')
    if ~ischar(medium.alpha_mode) || (~strcmp(medium.alpha_mode, 'no_absorption') && ~strcmp(medium.alpha_mode, 'no_dispersion'))
        error('medium.alpha_mode must be set to ''no_absorption'' or ''no_dispersion''');
    end
end

% extract which source fields have been given
source_case = 0;
if isfield(source, 'p0')
    source_case = source_case + 1;
end
if isfield(source, 'dp0dt')
    source_case = source_case + 2;
end

% check sensor fields
checkFieldNames(sensor, {'mask', 'record'});
enforceFields(sensor, {'mask'});

% check the sensor mask is binary
if sum(sensor.mask(:)) ~= numel(sensor.mask) - sum(sensor.mask(:) == 0)
    error('sensor.mask must be a binary grid (numeric values must be 0 or 1)');
end

% set default record options
record.p = true;
record.p_final = false;

% check for sensor.record and set usage flags - if no flags are given, the
% time history of the acoustic pressure is recorded by default
if isfield(sensor, 'record')

    % list of allowable record flags
    record.flags = {'p', 'p_final'};

    % check the contents of the cell array are valid inputs
    for record_index = 1:length(sensor.record)
        if ~ismember(sensor.record(record_index), record.flags);
            error(['''' sensor.record{record_index} ''' is not a valid input for sensor.record']);
        end
    end

    % set the usage flags depending on the cell array
    if ismember('p', sensor.record)
        record.p = true;
    else
        % set record.p to false if a user input for sensor.record is given
        % and 'p' is not set
        record.p = false;
    end
    if ismember('p_final', sensor.record)
        record.p_final = true;
    end

end

% =========================================================================
% CHECK OPTIONAL INPUTS
% =========================================================================

% shorten commonly used field names
t_array = kgrid.t_array;
c = medium.sound_speed;
if isfield(medium, 'alpha_coeff')
    a0 = medium.alpha_coeff;
    y = medium.alpha_power;
else
    a0 = 0;
end

% assign default input parameters
expand_grid = EXPAND_GRID_DEF;
mesh_plot = MESH_PLOT_DEF;
plot_frames = PLOT_FRAMES_DEF;
plot_freq = PLOT_FREQ_DEF;
plot_scale = PLOT_SCALE_DEF;
plot_sim = PLOT_SIM_DEF;
smooth_p0 = SMOOTH_P0_DEF;

% replace defaults with user defined values if provided and check inputs
if nargin < NUM_REQ_INPUT_VARIABLES
    error('Not enough input parameters');
elseif rem(nargin - NUM_REQ_INPUT_VARIABLES, 2) %**
    error('Optional input parameters must be given as param, value pairs');
elseif ~isempty(varargin)
    for input_index = 1:2:length(varargin)
        switch varargin{input_index}
            case 'ExpandGrid'
                expand_grid = varargin{input_index + 1};
                if ~islogical(expand_grid)
                    error('Optional input ExpandGrid must be Boolean');
                end
            case 'MeshPlot'
                mesh_plot = varargin{input_index + 1};
                if ~islogical(mesh_plot)
                    error('Optional input MeshPlot must be Boolean');
                end
            case 'PlotFrames'
                plot_frames = varargin{input_index + 1};
                if ~islogical(plot_frames)
                    error('Optional input PlotFrames must be Boolean');
                end
            case 'PlotFreq'
                plot_freq = varargin{input_index + 1};
                if ~(numel(plot_freq) == 1 && isnumeric(plot_freq))
                    error('Optional input PlotFreq must be a single numerical value');
                end
            case 'PlotScale'
                plot_scale = varargin{input_index + 1};
                if ~(numel(plot_scale) == 2 && isnumeric(plot_scale))
                    error('Optional input PlotScale must be a 2 element numerical array');
                end
            case 'PlotSim'
                plot_sim = varargin{input_index + 1};
                if ~islogical(plot_sim)
                    error('Optional input PlotSim must be Boolean');
                end
            case 'Smooth'
                smooth_p0 = varargin{input_index + 1};
                if ~islogical(smooth_p0)
                    error('Optional input Smooth must be Boolean');
                end
            otherwise
                error('Unknown optional input');
        end
    end
end

% cleanup unused variables
clear *_DEF NUM_REQ_INPUT_VARIABLES time;

% update command line status
disp('Running k-space simulation...');

% =========================================================================
% UPDATE COMMAND LINE STATUS
% =========================================================================

% update command line status
disp(['  time steps: ' num2str(length(t_array))]);
switch numDim(kgrid.k)
    case 1
        disp(['  input grid size: ' num2str(kgrid.Nx) ' grid points (' scaleSI(kgrid.x_size) 'm)']);
    case 2
        [x_sc, scale, prefix] = scaleSI(min(kgrid.x_size, kgrid.y_size)); %#ok<*ASGLU>
        disp(['  input grid size: ' num2str(kgrid.Nx) ' by ' num2str(kgrid.Ny) ' grid points (' num2str(kgrid.x_size*scale) ' by ' num2str(kgrid.y_size*scale) prefix 'm)']);
    case 3
        [x_sc, scale, prefix] = scaleSI(min([kgrid.x_size, kgrid.y_size, kgrid.z_size])); %#ok<ASGLU>
        disp(['  input grid size: ' num2str(kgrid.Nx) ' by ' num2str(kgrid.Ny) ' by ' num2str(kgrid.Nz) ' grid points (' num2str(kgrid.x_size*scale) ' by ' num2str(kgrid.y_size*scale) ' by ' num2str(kgrid.z_size*scale) prefix 'm)']);
end
disp(['  maximum supported frequency: ' scaleSI( kgrid.k_max * min(c(:)) / (2*pi) ) 'Hz']);

% =========================================================================
% CALCULATE FFT OF SOURCE FUNCTIONS AND EXPAND GRID
% =========================================================================

% smooth p0 distribution if required restoring the maximum magnitude
if smooth_p0
    if source_case ~= 2
        disp('  smoothing source.p0 distribution...');
        source.p0 = smooth(kgrid, source.p0, true);
    end
    if source_case > 1
        disp('  smoothing source.dp0dt distribution...');
        source.dp0dt = smooth(kgrid, source.dp0dt, true);
    end
end

if ~expand_grid
    % extract wavenumber matrix from kgrid
    k = kgrid.k;

    % extract the location of the sensor points
    sensor_index = (sensor.mask == 1);

    % compute FFT of source functions
    if source_case ~= 2
        p0_k = fftn(source.p0);
    end
    if source_case > 1
        dp0dt_k = fftn(source.dp0dt);
    end
else
    switch numDim(kgrid.k)
        case 1

            % create a larger kgrid to prevent wave wrapping
            kgrid_APE = makeGrid(kgrid.Nx*2, kgrid.dx);

            % extract wavenumber matrix from kgrid
            k = kgrid_APE.k;

            % expand the sensor mask
            sensor_mask = [sensor.mask, zeros(size(sensor.mask))];

            % expand the source functions
            if source_case ~= 2
                p0 = zeros(size(k));
                p0(1:kgrid.Nx) = source.p0;
            end
            if source_case > 1
                dp0dt = zeros(size(k));
                dp0dt(1:kgrid.Nx) = source.dp0dt;
            end

        case 2

            % create a larger kgrid to prevent wave wrapping
            kgrid_APE = makeGrid(kgrid.Nx*2, kgrid.dx, kgrid.Ny*2, kgrid.dy);

            % extract wavenumber matrix from kgrid
            k = kgrid_APE.k;

            % expand the sensor mask
            sensor_mask = zeros(size(k));
            sensor_mask(1:kgrid.Nx, 1:kgrid.Ny) = sensor.mask;

            % expand the source functions
            if source_case ~= 2
                p0 = zeros(size(k));
                p0(1:kgrid.Nx, 1:kgrid.Ny) = source.p0;
            end
            if source_case > 1
                dp0dt = zeros(size(k));
                dp0dt(1:kgrid.Nx, 1:kgrid.Ny) = source.dp0dt;
            end

        case 3

            % create a larger kgrid to prevent wave wrapping
            kgrid_APE = makeGrid(kgrid.Nx*2, kgrid.dx, kgrid.Ny*2, kgrid.dy, kgrid.Nz*2, kgrid.dz);

            % extract wavenumber matrix from kgrid
            k = kgrid_APE.k;

            % extract the location of the sensor points within the enlarged grid
            sensor_mask = zeros(size(k));
            sensor_mask(1:kgrid.Nx, 1:kgrid.Ny, 1:kgrid.Nz) = sensor.mask;

            % expand the source functions
            if source_case ~= 2
                p0 = zeros(size(k));
                p0(1:kgrid.Nx, 1:kgrid.Ny, 1:kgrid.Nz) = source.p0;
            end
            if source_case > 1
                dp0dt = zeros(size(k));
                dp0dt(1:kgrid.Nx, 1:kgrid.Ny, 1:kgrid.Nz) = source.dp0dt;
            end

    end

    % extract the location of the sensor points within the enlarged grid
    sensor_index = (sensor_mask == 1);

    % compute FFT of source functions
    if source_case ~= 2
        p0_k = fftn(p0);
    end
    if source_case > 1
        dp0dt_k = fftn(dp0dt);
    end

    % delete unused variables
    clear kgrid_APE sensor_mask p0 dp0dt source;
end

% =========================================================================
% DEFINE TIME PROPAGATORS
% =========================================================================

% shift the wavenumbers
k = ifftshift(k);

% find the index of the zero wave numbers
if source_case > 1
    k_0_index = find(k == 0);
end

% define time propagation function
if a0 ~= 0

    % convert attenuation to nepers
    a0 = db2neper(a0, y);

    % define coefficients
    if y == 2 || (isfield(medium, 'alpha_mode') && strcmp('no_dispersion', medium.alpha_mode))

        % compute Upsilon
        Upsilon = sqrt(1 - a0^2*c^(2*y)*k.^(2*y-2));

        % reset first value to avoid possible INF when k = 0
        Upsilon(1) = 1;

        % precalculate coefficients to speed up computation
        ckU = c*k.*Upsilon;
        coeff_1 = a0*c^(y + 1)*k.^(y);
        coeff_2 = 0;

    elseif (isfield(medium, 'alpha_mode') && strcmp('no_absorption', medium.alpha_mode))

        % compute Upsilon
        Upsilon = sqrt(1 - 2*a0*c^y*k.^(y-1)*tan(pi*y/2));

        % reset first value to avoid possible INF when k = 0
        Upsilon(1) = 1;

        % precalculate coefficients to speed up computation
        ckU = c*k.*Upsilon;
        coeff_1 = 0;
        coeff_2 = a0*c^y*k.^(y-1) ./ Upsilon;

        % reset first value to avoid possible INF when k = 0
        coeff_2(1) = 0;

    else

        % compute Upsilon
        Upsilon = sqrt(1 - a0^2*c^(2*y)*k.^(2*y-2) - 2*a0*c^y*k.^(y-1)*tan(pi*y/2));

        % reset first value to avoid possible INF when k = 0
        Upsilon(1) = 1;

        % precalculate coefficients to speed up computation
        ckU = c*k.*Upsilon;
        coeff_1 = a0*c^(y + 1)*k.^(y);
        coeff_2 = a0*c^y*k.^(y-1) ./ Upsilon;

        % reset first value to avoid possible INF when k = 0
        coeff_2(1) = 0;

    end

else
    % precalculate coefficients to speed up computation
    ck = c*k;
end

% =========================================================================
% PREPARE VISUALISATIONS AND STORAGE VARIABLES
% =========================================================================

% preallocate storage variables
if record.p
    sensor_data.p = zeros(sum(sensor.mask(:)), length(t_array));
end

% pre-compute suitable axes scaling factor
if plot_sim
    switch numDim(kgrid.k)
        case 1
            [x_sc, scale, prefix] = scaleSI(max(kgrid.x_vec));
        case 2
            [x_sc, scale, prefix] = scaleSI(max([kgrid.x_vec; kgrid.y_vec]));
        case 3
            [x_sc, scale, prefix] = scaleSI(max([kgrid.x_vec; kgrid.y_vec; kgrid.z_vec]));
    end
end

% initialise the figures used for animation if 'PlotSim' is set to 'true'
if plot_sim
    % create empty figure
    img = figure;

    % create waitbar and shift position so it doesn't overlap the figure window
    pbar = waitbar(0, 'Computing Pressure Field', 'Visible', 'off');
    posn_pbar = get(pbar, 'OuterPosition');
    posn_img = get(img, 'OuterPosition');
    posn_pbar(2) = max(min(posn_pbar(2) - posn_pbar(4), posn_img(2) - posn_pbar(4) - 10), 0);
    set(pbar, 'OuterPosition', posn_pbar, 'Visible', 'on');
end

% update command line status
disp(['  precomputation completed in ' scaleTime(toc)]);
disp('  starting time loop...');

% =========================================================================
% LOOP THROUGH TIME STEPS
% =========================================================================

% loop through each value of t
for t_index = 1:length(t_array)

    % extract the time point
    t = t_array(t_index);

    % compute pressure field
    switch (source_case + 3*(a0 ~= 0))
        case 1
            TP_p0 = cos(ck*t);
            p = ifftn(  p0_k.*TP_p0  );
        case 2
            TP_dp0dt = sin(ck*t) ./ (ck);
            TP_dp0dt(k_0_index) = t; %#ok<*FNDSB>
            p = ifftn(  dp0dt_k.*TP_dp0dt  );
        case 3
            TP_p0 = cos(ck*t);
            TP_dp0dt = sin(ck*t) ./ (ck);
            TP_dp0dt(k_0_index) = t;
            p = ifftn(  p0_k.*TP_p0 + dp0dt_k.*TP_dp0dt  );
        case 4
            TP_p0 = exp(-coeff_1*t) .* ( cos(ckU*t) - coeff_2 .* sin(ckU*t) );
            p = ifftn(  p0_k.*TP_p0  );
        case 5
            TP_dp0dt = sin(ckU*t) .* exp(-coeff_1*t) ./ ckU;
            TP_dp0dt(k_0_index) = t.* exp(-coeff_1(k_0_index)*t);
            p = ifftn(  dp0dt_k.*TP_dp0dt  );
        case 6
            TP_p0 = exp(-coeff_1*t) .* ( cos(ckU*t) - coeff_2 .* sin(ckU*t) );
            TP_dp0dt = sin(ckU*t) .* exp(-coeff_1*t) ./ ckU;
            TP_dp0dt(k_0_index) = t.* exp(-coeff_1(k_0_index)*t);
            p = ifftn(  p0_k.*TP_p0 + dp0dt_k.*TP_dp0dt  );
    end

    % extract required data
    if record.p
        sensor_data.p(:, t_index) = real(p(sensor_index));
    end

    % plot data if required
    if plot_frames || ( plot_sim && rem(t_index, plot_freq) == 0 )

        % update progress bar
        waitbar(t_index/length(t_array), pbar);
        drawnow;

        % update plot
        if plot_frames
            figure;
        end

        p = real(p);

        switch numDim(kgrid.k)
            case 1
                % extract the required pressure field from the enlarged grid
                if expand_grid
                    p = p(1:kgrid.Nx);
                end
                plot(kgrid.x*scale, p);
                xlabel(['x-position [' prefix 'm]']);
                set(gca, 'YLim', plot_scale);

            case 2
                % extract the required pressure field from the enlarged grid
                if expand_grid
                    p = p(1:kgrid.Nx, 1:kgrid.Ny);
                end

                if mesh_plot
                    mesh(kgrid.x_vec*scale, kgrid.y_vec*scale, p, 'EdgeColor', 'Black');
                    axis image;
                    set(gca, 'ZLim', plot_scale);
                    axis off;
                else
                    % add sensor mask onto plot
                    p(sensor.mask == 1) = plot_scale(2);

                    imagesc(kgrid.y_vec*scale, kgrid.x_vec*scale, p, plot_scale);
                    colormap(COLOR_MAP);
                    axis image;
                    ylabel(['x-position [' prefix 'm]']);
                    xlabel(['y-position [' prefix 'm]']);
                end

            case 3
                % extract the required pressure field from the enlarged grid
                if expand_grid
                    p = p(1:kgrid.Nx, 1:kgrid.Ny, 1:kgrid.Nz);
                end
                planeplot(kgrid, p, '', plot_scale, scale, prefix, COLOR_MAP);
        end
    end
end

% =========================================================================
% CLEAN UP
% =========================================================================

% assign output data
if ~isfield(sensor, 'record')
    sensor_data = sensor_data.p;
elseif record.p_final
    switch numDim(kgrid.k)
        case 1
            sensor_data.p_final = p(1:kgrid.Nx);
        case 2
            sensor_data.p_final = p(1:kgrid.Nx, 1:kgrid.Ny);
        case 3
            sensor_data.p_final = p(1:kgrid.Nx, 1:kgrid.Ny, 1:kgrid.Nz);
    end
end

% clean up used figures
if plot_sim
    close(img);
    close(pbar);
end

% update command line status
disp(['  computation completed in ' scaleTime(toc)]);

function planeplot(kgrid, data, data_title, plot_scale, scale, prefix, color_map)
% Subfunction to produce a plot of a three-dimensional matrix through the
% three central planes

subplot(2, 2, 1), imagesc(kgrid.y_vec*scale, kgrid.x_vec*scale, squeeze(data(:, :, kgrid.Nz/2)), plot_scale);
title([data_title 'x-y plane']);
axis image;
subplot(2, 2, 2), imagesc(kgrid.z_vec*scale, kgrid.x_vec*scale, squeeze(data(:, kgrid.Ny/2, :)), plot_scale);
title('x-z plane');
axis image;
xlabel(['(All axes in ' prefix 'm)']);
subplot(2, 2, 3), imagesc(kgrid.z_vec*scale, kgrid.y_vec*scale, squeeze(data(kgrid.Nx/2, :, :)), plot_scale);
title('y-z plane');
axis image;
colormap(color_map);
drawnow;


text:
function [params] = sv_calcMc(params)
    % Calculation of magnitude of completeness
% [params] = sv_calcMc(params)
% -------------------------------------
% Calculation of magnitude of completeness
%
% Input parameters:
%   params.mCatalog           Earthquake catalog
%   params.mPolygon           Polygon (defined by ex_selectgrid)
%   params.vX                 X-vector (defined by ex_selectgrid)
%   params.vY                 Y-vector (defined by ex_selectgrid)
%   params.vUsedNodes         Used nodes vX * vY defining the mPolygon (defined by ex_selectgrid)
%   params.bRandom            Perform random simulation (true) or real calculation (false)
%   params.nCalculation       Number of random simulations
%   params.bMap               Calculate a map (true) or a cross-section (false)
%   params.bNumber            Use constant number (true) or constant radius (false)
%   params.nNumberEvents      Number of earthquakes if bNumber is true
%   params.fMaxRadius         Maximum Radius using a constant number of events; works only when bNumber is true
%   params.fRadius            Radius of gridnode if bNumber is false
%   params.nMinimumNumber     Minimum number of earthquakes per node for determining a b-value
%   params.fMinMag            Lower limit of magnitude range for testing
%   params.fMaxMag            Upper limit of magnitude range for testing
%   params.bTimePeriod        Calculate seismicity difference for 2 periods (0) until start and end of catalog or
%                             a specific time period before and after fSplitTime (1)
%   params.fTimePeriodDays        Length of time periods
%   params.bTstart            Check for starting time of temporal mapping
%   params.fTstart            Starting time for temporal mapping
%   params.bBstnum            Check for boostrap sampling
%   params.fBstnum            Number of bootstrap samples
%   params.fBinning           Bin size for magnitude binning
%   params.sComment           Comment on calculation

% Output parameters:
%   Same as input parameters including
%   params.mValueGrid         Matrix of calculated values
%   params.vcsGridNames       Names of parameters calculated
%   Check sv_NodeCalcMc.m for a list of variables!!
%
% J. Woessner; woessner@seismo.ifg.ethz.ch
% updated: 11.06.03

report_this_filefun();

% Initialize
vResults = [];
params.sComment = [];
if isempty(params.fBinning)
    params.fBinning = 0.1;
end

% Determine time period of catalog
params.fTminCat = min(params.mCatalog.Date);
params.fTmaxCat = max(params.mCatalog.Date);

% Init result matrix
mValueGrid_ = [];

% Temporary saving the original catalog
mCatalog = params.mCatalog;

% Check for bootstrapping or not
% ------------------------------
% Case of calculations with bootstrapping
if params.bBstnum
    % Loop over time
    fTstart = params.fTstart;
    while fTstart < params.fTmaxCat
        mValueGrid_ = [];
        params.mCatalog = mCatalog;
        % Create Indices to catalog and select quakes in time period
        vSel = (fTstart <= params.mCatalog.Date & params.mCatalog.Date < fTstart+fTimePeriodDays);
        params.mCatalog = params.mCatalog.subset(vSel);
        [params.caNodeIndices] = ex_CreateIndexCatalog(params.mCatalog, params.mPolygon, params.bMap, params.nGriddingMode, ...
         params.nNumberEvents, params.fRadius, params.fSizeRectHorizontal, params.fSizeRectDepth);
        % Loop over all grid nodes
        hWaitbar1 = waitbar(0,'Calculating nodes...');
        set(hWaitbar1,'Numbertitle','off','Name','Node percentage');
        for nNode_ = 1:length(params.mPolygon(:,1))
            % Create node catalog
            mNodeCatalog_ = params.mCatalog(params.caNodeIndices{nNode_}, :);
            % Check for constant number of events calculations
            if (params.nGriddingMode == 0)
                [mNodeCatalog_] = ex_CheckMaxRadius(mNodeCatalog_, params.mPolygon, nNode_, params.caNodeIndices, params.fMaxRadius, params.nNumberEvents, params.bMap);
            end
            [nX,nY] = size(mNodeCatalog_);
            if (nX < params.nMinimumNumber)
               mValueGrid_= [mValueGrid_; NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN nX NaN NaN NaN NaN];
            else
                [rCalcNodeResult_] = sv_NodeCalcMc(params,mNodeCatalog_);
                % Store the results
                mValueGrid_= [mValueGrid_; rCalcNodeResult_.fMc_max rCalcNodeResult_.fMc_90 rCalcNodeResult_.fMc_95 rCalcNodeResult_.fMc_com...
                        rCalcNodeResult_.fMc_EMR rCalcNodeResult_.fMc_shi rCalcNodeResult_.fMc_Bst...
                        rCalcNodeResult_.fStd_Mc rCalcNodeResult_.fBvalue_Bst rCalcNodeResult_.fStd_B...
                        rCalcNodeResult_.fAvalue_Bst rCalcNodeResult_.fStd_A nX...
                        rCalcNodeResult_.bH rCalcNodeResult_.fPval rCalcNodeResult_.bH_Bst rCalcNodeResult_.fPval_Bst];
            end; % End of if on length(mNodeCatalog_(:,1)
            if rem(nNode_,500) == 0
                waitbar(nNode_/length(params.mPolygon(:,1)))
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                % Temporary saving
                params.vcsGridNames = cellstr(char('Mc max. curvature' , 'Mc 90% goodness of fit' , 'Mc 95% goodness of fit',...
                    'Mc best combination', 'Mc EMR-method', 'Mc(Shi-b-uncertainty)', 'Mc(Bst-mean)', 'Mc(Bst-2nd-moment)',...
                    'Mc(Bst-b)', 'Mc(b_2nd-moment)','Mc(Bst-a)', 'Mc(a_2nd-moment)','Number of events',...
                    'H(KST)','P(KST)','H(KST_Bst)','P(KST_Bst)'));
                params.mValueGrid = mValueGrid_;
                % Add parameter to params.sComment
                if  params.nGriddingMode == 0;   % Constant number
                    params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                            ' Time period (days) ' num2str(params.fTimePeriodDays) ' d, Constant number: ' num2str(params.nNumberEvents) ', MaxRadius: '...
                            num2str(params.fMaxRadius) ' km'];
                    vResults = params;
                    save(['result_Time' num2str(fTstart) '_Constnum_' num2str(params.nNumberEvents) '_MaxRad_' num2str(params.fMaxRadius)...
                            '_Nmin_' num2str(params.nMinimumNumber) '_Node' num2str(nNode_) '.mat'], 'vResults');
                elseif params.nGriddingMode == 1;   % Constant radius
                    params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                            ' Time period (days) ' num2str(params.fTimePeriodDays) ' d, Radius: ' num2str(params.fRadius) ' km, Nmin: ' num2str(params.nMinimumNumber)];
                    vResults = params;
                    save(['result_Time' num2str(fTstart) '_Rad_' num2str(params.fRadius) '_Nmin_' num2str(params.nMinimumNumber)  '_Node' num2str(nNode_) '.mat'], 'vResults');
                else  % Rectangle mode
                    params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                            ' Time period (days) ' num2str(params.fTimePeriodDays) ' d, Rect. X: ' num2str(params.fSizeRectHorizontal) ' km, Rect. Y: ' num2str(params.fSizeRectDepth)...
                            ' km, Nmin: ' num2str(params.nMinimumNumber)];
                    vResults = params;
                    save(['result_Time' num2str(fTstart) '_RectX_' num2str(params.fSizeRectHorizontal) '_RectY_' num2str(params.fSizeRectDepth)...
                            '_Nmin_' num2str(params.nMinimumNumber)  '_Node' num2str(nNode_) '.mat'], 'vResults');
                end; % END of params.nGriddingmode
                vResults =[];
            end; % End updating waitbar
        end; % for nNode
        close(hWaitbar1);
        % Parameter description
        params.vcsGridNames = cellstr(char('Mc max. curvature' , 'Mc 90% goodness of fit' , 'Mc 95% goodness of fit',...
                    'Mc best combination', 'Mc EMR-method', 'Mc(Shi-b-uncertainty)', 'Mc(Bst-mean)', 'Mc(Bst-2nd-moment)',...
                    'Mc(Bst-b)', 'Mc(b_2nd-moment)','Mc(Bst-a)', 'Mc(a_2nd-moment)','Number of events',...
                    'H(KST)','P(KST)','H(KST_Bst)','P(KST_Bst)'));
        params.mValueGrid = mValueGrid_;
        % Add parameter to params.sComment
        if  params.nGriddingMode == 0;   % Constant number
            params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                    ' Time period (days) ' num2str(params.fTimePeriodDays) ' d, Constant number: ' num2str(params.nNumberEvents) ', MaxRadius: '...
                    num2str(params.fMaxRadius) ' km'];
            vResults = params;
            save(['result_Time' num2str(fTstart) '_Constnum_' num2str(params.nNumberEvents) '_MaxRad_' num2str(params.fMaxRadius)...
                    '_Nmin_' num2str(params.nMinimumNumber) '.mat'], 'vResults');
        elseif params.nGriddingMode == 1;   % Constant radius
            params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                    ' Time period (days) ' num2str(params.fTimePeriodDays) ' d, Radius: ' num2str(params.fRadius) ' km, Nmin: ' num2str(params.nMinimumNumber)];
            vResults = params;
            save(['result_Time' num2str(fTstart) '_Rad_' num2str(params.fRadius) '_Nmin_' num2str(params.nMinimumNumber) '.mat'], 'vResults');
        else  % Rectangle mode
            params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                    ' Time period (days) ' num2str(params.fTimePeriodDays) ' d, Rect. X: ' num2str(params.fSizeRectHorizontal) ' km, Rect. Y: ' num2str(params.fSizeRectDepth)...
                    ' km, Nmin: ' num2str(params.nMinimumNumber)];
            vResults = params;
            save(['result_Time' num2str(fTstart) '_RectX_' num2str(params.fSizeRectHorizontal) '_RectY_' num2str(params.fSizeRectDepth)...
                    '_Nmin_' num2str(params.nMinimumNumber) '.mat'], 'vResults');
        end
        vResults =[];
        fTstart = fTstart+params.fTimePeriodDays;
    end; % End of while fTstart

    % Case of no bootstrapping
else
    % Loop over time
    fTstart = params.fTstart;
    while fTstart < params.fTmaxCat
        mValueGrid_ = [];
        params.mCatalog = mCatalog;
        % Create Indices to catalog and select quakes in time period
        vSel = (fTstart <= params.mCatalog.Date & params.mCatalog.Date < fTstart+params.fTimePeriodDays);
        params.mCatalog = params.mCatalog.subset(vSel);
        [params.caNodeIndices] = ex_CreateIndexCatalog(params.mCatalog, params.mPolygon, params.bMap, params.nGriddingMode, ...
         params.nNumberEvents, params.fRadius, params.fSizeRectHorizontal, params.fSizeRectDepth);
        % Loop over all grid nodes
        hWaitbar1 = waitbar(0,'Calculating nodes...');
        set(hWaitbar1,'Numbertitle','off','Name','Node percentage');
        for nNode_ = 1:length(params.mPolygon(:,1))
            % Create node catalog
            mNodeCatalog_ = params.mCatalog(params.caNodeIndices{nNode_}, :);
            % Check for constant number of events calculations
            if (params.nGriddingMode == 0)
                [mNodeCatalog_] = ex_CheckMaxRadius(mNodeCatalog_, params.mPolygon, nNode_, params.caNodeIndices, params.fMaxRadius, params.nNumberEvents, params.bMap);
            end
            [nX,nY] = size(mNodeCatalog_);
            if (nX < params.nMinimumNumber)
                mValueGrid_= [mValueGrid_; NaN NaN NaN NaN NaN NaN nX];
            else
                [rCalcNodeResult_] = sv_NodeCalcMc(params,mNodeCatalog_);
                mValueGrid_= [mValueGrid_; rCalcNodeResult_.fMc_max rCalcNodeResult_.fMc_90 rCalcNodeResult_.fMc_95 rCalcNodeResult_.fMc_com...
                        rCalcNodeResult_.fMc_EMR rCalcNodeResult_.fMc_shi nX];
            end; % End of if on length(mNodeCatalog_(:,1)
            if rem(nNode_,500) == 0
                waitbar(nNode_/length(params.mPolygon(:,1)))
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                % Temporary saving
                % Witout Lognormal, but mean and median for Norm Mc
                params.vcsGridNames = cellstr(char('Mc max. curvature' , 'Mc 90% goodness of fit' , 'Mc 95% goodness of fit',...
                    'Mc best combination', 'Mc EMR-method', 'Mc(Shi-b-uncertainty)','Number of events'));
                params.mValueGrid = mValueGrid_;
                % Add parameter to params.sComment
                if  params.nGriddingMode == 0;   % Constant number
                    params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                            ' Time period (days) ' num2str(params.fTimePeriodDays) ' d, Constant number: ' num2str(params.nNumberEvents) ', MaxRadius: '...
                            num2str(params.fMaxRadius) ' km'];
                    vResults = params;
                    save(['result_Time' num2str(fTstart) '_Constnum_' num2str(params.nNumberEvents) '_MaxRad_' num2str(params.fMaxRadius)...
                            '_Nmin_' num2str(params.nMinimumNumber) '_Node' num2str(nNode_) '.mat'], 'vResults');
                elseif params.nGriddingMode == 1;   % Constant radius
                    params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                            ' Time period (days) ' num2str(params.fTimePeriodDays) ' d, Radius: ' num2str(params.fRadius) ' km, Nmin: ' num2str(params.nMinimumNumber)];
                    vResults = params;
                    save(['result_Time' num2str(fTstart) '_Rad_' num2str(params.fRadius) '_Nmin_' num2str(params.nMinimumNumber)  '_Node' num2str(nNode_) '.mat'], 'vResults');
                else  % Rectangle mode
                    params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                            ' Time period (days) ' num2str(params.fTimePeriodDays) ' d, Rect. X: ' num2str(params.fSizeRectHorizontal) ' km, Rect. Y: ' num2str(params.fSizeRectDepth)...
                            ' km, Nmin: ' num2str(params.nMinimumNumber)];
                    vResults = params;
                    save(['result_Time' num2str(fTstart) '_RectX_' num2str(params.fSizeRectHorizontal) '_RectY_' num2str(params.fSizeRectDepth)...
                            '_Nmin_' num2str(params.nMinimumNumber)  '_Node' num2str(nNode_) '.mat'], 'vResults');
                end; % END of params.nGriddingmode
                vResults =[];
            end; % End updating waitbar
        end; % for nNode
        close(hWaitbar1);
        % Parameter description
        params.vcsGridNames = cellstr(char('Mc max. curvature' , 'Mc 90% goodness of fit' , 'Mc 95% goodness of fit',...
            'Mc best combination', 'Mc EMR-method', 'Mc(Shi-b-uncertainty)','Number of events'));
        params.mValueGrid = mValueGrid_;
        % Add parameter to params.sComment
        if  params.nGriddingMode == 0;   % Constant number
            params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                    ' Time period (days) ' num2str(params.fTimePeriodDays) ' d, Constant number: ' num2str(params.nNumberEvents) ', MaxRadius: '...
                    num2str(params.fMaxRadius) ' km'];
            vResults = params;
            save(['result_Time' num2str(fTstart) '_Constnum_' num2str(params.nNumberEvents) '_MaxRad_' num2str(params.fMaxRadius)...
                    '_Nmin_' num2str(params.nMinimumNumber) '.mat'], 'vResults');
        elseif params.nGriddingMode == 1;   % Constant radius
            params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                    ' Time period (days) ' num2str(params.fTimePeriodDays) ' d, Radius: ' num2str(params.fRadius) ' km, Nmin: ' num2str(params.nMinimumNumber)];
            vResults = params;
            save(['result_Time' num2str(fTstart) '_Rad_' num2str(params.fRadius) '_Nmin_' num2str(params.nMinimumNumber) '.mat'], 'vResults');
        else  % Rectangle mode
            params.sComment = ['Starttime ' num2str(fTstart) ' Spacing ' num2str(params.fSpacingHorizontal) ' deg.,'...
                    ' Time period (days) ' num2str(params.fTimePeriodDays) ' d, Rect. X: ' num2str(params.fSizeRectHorizontal) ' km, Rect. Y: ' num2str(params.fSizeRectDepth)...
                    ' km, Nmin: ' num2str(params.nMinimumNumber)];
            vResults = params;
            save(['result_Time' num2str(fTstart) '_RectX_' num2str(params.fSizeRectHorizontal) '_RectY_' num2str(params.fSizeRectDepth)...
                    '_Nmin_' num2str(params.nMinimumNumber) '.mat'], 'vResults');
        end
        vResults =[];
        fTstart = fTstart+params.fTimePeriodDays;
    end; % End of while fTstart
end; % END of if params.bBst

author: CelsoReyes


text:
function varargout = homkermap(varargin)
% VL_HOMKERMAP Homogeneous kernel map
%   V = VL_HOMKERMAP(X, N) computes a 2*N+1 dimensional approximated
%   kernel map for the Chi2 kernel. X is an array of data points. Each
%   point is expanded into a vector of dimension 2*N+1 and saved to
%   the output V. The expanded feature vectors are stacked along the
%   first dimension, so that the output array V has the same
%   dimensions of the input array X except for the first one, which is
%   2*N+1 times larger.
%
%   The function accepts the following options:
%
%   Kernel:: KCHI2
%     One of KCHI2 (Chi2 kernel), KINTERS (intersection kernel), KJS
%     (Jensen-Shannon kernel). The 'Kernel' option name can be omitted,
%     i.e. VL_HOMKERMAP(..., 'kernel', 'kchi2') has the same effect of
%     VL_HOMKERMAP(..., 'kchi2').
%
%   Period:: [automatically tuned]
%     Set the period of the kernel specturm. The approximation is
%     based on periodicizing the kernel specturm. If not specified,
%     the period is automatically set based on the heuristic described
%     in [2].
%
%   Window:: [RECTANGULAR]
%     Set the window used to truncate the spectrum before The window
%     can be either RECTANGULAR or UNIFORM window. See [2] and the API
%     documentation for details.
%
%   Gamma:: [1]
%     Set the homogeneity degree of the kernel. The standard kernels
%     are 1-homogeneous, but sometimes smaller values perform better
%     in applications. See [2] for details.
%
%   Example::
%     The following code results in approximatively the same
%     similarities matrices between points X and Y:
%
%       x = rand(10,1) ;
%       y = rand(10,100) ;
%       psix = vl_homkermap(x, 3) ;
%       psiy = vl_homkermap(y, 3) ;
%       figure(1) ; clf ;
%       ker = vl_alldist(x, y, 'kchi2') ;
%       ker_ = psix' * psiy ;
%       plot([ker ; ker_]') ;
%
%   Note::
%     The homogeneous kernels K(X,Y) are normally defined for
%     non-negative data only. VL_HOMKERMAP defines them for both
%     positive and negative data by using the definition
%     SIGN(X)SIGN(Y)K(ABS(X),ABS(Y)) -- note that other extensions are
%     possible as well (see [2]).
%
%   REFERENCES::
%     [1] A. Vedaldi and A. Zisserman
%     `Efficient Additive Kernels via Explicit Feature Maps',
%     Proc. CVPR, 2010.
%
%     [2] A. Vedaldi and A. Zisserman
%     `Efficient Additive Kernels via Explicit Feature Maps',
%     PAMI, 2011 (submitted).
%
%   See also: VL_HELP().
[varargout{1:nargout}] = vl_homkermap(varargin{:});

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function [PA,PO,dX] = marginalise(Xdata,ydata,Xtest,GP,ns)
%                       Xdata,Ydata,GPvar,RBF_GP_var,20,20
%
%   input -----------------------------------------------------------------
%
%       ------ Sample points ----------------------
%
%       o Xtest: (D x N), points to compute the joint distribution
%
%       ------ Gaussian Process parameters --------
%
%       o X: (D x N), set of stat points
%
%       o y: (1 x N), value at the each point
%
%       -------------------------------------------
%
%       o ns: (1 x D), sizes of joint dimensions nx = ns(1), ny = ns(2),...
%
%   output ----------------------------------------------------------------
%
%       o Mu: (1 x M), Gaussian Process smoothed y values at query points
%
%       o hy: (1 x M), marginal values at chosen discrete points
%

[D,N] = size(Xtest);

max_x = max(Xtest')';
min_x = min(Xtest')';

dX = abs(max_x - min_x)./ns';

PAO = GP(Xdata,ydata,Xtest);
PAO = reshape(PAO,ns);

if D == 2

    PA = sum(PAO.*dX(2),1);
    PO = sum(PAO.*dX(1),2);
    PA = PA(:);
    PO = PO(:);

else

    PA = sum(sum(PAO .* dX(3),4) .* dX(4),3);
    PO = sum(sum(PAO .* dX(2),1) .* dX(1),2);
    PA = squeeze(PA);
    PO = squeeze(PO);

end

end

author: epfl-lasa


text:
% VL_IKMEANSPUSH  Project data on integer K-means paritions
%   I = VL_IKMEANSPUSH(X,C) projects the data X to the integer K-meanns
%   clusters of centers C returning the cluster indeces I.
%
%   See also: VL_IKMEANS(), VL_HELP().

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).


text:
function linplus_test0195 ( )

%*****************************************************************************80
%
%% TEST0195 tests R8VEC_TO_R8GB, R8GB_TO_R8VEC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  m = 5;
  n = 8;
  ml = 2;
  mu = 1;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0195\n' );
  fprintf ( 1, '  For a general banded matrix,\n' );
  fprintf ( 1, '  R8VEC_TO_R8GB converts a real vector to a R8GB matrix.\n' );
  fprintf ( 1, '  R8GB_TO_R8VEC converts a R8GB matrix to a real vector.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix rows M =    %d\n', m );
  fprintf ( 1, '  Matrix columns N = %d\n', n );
  fprintf ( 1, '  Lower bandwidth ML  = %d\n', ml );
  fprintf ( 1, '  Upper bandwidth MU  = %d\n', mu );

  a = r8gb_indicator ( m, n, ml, mu );

  r8gb_print ( m, n, ml, mu, a, '  The R8GB indicator matrix:' );

  x = r8gb_to_r8vec ( m, n, ml, mu, a );

  k = 0;
  for j = 1 : n
    for i = 1 : 2*ml+mu+1
      k = k + 1;
      fprintf ( 1, '%4d  %4d  %4d  %14f\n', i, j, k, x(k) );
    end
  end

  a = r8vec_to_r8gb ( m, n, ml, mu, x );

  r8gb_print ( m, n, ml, mu, a, '  The recovered R8GB indicator matrix:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function quadrule_test16 ( )

%*****************************************************************************80
%
%% TEST16 tests GEN_LAGUERRE_EK_COMPUTE and LAGUERRE_SUM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2011
%
%  Author:
%
%    John Burkardt
%
  order_max = 10;

  nfunc = func_set ( 'COUNT', 'DUMMY' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST16\n' );
  fprintf ( 1, '  GEN_LAGUERRE_EK_COMPUTE computes a Gauss-Laguerre rule;\n' );
  fprintf ( 1, '  LAGUERRE_SUM carries it out.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Quadrature order will vary.\n' );
  fprintf ( 1, '  Integrand will vary.\n' );
  fprintf ( 1, '  The weight function is EXP ( - X ) * X^ALPHA.\n' );
  fprintf ( 1, '\n' );

  a = 0.0;
  alpha = 2.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The integration interval is [ %f, +oo ).\n', a );
  fprintf ( 1, '  ALPHA = %f\n', alpha );
  fprintf ( 1, '\n' );

  for ilo = 1 : 5 : nfunc

    ihi = min ( ilo + 4, nfunc );

    fprintf ( 1, '\n' );
    fprintf ( 1, '    ' );
    for i = ilo : ihi
      fprintf ( '%14s', fname(i) );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '\n' );

    for norder = 1 : order_max

      fprintf ( 1, '  %2d', norder );

      for i = ilo : ihi

        func_set ( 'SET', i );

        [ xtab, weight ] = gen_laguerre_ek_compute ( norder, alpha );

        result(i) = laguerre_sum ( @func, a, norder, xtab, weight );

        fprintf ( 1, '  %12f', result(i) );

      end

      fprintf ( 1, '\n' );

    end

  end

  return
end


text:
function [hPa] = psi2hPa(psi)
% Convert units of pressure from pounds per square inch to hectopascals.
% Chad Greene 2012
hPa = psi*68.9476;

author: Sable

repo: mcbench-benchmarks


text:
%%% Loading - Step 1.
% Summary:
% load data from 'cell_resp.stackf' etc, and saved into
% 'FishX_direct_load.mat' (X = number) in the same directory.
% Stimulus parsing is done with function 'StimulusKeyPreprocessing', but
% one is advised to also inspect the results manually.

% Validation:
% There is an option 'isDiscard50' to discard the noisier 50%, disabled now.
% There is also a section to manually select (draw polygons) regions of
% cells to discard that are anatomical outliers.

% Format:
% Cell-responses are converted into z-scores and kept that way through out.
% CInfo (cell info for the valid cells) stores the anatomical positions.
% For the rest of parameters saved, see the end of this file.

% Next:
% For the following Loading steps, only this 'FishX_direct_load.mat' is
% needed, the rest will be cleared from working memory.

%%
clear all;close all;clc

%% Set Manaully!

i_fish = 1;

M_dir = GetFishDirectories();

fpsec = 1.97; % Hz

%% load data
disp(['load data: fish ' num2str(i_fish)]);
datadir = M_dir{i_fish};
if exist(fullfile(datadir,'cell_resp_dim_lowcut.mat'), 'file'),
  load(fullfile(datadir,'cell_resp_dim_lowcut.mat'));
elseif exist(fullfile(datadir,'cell_resp_dim.mat'), 'file'),
    load(fullfile(datadir,'cell_resp_dim.mat'));
else
    errordlg('find data to load!');
end

load(fullfile(datadir,'cell_info.mat'));
if exist(fullfile(datadir,'cell_resp_lowcut.stackf'), 'file'),
  cell_resp_full = read_LSstack_fast_float(fullfile(datadir,'cell_resp_lowcut.stackf'),cell_resp_dim);
elseif exist(fullfile(datadir,'cell_resp.stackf'), 'file'),
    cell_resp_full = read_LSstack_fast_float(fullfile(datadir,'cell_resp.stackf'),cell_resp_dim);
else
    errordlg('find data to load!');
end

load(fullfile(datadir,'frame_turn.mat'));

% parse stimulus - should check manually to be sure!
[stimset,stim] = StimulusKeyPreprocessing(frame_turn,i_fish);

%% load anatomy
tiffname = fullfile(datadir,'ave.tif');
info = imfinfo(tiffname,'tiff');
nPlanes = length(info);
s1 = info(1).Height;
s2 = info(1).Width;
ave_stack = zeros(s1,s2,nPlanes);
for i=1:nPlanes,
    ave_stack(:,:,i) = imread(fullfile(datadir,'ave.tif'),i);
end

% x-y view
im = max(ave_stack,[],3);
out=imNormalize99(im);
anat_yx = repmat(out,[1 1 3]);

% y-z view
im = squeeze(max(ave_stack,[],2));
out=imNormalize99(im);
anat_yz = repmat(out,[1 1 3]);

% x-z view
im = squeeze(max(ave_stack,[],1));
out=imNormalize99(im);
out = flipud(out'); %%%% empirically necessary...
anat_zx = repmat(out,[1 1 3]);

dimv_yx = size(anat_yx);
dimv_yz = size(anat_yz);
dimv_zx = size(anat_zx);

%% NEW: fix the left-right flip in the anatomy stack and subsequent cell_info
ave_stack = fliplr(ave_stack);
anat_yx = fliplr(anat_yx);
anat_zx = fliplr(anat_zx);

[s1,s2,~] = size(ave_stack);
tic
for i_cell = 1:length(cell_info),
    % fix '.center'
    cell_info(i_cell).center(2) = s2-cell_info(i_cell).center(2)+1;
    % fix '.inds'
    IX = cell_info(i_cell).inds;
    [I,J] = ind2sub([s1,s2],IX);
    J = s2-J+1;
    cell_info(i_cell).inds = sub2ind([s1,s2],I,J);
    % fix '.x_minmax'
    cell_info(i_cell).x_minmax(1) = s2-cell_info(i_cell).x_minmax(1)+1;
    cell_info(i_cell).x_minmax(2) = s2-cell_info(i_cell).x_minmax(2)+1;
end
toc

%% compute z-score
disp('compute z-score')
cell_resp_z_full = zscore(cell_resp_full')';
nCells = cell_resp_dim(1);

%% Validating all cells

isDiscard50 = false; % option to pre-screen
if isDiscard50,
    %% Round 1: discard 50% noisy cells based on std of zscore of baseline

    % Compute std of dimest 10% of frames for each cell.
    prc = prctile(cell_resp_z_full,10,2);
    STD_full = zeros(nCells,1);
    for i = 1:nCells,
        ix = find(cell_resp_z_full(i,:)<prc(i));
        STD_full(i) = std(cell_resp_z_full(i,ix));
    end
    % Set threshold at 50%, i.e. discard 50% of all cells
    thr = prctile(STD_full,50);

    temp = find(STD_full>thr);
    [~,I] = sort(STD_full(temp));
    cIX = temp(I);

    % visualize cells to discard
    % gIX = (round((1:length(cIX))/1000)+1)'; % option: view sorted index
    gIX = round(cIX/1000)+1; % option: view anatomical z index
    M = cell_resp_z_full(cIX,1:1000);
    BasicPlotMaps(cIX,gIX,M,cell_info,stim,anat_yx,anat_yz,anat_zx);

    % I_v: index of valid cells
    I_v_Holder = ones(1,nCells);
    I_v_Holder(cIX) = 0;
    I_v = find(I_v_Holder);

    cell_resp_z = cell_resp_z_full(I_v,:);
    STD = STD_full(I_v);
    nCells = length(I_v);
    CInfo_0 = cell_info(I_v);

    %% visualize valid cells, sorted by noise
    [~,I] = sort(STD);
    cIX = I;
    gIX = (round((1:length(cIX))/1000)+1)'; % sorted index
    M = cell_resp_z(I,1:1000);
    BasicPlotMaps(cIX,gIX,M,CInfo_0,stim,anat_yx,anat_yz,anat_zx);

else
    %% Round 1 alternative: full set
    cell_resp_z = cell_resp_z_full;
    I_v_Holder = ones(1,nCells);
    I_v = (1:nCells)'; % valid Index. here not actually screened, I_v is fullset

    cIX = I_v; % 'cell-Index'
    gIX = (round((1:length(cIX))/1000)+1)'; % sorted 'group-Index'
    M = cell_resp_z_full(cIX,1:1000);
    CInfo_0 = cell_info;
    BasicPlotMaps(cIX,gIX,M,CInfo_0,stim,anat_yx,anat_yz,anat_zx);
end

%% %%%% Round 2: delete anatomically out-of-bound cells by hand
% only execute once (manual choice 1):
cHolder_Anat = []; % collect cIX of all out-of-bound cells

% draw cells first
% (choose start and stop index to draw; low numbers ~ ventral)

%% (optional) step 1: set limit to only plot ventral layer, helps to find outliers within that layer
I_start = 1;
I_stop = round(length(I_v)/10);

% plot those cells
cIX = I_start:I_stop;
gIX = round(cIX/1000)+1;
figure('Position',[100 0 1300 900]);
numK = round(cIX(end)/1000)+1;
[tot_image, dim_totimage] = DrawClustersOnMap_LSh(CInfo_0,cIX,gIX,numK,anat_yx,anat_yz,anat_zx,'hsv','full');

% Here manually draw polygons around cells to discard.
% Double click to connect last vertex to first vertex, then double click again within polygon to fix.
% then click again to start drawing the next one
% and when finished, break loop with Ctrl+C...
MaskArray = zeros(dim_totimage(1), dim_totimage(2));
k_zres = 20;
% draw polygon around outliers on Y-X view
for i = 1:100, % NOMINAL LOOP, break manually (with Ctrl+C, somehow the design didn't work)
    h_poly_yx = impoly;
    wait(h_poly_yx); % double click to finalize position!
    % update finalized polygon in bright color
    setColor(h_poly_yx,[0 1 1]);

    IJs = reshape([CInfo_0(I_start:I_stop).center],2,[])';
    A = sub2ind(dimv_yx(1:2),IJs(:,1),IJs(:,2));
    MaskArray = createMask(h_poly_yx);
    MaskArray(1:dimv_zx*k_zres+10,:) = [];
    B = find(MaskArray); % find indices of pixels within ROI
    cIX2 = find(ismember(A,B));
    cHolder_Anat = union(cHolder_Anat,cIX2);
    w = waitforbuttonpress;
    if w == 1,
        break; % this doesn't work ?!
    end
end

%% step 2: plot all cells (including the ones in the ventral layer, don't
I_start = 1;
I_stop = length(I_v);

% ...here until the end of the cell is the exact dupliate of the last cell...
cIX = I_start:I_stop;
gIX = round(cIX/1000)+1;
figure('Position',[100 0 1300 900]);
numK = round(cIX(end)/1000)+1;
[tot_image, dim_totimage] = DrawClustersOnMap_LSh(CInfo_0,cIX,gIX,numK,anat_yx,anat_yz,anat_zx,'hsv','full');

% Here manually draw polygons around cells to discard.
% Double click to connect last vertex to first vertex, then double click again within polygon to fix.
% then click again to start drawing the next one
% and when finished, break loop with Ctrl+C...
MaskArray = zeros(dim_totimage(1), dim_totimage(2));
k_zres = 20;
% draw polygon around outliers on Y-X view
for i = 1:100, % NOMINAL LOOP, break manually (with Ctrl+C, somehow the design didn't work)
    h_poly_yx = impoly;
    wait(h_poly_yx); % double click to finalize position!
    % update finalized polygon in bright color
    setColor(h_poly_yx,[0 1 1]);

    IJs = reshape([CInfo_0(I_start:I_stop).center],2,[])';
    A = sub2ind(dimv_yx(1:2),IJs(:,1),IJs(:,2));
    MaskArray = createMask(h_poly_yx);
    MaskArray(1:dimv_zx*k_zres+10,:) = [];
    B = find(MaskArray); % find indices of pixels within ROI
    cIX2 = find(ismember(A,B));
    cHolder_Anat = union(cHolder_Anat,cIX2);
    w = waitforbuttonpress;
    if w == 1,
        break; % this doesn't work ?!
    end
end

%% find extra outliers on top/bottom border of image (can't always get with polygon)
temp = [CInfo_0(:).center];
XY = reshape(temp',[],nCells)';
IX = find(XY(:,1)<15 | XY(:,1)> s1-15 ...
        | XY(:,2)<15 | XY(:,2)> s2-15);
cHolder_Anat = union(cHolder_Anat,IX);

%% test Plot: all antomy outliers
cIX = cHolder_Anat;
gIX = (1:length(cIX))';
figure;
DrawClustersOnMap_LSh(CInfo_0,cIX,gIX,numK,anat_yx,anat_yz,anat_zx,'hsv','full');

%% Clean up anatomy outliers
cIX_Invalid_Anat = I_v(cHolder_Anat); % convert back to index for original full set
I_v_Holder(cIX_Invalid_Anat) = 0;
I_v2 = find(I_v_Holder);

CR_raw = cell_resp_z_full(I_v2,:);
nCells = length(I_v2);
CInfo = cell_info(I_v2);

%% Plot all valid cells to save
cIX = (1:nCells)';
gIX = round(cIX/1000)+1;
M = CR_raw(cIX,1:1000);
figure;
DrawClustersOnMap_LSh(CInfo,cIX,gIX,numK,anat_yx,anat_yz,anat_zx,'hsv','full');

%% Save mat files
CInfo_full = cell_info;
tic
temp = fullfile(datadir,['Fish' num2str(i_fish) '_full_extrainfo.mat']);
save(temp,'cHolder_Anat','cIX_Invalid_Anat','I_v','I_v2','CInfo_full','-v7.3'); % CInfo_full not saved in mats before 4/24/15!! %% 'STD_full' saved before

temp = fullfile(datadir,['Fish' num2str(i_fish) '_direct_load_nodiscard.mat']);
varList = {'CR_raw','nCells','CInfo','anat_yx','anat_yz','anat_zx','ave_stack','fpsec','frame_turn'};
save(temp,varList{:},'-v7.3');
toc

%% (if needed)
% save(temp,'somethingelse','-append');

%% next run step 2: detrend

author: xiuyechen

repo: FishExplorer


text:
% Euclidean projection on a hyperplane
% Section 8.1.1, Boyd & Vandenberghe "Convex Optimization"
% Joelle Skaf - 10/04/05
%
% The projection of x0 on a hyperplane C = {x | a'*x = b} is given by
%           minimize || x - x0 ||^2
%               s.t.    a'*x = b
% It is also given by P_C(x0) = x0 + (b - a'*x0)*a/||a||^2

% Input data
randn('seed',0);
n  = 10;
a  = randn(n,1);
b  = randn(1);
x0 = randn(n,1);

% Analytical solution
fprintf(1,'Computing the analytical solution ...');
pc_x0 = x0 + (b - a'*x0)*a/norm(a)^2;
fprintf(1,'Done! \n');

% Solution via QP
fprintf(1,'Computing the optimal solution by solving a QP ...');

cvx_begin quiet
    variable x(n)
    minimize ( norm(x-x0)^2 )
    a'*x == b; %#ok
cvx_end

fprintf(1,'Done! \n');

% Verification
disp('--------------------------------------------------------------------------------');
disp('Verifying that p_C(x0) and x_star belong to the hyperplane C: ');
disp(['a^T*p_C(x0) - b = ' num2str(a'*pc_x0 - b)]);
disp(['a^T*x_star - b  = ' num2str(a'*x - b)]);
disp('Computing the distance between x0 and the hyperplane in each case');
disp(['||x0 - p_C(x0)|| = ' num2str(norm(x0 - pc_x0))]);
disp(['||x0 - x_star || = ' num2str(norm(x0 - x))]);
disp('Verifying that the analytical solution and the solution obtained via QP are equal: ');
disp([pc_x0 x])

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
% This file is part of the project NILM-Eval (https://github.com/beckel/nilm-eval).
% Licence: GPL 2.0 (http://www.gnu.org/licenses/gpl-2.0.html)
% Copyright: ETH Zurich, 2014
% Author: Christian Beckel;

figure_path = 'projects/buildsys/images/pie_chart/';
fontsize = 9;
% width = 4.8;
% height = 4.8;
households = [ 5 ];
if ~exist(figure_path, 'dir')
    mkdir(figure_path);
end

for household = households
    %params
    dataset = 'eco';
    % evalDays_type = 'completeSM_first90';
    evalDays_type = 'plug_statistics';
    granularity = 1;

    path_to_evalDays = strcat(pwd, '/input/autogen/evaluation_days/', dataset, '/', evalDays_type, '/', num2str(household, '%02d'), '.mat');
    load(path_to_evalDays); % evalDays

    % total consumption
    fprintf('Processing household %d, smart meter\n', household);
    smartmeter_consumption = read_smartmeter_data(dataset, household, evalDays, granularity, 'powerallphases');
    missing_values_idx_sm = smartmeter_consumption == -1;
    total_consumption = sum(smartmeter_consumption(~missing_values_idx_sm));

    % consumption of each appliance
    appliances = findAppliances(household, 'eco');
    % power_pct = zeros(length(appliances)+1, 1);
    avg_cons = zeros(length(appliances)+1, 1);
    for i = 1:length(appliances)
       fprintf('Processing household %d, plug %d/%d\n', household, i, length(appliances));
       plug_consumption = read_plug_data(dataset, household, appliances(i), evalDays, granularity);
       missing_values_idx_plug = plug_consumption == -1;
       sum_plug_consumption = sum(plug_consumption(~missing_values_idx_plug));
    %    ratio = nnz(~missing_values_idx_sm) / nnz(~missing_values_idx_plug);
    %    extrapolated_plug_consumption = ratio*sum_plug_consumption;
    %    power_pct(i) = extrapolated_plug_consumption / total_consumption;
        avg_cons(i) = sum_plug_consumption / nnz(~missing_values_idx_plug);
    end

    % power_pct(end) = 1 - sum(power_pct(~isnan(power_pct)));
    avg_cons(end) = (total_consumption / nnz(~missing_values_idx_sm)) - sum(avg_cons);

%     % plot pie chart
%     if household == 5
%         names{end+1} = '';
%         avg_cons(end+1) = avg_cons(end);
%         avg_cons(end-1) = 0;
%     end
%
    names = getApplianceNames(appliances);
    names{end+1} = 'Other';
    fig = figure();
    %title(strcat('Household: ', num2str(household)));
    % pie_handle = pie(power_pct);
    pie_handle = pie(avg_cons);
    % set(pie_handle(2:2:end));
    handle_array = findobj(pie_handle, 'Type', 'text');
%     percentages = get(handle_array, 'String');
%     new_labels = strcat(names',{': '}, percentages);
    new_labels = strcat(names', {': '}, num2str(avg_cons * 24/1000*30, 3));
    set(handle_array, {'String'}, new_labels);

%     fig = make_report_ready(fig, 'size', [width height], 'fontsize', fontsize);
    fig = make_report_ready(fig, 'fontsize', fontsize);

    filename = ['household_', num2str(household)];
    print('-depsc2', '-cmyk', '-r600', [figure_path, filename, '.eps']);
  	% saveas(fig_h, [figure_path, filename, '.png'], 'png');
    pause(1);
  	close(fig);

end

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function out=janaf(prop, spec,T);
% ---------------------------------------------------------------------
% function out=janaf(prop, spec, T)                 |   Version 1.01
% ---------------------------------------------------------------------
% Calculates JANAF curve fit according to JANAF virial equation
% Output is calculated in SI-units
%
% prop = 'c' for standard state specific heat
%      = 'h' for standard state enthalpy
%      = 's' for standard state entropy
% spec = 'CO2', 'H2O', 'CO', 'H2', 'O2', 'N2'
% T    = Temperature, vector allowed
%
% JANAF.mat required (contains the coefficients)
% ---------------------------------------------------------------------
% Last Change: 2003-07-18           |   (c)2003, Stefan Billig, Delphi

% check for correct syntax
if nargin~=3
    help janaf
    % end function
    return
end

% load coefficient table
load JANAF;
z=1;
% determine molecular weight from table
MWeight=eval(['MolWeight.' spec]);

for i=1:length(T)
    % choose temperature range vector
    if (T(i)>1000 & T(i)<=5000)
        eval(['ai=' spec '(1,:);'])
        out(z)=calc(ai, T(i), prop, MWeight);
        z=z+1;
    elseif (T(i)>=300 & T(i)<=1000)
        eval(['ai=' spec '(2,:);'])
        out(z)=calc(ai, T(i), prop, MWeight);
        z=z+1;
    else
        sprintf(['Temperature ' num2str(T(i)) 'K not between 300K and 5000K!'])
    end
end

%----------------------------------------------------------------------
function out=calc(ai, T, prop, MWeight)

R=8.314;
% calculate standard state value
switch prop
    case 'c'
        out=(ai(1)+ai(2)*T+ai(3)*T.^2+ai(4)*T.^3+ai(5)*T.^4)*R/MWeight;
    case 'h'
        out=(ai(1)+ai(2)/2*T+ai(3)/3*T.^2+ai(4)/4*T.^3+ai(5)/5*T.^4+ai(6)/T).*T*R/MWeight;
    case 's'
        out=(ai(1)*ln(T)+ai(2)*T+ai(3)/2*T.^2+ai(4)/3*T.^3+ai(5)/4*T.^4+ai(7))*R/MWeight;
end


text:
% This script test under 10-fold Cross Validation, divided by ID, how dtfv works
function [AUC, tt, tf, ft, ff] = auc_fun(label, score, pred)
lie_id = find(score<0);
if ~isempty(lie_id)
    if pred(lie_id(1)) == 0
        isign = -1;
    else
        isign = 1;
    end
else
    if pred(1) == 1
        isign = -1;
    else
        isign = 1;
    end
end
[~,~,~,AUC] = perfcurve(label,isign*score,0);
tt = sum(pred(label==1)==1);
tf = sum(pred(label==1)==0);
ft = sum(pred(label==0)==1);
ff = sum(pred(label==0)==0);

%fprintf('Partation:%d\tacc is %f, auc is %f, tt:%d, tf:%d, ft:%d, ff:%d\n',part, acc(1),AUC,tt,tf,ft,ff);
end

author: Doubaibai

repo: DARE

save_path: github-repos/MATLAB/Doubaibai-DARE


text:
function varargout = createTetrakaidecahedron()
%CREATETETRAKAIDECAHEDRON Create a 3D mesh representing a tetrakaidecahedron.
%
%   [V, E, F] = createTetrakaidecahedron;
%   Create a mesh structure representing a tetrakaidecahedron, composed of
%   both square and hexagonal faces. Tetrakaidecahedron can be used to tile
%   the 3D Euclidean space.
%
%   V is a 24-by-3 array with vertex coordinates,
%   E is a 36-by-2 array containing indices of neighbour vertices,
%   F is a 14-by-1 cell array containing vertex indices array of each face.
%
%   [V, F] = createTetrakaidecahedron;
%   Returns only the vertices and the face vertex indices.
%
%   MESH = createTetrakaidecahedron;
%   Returns the data as a mesh structure, with fields 'vertices', 'edges'
%   and 'faces'.
%
%   Example
%   [n, e, f] = createTetrakaidecahedron;
%   drawMesh(n, f);
%
%   See also
%   meshes3d, drawMesh

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2005-02-10
% Copyright 2005-2022 INRA - TPV URPOI - BIA IMASTE

nodes = [...
    1 0 2;0 1 2;-1 0 2;0 -1 2;...
    2 0 1;0 2 1;-2 0 1;0 -2 1;...
    2 1 0;1 2 0;-1 2 0;-2 1 0;-2 -1 0;-1 -2 0;1 -2 0;2 -1 0;...
    2 0 -1;0 2 -1;-2 0 -1;0 -2 -1;...
    1 0 -2;0 1 -2;-1 0 -2;0 -1 -2];

edges = [...
    1 2;1 4;1 5;2 3;2 6;3 4;3 7;4 8;...
    5 9;5 16;6 10;6 11;7 12;7 13;8 14;8 15;...
    9 10;9 17;10 18;11 12;11 18;12 19;13 14;13 19;14 20;15 16;15 20;16 17;....
    17 21;18 22;19 23;20 24;21 22;21 24;22 23;23 24];

faces = {...
    [1 2 3 4], ...
    [1 4 8 15 16 5], [1 5 9 10 6 2], [2 6 11 12 7 3], [3 7 13 14 8 4],...
    [5 16 17 9], [6 10 18 11], [7 12 19 13], [8 14 20 15],...
    [9 17 21 22 18 10], [11 18 22 23 19 12], [13 19 23 24 20 14], [15 20 24 21 17 16], ...
    [21 24 23 22]};
faces = faces';

% format output
varargout = formatMeshOutput(nargout, nodes, edges, faces);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function g=firkaiser(L,beta,varargin)
%FIRKAISER  Kaiser-Bessel window
%   Usage:  g=firkaiser(L,beta);
%           g=firkaiser(L,beta,...);
%
%   `firkaiser(L,beta)` computes the Kaiser-Bessel window of length *L* with
%   parameter *beta*. The smallest element of the window is set to zero when
%   the window has an even length. This gives the window perfect whole-point
%   even symmetry, and makes it possible to use the window for a Wilson
%   basis.
%
%   `firkaiser` takes the following flags at the end of the input arguments:
%
%     'normal'   Normal Kaiser-Bessel window. This is the default.
%
%     'derived'  Derived Kaiser-Bessel window.
%
%     'wp'       Generate a whole point even window. This is the default.
%
%     'hp'       Generate half point even window.
%
%   Additionally, `firkaiser` accepts flags to normalize the output. Please
%   see the help of |setnorm|. Default is to use `'null'` normalization.
%
%   Note that odd-length Derived Kaiser-Bessel windows are not
%   mathematically defined, yet they are supported by this code.
%
%   See also: firwin, setnorm
%
%   References: opsc89

%   AUTHOR: unknown. Additions by Clara Hollomey

if nargin<2
  error('Too few input arguments.');
end;

if numel(beta)>1
  error('beta must be a scalar.');
end;

% Define initial value for flags and key/value pairs.
definput.import={'setnorm'};
definput.importdefaults={'null'};
definput.flags.centering={'wp','hp'};
definput.flags.stype={'normal','derived'};

[flags,keyvals]=ltfatarghelper({},definput,varargin);

cent=0;
if flags.do_hp
  cent=.5;
end;

if flags.do_normal

  if (L == 1)
    g = 1;
  else
    m = L - 1;
    k = (0:m)'+rem(L,2)/2-.5+cent;
    k = 2*beta/(m)*sqrt(k.*(m-k));
    g = besseli(0,k)/besseli(0,beta);
  end;

  g=ifftshift(g);

  if ((flags.do_wp && rem(L,2)==0) || ...
      (flags.do_hp && rem(L,2)==1))

    % Explicitly zero last element. This is done to get the right
    % symmetry, and because that element sometimes turns negative.
    g(floor(L/2)+1)=0;
  end;

else

  %if rem(L,2)==1
  %  error('The length of the choosen window must be even.');
  %end;

  if flags.do_wp
    %if rem(L,4)==0
    %  L2=L/2+2;
    %else
      L2=floor(L/2+1);
    %end;
  else
    L2=floor((L+1)/2);
  end;

  % Compute a normal Kaiser window
  g_normal=fftshift(firkaiser(L2,beta,flags.centering));

  g1=sqrt(cumsum(g_normal(1:L2))./sum(g_normal(1:L2)));

  if flags.do_wp
    if rem(L,2)==0
      g=[flipud(g1);...
         g1(2:L/2)];
    else
      g=[flipud(g1);...
         g1(1:floor(L/2))];
    end;
  else
      if rem(L,2)==0
         g=[flipud(g1);0;...
         g1(1:end-1)];
      else
         g=[flipud(g1);...
         g1(1:end-1)];
      end
  end;

  if ((flags.do_wp && rem(L,2)==0)) %|| ...
      %(flags.do_hp && rem(L,2)==1))

    % Explicitly zero last element. This is done to get the right
    % symmetry, and because that element sometimes turn negative.
    g(floor(L/2)+1)=0;
  end;

end;

% The besseli computation sometimes generates a zero imaginary component.
g=real(g);

g=setnorm(g,flags.norm);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
%% Plane
% Concrete subclass of <ZeroVolShape.html |Shape|> representing a plane.

%%% Description
% |Plane| represents the shape of a plane.  It does not have a volume, but it is
% used to force a user-defined primary grid plane in <maxwell_run.html
% |maxwell_run|>.

%%% Construction
%  shape = Plane(normal_axis, intercept)
%  shape = Plane(normal_axis, intercept, dl_max)

% *Input Arguments*
%
% * |normal_axis|: axis normal to the plane.  It should be one of |Axis.x|,
% |Axis.y|, |Axis.z|.
% * |intercept|: location of the plane in the |normal_axis| direction.
% * |dl_max|: maximum grid size allowed in the plane.  It can be either |[dx dy
% dz]| or a single real number |dl| for |dx = dy = dz|.  If unassigned, |dl_max
% = Inf| is used.  In the |normal_axis| direction, |dl_max| is meaningless.

%%% Example
%   % Create an instance of Plane.
%   shape = Plane(Axis.y, 100);
%
%   % Use the constructed shape in maxwell_run().
%   [E, H] = maxwell_run({INITIAL ARGUMENTS}, 'OBJ', {'vacuum', 'none', 1.0}, shape, {REMAINING ARGUMENTS});

%%% See Also
% <Rectangle.html |Rectangle|>, <Line.html |Line|>, <Point.html |Point|>,
% <ZeroVolShape.html |ZeroVolShape|>, <maxwell_run.html |maxwell_run|>

classdef Plane < ZeroVolShape

	properties
		normal_axis
		intercept
	end

	methods
        function this = Plane(normal_axis, intercept, dl_max)
			chkarg(istypesizeof(normal_axis, 'Axis'), '"normal_axis" should be instance of Axis.');
			chkarg(istypesizeof(intercept, 'real'), '"intercept" should be real.');

			function level = lsf(x, y, z)
				chkarg(istypeof(x, 'real'), '"x" should be array with real elements.');
				chkarg(istypeof(y, 'real'), '"y" should be array with real elements.');
				chkarg(istypeof(z, 'real'), '"z" should be array with real elements.');
				chkarg(isequal(size(x), size(y), size(z)), '"x", "y", "z" should have same size.');

				loc = {x, y, z};
				level = -abs(loc{normal_axis} - intercept);
			end

			lprim = cell(1, Axis.count);
			for w = Axis.elems
				if w == normal_axis
					lprim{w} = intercept;
				else
					lprim{w} = [-Inf Inf];
				end
			end

			if nargin < 3  % no dl_max
				super_args = {lprim, @lsf};
			else
				dl_max = expand2row(dl_max, Axis.count);
				dl_max(normal_axis) = Inf;  % dl_max is meaningless in normal direction
				super_args = {lprim, @lsf, dl_max};
			end

			this = this@ZeroVolShape(super_args{:});
			this.normal_axis = normal_axis;
			this.intercept = intercept;
		end
	end
end

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


function price = PROJ_GMXB_Surrender(T, M, gmdb_params, modelInputs, N, alph)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Description: Price Gauranteed Minimum death benefits with period fees and early surrender in Levy Models
%              using the PROJ method
%
% Author:      Justin Lars Kirkby
% References:  (1) Valuation and optimal surrender of variable annuities
%                   with guaranteed minimum benefits and periodic fees,
%                   Kirkby and Aguilar 2022, Scandinavian Actuarial Journal
%
%              (2) Efficient Option Pricing By Frame Duality with The Fast
%              Fourier Transform, SIAM J. Financial Math., 2015
%
% Models Supported: Levy Processes, including jump diffusions and Black-Scholes model
% Returns: price of contract
%
% ----------------------
% Contract/Model Params
% ----------------------
% T   = time remaining until maturity
% M   = number of subintervals of [0,T] (total of M+1 monitoring points in time grid, including S_0)
%
% modelInputs
% ------------------------
% r     = interest rate
% q     = dividend yield
% rnCHF = risk netural characteristic function (function handle with single argument)
%
% gmdb_params
% ------------------------
% F_0       = initial fund value
% alpha_fee = period fee rate
% gamma     = surrender penalty rate, 1.0 = 100% fund lost upon surrender, 0.0 = 0% lost, ie no penalty
% g         = floor on growth
% c         = cap on growth
% death_prob_cond = conditional probability of death (mortality) table, based on age
%                   To construct see function make_mortality_table_pmf, set conditional=1
%
% ----------------------
% Numerical (PROJ) Params
% ----------------------
% alph =  grid with is 2*alph
% N  = number of grid points (power of 2, e.g. 2^12), resolution = 2*alph/(N-1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dt   = T/M;

r = modelInputs.r;
q = modelInputs.q;

rnCHF = modelInputs.rnCHF;
rnCHF = @(u)rnCHF(u).*exp(-1i*(r-q)*u*dt);  % undo the drift

%%%%%%%%%%%%%%%%%%%%%%%%%%
F_0 = gmdb_params.F_0;

g = gmdb_params.g;
c = gmdb_params.c;  % Cap on growth
c_S = c;  % Cap on growth for surrender
gamma = gmdb_params.gamma;

alpha_fee = gmdb_params.alpha_fee;
death_prob = gmdb_params.death_prob_cond;

%%% NOTE: assume constant fee and r,q
lambda_m = log(1-alpha_fee) + (r-q)*dt;

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

K = N/2;
dx = 2*alph/(N-1); a = 1/dx;

nnot = K/2;  % TODO: make sure 0 is on grid

%%%% Populate Beta coefficients for orthogonal projection
xmin = (1-K/2)*dx;

a2    = a^2;
Cons2 = 24*a2*exp(-r*dt)/N;
zmin  = (1 - K)*dx;  %Kbar corresponds to zero

dw    = 2*pi*a/N;
grand = (dw: dw: (N-1)*dw);
grand = exp(-1i*zmin*grand).*rnCHF(grand).*(sin(grand/(2*a))./grand).^2./(2+cos(grand/a));
beta  = Cons2*real(fft([1/(24*a2) grand]));   %%%%  NOTE: all toep matrices incorporate exp(-r*dt)

toepM = [beta(K:-1:1)'; 0 ; beta(2*K-1:-1:K +1)'];
toepM = fft(toepM);

%%%% Initial terminal payoff coefficients (recursion proceeds backwards in time)
grid = xmin + dx*(0:K-1);
expGrid = exp(grid);

% Initialize Value function at maturity
Vals = F_0 * max(exp(g*T), min(exp(c*T), expGrid));

Thet = zeros(K,1);

cusum_right = cumsum(beta(2*K:-1:K +1))';
cusum_left = [ fliplr(cumsum(beta(1:1:K-1)))';0];

for m=M-1:-1:0

    % Define the value function at next time
    Vfunc = @(x)spline(grid, Vals, x);

    % Define maturity/death benefit function
    Mfunc = @(x) F_0 * max(exp(g*(m+1)*dt), min(exp(c*(m+1)*dt), exp(x)));

    % Combine the two functions, weight by death prob, and shift by lambda
    pw = death_prob(m+1);
    H = @(x) pw*Mfunc(lambda_m + x) + (1-pw)*Vfunc(lambda_m + x);

    Vals = H(grid);

    Thet(2:K -1) = (Vals(1:K-2)+10*Vals(2:K-1)+Vals(3:K))/12;
    Thet(1)      = (13*Vals(1)+15*Vals(2)-5*Vals(3)+Vals(4))/48;
    Thet(K)      = 2*(13*Vals(K)+15*Vals(K-1)-5*Vals(K-2)+Vals(K-3))/48;

    p = ifft(toepM.*fft([Thet(1:K);zeros(K,1)]));

    % update values with augmentation on the boundaries
    Vals(1:K) = p(1:K) + cusum_left*Vals(1) + cusum_right*Vals(end);

    if m > 0 && gamma < 1  %%% Cant surrender in first period
        surrender = (1-gamma)*F_0*min(exp(c_S*m*dt), expGrid);
        Vals = max(Vals, surrender);
    end
end

price = Vals(nnot);

end


text:
function value = r8_exprel ( x )

%*****************************************************************************80
%
%% R8_EXPREL evaluates the exponential relative error term of an R8 argument.
%
%  Discussion:
%
%    The relative error term is ( exp ( x ) - 1 ) / x.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the exponential relative error term
%    at X.
%
  persistent nterms
  persistent xbnd

  if ( isempty ( nterms ) )
    alneps = log ( r8_mach ( 3 ) );
    xn = 3.72 - 0.3 * alneps;
    xln = log ( ( xn + 1.0 ) / 1.36 );
    nterms = r8_aint ( xn - ( xn * xln + alneps ) / ( xln + 1.36 ) + 1.5 );
    xbnd = r8_mach ( 3 );
  end

  absx = abs ( x );

  if ( absx < xbnd )
    value = 1.0;
  elseif ( absx <= 0.5 )
    value = 0.0;
    for i = 1 : nterms
      value = 1.0 + value * x / ( nterms + 2 - i );
    end
  else
    value = ( exp ( x ) - 1.0 ) / x;
  end

  return
end


text:
function J=spherAngGradient(xG,systemType,lRx,M)
%%SPHERANGGRADIENT Determine the gradient of a 3D spherical azimuth and
%          elevation measurement with respect to 3D position. Relativity
%          and atmospheric effects are not taken into account.
%
%INPUTS: xG The 3XN position vectors in the global coordinate system, each
%          with [x;y;z] components.
% systemType An optional parameter specifying the axes from which the
%          angles are measured in radians. Possible vaues are
%          0 (The default if omitted) Azimuth is measured counterclockwise
%            from the x-axis in the x-y plane. Elevation is measured up
%            from the x-y plane (towards the z-axis). This is consistent
%            with common spherical coordinate systems for specifying
%            longitude (azimuth) and geocentric latitude (elevation).
%          1 Azimuth is measured counterclockwise from the z-axis in the
%            z-x plane. Elevation is measured up from the z-x plane
%            (towards the y-axis). This is consistent with some spherical
%            coordinate systems that use the z axis as the boresight
%            direction of the radar.
%          2 This is the same as 0 except instead of being given
%            elevation, one desires the angle away from the z-axis, which
%            is (pi/2-elevation).
%          3 This is the same as 0 except azimuth is measured clockwise
%            from the y-axis in the x-y plane instead of counterclockwise
%            from the x-axis. This coordinate system often arises when
%            given "bearings" in a local East-North-Up coordinate system,
%            where the bearing directions are measured East of North.
%      lRx The 3X1 position vector of the receiver. If omitted, the
%          receiver is placed at the origin. All vectors in x are assumed
%          to be from the same receiver.
%        M A 3X3 rotation matrix from the global Coordinate system to the
%          orientation of the coordinate system at the receiver. If
%          omitted, it is assumed to be the identity matrix. All vectors in
%          x are assumed to have the same rotation matrix.
%
%OUTPUTS: J A 2X3XN set of N Jacobian matrices where the rows are
%           [azimuth;elevation] in that order and the columns take the
%           partial derivative of the row components with respect to
%           [x,y,z] in that order.
%
%The derivatives can be computed in a straightforward manner from
%the basic relation between spherical and Cartesian coordinates, which is
%given in Chapter 14.4.4.1 of [1], among other sources, for systemType=0.
%
%Note that singularities exist at the poles; that is when the elevation is
%+/-(pi/2).
%
%EXAMPLE:
%The results can be verified by comparing the value from numerical
%differentiation to the value returned by this function.
% x=[100;10;-80];
% systemType=0;
% f=@(x)getSpherAngle(x,systemType);
% fDim=2;
% numDiffVal=numDiff(x,f,fDim)
% funcVal=spherAngGradient(x,systemType)
%One will see that within reasonable finite precision bounds, the two
%values are equal.
%
%REFERENCES:
%[1] R. L. Duncombe, "Computational techniques," in Explanatory Supplement
%    to the Astronomical Almanac, 3rd ed., S. E. Urban and P. K.
%    Seidelmann, Eds. Mill Valley, CA: University Science Books, 2013.
%
%February 2017 David F.Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<4||isempty(M))
   M=eye(3,3);
end

if(nargin<3||isempty(lRx))
   lRx=zeros(3,1);
end

if(nargin<2||isempty(systemType))
   systemType=0;
end

N=size(xG,2);

%Transform from global coordinates to local coordinates.
xLocal=M*bsxfun(@minus,xG(1:3,:),lRx(1:3));

JTotal=zeros(2,3,N);
for curPoint=1:N
    x=xLocal(1,curPoint);
    y=xLocal(2,curPoint);
    z=xLocal(3,curPoint);

    r=norm(xLocal(:,curPoint));
    J=zeros(2,3);
    switch(systemType)
        case 0
            %Derivatives with respect to x.
            J(1,1)=-y/(x^2+y^2);
            J(2,1)=-x*z/(r^2*sqrt(x^2+y^2));

            %Derivatives with respect to y.
            J(1,2)=x/(x^2+y^2);
            J(2,2)=-y*z/(r^2*sqrt(x^2+y^2));

            %Derivatives with respect to z.
            J(1,3)=0;
            J(2,3)=sqrt(x^2+y^2)/r^2;
        case 1
            %Derivatives with respect to x.
            J(1,1)=z/(z^2+x^2);
            J(2,1)=-x*y/(r^2*sqrt(z^2+x^2));

            %Derivatives with respect to y.
            J(1,2)=0;
            J(2,2)=sqrt(z^2+x^2)/r^2;

            %Derivatives with respect to z.
            J(1,3)=-x/(z^2+x^2);
            J(2,3)=-z*y/(r^2*sqrt(z^2+x^2));
        case 2
            %Derivatives with respect to x.
            J(1,1)=-y/(x^2+y^2);
            J(2,1)=x*z/(r^2*sqrt(x^2+y^2));

            %Derivatives with respect to y.
            J(1,2)=x/(x^2+y^2);
            J(2,2)=y*z/(r^2*sqrt(x^2+y^2));

            %Derivatives with respect to z.
            J(1,3)=0;
            J(2,3)=-sqrt(x^2+y^2)/r^2;
        case 3
            %Derivatives with respect to x.
            J(1,1)=y/(x^2+y^2);
            J(2,1)=-x*z/(r^2*sqrt(x^2+y^2));

            %Derivatives with respect to y.
            J(1,2)=-x/(x^2+y^2);
            J(2,2)=-y*z/(r^2*sqrt(x^2+y^2));

            %Derivatives with respect to z.
            J(1,3)=0;
            J(2,3)=sqrt(x^2+y^2)/r^2;
        otherwise
            error('Invalid system type specified.')
    end

    %Rotate from local back to global coordinates.
    JTotal(:,:,curPoint)=J*M;
end
J=JTotal;
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function varargout = clipMeshVertices(v, f, b, varargin)
%CLIPMESHVERTICES Clip vertices of a surfacic mesh and remove outer faces
%
%   [V2, F2] = clipMeshVertices(V, F, B)
%   Clip a mesh represented by vertex array V and face array F, with the
%   box represented by B. The result is the set of vertices contained in
%   the box, and a new set of faces corresponding to original faces with
%   all vertices within the box.
%
%   [V2, F2] = clipMeshVertices(..., 'shape','sphere') Specify the shape.
%   Default is 'box'. But it's also possible to use a 'sphere'.
%
%   [V2, F2] = clipMeshVertices(..., 'inside',false) removes the inner
%   faces instead of the outer faces.
%
%   Example
%     [v, f] = createSoccerBall;
%     f = triangulateFaces(f);
%     box = [0 2 -1 2 -.5 2];
%     [v2, f2] = clipMeshVertices(v, f, box, 'inside', false);
%     figure('color','w'); view(3); axis equal
%     drawMesh(v, f, 'faceColor', 'none', 'faceAlpha', .2);
%     drawBox3d(box)
%     drawMesh(v2, f2, 'faceAlpha', .7);
%
%   See also
%   meshes3d, clipPoints3d
%
% ------
% Author: David Legland, oqilipo
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-04-07,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

% if input is given as a structure, parse fields
if isstruct(v)
    if nargin > 2; varargin = [b, varargin]; end
    b = f;
    f = v.faces;
    v = v.vertices;
end

parser = inputParser;
validStrings = {'box','sphere'};
addParameter(parser,'shape','box',@(x) any(validatestring(x, validStrings)));
addParameter(parser,'inside',true,@islogical);
parse(parser,varargin{:});

% clip the vertices
[v2, indVertices] = clipPoints3d(v, b,...
    'shape', parser.Results.shape, 'inside', parser.Results.inside);

% create index array for face indices relabeling
refInds = zeros(size(indVertices));
for i = 1:length(indVertices)
    refInds(indVertices(i)) = i;
end

% select the faces with all vertices within the box
if isnumeric(f)
    % Faces given as numeric array
    indFaces = sum(~ismember(f, indVertices), 2) == 0;
    f2 = refInds(f(indFaces, :));

elseif iscell(f)
    % Faces given as cell array
    nFaces = length(f);
    indFaces = false(nFaces, 1);
    for i = 1:nFaces
        indFaces(i) = sum(~ismember(f{i}, indVertices), 2) == 0;
    end
    f2 = f(indFaces, :);

    % re-label indices of face vertices (keeping horizontal index array)
    for i = 1:length(f2)
        f2{i} = refInds(f2{i})';
    end
end

switch nargout
    case 1
        mesh2.vertices=v2;
        mesh2.faces=f2;
        varargout{1}=mesh2;
    case 2
        varargout{1}=v2;
        varargout{2}=f2;
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips