text:
tic;
ssub = neuron_raw.options.ssub;
tsub = neuron_raw.options.tsub;

sframe = str2double(get(edit_begin, 'string'));
eframe = str2double(get(edit_end, 'string'));
num2read = str2double(get(edit_total, 'string'));

if and(ssub==1, tsub==1)
    neuron = neuron_raw;
    Y = double(data.Y(:, :, sframe+(1:num2read)-1));
    [d1s,d2s, T] = size(Y);
    fprintf('\nThe data has been loaded into RAM. It has %d X %d pixels X %d frames. \nLoading all data requires %.2f GB RAM\n\n', d1s, d2s, T, d1s*d2s*T*8/(2^30));
else
    [Y, neuron] = neuron_raw.load_data(nam_mat, sframe, num2read);
    [d1s,d2s, T] = size(Y);
    fprintf('\nThe data has been downsampled and loaded into RAM. It has %d X %d pixels X %d frames. \nLoading all data requires %.2f GB RAM\n\n', d1s, d2s, T, d1s*d2s*T*8/(2^30));
end
Y = neuron.reshape(Y, 1);
neuron_raw.P.p = 2;      %order of AR model

fprintf('Time cost in downsapling data:     %.2f seconds\n', toc);

author: zhoupc


text:
function n=countTestNodes(obj)
% Return the number of test nodes
%
% c=countTestCases(obj)
%
% Inputs:
%   obj             MoxUnitTestSuite instance
%
% Output:
%   c               Number of test nodes in obj
%
% NNO 2015

    n=obj.test_count;

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function M = Qdec2num(Qdec)
% M = Qdec2num(Qdec)
%
% Attempts to convert the two dimensional cell string array Qdec to numeric
% matrix M.
%
% Input
% Qdec: Two dimensional cell string array of Qdec data (eg. read with
% fReadQdec).
%
% Output
% M: Numeric matrix (without the name of the variables).
%
% $Revision: 1.2 $  $Date: 2015/01/06 17:14:53 $
% Original Author: Jorge Luis Bernal Rusiel
% CVS Revision Info:
%    $Author: mreuter $
%    $Date: 2015/01/06 17:14:53 $
%    $Revision: 1.2 $
%
Qdec = Qdec(2:end,:);
M = str2double(Qdec);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


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

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

% Last Modified by GUIDE v2.5 20-Apr-2004 10:07:56

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

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

% --- Executes just before TabPanel is made visible.
function TabPanel_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to TabPanel (see VARARGIN)

% Choose default command line output for TabPanel
handles.output = hObject;

% uicontrols handles list for tabA and tabB respectively
handles.tabA = [ handles.txNullA, handles.pbNullA, handles.axDisplay ];
handles.tabB = [ handles.txNullB, handles.pbNullB, handles.pbTest, ...
                 handles.chTest,  handles.lbTest, handles.rbTest ];

% Update handles structure
guidata(hObject, handles);

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

% --- Outputs from this function are returned to the command line.
function varargout = TabPanel_OutputFcn(hObject, eventdata, handles)
% Get default command line output from handles structure
varargout{1} = handles.output;

% --- If Enable == 'on', executes on mouse press in 5 pixel border.
% --- Otherwise, executes on mouse press in 5 pixel border or over tabA.
function tabA_ButtonDownFcn(hObject, eventdata, handles)
%
set(handles.tabA, 'Visible', 'on');
set(handles.tabB, 'Visible', 'off');

% --- If Enable == 'on', executes on mouse press in 5 pixel border.
% --- Otherwise, executes on mouse press in 5 pixel border or over tabB.
function tabB_ButtonDownFcn(hObject, eventdata, handles)
%

set(handles.tabA, 'Visible', 'off');
set(handles.tabB, 'Visible', 'on');

% --- Executes on selection change in lbTest.
function lbTest_Callback(hObject, eventdata, handles)

% --- Executes on button press in chTest.
function chTest_Callback(hObject, eventdata, handles)

% --- Executes on button press in pbTest.
function pbTest_Callback(hObject, eventdata, handles)

oStr = get( handles.lbTest, 'String');
nStr = 'Thanx. I needed it!';
lbStr = strvcat(oStr, nStr);
set( handles.lbTest, 'String', lbStr);

% --- Executes on button press in rbTest.
function rbTest_Callback(hObject, eventdata, handles)

% --- Executes on button press in pbNullB.
function pbNullB_Callback(hObject, eventdata, handles)

author: Sable

repo: mcbench-benchmarks


text:
function [HF,psi,x,y] = feed_dist(F,z,TF,PF)

parameters_dist
c = length(z); % c includes all species
options = optimset('Display','iter','Diagnostics','on','Algorithm','interior-point');

% Bubble point calculation
xB = z'; % Bubble point requirement (together with PF)
yB = xB; % Initial guess
TB = TF; % Initial guess
x0 = [yB; TB];
r  = fmincon('1',x0,[],[],[],[],[],[],@gB,options,xB,PF,c);
yB = r(1:c);
TB = r(end);

% Dew point calculation
yD = z'; % Dew point requirement (together with PF)
xD = yD; % Initial guess
TD = TF; % Initial guess
x0 = [xD; TD];
r  = fmincon('1',x0,[],[],[],[],[],[],@gD,options,yD,PF,c);
xD = r(1:c);
TD = r(end);

if TF <= TB
    psi = 0;
    x   = xB;
    y   = yB;
    [~,~,~,~,~,HF,~,~] = prsrk(x,PF*1e5,TF+273.15,pc,Tc,w,k,cpig,DHf,DGf,'L',Thermo);
elseif TF >= TD
    psi = 1;
    x   = xD;
    y   = yD;
    [~,~,~,~,~,HF,~,~] = prsrk(y,PF*1e5,TF+273.15,pc,Tc,w,k,cpig,DHf,DGf,'V',Thermo);
else
    % Two-phase calculation
    x  = z';    % Initial guess
    y  = x;     % Initial guess
    V  = 0.5*F; % Initial guess
    L  = V;     % Initial guess
    x0 = [x; y; V; L]; d = [F; TF; PF; z'];
    r = fmincon('1',x0,[],[],[],[],[],[],@gF,options,d,c);
    x = r(1:c); y = r(c+1:2*c);
    [~,~,~,~,~,HL,~,~] = prsrk(x,PF*1e5,TF+273.15,pc,Tc,w,k,cpig,DHf,DGf,'L',Thermo);
    [~,~,~,~,~,HV,~,~] = prsrk(y,PF*1e5,TF+273.15,pc,Tc,w,k,cpig,DHf,DGf,'V',Thermo);
    V = r(end-1); L = r(end);
    HF = 1/F*(HV*V + HL*L); psi = r(end-1)/F;
end

function [c,ceq] = gB(x0,x,PF,c)

parameters_dist
y = x0(1:c);
T = x0(end);

[~,~,phiL,~,~,~,~,~] = prsrk(x,PF*1e5,T+273.15,pc,Tc,w,k,cpig,DHf,DGf,'L',Thermo);
[~,~,phiV,~,~,~,~,~] = prsrk(y,PF*1e5,T+273.15,pc,Tc,w,k,cpig,DHf,DGf,'V',Thermo);
K = phiL./phiV;

ceq = [y - K.*x; sum(y) - 1];
c   = [];

function [c,ceq] = gD(x0,y,PF,c)

parameters_dist
x = x0(1:c);
T = x0(end);

[~,~,phiL,~,~,~,~,~] = prsrk(x,PF*1e5,T+273.15,pc,Tc,w,k,cpig,DHf,DGf,'L',Thermo);
[~,~,phiV,~,~,~,~,~] = prsrk(y,PF*1e5,T+273.15,pc,Tc,w,k,cpig,DHf,DGf,'V',Thermo);
K = phiL./phiV;

ceq = [y - K.*x; sum(x) - 1];
c   = [];

function [c,ceq] = gF(x0,d,c)

parameters_dist
x = x0(1:c);
y = x0(c+1:2*c);
V = x0(end-1);
L = x0(end);

F  = d(1);
TF = d(2);
PF = d(3);
z  = d(4:end);

[~,~,phiL,~,~,~,~,~] = prsrk(x,PF*1e5,TF+273.15,pc,Tc,w,k,cpig,DHf,DGf,'L',Thermo);
[~,~,phiV,~,~,~,~,~] = prsrk(y,PF*1e5,TF+273.15,pc,Tc,w,k,cpig,DHf,DGf,'V',Thermo);
K = phiL./phiV;

ceq = [y - K.*x; F*z - V*y - L*x; sum(x) - 1; sum(y) - 1];
c   = [];


text:
function numgrad = computeNumericalGradient(J, theta)
% numgrad = computeNumericalGradient(J, theta)
% theta: a vector of parameters
% J: a function that outputs a real-number. Calling y = J(theta) will return the
% function value at theta.

% Initialize numgrad with zeros   2x1  这是输入X的行数，其实也就是参数W的个数
numgrad = zeros(size(theta));

%% ---------- YOUR CODE HERE --------------------------------------
% Instructions:
% Implement numerical gradient checking, and return the result in numgrad.
% (See Section 2.3 of the lecture notes.)
% You should write code so that numgrad(i) is (the numerical approximation to) the
% partial derivative of J with respect to the i-th input argument, evaluated at theta.
% I.e., numgrad(i) should be the (approximately) the partial derivative of J with
% respect to theta(i).
%
% Hint: You will probably want to compute the elements of numgrad one at a time.

epsilon = 0.000001;
% 2 x2  得到参数的个数
numW = size(theta,1);
I = eye(numW);
I =  I * epsilon;
for j = 1:numW
    numgrad(j) = (J(theta+I(:,j)) - J(theta-I(:,j))) / (2 * epsilon);
end

%% ---------------------------------------------------------------
end

repo: MachineLearning


text:
%LandmarkMap Map of planar point landmarks
%
% A LandmarkMap object represents a square 2D environment with a number of landmark
% landmark points.
%
% Methods::
%   plot      Plot the landmark map
%   landmark   Return a specified map landmark
%   display   Display map parameters in human readable form
%   char      Convert map parameters to human readable string
%
% Properties::
%   map         Matrix of map landmark coordinates 2xN
%   dim         The dimensions of the map region x,y in [-dim,dim]
%   nlandmarks   The number of map landmarks N
%
% Examples::
%
% To create a map for an area where X and Y are in the range -10 to +10 metres
% and with 50 random landmark points
%        map = LandmarkMap(50, 10);
% which can be displayed by
%        map.plot();
%
% Reference::
%
%   Robotics, Vision & Control, Chap 6,
%   Peter Corke,
%   Springer 2011
%
% See also RangeBearingSensor, EKF.

% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB 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.
%
% RTB 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 Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

classdef LandmarkMap < handle
% TODO:
% add a name property, show in char()

    properties
        map    % map landmarks
        dim         % map dimension
        nlandmarks   % number of landmarks in map

        verbose
    end

    methods

        % constructor
        function map = LandmarkMap(nlandmarks, varargin)
        %LandmarkMap.LandmarkMap Create a map of point landmark landmarks
        %
        % M = LandmarkMap(N, DIM, OPTIONS) is a LandmarkMap object that represents N random point landmarks
        % in a planar region bounded by +/-DIM in the x- and y-directions.
        %
        % Options::
        % 'verbose'    Be verbose

        %% TODO: dim can be a 4-vector

            opt = [];
            [opt,args] = tb_optparse(opt, varargin);
            map.verbose = opt.verbose;

            if ~isempty(args) && isnumeric(args{1})
                dim = args{1};
            else
                dim = 10;
            end
            map.dim = dim;
            map.nlandmarks = nlandmarks;
            map.map = dim * (2*rand(2, nlandmarks)-1);
            map.verbose = false;
        end

        function f = landmark(map, k)
            %LandmarkMap.landmark Get landmarks from map
            %
            % F = M.landmark(K) is the coordinate (2x1) of the K'th landmark (landmark).
            f = map.map(:,k);
        end

        function plot(map, varargin)
            %LandmarkMap.plot Plot the map
            %
            % M.plot() plots the landmark map in the current figure, as a square
            % region with dimensions given by the M.dim property.  Each landmark
            % is marked by a black diamond.
            %
            % M.plot(LS) as above, but the arguments LS
            % are passed to plot and override the default marker style.
            %
            % Notes::
            % - The plot is left with HOLD ON.
            clf
            d = map.dim;
            axis equal
            axis([-d d -d d]);
            xlabel('x');
            ylabel('y');

            if nargin == 1
                args = {'kh'};
            else
                args = varargin;
            end
            h = plot(map.map(1,:)', map.map(2,:)', args{:});
            set(h, 'Tag', 'map');
            grid on
            hold on
        end

        function show(map, varargin)
        %map.SHOW Show the landmark map
        %
        % Notes::
        % - Deprecated, use plot method.
            warning('show method is deprecated, use plot() instead');
            map.plot(varargin{:});
        end

        function verbosity(map, v)
            %map.verbosity Set verbosity
            %
            % M.verbosity(V) set verbosity to V, where 0 is silent and greater
            % values display more information.
            map.verbose = v;
        end

        function display(map)
            %map.display Display map parameters
            %
            % M.display() displays map parameters in a compact
            % human readable form.
            %
            % Notes::
            % - This method is invoked implicitly at the command line when the result
            %   of an expression is a LandmarkMap object and the command has no trailing
            %   semicolon.
            %
            % See also map.char.
            loose = strcmp( get(0, 'FormatSpacing'), 'loose');
            if loose
                disp(' ');
            end
            disp([inputname(1), ' = '])
            disp( char(map) );
        end % display()

        function s = char(map)
        %map.char Convert map parameters to a string
        %
        % s = M.char() is a string showing map parameters in
        % a compact human readable format.
            s = 'LandmarkMap object';
            s = char(s, sprintf('  %d landmarks', map.nlandmarks));
            s = char(s, sprintf('  dimension %.1f', map.dim));
        end

    end % method

end % classdef

author: petercorke

repo: robotics-toolbox-matlab


text:
function model = gmm_estimate(data)
% gmm_estimate Estimates a GMM on a set of points
%
% Originally a part of: Maggot (developed within EU project CogX)
% Original author: Matej Kristan, 2009
%
% Input:
% - data (matrix): Points for which to estimate a model
%
% Output:
% - model (vector): values for corresponding points
%

N = length(data);
model.Mu = data;
model.Cov{N} = cell(N, 1);
for i = 1:N
    model.Cov{i} = 0;
end
model.w = ones(1,N) / N;

pdf0 = model;

% first we'll spherize the distribution
[Mu, C] = spherize(pdf0.Mu, pdf0.Cov, pdf0.w) ;
[U, S, V] = svd(C) ;
T = (diag(1./sqrt(diag(S))))*U' ;

pdf0.Mu = bsxfun(@minus, pdf0.Mu, Mu) ;
pdf0.Mu = T * pdf0.Mu ;
for i = 1 : length(pdf0.w)
    pdf0.Cov{i} = T*pdf0.Cov{i}*T' ;
end

C = T*C*T' ;
% calculate the optimal bandwidth by Kristan's estimator
H = optimal_bandwidth( pdf0.Mu, pdf0.Cov, pdf0.w, C, length(pdf0.w)) ;

pdf.Mu = Mu ;
pdf.Cov = {H} ;
pdf.w = 1 ;

iT = inv(T) ;
pdf.Mu = iT * pdf.Mu ;
pdf.Mu = bsxfun(@minus, pdf.Mu, -Mu) ;
for i = 1 : length(pdf.w)
    pdf.Cov{i} = iT * pdf.Cov{i} * iT';
end
H = pdf.Cov{1} ;

for i = 1:N
    model.Cov{i} = H ;
end

end

function [new_mu, new_Cov] = spherize(Mu, Cov, w)

if length(w)==1
    new_mu = Mu ;
    if ~isempty(Cov)
        new_Cov = Cov{1} ;
    else
        new_Cov = zeros(size(new_mu,1),size(new_mu,1)) ;
    end
    return ;
end

sumw = sum(w) ;
w = w / sumw ;
new_mu =  sum(bsxfun(@times,Mu, w),2) ;

n = size(new_mu,1) ;

if ~isempty(Cov)
    new_Cov = zeros(n,n) ;
    if n==1
        new_Cov = sum(w.*(cell2mat(Cov) + Mu.*Mu)) ;
    else
        for j=1:length(w)
            new_Cov = new_Cov + w(j)*( Cov{j} + Mu(:,j)*Mu(:,j)') ;
        end
    end
    new_Cov = new_Cov - new_mu*new_mu' ;
else
    new_Cov = zeros(n,n) ;
    if n==1
        new_Cov = sum(w.*(Mu.*Mu)) ;
    else
        for j=1:length(w)
            new_Cov = new_Cov + w(j)*(   Mu(:,j)*Mu(:,j)') ;
        end
    end
    new_Cov = new_Cov - new_mu*new_mu' ;

end

end

function [H] = optimal_bandwidth( Mu, Cov, w, Cov_smp, N_eff )

d = size(Mu,1) ;
G = (Cov_smp *(4 / ((d + 2) * N_eff))^(2 / (d + 4)));

alpha_scale = 1 ;
F = Cov_smp * alpha_scale; % for numerical stability. it could have been: F = Cov_smp;
% could also constrain to say that F = identity!
Rf2 = integral_squared_hessian(Mu, w, Cov, F, G);

h_amise = (N_eff ^ (-1) * det(F)^(-1 / 2) /( sqrt(4 * pi) ^ d * Rf2 * d ))^(1 / (d + 4)) ;
H = (F * h_amise ^ 2) * alpha_scale ;

end

function I = integral_squared_hessian(Mu, w, Cov, F, G)
% Calculates an integral over the squared Hessian of a Gaussian mixture
% model.
% Follows Wand and Jones "Kernel Smoothing", page 101., assuming H=h*F.

I = NaN ;
if ( isempty(Mu) )
    return;
end
% read dimension and number of components
[ d, N ]= size(Mu) ;

% precompute normalizer constNorm = ( 1 / 2pi)^(d/2)
constNorm = (1/(2*pi))^(d/2) ;
I = 0 ;

% test if F is identity for speedup
delta_F = sum(sum(abs(F-eye(size(F))))) ;
if delta_F < 1e-3
    % generate a summation over the nonsymmetric matrix
    for l1 = 1 : N
        S1 = Cov{l1}  + G ;
        Mu1 = Mu(:,l1) ;
        w1 = w(l1) ;
        for l2 = l1 : N
            S2 = Cov{l2};
            Mu2 = Mu(:,l2) ;
            w2 = w(l2) ;
            A = inv(S1 + S2) ;
            dm = (Mu1 - Mu2) ;
            m = dm'*A*dm ;
            f_t = constNorm*sqrt(det(A))*exp(-0.5*m) ;
            c = 2*sum(sum(A.*A'))*(1-2*m) + (1-m)^2 *trace(A)^2 ;

            % determine the weight of the term current
            if ( l1 == l2 )
                eta = 1 ;
            else
                eta = 2 ;
            end
            I = I + f_t*c*w2*w1*eta ;
        end
    end
else
    % generate a summation over the nonsymmetric matrix
    for l1 = 1 : N
        S1 = Cov{l1} ;
        Mu1 = Mu(:,l1) ;
        w1 = w(l1) ;
        for l2 = l1 : N
            S2 = Cov{l2} + G;
            Mu2 = Mu(:,l2) ;
            w2 = w(l2) ;
            A = inv(S1 + S2) ;
            dm = (Mu1 - Mu2) ;
            ds = dm'*A ;
            b = ds'*ds ;
            B = A - 2*b ;
            C = A - b ;

            f_t = constNorm*sqrt(det(A))*exp(-0.5*ds*dm) ;
            c = 2*trace(F*A*F*B) + trace(F*C)^2 ;

            % determine the weight of the term current
            if ( l1 == l2 )
                eta = 1 ;
            else
                eta = 2 ;
            end
            I = I + f_t*c*w2*w1*eta ;
        end
    end
end

end

author: votchallenge

repo: toolkit-legacy


text:
% Filename: ora_foc.m
% Oustaloup-Recursive-Approximation for fractional order differentiator
% Reference
% [1] Oustaloup, A.; Levron, F.; Mathieu, B.; Nanot, F.M.; "Frequency-band
%     complex noninteger differentiator: characterization and synthesis".
%     EEE Transactions on Circuits and Systems I: Fundamental Theory and
%     Applications, I , Volume: 47 Issue: 1, Jan. 2000, Page(s): 25 -39
% [2] D. Xue, Y.Q. Chen and D. Atherton. "Linear feedback control -
%     analysis and design with Matlab 6.5". Textbook draft. Chapter 9
%     "Fractional order control - An introduction". PDF available upon
%     request. Send request email to "yqchen@ieee.org".
% by YangQuan Chen. Nov. 2001.
% Utah State University. http://www.csois.usu.edu/people/yqchen/
% FOC web http://mechatronics.ece.usu.edu/foc/
%
% Input variables:
%       r: the fractional order as in s^r, r is a real number
%       N: order of the finite TF approximation for both (num/den)
%           (Note: 2N+1 recursive z-p pairs)
%       w_L: low frequency limit of the range of the frequency of interest
%       w_H: upper freq. limit of the range of the frequency of interest
% Output:
%       sys_foc: continuous time transfer function system object (TF)
% Sample values: w_L=0.1;w_H=1000;  r=0.5;  N=4;
% Existing problem: Be careful when doing "c2d", I met some problems.
function [sys_foc]=ora_foc(r,N,w_L,w_H)
w_L=w_L*0.1;w_H=w_H*10; % enlarge the freq. range of interest for goodness
mu=w_H/w_L;                %
w_u=sqrt(w_H*w_L);
alpha=(mu)^(r/(2*N+1));
eta=(mu)^((1-r)/(2*N+1));
k=-N:N;
w_kp=(mu).^( (k+N+0.5-0.5*r)/(2*N+1) )*w_L;
w_k=(mu).^( (k+N+0.5+0.5*r)/(2*N+1) )*w_L;
D_N_K=(w_u/w_H)^r * prod(w_k) / prod(w_kp);
D_N_P=-w_k;D_N_Z=-w_kp;
[num,den]=zp2tf(D_N_Z',D_N_P',D_N_K);
sys_foc=tf(num,den);


text:
function im = ibwt(decomp)
% function im = ibwt(decomp)
%
% BWT reconstruction
% Version 1.2
%
% Arguments:
%  decomp: A square BWT decomposition, as produced by bwt.m
%
% Result:
%  im: An image
%
% Citation:
%  Willmore B, Prenger RJ, Wu MC and Gallant JL (2008). The Berkeley
%  Wavelet Transform: A biologically-inspired orthogonal wavelet transform.
%  Neural Computation 20:6, 1537-1564
%
% The article is available at:
%  <http://dx.doi.org/10.1162/neco.2007.05-07-513>
%
% Copyright (c) 2008 Ben Willmore
%
% 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.

sz = size(decomp);

if (length(sz) ~= 2) || (sz(1) ~= sz(2))
  disp('Input must be square');
  im = nan;
  return;
end

sz = sz(1);

numlevels = log(sz)/log(3);

if ( (numlevels-floor(numlevels)) > abs(numlevels)*eps )
  disp('Input side length must be a power of 3');
  im = nan;
  return;
end

im = zeros(sz);

for level = 1:numlevels
  ssz = size(decomp,1);
  decomp_thislevel = decomp;

  if (ssz>3)
    decomp_thislevel(end-ssz/3+1:end,1:ssz/3) = 0;
  end

  im_thislevel = ibwt_onelevel(decomp_thislevel)/sz*ssz*((3^level)^2);
  idx = ceil((1:sz)*ssz/sz);
  im = im + im_thislevel(idx,idx); % upsample by sz/ssz times

  if (ssz>1)
    decomp = decomp(end-ssz/3+1:end,1:ssz/3);
  end
end


text:
function [vertices, faces] = curveToMesh(curve, varargin)
%CURVETOMESH  Create a mesh surrounding a 3D curve
%
%   [V, F] = curveToMesh(CURVE)
%   Computes the vertices and the faces of the mesh surrounding the
%   specified 3D curve.
%
%   [V, F] = curveToMesh(CURVE, THICKNESS)
%   Specifies the thickness of the mesh (distance between mesh vertices and
%   curve vertices). Default is 0.5.
%
%   [V, F] = curveToMesh(CURVE, THICKNESS, NCORNERS)
%   Also specifies the number of mesh vertices around each curve vertex.
%   Default is 8.
%
%
%   Example
%     % Creates a tubular mesh around a trefoil knot curve
%     t = linspace(0, 2*pi, 200)';
%     x = sin(t) + 2 * sin(2 * t);
%     y = cos(t) - 2 * cos(2 * t);
%     z = -sin(3 * t);
%     curve = [x, y, z];
%     [v2, f2] = curveToMesh(curve, .5, 16);
%     figure;
%     drawMesh(v2, f2);
%     axis equal; view(3);
%     axis([-4 4 -4 4 -2 2]);
%
%   See also
%     meshes3d, torusMesh, surfToMesh

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2015-01-07,    using Matlab 8.4.0.150421 (R2014b)
% Copyright 2015 INRA - Cepia Software Platform.

radius = .1;
if nargin > 1
    radius = varargin{1};
end

nCorners = 8;
if nargin > 2
    nCorners = varargin{2};
end

nNodes = size(curve, 1);
nVerts = nNodes * nCorners;

vertices = zeros(nVerts, 3);

% create reference corners, that will be rotated and translated
t = linspace(0, 2*pi, nCorners + 1)';
t(end) = [];
baseCorners = radius * [cos(t) sin(t) zeros(size(t))];

for iNode = 1:nNodes
    % coordinate of current node
    node = curve(iNode, :);

    % compute local tangent vector
    iNext = mod(iNode, nNodes) + 1;
    tangentVector = normalizeVector3d(curve(iNext, :) - node);

    % convert to spherical coordinates
    [theta, phi, rho] = cart2sph2(tangentVector); %#ok<ASGLU>

    % apply transformation to place corners around current node
    rotY = createRotationOy(theta);
    rotZ = createRotationOz(phi);
    trans = createTranslation3d(node);
    transformMatrix = trans * rotZ * rotY;
    corners = transformPoint3d(baseCorners, transformMatrix);

    % concatenate with other corners
    vertices( (1:nCorners) + (iNode - 1) * nCorners, :) = corners;
end

% indices of vertices
inds = (1:nVerts)';
add1 = repmat([ones(nCorners-1, 1) ; 1-nCorners], nNodes, 1);

% generate faces
faces = [inds ...
    mod(inds + add1 - 1, nVerts) + 1 ...
    mod(inds + nCorners + add1 - 1, nVerts) + 1 ...
    mod(inds + nCorners - 1, nVerts) + 1];

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
% en_softth - soft threshold function for the Elastic-net regularization
%
% Copyright(c) 2009 Ryota Tomioka
% This software is distributed under the MIT license. See license.txt

function [vv,ss]=en_softth(vv,lambda,info)

n = size(vv,1);
theta  = info.theta;

if theta<1
  I=find(abs(vv)>lambda*theta);
  vv=sparse(I,1,(abs(vv(I))-lambda*theta).*sign(vv(I))/(1+lambda*(1-theta)),n,1);
else
  vv=l1_softth(vv,lambda,info);
end

ss=en_spec(vv,theta);

repo: REST


text:
function dataout = ft_annotate(cfg, datain)

% FT_ANNOTATE returns the same output data as the user has provided as input, but allows
% to add comments to that data structure. These comments are stored along with the other
% provenance information and can be displayed with FT_ANALYSISPIPELINE. Adding comments
% is especially useful if you have manually (i.e. in plain MATLAB) modified the data
% structure, whereby some provenance information is missing.
%
% Use as
%   outdata = ft_annotate(cfg, indata)
% where the input data structure can be any of the FieldTrip data structures and
% the configuration structure should contain
%   cfg.comment    = string
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output structure.
%
% See also FT_ANALYSISPIPELINE, FT_MATH

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

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar datain
ft_preamble provenance datain

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% ensure that the required options are present
cfg = ft_checkconfig(cfg, 'required', 'comment');

% simply copy the input data to the output data
dataout = datain;

% this line is meant to provide some information
fprintf('adding the comment: %s\n', cfg.comment);

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous datain
ft_postamble provenance dataout
ft_postamble history dataout
ft_postamble savevar dataout

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function test23
%TEST23 test cs_dmspy
%
% Example:
%   test23
% See also: testall

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

clear functions

randn ('state', 0) ;
rand ('state', 0) ;

clf

for trials = 1:1000

    % m = fix (100 * rand (1)) ;
    n = fix (100 * rand (1)) ;
    m = n ;
    % d = 0.1 * rand (1) ;
    d = rand (1) * 4 * max (m,n) / max (m*n,1) ;
    A = sprandn (m,n,d) ;
    % S = sprandn (m,m,d) + speye (m) ;

    if (~ispc)
        if (rand ( ) > .5)
            A = A + 1i * sprand (A) ;
        end
    end

    cs_dmspy (A) ;
    drawnow

    % pause
end


text:
classdef PTKSaveAirwayMesh < PTKPlugin
    % PTKSaveAirwayMesh. Plugin for saving an STL file containing the surface
    %     mesh of the airways, pruned to the segmental level
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Save airway<br>mesh'
        ToolTip = 'Segments the airway lumen down to the segmental bronchi, and then creates and saves an STL mesh of the lumen surface'
        Category = 'Airways'

        AllowResultsToBeCached = true
        AlwaysRunPlugin = true
        PluginType = 'ReplaceOverlay'
        HidePluginInDisplay = false
        FlattenPreviewImage = true
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 2
        GeneratePreview = false
    end

    methods (Static)
        function results = RunPlugin(dataset, reporting)
            airways_by_segment_pruned = dataset.GetResult('PTKAirwaysPrunedBySegment');
            results = airways_by_segment_pruned.PrunedSegmentsByLobeImage;
            results.ChangeRawImage(results.RawImage > 0);

            smoothing_size = 1;
            filename = 'PrunedAirwaysSurfaceMesh.stl';
            results_upsampled = results.Copy;
            results_upsampled.AddBorder(6);
            results_upsampled.DownsampleImage(0.25);

            % Specifies the coordinate system to use when saving out files
            coordinate_system = MimCoordinateSystem.DicomUntranslated;

            template_image = results;

            small_structures = true;
            dataset.SaveSurfaceMesh('PTKSaveAirwayMesh', 'Airway mesh', filename, 'Surface mesh of the segmented airways down to the level of the pulmonary segments' , results, smoothing_size, small_structures, coordinate_system, template_image);
        end

        function results = GenerateImageFromResults(results, image_templates, reporting)
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [model] = tapas_sem_hier_prepare_model(data, model, inference)
%%
%
% Input
%
% Output
%

% aponteeduardo@gmail.com
% copyright (C) 2016
%

model.graph{1}.htheta.model = ...
    tapas_sem_hier_prepare_ptheta(model.graph{1}.htheta.model);

nc = size(model.graph{2}.htheta.T, 2);
ns = numel(data);

if ~isstruct(model.graph{4}.htheta.y)
    model.graph{4}.htheta.y = struct();
end

if ~isfield(model.graph{4}.htheta.y, 'mu')
    % The prior is added after wards
    model.graph{4}.htheta.y.mu = 0 * model.graph{1}.htheta.model.njm' * ...
        model.graph{1}.htheta.model.mu;
end

np = numel(model.graph{4}.htheta.y.mu);

% Mean prior variance
mu_pe = model.graph{1}.htheta.model.njm' * model.graph{1}.htheta.model.pm;
% Assume that the variance of the variance is only max 50
sg_pe = 0.5 * mu_pe;

[alpha, beta] = tapas_gamma_moments_to_ab(mu_pe, sg_pe);

if ~isfield(model.graph{4}.htheta.y, 'alpha')
    model.graph{4}.htheta.y.alpha =  alpha; % Wide prior
end

if ~isfield(model.graph{4}.htheta.y, 'beta')
    model.graph{4}.htheta.y.beta = beta;
end

model.graph{4}.htheta.y.eta = ones(size(mu_pe));

end

author: translationalneuromodeling

repo: tapas


text:
function val = viewGetFlat(vw,param,varargin)
% Get data from various view structures
%
% This function is wrapped by viewGet. It should not be called by anything
% else other than viewGet.
%
% This function retrieves information from the view that relates to a
% specific component of the application.
%
% We assume that input comes to us already fixed and does not need to be
% formatted again.

if notDefined('vw'), vw = getCurView; end
if notDefined('param'), error('No parameter defined'); end

mrGlobals;
val = [];

switch param

    case 'graycoords'
        % 'graycoords' is also an alias for coords in volume/gray views:
        if ~isequal(vw.viewType, 'Flat')
            val = viewGet(vw, 'Coords');
            return
        end

        % Example usage for FLAT view:
        % val = viewGet(vw{1},'graycoords','left');
        if length(varargin) ~= 1,
            error('You must specify which hemisphere.');
        end
        hname = varargin{1};
        switch hname
            case 'left'
                val = vw.grayCoords{1};
            case 'right'
                val = vw.grayCoords{2};
            otherwise
                error('Bad hemisphere name');
        end
    case 'leftpath'
        if checkfields(vw,'leftPath'), val = vw.leftPath; end
    case 'rightpath'
        if checkfields(vw,'rightPath'), val = vw.rightPath; end
    case 'fliplr'
        if checkfields(vw,'flipLR'), val = vw.flipLR; end
    case 'imagerotation'
        if checkfields(vw,'rotateImageDegrees'),val = vw.rotateImageDegrees; end
    case 'hemifromcoords'
        if ~exist('varargin', 'var') || isempty(varargin)
            val = [];
            warning('vista:viewError','Need coords to determine hemisphere');
            return;
        else
            % get left and right nodes to compare to the coords
            l = viewGet(vw, 'allLeftNodes');
            r = viewGet(vw, 'allRightNodes');
            l = l(1:3, :);
            r = r(1:3, :);
            % get coords in proper orientation
            coords = varargin{1};
            if ~isequal(size(coords, 1), 3) && isequal(size(coords, 2), 3)
                coords = coords';
            end
            % nodes are [cor axi sag] but coords are [axi cor sag]
            x = coords(2, :);
            y = coords(1, :);
            z = coords(3, :);
            [~, right] =  intersectCols(single([x; y; z]), r);
            [~, left]  =  intersectCols(single([x; y; z]), l);
            val = nan(1,length(x));
            val(left) = 1;
            val(right) = 2;
            val(intersect(left , right)) = nan;
        end
    case 'roihemi'
        if ~exist('varargin', 'var') || isempty(varargin)
            coords = viewGet(vw, 'roiCoords');
        else
            coords = varargin{1};
        end
        hemi = viewGet(vw,  'hemiFromCoords', coords);
        hemi = round(nanmean(hemi));
        if hemi == 1, val = 'left' ;
        elseif hemi == 2, val = 'right';
        else val = []; warning('vista:viewError','Could not determine hemifield');
        end

    otherwise
        error('Unknown viewGet parameter');

end

return


text:
% PURPOSE: Gets the halfamp cutoff frequency value using filtfilt.m
%
% WARNING: For working with filtfilt.m function!
%
% filtfilt result has the following characteristics:
%       a) Zero-phase distortion
%       b) A filter transfer function, which equals the squared magnitude of the original filter transfer function
%           WARNING by JLC: THIS IMPLIES THAT THE CUTOFF FREQUENCY IS NOT AT -3dB ANYMORE, IT IS AT -6dB INSTEAD!!!
%       c) A filter order that is double the order of the filter specified by b and a
%
% *** This function is part of ERPLAB Toolbox ***
% Author: Javier Lopez-Calderon
% Center for Mind and Brain
% University of California, Davis,
% Davis, CA
% 2009

function frec3dB = halfamp(b, a, fs)
try
        [hf,f1] = freqz(b,a,10000,fs);
        hf2 = hf^2; % filtfilt has a transfer function, which equals the squared magnitude of the original filter transfer function.
        [v loc] = min(abs(0.707-abs(hf2))); % frequency at gain 70.7%
        frec3dB = f1(loc);
catch
        frec3dB = [];
end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function suborder_num = tetrahedron_nco_suborder_num ( rule )

%*****************************************************************************80
%
%% TETRAHEDRON_NCO_SUBORDER_NUM returns the number of suborders for an NCO rule.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Peter Silvester,
%    Symmetric Quadrature Formulae for Simplexes,
%    Mathematics of Computation,
%    Volume 24, Number 109, January 1970, pages 95-100.
%
%  Parameters:
%
%    Input, integer RULE, the index of the rule.
%
%    Output, integer SUBORDER_NUM, the number of suborders of the rule.
%
  if ( rule == 1 )
    suborder_num = 1;
  elseif ( rule == 2 )
    suborder_num = 1;
  elseif ( rule == 3 )
    suborder_num = 2;
  elseif ( rule == 4 )
    suborder_num = 3;
  elseif ( rule == 5 )
    suborder_num = 5;
  elseif ( rule == 6 )
    suborder_num = 6;
  elseif ( rule == 7 )
    suborder_num = 9;

  else

    suborder_num = -1;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TETRAHEDRON_NCO_SUBORDER_NUM - Fatal error!\n' );
    fprintf ( 1, '  Illegal RULE = %d\n', rule );
    error ( 'TETRAHEDRON_NCO_SUBORDER_NUM - Fatal error!\n' )

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function pass = test_std(pref)
% Test the chebfun3/max command.

if ( nargin < 1 )
    pref = chebfunpref;
end
tol = 1e3*pref.cheb3Prefs.chebfun3eps;

% Note that std(chebfun(@(x) x.^2)) = sqrt(4/45).

f1 = chebfun3(@(x,y,z) x.^2 + y.*z); % f1 has x-std equal to sqrt(4/45).
f2 = chebfun3(@(x,y,z) y.^2 + x.*z); % f2 has y-std equal to sqrt(4/45).
f3 = chebfun3(@(x,y,z) z.^2 + x.*y); % f3 has z-std equal to sqrt(4/45).

exact = chebfun2(@(x,y) sqrt(4/45)); % x-std of y.*z equals 0.

h1 = std(f1);
h2 = std(f1,[],1);
h3 = std(f2,[],2);
h4 = std(f3,[],3);

pass(1) = norm(h1 - exact) < tol;
pass(2) = norm(h2 - exact) < tol;
pass(3) = norm(h3 - exact) < tol;
pass(4) = norm(h4 - exact) < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
classdef trigdouble < chebdouble
%FOURDOUBLE   Fourier double class.
%
%   See the CHEBDOUBLE class for details.
%
%   This class in intended solely as a worker-class for PDESOLVE.

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

    methods

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS CONSTRUCTOR:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        function obj = trigdouble(varargin)

            % Call the CHEBDOUBLE constructor:
            obj = obj@chebdouble(varargin{:});

        end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  DIFF  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        function u = diff(u, k)
            %DIFF   Compute the k-th derivative of u using Fourier
            % differentiation matrices defined by diffmat.

            % Store the diffmat D as a persistent variable to allow speeding up
            % if we work with the same discretization at multiple time steps.
            % Note that the matrix D is independent of the domain, since it is
            % scaled separately below.
            persistent D

            % Assume first-order derivative:
            if ( nargin == 1 )
                k = 1;
            end

            N = length(u.values);

            % Construct D if we don't match a previous discretization:
            if ( isempty(D) || numel(D) < k || size(D{k}, 1) ~= N )
                D{k} = trigcolloc.diffmat(N, k); % Diffmat
            end

            % Interval scaling. (Note: trigtech.diffmat is defined on [0, 2*pi))
            c = 2/diff(u.domain);

            % Muliplying by the kth-order differentiation matrix:
            u.values = c^k*(D{k}*u.values);

            % Update the difforder:
            u.diffOrder = u.diffOrder + k;

        end

        %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  SUM  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % The differential operators
        function I = sum(u, a, b)
            %SUM  Compute the integral of u.

            persistent W

            if ( nargin > 1 )
                % TODO: Add support.
                error('CHEBFUN:FOURDOUBLE:sum:notImplemented', ...
                    'Partial integrals not implemented yet.')
            end

            N = length(u.values);

            % Retrieve or compute weights:
            if ( N > 5 && numel(W) >= N && ~isempty(W{N}) )
                % Weights are already in storage!
            else
                c = diff(u.domain)/2; % Interval scaling.
                W{N} = c*trigtech.quadwts(N);
            end

            % Find the sum by muliplying by the weights vector:
            I = W{N}*u.values;

        end

        function I = integral(varargin)
            I = sum(varargin{:});
        end

        %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  CUMSUM  %%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % The differential operators
        function u = cumsum(u)
            %CUMSUM   Compute the indefinite integral of the Chebyshev
            %         interpolant to u.

            % TODO: Add support.
            error('CHEBFUN:FOURDOUBLE:cumsum:notImplemented', ...
                'CUMSUM not implemented yet.')

        end

        %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  FRED  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        function u = fred(K, u)
            %FRED  Fredholm operator with kernel K.
            %   FRED(K, U) computes the action of the Fredholm operator with
            %   kernel K on the Chebyshev interpolant to the points in the
            %   vector U.

            % TODO: Add support.
            error('CHEBFUN:FOURDOUBLE:fred:notImplemented', ...
                'FRED not implemented yet.')

        end

        %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  VOLT  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        function u = volt(K, u)
            %VOLT  Volterra operator with kernel K.
            %   VOLT(K, U) computes the action of the Volterra operator with
            %   kernel K on the Chebyshev interpolant to the points in the
            %   vector U.

            % TODO: Add support.
            error('CHEBFUN:FOURDOUBLE:volt:notImplemented', ...
                'VOLT not implemented yet.')

        end

        %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  FEVAL  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        function out = feval(u, y)
            %FEVAL  Evaluate polynomial interpolant of data {X_four, U} at a
            % point y using barycentric interpolation.

            persistent dom x

            n = length(u.values);
            udom = u.domain;

            if ( length(x) ~= n || isempty(dom) || ~all(dom == udom) )
                x = trigpts(n, dom);
                dom = udom;
            end

            out = trigBary(y, u.values, x, udom);

        end

    end

end


text:
function n = poisson_fixed_time ( lambda, time )

%*****************************************************************************80
%
%% POISSON_FIXED_TIME counts the Poisson events in a fied time.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real LAMBDA, the average number of events per unit time.
%
%    Input, real TIME, the amount of time to observe.
%
%    Output, integer N, the number of Poisson events observed.
%
  n = 0;
  t = 0.0;

  while ( t < time )
    dt = - log ( rand ( 1, 1 ) ) / lambda;
    n = n + 1;
    t = t + dt;
  end

  return
end


text:
% testThetaChange
cnnConfig = config();
[Theta, meta] = cnnInitParams(cnnConfig);
newTheta = thetaChange(Theta, meta, 'vec2stack', cnnConfig);
oldTheta = thetaChange(newTheta, meta, 'stack2vec', cnnConfig);
newOldTheta = thetaChange(oldTheta, meta, 'vec2stack', cnnConfig);

if isequal(newTheta,newOldTheta);
    fprintf('The thetaConfig function works well, congratulations\n');
else
    fprintf('The thetaConfig function does not work well, please check!\n');
end

if isequal(Theta,oldTheta);
    fprintf('The thetaConfig function works well, congratulations\n');
else
    fprintf('The thetaConfig function does not work well, please check!\n');
end

author: xuzhenqi

repo: cnn

save_path: github-repos/MATLAB/xuzhenqi-cnn


text:
% Main Script
close all
PsiPsi = [];
PsiU = [];
Phi = [];
CostQ = [];
x_save = []; % Keep track of the states
t_save = []; % Keep track of the time

currentTime = 0; % Indicating the current time throught the simulation

PhiLength = numel(Phi_fun(zeros(1,4)));
PsiLength = numel(Psi_fun(zeros(1,4)));

r = 1;  % Weight on u

N = min(80,2*(PsiLength+PsiLength));
IterMax = 20;  % Number of iterations

X = [0.1,-5,0.2,2,zeros(1,PsiLength^2+PsiLength+1)];

T = 0.01; % Length for each time interval

%% Online Data Collection
for i = 1:N
	[t,X] = ode45(@adpSysWrapper, ...
		[i-1,i]*T, ...
		[X(end,1:4) zeros(1, PsiLength^2 + PsiLength + 1)]);
	Phi    = [Phi;
		Phi_fun(X(end,1:4))-Phi_fun(X(1,1:4))];        %#ok<AGROW>
	PsiPsi = [PsiPsi;
		X(end,4 + (1:PsiLength^2))];                   %#ok<AGROW>
	PsiU   = [PsiU;
		X(end,4 + PsiLength^2 + (1:PsiLength))];       %#ok<AGROW>
	CostQ = [CostQ;
		X(end,end)];			                         %#ok<AGROW>
	t_save = [t_save;
		t(:)];                                         %#ok<AGROW>
	x_save = [x_save;
		X(:,1:4)];                                     %#ok<AGROW>
	currentTime = t_save(end);
end

% Off-Policy Learning. Solve the matrix A*pw = B

w = zeros(PsiLength,1);

for i = 1:IterMax
	A = [Phi -2*r*PsiU-2*r*PsiPsi*kron(w,eye(PsiLength))];
	% Note: To be consistent with the notations in (Y Jiang &  ZP Jiang,
	% TNNLS 2014), the above line should be replaced with
	% A = [Phi 2*r*PsiU-2*r*PsiPsi*kron(w,eye(PsiLength))];
	% i.e., consider u = w'*phi, not u = -w'*psi
	% In that case, the simpleSysWrapper should be called as
	% >> simpleSysWrapper(t,x,-w)
	B = -(CostQ + PsiPsi*kron(w,w));
	pw = A\B;
	p = pw(1:PhiLength);
	w = pw(PhiLength+1:end);
	if i == 1
		p0 = p; % save the initial value function
	else
		if norm(p-pp) <= 0.01
			Iter = i;
			break
		end
	end
	pp = p; %save the previous p for convergence check
end

%% Post learning
% Terminate exploration noise but keep appying the inital gains until the
% states enters the region of attraction

% Compute the region of attraction
D = getRegionOfAttraction(p);
% Keep checking if state is in the region of attraction
currentStates = x_save(end,:);
while  ~isInRegionOfAttraction(currentStates,p,D);
	[t,y] = ode45(@(t,x) simpleSysWrapper(t,x,w*0), ...
		currentTime+[0,0.1], ...
		currentStates);
	t_save = [t_save; t];
	x_save = [x_save; y];
	currentStates = x_save(end,:);
	currentTime = t_save(end);
end

%% Update controller and finish the rest of the simulation
[t,y] = ode45(@(t,x) simpleSysWrapper(t,x,-w), ...
	currentTime+[0,5], ...
	currentStates);
t_save = [t_save
	t(:)];
x_save = [x_save
	y(:,1:4)];

% Also compare with unlearned performance
[t0,y0] = ode45(@(t,x) simpleSysWrapper(t,x,w*0), ...
	currentTime+[0,5], ...
	currentStates);

%% Plotting results
% Create figure folder
if exist('simFigures','dir') == 0
	mkdir('simFigures');
end

%% Figure 1 Time course
figure(1)
% subplot(221)
plot(t_save,x_save(:,1),t0, y0(:,1), 'r--', 'LineWidth', 2)
xlabel('Time (sec)')
legend1 = legend('x_1 (Under ADP)', 'x_1 (Unlearned)');
set(legend1,'FontSize',12);
xlim([0,5])
% %
% subplot(222)
% plot(t_save,x_save(:,2),t0, y0(:,2), 'r--', 'LineWidth', 2)
% xlabel('Time (sec)')
% legend('x_b (Under ADP)', 'd/dt x_b (Unlearned)')
% %
% subplot(223)
% plot(t_save,x_save(:,3),t0, y0(:,3), 'r--', 'LineWidth', 2)
% xlabel('Time (sec)')
% legend('x_b (Under ADP)', 'x_w (Unlearned)')
% %
% subplot(224)
% plot(t_save,x_save(:,4),t0, y0(:,4), 'r--', 'LineWidth', 2)
% xlabel('Time (sec)')
% legend('x_b (Under ADP)', 'd/dt x_w (Unlearned)')

% Create textarrow
annotation('textarrow',[0.366326530612245 0.268112244897959],...
	[0.681349206349206 0.743253968253969],'String',{'Controller Updated'},...
	'FontSize',12);

print('.\simFigures\Ch3_ex1_fig1_x','-depsc');
%% Figure 2
% Plot the value function surfaces and compare the inital one and the
% optimized one
xxb = -0.5:0.05:0.5;
xxw = -0.2:0.02:0.2;
[XX,YY] = meshgrid(xxb, xxw);
VV = zeros(size(XX));
VV0 = VV;

for i = 1:numel(XX)
	VV0(i) = p0'*Phi_fun([XX(i),0,YY(i),0])';
	VV(i) = p'*Phi_fun([XX(i),0,YY(i),0])';
end

figure(2)
surf(XX,YY,VV0)
hold on
surf(XX,YY,VV)
hold off

xlabel('x_b', 'FontSize', 14);
ylabel('x_w', 'FontSize', 14);

% Create textarrow
annotation('textarrow',[0.821428571428571 0.889285714285713],...
	[0.821759259259259 0.652380952380954],'String','V_0(x_b,0,x_w,0)',...
	'FontSize',12);

% Create textarrow
annotation('textarrow',[0.833928571428571 0.868148148148148],...
	[0.157142857142857 0.336798336798337],'String', ...
	['V_' num2str(Iter-1) '(x_b,0,x_w,0)'],...
	'FontSize',12);

print('.\simFigures\Ch3_ex1_fig2_v','-depsc')

%% Figure 3 -- Region of Attraction
% Note: This section will be skipped if you are running an order version of
% MATLAB which does not have the function delaunayTriangulation

SkipSection = false;
try
	delaunayTriangulation;
catch e
	if strcmp(e.identifier, 'MATLAB:UndefinedFunction');
		SkipSection = true;
		warning(['The Region of attraction will not be plotted because', ...
			' MATLAB version is too old'])
	end
end

if SkipSection == false
	indexToRemove = [];

	xxb = linspace(-0.5,0.5,30);
	xxw = linspace(-0.2,0.2,30);
	xxdb = linspace(-5,5,30);

	[XX,YY,ZZ] = meshgrid(xxb, xxdb, xxw);

	for i = 1:numel(XX)
		if ~isInRegionOfAttraction([XX(i) YY(i) ZZ(i) 0],p,D)
			indexToRemove = [indexToRemove i];
		end
	end
	XX(indexToRemove) = [];
	YY(indexToRemove) = [];
	ZZ(indexToRemove) = [];

	DT = delaunayTriangulation(XX(:),YY(:),ZZ(:));
	figure(3)
	k = convexHull(DT);
	faceColor  = [0.6875 0.8750 0.8984];
	trisurf(k,DT.Points(:,1),DT.Points(:,2),DT.Points(:,3), ...
		'FaceColor', faceColor, ...   %'EdgeAlpha', 0.2, ...
		'FaceAlpha', 0.3)
	xlim([-0.5 0.5])
	ylim([-5 5])
	zlim([-0.2 0.2])
	% Create textbox
	annotation('textbox',...
		[0.453 0.168 0.07 0.017],...
		'String','\Omega',...
		'LineStyle','none',...
		'FontSize',12,...
		'FitBoxToText','off');
	% Create textarrow
	annotation('textarrow',[0.735714285714286 0.68],...
		[0.85952380952381 0.72],'String',{ '\Omega_i estimate'},'FontSize',12);
	xlabel('x_b', 'FontSize', 14);
	ylabel('x''_b', 'FontSize', 14);
	zlabel('x_w', 'FontSize', 14);

	print('.\simFigures\Ch3_ex1_fig3_Omega','-depsc')
end

%% Save Numerical Results and Clean Up
% Uncomment the following code to save data

% if exist('simResults','dir') == 0
% 	mkdir('simResults')
% end
% save .\simResults\simResults.mat

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function W = Fold(W, dim, i)
    dim = circshift(dim, [1-i, 1-i]);
    W = shiftdim(reshape(W, dim), length(dim)+1-i);
end

author: andrewssobral


text:
function DisplayLocation(city, lat, lng)
%#event
    disp(['Latitude = ' num2str(lat) '; Longitude = ' num2str(lng)]);
    disp(['Grid center near: ' city]);


text:
function test_pull1751

% WALLTIME 00:20:00
% MEM 3gb
% DEPENDENCY ft_nirs_signalqualityindex

% load sample set of data
fname = dccnpath('/home/common/matlab/fieldtrip/data/test/original/nirs/artinis/Sofia/sample_data.oxy4');

cfg = [];
cfg.dataset = fname;
cfg.channel = 'nirs';
data_test = ft_preprocessing(cfg);

% load expected data output
load('/home/common/matlab/fieldtrip/data/test/sqi/output_sqi.mat');

% compute sqi
cfg = [];
data_out = ft_nirs_signalqualityindex(cfg, data_test);

% check whether it matches the expected output
assert(isequaln(data_out.trial, output_trials))

% check whether the size of the output matrix is consistent with the input one
assert(isequal(size(output_trials{1}),size(data_test.trial{1})))

%%%% check whether it raises error when there are problems with the labels
% case 1: channel labels do not match the expected order
data_non_ordered = data_test;
labels = data_non_ordered.label;
idx_rand = randperm(numel(labels));
labels_new = labels(idx_rand);

trials = data_non_ordered.trial;
trials_new{1, 1} = trials{1, 1}(idx_rand,:);

data_non_ordered.trial = trials_new;
data_non_ordered.label = labels_new;

try
    ft_nirs_signalqualityindex(cfg, data_non_ordered);
    raises_error = 0;
catch
    raises_error = 1;
end

if ~raises_error
    error('Does not raise error when labels do not match expected order.')
end

% case 2: more than one label associated to the same channel
data_wrong_num_labels = data_test;
labels = data_wrong_num_labels.label;
labels{3} = 'Rx1-Tx1 [900nm]';
data_wrong_num_labels.label = labels;

try
    ft_nirs_signalqualityindex(cfg, data_wrong_num_labels);
    raises_error = 0;
catch
    raises_error = 1;
end

if ~raises_error
    error('Does not raise error when channel combination is not associated with a number of wavelengths different than 2.')
end


text:
function [model, activation] = bp_forward(model, batch)
% backprop feed forward used for discriminative finetuning

global kConv_forward2;
global kConv_forward_c;

activation = cell(model.numLayer, 1);
for l = 2 : model.numLayer
    activation{l-1} = batch;
    if l == 2
        stride = model.layers{l}.stride;
        hidden_presigmoid = myConvolve2(kConv_forward2, batch, model.layers{l}.w, stride, 'forward');
        hidden_presigmoid = bsxfun(@plus, hidden_presigmoid, permute(model.layers{l}.c, [2,3,4,5,1]));
        batch = 1 ./ (1 + exp(-hidden_presigmoid));
    elseif strcmp(model.layers{l}.type, 'convolution')
        stride = model.layers{l}.stride;
        hidden_presigmoid = myConvolve(kConv_forward_c, batch, model.layers{l}.w, stride, 'forward');
        hidden_presigmoid = bsxfun(@plus, hidden_presigmoid, permute(model.layers{l}.c, [2,3,4,5,1]));
        batch = 1 ./ (1 + exp(-hidden_presigmoid));
    elseif strcmp(model.layers{l}.type, 'fullconnected') && l < model.numLayer
        batch_size = size(batch,1);
        batch = reshape(batch, batch_size, []);
        hidden_presigmoid = bsxfun(@plus, double(batch) * double(model.layers{l}.w), double(model.layers{l}.c));
        batch = 1 ./ (1 + exp(-hidden_presigmoid));
    else
        batch_size = size(batch,1);
        batch = reshape(batch, batch_size, []);
        temp = bsxfun(@plus, batch * model.layers{l}.w, model.layers{l}.c);
        energy = exp( bsxfun (@minus, temp, max(temp, [], 2)));
        batch = bsxfun(@rdivide, energy, sum(energy,2));
    end
end
activation{end} = batch;

author: zhirongw


text:
% this function write contents of "txt" on an image and then returns its image.
% "h" is the texyt height in pixels
% "figNum" is the number for figure.
% note that you may not interacti with matlab while running this command since the figure window should be on the top.

function bw = text_to_image(txt, h, figNum)
f = h/1.2;
sz1 = round(f*4);
sz2 = round(length(txt)*f*1.5);
r = 1;
if sz2 > 800
  r = sz2/800;
  sz2 = round(sz2/r);
  f = round(f/r);
end

im1 = ones(sz1, sz2, 3);
flag1 = 0;
if ~exist('figNum')
  figNum = ceil((rand*1000)+1);
  flag1 = 1;
end

figure(figNum);
imshow(im1);
text(10,f*3, txt, 'fontsize', f,'color', 'b');
im2 = getframe;
if flag1
  close(figNum);
end

im2 = im2.cdata;
bw2 = im2bw(im2);
f1 = find(sum(~bw2));
f2 = find(sum(~bw2'));
x1 = f1(1);
y1 = f2(1);
x2 = f1(end);
y2 = f2(end);
bw = ~bw2(y1:y2, x1:x2);

if r~=1
  bw = im2bw(imresize(double(bw), h/size(bw,1)));
end

author: VisDrone

repo: DroneCrowd

save_path: github-repos/MATLAB/VisDrone-DroneCrowd


text:
function d = phasecorr(A,B)
%PHASECORR Compute phase correlation matrix
%
%   D = PHASECORR(A, B) computes the phase correlation matrix from input
%   2-D images A and B. PHASECORR returns the 'full' correlation matrix
%   similar to normxcorr2.

%   Copyright 2013 The MathWorks, Inc.

size_A  = size(A);
size_B  = size(B);

% Let fft2 zero pad time domain signals such that we form the 'full'
% correlation matrix analogous to the result produced by normxcorr2.
outSize = size_A + size_B - 1;

% Form 2-D spectra of moving and fixed.
% Window each signal prior to taking FFT to help reduce ringing effects in
% frequency domain due to finite image.
A = fft2(A,outSize(1),outSize(2));
B = fft2(B,outSize(1),outSize(2));

% Form phase correlation matrix, d
% Use 'symmetric' option as performance optimization. We expect that the
% input moving images A and B are real valued, so d should always be real
% valued.
ABConj = A .* conj(B);
d = ifft2(ABConj ./ abs(eps+ABConj),'symmetric');

author: zhouyuanzxcv

repo: Hyperspectral


text:
function F = power( F, G )
%.^ Componentwise power for DISKFUNV.
%   F.^G where F is a DISKFUNV and G is a double returns the result from
%   componentwise powers.
%
%   F.^G where F is a double and G is a CHEBFUN2 returns from componentwise
%   powers.

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

% Empty check:
if ( isempty( F ) || isempty( G ) )
    F = diskfunv;
    return
end

if ( isa( F, 'double' ) )       % scalar . ^ DISKFUNV

    if ( numel(F) == 1 )

        scalar = F;
        F = G;
        F.components{1} = power( scalar, G.components{1} );
        F.components{2} = power( scalar, G.components{2} );

    else

        error('CHEBFUN:DISKFUNV:power:double', 'Dimension mismatch.');

    end

elseif ( isa(G, 'double') )      % DISKFUNV . ^ scalar

    if ( numel(G) == 1 )

        scalar = G;
        F.components{1} = power( F.components{1}, scalar );
        F.components{2} = power( F.components{2}, scalar );

    else

        error('CHEBFUN:DISKFUNV:power:double',...
                     'DISKFUNV and scalar dimension mismatch.');

    end

elseif ( isa(F,'diskfunv') && isa(G,'diskfunv') )  % DISKFUNV.^DISKFUNV

    error('CHEBFUN:DISKFUNV:power:size', 'DISKFUNV dimension mismatch.');

else  % error

    error('CHEBFUN:DISKFUNV:power:inputs', 'Unrecognized input arguments.');

end

end


text:
function Y=sum(Y)
%SUM (overloaded)

Y.basis = sum(Y.basis,1);
Y.dim(1) = 1;
Y.dim(2) = 1;
% Reset info about conic terms
Y.conicinfo = [0 0];
Y = clean(Y);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
clear all; close all; clc

% Connect to Arduino
tclab;

% Run time in minutes
run_time = 10;

% Number of cycles (1 cycle per 3 seconds)
loops = round(20*run_time);

% heater input
Q1 = zeros(1,loops);
Q2 = zeros(1,loops);
Q1(3:end) = 100.0;
Q1(50:end) = 0.0;
Q1(100:end) = 80.0;

Q2(25:end) = 60.0;
Q2(75:end) = 100.0;
Q2(125:end) = 25.0;

for i = 130:180
    if mod(i,10)==0
        Q1(i:i+10) = rand(1) * 100;
    end
    if mod(i+5,10)==0
        Q2(i:i+10) = rand(1) * 100;
    end
end

% Temperature (degC)
T1 = ones(1,loops) * T1C(); % measured T
T2 = ones(1,loops) * T2C(); % measured T
time = zeros(1,loops);

start_time = clock;
prev_time = start_time;

% dynamic plot (note: subplots needs to be declared here first)
figure(1)
subplot(2,1,1)
hold on, grid on
anexp1 = animatedline('LineStyle','-', 'Color', 'k', 'LineWidth', 2);
anexp2 = animatedline('LineStyle','-', 'Color', 'b', 'LineWidth', 2);
ylabel('Temperature \circC')
legend('T_1 Measured', 'T_2 Measured', ...
    'Location', 'northwest')
subplot(2,1,2)
hold on, grid on
anQ1 = animatedline('LineStyle','-', 'Color', 'k', 'LineWidth', 2);
anQ2 = animatedline('LineStyle','--', 'Color', 'b', 'LineWidth', 2);
ylabel('Power Level Q (%)')
legend('Q_1', 'Q_2', 'Location', 'northwest')
xlabel('Time (sec)')

for ii = 1:loops
    % adjust power level
    h1(Q1(ii));
    h2(Q2(ii));

    % Pause Sleep time
    pause_max = 3.0;
    pause_time = pause_max - etime(clock,prev_time);
    if pause_time >= 0.0
        pause(pause_time - 0.01)
    else
        pause(0.01)
    end

    % Record time and change in time
    t = clock;
    dt = etime(t,prev_time);
    if ii>=2
        time(ii) = time(ii-1) + dt;
    end
    prev_time = t;

    % read and record from temperature controller
    T1(ii) = T1C();
    T2(ii) = T2C();

    % plot
    addpoints(anexp1,time(ii),T1(ii))
    addpoints(anexp2,time(ii),T2(ii))
    addpoints(anQ1,time(ii),Q1(ii))
    addpoints(anQ2,time(ii),Q2(ii))
    drawnow
end

h1(0);
h2(0);
disp('Heaters off')
% turn off heater but keep LED on if T > 50
if (T1C() || T2C()) > 50
    led(1)
    disp(['Warning, heater temperature 1 =', num2str(T1C())])
    disp(['Warning, heater temperature 2 =', num2str(T2C())])
else
    led(0)
end

% save txt file with data
data = [time',Q1',Q2',T1',T2'];
csvwrite('data.txt',data);

% save as csv file with headers
fid = fopen('data.csv','w');
fid2 = fopen('data.txt','r');
fprintf(fid,'%s\n',['time,Q1,Q2,TC1,TC2']);
while ~feof(fid2)
    line = fgetl(fid2);
    fprintf(fid,'%s\n',line);
end
fclose(fid);
fclose(fid2);

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino