text:
% psf_mismatch_example2.m
%
% 2D example showing the effect of PSF mismatch on ML-EM algorithm
%
% Copyright 2001-8-24, Jeff Fessler, The University of Michigan

if ~has_aspire, return, end

%
% generate data
%
if ~isvar('yi'), printm 'setup psf_mismatch_example2'
	f.dir = test_dir;
	f.dsc0 = [f.dir 't0.dsc'];
	f.dsc1 = [f.dir 't1.dsc'];
	f.dsc2 = [f.dir 't2.dsc'];
	f.wtf0 = strrep(f.dsc0, 'dsc', 'wtf');
	f.wtf1 = strrep(f.dsc1, 'dsc', 'wtf');
	f.wtf2 = strrep(f.dsc2, 'dsc', 'wtf');
	if 1
		Fwhm0 = 5;
		os_run(['wt -chat 0 dsc 12 fwhm_detector 5 >! ' f.dsc0])
		os_run(['echo y | wt -chat 0 gen ' f.dsc0])
		Fwhm1 = 7;
		os_run(['wt -chat 0 dsc 12 fwhm_detector 7 >! ' f.dsc1])
		os_run(['echo y | wt -chat 0 gen ' f.dsc1])
		Fwhm2 = 2;
		os_run(['wt -chat 0 dsc 12 fwhm_detector 2 >! ' f.dsc2])
		os_run(['echo y | wt -chat 0 gen ' f.dsc2])
	end

%	G0 = Gtomo2_wtmex('f.wtf0'); % cannot use with multiple .wtf's

	[t nx ny nb na] = wtf_read(f.wtf0);
	ig = image_geom('nx', nx, 'ny', ny, 'dx', 1);
	ig.mask = reshape(sum(t) ~= 0, nx, ny);

	G0 = Gsparse(f.wtf0, 'mask', ig.mask);
	G1 = Gsparse(f.wtf1, 'mask', ig.mask);
	G2 = Gsparse(f.wtf2, 'mask', ig.mask);

	sg = sino_geom('par', 'nb', nb, 'na', na, 'dr', 1);
	xtrue = ellipse_im(ig, [6.5 -0.5 5 5 0 100], 'oversample', 4);

	ri = 1;
	yi = sg.shape(G0 * xtrue(ig.mask)) + ri;
	im(yi, 'yi'), cbar
prompt
end

if ~isvar('Gb2'), printm 'Gbs'
	f.nblock = 6;
	Gb0 = Gblock(G0, f.nblock);
	Gb1 = Gblock(G1, f.nblock);
	Gb2 = Gblock(G2, f.nblock);
prompt
end

% uniform initial image
xinit = ones(sum(ig.mask(:)),1);

% FBP
if 0 && ~isvar('xfbp'), printm 'fbp'
	xfbp = em_fbp(sg, ig, yi, 1, ri);
	im(xfbp), cbar
prompt
end

if ~isvar('x2') && 0
	x0 = eml_em(xinit, G0, yi(:), 1, ri, 'isave', 'all', 'niter', f.niter);
	x1 = eml_em(xinit, G1, yi(:), 1, ri, 'isave', 'all', 'niter', f.niter);
	x2 = eml_em(xinit, G2, yi(:), 1, ri, 'isave', 'all', 'niter', f.niter);
	x0 = ig.embed(x0);
	x1 = ig.embed(x1);
	x2 = ig.embed(x2);
end

%
% OS-EM iterations
%
if ~isvar('xo2'), printm 'run os-em'
	f.niter = 40;
	rri = ri * ones(nb,na);
	osem = @(G) eml_osem(xinit, G, ...
		yi, [], rri, 'niter', f.niter, 'isave', 'all');
	xo0 = osem(Gb0);
	xo1 = osem(Gb1);
	xo2 = osem(Gb2);
	xo0 = ig.embed(xo0);
	xo1 = ig.embed(xo1);
	xo2 = ig.embed(xo2);
	im(xo2)
prompt
end

if ~isvar('fw2'), printm 'find widths'
	ix = 39+[-9:9];
	iy = 33+[-9:9];

	xx0 = xo0(ix,iy,2:end);
	xx1 = xo1(ix,iy,2:end);
	xx2 = xo2(ix,iy,2:end);
	im(xx0)
	clear fw0 fw1 fw2
	for ii=1:f.niter
		fw0(ii) = fwhm2(xx0(:,:,ii));
		fw1(ii) = fwhm2(xx1(:,:,ii));
		fw2(ii) = fwhm2(xx2(:,:,ii));
	end

	xx = xtrue(ix,iy);
	fx = fwhm2(xx);
end

if 1 && im
	im clf
	ii = 1:f.niter;
	plot(ii(1:4:end), fw0(1:4:end), 'yo', ...
		ii(1:4:end), fw1(1:4:end), 'cx', ...
		ii(1:4:end), fw2(1:4:end), 'g+', ...
		ii, fx + 0*ii, 'b-', ...
		ii, fw0, 'y-', ...
		ii, fw1, 'c-', ...
		ii, fw2, 'g-')
	axisy(3,11)
	legend(	sprintf('True system PSF, FWHM=%g', Fwhm0), ...
		sprintf('Model PSF too big, FWHM=%g', Fwhm1), ...
		sprintf('Model PSF too small, FWHM=%g', Fwhm2), ...
		'ideal')
	xlabel 'Iteration of OS-EM algorithm'
	ylabel 'Width of reconstructed circle'
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function pde = elli3DorthocircIntf(bm,bp,rx,ry,rz)
%% USAGE: polynomial solution for Poisson equation
%  Last Modified: 02/21/2020 by Xu Zhang

%% PDE Structure
pde = struct('intf',@intf,'f',@f,'fm',@fm,'fp',@fp,'exactu',@exactu,...
    'um',@um,'up',@up,'Dxu',@Dxu,'Dxum',@Dxum,'Dxup',@Dxup,'Dyu',@Dyu,...
    'Dyum',@Dyum,'Dyup',@Dyup,'Dzu',@Dzu,'Dzum',@Dzum,'Dzup',@Dzup,...
    'A',@A,'Am',@Am,'Ap',@Ap,'one',@one,'gD',@gD);

pde.bm = bm;
pde.bp = bp;
%% interface function
    function u = intf(x,y,z)
        u = log(F1(x,y,z).*F2(x,y,z).*F3(x,y,z)-r(x,y,z)+1);
    end

%% exact solution
    function u = exactu(x,y,z)
        u = um(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = up(x(id),y(id),z(id));
    end
    function u = um(x,y,z)
        u = intf(x,y,z)/bm;
    end
    function u = up(x,y,z)
        u = intf(x,y,z)/bp;
    end
%% Boundary Function
    function u = gD(x,y,z)
        u = exactu(x,y,z);
    end
%% Derivative of the exact solution
    function u = Dxu(x,y,z)
        u = Dxum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Dxup(x(id),y(id),z(id));
    end
    function u = Dxum(x,y,z)
        fx = (F1x(x,y,z).*F2(x,y,z).*F3(x,y,z)+F1(x,y,z).*F2x(x,y,z).*F3(x,y,z)+...
            F1(x,y,z).*F2(x,y,z).*F3x(x,y,z)) - Dxr(x,y,z);
        u = (F1(x,y,z).*F2(x,y,z).*F3(x,y,z)-r(x,y,z)+1).^(-1).*fx/bm;
    end
    function u = Dxup(x,y,z)
        fx = (F1x(x,y,z).*F2(x,y,z).*F3(x,y,z)+F1(x,y,z).*F2x(x,y,z).*F3(x,y,z)+...
            F1(x,y,z).*F2(x,y,z).*F3x(x,y,z)) - Dxr(x,y,z);
        u = (F1(x,y,z).*F2(x,y,z).*F3(x,y,z)-r(x,y,z)+1).^(-1).*fx/bp;
    end

    function u = Dyu(x,y,z)
        u = Dyum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Dyup(x(id),y(id),z(id));
    end
    function u = Dyum(x,y,z)
        fy = (F1y(x,y,z).*F2(x,y,z).*F3(x,y,z)+F1(x,y,z).*F2y(x,y,z).*F3(x,y,z)+...
            F1(x,y,z).*F2(x,y,z).*F3y(x,y,z)) - Dyr(x,y,z);
        u = (F1(x,y,z).*F2(x,y,z).*F3(x,y,z)-r(x,y,z)+1).^(-1).*fy/bm;
    end
    function u = Dyup(x,y,z)
        fy = (F1y(x,y,z).*F2(x,y,z).*F3(x,y,z)+F1(x,y,z).*F2y(x,y,z).*F3(x,y,z)+...
            F1(x,y,z).*F2(x,y,z).*F3y(x,y,z)) - Dyr(x,y,z);
        u = (F1(x,y,z).*F2(x,y,z).*F3(x,y,z)-r(x,y,z)+1).^(-1).*fy/bp;
    end

    function u = Dzu(x,y,z)
        u = Dzum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Dzup(x(id),y(id),z(id));
    end
    function u = Dzum(x,y,z)
        fz = (F1z(x,y,z).*F2(x,y,z).*F3(x,y,z)+F1(x,y,z).*F2z(x,y,z).*F3(x,y,z)+...
            F1(x,y,z).*F2(x,y,z).*F3z(x,y,z)) - Dzr(x,y,z);
        u = (F1(x,y,z).*F2(x,y,z).*F3(x,y,z)-r(x,y,z)+1).^(-1).*fz/bm;
    end
    function u = Dzup(x,y,z)
        fz = (F1z(x,y,z).*F2(x,y,z).*F3(x,y,z)+F1(x,y,z).*F2z(x,y,z).*F3(x,y,z)+...
            F1(x,y,z).*F2(x,y,z).*F3z(x,y,z)) - Dzr(x,y,z);
        u = (F1(x,y,z).*F2(x,y,z).*F3(x,y,z)-r(x,y,z)+1).^(-1).*fz/bp;
    end

%% right hand side function
    function u = f(x,y,z)
        u = fm(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = fp(x(id),y(id),z(id));
    end
    function u = fm(x,y,z)
        r_xx = ry^2*(rz*2);
        r_yy = ry^2*(rz*2);
        r_zz = ry^2*(rz*2);

        fx = (F1x(x,y,z).*F2(x,y,z).*F3(x,y,z)+F1(x,y,z).*F2x(x,y,z).*F3(x,y,z)+...
            F1(x,y,z).*F2(x,y,z).*F3x(x,y,z)) - Dxr(x,y,z);
        fy = (F1y(x,y,z).*F2(x,y,z).*F3(x,y,z)+F1(x,y,z).*F2y(x,y,z).*F3(x,y,z)+...
            F1(x,y,z).*F2(x,y,z).*F3y(x,y,z)) - Dyr(x,y,z);
        fz = (F1z(x,y,z).*F2(x,y,z).*F3(x,y,z)+F1(x,y,z).*F2z(x,y,z).*F3(x,y,z)+...
            F1(x,y,z).*F2(x,y,z).*F3z(x,y,z)) - Dzr(x,y,z);

        fxx = F1xx(x,y,z).*F2(x,y,z).*F3(x,y,z) + F1x(x,y,z).*F2x(x,y,z).*F3(x,y,z) + F1x(x,y,z).*F2(x,y,z).*F3x(x,y,z) +...
            F1(x,y,z).*F2xx(x,y,z).*F3(x,y,z) + F1x(x,y,z).*F2x(x,y,z).*F3(x,y,z) + F1(x,y,z).*F2x(x,y,z).*F3x(x,y,z) +...
            F1(x,y,z).*F2(x,y,z).*F3xx(x,y,z) + F1(x,y,z).*F2x(x,y,z).*F3x(x,y,z) + F1x(x,y,z).*F2(x,y,z).*F3x(x,y,z) - r_xx;

        fyy = F1yy(x,y,z).*F2(x,y,z).*F3(x,y,z) + F1y(x,y,z).*F2y(x,y,z).*F3(x,y,z) + F1y(x,y,z).*F2(x,y,z).*F3y(x,y,z) +...
            F1(x,y,z).*F2yy(x,y,z).*F3(x,y,z) + F1y(x,y,z).*F2y(x,y,z).*F3(x,y,z) + F1(x,y,z).*F2y(x,y,z).*F3y(x,y,z) +...
            F1(x,y,z).*F2(x,y,z).*F3yy(x,y,z) + F1(x,y,z).*F2y(x,y,z).*F3y(x,y,z) + F1y(x,y,z).*F2(x,y,z).*F3y(x,y,z) - r_yy;

        fzz = F1zz(x,y,z).*F2(x,y,z).*F3(x,y,z) + F1z(x,y,z).*F2z(x,y,z).*F3(x,y,z) + F1z(x,y,z).*F2(x,y,z).*F3z(x,y,z) +...
            F1(x,y,z).*F2zz(x,y,z).*F3(x,y,z) + F1z(x,y,z).*F2z(x,y,z).*F3(x,y,z) + F1(x,y,z).*F2z(x,y,z).*F3z(x,y,z) +...
            F1(x,y,z).*F2(x,y,z).*F3zz(x,y,z) + F1(x,y,z).*F2z(x,y,z).*F3z(x,y,z) + F1z(x,y,z).*F2(x,y,z).*F3z(x,y,z) - r_zz;

        u = -(F1(x,y,z).*F2(x,y,z).*F3(x,y,z)-r(x,y,z)+1).^(-1).*(fxx+fyy+fzz) +...
            (F1(x,y,z).*F2(x,y,z).*F3(x,y,z)-r(x,y,z)+1).^(-2).*(fx.^2+fy.^2+fz.^2);
    end
    function u = fp(x,y,z)
        u = fm(x,y,z);
    end

%% Diffusion coefficient function
    function u = A(x,y,z)
        u = Am(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Ap(x(id),y(id),z(id));
    end
    function u = Am(x,y,z)
        u = bm*ones(size(x));
    end
    function u = Ap(x,y,z)
        u = bp*ones(size(x));
    end

%% Other function
    function u = one(x,y,z)
        u = ones(size(x));
    end
    function u = r(x,y,z)
        u = ry^2*(1+rz*(x.^2 + y.^2 + z.^2));
    end
    function u = Dxr(x,y,z)
        u = ry^2*(rz*2*x);
    end
    function u = Dyr(x,y,z)
        u = ry^2*(rz*2*y);
    end
    function u = Dzr(x,y,z)
        u = ry^2*(rz*2*z);
    end
    function u = F1(x,y,z)
        u = (x.^2 + y.^2 -rx^2).^2 + z.^2;
    end
    function u = F2(x,y,z)
        u = (x.^2 + z.^2 -rx^2).^2 + y.^2;
    end
    function u = F3(x,y,z)
        u = (y.^2 + z.^2 -rx^2).^2 + x.^2;
    end
    function u = F1x(x,y,z)
        u = 2*(x.^2 + y.^2 -rx^2).*(2*x);
    end
    function u = F1xx(x,y,z)
        u = 2*(x.^2 + y.^2 -rx^2).*(2) + 8*x.^2;
    end
    function u = F2x(x,y,z)
        u = 2*(x.^2 + z.^2 -rx^2).*(2*x);
    end
    function u = F2xx(x,y,z)
        u = 2*(x.^2 + z.^2 -rx^2).*(2) + 8*x.^2;
    end
    function u = F3x(x,y,z)
        u = 2*x;
    end
    function u = F3xx(x,y,z)
        u = 2*ones(size(x));
    end
    function u = F1y(x,y,z)
        u = 2*(x.^2 + y.^2 -rx^2).*(2*y);
    end
    function u = F1yy(x,y,z)
        u = 2*(x.^2 + y.^2 -rx^2).*(2) + 8*y.^2;
    end
    function u = F2y(x,y,z)
        u = 2*y;
    end
    function u = F2yy(x,y,z)
        u = 2*ones(size(y));
    end
    function u = F3y(x,y,z)
        u = 2*(y.^2 + z.^2 -rx^2).*(2*y);
    end
    function u = F3yy(x,y,z)
        u = 2*(y.^2 + z.^2 -rx^2)*(2) + 8*y.^2;
    end
    function u = F1z(x,y,z)
        u = 2*z;
    end
    function u = F1zz(x,y,z)
        u = 2*ones(size(z));
    end
    function u = F2z(x,y,z)
        u = 2*(x.^2 + z.^2 -rx^2).*(2*z);
    end
    function u = F2zz(x,y,z)
        u = 2*(x.^2 + z.^2 -rx^2).*(2) + 8*z.^2;
    end
    function u = F3z(x,y,z)
        u = 2*(y.^2 + z.^2 -rx^2).*(2*z);
    end
    function u = F3zz(x,y,z)
        u = 2*(y.^2 + z.^2 -rx^2).*(2) + 8*z.^2;
    end
end

repo: ifem


text:
function X=brutepersp(X,t,y)

X.basis = [zeros(size(X.basis,1),1) X.basis];
X.lmi_variables = [t y];
[i,j] = sort(X.lmi_variables);
X.basis = X.basis(:,[1 1+j]);
X.lmi_variables = i;
X.conicinfo = [0 0];

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function ef(varargin)
% Command-style shortcut for exp_fullform(varargin)
%
% See also:
%   exp_fullform
%
% Example:
% >> ef var1 var2
%    var1 =
%      ...
%    var2 =
%      ...
%

% Copyright (C) Christian Kothe, SCCN, 2011, christian@sccn.ucsd.edu
%
% 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 2 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,
% write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
% USA

if ~iscellstr(varargin)
    error('ef is to be used as a command, i.e. without brackets.'); end

for i=1:length(varargin)
    disp([varargin{i} ' = ']);
    disp(['  ' exp_fullform(evalin('caller',varargin{i}))]);
    disp('');
end

repo: REST


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%    gb 05/05/05
%
% This script has to be executed directly in the command line
% It plots the MSE graphs for all the sujects that are in the repository
% folder : /Snarp/u1/data/reading_longitude/fmri
%
% After a graph is displayed, type Enter to plot the next one.
%
% Type Ctrl + c in the command line to terminate the program
%
% Just make sure the VISTASOFT directory is in the path and type :
% scriptReadMSE
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all
clear all

if isunix
    networkPath = '/snarp1/u1';
else
    networkPath = '\\snarp\u1\';
end

rootDir = fullfile(networkPath,'data','reading_longitude','fmri');

cd(rootDir);
[fileNum dirName] = countDirs(pwd);

for count = 3:fileNum
    try
        currentDir = dirName{count};

        cd(fullfile(rootDir,currentDir));
        motionCompPlotMSE(currentDir);
        pause;

    end

    close all;

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function x = subsasgn( x, S, y )

%   Disciplined convex/geometric programming information for SUBSASGN:
%      Subscripting can be used to change values of elements or slices
%      of CVX variables in the same manner as with numeric arrays. All
%      conventions are preserved, including the colon ':' and 'end'
%      notation, as well as the ability to "expand" a CVX variable by
%      assigning a value to a location outside of its current dimensions
%      (e.g., X(end+1)=0).
%
%      One notable exception is this: if the right-hand side is a CVX
%      expression, then the left-hand side must be as well. So, for
%      example, the following will fail:
%         variable x;
%         y = ones(3,1);
%         y(2) = x;
%      This is because MATLAB does not know how to automatically promote
%      'y' to a CVX variable so that it can accept 'x' as an element. If
%     you want
%      to accomplish something like this, you must manually convert y into a
%      CVX variable first, as follows:
%         variable x;
%         y = cvx(ones(3,1));
%         y(2) = x;

error( nargchk( 3, 3, nargin ) );

%
% Test subscripts
%

szx = size( x );
szy = size( y );
nlx = prod( szx );
try
    temp = reshape( 1 : nlx, szx );
    ndx_x = builtin( 'subsasgn', temp, S, zeros( szy ) );
catch
    error( lasterr );
end
szx_n = size( ndx_x );

%
% Assign data
%

x = cvx( x );
bx = x.basis_;
if any( szx_n < szx ),
    bx = bx( :, ndx_x );
else
    if any( szx_n > szx ),
        bx( :, end + 1 ) = 0;
        ndx_x( ndx_x == 0 ) = size( bx, 2 );
        bx = bx( :, ndx_x );
        temp = reshape( 1 : prod( szx_n ), szx_n );
    end
    ndx_x = builtin( 'subsref', temp, S );
    ndx_x = ndx_x( : );
    nlz = length( ndx_x );
    y = cvx( y );
    by = y.basis_;
    nx = size( bx, 1 );
    [ ny, my ] = size( by );
    if nx < ny,
        if issparse( by ) && ~issparse( bx ), bx = sparse( bx ); end
        bx( ny, : ) = 0;
    elseif nx > ny,
        if issparse( bx ) && ~issparse( by ), by = sparse( by ); end
        by( nx, : ) = 0;
    end
    if my < nlz,
        by = by( :, ones( 1, nlz ) );
    end
    bx( :, ndx_x ) = by;
end

%
% Create the new object
%

x = cvx( szx_n, bx );

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.


text:
function rf=lpcla2rf(la)
%LPCLA2RF Convert log areas to reflection coefficients RF=(LA)

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: lpcla2rf.m,v 1.4 2007/05/04 07:01:38 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 2 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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nf,p2]=size(la);
rf=-tanh((la(:,1:p2-1)-la(:,2:p2))/2);

author: decouples


text:
function varargout = prtRvUtilPlotMvnEllipse(mu,Sigma,nStds,n)
% prtRvUtilPlotMvnEllipse plots the confidence region of a multivariate Gaussian
%
% Syntax: [ellipseHandle, contourXY] = prtRvUtilPlotMvnEllipse(mu,Sigma,nStds)
%
% Adapted From: http://www.mathworks.com/matlabcentral/fileexchange/8793
% Then re-adpated from Beals Mixture of Factor Analyzes Code

%assert(ismember(numel(mu),[2 3]),'prtRvUtilPlotMvnEllipse only works in 2 or 3 /lr   Bv  5,plot(sigmaEigVectors(1,iVec)*[-1 1]*sigmaEigValues(iVec,iVec)*nStd + cMean(1),sigmaEigVectors(2,iVec)*[-1 1]*sigmaEigValues(iVec,iVec)*nStd + cMean(2),'color',cmap(max(ceil(cPi*size(cmap,1)),1),:));dimensions.');

if nargin < 3 || isempty(nStds)
    nStds = 1;
end
if nargin < 4 || isempty(n)
    n = 20;
end
mu = mu(:);

switch numel(mu)
    case 2
        theta = (0:1:n-1)/(n-1)*2*pi;

        epoints = nStds*sqrtm(Sigma) * [cos(theta); sin(theta)]*1   + mu*ones(1,n);

        ellipseHandle = plot(epoints(1,:),epoints(2,:),'LineWidth',3);

    case 3
        theta = (0:1:n-1)'/(n-1)*pi;
        phi = (0:1:n-1)/(n-1)*2*pi;

        sx = sin(theta)*cos(phi);
        sy = sin(theta)*sin(phi);
        sz = cos(theta)*ones(1,n);

        svect = [reshape(sx,1,n*n); reshape(sy,1,n*n); reshape(sz,1,n*n)];
        epoints = nStds*sqrtm(Sigma) * svect  + mu*ones(1,n*n);

        ex = reshape(epoints(1,:),n,n);
        ey = reshape(epoints(2,:),n,n);
        ez = reshape(epoints(3,:),n,n);

        ellipseHandle = mesh(ex,ey,ez);
end
% [Ev,D] = eig(Sigma);
% iV = inv(Sigma);
%
% % Find the l    arger projection
% P = [1,0;0,0];  % X-axis projection operator
% P1 = P * 2*sqrt(D(1,1)) * Ev(:,1);
% P2 = P * 2*sqrt(D(2,2)) * Ev(:,2);
% if abs(P1(1)) >= abs(P2(1)),
%     Plen = P1(1);
% else
%     Plen = P2(1);
% end
% count = 1;
% step = 0.001*Plen;
% Contour1 = zeros(2001,2);
% Contour2 = zeros(2001,2);
% for x = -Plen:step:Plen,
%     a = iV(2,2);
%     b = x * (iV(1,2)+iV(2,1));
%     c = (x^2) * iV(1,1) - 1;
%     Root1 = (-b + sqrt(b^2 - 4*a*c))/(2*a);
%     Root2 = (-b - sqrt(b^2 - 4*a*c))/(2*a);
%     if isreal(Root1),
%         Contour1(count,:) = [x,Root1] + mu(:)';
%         Contour2(count,:) = [x,Root2] + mu(:)';
%         count = count + 1;
%     end
% end
% Contour1 = Contour1(1:count-1,:);
% Contour2 = [Contour1(1,:);Contour2(1:count-1,:);Contour1(count-1,:)];
%
% contourXY = cat(1,Contour1,flipud(Contour2));
%
% ellipseHandle = plot(contourXY(:,1),contourXY(:,2));

varargout = {};
if nargout
    varargout = {ellipseHandle};
end

end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
%

% problem 7 - convolution of x[n] and h[n]

n=0:10;
u=ones(size(n));
n1=0:3;
u4_1=zeros(size(n1));
n2=4:10;
u4_2=ones(size(n2));
u4=[u4_1 u4_2];
x=u-u4;
h=0.7.^n;
y=conv(x,h);
stem(0:20,y);
xlim([-1 21]);
legend('y[n]')

author: Sable

repo: mcbench-benchmarks


text:
%--------------------------------------------------------------------------
%   某个AD测试复制过来的ADC性能分析代码，还没研究注释其测试效果
%   未完成->打算与函数ad_analyzer整合
%--------------------------------------------------------------------------
function [ENOB, SINAD, SFDR, SNR] = LinsiFFT(Dout,fclk)

%   Recenter the digital sine wave
    %Dout=Dout-(2^numbit-1)/2;
    Dout=Dout-mean(Dout);
    numpt=length(Dout);
    har_N=9;
 % Plot results in the time domain

    %If no window function is used, the input tone must be chosen to be unique and with
    %regard to the sampling frequency. To achieve this prime numbers are introduced and the
    %input tone is determined by fIN = fSAMPLE * (Prime Number / Data Record Size).
    %To relax this requirement, window functions such as HANNING and HAMING (see below) can
    %be introduced, however the fundamental in the resulting FFT spectrum appears 'sharper'
    %without the use of window functions.

    Doutw=Dout.*blackmanharris(numpt);

    %Performing the Fast Fourier Transform
    Dout_spect=fft(Doutw);
    Dout_spect(1:3)=0;
   % Dout_spect(numpt/2-4:numpt/2)=0;
    %Recalculate to dB
    Dout_dB=20*log10(abs(Dout_spect));
    maxdB=max(Dout_dB(1:numpt/2));
  %Span of the input frequency on each side
    span=max(round(numpt/200),5);    % default (numpt/200,5)
    span=5;
    %Calculate SNR, SINAD, THD and SFDR values
    %Find the signal bin number, DC = bin 1
    fin=find(Dout_dB(1:numpt/2)==maxdB);
    Dout_dB_tmp=[Dout_dB(1:fin-span); Dout_dB(fin+span:numpt/2)];
    Sec_maxdB=max(Dout_dB_tmp(1:end));

    %Determine power spectrum
    spectP=(abs(Dout_spect)).*(abs(Dout_spect));

    %Find DC offset power
    Pdc=sum(spectP(1:span));
    %Extract overall signal power
    Ps=sum(spectP(fin-span:fin+span));
    %Vector/matrix to store both frequency and power of signal and harmonics
    Fh=[];
    %The 1st element in the vector/matrix represents the signal, the next element represents
    %the 2nd harmonic, etc.
    Ph=[];
  %Approximate search span for harmonics on each side
    spanh=5;
    %Find harmonic frequencies and power components in the FFT spectrum
    for har_num=1:har_N
    %Input tones greater than fSAMPLE are aliased back into the spectrum
    tone=rem((har_num*(fin-1)+1)/numpt,1);
    if tone>0.5
    %Input tones greater than 0.5*fSAMPLE (after aliasing) are reflected
    tone=1-tone;
    end
    Fh=[Fh tone];
    %For this procedure to work, ensure the folded back high order harmonics do not overlap
    %with DC or signal or lower order harmonics
     index_tmp=round(tone*numpt)-spanh;
    if index_tmp<1
        index_tmp=1;
    end
    har_peak=max(spectP(index_tmp:round(tone*numpt)+spanh));
    har_bin=find(spectP(index_tmp:round(tone*numpt)+spanh)==har_peak);
    har_bin=har_bin+round(tone*numpt)-spanh-1;
    if har_bin<2
        har_bin=2;
    end
    Ph=[Ph sum(spectP(har_bin-1:har_bin+1))];
    end
    %Determine the total distortion power
    Pd=sum(Ph(2:har_N));
    %Determine the noise power
    Pn=sum(spectP(1:numpt/2))-Ps-Pd;

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    format;
    A=(max(Dout)-min(Dout))/2;
    SINAD=10*log10(Ps/(Pn+Pd));
    ENOB=(SINAD-1.76)/6.02;
    SNR=10*log10(Ps/Pn);
    %THD is calculated from 2nd through 5th order harmonics
    THD=10*log10(Pd/Ph(1));
    SFDR=maxdB-Sec_maxdB;
    %Signal & Harmonic Power Components
    HD=10*log10(Ph(1:har_N)/Ph(1));

    %Display the results in the frequency domain with an FFT plot

    fclk=fclk/1e6; % change Mhz to Khz
    FFin=fin/numpt*fclk;
    fin_dis=round(FFin*10)/10;
    SINAD_dis=round(SINAD*100)/100;
    ENOB_dis=round(ENOB*100)/100;
    SFDR_dis=round(SFDR*100)/100;
    SNR_dis=round(SNR*100)/100;
    %For TTIMD, use the following short routine, normalized to 锟?.5dB full-scale.

    set (gcf,'Position',[10,10,650,400]);
    set(gca,'fontsize', 14);
    %plot([1:numpt/2].*fclk/numpt,Dout_dB(1:numpt/2)-maxdB);
    stem([0:numpt/2-1].*fclk/numpt,Dout_dB(1:numpt/2)-maxdB,'LineWidth',1.5,'BaseValue',-120);
    set(gca,'FontSize',14,'FontWeight','bold');
    grid on;
    %title('FFT PLOT');
    xlabel('Frequency [MHz]', 'fontsize', 16,'fontweight', 'bold');
    ylabel('[dB]', 'fontsize', 16,'fontweight', 'bold');
    a1=axis;
    axis([a1(1) fclk/2 -120 a1(4)]);
    %text(10,-10,['SINAD=',num2str(SINAD) ' ','ENOB=' num2str(ENOB) ' ','SFDR=' num2str(SFDR)  '    ','FFT--' num2str(numpt)], 'fontsize', 12, 'fontweight', 'bold');
    %text(fclk/4-14*fclk/100,-6,[num2str(numbit) '-bit' '  mode' ], 'fontsize', 14, 'fontweight', 'bold');
    text(fclk/3-14*fclk/100,-6,['fs=',num2str(fclk) 'MS/s' ], 'fontsize', 16, 'fontweight', 'bold');
    text(fclk/3-14*fclk/100,-14,['fin=',num2str(fin_dis) 'MHz' ], 'fontsize', 16, 'fontweight', 'bold');
    %text(fclk/4-14*fclk/100,-42,['SNR=' num2str(SNR_dis) 'dB'], 'fontsize', 14, 'fontweight', 'bold');
    text(fclk/3-14*fclk/100,-22,['SNDR=' num2str(SINAD_dis) 'dB' ], 'fontsize', 16, 'fontweight', 'bold');
    text(fclk/3-14*fclk/100,-30,['ENOB=' num2str(ENOB_dis) 'bit'], 'fontsize', 16, 'fontweight', 'bold');
    text(fclk/3-14*fclk/100,-38,['SFDR=' num2str(SFDR_dis) 'dB'], 'fontsize', 16, 'fontweight', 'bold');
    text(fclk/3-14*fclk/100,-46,['FFT--' num2str(numpt)], 'fontsize', 16, 'fontweight', 'bold')

    disp('----LX08D1000 ADC Testing Results----');
    dispstr=strcat('ENOB =',char(20),num2str(ENOB),' bits');
    disp(dispstr);
    dispstr=strcat('SINAD =',char(20),num2str(SINAD),' dB');
    disp(dispstr);
    dispstr=strcat('SNR =',char(20),num2str(SNR),' dB');
    disp(dispstr);
    dispstr=strcat('SFDR =',char(20),num2str(SFDR),' dB');
    disp(dispstr);
    disp('THD is calculated from 2nd through 5th order harmonics');
    dispstr=strcat('THD =',char(20),num2str(THD),' dB');
    disp(dispstr);
    dispstr=strcat('HD(1st~9th) in dB =',char(20),num2str(HD));
    disp(dispstr);

repo: SP


text:
%GM_PHD_Construct_Update_Components
%Last modified 12th September 2013
%Matlab code by Bryan Clarke b.clarke@acfr.usyd.edu.au

%This file creates the components needed for performing a Kalman filter update on the
%targets using the measurement.
s = sprintf('Step 3: Constructing update components for all targets, new and existing.');
disp(s);

%We need to clear the data structures each iteration
eta = [];
S = [];
K = [];
P_k_k = [];

for j = 1:numTargets_Jk_k_minus_1
    m_j = mk_k_minus_1(:,j);
    eta_j = H2 * m_j;%Observation model. Assume we see position AND velocity of the target.

    P_range = calculateDataRange4(j); %4x4 array

    PHt = Pk_k_minus_1(:,P_range) * H2'; %Taken from Tim Bailey's EKF code. 4x4 array

    %Calculate K via Tim Bailey's method.
    S_j = R2 + H2 * PHt;
    %At this point, Tim Bailey's code makes S_j symmetric. In this case, it leads to the matrix being non-positive definite a lot of the time and chol() crashes.
    %So we won't do that.
    SChol= chol(S_j);

    SCholInv= SChol \ eye(size(SChol)); % triangular matrix, invert via left division
    W1 = PHt * SCholInv;

    K_j = W1 * SCholInv';

    P_j = Pk_k_minus_1(:,P_range) - W1*W1';%4x4 array
    %End Tim Bailey's code.

    eta = [eta, eta_j];
    S = [S, S_j];
    K = [K, K_j];
    P_k_k = [P_k_k, P_j];
end


text:
function [varargout]=uniqueIntegerRow(F)

Fs=sort(F,2); %Sort so faces with same nodes have the same rows
[~,IND_F,IND_F_2]=unique(Fs,'rows');
F_uni=F(IND_F,:);

%%

varargout{1}=F_uni;
varargout{2}=IND_F;
varargout{3}=IND_F_2;

%%
% _*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 [uOutput] = scecdcimp(nFunction, sFilename)

% Filter function switchyard
if nFunction == 0     % Return info about filter
  uOutput = 'NIED (Kanto-Tokai) catalog with focal mechanism';
elseif nFunction == 2 % Return filename of help-file
  uOutput = 'import_NIED_doc.html';
elseif nFunction == 1 % Import and return catalog
  % Read formated data
  mData = textread(sFilename, '%s', 'delimiter', '\n', 'whitespace', '');
  % Create empty catalog
  uOutput = zeros(length(mData), 13);
  % Loop thru all lines of catalog and convert them
  for i = 1:length(mData)
    if rem(i,100) == 0 ; disp([ num2str(i) ' of ' num2str(length(mData)) ' events processed ']); end
    try
      l = find(mData{i} == ' ');
      mData{i}(l) = '0';
      uOutput(i,1) = str2num(mData{i}(37:44));  % Longitude
      uOutput(i,2) =  str2num(mData{i}(29:35)); % Latitude
      uOutput(i,3) = str2num(mData{i}(8:11));   % Year
      uOutput(i,4) = str2num(mData{i}(12:13));  % Month
      uOutput(i,5) = str2num(mData{i}(14:15));  % Day
      uOutput(i,6) = str2num(mData{i}(52:55));  % Magnitude
      uOutput(i,7) = str2num(mData{i}(46:50));  % Depth
      uOutput(i,8) = str2num(mData{i}(17:18));  % Hour
      uOutput(i,9) = str2num(mData{i}(20:21));  % Minute

      uOutput(i,10) = str2num(mData{i}(113:120));  % strike of nodal plane 1
      uOutput(i,11) = str2num(mData{i}(122:129));  % dip of nodal plane
      uOutput(i,12) = str2num(mData{i}(131:138));  % rake of nodal plane 1

      uOutput(i,13) = str2num(mData{i}(107:111));  % quality

      uOutput(i,14) = str2num(mData{i}(140:147));  % strike of nodal plane 2
      uOutput(i,15) = str2num(mData{i}(149:156));  % dip of nodal plane 2
      uOutput(i,16) = str2num(mData{i}(158:165));  % rake of nodal plane 2

      %Create decimal year
      uOutput(i,3) = decyear([uOutput(i,3) uOutput(i,4) uOutput(i,5) uOutput(i,8) uOutput(i,9)]);
    catch
      disp(['Import: Problem in line ' num2str(i) ' of ' sFilename '. Line ignored.']);
      uOutput(i,:) = uOutput(i,:)*nan;
    end
  end
  l = isnan(uOutput(:,1));
  uOutput(l,:) = [];
end

author: CelsoReyes


text:
function output = rdivide(input1, input2)

if ~iscell(input1)
  input1 = repmat({input1}, size(input2));
end
if ~iscell(input2)
  input2 = repmat({input2}, size(input1));
end
output = cellfun(@rdivide, input1, input2, 'uniformoutput', false);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function W = compute_mesh_weight(vertex,face,type,options)

% compute_mesh_weight - compute a weight matrix
%
%   W = compute_mesh_weight(vertex,face,type,options);
%
%   W is sparse weight matrix and W(i,j)=0 is vertex i and vertex j are not
%   connected in the mesh.
%
%   type is either
%       'combinatorial': W(i,j)=1 is vertex i is conntected to vertex j.
%       'distance': W(i,j) = 1/d_ij^2 where d_ij is distance between vertex
%           i and j.
%       'conformal': W(i,j) = cot(alpha_ij)+cot(beta_ij) where alpha_ij and
%           beta_ij are the adjacent angle to edge (i,j)
%
%   If options.normalize=1, the the rows of W are normalize to sum to 1.
%
%   Copyright (c) 2007 Gabriel Peyre

options.null = 0;
[vertex,face] = check_face_vertex(vertex,face);

nface = size(face,1);
n = max(max(face));

verb = getoptions(options, 'verb', n>5000);

if nargin<3
    type = 'conformal';
end

switch lower(type)
    case 'combinatorial'
        W = triangulation2adjacency(face);
    case 'distance'
        W = my_euclidean_distance(triangulation2adjacency(face),vertex);
        W(W>0) = 1./W(W>0);
        W = (W+W')/2;
    case 'conformal'
        % conformal laplacian
        W = sparse(n,n);
        ring = compute_vertex_face_ring(face);
        for i = 1:n
            if verb
                progressbar(i,n);
            end
            for b = ring{i}
                % b is a face adjacent to a
                bf = face(:,b);
                % compute complementary vertices
                if bf(1)==i
                    v = bf(2:3);
                elseif bf(2)==i
                    v = bf([1 3]);
                elseif bf(3)==i
                    v = bf(1:2);
                else
                    error('Problem in face ring.');
                end
                j = v(1); k = v(2);
                vi = vertex(:,i);
                vj = vertex(:,j);
                vk = vertex(:,k);
                % angles
                alpha = myangle(vk-vi,vk-vj);
                beta = myangle(vj-vi,vj-vk);
                % add weight
                W(i,j) = W(i,j) + cot( alpha );
                W(i,k) = W(i,k) + cot( beta );
            end
        end
    otherwise
        error('Unknown type.')
end

if isfield(options, 'normalize') && options.normalize==1
    W = diag(sum(W,2).^(-1)) * W;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function beta = myangle(u,v);

du = sqrt( sum(u.^2) );
dv = sqrt( sum(v.^2) );
du = max(du,eps); dv = max(dv,eps);
beta = acos( sum(u.*v) / (du*dv) );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function W = my_euclidean_distance(A,vertex)

if size(vertex,1)<size(vertex,2)
    vertex = vertex';
end

[i,j,s] = find(sparse(A));
d = sum( (vertex(i,:) - vertex(j,:)).^2, 2);
W = sparse(i,j,d);

author: gpeyre

repo: matlab-toolboxes


text:
function out = proj_hyperplane_box(x,a,b,l,u)
%PROJ_HYPERPLANE_BOX computes the orthogonal projection onto the intersection of a hyperplane and a box {x:<a,x>=b,l<=x<=u}
%
%  Usage:
%  out = PROJ_HYPERPLANE_BOX(x,a,b,[l],[u])
%  ===========================================
%  Input:
%  x - point to be projected (vector/matrix)
%  a - vector/matrix
%  b - scalar
%  l - lower bound (vector/matrix/scalar) [default: -inf]
%  u - upper bound (vector/matrix/scalar) [default: inf]
%  ===========================================
%  Assumptions:
%  The intersection of the hyperplane and the box is nonempty
%  ===========================================
%  Output:
%  out - projection vector

% This file is part of the FOM package - a collection of first order methods for solving convex optimization problems
% Copyright (C) 2017 Amir and Nili Beck
%
% 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/>.

if (nargin < 3)
    error ('usage: proj_hyperplane_box(x,a,b,[l],[u])') ;
end

if (nargin < 5)
    %upper bound is not given, setting to defalut value :inf
    u = inf ;
end
if ((nargin < 4) || (isempty( l)))
    %lower bound is not given, setting to defalut value :-inf
    l = -inf;
end

%checking that sum (l) < b < sum (u)

sumlb = trace(a'*((l .* (sign(a)>0)) + (u .* (sign(a)<0)))) ;
sumub = trace(a'*((u .* (sign(a)>0)) + (l .* (sign(a)<0)))) ;

if ((sumlb > b) || (any(any(l > u))) ||  (sumub < b))
    error('Set is infeasible') ;
end

%solve with equality
%defining f on lambda - to be used by the bisetion

eps = 1e-10 ;

f= @(lam)   trace(a'*min(max(x-lam*a,l),u))-b;
lambda_min = -1;
while(f(lambda_min)<0)
    lambda_min = lambda_min *2  ;
end

lambda_max = 1;
while(f(lambda_max)>0)
    lambda_max = lambda_max *2  ;
end

final_lam = bisection(f,lambda_min,lambda_max,eps) ;
out=  min(max(x-final_lam*a,l),u) ;

author: hiroyuki-kasai

repo: SGDLibrary


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

% ADD2DATE Add a yearfraction quantity V to a Date D
%   R = ADD2DATE(D,V) will add a double valued quantity V to
%   a date string or serial date number D, and return the new date string R
%
%   INPUT PARAMETERS:
%   D:  string or serial number, defining the date value
%   V:  double scalar, defining the quantity to add to D
%
%   RETURN PARAMETERS:
%   R:  returning the new date number R

%   ASSUMPTIONS:
%       V is given in terms of yearfractions
%       1 year equals 12 month
%       1 month equals 30 days
%       yearfractions < 1/30 are rounded to plus infinity

function R = add2date(D,V)

if ischar(D)
    dateval = datenum(D);
else
    dateval = D;
end

nr_month = floor(V*12);
nr_days = ceil((V*12-nr_month)*30);

tmpVec = zeros(size(V));
for i = 1:length(V)
   tmpVec(i) = addtodate(dateval,nr_month(i),'month');
end
R = daysadd(tmpVec,nr_days);

end


text:
function logical_mask = build_logical_mask(nt, maxNt, Ns)

%   Copyright (C) 2009 Cesare Magri
%   Version: 1.0.4

% -------
% LICENSE
% -------
% This software is distributed free under the condition that:
%
% 1. it shall not be incorporated in software that is subsequently sold;
%
% 2. the authorship of the software shall be acknowledged and the following
%    article shall be properly cited in any publication that uses results
%    generated by the software:
%
%      Magri C, Whittingstall K, Singh V, Logothetis NK, Panzeri S: A
%      toolbox for the fast information analysis of multiple-site LFP, EEG
%      and spike train recordings. BMC Neuroscience 2009 10(1):81;
%
% 3.  this notice shall remain in place in each source file.

column_increment = 0:maxNt:(Ns*maxNt)-1;

mask = zeros(maxNt, Ns);
mask(column_increment + nt(:).') = 1;
mask(end:-1:1,:) = cumsum(mask(end:-1:1,:));
% Converting 1 to logical "true":
logical_mask = mask>0;


text:
function varargout = eventlist2erpGUI(varargin)
% EVENTLIST2ERPGUI M-file for eventlist2erpGUI.fig
%      EVENTLIST2ERPGUI, by itself, creates a new EVENTLIST2ERPGUI or raises the existing
%      singleton*.
%
%      H = EVENTLIST2ERPGUI returns the handle to a new EVENTLIST2ERPGUI or the handle to
%      the existing singleton*.
%
%      EVENTLIST2ERPGUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in EVENTLIST2ERPGUI.M with the given input arguments.
%
%      EVENTLIST2ERPGUI('Property','Value',...) creates a new EVENTLIST2ERPGUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before eventlist2erpGUI_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to eventlist2erpGUI_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help eventlist2erpGUI

% Last Modified by GUIDE v2.5 26-Aug-2013 12:59:22

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
        'gui_Singleton',  gui_Singleton, ...
        'gui_OpeningFcn', @eventlist2erpGUI_OpeningFcn, ...
        'gui_OutputFcn',  @eventlist2erpGUI_OutputFcn, ...
        'gui_LayoutFcn',  [] , ...
        'gui_Callback',   []);
if nargin && ischar(varargin{1})
        gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
        [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
        gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

%--------------------------------------------------------------------------
function eventlist2erpGUI_OpeningFcn(hObject, eventdata, handles, varargin)
% Choose default command line output for eventlist2erpGUI
handles.output = hObject;
ERP = varargin{1};

if isfield(ERP, 'EVENTLIST')
        if ~isempty(ERP.EVENTLIST)
                nel = length(ERP.EVENTLIST);
        else
                title      = 'ERPLAB: eventlist2erpGUI()';
                question{1} = 'EVENTLIST structure is empty in this ERPset';
                question{2} = 'Do you want to continue anyway?';
                button = askquest(question, title);

                if ~strcmpi(button,'yes')
                        disp('User selected Cancel')
                        return
                end
                nel = 0;
        end
else
        title      = 'ERPLAB: eventlist2erpGUI()';
        question{1} = 'EVENTLIST structure was not found in this ERPset';
        question{2} = 'Do you want to continue anyway?';
        button = askquest(question, title);

        if ~strcmpi(button,'yes')
                disp('User selected Cancel')
                return
        end
        nel = 0;
end

handles.nel = nel;

if nel==0
        text_message1 = sprintf('This ERPset contains %g EVENTLIST', nel);
        text_message2 = '';
        set(handles.text_message1, 'String', text_message1)
        set(handles.radiobutton_make1, 'Enable', 'off')
        set(handles.radiobutton_replace, 'Enable', 'off')
        set(handles.edit_replace, 'Enable', 'off')
        set(handles.text_replace, 'Enable', 'off')
elseif nel==1
        text_message1 = sprintf('This ERPset contains %g EVENTLIST', nel);
        text_message2 = 'What would you like to do with the imported EVENTLIST?';
        set(handles.text_message1, 'String', text_message1)
        set(handles.radiobutton_replace, 'Enable', 'off')
        set(handles.edit_replace, 'Enable', 'off')
        set(handles.text_replace, 'Enable', 'off')
else
        text_message1 = sprintf('This ERPset contains %g EVENTLISTs (1-%g)', nel, nel);
        text_message2 = 'What would you like to do with the imported EVENTLIST?';
        set(handles.text_message1, 'String', text_message1)
end

set(handles.text_message2, 'String', text_message2)
text_message3 = sprintf('Append imported EVENTLIST as #%g', nel+1);
set(handles.radiobutton_append, 'String', text_message3)
set(handles.radiobutton_append, 'Value', 1)

%
% Name & version
%
version = geterplabversion;
set(handles.gui_chassis,'Name', ['ERPLAB ' version '   -   Import EVENTLIST from Text to ERP GUI'])

%
% Color GUI
%
handles = painterplab(handles);

%
% Set font size
%
handles = setfonterplab(handles);

% Update handles structure
guidata(hObject, handles);

% help
% helpbutton

% UIWAIT makes eventlist2erpGUI wait for user response (see UIRESUME)
uiwait(handles.gui_chassis);

%--------------------------------------------------------------------------
function varargout = eventlist2erpGUI_OutputFcn(hObject, eventdata, handles)

varargout{1} = handles.output;

% The figure can be deleted now
delete(handles.gui_chassis);
pause(0.5)

%--------------------------------------------------------------------------
function radiobutton_make1_Callback(hObject, eventdata, handles)

if get(hObject, 'Value')
        set(handles.radiobutton_replace, 'Value', 0)
        set(handles.radiobutton_append, 'Value', 0)
else
        set(handles.radiobutton_make1, 'Value', 1)
end

%--------------------------------------------------------------------------
function radiobutton_replace_Callback(hObject, eventdata, handles)

if get(hObject, 'Value')
        set(handles.radiobutton_make1, 'Value', 0)
        set(handles.radiobutton_append, 'Value', 0)
else
        set(handles.radiobutton_replace, 'Value', 1)
end

%--------------------------------------------------------------------------
function radiobutton_append_Callback(hObject, eventdata, handles)

if get(hObject, 'Value')
        set(handles.radiobutton_replace, 'Value', 0)
        set(handles.radiobutton_make1, 'Value', 0)
else
        set(handles.radiobutton_append, 'Value', 1)
end

%--------------------------------------------------------------------------
function pushbutton_cancel_Callback(hObject, eventdata, handles)

handles.output = [];

% Update handles structure
guidata(hObject, handles);
uiresume(handles.gui_chassis);

%--------------------------------------------------------------------------
function pushbutton_OK_Callback(hObject, eventdata, handles)

eventlistname = get(handles.edit_eventlistname, 'String');
nel = handles.nel;

if isempty(strtrim(eventlistname))
        msgboxText{1} =  'Error: You must choose an EVENTLIST text file.';
        title = 'ERPLAB: eventlist2erpGUI';
        errorfound(msgboxText, title);
        return
end

if get(handles.radiobutton_make1, 'Value')
        option = 0;
        index  = [];

elseif get(handles.radiobutton_replace, 'Value')

        option = 1;
        index = str2num(get(handles.edit_replace, 'String'));

        if isempty(index) || index<1 || index>nel
                msgboxText{1} =  'Error: invalid EVENTLIST index.';
                title = 'ERPLAB: eventlist2erpGUI';
                errorfound(msgboxText, title);
                return
        end

elseif get(handles.radiobutton_append, 'Value')
        option = 2;
        index  = [];

else
        return
end

outputarray = {eventlistname, option, index};
handles.output = outputarray;

% Update handles structure
guidata(hObject, handles);
uiresume(handles.gui_chassis);

%--------------------------------------------------------------------------
function edit_replace_Callback(hObject, eventdata, handles)

%--------------------------------------------------------------------------
function edit_replace_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
end

%--------------------------------------------------------------------------
function edit_eventlistname_Callback(hObject, eventdata, handles)

%--------------------------------------------------------------------------
function edit_eventlistname_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
end

%--------------------------------------------------------------------------
function pushbutton_browse_Callback(hObject, eventdata, handles)

[bdfilename,bdfpathname] = uigetfile({'*.txt';'*.dat';'*.*'},'Select an EVENTLIST text file');

if isequal(bdfilename,0)
        disp('User selected Cancel')
        return
else
        set(handles.edit_eventlistname, 'String', fullfile(bdfpathname, bdfilename));

        % Update handles structure
        guidata(hObject, handles);

        flinkname = fullfile(bdfpathname, bdfilename);
        disp(['For EVENTLIST text file, user selected ' flinkname ''])
end

%--------------------------------------------------------------------------
function gui_chassis_CloseRequestFcn(hObject, eventdata, handles)

if isequal(get(handles.gui_chassis, 'waitstatus'), 'waiting')

        handles.output = [];
        %Update handles structure
        guidata(hObject, handles);
        uiresume(handles.gui_chassis);
else
        % The GUI is no longer waiting, just close it
        delete(handles.gui_chassis);
end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function value = r8mat_sum ( m, n, a )

%*****************************************************************************80
%
%% R8MAT_SUM returns the sum of the entries of an R8MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, real A(M,N), the array.
%
%    Output, real VALUE, the sum of the entries.
%
  value = sum ( sum ( a(1:m,1:n) ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function a = r8vec_permute ( n, a, p )

%*****************************************************************************80
%
%% R8VEC_PERMUTE permutes an R8VEC in place.
%
%  Discussion:
%
%    This routine permutes an array of real "objects", but the same
%    logic can be used to permute an array of objects of any arithmetic
%    type, or an array of objects of any complexity.  The only temporary
%    storage required is enough to store a single object.  The number
%    of data movements made is N + the number of cycles of order 2 or more,
%    which is never more than N + N/2.
%
%    P(I) = J means that the I-th element of the output array should be
%    the J-th element of the input array.  P must be a legal permutation
%    of the integers from 1 to N, otherwise the algorithm will
%    fail catastrophically.
%
%  Example:
%
%    Input:
%
%      N = 5
%      P = (   2,   4,   5,   1,   3 )
%      A = ( 1.0, 2.0, 3.0, 4.0, 5.0 )
%
%    Output:
%
%      A    = ( 2.0, 4.0, 5.0, 1.0, 3.0 ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 May 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of objects.
%
%    Input, real A(N), the array to be permuted.
%
%    Input, integer P(N), the permutation.
%
%    Output, real A(N), the permuted array.
%
  ierror = perm1_check ( n, p );

  if ( ierror ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8VEC_PERMUTE - Fatal error!\n' );
    fprintf ( 1, '  PERM1_CHECK finds the permuation illegal.\n' );
    error ( 'R8VEC_PERMUTE - Fatal error!' );
  end
%
%  Search for the next element of the permutation that has not been used.
%
  for istart = 1 : n

    if ( p(istart) < 0 )

      continue

    elseif ( p(istart) == istart )

      p(istart) = -p(istart);
      continue

    else

      a_temp = a(istart);
      iget = istart;
%
%  Copy the new value into the vacated entry.
%
      while ( 1 )

        iput = iget;
        iget = p(iget);

        p(iput) = -p(iput);

        if ( iget < 1 || n < iget )
          fprintf ( 1, '\n' );
          fprintf ( 1, 'R8VEC_PERMUTE - Fatal error!\n' );
          fprintf ( 1, '  A permutation index is out of range.\n' );
          fprintf ( 1, '  P(%d) = %d\n', iput, iget );
          error ( 'R8VEC_PERMUTE - Fatal error!' );
        end

        if ( iget == istart )
          a(iput) = a_temp;
          break
        end

        a(iput) = a(iget);

      end

    end

  end
%
%  Restore the signs of the entries.
%
% p(1:n) = -p(1:n);

  return
end


text:
function fiff_end_file(fid)
%
% fiff_end_file(fid)
%
% Writes the closing tags to a fif file and closes the file
%
%     fid           An open fif file descriptor
%

%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%   Revision 1.3  2006/04/23 15:29:40  msh
%   Added MGH to the copyright
%
%   Revision 1.2  2006/04/10 23:26:54  msh
%   Added fiff reading routines
%
%   Revision 1.1  2005/12/05 16:01:04  msh
%   Added an initial set of fiff writing routines.
%

global FIFF;
if isempty(FIFF)
    FIFF = fiff_define_constants();
end

me = 'MNE:fiff_end_file';

if nargin ~= 1
    error(me, 'An open file id required as an argument');
end

datasize = 0;
count = fwrite(fid, int32(FIFF.FIFF_NOP), 'int32');
if count ~= 1
    error(me, 'write failed');
end
count = fwrite(fid, int32(FIFF.FIFFT_VOID), 'int32');
if count ~= 1
    error(me, 'write failed');
end
count = fwrite(fid, int32(datasize), 'int32');
if count ~= 1
    error(me, 'write failed');
end
count = fwrite(fid, int32(FIFF.FIFFV_NEXT_NONE), 'int32');
if count ~= 1
    error(me, 'write failed');
end
fclose(fid);

return


text:
%%Ex. 8 Extracting an individual element of an array

a = [3 6 7];
b = [1 9 4 5];
c = a(2) + b(4)

%Output:  c = 11

author: TheAlgorithms

repo: MATLAB-Octave


text:
function table = l4vec_data_read ( input_filename,  n )

%*****************************************************************************80
%
%% L4VEC_DATA_READ reads data from an L4VEC file.
%
%  Discussion:
%
%    An L4VEC is a vector of L4's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer N, the number of columns in the data.
%
%    Output, integer TABLE(N), the data.
%
  table = zeros ( n, 1 );

  i = 0;

  while ( i < n )

    line = fgets ( input_unit );

    if ( line == -1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'L4VEC_DATA_READ - Error!\n' );
      fprintf ( 1, '  End of input while reading data.\n' );
      error ( 'L4VEC_DATA_READ - Error!' );
    end

    if ( line(1) == '#' )

    elseif ( s_len_trim ( line ) == 0 )

    else

      i = i + 1;
      table(i) = s_to_l4 ( line );

    end

  end

  fclose ( input_unit );

  return
end


text:
function paths = zero_jump_link(frames, downweight_notrack)
% AUTORIGHTS
% ---------------------------------------------------------
% Copyright (c) 2015, Georgia Gkioxari
%
% This file is part of the Action Tubes code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.

% at each iteration in the while loop, a new path is generated and the boxes, scores and box boxes_idx
% associated with this path is removed so that when the next path will be generated it will consider the remaining
% boxes. e.g. from frame 1 , out of 5 boxes, box no 2 is used to generate path 1, so we
% remove the box from the list, so when we will generate path 2, we will not use this box 2 again

% This version is corrected to keep track of the correct boxes indices in each
% path , as the boxes are removed after generating each path , to keep track of the indices
% of the boxes following lines are added :

% index = frames(V(j)).boxes_idx(id);
% index =  [index;  frames(V(j+1)).boxes_idx(id)]; %% newly added
% frames(V(j)).boxes_idx(id)  = []; %% newly added
% ---------------------------------------------------------

% set of vertices
num_frames = length(frames);
V = 1:num_frames;
isempty_vertex = false(length(V),1);
if nargin < 2
  downweight_notrack = false;
end
% CHECK: make sure data is not empty
for i=1:num_frames
    if isempty(frames(i).boxes)
        disp(i);
        error('Empty frame');
    end
end

%%--- dymamic programming (K-paths)
path_counter = 0;

while ~any(isempty_vertex)

    %fprintf('# of iteration = %d\n',path_counter);
    % step1: initialize data structures
    T = length(V);
    for i=1:T
        num_states = size(frames(V(i)).boxes,1);
        data(i).scores = zeros(num_states,1);
        data(i).index  = nan(num_states,1);
        data(i).track_scores = nan(num_states,1);

    end

    % step2: solve viterbi -- forward pass --
    % we do the forward pass from the frame before the last frame (T-1),
    % computing the edge scores between the boxes at last frame (T) and the frame before the last frame (T-1), and continue this till we
    % reach the first frame,
    % at the first iteration when we compute the edge scores between frames
    % T-1 and frame T and add it with the data(i+1).scores', the
    % data(i+1).scores' values are all 0, this is the case only at first
    % iteration, and at the end of the first iteration data(i+1).scores'
    % are initialised with the max edge scores for each boxes in the frame
    % T-1, in the subsequent iteration, we add the max edge scores of boxes from i+1-th frame to the edge scores of boxes of i-th frame
    % and then we compute the max_edge_scores for i-th frame and store them in data(i).scores and the indices of the boxes havining max_edge_scores
    % with the i+1 frame boxes are stored in data(i).index
    for i=T-1:-1:1
        [edge_score, data(i).track_scores]  = score_of_edge(frames(V(i)),frames(V(i+1)), downweight_notrack);
        edge_score = bsxfun(@plus,edge_score,data(i+1).scores');
        [data(i).scores, data(i).index] = max(edge_score,[],2);
    end
    % step3: decode -- backward pass of Viterbi - back tracing
    path_counter = path_counter+1;
    [s, si] = sort(data(1).scores,'descend'); % sort the max_edge_scores of boxes in decreasing order in the first frame
    id = si(1); % index of the box having maximum edge score in the first frame
    score = data(1).scores(id);
    %index = id;
    index = frames(V(1)).boxes_idx(id);
    boxes = frames(V(1)).boxes(id,:);
    scores = frames(V(1)).scores(id,1);
    for j=1:T-1
        id = data(j).index(id);
        index =  [index;  frames(V(j+1)).boxes_idx(id)]; %% newly added
        boxes =  [boxes;  frames(V(j+1)).boxes(id,:)];
        scores = [scores; frames(V(j+1)).scores(id,:)];
    end
    paths(path_counter).total_score = score/num_frames;
    paths(path_counter).idx   = index;
    paths(path_counter).boxes = [boxes scores];

    [top_scores] = sort(scores,'descend'); % sort the max_edge_scores of boxes in decreasing order in the first frame
    mean_top_scrs = mean(top_scores(1:ceil(numel(top_scores)*0.5)));
    g_filter = [1 4 6 4 1]/16;
    smooth_scores = imfilter(scores(:),g_filter(:),'symmetric');
    paths(path_counter).smooth_scores = smooth_scores ;

    scores = scores +mean_top_scrs;
    paths(path_counter).scores = scores ;

    % step4: remove covered boxes
    for j=1:T
        id = paths(path_counter).idx(j);
        id11 = id==frames(V(j)).boxes_idx; % box id 2 remove, because it is included in this path
        id2Rem = find(id11); % % box id 2 remove
        frames(V(j)).boxes(id2Rem,:) = [];
        %frames(V(j)).feat(id2Rem,:)  = [];
        frames(V(j)).scores(id2Rem)  = [];
        frames(V(j)).boxes_idx(id2Rem)  = []; %% newly added
        isempty_vertex(j) = isempty(frames(V(j)).boxes);
    end

end

% -------------------------------------------------------------------------
function [score, track_score] = score_of_edge(v1,v2, downweight_notrack)
% -------------------------------------------------------------------------

N1 = size(v1.boxes,1);
N2 = size(v2.boxes,1);
score = nan(N1,N2);
track_score = nan(N1,N2);

for i1=1:N1
    % scores
    scores2 = v2.scores;
    scores1 = v1.scores(i1);
    score(i1,:) = scores1+scores2';
end
if isfield(v2, 'trackedboxes') && ~isempty(v1.trackedboxes)
  if numel(v2.trackedboxes) < 3
  overlap_ratio_1 = get_overlap_MtoN(v1.boxes, v1.trackedboxes{1,1});
  overlap_ratio_2 = get_overlap_MtoN(v2.boxes, v1.trackedboxes{1,2});
  track_score = round(overlap_ratio_1) * round(overlap_ratio_2');
  score(track_score > 0) = score(track_score > 0)  + 1 ;

  if downweight_notrack
      score(~track_score) = score(~track_score) -1;
  end
  track_score = double(track_score > 0);

  else
    if ~isempty(v1.trackedboxes)
      overlap_ratio_a1 = get_overlap_MtoN(v1.boxes, v1.trackedboxes{1,2});
      overlap_ratio_a2 = get_overlap_MtoN(v2.boxes, v1.trackedboxes{1,3});
      track_12 = round(overlap_ratio_a1) * round(overlap_ratio_a2');
    else
      track_12 = zeros(size(score));
    end
    overlap_ratio_b1 = get_overlap_MtoN(v1.boxes, v2.trackedboxes{1,1});
    overlap_ratio_b2 = get_overlap_MtoN(v2.boxes, v2.trackedboxes{1,2});
    track_21 = round(overlap_ratio_b1) * round(overlap_ratio_b2');
  score(track_12 > 0) = score(track_12 > 0)  + .5 ;
  score(track_21 > 0) = score(track_21 > 0)  + .5 ;

  track_score = double(track_12 > 0) + double(track_21 > 0);
  end
end

% -------------------------------------------------------------------------
function iou = inters_union(bounds1,bounds2)
% -------------------------------------------------------------------------

inters = rectint(bounds1,bounds2);
ar1 = bounds1(:,3).*bounds1(:,4);
ar2 = bounds2(:,3).*bounds2(:,4);
union = bsxfun(@plus,ar1,ar2')-inters;

iou = inters./(union+eps);

author: feichtenhofer

repo: Detect-Track


text:
%DEMO_MODELASSESMENT2  Demonstration for model assessment when the observation
%                      model is non-Gaussian
%
%  Description
%    We will consider the classification problem in demo_classific.
%    The analysis is conducted with full Gaussian process using
%    both probit and logit likelihood. The performance of these two
%    models are compared by evaluating the ten-fold
%    cross-validation, leave-one-out cross-validation, WAIC, DIC
%    and the effective number of parameters The inference will be
%    conducted using maximum a posterior (MAP) estimate for the
%    parameters using EP and Laplace approximation, via full Markov
%    chain Monte Carlo (MCMC) and with an integration approximation
%    (IA) for the parameters.
%
%    This demo is organised in two parts:
%     1) data analysis with with probit likelihood
%     2) data analysis with with logit likelihood
%
%  See also
%    DEMO_CLASSIFIC1, DEMO_MODELASSESMENT1
%
% Copyright (c) 2009-2010 Jarno Vanhatalo
% Copyright (c) 2010 Aki Vehtari

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

% =====================================
% 1) data analysis with probit likelihood
% =====================================
disp('Data analysis with probit likelihood')
S = which('demo_classific');
L = strrep(S,'demo_classific.m','demodata/synth.tr');
x=load(L);
y=x(:,end);
y = 2.*y-1;
x(:,end)=[];
[n, nin] = size(x);

DIC=repmat(NaN,1,8);DIC2=repmat(NaN,1,8);DIC_latent=repmat(NaN,1,8);
p_eff=repmat(NaN,1,8);p_eff2=repmat(NaN,1,8);p_eff_latent=repmat(NaN,1,8);p_eff_latent2=repmat(NaN,1,8);

% Create covariance functions
gpcf = gpcf_sexp('lengthScale', [0.9 0.9], 'magnSigma2', 2);

% Set the prior for the parameters of covariance functions
pl = prior_logunif();
gpcf = gpcf_sexp(gpcf, 'lengthScale_prior', pl,'magnSigma2_prior', pl); %

% Create the GP structure
gp = gp_set('lik', lik_probit, 'cf', gpcf, 'jitterSigma2', 1e-4);

% ------- Laplace approximation --------
disp(['Probit with Laplace integration over the latent values '; ...
      'and MAP estimate for the parameters                    '])

% Set the approximate inference method
gp = gp_set(gp, 'latent_method', 'Laplace');

n=length(y);

% Set the options for the optimization
opt=optimset('TolFun',1e-3,'TolX',1e-3);
% Optimize with the quasi-Newton method
gp=gp_optim(gp,x,y,'opt',opt,'optimf',@fminlbfgs);

% Evaluate the effective number of parameters and DIC with focus on
% latent variables.
models{1} = 'pr_Laplace';
p_eff_latent(1) = gp_peff(gp, x, y);
[DIC_latent(1), p_eff_latent2(1)] = gp_dic(gp, x, y);
WAIC(1) = gp_waic(gp,x,y);

% Evaluate the 10-fold cross-validation results.
cvres = gp_kfcv(gp, x, y, 'display', 'fold');
mlpd_cv(1) = cvres.mlpd_cv;

% Evaluate the leave-one-out cross-validation results.
[Ef,Varf,lpy] =  gp_loopred(gp, x, y);
mlpd_loo(1) = mean(lpy);

% ------- Expectation propagation --------
disp(['Probit with EP integration over the latent values and MAP '; ...
      'estimate for the parameters                               '])

% Set the approximate inference method
gp = gp_set(gp, 'latent_method', 'EP');

% Set the options for the optimization
opt=optimset('TolFun',1e-3,'TolX',1e-3);
% Optimize with the BFGS quasi-Newton method
gp=gp_optim(gp,x,y,'opt',opt,'optimf',@fminlbfgs);

% Evaluate the effective number of parameters and DIC with focus on
% latent variables.
models{2} = 'pr_EP';
p_eff_latent(2) = gp_peff(gp, x, y) ;
[DIC_latent(2), p_eff_latent2(2)] = gp_dic(gp, x, y);
WAIC(2) = gp_waic(gp,x,y);

% Evaluate the 10-fold cross-validation results.
cvres = gp_kfcv(gp, x, y, 'display', 'fold');
mlpd_cv(2) = cvres.mlpd_cv;

% Evaluate the leave-one-out cross-validation results.
[Ef,Varf,lpy] =  gp_loopred(gp, x, y);
mlpd_loo(2) = mean(lpy);

% ------- MCMC ---------------
disp(['Probit with MCMC integration over the latent values and '; ...
      'the parameters                                          '])

% Set the approximate inference method
gp = gp_set(gp, 'latent_method', 'MCMC');

% Sample
mcopt.nsamples=220;mcopt.display=20;
[rgp,gp]=gp_mc(gp, x, y, mcopt);
rgp=thin(rgp, 21, 2);

% Evaluate the effective number of parameters and DIC with focus on
% latent variables.
models{3} = 'pr_MCMC';
[DIC(3), p_eff(3)] =  gp_dic(rgp, x, y, 'focus', 'hyper');
[DIC2(3), p_eff2(3)] =  gp_dic(rgp, x, y);
WAIC(3) = gp_waic(rgp,x,y);

% Evaluate the 10-fold cross-validation results.
mcopt.nsamples=50;mcopt.display=20;
cvres = gp_kfcv(gp, x, y, 'inf_method', 'MCMC', 'opt', mcopt, 'display', 'fold');
mlpd_cv(3) = cvres.mlpd_cv;

% Evaluate the leave-one-out cross-validation results.
[Ef,Varf,lpy] =  gp_loopred(rgp, x, y);
mlpd_loo(3) = mean(lpy);

% --- Integration approximation approach ---
disp(['Probit with EP integration over the latent values and '; ...
      'grid integration over the parameters                  '])

% Use EP
gp = gp_set(gp, 'latent_method', 'EP');

% Set the options for the optimization
opt=optimset('TolFun',1e-3,'TolX',1e-3);
% Optimize with the BFGS quasi-Newton method
gp=gp_optim(gp,x,y,'opt',opt,'optimf',@fminlbfgs);

% now perform the integration
clear opt
opt.int_method = 'grid';
opt.step_size = 2;
gp_array = gp_ia(gp, x, y, opt);

models{4} = 'pr_IA';
[DIC(4), p_eff(4)] =  gp_dic(gp_array, x, y, 'focus', 'hyper');
[DIC2(4), p_eff2(4)] =  gp_dic(gp_array, x, y);
WAIC(4) = gp_waic(gp_array,x,y);

% Then the 10 fold cross-validation.
cvres = gp_kfcv(gp, x, y, 'inf_method', 'IA', 'opt', opt, 'display', 'fold');
mlpd_cv(4) = cvres.mlpd_cv;

% Evaluate the leave-one-out cross-validation results.
[Ef,Varf,lpy] =  gp_loopred(gp_array, x, y);
mlpd_loo(4) = mean(lpy);

% =====================================
% 2) data analysis with logit likelihood
% =====================================
disp('Data analysis with logit likelihood')

S = which('demo_classific');
L = strrep(S,'demo_classific.m','demodata/synth.tr');
x=load(L);
y=x(:,end);
y = 2.*y-1;
x(:,end)=[];
[n, nin] = size(x);

% Create covariance functions
gpcf = gpcf_sexp('lengthScale', [0.9 0.9], 'magnSigma2', 2);

% Set the prior for the parameters of covariance functions
pl = prior_logunif();
gpcf = gpcf_sexp(gpcf, 'lengthScale_prior', pl,'magnSigma2_prior', pl); %

% Create the likelihood structure
lik = ('init');

% Create the GP structure
gp = gp_set('lik', lik_logit, 'cf', gpcf, 'jitterSigma2', 1e-4);

% ------- Laplace approximation --------
disp(['Logit with Laplace integration over the latent values and '; ...
      'MAP estimate for the parameters                           '])

% Set the approximate inference method
gp = gp_set(gp, 'latent_method', 'Laplace');

% Set the options for the optimization
opt=optimset('TolFun',1e-3,'TolX',1e-3);
% Optimize with the BFGS quasi-Newton method
gp=gp_optim(gp,x,y,'opt',opt,'optimf',@fminlbfgs);

% Evaluate the effective number of parameters and DIC with focus on
% latent variables.
models{5} = 'lo_Laplace';
p_eff_latent(5) = gp_peff(gp, x, y);
[DIC_latent(5), p_eff_latent2(5)] = gp_dic(gp, x, y);
WAIC(5) = gp_waic(gp,x,y);

% Evaluate the 10-fold cross-validation results.
cvres = gp_kfcv(gp, x, y, 'display', 'fold');
mlpd_cv(5) = cvres.mlpd_cv;

% Evaluate the leave-one-out cross-validation results.
[Ef,Varf,lpy] =  gp_loopred(gp, x, y);
mlpd_loo(5) = mean(lpy);

% ------- Expectation propagation --------
disp(['Logit with EP integration over the latent values and MAP'; ...
      'estimate for the parameters                             '])

% Set the approximate inference method
gp = gp_set(gp, 'latent_method', 'EP');

% Set the options for the optimization
opt=optimset('TolFun',1e-3,'TolX',1e-3);
% Optimize with the BFGS quasi-Newton method
gp=gp_optim(gp,x,y,'opt',opt,'optimf',@fminlbfgs);

% Evaluate the effective number of parameters and DIC with focus on
% latent variables.
models{6} = 'lo_EP';
p_eff_latent(6) = gp_peff(gp, x, y) ;
[DIC_latent(6), p_eff_latent2(6)] = gp_dic(gp, x, y);
WAIC(6) = gp_waic(gp,x,y);

% Evaluate the 10-fold cross-validation results.
cvres = gp_kfcv(gp, x, y, 'display', 'fold');
mlpd_cv(6) = cvres.mlpd_cv;

% Evaluate the leave-one-out cross-validation results.
[Ef,Varf,lpy] =  gp_loopred(gp, x, y);
mlpd_loo(6) = mean(lpy);

% ------- MCMC ---------------
disp(['Logit with MCMC integration over the latent values and '; ...
      'the parameters                                         '])

% Set the approximate inference method
gp = gp_set(gp, 'latent_method', 'MCMC');

% Sample
mcopt.nsamples=200;mcopt.display=20;
[rgp,gp] = gp_mc(gp, x, y, mcopt);
rgp=thin(rgp, 21, 2);

% Evaluate the effective number of parameters and DIC with focus on latent variables.
models{7} = 'lo_MCMC';
[DIC(7), p_eff(7)] =  gp_dic(rgp, x, y, 'focus', 'hyper');
[DIC2(7), p_eff2(7)] =  gp_dic(rgp, x, y);
WAIC(7) = gp_waic(rgp,x,y);

% Evaluate the 10-fold cross-validation results.
mcopt.nsamples=50;mcopt.display=20;
cvres = gp_kfcv(gp, x, y, 'inf_method', 'MCMC', 'opt', mcopt, 'display', 'fold');
mlpd_cv(7) = cvres.mlpd_cv;

% Evaluate the leave-one-out cross-validation results.
[Ef,Varf,lpy] =  gp_loopred(rgp, x, y);
mlpd_loo(7) = mean(lpy);

% --- Integration approximation approach ---
disp(['Logit with EP integration over the latent values and grid '; ...
      'integration over the parameters                           '])

% Use EP
gp = gp_set(gp, 'latent_method', 'EP');

% Set the options for the optimization
opt=optimset('TolFun',1e-3,'TolX',1e-3);
% Optimize with the BFGS quasi-Newton method
gp=gp_optim(gp,x,y,'opt',opt,'optimf',@fminlbfgs);

% now perform the integration
clear opt
opt.int_method = 'grid';
opt.step_size = 2;
gp_array = gp_ia(gp, x, y, opt);

models{8} = 'lo_IA';
[DIC(8), p_eff(8)] =  gp_dic(gp_array, x, y, 'focus', 'hyper');
[DIC2(8), p_eff2(8)] =  gp_dic(gp_array, x, y);
WAIC(8) = gp_waic(gp_array,x,y);

% Then the 10 fold cross-validation.
cvres = gp_kfcv(gp, x, y, 'inf_method', 'IA', 'opt', opt, 'display', 'fold');
mlpd_cv(8) = cvres.mlpd_cv;

% Evaluate the leave-one-out cross-validation results.
[Ef,Varf,lpy] =  gp_loopred(gp_array, x, y);
mlpd_loo(8) = mean(lpy);

%========================================================
% PART 4 Print the results
%========================================================
disp('Summary of the results')

S = '       ';
for i = 1:length(models)
    S = [S '  ' models{i}];
end

S = sprintf([S '\n CV-mlpd  %6.2f    %6.2f  %6.2f  %6.2f   %6.2f    %6.2f  %6.2f  %6.2f'], mlpd_cv);
S = sprintf([S '\n LOO-mlpd %6.2f    %6.2f  %6.2f  %6.2f   %6.2f    %6.2f  %6.2f  %6.2f'], mlpd_loo);
S = sprintf([S '\n ']);
S = sprintf([S '\n WAIC     %6.2f    %6.2f  %6.2f  %6.2f   %6.2f    %6.2f  %6.2f  %6.2f'], WAIC);
S = sprintf([S '\n ']);
S = sprintf([S '\n DIC_h    %6.2f    %6.2f  %6.2f  %6.2f   %6.2f    %6.2f  %6.2f  %6.2f'], DIC);
S = sprintf([S '\n DIC_a    %6.2f    %6.2f  %6.2f  %6.2f   %6.2f    %6.2f  %6.2f  %6.2f'], DIC2);
S = sprintf([S '\n DIC_l    %6.2f    %6.2f  %6.2f  %6.2f   %6.2f    %6.2f  %6.2f  %6.2f'], DIC_latent);
S = sprintf([S '\n peff_h   %6.2f    %6.2f  %6.2f  %6.2f   %6.2f    %6.2f  %6.2f  %6.2f'], p_eff);
S = sprintf([S '\n peff_a   %6.2f    %6.2f  %6.2f  %6.2f   %6.2f    %6.2f  %6.2f  %6.2f'], p_eff2);
S = sprintf([S '\n peff_l   %6.2f    %6.2f  %6.2f  %6.2f   %6.2f    %6.2f  %6.2f  %6.2f'], p_eff_latent);
S = sprintf([S '\n peff_l2  %6.2f    %6.2f  %6.2f  %6.2f   %6.2f    %6.2f  %6.2f  %6.2f'], p_eff_latent2);
S = sprintf([S '\n ']);
S = sprintf([S '\n The notation is as follows:']);
S = sprintf([S '\n pr_*     = probit likelihood and inference method']);
S = sprintf([S '\n lo_*     = logit likelihood and inference method']);
S = sprintf([S '\n CV-mlpd  = mean log predictive density from the 10-fold CV. ']);
S = sprintf([S '\n LOO-mlpd = mean log predictive density from the 10-fold CV. ']);
S = sprintf([S '\n WAIC     = Widely applicable information criterion. ']);
S = sprintf([S '\n DIC_h    = DIC with focus on parameters. ']);
S = sprintf([S '\n DIC_a    = DIC with focus on parameters and laten variables (all). ']);
S = sprintf([S '\n DIC_l    = DIC with focus on latent variables. ']);
S = sprintf([S '\n peff_h   = effective number of parameters (latent variables marginalized). ']);
S = sprintf([S '\n peff_a   = effective number of parameters and latent variables. ']);
S = sprintf([S '\n peff_l   = effective number of latent variables evaluated with gp_peff. ']);
S = sprintf([S '\n peff_l2  = effective number of latent variables evaluated with gp_dic. ']);
S = sprintf([S '\n '])

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
% The COBRAToolbox: testWriteCbModel.m
%
% Purpose:
%     - test the writeCbModel function
%
% Authors:
%     - Laurent Heirendt

global CBTDIR

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testWriteCbModel'));
cd(fileDir);

% Note: test is only compatible with Matlab R2015b or later
if ~verLessThan('matlab', '8.6')
    % define tolerance
    tol = 1e-6;

    % load the model
    model = getDistributedModel('ecoli_core_model.mat');

    % write the model to an xls file
    writeCbModel(model, 'xlsx', 'testData');

    % read in the xls model file
    modelIn = xls2model('testData.xlsx');

    % convert an old style model
    model = convertOldStyleModel(model);

    % test
    assert(isequal(model.lb, modelIn.lb))
    assert(isequal(model.b, modelIn.b))
    assert(isequal(model.ub, modelIn.ub))
    assert(isequal(sort(model.mets), sort(modelIn.mets)))
    assert(isequal(model.csense, modelIn.csense))
    assert(isequal(model.osenseStr, modelIn.osenseStr))
    assert(isequal(model.rxns, modelIn.rxns))
    assert(isequal(sort(model.genes), sort(modelIn.genes)))
    assert(isequal(model.c, modelIn.c))

    % NOTE: model.rules and model.S are different from modelIn.S and modelIn.rules
    %       as the metabolites are not ordered in the same way.

    solverOK = changeCobraSolver('glpk');

    if solverOK
        % run an LP and compare the solutions
        solModel = optimizeCbModel(model);
        solModelIn = optimizeCbModel(modelIn);

        assert(abs(solModel.f - solModelIn.f) < tol)
        assert(solModel.stat == solModelIn.stat)
    end

    % remove the generated file
    delete('testData.xlsx');
else
    fprintf('\ntestWriteCbModel is not compatible with this version of MATLAB. Please upgrade your version of MATLAB.\n\n');
end

% test varargin

% load the model
model = getDistributedModel('ecoli_core_model.mat');

% write out using varargin
outmodel1 = writeCbModel(model, 'format', 'mat', 'fileName', 'testModel1');
outmodel2 = writeCbModel(model, 'mat', 'testModel2.mat');

% read in the testModel and testModel2
testModel1 = readCbModel('testModel1.mat');
testModel2 = readCbModel('testModel2.mat');
testModel2.description = '';
testModel1.description = '';
assert(isequal(testModel1, testModel2));

% test the legacy signature with more input arguments
[compSymbols, compNames] = getDefaultCompartmentSymbols();
outmodel3 = writeCbModel(model, 'mat', 'testModel3.mat', compNames, compSymbols);
testModel3 = readCbModel('testModel3.mat');
testModel3.description = '';
assert(isequal(testModel1, testModel3));

% test new signature
outmodel4 = writeCbModel(model, 'format', 'mat', 'fileName', 'testModel4.mat', 'compNames', compNames, 'compSymbols', compSymbols);
outmodel5 = writeCbModel(model, 'format', 'mat', 'compNames', compNames, 'compSymbols', compSymbols,'fileName', 'testModel5.mat');
testModel4 = readCbModel('testModel4.mat');
testModel5 = readCbModel('testModel5.mat');
testModel4.description = '';
testModel5.description = '';
assert(isequal(testModel4, testModel5));

% remove generate files during testing
for i = 1:5
    delete(['testModel', num2str(i),'.mat']);
end

% change to old directory
cd(currentDir);

author: opencobra

repo: cobratoolbox


text:
function obj = coeftq1_hand_eye_func_new(in1)
coef_J11 = in1(:,11);
coef_J21 = in1(:,21);
coef_J28 = in1(:,28);
coef_J33 = in1(:,33);
coef_J43 = in1(:,43);
coef_J50 = in1(:,50);
coef_J55 = in1(:,55);
coef_J62 = in1(:,62);
coef_J67 = in1(:,67);
coef_J72 = in1(:,72);
coef_J79 = in1(:,79);
obj = [coef_J11,coef_J21,coef_J28,coef_J33,coef_J43,coef_J50,coef_J55,coef_J62,coef_J67,coef_J72,coef_J79];

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
addpath(genpath('/dartfs-hpc/rc/home/m/f0042vm/software/fdaM'))

% create a basis set

dt = 0.0288;
order = 3;
l = 25;
degree = 8;

basis = create_bspline_basis([0,l/dt], degree+4, order);
bf = full(eval_basis((1:l/dt),basis));
bf = bf(:,3:end-2);

obf = spm_orth(bf); % orthogonalize bf

figure(1);
subplot(3,1,1);
cla
h1 = plot(bf,'r');
title('Spline Basis Set');
subplot(3,1,2);
h2 = plot(obf,'b');
title('SPM orthogonalized splines');

% compute inverse orthogonal transform

iO = (obf'*obf)\obf'*bf;
iO(iO < 0) = 0;

% generate a random walk;

rand_walk = sum(triu(repmat(randn(length(bf),1),1,length(bf))));

figure(1);
subplot(3,1,3);
cla
h1=plot(rand_walk,'bla');
hold on;

% fit both basis fnuctions to the random walk

B = (bf'*bf)\bf'*rand_walk';
orth_B = (obf'*obf)\obf'*rand_walk';

h2=plot(obf*orth_B,'*b');
h3=plot(bf*B,'r');
h4=plot(bf.*repmat(B',size(bf,1),1),'g')

legend([h1(1),h2(1),h3(1),h4(1)],'Random Walk','spm spline model','spline model','weighted splines','location','southeast');
title('Spline model can be obtained with either basis set');

% so now we have the following
% splineModel = bf*B
% splineModel = obf*orth_B
% obf*iO = bf
% so we can rearrange
% obf*iO*b = obf*orth_B
% iO*b = orth_B
% b = inv(iO)*orth_B

i0 = (bf'*obf)'\obf'*obf;
iO(iO<0) = 0;

sum(inv(iO)*orth_B)
sum(B)

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [feasible,feaslistLMI] = isfeasible(F,tol)

if nargin == 1
    tol = 0;
end
feaslistLMI = check(F);
feasible = all(feaslistLMI >= -tol);

author: yalmip

repo: YALMIP


text:
function x = nsctrec(y, dfilt, pfilt)
% NSCTREC   Nonsubsampled Contourlet Reconstruction
%
%	x = nsscrec(y, [dfilt, pfilt] )
%
% INPUT:
%   y:  a cell vector of length length(nlevs) + 1, where except y{1} is
%       the lowpass subband, each cell corresponds to one pyramidal
%       level and is a cell vector that contains bandpass directional
%       subbands from the DFB at that level.
%   dfilt:
%       a string, filter name for the directional decomposition step.
%       It is optional with default value 'dmaxflat7'. See dfilters.m for all
%       available filters.
%   pfilt:
%       a string, filter name for the pyramidal decomposition step.
%       It is optional with default value 'maxflat'. See atrousfilters.m for
%       all available filters.
%
% OUTPUT:
%   x:
%       a matrix, reconstructed image
%
% See also: ATROUSFILTERS, DFILTERS, NSCTDEC, NSFBREC, NSDFBREC

%
% HISTORY:
%     02/17/04 Created by Jianping Zhou.
%     08/07/04 Modified by Jianping Zhou. Incorporte the fast implementation
%              of convolution algorithm by Jason Laska.
%     08/30/04 Modified by Arthur. L. Cunha, added more pyramid filters.
%     10/17/04 Modified by Arthur. L. Cunha, replaced periodic with symmetric extension
%     01/24/05 Modified by Jianping Zhou, changed function names and corrected a bug.
%     10/31/05 Modified by Arthur L Cunha, corrected a bug in the pyramid decomposition (a wrong index)

% Check input
if ~exist('dfilt', 'var')
    dfilt = 'dmaxflat7' ;
end;

if ~exist('pfilt', 'var')
    pfilt = 'maxflat' ;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get fan filters, parallelogram filters, and pyramid filters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get the directional filters for the critically sampled DFB.
filters = cell(4) ;
[h1, h2] = dfilters(dfilt, 'r');
% A scale is required for the nonsubsampled case.
h1 = h1./sqrt(2) ;
h2 = h2./sqrt(2) ;

% Generate the first-level fan filters by modulations.
filters{1} = modulate2(h1, 'c');
filters{2} = modulate2(h2, 'c');

% Obtain the parallelogram filters from the diamond filters
[filters{3}, filters{4}] = parafilters( h1, h2 ) ;

% Currently only one filter by Arthur Cunha
% It has been normalized.
[h1, h2, g1, g2] = atrousfilters(pfilt) ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Nonsubsampled Contourlet transform with tree structure filter banks
% Nonsubsampled pyramids make multiresolution decomposition.
% Nonsubsampled directional filter banks make directional decomposition.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n = length(y) - 1;
xlo = y{1};
nIndex = n-1;
for i=1:n

    % Process the detail subbands
    if iscell( y{i+1} )
        % Nonsubsampled DFB reconstruction
        xhi = nsdfbrec( y{i+1}, filters );
    else
        % No DFB decomposition, copy directly
        xhi = y{i+1};
    end

    % Nonsubsampled Pyramid reconstruction
    x = nsfbrec(xlo, xhi, g1, g2, nIndex);

    % Prepare for the next level
    xlo = x ;
    nIndex = nIndex -1;
end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function title = p01_title ( )

%*****************************************************************************80
%
%% P01_TITLE returns a title for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, a title for the problem.
%
  title = 'The constant function';

  return
end


text:
 function [w0, phi] = adw_fan_new(sg, ig, wi)
%function [w0, phi] = adw_fan_new(sg, ig, wi)
%
% Compute the angular-dependent weighting for fan-beam geometry.
% w0(\Phi) = w(x0, y0, \Phi)
% For fully corrected penalty:
%	w(s',\beta') * J(s') | \phi'=\Phi + ...
%	w(s',\beta') * J(s') | \phi'=\Phi-pi
% See fessler chapter 3.5.3 (?).
% Useful for regularization design and variance prediction.
%
% in:
%	sg	sino_geom()
%	ig	image_geom()
%	wi	[nb,na] var(yi)
% out:
%	w0	[np,na] angular-dependent weighting, np = sum(ig.mask(:))
%
% Copyright 2005-4-07, Yingying Zhang & Jeff Fessler, The University of Michigan

if nargin == 1 && streq(sg, 'test'), adw_fan_test, return, end
if nargin < 3, ir_usage, end

[x y] = ndgrid(ig.x, ig.y); % image domain pixel locations

r0 = sqrt(x.^2 + y.^2);
phi0 = atan2(y,x); % use it for our purpose since atan will have pi flip
% [phi0 r0] = cart2pol(x,y);
% image domain dummy variables for polar coordinates

% dummy variables in freq. domain
phi = sg.ar;

% sinogram domain locations
beta = sg.ar; % angles in radians
gamma = sg.s / sg.dsd;

% loop over \phi_i
Dcf = sg.dsd + sg.dfs;
ratio_fcf = sg.dfs / Dcf;
% w0 = zeros(ig.nx, ig.ny, sg.na);
w0 = zeros([size(r0(ig.mask),1) sg.na]);

ticker reset
for ia = 1:sg.na
	ticker(mfilename, ia, sg.na)

	% -------------------------
	% for \phi' = \phi
	% -------------------------

	% compute new variables in (3.5.2)
	r_p = r0 .* cos(phi(ia) - phi0); % use beta, not outer_sum(gamma, beta)
%	since evaluate at \phi_p = \phi(view angle)
%	r_p = x .* cos(phi(ia)) + y .* sin(phi(ia));
	s_p = Dcf .* (asin(r_p / sg.dso) - asin(ratio_fcf .* r_p / sg.dso));
	beta_p = phi(ia) - asin(r_p / sg.dso);

	% compute Jacobian determinant (2.7.9)
	Jacob_p = abs(sg.dso * cos(s_p / sg.dsd) - ...
	sg.offset * sin(s_p / sg.dsd)) / sg.dsd;

	% find corresponding wi: nearest neighbor method
	s_cen = sg.w;
	% caution: this looks wrong!!! (sg.nb + 1)/2 + ob.offset_s; % for nb = even % <-----(yy) plus:(to me account for the right index) or (hugo) minus?
	s_loc = round(s_p / sg.d + s_cen);
	s_loc = min(s_loc, sg.nb);
	s_loc = max(s_loc, 1);
	beta_loc = round((mod(beta_p,2*pi)-deg2rad(sg.orbit_start)) / (2*pi/sg.na)) + 1;
	beta_loc = min(beta_loc, sg.na);
	beta_loc = max(beta_loc, 1);

	% find the corresponding index in column-stack wi
	loc = 1 + (s_loc(:)-1) + (beta_loc(:)-1)*sg.nb;

	% compute w0 at \phi' = \phi
%	w0(:,:,ia) = w0(:,:,ia) + reshape(wi(loc(:)), [ig.nx ig.ny]) .* (1 ./ Jacob_p);
	wi_temp = wi(loc(:));
	wi_p = wi_temp(ig.mask);
	Jacob = 1 ./ Jacob_p(ig.mask);
	w0(:,ia) = wi_p .* Jacob;

	clear s_loc beta_loc loc

	% -------------------------
	% for \phi' = \phi - pi
	% -------------------------
	phi_pi = phi(ia) - pi;
	r_p = r0 .* cos(phi_pi - phi0);
%	r_p = x .* cos(phi(ia)) + y .* sin(phi(ia));
	s_p = Dcf .* (asin(r_p / sg.dso) - asin(ratio_fcf .* r_p / sg.dso));
	beta_p = phi_pi - asin(r_p / sg.dso);

	% Jacob_p again: not needed since same
%	Jacob_p = abs(sg.dso*cos(s_p / sg.dsd) - ...
%	sg.offset * sin(s_p / sg.dsd)) / sg.dsd;

	% find corresponding wi: nearest neighbor method
	s_loc = round(s_p / sg.d + s_cen);
	s_loc = min(s_loc, sg.nb);
	s_loc = max(s_loc, 1);

	beta_loc = round((mod(beta_p,2*pi)-deg2rad(sg.orbit_start)) / (2*pi/sg.na)) + 1;
	beta_loc = min(beta_loc, sg.na);
	beta_loc = max(beta_loc, 1);

	% find the corresponding index in column-stack wi
	loc = 1 + (s_loc(:)-1) + (beta_loc(:)-1)*sg.nb;

	% compute w0 at \phi' = \phi - pi and sum with \phi' = \phi
%	w0(:,:,ia) = w0(:,:,ia) + reshape(wi(loc(:)), [ig.nx ig.ny]) .* (1 ./ Jacob_p);
	wi_temp = wi(loc(:));
	wi_p = wi_temp(ig.mask);
	w0(:,ia) = w0(:,ia) + wi_p .* Jacob;

	clear s_loc beta_loc loc

	% angle_dependent blur
%	b0(:,:,ia) = (sg.ds / Dcf) .* sqrt(sg.dso^2 + r0^2 + 2 .* sg.dso .* Dcf .* cos(phi(ia)-phi0));
end

function adw_fan_test
down = 8;
sg = sino_geom('ge1', 'down', down);
ig = image_geom('nx', 512, 'fov', 500, 'down', down);
wi = sg.ones;
w0 = adw_fan_new(sg, ig, wi);
w0 = ig.embed(w0);
im clf, im(w0(:,:,1:10:end)), cbar

author: JeffFessler

repo: mirt


text:
function [ y, j, f, ierror ] = yjf_check_english ( y, j, f )

%*****************************************************************************80
%
%% YJF_CHECK_ENGLISH normalizes an English YJF date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, J, real F, the YJF date.
%
%    Output, integer IERROR, nonzero if there was an error.
%
  [ y, j, ierror ] = yj_check_english ( y, j );

  if ( ierror ~= 0 )
    return
  end
%
%  Force the fraction to lie between 0 and 1.
%
  while ( f < 0.0 )
    f = f + 1.0
    j = j - 1;
  end

  while ( 1.0 <= f )
    f = f - 1.0;
    j = j + 1;
  end

  return
end


text:
function rvec = s_to_r8vec ( s, n )

%*****************************************************************************80
%
%% S_TO_R8VEC reads an R8VEC from a string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be read.
%
%    Input, integer N, the number of values expected.
%
%    Output, real RVEC(N), the values read from the string.
%
  i = 0;
  rvec = zeros ( n, 1 );

  ilo = 1;

  while ( i < n )

    i = i + 1;

    [ rvec(i), lchar ] = s_to_r8 ( s(ilo:end) );

    ilo = ilo + lchar;

  end

  return
end


text:
  function x = flipdims(x, varargin)
%|function x = flipdims(x, varargin)
%|
%| Generalization of flipdim() that flips all dimensions.
%| This is useful for the adjoint of convolution operator.
%|
%| option
%|	'odd'	0|1	if 1, make each dimension odd sized after flipping
%|			by appending a zero.  default: 0
%|
%| Jeff Fessler

if nargin < 1, ir_usage, end

arg.odd = 0;
arg = vararg_pair(arg, varargin);

for ii = 1:ndims(x)

	if arg.odd && ~rem(size(x,ii), 2)
		sz = size(x);
		sz(ii) = 1;
		try
			tmp = zeros(sz, 'like', x);
		catch % for octave
			tmp = zeros(sz, class(x));
		end
		x = cat(ii, x, tmp); % append one zero along dim ii
	end

	x = flip(x, ii);
end


text:
function title = p03_title ( )

%*****************************************************************************80
%
%% P03_TITLE returns the title of problem 3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = 'The pulse data, 0 0 0 0 0 1 0 0 0 0 0';

  return
end


text:
% Calculate feature vector for a triangle
%

% Feature vector is:
%   centroid_x
%   centroid_y
%   orientation_1
%   orientation_2
%   orientation_3
%   area
function features = triangleFeatures(points)
    features = zeros(6, 1);

    fCentrX = 1;
    fCentrY = 2;
    fOr1 = 3;
    fOr2 = 4;
    fOr3 = 5;
    fArea = 6;

    centr = centroid(points);
    features(fCentrX) = centr(1);
    features(fCentrY) = centr(2);

    line = createLine(centr, points(1, :));
    features(fOr1) = rad2deg(lineAngle(line));

    line = createLine(centr, points(2, :));
    features(fOr2) = rad2deg(lineAngle(line));

    line = createLine(centr, points(3, :));
    features(fOr3) = rad2deg(lineAngle(line));

%     features(fArea) = triangleArea(points(1, :), points(2, :), points(3, :));

    features(6) = points(1, 1);
    features(7) = points(1, 2);
    features(8) = points(2, 1);
    features(9) = points(2, 2);
    features(10) = points(3, 1);
    features(11) = points(3, 2);
end

author: kavli-ntnu