text:
function [neg, zer, nd, pos] = convertGene2PathwayInteractions(E, epSys, uSys)
% This function was written to find out the type of epistatic interactions
% that either exist between two different or within same subsystem
%
% USAGE:
%
%     [neg,zer,nd,pos] = convertGene2PathwayInteractions(E,epSys,uSys)
%
% INPUT:
%    E:     A square epistatic interaction matrix (or genes or reactions)
%    epSys: Subsystems belonging to the gene at that index (a second
%           order cell array) Each cell array may contain one or more
%           subsystems.
%           e.g. for first gene or reaction in E:
%           epCmpt{1,1} = {'Glycolysis';'TCA cycle'}
%    uSys:  Unique Subsystems in the epCmpt (use this arguments if
%           interested in interactions between selected subsystems in the epSys)
%
% OUTPUTS:
%    neg:            matrix of number of aggravating (negative) pathway-pair interactions
%    zer:            matrix of number of no-epistatic interactions pathway-pair interactions
%    nd:             matrix of number of non-decisive interactions pathway-pair interactions
%    pos:            matrix of number of buffering (positive) interactions pathway-pair interactions
%
% NOTE:
%    See figures in following publication:
%    Joshi CJ and Prasad A, 2014, "Epistatic interactions among metabolic genes
%    depend upon environmental conditions", Mol. BioSyst., 10, 2578-2589.
%
% .. Authors:
%     - Chintan Joshi 10/26/2018

if (nargin < 3)
    uSys = unique(convertMyCell2List(epSys, 2));
end

negE = (E < -0.25);
zerE = (E >= -0.25) & (E < 0.25);
ndE = (E >= 0.25) & (E < 0.85);
posE = (E >= 0.85);
g1 = 0;
g2 = 0;
pos = zeros(length(uSys), length(uSys)); neg = pos; nd = pos; zer = pos;
for i = 1:length(uSys)
    index1 = cellfun(@strcmp, epSys, repmat(uSys(i), length(epSys), 1));
    for j = 1:length(uSys)
        index2 = cellfun(@strcmp, epSys, repmat(uSys(j), length(epSys), 1));
        pos(i, j) = sum(sum(posE(index1, index2)));
        neg(i, j) = sum(sum(negE(index1, index2)));
        nd(i, j) = sum(sum(ndE(index1, index2)));
        zer(i, j) = sum(sum(zerE(index1, index2)));
    end
end

function B = convertMyCell2List(A, dimSense)

% This function linearizes a cell array if some of the cells in the array are another embbedded cell arrays.
% (hence, two degrees of cell array)

% This will only work for atmost two degrees of cell array.

% A=cell array to be linearized.
% dimSense=determines if inner cells are row (dimSense=1) vectors or columns (dimSense=2, default)

cnt = 0;
if nargin < 2
    dimSense = 2;
end
for i = 1:length(A)
    for j = 1:length(A{i, 1})
        cnt = cnt + 1;
        if dimSense == 1
            B{cnt, 1} = A{i, 1}{1, j};
        elseif dimSense == 2
            B{cnt, 1} = A{i, 1}{j, 1};
        end
    end
end

author: opencobra

repo: cobratoolbox


text:
function [y] = spm_fx_fmri_linear(x,u,P,M)
% state equation for a dynamic model of fMRI (linear version)
% responses
% FORMAT [y] = spm_fx_fmri_linear(x,u,P,M)
% x      - state vector
%   x(:,1) - excitatory neuronal activity             ue
%   x(:,2) - vascular signal                          s
%   x(:,3) - rCBF                                  ln(f)
%   x(:,4) - venous volume                         ln(v)
%   x(:,5) - deoyxHb                               ln(q)
%  [x(:,6) - inhibitory neuronal activity             ui]
%
% y      - dx/dt
%
%___________________________________________________________________________
%
% References for hemodynamic & neuronal state equations:
% 1. Buxton RB, Wong EC & Frank LR. Dynamics of blood flow and oxygenation
%    changes during brain activation: The Balloon model. MRM 39:855-864,
%    1998.
% 2. Friston KJ, Mechelli A, Turner R, Price CJ. Nonlinear responses in
%    fMRI: the Balloon model, Volterra kernels, and other hemodynamics.
%    Neuroimage 12:466-477, 2000.
% 3. Stephan KE, Kasper L, Harrison LM, Daunizeau J, den Ouden HE,
%    Breakspear M, Friston KJ. Nonlinear dynamic causal models for fMRI.
%    Neuroimage 42:649-662, 2008.
% 4. Marreiros AC, Kiebel SJ, Friston KJ. Dynamic causal modelling for
%    fMRI: a two-state model.
%    Neuroimage. 2008 Jan 1;39(1):269-78.
%__________________________________________________________________________

% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston & Klaas Enno Stephan
% $Id: spm_fx_fmri_linear.m 4052 2010-08-27 19:22:44Z karl $

% Neuronal motion
%==========================================================================
P.B   = full(P.B);                       % bi-linear parameters
P.C   = P.C/16;                          % exogenous parameters
P.D   = full(P.D);                       % nonlinear parameters

% excitatory connections
%--------------------------------------------------------------------------
for i = 1:size(P.B,3)
    P.A = P.A + u(i)*P.B(:,:,i);
end

% and nonlinear (state) terms
%--------------------------------------------------------------------------
for i = 1:size(P.D,3)
    P.A = P.A + x(i,1)*P.D(:,:,i);
end

% implement differential state equation y = dx/dt (neuronal)
%--------------------------------------------------------------------------
y    = x;
if size(x,2) == 5

    % one neuronal state per region
    %----------------------------------------------------------------------
    y(:,1) = P.A*x(:,1) + P.C*u(:);

else

    % extrinsic (two neuronal states)
    %----------------------------------------------------------------------
    A      = exp(P.A)/8;             % enforce positivity
    IE     = diag(diag(A));          % inhibitory to excitatory
    EE     = A - IE;                 % excitatory to excitatory
    EI     = 1;                      % excitatory to inhibitory
    SE     = 1;                      % self-inhibition (excitatory)
    SI     = 2;                      % self-inhibition (inhibitory)

    % motion - excitatory and inhibitory: y = dx/dt
    %----------------------------------------------------------------------
    y(:,1) = EE*x(:,1) - SE*x(:,1) - IE*x(:,6) + P.C*u(:);
    y(:,6) = EI*x(:,1) - SI*x(:,6);

end

% Hemodynamic motion
%==========================================================================

% hemodynamic parameters
%--------------------------------------------------------------------------
%   H(1) - signal decay                                   d(ds/dt)/ds)
%   H(2) - autoregulation                                 d(ds/dt)/df)
%   H(3) - transit time                                   (t0)
%   H(4) - exponent for Fout(v)                           (alpha)
%   H(5) - resting oxygen extraction                      (E0)
%   H(6) - ratio of intra- to extra-vascular components   (epsilon)
%          of the gradient echo signal
%--------------------------------------------------------------------------
H      = [0.65 0.41 2.00 0.32 0.34];
H      = [0.64 0.32 2.00 0.32 0.32];

% signal decay
%--------------------------------------------------------------------------
sd     = H(1)*exp(P.decay);

% transit time
%--------------------------------------------------------------------------
tt     = H(3)*exp(P.transit);

% Fout = f(v) - outflow
%--------------------------------------------------------------------------
fv     = 1 + x(:,4)/H(4);

% e = f(f) - oxygen extraction x rCBF
%--------------------------------------------------------------------------
ff     = (1 - log(H(5))).*x(:,3);

% implement differential state equation y = dx/dt (hemodynamic)
%--------------------------------------------------------------------------
y(:,2) = x(:,1) - sd.*x(:,2) - H(2)*x(:,3);
y(:,3) = x(:,2);
y(:,4) = (1  + x(:,3) - fv)./tt;
y(:,5) = (ff - x(:,5) + x(:,4))./tt;
y      = y(:);

author: spm


text:
function [u1, u2] = projectSVD(U1, u)
% Assumes that `U1` has orthonormal columns
% so that :math:`U1 U1^T` is a projector.
% It returns the projections
%
% .. math::
%    u1 &= U1 U1^T u \\
%    u2 &= (I - U1 U1^T) u
%
% USAGE:
%
%    [u1, u2] = projectSVD(U1, u)
%
% EXAMPLE:
%
%    [U1, D1, V1, r] = subspaceSVD(S);
%    [uC, uL] = projectSVD(U1, u);
%    [vfR, vfN] = projectSVD(V1, vf);
%    [vrR, vrN] = projectSVD(V1, vr);
%
% .. Author: - Michael Saunders 29 Jul 2009 First version of projectSVD.m written as alternative to Ronan's projectOntoSubspace.m.

  u1 = U1*(U1'*u);
  u2 = u - u1;


text:
function [ a, b ] = p22_lim ( )

%*****************************************************************************80
%
%% P22_LIM returns the integration limits for problem 22.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A, B, the limits of integration.
%
  a = 0.0;
  b = 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function yp=predict(tb,yb,n,f_unk,B_un,mass);

% functions of the generalized state equations of a Linear System

g(1:n,1)=yb(n+1:2*n);
[stiff,damp,xkp,xcp]=kcm(n,yb);
g(n+1:2*n)=B_un*f_unk+inv(mass)*(-stiff*yb(1:n)-damp*yb(n+1:2*n));
g(2*n+1:3*n+2)=zeros(n+2,1);

yp =g;

author: Sable

repo: mcbench-benchmarks


text:
function movie_sound_field_mono_wfs(X,Y,Z,xs,src,f,outfile,conf)
%MOVIE_SOUND_FIELD_MONO_WFS movie of a WFS sound field
%
%   Usage: movie_sound_field_mono_wfs(X,Y,Z,xs,src,f,outfile,conf)
%
%   Input parameters:
%       X           - x-axis / m; single value or [xmin,xmax]
%       Y           - y-axis / m; single value or [ymin,ymax]
%       Z           - z-axis / m; single value or [zmin,zmax]
%       xs          - position of point source / m
%       src         - sourcetype of the virtual source:
%                         'pw' - plane wave (xs, ys are the direction of the
%                                plane wave in this case)
%                         'ps' - point source
%                         'fs' - focused source
%       f           - monochromatic frequency / Hz
%       outfile     - name for the movie file
%       conf        - configuration struct (see SFS_config)
%
%   MOVIE_SOUND_FIELD_MONO_WFS(X,Y,Z,xs,src,f,L,outfile,conf) generates a
%   movie of simulations of a sound field of the given source positioned at xs
%   using a WFS driving function in the temporal domain with different phase.
%
%   See also: sound_field_mono_wfs, plot_sound_field

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 8;
nargmax = 8;
narginchk(nargmin,nargmax);
isargvector(X,Y,Z);
isargxs(xs);
isargpositivescalar(f);
isargchar(src,outfile);
isargstruct(conf);

%% ===== Configuration ==================================================
% Plotting
useplot = conf.plot.useplot;
% Temporary dir
tmpdir = conf.tmpdir;

%% ===== Simulation =====================================================
phase = linspace(2*pi,0,25);
% Generate a random number string for the tmp files
rn = sprintf('%04.0f',10000*rand);
% Simulate the time by different phase values
for ii = 1:length(phase)-1
    conf.phase = phase(ii);
    conf.plot.useplot = 0;
    % Calculate sound field for the given phase
    [P,x,y,z,x0] = sound_field_mono_wfs(X,Y,Z,xs,src,f,conf);

    % === Save temporary data ===
    if ~exist(tmpdir,'dir')
        mkdir(tmpdir);
    end
    conf.plot.file = sprintf('%s/%s_%04.0f.png',tmpdir,rn,ii);
    plot_sound_field(P,x,y,z,x0,conf);
end

%% ===== Create movie ====================================================
conf.plot.useplot = useplot;
generate_movie(outfile,tmpdir,rn);

% Clean up tmp files
delete([tmpdir,'/',rn,'*.png']);

%% ===== Show movie ======================================================
if useplot
    status = system('which mplayer');
    if status
        error('%s: mplayer is needed to show this movie.',upper(mfilename));
    else
        cmd = sprintf('mplayer %s -loop 0',outfile);
        system(cmd);
    end
end

author: sfstoolbox

repo: sfs-matlab


text:
% Create a spectrum with a linear phase, a given delay in samples
%
% Description
%  If the spectrum length is even, the Nyquist's phase is managed as follow:
%   sign(real(exp((delay*1i*pi))))
%
% Inputs
%  delay    : [samples]
%  fftlen   : length of the fft
%
% Outputs
%  shift    : the delay-spectrum
%
% Example
%  S = fft(...);
%  D = delay2spec(12, 1024);
%  S = S.*D;  % shift the time signal by 12 samples to the left
%
% Copyright (c) 2008 Ircam/CNRS-UMR9912-STMS
%
% License
%  This file is under the LGPL license,  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. This file 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.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Gilles Degottex <gilles.degottex@ircam.fr>
%

function shift = delay2spec(delay, dftlen)
    if delay==0
        shift = ones(1,dftlen);
    else
        % odd length
        if mod(dftlen,2)==1
            shift = exp((delay*2i*pi./dftlen).*(1:(dftlen-1)/2));
            shift = [1, shift(1:end), conj(shift(end:-1:1))];

        % even length
        else
            shift = exp((delay*2i*pi./dftlen).*(1:dftlen/2));
            shift = [1, shift(1:end-1), sign(real(shift(end))), conj(shift(end-1:-1:1))];
        end
    end
return

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
%  Figure 6.60      Feedback Control of Dynamic Systems, 5e
%                   Franklin, Powell, Emami
%

clear all
close all

num=10;
den=conv([1 0],[1/2.5 1]);
den=conv(den,[1/6 1]);
w=logspace(-1,2,100);
[mag,phas]=bode(num,den,w);
[OLgm,OLpm,OLwcg,OLwcp]=margin(mag,phas,w)
%Lead compensator
numl=10*[1 2];
denl=[1 20];
num1=conv(num,numl);
den1=conv(den,denl);
[magcl,phascl]=bode(num1,den1,w);
[D1gm,D1pm,D1wcg,D1wcp]=margin(magcl,phascl,w)
numll=10*conv(numl,[1 4]);
denll=conv(denl,[1 40]);
num2=conv(num,numll);
den2=conv(den,denll);
[magcll,phascll]=bode(num2,den2,w);
[D2gm,D2pm,D2wcg,D2wcp]=margin(magcll,phascll,w)
subplot(2,1,1)
loglog(w,mag,'-',w,magcl,'--',w,magcll,'-.',w,ones(100,1),'-');
axis([.1 100 .01 10])
grid;
ylabel('Magnitude');
title('Fig. 6.60 Bode plot for Example 6.16.   (a) magnitude');
subplot(2,1,2)
semilogx(w,phas,'-',w,phascl,'--',w,phascll,'-.',w,-180*ones(100,1),'-');
axis([.1 100 -250 -50])
grid;
xlabel('\omega (rad/sec)');
ylabel('phase (deg)');
title('Fig. 6.60 (b) phase.');


text:
function p = rad(p);
%RAD          Radius of (interval) polynomial (same as p.rad)
%
%   r = rad(p)
%

% written  10/04/02     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  p.c = rad(p.c);
  p = normalize(p);

repo: OpenMAS


text:
function nameout = cleanvarname(namein)

% Routine to remove unallowed characters from strings
% nameout can be use as a variable or field in a structure

% custom change
if strcmp(namein,'#')
    namein = 'nb';
end

% 1st char must be a letter
for l=1:length(namein)
    if isletter(namein(l))
        nameout = namein(l:end);
        break % exist for loop
    end
end

% remove usual suspects
stringcheck = [strfind(namein,'('), ...
    strfind(namein,')') ...
    strfind(namein,'[') ...
    strfind(namein,']') ...
    strfind(namein,'{') ...
    strfind(namein,'}') ...
    strfind(namein,'-') ...
    strfind(namein,'+') ...
    strfind(namein,'*') ...
    strfind(namein,'/') ...
    strfind(namein,'#') ...
    strfind(namein,'%') ...
    strfind(namein,'&') ...
    strfind(namein,'@') ...
    ];

if ~isempty(stringcheck)
    nameout(stringcheck) = [];
end

% last check
if ~isvarname(nameout)
    error('the variable name to use is still invalid, check chars to remove')
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
%KWAVETRANSDUCER   Class definition for k-Wave transducer
%
% DESCRIPTION:
%       See makeTransducer for function arguments.
%
% ABOUT:
%       author          - Bradley Treeby
%       date            - 9th December 2010
%       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

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

% Define as a handle class - A handle class constructor returns a handle
% object that is a reference to the object created. You can assign the
% handle object to multiple variables or pass it to functions without
% causing MATLAB to make a copy of the original object. A function that
% modifies a handle object passed as an input argument does not need to
% return the object.
classdef kWaveTransducer < handle

    % define the properties of the transducer that cannot be modified by
    % the user after the transducer is initialised (these parameters are
    % stored). The numbers assigned here are the default values used if the
    % parameters are not set explicitly by the user.
    properties (GetAccess = 'public', SetAccess = 'private')

        % the total number of transducer elements
        number_elements = 128;

        % the width of each element in grid points
        element_width = 1;

        % the length of each element in grid points
        element_length = 20;

        % the spacing (kerf width) between the transducer elements in grid
        % points
        element_spacing = 0;

        % the position of the corner of the transducer in the grid
        position = [1, 1, 1];

        % the radius of curvature of the transducer [m]
        radius = inf;

    end

    % define the properties of the transducer that can be modified by the
    % user after the transducer is initialised (these parameters are
    % stored). The numbers assigned here are the default values used if the
    % parameters are not set explicitly by the user.
    properties (Access = 'public')

        % the transducer elements that are currently active elements
        active_elements;

        % the focus depth in the elevation direction [m]
        elevation_focus_distance = inf;

        % transmit apodization
        transmit_apodization = 'Rectangular';

        % receive apodization
        receive_apodization = 'Rectangular';

        % sound speed used to calculate beamforming delays [m/s]
        sound_speed = 1540;

        % focus distance used to calculate beamforming delays [m]
        focus_distance = inf;

        % time index to start recording if transducer is used as a sensor
        record_start_index = 1;

    end

    % define the dependent properties (these parameters are computed when
    % queried). Some of these properties act like normal user-modifiable
    % properties (e.g., steering_angle), but instead have explicit set and
    % get functions that change an internal stored value. This allows error
    % checking when the properties are changed.
    properties(Dependent = true)

        % binary mask of the active transducer elements
        active_elements_mask;

        % indexed mask of the active transducer elements, where the index
        % indicates the transducer element that each grid point belongs to.
        indexed_active_elements_mask;

        % binary mask of all the transducer elements (both active and
        % inactive)
        all_elements_mask;

        % indexed mask of all the transducer elements, where the index
        % indicates the transducer element that each grid point belongs to.
        indexed_elements_mask;

        % pseudonym for the active_elements_mask to maintain compatability
        % when the transducer is used in place of sensor or source.
        mask;

        % total width of the transducer in grid points
        transducer_width;

        % current number of active transducer elements
        number_active_elements;

        % user defined input signal with appended zeros dependent on the
        % beamforming (steering_angle, focus_distance) settings
        input_signal;

        % maximum steering angle - this sets the beamforming delay offset
        steering_angle_max;

        % number of zeros appended to input signal to allow a single time
        % series to be used within kspaceFirstOrder3D (either set to 'auto'
        % or based on the setting for steering_angle_max)
        appended_zeros;

        % offset used to make all the delays in the delay_mask positive
        % (either set to 'auto' or based on the setting for
        % steering_angle_max)
        beamforming_delays_offset;

        % steering angle used to calculate beamforming delays [deg]
        steering_angle;

    end

    % define the hidden properties - these cannot be seen or accessed
    % directly by the user, and are used to store internal properties
    properties (Hidden = true, Access = 'private')

        % indexed sensor mask from which the other masks are computed
        indexed_mask;

        % indexed mask for the voxels in each element used to compute
        % elevation focussing
        indexed_element_voxel_mask;

        % size of the grid in which the transducer is defined
        stored_grid_size;

        % corresponding grid spacing
        grid_spacing;

        % time spacing
        dt;

        % original copy of the user defined input signal
        stored_input_signal;

        % stored value of appended_zeros (accessed using get and
        % set methods). This is used to set the number of zeros that are
        % appended and prepended to the input signal.
        stored_appended_zeros = 'auto';

        % stored value of the minimum beamforming delay. This is used to
        % offset the delay mask so that all the delays are >= 0
        stored_beamforming_delays_offset = 'auto';

        % stored value of the steering_angle_max (accessed using get and
        % set methods). This can be set by the user and is used to derive
        % the two parameters above.
        stored_steering_angle_max = 'auto';

        % stored value of the steering_angle (accessed using get and
        % set methods)
        stored_steering_angle = 0;

    end

    % constructor function
    methods
        function transducer = kWaveTransducer(kgrid, transducer_properties)

            % allocate the grid size and spacing
            transducer.stored_grid_size = [kgrid.Nx, kgrid.Ny, kgrid.Nz];
            transducer.grid_spacing = [kgrid.dx, kgrid.dy, kgrid.dz];

            % allocate the temporal spacing
            if isnumeric(kgrid.dt)
                transducer.dt = kgrid.dt;
            elseif isnumeric(kgrid.t_array)
                transducer.dt = kgrid.t_array(2) - kgrid.t_array(1);
            else
                error('kgrid.dt or kgrid.t_array must be explicitly defined');
            end

            % check for properties input
            if nargin == 1
                transducer_properties = [];
            end

            % check the input fields
            checkFieldNames(transducer_properties, {...
                'number_elements', 'element_width', 'element_length', 'element_spacing',...
                'position', 'radius', ...
                'active_elements', 'elevation_focus_distance', 'focus_distance', 'steering_angle',...
                'receive_apodization', 'transmit_apodization',...
                'sound_speed', 'input_signal', 'steering_angle_max'});

            % replace default settings with user defined properties if
            % given
            % -------------------------------------------------------------

            if isfield(transducer_properties, 'number_elements')
                % force value to be a positive integer
                transducer.number_elements = round(abs(transducer_properties.number_elements));
            end

            if isfield(transducer_properties, 'element_width')
                % force value to be a positive integer
                transducer.element_width = round(abs(transducer_properties.element_width));
            end

            if isfield(transducer_properties, 'element_length')
                % force value to be a positive integer
                transducer.element_length = round(abs(transducer_properties.element_length));
            end

            if isfield(transducer_properties, 'element_spacing')
                % force value to be a positive integer
                transducer.element_spacing = round(abs(transducer_properties.element_spacing));
            end

            if isfield(transducer_properties, 'position')
                % force values to be positive integers
                transducer.position = round(abs(transducer_properties.position));
            end

            if isfield(transducer_properties, 'radius')
                transducer.radius = transducer_properties.radius;

                % only allow an infinite radius for now
                if ~isinf(transducer.radius)
                    error('Only a value of transducer.radius = inf is currently supported');
                end
            end

            if isfield(transducer_properties, 'active_elements')
                transducer.active_elements = transducer_properties.active_elements;
            else
                transducer.active_elements = ones(transducer_properties.number_elements, 1);
            end

            if isfield(transducer_properties, 'elevation_focus_distance')
                transducer.elevation_focus_distance = transducer_properties.elevation_focus_distance;
            end

            if isfield(transducer_properties, 'receive_apodization')
                % if a user defined apodization is given explicitly, check
                % the length of the input
                if isnumeric(transducer_properties.receive_apodization) ...
                        && (length(transducer_properties.receive_apodization) ~= transducer.number_active_elements)
                    error('The length of the receive apodization input must match the number of active elements');
                end

                % assign the input
                transducer.receive_apodization = transducer_properties.receive_apodization;
            end

            if isfield(transducer_properties, 'transmit_apodization')
                % if a user defined apodization is given explicitly, check
                % the length of the input
                if isnumeric(transducer_properties.transmit_apodization) ...
                        && (length(transducer_properties.transmit_apodization) ~= transducer.number_active_elements)
                    error('The length of the transmit apodization input must match the number of active elements');
                end

                % assign the input
                transducer.transmit_apodization = transducer_properties.transmit_apodization;
            end

            if isfield(transducer_properties, 'sound_speed')
                transducer.sound_speed = transducer_properties.sound_speed;

                % check to see the sound_speed is positive
                if ~(transducer.sound_speed> 0)
                    error('transducer.sound_speed must be greater than 0');
                end
            end

            if isfield(transducer_properties, 'focus_distance')
                transducer.focus_distance = transducer_properties.focus_distance;
            end

            if isfield(transducer_properties, 'input_signal')
                transducer.stored_input_signal = transducer_properties.input_signal;

                % force the input signal to be a column vector
                if numDim(transducer.input_signal) == 1
                    transducer.stored_input_signal = reshape(transducer.stored_input_signal, [], 1);
                else
                    error('transducer.input_signal must be a one-dimensional array');
                end
            end

            if isfield(transducer_properties, 'steering_angle_max')
                % set the maximum steering angle using the set method (this
                % avoids having to duplicate error checking)
                transducer.steering_angle_max = transducer_properties.steering_angle_max;
            end

            if isfield(transducer_properties, 'steering_angle')
                % set the maximum steering angle using the set method (this
                % avoids having to duplicate error checking)
                transducer.steering_angle = transducer_properties.steering_angle;
            end

            % -------------------------------------------------------------

            % check the transducer fits into the grid
            if sum(transducer.position(:) == 0)
                error('The defined transducer must be positioned within the grid');
            elseif (transducer.position(2) + transducer.number_elements*transducer.element_width + (transducer.number_elements - 1)*transducer.element_spacing) > transducer.stored_grid_size(2)
                error('The defined transducer is too large or positioned outside the grid in the y-direction');
            elseif (transducer.position(3) + transducer.element_length) > transducer.stored_grid_size(3)
                transducer.position(3)
                transducer.element_length
                transducer.stored_grid_size(3)
                error('The defined transducer is too large or positioned outside the grid in the z-direction');
            elseif transducer.position(1) > transducer.stored_grid_size(1)
                error('The defined transducer is positioned outside the grid in the x-direction');
            end

            % assign the data type for the transducer matrix based on the
            % number of different elements (uint8 supports 255 numbers so
            % most of the time this data type will be used)
            if transducer.number_elements < intmax('uint8');
                mask_type = 'uint8';
            elseif transducer.number_elements < intmax('uint16');
                mask_type = 'uint16';
            elseif transducer.number_elements < intmax('uint32');
                mask_type = 'uint32';
            else
                mask_type = 'uint64';
            end

            % create an empty transducer mask (the grid points within
            % element 'n' are all given the value 'n')
            transducer.indexed_mask = zeros(transducer.stored_grid_size, mask_type);

            % create a second empty mask used for the elevation beamforming
            % delays (the grid points across each element are numbered 1 to
            % M, where M is the number of grid points in the elevation
            % direction)
            transducer.indexed_element_voxel_mask = zeros(transducer.stored_grid_size, mask_type);

            % create the corresponding indexing variable 1:M
            element_voxel_index = repmat(1:transducer.element_length, transducer.element_width, 1);

            % for each transducer element, calculate the grid point indices
            for element_index = 1:transducer.number_elements

                % assign the current transducer position
                element_pos_x = transducer.position(1);
                element_pos_y = transducer.position(2) + (transducer.element_width + transducer.element_spacing)*(element_index - 1);
                element_pos_z = transducer.position(3);

                % assign the grid points within the current element the
                % index of the element
                transducer.indexed_mask(element_pos_x, element_pos_y:element_pos_y + transducer.element_width - 1, element_pos_z:element_pos_z + transducer.element_length - 1) = element_index;

                % assign the individual grid points an index corresponding
                % to their order across the element
                transducer.indexed_element_voxel_mask(element_pos_x, element_pos_y:element_pos_y + transducer.element_width - 1, element_pos_z:element_pos_z + transducer.element_length - 1) = element_voxel_index;

            end

            % double check the transducer fits within the desired grid size
            if any(size(transducer.indexed_elements_mask) ~= transducer.stored_grid_size)
                error('Desired transducer is larger than the input grid_size');
            end

        end
    end

    % set and get functions for dependent variables
    methods

        % return input signal
        function signal = get.input_signal(obj)

            % copy the stored input signal
            signal = obj.stored_input_signal;

            % check the signal is not empty
            if isempty(signal)
                error('Transducer input signal is not defined');
            end

            % automatically prepend and append zeros if the beamforming
            % delay offset is set

            % check if the beamforming delay offset is set. If so, use this
            % number to prepend and append this number of zeros to the
            % input signal. Otherwise, calculate how many zeros are needed
            % and prepend and append these.
            if ~strcmp(obj.stored_appended_zeros, 'auto')

                % use the current value of the beamforming offset to add
                % zeros to the input signal
                signal = [zeros(obj.stored_appended_zeros, 1); signal; zeros(obj.stored_appended_zeros, 1)];

            else

                % get the current delay beam forming
                delay_mask = obj.delay_mask;

                % find the maximum delay
                delay_max = max(delay_mask(:));

                % count the number of leading zeros in the input signal
                leading_zeros = find(signal ~= 0, 1) - 1;

                % count the number of trailing zeros in the input signal
                trailing_zeros = find(flipud(signal) ~= 0, 1) - 1;

                % check the number of leading zeros is sufficient given the
                % maximum delay
                if leading_zeros < delay_max + 1
                    disp(['  prepending transducer.input_signal with ' num2str(delay_max - leading_zeros + 1) ' leading zeros']);

                    % prepend extra leading zeros
                    signal = [zeros(delay_max - leading_zeros + 1, 1); signal];
                end

                % check the number of leading zeros is sufficient given the
                % maximum delay
                if trailing_zeros < delay_max + 1
                    disp(['  appending transducer.input_signal with ' num2str(delay_max - trailing_zeros + 1) ' trailing zeros']);

                    % append extra trailing zeros
                    signal = [signal; zeros(delay_max - trailing_zeros + 1, 1)];
                end
            end

        end

        function transducer_width = get.transducer_width(obj)
            % return the overall length of the transducer
            transducer_width = obj.number_elements*obj.element_width + (obj.number_elements - 1)*obj.element_spacing;
        end

        function number_active_elements = get.number_active_elements(obj)
            % return the current number of active elements
            number_active_elements = sum(obj.active_elements(:));
        end

        % allow mask query to allow compatability with regular sensor
        % structure - return the active sensor mask
        function mask = get.mask(obj)
            mask = obj.active_elements_mask;
        end

        % return a binary mask showing the locations of the active elements
        function mask = get.active_elements_mask(obj)

            % copy the indexed elements mask
            mask = obj.indexed_mask;

            % remove inactive elements from the mask
            for element_index = 1:obj.number_elements
                mask(mask == element_index) = obj.active_elements(element_index);
            end

            % convert remaining mask to binary
            mask(mask ~= 0) = 1;

        end

        % return a binary mask showing the locations of all the elements
        % (both active and inactive)
        function mask = get.all_elements_mask(obj)
            mask = obj.indexed_mask;
            mask(mask ~= 0) = 1;
        end

        % return a copy of the indexed elements mask
        function mask = get.indexed_elements_mask(obj)
            mask = obj.indexed_mask;
        end

        % return a copy of the indexed elements mask only for the active
        % elements
        function mask = get.indexed_active_elements_mask(obj)

            % copy the indexed elements mask
            mask = obj.indexed_mask;

            % remove inactive elements from the mask
            for element_index = 1:obj.number_elements
                if ~(obj.active_elements(element_index))
                    mask(mask == element_index) = 0;
                end
            end

            % force the lowest element index to be 1
            lowest_active_element_index = find(obj.active_elements, 1);
            mask(mask ~= 0) = mask(mask ~= 0) - lowest_active_element_index + 1;

        end

        % return the stored value of the steering angle
        function steering_angle = get.steering_angle(obj)
            steering_angle = obj.stored_steering_angle;
        end

        % set the stored value of the steering angle
        function set.steering_angle(obj, steering_angle)

            % force to be scalar
            steering_angle = steering_angle(1);

            % check if the steering angle is between -90 and 90
            if ~(steering_angle < 90 && steering_angle > -90)
                error('Input for steering_angle must be betweeb -90 and 90 degrees');
            end

            % check if the steering angle is less than the maximum steering
            % angle
            if ~strcmp(obj.stored_steering_angle_max, 'auto') && (abs(steering_angle) > obj.stored_steering_angle_max)
                error('Input for steering_angle cannot be greater than steering_angle_max');
            end

            % update the stored value
            obj.stored_steering_angle = steering_angle;

        end

        % return the stored value of the maximum steering angle
        function steering_angle_max = get.steering_angle_max(obj)
            steering_angle_max = obj.stored_steering_angle_max;
        end

        % set the stored value of the maximum steering angle and
        % beamforming offset
        function set.steering_angle_max(obj, steering_angle_max)

            % force input to be scalar and positive
            steering_angle_max = steering_angle_max(1);

            % check the steering angle is within range
            if ~(steering_angle_max < 90 && steering_angle_max > -90)
                error('Input for steering_angle_max must be between -90 and 90');
            end

            % check the maximum steering angle is greater than the current
            % steering angle
            if ~strcmp(obj.stored_steering_angle_max, 'auto') && (abs(obj.stored_steering_angle) > steering_angle_max)
                error('Input for steering_angle_max cannot be less than the current steering_angle');
            end

            % overwrite the stored value
            obj.stored_steering_angle_max = steering_angle_max;

            % store a copy of the current value for the steering angle
            current_steering_angle = obj.stored_steering_angle;

            % overwrite with the user defined maximum value
            obj.stored_steering_angle = steering_angle_max;

            % set the beamforming delay offset to zero (this means the
            % delays will contain negative values which we can use to
            % derive the new values for the offset)
            obj.stored_appended_zeros = 0;
            obj.stored_beamforming_delays_offset = 0;

            % get the element beamforming delays and reverse
            delay_times = -obj.beamforming_delays;

            % get the maximum and minimum beamforming delays
            max_delay = max(delay_times);
            min_delay = min(delay_times);

            % add the maximum and minimum elevation delay if the elevation
            % focus is not set to infinite
            if ~isinf(obj.elevation_focus_distance)
                max_delay = max_delay + max(obj.elevation_beamforming_delays);
                min_delay = min_delay + min(obj.elevation_beamforming_delays);
            end

            % set the beamforming offset to the difference between the
            % maximum and minimum delay
            obj.stored_appended_zeros = max_delay - min_delay;

            % set the minimum beamforming delay (converting to a positive
            % number)
            obj.stored_beamforming_delays_offset = - min_delay;

            % reset the previous value of the steering angle
            obj.stored_steering_angle = current_steering_angle;

        end

        % return the stored value of appended zeros
        function appended_zeros = get.appended_zeros(obj)
            appended_zeros = obj.stored_appended_zeros;
        end

        % return the stored value of the offset used to force the values in
        % delay_mask to be >= 0
        function beamforming_delays_offset = get.beamforming_delays_offset(obj)
            beamforming_delays_offset = obj.stored_beamforming_delays_offset;
        end

    end

    % general class methods
    methods

        % calculate the beamforming delays based on the focus and steering
        % settings
        function delay_times = beamforming_delays(obj)

            % calculate the element pitch in [m]
            element_pitch = (obj.element_spacing + obj.element_width)*obj.grid_spacing(2);

            % create indexing variable
            element_index = -(obj.number_active_elements - 1)/2:(obj.number_active_elements - 1)/2;

            % check for focus depth
            if isinf(obj.focus_distance)
                % calculate time delays for a steered beam
                delay_times = element_pitch*element_index*sin(obj.steering_angle*pi/180)/(obj.sound_speed);        % [s]
            else
                % calculate time delays for a steered and focussed beam
                delay_times = obj.focus_distance/obj.sound_speed * (1 - sqrt(...
                    1 ...
                    + (element_index*element_pitch./obj.focus_distance).^2 ...
                    - 2 * (element_index*element_pitch./obj.focus_distance) * sin(obj.steering_angle*pi/180) ...
                    ));      % [s]
            end

            % convert the delays to be in units of time points
            delay_times = round(delay_times/obj.dt);

        end

        % calculate the elevation beamforming delays based on the focus
        % setting
        function delay_times = elevation_beamforming_delays(obj)

            if ~isinf(obj.elevation_focus_distance)

                % create indexing variable
                element_index = -(obj.element_length - 1)/2:(obj.element_length - 1)/2;

                % calculate time delays for a focussed beam
                delay_times = (obj.elevation_focus_distance - sqrt( (element_index*obj.grid_spacing(3)).^2 + obj.elevation_focus_distance^2 ))./obj.sound_speed;

                % convert the delays to be in units of time points and then
                % reverse
                delay_times = -round(delay_times/obj.dt);

            else

                % create an empty array
                delay_times = zeros(1, obj.element_length);

            end

        end

        % function to create a scan line based on the input sensor data and
        % the current apodization and beamforming setting
        function line = scan_line(obj, sensor_data)

            % get the current apodization setting
            apodization = obj.get_receive_apodization;

            % get the current beamforming weights and reverse
            delays = -obj.beamforming_delays;

            % offset the received sensor_data by the beamforming delays and
            % apply receive apodization
            for element_index = 1:obj.number_active_elements
                if delays(element_index) > 0
                    % shift element data forwards
                    sensor_data(element_index, :) = apodization(element_index).*[sensor_data(element_index, 1 + delays(element_index):end), zeros(1, delays(element_index))];
                elseif delays(element_index) < 0
                    % shift element data backwards
                    sensor_data(element_index, :) = apodization(element_index).*[zeros(1, -delays(element_index)), sensor_data(element_index, 1:end + delays(element_index))];
                end
            end

            % form the a-line summing across the elements
            line = sum(sensor_data);

        end

        % function to average the individual time series recorded at each
        % grid point of each transducer element (as returned by the C++
        % code) and return a single time series per active transducer
        % element (as returned by the MATLAB code)
        function sensor_data_sum = combine_sensor_data(obj, sensor_data)

            % check the data is the correct size
            if size(sensor_data, 1) ~= (obj.number_active_elements * obj.element_width * obj.element_length)
                error('The number of time series in the input sensor_data must match the number of grid points in the active tranducer elements.');
            end

            % get index of which element each time series belongs to
            ind = obj.indexed_active_elements_mask(obj.indexed_active_elements_mask > 0);

            % create empty output
            sensor_data_sum = zeros(obj.number_active_elements, size(sensor_data, 2));

            % check if an elevation focus is set
            if isinf(obj.elevation_focus_distance)

                % loop over time series and sum
                for ii = 1:length(ind)
                    sensor_data_sum(ind(ii), :) = sensor_data_sum(ind(ii), :) + sensor_data(ii, :);
                end

            else

                % get the elevation delay for each grid point of each
                % active transducer element (this is given in units of grid
                % points)
                dm = obj.delay_mask(2);
                dm = dm(obj.active_elements_mask ~= 0);

                % loop over time series, shift and sum
                for ii = 1:length(ind)
                    sensor_data_sum(ind(ii), 1:end - dm(ii)) = sensor_data_sum(ind(ii), 1:end - dm(ii)) + sensor_data(ii, 1 + dm(ii):end);
                end

            end

            % divide by number of time series in each element
            sensor_data_sum = sensor_data_sum .* (1/ (obj.element_width * obj.element_length));

        end

        % plot the transducer using voxelPlot
        function plot(obj)
            voxelPlot(double(obj.all_elements_mask));
        end

        % print out the transducer properties
        function properties(obj)
            disp(' ');
            disp('k-Wave Transducer Properties');
            disp(['  transducer position:       [' num2str(obj.position) ']']);
            disp(['  transducer width:          ' scaleSI(obj.transducer_width*obj.grid_spacing(2)) 'm (' num2str(obj.transducer_width) ' grid points)']);
            disp(['  number of elements:        ' num2str(obj.number_elements)]);
            disp(['  number of active elements: ' num2str(obj.number_active_elements) ' (elements ' num2str(find(obj.active_elements, 1, 'first')) ' to ' num2str(find(obj.active_elements, 1, 'last')) ')']);
            disp(['  element width:             ' scaleSI(obj.element_width*obj.grid_spacing(2)) 'm (' num2str(obj.element_width) ' grid points)']);
            disp(['  element spacing (kerf):    ' scaleSI(obj.element_spacing*obj.grid_spacing(2)) 'm (' num2str(obj.element_spacing) ' grid points)']);
            disp(['  element pitch:             ' scaleSI((obj.element_spacing + obj.element_width)*obj.grid_spacing(2)) 'm (' num2str((obj.element_spacing + obj.element_width)) ' grid points)']);
            disp(['  element length:            ' scaleSI(obj.element_length*obj.grid_spacing(3)) 'm (' num2str(obj.element_length) ' grid points)']);
            disp(['  sound speed:               ' num2str(obj.sound_speed) 'm/s']);
            if isinf(obj.focus_distance)
                disp('  focus distance:             infinite');
            else
                disp(['  focus distance:            ' scaleSI(obj.focus_distance) 'm']);
            end
            if isinf(obj.elevation_focus_distance)
                disp('  elevation focus distance:  infinite');
            else
                disp(['  elevation focus distance:  ' scaleSI(obj.elevation_focus_distance) 'm']);
            end
            disp(['  steering angle:            ' num2str(obj.steering_angle) ' degrees']);
            disp(['  steering angle max:        ' num2str(obj.steering_angle_max) ' degrees']);
        end

    end

    % hidden class methods used with kspaceFirstOrder3D
    methods (Hidden = true, Access = 'public')

        % return a delay mask based on the transmit beamforming delays.
        % This is used if the transducer is a source. The delays are set
        % for each grid point of the active_elements_mask and define which
        % element of the input_signal should be used at each source point.
        % This mask is incremented within kspaceFirstOrder3D after each
        % time step. The mask is indexed from 0 and is incremented by 1
        % within kspaceFirstOrder3D .
        function mask = delay_mask(obj, mode)
            % mode == 1: both delays
            % mode == 2: elevation only
            % mode == 3: azimuth only

            % assign the delays to a new mask using the indexed_element_mask
            indexed_active_elements_mask = obj.indexed_active_elements_mask;
            active_elements_index = find(indexed_active_elements_mask ~= 0);
            mask = zeros(obj.stored_grid_size);

            % calculate azimuth focus delay times provided they are not all
            % zero
            if (~isinf(obj.focus_distance) || (obj.steering_angle ~= 0)) && (nargin == 1 || mode ~= 2)

                % get the element beamforming delays and reverse
                delay_times = -obj.beamforming_delays;

                % add delay times
                mask(active_elements_index) = delay_times(indexed_active_elements_mask(active_elements_index));

            end

            % calculate elevation focus time delays provided each element
            % is longer than one grid point
            if ~isinf(obj.elevation_focus_distance) && (obj.element_length > 1) && (nargin == 1 || mode ~= 3)

                % get elevation beamforming delays
                elevation_delay_times = obj.elevation_beamforming_delays;

                % get current mask
                element_voxel_mask = obj.indexed_element_voxel_mask;

                % add delay times
                mask(active_elements_index) = mask(active_elements_index) + elevation_delay_times(element_voxel_mask(active_elements_index)).';

            end

            % shift delay times (these should all be >= 0, where a value of
            % 0 means no delay)
            if strcmp(obj.stored_appended_zeros, 'auto')
                mask(active_elements_index) = mask(active_elements_index) - min(min(min(mask(active_elements_index))));
            else
                mask(active_elements_index) = mask(active_elements_index) + obj.stored_beamforming_delays_offset;
            end

        end

        % return a binary mask based on the elevation focus. This sets the
        % elements to extract at each timestep from the sensor_data_buffer
        % used in kspaceFirstOrder3D. This is used if the transducer is a
        % sensor and the elevation focus is set. The mask is of the form:
        %
        %   1 0 0 0
        %   0 1 0 0
        %   0 0 1 0
        %   0 0 1 0
        %   0 0 0 1
        %   0 0 0 1
        %   0 0 0 1
        %   0 0 1 0
        %   0 0 1 0
        %   0 1 0 0
        %   1 0 0 0
        %
        % The mask remains constant, but extracts new elements at each time
        % step as the FIFO buffer sensor_data_buffer is updated
        function mask = elevation_beamforming_mask(obj)

            % get elevation beamforming mask
            delay_mask = obj.delay_mask(2);

            % extract the active elements
            delay_mask = delay_mask(obj.active_elements_mask ~= 0);

            % force delays to start from zero
            delay_mask = delay_mask - min(delay_mask(:));

            % create an empty output mask
            mask = zeros(length(delay_mask), max(delay_mask) + 1);

            % populate the mask by setting 1's at the index given by the
            % delay time
            for index = 1:length(delay_mask)
                mask(index, delay_mask(index) + 1) = 1;
            end

            % flip the mask so the shortest delays are at the right
            mask = fliplr(mask);
        end

        % convert the transmit apodization into the form of a element mask,
        % where the apodization values are placed at the grid points
        % belonging to the active transducer elements. These values are
        % then extracted in the correct order within
        % kspaceFirstOrder_inputChecking using apodization =
        % transmit_apodization_mask(active_elements_mask ~= 0)
        function mask = transmit_apodization_mask(obj)

            % check if a user defined apodization is given and whether this
            % is still the correct size (in case the number of active
            % elements has changed)
            if isnumeric(obj.transmit_apodization)
                if (length(obj.transmit_apodization) ~= obj.number_active_elements)
                    error('The length of the transmit apodization input must match the number of active elements');
                else
                    % assign apodization
                    apodization = obj.transmit_apodization;
                end
            else
                % create apodization using getWin
                apodization = getWin(obj.number_active_elements, obj.transmit_apodization);
            end

            % create an empty mask;
            mask = zeros(obj.stored_grid_size);

            % assign the apodization values to every grid point in the
            % transducer
            mask_index = obj.indexed_active_elements_mask;
            mask_index = mask_index(mask_index ~= 0);
            mask(obj.active_elements_mask == 1) = apodization(mask_index);

        end

        % allow the mask to be resized by the simulation functions (this is
        % used when 'PMLInside' is set to false)
        function expand_grid(obj, expand_size)
            obj.indexed_mask = expandMatrix(obj.indexed_mask, expand_size, 0);
        end

        % allow the mask to be resized by the simulation functions (this is
        % used when 'PMLInside' is set to false)
        function retract_grid(obj, retract_size)
            obj.indexed_mask = obj.indexed_mask(1 + retract_size(1):end - retract_size(1), 1 + retract_size(2):end - retract_size(2), 1 + retract_size(3):end - retract_size(3));
        end

        % return the size of the grid
        function grid_size = grid_size(obj)
            grid_size = obj.stored_grid_size;
        end

    end

    % internal class methods only accessible by other functions within
    % kWaveTransducer
    methods (Hidden = true, Access = 'private')

        % return the transmit apodization, converting strings of window
        % type to actual numbers using getWin
        function apodization = get_transmit_apodization(obj)

            % check if a user defined apodization is given and whether this
            % is still the correct size (in case the number of active
            % elements has changed)
            if isnumeric(obj.transmit_apodization)
                if (length(obj.transmit_apodization) ~= obj.number_active_elements)
                    error('The length of the transmit apodization input must match the number of active elements');
                else
                    % assign apodization
                    apodization = obj.transmit_apodization;
                end
            else
                % create apodization using getWin
                apodization = getWin(obj.number_active_elements, obj.transmit_apodization);
            end
        end

        % return the receive apodization, converting strings of window
        % type to actual numbers using getWin
        function apodization = get_receive_apodization(obj)

            % check if a user defined apodization is given and whether this
            % is still the correct size (in case the number of active
            % elements has changed)
            if isnumeric(obj.receive_apodization)
                if (length(obj.receive_apodization) ~= obj.number_active_elements)
                    error('The length of the receive apodization input must match the number of active elements');
                else
                    % assign apodization
                    apodization = obj.receive_apodization;
                end
            else
                % create apodization using getWin
                apodization = getWin(obj.number_active_elements, obj.receive_apodization);
            end
        end

    end
end

repo: Aero-matlab


text:
function x = project_exp(v)
% PROJECT_EXP    Project points onto the exponential cone.
%
%   When v is a 1 x 3 vector, project_exp(v) is the projection
%   of v onto the exponential cone. When v is an n x 3 vector,
%   project_exp(v) projects each row of v onto the cone
%   in a vectorized fashion.
%
%   For reference, the exponential cone and its dual are given by
%     Kexp   = { (x,y,z) | ye^(x/y) <= z, y > 0 }
%     Kexp^* = { (u,v,w) | u < 0, -ue^(v/u) <= ew } cup { (0,v,w) | v,w >= 0 }

    r = v(:,1); s = v(:,2); t = v(:,3);
    x = nan(size(v));

    % v in cl(Kexp)
    idx = ( (s.*exp(r./s) <= t & s > 0) | (r <= 0 & s == 0 & t >= 0) );
    x(idx,:) = v(idx,:);

    % -v in Kexp^*
    idx = ( (-r < 0 & r.*exp(s./r) <= -exp(1).*t) | (-r == 0 & -s >= 0 & -t >= 0) );
    x(idx,:) = 0;

    % special case with analytical solution
    idx = (r < 0 & s < 0);
    x(idx,:) = v(idx,:);
    x(idx,2) = max(x(idx,2),0);
    x(idx,3) = max(x(idx,3),0);

    % minimize ||x - v||^2 subject to se^{r/s} = t via primal-dual Newton method
    % these components are computed serially, so much slower
    idx = find(isnan(x(:,1)));

    g     = @(w)   w(2)*exp(w(1)/w(2)) - w(3);
    gradg = @(w)   [ exp(w(1)/w(2)); exp(w(1)/w(2))*(1 - w(1)/w(2)); -1 ];

    alpha = 0.001; beta = 0.5;

    for i = 1:length(idx)
        disp('newton')
        u = v(idx(i),:)';
        u(2) = max(u(2),1); u(3) = max(u(3),1);
        y = 1; % dual variable

        r = @(w,z) [ w - v(idx(i),:)' + z*gradg(w); g(w) ];

        for iter = 1:100
            KKT = [ eye(3)+y*hessg(u), gradg(u) ; gradg(u)', 0 ];
            z = KKT \ -r(u,y);
            du = z(1:3);
            dy = z(4);

            % backtracking line search
            t = 1;
            ustep = u + t*du; ystep = y + t*dy;
            while ustep(2) < 0 || (norm(r(ustep, ystep)) > (1 - alpha*t)*norm(r(u, y)))
                t = beta*t;
                ustep = u + t*du; ystep = y + t*dy;
            end

            u = ustep; y = ystep;

            if abs(g(u)) < 1e-8 && norm(r(u,y)) <= 1e-8
                x(idx(i),:) = u;
                break;
            end
        end
    end
end

function h = hessg(w)
    r = w(1); s = w(2); t = w(2);
    h = exp(r/s)*[ 1/s,    -r/s^2,   0;
                   -r/s^2, r^2/s^3,  0;
                   0,      0,        0 ];
end

author: cvxgrp

repo: proximal

save_path: github-repos/MATLAB/cvxgrp-proximal


text:
%% Autoregressive Conditional Mean, Variance and Kurtosis
% Allows the estimation of the Autoregressive Conditional Kurtosis Model
% presented in Brooks, C., Burke, S., P., and Persand, G., (2005),
% "Autoregressive Conditional Kurtosis Model", Journal of Financial
% Econometrics, 3(3),339-421.
%
%% *_Mean Models_*
%
% $$ARMAX(AR, MA, X): r_t = a_0 + {\sum_{i=1}^n}{a_1}{r_{t-i}} + {\sum_{j=1}^k}{a_2}{\varepsilon}_{t-j} + {\sum_{l=1}^m}{a_3}{X_l} + {\varepsilon}_t$
%
%% *_Variance Models_*
%
% $$GARCH(P,Q,Y): {\sigma}_t^2 = b_0 + {\sum_{i=1}^p}b_{1,i}{\varepsilon}_{t-i}^2 + {\sum_{j=1}^q}b_{2,j}{\sigma}_{t-q}^2 + {\sum_{l=1}^m}{b_3}{Y_l}$
% $$GJR-GARCH(P,Q,Y): {\sigma}_t^2 = b_0 + {\sum_{i=1}^p}b_{1,i}{\varepsilon}_{t-i}^2 + {\sum_{j=1}^q}b_{2,j}{\sigma}_{t-j}^2 + {\sum_{i=1}^p}b_{3,i}{\varepsilon}_{t-i}^2*I_{t-i} + {\sum_{l=1}^m}{b_3}{Y_l}$
% $$AGARCH(P,Q,Y): {\sigma}_t^2 = b_0 + {\sum_{i=1}^p}b_{1,i}({\varepsilon}_{t-i} + {\gamma_{t-p}}))^2 + {\sum_{j=1}^q}b_{2,j}{\sigma}_{t-j} + {\sum_{l=1}^m}{b_3}{Y_l}$
% $$NAGARCH(P,Q,Y): {\sigma}_t^2 = b_0 + {\sum_{i=1}^p}b_{1,i}({\varepsilon}(t-i)/{\sqrt{{\sigma}_{t-i}^2}} + {\sum_{i=1}^p}{\gamma_{t-i}}^2 +  {\sum_{j=1}^q}b_{2,j}{\sigma}_{t-j}^2 + {\sum_{l=1}^m}{b_3}{Y_l}$
%
%% *_Kurtosis Models_*
% $$GARCH-K(P,Q): k_t = d_0 + {\sum_{i=1}^p}d_{1,i}{\varepsilon}_{t-i}^4/{\sigma}_{t-i}^2 +{\sum_{j=1}^q}d_{2,j}k_{t-q}$
%
%% *_Distribution_*
%
% $$f(x) = \frac{{\Gamma}\left(\frac{{\nu_t}+1}{2} \right)}{\sqrt{{\nu_t}{\pi}}{\Gamma} \left( \frac{{\nu_t}}{2} \right)}\left(1+\frac{\epsilon_t^2}{\nu_t} \right)^{-\frac{\nu_t+1}{2}}$
%
% where the degrees of freedom can be expressed as a function of conditional kurtosis
%
% $$\nu_t = \frac{4k_t - 6}{k_t - 3}$
%
% <..\readme\readme.html Return to Main>


text:
function N = my_per_vertex_normals(V,F,varargin)
  % MY_PER_VERTEX_NORMALS  Compute per-vertex (area-weighted) normals over a mesh % (V,F)
  %
  % N = per_vertex_normals(V,F)
  %
  % Inputs:
  %   V  #V by 3 list of vertex positions
  %   F  #F by 3 list of triangle indices
  % Outputs:
  %   N  #V by 3 list of vertex normals, area-weighted
  %

  %Compute per-face normals.
  FN = normalizerow(normals(V,F));

  %Average to compute per-vertex normals.
  N = ...

end

author: odedstein

repo: sgi-introduction-course


text:
function view = computeSpatialGradient(view)
%
% view = computeSpatialGradient(view)
%
% Loads the meanMap, then uses estFilIntGrad to compute the spatial gradient.
%
% To use this spatial gradient for inhomogeneity correction,
% - set dataTYPES(*).blockedAnalysis(*).inhomoCorrect = 2
% - load spatial gradient (from the File menu)
% - compute corAnal (from the Analysis menu)
%
% djh, 7/11/01
% Ress, 04/05 Modified to allow NaNs in meanMap.

% Tried various alternatives for computing spatial gradient:
% - more blurring in estFilIntGrad
%   [not much difference]
% - aniso3
%   [looks essentially like no blurring]
% - do nothing
%   [looks a lot like estFilIntGrad, except for a scale factor]

% Load meanMap from spatialGradMap.mat, if it exists
view = loadMeanMap(view);

% Initialize to empty cell array
nScans = numScans(view);
map = cell(1,nScans);

%% Compute the robust estimate of spatial gradient from meanMap
% put up a mrvWaitbar if needed
verbose = prefsVerboseCheck;
if verbose
    waitHandle = mrvWaitbar(0,'Computing spatial gradient from mean images.  Please wait...');
end

% main loop: across scans
for iScan = 1:nScans
    map1 = view.map{iScan};
    map1(~isfinite(map1)) = min(map1(:));
    map{iScan} = estFilIntGrad(map1);
    %map{iScan} = estPolIntGrad(view.map{iScan},[3,3,3],1);
    %map{iScan} = aniso3(view.map{iScan},'tukeyPsi',100,100);

    % a 2nd check for non-finite values is necessary (e.g., for motioncomp
    % data)
    map{iScan}( isnan(map{iScan}) | isinf(map{iScan}) ) = min(map1(:));

    if verbose, mrvWaitbar(iScan/nScans); end
end

if verbose
    close(waitHandle);
end

% Initial save; so if setParameterMap fails we're not screwed
mapName = 'spatialGrad';
save( fullfile(dataDir(view), 'spatialGrad.mat'), 'map', 'mapName' );

% Set parameter map
view = setParameterMap(view,map,'spatialGrad');

% Save file (again)
saveParameterMap(view);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [x, c, funVal, ValueL]=LogisticR(A, y, z, opts)
%
%%
% Function LogisticR
%      Logistic Loss with the L1-norm Regularization
%
%% Problem
%
%  min  f(x,c) = - weight_i * log (p_i) + 1/2 rsL2 * ||x||_2^2
%                 + z * \|x\|_1
%
%  By default, rsL2=0.
%
%  a_i denotes a training sample,
%      and a_i' corresponds to the i-th row of the data matrix A
%
%  y_i (either 1 or -1) is the response
%
%  p_i= 1/ (1+ exp(-y_i (x' * a_i + c) ) ) denotes the probability
%
%  weight_i denotes the weight for the i-th sample
%
%% Input parameters:
%
%  A-         Matrix of size m x n
%                A can be a dense matrix
%                         a sparse matrix
%                         or a DCT matrix
%  y -        Response vector (of size mx1)
%  z -        Lq/L1 norm regularization parameter (z >=0)
%  opts-      Optional inputs (default value: opts=[])
%
%% Output parameters:
%  x-         The obtained weight of size n x 1
%  c-         The obtained intercept (scalar)
%  funVal-    Function value during iterations
%
%% Copyright (C) 2009-2010 Jun Liu, and Jieping Ye
%
% You are suggested to first read the Manual.
%
% For any problem, please contact with Jun Liu via j.liu@asu.edu
%
% Last modified on February 18, 2010.
%
%% Related papers
%
% [1]  Jun Liu, Jianhui Chen, and Jieping Ye, Large-Scale Sparse Logistic
%      Regression, KDD09
%
% [2]  Jun Liu and Jieping Ye, Efficient Euclidean Projections
%      in Linear Time, ICML 2009.
%
%% Related functions
%
%  sll_opts, initFactor, pathSolutionLeast
%  LogisticC, nnLogisticR, nnLogisticC
%
%%

%% Verify and initialize the parameters
%%
if (nargin <3)
    error('\n Inputs: A, y, and z should be specified!\n');
end

[m,n]=size(A);

if (length(y) ~=m)
    error('\n Check the length of y!\n');
end

if (z<=0)
    error('\n z should be positive!\n');
end

opts=sll_opts(opts); % run sll_opts to set default values (flags)

%% Detailed initialization
%% Normalization

% Please refer to sll_opts for the definitions of mu, nu and nFlag
%
% If .nFlag =1, the input matrix A is normalized to
%                     A= ( A- repmat(mu, m,1) ) * diag(nu)^{-1}
%
% If .nFlag =2, the input matrix A is normalized to
%                     A= diag(nu)^{-1} * ( A- repmat(mu, m,1) )
%
% Such normalization is done implicitly
%     This implicit normalization is suggested for the sparse matrix
%                                    but not for the dense matrix
%

if (opts.nFlag~=0)
    if (isfield(opts,'mu'))
        mu=opts.mu;
        if(size(mu,2)~=n)
            error('\n Check the input .mu');
        end
    else
        mu=mean(A,1);
    end

    if (opts.nFlag==1)
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=n)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,1)/m).^(0.5); nu=nu';
        end
    else % .nFlag=2
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=m)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,2)/n).^(0.5);
        end
    end

    ind_zero=find(abs(nu)<= 1e-10);    nu(ind_zero)=1;
    % If some values in nu is typically small, it might be that,
    % the entries in a given row or column in A are all close to zero.
    % For numerical stability, we set the corresponding value to 1.
end

if (~issparse(A)) && (opts.nFlag~=0)
    fprintf('\n -----------------------------------------------------');
    fprintf('\n The data is not sparse or not stored in sparse format');
    fprintf('\n The code still works.');
    fprintf('\n But we suggest you to normalize the data directly,');
    fprintf('\n for achieving better efficiency.');
    fprintf('\n -----------------------------------------------------');
end

%% Group & others

% The parameter 'weight' contains the weight for each training sample.
% See the definition of the problem above.
% The summation of the weights for all the samples equals to 1.
if (isfield(opts,'sWeight'))
    sWeight=opts.sWeight;

    if ( length(sWeight)~=2 || sWeight(1) <=0 || sWeight(2) <= 0)
        error('\n Check opts.sWeight, which contains two positive values');
    end

    % we process the weight, so that the summation of the weights for all
    % the samples is 1.

    p_flag=(y==1);                  % the indices of the postive samples
    m1=sum(p_flag) * sWeight(1);    % the total weight for the positive samples
    m2=sum(~p_flag) * sWeight(2);   % the total weight for the positive samples

    weight(p_flag,1)=sWeight(1)/(m1+m2);
    weight(~p_flag,1)=sWeight(2)/(m1+m2);
else
    weight=ones(m,1)/m;             % if not specified, we apply equal weight
end

%% Starting point initialization

% process the regularization parameter

% L2 norm regularization
if isfield(opts,'rsL2')
    rsL2=opts.rsL2;
    if (rsL2<0)
        error('\n opts.rsL2 should be nonnegative!');
    end
else
    rsL2=0;
end

p_flag=(y==1);                  % the indices of the postive samples
m1=sum(weight(p_flag));         % the total weight for the positive samples
m2=1-m1;                        % the total weight for the positive samples

% L1 norm regularization
if (opts.rFlag==0)
    lambda=z;
else % z here is the scaling factor lying in [0,1]
    if (z<0 || z>1)
        error('\n opts.rFlag=1, and z should be in [0,1]');
    end

    % we compute ATb for computing lambda_max, when the input z is a ratio

    b(p_flag,1)=m2;  b(~p_flag,1)=-m1;
    b=b.*weight;

    % compute AT b
    if (opts.nFlag==0)
        ATb =A'*b;
    elseif (opts.nFlag==1)
        ATb= A'*b - sum(b) * mu';  ATb=ATb./nu;
    else
        invNu=b./nu;               ATb=A'*invNu-sum(invNu)*mu';
    end

    lambda_max=max(abs(ATb));
    lambda=z*lambda_max;

    rsL2=rsL2*lambda_max; % the input rsL2 is a ratio of lambda_max
end

% initialize a starting point
if opts.init==2
    x=zeros(n,1); c=log(m1/m2);
else
    if isfield(opts,'x0')
        x=opts.x0;
        if (length(x)~=n)
            error('\n Check the input .x0');
        end
    else
        x=zeros(n,1);
    end

    if isfield(opts,'c0')
        c=opts.c0;
    else
        c=log(m1/m2);
    end
end

% compute A x
if (opts.nFlag==0)
    Ax=A* x;
elseif (opts.nFlag==1)
    invNu=x./nu; mu_invNu=mu * invNu;
    Ax=A*invNu -repmat(mu_invNu, m, 1);
else
    Ax=A*x-repmat(mu*x, m, 1);    Ax=Ax./nu;
end

%% The main program

%% The Armijo Goldstein line search scheme + accelearted gradient descent

if (opts.mFlag==0 && opts.lFlag==0)

    bFlag=0; % this flag tests whether the gradient step only changes a little

    L=1/m + rsL2; % the intial guess of the Lipschitz continuous gradient

    weighty=weight.*y;
    % the product between weight and y

    % assign xp with x, and Axp with Ax
    xp=x; Axp=Ax; xxp=zeros(n,1);
    cp=c;         ccp=0;

    alphap=0; alpha=1;

    for iterStep=1:opts.maxIter
        % --------------------------- step 1 ---------------------------
        % compute search point s based on xp and x (with beta)
        beta=(alphap-1)/alpha;    s=x + beta* xxp;  sc=c + beta* ccp;

        % --------------------------- step 2 ---------------------------
        % line search for L and compute the new approximate solution x

        % compute As=A*s
        As=Ax + beta* (Ax-Axp);

        % aa= - diag(y) * (A * s + sc)
        aa=- y.*(As+ sc);

        % fun_s is the logistic loss at the search point
        bb=max(aa,0);
        fun_s= weight' * ( log( exp(-bb) +  exp(aa-bb) ) + bb )+...
            rsL2/2 * s'*s;

        % compute prob=[p_1;p_2;...;p_m]
        prob=1./( 1+ exp(aa) );

        % b= - diag(y.* weight) * (1 - prob)
        b= -weighty.*(1-prob);

        gc=sum(b); % the gradient of c

        % compute g= AT b, the gradient of x
        if (opts.nFlag==0)
            g=A'*b;
        elseif (opts.nFlag==1)
            g=A'*b - sum(b) * mu';  g=g./nu;
        else
            invNu=b./nu;              g=A'*invNu-sum(invNu)*mu';
        end

        g=g+ rsL2 * s;  % add the squared L2 norm regularization term

        % copy x and Ax to xp and Axp
        xp=x;    Axp=Ax;
        cp=c;

        while (1)
            % let s walk in a step in the antigradient of s to get v
            % and then do the Lq/L1-norm regularized projection
            v=s-g/L;  c= sc- gc/L;

            % L1-norm regularized projection
            x=sign(v).*max(abs(v)-lambda / L,0);

            v=x-s;  % the difference between the new approximate solution x
            % and the search point s

            % compute A x
            if (opts.nFlag==0)
                Ax=A* x;
            elseif (opts.nFlag==1)
                invNu=x./nu; mu_invNu=mu * invNu;
                Ax=A*invNu -repmat(mu_invNu, m, 1);
            else
                Ax=A*x-repmat(mu*x, m, 1);     Ax=Ax./nu;
            end

            % aa= - diag(y) * (A * x + c)
            aa=- y.*(Ax+ c);

            % fun_x is the logistic loss at the new approximate solution
            bb=max(aa,0);
            fun_x= weight'* ( log( exp(-bb) +  exp(aa-bb) ) + bb ) +...
                rsL2/2 * x'*x;

            r_sum= (v'*v + (c-sc)^2) / 2;
            l_sum=fun_x - fun_s - v'* g - (c-sc)* gc;

            if (r_sum <=1e-20)
                bFlag=1; % this shows that, the gradient step makes little improvement
                break;
            end

            % the condition is fun_x <= fun_s + v'* g + c * gc
            %                           + L/2 * (v'*v + (c-sc)^2 )
            if(l_sum <= r_sum * L)
                break;
            else
                L=max(2*L, l_sum/r_sum);
                %fprintf('\n L=%e, r_sum=%e',L, r_sum);
            end
        end

        % --------------------------- step 3 ---------------------------
        % update alpha and alphap, and check whether converge
        alphap=alpha; alpha= (1+ sqrt(4*alpha*alpha +1))/2;

        ValueL(iterStep)=L;
        % store values for L

        xxp=x-xp;    ccp=c-cp;
        funVal(iterStep)=fun_x + lambda* sum(abs(x));

        if (bFlag)
            % fprintf('\n The program terminates as the gradient step changes the solution very small.');
            break;
        end

        switch(opts.tFlag)
            case 0
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                        break;
                    end
                end
            case 1
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                            opts.tol* funVal(iterStep-1))
                        break;
                    end
                end
            case 2
                if ( funVal(iterStep)<= opts.tol)
                    break;
                end
            case 3
                norm_xxp=sqrt(xxp'*xxp);
                if ( norm_xxp <=opts.tol)
                    break;
                end
            case 4
                norm_xp=sqrt(xp'*xp);    norm_xxp=sqrt(xxp'*xxp);
                if ( norm_xxp <=opts.tol * max(norm_xp,1))
                    break;
                end
            case 5
                if iterStep>=opts.maxIter
                    break;
                end
        end
    end
end

%% Reformulated problem + Nemirovski's line search scheme

% Original problem:           min  f(x,c) = - weight_i * log (p_i) + 1/2 rsL2 * ||x||_2^2
%                                                      + z * \|x\|_1
%
% Reformulated problem:   min  f(x,c) = - weight_i * log (p_i) + 1/2 rsL2 * ||x||_2^2
%                                                      + z * 1^T t
%                                   s.t., |x| <=t

if (opts.mFlag==1 && opts.lFlag==0)

    bFlag=0; % this flag tests whether the gradient step only changes a little

    L=1/m + rsL2; % the intial guess of the Lipschitz continuous gradient

    weighty=weight.*y;
    % the product between weight and y

    % assign xp with x, and Axp with Ax
    xp=x;         Axp=Ax;     xxp=zeros(n,1);
    cp=c;         ccp=0;

    t=abs(x); tp=t;
    % t is the upper bound of absolute value of x

    alphap=0; alpha=1;

    for iterStep=1:opts.maxIter
        % --------------------------- step 1 ---------------------------
        % compute search point s based on xp and x (with beta)
        beta=(alphap-1)/alpha;    s=x + beta* xxp;  sc=c + beta* ccp; s_t= t + beta * (t -tp);

        % --------------------------- step 2 ---------------------------
        % line search for L and compute the new approximate solution x

        % compute As=A*s
        As=Ax + beta* (Ax-Axp);

        % aa= - diag(y) * (A * s + sc)
        aa=- y.*(As+ sc);

        % fun_s is the logistic loss at the search point
        bb=max(aa,0);
        fun_s= weight' * ( log( exp(-bb) +  exp(aa-bb) ) + bb )+...
            rsL2/2 * s'*s;

        % compute prob=[p_1;p_2;...;p_m]
        prob=1./( 1+ exp(aa) );

        % b= - diag(y.* weight) * (1 - prob)
        b= -weighty.*(1-prob);

        gc=sum(b); % the gradient of c

        % compute g= AT b, the gradient of x
        if (opts.nFlag==0)
            g=A'*b;
        elseif (opts.nFlag==1)
            g=A'*b - sum(b) * mu';  g=g./nu;
        else
            invNu=b./nu;              g=A'*invNu-sum(invNu)*mu';
        end

        g=g+ rsL2 * s;  % add the squared L2 norm regularization term

        % copy x and Ax to xp and Axp
        xp=x;    Axp=Ax;
        cp=c;    tp=t;

        while (1)
            % let s walk in a step in the antigradient of s to get v
            % and then do the Lq/L1-norm regularized projection
            c= sc- gc/L;
            u=s-g/L;
            v= s_t - lambda / L;

            % projection
            [x, t]=ep1R(u,v,n);

            v=x-s;  % the difference between the new approximate solution x
            % and the search point s
            v_t=t-s_t;

            % compute A x
            if (opts.nFlag==0)
                Ax=A* x;
            elseif (opts.nFlag==1)
                invNu=x./nu; mu_invNu=mu * invNu;
                Ax=A*invNu -repmat(mu_invNu, m, 1);
            else
                Ax=A*x-repmat(mu*x, m, 1);     Ax=Ax./nu;
            end

            % aa= - diag(y) * (A * x + c)
            aa=- y.*(Ax+ c);

            % fun_x is the logistic loss at the new approximate solution
            bb=max(aa,0);
            fun_x= weight'* ( log( exp(-bb) +  exp(aa-bb) ) + bb ) +...
                rsL2/2 * x'*x;

            r_sum= (v'*v + (c-sc)^2 + v_t'*v_t) / 2;
            l_sum=fun_x - fun_s - v'* g - (c-sc)* gc;

            if (r_sum <=1e-20)
                bFlag=1; % this shows that, the gradient step makes little improvement
                break;
            end

            % the condition is fun_x <= fun_s + v'* g + c * gc
            %                           + L/2 * (v'*v + (c-sc)^2 )
            if(l_sum <= r_sum * L)
                break;
            else
                L=max(2*L, l_sum/r_sum);
                %fprintf('\n L=%e, r_sum=%e',L, r_sum);
            end
        end

        % --------------------------- step 3 ---------------------------
        % update alpha and alphap, and check whether converge
        alphap=alpha; alpha= (1+ sqrt(4*alpha*alpha +1))/2;

        ValueL(iterStep)=L;
        % store values for L

        xxp=x-xp;    ccp=c-cp;
        funVal(iterStep)=fun_x + lambda* sum(t);

        if (bFlag)
            % fprintf('\n The program terminates as the gradient step changes the solution very small.');
            break;
        end

        switch(opts.tFlag)
            case 0
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                        break;
                    end
                end
            case 1
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                            opts.tol* funVal(iterStep-1))
                        break;
                    end
                end
            case 2
                if ( funVal(iterStep)<= opts.tol)
                    break;
                end
            case 3
                norm_xxp=sqrt(xxp'*xxp + norm(t-tp)^2 + (c-cp)^2);
                if ( norm_xxp <=opts.tol)
                    break;
                end
            case 4
                norm_xp=sqrt(xp'*xp + cp^2 + tp'*tp);    norm_xxp=sqrt(xxp'*xxp + norm(t-tp)^2 + (c-cp)^2);
                if ( norm_xxp <=opts.tol * max(norm_xp,1))
                    break;
                end
            case 5
                if iterStep>=opts.maxIter
                    break;
                end
        end
    end
end

%% Reformulated problem + adaptive line search scheme

% Original problem:           min  f(x,c) = - weight_i * log (p_i) + 1/2 rsL2 * ||x||_2^2
%                                                      + z * \|x\|_1
%
% Reformulated problem:   min  f(x,c) = - weight_i * log (p_i) + 1/2 rsL2 * ||x||_2^2
%                                                      + z * 1^T t
%                                   s.t., |x| <=t

if (opts.mFlag==1 && opts.lFlag==1)

    bFlag=0; % this flag tests whether the gradient step only changes a little

    L=1/m + rsL2; % the intial guess of the Lipschitz continuous gradient

    weighty=weight.*y;
    % the product between weight and y

    gamma=1;
    % we shall set the value of gamma = L,
    % and L is appropriate for the starting point

    % assign xp with x, and Axp with Ax
    xp=x; Axp=Ax; xxp=zeros(n,1);
    cp=c;         ccp=0;
    t=abs(x); tp=t;
    % t is the upper bound of absolute value of x

    for iterStep=1:opts.maxIter
        while (1)
            if (iterStep~=1)
                alpha= (-gamma+ sqrt(gamma*gamma + 4* L * gamma)) / (2*L);
                beta= (gamma - gamma* alphap) / (alphap * gamma + alphap* L * alpha);
                % beta is the coefficient for generating search point s

                s=x + beta* xxp;   s_t= t + beta * (t -tp); sc=c + beta* ccp;
                As=Ax + beta* (Ax-Axp);
            else
                alpha= (-1+ sqrt(5)) / 2;  beta=0;
                s=x; s_t=t; sc= c;
                As=Ax;
            end

            % aa= - diag(y) * (A * s + sc)
            aa=- y.*(As+ sc);

            % fun_s is the logistic loss at the search point
            bb=max(aa,0);
            fun_s= weight' * ( log( exp(-bb) +  exp(aa-bb) ) + bb )+...
                rsL2/2 * s'*s;

            % compute prob=[p_1;p_2;...;p_m]
            prob=1./( 1+ exp(aa) );

            % b= - diag(y.* weight) * (1 - prob)
            b= -weighty.*(1-prob);

            gc=sum(b); % the gradient of c

            % compute g= AT b, the gradient of x
            if (opts.nFlag==0)
                g=A'*b;
            elseif (opts.nFlag==1)
                g=A'*b - sum(b) * mu';  g=g./nu;
            else
                invNu=b./nu;              g=A'*invNu-sum(invNu)*mu';
            end

            g=g+ rsL2 * s;  % add the squared L2 norm regularization term

            % let s walk in a step in the antigradient of s to get v
            % and then do the Lq/L1-norm regularized projection
            cnew= sc- gc/L;
            u=s-g/L;
            v= s_t - lambda / L;

            % projection
            [xnew, tnew]=ep1R(u,v,n);

            v=xnew-s;  % the difference between the new approximate solution x
            % and the search point s
            v_t=tnew-s_t;

            % compute A xnew
            if (opts.nFlag==0)
                Axnew=A* xnew;
            elseif (opts.nFlag==1)
                invNu=xnew./nu; mu_invNu=mu * invNu;
                Axnew=A*invNu -repmat(mu_invNu, m, 1);
            else
                Axnew=A*xnew-repmat(mu*xnew, m, 1);     Axnew=Axnew./nu;
            end

            % aa= - diag(y) * (A * x + c)
            aa=- y.*(Axnew+ cnew);

            % fun_x is the logistic loss at the new approximate solution
            bb=max(aa,0);
            fun_x= weight'* ( log( exp(-bb) +  exp(aa-bb) ) + bb ) +...
                rsL2/2 * xnew'*xnew;

            r_sum= (v'*v + (cnew-sc)^2 + v_t'*v_t) / 2;
            l_sum=fun_x - fun_s - v'* g - (cnew-sc)* gc;

            if (r_sum <=1e-20)
                bFlag=1; % this shows that, the gradient step makes little improvement
                break;
            end

            % the condition is fun_x <= fun_s + v'* g + c * gc
            %                           + L/2 * (v'*v + (c-sc)^2 )
            if(l_sum <= r_sum * L)
                break;
            else
                L=max(2*L, l_sum/r_sum);
                %fprintf('\n L=%e, r_sum=%e',L, r_sum);
            end
        end

        gamma=L* alpha* alpha;    alphap=alpha;
        % update gamma, and alphap

        ValueL(iterStep)=L;
        % store values for L

        tao=L * r_sum / l_sum;
        if (tao >=5)
            L=L*0.8;
        end
        % decrease the value of L

        xp=x;  x=xnew; xxp=x-xp;
        Axp=Ax; Ax=Axnew;
        tp=t; t=tnew;
        cp=c; c=cnew; ccp=c-cp;
        % update x, t, c, and Ax

        funVal(iterStep)=fun_x + lambda* sum(t);

        if (bFlag)
            % fprintf('\n The program terminates as the gradient step changes the solution very small.');
            break;
        end

        switch(opts.tFlag)
            case 0
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                        break;
                    end
                end
            case 1
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                            opts.tol* funVal(iterStep-1))
                        break;
                    end
                end
            case 2
                if ( funVal(iterStep)<= opts.tol)
                    break;
                end
            case 3
                norm_xxp=sqrt(xxp'*xxp + norm(t-tp)^2 + (c-cp)^2);
                if ( norm_xxp <=opts.tol)
                    break;
                end
            case 4
                norm_xp=sqrt(xp'*xp + cp^2 + tp'*tp);    norm_xxp=sqrt(xxp'*xxp + norm(t-tp)^2 + (c-cp)^2);
                if ( norm_xxp <=opts.tol * max(norm_xp,1))
                    break;
                end
            case 5
                if iterStep>=opts.maxIter
                    break;
                end
        end
    end
end

%%
if(opts.mFlag==0 && opts.lFlag==1)
    error('\n The function does not support opts.mFlag=0 & opts.lFlag=1!');
end

author: thomaskuestner


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

nsteps = 10;
Vfrac_final = 0.15;
Perimeter_target=3.5;
optimality_final =1e-4;
constr_final =1e-4;

BCscale_factor = 0.3;
HJiter0 = 1;
e2 = 1;
N_holes = [17 9 9];
R_holes = 0.2;
phase_holes = [0 0 0];

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

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

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

plotting = 1;
printing = 0;
monitoring = 1;
monitoring_interval = 1;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function[c]=materncfun(al,ga)
%MATERNCFUN Returns the normalization or C-function for a Matern process.
%
%   MATERNCFUN is a low-level function used in the JMATERN package.
%
%   MATERNCFUN(ALPHA) returns the Matern C-function defined in Eqn. (44) of
%   Lilly et al. (2017).
%
%   MATERNCFUN(ALPHA,GAMMA) returns C-function for the two-parameter
%   generalized Matern process.  See MATERNSPEC for details.
%
%   Usage: c=materncfun(alpha);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2015--2020 J.M. Lilly --- type 'help jlab_license' for details

%c=frac(1,2*pi)*frac(gamma(1/2)*gamma(alpha-1/2),gamma(alpha));
if nargin==1
    if al>1/2
        c=frac(1,2*sqrt(pi))*frac(gamma(al-1/2),gamma(al));
    else
        c=frac(2.^(abs(al-1/2)-abs(al)),sqrt(2*pi)).*...
            frac(gamma(abs(al-1/2)),gamma(abs(al)));
    end
else
    c=frac(1,2*pi*ga).*beta(frac(1,2*ga),frac(1,ga).*(al-1/2));
end
%Verified that these two forms are the same for al>1/2
%al1=[1/2:.01:10]';
%gamma1=[0:.01:10]';
%[al,gamma]=meshgrid(al1,gamma1);
%c=materncfun(al,gamma);
%jpcolor(al1,gamma1,log10(c)),caxis([-1 0])

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function [ Ycap ] = lraNTF( Y,opts )
%% Low-rank approximation based Non-Negative Tensor(CP) Factorization
%   Usage: Ycap=lraNCP(Y,opts);
%   Input: Y can be ktensor/ttensor/tensor/double
%  Output: Ycap is a nonnegative ktensor
%
%  Parameters:
%   opts.NumOfComp:  rank of Ycap
%   opts.lra_rank:   rank of low-rank approximation of Y
%   opts.CPalg:      a -mat file including the configuration of
%                           unconstrained CPD
%   opts.maxiter:    max. of global iterations    [default:50]
%   opts.maxiniter: max. of inner iterations. [default:20]
%   opts.nlssolver: [HALS|APG|MU] the algorithm used for solving NLS
%                      sub-problem
%   opts.initmode:  [cp|random]
%   opts.tol:        stop threshold. [default:1e-6]
%   opts.sparsity  an N-by-1 vector in [0,1) to specify the sparsity of
%                   results (l1-norm based sparsity).
%
%
%
% Cite this work currently:
%  Guoxu Zhou; Cichocki, A.; Shengli Xie; , "Fast Nonnegative Matrix/Tensor Factorization Based on Low-Rank Approximation,"
%  IEEE Transactions on Signal Processing, vol.60, no.6, pp.2928-2940, June 2012
%  doi: 10.1109/TSP.2012.2190410
%  URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=6166354&isnumber=6198804
%
% By Guoxu Zhou
% E-Mail: zhouguoxu@gmail.com
% Last updated: 2013-09-09

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function [ value_min, value_mean, value_max, value_var ] = ...
  tet_mesh_quality2 ( node_num, node_xyz, tetra_order, tetra_num, ...
  tetra_node )

%*****************************************************************************80
%
%% TET_MESH_QUALITY2 returns a tet mesh quality factor.
%
%  Discussion:
%
%    The tet mesh quality measure is the minimum of the
%    corresponding tetrahedron quality measure, over all tetrahedrons in the
%    tet mesh.
%
%    This routine is designed for an order-4 tet mesh.  Order 10 tet meshes
%    may be input, but the extra nodes are ignored.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_XYZ(3,NODE_NUM), the nodes.
%
%    Input, integer TETRA_ORDER, the order of the tetrahedrons.
%
%    Input, integer TETRA_NUM, the number of tetrahedrons.
%
%    Input, integer TETRA_NODE(TETRA_ORDER,TETRA_NUM), the indices of the nodes
%    that make up the tetrahedrons.
%
%    Output, real VALUE_MIN, VALUE_MEAN, VALUE_MAX, VALUE_VAR,
%    the minimum, mean, maximum and variance of the quality measure.
%
  dim_num = 3;

  for tetra = 1 : tetra_num

    tetrahedron(1:dim_num,1:4) = node_xyz(1:dim_num,tetra_node(1:4,tetra));

    tetrahedron_quality(tetra) = tetrahedron_quality2_3d ( tetrahedron );

  end

  value_max  = max  ( tetrahedron_quality(1:tetra_num) );
  value_min  = min  ( tetrahedron_quality(1:tetra_num) );
  value_mean = mean ( tetrahedron_quality(1:tetra_num) );
  value_var  = var  ( tetrahedron_quality(1:tetra_num) );

  return
end


text:
function [bestEpoch, bestNet]= pickBestNet(sessionID, N, verbose)
    if nargin<2, N= 5; end
    if nargin<3, verbose= 2; end

    paths= localPaths();
    outFnLatest= sprintf('%s%s_latest.mat', paths.outPrefix, sessionID);
    res= load( outFnLatest, 'obj', 'opts', 'auxData');

    if res.opts.epochTestFrequency~=1
        error('This code assumes epochTestFrequency==1 (it is %d)', res.opts.epochTestFrequency);
    end
    bestEpoch= getBestEpoch(res.obj.val.recall, res.opts.recallNs, N);
    assert(~isempty(bestEpoch));

    if verbose>0
        whichRecallInds= 1:6;
        hline= repmat('=', 1, length(whichRecallInds)*5-1 + 14);

        relja_display('%s Best epoch: %d (out of %d)', sessionID, bestEpoch, size(res.obj.val.recall,2) );

        if isfield(res.obj, 'pretrain')
            offtheshelfValRecs= res.obj.pretrain.val.recall( whichRecallInds );
        end
        bestValRecs= res.obj.val.recall( whichRecallInds, bestEpoch );

        relja_display('%s', hline);
        recallStr= sprintf('%04d ', res.opts.recallNs(whichRecallInds) );
        relja_display('Recall@N      %s', recallStr);
        relja_display('%s', hline);

        if exist('offtheshelfValRecs', 'var')
            offtheshelfStr= sprintf('%.2f ', offtheshelfValRecs);
        end
        trainedStr= sprintf('%.2f ', bestValRecs);

        relja_display('off-the-shelf %s', offtheshelfStr);
        relja_display('our trained   %s', trainedStr);

        if exist('offtheshelfValRecs', 'var')
            relImpStr= sprintf('%.2f ', bestValRecs./offtheshelfValRecs);
            relja_display('trained/shelf %s', relImpStr);
        end

        if verbose>1

            figure; plotResults(res.obj, res.opts, res.auxData);

            figure;
            plot( res.opts.recallNs, res.obj.pretrain.val.recall, 'bx-' );
            hold on;
            plot( res.opts.recallNs, res.obj.val.recall(:, bestEpoch), 'ro-' );
            xlabel('N');
            ylabel('Recall@N');
            xlim([0,50]);
            grid on;
            title( sprintf('%s %s %s %s %s', sessionID, res.opts.netID, res.opts.layerName, res.opts.dbValName, res.opts.method), 'Interpreter', 'none' );
            legend( 'off-the-shelf', 'our trained', 'Location', 'SouthEast');
        end
    end

    if nargout>1
        outFnBest= sprintf('%s%s_ep%06d_latest.mat', paths.outPrefix, sessionID, bestEpoch);
        load(outFnBest, 'net');
        bestNet= net; clear net;
        bestNet.meta.sessionID= sessionID;
        bestNet.meta.epoch= bestEpoch;

        % remove unneeded momentum
        for iLayer= 1:length(bestNet.layers)
            if isfield(bestNet.layers{iLayer}, 'momentum')
                bestNet.layers{iLayer}= rmfield(bestNet.layers{iLayer}, 'momentum');
            elseif isprop(bestNet.layers{iLayer}, 'momentum')
                bestNet.layers{iLayer}.momentum= [];
            end
        end
    end
end

author: Relja

repo: netvlad

save_path: github-repos/MATLAB/Relja-netvlad


text:
%% fmin_symbolic_form
% Below is a demonstration of the features of the |fmin_symbolic_form| function

%%
clear; close all; clc;

%% Syntax
% |Fmin=fmin_symbolic_form(Pn,OPT_struct,F);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function demo4surf
% Demonstration of a bilinear surface.
%

% D.M. Spink
% Copyright (c) 2000

srf = nrb4surf([0.0 0.0 0.5],[1.0 0.0 -0.5],[0.0 1.0 -0.5],[1.0 1.0 0.5]);
nrbplot(srf,[10,10]);
title('Construction of a bilinear surface.');