text:
function [ fea, out ] = ex_convdiff4( varargin )
%EX_CONVDIFF4 1D Burgers equation (convection and diffusion) example.
%
%   [ FEA, OUT ] = EX_CONVDIFF4( VARARGIN ) 1D Burgers equation with steady solution,
%   u_t + (b*u-c)*u_x - nu*u_xx = 0 with exact solution c/b*(1-tanh(c/(2*nu)*(x-x0))).
%   Tests both time dependent and steady solvers. Accepts the following property/value pairs.
%
%       Input       Value/{Default}        Description
%       -----------------------------------------------------------------------------------
%       b           scalar {1.0}           Strength of nonlinearity
%       c           scalar {0.13}          Convection velocity
%       nu          scalar {0.01}          Diffusion coefficient
%       x0          scalar {0.5}           Posision of smooth shock
%       hmax        scalar {1/25}          Max grid cell size
%       ischeme     scalar {-1}            Solver scheme (<0 = stationary)
%       nsolve      scalar {2}             Nonlinear solver (when ischeme<0)
%       dt          scalar {0.1}           Time step size
%       sfun        string {sflag1}        Shape function
%       iplot       scalar 0/{1}           Plot solution (=1)
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       fea         struct                 Problem definition struct
%       out         struct                 Output struct

% Copyright 2013-2022 Precise Simulation, Ltd.

cOptDef = { ...
            'b',        1.0; ...
            'c',        0.13; ...
            'nu',       0.01; ...
            'x0',       0.5; ...
            'hmax',     1/20; ...
            'ischeme'   -1; ...
            'nsolve'    2; ...
            'dt'        0.1; ...
            'sfun',     'sflag1'; ...
            'iplot',    1; ...
            'tol',      1e-2; ...
            'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;

b  = opt.b;
c  = opt.c;
nu = opt.nu;
x0 = opt.x0;
refsol = [num2str(c/b),'*(1-tanh(',num2str(c/(2*nu)),'*(x-',num2str(x0),')))'];
bcd1   = c/b*(1-tanh((c/(2*nu)*(0-x0))));
bcd2   = c/b*(1-tanh((c/(2*nu)*(1-x0))));

% Grid generation.
fea.grid = linegrid( 1/opt.hmax, 0, 1 );

% Problem definition.
fea.sdim  = { 'x' };
fea = addphys( fea, @convectiondiffusion );
fea.phys.cd.sfun = { opt.sfun };
fea.phys.cd.eqn.coef{2,4} = { opt.nu };
fea.phys.cd.eqn.coef{3,4} = { [num2str(b),'*c-',num2str(c)] };
fea = parsephys(fea);

% Parse and solve problem.
fea = parseprob( fea );
fea.bdr.d{1} = bcd1;
fea.bdr.d{2} = bcd2;
fea.bdr.n    = cell(1,2);
x = fea.grid.p';
if( strcmp( opt.sfun,'sflag2' ) )
  x = [ x; (x(2:end)+x(1:end-1))/2 ];
end
if( opt.ischeme<0 )
  init = 0;
  jac.form  = {[1;1]};
  jac.coef  = {[num2str(b),'*cx']};
  fea.sol.u = solvestat( fea, 'fid', fid, 'init', init, 'maxnit', 1000, 'nsolve', 2, 'jac', jac, 'nsolve', opt.nsolve );
else
  init = [num2str(bcd1),'+x*',num2str(bcd2-bcd1)];
  [fea.sol.u,tlist] = solvetime( fea, 'fid', fid, 'init', init, 'ischeme', opt.ischeme, 'tstep', opt.dt, 'tmax', 10 );
end

% Postprocessing.
if( opt.iplot>0 )
  figure
  i_sol = size(fea.sol.u,2);
  [~,ix] = sort( x );
  postplot( fea, 'surfexpr', 'c', 'solnum', i_sol );
  hold on
  u_r = real( eval( refsol ) );
  plot( sort(x), u_r(ix), 'r--' );
  title( 'Steady solution' )
  xlabel( 'x' )
  drawnow
end

% Error checking.
i_sol = size(fea.sol.u,2);
u_i   = fea.sol.u(:,i_sol);
u_r   = real( eval( refsol ) );
errnm = norm( u_i - u_r )/norm( u_r );
out.err  = errnm;
out.pass = all( errnm<opt.tol );

if ( nargout==0 )
  clear fea out
end

author: precise-simulation

repo: featool-multiphysics


text:
% Copyright 2011 Zdenek Kalal
%
% This file is part of TLD.
%
% TLD 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.
%
% TLD 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 TLD.  If not, see <http://www.gnu.org/licenses/>.

function flag = bb_isin(bb,imsize)

flag = bb(1,:) >= 1 & ...
bb(2,:) >= 1 & ...
bb(3,:) <= imsize(2) & ...
bb(4,:) <= imsize(1);

repo: OpenTLD


text:
function surf = imJointSurface(img, L1, L2, varargin)
%IMJOINTSURFACE Surface area of the interface between two labels
%
%   Deprecated, use 'imJointSurfaceArea' instead.
%
%   S = imJointSurface(LBL, L1, L2)
%   Estimates the joint surface area between the two labels L1 and L2 in
%   the label image LBL.
%
%   S = imJointSurface(LBL, L1, L2, NDIRS)
%   Specifies the number of directions used for estimating surface area.
%   NDIRS can be either 3 or 13 (the default).
%
%   S = imJointSurface(..., RESOL)
%   Specifies image resolution. RESOL is a 1-by-3 row vector containing
%   resolution in the X, Y and Z direction (in that order).
%
%
%   Example
%     % generate a demo image
%     img = discreteBall(1:10, 1:100, 1:100, [50.12 50.23 50.34 40]);
%     % convert to image with two different labels
%     img2 = uint8(img + 1);
%     % compute joint surface area
%     imJointSurface(img2, 1, 2)
%     ans =
%         2.0102e+004
%
%   See also
%     imJointSurfaceArea, imSurface

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2010-07-26,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

warning('MatImage:deprecated', ...
    'function imJointSurface is obsolete, use imJointSurfaceArea instead');

% check image dimension and type
if ndims(img) ~= 3 || islogical(img)
    error('first argument should be a 3D image');
end

%% Process input arguments

% default number of directions
ndir = 13;

% default image resolution
delta = [1 1 1];

% Process user input arguments
while ~isempty(varargin)
    var = varargin{1};
    if ~isnumeric(var)
        error('option should be numeric');
    end

    % option is either connectivity or resolution
    if isscalar(var)
        ndir = var;
    else
        delta = var;
    end
    varargin(1) = [];
end

%% Initialisations

% distances between a pixel and its neighbours.
d1  = delta(1);
d2  = delta(2);
d3  = delta(3);

% volume of a voxel (used for computing line densities)
vol = d1 * d2 * d3;

%% Main processing for 3 directions

% number of transitions along the 3 main directions
n1a = sum(sum(sum(img(:,1:end-1,:)==L1 & img(:,2:end,:)==L2)));
n1b = sum(sum(sum(img(:,1:end-1,:)==L2 & img(:,2:end,:)==L1)));
n2a = sum(sum(sum(img(1:end-1,:,:)==L1 & img(2:end,:,:)==L2)));
n2b = sum(sum(sum(img(1:end-1,:,:)==L2 & img(2:end,:,:)==L1)));
n3a = sum(sum(sum(img(:,:,1:end-1)==L1 & img(:,:,2:end)==L2)));
n3b = sum(sum(sum(img(:,:,1:end-1)==L2 & img(:,:,2:end)==L1)));

if ndir == 3
    % compute surface area by averaging over the 3 main directions
    surf = 4/3 * ((n1a+n1b)/d1 + (n2a+n2b)/d2 + (n3a+n3b)/d3) / 2 * vol;
    return;
end

%% Additional processing for 13 directions

% Number of connected components along diagonals contained in the three
% main planes
n4a = sum(sum(sum(img(2:end,1:end-1,:)==L1   & img(1:end-1,2:end,:)==L2)));
n4b = sum(sum(sum(img(2:end,1:end-1,:)==L2   & img(1:end-1,2:end,:)==L1)));
n5a = sum(sum(sum(img(1:end-1,1:end-1,:)==L1 & img(2:end,2:end,:)==L2)));
n5b = sum(sum(sum(img(1:end-1,1:end-1,:)==L2 & img(2:end,2:end,:)==L1)));
n6a = sum(sum(sum(img(:,2:end,1:end-1)==L1   & img(:,1:end-1,2:end)==L2)));
n6b = sum(sum(sum(img(:,2:end,1:end-1)==L2   & img(:,1:end-1,2:end)==L1)));
n7a = sum(sum(sum(img(:,1:end-1,1:end-1)==L1 & img(:,2:end,2:end)==L2)));
n7b = sum(sum(sum(img(:,1:end-1,1:end-1)==L2 & img(:,2:end,2:end)==L1)));
n8a = sum(sum(sum(img(2:end,:,1:end-1)==L1   & img(1:end-1,:,2:end)==L2)));
n8b = sum(sum(sum(img(2:end,:,1:end-1)==L2   & img(1:end-1,:,2:end)==L1)));
n9a = sum(sum(sum(img(1:end-1,:,1:end-1)==L1 & img(2:end,:,2:end)==L2)));
n9b = sum(sum(sum(img(1:end-1,:,1:end-1)==L2 & img(2:end,:,2:end)==L1)));

%TODO: add the case of 9 directions ?

% Number of connected components along lines corresponding to diagonals of
% the unit cube
n10a = sum(sum(sum(img(1:end-1,1:end-1,1:end-1)==L1 & img(2:end,2:end,2:end)==L2)));
n10b = sum(sum(sum(img(1:end-1,1:end-1,1:end-1)==L2 & img(2:end,2:end,2:end)==L1)));
n11a = sum(sum(sum(img(2:end,1:end-1,1:end-1)==L1 & img(1:end-1,2:end,2:end)==L2)));
n11b = sum(sum(sum(img(2:end,1:end-1,1:end-1)==L2 & img(1:end-1,2:end,2:end)==L1)));
n12a = sum(sum(sum(img(1:end-1,2:end,1:end-1)==L1 & img(2:end,1:end-1,2:end)==L2)));
n12b = sum(sum(sum(img(1:end-1,2:end,1:end-1)==L2 & img(2:end,1:end-1,2:end)==L1)));
n13a = sum(sum(sum(img(2:end,2:end,1:end-1)==L1 & img(1:end-1,1:end-1,2:end)==L2)));
n13b = sum(sum(sum(img(2:end,2:end,1:end-1)==L2 & img(1:end-1,1:end-1,2:end)==L1)));

% space between 2 voxels in each direction
d12  = hypot(d1, d2);
d13  = hypot(d1, d3);
d23  = hypot(d2, d3);
d123 = sqrt(d1^2 + d2^2 + d3^2);

% Compute weights corresponding to surface fraction of spherical caps
% For isotropic case, weights correspond to:
% c1 = 0.04577789120476 * 2;  % Ox
% c2 = 0.04577789120476 * 2;  % Oy
% c3 = 0.04577789120476 * 2;  % Oz
% c4 = 0.03698062787608 * 2;  % Oxy
% c6 = 0.03698062787608 * 2;  % Oxz
% c8 = 0.03698062787608 * 2;  % Oyz
% c10 = 0.03519563978232 * 2;  % Oxyz
c = computeDirectionWeights3d13(delta);

% compute the weighted sum of each direction
% intersection count * direction weight / line density
surf = 4 * vol * (...
    (n1a+n1b)*c(1)/d1 + (n2a+n2b)*c(2)/d2 + (n3a+n3b)*c(3)/d3 + ...
    (n4a+n4b+n5a+n5b)*c(4)/d12 + (n6a+n6b+n7a+n7b)*c(6)/d13 + ...
    (n8a+n8b+n9a+n9b)*c(8)/d23 + ...
    (n10a + n10b + n11a + n11b + n12a + n12b + n13a + n13b) * c(10) / d123) / 2;

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [features, sp2reg] = computeFeatures(imName, paths, typ, param, data)

	clusters = data.clusters;
	superpixels = data.superpixels;
	nSP = max(superpixels(:));

	%For each region set compute the regions features
	for i = 1:size(clusters, 2),
		%calclate the sp2reg for
		cluster = clusters(:,i);
		nR = max(cluster);
		sp2reg{i} = false(nSP,nR);
		sp2reg{i}(sub2ind(size(sp2reg{i}),[1:nSP]',cluster)) = true;
		if(i > 1)
			[TF LOC] = ismember(sp2reg{i}', sp2reg{i-1}', 'rows');
			features{i}(:,find(TF)) = features{i-1}(:,LOC(find(TF)));
			if(~isempty(find(~TF)))
				features{i}(:,find(~TF)) = param.fName(superpixels, sp2reg{i}(:,find(~TF)), param, data);
			end
		else
			%Get features.
			features{i} = param.fName(superpixels, sp2reg{i}, param, data);
		end
		fprintf('.');
	end

	features = features;
	sp2reg = sp2reg;

	%Save the things
	fileName = fullfile(paths.featuresDir, typ, strcat(imName, '.mat'));
	save(fileName, 'clusters', 'superpixels', 'sp2reg', 'features');
end

author: s-gupta

repo: rgbd

save_path: github-repos/MATLAB/s-gupta-rgbd


text:
function cluster_orthviews_overlap_3colors(mask1, mask2, mask3, varargin)
% Plot three clusters on the orthviews, and their intersections for
% activations and deactivations in intermediate colors
%
% :Usage:
% ::
%
%    cluster_orthviews_overlap_3colors(mask1, mask2, mask3, ['colors', colors cell], ['surface'] )
%
% Positive effects only!
%
% :Inputs:
%
%   mask1 mask2 mask3 are either image names (preferred!) or clusters
%   structures.  for clusters structures, need to add .dim field
%   and you need the 2010 object-oriented code in the canlab repository.
%
% ..
%    tor wager, june 2010
%    A good function, but not polished yet.  needs debugging for various
%    minor issues.
% ..

colors = {[0 0 1] [1 0 0] [1 1 0]};
dosurface = 0;

for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}
            % reserved keywords
            case 'betas'
            case 'design'

                % functional commands
            case 'colors', colors = varargin{i+1};
            case 'surface', dosurface = 1;

            otherwise, warning(['Unknown input string option:' varargin{i}]);
        end
    end
end

if isstruct(mask1) && isstruct(mask2) && isstruct(mask3)
    loadfromfile = 0;
elseif ischar(mask1) && ischar(mask2) && ischar(mask3)
    loadfromfile = 1;
else
    error('Inputs must be all cl structs or mask filenames');
end

if loadfromfile

    [maskInfo1, dat1] = iimg_read_img(mask1, 2);
    [maskInfo2, dat2] = iimg_read_img(mask2, 2);
    [maskInfo3, dat3] = iimg_read_img(mask3, 2);

    % positive only
    dat1 =  dat1 > 0;
    dat2 =  dat2 > 0;
    dat3 =  dat3 > 0;

    int_dat = iimg_intersection(mask1, mask2, 'posneg'); % intersections
    pos12 = int_dat(:, 1);

    int_dat = iimg_intersection(mask1, mask3, 'posneg'); % intersections
    pos13 = int_dat(:, 1);

    int_dat = iimg_intersection(mask2, mask3, 'posneg'); % intersections
    pos23 = int_dat(:, 1);

    pos123 = pos12 & pos13;

    intersect_all = any([pos12 pos13 pos23 pos123], 2);

    %remove higher-order intersections from  mask dat
    dat1(intersect_all) = 0;
    dat2(intersect_all) = 0;
    dat3(intersect_all) = 0;

    pos12(pos123) = 0;
    pos23(pos123) = 0;
    pos13(pos123) = 0;

    % make clusters
    clear cl

    cl{1} =  iimg_indx2clusters(dat1, maskInfo1);
    cl{2} =  iimg_indx2clusters(dat2, maskInfo2);
    cl{3} =  iimg_indx2clusters(dat3, maskInfo3);

    cl{4} =  iimg_indx2clusters(pos12, maskInfo1);
    cl{5} =  iimg_indx2clusters(pos13, maskInfo2);
    cl{6} =  iimg_indx2clusters(pos23, maskInfo3);

    cl{7} =  iimg_indx2clusters(pos123, maskInfo1);

else
    % Work from clusters directly

    % convert to masks

    if ~isfield(mask1, 'dim') || isempty(mask1(1).dim) || ...
            ~isfield(mask2, 'dim') || isempty(mask2(1).dim) || ...
            ~isfield(mask3, 'dim') || isempty(mask3(1).dim)

        error('Enter .dim (image dimension in voxels) field in all cl structs. e.g., [91 109 91]');
    end

    r = cluster2region(mask1);
    m1 = region2imagevec(r);

    r = cluster2region(mask2);
    m2 = region2imagevec(r);

    r = cluster2region(mask3);
    m3 = region2imagevec(r);

    m1 = replace_empty(m1);
    m2 = replace_empty(m2);
    m3 = replace_empty(m3);

    dat1 =  m1.dat > 0;
    dat2 =  m2.dat > 0;
    dat3 =  m3.dat > 0;

    pos123 = all([dat1 dat2 dat3], 2);
    pos12 = all([dat1 dat2], 2);
    pos23 = all([dat2 dat3], 2);
    pos13 = all([dat1 dat3], 2);

    intersect_all = any([pos12 pos13 pos23 pos123], 2);

    %remove higher-order intersections from  mask dat
    dat1(intersect_all) = 0;
    dat2(intersect_all) = 0;
    dat3(intersect_all) = 0;

    pos12(pos123) = 0;
    pos23(pos123) = 0;
    pos13(pos123) = 0;

    % make clusters
    clear cl

    cl{1} =  iimg_indx2clusters(dat1, m1.volInfo);
    cl{2} =  iimg_indx2clusters(dat2, m2.volInfo);
    cl{3} =  iimg_indx2clusters(dat3, m3.volInfo);

    cl{4} =  iimg_indx2clusters(pos12, m1.volInfo);
    cl{5} =  iimg_indx2clusters(pos13, m1.volInfo);
    cl{6} =  iimg_indx2clusters(pos23, m1.volInfo);

    cl{7} =  iimg_indx2clusters(pos123, m1.volInfo);

end

% set colors for overlap
% colors{end+1} = {mean([colors{1}; colors{2}])};
% colors{end+1} = {mean([colors{1}; colors{3}])};
% colors{end+1} = {mean([colors{2}; colors{3}])};
%
% colors{end+1} = {mean([colors{1}; colors{2}; colors{3}])};

colors(end+1) = {2*mean([colors{1}; colors{2}])}; colors{end}(colors{end} > 1) = 1;
colors(end+1) = {2*mean([colors{1}; colors{3}])}; colors{end}(colors{end} > 1) = 1;
colors(end+1) = {2*mean([colors{2}; colors{3}])}; colors{end}(colors{end} > 1) = 1;

colors(end+1) = {3*mean([colors{1}; colors{2}; colors{3}])};

cluster_orthviews();
for i = 1:7
    if ~isempty(cl{i}), cluster_orthviews(cl{i}, colors(i), 'add', 'solid'); end
end

cluster_orthviews_montage(10, 'axial', [], 'onerow', 'range', [-35 55]);
h = findobj(gcf,'Type', 'Text'); set(h, 'FontSize', 18) %delete(h);

cluster_orthviews_montage(10, 'sagittal', [], 'onerow', 'range', [-40 40]);
h = findobj(gcf,'Type', 'Text'); set(h, 'FontSize', 18) %delete(h);

if dosurface

    create_figure;
    s = addbrain('hires left'); set(s, 'FaceColor', [.5 .5 .5], 'FaceAlpha', 1);
    for i = 1:7
        cluster_surf(cl{i}, colors(i), 2, s);
    end
    for i = 4:7
        cluster_surf(cl{i}, colors{i}(1), 2, s);
    end

    scn_export_papersetup; saveas(gcf, 'cluster_overlap_3colors_Lmedial.png');
    view(270, 0); lightFollowView
    scn_export_papersetup; saveas(gcf, 'cluster_overlap_3colors_Lateral.png');

    create_figure('s2');
    s = addbrain('hires right'); set(s, 'FaceColor', [.5 .5 .5], 'FaceAlpha', 1);
    for i = 1:3
        cluster_surf(cl{i}, colors(i), 2, s);
    end
    for i = 4:7
        cluster_surf(cl{i}, colors{i}(1), 2, s);
    end
    lightFollowView
    scn_export_papersetup; saveas(gcf, 'cluster_overlap_3colors_Rmedial.png');
    view(90, 0); lightFollowView
    scn_export_papersetup; saveas(gcf, 'cluster_overlap_3colors_Rateral.png');

    create_figure('s3');
    s = addbrain('limbic');
    s = [s addbrain('brainstem')];
    axis auto
    set(s, 'FaceColor', [.5 .5 .5]);
    set(s(end-1), 'FaceAlpha', .15);
    view(135, 10)
    for i = 1:3
        cluster_surf(cl{i}, colors(i), 2, s);
    end
    for i = 4:7
        cluster_surf(cl{i}, colors{i}(1), 2, s);
    end
    scn_export_papersetup; saveas(gcf, 'cluster_overlap_3colors_Lmedial.png');

end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [D, Y, optinf] = cbpdndl(D0, S, lambda, opt)

% cbpdndl -- Convolutional BPDN Dictionary Learning
%
%         argmin_{x_m,d_m} (1/2) \sum_k ||\sum_m d_m * x_k,m - s_k||_2^2 +
%                           lambda \sum_k \sum_m ||x_k,m||_1
%
%         The solution is computed using Augmented Lagrangian methods
%         (see boyd-2010-distributed) with efficient solution of the
%         main linear systems (see wohlberg-2016-efficient).
%
% Usage:
%       [D, Y, optinf] = cbpdndl(D0, S, lambda, opt)
%
% Input:
%       D0          Initial dictionary
%       S           Input images
%       lambda      Regularization parameter
%       opt         Options/algorithm parameters structure (see below)
%
% Output:
%       D           Dictionary filter set (3D array)
%       X           Coefficient maps (4D array)
%       optinf      Details of optimisation
%
%
% Options structure fields:
%   Verbose          Flag determining whether iteration status is displayed.
%                    Fields are iteration number, functional value,
%                    data fidelity term, l1 regularisation term, and
%                    primal and dual residuals (see Sec. 3.3 of
%                    boyd-2010-distributed). The values of rho and sigma
%                    are also displayed if options request that they are
%                    automatically adjusted.
%   MaxMainIter      Maximum main iterations
%   AbsStopTol       Absolute convergence tolerance (see Sec. 3.3.1 of
%                    boyd-2010-distributed)
%   RelStopTol       Relative convergence tolerance (see Sec. 3.3.1 of
%                    boyd-2010-distributed)
%   L1Weight         Weight array for L1 norm
%   Y0               Initial value for Y
%   U0               Initial value for U
%   G0               Initial value for G (overrides D0 if specified)
%   H0               Initial value for H
%   rho              Augmented Lagrangian penalty parameter
%   AutoRho          Flag determining whether rho is automatically updated
%                    (see Sec. 3.4.1 of boyd-2010-distributed)
%   AutoRhoPeriod    Iteration period on which rho is updated
%   RhoRsdlRatio     Primal/dual residual ratio in rho update test
%   RhoScaling       Multiplier applied to rho when updated
%   AutoRhoScaling   Flag determining whether RhoScaling value is
%                    adaptively determined (see wohlberg-2015-adaptive). If
%                    enabled, RhoScaling specifies a maximum allowed
%                    multiplier instead of a fixed multiplier
%   sigma            Augmented Lagrangian penalty parameter
%   AutoSigma        Flag determining whether sigma is automatically
%                    updated (see Sec. 3.4.1 of boyd-2010-distributed)
%   AutoSigmaPeriod  Iteration period on which sigma is updated
%   SigmaRsdlRatio   Primal/dual residual ratio in sigma update test
%   SigmaScaling     Multiplier applied to sigma when updated
%   AutoSigmaScaling Flag determining whether SigmaScaling value is
%                    adaptively determined (see wohlberg-2015-adaptive). If
%                    enabled, SigmaScaling specifies a maximum allowed
%                    multiplier instead of a fixed multiplier.
%   StdResiduals     Flag determining whether standard residual definitions
%                    (see Sec 3.3 of boyd-2010-distributed) are used instead
%                    of normalised residuals (see wohlberg-2015-adaptive)
%   XRelaxParam      Relaxation parameter (see Sec. 3.4.3 of
%                    boyd-2010-distributed) for X update
%   DRelaxParam      Relaxation parameter (see Sec. 3.4.3 of
%                    boyd-2010-distributed) for D update
%   LinSolve         Linear solver for main problem: 'SM' or 'CG'
%   MaxCGIter        Maximum CG iterations when using CG solver
%   CGTol            CG tolerance when using CG solver
%   CGTolAuto        Flag determining use of automatic CG tolerance
%   CGTolFactor      Factor by which primal residual is divided to obtain CG
%                    tolerance, when automatic tolerance is active
%   NoBndryCross     Flag indicating whether all solution coefficients
%                    corresponding to filters crossing the image boundary
%                    should be forced to zero.
%   DictFilterSizes  Array of size 2 x M where each column specifies the
%                    filter size (rows x columns) of the corresponding
%                    dictionary filter
%   NonNegCoef       Flag indicating whether solution should be forced to
%                    be non-negative
%   ZeroMean         Force learned dictionary entries to be zero-mean
%
%
% Author: Brendt Wohlberg <brendt@lanl.gov>  Modified: 2015-12-18
%
% This file is part of the SPORCO library. Details of the copyright
% and user license can be found in the 'License' file distributed with
% the library.

if nargin < 4,
  opt = [];
end
checkopt(opt, defaultopts([]));
opt = defaultopts(opt);

% Set up status display for verbose operation
hstr = ['Itn   Fnc       DFid      l1        Cnstr     '...
        'r(X)      s(X)      r(D)      s(D) '];
sfms = '%4d %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e';
nsep = 84;
if opt.AutoRho,
  hstr = [hstr '     rho  '];
  sfms = [sfms ' %9.2e'];
  nsep = nsep + 10;
end
if opt.AutoSigma,
  hstr = [hstr '     sigma  '];
  sfms = [sfms ' %9.2e'];
  nsep = nsep + 10;
end
if opt.Verbose && opt.MaxMainIter > 0,
  disp(hstr);
  disp(char('-' * ones(1,nsep)));
end

% Collapsing of trailing singleton dimensions greatly complicates
% handling of both SMV and MMV cases. The simplest approach would be
% if S could always be reshaped to 4d, with dimensions consisting of
% image rows, image cols, a single dimensional placeholder for number
% of filters, and number of measurements, but in the single
% measurement case the third dimension is collapsed so that the array
% is only 3d.
if size(S,3) > 1,
  xsz = [size(S,1) size(S,2) size(D0,3) size(S,3)];
  % Insert singleton 3rd dimension (for number of filters) so that
  % 4th dimension is number of images in input s volume
  S = reshape(S, [size(S,1) size(S,2) 1 size(S,3)]);
else
  xsz = [size(S,1) size(S,2) size(D0,3) 1];
end
Nx = prod(xsz);
Nd = prod(xsz(1:2))*size(D0,3);
cgt = opt.CGTol;

% Dictionary size may be specified when learning multiscale
% dictionary
if isempty(opt.DictFilterSizes),
  dsz = [size(D0,1) size(D0,2)];
else
  dsz = opt.DictFilterSizes;
end

% Mean removal and normalisation projections
Pzmn = @(x) bsxfun(@minus, x, mean(mean(x,1),2));
Pnrm = @(x) bsxfun(@rdivide, x, sqrt(sum(sum(x.^2, 1), 2)));

% Projection of filter to full image size and its transpose
% (zero-pad and crop respectively)
Pzp = @(x) zpad(x, xsz(1:2));
PzpT = @(x) bcrop(x, dsz);

% Projection of dictionary filters onto constraint set
if opt.ZeroMean,
  Pcn = @(x) Pnrm(Pzp(Pzmn(PzpT(x))));
else
  Pcn = @(x) Pnrm(Pzp(PzpT(x)));
end

% Start timer
tstart = tic;

% Project initial dictionary onto constraint set
D = Pnrm(D0);

% Compute signal in DFT domain
Sf = fft2(S);

% Set up algorithm parameters and initialise variables
rho = opt.rho;
if isempty(rho), rho = 50*lambda+1; end;
if opt.AutoRho,
  asgr = opt.RhoRsdlRatio;
  asgm = opt.RhoScaling;
end
sigma = opt.sigma;
if isempty(sigma), sigma = size(S,3); end;
if opt.AutoSigma,
  asdr = opt.SigmaRsdlRatio;
  asdm = opt.SigmaScaling;
end
optinf = struct('itstat', [], 'opt', opt);
rx = Inf;
sx = Inf;
rd = Inf;
sd = Inf;
eprix = 0;
eduax = 0;
eprid = 0;
eduad = 0;

% Initialise main working variables
X = [];
if isempty(opt.Y0),
  Y = zeros(xsz, class(S));
else
  Y = opt.Y0;
end
Yprv = Y;
if isempty(opt.U0),
  if isempty(opt.Y0),
    U = zeros(xsz, class(S));
  else
    U = (lambda/rho)*sign(Y);
  end
else
  U = opt.U0;
end
Df = [];
if isempty(opt.G0),
  G = Pzp(D);
else
  G = opt.G0;
end
Gprv = G;
if isempty(opt.H0),
  if isempty(opt.G0),
    H = zeros(size(G), class(S));
  else
    H = G;
  end
else
  H = opt.H0;
end
Gf = fft2(G, size(S,1), size(S,2));
GSf = bsxfun(@times, conj(Gf), Sf);

% Main loop
k = 1;
while k <= opt.MaxMainIter && (rx > eprix|sx > eduax|rd > eprid|sd >eduad),

  % Solve X subproblem. It would be simpler and more efficient (since the
  % DFT is already available) to solve for X using the main dictionary
  % variable D as the dictionary, but this appears to be unstable. Instead,
  % use the projected dictionary variable G
  Xf = solvedbi_sm(Gf, rho, GSf + rho*fft2(Y - U));
  X = ifft2(Xf, 'symmetric');
  clear Xf Gf GSf;

  % See pg. 21 of boyd-2010-distributed
  if opt.XRelaxParam == 1,
    Xr = X;
  else
    Xr = opt.XRelaxParam*X + (1-opt.XRelaxParam)*Y;
  end

  % Solve Y subproblem
  Y = shrink(Xr + U, (lambda/rho)*opt.L1Weight);
  if opt.NonNegCoef,
    Y(Y < 0) = 0;
  end
  if opt.NoBndryCross,
    %Y((end-max(dsz(1,:))+2):end,:,:,:) = 0;
    Y((end-size(D0,1)+2):end,:,:,:) = 0;
    %Y(:,(end-max(dsz(2,:))+2):end,:,:) = 0;
    Y(:,(end-size(D0,2)+2):end,:,:) = 0;
  end
  Yf = fft2(Y);
  YSf = sum(bsxfun(@times, conj(Yf), Sf), 4);

  % Update dual variable corresponding to X, Y
  U = U + Xr - Y;
  clear Xr;

  % Compute primal and dual residuals and stopping thresholds for X update
  nX = norm(X(:)); nY = norm(Y(:)); nU = norm(U(:));
  if opt.StdResiduals,
    % See pp. 19-20 of boyd-2010-distributed
    rx = norm(vec(X - Y));
    sx = norm(vec(rho*(Yprv - Y)));
    eprix = sqrt(Nx)*opt.AbsStopTol+max(nX,nY)*opt.RelStopTol;
    eduax = sqrt(Nx)*opt.AbsStopTol+rho*nU*opt.RelStopTol;
  else
    % See wohlberg-2015-adaptive
    rx = norm(vec(X - Y))/max(nX,nY);
    sx = norm(vec(Yprv - Y))/nU;
    eprix = sqrt(Nx)*opt.AbsStopTol/max(nX,nY)+opt.RelStopTol;
    eduax = sqrt(Nx)*opt.AbsStopTol/(rho*nU)+opt.RelStopTol;
  end
  clear X;

  % Compute l1 norm of Y
  Jl1 = sum(abs(vec(bsxfun(@times, opt.L1Weight, Y))));

  % Update record of previous step Y
  Yprv = Y;

  % Solve D subproblem. Similarly, it would be simpler and more efficient to
  % solve for D using the main coefficient variable X as the coefficients,
  % but it appears to be more stable to use the shrunk coefficient variable Y
  if strcmp(opt.LinSolve, 'SM'),
    Df = solvemdbi_ism(Yf, sigma, YSf + sigma*fft2(G - H));
  else
    [Df, cgst] = solvemdbi_cg(Yf, sigma, YSf + sigma*fft2(G - H), ...
                              cgt, opt.MaxCGIter, Df(:));
  end
  clear YSf;
  D = ifft2(Df, 'symmetric');
  if strcmp(opt.LinSolve, 'SM'), clear Df; end

  % See pg. 21 of boyd-2010-distributed
  if opt.DRelaxParam == 1,
    Dr = D;
  else
    Dr = opt.DRelaxParam*D + (1-opt.DRelaxParam)*G;
  end

  % Solve G subproblem
  G = Pcn(Dr + H);
  Gf = fft2(G);
  GSf = bsxfun(@times, conj(Gf), Sf);

  % Update dual variable corresponding to D, G
  H = H + Dr - G;
  clear Dr;

  % Compute primal and dual residuals and stopping thresholds for D update
  nD = norm(D(:)); nG = norm(G(:)); nH = norm(H(:));
  if opt.StdResiduals,
    % See pp. 19-20 of boyd-2010-distributed
    rd = norm(vec(D - G));
    sd = norm(vec(sigma*(Gprv - G)));
    eprid = sqrt(Nd)*opt.AbsStopTol+max(nD,nG)*opt.RelStopTol;
    eduad = sqrt(Nd)*opt.AbsStopTol+sigma*nH*opt.RelStopTol;
  else
    % See wohlberg-2015-adaptive
    rd = norm(vec(D - G))/max(nD,nG);
    sd = norm(vec(Gprv - G))/nH;
    eprid = sqrt(Nd)*opt.AbsStopTol/max(nD,nG)+opt.RelStopTol;
    eduad = sqrt(Nd)*opt.AbsStopTol/(sigma*nH)+opt.RelStopTol;
  end

  % Apply CG auto tolerance policy if enabled
  if opt.CGTolAuto && (rd/opt.CGTolFactor) < cgt,
    cgt = rd/opt.CGTolFactor;
  end

  % Compute measure of D constraint violation
  Jcn = norm(vec(Pcn(D) - D));
  clear D;

  % Update record of previous step G
  Gprv = G;

  % Compute data fidelity term in Fourier domain (note normalisation)
  Jdf = sum(vec(abs(sum(bsxfun(@times,Gf,Yf),3)-Sf).^2))/(2*xsz(1)*xsz(2));
  clear Yf;
  Jfn = Jdf + lambda*Jl1;

  % Record and display iteration details
  tk = toc(tstart);
  optinf.itstat = [optinf.itstat;...
       [k Jfn Jdf Jl1 rx sx rd sd eprix eduax eprid eduad rho sigma tk]];
  if opt.Verbose,
    dvc = [k, Jfn, Jdf, Jl1, Jcn, rx, sx, rd, sd];
    if opt.AutoRho,
      dvc = [dvc rho];
    end
    if opt.AutoSigma,
      dvc = [dvc sigma];
    end
    disp(sprintf(sfms, dvc));
  end

  % See wohlberg-2015-adaptive and pp. 20-21 of boyd-2010-distributed
  if opt.AutoRho,
    if k ~= 1 && mod(k, opt.AutoRhoPeriod) == 0,
      if opt.AutoRhoScaling,
        rhomlt = sqrt(rx/sx);
        if rhomlt < 1, rhomlt = 1/rhomlt; end
        if rhomlt > opt.RhoScaling, rhomlt = opt.RhoScaling; end
      else
        rhomlt = opt.RhoScaling;
      end
      rsf = 1;
      if rx > opt.RhoRsdlRatio*sx, rsf = rhomlt; end
      if sx > opt.RhoRsdlRatio*rx, rsf = 1/rhomlt; end
      rho = rsf*rho;
      U = U/rsf;
    end
  end
  if opt.AutoSigma,
    if k ~= 1 && mod(k, opt.AutoSigmaPeriod) == 0,
      if opt.AutoSigmaScaling,
        sigmlt = sqrt(rd/sd);
        if sigmlt < 1, sigmlt = 1/sigmlt; end
        if sigmlt > opt.SigmaScaling, sigmlt = opt.SigmaScaling; end
      else
        sigmlt = opt.SigmaScaling;
      end
      ssf = 1;
      if rd > opt.SigmaRsdlRatio*sd, ssf = sigmlt; end
      if sd > opt.SigmaRsdlRatio*rd, ssf = 1/sigmlt; end
      sigma = ssf*sigma;
      H = H/ssf;
    end
  end

  k = k + 1;

end

D = PzpT(G);

% Record run time and working variables
optinf.runtime = toc(tstart);
optinf.Y = Y;
optinf.U = U;
optinf.G = G;
optinf.H = H;
optinf.lambda = lambda;
optinf.rho = rho;
optinf.sigma = sigma;
optinf.cgt = cgt;
if exist('cgst'), optinf.cgst = cgst; end

if opt.Verbose && opt.MaxMainIter > 0,
  disp(char('-' * ones(1,nsep)));
end

return

function u = vec(v)

  u = v(:);

return

function u = shrink(v, lambda)

  if isscalar(lambda),
    u = sign(v).*max(0, abs(v) - lambda);
  else
    u = sign(v).*max(0, bsxfun(@minus, abs(v), lambda));
  end

return

function u = zpad(v, sz)

  u = zeros(sz(1), sz(2), size(v,3), size(v,4), class(v));
  u(1:size(v,1), 1:size(v,2),:,:) = v;

return

function u = bcrop(v, sz)

  if numel(sz) <= 2,
    if numel(sz) == 1
      cs = [sz sz];
    else
      cs = sz;
    end
    u = v(1:cs(1), 1:cs(2), :);
  else
    cs = max(sz,[],2);
    u = zeros(cs(1), cs(2), size(v,3), class(v));
    for k = 1:size(v,3),
      u(1:sz(1,k), 1:sz(2,k), k) = v(1:sz(1,k), 1:sz(2,k), k);
    end
  end

return

function opt = defaultopts(opt)

  if ~isfield(opt,'Verbose'),
    opt.Verbose = 0;
  end
  if ~isfield(opt,'MaxMainIter'),
    opt.MaxMainIter = 1000;
  end
  if ~isfield(opt,'AbsStopTol'),
    opt.AbsStopTol = 1e-6;
  end
  if ~isfield(opt,'RelStopTol'),
    opt.RelStopTol = 1e-4;
  end
  if ~isfield(opt,'L1Weight'),
    opt.L1Weight = 1;
  end
  if ~isfield(opt,'Y0'),
    opt.Y0 = [];
  end
  if ~isfield(opt,'U0'),
    opt.U0 = [];
  end
  if ~isfield(opt,'G0'),
    opt.G0 = [];
  end
  if ~isfield(opt,'H0'),
    opt.H0 = [];
  end
  if ~isfield(opt,'rho'),
    opt.rho = [];
  end
  if ~isfield(opt,'AutoRho'),
    opt.AutoRho = 0;
  end
  if ~isfield(opt,'AutoRhoPeriod'),
    opt.AutoRhoPeriod = 10;
  end
  if ~isfield(opt,'RhoRsdlRatio'),
    opt.RhoRsdlRatio = 10;
  end
  if ~isfield(opt,'RhoScaling'),
    opt.RhoScaling = 2;
  end
  if ~isfield(opt,'AutoRhoScaling'),
    opt.AutoRhoScaling = 0;
  end
  if ~isfield(opt,'sigma'),
    opt.sigma = [];
  end
  if ~isfield(opt,'AutoSigma'),
    opt.AutoSigma = 0;
  end
  if ~isfield(opt,'AutoSigmaPeriod'),
    opt.AutoSigmaPeriod = 10;
  end
  if ~isfield(opt,'SigmaRsdlRatio'),
    opt.SigmaRsdlRatio = 10;
  end
  if ~isfield(opt,'SigmaScaling'),
    opt.SigmaScaling = 2;
  end
  if ~isfield(opt,'AutoSigmaScaling'),
    opt.AutoSigmaScaling = 0;
  end
  if ~isfield(opt,'StdResiduals'),
    opt.StdResiduals = 0;
  end
  if ~isfield(opt,'XRelaxParam'),
    opt.XRelaxParam = 1;
  end
  if ~isfield(opt,'DRelaxParam'),
    opt.DRelaxParam = 1;
  end
  if ~isfield(opt,'LinSolve'),
    opt.LinSolve = 'SM';
  end
  if ~isfield(opt,'MaxCGIter'),
    opt.MaxCGIter = 1000;
  end
  if ~isfield(opt,'CGTol'),
    opt.CGTol = 1e-3;
  end
  if ~isfield(opt,'CGTolAuto'),
    opt.CGTolAuto = 0;
  end
  if ~isfield(opt,'CGTolAutoFactor'),
    opt.CGTolFactor = 50;
  end
  if ~isfield(opt,'NoBndryCross'),
    opt.NoBndryCross = 0;
  end
  if ~isfield(opt,'DictFilterSizes'),
    opt.DictFilterSizes = [];
  end
  if ~isfield(opt,'NonNegCoef'),
    opt.NonNegCoef = 0;
  end
  if ~isfield(opt,'ZeroMean'),
    opt.ZeroMean = 0;
  end

return

author: thfylsty


text:
function truefalse = ch_is_digit ( c )

%*****************************************************************************80
%
% CH_IS_DIGIT returns TRUE if the character C is a digit.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 December 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character C, a character.
%
%    Output, integer TRUEFALSE, is TRUE (1) if C is a digit, FALSE (0) otherwise.
%
  TRUE = 1;
  FALSE = 0;

  if ( '0' <= c & c <= '9' )
    truefalse = TRUE;
  else
    truefalse = FALSE;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function VectorFitting_Sparam ()
clear; close all; FONTSIZE=20;

file='../DC_ringmod_type1_R=10,gap=180,Lc=0,wg=500,lambda=1550,mesh=2,angle=30.mat';
load (file)

c=3e8; wavelength=c./f*1e6;
fSCALING = 1e14; f=f/fSCALING; s=1i*2*pi*f;  Np=length(f);
hwait = waitbar(0,'Please wait...');

% average the S parameters for the symmetric parameters that were simulated twice,
% by considering amplitude and phase separately
S1221 = (abs(S12)+abs(S21))/2 .* exp ( 1i * (unwrap(angle(S12)) + unwrap(angle(S21)) ) /2 );
S12=S1221; S21 = S1221;
S4132 = (abs(S41)+abs(S32))/2 .* exp ( 1i * (unwrap(angle(S41)) + unwrap(angle(S32)) ) /2 );
S41=S4132; S32 = S4132;
S13=S31; S23=S41; S33=S11; S43=S21; S14=S32; S24=S42; S34=S12; S44=S22;

for i=1:Np
  Sparam = [ [ S11(i),S12(i),S13(i),S14(i)];
    [ S21(i),S22(i),S23(i),S24(i)];
    [ S31(i),S32(i),S33(i),S34(i)];
    [ S41(i),S42(i),S43(i),S44(i)] ] ;
  Test1(i) = norm(Sparam);
  Sparam_w (:,:,i)=Sparam;
end

% Check if S-Parameters are already passive; if not, perform Vector Fit,
% otherwise, export.
if ~isempty(find(Test1>1))

  % rational fit, sweep number of parameters, N:
  optsN=20:1:100;
  opts.poletype='lincmplx'; opts.parametertype='S'; opts.stable=0; opts.Niter_out=10;
  Npassive=[]; rms3passive=[];
  for i=1:length(optsN)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    opts.N=optsN(i);    % Vector Fit:
    [SER,rmserr,Hfit,opts2]=VFdriver(Sparam_w,s,[],opts);

    %%%%%%%%%%%%%%%%%%%%%%  Enforce Passivity:
    [SER,H_passive,opts3,wintervals]=RPdriver(SER,s,opts);
    % Note: added output parameter wintervals to RPdriver.

    % Find rms error:
    tell=0; Nc=length(SER.D);
    for col=1:Nc
      for row=col:Nc % makes assumption that S = S'
        tell=tell+1; % make a single vector:
        Sparam10(tell,:)= squeeze(Sparam_w(row,col,:)).';
        fit10(tell,:)= squeeze(Hfit(row,col,:)).';
        fitP(tell,:)= squeeze(H_passive(row,col,:)).';
      end
    end
    rms2(i) = sqrt(sum(sum(abs((Sparam10-fit10).^2))))/sqrt(4*Np);
    rms3(i) = sqrt(sum(sum(abs((Sparam10-fitP).^2))))/sqrt(4*Np);

    % Determine if the Passivity Enforcement was successful.
    if (rms3(i) < 1) && isempty(wintervals)
      Npassive=[Npassive optsN(i)]; rms3passive = [rms3passive rms3(i)];
    end

    waitbar (i/length(optsN), hwait);
    if (rms3(i) < 1e-4) && isempty(wintervals) ; break; end
  end
  close (hwait);

  % Plot error versus number of fitting parameters:
  figure;
  semilogy(optsN(1:i),rms2,'ro-', 'LineWidth',2, 'MarkerSize',7); hold all;
  labels={}; labels{end+1} = 'Vector Fit, rms';
  plot (Npassive, rms3passive, 'kx', 'MarkerSize',14, 'LineWidth',3);
  labels{end+1} = 'Passivity-Enforced Fit, rms';
  xlabel('Fitting order, N');   ylabel('rms error');
  for ii=1:length(labels); labels{ii}=[labels{ii} ' ' char(31) ]; end;
  legend (labels,'Location','NorthEast');
  printfig(file,'convergence');

  % Plot residue values:
  figure;
  residues=sort(squeeze(prod(prod(abs(SER.R),1),2).^(1/4)), 1, 'descend');
  semilogy(residues,'-s', 'LineWidth',3, 'MarkerSize',8);
  xlabel('Residue index'); ylabel('Residue magnitude');
  printfig(file,'residues');

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Plot S-Parameters + fit functions:
  figure;
  Nc=length(SER.D);
  for row=1:Nc
    for col= row:Nc
      dum1=squeeze(Sparam_w(row,col,:));
      dum2=squeeze(H_passive(row,col,:));
      h1=semilogy(wavelength,abs(dum1),'b','LineWidth',3); hold on
      h2=semilogy(wavelength,abs(dum2),'r-.','LineWidth',4);
      h3=semilogy(wavelength,abs(dum2-dum1),'g--','LineWidth',3);
    end
  end
  hold off
  xlabel('Wavelength'); ylabel('Amplitude [S]');
  axis tight; Yl =ylim; ylim ([Yl(1)*10 Yl(2)*2]);
  labels={}; labels{end+1} = 'FDTD S-Parameters';
  labels{end+1} = 'Passivity-enforced Fit'; labels{end+1}='Deviation';
  for i=1:length(labels); labels{i}=[labels{i} ' ' char(31)]; end;
  legend (labels,'Location','Best');
  printfig(file,'VF2a');

  figure;
  Nc=length(SER.D);
  for row= 1:Nc  %1
    for col= row:Nc  %3
      dum1=squeeze(Sparam_w(row,col,:));
      dum2=squeeze(H_passive(row,col,:));
      h1=semilogy(wavelength,unwrap(angle(dum1)),'b','LineWidth',3); hold on
      h2=plot(wavelength,unwrap(angle(dum2)),'r-.','LineWidth',4);
      h3=plot(wavelength,abs(unwrap(angle(dum2))-unwrap(angle(dum1))),'g--','LineWidth',3);
    end
  end
  hold off
  xlabel('Wavelength'); ylabel('Phase [S]');
  axis tight; Yl =ylim; ylim ([Yl(1)*10 Yl(2)*2]);
  labels={}; labels{end+1} = 'FDTD S-Parameters';
  labels{end+1} = 'Passivity-enforced Fit'; labels{end+1}='Deviation';
  for i=1:length(labels); labels{i}=[labels{i} ' ' char(31)]; end;
  legend (labels,'Location','Best');
  printfig(file,'VF2p');

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Passivity test results:
  for i=1:length(f)
    Test0(i) = norm(Sparam_w(:,:,i));
    Test1(i) = norm(Hfit(:,:,i));
    Test2(i) = norm(H_passive(:,:,i));
  end
  figure;
  plot (wavelength,Test0,'LineWidth',2); hold all;
  plot (wavelength,Test1,'--','LineWidth',2)
  plot (wavelength,Test2,'LineWidth',4)
  plot(wavelength,ones(length(f),1),'--','LineWidth',3);
  labels={}; labels{1} = 'FDTD S-Parameters'; labels{2} = 'Rational Fit'; labels{3} = 'Passivity-enforced Fit'; labels{4}='Passivity limit';
  for i=1:length(labels); labels{i}=[labels{i} ' ' char(31)]; end; legend (labels,'Location','Best');
  axis tight
  xlabel ('Wavelength [\mum]');
  ylabel ('S-Parameter Passivity Test');
  printfig(file,'passivitytest3');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% export S parameters to INTERCONNECT
fid = fopen([file '.sparam'],'w');  Nc=4;
for row= 1:Nc
  for col= 1:Nc
    fprintf(fid,'%s\n',[ '(''port ' num2str(col) ''',''TE'',1,''port ' num2str(row) ''',1,''transmission'')' ]  );
    fprintf(fid,'%s\n',[ '(' num2str(Np) ',3)' ]  );
    dum2=squeeze(H_passive(row,col,:));
    mag = abs(dum2);
    phase = unwrap (angle(dum2)); % figure; plot (wavelength, phase);
    for i=1:Np
      fprintf(fid,'%g	%g	%g\n', f(i)*1e14, mag(i), phase(i));
    end
  end
end
fclose(fid);

function printfig (file, b)
global PRINT_titles;
PRINT_titles=0;
FONTSIZE=20;
set(get(gca,'xlabel'),'FontSize',FONTSIZE);
set(get(gca,'ylabel'),'FontSize',FONTSIZE);
set(get(gca,'title'),'FontSize',FONTSIZE-5);
set(gca,'FontSize',FONTSIZE-2);
if PRINT_titles==0
  delete(get(gca,'title'))
end
%a=strfind(file,'.'); file(a)=',';
pdf = [file(1:end-4) '_' b '.pdf'];
print ('-dpdf','-r300', pdf);
system([ 'pdfcrop ' pdf ' ' pdf ' &' ]);
% system(['acroread ' pdf '.pdf &']);

author: lukasc-ubc

repo: SiliconPhotonicsDesign


text:
function p = p25_p_get ( )

%*****************************************************************************80
%
%% P25_P_GET gets the value of a parameter for problem 25.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real P, the value of the parameter that multiplies the
%    Gaussian noise.
%
  p = 0.0;
  p = p25_p_val ( 'GET', p );

  return
end


text:
function a = zero_eigen_right ( n )

%*****************************************************************************80
%
%% ZERO_EIGEN_RIGHT returns the right eigenvectors of the ZERO matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, N, the order of the matrix.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for i = 1 : n
    for j = 1 : n
      if ( i == j )
        a(i,j) = 1.0;
      else
        a(i,j) = 0.0;
      end
    end
  end

  return
end


text:
function model = setupComponentContribution(model,molFileDir,cid,printLevel)
% Estimates standard transformed reaction Gibbs energy and directionality
% at in vivo conditions in multicompartmental metabolic reconstructions.
% Has external dependencies on the COBRA toolbox, the component
% contribution method, Python (with numpy and Open Babel bindings),
% ChemAxon's Calculator Plugins, and Open Babel. See details on
% availability at the end of help text.
%
% modelT = setupThermoModel(model,molfileDir,cid,T,cellCompartments,ph,...
%                           is,chi,concMin,concMax,confidenceLevel)
%
% INPUTS
% model             Model structure with following fields:
% .S                m x n stoichiometric matrix.
% .mets             m x 1 array of metabolite identifiers.
% .rxns             n x 1 array of reaction identifiers.
% .metFormulas      m x 1 cell array of metabolite formulas. Formulas for
%                   protons should be H, and formulas for water should be
%                   H2O.
% .metCharges       m x 1 numerical array of metabolite charges.
%
% OPTIONAL INPUTS
% molFileDir                Path to a directory containing molfiles for the
%                           major tautomer of the major microspecies of
%                           each metabolite at pH 7. Molfiles should be
%                           named with the metabolite identifiers in
%                           model.mets (without compartment assignments).
%                           Not required if cid are specified.
% cid                       m x 1 cell array of KEGG Compound identifiers.
%                           Not required if molfiledir is specified.
% model.metCompartments     m x 1 array of metabolite compartment
%                           assignments. Not required if metabolite
%                           identifiers are strings of the format ID[*]
%                           where * is the appropriate compartment
%                           identifier.
%
% OUTPUTS
% model                 Model structure with following additional fields:
% .inchi                Structure containing four m x 1 cell array's of
%                       IUPAC InChI strings for metabolites, with varying
%                       levels of structural detail.
% .pKa                  m x 1 structure containing metabolite pKa values
%                       estimated with ChemAxon's Calculator Plugins.
% .pseudoisomers        p x 4 matrix with the following columns:
%                       1. Metabolite index.
%                       2. Estimated pseudoisomer standard Gibbs energy.
%                       3. Number of hydrogen atoms in pseudoisomer
%                       chemical formula.
%                       4. Charge on pseudoisomer.
%
% WRITTEN OUTPUTS
% MetStructures.sdf     An SDF containing all structures input to the
%                       component contribution method for estimation of
%                       standard Gibbs energies.
%
%
% Ronan M. T. Fleming, Sept. 2012   Version 1.0
% Hulda S. H., Dec. 2012            Version 2.0

%% Configure inputs
% Retreive molfiles from KEGG if KEGG ID are given. Otherwise use molfiles
% in molfileDir.
if ~exist('cid','var')
    cid = [];
end
if ~exist('printLevel','var')
    printLevel = 1;
end
%% Get metabolite structures
if ~isempty(cid)
    molFileDir = 'molfilesFromKegg';
    fprintf('\nRetreiving molfiles from KEGG.\n');
    takeMajorMS = true; % Convert molfile from KEGG to major tautomer of major microspecies at pH 7
    pH = 7;
    takeMajorTaut = true;
    kegg2mol(cid,molFileDir,model.mets,takeMajorMS,pH,takeMajorTaut); % Retreive mol files
end

if printLevel>0
fprintf('Creating MetStructures.sdf from molfiles.\n')
end

sdfFileName = 'MetStructures.sdf';
includeRs = 0; % Do not include structures with R groups in SDF
[sdfMetList,noMolMetList] = mol2sdf(model.mets,molFileDir,sdfFileName,includeRs);

if printLevel>0
fprintf('Converting SDF to InChI strings.\n')
end
model.inchi = createInChIStruct(model.mets,sdfFileName);
compositeBool = ~cellfun('isempty',regexp(model.inchi.nonstandard,'\.')); % Remove InChI for composite compounds as they cause problems later.
model.inchi.standard(compositeBool) = cell(sum(compositeBool),1);
model.inchi.standardWithStereo(compositeBool) = cell(sum(compositeBool),1);
model.inchi.standardWithStereoAndCharge(compositeBool) = cell(sum(compositeBool),1);
model.inchi.nonstandard(compositeBool) = cell(sum(compositeBool),1);

%% Estimate metabolite pKa values with ChemAxon calculator plugins and determine all relevant pseudoisomers.
if printLevel>0
fprintf('Estimating metabolite pKa values.\n');
end
npKas = 20; % Number of acidic and basic pKa values to estimate
takeMajorTaut = false; % Estimate pKa for input tautomer. Input tautomer is assumed to be the major tautomer for the major microspecies at pH 7.
model.pseudoisomers = estimate_pKa(model.mets,model.inchi.nonstandard,npKas,takeMajorTaut); % Estimate pKa and determine pseudoisomers
model.pseudoisomers = rmfield(model.pseudoisomers,'met');

% Add number of hydrogens and charge for metabolites with no InChI
if any(~[model.pseudoisomers.success]) && printLevel>0
    fprintf('Assuming that metabolite species in model.metFormulas are representative for metabolites where pKa could not be estimated.\n');
end
nonphysicalMetSpecies = {};
for i = 1:length(model.mets)
    model_z = model.metCharges(i); % Get charge from model
    model_nH = numAtomsOfElementInFormula(model.metFormulas{i},'H'); % Get number of hydrogens from metabolite formula in model
    if ~model.pseudoisomers(i).success
        model.pseudoisomers(i).zs = model_z;
        model.pseudoisomers(i).nHs = model_nH;
        model.pseudoisomers(i).majorMSpH7 = true; % Assume species in model is the major (and only) metabolite species %RF: this seems dubious
    end
    if ~any(model.pseudoisomers(i).nHs == model_nH)
        nonphysicalMetSpecies = [nonphysicalMetSpecies; model.mets(i)];
    end
end
if ~isempty(nonphysicalMetSpecies)
    nonphysicalMetSpecies = unique(regexprep(nonphysicalMetSpecies,'\[\w\]',''));
    if printLevel>1
        fprintf('%s\n','#H in model.metFormulas does not match any of the species calculated mol file for metabolites:')
        for n=1:length(nonphysicalMetSpecies)
            fprintf('%s\n',nonphysicalMetSpecies{n});%,model.metFormulas{m});
        end
    end
end

author: opencobra

repo: cobratoolbox


text:
%finds best source and distance to the best source
% if distance is negative, the best source cannot be found (for example, because the propagation was stopped before it reached this point)
% Danil Kirsanov, 09/2007

function [source_id, distance] = geodesic_best_source(algorithm, destination)

global geodesic_library;

if nargin == 2
    d = geodesic_convert_surface_points({destination});

    tmp = 1;
    [source_id, tmp1, distance] = calllib(geodesic_library, 'distance_and_source', algorithm.id, d, tmp);
    source_id = source_id + 1;
else                                    %return distances and sources for all vertices
    tmp = libpointer('doublePtrPtr');
    tmp1 = libpointer('int32PtrPtr');

    [num_vertices, d, s] = calllib(geodesic_library, 'distance_and_source_for_all_vertices', algorithm.id, tmp, tmp1);

    setdatatype(d, 'doublePtr', num_vertices);
    distance = d.Value;
    setdatatype(s, 'int32Ptr', num_vertices);
    source_id = s.Value + 1;
end

author: Sable

repo: mcbench-benchmarks


text:
function Index = SVM_AFdetection_withoutTrainingModel(test_feature,test_class)
%******************************************************
% $ This function is used for calculating the SVM-based AF detection
% indices using the trained 9-fold models from MIT AF database. The
% determination of the AF is based on the voting resuts from the 9-fold
% models
%
% $ Reference: Q. Li, C. Y. Liu, J. Oster and G. D. Clifford. Chapter:
% Signal processing and feature selection preprocessing for classification
% in noisy healthcare data. In book: Machine Learning for Healthcare
% Technologies, Edition: 1st, Publisher: IET, Editors: David A. Clifton,
% 2016.
%
% Roberta Collaca, Julian Oster, Clifford
%
% $ Variable declaration: Input:
%   test_feature: the feature metrix from AF detection, you need to run
%       AF_feature function first to obtain the AF features, each row is a
%       feature vector from an RR interval time series
%   test_class: the labels of AF for the tested RR ineterval time series,
%       N*1 metrix, 1 for AF and 0 for non-AF
% Output:
%   Index: AF detection indices output, with the defined order of indices
%   of TP FN FP TN Se Sp Acc PPV NPV J in turn.
%
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information

if nargin<2
    test_class = ones(length(test_feature),1);
end

% if strcmp(computer,'GLNXA64')
%     addpath('/home/adriana/libsvm-320/matlab');
% elseif strcmp(computer, 'MACI64')
%     addpath('/Users/adri/Documents/MATLAB/libsvm-3.20/libsvm-3.20/matlab/');
% end

model_all = load('models.mat');
meanx_all = load('means.mat');
stdx_all = load('stds.mat');

select=[2 3 4 5 6 7 12 14];
test_feature=test_feature(:,select);

predict_all=[];
for k = 1:9
    model=model_all.model_all{k};
    meanx=meanx_all.meanx_all{k};
    stdx=stdx_all.stdx_all{k};
    xtest=bsxfun(@minus, test_feature, meanx);
    xtest=bsxfun(@rdivide, xtest,stdx);
    [predict_all(:,k), accuracy_test, ytest_out] = svmpredict(test_class, xtest, model, '-b 1');
end

vote=sum(predict_all');
vote(find(vote<5))=0;
vote(find(vote>=5))=1;
vote=vote';

TP=length(find(test_class==1 & vote==1));
FP=length(find(test_class==0 & vote==1));
FN=length(find(test_class==1 & vote==0));
TN=length(find(test_class==0 & vote==0));
Se=TP/(TP+FN);
Sp=TN/(TN+FP);
Acc=(TP+TN)/length(test_class);
PPV=TP/(TP+FP);
NPV=TN/(TN+FN);
J=Se+Sp-1;
% Index=[TP FN FP TN Se Sp Acc PPV NPV J];

% Edited by Adriana N Vest
% Changed output Index to include only AF prediction, not accuracy
Index = [vote];
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function [X,Y,Out] = TMac(data,known,Nway,coreNway,opts)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% tensor completion by parallel matrix factorization
% ==============================================================================
% solve
%  min_{X,Y} sum_{n=1}^N 0.5*alpha(n)*||Proj_{known}(Xn*Yn-Unfold_n(M))||_F^2
% ===============================================================================
%
% Input:
%       data: observed entries of the underlying tensor
%       known: indices of observed entries
%       Nway: the dimension of the underlying tensor
%       coreNway: estimated ranks of all mode matricizations
%       opts.
%           maxit: maximum number of iterations (default: 500)
%           tol: stopping tolerance (default: 1e-4)
%           maxT: maximum running time (sec) (default: 1e6)
%           alpha: weights in the model (default: alpha(n) = 1/N, any n)
%           alpha_adj: determine whether dynamically update alpha
%                       (default: 1 (yes))
%           rank_adj: determine rank-adjusting strategy
%                       (1: increase; -1: decrease; 0: fix)
%           rank_inc: rank increment if rank-increasing strategy is used
%           rank_min: minimum rank estimation
%           rank_max: maximum rank estimation
%
% Output:
%       X,Y: cell structs
%       Out.
%           iter: number of iterations
%           relerr1: relative change array of total fitting
%           relerr2: array of total fitting
%           alpha: final weights alpha (may be different from input)

N = length(Nway); coNway = zeros(1,N);
for n = 1:N
    coNway(n) = prod(Nway)/Nway(n);
end
nrmb = norm(data);
%% Parameters and defaults
if isfield(opts,'maxit')      maxit = opts.maxit;     else maxit = 500;             end
if isfield(opts,'tol')        tol = opts.tol;         else tol = 1e-4;              end
if isfield(opts,'maxT')       maxT = opts.maxT;       else maxT = 1e6;              end
if isfield(opts,'alpha')      alpha = opts.alpha;     else alpha = ones(1,N)/N;     end
if isfield(opts,'alpha_adj')
    alpha_adj = opts.alpha_adj;
else
    alpha_adj = 1;
end

if isfield(opts,'rank_adj') rank_adj = opts.rank_adj; else rank_adj = zeros(1,N);   end
if isfield(opts,'rank_inc') rank_inc = opts.rank_inc; else rank_inc = ones(1,N);    end
if isfield(opts,'rank_min') rank_min = opts.rank_min; else rank_min = ones(1,N);    end
if isfield(opts,'rank_max') rank_max = opts.rank_max; else rank_max = 50*ones(1,N); end

%% Data preprocessing and initialization
if isfield(opts,'X0')
    X = opts.X0;
else
    X = cell(1,N);
    for n = 1:N
        X{n} = randn(Nway(n),coreNway(n));
    end
end

if isfield(opts,'Y0')
    Y = opts.Y0;
else
    Y = cell(1,N);
    for n = 1:N
        Y{n} = randn(coreNway(n),coNway(n));
    end
end

% rescale the initial point based on number of elements
estMnrm = sqrt(nrmb^2*(prod(Nway)/length(known)));

for n = 1:N
    X{n} = X{n}/norm(X{n},'fro')*estMnrm^(Nway(n)/(Nway(n)+coNway(n)));
    Y{n} = Y{n}/norm(Y{n},'fro')*estMnrm^(coNway(n)/(Nway(n)+coNway(n)));
end

X0 = X; Y0 = Y;

[known,id] = sort(known); data = data(id);

M = zeros(Nway); M(known) = data;

sx = cell(1,N);
reschg = ones(1,N);
reschg_tol = max(1e-2,10*tol);
rank_inc_num = sum(rank_adj==1);

res0 = zeros(1,N); TotalRes = 0;
res = res0;
for n = 1:N
    Mn = Fold(X{n}*Y{n},Nway,n);
    res0(n) = norm(Mn(known)-data);
    TotalRes = TotalRes+res0(n);
end
solX = ones(1,N);
Xsq = cell(1,N); Yt = cell(1,N); spI = cell(1,N);
for n = 1:N
    Yt{n} = Y{n}';
end

Out.rank = coreNway;

start_time = tic;

alpha = alpha/sum(alpha);

fprintf('Iteration:     ');

for k = 1:maxit
    fprintf('\b\b\b\b\b%5i',k);

    % update (X,Y)
    for n = 1:N
        if alpha(n) > 0
            Mn = Unfold(M,Nway,n);
            if solX(n)
                X{n} = Mn*Yt{n};
            end
            solX(n) = 1;

            Xsq{n} = X{n}'*X{n};
            Y{n} = pinv(Xsq{n})*X{n}'*Mn;
            Yt{n} = Y{n}';

            if rank_adj(n) == -1
                %rank_dec_adaptive();
                rank_dec_adaptive(Xsq,n,coreNway,rank_min,rank_adj,X,X0,Y,Y0);
            end
        end
    end

    % update M
    Mn = Fold(X{1}*Y{1},Nway,1);
    res(1) = norm(Mn(known)-data);
    M = alpha(1)*Mn;
    for n = 2:N
        if alpha(n) > 0
        Mn = Fold(X{n}*Y{n},Nway,n);
        res(n) = norm(Mn(known)-data);
        M = M+alpha(n)*Mn;
        end
    end

    % pass the true tensor M for evaluation
    if isfield(opts,'Mtr')
        Out.truerel(k) = norm(M(:)-opts.Mtr(:))/norm(opts.Mtr(:));
    end

    M(known) = data;

    TotalRes0 = TotalRes;
    TotalRes = 0;
    for n = 1:N
        if alpha(n) > 0
            TotalRes = TotalRes+res(n)^2;
        end
    end
    ratio = res./res0;   reschg = abs(1-ratio);

    if rank_inc_num > 0
        for n = 1:N
            if alpha(n) > 0
            if coreNway(n) < rank_max(n) && reschg(n) < reschg_tol
                %rank_inc_adaptive();
                rank_inc_adaptive(M,Y,n,Yt,Y0,coreNway,rank_max,rank_inc,rank_inc_num,nstall,X0,solX,Nway,coNway);
            end
            end
        end
    end

    % adaptively update weight

    if alpha_adj ~= 0
        alpha = 1./(res.^2); alpha = alpha/sum(alpha);
    end

    % record how the rank estimates are updated
    Out.rank = [Out.rank; coreNway];

    % --- diagnostics, reporting, stopping checks ---
    relerr1 = abs(TotalRes-TotalRes0)/(TotalRes0+1);
    relerr2 = sum(alpha.*res)/nrmb;

    % reporting
    Out.hist_rel(1,k) = relerr1;
    Out.hist_rel(2,k) = relerr2;

    % check stopping criterion
    crit = relerr1<tol;
    if crit; nstall = nstall+1; else nstall = 0; end
    if nstall>=3 || relerr2<tol break; end
    if toc(start_time)>maxT; break; end;

    X0 = X; Y0 = Y; M0 = M;
    res0 = res;
end % main
fprintf('\n'); Out.iter = k;
Out.alpha = alpha;

end

author: andrewssobral


text:
function tests = test_imBoxFilter(varargin)
% Test suite for function imBoxFilter.
%
%   Test case for the file imBoxFilter
%
%   Example
%   test_imBoxFilter
%
%   See also
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2018-11-22,    using Matlab 8.6.0.267246 (R2015b)
% Copyright 2018 INRA - Cepia Software Platform.

tests = functiontests(localfunctions);

function test_SquareUInt8_5x5(testCase) %#ok<*DEFNU>
% Test call of function without argument

img = zeros([20, 20], 'uint8');
img(6:15, 6:15) = 255;
dims = [5 5];

res = imBoxFilter(img, dims);

assertEqual(testCase, size(res), size(img));
assertEqual(testCase, class(res), class(img));

function test_SquareFloat_5x5(testCase) %#ok<*DEFNU>
% Test call of function without argument

img = zeros([20, 20], 'uint8');
img(6:15, 6:15) = 255;
dims = [5 5];

resRef = imBoxFilter(img, dims);
res = uint8(imBoxFilter(double(img), dims));

assertEqual(testCase, res, resRef);

function test_CubeFloat_7x5x3(testCase) %#ok<*DEFNU>
% Test call of function without argument

img = zeros([40 40 40], 'uint8');
img(11:30, 11:30, 11:30) = 255;
dims = [7 5 3];

resRef = imBoxFilter(img, dims);
res = uint8(imBoxFilter(double(img), dims));

assertEqual(testCase, res, resRef);


text:
function value = p13_f ( dim_num, point_num, x )

%*****************************************************************************80
%
%% P13_F evaluates the integrand for problem 13.
%
%  Dimension:
%
%    DIM_NUM arbitrary.
%
%  Region:
%
%    0 <= X(1:DIM_NUM) <= 1
%
%  Integrand:
%
%    product ( 1 <= i <= dim_num ) t(n(i))(2*x(i)-1)
%
%    where T(N)(X) is the Chebyshev polynomial of order N,
%    and N(I) = mod ( i, 4 ) + 1.
%
%  Exact Integral:
%
%    0
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Harald Niederreiter,
%    Implementation and Tests of Low-Discrepancy Sequences,
%    ACM Transactions on Modeling and Computer Simulation,
%    Volume 2, Number 3, July 1992, pages 195-213.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the argument.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, real X(DIM_NUM,POINT_NUM), the evaluation points.
%
%    Output, real VALUE(POINT_NUM), the integrand values.
%
  value(1:point_num) = 1.0;

  for point = 1 : point_num

    for dim = 1 : dim_num

      t = 2.0 * x(dim,point) - 1.0;
      k = mod ( dim, 4 );

      if ( k == 1 )
        factor = t;
      elseif ( k == 2 )
        factor = 2.0 * t - 1.0;
      elseif ( k == 3 )
        factor = ( 4.0 * t - 3.0 ) * t;
      elseif ( k == 4 )
        factor = ( 8.0 * t - 8.0 * t + 1.0 );
      end

      value(point) = value(point) * factor;

    end

  end

  p13_i4 ( 'I', '#', point_num );

  return
end


text:
function [scanOut3M,bbox3M] = cropAroundCOM(scanNum,strName,outputImgSizeV,...
                              bkgVal,planC)
% cropAroundCOM.m
% Crop input scan to specified dimensions around center of mass.
%--------------------------------------------------------------------------
% INPUTS
% scanNum            : Input scan no.
% strName            : Structure to crop around
% outputImgSizeV     : Output dimensions [rows, cols, slcs]
% bkgVal             : User-input intensity assigned to voxels outside
%                      strName. Leave empty to skip.
% planC
%--------------------------------------------------------------------------
% AI 04/22/22

%% Get input scan array
indexS = planC{end};
scan3M = double(getScanArray(scanNum,planC));
CToffset = double(planC{indexS.scan}(scanNum).scanInfo(1).CTOffset);
scan3M = scan3M - CToffset;
origSizV = size(scan3M);

if length(outputImgSizeV)==2
    outputImgSizeV(3) = size(scan3M,3);
end

%% Get structure mask
strNum = getMatchingIndex(strName,...
    {planC{indexS.structures}.structureName},'exact');
assocScanV = getStructureAssociatedScan(strNum,planC);

% Check for associated scan num
if isempty(strNum)
    error('Structure ''%s'' not found.',strName)
else
    strNum = strNum(assocScanV == scanNum);
    if isempty(strNum)
        error('Structure ''%s'' not found.',strName);
    end
end
cropStr3M = getStrMask(strNum,planC);

%Assign bkg intensity
if ~isempty(bkgVal)
   scan3M(~cropStr3M) = bkgVal;
end

%% Crop around COM

%Compute COM (3D)
bbox3M = false(origSizV);
%[x3M, y3M, z3M] = ndgrid(1:origSizV(1), 1:origSizV(2), 1:origSizV(3));
%COMv = round(mean([x3M(cropStr3M), y3M(cropStr3M), z3M(cropStr3M)]));
[rV, cV, sV] = find3d(cropStr3M);
COMv = round(mean([rV; cV; sV],2));

%Compute crop extents around COM
rStart = max(1,COMv(1)-outputImgSizeV(1)/2);
cStart = max(1,COMv(2)-outputImgSizeV(2)/2);

rEnd = min(COMv(1)+outputImgSizeV(1)/2, origSizV(1));
cEnd = min(COMv(2)+outputImgSizeV(2)/2, origSizV(2));

if ~isequal(outputImgSizeV(3),origSizV(3))
    sEnd = min(COMv(3)+outputImgSizeV(3)/2,origSizV(3));
    sStart = max(1,COMv(3)-outputImgSizeV(3)/2);
else
    sStart = 1;
    sEnd = outputImgSizeV(3) + 1;
end

%% Calc. padding required for specified output dimensions where needed
cropSizV = [rEnd - rStart, cEnd - cStart, sEnd - sStart];
if cropSizV(1)<outputImgSizeV(1)
    xPad = floor(outputImgSizeV(1)/2 - cropSizV(1)/2);
else
    xPad = 1;
end
if cropSizV(2)<outputImgSizeV(2)
    yPad = floor(outputImgSizeV(2)/2 - cropSizV(2)/2);
else
    yPad = 1;
end
if cropSizV(3)<outputImgSizeV(3)
    zPad = floor(outputImgSizeV(3)/2 - cropSizV(3)/2);
else
    zPad = 1;
end

%% Populate output scan
% Initialize output scan array
cornerCube = scan3M(1:5,1:5,1:5);
bgMean = mean(cornerCube(:));
scanOut3M = bgMean*ones([outputImgSizeV(1:2),size(scan3M,3)]);

%Populate with cropped scan
bbox3M(rStart:rEnd-1,cStart:cEnd-1,sStart:sEnd-1) = true;
scanCropM = scan3M(rStart:rEnd-1,cStart:cEnd-1,sStart:sEnd-1);
scanOut3M(xPad:xPad+cropSizV(1)-1,yPad:yPad+cropSizV(2)-1,...
    zPad:zPad+cropSizV(3)-1) = scanCropM;

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function boundary = p03_boundary_nearest ( m, n, point )

%*****************************************************************************80
%
%% P03_BOUNDARY_NEAREST returns a nearest boundary point in problem 03.
%
%  Discussion:
%
%    The given input point need not be inside the region.
%
%    In some cases, more than one boundary point may be "nearest",
%    but only one will be returned.
%
%    31 August 2005: Thanks to Hua Fei for pointing out that a previous
%    version of this routine gave inaccurate results for points that were
%    significantly far from the box.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real POINT(M,N), the points.
%
%    Output, real BOUNDARY(M,N), points on the boundary
%    that are nearest to each point.
%
  center = [ 0.0, 0.0 ];
  r1 = 1.0;
  r2 = 0.4;
  x1 = -1.0;
  x2 = +1.0;
  y1 = -1.0;
  y2 = +1.0;

  for j = 1 : n

    x = point(1,j);
    y = point(2,j);
%
%  Special case: a point at the center of the box.
%  The closest point is ANY point on the circle.
%
    if ( point(1,j) == center(1) & point(2,j) == center(2) )
      boundary(1,j) = center(1) + r2;
      boundary(2,j) = center(2);
      continue
    end
%
%  Is the point to the left of the box?
%
    if ( x <= x1 )

      boundary(1,j) = x1;

      if ( y <= y1 )
        boundary(2,j) = y1;
      elseif ( y <= y2 )
        boundary(2,j) = y;
      elseif ( y2 <= y )
        boundary(2,j) = y2;
      end
%
%  To the right of the box?
%
    elseif ( x2 <= x )

      boundary(1,j) = x2;

      if ( y <= y1 )
        boundary(2,j) = y1;
      elseif ( y <= y2 )
        boundary(2,j) = y;
      elseif ( y2 <= y )
        boundary(2,j) = y2;
      end
%
%  Below the middle of the box?
%
    elseif ( y <= y1 )

      boundary(1,j) = x;
      boundary(2,j) = y1;
%
%  Above the middle of the box?
%
    elseif ( y2 <= y )

      boundary(1,j) = x;
      boundary(2,j) = y2;
%
%  Inside the box.
%  Figure out which side is closest by drawing the diagonal lines.
%
    else
%
%  Y is small.
%
      if (     y <= x &  y <= -x )
        boundary(1,j) = x;
        boundary(2,j) = y1;
%
%  X is big.
%
      elseif ( y <= x & -y <=  x )
        boundary(1,j) = x2;
        boundary(2,j) = y;
%
%  Y is big.
%
      elseif ( x <= y & -x <=  y )
        boundary(1,j) = x;
        boundary(2,j) = y2;
%
%  X is small.
%
      elseif ( x <= y &  x <= -y )
        boundary(1,j) = x1;
        boundary(2,j) = y;
      end
%
%  For points inside the box, the boundary of the circle might be closer than
%  the boundary of the box.
%
      r = sqrt ( sum ( ( point(1:m,j) - center(1:m)' ).^2 ) );

      dist_circle = abs ( r - r2 );
      dist_box = sqrt ( sum ( ( point(1:m,j) - boundary(1:m,j) ).^2 ) );

      if ( dist_circle <  dist_box )
        boundary(1:m,j) = center(1:m)' + r2 / r * ( point(1:m,j) - center(1:m)' );
      end

    end

  end

  return
end


text:
function vals = polys_vals(polys,ts,tt,r)
idx = 1;
N = length(tt);
vals = zeros(1,N);
for i = 1:N
    t = tt(i);
    if t<ts(idx)
        vals(i) = 0;
    else
        while idx<length(ts) && t>ts(idx+1)+0.0001
            idx = idx+1;
        end
        vals(i) = poly_val(polys(:,idx),t,r);
    end
end
end

author: symao


text:
function plot_vf(vf, M, options)

% plot_vf - plot a vector field with
%   an optional image in the background.
%
% plot_vf(vf, M, options);
%
%   WORKS ONLY FOR 2D VECTOR FIELDS
%
%   set options.display_streamlines=1 to display streamlines.
%
%   See also: plot_tensor_field.
%
%   Copyright (c) 2004 Gabriel Peyre

if nargin<2
    M = [];
end

if iscell(vf)
    nrows = 1;
    if length(vf)>3
        nrows = 2;
    end
    if length(vf)>6
        nrows = 3;
    end
    ncols = ceil(length(vf)/nrows);
    lgd = getoptions(options, 'lgd', []);
    clf;
    for i=1:length(vf)
        if iscell(M)
            Mi = M{i};
        else
            Mi = M;
        end
        subplot(nrows,ncols,i);
        plot_vf(vf{i}, Mi, options);
        if iscell(lgd)
            title(lgd{i});
        else
            title(lgd);
        end
        axis tight;
    end
    return;
end

options.null = 1;
is_oriented = getoptions(options, 'is_oriented', 1);
strech_factor = getoptions(options, 'strech_factor', .6);
reorient = getoptions(options, 'reorient', 0);
linestyle = getoptions(options, 'linestyle', 'b');
display_streamlines = getoptions(options, 'display_streamlines', 0);
streamline_density = getoptions(options, 'streamline_density', 8);
streamline_width = getoptions(options, 'streamline_width', 1);
line_width = getoptions(options, 'line_width', 1);
display_arrows = getoptions(options, 'display_arrows', 1);
subsampling = getoptions(options, 'subsampling', []);
normalize_flow = getoptions(options, 'normalize_flow', 0);

if display_streamlines && ~isfield(options, 'display_arrows')
    display_arrows = 0;
end

if display_arrows==1 && not(isempty(subsampling))
    vf = vf(1:subsampling:end,1:subsampling:end,:);
end

if size(vf,3)~=2
    warning('Dimension >2, cropping ...');
    vf = vf(:,:,1:2);
end

if reorient
    % reorient the vf to x>0
    epsi = sign(vf(:,:,1));
    I = find( epsi==0 );
    epsi(I) = 1;
    vf(:,:,1) = vf(:,:,1).*epsi;
    vf(:,:,2) = vf(:,:,2).*epsi;
end

if normalize_flow
    vf = perform_vf_normalization(vf);
end

n = size(vf,1);
p = size(vf,2);

x = 0:1/(n-1):1;
y = 0:1/(p-1):1;
[Y,X] = meshgrid(y,x);

hold on;

if display_arrows
    imagesc(x,y,M');
    if is_oriented
        h = quiver(X,Y,vf(:,:,1),vf(:,:,2), strech_factor, linestyle);
    else
        h = quiver(X,Y,vf(:,:,1),vf(:,:,2), strech_factor*0.7, linestyle);
        h = quiver(X,Y,-vf(:,:,1),-vf(:,:,2), strech_factor*0.7, linestyle);
    end
    axis xy;
    axis equal;
    set(h, 'LineWidth', line_width);
end

if display_streamlines
    if not(isempty(M))
        imagesc(M');
    end
    [X,Y] = meshgrid(1:n,1:p);
    [XY,tmp] = streamslice(X,Y, vf(:,:,2), vf(:,:,1) ,streamline_density);
    % reverse stream
    for i=1:length(XY)
        XY{i} = XY{i}(:,2:-1:1);
    end
    h = streamline(XY);
    set(h, 'LineWidth', streamline_width);
    axis ij;
end

axis off;
hold off;

author: gpeyre

repo: matlab-toolboxes


text:
function [ssarr, l_arr, cu_arr]=refparams_vecgsm(org,subands,M)

%This function computes the parameters of the reference image. This is
%called by vifvec.m.

for i=1:length(subands);
    sub=subands(i);
    y=org{sub};

    sizey=floor(size(y)./M)*M; % crop  to exact multiple size
    y=y(1:sizey(1),1:sizey(2));

    % Collect MxM blocks. Rearrange each block into an
    % M^2 dimensional vector and collect all such vectors.
    % Collece ALL possible MXM blocks (even those overlapping) from the subband
    temp=[];
    for j=1:M
        for k=1:M
            temp=cat(1,temp,reshape(y(k:end-(M-k), j:end-(M-j)),1,[]));
        end
    end

    % estimate mean and covariance
    mcu=mean(temp')';
    cu=((temp-repmat(mcu,1,size(temp,2)))*(temp-repmat(mcu,1,size(temp,2)))')./size(temp,2); % covariance matrix for U

    % Collect MxM blocks as above. Use ONLY non-overlapping blocks to
    % calculate the S field
    temp=[];
    for j=1:M
        for k=1:M
            temp=cat(1,temp,reshape(y(k:M:end, j:M:end),1,[]));
        end
    end

    % Calculate the S field
    ss=(inv(cu)*temp);
    ss=sum(ss.*temp)./(M*M);
    ss=reshape(ss,sizey/M);

    % Eigen-decomposition
    [v,d]=eig(cu);
    l_arr(sub,:)=diag(d)';

    % rearrange for output
    ssarr{sub}=ss;
    temp=0;
    d=diag(d);
    cu_arr{sub}=cu;
end

repo: SelfExSR


text:
function geometry_test02035 ( )

%*****************************************************************************80
%
%% TEST02035 tests CYLINDER_SAMPLE_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
  n = 20;

  p1 = [ 0.0; -2.0; 0.0 ];
  p2 = [ 0.0;  2.0; 0.0 ];
  r = 1.0;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02035\n' );
  fprintf ( 1, '  CYLINDER_SAMPLE_3D samples points in a cylinder.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Radius R = %f\n', r );
  fprintf ( 1, '  Center of bottom disk = %f  %f  %f\n', p1(1:3,1) );
  fprintf ( 1, '  Center of top disk =    %f  %f  %f\n', p2(1:3,1) );

  [ p, seed ] = cylinder_sample_3d ( p1, p2, r, n, seed );

  r8mat_transpose_print ( 3, n, p, '  Sample points:' );

  return
end


text:
function column_num = file_column_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_COLUMN_COUNT counts the columns in the first line of a file.
%
%  Discussion:
%
%    The file is assumed to be a simple text file.
%
%    Most lines of the file are presumed to consist of COLUMN_NUM words,
%    separated by spaces.  There may also be some blank lines, and some
%    comment lines, which have a "#" in column 1.
%
%    The routine tries to find the first non-comment non-blank line and
%    counts the number of words in that line.
%
%    If all lines are blanks or comments, it goes back and tries to analyze
%    a comment line.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the file.
%
%    Output, integer COLUMN_NUM, the number of columns in the file.
%
  FALSE = 0;
  TRUE = 1;
%
%  Open the file.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_COLUMN_COUNT - Error!' );
    column_num = -1;
    return;
  end
%
%  Read one line, but skip blank lines and comment lines.
%  Use FGETL so we drop the newline character!
%
  got_one = FALSE;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( s_len_trim ( line ) == 0 )

    elseif ( line(1) == '#' )

    else
      got_one = TRUE;
      break;
    end

  end

  fclose ( input_unit );

  if ( got_one == FALSE )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Warning!\n' );
    fprintf ( 1, '  The file does not seem to contain any data.\n' );
    column_num = -1;
    return;
  end

  column_num = s_word_count ( line );

  return
end


text:
function [v,g] = ctransform(x, u, y, p, epsilon)

% ctransform - c-transform and gradient of semi-discrete energy
%
%   [v,g] = ctransform(x,u, y, epsilon);
%
%   u is sampled at points x.
%   v is u^c sampled at points y.
%
%   x must be of size (dim,n)
%   y must be of size (dim,m)
%
% If epsilon==0
%   u^c(y) = min_x  |x-y|^p-u(x)
% else
%   u^c(y) = - epsilon * log( sum_x  exp( (-|x-y|^p + u(x))/epsilon )
%
%   The semi discrete energy between continuous measure a and discrete
%   measure b is
%       E(u) = int u^c(y) da(y) - <u,b>
%   g(:,j) is the gradient of u -> u^c(y_j).
%
%   Copyright (c) 2017 Gabriel Peyre

if nargin<5
    epsilon = 0;
end

n = size(x,2);
m = size(y,2);

c = distmat(x,y).^p; % cost
S = c -  repmat(u(:), [1 m]);
[v,I] = min(S);
if epsilon>0
    % use soft-min, stabilized
    if 1
        S = S - repmat(v, [size(S,1) 1]);
        v = - epsilon * log( sum(  exp( -S/epsilon ) ) ) + v;
    else
        v = - epsilon * log( sum(  exp( -S/epsilon ) ) );
    end
end
v = v(:);

if nargout>1
    % compute also the gradient
    if epsilon==0
        % simply indicator function of the Laguerre cells, indicated by the index I
        g = sparse( I, 1:m, ones(m,1),n, m );
    else
        % smoothed-indicator
        g = exp( -S/epsilon );
        g = g ./ repmat( sum(g), [n 1] );
    end

end

end

author: optimaltransport


text:
classdef CEC2020_F2 < PROBLEM
% <single> <real>
% Shifted and rotated Schwefel's function

%------------------------------- Reference --------------------------------
% C .T. Yue, K. V. Price, P. N. Suganthan, J. J. Liang, M. Z. Ali, B. Y.
% Qu, N. H. Awad, and P. P Biswas, Problem definitions and evaluation
% criteria for the CEC 2020 special session and competition on single
% objective bound constrained numerical optimization, Zhengzhou University,
% China and Nanyang Technological University, Singapore, 2019.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        O;      % Optimal decision vector
        Mat;	% Rotation matrix
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'CEC2020.mat'),'Data');
            obj.O = Data{2}.o;
            obj.M = 1;
            if isempty(obj.D) || obj.D < 10
                obj.D   = 5;
                obj.Mat = Data{2}.M_5;
            elseif obj.D < 15
                obj.D   = 10;
                obj.Mat = Data{2}.M_10;
            elseif obj.D < 20
                obj.D   = 15;
                obj.Mat = Data{2}.M_15;
            else
                obj.D   = 20;
                obj.Mat = Data{2}.M_20;
            end
            obj.lower    = zeros(1,obj.D) - 100;
            obj.upper    = zeros(1,obj.D) + 100;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            Z = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            Y = 10*Z*obj.Mat';
            Z = Y + 4.2097e2;
            g         = Z.*sin(sqrt(abs(Z)));
            temp      = 500 - mod(Z(Z>500),500);
            g(Z>500)  = temp.*sin(sqrt(abs(temp))) - (Z(Z>500)-500).^2/10000/obj.D;
            temp      = mod(abs(Z(Z<-500)),500) - 500;
            g(Z<-500) = temp.*sin(sqrt(abs(temp))) - (Z(Z<-500)-500)/10000/obj.D;
            PopObj = 1100 + 418.9829*obj.D - sum(g,2);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
classdef PrincipalDirectionComputerFactory < handle

    methods (Access = public, Static)

        function pdc = create(cParams)

            switch cParams.type
                case '2D'
                    pdc = PrincipalDirectionComputerIn2D(cParams);
                case '3D'
                    pdc = PrincipalDirectionComputerIn3D(cParams);
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function excelrecord4favarfcn(favar, IRFperiods, endo, IRFt, strctident, pref)
% script excelrecord4
% records the information contained in the worksheet 'IRF' of the excel spreadsheet 'results.xls'

% create the cell that will be saved on excel
irfcell={};

% build preliminary elements: space between the tables
horzspace=repmat({''},2,5*favar.IRF.npltXshck);
vertspace=repmat({''},IRFperiods+3,1);

%subsample
endo2={endo{favar.IRF.plotXshock_index,1}};
informationvariablestrings=favar.informationvariablestrings(1,favar.plotX_index);

% loop over variables (vertical dimension)
for ii=1:favar.npltX
tempcell={};
   % loop over shocks (horizontal dimension)
   for jj=1:favar.IRF.npltXshck %:nfactorshocks, insert again when all dimensions are considered in the irf favar estimation
   % create cell of IRF record for variable ii in response to shock jj
      % if a sign restriction identification scheme has been used, use the structural shock labels
      if IRFt==4
      temp=['response of ' informationvariablestrings{1,ii} ' to ' strctident.signreslabels{jj,1} ' shocks'];
      else
      temp=['response of ' informationvariablestrings{1,ii} ' to ' endo2{1,jj} ' shocks'];
      end
   irf_ij=[temp {''} {''} {''};{''} {''} {''} {''};{''} {'lw. bound'} {'median'} {'up. bound'};num2cell((1:IRFperiods)') num2cell((favar.IRF.favar_irf_estimates{ii,jj})')];
   tempcell=[tempcell irf_ij vertspace];
   end
irfcell=[irfcell;horzspace;tempcell];
end

% trim
irfcell=irfcell(3:end,1:end-1);

% write in excel
if pref.results==1
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),irfcell,'IRF FAVAR','B2');
end

author: european-central-bank

repo: BEAR-toolbox


text:
  function ob = Gembed(varargin)
%|function ob = Gembed([args])
%| Construct Gembed object that embes a short vector into a longer vector.
%| with dimensions [(Nd)].
%|
%| option1:
%|	'list'	[N]		list of array entries
%|	'odim'	ndims
%|
%| option2:
%|	'samp'	logical [(Nd)]	sampling pattern of array entries
%|
%| out
%|	ob	[M N]	fatrix2 object, where M = numel(samp), N = length(list)
%|
%| The two versions are related by list = find(samp), odim = size(samp).
%|
%| Copyright 2010-12-01, Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end
if nargin == 1 && streq(varargin{1}, 'test'), Gembed_test, return, end

arg.samp = [];
arg.list = [];
arg.odim = [];
arg = vararg_pair(arg, varargin);

if ~isempty(arg.samp) && isempty(arg.list) && isempty(arg.odim)
	ob = Gembed_samp_build(arg.samp);

elseif isempty(arg.samp) && ~isempty(arg.list) && ~isempty(arg.odim)
	ob = Gembed_list_build(arg.odim, arg.list);

else
	fail('must give exactly one of samp, or list and odim')
end

% Gembed_list_build()
function ob = Gembed_list_build(odim, list)
arg.odim = odim;
arg.list = list;

ob = fatrix2( ...
	'accept1d', true, ... % trick: OK because of forw(forw()) in mtimes2
	'idim', numel(arg.list), ...
	'odim', arg.odim, ...
	'arg', arg, ...
	'forw', @Gembed_list_forw, ...
	'back', @Gembed_list_back, ...
	'caller', [mfilename '(list)']);
%	'mask', true(numel(arg.list),1), ...

% Gembed_list_forw(): y = A * x
% in
%	x	[N]
% out
%	y	[(M)]
%
function y = Gembed_list_forw(arg, x)

y = zeros([arg.odim 1], class(x)); % [(M)]
y(arg.list) = x;

% Gembed_list_back(): x = A' * y
% in
%	y	[(M)]
% out
%	x	[N]
%
function x = Gembed_list_back(arg, y)

x = y(arg.list);

% Gembed_samp_build()
function ob = Gembed_samp_build(samp)
arg.samp = samp;
arg.odim = size(samp);
if arg.odim(end) == 1
	arg.odim = arg.odim(1:end-1); % remove trailing '1'
end
if ~islogical(arg.samp), error 'samp must be logical', end

ob = fatrix2( ...
	'idim', sum(arg.samp(:)), ...
	'odim', arg.odim, ...
	'arg', arg, ...
	'forw', @Gembed_samp_forw, ...
	'back', @Gembed_samp_back, ...
	'caller', [mfilename '(list)']);

% Gembed_samp_forw()
function y = Gembed_samp_forw(arg, x)
y = zeros([arg.odim 1], class(x)); % [(M)]
y(arg.samp) = x;

% Gembed_samp_back()
function x = Gembed_samp_back(arg, y)
x = y(arg.samp);

% Gembed_test
function Gembed_test

odim = [5 2];
list = [2 1 9 7];
odim = prod(odim); % required because of fatrix2 full 1d ambiguity

for ii=1:2
	if ii==1 % type list
		A = Gembed('odim', odim, 'list', list);
	else % type samp
		samp = false([odim 1]);
		samp(list) = 1;
		A = Gembed('samp', samp);
	end

	fatrix2_tests(A)
	test_adjoint(A);
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
classdef CaffePooling < dagnn.Filter
%CAFFEPOOLING - This layer implements Caffe-style pooling, which applies
% the pooling kernel to the padding, as well as to the input (unlike matconvnet
% which does not include the padding itself in the kernel computation).
% To mimic this effect, we "pre-pad" the input using an empty convolution, then
% pass to the standard pooling layer.

  properties
    method = 'max'
    poolSize = [1 1]
    opts = {'cuDNN'}
  end

  methods
    function outputs = forward(self, inputs, params)
      outputs{1} = vl_nncaffepool(inputs{1}, self.poolSize, ...
                                  'pad', self.pad, ...
                                  'method', self.method, ...
                                  'stride', self.stride, ...
                                  'extraArgs', self.opts) ;
    end

    function [derInputs, derParams] = backward(self, inputs, params, derOutputs)
      derInputs{1} = vl_nncaffepool(inputs{1}, self.poolSize, derOutputs{1}, ...
                                  'pad', self.pad, ...
                                  'method', self.method, ...
                                  'stride', self.stride, ...
                                  'extraArgs', self.opts) ;
      derParams = {} ;
    end

    function kernelSize = getKernelSize(obj)
      kernelSize = obj.poolSize ;
    end

    function outputSizes = getOutputSizes(obj, inputSizes)
      outputSizes = getOutputSizes@dagnn.Filter(obj, inputSizes) ;
      outputSizes{1}(3) = inputSizes{1}(3) ;
    end

    function obj = CaffePooling(varargin)
      obj.load(varargin) ;
    end
  end
end

repo: MFT


text:
function E = bryant2angular(rpy)
% Jacobian of Omega in terms of rolldot, pitchdot and yawdot
% "Bryant Angle Representation"
% http://www.u.arizona.edu/~pen/ame553/Hallway/CAAMS/CAAMS_15.pdf

roll = rpy(1); %x-axis
pitch = rpy(2); %y-axis
yaw = rpy(3); %z-axis

% paul's version
% E = [ cos(yaw)*cos(pitch),           -sin(yaw),                   0;
%       sin(yaw)*cos(pitch),            cos(yaw),                   0;
%               -sin(pitch),                   0,                   1];

% bryant angle version
E = [   cos(roll)*cos(yaw),  sin(yaw),  0;
      -cos(pitch)*sin(yaw),  cos(yaw),  0;
                sin(pitch),         0,  1];

end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
% [x, Covx, Covx_x, entropy] = rts_smoother_step(x, Covx, x_s, Covx_s, A, Q)

% Last modified 2011-10-19
% Copyright (c) Jaakko Luttinen (jaakko.luttinen@aalto.fi)

function [x, Covx, Covx_x, entropy] = rts_smoother_step(x, Covx, x_s, Covx_s, A, Q)
% Perform the RTS smoothing step

x_p = A*x;
Covx_p = A*Covx*A' + Q;

S = (Covx*A') / Covx_p;
x = x + S*(x_s-x_p);
if nargout >= 2
  Covx = Covx + S*(Covx_s-Covx_p)*S';
end
if nargout >= 3
  Covx_x = S*Covx_s;
end
if nargout >= 4
  % Compute entropy term:
  %
  % INT[ p(x_n, x_(n+1) | Y) log p(x_n | Y, x_(n+1)) dx_n dx_(n+1) ]
  Cov_joint = [Covx    Covx_x
               Covx_x' Covx_s];
  entropy = gaussian_entropy(logdet_cov(Cov_joint), size(Cov_joint,1)) ...
            - gaussian_entropy(logdet_cov(Covx_s), size(Covx_s,1));

end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function g = p41_g ( n, x )

%*****************************************************************************80
%
%% P41_G evaluates the gradient for problem 41.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the values of the variables.
%
%    Output, real G(N), the gradient of the objective function.
%
  g = zeros ( 3, 1 );

  g(1) = 400.0 * x(1)^3 - 400.0 * x(2) * x(1) ...
    + 2.0 * x(1) - 2.0;

  g(2) = -200.0 * x(1)^2 + 220.2 * x(2) + 19.8 * x(4) - 40.0;

  g(3) = -360.0 * x(3) * x(4) + 360.0 * x(3)^3 ...
    + 2.0 * x(3) - 2.0;

  g(4) = + 180.0 * x(4) - 180.0 * x(3)^2 + 20.2 * x(4) ...
    + 19.8 * x(2) - 40.0;

  return
end


text:
function output = createOutput
    output.iterInit    = 0;
    output.cost        = 0;
    output.KKTError       = 0;
    output.timeElapsed.searchDirection = 0;
    output.timeElapsed.lineSearch = 0;
    output.timeElapsed.KKTErrorCheck = 0;
    output.timeElapsed.total = 0;
    output.rho         = 0;
    output.iterTotal     = 0;
    output.exitflag    = 0;
end

author: deng-haoyang

repo: ParNMPC


text:
function angles = wrapMinusPiPifast(angles)

tol = 0.001;

angles(angles < -pi) = angles(angles < -pi) + 2*pi;
angles(angles < -pi) = angles(angles < -pi) + 2*pi;

angles(angles > pi) = angles(angles > pi) - 2*pi;
angles(angles > pi) = angles(angles > pi) - 2*pi;

angles(angles == pi) = angles(angles == pi) - tol;
angles(angles == -pi) = angles(angles == -pi) + tol;

assert(~any(angles > pi));
assert(~any(angles < -pi));

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut