text:
%% SectoralCylinder
% Concrete subclass of <GenericCylinder.html |GenericCylinder|> representing a
% cylinder with a sectoral cross section.

%%% Description
% |SectoralCylinder| represents the shape of a sectoral cylinder.  Its cross
% section is a circular sector (angular portion of a disk).  The axis of the
% cylinder should be aligned with one of the axes of the Cartesian coordinate
% system.

%%% Construction
%  shape = SectoralCylinder(normal_axis, height, center, radius, theta, d_theta)
%  shape = SectoralCylinder(normal_axis, height, center, radius, theta, d_theta, dl_max)
%
% *Input Arguments*
%
% * |normal_axis|: axis of the cylinder.  It should be one of |Axis.x|,
% |Axis.y|, |Axis.z|.
% * |height|: size of the cylinder along its axis.
% * |center|: center of the cylinder in the format of |[x y z]|.  For
% |normal_axis = Axis.z|, |(x, y)| is the coordinate of the center of the
% circle.
% * |radius|: radius of the circle
% * |theta|: beginning angle of the sector in radian
% * |d_theta|: angular width of the sector in radian between -2*pi and 2*pi.
% * |dl_max|: maximum grid size allowed in the cylinder.  It can be either |[dx
% dy dz]| or a single real number |dl| for |dx = dy = dz|.  If unassigned,
% |dl_max = Inf| is used.

%%% Example
%   % Create an instance of SectoralCylinder.
%   shape = SectoralCylinder(Axis.z, 100, [0 0 50], 50, pi/6, pi/3);
%
%   % Use the constructed shape in maxwell_run().
%   [E, H] = maxwell_run({INITIAL ARGUMENTS}, 'OBJ', {'vacuum', 'none', 1.0}, shape, {REMAINING ARGUMENTS});

%%% See Also
% <SectoralCylinder.html |CircularCylinder|>, <CircularCylinder.html
% |CircularCylinder|>, <CircularShellCylinder.html
% |CircularShellCylinder|>, <EllipticCylinder.html |EllipticCylinder|>,
% <PolyognalCylinder.html |PolygonalCylinder|>, <Shape.html |Shape|>,
% <maxwell_run.html |maxwell_run|>

classdef SectoralCylinder < GenericCylinder

	properties (SetAccess = immutable)
		lsf_th  % level set function for angles
	end

	methods
        function this = SectoralCylinder(normal_axis, height, center, radius, theta, d_theta, dl_max)
			chkarg(istypesizeof(normal_axis, 'Axis'), '"normal_axis" should be instance of Axis.');
			chkarg(istypesizeof(height, 'real') && height > 0, '"height" should be positive.');
			chkarg(istypesizeof(center, 'real', [1, Axis.count]), ...
				'"center" should be length-%d row vector with real elements.', Axis.count);
			chkarg(istypesizeof(radius, 'real') && radius > 0, '"radius" should be positive.');
			chkarg(istypesizeof(theta, 'real'), '"theta" should be real.');
			chkarg(istypesizeof(d_theta, 'real') && (d_theta <= 2*pi || d_theta >= -2*pi), '"d_theta" should be real between -pi and pi.');

			thetas = sort([theta, theta + d_theta]);
			cth = mean(thetas);  % center of theta range
			cth = mod(cth + pi, 2*pi) - pi;  % -pi <= c_th < pi (range of atan2(y,x))
			sth = diff(thetas) / 2;  % semiwidth of theta range

			cr = radius/2;  % center of radius range
			sr = radius/2;  % semiwidth of radius range

			function level = lsf_th(th)
				dth = th - cth;
				dth = mod(dth + pi, 2*pi) - pi;
				level = 1 - abs(dth./sth);
			end

			[h, v, n] = cycle(normal_axis);

			% lsf2d() can handle rho = [p q] with column vectors p and q.  The
			% level set function is the one for a rectangle defined in the
			% (theta, radius) domain.
			function level = lsf2d(p, q)
				chkarg(istypeof(p, 'real'), '"p" should be array with real elements.');
				chkarg(istypeof(q, 'real'), '"q" should be array with real elements.');
				chkarg(isequal(size(p), size(q)), '"p" and "q" should have same size.');

				lc = {p - center(h), q - center(v)};  % locations in center-of-mass coordinates

				theta_pt = atan2(lc{Dir.v}, lc{Dir.h});  % -pi <= theta_rho < pi
				r_pt = sqrt(lc{Dir.h}.^2 + lc{Dir.v}.^2);
				zero_r_pt = (r_pt==0);  % atan2(0,0) is not well-defined; handle such cases separately

				dr = r_pt - cr;
				level = min(lsf_th(theta_pt), 1 - abs(dr./sr));
				level(zero_r_pt) = 0;
			end

			lprim = cell(1, Axis.count);
			lprim{h} = [radius * cos(thetas) + center(h), center(h)];
			lprim{v} = [radius * sin(thetas) + center(v), center(v)];
			lprim{n} = [-height height]/2 + center(n);

			if lsf_th(0) > 0  % sector contains +x-direction from center
				lprim{h} = [lprim{h}, center(h) + radius];
			end
			if lsf_th(pi/2) > 0  % sector contains +y-direction from center
				lprim{v} = [lprim{v}, center(v) + radius];
			end
			if lsf_th(pi) > 0  % sector contains -x-direction from center
				lprim{h} = [lprim{h}, center(h) - radius];
			end
			if lsf_th(3*pi/2) > 0  % sector contains -y-direction from center
				lprim{v} = [lprim{v}, center(v) - radius];
			end

			if nargin < 7  % no dl_max
				super_args = {normal_axis, @lsf2d, lprim};
			else
				super_args = {normal_axis, @lsf2d, lprim, dl_max};
			end

			this = this@GenericCylinder(super_args{:});
			this.lsf_th = @lsf_th;
		end
	end
end

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
classdef AeroAnglesPolySteeringModel < AbstractAnglePolySteeringModel
    %AeroAnglesPolySteeringModel Summary of this class goes here
    %   Detailed explanation goes here

    properties
        bankModel(1,1) PolynominalModel = PolynominalModel(0,0,0,0);
        aoAModel(1,1) PolynominalModel = PolynominalModel(0,0,0,0);
        slipModel(1,1) PolynominalModel = PolynominalModel(0,0,0,0);

        bankContinuity(1,1) logical = true;
        aoAContinuity(1,1) logical = true;
        slipContinuity(1,1) logical = true;
    end

    methods
        function R_body_2_inertial = getBody2InertialDcmAtTime(obj, ut, rVect, vVect, bodyInfo)
            bankAng = obj.bankModel.getValueAtTime(ut);
            angOfAttack = obj.aoAModel.getValueAtTime(ut);
            angOfSideslip = obj.slipModel.getValueAtTime(ut);

%             [~, ~, ~, dcm] = computeBodyAxesFromAeroAngles(ut, rVect, vVect, bodyInfo, angOfAttack, angOfSideslip, bankAng);
            baseFrame = bodyInfo.getBodyFixedFrame();
            [~, ~, ~, R_body_2_inertial] = computeInertialBodyAxesFromFrameAeroAngles(ut, rVect, vVect, bodyInfo, bankAng, angOfAttack, angOfSideslip, baseFrame);
        end

        function [angleModel, continuity] = getAngleNModel(obj, n)
            angleModel = PolynominalModel.empty(1,0);

            switch n
                case 1
                    angleModel = obj.bankModel;
                    continuity = obj.bankContinuity;
                case 2
                    angleModel = obj.aoAModel;
                    continuity = obj.aoAContinuity;
                case 3
                    angleModel = obj.slipModel;
                    continuity = obj.slipContinuity;
            end
        end

        function t0 = getT0(obj)
            t0 = obj.slipModel.t0;
        end

        function setT0(obj, newT0)
            obj.bankModel.t0 = newT0;
            obj.aoAModel.t0 = newT0;
            obj.slipModel.t0 = newT0;
        end

        function setConstTerms(obj, bankConst, aoaConst, slipConst)
            obj.bankModel.constTerm = bankConst;
            obj.aoAModel.constTerm = aoaConst;
            obj.slipModel.constTerm = slipConst;
        end

        function setLinearTerms(obj, bank, aoa, slip)
            obj.bankModel.linearTerm = bank;
            obj.aoAModel.linearTerm = aoa;
            obj.slipModel.linearTerm = slip;
        end

        function setAccelTerms(obj, bank, aoa, slip)
            obj.bankModel.accelTerm = bank;
            obj.aoAModel.accelTerm = aoa;
            obj.slipModel.accelTerm = slip;
        end

        function setTimeOffsets(obj, timeOffset)
            obj.bankModel.tOffset = timeOffset;
            obj.aoAModel.tOffset = timeOffset;
            obj.slipModel.tOffset = timeOffset;
        end

        function [angle1Cont, angle2Cont, angle3Cont] = getContinuityTerms(obj)
            angle1Cont = obj.bankContinuity;
            angle2Cont = obj.aoAContinuity;
            angle3Cont = obj.slipContinuity;
        end

        function setContinuityTerms(obj, angle1Cont, angle2Cont, angle3Cont)
            obj.bankContinuity = angle1Cont;
            obj.aoAContinuity = angle2Cont;
            obj.slipContinuity = angle3Cont;
        end

        function setConstsFromDcmAndContinuitySettings(obj, dcm, ut, rVect, vVect, bodyInfo)
            if(obj.bankContinuity || obj.aoAContinuity || obj.slipContinuity)
                [bankAng,angOfAttack,angOfSideslip] = computeAeroAnglesFromBodyAxes(ut, rVect, vVect, bodyInfo, dcm(:,1), dcm(:,2), dcm(:,3));

                if(obj.bankContinuity)
                    obj.bankModel.constTerm = bankAng;
                end

                if(obj.aoAContinuity)
                    obj.aoAModel.constTerm = angOfAttack;
                end

                if(obj.slipContinuity)
                    obj.slipModel.constTerm = angOfSideslip;
                end
            end
        end

        function setInitialAttitudeFromState(obj, stateLogEntry, tOffsetDelta)
            t0 = stateLogEntry.time;
            obj.setT0(t0);

            obj.bankModel.tOffset = obj.bankModel.tOffset + tOffsetDelta;
            obj.aoAModel.tOffset = obj.aoAModel.tOffset + tOffsetDelta;
            obj.slipModel.tOffset = obj.slipModel.tOffset + tOffsetDelta;
        end

        function [angle1Name, angle2Name, angle3Name] = getAngleNames(~)
            angle1Name = 'Bank Angle';
            angle2Name = 'Angle of Attack';
            angle3Name = 'Side Slip Angle';
        end

        function newSteeringModel = deepCopy(obj)
            newSteeringModel = AeroAnglesPolySteeringModel(obj.bankModel.deepCopy(), obj.aoAModel.deepCopy(), obj.slipModel.deepCopy());
            newSteeringModel.bankContinuity = obj.bankContinuity;
            newSteeringModel.aoAContinuity = obj.aoAContinuity;
            newSteeringModel.slipContinuity = obj.slipContinuity;
        end

        function optVar = getNewOptVar(obj)
            optVar = SetAeroSteeringModelActionOptimVar(obj);
        end

        function optVar = getExistingOptVar(obj)
            optVar = obj.optVar;
        end

        function [addActionTf, steeringModel] = openEditSteeringModelUI(obj, lv, useContinuity)
            output = AppDesignerGUIOutput({false, obj});
            lvd_EditActionSetSteeringModelGUI_App(obj, lv, useContinuity, output);
            addActionTf = output.output{1};
            steeringModel = output.output{2};
        end
    end

    methods(Access=private)
        function obj = AeroAnglesPolySteeringModel(bankModel, aoAModel, slipModel)
            obj.bankModel = bankModel;
            obj.aoAModel = aoAModel;
            obj.slipModel = slipModel;
        end
    end

    methods(Static)
        function model = getDefaultSteeringModel()
            bankModel = PolynominalModel(0,0,0,0);
            aoAModel = PolynominalModel(0,0,0,0);
            slipModel = PolynominalModel(0,0,0,0);

            model = AeroAnglesPolySteeringModel(bankModel, aoAModel, slipModel);
        end

        function typeStr = getTypeNameStr()
            typeStr = SteeringModelEnum.AeroAnglesPoly.nameStr;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function coeff = fourier_comp(func,m,L)
%FOURIER_COMP Numerically evaluate Fourier series coefficient.
% coeff = fourier_comp(func,m,L) tries to approximate the function func
% from -L to L with m term Fourier series using quad (MATLAB functions).
% func is a function handle, should accept a vector argument x and return a
% vector result y
% This function return a structure with these fields:
% coeff.a0
% coeff.an
% coeff.bn
% The Fourier series for function f(x) is given below.
%                       inf
%                      -----
%                      \
%        f(x) = a0/2 +  \  an * sin(n*pi*x/L) + bn * con(n*pi*x/L)
%                       /
%                      /
%                      -----
%                      n = 1
%
%   Example:
%
%   f = @(x)x.*cos(x);
%   coeff = fourier_comp(f,3,pi)
%     coeff =
%     a0: 0
%     an: [-0.5000 1.3333 -0.7500]
%     bn: [1.4136e-016 0 7.0679e-017]
%
%   See also
%   fourier_gui
%
% Author: Amin Bashi
% Created: Jul 2009
% Copyright 2009

% func must be the function handle

for n = 1:m
    an(n) = quad(@fa,-L,L)/L;
    bn(n) = quad(@fb,-L,L)/L;
end
a0 = quad(func,-L,L)/L;
coeff.a0 = a0;
coeff.an = an;
coeff.bn = bn;
    function y = fa(t)
        y = func(t).*sin(n*pi*t/L);
    end
    function y = fb(t)
        y = func(t).*cos(n*pi*t/L);
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function set_device(device_id)
% set_device(device_id)
%   set Caffe's GPU device ID

CHECK(isscalar(device_id) && device_id >= 0, ...
  'device_id must be non-negative integer');
device_id = double(device_id);

caffe_('set_device', device_id);

end

author: aimerykong

repo: deepImageAestheticsAnalysis


text:
function [output,status] = urlwrite(urlChar,location,method,params)
%URLWRITE Save the contents of a URL to a file.
%   URLWRITE(URL,FILENAME) saves the contents of a URL to a file.  FILENAME
%   can specify the complete path to a file.  If it is just the name, it will
%   be created in the current directory.
%
%   F = URLWRITE(...) returns the path to the file.
%
%   F = URLWRITE(...,METHOD,PARAMS) passes information to the server as
%   part of the request.  The 'method' can be 'get', or 'post' and PARAMS is a
%   cell array of param/value pairs.
%
%   [F,STATUS] = URLWRITE(...) catches any errors and returns the error code.
%
%   Examples:
%   urlwrite('http://www.mathworks.com/',[tempname '.html'])
%   urlwrite('ftp://ftp.mathworks.com/README','readme.txt')
%   urlwrite(['file:///' fullfile(prefdir,'history.m')],'myhistory.m')
%
%   From behind a firewall, use the Preferences to set your proxy server.
%
%   See also URLREAD.

%   Matthew J. Simoneau, 13-Nov-2001
%   Copyright 1984-2011 The MathWorks, Inc.
%   $Revision: 1.4.4.14 $ $Date: 2011/09/03 22:43:01 $

% This function requires Java.
if ~usejava('jvm')
   error(message('MATLAB:urlwrite:NoJvm'));
end

import com.mathworks.mlwidgets.io.InterruptibleStreamCopier;

% Be sure the proxy settings are set.
com.mathworks.mlwidgets.html.HTMLPrefs.setProxySettings

% Check number of inputs and outputs.
if ~ischar(urlChar)
    error('MATLAB:urlwrite:InvalidInput','The first input, the URL, must be a character array.');
end
if ~ischar(location)
    error('MATLAB:urlwrite:InvalidInput','The second input, a filename, must be a character array.');
end
if (nargin > 2) && ~strcmpi(method,'get') && ~strcmpi(method,'post')
    error('MATLAB:urlwrite:InvalidInput','Second argument must be either "get" or "post".');
end

% Do we want to throw errors or catch them?
if nargout == 2
    catchErrors = true;
else
    catchErrors = false;
end

% Set default outputs.
output = '';
status = 0;

% GET method.  Tack param/value to end of URL.
if (nargin > 2) && strcmpi(method,'get')
    if mod(length(params),2) == 1
        error('MATLAB:urlwrite:InvalidInput','Invalid parameter/value pair arguments.');
    end
    for i=1:2:length(params)
        if (i == 1), separator = '?'; else separator = '&'; end
        param = char(java.net.URLEncoder.encode(params{i}));
        value = char(java.net.URLEncoder.encode(params{i+1}));
        urlChar = [urlChar separator param '=' value];
    end
end

% Create a urlConnection.
[urlConnection,errorid,errormsg] = plugin_urlreadwrite(mfilename,urlChar);
if isempty(urlConnection)
    if catchErrors, return
    else error(errorid,errormsg);
    end
end

pluginSize = 2*plugin_urlsize(urlChar)/1000000;
timeOut = max(round(pluginSize*1000), 5000);
urlConnection.setReadTimeout(timeOut); % timeout in 5 seconds

% POST method.  Write param/values to server.
if (nargin > 2) && strcmpi(method,'post')
    try
        urlConnection.setDoOutput(true);
        urlConnection.setRequestProperty( ...
            'Content-Type','application/x-www-form-urlencoded');
        printStream = java.io.PrintStream(urlConnection.getOutputStream);
        for i=1:2:length(params)
            if (i > 1), printStream.print('&'); end
            param = char(java.net.URLEncoder.encode(params{i}));
            value = char(java.net.URLEncoder.encode(params{i+1}));
            printStream.print([param '=' value]);
        end
        printStream.close;
    catch
        if catchErrors, return
        else error('MATLAB:urlwrite:ConnectionFailed','Could not POST to URL.');
        end
    end
end

% Specify the full path to the file so that getAbsolutePath will work when the
% current directory is not the startup directory and urlwrite is given a
% relative path.
file = java.io.File(location);
if ~file.isAbsolute
   location = fullfile(pwd,location);
   file = java.io.File(location);
end

% Make sure the path isn't nonsense.
try
   file = file.getCanonicalFile;
catch
   error('MATLAB:urlwrite:InvalidOutputLocation','Could not resolve file "%s".',char(file.getAbsolutePath));
end

% Open the output file.
try
    fileOutputStream = java.io.FileOutputStream(file);
catch
    error('MATLAB:urlwrite:InvalidOutputLocation','Could not open output file "%s".',char(file.getAbsolutePath));
end

% Read the data from the connection.
try
    inputStream = urlConnection.getInputStream;
    % This StreamCopier is unsupported and may change at any time.
    isc = InterruptibleStreamCopier.getInterruptibleStreamCopier;
    isc.copyStream(inputStream,fileOutputStream);
    inputStream.close;
    fileOutputStream.close;
    output = char(file.getAbsolutePath);
catch
    fileOutputStream.close;
    delete(file);
    if catchErrors, return
    else error('MATLAB:urlwrite:ConnectionFailed','Error downloading URL. Your network connection may be down or your proxy settings improperly configured.');
    end
end

status = 1;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% If all is well, all of these scripts should run without errors.

% bnets
cg1
cg2
discrete1
fa1
gaussian1
gaussian2
if exist('@gibbs_sampling_inf_engine/private/compute_posterior','file')
  % only exists if installC has been run
  gibbs_test1
end
learn1
lw1
mfa1
mixexp1
mixexp2
mixexp3
mog1
mpe1
mpe2
qmr1
qmr2
sample1
softev1
softmax1
sprinkler1

% belief propagation
belprop_polytree_discrete
belprop_polytree_gauss % alag
belprop_polytree_cg
belprop_loop1_discrete
belprop_loop1_gauss
belprop_loopy_discrete
belprop_loopy_gauss
belprop_loopy_cg % like cg1

% factor graphs
%fg1   failed since marginals were not exact

fg2
fg3
fg_mrf1
fg_mrf2

% Structure learning
bic1
cooper_yoo
k2demo1
mcmc1
model_select1
pc1
%pc2   failed due to numerical problems in KPMstats/cond_indep_fisher_z

% limids
asia_dt1
id1
oil1
pigs1

% dbns
arhmm1
bat1
bkff1
chmm1
dhmm1
filter_test1
ghmm1
kalman1
kjaerulff1
loopy_dbn1
mhmm1
mildew1
reveal1
viterbi1
water1

% HHMMs
abcd_hhmm
sample_square_hhmm_discrete
%learn_square_hhmm_cts
sample_motif_hhmm

%sparse jtree engine & ndx 2TBN engine
if exist('@jtree_sparse_inf_engine/init_pot','file')
  % only exists if installC has been run
  discrete2
  discrete3
  filter_test1
  water2
end

%find . -path '*.m' -exec wc -l {} \; | ~/count.pl

% we cannot use tic;toc to time test_BNT, since functions within this script
% reset the tic;toc timer. Hence we use the following:
%clock0=clock; cpu0 = cputime; test_BNT; cpu=cputime-cpu0; elapsed=etime(clock, clock0)

author: beckel

repo: nilm-eval

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


text:
function varargout = lap(varargin)
%LAP   Laplacian of a DISKFUN.
%   L = LAP(F) returns a DISKFUN representing the Laplacian of F.
%
%   This is shorthand for LAPLACIAN(F).
%
% See also DISKFUN/LAPLACIAN.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

[varargout{1:nargout}] = laplacian(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%% FOURCURL3P2TEST
%
% Use the quadraticNedelec elements to approximate the velocity u and
% the stream function w.
%
% We solve the following equations:
%
%  -w + curl curl u = 0
%  curl curl w + u  = f
%
% with Dirichlet boundary condition
%
% u\times n = (curl u )\times n  = 0  on \partial \Omega
%
% Please check fourCurl3doc for details.
%
% Lin Zhong, May, 2013. Clean up. Long Chen, Oct 21, 2018
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

close all;
clear variables;

%% Problem Setting
pde = fourCurl3data;
[node,elem] = cubemesh([-1,1,-1,1,-1,1],2);
bdFlag = setboundary3(node,elem,'Dirichlet');

%% Parameters
maxIt = 3;
errwIwhL2 = zeros(maxIt,1);
errwIwhHcurl = zeros(maxIt,1);
errwL2 = zeros(maxIt,1);
errwwhHcurl = zeros(maxIt,1);
erruIuhL2 = zeros(maxIt,1);
erruIuhHcurl = zeros(maxIt,1);
erruL2 = zeros(maxIt,1);
erruuhHcurl = zeros(maxIt,1);
assembleTime = zeros(maxIt,1);
solverTime = zeros(maxIt,1);
N = zeros(maxIt,1);
h = zeros(maxIt,1);
option = [];

for k = 1:maxIt
    % refine grid
    [node,elem,bdFlag] = uniformrefine3(node,elem,bdFlag);
    % solve the equation
    [w,u,eqn,info] = fourCurl3P2(node,elem,bdFlag,pde,option);
    % compute error
    uI = edgeinterpolate2(pde.exactu,node,eqn.edge,eqn.face,eqn.face2edge);
    wI = edgeinterpolate2(pde.curlcurlu,node,eqn.edge,eqn.face,eqn.face2edge);
    errwL2(k) = getL2error3ND2(node,elem,pde.curlcurlu,w);
    erruL2(k) = getL2error3ND2(node,elem,pde.exactu,u);
    errwwhHcurl(k) = getHcurlerror3ND2(node,elem,pde.curlcurlcurlu,w);
    erruuhHcurl(k) = getHcurlerror3ND2(node,elem,pde.curlu,u);
    errwIwhL2(k) = sqrt((w-wI)'*eqn.M*(w-wI));
    erruIuhL2(k) = sqrt((u-uI)'*eqn.M*(u-uI));
    errwIwhHcurl(k) = sqrt((w-wI)'*eqn.A*(w-wI));
    erruIuhHcurl(k) = sqrt((u-uI)'*eqn.A*(u-uI));
    % record information
    N(k) = length(u)+length(w);
    h(k) = 1./(size(node,1)^(1/3)-1);
    assembleTime(k) = info.assembleTime;
    solverTime(k) = info.solverTime;
end

%% Plot convergence rates
subplot(2,2,1)
showrateh2(h, errwIwhL2(1:maxIt), 1 , 'k-+', '||w_I-w_h||',...
           h, erruIuhL2(1:maxIt), 1 , 'r-*', '||u_I-u_h||' );
subplot(2,2,2)
showrateh2(h, errwL2(1:maxIt), 1 , 'k-+', '||w-w_h||',...
           h, erruL2(1:maxIt), 1 , 'r-*', '||u-u_h||' );
subplot(2,2,3)
showrateh2(h, errwIwhHcurl(1:maxIt), 1 , 'b-+', '||w_I-w_h||_1',...
           h, erruIuhHcurl(1:maxIt), 1 , 'g-*', '||u_I-u_h||_1' );
subplot(2,2,4)
showrateh2(h, errwwhHcurl(1:maxIt), 1 , 'b-+', '||w-w_h||_1',...
           h, erruuhHcurl(1:maxIt), 1 , 'g-*', '||u-u_h||_1' );

%% Output
err = struct('N', N, 'h', h, ...
             'wL2',errwL2(1:maxIt),'uL2',erruL2(1:maxIt),...
             'wIwhL2',errwIwhL2(1:maxIt), 'uIuhL2',erruIuhL2(1:maxIt),...
             'wIwhHcurl',errwIwhHcurl(1:maxIt),'uIuhHcurl',erruIuhHcurl(1:maxIt),...
             'wwhHcurl',errwwhHcurl(1:maxIt),'uuhHcurl',erruuhHcurl(1:maxIt));
time = struct('N',N,'assemble',assembleTime(1:maxIt),'solver',solverTime(1:maxIt));

%% Display error on screen
disp('Table: Error')
colname = {'#Dof','h','||w_I-w_h||','||w-w_h||','||w_I-w_h||_1','||w-w_h||_1'};
disptable(colname,err.N,[],err.h,'%0.3e',err.wIwhL2,'%0.5e',err.wL2,'%0.5e',...
              err.wIwhHcurl,'%0.5e',err.wwhHcurl,'%0.5e');

colname = {'#Dof','h','||u_I-u_h||','||u-u_h||','||u_I-u_h||_1','||u-u_h||_1'};
disptable(colname,err.N,[],err.h,'%0.3e',err.uIuhL2,'%0.5e',err.uL2,'%0.5e',...
              err.uIuhHcurl,'%0.5e',err.uuhHcurl,'%0.5e');

%     display('Table: CPU time')
colname = {'#Dof','Assemble','Solve'};
disptable(colname,time.N,[],time.assemble,'%0.2e',time.solver,'%0.2e');

repo: ifem


text:
function y=v_rhartley(x,n)
%V_RHARTLEY Calculate the Hartley transform of real data Y=(X,N)
% Data is truncated/padded to length N if specified.
% The inverse transformation is x=hartley(y,n)/n

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: v_rhartley.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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin < 2
  y=fft(real(x));
else
  y=fft(real(x),n);
end
y=real(y)-imag(y);

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function sucra()
    % This script evaluates the percentage of space time coevered by
    %alarms
    %
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals
    global iala
    re = [];
    % Stefan Wiemer    4/95

    report_this_filefun();

    abo = abo2;

    for tre2 = min(abo(:,4)):0.1:max(abo(:,4)-0.1)
        abo = abo2;
        abo(:,5) = abo(:,5)* days(ZG.bin_dur) + ZG.primeCatalog.Date(1);
        l = abo(:,4) >= tre2;
        abo = abo(l,:);
        l = abo(:,3) < ZG.tresh_km;
        abo = abo(l,:);
        set(gca,'NextPlot','add')

        % space time volume covered by alarms
        if isempty(abo)
            Va = 0;
        else
            Va = sum(pi*abo(:,3).^2)*iala;
        end

        % All space time
        [len, ncu] = size(cumuall);

        r = loc(3,:);
        %r = reshape(cumuall(len,:),length(gy),length(gx));
        %r=reshape(normlap2,length(yvect),length(xvect));
        l = r < ZG.tresh_km;
        V = sum(pi*r(l).^2*(teb-t0b));
        disp([' Zalarm = ' num2str(tre2)])
        disp([' =============================================='])
        disp([' Total space-time volume (R<Rmin):  ' num2str(V)])
        disp([' Space-time volume covered with alarms (R<Rmin):  ' num2str(Va)])
        disp([' Percent of total covered with alarms (R<Rmin):  ' num2str(Va/V*100) ' Percent' ])

        re = [re ; tre2 Va/V*100 ];
    end   % for tre2

    figure

    axis off

    uicontrol('Units','normal',...
        'Position',[.0 .65 .08 .06],'String','Save ',...
        'Callback',{@calSave9, re(:,1), re(:,2)})

    rect = [0.20,  0.10, 0.70, 0.60];
    axes('position',rect)
    set(gca,'NextPlot','add')
    pl = semilogy(re(:,1),re(:,2),'r');
    set(pl,'LineWidth',1.5)
    pl = semilogy(re(:,1),re(:,2),'ob');
    set(pl,'LineWidth',1.5,'MarkerSize',10)
    set(gca,'YScale','log')

    set(gca,'visible','on','FontSize',ZmapGlobal.Data.fontsz.m,'FontWeight','bold',...
        'FontWeight','bold','LineWidth',1.5,...
        'Box','on')
    grid

    ylabel('Va/Vtotal in %')
    xlabel('Zalarm ')
    watchoff

end

author: CelsoReyes


text:
 function y = mtimes(ob, x)
%function y = mtimes(ob, x)
%	y = G * x	or x = G' * y
%	Copyright 2002-2-20	Jeff Fessler	The University of Michigan

%
%	scalar * G
%
if isa(ob, 'double') & length(x) == 1 & isa(x, 'Gtomo2_dsc')
	y = x;
	y.scale = ob;
	return
end

if ob.apower ~= 1, error notdone, end

%
%	partial projection or backprojection (for ordered subsets)
%
if ob.is_subset
	%
	%	Gt(:,ii)' * x		partial forward projection
	%
	if ~ob.is_transpose
		if ob.is_masked
			x = embed(x, ob.mask);
		end
		y = wtfmex('dsc,proj', ob.arg', single(x), ob.mask, ...
			int32(ob.ia_start), int32(ob.ia_inc), ...
			int32(ob.nthread), ob.chat);
		y = y(:,(ob.ia_start+1):ob.ia_inc:ob.na);

	%
	%	Gt(:,ii) * y		partial back projection
	%
	else
		y = zeros(ob.nb,ob.na);
		ia = (ob.ia_start+1):ob.ia_inc:ob.na;
		y(:,ia) = reshape(x, ob.nb, length(ia));
		y = wtfmex('dsc,back', ob.arg', single(y), ob.mask, ...
			int32(ob.ia_start), int32(ob.ia_inc), ...
			int32(ob.nthread), ob.chat);
		if ob.is_masked
			y = y(ob.mask);
		end
	end

%
%	full projection
%
elseif ~ob.is_transpose
	if ob.is_masked
		x = embed(x, ob.mask);
	end
	y = wtfmex('dsc,proj', ob.arg', single(x), ob.mask, ...
		int32(0), int32(1), int32(ob.nthread), ob.chat);

%
%	full back-projection
%
else
	y = wtfmex('dsc,back', ob.arg', single(x), ob.mask, ...
		int32(0), int32(1), int32(ob.nthread), ob.chat);
	if ob.is_masked
		y = y(ob.mask);
	end
end

y = ob.scale * double(y(:));

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [l1, l2, l3]=dtiEigenvaluesFromWestinShapes(cl, cp, vol, method)
%V is the volume of original tensor
%solution
%Method specifies whether the Westin shapes aligned are computed with
%simple ("new") denominator, l1, or original definition (old) denominator,
%l1+l2+l3
if ~exist('method', 'var')
    method='westinShapes_l1';
end

switch method
    case 'westinShapes_l1'
%westin shapes are new(simplified) versions, NOT the ones computed by dtiComputeWestinShapes:
%cl=(l1-l2)/l1;
%cp=(l2-l3)/l1;
%cs=l3/l1;

l1_sol(:, 1)=-((-3/pi).^(1/3).*vol.^(1/3))./(2^(2/3).*((-1+cl).*(-1+cl+cp)).^(1/3));
l3_sol(:, 1)=(-1+cl+cp).*l1_sol(:, 1);
l1_sol(:, 2)=(3/pi)^(1/3)./(2^(2/3).*(((-1+cl).*(-1+cl+cp))./vol).^(1/3));
l3_sol(:, 2)=-(-1+cl+cp).*l1_sol(:, 2);
l1_sol(:, 3)=-((-1)^(2/3).*(3/pi)^(1/3))./(2^(2/3).*(((-1+cl).*(-1+cl+cp))./vol).^(1/3));
l3_sol(:, 3)=(-1+cl+cp).*l1_sol(:, 3);

    case 'westinShapes_lsum'
%westin shapes as those computed by dtiComputeWestinShapes:
%cl=(l1-l2)/(l1+l2+l3);
%cp=(l2-l3)/(l1+l2+l3);
%cs=l3/(l1+l2+l3);

%I am not sure whether the results produced by this method make sence -- at
%least when protted in barycentric coordinates. If you want to use it,
%check the code.

l1_sol(:, 1)=-((3/pi).^(1/3).*(-(2+4*cl+cp).^2.*vol).^(1/3))./(2*((-2+2*cl-cp).*(-1+cl+cp)).^(1/3));
l3_sol(:, 1)=-2*l1_sol(1).*(-1+cl+cp)./(2+4*cl+cp);
l1_sol(:, 2)=((3/pi).^(1/3))./(2*(((-2+2*cl-cp).*(-1+cl+cp))./((2+4*cl+cp).^2.*vol)).^(1/3));
l3_sol(:, 2)=-2*l1_sol(2).*(-1+cl+cp)./(2+4*cl+cp);
l1_sol(:, 3)=((-1).^(2/3).*(3/pi).^(1/3))./(2*(((-2+2*cl-cp).*(-1+cl+cp))./((2+4*cl+cp).^2.*vol)).^(1/3));
l3_sol(:, 3)=-2*l1_sol(3).*(-1+cl+cp)./(2+4*cl+cp);

    otherwise
        fprintf('Enter either "westinShapes_lsum" or "westinShapes_l1" for method'); return;
end

%The three solutions are only different in that some of them are not real! The second one is usually good enough.
solN=1;

while(~isreal(l1_sol(:, solN)) || ~isreal(l3_sol(:, solN)))
solN=solN+1;
end
l1=l1_sol(:, solN);
l3=l3_sol(:, solN);
l2=vol./(l1.*l3.*4.*pi./3);
l2(isnan(l2))=0;

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function row_num = file_row_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_ROW_COUNT counts the number of row records in a file.
%
%  Discussion:
%
%    Each input line is a "RECORD".
%
%    The records are divided into three groups:
%
%    * BLANK LINES (nothing but blanks)
%    * COMMENT LINES (begin with a '#')
%    * DATA RECORDS (anything else)
%
%    The value returned by the function is the number of data records.
%
%    By the way, if the MATLAB routine FGETS is used, instead of
%    FGETL, then the variable LINE will include line termination
%    characters, which means that a blank line would not actually
%    have zero characters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 December 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the input file.
%
%    Output, integer ROW_NUM, the number of rows found.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_ROW_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_ROW_COUNT - Error!' );
    row_num = -1;
    return;
  end

  blank_num = 0;
  comment_num = 0;
  row_num = 0;

  record_num = 0;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    record_num = record_num + 1;
    record_length = s_len_trim ( line );

    if ( record_length <= 0 )
      blank_num = blank_num + 1;
    elseif ( line(1) == '#' )
      comment_num = comment_num + 1;
    else
      row_num = row_num + 1;
    end

  end

  fclose ( input_unit );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function out=cellFilter(filterFunction,in)

keep=cellfun(filterFunction,in);
out=in(keep);

end

author: kristinbranson

repo: JAABA


text:
function ydot = lorenzeq(t,y)
%LORENZEQ Equation of the Lorenz chaotic attractor.
%   ydot = lorenzeq(t,y).
%   The differential equation is written in almost linear form.

SIGMA = 10.;
RHO = 28.;
BETA = 8./3.;

A = [ -BETA    0     y(2)
    0  -SIGMA   SIGMA
    -y(2)   RHO    -1  ];

ydot = A*y;


text:
%PLOT_HOMLINE Draw a line in homogeneous form
%
% PLOT_HOMLINE(L) draws a 2D line in the current plot defined in homogenous
% form ax + by + c = 0  where L (3x1) is L = [a b c].
% The current axis limits are used to determine the endpoints of
% the line.  If L (3xN) then N lines are drawn, one per column.
%
% PLOT_HOMLINE(L, LS) as above but the MATLAB line specification LS is given.
%
% H = PLOT_HOMLINE(...) as above but returns a vector of graphics handles for the lines.
%
% Notes::
% - The line(s) is added to the current plot.
% - The line(s) can be drawn in 3D axes but will always lie in the
%   xy-plane.
%
% Example::
%          L = homline([1 2]', [3 1]'); % homog line from (1,2) to (3,1)
%          plot_homline(L, 'k--'); % plot dashed black line
%
% See also PLOT_BOX, PLOT_POLY, HOMLINE.

% Copyright (C) 1993-2019 Peter I. Corke
%
% This file is part of The Spatial Math Toolbox for MATLAB (SMTB).
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
% of the Software, and to permit persons to whom the Software is furnished to do
% so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
% FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
% COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
% IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
% CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
%
% https://github.com/petercorke/spatial-math

function handles = plot_homline(lines, varargin)

	% get plot limits from current graph
	xlim = get(gca, 'XLim');
	ylim = get(gca, 'YLim');

    ish = ishold;
    hold on;

    if min(size(lines)) == 1
        lines = lines(:);
    end

    assert(numrows(lines) == 3, 'SMTB:plot_homline:badarg', 'Input must be a 3-vector or 3xN matrix');

	h = [];
	% for all input lines (columns
	for l=lines
        if abs(l(2)) > abs(l(1))
            y = (-l(3) - l(1)*xlim) / l(2);
            hh = plot(xlim, y, varargin{:});
        else
            x = (-l(3) - l(2)*ylim) / l(1);
            hh = plot(x, ylim, varargin{:});
        end
		h = [h; hh];
	end

    if ~ish
        hold off
    end

	if nargout > 0,
		handles = h;
	end

author: petercorke

repo: spatialmath-matlab


text:
function [x, funVal, ValueL]=LeastR(A, y, z, opts)
%
%%
% Function LeastR
%      Least Squares Loss with the L1-norm Regularization
%
%% Problem
%
%  min  1/2 || A x - y||^2 + 1/2 rsL2 * ||x||_2^2 + z * ||x||_1
%
%  By default, rsL2=0.
%  When rsL2 is nonzero, this corresponds to the well-know elastic net.
%
%% Input parameters:
%
%  A-         Matrix of size m x n
%                A can be a dense matrix
%                         a sparse matrix
%                         or a DCT matrix
%  y -        Response vector (of size mx1)
%  z -        L_1 norm regularization parameter (z >=0)
%  opts-      Optional inputs (default value: opts=[])
%
%% Output parameters:
%  x-         Solution
%  funVal-    Function value during iterations
%
%% Copyright (C) 2009-2010 Jun Liu, and Jieping Ye
%
% You are suggested to first read the Manual.
%
% For any problem, please contact with Jun Liu via j.liu@asu.edu
%
% Last modified on February 18, 2010.
%
%% Related papers
%
% [1]  Jun Liu and Jieping Ye, Efficient Euclidean Projections
%      in Linear Time, ICML 2009.
%
% [2]  Jun Liu and Jieping Ye, Sparse Learning with Efficient Euclidean
%      Projections onto the L1 Ball, Technical Report ASU, 2008.
%
%% Related functions
%
%  sll_opts, initFactor, pathSolutionLeast
%  LeastC, nnLeastR, nnLeastC,
%
%%

%% Verify and initialize the parameters
%%

% Verify the number of input parameters
if (nargin <3)
    error('\n Inputs: A, y and z should be specified!\n');
elseif (nargin==3)
    opts=[];
end

% Get the size of the matrix A
[m,n]=size(A);% m -- sample size,n -- feature dimension%******************************

% Verify the length of y
%if(size(y,1) ~=m)%********************************
if (length(y) ~=m)
    error('\n Check the length of y!\n');
end

% Verify the value of z
if (z<0)
    error('\n z should be nonnegative!\n');
end

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

%% Detailed initialization
%% Normalization

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

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

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

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

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

%% Starting point initialization

% compute AT y
if (opts.nFlag==0)
    ATy=A'*y; %A-->mxn;y-->mxk****************************
elseif (opts.nFlag==1)
    ATy=A'*y - sum(y) * mu';  ATy=ATy./nu;
else
    invNu=y./nu;              ATy=A'*invNu-sum(invNu)*mu';
end

% process the regularization parameter

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

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

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

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

% initialize a starting point
if opts.init==2
    %x=zeros(n,k);%*************************************
    x=zeros(n,1);
else
    if isfield(opts,'x0')
        x=opts.x0;
        if (length(x)~=n)
            error('\n Check the input .x0');
        end
    else
        x=ATy;  % if .x0 is not specified, we use ratio*ATy,
        % where ratio is a positive value
    end
end

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

if (opts.init==0) % If .init=0, we set x=ratio*x by "initFactor"
    % Please refer to the function initFactor for detail

    x_norm=sum(abs(x)); x_2norm=x'*x;%*****************************corresponding to L1norm,Fnorm
    if x_norm>=1e-6
        ratio=initFactor(x_norm, Ax, y, lambda,'LeastR', rsL2, x_2norm);
        x=ratio*x;    Ax=ratio*Ax;
    end
end

%% The main program

%% The Armijo Goldstein line search scheme + accelearted gradient descent
if (opts.mFlag==0 && opts.lFlag==0)

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

    L=1 + rsL2;
    % We assume that the maximum eigenvalue of A'A is over 1

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

    % alphap and alpha are used for computing the weight in forming search point
    alphap=0; alpha=1;

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

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

        % compute the gradient (g) at s
        As=Ax + beta* (Ax-Axp);

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

        % obtain the gradient g
        g=ATAs-ATy + rsL2 * s;%************************** +rsL2*s*inv(sigma)

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

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

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

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

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

            Av=Ax -As;
            r_sum=v'*v; l_sum=Av'*Av;

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

            % the condition is ||Av||_2^2 <= (L - rsL2) * ||v||_2^2
            if(l_sum <= r_sum * (L-rsL2))
                break;
            else
                L=max(2*L, l_sum/r_sum + rsL2);
                % fprintf('\n L=%5.6f',L);
            end
        end

        ValueL(iterStep)=L;

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

        xxp=x-xp;   Axy=Ax-y;
        funVal(iterStep)=Axy'* Axy/2 + rsL2/2 * x'*x + sum(abs(x)) * lambda;

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

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

%% Reformulated problem + Nemirovski's scheme

% .mFlag=1, and .lFlag=0
%  refomulate the problem as the constrained convex optimization
%  problem, and then apply Armijo Goldstein line search scheme

% Problem:
%    min  1/2 || A x - y||^2 + 1/2 rsL2 * ||x||_2^2 + z * t' * 1
%    s.t.   |x| <= t

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

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

    L=1 + rsL2;
    % We assume that the maximum eigenvalue of A'A is over 1

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

    % alphap and alpha are used for computing the weight in forming search point
    alphap=0; alpha=1;

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

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

        % compute the gradient (g) at s
        As=Ax + beta* (Ax-Axp);

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

        % obtain the gradient g
        g=ATAs-ATy + rsL2 * s;

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

        while (1)
            % let s walk in a step in the antigradient of s to get v
            % and then do the l1-norm regularized projection

            u=s-g/L;
            v= s_t - lambda / L;

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

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

            v_t=t-s_t;

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

            Av=Ax -As;
            r_sum=v'*v + v_t'*v_t; l_sum=Av'*Av + v'*v * rsL2;

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

            % the condition is ||Av||_2^2 + rsL2 * ||v||_2^2
            %                       <= L * (||v||_2^2 + ||v_t|| _2^2 )
            if(l_sum <= r_sum * L)
                break;
            else
                L=max(2*L, l_sum/r_sum);
                % fprintf('\n L=%5.6f',L);
            end
        end

        ValueL(iterStep)=L;

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

        xxp=x-xp;   Axy=Ax-y;
        funVal(iterStep)=Axy'* Axy/2 + rsL2/2 * x'*x + sum(t) * lambda;

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

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

end

%% adaptive line search

% .mFlag=1, and .lFlag=1
%  refomulate the problem as the constrained convex optimization
%  problem, and then apply adaptive line search scheme

% Problem:
%    min  1/2 || A x - y||^2 + 1/2 rsL2 * ||x||_2^2 + z * t' * 1
%    s.t.   |x| <= t

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

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

    L=1 + rsL2;
    % We assume that the maximum eigenvalue of A'A is over 1

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

    xp=x; Axp=Ax;
    % store x and Ax
    xxp=zeros(n,1);
    % the difference of x and xp
    t=abs(x); tp=t;
    % t is the upper bound of absolute value of x

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

    % We begin the adaptive line search in the following
    %
    % Note that, in the line search, L and beta are changing

    for iterStep=1:opts.maxIter

        ATAxp=ATAx;
        % store ATAx to ATAxp

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

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

                s=x + beta* xxp;   s_t= t + beta * (t -tp);
                As=Ax + beta* (Ax-Axp);
                ATAs=ATAx + beta * (ATAx- ATAxp);
                % compute the search point s, A * s, and A' * A * s
            else
                alpha= (-1+ sqrt(5)) / 2;
                beta=0; s=x; s_t=t; As=Ax; ATAs=ATAx;
            end

            g=ATAs-ATy + rsL2 * s;
            % compute the gradient g

            % let s walk in a step in the antigradient of s
            u=s-g/L;
            v= s_t - lambda / L;

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

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

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

            Av=Axnew -As;
            r_sum=v'*v + v_t'*v_t; l_sum=Av'*Av + v'*v * rsL2;

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

            % the condition is ||Av||_2^2 + rsL2 * ||v||_2^2
            %                       <= L * (||v||_2^2 + ||v_t|| _2^2 )
            if(l_sum <= r_sum * L)
                break;
            else
                L=max(2*L, l_sum/r_sum);
                % fprintf('\n L=%5.6f',L);
            end
        end
        %--------- Line Search for L ends

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

        ValueL(iterStep)=L;

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

        xp=x;    x=xnew; xxp=x-xp;
        Axp=Ax;  Ax=Axnew;
        % update x and Ax with xnew and Axnew
        tp=t; t=tnew;
        % update tp and t

        Axy=Ax-y;
        funVal(iterStep)=Axy' * Axy/2 + rsL2/2 * x'*x + lambda * sum(t);
        % compute function value

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

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

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

author: zzstefan

repo: BrainNetClass


text:
function [c,tc]=melcepst(s,fs,w,nc,p,n,inc,fl,fh)
%MELCEPST Calculate the mel cepstrum of a signal C=(S,FS,W,NC,P,N,INC,FL,FH)
%
%
% Simple use: (1) c=melcepst(s,fs)          % calculate mel cepstrum with 12 coefs, 256 sample frames
%			  (2) c=melcepst(s,fs,'E0dD')   % include log energy, 0th cepstral coef, delta and delta-delta coefs
%
% Inputs:
%     s	  speech signal
%     fs  sample rate in Hz (default 11025)
%     w   mode string (see below)
%     nc  number of cepstral coefficients excluding 0'th coefficient [default 12]
%     p   number of filters in filterbank [default: floor(3*log(fs)) =  approx 2.1 per ocatave]
%     n   length of frame in samples [default power of 2 < (0.03*fs)]
%     inc frame increment [default n/2]
%     fl  low end of the lowest filter as a fraction of fs [default = 0]
%     fh  high end of highest filter as a fraction of fs [default = 0.5]
%
%		w   any sensible combination of the following:
%
%               'R'  rectangular window in time domain
%				'N'	 Hanning window in time domain
%				'M'	 Hamming window in time domain (default)
%
%               't'  triangular shaped filters in mel domain (default)
%               'n'  hanning shaped filters in mel domain
%               'm'  hamming shaped filters in mel domain
%
%				'p'	 filters act in the power domain
%				'a'	 filters act in the absolute magnitude domain (default)
%
%               '0'  include 0'th order cepstral coefficient
%				'E'  include log energy
%				'd'	 include delta coefficients (dc/dt)
%				'D'	 include delta-delta coefficients (d^2c/dt^2)
%
%               'z'  highest and lowest filters taper down to zero (default)
%               'y'  lowest filter remains at 1 down to 0 frequency and
%			   	     highest filter remains at 1 up to nyquist freqency
%
%		       If 'ty' or 'ny' is specified, the total power in the fft is preserved.
%
% Outputs:	c     mel cepstrum output: one frame per row. Log energy, if requested, is the
%                 first element of each row followed by the delta and then the delta-delta
%                 coefficients.
%           tc    fractional time in samples at the centre of each frame
%                 with the first sample being 1.
%

% BUGS: (1) should have power limit as 1e-16 rather than 1e-6 (or possibly a better way of choosing this)
%           and put into VOICEBOX
%       (2) get rdct to change the data length (properly) instead of doing it explicitly (wrongly)

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: melcepst.m 4914 2014-07-24 08:44:26Z 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<2 fs=11025; end
if nargin<3 w='M'; end
if nargin<4 nc=12; end
if nargin<5 p=floor(3*log(fs)); end
if nargin<6 n=pow2(floor(log2(0.03*fs))); end
if nargin<9
   fh=0.5;
   if nargin<8
     fl=0;
     if nargin<7
        inc=floor(n/2);
     end
  end
end

if isempty(w)
   w='M';
end
if any(w=='R')
   [z,tc]=enframe(s,n,inc);
elseif any (w=='N')
   [z,tc]=enframe(s,hanning(n),inc);
else
   [z,tc]=enframe(s,hamming(n),inc);
end
f=rfft(z.');
[m,a,b]=melbankm(p,n,fs,fl,fh,w);
pw=f(a:b,:).*conj(f(a:b,:));
pth=max(pw(:))*1E-20;
if any(w=='p')
   y=log(max(m*pw,pth));
else
   ath=sqrt(pth);
   y=log(max(m*abs(f(a:b,:)),ath));
end
c=rdct(y).';
nf=size(c,1);
nc=nc+1;
if p>nc
   c(:,nc+1:end)=[];
elseif p<nc
   c=[c zeros(nf,nc-p)];
end
if ~any(w=='0')
   c(:,1)=[];
   nc=nc-1;
end
if any(w=='E')
   c=[log(max(sum(pw),pth)).' c];
   nc=nc+1;
end

% calculate derivative

if any(w=='D')
  vf=(4:-1:-4)/60;
  af=(1:-1:-1)/2;
  ww=ones(5,1);
  cx=[c(ww,:); c; c(nf*ww,:)];
  vx=reshape(filter(vf,1,cx(:)),nf+10,nc);
  vx(1:8,:)=[];
  ax=reshape(filter(af,1,vx(:)),nf+2,nc);
  ax(1:2,:)=[];
  vx([1 nf+2],:)=[];
  if any(w=='d')
     c=[c vx ax];
  else
     c=[c ax];
  end
elseif any(w=='d')
  vf=(4:-1:-4)/60;
  ww=ones(4,1);
  cx=[c(ww,:); c; c(nf*ww,:)];
  vx=reshape(filter(vf,1,cx(:)),nf+8,nc);
  vx(1:8,:)=[];
  c=[c vx];
end

if nargout<1
   [nf,nc]=size(c);
%    t=((0:nf-1)*inc+(n-1)/2)/fs;
   ci=(1:nc)-any(w=='0')-any(w=='E');
   imh = imagesc(tc/fs,ci,c.');
   axis('xy');
   xlabel('Time (s)');
   ylabel('Mel-cepstrum coefficient');
	map = (0:63)'/63;
	colormap([map map map]);
	colorbar;
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function simple_xy_rd ( )

%*****************************************************************************80
%
%% SIMPLE_XY_RD reads a 2 dimensional array from a NETCDF file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 December 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Russ Rew, Glenn Davis, Steve Emmerson, Harvey Davies, Ed Hartne,
%    The NetCDF User's Guide,
%    Unidata Program Center, March 2009.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SIMPLE_XY_RD:\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read a small 2D array from a NETCDF file.\n' );
%
%  Open the file.
%
  nc_filename = 'simple_xy.nc';
  ncid = netcdf.open ( nc_filename, 'NOWRITE' );
%
%  Get the varid of the data based on its name.
%
  varid = netcdf.inqVarID ( ncid, 'data_out' );
%
%  Read the data.
%
  data_out = netcdf.getVar ( ncid, varid );
%
%  Compare the data to the expected value.
%
  nj = 6;
  ni = 12;

  for j = 1 : nj
    for i = 1 : ni
      if ( data_out(i,j) ~= ( j - 1 ) * ni - ( i - 1 ) )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'SIMPLE_XY_RD - Fatal error!\n' );
        fprintf ( 1, '  Entry (%d,%d) of array seems to be wrong.\n', i, j );
      end
    end
  end
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The data is correct.\n' );
%
%  Close the file.
%
  netcdf.close ( ncid );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The data was read from the file "%s".\n', nc_filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SIMPLE_XY_RD:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end


text:
function r=rmse(data,estimate)
% Function to calculate root mean square error from a data vector or matrix
% and the corresponding estimates.
% Usage: r=rmse(data,estimate)
% Note: data and estimates have to be of same size
% Example: r=rmse(randn(100,100),randn(100,100));

% delete records with NaNs in both datasets first
I = ~isnan(data) & ~isnan(estimate);
data = data(I); estimate = estimate(I);

r=sqrt(sum((data(:)-estimate(:)).^2)/numel(data));


text:
%[2016]-"Jaya: A simple and new optimization algorithm for solving
%constrained and unconstrained optimization problems"

% (9/12/2020)

function JA = jJayaAlgorithm(feat,label,opts)
% Parameters
lb    = 0;
ub    = 1;
thres = 0.5;

if isfield(opts,'N'), N = opts.N; end
if isfield(opts,'T'), max_Iter = opts.T; end
if isfield(opts,'thres'), thres = opts.thres; end

% Objective function
fun = @jFitnessFunction;
% Number of dimensions
dim = size(feat,2);
% Initial (26)
X   = zeros(N,dim);
for i = 1:N
	for d = 1:dim
    X(i,d) = lb + (ub - lb) * rand();
  end
end
% Fitness
fit  = zeros(1,N);
fitG = inf;
for i = 1:N
  fit(i) = fun(feat,label,(X(i,:) > thres),opts);
  % Best
  if fit(i) < fitG
    fitG = fit(i);
    Xgb  = X(i,:);
  end
end
% Pre
Xnew = zeros(N,dim);

curve = zeros(1,max_Iter);
curve(1) = fitG;
t = 2;
% Iteration
while t <= max_Iter
  % Identify best & worst in population
  [~, idxB] = min(fit);
  Xbest     = X(idxB,:);
  [~, idxW] = max(fit);
  Xworst    = X(idxW,:);
  % Start
  for i = 1:N
    for d = 1:dim
      % Random numbers
      r1 = rand();
      r2 = rand();
      % Update (1)
      Xnew(i,d) = X(i,d) + r1 * (Xbest(d) - abs(X(i,d))) - ...
        r2 * (Xworst(d) - abs(X(i,d)));
    end
    % Boundary
    XB = Xnew(i,:); XB(XB > ub) = ub; XB(XB < lb) = lb;
    Xnew(i,:) = XB;
  end
  % Fitness
  for i = 1:N
    Fnew = fun(feat,label,(Xnew(i,:) > thres),opts);
    % Greedy selection
    if Fnew < fit(i)
      fit(i) = Fnew;
      X(i,:) = Xnew(i,:);
    end
    % Best
    if fit(i) < fitG
      fitG = fit(i);
      Xgb  = X(i,:);
    end
  end
  % Save
  curve(t) = fitG;
  fprintf('\nIteration %d Best (JA)= %f',t,curve(t))
  t = t + 1;
end
% Select features based on selected index
Pos   = 1:dim;
Sf    = Pos((Xgb > thres) == 1);
sFeat = feat(:,Sf);
% Store results
JA.sf = Sf;
JA.ff = sFeat;
JA.nf = length(Sf);
JA.c  = curve;
JA.f  = feat;
JA.l  = label;
end

author: JingweiToo

repo: Wrapper-Feature-Selection-Toolbox


text:
 function [mat, nexp] = filtmat(arg0, arg1, arg2, arg3, arg4)
%function mat = filtmat('1d', kernx, nx)
%			1d nonperiodic
%function mat = filtmat({'1d' 'per'}, kernx, nx)
%			1d periodic
%function mat = filtmat('2d,mul,sep', kernx, kerny, [nx ny])
%			2d multiplicatively separable nonperiodic
%function mat = filtmat('3d,mul,sep', kernx, kerny, kernz, [nx ny nz])
%			3d multiplicatively separable nonperiodic
%			kernel(x,y,z) = kernx(x) * kerny(y) * kernz(z)
%function [mat, nexp] = filtmat({'3d,mul,sep' 'expand'}, kernx, kerny, kernz, [nx ny nz])
%			3d multiplicatively separable nonperiodic, expanded
%function mat = filtmat('1d,causal', kernx, nx)
%			1d causal filter, ala convmtx
%
%	Form matrix to do 3d linear filtering - odd kernel length.
%	expanding is like what happens in 'conv' - output is longer

if ~nargin, help filtmat, error arg, end
flag_expand = any(strcmp(arg0, 'expand'));

if nargin == 3 && strcmp(arg0, '1d,causal')
	mat = filtmat1causal(arg1, arg2);
	return

elseif nargin == 3 && any(strcmp(arg0, '1d'))
	if any(strcmp(arg0, 'per'))
		mat = filtmat1per(arg1, arg2);
	else
		mat = filtmat1(arg1, arg2, flag_expand);
	end
	return

elseif nargin == 4 && any(strcmp(arg0, '2d,mul,sep'))
	mat = filtmat('3d,mul,sep', arg1, arg2, 1, [arg3 1]);
	return

elseif nargin == 5 && any(strcmp(arg0, '3d,mul,sep'))
	k.x = arg1;
	k.y = arg2;
	k.z = arg3;
	n.x = arg4(1);
	n.y = arg4(2);
	n.z = arg4(3);
else
	nargin, help filtmat, error arg
end

	if flag_expand

		[p.x, o.x] = filtmat1(k.x, n.x, flag_expand);
		[p.y, o.y] = filtmat1(k.y, n.y, flag_expand);
		[p.z, o.z] = filtmat1(k.z, n.z, 0);	% no z expansion

		p.x = kron(speye(n.y*n.z), p.x);
		p.y = kron(speye(n.z), kron(p.y, speye(o.x)));
		p.z = kron(p.z, speye(o.x*o.y));

		if 0
			n.xyz = n.x * n.y * n.z;
			t = (o.y-n.y)/2 * o.x;
			p.x = [sparse(t,n.xyz); p.x; sparse(t,n.xyz)];
			t = zeros(o.x,o.y);
			t([1:n.x]+(o.x-n.y)/2-1,:) = 1;
			tt = sparse(o.x*o.y, n.x*n.y);
			tt(find(t(:)),:) = p.y;
			p.y = tt;
			p.z = sparse(1,1);
			if n.z ~= 1, error notdone, end
		end
	else
		p.x = filtmat1(k.x, n.x, flag_expand);
		p.y = filtmat1(k.y, n.y, flag_expand);
		p.z = filtmat1(k.z, n.z, flag_expand);

		p.x = kron(speye(n.y*n.z), p.x);
		p.y = kron(speye(n.z), kron(p.y, speye(n.x)));
		p.z = kron(p.z, speye(n.x*n.y));
		o = n;
	end

%	mat = p.x + p.y + p.z;
	if flag_expand, warning('fix: expand may not work for multiplicative?'), end
	mat = p.z * p.y * p.x;
	nexp = [o.x o.y o.z];

function mat = filtmat1causal(kern, n)
	nk = length(kern);
	if nk < n, error n, end
	if n ~= nk, error 'not done', end
	t = kern(:,ones(1,n));
	t = [t; zeros(nk,nk)];
	t = t(:);
	t((end-nk+1):end) = [];
	mat = reshape(t, 2*nk-1,nk);
	mat = mat(1:nk,:);

function [mat, nexp] = filtmat1(kern, n, flag_expand)
%	1d matrix that does linear filtering - odd kernel length.
%	flag_expand means output has more entries than input (n+nk-1)

	nk = length(kern);
	ic = (nk+1)/2;		% center index
	if round(ic) ~= ic,	error('odd kernel only'), end

	if flag_expand
		n = n + nk-1;
	end
	nexp = n;

	mat = diag(ones(n,1)*kern(ic));
	for iv=1:(ic-1)
		mat = mat + diag(ones(n-iv,1) * kern(ic-iv), -iv) ...
			  + diag(ones(n-iv,1) * kern(ic+iv), iv);
	end

	if flag_expand
		mat = mat(:,ic:(end-ic+1));
	end
	mat = sparse(mat);

function mat = filtmat1per(kern, nn)
%	1d matrix that does periodic linear filtering.  odd or even kernel ok.

	nk = length(kern);
	ic = (nk+1)/2;

	mat = spdiag(ones(nn,1)*kern(ic));
	for kk=1:(ic-1)
		mat = mat ...
			+ diag(ones(nn-kk,1) * kern(ic-kk),	-kk) ...
			+ diag(ones(kk   ,1) * kern(ic-kk),	nn-kk);
	end

	for kk=1:floor(nk/2)
		mat = mat ...
			+ diag(ones(nn-kk,1) * kern(ic+kk),	kk) ...
			+ diag(ones(kk   ,1) * kern(ic+kk),	kk-nn);
	end
	mat = sparse(mat);


text:
function [ vid_locs, bboxes, gts_all, invalid_frames ] = CollectTestData( db_root, bb_root, extra_dir )
%COLLECTTESTDATA Summary of this function goes here
%   Detailed explanation goes here

    cat_1 = [ 114, 124, 125, 126, 150, 158, 401, 402, 505, 506, 507, 508, 509, 510, 511, 514, 515, 518, 519, 520, 521, 522, 524, 525, 537, 538, 540, 541, 546, 547, 548];
    cat_2 = [203, 208, 211, 212, 213, 214, 218, 224, 403, 404, 405, 406, 407, 408, 409, 412, 550, 551, 553];
    cat_3 = [410, 411, 516, 517, 526, 528, 529, 530, 531, 533, 557, 558, 559, 562];

    all_test = cat(2, cat_1, cat_2, cat_3);

    vid_locs = cell(numel(all_test),1);
    bboxes = cell(numel(all_test),1);
    gts_all = cell(numel(all_test),1);
    invalid_frames = cell(numel(all_test),1);
    for i=1:numel(all_test)
        vid_locs{i} = [db_root '/', num2str(all_test(i)), '/vid.avi'];
        bboxes{i} = dlmread([bb_root,  num2str(all_test(i)), '_dets.txt'], ',');

        %% Grab the ground truth
        fps_all = dir([db_root, '/', num2str(all_test(i)), '/annot/*.pts']);
        gt_landmarks = zeros([68, 2, size(fps_all)]);
        for k = 1:size(fps_all)
            gt_landmarks_frame = dlmread([db_root, '/', num2str(all_test(i)), '/annot/', fps_all(k).name], ' ', 'A4..B71');
            gt_landmarks(:,:,k) = gt_landmarks_frame;
        end
        % Remove unreliable frames
        if(exist([extra_dir, '/', num2str(all_test(i)), '.mat'], 'file'))
            inv_frames = load([extra_dir, '/', num2str(all_test(i)), '.mat']);

            gt_landmarks(:,:,int32(inv_frames.error)) = [];
            invalid_frames{i} = inv_frames.error;
        end
        gts_all{i} = gt_landmarks;
    end

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function copyBinFile(ops1)

fid    = cell(size(ops1));
for i = 1:size(ops1,1)
    for j = 1:size(ops1,2)
        fid{i,j}           = fopen(ops1{i,j}.RegFile, 'r');
        folder = fullfile(ops1{i,j}.RegFileBinLocation, ops1{i,j}.mouse_name, ...
            ops1{i,j}.date, ops1{i,j}.CharSubDirs);
        if ~exist(folder, 'dir')
            mkdir(folder)
        end
        fidCopy = fopen(fullfile(folder, ...
            sprintf('%s_%s_%s_plane%d.bin', ops1{i,j}.mouse_name, ops1{i,j}.date, ...
            ops1{i,j}.CharSubDirs, i)), 'w');
        sz = ops1{i,j}.Lx * ops1{i,j}.Ly;
        parts = ceil(sum(ops1{i,j}.Nframes) / 2000);
        for p = 1:parts
            toRead = 2000;
            if p == parts
                toRead = sum(ops1{i,j}.Nframes) - 2000 * (parts-1);
            end
            data = fread(fid{i,j},  sz*toRead, '*int16');
            fwrite(fidCopy, data, class(data));
        end
        fclose(fidCopy);
        fclose(fid{i,j});
    end
end

author: cortex-lab


text:
function exact = p05_exact ( )

%*****************************************************************************80
%
%% P05_EXACT returns the exact integral for problem 5.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real EXACT, the value of the integral.
%
  exact = 4.0 / 9.0;

  return
end


text:
function [acc pred softpred]  = testBNB(model, allSNumBi, labels, params)

nbsenttst = allSNumBi;
nblbltst = labels;

pred = zeros(size(labels));

for i = 1:length(nbsenttst)
    %words(nbsenttst{i})
    nbsentraw = nbsenttst{i};
    nbsentraw=nbsentraw(nbsentraw>0);
    ind = zeros(params.dictsize, 1);
%     updates = unique(nbsentraw);
%     counts = power(histc(nbsentraw, updates)', params.testp);

    ind(nbsentraw) = 1;
    ind = logical(ind);
    neglog = sum(model.lpwneg(ind)) + sum(model.mlpwneg(~ind));
    poslog = sum(model.lpwpos(ind)) + sum(model.mlpwpos(~ind));
    if neglog + model.pneg  > poslog + model.ppos
        pred(i) = 0;
    else
        pred(i) = 1;
    end
    softpred(i) = poslog - neglog;
end

acc = sum(pred == labels) / length(labels);
acc = acc*100

author: sidaw

repo: nbsvm

save_path: github-repos/MATLAB/sidaw-nbsvm


text:
function [f, e] = fv(B, nk, thmax, noplot)
%FV     Field of values (or numerical range).
%       FV(A, NK, THMAX) evaluates and plots the field of values of the
%       NK largest leading principal submatrices of A, using THMAX
%       equally spaced angles in the complex plane.
%       The defaults are NK = 1 and THMAX = 16.
%       (For a `publication quality' picture, set THMAX higher, say 32.)
%       The eigenvalues of A are displayed as `x'.
%       Alternative usage: [F, E] = FV(A, NK, THMAX, 1) suppresses the
%       plot and returns the field of values plot data in F, with A's
%       eigenvalues in E.   Note that NORM(F,INF) approximates the
%       numerical radius,
%                 max {abs(z): z is in the field of values of A}.

%       Theory:
%       Field of values FV(A) = set of all Rayleigh quotients. FV(A) is a
%       convex set containing the eigenvalues of A.  When A is normal FV(A) is
%       the convex hull of the eigenvalues of A (but not vice versa).
%               z = x'Ax/(x'x),  z' = x'A'x/(x'x)
%               => REAL(z) = x'Hx/(x'x),   H = (A+A')/2
%       so      MIN(EIG(H)) <= REAL(z) <= MAX(EIG(H)),
%       with equality for x = corresponding eigenvectors of H.  For these x,
%       RQ(A,x) is on the boundary of FV(A).
%
%       Based on an original routine by A. Ruhe.
%
%       References:
%       R. A. Horn and C. R. Johnson, Topics in Matrix Analysis, Cambridge
%            University Press, 1991; sec. 1.5.
%       A. S. Householder, The Theory of Matrices in Numerical Analysis,
%            Blaisdell, New York, 1964; sec. 3.3.
%       C. R. Johnson, Numerical determination of the field of values of a
%            general complex matrix, SIAM J. Numer. Anal., 15 (1978),
%            pp. 595-602.

if nargin < 2 | isempty(nk), nk = 1; end
if nargin < 3 | isempty(thmax), thmax = 16; end
thmax = thmax - 1;  % Because code below uses thmax + 1 angles.

iu = sqrt(-1);
[n, p] = size(B);
if n ~= p, error('Matrix must be square.'), end
f = [];
z = zeros(2*thmax+1,1);
e = eig(B);

% Filter out cases where B is Hermitian or skew-Hermitian, for efficiency.
if isequal(B,B')

   f = [min(e) max(e)];

elseif isequal(B,-B')

   e = imag(e);
   f = [min(e) max(e)];
   e = iu*e; f = iu*f;

else

for m = 1:nk

   ns = n+1-m;
   A = B(1:ns, 1:ns);

   for i = 0:thmax
      th = i/thmax*pi;
      Ath = exp(iu*th)*A;               % Rotate A through angle th.
      H = 0.5*(Ath + Ath');             % Hermitian part of rotated A.
      [X, D] = eig(H);
      [lmbh, k] = sort(real(diag(D)));
      z(1+i) = rq(A,X(:,k(1)));         % RQ's of A corr. to eigenvalues of H
      z(1+i+thmax) = rq(A,X(:,k(ns)));  % with smallest/largest real part.
   end

   f = [f; z];

end
% Next line ensures boundary is `joined up' (needed for orthogonal matrices).
f = [f; f(1,:)];

end
if thmax == 0; f = e; end

if nargin < 4

   ax = cpltaxes(f);
   plot(real(f), imag(f))      % Plot the field of values
   axis(ax);
   axis('square');

   hold on
   plot(real(e), imag(e), 'x')    % Plot the eigenvalues too.
   hold off

end

function z = rq(A,x)
%RQ      Rayleigh quotient.
%        RQ(A,x) is the Rayleigh quotient of A and x, x'*A*x/(x'*x).

z = x'*A*x/(x'*x);

author: opencobra

repo: cobratoolbox


text:
%FIX_LINES  Improves the line style of eps files generated by print
%
% Examples:
%   fix_lines fname
%   fix_lines fname fname2
%   fstrm_out = fixlines(fstrm_in)
%
% This function improves the style of lines in eps files generated by
% MATLAB's print function, making them more similar to those seen on
% screen. Grid lines are also changed from a dashed style to a dotted
% style, for greater differentiation from dashed lines.
%
% The function also places embedded fonts after the postscript header, in
% versions of MATLAB which place the fonts first (R2006b and earlier), in
% order to allow programs such as Ghostscript to find the bounding box
% information.
%
%IN:
%   fname - Name or path of source eps file.
%   fname2 - Name or path of destination eps file. Default: same as fname.
%   fstrm_in - File contents of a MATLAB-generated eps file.
%
%OUT:
%   fstrm_out - Contents of the eps file with line styles fixed.

% Copyright: (C) Oliver Woodford, 2008-2014

% The idea of editing the EPS file to change line styles comes from Jiro
% Doke's FIXPSLINESTYLE (fex id: 17928)
% The idea of changing dash length with line width came from comments on
% fex id: 5743, but the implementation is mine :)

% Thank you to Sylvain Favrot for bringing the embedded font/bounding box
% interaction in older versions of MATLAB to my attention.
% Thank you to D Ko for bringing an error with eps files with tiff previews
% to my attention.
% Thank you to Laurence K for suggesting the check to see if the file was
% opened.

% 01/03/15: Issue #20: warn users if using this function in HG2 (R2014b+)
% 27/03/15: Fixed out of memory issue with enormous EPS files (generated by print() with OpenGL renderer), related to issue #39

function fstrm = fix_lines(fstrm, fname2)

% Issue #20: warn users if using this function in HG2 (R2014b+)
if using_hg2
    warning('export_fig:hg2','The fix_lines function should not be used in this Matlab version.');
end

if nargout == 0 || nargin > 1
    if nargin < 2
        % Overwrite the input file
        fname2 = fstrm;
    end
    % Read in the file
    fstrm = read_write_entire_textfile(fstrm);
end

% Move any embedded fonts after the postscript header
if strcmp(fstrm(1:15), '%!PS-AdobeFont-')
    % Find the start and end of the header
    ind = regexp(fstrm, '[\n\r]%!PS-Adobe-');
    [ind2, ind2] = regexp(fstrm, '[\n\r]%%EndComments[\n\r]+');
    % Put the header first
    if ~isempty(ind) && ~isempty(ind2) && ind(1) < ind2(1)
        fstrm = fstrm([ind(1)+1:ind2(1) 1:ind(1) ind2(1)+1:end]);
    end
end

% Make sure all line width commands come before the line style definitions,
% so that dash lengths can be based on the correct widths
% Find all line style sections
ind = [regexp(fstrm, '[\n\r]SO[\n\r]'),... % This needs to be here even though it doesn't have dots/dashes!
       regexp(fstrm, '[\n\r]DO[\n\r]'),...
       regexp(fstrm, '[\n\r]DA[\n\r]'),...
       regexp(fstrm, '[\n\r]DD[\n\r]')];
ind = sort(ind);
% Find line width commands
[ind2, ind3] = regexp(fstrm, '[\n\r]\d* w[\n\r]');
% Go through each line style section and swap with any line width commands
% near by
b = 1;
m = numel(ind);
n = numel(ind2);
for a = 1:m
    % Go forwards width commands until we pass the current line style
    while b <= n && ind2(b) < ind(a)
        b = b + 1;
    end
    if b > n
        % No more width commands
        break;
    end
    % Check we haven't gone past another line style (including SO!)
    if a < m && ind2(b) > ind(a+1)
        continue;
    end
    % Are the commands close enough to be confident we can swap them?
    if (ind2(b) - ind(a)) > 8
        continue;
    end
    % Move the line style command below the line width command
    fstrm(ind(a)+1:ind3(b)) = [fstrm(ind(a)+4:ind3(b)) fstrm(ind(a)+1:ind(a)+3)];
    b = b + 1;
end

% Find any grid line definitions and change to GR format
% Find the DO sections again as they may have moved
ind = int32(regexp(fstrm, '[\n\r]DO[\n\r]'));
if ~isempty(ind)
    % Find all occurrences of what are believed to be axes and grid lines
    ind2 = int32(regexp(fstrm, '[\n\r] *\d* *\d* *mt *\d* *\d* *L[\n\r]'));
    if ~isempty(ind2)
        % Now see which DO sections come just before axes and grid lines
        ind2 = repmat(ind2', [1 numel(ind)]) - repmat(ind, [numel(ind2) 1]);
        ind2 = any(ind2 > 0 & ind2 < 12); % 12 chars seems about right
        ind = ind(ind2);
        % Change any regions we believe to be grid lines to GR
        fstrm(ind+1) = 'G';
        fstrm(ind+2) = 'R';
    end
end

% Define the new styles, including the new GR format
% Dot and dash lengths have two parts: a constant amount plus a line width
% variable amount. The constant amount comes after dpi2point, and the
% variable amount comes after currentlinewidth. If you want to change
% dot/dash lengths for a one particular line style only, edit the numbers
% in the /DO (dotted lines), /DA (dashed lines), /DD (dot dash lines) and
% /GR (grid lines) lines for the style you want to change.
new_style = {'/dom { dpi2point 1 currentlinewidth 0.08 mul add mul mul } bdef',... % Dot length macro based on line width
             '/dam { dpi2point 2 currentlinewidth 0.04 mul add mul mul } bdef',... % Dash length macro based on line width
             '/SO { [] 0 setdash 0 setlinecap } bdef',... % Solid lines
             '/DO { [1 dom 1.2 dom] 0 setdash 0 setlinecap } bdef',... % Dotted lines
             '/DA { [4 dam 1.5 dam] 0 setdash 0 setlinecap } bdef',... % Dashed lines
             '/DD { [1 dom 1.2 dom 4 dam 1.2 dom] 0 setdash 0 setlinecap } bdef',... % Dot dash lines
             '/GR { [0 dpi2point mul 4 dpi2point mul] 0 setdash 1 setlinecap } bdef'}; % Grid lines - dot spacing remains constant

% Construct the output
% This is the original (memory-intensive) code:
%first_sec = strfind(fstrm, '% line types:'); % Isolate line style definition section
%[second_sec, remaining] = strtok(fstrm(first_sec+1:end), '/');
%[remaining, remaining] = strtok(remaining, '%');
%fstrm = [fstrm(1:first_sec) second_sec sprintf('%s\r', new_style{:}) remaining];
fstrm = regexprep(fstrm,'(% line types:.+?)/.+?%',['$1',sprintf('%s\r',new_style{:}),'%']);

% Write the output file
if nargout == 0 || nargin > 1
    read_write_entire_textfile(fname2, fstrm);
end
end

author: altmany


text:
function verbose = prefsFormatCheck
% checks 'VISTA' preference groupfor a file format setting
%
%   verbose = prefsFormatCheck;
%
% This variable determines how many data files (including t1-weighted
% anatomies, segmentation files, and parameter maps) are saved in mrVista.
%
% There are two possible values allowed:
%	'default': This indicates that the pre-existing, specialized formats
%	adopted by mrLoadRet/mrVista will be used:
%		volume anatomies will use the mrGray *.dat format;
%		inplane anatomies will be saved as Inplane/anat.mat in each
%		session;
%		segmentation files will be mrGray *.class / *.gray files;
%		parameter maps will be saved as .mat files.
%
%	'nifti': these formats will all employ NIFTI files. The segmentation
%	files (normally kept in .class files for each hemisphere) should be
%	saved as a labeled NIFTI file from ITKGray; anatomies and parameter
%	maps will be loaded/saved as compressed NIFTI. ***as of 6/2008, this
%	functionality is not fully implemented.***
%		* The 'grow_gray' MEX file needs to be tested and stable on all
%		machines: whereas before mrGray was a Windows-only bottleneck, once
%		the gray matter graph was saved, it could be read from any system.
%		The grow_gray would need to be employed everytime a segmentation is
%		installed, and is currently not stable enough to rely on it;
%		* The inplane anatomy/paramter map code needs to be updated to
%		allow for NIFTI files.
%
% ras, 09/2006.

if ~ispref('VISTA', 'fileFormat') | ~ismember( getpref('VISTA', 'fileFormat'), {'default' 'nifti'} )
    setpref('VISTA', 'fileFormat', 'default');
    fprintf('[%s]: ', mfilename);
    fprintf('Initializing VISTA preference ''fileFormat'' to ''default''. \n');
    fprintf('This will cause files to be saved and loaded in the \n');
	fprintf('legacy, specialized mrVista format: \n');
	fprintf('\tAnatomy files are in mrGray *.dat format;\n');
	fprintf('\tSegmentation files are in mrGray *.class and *.gray formats;\n');
	fprintf('\tInplane anatomies are saved as Inplane/anat.mat;\n');
	fprintf('\tParameter maps are saved as .mat files.\n');
	fprintf('To save these files in NIFTI format, run: \n');
	fprintf('setpref(''VISTA'', ''fileFormat'', ''nifti''). \n');
end

verbose = getpref('VISTA', 'fileFormat');

return


text:
function [restr,tpl]=create_restrictions_and_markov_chains1(tpl)
% create_restrictions_and_markov_chains1 -- creates linear restrictions and
% markov chains for the SVAR model in which coefficients are switching
% across all equations (synchronized case)
%
% ::
%
%
%   [restr,tpl]=create_restrictions_and_markov_chains1(tpl)
%
% Args:
%
%    - **tpl** [struct]: template created for SVAR objects
%
% Returns:
%    :
%
%    - **restr** [cell]: two column-cell (see below). The first column
%    contains COEF objects or linear combinations of COEF objects, which are
%    themselves COEF objects.
%
%    - **tpl** [struct]: modified template
%
% Note:
%
%    - The syntax to construct an advanced COEF object is
%    a=coef(eqtn,vname,lag,chain_name,state)
%      - **eqtn** [integer|char]: integer or variable name
%      - **vname** [integer|char]: integer or variable name
%      - **lag** [integer]: integer or variable name
%      - **chain_name** [char]: name of the markov chain
%      - **state** [integer]: state number
%
%    - RISE sorts the endogenous variables alphabetically and use this order
%    to tag each equation in SVAR and RFVAR models.
%
%    - The lag coefficients are labelled a0, a1, a2,...,ak, for a model with k
%    lags. Obviously, a0 denotes the contemporaneous coefficients.
%
%    - The constant terms labelled c_1_1, c_2_2,...,c_n_n, for a model with n
%    endogenous variables.
%
%    - The standard deviations labelled sig_1_1, sig_2_2,...,sig_n_n, for a
%    model with n endogenous variables.
%
% Example:
%    coef('pi','ygap',0,'policy',1)
%
%    coef(2,3,0,'policy',1)
%
%    coef(2,'ygap',0,'policy',1)
%
%    coef('pi',3,0,'policy',1)
%
%    See also:

% We add a Markov chain to the template
%----------------------------------------
% N.B: The chain controls the coefficients of all equations but not the
% variance, which remains constant.
last=numel(tpl.markov_chains);
tpl.markov_chains(last+1)=struct('name','syncoef',...
    'states_expected_duration',[3+1i,3+1i],...
    'controlled_parameters',{{'c','a0','a1','a2'}});

% syntax is coef(eqtn,vname,lag,chain_name,state)
%------------------------------------------------
restr=cell(0,2);

numberOfStates=numel(tpl.markov_chains(end).states_expected_duration);

for istate=1:numberOfStates
    restr=[restr
        {
        % first equation or "FFR" equation:
        %----------------------------------
        coef(1,'pi',1,'syncoef',istate),0
        coef(1,'pi',2,'syncoef',istate),0
        coef(1,'ygap',1,'syncoef',istate),0
        coef(1,'ygap',2,'syncoef',istate),0
        coef(1,'FFR',2,'syncoef',istate),0
        % second equation or "pi" equation
        %----------------------------------
        coef(2,'FFR',0,'syncoef',istate),0
        coef(2,'FFR',1,'syncoef',istate),0
        coef(2,'FFR',2,'syncoef',istate),0
        coef(2,'ygap',1,'syncoef',istate),0
        coef(2,'ygap',2,'syncoef',istate),0
        % third equation or "ygap" equation
        %-----------------------------------
        coef(3,'FFR',1,'syncoef',istate),0
        coef(3,'FFR',2,'syncoef',istate),0
        coef(3,'pi',1,'syncoef',istate),0
        coef(3,'pi',2,'syncoef',istate),0
        coef(3,'pi',0,'syncoef',istate)+coef(3,'FFR',0,'syncoef',istate),0
        }
        ]; %#ok<AGROW>
end

end

author: jmaih


text:
function [w]=obtencion_w(R)
    j=1;
    x=1;
    i=1;

    %Desmontamos la matriz
    while j<100
        r11(j)=R{j}(1,1);
        r12(j)=R{j}(1,2);
        r13(j)=R{j}(1,3);

        r21(j)=R{j}(2,1);
        r22(j)=R{j}(2,2);
        r23(j)=R{j}(2,3);

        r31(j)=R{j}(3,1);
        r32(j)=R{j}(3,2);
        r33(j)=R{j}(3,3);

        j=j+1;
    end

    %Derivamos los vectores
    rd11=diff(r11);
    rd12=diff(r12);
    rd13=diff(r13);

    rd21=diff(r21);
    rd22=diff(r22);
    rd23=diff(r23);

    rd31=diff(r31);
    rd32=diff(r32);
    rd33=diff(r33);

    %Motamos la matriz derivada
    while x<99
        dR{x}(1,1)=rd11(x);
        dR{x}(1,2)=rd12(x);
        dR{x}(1,3)=rd13(x);

        dR{x}(2,1)=rd21(x);
        dR{x}(2,2)=rd22(x);
        dR{x}(2,3)=rd23(x);

        dR{x}(3,1)=rd31(x);
        dR{x}(3,2)=rd32(x);
        dR{x}(3,3)=rd33(x);

        x=x+1;
    end

    %Multiplicamos dR por R(transpuesta) para sacar w
    while i<99
        w{i}=dR{i}*(R{i}.');

        i=i+1;
    end
end

repo: ARTE


text:
function forest = SCiForest(X, t, psi, num_attributes, num_hyperplanes)
%iForest Construct a "forest" (set) of isolation trees
%   X: the data, each row is an instance
%   t: the number of trees to construct
%   psi: the sub-sampling size

height_limit = ceil(log2(psi));
[numdata ~] = size(X);
if nargin < 5
    num_hyperplanes = 10;
end
if nargin < 4
    num_attributes = 2;
end

for i=1:t
    sampleind = randsample(1:numdata, psi);
    X_sample = X(sampleind, :);
    forest(i) = SCiTree(X_sample, num_attributes, num_hyperplanes, 0, height_limit);
end

end

author: dsmi-lab-ntust

repo: AnomalyDetectionToolbox


text:
% Focussed Detector in 3D Example
%
% This example shows how k-Wave can be used to model the output of a
% focussed bowl detector in 3D where the directionality arises from
% spatially averaging across the detector surface. It builds on the
% Focussed Detector in 2D example.
%
% author: Ben Cox
% date: 29th October 2010
% last update: 18th September 2012
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
%
% k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

clear all;

% =========================================================================
% SIMULATION
% =========================================================================

% create the computational grid
Nx = 64;            % number of grid points in the x direction
Ny = 64;            % number of grid points in the y direction
Nz = 64;            % number of grid points in the z direction
dx = 100e-3/Nx;     % grid point spacing in the x direction [m]
dy = dx;            % grid point spacing in the y direction [m]
dz = dx;            % grid point spacing in the z direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);

% define the properties of the propagation medium
medium.sound_speed = 1500;	% [m/s]

% create the time array
[kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed);
Nt = length(kgrid.t_array);

% create a concave sensor
radius = Nx/4-1;
height = 10;
ss = makeSphericalSection(radius, height);

% add it to a mask of the correct size
sensor.mask = zeros(Nx, Ny, Nz);
[a,b,c] = size(ss);
sphere_offset = 10;
sensor.mask(sphere_offset+(1:a),(Ny-b+1)/2+(1:b),(Nz-c+1)/2+(1:c)) = ss;

% define a time varying sinusoidal source
source_freq = 0.25e6;
source_mag = 1;
source.p = source_mag*sin(2*pi*source_freq*kgrid.t_array);
source.p = filterTimeSeries(kgrid, medium, source.p);

% place the first point source near the focus of the detector
source1 = zeros(Nx, Ny, Nz);
source1(sphere_offset+radius, Ny/2+1, Nz/2+1) = 1;

% place the second point source off axis
source2 = zeros(Nx, Ny, Nz);
source2(sphere_offset+radius, Ny/2+6, Nz/2+6) = 1;

% run the first simulation
source.p_mask = source1;
input_args = {'PMLSize', 10, 'DataCast', 'single', 'PlotSim', false};
sensor_data1 = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});

% average the data recorded at each grid point to simulate the measured
% signal from a single element focussed detector
sensor_data1 = sum(sensor_data1, 1);

% run the second simulation
source.p_mask = source2;
sensor_data2 = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});

% average the data recorded at each grid point to simulate the measured
% signal from a single element focussed detector
sensor_data2 = sum(sensor_data2, 1);

% =========================================================================
% VISUALISATION
% =========================================================================

% plot the detector and on-axis and off-axis point sources
voxelPlot(sensor.mask + source1 + source2);
view([99, 18]);

% plot the time series corresponding to the on-axis and off-axis sources
figure
[t_sc, t_scale, t_prefix] = scaleSI(kgrid.t_array(end));
plot(kgrid.t_array.*t_scale, sensor_data1, '-');
hold on
plot(kgrid.t_array.*t_scale, sensor_data2, 'r-');
xlabel(['Time [' t_prefix 's]']);
ylabel('Average Pressure Measured By Focussed Detector [Pa]');
legend('Source on axis', 'Source off axis');

repo: Aero-matlab


text:
function [R T] = fX4_to_RT(Tm)

R = Tm(1:3,1:3);
T = Tm(1:3,4)';

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
% @Author: amishkin
% @Date:   2018-07-10T14:56:18-07:00
% @Email:  amishkin@cs.ubc.ca
% @Last modified by:   aaronmishkin
% @Last modified time: 2018-07-25T13:29:13-07:00

% ###############################################
% ############# Implements mf-exact #############
% ###############################################

function [nlz, log_loss, Sigma, mu] = mf_exact(method_name, y, X, gamma, y_te, X_te, options, mu_start, sigma_start)

    fprintf('%s\n',method_name);
    [N,D] = size(X);

    % set default options
    [max_iters, lowerBoundTol, display, num_samples, beta_start, alpha_start, decay_rate, mini_batch_size] = myProcessOptions(options, ...
    'max_iters', 2000, 'lowerBoundTol',1e-4, 'display', 1, 'num_samples', 1, 'beta', 0.1, 'alpha', 0.8, 'decay_rate', 0.55, 'mini_batch_size', N);

    % minfunc options
    optMinFunc = struct('display', display,...
        'Method', 'lbfgs',...
        'DerivativeCheck', 'off',...
        'LS', 2,...
        'recordPath', 1, ...
        'recordPathIters', 1, ...
        'MaxIter', max_iters+1,...
        'MaxFunEvals', max_iters+1,...
        'TolFun', lowerBoundTol,......
        'TolX', lowerBoundTol);

    V = sigma_start;
    m = mu_start;

    v0 = [m; sqrt(diag(V))];
    funObj = @funObj_mfvi_exact;

    % compute loss at iter =0
    post_dist.mean = m;
    post_dist.covMat =V;
    iter = 0;
    [pred, log_lik]=get_loss(iter, post_dist, X, y, gamma, X_te, y_te);
    nlz(iter+1)=-log_lik;%nlz0
    log_loss(iter+1)=pred.log_loss;%log_loss0

    % optimize using minfunc
    [v, f, exitflag, inform] = minFunc(funObj, v0, optMinFunc, y, X, gamma);
    v_all = inform.trace.x;

    % compute loss for iter>0
    for ii=1:size(v_all,2)
        vi = v_all(:,ii);
        post_dist.mean = vi(1:D);
        U = diag(vi(D+1:end));
        post_dist.covMat = U'*U;
        if sum(eig(post_dist.covMat) <= 1e-8) > 0
            ALERT = 'NOT PD'
            nlz(ii) = NaN;
            log_loss(ii) = NaN;
            continue
        end
    end
    [pred, log_lik] = get_loss(ii-1, post_dist, X, y, gamma, X_te, y_te);
    nlz =-log_lik;
    log_loss = pred.log_loss;

    Sigma = post_dist.covMat;
    mu = post_dist.mean;
end

author: emtiyaz

repo: vadam

save_path: github-repos/MATLAB/emtiyaz-vadam