text:
%% LineSrc
% Concrete subclass of <Source.html |Source|> representing a constant electric
% dipole distribution along a line.

%%% Description
% |LineSrc| is to generate a cylindrical wave in a homogeneous medium.  It can
% creates a cylindrical wave propagating not only in the radial direction, but
% also in the direction slanted toward the line axis.  It also supports the
% polarization normal to the line axis.

%%% Construction
%  src = LineSrc(axis, intercept)
%  src = LineSrc(axis, intercept, polarization)
%  src = LineSrc(axis, intercept, polarization, IorK)
%  src = LineSrc(axis, intercept, polarization, IorK, ka_Bloch)
%  src = LineSrc(axis, intercept, polarization, IorK, theta, wvlen)
%
% *Input Arguments*
%
% * |axis|: direction of the line.  It should be one of |Axis.x|, |Axis.y|,
% |Axis.z|.
% * |intercept|: location of the line in the plane normal to the line.  For
% |axis == Axis.y|, it is in the format of |[z x]|.
% * |polarization|: direction of the dipoles distributed on the plane.  It
% should be one of |Axis.x|, |Axis.y|, |Axis.z|.  If unassigned, |polarization =
% axis| is used.
% * |IorK|: for |polarization == axis|, complex amplitude of the current that
% the dipoles drive; for |polarization ~= axis|, complex amplitude of the sheet
% current density that the dipoles drive.  If unassigned, |IorK = 1.0| is used.
% * |ka_Bloch|: Bloch wavevector in the axis of this line.  If unassigned,
% |ka_Bloch = 0| is used.
% * |theta|: oblique incidence angle.  |abs(theta)| should not exceed |pi/2|. If
% unassigned, |theta = 0.0| is used.
% * |wvlen|: wavelength of the cylindrical wave in the background medium.  It is
% required to set up the Bloch boundary condition for oblique incidence. |wvlen|
% is not the vacuum wavelength used in the frequency-domain Maxwell's equations,
% but the wavelength in the medium where this |LineSrc| is located.

%%% Note
% In the finite-difference grid, the dipoles in |LineSrc| are located at the
% _E_-field points.  This poses a condition on |intercept| argument in the
% constructor: if |polarization == axis|, both elements of |intercept| should be
% at dual grid points; if |polarization ~= axis|, the element of |intercept|
% corresponding to the |polarization| component should be at a primary grid
% point, and the other element of |intercept| should be at a dual grid point.
% Therefore, make sure that the elements of |intercept| at primary grid points
% do not overlap with the locations of the vertices of <Shape.html |Shape|>;
% otherwise dynamic grid generation in <moxwell_run.html |maxwell_run|> will
% fail.

%%% Example
%   % Create an instance of LineSrc.
%   src = LineSrc(Axis.z, [0 0]);
%
%   % Use the constructed src in maxwell_run().
%   [E, H] = maxwell_run({INITIAL ARGUMENTS}, 'SRCJ', src);

%%% See Also
% <PointSrc.html |PointSrc|>, <PlaneSrc.html |PlaneSrc|>, <maxwell_run.html
% |maxwell_run|>

classdef LineSrc < Source & WithBloch

	properties (SetAccess = immutable)
		axis  % axis of line: one of Axis.x, Axis.y, Axis.z
		intercept  % [h, v]: location of intercept in transverse plane
		polarization  % one of Axis.x, Axis.y, Axis.z
		kBloch  % Bloch k-vector [kx, ky, kz]
		IorK  % current I if polarization == axis; sheet current density K otherwise
	end

	methods
		function this = LineSrc(axis, intercept, polarization, IorK, theta, wvlen)
			chkarg(istypesizeof(axis, 'Axis'), '"axis" should be instance of Axis.');
			chkarg(istypesizeof(intercept, 'real', [1 Dir.count]), ...
				'"intercept" should be length-%d row vector with real elements.', Dir.count);

			if nargin < 3  % no polarization
				polarization = axis;
			end
			chkarg(istypesizeof(polarization, 'Axis'), '"polarization" should be instance of Axis.');

			if nargin < 4  % no IorK
				IorK = 1.0;
			end
			chkarg(istypesizeof(IorK, 'complex'), '"IorK" should be complex.');

			if nargin < 5  % no theta (or ka_Bloch)
				ka_Bloch = 0;
			elseif nargin < 6  % no wvlen
				ka_Bloch = theta;  % theta is in fact ka_Bloch
				chkarg(istypesizeof(ka_Bloch, 'real'), '"ka_Bloch" should be real.');
			else  % with wvlen
				assert(nargin == 6);
				chkarg(istypesizeof(theta, 'real') && theta >= -pi/2 && theta <= pi/2, ...
					'"theta" should be polar angle in radian between -pi/2 and pi/2.');
				chkarg(istypesizeof(wvlen, 'real') && wvlen > 0, '"wvlen" should be positive.');
				ka_Bloch = (2*pi/wvlen) * sin(theta);
			end

			lgrid = cell(1, Axis.count);
			laltgrid = cell(1, Axis.count);
			[h, v] = cycle(axis);  % h, v ~= axis
			if polarization == h
				laltgrid{h} = intercept(Dir.h);
				lgrid{v} = intercept(Dir.v);
			elseif polarization == v
				lgrid{h} = intercept(Dir.h);
				laltgrid{v} = intercept(Dir.v);
			else  % polarization == axis
				lgrid{h} = intercept(Dir.h);
				lgrid{v} = intercept(Dir.v);
			end

			line = Line(axis, intercept);
			this = this@Source(lgrid, laltgrid, line);

			this.axis = axis;
			this.intercept = intercept;
			this.polarization = polarization;
			this.IorK = IorK;

			this.kBloch = [0 0 0];
			this.kBloch(this.axis) = ka_Bloch;
		end

		function [index_cell, JMw_patch] = generate_kernel(this, w_axis, grid3d)
			grid3d.set_kBloch(this);
			[h, v, a] = cycle(this.axis);  % h, v: axes in plane normal to this line
			p = this.polarization;
			if w_axis ~= p
				JMw_patch = [];
				index_cell = cell(1, Axis.count);
			else  % w_axis == p
				index_cell = {':', ':', ':'};
				if p == a
					% Set index_cell{h} and index_cell{v}.
					g = this.gt;
					axes = [h v];
					for d = Dir.elems
						ind = ind_for_loc(this.intercept(d), axes(d), g, grid3d);
						index_cell{axes(d)} = ind;
					end

					% Set Jw.
					dh = grid3d.dl{h,g}(index_cell{h});
					dv = grid3d.dl{v,g}(index_cell{v});
					I = this.IorK;
					JMw = I / (dh * dv);  % I = J * (area)
					ga = alter(this.gt);
				else  % p ~= a
					% Set index_cell{h}.
					if h == p
						g = alter(this.gt);
					else
						g = this.gt;
					end
					ind = ind_for_loc(this.intercept(Dir.h), h, g, grid3d);
					index_cell{h} = ind;

					% Set index_cell{v}
					g = alter(g);
					ind = ind_for_loc(this.intercept(Dir.v), v, g, grid3d);
					index_cell{v} = ind;

					% Set Jw.
					q = setdiff([h v], p);
					dq = grid3d.dl{q, this.gt}(index_cell{q});
					K = this.IorK;
					JMw = K / dq;
					ga = this.gt;
				end

				% Set Jw_patch.
				la = grid3d.l{a, ga};
				ka = this.kBloch(a);

				JMw_patch = JMw .* exp(-1i * (ka .* la)).';
				JMw_patch = ipermute(JMw_patch, int([a h v]));
			end
		end
	end
end

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function [B] = bit2B(bit)
% Convert computery things from bits to bytes.
% Chad A. Greene 2012
B = bit*2^-3;

author: Sable

repo: mcbench-benchmarks


text:
function p=v_paramsetch(d,q,m,c,t)
%V_PARAMSETCH update and check parameter values p=(d,q,m,c,t)
% Usage: (1) function x=func(y,q)
%            d=struct('a',1,'b',2,'c',3); % default parameters
%            p=v_paramsetch(d,q); % update selected parameters
%
%        (2) function x=func(y,q)
%            d=struct('a',1,'b',2,'c',3); % default parameters
%            c={'p.a>0 && p.a<5','p.b>0'};
%            p=v_paramsetch(d,q,'E',c); % check parameter ranges
%
%        (3) t={'a','description of parameter a';'c','and of parameter c'}
%            p=v_paramsetch(d,q,'l',c,t); % list values with optional descritions
%                                       % '-','*','+' indicates default, updated and new fields
%
%  Inputs:
%       d  default parameter structure
%       q  new parameter values either a struct or alternatively matrix with
%          each row a different variable in the same order as the fields of d
%       m  mode string: any combination of the following
%           'a' ignore additional fields in q
%           'A' additional fields in q constitute an error
%           'e' print errors but don't exit
%           'E' print errors and exit
%           'l' list fields and their values (default if no output)
%       c  cell array with parameter checking conditions e.g. 'p.a>3' (use p for structure name)
%       t  cell array with descriptive text for each field in a new row. Either in
%          the form t(:,*)={'field' 'description'} or a single column of
%          descriptions in the same order as the fields of d
%
% Outputs:
%       p  output parameter structure
%

%      Copyright (C) Mike Brookes 2017
%      Version: $Id: v_paramsetch.m 10865 2018-09-21 17:22:45Z dmb $
%
%   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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p=d;                % initialize to the default values
numerr=0;           % initialize error count
dn=fieldnames(d);   % default list of parameter fields
ndn=length(dn);     % number of default parameter fields
dup=zeros(ndn,1);   % update flags for default fields
% sort out input arguments
if nargin<5
    t={'' ''};      % define an empty description array
    if nargin<4
        c=cell(0);  % define an empty check condition array
        if nargin<3
            m='';   % define an empty mode string
        end
    end
end
% now update the selected fields
if nargin>1 && numel(q)>0   % if update argument exists
    if isstruct(q)          % if update argument is a structure
        qn=fieldnames(q);   % field names to update
        addnew=~any(lower(m)=='a'); % new fields should be added into p
        adderr=any(m=='A');         % new fields constitute an error
        for i=1:length(qn)
            fi=qn{i};
            old=isfield(p,fi); % is this an existing field ?
            if addnew || old
                p.(fi)=q.(fi);
            end
            if old
                dup(find(strcmp(fi,dn),1))=1; % indicate field is updated
            end
            if adderr && ~old
                fprintf(2,'%s is an unknown parameter field\n',fi);
                numerr=numerr+1; % increment error count
            end
        end
    else                            % else update argument is a matrix
        nq=min(size(q,1),ndn);      % number of fields to update
        dup(1:nq)=1;                % indicate which fields are updated
        for i=1:nq
            p.(dn{i})=q(i,:);
        end
        if size(q,1)>nq && any(lower(m)=='e')
            fprintf(2,'More than %d parameters specified\n',nq);
            numerr=numerr+1;
        end
    end
end
% Apply parameter checks
if any(lower(m)=='e') && numel(c)>0
    for i=1:numel(c)
        if any(~eval(c{i}))
            numerr=numerr+1;
            fprintf(2,'Parameter check failed: %s\n',c{i});
        end
    end
end
% print out a list of the parameters if requested
if ~nargout || any(m=='l')
    pn=fieldnames(p);
    nf=length(pn);
    st=size(t);
    for i=1:nf
        fi=pn{i};
        vi=p.(fi);
        if i>ndn
            cat='+';
        else
            cat='-'+('*'-'-')*dup(i);
        end
        if st(2)>1
            jti=find(strcmp(fi,t(:,1)),1);
            if ~isempty(jti)
                jti=t{jti,2}; % description string
            end
        elseif i<=st(1)
            jti=t{i,1}; % description string
        else
            jti=[];
        end
        if isnumeric(vi) && length(vi)==numel(vi) && isreal(vi) % can print on one line
            fit=fi;
            if size(vi,1)>1
                fit=[fi ''''];
            end
            fprintf('%3d%c %s =',i,cat,fit);
            fprintf(' %g',vi);
            if isempty(jti)
                fprintf('\n');
            else
                fprintf(' = %s\n',jti);
            end
        else
            fprintf('%3d%c %s =',i,cat,fi);
            if isempty(jti)
                fprintf('\n');
            else
                fprintf(' %s =\n',jti);
            end
            disp(vi);
        end
    end
end
if numerr>0 && any(m=='E')
    error('%d error%c in parameter specification',numerr,(' '+(numerr>1)*('s'-' ')));
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
%DEMO_QUANTILEGP Demonstration of Quantile GP regression
%
%  Description
%    This demo demonstrates Quantile-GP regression with toy data. QGP can
%    be used for estimating quantiles of interest of the response variable
%    with respect to input variables.
%
%   Reference
%     Boukouvalas et al. (2012). Direct Gaussian Process Quantile
%     Regression Using Expectation Propagation. In Proceedings of the
%     29th International Conference on Machine Learning, Edinburgh,
%     Scotland, UK, 2012.
%
%  See also  LIK_QGP, DEMO_*
%
% Copyright (c) 2013 Ville Tolvanen

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

prevstream=setrandstream(0);
% Toy data
n=200;
x=linspace(-3,2,n)';
y=((x.^3 + 3.*x.^2-6.*x-8)/4);
x=(x-mean(x))/std(x);
y=(y-mean(y))/std(y);
xt=x;
yt=y;

% add heteroscedastic noise
noise=0.2.*randn(n,1);
noise=noise + 0.2.*exp(-1.5.*x).*randn(n,1);
y=y+noise;

% Create covariance function
gpcf=gpcf_sexp('lengthScale', 1, 'magnSigma2', 1);

% Create quantile-GP likelihood
lik=lik_qgp('sigma2', 0.5, 'sigma2_prior', prior_t('s2',10,'nu',1));

% We want to infer the 0.25, 0.5 and 0.75 quantiles
lik=lik_qgp(lik, 'quantile', 0.05);
lik2=lik_qgp(lik, 'quantile', 0.5);
lik3=lik_qgp(lik, 'quantile', 0.95);

lik4=lik_qgp(lik, 'quantile', 0.25);
lik5=lik_qgp(lik, 'quantile', 0.75);

% Create gp structure. Quantile-GP only works with EP or MCMC. Lets use
% MCMC for 5%, 50% and 95% quantiles.
gp = gp_set('lik', lik, 'cf', gpcf, 'jitterSigma2', 1e-6, 'latent_method', 'MCMC');
gp2 = gp_set('lik', lik2, 'cf', gpcf, 'jitterSigma2', 1e-6, 'latent_method', 'MCMC');
gp3 = gp_set('lik', lik3, 'cf', gpcf, 'jitterSigma2', 1e-6, 'latent_method', 'MCMC');

% EP for 25% and 75% quantiles (with low number of data points, EP is
% unstable for very low or high quantiles)
gp4 = gp_set('lik', lik4, 'cf', gpcf, 'jitterSigma2', 1e-6, 'latent_method', 'EP');
gp5 = gp_set('lik', lik5, 'cf', gpcf, 'jitterSigma2', 1e-6, 'latent_method', 'EP');

% Set options for optimization
opt=optimset('TolX',1e-4, 'TolFun', 1e-4, 'Display', 'iter', 'derivativecheck', 'off');

% Sample values
ssls_opt.latent_opt.repeat=100;
rgp=gp_mc(gp, x, y, 'nsamples',200,'ssls_opt', ssls_opt, 'display', 10);
rgp=thin(rgp,10);
rgp2=gp_mc(gp2, x, y, 'nsamples', 200, 'ssls_opt', ssls_opt, 'display', 10);
rgp2=thin(rgp2,10);
rgp3=gp_mc(gp3, x, y, 'nsamples', 200, 'ssls_opt', ssls_opt, 'display', 10);
rgp3=thin(rgp3,10);

% Optimize hyperparameters
gp4=gp_optim(gp4, x, y, 'opt', opt, 'optimf', @fminlbfgs);
gp5=gp_optim(gp5, x, y, 'opt', opt, 'optimf', @fminlbfgs);

% Do predictions for training points
[Ef1,Varf1]=gp_pred(rgp, x, y, xt);
[Ef2,Varf2]=gp_pred(rgp2, x, y, xt);
[Ef3,Varf3]=gp_pred(rgp3, x, y, xt);
[Ef4,Varf4]=gp_pred(gp4, x, y, xt);
[Ef5,Varf5]=gp_pred(gp5, x, y, xt);

% Visualize
figure(1);
plot(xt, Ef1, '-r',xt, Ef4,'-m', xt, Ef2, '-k', xt, Ef3, '-r', xt, y, '.',xt, Ef5,'-m');
legend('5% / 95%', '25% / 75%', '50%');
setrandstream(prevstream);

author: gpstuff-dev

repo: gpstuff

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


text:
function [K0,K1,K2,H1] = spm_kernels(varargin)
% returns global Volterra kernels for a MIMO Bilinear system
% FORMAT [K0,K1,K2] = spm_kernels(M,P,N,dt)            - output kernels
% FORMAT [K0,K1,K2] = spm_kernels(M0,M1,N,dt)          - state  kernels
% FORMAT [K0,K1,K2] = spm_kernels(M0,M1,L1,N,dt)       - output kernels (1st)
% FORMAT [K0,K1,K2] = spm_kernels(M0,M1,L1,L2,N,dt)    - output kernels (2nd)
%
% M,P   - model structure and parameters;
%         or its bilinear reduction:
%
% M0    - (n x n)     df(q(0),0)/dq                    - n states
% M1    - {m}(n x n)  d2f(q(0),0)/dqdu                 - m inputs
% L1    - (l x n)     dldq                             - l outputs
% L2    - {m}(n x n)  dl2dqq
%
% N     - kernel depth       {intervals}
% dt    - interval           {seconds}
%
% Volterra kernels:
%---------------------------------------------------------------------------
% K0    - (1 x l)             = K0(t)         = y(t)
% K1    - (N x l x m)         = K1i(t,s1)     = dy(t)/dui(t - s1)
% K2    - (N x N x l x m x m) = K2ij(t,s1,s2) = d2y(t)/dui(t - s1)duj(t - s2)
%
%___________________________________________________________________________
% Returns Volterra kernels for bilinear systems of the form
%
%         dq/dt   = f(q,u) = M0*q + M1{1}*q*u1 + ... M1{m}*q*um
%            y(i) = L1(i,:)*q + q'*L2{i}*q
%
% where q = [1 x(t)] are the states augmented with a constant term
%
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_kernels.m 6937 2016-11-20 12:30:40Z karl $

% assign inputs
%--------------------------------------------------------------------------
if nargin == 4

    M0 = varargin{1};
    M1 = varargin{2};
    N  = varargin{3};
    dt = varargin{4};

elseif nargin == 5

    M0 = varargin{1};
    M1 = varargin{2};
    L1 = varargin{3};
    N  = varargin{4};
    dt = varargin{5};

elseif nargin == 6

    M0 = varargin{1};
    M1 = varargin{2};
    L1 = varargin{3};
    L2 = varargin{4};
    N  = varargin{5};
    dt = varargin{6};
end

% bilinear reduction if necessary
%--------------------------------------------------------------------------
if isstruct(M0)
    [M0,M1,L1,L2] = spm_bireduce(M0,M1);
end

% Volterra kernels for bilinear systems
%==========================================================================

% make states the outputs (i.e. remove constant) if L1 is not specified
%--------------------------------------------------------------------------
try
    L1;
catch
    L1 = speye(size(M0));
    L1 = L1(2:end,:);
end
try
    L2;
catch
    L2 = [];
end

% parameters
%--------------------------------------------------------------------------
N     = fix(N);                     % kernel depth
n     = size(M0,1);                 % state variables
m     = size(M1,1);                 % inputs
l     = size(L1,1);                 % outputs
H1    = zeros(N,n,m);
K1    = zeros(N,l,m);
K2    = zeros(N,N,l,m,m);
M0    = full(M0);

% pre-compute matrix exponentials
%--------------------------------------------------------------------------
e1    = expm( dt*M0);
e2    = expm(-dt*M0);
for p = 1:m
    M{1,p} = e1*M1{p}*e2;
end
ei    = e1;
for i = 2:N
    ei    = e1*ei;
    for p = 1:m
        M{i,p} = e1*M{i - 1,p}*e2;
    end
end

% 0th order kernel
%--------------------------------------------------------------------------
X0    = sparse(1,1,1,n,1);
if nargout > 0
    H0    = ei*X0;
    K0    = L1*H0;
end

% 1st order kernel
%--------------------------------------------------------------------------
if nargout > 1
    for p = 1:m
        for i = 1:N
            H1(i,:,p) = M{i,p}*H0;
            K1(i,:,p) = H1(i,:,p)*L1';
        end
    end
end

% 2nd order kernels
%--------------------------------------------------------------------------
if nargout > 2
    for p = 1:m
        for q = 1:m
            for j = 1:N
                H  = L1*M{j,q}*H1((j:N),:,p)';
                K2(j,j:N,:,q,p) = H';
                K2(j:N,j,:,p,q) = H';

            end
        end
    end

    if isempty(L2), return, end

    % add output nonlinearity
    %----------------------------------------------------------------------
    for i = 1:m
        for j = 1:m
            for p = 1:l
                K2(:,:,p,i,j) = K2(:,:,p,i,j) + H1(:,:,i)*L2{p}*H1(:,:,j)';
            end
        end
    end
end

author: spm


text:
function res = imGaussianFilter(img, kernelSize, sigma, varargin)
%IMGAUSSIANFILTER Apply gaussian filter to an image, using separability
%
%   IMGF = imGaussianFilter(IMG, SIZE, SIGMA)
%   Applies gaussian filtering on input image IMG with kernel size defined
%   by SIZE and SIGMA.
%   IMG is the input image,
%   SIZE is the size of the convolution kernel, either as a scalar, or as a
%   1-by-ND row vector containging size in the X, Y, and eventually z
%   direction. Default is 3.
%   SIGMA is the width of the kernel, either as a scalar (the same sigma
%   will be used in each direction), or as a row vector containing sigmax,
%   sigmay, and eventually sigmaz. Default is SIZE/2.
%
%   IMGF = imGaussianFilter(IMG, SIZE, SIGMA, OPTIONS)
%   Apply the same kind of options than for imfilter.
%
%   The function works for 2D or 3D images, for grayscale or color images.
%   In case of color images, the filtering is repeated for each channel of
%   the image.
%
%
%   Example
%     % Gaussian filter of a grayscale image
%     img = imread('cameraman.tif');
%     imgf = imGaussianFilter(img, 11, 4);
%     % is equivalent, but is in general faster, that:
%     imgf2 = imfilter(img, fspecial('gaussian', 11, 4));
%
%     % Using anisotropic filtering
%     img = imread('cameraman.tif');
%     imgf = imGaussianFilter(img, [13 5], [4 2]);
%     figure; subplot(121); imshow(img); subplot(122); imshow(imgf);
%
%     % Gaussian filtering of a color image
%     img = imread('peppers.png');
%     imgf = imGaussianFilter(img, [5 5], [2 2]);
%     imshow(imgf)
%
%   Note that there can be slight differences due to rounding effects. To
%   minimize them, it is possible to use something like:
%   imgf3 = uint8(imGaussianFilter(single(img), 11, 4));
%
%   See also
%   gaussianKernel3d, imfilter, fspecial
%
%   Requires
%   imfilter in the image processing toolbox

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

%   HISTORY
%   2017.08.02 use replicate option as default

%% Dispatch processing in case of color images

% compute image dimension
nd = ndims(img);

% case of color images

if isColorImage(img)
    % allocate result for memory by duplicating original image
    res = img;
    if nd == 3
        for i = 1:3
            % filter each channel of the 2D image
            res(:,:,i) = imGaussianFilter(img(:,:,i), ...
                kernelSize, sigma, varargin{:});
        end
    else
        for i = 1:3
            % filter each channel of the 3D image
            res(:,:,i,:) = imGaussianFilter(img(:,:,i,:), ...
                kernelSize, sigma, varargin{:});
        end
    end
    return;

end

%% Process input arguments

% process kernel size
if nargin < 2
    kernelSize = 3;
end
if length(kernelSize) == 1
    kernelSize = repmat(kernelSize, 1, nd);
end

% process filter sigma
if nargin < 3
    sigma = kernelSize / 2;
end
if length(sigma) == 1
    sigma = repmat(sigma, 1, nd);
end

% default options
if isempty(varargin)
    varargin = {'replicate'};
end

%% Main processing

% init result
res = img;

% apply a 1D filter in each direction
for i = 1:nd
    % compute spatial reference
    refSize = (kernelSize(i) - 1) / 2;
    s0 = floor(refSize);
    s1 = ceil(refSize);
    lx = -s0:s1;

    % compute normalized kernel
    sigma2 = 2 * sigma(i) .^ 2;
    h = exp(-(lx.^2 / sigma2));
    h = h / sum(h);

    % reshape
    newDim = [ones(1, i-1) kernelSize(i) ones(1, nd-i)];
    newDim = newDim([2 1 3:nd]);

    h = reshape(h, newDim);

    % apply filtering along one direction
    res = imfilter(res, h, varargin{:});
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function out = loggamma(x, UBorLB)
%natural logarithm of the gamma function
%UBorLB = 'ub' or 'lb', upper or lower bound

%
%   Created in 2013 by Victoria Kostina (vkostina@caltech.edu)
%

out = gamma(x);
if isfinite(out)
    out = log(out);
    return;
end
%  Sharp upper and lower bounds for the gamma function
% Horst Alzer
beta = 1/1620;
y = x - 1;
switch lower(UBorLB)
    case 'ub'
        out = .5*log(2*pi*y) + y.*log(y/exp(1)) + y/2.*log(y.*sinh(1/y)).*(1 + beta/y^5); %upper bound
    case 'lb'
        out = .5*log(2*pi*x) + y.*log(y/exp(1)) + y/2.*log(y.*sinh(1/y)); %lower bound
    otherwise
        disp('gammafun(x, UBorLB): UB or LB?');
        out = NaN;
end
end

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
function dist=dist(x,D)
n=size(x,2);
dist=0;
for i=1:n-1
    dist=dist+D(x(i),x(i+1));
end
dist=dist+D(x(1),x(n));

author: HuangCongQing


text:
function obj=MOxUnitErroredTestOutcome(test_,duration,error_)
% Errored test outcome
%
% Inputs:
%   test_               MOxUnitTestCase of the test that errored
%   duration            time in seconds it took to run the test
%   error_              struct represting the contents of lasterror() after
%                       the test failed
%
% Returns:
%   obj                 MOxUnitErroredTestOutcome instance, which is a
%                       subclass of MoxUnitTestOutcome
%
% See also: lasterror, MoxUnitTestOutcome

    s=struct();
    s.error=error_;
    obj=class(s,'MOxUnitErroredTestOutcome',...
                MOxUnitTestOutcome(test_,duration));

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function cov = ma2covv(ma,Peps,ncov)

%function cov = ma2covv(ma,Peps,ncov)
%   Transforms Moving Average parameters MA into a covariancefunction cov
%   ncov is the number of covariances
%   For a zero-mean process, the covariance function has been defined as
%   R(s)(.) = E{X[n+s]<X[n],.>},
%   or, in old-fashioned notation:
%   R(s) = E{X[n+s]X[n]'} where X[n] is a column vector.

%S. de Waele, March 2003.

s = kingsize(ma);
order = s(3)-1;
dim = s(1); I = eye(dim);

if nargin <  3, ncov = order; end
Peps_mapl = timesv(Peps,transposeLTI(ma)); % = Peps*MA'
cov = convv(ma,Peps_mapl);                 % = MA*Peps*MA'
cov = cov(:,:,1+order:end);
%Create covariance function with proper length
if ncov>order
   cov = cat(3,cov,zeros(dim,dim,ncov-order));
end
if ncov<order
	cov = cov(:,:,1+order);
end


text:
function[v,normH10_ab,scal01] = scalProds(a,kernel,rho)
% This function computes the vector which entries are
% v(i) = 2 \pi\int_{a(1)}^{a(2)}r(f'(r)ep'(r))
% where f' is the der property of kernel
% and e_p(r) = Cp(rho) J_0(rho(p) r)r
% Where Cp is a normalization constant in $H^1_0$. (See function 'Cp').
% The function uses oscillatory integrals through the matlab function
% "integral" (which is vectorized, but quite slow !).
% For big applications, we recommend using a custom Kernel object where the
% scalFunc is optimized (explicit if possible).

warning off MATLAB:integral:NonFiniteValue

if nargin==0
    run('scalProdsVal.m')
    % Unitary Test
else
    scalFunc = kernel.scalFunc;
    normFunc = kernel.normFunc;

    if length(a) == 1
        b = 1;
    else
        b = a(2);
        a = a(1);
    end
    % Vector valued function
    % (i.e. when x is scalar, fun(x) is a vector (here of size
    % [1,length(rho)])
    v = scalFunc(a,b,rho);

    if a==0
        leftEdge = 0;
    else
        if kernel.singular
            leftEdge = Inf;
        else
            leftEdge = scalFunc(0,a,rho);
        end
    end

    if b==1
        rightEdge = 0;
    else
        rightEdge = scalFunc(b,1,rho);
    end

    scal01 = (leftEdge + v + rightEdge);
    normH10_ab = normFunc(a,b);

end

v = v(:);
scal01 = scal01(:);

warning on MATLAB:integral:NonFiniteValue

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function y = lt(X,Y)
%LE (overloaded)

if isa(X,'blkvar')
    X = sdpvar(X);
end

if isa(Y,'blkvar')
    Y = sdpvar(Y);
end

try
    y = constraint(X,'<=',Y);
catch
    error(lasterr)
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% Authors: Fabrizio Smeraldi: http://www.eecs.qmul.ac.uk/~fabri
%          George Azzopardi : http://www.cs.rug.nl/~george
% Date: April 2009  -  URL  : http://www.ranklets.net
% Free for non-commercial use, all other rights reserved
%
% Efficiently Compute Intensity Ranklets using Distribution Count
% sorting
%
% See: F. Smeraldi, Fast algorithms for the computation of Ranklets,
% Proceedings of ICIP, pp 3969-3972, 2009
%
% Inputs:
%    1. imgFilename - The filename of the image to be processed
%    2. xhalfw      - The half width of the kernel; e.g. 5
%    3. yhalfw      - The half height of the kernel; e.g. 5
% Outputs:
%    1. vertical    - The response for the vertically tuned filter
%    2. horizontal  - The response for the horizontally tuned filter
%    3. diagonal    - The response for the diagonally tuned filter

function [vertical horizontal diagonal] = RankletFilter(imgFilename, xhalfw, yhalfw)

% Read image file and determine its size
if isstring(imgFilename)
    img = imread(imgFilename);
else
    img = imgFilename;
    [xsize,ysize] = size(img);
end

% Initialize the responses
horizontal = zeros(size(img));
vertical   = zeros(size(img));
diagonal   = zeros(size(img));

% Filter the image with 3 filter orientations
for x = xhalfw+1:xsize - xhalfw
  for y = yhalfw+1:ysize - yhalfw
    % [vertical(y,x) horizontal(y,x) diagonal(y,x)] = ranklet(x,y,xhalfw,yhalfw,img);
    [vertical(x,y) horizontal(x,y) diagonal(x,y)] = ranklet(x,y,xhalfw,yhalfw,img);
  end
end

% Computes the ranklet value for the given window for 3 orientations
% Inputs:
%    1. x0, y0 - The coordinate of the middle pixel in the current window
%    2. xhalfw - The half width of the kernel
%    3. yhalfw - The half height of the kernel
%    4. img    - The image
% Outputs:
%    1. vertrk  - The ranklet value for the vertically tuned filter
%    2. horizrk - The ranklet value for the horizontally tuned filter
%    3. diagrk  - The ranklet value for the diagonally tuned filter
function [vertrk, horizrk, diagrk] = ranklet(x0, y0, xhalfw, yhalfw, img)

% entire image window
window = img(x0 - xhalfw:x0 + xhalfw - 1, y0 - yhalfw:y0 + yhalfw - 1);

% compute histogram between min and max intensity values only
imin     = min(min(window));
imax     = max(max(window));
if imin == imax
    vertrk = NaN;
    horizrk = NaN;
    diagrk = NaN;
    return;
end
dynrange = double(imin:imax);
histo    = hist(window(:), dynrange);

% discard zeros in histogram and compact it
support     = histo > 0;
histo       = histo(support);
nonzerobins = sum(support);

% Midranks calculation matrix - this has got elements
%   [ 1/2  0    0    0 ...]
%   [ 1/2  1/2  0    0 ...]
%   [ 0    1/2  1/2  0 ...]
% and so on
midm = 1/2 * (diag(ones(nonzerobins,1)) + diag(ones(nonzerobins - 1,1),-1));
cumhisto = cumsum(histo)'; % Vertical vector
midranks = midm(1:nonzerobins,1:nonzerobins) * cumhisto + 0.5;

% Keep only nonzero bins for next histograms
range=dynrange(support);

% Compute T (Treatment) sets:

% Vertical
window = img(x0-xhalfw:x0-1, y0-yhalfw:y0+yhalfw-1);
vt     = hist(window(:), range);

% Horizontal
window = img(x0-xhalfw:x0+xhalfw-1, y0-yhalfw:y0-1);
ht     = hist(window(:), range);

% Diagonal
window  = img(x0-xhalfw:x0-1, y0-yhalfw:y0-1);
window2 = img (x0:x0+xhalfw-1, y0:y0+yhalfw-1);
dt      = hist(window(:), range) + hist(window2(:), range);

% Scalar product of histogram of test set
% with midranks computes Wilcoxon statistics
vertrk  = vt * midranks;
horizrk = ht * midranks;
diagrk  = dt * midranks;

% Normalize
mean    = xhalfw*yhalfw*(4*xhalfw*yhalfw+1);
scaling = 2* (xhalfw*yhalfw) * (xhalfw * yhalfw);
vertrk  = (vertrk  - mean) / scaling;
horizrk = (horizrk - mean) / scaling;
diagrk  = (diagrk  - mean) / scaling;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [Y, cardiCandi] = SparseCodingwithBasis(U, fea, options)
% SparseCodingwithBasis: Sparse Coding with the Provided Basis
%
%     [Y, cardiCandi] = SparseCodingwithBasis(U, fea, options)
%
%     Input:
%
%                 U     - basis matrix. Each column of U is a basis vector.
%               fea     - data matrix. Each column of fea is a sample.
%           options     - Struct value in Matlab. The fields in options
%                         that can be set:
%
%              ReguParaType    -  'LARs': use LARs to solve the LASSO
%                                 problem. You need to specify the
%                                 cardinality requirement in Cardi.
%
%                                 'SLEP': use SLEP to solve the LASSO
%                                 problem. You need to specify the 'ReguAlpha'
%                                 Please see http://www.public.asu.edu/~jye02/Software/SLEP/
%                                  for details on SLEP. (The Default)
%
%                      Cardi   -  An array to specify the cadinality requirement.
%                                 Default [10:10:50]
%
%                  ReguAlpha   -  regularization paramter for L1-norm regularizer
%                                 Default 0.05
%
%
%     Output:
%
%                 Y     - A cell.  matrix. each element in the cell will be
%                         a sparse representation of the input fea.
%        cardiCandi     - An array. The length of 'cardiCandi' is the same
%                         as the length of 'Y'.
%            elapse     - Time on SparseCodingwithBasis
%
%
%    Examples:
%
% See SCC
%
%Reference:
%
%   [1] Deng Cai, Hujun Bao, Xiaofei He, "Sparse Concept Coding for Visual
%   Analysis," The 24th IEEE Conference on Computer Vision and Pattern
%   Recognition (CVPR 2011), pp. 2905-2910, Colorado Springs, CO, USA,
%   20-25 June 2011.
%
%   version 2.0 --Jan/2012
%   version 1.0 --Aug/2009
%
%   Written by Deng Cai (dengcai AT gmail.com)
%

BasisNum = size(U,2);

if strcmpi(options.ReguParaType,'LARs')
    cardiCandi = options.Cardi;
    cardiCandi(cardiCandi > BasisNum) = [];
    cardiCandi = unique(cardiCandi);
else
    cardiCandi = 1;
end

nSmp = size(fea,2);

Y = cell(length(cardiCandi),1);
for i = 1:length(Y)
    Y{i} = zeros(nSmp,BasisNum);
end

if strcmpi(options.ReguParaType,'LARs')
    Gram = U'*U;
    Gram = max(Gram,Gram');

    for i = 1:nSmp
        eigvector_T = lars(U, fea(:,i),'lasso', -BasisNum,1,Gram,cardiCandi,0);

        [tm,tn] = size(eigvector_T);
        tCar = zeros(tn,1);
        for k = 1:tn
            tCar(k) = length(find(eigvector_T(:,k)));
        end

        for cardidx = 1:length(cardiCandi)
            ratio = cardiCandi(cardidx);
            iMin = find(tCar == ratio);
            if isempty(iMin)
                error('Card dose not exist!');
            end
            Y{cardidx}(i,:) = eigvector_T(:,iMin(end));
        end
    end

    for i = 1:length(Y)
        if cardiCandi(i)/BasisNum < 0.6
            Y{i} = sparse(Y{i});
        end
    end
elseif strcmpi(options.ReguParaType,'SLEP')
    opts=[];
    opts.rFlag=1;       % the input parameter 'ReguAlpha' is a ratio in (0, 1)
    for i = 1:nSmp
        x = LeastR(U, fea(:,i), options.ReguAlpha, opts);
        Y{1}(i,:) = x';
    end
else
    error('Method does not exist!');
end

author: ZJULearning

repo: MatlabFunc


text:
% Copyright 2016 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
% http ://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

%% Local Laplacian on Parrot.
runOnFilenames('../../images/', ...
    'source_low.png', 'target_low.png', ...
    'source.png', 'target.png');

% fprintf('Program paused, press Enter to continue.\n');
% pause;
% close all;
%
% %% Local Laplacian on Parrot (grayscale).
% runOnFilenames('../../images/', ...
%     'low_res_in.png', 'low_res_out_gray.png', ...
%     'high_res_in.png', 'high_res_ground_truth_out_gray.png');
%
% fprintf('Program paused, press Enter to continue.\n');
% pause;
% close all;

author: mahmoudnafifi


text:
function coef=comp_dgtreal_ola(f,g,a,M,Lb,phasetype)
%
%  This function implements periodic convolution using overlap-add. The
%  window g is supposed to be extended by fir2iir.

[L W]=size(f);
gl=length(g);

N=L/a;
M2=floor(M/2)+1;

% Length of extended block and padded g
Lext=Lb+gl;

% Number of blocks
Nb=L/Lb;

% Number of time positions per block
Nblock   = Lb/a;

if rem(Nblock,1)~=0
  error('The length of the time shift must devide the block length.');
end;

% Number of time positions in half extension
b2 = gl/2/a;

if rem(b2,1)~=0
  error(['The length of the time shift must devide the window length by ' ...
         'an even number.'])
end;

% Extend window to length of extended block.
gpad=fir2long(g,Lext);

coef=zeros(M2,N,W,assert_classname(f,g));

for ii=0:Nb-1

  block=comp_sepdgtreal(postpad(f(ii*Lb+1:(ii+1)*Lb,:),Lext),gpad,a,M,phasetype);

  % Large block
  coef(:,ii*Nblock+1:(ii+1)*Nblock,:) = coef(:,ii*Nblock+1:(ii+1)*Nblock,:)+block(:,1:Nblock,:);

  % Small block +
  s_ii=mod(ii+1,Nb);
  coef(:,s_ii*Nblock+1   :s_ii*Nblock+b2,:) = coef(:,s_ii*Nblock+1 ...
                                                   :s_ii*Nblock+b2,:)+ block(:,Nblock+1:Nblock+b2,:);

  % Small block -
  s_ii=mod(ii-1,Nb)+1;
  coef(:,s_ii*Nblock-b2+1:s_ii*Nblock,:) =coef(:,s_ii*Nblock-b2+1:s_ii*Nblock,:)+ block(:,Nblock+b2+1:Nblock+2*b2,:);

end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function X = regularize(X, reg)

%error('Do not use this function');

% $$$ normX = normest(X, 1e-1);

n = length(X);
I = diag(sparse(ones(n,1)));

if nargin < 2
  reg = 1;
end
normX = reg;
coeff = 1e-6;
X = X + coeff*normX*I;

% $$$ while ~success
% $$$   try
% $$$     U = chol(X + 1e-8*normX*I, varargin{:});
% $$$     success = true;
% $$$   catch
% $$$     times = times + 1;
% $$$     warning(sprintf(['Cholesky decomposition failed %d-th time, increasing ' ...
% $$$                      'regularization'], times))
% $$$     normX = normX*100;
% $$$   end
% $$$ end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function D = cosine(A,B)
% {cosine} computes the cosine distance.
%
%      D = cosine(A,B)
%
%      A: M-by-P matrix of M P-dimensional vectors
%      B: N-by-P matrix of M P-dimensional vectors
%
%      D: M-by-N distance matrix
%
% Author: Stefano Melacci (2009)
%         mela@dii.unisi.it
%         * based on the code of Vikas Sindhwani, vikas.sindhwani@gmail.com

if (size(A,2) ~= size(B,2))
    error('A and B must be of the same dimensionality.');
end

if (size(A,2) == 1) % if dim = 1...
    A = [A, zeros(size(A,1),1)];
    B = [B, zeros(size(B,1),1)];
end

aa=sum(A.*A,2);
bb=sum(B.*B,2);
ab=A*B';

% to avoid NaN for zero norms
aa((aa==0))=1;
bb((bb==0))=1;

D = real(ones(size(A,1),size(B,1))-(1./sqrt(kron(aa,bb'))).*ab);

author: viggin

repo: domain-adaptation-toolbox


text:
function xy=solveBivarMultilinSys(a,c)
%%SOLVEBIVARMULTILINSYS Given bivariate multilinear equations of the form
%       c(i)=a(1,i)+a(2,i)*x+a(3,i)*y+a(4,i)*x*y
%       for i=1,2, solve for x and y.
%
%INPUTS: a A 4X2 vector of coefficients as defined above.
%        c A length 2 vector of the values of the multilinear equations.
%
%OUTPUTS: xy A 2X2 vector of [x;y] solutions. xy(:,i) is the ith solution.
%
%EXAMPLE:
%We generate a random system and obtain both solutions. One solutions
%agrees with the [x;y] values used to generate the original c values and
%the other solution does not, though it does agree with c values.
% a=100*randn(4,2);
% x=1;
% y=-12;
% c=zeros(2,1);
% for k=1:2
%     c(k)=a(1,k)+a(2,k)*x+a(3,k)*y+a(4,k)*x*y;
% end
% xy=solveBivarMultilinSys(a,c)%One solution will be [1;-12]
%
%April 2021 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

c1=c(1);
c2=c(2);
a11=a(1,1);
a21=a(2,1);
a31=a(3,1);
a41=a(4,1);
a12=a(1,2);
a22=a(2,2);
a32=a(3,2);
a42=a(4,2);

commonTerm=-(a12*a41-a11*a42+a42*c1-a41*c2);
addTermX=a21*a32-a22*a31+commonTerm;
addTermY=a22*a31-a21*a32+commonTerm;
radTerm=sqrt(-4*(a22*a41-a21*a42)*(a12*a31-a11*a32+a32*c1-a31*c2)+(a22*a31-a21*a32+a12*a41-a11*a42+a42*c1-a41*c2)^2);

xDenom=2*(a22*a41-a21*a42);
yDenom=2*(a32*a41-a31*a42);

x1=(addTermX-radTerm)/xDenom;
y1=(addTermY+radTerm)/yDenom;
x2=(addTermX+radTerm)/xDenom;
y2=(addTermY-radTerm)/yDenom;

xy=[x1,x2;
    y1,y2];

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function test_bug1227

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_sourcewrite

[pnt, tri] = mesh_sphere(162);

source = [];
source.pos = pnt;
source.tri = tri;
source.pow = randn(162,100);
source.dimord = 'pos_time';
source.time = (0:99)/100;

cfg = [];
cfg.filename  = [tempname,'.gii'];
cfg.parameter = 'pow';
ft_sourcewrite(cfg, source);
if exist(cfg.filename, 'file')
  fprintf('gifti file created\n');
  mrig = ft_read_headshape(cfg.filename);
  delete(cfg.filename);
  fprintf('gifti file deleted\n');
end

[x, y, z]    = ndgrid(1:4,1:5,1:6);
source2      = [];
source2.pos  = [x(:) y(:) z(:)];
source2.dim  = [4 5 6];
source2.time = (0:99)/100;
source2.pow  = randn(120,100);

cfg.filename = [tempname,'.nii'];
ft_sourcewrite(cfg, source2);
if exist(cfg.filename, 'file')
  fprintf('nifti file created\n');
  mrin = ft_read_mri(cfg.filename);
  delete(cfg.filename);
  fprintf('nifti file deleted\n');
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [net,stats] = cnn_train_dag_mcn2(net, imdb, getBatch, varargin)
%CNN_TRAIN_DAG Demonstrates training a CNN using the DagNN wrapper
%    CNN_TRAIN_DAG() is similar to CNN_TRAIN(), but works with
%    the DagNN wrapper instead of the SimpleNN wrapper.

% Copyright (C) 2014-16 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).
addpath(fullfile(vl_rootnn, 'examples'));

opts.expDir = fullfile('data','exp') ;
opts.continue = true ;
opts.batchSize = 256 ;
opts.numSubBatches = 1 ;
opts.train = [] ;
opts.val = [] ;
opts.gpus = [] ;
opts.prefetch = false ;
opts.epochSize = inf;
opts.numEpochs = 300 ;
opts.learningRate = 0.001 ;
opts.weightDecay = 0.0005 ;

opts.solver = [] ;  % Empty array means use the default SGD solver
[opts, varargin] = vl_argparse(opts, varargin) ;
if ~isempty(opts.solver)
  assert(isa(opts.solver, 'function_handle') && nargout(opts.solver) == 2,...
    'Invalid solver; expected a function handle with two outputs.') ;
  % Call without input arguments, to get default options
  opts.solverOpts = opts.solver() ;
end

opts.momentum = 0.9 ;
opts.saveSolverState = true ;
opts.nesterovUpdate = false ;
opts.randomSeed = 0 ;
opts.profile = false ;

opts.parameterServer.method = 'tmove'; %'tmove' or 'mmap';
opts.parameterServer.prefix = ['matconvnet' num2str(feature('getpid'))] ;
opts.parameterServer.memoryMapFile = fullfile(tempdir, 'ramdisk', opts.parameterServer.prefix) ;

opts.derOutputs = {'objective', 1} ;
opts.extractStatsFn = @extractStats ;
opts.plotStatistics = true;
opts.postEpochFn = [] ;  % postEpochFn(net,params,state) called after each epoch; can return a new learning rate, 0 to stop, [] for no change

opts.valmode = '30samples';
opts.temporalStride = 1;
opts.backpropDepth = [];
opts.numValFrames = 3;
opts.nFramesPerVid = 5;
opts.saveAllPredScores = false;
opts.denseEval = 0;
opts.cudnnWorkspaceLimit = [];
opts.plotDiagnostics = false;
opts.augmentation = '';
opts.temporalFullConvTest = true;
[opts, varargin] = vl_argparse(opts, varargin) ;

if ~exist(opts.expDir, 'dir'), mkdir(opts.expDir) ; end
if isempty(opts.train), opts.train = find(imdb.images.set==1) ; end
if isempty(opts.val), opts.val = find(imdb.images.set==2) ; end
if isscalar(opts.train) && isnumeric(opts.train) && isnan(opts.train)
  opts.train = [] ;
end
if isscalar(opts.val) && isnumeric(opts.val) && isnan(opts.val)
  opts.val = [] ;
end

% -------------------------------------------------------------------------
%                                                            Initialization
% -------------------------------------------------------------------------

evaluateMode = isempty(opts.train) ;
if ~evaluateMode
  if isempty(opts.derOutputs)
    error('DEROUTPUTS must be specified when training.\n') ;
  end
end

% -------------------------------------------------------------------------
%                                                        Train and validate
% -------------------------------------------------------------------------

modelPath = @(ep) fullfile(opts.expDir, sprintf('net-epoch-%d.mat', ep));
modelFigPath = fullfile(opts.expDir, 'net-train.pdf') ;

start = opts.continue * findLastCheckpoint(opts.expDir) ;
if start >= 1
  fprintf('%s: resuming by loading epoch %d\n', mfilename, start) ;
  [net, state, stats] = loadState(modelPath(start)) ;
else
  state = [] ;
end

for epoch=start+1:opts.numEpochs

  % Set the random seed based on the epoch and opts.randomSeed.
  % This is important for reproducibility, including when training
  % is restarted from a checkpoint.

  rng(epoch + opts.randomSeed) ;
  prepareGPUs(opts, epoch == start+1) ;

  % Train for one epoch.
  params = opts ;
  params.epoch = epoch ;
  params.learningRate = opts.learningRate(min(epoch, numel(opts.learningRate))) ;
  params.train = opts.train(randperm(numel(opts.train))) ; % shuffle
  params.train = params.train(1:min(opts.epochSize, numel(opts.train)));
  params.val = opts.val ;
  params.imdb = imdb ;
  params.getBatch = getBatch ;
  if iscell(opts.backpropDepth)
    if epoch <= numel(opts.backpropDepth)
      params.backpropDepth = opts.backpropDepth{epoch};
    else
      params.backpropDepth = [];
    end
  else
    params.backpropDepth = opts.backpropDepth;
  end

  if numel(opts.gpus) <= 1
    [net, state] = processEpoch(net, state, params, 'train') ;
    [net, state] = processEpoch(net, state, params, 'val') ;
    if ~evaluateMode
      saveState(modelPath(epoch), net, state) ;
    end
    lastStats = state.stats ;
  else
    assert(isempty(params.backpropDepth), 'Cannot use backpropDepth in multi GPU mode.')
    spmd
      [net, state] = processEpoch(net, state, params, 'train') ;
      [net, state] = processEpoch(net, state, params, 'val') ;
      if labindex == 1 && ~evaluateMode
        saveState(modelPath(epoch), net, state) ;
      end
      lastStats = state.stats ;
    end
    lastStats = accumulateStats(lastStats) ;
  end

  stats.train(epoch) = lastStats.train ;
  stats.val(epoch) = lastStats.val ;
  clear lastStats ;
  saveStats(modelPath(epoch), stats, opts) ;

  if opts.plotStatistics
    switchFigure(1) ; clf ;
    values = [] ;   values_loss = [] ;
    leg = {} ;   leg_loss = {} ;
    for s = {'train', 'val'}
      s = char(s) ;
      for f = setdiff(fieldnames(stats.(s))', {'num', 'time','scores', 'allScores'})
        f = char(f) ;
        if isempty(strfind(f,'err'))
          leg_loss{end+1} = sprintf('%s (%s)', f, s) ;
          tmp = [stats.(s).(f)] ;
          values_loss(end+1,:) = tmp(1,:)' ;
        else
          leg{end+1} = sprintf('%s (%s)', f, s) ;
          tmp = [stats.(s).(f)] ;
          values(end+1,:) = tmp(1,:)' ;
        end
        tmp = [stats.(s).(f)];
        fprintf('%s (%s):%.3f\n', f, s, tmp(end))
      end
    end

    if ~isempty(values_loss)
      subplot(1,2,1) ; plot(1:epoch, values_loss','o-') ;
      legend(leg_loss{:},'Location', 'northoutside'); xlabel('epoch') ; ylabel('objective') ;
      subplot(1,2,2) ; plot(1:epoch, values','o-') ; ylim([0 1])
      legend(leg{:},'Location', 'northoutside') ; xlabel('epoch') ; ylabel('error') ;
      grid on ;
      drawnow ;
      print(1, modelFigPath, '-dpdf') ;
    end
  end

  if ~isempty(opts.postEpochFn)
    if nargout(opts.postEpochFn) == 0
      opts.postEpochFn(net, params, state) ;
    else
      lr = opts.postEpochFn(net, params, state) ;
      if ~isempty(lr), opts.learningRate = lr; end
      if opts.learningRate == 0, break; end
    end
  end
end

% With multiple GPUs, return one copy
if isa(net, 'Composite'), net = net{1} ; end

% -------------------------------------------------------------------------
function [net, state] = processEpoch(net, state, params, mode)
% -------------------------------------------------------------------------
% Note that net is not strictly needed as an output argument as net
% is a handle class. However, this fixes some aliasing issue in the
% spmd caller.

% initialize with momentum 0
if isfield(state, 'momentum'), state.solverState = state.momentum; state.momentum = []; end
if isempty(state) || isempty(state.solverState)
  state.solverState = cell(1, numel(net.params)) ;
  state.solverState(:) = {0} ;
end

% move CNN  to GPU as needed
numGpus = numel(params.gpus) ;
if numGpus >= 1
  net.move('gpu') ;
  for i = 1:numel(state.solverState)
    s = state.solverState{i} ;
    if isnumeric(s)
      state.solverState{i} = gpuArray(s) ;
    elseif isstruct(s)
      state.solverState{i} = structfun(@gpuArray, s, 'UniformOutput', false) ;
    end
  end
end
if numGpus > 1
  parserv = ParameterServer(params.parameterServer) ;
  net.setParameterServer(parserv) ;
else
  parserv = [] ;
end

% profile
if params.profile
  if numGpus <= 1
    profile clear ;
    profile on ;
  else
    mpiprofile reset ;
    mpiprofile on ;
  end
end

num = 0 ;
epoch = params.epoch ;
subset = params.(mode) ;
adjustTime = 0 ;
moreopts = [];

stats.num = 0 ; % return something even if subset = []
stats.time = 0 ;

stats.scores = [] ;
stats2.err1 = 0;
stats2.err5 = 0;
net.backpropDepth = params.backpropDepth;

if ~strcmp(mode,'train')
  net.mode = 'test';
  dataset =  ceil(params.imdb.images.set(subset(1))/2);
  nClasses = numel(params.imdb.classes.name);
  if nClasses < 5
    nClasses = numel(params.imdb.classes.name{dataset});
  end
  stats2.scores = zeros(nClasses, numel(subset));

  moreopts.frameSample = 'uniformly';
  moreopts.augmentation = 'uniform';
  moreopts.keepFramesDim = true; % make getBatch output 5 dimensional

  if strcmp(params.valmode,'30samples')
    % sample less frames and crops:
    moreopts.numAugments = 6;
    moreopts.nFramesPerVid = 5;
  elseif strcmp(params.valmode,'centreSamplesFast')
    % sample less frames and crops:
    moreopts.numAugments = 2;
    moreopts.nFramesPerVid = 3;
  elseif strcmp(params.valmode,'250samples') ,
    moreopts.numAugments = 10;
    moreopts.nFramesPerVid = 25;
  elseif strcmp(params.valmode,'dense')
    moreopts.augmentation = 'none';
    moreopts.numAugments = 0;
    moreopts.nFramesPerVid = 25;
    moreopts.keepFramesDim = true;
    params.batchSize = numlabs;
    params.numSubBatches = numlabs;
  elseif strcmp(params.valmode,'temporalStrideRandom')
    moreopts.nFrameStack = params.nFramesPerVid;
    moreopts.temporalStride = ceil(median(params.temporalStride));
    moreopts.temporalStride = max(params.temporalStride);
    params.batchSize =  32*numlabs ;
    params.numSubBatches = params.batchSize; % has to be
    moreopts.nFramesPerVid = params.numValFrames;
  end

  if params.denseEval
    moreopts.augmentation = 'none';
    moreopts.numAugments = 2;
  end
  if params.temporalFullConvTest
    params.nFramesPerVid = params.numValFrames;
  end
  pred_layers = [];
  for l=1:numel(net.layers)
    if isempty( net.layers(l).params ), continue; end;
    if size(net.params(net.getParamIndex(net.layers(l).params{1})).value,4) == nClasses || ...
        size(net.params(net.getParamIndex(net.layers(l).params{1})).value,5) == nClasses % 3D FC layer
          pred_layers = [pred_layers net.layers(l).outputIndexes];
          net.vars(net.layers(l).outputIndexes).precious = 1;
    end
    if isa(net.layers(l).block, 'dagnn.LSTM')
      pred_layers = [pred_layers net.layers(l).outputIndexes(1)];
      net.vars(net.layers(l).outputIndexes(1)).precious = 1 ;
    end
  end

  if params.saveAllPredScores
    stats2.allScores = zeros(numel(pred_layers),moreopts.numAugments* moreopts.nFramesPerVid/params.nFramesPerVid, nClasses,  numel(subset));
  end
else
  net.mode = 'normal';
  moreopts = [];
end

start = tic ;
for t=1:params.batchSize:numel(subset)
  fprintf('%s: epoch %02d: %3d/%3d:', mode, epoch, ...
          fix((t-1)/params.batchSize)+1, ceil(numel(subset)/params.batchSize)) ;
  batchSize = min(params.batchSize, numel(subset) - t + 1) ;

  for s=1:params.numSubBatches
    % get this image batch and prefetch the next
    batchStart = t + (labindex-1) + (s-1) * numlabs ;
    batchEnd = min(t+params.batchSize-1, numel(subset)) ;
    batch = subset(batchStart : params.numSubBatches * numlabs : batchEnd) ;
    num = num + numel(batch) ;
    if numel(batch) == 0, continue ; end

    inputs = params.getBatch(params.imdb, batch, moreopts) ;

    moreopts.frameList = [];
    if strcmp(net.mode, 'test') && strcmp(params.valmode,'temporalStrideRandom')
      for i = 2:4:numel(inputs)
        sz = size(inputs{i});
        inputs{i} = gather(inputs{i});
        nFramesPerVid = sz(5)/moreopts.numAugments;
        chunks = ceil(nFramesPerVid / params.nFramesPerVid);
        inputs{i} = reshape(inputs{i}, sz(1), sz(2), sz(3),   [], params.nFramesPerVid);
        inputs{i} = permute(inputs{i} , [1 2 3 5 4]);
      end
    end
    net.meta.curNumFrames = repmat(size(inputs{2},4) / numel(inputs{4}),1,numel(net.layers)); % nFrames = instances/labels

    net.meta.curBatchSize = numel(batch);
    inputs{end+1} = 'inputSet'; inputs{end+1} = ceil(params.imdb.images.set(batch)/2); % dataset

    if params.prefetch
      if s == params.numSubBatches
        batchStartNext = t + (labindex-1) + params.batchSize ;
        batchEndNext = min(t+2*params.batchSize-1, numel(subset)) ;
      else
        batchStartNext = batchStart + numlabs ; batchEndNext = batchEnd;
      end
      nextBatch = subset(batchStartNext : params.numSubBatches * numlabs : batchEndNext) ;
      if ~isempty(nextBatch)
        moreopts.frameList = params.getBatch(params.imdb, nextBatch, moreopts) ;
      else
        moreopts.frameList = NaN ;
      end
    end
    if ndims(inputs{2})>4  % average over frames
      dataset = inputs{end}(1);
      nClasses = numel(params.imdb.classes.name);
      if nClasses < 5
        nClasses = numel(params.imdb.classes.name{dataset});
      end

      frame_predictions = cell(numel(pred_layers),size(inputs{2},5));

      for fr = 1:size(inputs{2},5)
        frame_inputs = inputs;
        net.meta.curNumFrames = repmat(size(inputs{2},4) / numel(inputs{4}),1,numel(net.layers)); % nFrames = instances/labels
        for i = 2:4:numel(inputs)
          if size(frame_inputs{i},5) > 1
            frame_inputs{i}=frame_inputs{i}(:,:,:,:,fr);
          end
        end
        if strcmp(mode, 'train')
          net.accumulateParamDers = (s ~= 1) ;
          net.eval(frame_inputs, derOutputs) ;
        else
          net.eval(frame_inputs) ;
        end
        [frame_predictions{:,fr}] = deal(net.vars(pred_layers).value) ;
      end

      tmp = [];
      for k = 1:numel(pred_layers)
        frame_predictions(k,:)= cellfun(@(x) mean(mean(x,1),2), frame_predictions(k,:), 'UniformOutput', false);
        if numel(batch) == 1
          frame_predictions(k,:)= cellfun(@(x) mean(x,4), frame_predictions(k,:), 'UniformOutput', false);
        end
        tmp = [tmp; frame_predictions{k,:}];
      end
      frame_predictions = tmp;

      if min(net.meta.curNumFrames) > 1
        frame_predictions = mean(frame_predictions,4);
      end
      if  params.saveAllPredScores
        stats2.allScores(:,:,:,batchStart : params.numSubBatches * numlabs : batchEnd) = gather(frame_predictions);
      end
      % average over time (dim+1) and spatial locations and batches
      frame_predictions = mean(mean(mean(frame_predictions),1),2);
      [err1, err5] = error_multiclass(params, inputs{4}, gather(frame_predictions));
      stats2.err1 = (stats2.err1 + err1);
      stats2.err5 = (stats2.err5 + err5);
    else  % inputs four dimensional
      if strcmp(mode, 'train')
        net.mode = 'normal' ;
        net.accumulateParamDers = (s ~= 1) ;
        net.eval(inputs, params.derOutputs, 'holdOn', s < params.numSubBatches) ;
      else
        net.mode = 'test' ;
        net.eval(inputs) ;
      end
    end
    if strcmp(mode, 'val') && ndims(inputs{2})>4
      stats2.scores(:, batchStart : params.numSubBatches * numlabs : batchEnd) = squeeze(gather(frame_predictions));
    end
  end

  % Accumulate gradient.
  if strcmp(mode, 'train')
    if ~isempty(parserv), parserv.sync() ; end
    state = accumulateGradients(net, state, params, batchSize, parserv) ;
  end

  % Get statistics.
  time = toc(start) + adjustTime ;
  batchTime = time - stats.time ;
  stats.num = num ;
  stats.time = time ;
  stats = params.extractStatsFn(stats,net) ;
  if ndims(inputs{2})>4  % average over frames
    for f = fieldnames(stats2)'
        f = char(f) ;  stats.(f) = stats2.(f);
    end
  end
  currentSpeed = batchSize / batchTime ;
  averageSpeed = (t + batchSize - 1) / time ;
  if t == 3*params.batchSize + 1
    % compensate for the first three iterations, which are outliers
    adjustTime = 4*batchTime - time ;
    stats.time = time + adjustTime ;
  end

  fprintf('lr: %.0e, %.1f (%.1f) Hz',params.learningRate, averageSpeed, currentSpeed) ;
  for f = setdiff(fieldnames(stats)',  {'num', 'time','scores', 'allScores'})
    f = char(f) ;
    if ndims(inputs{2})>4  && any(strcmp(f, {'err1', 'err5'}))
      n = (t + batchSize - 1) / max(1,numlabs) ;
      stats.(f) = stats.(f) / n;
    end
      fprintf(' %s:%.3f', f, stats.(f)) ;
  end
  fprintf('\n') ;

  %  debug info
  if params.plotDiagnostics && numGpus <= 1
    figure(2) ; net.diagnose('Vars',1,'Params',1,'Time',1) ; drawnow ;
  end
end

% Save back to state.
state.stats.(mode) = stats ;
if params.profile
  if numGpus <= 1
    state.prof.(mode) = profile('info') ;
    profile off ;
  else
    state.prof.(mode) = mpiprofile('info');
    mpiprofile off ;
  end
end
if ~params.saveSolverState
  state.solverState = [] ;
else
  for i = 1:numel(state.solverState)
    s = state.solverState{i} ;
    if isnumeric(s)
      state.solverState{i} = gather(s) ;
    elseif isstruct(s)
      state.solverState{i} = structfun(@gather, s, 'UniformOutput', false) ;
    end
  end
end

net.reset() ;
net.move('cpu') ;

% -------------------------------------------------------------------------
function state = accumulateGradients(net, state, params, batchSize, parserv)
% -------------------------------------------------------------------------
numGpus = numel(params.gpus) ;
otherGpus = setdiff(1:numGpus, labindex) ;

for p=1:numel(net.params)

  if ~isempty(parserv)
    parDer = parserv.pullWithIndex(p) ;
  else
    parDer = net.params(p).der ;
  end
  if isempty(parDer)
    % fprintf('empty param for: %s\n', net.params(p).name) ;
    continue;
  end
  switch net.params(p).trainMethod

    case 'average' % mainly for batch normalization
      thisLR = net.params(p).learningRate ;
      net.params(p).value = vl_taccum(...
          1 - thisLR, net.params(p).value, ...
          (thisLR/batchSize/net.params(p).fanout),  parDer) ;

    case 'gradient'
      thisDecay = params.weightDecay * net.params(p).weightDecay ;
      thisLR = params.learningRate * net.params(p).learningRate ;

      if thisLR>0 || thisDecay>0
        % Normalize gradient and incorporate weight decay.
        parDer = vl_taccum(1/batchSize, parDer, ...
                           thisDecay, net.params(p).value) ;

        if isempty(params.solver)
          % Default solver is the optimised SGD.
          % Update momentum.
          state.solverState{p} = vl_taccum(...
            params.momentum, state.solverState{p}, ...
            -1, parDer) ;

          % Nesterov update (aka one step ahead).
          if params.nesterovUpdate
            delta = params.momentum * state.solverState{p} - parDer ;
          else
            delta = state.solverState{p} ;
          end

          % Update parameters.
          net.params(p).value = vl_taccum(...
            1,  net.params(p).value, thisLR, delta) ;

        else
          % call solver function to update weights
          [net.params(p).value, state.solverState{p}] = ...
            params.solver(net.params(p).value, state.solverState{p}, ...
            parDer, params.solverOpts, thisLR) ;
        end
      end
    otherwise
      error('Unknown training method ''%s'' for parameter ''%s''.', ...
        net.params(p).trainMethod, ...
        net.params(p).name) ;
  end
end

% -------------------------------------------------------------------------
function stats = accumulateStats(stats_)
% -------------------------------------------------------------------------

for s = {'train', 'val'}
  s = char(s) ;
  total = 0 ;

  % initialize stats stucture with same fields and same order as
  % stats_{1}
  stats__ = stats_{1} ;
  names = fieldnames(stats__.(s))' ;
  values = zeros(1, numel(names)) ;
  fields = cat(1, names, num2cell(values)) ;
  stats.(s) = struct(fields{:}) ;

  for g = 1:numel(stats_)
    stats__ = stats_{g} ;
    num__ = stats__.(s).num ;
    total = total + num__ ;

    for f = setdiff(fieldnames(stats__.(s))', 'num')
      f = char(f) ;
      stats.(s).(f) = stats.(s).(f) + stats__.(s).(f) * num__ ;

      if g == numel(stats_)
        stats.(s).(f) = stats.(s).(f) / total ;
      end
    end
  end
  stats.(s).num = total ;
end

% -------------------------------------------------------------------------
function stats = extractStats(stats, net)
% -------------------------------------------------------------------------
sel = find(cellfun(@(x) isa(x,'dagnn.Loss')||isa(x,'dagnn.Loss'), {net.layers.block})) ;
for i = 1:numel(sel)
  if net.layers(sel(i)).block.ignoreAverage, continue; end
  stats.(net.layers(sel(i)).name) = net.layers(sel(i)).block.average ;
end

% -------------------------------------------------------------------------
function saveState(fileName, net_, state)
% -------------------------------------------------------------------------
net = net_.saveobj() ;
state.stats = [];
save(fileName, 'net', 'state') ;

% -------------------------------------------------------------------------
function saveStats(fileName, stats, opts)
% -------------------------------------------------------------------------
if exist(fileName)
  save(fileName, 'stats', 'opts', '-append') ;
else
  save(fileName, 'stats', 'opts') ;
end

% -------------------------------------------------------------------------
function [net, state, stats] = loadState(fileName)
% -------------------------------------------------------------------------
state = [];
load(fileName, 'net', 'state', 'stats') ;
net = dagnn.DagNN.loadobj(net) ;
if isempty(whos('stats'))
  if isfield(state, 'stats')
    stats = state.stats;
  else
    error('Epoch ''%s'' was only partially saved. Delete this file and try again.', ...
        fileName) ;
  end
end

% -------------------------------------------------------------------------
function epoch = findLastCheckpoint(modelDir)
% -------------------------------------------------------------------------
list = dir(fullfile(modelDir, 'net-epoch-*.mat')) ;
tokens = regexp({list.name}, 'net-epoch-([\d]+).mat', 'tokens') ;
epoch = cellfun(@(x) sscanf(x{1}{1}, '%d'), tokens) ;
epoch = max([epoch 0]) ;

% -------------------------------------------------------------------------
function switchFigure(n)
% -------------------------------------------------------------------------
if get(0,'CurrentFigure') ~= n
  try
    set(0,'CurrentFigure',n) ;
  catch
    figure(n) ;
  end
end

% -------------------------------------------------------------------------
function clearMex()
% -------------------------------------------------------------------------
clear vl_tmove vl_imreadjpeg ;

% -------------------------------------------------------------------------
function prepareGPUs(opts, cold)
% -------------------------------------------------------------------------
numGpus = numel(opts.gpus) ;
if numGpus > 1
  % check parallel pool integrity as it could have timed out
  pool = gcp('nocreate') ;
  if ~isempty(pool) && pool.NumWorkers ~= numGpus
    delete(pool) ;
  end
  pool = gcp('nocreate') ;
  if isempty(pool)
    parpool('local', numGpus) ;
    cold = true ;
  end

end
if numGpus >= 1 && cold
  fprintf('%s: resetting GPU\n', mfilename)
  clearMex() ;
  if numGpus == 1
    gpuDevice(opts.gpus)
  else
    spmd
      clearMex() ;
      gpuDevice(opts.gpus(labindex))
    end
  end
end

function [err1, err5] = error_multiclass(opts, labels, predictions)
% -------------------------------------------------------------------------
[~,predictions] = sort(predictions, 3, 'descend') ;
error = ~bsxfun(@eq, predictions, reshape(labels, 1, 1, 1, [])) ;
err1 = sum(sum(sum(error(:,:,1,:)))) ;
err5 = sum(sum(sum(min(error(:,:,1:5,:),[],3)))) ;

author: feichtenhofer

repo: st-resnet


text:
function asa121_test01 ( )

%*****************************************************************************80
%
%% TEST01 demonstrates the use of TRIGAMMA.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01:\n' );
  fprintf ( 1, '  TRIGAMMA computes the Trigama function.\n' );
  fprintf ( 1, '  Compare the result to tabulated values.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '          X         ' );
  fprintf ( 1, 'FX                        FX2\n' );
  fprintf ( 1, '                    ' );
  fprintf ( 1, '(Tabulated)               (TRIGAMMA)              DIFF\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = trigamma_values ( n_data );

    if ( n_data == 0 )
      break
    end

    [ fx2, ifault ] = trigamma ( x );

    fprintf ( 1, '  %12.4f  %24.16e  %24.16e  %10.4e\n', ...
    x, fx, fx2, abs ( fx - fx2 ) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% AMG TEST II: Different boundary conditions
%
% We consider the linear finite element discretization of Poisson equation
% on the unstructured mesh with Dirichlet or Neumann boundary conditions.

%%
load lakemesh
%% Dirichlet problem
close all;
[N,itStep,time,err] = amgtest(node,elem,1);
%%
colHeaders = {'Unknowns','Iterations','Time (sec)','Error'};
makeHtmlTable([N itStep time err],[],[],colHeaders,[],6);
%%
r = showrate(N,time,2);
xlabel('N'); ylabel('Time');
title(['Complexity is N^{' num2str(r) '}'] ,'Fontsize', 14);

%% Neumann problem
close all;
[N,itStep,time,err] = amgtest(node,elem,0);
%%
colHeaders = {'Unknowns','Iterations','Time (sec)','Error'};
makeHtmlTable([N itStep time err],[],[],colHeaders,[],6);
%%
r = showrate(N,time,2);
xlabel('N'); ylabel('Time');
title(['Complexity is N^{' num2str(r) '}'] ,'Fontsize', 14);

repo: Aero-matlab


text:
                    function lagre1
% ***************************************************************
%                   P r o g r a m   LAGRE1
% ***************************************************************
%
%  PURPOSE:
%    Derive differential equation of motion of a mechanical system
%    with one degree of freedom by means of Lagrange equation
%          d/dt(dL/dqt) - dL/dq = QN(t,q,qt)
%    Resolve the equation numerically and plots the graphics
%    of the coordinate, velocity and phase plane.
%    If possible, the program could solve the problem analytically.
%
%  INPUT DATA:
%    L    - expression of the Lagrangian L = L(t,q,qt);
%    QN   - generalized non potential force QN = QN(t,q,qt);
%    q0   - initial value of the coordinate;
%    qt0  - initial value of the velocity;
%    Tend - upper bound of the integration;
%    eps  - precision of the calculations;
%    np   - number of parameters .
%    P{1}, P{2}, ..., P{np} - names of the parameters (array of cells);
%
%  NOTES:
%   1. The coordinate is designed by the symbol 'q' and velocity by 'qt';
%   2. The physical names of the parameters are assigned to the
%      cells of the array P like this: P{1}='m', P{2}='c',...;
%   3. For analytical solution the values of Tend, eps, np and P are not
%      needed.
%   4. Initial values q0, qt0 must to be entered as strings, even though
%      they represent numbers!
%   5. All the data can be input from file or in interactive mode.
%
%  EXAMPLE of DATA FILE:
%   % Data for problem ...
%    L    = '1/2*a*qt^2 + 9.81*l*cos(q)';
%    QN   = '-k*qt';
%    q0   = '0.5'; ( or q0 = 'q0';)
%    qt0  = '10';  ( or qt0 = 'qt0';)
%    Tend = 20;
%    eps  = 1.e-7;
%    np   = 3;
%    P{1} = 'a';
%    P{2} = 'l';
%    P{3} = 'k';

% ---------------------------------------------------------
%                DATA INPUT OF THE PROBLEM
% =========================================================

 clear
 disp(' ');
 disp(' How will you input the data ?    ');
 disp('     1. From a data file;         ');
 disp('     2. In interactive mode.      ');
 ans = input(' Number of Your choice : '  );
 flag = 0;
 if ans == 1
    while 1
       disp(' ');
       indat = input(' Input the name of the data file :', 's');
       if exist([cd,'\',indat]) % Search only in current directory
          eval(indat);
          flag = 1; break  % Successful
        else               % Unsuccessful
          disp(' ');
          disp([' File ',indat,' not exist!'])
          disp(' You have to:')
          disp(' 1. Enter another DATA file name, or')
          disp(' 2. Input the DATA interactively !')
          ans2 = input(' Your choice, please: ');
          if ans2 == 2, break , end
       end
    end
 end
if flag == 0
    %Input of the data in-line mode
    L    = input(' Expression of the Lagrangian  L : ','s' );
    QN   = input(' Generalized unpotential force QN : ','s');
    q0   = input(' Initial value of the coordinate q0 : '  );
    qt0  = input(' Initial value of the velocity qt0 : '   );
    Tend = input(' Upper bound of the integration Tend : ' );
    eps  = input(' Precision of the calculations eps : '   );
    np   = input(' Number of parameters  np : '            );
    % Asigning names of the parameters
    if np > 0
       disp(' ');
       disp(' Enter names of the parameters:')
       for i = 1:np
           ii = num2str(i);
           P{i} = input([' Name of parameter P',ii,': '],'s');
       end
    end
 end

% ---------------------------------------------------------
%       Deriving Differential Equation of Motion
% =========================================================

 syms t q qt qtt Dq D2q
 Lqt = diff(L,qt);
 deq = diff(Lqt,q)*qt + diff(Lqt,qt)*qtt + ...
       diff(Lqt,t) - diff(L,q) - QN;
 deq = subs(deq,{qt,qtt},{Dq,D2q});
 deq = simple(deq);
 disp(' ');
 disp('    D i f f e r e n t i a l   E q u a t i o n ')
 disp(' ***********************************************')
 disp(' ');
 disp([' DEq: ', char(deq),' = 0'])
%
% ---------------------------------------------------------
%                  ANALYTICAL SOLUTION
% =========================================================

 disp(' ');
 ans = input(' Would you like analitical solution? (Y/N): ','s');
 if ans =='Y' | ans == 'y'
    if ~isstr(q0) , q0  = num2str(q0) ; end % Repairing user
    if ~isstr(qt0), qt0 = num2str(qt0); end % input errors!
    inicond = ['q(0)=',q0,',Dq(0)=',qt0];
    q = dsolve(char(deq), inicond, 't');
    if ~isempty(q)
        disp(' ');
        disp('   Low of Motion   ');
        disp(' ***************** ');
        disp(' ');
        disp('q = '); pretty(q)
        disp(' ');
        fname = input(' Name of file to write solution: ','s');
        save(fname, 'q');
    end
       disp(' ');
       ans = input(' Would you like numerical solution? (Y/N): ','s');
       if ans == 'N' | ans == 'n', return, end
       q = 'q'; % Clear analytical solution from q !
 end

% ---------------------------------------------------------
%                  NUMERICAL SOLUTION
% =========================================================

% Input the name of the file-function
disp(' ');
fname = input(' Name of the File-function to be generated: ','s');
flag1 = 'Y';
if exist([cd,'\',fname]) % Search only in curent directory!
    disp(' ');
    disp([' A file-function with name ',fname,' already exist !'])
    flag1 = input(' Overwrite it ? (Y/N): ', 's');
end

% ---------------------------------------------------------
%              GENERATING THE FILE-FUNCTION
% ---------------------------------------------------------

if ( flag1 == 'Y' | flag1 == 'y' )
   qtt = solve(deq, 'D2q'); % Second derivative
   qtt = subs(qtt, {q,Dq},{'y(1)','y(2)'}); % Canonization
% Opening the file to write file-function
   [Fid,mes] = fopen([fname,'.m'],'wt');
% Generating the string with physical parameters: m, c ...
   strpar = '';
   for j = 1:np
      strpar = [strpar,',',P{j}];
   end
   disp(' ');
   titl = input(' Denomination of the Problem: ','s');
%       Writing the headline of the File-function
   fprintf(Fid,['function yt = ',fname,'(t,y',strpar,')\n']);
   fprintf(Fid,['%% ',titl]);
%       Writing the first derivatives
   fprintf(Fid,'\n%% The first derivatives\n');
   fprintf(Fid, '  yt(1) = y(2); \n');
   %fprintf(Fid,['  yt(2) = ',char(b/a),'; \n']);
   fprintf(Fid,['  yt(2) = ',char(qtt),'; \n']);
   fprintf(Fid,'  yt = yt''; \n');
   fprintf(Fid,['%% *** End of File-function ',fname,' ***']);
   fclose(Fid);
   edit(fname)
end

% ---------------------------------------------------------
%        INTEGRATION AND VISUALIZATION OF THE REZULTS
% ---------------------------------------------------------

flag2 = 0;
% Initial entering values of the parameters and generating
% the string with parameters 'P{1}, P{2}, ..., P{np}' to be
% passed to the File-function as actual arguments
if np > 0
   PP = P; % Saving the physical names of the parameters in PP
   parameters = ' ';
   disp(' ');
   disp(' Input the numerical values of the parameters: ')
   for i = 1:np
       i = num2str(i);
       eval(['P{',i,'}=input([''   '',P{',i,'},'' = '']);']);
       parameters = [parameters,',P{',i,'}'];
   end
 else
   parameters = [];
end
% Check-up type of q0 and qt0 and correct
% it if needed
if ischar(q0)
    q0 = str2num(q0);
    if isempty(q0), q0 = input(' q0 = '); end
end
if ischar(qt0)
    qt0 = str2num(qt0);
    if isempty(qt0), qt0 = input(' qt0 = '); end
end
while 1
    if flag2 == 1
        disp(' ');
        eps  = input(' Precision of the computations eps: ');
        Tend = input(' Upper bound of the integration Tend: ');
        q0   = input(' Initial coordinate q0: ');
        qt0  = input(' Initial velocity qt0: ');
        if np > 0
          P = PP; % Restoring the names of the parameters !
          disp(' ');
          disp(' Input the numerical values of the parameters: ')
          for i = 1:np
              i = num2str(i);
              eval(['P{',i,'}=input([''   '',P{',i,'},'' = '']);']);
          end
        end
    end
    y0 = [q0 qt0]; % initial conditions
    options = odeset('AbsTol',eps,'RelTol',100*eps);
    % Choosing of the Solver
    disp('                                        ');
    disp('      Choose the proper Solver:         ');
    disp('  -------------------------------       ');
    disp(' A. Non stiff differential equations    ');
    disp('   1. ode45   - middle precision;       ');
    disp('   2. ode23   - low precision;          ');
    disp('   3. ode113  - from low to upper.      ');
    disp('                                        ');
    disp(' B. Stiff differential equations        ');
    disp('   1. ode15s  - from low to upper;      ');
    disp('   2. ode23s  - low precision;          ');
    disp('   3. ode23t  - middle precision;       ');
    disp('   4. ode23tb - low precision.          ');
    disp('                                        ');
    solver = input(' The name of the Solver: ','s');

    % Integration of the Differential Equations

    eval(['[t,y] = feval(solver,eval([''@'',fname]),',...
                  '[0 Tend],y0,options',parameters,');']);
    % Plotting graphs
    tmin  = min(t); tmax = max(t);
    y1min = min(y(:,1));
    y1max = max(y(:,1));
    y2min = min(y(:,2));
    y2max = max(y(:,2));
    dy1   = y1max - y1min;
    dy2   = y2max - y2min;
    xmin  = y1min - 0.1*dy1;
    xmax  = y1max + 0.1*dy1;
    ymin  = y2min - 0.1*dy2;
    ymax  = y2max + 0.1*dy2;
    % Coordinate q = q(t)
    figure % 1
    comet(t,y(:,1))
    plot(t,y(:,1),[tmin tmax],[0 0],'k'), grid on
    axis([tmin, tmax, xmin, xmax]);
    set(gca,'FontName','Arial Cyr','FontSize',12);
    title('Low of motion {\itq} = {\itq}({\itt})')
    xlabel('{\itt}'); ylabel('{\itq}'); pause
    % Velocity qt = qt(t)
    figure % 2
    comet(t,y(:,2))
    plot(t,y(:,2),[tmin tmax],[0 0],'k'), grid on
    axis([tmin, tmax, ymin, ymax]);
    set(gca,'FontName','Arial','FontSize',12);
    title('Velocity {\it qt} = {\it qt}({\itt})')
    xlabel('{\itt}'); ylabel('{\it qt}'); pause
    % Coordinate and Velocity
    figure % 3
    subplot(2,1,1), plot(t,y(:,1),[tmin tmax],[0 0],'k')
    grid on, axis([tmin, tmax, xmin, xmax]);
    set(gca,'FontName','Arial','FontSize',12);
    title('Low of motion {\itq} = {\itq}({\itt})')
    subplot(2,1,2), plot(t,y(:,2),[tmin tmax],[0 0],'k')
    grid on, axis([tmin, tmax, ymin, ymax]);
    set(gca,'FontName','Arial','FontSize',12);
    title('Velocity {\it qt} = {\it qt}({\itt})'), pause
    % Phase Plane
    figure % 4
    subplot(1,1,1)
    comet(y(:,1),y(:,2))
    plot(y(:,1),y(:,2), [xmin, xmax],[0 0],'k',...
                  [0 0],[ymin, ymax],'k'), grid on
    axis([xmin, xmax, ymin, ymax]);
    set(gca,'FontName','Arial Cyr','FontSize',12);
    title(' Phase Plane {\it qt} = {\it qt}({\itq})')
    xlabel('{\itq}'), ylabel('{\it qt}'), pause
    flag2 = 1;
    close all
    disp(' ');
    ans = input(' Would you like to continue? (Y/N): ','s');
    if ans == 'n' | ans == 'N', break, end
end

%  *************** End of Program LAGRE1 ******************

author: Sable

repo: mcbench-benchmarks


text:
function correct_rate = test_svm_classifier(label_ind)
% Train an SVM classifier with the reduced TIDIGITS set

TRAIN_DIR = 'temp/train';
TEST_DIR = 'temp/test';
LABEL_IND = label_ind;

USE_PCA = false;
NUM_OF_COEFFS = 5;

[train_obj, train_features] = get_files_and_mfcc_features(TRAIN_DIR, LABEL_IND, '*.wav');
[test_obj, test_features] = get_files_and_mfcc_features(TEST_DIR, LABEL_IND, '*.wav');

% Remove NaNs from cells
%test_features.mfcc_delta.Std = remove_nans_from_cell(test_features.mfcc_delta.Std);
%test_features.mfcc_delta.Mean = remove_nans_from_cell(test_features.mfcc_delta.Mean);

train = [train_features.mfcc.Mean' train_features.mfcc.Std' train_features.mfcc_delta.Mean' train_features.mfcc_delta.Std' ...
         train_features.centroid.Mean' train_features.centroid.Std' ...
         train_features.rms.Mean' train_features.rms.Std'];
train_labels = get(train_obj, 'Label')';

test = [test_features.mfcc.Mean' test_features.mfcc.Std' test_features.mfcc_delta.Mean' test_features.mfcc_delta.Std' ...
        test_features.centroid.Mean' test_features.centroid.Std' ...
        test_features.rms.Mean' test_features.rms.Std'];
test_labels = get(test_obj, 'Label')';

% prewhiten requires having drtoolbox added to Matlab path.
% Prewhitening basically gets rid of correlated features, or features with low
% variance.
train = prewhiten(train);
test = prewhiten(test);

% save features train train_labels test test_labels;
% load features;

if USE_PCA
    [~, train_pcvec] = pca(train);
    train_pcvec = train_pcvec(1:NUM_OF_COEFFS, :);
    train = train * train_pcvec';
    test = test * train_pcvec';
end

% s = svmtrain(train, train_labels);
% c = svmclassify(s, test);
c = multisvm(train, train_labels, test, 'tolkkt', 1e-2, 'kktviolationlevel', 0.1);
u = unique(train_labels);
display(['Number of unique labels: ' num2str(length(u))]);
correct = strcmp(c, test_labels);
confusionmat(test_labels, c)
classperf(test_labels, c)
svm_correct_rate = sum(correct)/length(correct);
display(svm_correct_rate);

% classify using GMM
NUM_OF_GAUSSIANS = 10;
NUM_OF_ITERATIONS = 20;
[mu_train, sigma_train, c_train] = ...
    GMM.gmm_training(train', train_labels, NUM_OF_GAUSSIANS, ...
                     NUM_OF_ITERATIONS);
c = GMM.gmm_classification(test', mu_train, sigma_train, c_train);
correct = strcmp(c, test_labels);
gmm_correct_rate = sum(correct)/length(correct);
confusionmat(test_labels, u(c))
classperf(test_labels, u(c))
display(gmm_correct_rate);

% classify using K-NN
K = 3; % number of neighbors to use
c = knnclassify(test, train, train_labels, K);
correct = strcmp(c, test_labels);
knn_correct_rate = sum(correct)/length(correct);
confusionmat(test_labels, c)
classperf(test_labels, c)
display(knn_correct_rate);

% classify using MIR
result = mirclassify(train_obj, {train_features.mfcc.Mean, train_features.mfcc.Std, train_features.centroid.Mean, train_features.centroid.Std}, ...
    test_obj, {test_features.mfcc.Mean, test_features.mfcc.Std, test_features.centroid.Mean, test_features.centroid.Std});
mir_correct_rate = get(result, 'Correct');
display(mir_correct_rate);

correct_rate = [svm_correct_rate, gmm_correct_rate, knn_correct_rate, ...
    mir_correct_rate];

end

function without_nans = remove_nans_from_cell(c)
    without_nans = c;
    len = length(c{1});
    without_nans(cellfun(@(x) all(isnan(x)), without_nans)) = mat2cell(zeros(len, 1));
    without_nans = cell2mat(without_nans);
end

function [svm_struct, groups] = train_svm(features, labels)
    groups = unique(labels);
    num_of_groups = length(groups);
%     num_of_samples = length(labels);
    svm_struct = cell(num_of_groups);

    for i = 1:num_of_groups
        for j = 1:num_of_groups
            i_features = features(cell2mat(labels) == groups{i}, :);

            if (j < i)
                j_features = features(cell2mat(labels) == groups{j}, :);
            else
                j_features = features(cell2mat(labels) ~= groups{i}, :);
            end;

            ij_features = [i_features; j_features];
            binary_labels = [ones(size(i_features,1), 1); zeros(size(j_features,1), 1)];
            svm_struct{i,j} = svmtrain(ij_features, binary_labels);
        end
    end
end

function class = classify_svm(svm_struct, features, groups)
    num_of_groups = length(svm_struct);
    num_of_samples = size(features, 1);
    labels = zeros(num_of_groups, num_of_groups, num_of_samples);

    for i = 1:num_of_groups
        for j = 1:num_of_groups
            labels(i, j, :) = svmclassify(svm_struct{i, j}, features);
        end
    end

    win_times = squeeze(sum(labels, 1));
    [~, class] = max(win_times);
    class = groups(class);
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
%  computing the Householder transformation  w = v*H  with
%            H = I - 2*u*u'
%
%  syntax  >> w = HouseholderTransformRight(u,v)


text: unfittedType = 'BOUNDARY';


text:
% Set Manuel's color preferences for programming with Matlab and Fortran

%% Setup the color pallet

% divide by 256 for Matlab
sol.base03  = [ 59  59  59] / 256;
sol.base02  = [  7  54  66] / 256;
sol.base01  = [ 88 110 117] / 256;
sol.base00  = [101 123 131] / 256;
sol.base0   = [128 255 128] / 256;
sol.base1   = [147 161 161] / 256;
sol.base2   = [238 232 213] / 256;
sol.base3   = [253 246 227] / 256;
sol.yellow  = [181 137   0] / 256;
sol.orange  = [255,128,  0] / 256;
sol.red     = [220  50  47] / 256;
sol.magenta = [211  54 130] / 256;
sol.violet  = [108 113 196] / 256;
sol.blue    = [  0 128 255] / 256;
sol.cyan    = [ 42 161 152] / 256;
sol.green   = [133 153   0] / 256;
sol.gray    = [160,160,160] / 256;

%% Setup Color Scheme for Light/Dark Options
    sycb = 0;           % Set 'Use system colors' checkbox to False
    txc  = sol.base0;   % Set 'Text' color
    bgc  = sol.base03;  % Set 'Background' color
    kwd  = sol.orange;  % Set 'Keywords' color
    cmt  = sol.gray;    % Set 'Comments' color - spec calls for base1, not sure about that...
    str  = sol.blue;    % Set 'Strings' color
    ustr = sol.orange;  % Set 'Unterminated strings' color
    scmd = sol.yellow;  % Set 'System commands' color
    errs = sol.red;     % Set 'Errors' color
    hyp  = sol.blue;    % Set 'Hyperlinks' color
    warn = sol.orange;  % Set 'Warnings' color
    afhb = 0;           % Set 'Autofix highlight' checkbox to False
    afh  = sol.base00;  % Set 'Autofix highlight' color
%     ahib = 0;           % Set 'Automatically highlight' checkbox to False
%     ahi  = sol.violet;  % Set 'Automatically highlight' color
    vwss = sol.green;   % Set 'Variables with shared scope' color
    clhb = 0;           % Set 'Highlight cells' checkbox to False
    hclb = 1;           % Set 'Highlight current line' checkbox to True
    hcl  = sol.base02;  % Set 'Highlight current line' color
    slnb = 1;           % Set 'Show line numbers' checkbox to True
    rtlb = 1;           % Set 'Show line' checkbox in Right-hand text limit to True

%% Reset colors
%     sycb = 1;           % Set 'Use system colors' checkbox to True
% Need to get colors setup to return to default

%% Desktop tool colors
com.mathworks.services.Prefs.setBooleanPref('ColorsUseSystem',sycb); clear('sycb')

com.mathworks.services.Prefs.setColorPref('ColorsText',java.awt.Color(txc(1), txc(2), txc(3)));
com.mathworks.services.ColorPrefs.notifyColorListeners('ColorsText'); clear('txc')

com.mathworks.services.Prefs.setColorPref('ColorsBackground',java.awt.Color(bgc(1), bgc(2), bgc(3)));
com.mathworks.services.ColorPrefs.notifyColorListeners('ColorsBackground'); clear('bgc')

%% Setup MATLAB syntax highlighting colors
com.mathworks.services.Prefs.setColorPref('Colors_M_Keywords',java.awt.Color(kwd(1), kwd(2), kwd(3)));
com.mathworks.services.ColorPrefs.notifyColorListeners('Colors_M_Keywords'); clear('kwd')

com.mathworks.services.Prefs.setColorPref('Colors_M_Comments',java.awt.Color(cmt(1), cmt(2), cmt(3)));
com.mathworks.services.ColorPrefs.notifyColorListeners('Colors_M_Comments'); clear('cmt')

com.mathworks.services.Prefs.setColorPref('Colors_M_Strings',java.awt.Color(str(1), str(2), str(3)));
com.mathworks.services.ColorPrefs.notifyColorListeners('Colors_M_Strings'); clear('str')

com.mathworks.services.Prefs.setColorPref('Colors_M_UnterminatedStrings',java.awt.Color(ustr(1), ustr(2), ustr(3)));
com.mathworks.services.ColorPrefs.notifyColorListeners('Colors_M_UnterminatedStrings'); clear('ustr')

com.mathworks.services.Prefs.setColorPref('Colors_M_SystemCommands',java.awt.Color(scmd(1), scmd(2), scmd(3)));
com.mathworks.services.ColorPrefs.notifyColorListeners('Colors_M_SystemCommands'); clear('scmd')

com.mathworks.services.Prefs.setColorPref('Colors_M_Errors',java.awt.Color(errs(1), errs(2), errs(3)));
com.mathworks.services.ColorPrefs.notifyColorListeners('Colors_M_Errors');clear('errs')

%% Other colors
com.mathworks.services.Prefs.setColorPref('Colors_HTML_HTMLLinks',java.awt.Color(hyp(1), hyp(2), hyp(3)));
com.mathworks.services.ColorPrefs.notifyColorListeners('Colors_HTML_HTMLLinks'); clear('hyp')

%% Code analyzer colors
com.mathworks.services.Prefs.setColorPref('Colors_M_Warnings',java.awt.Color(warn(1), warn(2), warn(3)));
com.mathworks.services.ColorPrefs.notifyColorListeners('Colors_M_Warnings'); clear('warn')

com.mathworks.services.Prefs.setBooleanPref('ColorsUseMLintAutoFixBackground',afhb);
com.mathworks.services.Prefs.setColorPref('ColorsMLintAutoFixBackground',java.awt.Color(afh(1), afh(2), afh(3)));
com.mathworks.services.ColorPrefs.notifyColorListeners('ColorsMLintAutoFixBackground'); clear('afhb','afh')

%% Variable and function colors
% com.mathworks.services.Prefs.setBooleanPref('Editor.VariableHighlighting.Automatic',ahib);
% com.mathworks.services.Prefs.setColorPref('Editor.VariableHighlighting.Color',java.awt.Color(ahi(1), ahi(2), ahi(3)));
% com.mathworks.services.ColorPrefs.notifyColorListeners('Editor.VariableHighlighting.Color'); clear('ahib','ahi')

% *** need to add checkbox option for variables with shared scope here
com.mathworks.services.Prefs.setColorPref('Editor.NonlocalVariableHighlighting.TextColor',java.awt.Color(vwss(1), vwss(2), vwss(3)));
com.mathworks.services.ColorPrefs.notifyColorListeners('Editor.NonlocalVariableHighlighting.TextColor');clear('vwss')

%% Cell display options
com.mathworks.services.Prefs.setBooleanPref('EditorCodepadHighVisible',clhb); clear('clhb')
% Highlight cells color is 'Editorhighlight-lines'

%% Editor/Debugger General display options
com.mathworks.services.Prefs.setBooleanPref('Editorhighlight-caret-row-boolean',hclb);
com.mathworks.services.Prefs.setColorPref('Editorhighlight-caret-row-boolean-color',java.awt.Color(hcl(1), hcl(2), hcl(3)));
com.mathworks.services.ColorPrefs.notifyColorListeners('Editorhighlight-caret-row-boolean-color'); clear('hclb','hcl')
com.mathworks.services.Prefs.setBooleanPref('EditorShowLineNumbers',slnb); clear('slnb')
com.mathworks.services.Prefs.setBooleanPref('EditorRightTextLineVisible',rtlb);
%EditorRightTextLineLimit=I80
%EditorRightTextLimitLineWidth=I1

%% Cleanup
clear('sol'); clear('rtlb')


text:
function entry = findindict(c,dcode)
% Look up an entry in the dictionary
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id: findindict.m 1143 2008-02-07 19:33:33Z spm $

entry = [];
d = getdict;
d = d.(dcode);
if ischar(c)
    for i=1:length(d),
        if strcmpi(d(i).label,c),
            entry = d(i);
            break;
        end;
    end;
elseif isnumeric(c) && numel(c)==1
    for i=1:length(d),
        if d(i).code==c,
            entry = d(i);
            break;
        end;
    end;
else
    error(['Inappropriate code for ' dcode '.']);
end;
if isempty(entry)
    fprintf('\nThis is not an option.  Try one of these:\n');
    for i=1:length(d)
        fprintf('%5d) %s\n', d(i).code, d(i).label);
    end;
    %fprintf('\nNO CHANGES MADE\n');
end;

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
%% Search the peaks of a spherical function represented by SH coefficients
%% using gradient-based numerical optimization technique (Quasi-Newton)
%%
%%
%% Example:
%%   [angles] = find_peaks(coeff,init_pos);
%%  where 'init_pos' is an nx2 array giving n initial directions in
%%  pair of [polar, azimuth], e.g. [80 10; 80 110; 80 150]

function [angles] = find_peaks(coeff, init_pos)
%%=============================================================
%% Project:   Spherical Harmonics
%% Module:    $RCSfile: find_peaks.m,v $
%% Language:  MATLAB
%% Author:    $Author: bjian $
%% Date:      $Date: 2007/12/27 06:23:35 $
%% Version:   $Revision: 1.3 $
%%=============================================================

n_try = size(init_pos,1);
[n,k] = size(coeff);
angles = zeros(n, n_try*2);
%assume dl=2;
degree = (sqrt((8*k+1))-3)/2;

options = optimset('GradObj','on','Display', 'off','LargeScale','off', 'MaxIter', 100,  'TolFun',1e-06, 'TolX',1e-08);
tic
for i=1:n
     %disp(sprintf(' i = %d',i));
     angles_i = zeros(1, n_try*2);
     for ii=1:n_try
			 [peak_angle] = fminunc(@(x)real_spherical_harmonics(x,coeff(i,:), degree, 2), init_pos(ii,1:2)*pi/180,  options);
			 angles_i(2*ii-1:2*ii) = peak_angle;
     end
     %angles_i*180/pi
     angles(i,:) = angles_i;
end
toc
angles = angles*180/pi;


text:
function F = ge(X,Y)
% Internal class for constraint lists

try
    % Try to evaluate
    if isa(X,'constraint')
        % (z > w) > y
        try
            Z = X.List{end} - Y;
        catch
            Y = reshape(Y,[],1);
            Z = X.List{end} - Y;
        end
        F = X;
        F.List{end+1} = '>=';
        F.List{end+1} = Y;
        F.Evaluated{end+1} = Z;
        F.ConstraintID(end+1) = yalmip('ConstraintID');
        F.strict(end+1) = 0;
    else
        % z > (w > y)
        try
            Z = X - Y.List{1};
        catch
            X = reshape(X,[],1);
            Z = X - Y.List{1};
        end
        F = Y;
        F.List = {X,'>=',F.List{:}};
        F.Evaluated = {Z,F.Evaluated{:}};
        F.ConstraintID = [yalmip('ConstraintID') F.ConstraintID];
        F.strict = [1 F.strict];
    end
catch
    error(lasterr);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [reducedY, cl, cluster_indx] = bayes_meta_feature_abstract(Y, thresh, shrink, bayes_model, volInfo, doplot)
% Purpose: Instead of original matrix of activations in studies x voxels,
% we may want to work with a data matrix of activations in contiguous
% regions x voxels
%
% This function calculates Y = 0,1 for Y = 1 in any voxel in each
% contiguous cluster
%
% :Examples:
% ::
%
%    reducedY = bayes_meta_feature_abstract(Y, .1, bayes_model, volInfo);
%    xval = classify_naive_bayes('xval', reducedY, Xi, 0, 0, bestk, bestg );
%    xval.prop_correct_by_class
%
% needs bayes_model.pa1_given_t, nclasses, params.k

fprintf(1, 'Getting likelihood. ');

nfeatures_orig = length(bayes_model.whkeep);

pa1_given_t = zeros(nfeatures_orig, bayes_model.nclasses);
pa1_given_t(bayes_model.whkeep,:) = bayes_model.pa1_given_t;

evidence_yes = (sum(Y(:, bayes_model.whkeep_from_notempty)) ./ bayes_model.nobs)';
evidence_yes = evidence_yes(:, ones(1, bayes_model.nclasses));
tmp = zeros(nfeatures_orig, bayes_model.nclasses);
tmp(bayes_model.whkeep,:) = evidence_yes;
evidence_yes = tmp;

lr = ( (pa1_given_t + shrink) ./ (evidence_yes + shrink) ) - 1;

% % pt = zeros(nfeatures_orig, bayes_model.nclasses);
% % pt(bayes_model.whkeep,:) = bayes_model.pt_given_act1;

lr(~bayes_model.whkeep, :) = 0;

if doplot, create_figure('Hist', 1, 2); hist(lr); plot_vertical_line(thresh); plot_vertical_line(-thresh); end

lr(lr < thresh & lr > -thresh) = 0;

if doplot, subplot(1, 2, 2);
    tmp = lr(:);
    tmp(tmp == 0) = [];
    hist(tmp, 100); plot_vertical_line(thresh); plot_vertical_line(-thresh);
    sum(lr > 0)
    sum(lr < 0)
end

% -------------------------------------------------------
fprintf(1, 'Getting clusters. ');

% get all clusters
% % % cl = iimg_indx2clusters(lr(:, 1), volInfo);
% % %
% % % for i = 2:size(lr, 2)
% % %     cl = [cl iimg_indx2clusters(lr(:, i), volInfo)];
% % % end
% % %
% % % % reduce to remove overlap
% % % clu = clusters2CLU(cl);
% % % cl = tor_extract_rois([], clu, clu);
% this code does the same as the code above:
[cl, cluster_indx] = iimg_indx2clusters(any(lr')', volInfo);

fprintf(1, '%3.0f regions. ', length(cl))
%
%cl = classify_naive_bayes('plot', bayes_model, 'lr', thresh, {[1 .7 0] [0 0 1]});

fprintf(1, 'Getting data in clusters. \n');

reducedY = Meta_cluster_tools('getdata',cl, Y',volInfo);

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


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.

function price = Price_KnockIn(S,K,U,C,r,T)

minS = min(S(:,2:end),[],2);
Ind = find(minS<U);
payoff = zeros(size(S,1),1);

if C == 1
    %payoff(Ind) = rebate ;
    payoff(Ind) = max(S(Ind,end) - K,0);
else
    %payoff(Ind) = rebate;
    payoff(Ind) = max(K-S(Ind,end) ,0);
end
    price = mean(payoff) * exp(-r*T);
end