text:
%The purpose of this program is to implement Galerkin method over
% the entire domain for solving the following general 2nd order,
% homogeneous, Boundary Value problem (BVP) with constant coefficients, and
% then comparing the answer with the exact solution.
%==========================================================================
%               ax"(t)+bx'(t)+cx(t)=0 for t1<=t<=t2
%                   BC: x(t1)=x1 and x(t2)=x2
% >> BVP_Galerkin1(a,b,c,t1,t2,x1,x2,n)
% where "n" is the number of trial functions
% The output of this program is
% 1- The approximated solution of x(t) vs. exact solution of x(t)
% 2- The approximated x'(t) vs. exact x'(t)
% 3- The approximated x"(t) vs. exact x"(t)
% ======Example============================================================
% Equation                  x"(t)+ x'(t)+ x(t)=0
% Boundary values           x(1)=2,     x(10)=0;
% Solution: We have:
%                           a=1;    b=1;   c=1;
%                           t1=1;       t2=10;
%                           x1=2;       x2=0;
% Using n=8 tiral functions,
% >>BVP_Galerkin(1,2,3,1,10,2,0,8)
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%                   October.26, 2010, Ramin Shamshiri
%                       ramin.sh@ufl.edu
%            Doctoral student at the University of Florida
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

function BVP_Galerkin1(a,b,c,t1,t2,x1,x2,n) % Declare function
dt=0.001; % increment
%% Begin Approximate solution via Galerkin method over the entire domain

for i=1:n
    for j=1:n
    K1(i,j)=(i*j)*(t2^(j+i-1)-t1^(j+i-1))/(i+j-1);          % K1 matrix corresponding to the x"(t)
    K2(i,j)=j*((t2^(i+j))-(t1^(i+j)))/(i+j);                % K2 matrix corresponding to the x'(t)
    K3(i,j)=(t2^(i+j+1)-(t1^(i+j+1)))/(i+j+1);              % K3 matrix corresponding to the x(t)
    end
end
% Overall Matrix
K=-a*K1+b*K2+c*K3;

%Applying boundary conditions and Solving nxn system of equations
for i=1:n-1
     for j=1:n
         A(i,j)=K(i,j)+((j*a*t2^(j+i-1))-(j*a*t1^(j+i-1)));

     end
end

%Applying boundary conditions
for j=1:n
    A(n-1,j)=(t1^j);
    A(n,j)=(t2^j);
end

for i=1:n-1
    B(i,1)=0;
end
B(n-1,1)=x1;
B(n,1)=x2;

C=inv(A)*B; % Ci: Coefficients of terms

t=t1:dt:t2;
t=t';
% Creating vertical axis values, x(t)
x=  zeros((t2-t1)/dt+1,1);
for i=1:(t2-t1)/dt+1
    for j=1:n
        x(i,1)=x(i,1)+ C(j,1)*(t(i,1)^j);
    end
end

% calculating x'(t)
dx=zeros((t2-t1)/dt+1,1);
for i=2:(t2-t1)/dt+1
    dx(i,1)=(x(i)-x(i-1))/dt;
end
% calculating x"(t)
d2x=zeros((t2-t1)/dt+1,1);
for i=3:(t2-t1)/dt+1
    d2x(i,1)=(dx(i)-dx(i-1))/dt;
end

 figure1 = figure('Color',[1 1 1]);
% Plotting approximate solution via Galerkin method over the entire domain
subplot(3,1,1,'Parent',figure1);
% Plotting the Approximate solution in green (Galerkin method over the entire domain)
plot(t,x,'-.', 'Color','r', 'LineWidth',2), grid on; hold on;
% Plotting the approximate dx
subplot(3,1,2,'Parent',figure1);
plot(t,dx, '-.', 'Color','r', 'LineWidth',2), grid on; hold on;
% Plotting the approximate d2u
subplot(3,1,3,'Parent',figure1);
plot(t,d2x, '-.', 'Color','r', 'LineWidth',2), grid on; hold on;
%% End of approximate solution via Galerkin method over the entire domain

%% ==============Begin Calculating Exact Soution==========================
r=roots([a b c]);   % Roots of auxiliary equation
r1=r(1);  r2=r(2);

% Bulding t-axis values
t=t1:dt:t2; t=t';
s=size(t);

% Bulding x(t) axis values
x=zeros(s(1),1);
if r1==r2   % For critically damped case

    C=inv([exp(r1*t1)  t1*exp(r2*t1);exp(r1*t2)  t2*exp(r2*t2)])*([x1;x2]);
    C1=C(1);
    C2=C(2);

    for i=1:s(1)
        x(i)=C1*exp(r1*t(i))+C2*t(i)*exp(r2*t(i));
    end

else        % For Under-damped or Over-damped case

    C=inv([exp(r1*t1) exp(r2*t1);exp(r1*t2) exp(r2*t2)])*([x1;x2]);
    C1=C(1);
    C2=C(2);

    for i=1:s(1)
        x(i)=C1*exp(r1*t(i))+C2*exp(r2*t(i));
    end
end

% calculating x'(t)
dx=zeros(s(1),1);
for i=2:s(1)
    dx(i,1)=(x(i)-x(i-1))/dt;
end

% calculating x"(t)
d2x=zeros(s(1),1);
for i=3:s(1)
    d2x(i,1)=(dx(i)-dx(i-1))/dt;
end

% Plotting Exact solution
subplot(3,1,1,'Parent',figure1);
% Plotting the exact solution x(t) in blue
plot(t,x, 'Color','b', 'LineWidth',2);
xlabel('Time (sec)','FontSize',12);
ylabel('x(t) (m)','FontSize',12);
legend1 = legend(subplot(3,1,1),'show');
legend('Approximate x(t)','Exact x(t)')

% Plotting the exact solution x'(t) in blue
subplot(3,1,2,'Parent',figure1);
plot(t,dx, 'Color','b', 'LineWidth',2);
xlabel('Time (sec)','FontSize',12);
ylabel('dx/dt (m/s)','FontSize',12);
legend1 = legend(subplot(3,1,2),'show');
legend('Approximate dx(t)','Exact dx(t)')

% Plotting the exact solution x"(t) in blue
subplot(3,1,3,'Parent',figure1);
plot(t,d2x, 'Color','b', 'LineWidth',2);
xlabel('Time (sec)','FontSize',12);
ylabel('d^2x/dt^2 m/s^2','FontSize',12);
legend1 = legend(subplot(3,1,3),'show');
legend('Approximate d2x(t)','Exact dx(t)')
clc
%% =====================End of exact solution==============================


text:
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% /**
%  * Copyright (C) 2016 CoDyCo
%  * @author: Daniele Pucci, Gabriele Nava
%  * Permission is granted to copy, distribute, and/or modify this program
%  * under the terms of the GNU General Public License, version 2 or any
%  * later version published by the Free Software Foundation.
%  *
%  * A copy of the license can be found at
%  * http://www.robotcub.org/icub/license/gpl.txt
%  *
%  * 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
%  */
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% NOTE: THIS SCRIPT IS RUN AUTOMATICALLY WHEN THE USER STARTS THE ASSOCIATED
% SIMULINK MODEL. NO NEED TO RUN THIS SCRIPT EVERY TIME.
clearvars -except sl_synch_handles simulinkStaticGUI
clc

%% GENERAL SIMULATION INFO
%
% If you are simulating the robot with Gazebo, remember that it is required
% to launch Gazebo as follows:
%
%     gazebo -slibgazebo_yarp_clock.so
%
% and properly set the environmental variable YARP_ROBOT_NAME in the .bashrc.

% Simulation time in seconds. For long simulations, put an high number
% (not inf) for allowing automatic code generation
Config.SIMULATION_TIME = 600000;

% Controller period [s]
Config.tStep           = 0.01;

% If Config.SAVE_WORKSPACE = True, every time the simulink model is run the
% workspace is saved after stopping the simulation
Config.SAVE_WORKSPACE  = false;

% if TRUE, the controller will STOP if the joints hit the joints limits
% and/or if the (unsigned) difference between two consecutive joints
% encoders measurements is greater than a given threshold.
Config.EMERGENCY_STOP_WITH_JOINTS_LIMITS  = false;
Config.EMERGENCY_STOP_WITH_ENCODER_SPIKES = true;

% Verify that the integration time has been respected during the simulation
Config.CHECK_INTEGRATION_TIME = true;
Config.PLOT_INTEGRATION_TIME  = false;

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CONFIGURATIONS COMPLETED: loading gains and parameters for the specific robot
run(strcat('app/robots/',getenv('YARP_ROBOT_NAME'),'/configRobot.m'));
run(strcat('app/robots/',getenv('YARP_ROBOT_NAME'),'/configJointsControl.m'));
run(strcat('app/robots/',getenv('YARP_ROBOT_NAME'),'/gainsAndReferences.m'));

author: robotology

repo: whole-body-controllers


text:
function [Fconv,info] = nonlineartocone(F)
%NONLINEARTOCONE Convert nonlinear constraint to second order cone

if islinear(F)
    Fconv = F;
    info = 0;
    return
end

[n,m]=size(F);
if (n*m==1)
    Fconv = [];
    % Involved in polynomial constraints
    sqrList = yalmip('nonlinearvariables');
    h = sort(unique(sqrList));
    % % SDP relaxation
    for i = 1:size(sqrList,1)
        left_index = find(h==sqrList(i,1));
        x1_index = find(h==sqrList(i,2));
        x2_index = find(h==sqrList(i,3));
        Z{i}=zeros(length(h),length(h));
        if x1_index==x2_index
            Z{i}(x2_index,x1_index) = 1;
        else
            Z{i}(x1_index,x2_index) = 0.5;
            Z{i}(x2_index,x1_index) = 0.5;
        end
    end
    vars = getvariables(F);
    base = getbase(F);
    nonlinvars = find(ismember(vars,sqrList(:,1)));
    linvars = find(~ismember(vars,sqrList(:,1)));
    % Construct quadratic
    Q = zeros(length(h));
    for i = 1:length(nonlinvars)
        indexinlist = find(vars(nonlinvars(i))==sqrList(:,1));
        indexinh = find(vars(nonlinvars(i))==h);
        Q = Q-Z{indexinlist}*base(1+find(h(indexinh)==vars));
    end
    used = find(any(Q));
    [B,r] = chol(Q(used,used));
    if r==0
        linear = base(1);
        if ~isempty(linvars)
            linear = linear+base(1+find(ismember(linvars,vars)))*recover(linvars);
        end
        ConesAxb=[2*B*recover(h(used));1-linear];
        Conescxd=1+linear;
        Fconv = cone(ConesAxb,Conescxd);
        info = 0;
    else
        Fconv = F;
        if nargout == 2
            info = 1;
        else
            warning('Cannot re-write to second order cone constraint');
        end
    end
else
    if nargout == 2
        Fconv = F;
        info = 1;
    else
        error('nonlineartocone can only be applied to scalar inequalities')
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [featuresM,volToEval,mask3M,minr,maxr,minc,maxc,mins,maxs] = getLawsAndHaralickFeatures(structNum,...
    rowMargin,colMargin,slcMargin,minIntensity,maxIntensity,planC,varargin)

% --- Optional inputs ---
% varargin{1} : haralOnlyFlag (default : 1)
% varargin{2} : Vector of flags indicating haralick textures to be computed
%               (default : all ones)
%               flagV = [energyFlg,entropyFlag,sumAvgFlg,homogFlg,...
%               contrastFlg,corrFlg,clustShadFlg,clustPromFlg,haralCorrFlg];
% varargin{3} : No. gray levels   (default:32 )
% varargin{4} : patchRadius       (default:[1,2])
% ------------------------
% AI 9/1/17 Also returns volToEval, mask3M, boundingbox coords
% AI 9/7/17 Added optional inputs

if ~exist('planC','var')
    global planC
end
indexS = planC{end};

% Get optional inputs
minargs = 7;
maxargs = 8;
nOptArgs = nargin - minargs;
optC = {1,ones(1,9),32,[1,2]};
[optC{1:nOptArgs}] = varargin{:};
[haralOnlyFlag,flagV,nLevel,patchRadius] = optC{:};

% get the region of interest
if ~iscell(structNum)
    [volToEval,maskBoundingBox3M] = ...
        getROI(structNum,rowMargin,colMargin,slcMargin,planC);
else
    volToEval = structNum{1};
    maskBoundingBox3M = structNum{2};
end

nanIntenityV = volToEval < -400;

% Generate Law's texture
if ~haralOnlyFlag
    % padd with mean intensities
    meanVol = nanmean(volToEval(:));
    if exist('padarray.m','file')
        paddedVolM = padarray(volToEval,[5 5 5],meanVol,'both');
    else
        paddedVolM = padarray_oct(volToEval,[5 5 5],meanVol,'both');
    end
    lawsMasksS = getLawsMasks();

    fieldNamesC = fieldnames(lawsMasksS);
    numFeatures = length(fieldNamesC);
    % initialize features matrix
    featuresM = zeros(sum(maskBoundingBox3M(~nanIntenityV)),numFeatures);
    for i = 1:numFeatures
        disp(i)
        text3M = convn(paddedVolM,lawsMasksS.(fieldNamesC{i}),'same');
        text3M = text3M(6:end-5,6:end-5,6:end-5);
        % featuresM(:,i) = text3M(maskBoundingBox3M); % for non cubic roi
        featuresM(:,i) = text3M(~nanIntenityV); % for the entire cubic roi
    end

else
    featuresM = [];
end

% Intensity as a feature
featuresM(:,end+1) = volToEval(~nanIntenityV);

% Filter intensities
nanIntenityV = volToEval < -400;
%nanIntenityV = false(size(nanIntenityV));
volToEval(nanIntenityV) = NaN;
% volToEval(volToEval > maxIntensity) = NaN;

% Genarate Haralick textures
separateDirnFlag = 1;
nOff = 1;
if separateDirnFlag
    offsetsM = getOffsets(2);
    nOff = size(offsetsM,1);
    numVox = sum(~nanIntenityV(:));
    nanIntenityV = repmat(nanIntenityV,[1 1 1 nOff]);
end

%featuresM = zeros(sum(maskBoundingBox3M(:)),0);
for patchSiz = patchRadius

    patchSizeV  = [patchSiz patchSiz 0];

    for numGrLevels = nLevel
        offsetsM = getOffsets(2);

        % Haralick texture
        waitH = NaN;

        [energy3M,entropy3M,sumAvg3M,corr3M,invDiffMom3M,contrast3M,...
            clustShade3M,clustPromin3M,haralCorr3M] = gpuTextureByPatch(volToEval,...
            numGrLevels,patchSizeV,offsetsM,flagV,waitH,minIntensity,maxIntensity,separateDirnFlag);

        tempC = {energy3M,entropy3M,sumAvg3M,corr3M,invDiffMom3M,contrast3M,...
            clustShade3M,clustPromin3M,haralCorr3M};
        idx = find(flagV);
        for n = 1:length(idx)
            featuresM(:,end+1:end+nOff) = reshape(tempC{idx(n)}(~nanIntenityV),numVox,nOff);
        end

    end % gray level
end % patch radius

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function wavelet_test09 ( )

%*****************************************************************************80
%
%% WAVELET_TEST09 tests DAUB18_TRANSFORM and DAUB18_TRANSFORM_INVERSE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 July 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'WAVELET_TEST09\n' );
  fprintf ( 1, '  DAUB18_TRANSFORM computes the DAUB18 transform of a vector.\n' );
  fprintf ( 1, '  DAUB18_TRANSFORM_INVERSE inverts it.\n' );
%
%  Random data.
%
  n = 16;
  seed = 123456789;
  [ u, seed ] = r8vec_uniform_01 ( n, seed );

  v = daub18_transform ( n, u );

  w = daub18_transform_inverse ( n, v );

  fprintf ( 1, '\n' );
  fprintf ( 1, '   i      U(i)       D18(U)(i)  D18inv(D18(U))(i)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %2d  %10.4f  %10.4f  %10.4f\n', i, u(i), v(i), w(i) );
  end
%
%  Constant signal.
%
  n = 8;
  u(1:n) = 1.0;

  v = daub18_transform ( n, u );

  w = daub18_transform_inverse ( n, v );

  fprintf ( 1, '\n' );
  fprintf ( 1, '   i      U(i)        D18(U)(i)  D18inv(D18(U))(i)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %2d  %10.4f  %10.4f  %10.4f\n', i, u(i), v(i), w(i) );
  end
%
%  Linear signal.
%
  n = 16;
  a_first = 1.0;
  a_last = n;
  u = linspace ( a_first, a_last, n );

  v = daub18_transform ( n, u );

  w = daub18_transform_inverse ( n, v );

  fprintf ( 1, '\n' );
  fprintf ( 1, '   i      U(i)        D18(U)(i)  D18inv(D18(U))(i)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %2d  %10.4f  %10.4f  %10.4f\n', i, u(i), v(i), w(i) );
  end
%
%  Quadratic data.
%
  n = 8;
  u(1) = 25.0;
  u(2) = 16.0;
  u(3) = 9.0;
  u(4) = 4.0;
  u(5) = 1.0;
  u(6) = 0.0;
  u(7) = 1.0;
  u(8) = 4.0;

  v = daub18_transform ( n, u );

  w = daub18_transform_inverse ( n, v );

  fprintf ( 1, '\n' );
  fprintf ( 1, '   i      U(i)        D18(U)(i)  D18inv(D18(U))(i)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %2d  %10.4f  %10.4f  %10.4f\n', i, u(i), v(i), w(i) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function odf = FourierODF(C,varargin)
% defines an ODF by its Fourier coefficients
%
% Syntax
%   odf = FourierODF(C,CS,SS)
%
% Input
%  C      - Fourier coefficients / C coefficients
%  CS, SS - crystal, specimen @symmetry
%
% Output
%  odf - @SO3Fun
%
% See also
% uniformODF unimodalODF BinghamODF fibreODF

odf = SO3FunHarmonic(C,varargin{:});

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function varargout = fillMeshFaces(varargin)
%FILLMESHFACES Fill the faces of a mesh with the specified colors.
%
%   fillMeshFaces(V, F, VERTEXCOLORS)
%   Colorizes a mesh by filling faces with an array of values. The colors
%   can be a NV-by-1 array of values, or a NV-by-3 array of values.
%   Face filling uses 'interp' coloring mode.
%
%   fillMeshFaces(V, F, FACECOLORS)
%   Colorizes the mesh by specifying the value or the color associated to
%   each face. Face filling uses 'flat' coloring mode.
%
%   fillMeshFaces(..., PNAME, PVALUE)
%   Specifies additional parameters that will be passed to the 'patch'
%   function.
%
%   Example
%     % Colorize mesh based on z-coordinate of vertices.
%     [v, f] = createIcosahedron;
%     values = v(:,3);
%     figure; axis equal; view(3);
%     fillMeshFaces(v, f, values);
%
%     % Colorize mesh using specific color for each face
%     [v, f] = createIcosahedron;
%     colors = jet(20);
%     figure; axis equal; view(3);
%     fillMeshFaces(v, f, colors);
%
%   See also
%     drawMesh

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2020-04-16, using Matlab 9.7.0.1247435 (R2019b) Update 2
% Copyright 2020-2022 INRAE - BIA Research Unit - BIBS Platform (Nantes)

%% Parse input arguments

% extract first argument
var1 = varargin{1};
varargin(1) = [];

% Check if first input argument is an axes handle
if isAxisHandle(var1)
    ax = var1;
    var1 = varargin{1};
    varargin(1) = [];
else
    ax = gca;
end

% Check if the input is a mesh structure
if isstruct(var1)
    % extract data to display
    vertices = var1.vertices;
    faces = var1.faces;
else
    % assumes input is given with vertices+faces arrays
    vertices = var1;
    faces = varargin{1};
    varargin(1) = [];
end

% next argument is face color
colors = varargin{1};
varargin(1) = [];

% adapt the face color key value depending on the size of the "color" input
% argument
faceColorMode = 'interp';
if size(colors, 1) == size(faces, 1)
    faceColorMode = 'flat';
end

% array FACES is a NF-by-NV indices array, with NV number of vertices of
% each face, and NF number of faces
h = patch('Parent', ax, ...
    'vertices', vertices, 'faces', faces, 'FaceVertexCData', colors, ...
    'FaceColor', faceColorMode, varargin{:});

%% Process output arguments

% format output parameters
if nargout > 0
    varargout = {h};
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
% HARRIS - Harris corner detector
%
% Usage:                 cim = harris(im, sigma)
%                [cim, r, c] = harris(im, sigma, thresh, radius, disp)
%  [cim, r, c, rsubp, csubp] = harris(im, sigma, thresh, radius, disp)
%
% Arguments:
%            im     - image to be processed.
%            sigma  - standard deviation of smoothing Gaussian. Typical
%                     values to use might be 1-3.
%            thresh - threshold (optional). Try a value ~1000.
%            radius - radius of region considered in non-maximal
%                     suppression (optional). Typical values to use might
%                     be 1-3.
%            disp   - optional flag (0 or 1) indicating whether you want
%                     to display corners overlayed on the original
%                     image. This can be useful for parameter tuning. This
%                     defaults to 0
%
% Returns:
%            cim    - binary image marking corners.
%            r      - row coordinates of corner points.
%            c      - column coordinates of corner points.
%            rsubp  - If five return values are requested sub-pixel
%            csubp  - localization of feature points is attempted and
%                     returned as an additional set of floating point
%                     coords. Note that you may still want to use the integer
%                     valued coords to specify centres of correlation windows
%                     for feature matching.
%
% If thresh and radius are omitted from the argument list only 'cim' is returned
% as a raw corner strength image.  You may then want to look at the values
% within 'cim' to determine the appropriate threshold value to use. Note that
% the Harris corner strength varies with the intensity gradient raised to the
% 4th power.  Small changes in input image contrast result in huge changes in
% the appropriate threshold.
%
% Note that this code computes Noble's version of the detector which does not
% require the parameter 'k'.  See comments in code if you wish to use Harris'
% original measure.
%
% See also: NONMAXSUPPTS, DERIVATIVE5

% References:
% C.G. Harris and M.J. Stephens. "A combined corner and edge detector",
% Proceedings Fourth Alvey Vision Conference, Manchester.
% pp 147-151, 1988.
%
% Alison Noble, "Descriptions of Image Surfaces", PhD thesis, Department
% of Engineering Science, Oxford University 1989, p45.

% Copyright (c) 2002-2010 Peter Kovesi
% Centre for Exploration Targeting
% The University of Western Australia
% http://www.csse.uwa.edu.au/~pk/research/matlabfns/
%
% 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, 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.

% March    2002 - Original version
% December 2002 - Updated comments
% August   2005 - Changed so that code calls nonmaxsuppts
% August   2010 - Changed to use Farid and Simoncelli's derivative filters

function [cim, r, c, rsubp, csubp] = harris(im, sigma, thresh, radius, disp)

    error(nargchk(2,5,nargin));
    if nargin == 4
	disp = 0;
    end

    if ~isa(im,'double')
	im = double(im);
    end

    subpixel = nargout == 5;

    % Compute derivatives and elements of the structure tensor.
    [Ix, Iy] = derivative5(im, 'x', 'y');
    Ix2 = gaussfilt(Ix.^2,  sigma);
    Iy2 = gaussfilt(Iy.^2,  sigma);
    Ixy = gaussfilt(Ix.*Iy, sigma);

    % Compute the Harris corner measure. Note that there are two measures
    % that can be calculated.  I prefer the first one below as given by
    % Nobel in her thesis (reference above).  The second one (commented out)
    % requires setting a parameter, it is commonly suggested that k=0.04 - I
    % find this a bit arbitrary and unsatisfactory.

    cim = (Ix2.*Iy2 - Ixy.^2)./(Ix2 + Iy2 + eps); % My preferred  measure.
%    k = 0.04;
%    cim = (Ix2.*Iy2 - Ixy.^2) - k*(Ix2 + Iy2).^2; % Original Harris measure.

    if nargin > 2   % We should perform nonmaximal suppression and threshold

	if disp  % Call nonmaxsuppts to so that image is displayed
	    if subpixel
		[r,c,rsubp,csubp] = nonmaxsuppts(cim, radius, thresh, im);
	    else
		[r,c] = nonmaxsuppts(cim, radius, thresh, im);
	    end
	else     % Just do the nonmaximal suppression
	    if subpixel
		[r,c,rsubp,csubp] = nonmaxsuppts(cim, radius, thresh);
	    else
		[r,c] = nonmaxsuppts(cim, radius, thresh);
	    end
	end
    end

author: jianxiongxiao

repo: ProfXkit


text:
function [A, AH, normA] = MakeTransforms(type, N, params)

% [A, AH] = MakeTransforms(type, N, params)
% with A AH = I
%
% INPUT
%    N : signal length
%    type : type of transform can be 'DFT' or 'STFT'
%    params : parameters for transform
%
%    For DFT, params = [Nfft]
%       Nfft = number of DFT coefficients, Nfft >= N
%       Nfft/N is over-sampling factor
%
%   For STFT, params = [R, Nfft]
%       R = length of frame
%       Nfft = number of DFT coefficients, Nfft >= R
%
%    For DCT, params = [Ndct]
%       Ndct = number of DCT coefficients, Ndct >= N
%       Ndct/N is over-sampling factor
%
%   For STDCT, params = [R, Ndct]
%       R = length of frame
%       Ndct = number of DCT coefficients, Ndct >= R
%
% OUTPUT
%    A, AH : function handles for transform
%       AH is the conjugate transpose of A
%       Note: input to AH must be a row vector of length N

% DFT : Discrete Fourier transform with zero-padding
% STFT : Short-time Fourier transform
% DCT : Discrete cosine transform with zero-padding
% STDCT : Short-time DCT

switch type
    case '2D-CDWT'
        [Faf, Fsf] = FSfarras;
        [af, sf] = dualfilt1;
        AH = @(x) cplxdual2D(x,params(1),Faf,af);
        A = @(x) icplxdual2D(x,params(1),Fsf,sf);

        normA = 1;

    case 'DWT'
        [h0, h1, g0, g1] = daubf(params(2));    % Get filters
        AH = @(x) udwt(x, params(1), h0, h1);
        A = @(x) iudwt(x, params(1), g0, g1);
        normA = 1;
    case '2D-DWT'
        [Faf, Fsf] = FSfarras;
        [af, sf] = dualfilt1;
        AH = @(x) dualtree2D(x,params(1), Faf, af);
        A = @(x) idualtree2D(x,params(1), Fsf, sf);

        normA = 1;
    case 'DFT'

        % over-sampled DFT (oversampling factor = Nfft/N)
        Nfft = params(1);
        truncate = @(x, N) x(1:N);
        AH = @(x) fft(x, Nfft)/sqrt(Nfft);
        A = @(X) truncate(ifft(X), N)*sqrt(Nfft);

        normA = sqrt(N/Nfft);

    case 'STFT'

        R = params(1);
        M = params(2);
        K = params(3);
        Nfft = params(4);
        AH = @(x) pSTFT2(x, R, M, K, Nfft);
        A = @(x) ipSTFT2(x,R,M,K,N);

        normA = sqrt(R/(M*Nfft));

    case 'DCT'

        % over-sampled DFT (oversampling factor = Nfft/N)
        Ndct = params(1);
        truncate = @(x, N) x(1:N);
        AH = @(x) dct(x, Ndct);
        A = @(X) truncate(idct(X), N);

        normA = sqrt(N/Ndct);

    case 'STDCT'

        R = params(1);
        M = params(2);
        K = params(3);
        Ndct = params(4);
        AH = @(x) stdct(x, R,M,K, Ndct);
        A = @(x) istdct(x, R,M,K, N);

        normA = sqrt(R/(2*Ndct));

end

% x = randn(1, N);         % assuming row vector ! ! !
%
% % Verify perfect reconstruction property of transform 1 (A*AH = I)
% c = AH(x);
% err = x - A(c)';
% re = max(abs(err(:)));          % should be zero
% fprintf('Reconstruction error = %f\n', re)
%
% % Verify Parseval energy identity
% E = sum(abs(x(:)).^2);
% E1 = sum(abs(c(:)).^2);         % should be unity
% fprintf('Energy ratio = %f\n', E/E1)

author: thomaskuestner


text:
function err ( )

%*****************************************************************************80
%
%% ERR computes the L2 and H1 errors using a 2 point Gauss quadrature rule.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 April 2006
%

%
%  Set up the quadrature information.
%
  [ xq, wq ] = gauss3pt ( );

  el2 = 0.0;
  eh1 = 0.0;

  for it = 1 : nel

    for iq = 1 : nq
      ar = area(it) * wq(iq);
      x = xc(node(it,1)) + area(it) / 2.0 * ( xq(iq) + 1.0 );
      uh = 0.0;
      uhx = 0.0;
      for in = 1 : nnodes
        ip = node(it,in);
        [ bb, bx ] = quadbf ( x, it, in, xc, node );
        i = indx(ip);
        if ( i ~= 0 )
          uh = uh + bb * f(i);
          uhx = uhx + bx * f(i);
        end
      end
      [ uex, uexx ] = exact ( x );
      el2 = el2 + ar * ( uh - uex )^2;
      eh1 = eh1 + ar * ( uhx - uexx )^2 );
    end
  end

  el2 = sqrt ( el2 )
  eh1 = sqrt ( eh1 )

  return
end


text:
function [ts] = PL2TsBySource( filename, source, channelInSource, unit )
% PL2TsBySource( filename, source, channelInSource, unit ): read spike timestamps from a .pl2 file
%
% ts = PL2TsBySource( filename, source, channelInSource, unit );
% ts = PL2TsBySource( filename, 'SPK', 1, 2 );
% ts = PL2TsBySource( filename, 7, 1, 3 );
%
% INPUT:
%   filename - if empty string, will use File Open dialog
%   source - source name or source numeric id
%   channel -  1-based channel number within source
%           To print the list of spike channels, use these commands:
%
%           pl2 = PL2GetFileIndex(filename);
%           PL2Print(pl2.SpikeChannels);
%
%   unit  - unit number (0 - unsorted, 1-255 sorted units)
%
% OUTPUT:
%   ts - array of timestamps (in seconds)

n = 0;
ts = [];

if nargin ~= 4
    error 'expected 4 input arguments';
end

filename = internalPL2ResolveFilename(filename);
pl2 = PL2GetFileIndex(filename);
channelNumber = internalPL2ResolveChannelBySource(pl2.SpikeChannels, source, channelInSource);
ts = Pl2Ts(filename, channelNumber, unit);

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% Superclass for linear filters.
%
% This file is part of the Kernel Adaptive Filtering Toolbox for Matlab.
% https://github.com/steven2358/kafbox/

classdef linear_filter < base_estimator

end

repo: kafbox


text:
function [Ydt, X, R] = detrend_data(Y,nk, method)
%% detrend fluorescence signals with B-spline basis
%% Inputs:
%   Y: d X T matrix, video data
%   nk: scalar, number of knots
%   method: {'local_min', 'spline'} method for removing the trend
% Outputs:
%   Ydt: d X T matrix, detrended data
%   X: T*M matrix, each column is one basis
%   R: d*M matrix, coefficients of all basis for all pixels

%% create basis
[~, T] = size(Y);

if ~exist('nk', 'var')
    nk = 5;
end
if ~exist('method', 'var') || isempty(method)
    method = 'spline';
end

if strcmpi(method, 'spline')
    X = bsplineM((1:T)', linspace(1, T, nk), 4);

    %% compute coefficients of all spline basis
    R = (Y*X)/(X'*X);

    %% compute detrended data
    Ydt = Y-R*X';
else
    k = ceil(T/nk);
    [d, T] = size(Y);
    Tnew = ceil(T/k)*k;
    if T~=Tnew
        Y(:, (T+1):Tnew) = repmat(Y(:, T), [1, Tnew-T]);
    end
    Y = reshape(Y, d, k, []);
    Ydt = reshape(bsxfun(@minus, Y, min(Y,[], 2)), d, []);
    Ydt = Ydt(:, 1:T);
    X = [];
    R = [];
end

author: zhoupc


text:
%% SeDuMi Interface Guide

% See the following Wiki page:
% http://www.i2c2.aut.ac.nz/Wiki/OPTI/index.php/Solvers/SEDUMI

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function result = compare(c1,c2, tolerance)

%tstoolbox/@core/compare
%   Syntax:
%     * compare(c1,c2, tolerance)
%
%   Input Arguments:
%     * c1,c2 core object of two signals
%     * tolerance tolerance of the signals's RMS value (default
%       tolerance=1e-6)
%
%   compare compare two signals whether they have equal values slight
%   differences due to rounding errors are ignored depending on the value
%   of tolerance when signals are found to be not equal, a zero is
%   returned.
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

if nargin < 2
	help(mfilename)
	return
end

if nargin < 3
	tolerance = 1e-6;
end

if dlens(c1) ~= dlens(c2)
	result = 0;
else
	ref = data(rms(c1));
	diff = data(rms(c1-c2));
	if any((diff ./ ref) > tolerance)
		result = 0;
	else
		result = 1;
	end
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function a = daub10 ( n )

%*****************************************************************************80
%
%% DAUB10 returns the DAUB10 matrix.
%
%  Discussion:
%
%    The DAUB10 matrix is the Daubechies wavelet transformation matrix
%    with 10 coefficients.
%
%    Note that in the reference, the coefficient 0.0775714938400459
%    is given incorrectly, with the "8" misrepresented as a "0".
%
%  Properties:
%
%    The family of matrices is nested as a function of N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 July 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Gilbert Strang, Truong Nguyen,
%    Wavelets and Filter Banks,
%    Wellesley-Cambridge Press, 1997,
%    ISBN: 0-9614088-7-1,
%    LC: TK7872.F5S79 / QA403.3.S87
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be at least 10, and a multiple of 2.
%
%    Output, real A(N,N), the matrix.
%
  c = [ ...
    0.1601023979741929, ...
    0.6038292697971895, ...
    0.7243085284377726, ...
    0.1384281459013203, ...
   -0.2422948870663823, ...
   -0.0322448695846381, ...
    0.0775714938400459, ...
   -0.0062414902127983, ...
   -0.0125807519990820, ...
    0.0033357252854738 ];

  if ( n < 10 || mod ( n, 2 ) ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'DAUB10 - Fatal error!\n' );
    fprintf ( 1, '  N must be at least 10 and a multiple of 2.\n' );
    error ( 'DAUB10 - Fatal error!' );
  end

  a = zeros ( n, n );

  for i = 1 : 2 : n - 1

    a(i,i)                  =    c(1);
    a(i,i+1)                =    c(2);
    a(i,i4_wrap(i+2,1,n))   =    c(3);
    a(i,i4_wrap(i+3,1,n))   =    c(4);
    a(i,i4_wrap(i+4,1,n))   =    c(5);
    a(i,i4_wrap(i+5,1,n))   =    c(6);
    a(i,i4_wrap(i+6,1,n))   =    c(7);
    a(i,i4_wrap(i+7,1,n))   =    c(8);
    a(i,i4_wrap(i+8,1,n))   =    c(9);
    a(i,i4_wrap(i+9,1,n))   =   c(10);

    a(i+1,i)                =   c(10);
    a(i+1,i+1)              =  - c(9);
    a(i+1,i4_wrap(i+2,1,n)) =    c(8);
    a(i+1,i4_wrap(i+3,1,n)) =  - c(7);
    a(i+1,i4_wrap(i+4,1,n)) =    c(6);
    a(i+1,i4_wrap(i+5,1,n)) =  - c(5);
    a(i+1,i4_wrap(i+6,1,n)) =    c(4);
    a(i+1,i4_wrap(i+7,1,n)) =  - c(3);
    a(i+1,i4_wrap(i+8,1,n)) =    c(2);
    a(i+1,i4_wrap(i+9,1,n)) =  - c(1);

  end

  return
end


text:
% [INPUT]
% p = A vector of floats [0,Inf) of length t representing the prices.
% r = A vector of floats (-Inf,Inf) of length t representing the logarithmic returns.
% v = A vector of floats [0,Inf) of length t representing the trading volumes.
% cp = A vector of floats [0,Inf) of length t representing the market capitalization.
% bwl = An integer [90,252] representing the dimension of the long bandwidth.
% bwm = An integer [21,90) representing the dimension of the medium bandwidth.
% bws = An integer [5,21) representing the dimension of the short bandwidth.
%
% [OUTPUT]
% hhlr = A column vector of floats [0,Inf) of length t representing the Hui-Heubel Liquidity Ratio.
% tr = A column vector of floats [0,Inf) of length t representing the Turnover Ratio.
% vr = A column vector of floats [0,Inf) of length t representing the Variance Ratio.

function [hhlr,tr,vr] = liquidity_metrics(varargin)

    persistent ip;

    if (isempty(ip))
        ip = inputParser();
        ip.addRequired('p',@(x)validateattributes(x,{'double'},{'real' 'finite' 'nonnegative' 'vector' 'nonempty'}));
        ip.addRequired('r',@(x)validateattributes(x,{'double'},{'real' 'finite' 'vector' 'nonempty'}));
        ip.addRequired('v',@(x)validateattributes(x,{'double'},{'real' 'finite' 'nonnegative' 'vector' 'nonempty'}));
        ip.addRequired('cp',@(x)validateattributes(x,{'double'},{'real' 'finite' 'nonnegative' 'vector' 'nonempty'}));
        ip.addRequired('bwl',@(x)validateattributes(x,{'double'},{'real' 'finite' 'integer' '>=' 90 '<=' 252 'scalar'}));
        ip.addRequired('bwm',@(x)validateattributes(x,{'double'},{'real' 'finite' 'integer' '>=' 21 '<' 90 'scalar'}));
        ip.addRequired('bws',@(x)validateattributes(x,{'double'},{'real' 'finite' 'integer' '>=' 5 '<' 21 'scalar'}));
    end

    ip.parse(varargin{:});

    ipr = ip.Results;
    [p,r,v,cp,bwl,bwm,bws] = validate_input(ipr.p,ipr.r,ipr.v,ipr.cp,ipr.bwl,ipr.bwm,ipr.bws);

    nargoutchk(3,3);

    [hhlr,tr,vr] = liquidity_metrics_internal(p,r,v,cp,bwl,bwm,bws);

end

function [hhlr,tr,vr] = liquidity_metrics_internal(p,r,v,cp,bwl,bwm,bws)

    hhlr = calculate_hhlr(p,v,cp,bwl,bws);
    tr = calculate_tr(v,cp,bwl);
    vr = calculate_vr(r,bwl,bwm);

end

function hhlr = calculate_hhlr(p,v,cp,bwl,bws)

    tr = v ./ cp;
    tr(~isfinite(tr)) = 0;

    windows_p = extract_rolling_windows(p,bws);
    dp = cellfun(@(x)(max(x) - min(x)) / min(x),windows_p);

    alpha = 2 / (bwl + 1);

    hhlr = dp ./ tr;
    hhlr(~isfinite(hhlr)) = 0;
    hhlr(1:bws) = mean(hhlr(bws+1:bws*2+1));
    hhlr = [hhlr(1); filter(alpha,[1 (alpha - 1)],hhlr(2:end),(1 - alpha) * hhlr(1))];
    hhlr = (hhlr - min(hhlr)) ./ (max(hhlr) - min(hhlr));

end

function tr = calculate_tr(v,cp,bwl)

    alpha = 2 / (bwl + 1);

    tr = v ./ cp;
    tr(~isfinite(tr)) = 0;
    tr = [tr(1); filter(alpha,[1 (alpha - 1)],tr(2:end),(1 - alpha) * tr(1))];
    tr = (tr - min(tr)) ./ (max(tr) - min(tr));

end

function vr = calculate_vr(r,bwl,bwm)

    alpha = 2 / (bwl + 1);
    t = bwl / bwm;

    windows_long = extract_rolling_windows(r,bwl);
    var_long = cellfun(@var,windows_long);

    windows_short = extract_rolling_windows(r,bwm);
    var_short = cellfun(@var,windows_short);

    vr = var_long ./ (t .* var_short);
    vr(~isfinite(vr)) = 0;
    vr(1:bwm) = mean(vr(bwm+1:bwm*2+1));
    vr = [vr(1); filter(alpha,[1 (alpha - 1)],vr(2:end),(1 - alpha) * vr(1))];

end

function [p,r,v,cp,bwl,bwm,bws] = validate_input(p,r,v,cp,bwl,bwm,bws)

    data = {p(:) r(:) v(:) cp(:)};

    l = unique(cellfun(@numel,data));

    if (numel(l) ~= 1)
        error('The number of elements of ''p'', ''r'' and ''v'' must be equal.');
    end

    if (l < 5)
        error('The value of ''p'', ''r'' and ''v'' is invalid. Expected inputs to be vectors containing at least 5 elements.');
    end

    [p,r,v,cp] = deal(data{:});

    if (bwl < (bwm * 2))
        error(['The long bandwidth (' num2str(bwl) ') must be at least twice the medium bandwidth (' num2str(bwm) ').']);
    end

    if (bwm < (bws * 2))
        error(['The medium bandwidth (' num2str(bwm) ') must be at least twice the short bandwidth (' num2str(bws) ').']);
    end

end

author: TommasoBelluzzo

repo: SystemicRisk


text:
function linplus_test23 ( )

%*****************************************************************************80
%
%% TEST23 tests R8GB_FA, R8GB_TRF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  m = 5;
  n = m;
  ml = 1;
  mu = 1;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST23\n' );
  fprintf ( 1, '  For a general banded matrix,\n' );
  fprintf ( 1, '  R8GB_FA factors, using LINPACK conventions;\n' );
  fprintf ( 1, '  R8GB_TRF factors, using LAPACK conventions;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix rows M              = %d\n', m );
  fprintf ( 1, '  Matrix columns N           = %d\n', n );
  fprintf ( 1, '  Lower bandwidth ML         = %d\n', ml );
  fprintf ( 1, '  Upper bandwidth MU         = %d\n', mu );
%
%  Set the matrix.
%
  [ a, seed ] = r8gb_random ( m, n, ml, mu, seed );
%
%  Factor the matrix.
%
  [ a_lu, pivot, info ] = r8gb_fa ( n, ml, mu, a );

  r8gb_print ( m, n, ml, mu, a_lu, '  The R8GB_FA factors:' );
%
%  Set the matrix.
%
  seed = 123456789;
  [ a, seed ] = r8gb_random ( m, n, ml, mu, seed );
%
%  Factor the matrix.
%
  [ a_lu, pivot, info ] = r8gb_trf ( m, n, ml, mu, a );

  r8gb_print ( m, n, ml, mu, a_lu, '  The R8GB_TRF factors:');

  return
end


text:
function COV = UnTangent_space(T,C)
NTrial = size(T,2);
N_elec = (sqrt(1+8*size(T,1))-1)/2;
COV = zeros(N_elec,N_elec,NTrial);

if nargin<2
    C = riemann_mean(COV);
end

index = reshape(triu(ones(N_elec)),N_elec*N_elec,1)==0;

Out = zeros(N_elec*N_elec,NTrial);

Out(not(index),:) = T;
P = C^0.5;
for i=1:NTrial
  tmp = reshape(Out(:,i),N_elec,N_elec,[]);
  tmp =  diag(diag(tmp))+triu(tmp,1)/sqrt(2) + triu(tmp,1)'/sqrt(2);
  tmp = P*tmp*P;
  COV(:,:,i) = RiemannExpMap(C,tmp);
end

author: alexandrebarachant

repo: covariancetoolbox


text:
%................................................................

function stiffness=formStiffness2D(GDof,numberElements,...
    elementNodes,numberNodes,nodeCoordinates,D,thickness)

% compute stiffness matrix
% for plane stress Q4 elements

stiffness=zeros(GDof);

% 2 by 2 quadrature
[gaussWeights,gaussLocations]=gauss2d('2x2');

for e=1:numberElements
  numNodePerElement = length(elementNodes(e,:));
  numEDOF = 2*numNodePerElement;
  elementDof=zeros(1,numEDOF);
  for i = 1:numNodePerElement
      elementDof(2*i-1)=2*elementNodes(e,i)-1;
      elementDof(2*i)=2*elementNodes(e,i);
  end

  % cycle for Gauss point
  for q=1:size(gaussWeights,1)
    GaussPoint=gaussLocations(q,:);
    xi=GaussPoint(1);
    eta=GaussPoint(2);

% shape functions and derivatives
    [shapeFunction,naturalDerivatives]=shapeFunctionQ4(xi,eta);

% Jacobian matrix, inverse of Jacobian,
% derivatives w.r.t. x,y
    [Jacob,invJacobian,XYderivatives]=...
        Jacobian(nodeCoordinates(elementNodes(e,:),:),naturalDerivatives);

%  B matrix
    B=zeros(3,numEDOF);
    B(1,1:2:numEDOF)       = XYderivatives(:,1)';
    B(2,2:2:numEDOF)  = XYderivatives(:,2)';
    B(3,1:2:numEDOF)       = XYderivatives(:,2)';
    B(3,2:2:numEDOF)  = XYderivatives(:,1)';

% stiffness matrix
    stiffness(elementDof,elementDof)=...
        stiffness(elementDof,elementDof)+...
        B'*D*thickness*B*gaussWeights(q)*det(Jacob);
  end
end

repo: Aero-matlab


text:
function err = RotateMatrixCriterion(R,U)

dx = R*U(:,1)-[1;0;0;0];
dy = R*U(:,2)-[0;1;0;0];
err = dx'*dx + dy'*dy;

author: kristinbranson

repo: JAABA


text:
function x = emailFeatures(word_indices)
%EMAILFEATURES takes in a word_indices vector and produces a feature vector
%from the word indices
%   x = EMAILFEATURES(word_indices) takes in a word_indices vector and
%   produces a feature vector from the word indices.

% Total number of words in the dictionary
n = 1899;

% You need to return the following variables correctly.
x = zeros(n, 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return a feature vector for the
%               given email (word_indices). To help make it easier to
%               process the emails, we have have already pre-processed each
%               email and converted each word in the email into an index in
%               a fixed dictionary (of 1899 words). The variable
%               word_indices contains the list of indices of the words
%               which occur in one email.
%
%               Concretely, if an email has the text:
%
%                  The quick brown fox jumped over the lazy dog.
%
%               Then, the word_indices vector for this text might look
%               like:
%
%                   60  100   33   44   10     53  60  58   5
%
%               where, we have mapped each word onto a number, for example:
%
%                   the   -- 60
%                   quick -- 100
%                   ...
%
%              (note: the above numbers are just an example and are not the
%               actual mappings).
%
%              Your task is take one such word_indices vector and construct
%              a binary feature vector that indicates whether a particular
%              word occurs in the email. That is, x(i) = 1 when word i
%              is present in the email. Concretely, if the word 'the' (say,
%              index 60) appears in the email, then x(60) = 1. The feature
%              vector should look like:
%
%              x = [ 0 0 0 0 1 0 0 0 ... 0 0 0 0 1 ... 0 0 0 1 0 ..];
%
%

for i = 1:length(word_indices)
    x(word_indices(i)) = 1
end

% =========================================================================

end

author: Ayatans

repo: Machine-Learning-homework


text:
function PlotTempData(city, xData, localAnomData, fitval, latStr, lngStr)
%#event
xlabelData = xData{1};
xticks = xData{2};
xtickValues = xticks{1};
xtickLabels = xticks{2};

% Draw a pretty picture
figure;
hold on
plot(xlabelData, localAnomData, 'b.', 'MarkerSize', 4);
line(xlabelData, fitval, 'Color', 'red', 'LineWidth', 2);
hold off
set(gca, 'XLim', [xlabelData(1), xlabelData(end)]);
set(gca, 'XTick', xtickValues, 'XTickLabel', xtickLabels);

% Annotate the plot
titleStr = sprintf(['Monthly Local Temperature Anomaly\nLatitude %s   ' ...
                    'Longitude %s\nGrid Center: %s'], latStr, lngStr, ...
                    city);
title(titleStr, 'FontWeight', 'Bold');

xlabel('Date', 'FontWeight', 'Bold');
ylabel('Variance Degrees Centigrade', 'FontWeight', 'Bold');

author: Sable

repo: mcbench-benchmarks


text:
function bool = index_selected(iter, indeces)

if islogical(indeces)
  if isscalar(indeces)
    bool = indeces;
  else
    error('Logical values must not be matrices');
  end
elseif isscalar(indeces)
  if indeces > 0
    bool = iter >= indeces;
  else
    error('Index must be positive integer.');
  end
else
  bool = any(iter==indeces);
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [mri] = align_fsaverage2mni(mri)

% ALIGN_FSAVERAGE2MNI performs an affine alignment of the anatomical volume from
% FSAVERAGE towards MNI coordinates. Only the homogeneous transformation matrix is
% modified and the coordsys-field is updated.
%
% Use as
%   mri = align_fsaverage2mni(mri)
% where the first input argument is a FieldTrip MRI-structure.
%
% with fsaverage we mean MNI305
% with mni       we mean MNI152, i.e. the template used in SPM
%
% See http://freesurfer.net/fswiki/CoordinateSystems
%
% See also ALIGN_CTF2ACPC, ALIGN_NEUROMAG2ACPC

fsaverage2mni = [
   0.9975   -0.0073    0.0176   -0.0429
   0.0146    1.0009   -0.0024    1.5496
  -0.0130   -0.0093    0.9971    1.1840
  ];

assert(strcmp(mri.coordsys, 'fsaverage'), 'incorrect input coordinate system ''%s''', mri.coordsys);
mri.transform = fsaverage2mni * mri.transform;
mri.coordsys = 'mni';

author: spm


text:
function [dt6, xformToAcpc, mmPerVoxel, fileName, desc, intentName] = dtiLoadTensorsFromNifti(ni)
%Load a dt6 (tensor) data file in mrDiffusion format from NIFTI format.
%
%  [dt6, xformToAcpc, mmPerVox, fileName, desc, intentName] = dtiLoadTensorsFromNifti(niFile)
%
% The six entries are the diffusion tensor values, derived from the raw
% data.  The raw data can be in many different directions.  The diffusion
% tensor is a 3x3 positive-definite matrix, D.  The entries in the matrix
% are stored in a vector:
%
%    Dxx Dyy Dzz Dxy Dxz Dyz
%
% or equivalently
%
%    D(1,1), D(2,2), D(3,3), D(1,2), D(1,3), D(2,3)
%
% HISTORY:
%  2007.10.04 AJS: Wrote it.
%
% (c) Stanford VISTA Team

if(~exist('ni','var')||isempty(ni))
    [f,p] = uigetfile({'*.nii.gz';'*.*'},'Select a NIFTI tensor file...');
    if(isnumeric(f)), disp('User canceled.'); return; end
    ni = fullfile(p,f);
end

if(ischar(ni))
    if exist(ni,'file'), ni = niftiRead(ni);
    else error('Can not find file %s\n',ni);
    end
end

fileName = ni.fname;

% We convert from the 5d, lower-tri row order NIFTI tensor format used by
% other groups, such as FSL, which stores (x,y,z,1,directions)
%    direction ordering:  (Dxx Dxy Dyy Dxz Dyz Dzz)
%
% to our 4d tensor format, (x,y,z,directions)
%    direction ordering: (Dxx Dyy Dzz Dxy Dxz Dyz).
%
dt6 = double(squeeze(ni.data(:,:,:,1,[1 3 6 2 4 5])));
xformToAcpc = ni.qto_xyz;
mmPerVoxel = ni.pixdim(1:3);
desc = ni.descrip;
intentName = ni.intent_name;

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%% housekeeping
close all
clear
clc

%% It just amounts to adding a keyword in the parameter list

he=rise('hybrid_expect');

%%
he=solve(he);

author: jmaih


text:
function [ raw_label ] = frame2rawlabel(label, win_len, win_step)

num_frame = length(label);

total_len = (num_frame-1)*win_step + win_len;

raw_label = zeros(1, total_len);
% temp_label = zeros(1,win_len);
start_indx = 0;
i = 1;
while(1)
    if(start_indx+win_len>total_len)
        break;
    end
    if i==1
        raw_label(start_indx+1 : start_indx + win_len) = label(i);
    else
        temp_label(:) = label(i);
        raw_label(start_indx+1 : start_indx + win_len) = raw_label(start_indx+1 : start_indx + win_len) + temp_label;
%         temp_label = zeros(1,total_len);
    end
    i = i + 1;
    start_indx = start_indx + win_step;
end

raw_label = raw_label>= 1;

end

author: jtkim-kaist

repo: VAD

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


text:
function n = navier_q2(xy,mv,flowsol)
%navier_q2  Q2 convection matrix
%   N = navier_q2(xy,ev,flowsol);
%   input
%          xy         Q2 nodal coordinate vector
%          ev         element mapping matrix
%          flowsol    Q2-Q1 or Q2-P1 flow solution
%   output
%          N          Q2 scalar convection matrix
%
%   Natural boundary conditions apply.
%   IFISS function: DJS; 11 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
nngpt=9;
x=xy(:,1); y=xy(:,2);
nvtx=length(x);   nel=length(mv(:,1));
usol=flowsol(1:nvtx); vsol=flowsol(nvtx+1:2*nvtx);
fprintf('setting up Q2 convection matrix...  ')
%
% initialise global matrices
      n = sparse(nvtx,nvtx);
%
%
% set up 3x3 Gauss points
      gpt=sqrt(0.6);
      s(1) = -gpt; t(1) = -gpt; wt(1)=25/81;
      s(2) =  gpt; t(2) = -gpt; wt(2)=25/81;
      s(3) =  gpt; t(3) =  gpt; wt(3)=25/81;
      s(4) = -gpt; t(4) =  gpt; wt(4)=25/81;
      s(5) =  0.0; t(5) = -gpt; wt(5)=40/81;
      s(6) =  gpt; t(6) =  0.0; wt(6)=40/81;
      s(7) =  0.0; t(7) =  gpt; wt(7)=40/81;
      s(8) = -gpt; t(8) =  0.0; wt(8)=40/81;
      s(9) =  0.0; t(9) =  0.0; wt(9)=64/81;
%
% inner loop over elements
        for ivtx = 1:4
        xl_v(:,ivtx) = x(mv(:,ivtx));
        yl_v(:,ivtx) = y(mv(:,ivtx));
        end
        for idx = 1:9
		xsl(:,idx) = usol(mv(:,idx));
		ysl(:,idx) = vsol(mv(:,idx));
		end
      ne = zeros(nel,9,9);
%
% loop over Gauss points
         for igpt = 1:nngpt
         sigpt=s(igpt);
         tigpt=t(igpt);
         wght=wt(igpt);
%  evaluate derivatives etc
         [jac,invjac,phi,dphidx,dphidy] = deriv(sigpt,tigpt,xl_v,yl_v);
         [psi,dpsidx,dpsidy] = qderiv(sigpt,tigpt,xl_v,yl_v);
         u_x = zeros(nel,1); u_y=zeros(nel,1);
            for k=1:9
		    u_x(:) = u_x(:) + xsl(:,k) .* psi(:,k);
		    u_y(:) = u_y(:) + ysl(:,k) .* psi(:,k);
		    end
		for j = 1:9
            for i = 1:9
				ne(:,i,j)  = ne(:,i,j)  + wght*u_x(:).*psi(:,i).*dpsidx(:,j);
                ne(:,i,j)  = ne(:,i,j)  + wght*u_y(:).*psi(:,i).*dpsidy(:,j);
             end
	    end
%
% end of Gauss point loop
         end
%
%%  element assembly into global matrix
      for krow=1:9
	  nrow=mv(:,krow);
          for kcol=1:9
		  ncol=mv(:,kcol);
          n = n + sparse(nrow,ncol,ne(:,krow,kcol),nvtx,nvtx);
	      end
       end
%
fprintf('done.\n')
return

author: johannesgerer

repo: jburkardt-m


text:
function [Y, sel] = vl_colsubset(X,n,varargin)
% VL_COLSUBSET Select a given number of columns
%   Y = VL_COLSUBSET(X, N) returns a random subset Y of N columns of
%   X. The selection is order-preserving and without replacement. If N
%   is larger or equal to the number of columns of X (e.g. N = Inf),
%   then the function returns all the columns (i.e., Y = X).
%
%   If 0 < N < 1, then the function returns a fraction N of the
%   columns (rounded to the closest integer).
%
%   [Y, SEL] = VL_COLSUBSET(...) returns the indexes SEL of the
%   selected columns.
%
%   The function accepts the following options:
%
%   Beginning::
%     Returns the fist N columns.
%
%   Ending::
%     Returns the last N columns.
%
%   Random:: [default]
%     Returns N columns selected at random (using RANDPERM()).
%
%   Uniform::
%     Returns N uniformly spaced columns.
%
%   Largest::
%     Returns the N largest columns (using SORTROWS()).
%
%   Smallest::
%     Returns the N smallest columns (using SORTROWS()).
%
%  See also: VL_HELP().

% Authors: Andrea Vedaldi

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

if nargin < 2, n = 1 ; end

mode = 'random' ;
i = 1 ;
while i <= length(varargin)
  switch lower(varargin{i})
    case {'beginning', ...
          'ending', ...
          'random', ...
          'uniform', ...
          'largest', ...
          'smallest'}
      mode = lower(varargin{1}) ;
      i = i + 1 ;
    otherwise
      error('Unknown option ''%s''.', varargin{i}) ;
  end
end

m = size(X,2) ;

if n < 0, error('N must not be smaller than 0.') ; end
if n ~= round(n)
  if n > 1
    error('N must be a natural number, +inf, or a fraction in 0 and 1.') ;
  end
  n = round(m * n) ;
end

n = min(m,n) ;

switch mode
  case 'random'
    perm = randperm(m) ;
    sel  = sort(perm(1:n)) ;
  case 'beginning'
    perm = 1:m ;
    sel  = sort(perm(1:n)) ;
  case 'ending'
    perm = m:-1:1 ;
    sel  = sort(perm(1:n)) ;
  case 'uniform'
    if n < 1
      sel = [] ;
    else
      sel = round(linspace(1, m, min(m,n))) ;
    end
  case 'largest'
    [drop, perm] = sortrows(X') ;
    sel = sort(perm(end-n+1:end)) ;
  case 'smallest'
    [drop, perm] = sortrows(X') ;
    sel = sort(perm(1:n)) ;
end

Y = X(:, sel) ;

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function y = vl_nndropout_wrapper(x, mask, test, dzdy)
%VL_NNDROPOUT_WRAPPER AutoNN wrapper for MatConvNet's vl_nndropout
%   This wrapper bypasses the dropout operation when in test mode.

% Copyright (C) 2017 Joao F. Henriques.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

  if test
    % pass-through input unchanged, when in test-mode
    if nargin < 4
      y = x ;
    else  % test-mode backward pass
      y = dzdy ;
    end
  else
    if nargin < 4
      y = vl_nndropout(x, 'mask', mask) ;
    else  % backward pass
      y = vl_nndropout(x, dzdy, 'mask', mask) ;
    end
  end
end

repo: MFT


text:
classdef dl_list < handle

% Copyright (C) 2012 Jean-Pierre de la Croix
% see the LICENSE file included with this software

    properties
        head_
        tail_
        size_
        iterator_
    end

    methods
        function obj = dl_list()
            obj.head_ = [];
            obj.tail_ = obj.head_;
            obj.size_ = 0;
        end

        function insert_key(obj, key, index)
            n = mcodekit.list.dl_list_node(key);
            index = max(min(index, obj.size_+1), 1);
            if (obj.size_ == 0)
                obj.head_ = n;
                obj.tail_ = obj.head_;
            else
                if (index == 1)
                    n.next_ = obj.head_;
                    obj.head_.prev_ = n;
                    obj.head_ = n;
                elseif (index == obj.size_+1)
                    n.prev_ = obj.tail_;
                    obj.tail_.next_ = n;
                    obj.tail_ = n;
                else
                    m = obj.head_;
                    for i = 2:index
                        m = m.next_;
                    end
                    n.prev_ = m.prev_;
                    m.prev_.next_ = n;
                    n.next_ = m;
                    m.prev_ = n;
                end
            end
            obj.size_ = obj.size_+1;
        end

        function key = remove_key(obj, index)
            index = max(min(index, obj.size_), 1);
            if (obj.size_ == 0)
                key = [];
            elseif (obj.size_ == 1)
                key = obj.tail_.key_;
                delete(obj.tail_);
                obj.tail_ = [];
                obj.head_ = obj.tail_;
                obj.size_ = 0;
            elseif (index == 1)
                key = obj.head_.key_;
                n = obj.head_;
                obj.head_ = n.next_;
                obj.head_.prev_ = [];
                delete(n);
                obj.size_ = obj.size_-1;
            elseif (index == obj.size_)
                key = obj.tail_.key_;
                n = obj.tail_;
                obj.tail_ = n.prev_;
                obj.tail_.next_ = [];
                delete(n);
                obj.size_ = obj.size_-1;
            else
                m = obj.head_;
                for i = 2:index
                    m = m.next_;
                end
                key = m.key_;
                m.prev_.next_ = m.next_;
                m.next_.prev_ = m.prev_;
                delete(m);
                obj.size_ = obj.size_-1;
            end
        end

        function append_key(obj, key)
            obj.insert_key(key, obj.size_+1);
        end

        function key = get_key(obj, index)
            assert(index > 0, 'index has to be positive');
            index = max(min(index, obj.size_), 1);
            n = obj.head_;
            for i = 2:index
                n = n.next_;
            end
            key = n.key_;
        end

        function key = get_first_key(obj)
            assert(obj.size_ > 0, 'list is empty');
            key = obj.head_.key_;
        end

        function key = get_last_key(obj)
            assert(obj.size_ > 0, 'list is empty');
            key = obj.tail_.key_;
        end

        function iterator = get_iterator(obj)
            iterator = mcodekit.list.dl_list_iterator(obj);
        end
    end
end

author: jdelacroix

repo: simiam

save_path: github-repos/MATLAB/jdelacroix-simiam


text:
function overlap_ratio = get_overlap_1toN(rect1, rect2)
    rect1 = repmat(rect1, [size(rect2, 1) 1]);

    area1 = (rect1(:, 3) - rect1(:, 1)) .* (rect1(:, 4) - rect1(:, 2));
    area2 = (rect2(:, 3) - rect2(:, 1)) .* (rect2(:, 4) - rect2(:, 2));

    l = max(rect1(:, 1), rect2(:, 1));
    r = min(rect1(:, 3), rect2(:, 3));
    t = max(rect1(:, 2), rect2(:, 2));
    b = min(rect1(:, 4), rect2(:, 4));

    w = r - l;
    h = b - t;
    overlap = w .* h;
    overlap(w < 0 | h < 0) = 0;

    overlap_ratio = overlap ./ (area1 + area2 - overlap);
end

author: feichtenhofer

repo: Detect-Track


text:
% demonstrate usage of classification
%
% Copyright (c) by Carl Edward Rasmussen and Hannes Nickisch 2013-10-16.
clear all, close all

%% SAY WHICH CODE WE WISH TO EXERCISE
id = [1,1]; % use Gauss/Exact
id = [1,2; 2,2; 3,2]; % compare Laplace
id = [1,3; 2,3; 3,3]; % study EP
id = [1,5; 2,5]; % look into KL (takes quite a while)
id = [1,4; 2,4; 3,4]; % deal with VB

seed = 943; randn('seed',seed), rand('seed',seed)

ntr = 50; nte = 1e4;                        % number of training and test points
xtr = 10*sort(rand(ntr,1));                                     % sample dataset
p = @(x) 1./(1+exp(-5*sin(x)));                  % "true" underlying probability
ytr = 2*(p(xtr)>rand(ntr,1))-1;                              % draw labels +1/-1
i = randperm(ntr); nout = 3;                                      % add outliers
ytr(i(1:nout)) = -ytr(i(1:nout));
xte = linspace(0,10,1e4)';                    % support, we test our function on
cov = {@covSEiso}; sf = 1; ell = 0.7;                             % setup the GP
hyp0.cov  = log([ell;sf]);
mean = {@meanZero};                                                   % m(x) = 0
hyp0.mean = [];
lik_list = {'likGauss','likErf','likLogistic'};          % allowable likelihoods
inf_list = {'infExact','infLaplace','infEP','infVB','infKL'};   % inference algs

Ncg = 50;                                   % number of conjugate gradient steps
sdscale = 0.5;                  % how many sd wide should the error bars become?
col = {'k',[.8,0,0],[0,.5,0],'b',[0,.75,.75],[.7,0,.5]};                % colors
ymu{1} = 2*p(xte)-1; ys2{1} = 0;
for i=1:size(id,1)
  lik = lik_list(id(i,1));                                % setup the likelihood
  if strcmp(lik,'likGauss')
    sn = .2; hyp0.lik = log(sn);
  else
    hyp0.lik = [];
  end
  inf = inf_list{id(i,2)};
  fprintf('OPT: %s/%s\n',lik_list{id(i,1)},inf_list{id(i,2)})
  hyp = minimize(hyp0,'gp', -Ncg, inf, mean, cov, lik, xtr, ytr);   % opt hypers
  [ymu{i+1}, ys2{i+1}] = gp(hyp, inf, mean, cov, lik, xtr, ytr, xte);  % predict
  [nlZ(i+1)] = gp(hyp, inf, mean, cov, lik, xtr, ytr);
end

figure, hold on
for i=1:size(id,1)+1
  plot(xte,ymu{i},'Color',col{i},'LineWidth',2)
  if i==1
    leg = {'function'};
  else
    leg{end+1} = sprintf('%s/%s -lZ=%1.2f',...
                                lik_list{id(i-1,1)},inf_list{id(i-1,2)},nlZ(i));
  end
end
for i=1:size(id,1)+1
  ysd = sdscale*sqrt(ys2{i});
  fill([xte;flipud(xte)],[ymu{i}+ysd;flipud(ymu{i}-ysd)],...
       col{i},'EdgeColor',col{i},'FaceAlpha',0.1,'EdgeAlpha',0.3);
end
for i=1:size(id,1)+1, plot(xte,ymu{i},'Color',col{i},'LineWidth',2), end
plot(xtr,ytr,'k+'), plot(xtr,ytr,'ko'), legend(leg)


text:
classdef CEC2017_F3 < PROBLEM
% <single> <real> <constrained>
% CEC'2017 constrained optimization benchmark problem

%------------------------------- Reference --------------------------------
% G. Wu, R. Mallipeddi, and P. N. Suganthan, Problem definitions and
% evaluation criteria for the CEC 2017 competition on constrained real-
% parameter optimization, National University of Defense Technology, China,
% 2016.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        O;  % Optimal decision vector
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'CEC2017.mat'),'Data');
            obj.O = Data{3};
            obj.M = 1;
            if isempty(obj.D); obj.D = 10; end
            obj.D = min(obj.D,length(obj.O));
            obj.lower    = zeros(1,obj.D) - 100;
            obj.upper    = zeros(1,obj.D) + 100;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            Z = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            PopObj = sum(cumsum(Z,2).^2,2);
        end
        %% Calculate constraint violations
        function PopCon = CalCon(obj,PopDec)
            Z = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            PopCon(:,1) = sum(Z.^2-5000*cos(0.1*pi*Z)-4000,2);
            PopCon(:,2) = abs(sum(Z.*sin(0.1*pi*Z),2)) - 1e-4;
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [power] = VBA_power(type, alpha, Estat, dof)
% // VBA toolbox //////////////////////////////////////////////////////////
%
% [power] = VBA_power(type, alpha, Estat, dof)
% Statistical power of a t- or F-test
%
% IN:
%   - type: type of test. Can be set to 't' or 'F'.
%   - alpha: statistical significance level
%   - Estat: expected effect size (z-score for a t-test, F value for an
%            F-test)
%   - dof: degrees of freedom
%
% OUT:
%   - power: ensuing statistical power
%
% /////////////////////////////////////////////////////////////////////////

%  ensure numerical stability
dof(isinf(dof)) = 1e8;
dof(dof==0) = 1e-8;

switch type
    case 't'
        tcritic = VBA_spm_invTcdf (1 - alpha, dof);
        power = 1 - VBA_spm_Tcdf (tcritic - Estat, dof);
    case 'F'
        Fcritic = VBA_spm_invFcdf (1 - alpha, dof(1), dof(2));
        power = 1-VBA_ncfcdf (Fcritic, dof(1), dof(2), Estat);
    otherwise
        disp('*** VBA_power: ''type'' argument can only be ''t'' or ''F''')
        power = [];
        return
end

author: MBB-team

repo: VBA-toolbox


text:
%% freeBoundaryPatch
% Below is a demonstration of the features of the |freeBoundaryPatch| function

%%
clear; close all; clc;

%% Syntax
% |[indFree]=freeBoundaryPatch(F);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [ daten f107 ] = read_solarflux( filename )

% INPUT
% 'filename'  : full path of solar flux observation file
% i.e. ...
% \STP\SOLAR_DATA\SOLAR_RADIO\FLUX\Penticton_Observed\daily\DAILYPLT.OBS

% OUTPUTS
% 'daten'     : array of MATLAB datenum dates
% 'f107'      : corresponding solar flux values

% Reads solar flux data files from Penticton Observatory
% daily F10.7 flux values, downloaded from ftp.ngdc.noaa.gov

% Author : John A. Smith, CIRES/NOAA, Univ. of Colorado at Boulder

fid = fopen(filename,'r','l','US-ASCII'); % open file

% check file open successful
if fid==-1
    error('Could not open file')
end

dat=textscan(fid,'%4d %2d %2d %f','treatasempty','.'); % read file contents

fclose(fid); % close file

% assign year, month and day
year=dat{1};
month=dat{2};
day=dat{3};

daten=datenum(double([year month day])); % convert to MATLAB datenum

f107=dat{4}; % assign f10.7 data

f107(f107==0)=NaN; % NaN empty fields

end


text:
% area of ellipse
function [data,units] = compute_area(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);
  data{i} = (2*trx(fly).a_mm).*(2*trx(fly).b_mm)*pi;
end
units = parseunits('mm^2');


text:
function [ r, s, area ] = node_reference ( code )

%*****************************************************************************80
%
%% NODE_REFERENCE returns the basis nodes for any available element.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character CODE(*), identifies the element desired.
%    Legal values include 'Q4', 'Q8', 'Q9', 'Q12', 'Q16', 'QL',
%    'T3', 'T4', 'T6' and 'T10'.
%
%    Output, real R(N), S(N), the coordinates of the basis nodes.
%
%    Output, real AREA, the area of the element.
%
  if ( s_eqi ( code, 'Q4' ) )
    [ r, s, area ] = node_reference_q4 ( );
  elseif ( s_eqi ( code, 'Q8' ) )
    [ r, s, area ] = node_reference_q8 ( );
  elseif ( s_eqi ( code, 'Q9' ) )
    [ r, s, area ] = node_reference_q9 ( );
  elseif ( s_eqi ( code, 'Q12' ) )
    [ r, s, area ] = node_reference_q12 ( );
  elseif ( s_eqi ( code, 'Q16' ) )
    [ r, s, area ] = node_reference_q16 ( );
  elseif ( s_eqi ( code, 'QL' ) )
    [ r, s, area ] = node_reference_ql ( );
  elseif ( s_eqi ( code, 'T3' ) )
    [ r, s, area ] = node_reference_t3 ( );
  elseif ( s_eqi ( code, 'T4' ) )
    [ r, s, area ] = node_reference_t4 ( );
  elseif ( s_eqi ( code, 'T6' ) )
    [ r, s, area ] = node_reference_t6 ( );
  elseif ( s_eqi ( code, 'T10' ) )
    [ r, s, area ] = node_reference_t10 ( );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'NODE_REFERENCE - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of CODE = "%s"\n', code );
    error ( 'NODE_REFERENCE - Fatal error!' );
  end

  return
end


text:
function [ y2, m2, d2, f2 ] = yjf_to_ymdf_common ( y1, j1, f1 )

%*****************************************************************************80
%
%% YJF_TO_YMDF_COMMON converts a Common date from YJF to YMDF format.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y1, J1, real F1, the YJF date.
%
%    Output, integer Y2, M2, D2, real F2,
%    the YMDF date.
%

%
%  Copy the input.
%
  y2 = y1;
  j2 = j1;
  f2 = f1;
%
%  Check the input.
%
  [ y2, j2, f2, ierror ] = yjf_check_common ( y2, j2, f2 );

  if ( ierror ~= 0 )
    y2 = 0;
    m2 = 0;
    d2 = 0;
    f2 = 0.0;
    return
  end
%
%  Convert the input.
%
  d2 = j2;
  m2 = 1;

  [ y2, m2, d2 ] = day_borrow_common ( y2, m2, d2 );

  [ y2, m2, d2 ] = day_carry_common ( y2, m2, d2 );

  return
end


text:
% mff_fileio_read_header() - import MFF header
%
% Usage:
%   >> header = read_eeglabheader(filename);
%
% Inputs:
%   filename - [string] file name
%
% Outputs:
%   header   - FILEIO toolbox type structure

% This file is part of mffmatlabio.
%
% mffmatlabio 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.
%
% mffmatlabio 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 mffmatlabio.  If not, see <https://www.gnu.org/licenses/>.

function header = mff_fileio_read_header(mffFile)

if nargin < 1
    help mff_fileio_read_header;
    return;
end

% import basic information
tmp = mff_import(mffFile);
curPath = pwd;
p = fileparts(which('ft_read_header'));
cd(fullfile(p, 'private'));
header = read_eeglabheader( tmp );
cd(curPath);
header.orig = tmp;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ox, otheta, ollh, olpp, v] = tapas_ti_step_mh(y, u, ox, otheta, ...
    ollh, olpp, ptheta, htheta, pars)
%% Estimates the posterior probability of the parameters using MCMC combined
% with path sampling.
%
% Input
%   y       Saccade data.
%   u       Experimental input.
%   ptheta  Priors of the model
%   htheta  Kernel transformation.
%   pars    Parameters for the MCMC.
%           pars.verbose   True or false. Defs. False.
%           pars.mc3it      Number of iterations for MC3. Defs. 0.
%           pars.kup        Number of steps for the kernel update.
%                           Defs. 500.
%           pars.seed       Seed for the random generation. If zero use
%                           rng('shuffle'). Defaults to zero.
%           pars.samples    If true, stores the samples from the chain at
%                           lowest temperature. Defaults to zero.
%
% Output
%   ps      Samples from the posterior distribuion
%   fe      Free energy estimate.
%
% Uses a standard MCMC on a population and applies an exchange operator
% to improve the mixing.
%
%
% Real parameter evolutionary Monte Carlo with applications to Bayes Mixture
% Models, Journal of American Statistica Association, Liang & Wong 2001.

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

if isfield(ptheta, 'mh_propose_sample')
    propose_sample = htheta.mh_propose_sample;
else
    propose_sample = @tapas_ti_propose_gaussian;
end

[ntheta, ratio] = propose_sample(otheta, ptheta, htheta);
[nllh, nx] = ptheta.llh(y, [], u, ntheta, ptheta);

T = ptheta.T;

nllh = sum(nllh, 1);
nlpp = sum(ptheta.lpp(y, nx, u, ntheta, ptheta), 1);

v = nllh .* T + nlpp - (ollh .* T + olpp) + ratio;
nansv = isnan(v);
v(nansv) = -inf;

v = rand(size(v)) < exp(bsxfun(@min, v, 0));

ollh(v) = nllh(v);
olpp(v) = nlpp(v);
otheta(:, v) = ntheta(:, v);
ox(:, v) = nx(:, v);

assert(all(-inf < ollh), 'mpdcm:ti:mh', '-inf value in the new samples');

end

author: translationalneuromodeling

repo: tapas


text:
function S= bbciutil_transformProc2FcnParam(S)
%TRANSFORMPROC2FCNPARAM - Transform 'proc' field to fcn/param fields
%
%This function just transforms the way a sequence of processing steps
%is given. This is best exmplained by an example:
% S.proc=  {{'proc_baseline', [-200 0]}, ...
%           {'proc_jumpingMeans', 5, 10}};
% T.fcn= {'proc_baseline', 'proc_jumpingMeans'}
% T.param= {{[-200 0]}, {5, 10}};
%The 'proc' representation is more intuitive for definition, but the
%fcn/param representation is better for performance. There the representation
%is transform initially in bbci_apply.
%
%Synopsis:
%  T= bbciutil_transformProc2FcnParam(S)
%
%Arguments:
%  S - Structure which has the proc
%
%Output:
%  T - Updated structure

% 02-2011 Benjamin Blankertz

if ~isfield(S, 'proc'),
  for k= 1:length(S),
    if (~isfield(S,'param') || isempty(S(k).param)) && ...
          isfield(S(k), 'fcn') && ~isempty(S(k).fcn),
      S(k).param= repmat({{}}, [1 length(S(k).fcn)]);
    end
  end
  return;
end

for k= 1:length(S),
  if isfield(S(k),'fcn') && ~isempty(S(k).fcn),
%    warning('struct has already nonempty field fcn');
    continue;
  end
  nProcs= length(S(k).proc);
  S(k).fcn= cell(1, nProcs);
  S(k).param= cell(1, nProcs);
  for p= 1:nProcs,
    [S(k).fcn{p}, S(k).param{p}]= misc_getFuncParam(S(k).proc{p});
  end
end
S= rmfield(S, 'proc');

author: bbci


text:
function clsts= getClusters(net, opts, clstFn, k, dbTrain, trainDescFn)

    if ~exist(clstFn, 'file')

        if ~exist(trainDescFn, 'file')

            simpleNnOpts= {'conserveMemory', true, 'mode', 'test'};

            if opts.useGPU
                net= relja_simplenn_move(net, 'gpu');
            end

            % ---------- extract training descriptors

            relja_display('Computing training descriptors');

            nTrain= 50000;
            nPerImage= 100;
            nIm= ceil(nTrain/nPerImage);

            rng(43);
            trainIDs= randsample(dbTrain.numImages, nIm);

            nTotal= 0;

            prog= tic;

            for iIm= 1:nIm
                relja_progress(iIm, nIm, 'extract train descs', prog);

                % --- extract descriptors

                % didn't want to complicate with batches here as it's only done once (per network and training set)

                im= vl_imreadjpeg({[dbTrain.dbPath, dbTrain.dbImageFns{iIm}]});
                im= im{1};

                % fix non-colour images
                if size(im,3)==1
                    im= cat(3,im,im,im);
                end

                im(:,:,1)= im(:,:,1) - net.meta.normalization.averageImage(1,1,1);
                im(:,:,2)= im(:,:,2) - net.meta.normalization.averageImage(1,1,2);
                im(:,:,3)= im(:,:,3) - net.meta.normalization.averageImage(1,1,3);

                if opts.useGPU
                    im= gpuArray(im);
                end

                res= vl_simplenn(net, im, [], [], simpleNnOpts{:});
                descs= gather(res(end).x);
                descs= reshape( descs, [], size(descs,3) )';

                % --- sample descriptors

                nThis= min( min(nPerImage, size(descs,2)), nTrain - nTotal );
                descs= descs(:, randsample( size(descs,2), nThis ) );

                if iIm==1
                    trainDescs= zeros( size(descs,1), nTrain, 'single' );
                end

                trainDescs(:, nTotal+[1:nThis])= descs;
                nTotal= nTotal+nThis;
            end

            trainDescs= trainDescs(:, 1:nTotal);

            % move back to CPU addLayers() assumes it
            if opts.useGPU
                net= relja_simplenn_move(net, 'cpu');
            end

            save(trainDescFn, 'trainDescs');
        else
            relja_display('Loading training descriptors');
            load(trainDescFn, 'trainDescs');
        end

        % ---------- Cluster descriptors

        relja_display('Computing clusters');
        clsts= yael_kmeans(trainDescs, k, 'niter', 100, 'verbose', 0, 'seed', 43);
        clear trainDescs;

        save(clstFn, 'clsts');
    else
        relja_display('Loading clusters');
        load(clstFn, 'clsts');
        assert(size(clsts, 2)==k);
    end

end

author: Relja

repo: netvlad

save_path: github-repos/MATLAB/Relja-netvlad


text:
function model = ml_trainlda(varargin)
% Learn a linear predictive model by (regularized) Linear Discriminant Analysis.
% Model = ml_trainlda(Trials, Targets, Lambda, Options...)
%
% LDA is one of the simplest and oldest learning algorithms, originally introduced by Fisher in [1].
% Its basic assumption is that the data originates from two classes (extended to more via 1-vs-1
% voting), where the data in each class is distributed in the feature space according to a normal
% distribution. The shape of the distribution is assumed to be identical for both classes, and the
% relative weight/prior probability of each class is (by default) assumed to be identical, as well.
% An example would be data with features generated as a weighted sum of some latent variables which
% take on different values between conditions (but identical values within each condition), and
% which are superimposed with gaussian noise (which is identically distributed across conditions,
% e.g. from a large sum of independent random variables).
%
% Despite its simplicity, LDA assumes more structure in the data than is usually necessary (namely a
% certain normal distribution per class), and sometimes, more than what can be satisfactorily
% learned from the data, so that, even when the assumptions are fulfilled, the method is not
% guaranteed to give the optimal result. The estimation of the per-class covariance matrix is a
% notoriously data-hungry (and especially outlier-sensitive) step, and the main weakness of standard
% LDA.
%
% ml_trainlda offers three advanced covariance estimators, with different trade-offs to mitigate
% these problems. The variants 'shrinkage' and 'independence' each introduce a regularization
% parameter [2] which controls the complexity of the estimated matrix, and which need to be selected
% in a parameter search (which is orders of magnitude more time-consuming). The 'auto' variant
% computes the degree of regularization analyically in closed form, and is therefore both fast and
% (in some sense) optimal (but more restricted than 'independence') [3]. When enough trials are
% available, full covariance matrices are learned, but the less trials are given, the more the
% covariance estimates degrade to spherical (though well-formed) ones. Auto is the default for lda.
%
% While regularization can automatically control the complexity of a classifier, it is not a panacea
% which allows to add arbitrarily many features, since with each additional feature, the amount of
% structure (here:correlation) that can be captured in the remaining ones gets reduced.
%
% Within the toolbox, LDA is one of the bread-and-butter classifiers, and is worth trying in every
% reasonably simple classification task, especially for its speed. The major problem with LDA
% compared to other available classifiers is that it easier to break it with outliers than, for
% example, support vector machines or logistic regression. Another weakness is that the outputs are
% relatively primitive probability estimates, in contrast to, for example, logistic regression
% (somewhat better) or relevance vector machines (clearly better).
%
% In:
%   Trials       : training data, as in ml_train
%
%   Targets      : target variable, as in ml_train
%
%   Lambda       : optional regularization parameter, reasonable range: 0:0.1:1, greater is stronger
%                  requires that the regularization mode is set to either 'shrinkage' or 'independence' (default: [])
%
%   Options  : optional name-value parameters to control the training details:
%              'regularization' -> 'shrinkage': covariance shrinkage, depends on plambda
%                                  'independence': feature independence, depends on plambda
%                                  'auto': analytical covariance shrinkage, plambda is ignored (default)
%              'weight_bias' -> 0/1, take unequal class priors into account for bias calculation
%                               default: 0
%              'weight_cov' -> 0/1, take unequal class priors into account for covariance calculation
%                              default: 0
% Out:
%   Model   : a linear model; w is the linear weights, b is the bias; classes indicates the class labels which the model predicts
%
% Examples:
%   % learn a standard shrinkage LDA model
%   model = ml_trainlda(trials,targets);
%
%   % take unequal class priors into account for both the bias and the covariance matrix
%   model = ml_trainlda(trials,targets,[],'weight_bias',1,'weight_cov',1);
%
%   % use a different type of regularization, which controls feature independence and requires cross-validation
%   model = utl_searchmodel({trials,target},'args',{{'lda',search(0:0.1:1),'regularization','independence'}})
%
%
% See also:
%   ml_predictlda
%
% References:
%   [1] Fisher, R. "The use of multiple measurements in taxonomic problems."
%       Annals Eugen. 7 (1936), 188, 179.
%   [2] Friedman, J. "Regularized discriminant analysis."
%       Journal of the American Statistical Association 84, 405 (1989), 175, 165.
%   [3] O. Ledoit and M. Wolf, "A well-conditioned estimator for large-dimensional covariance matrices"
%       J Multivar Anal, 88(2): 365-411, 2004.
%
%                           Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                           2010-04-03
dp;

arg_define([0 3],varargin, ...
    arg_norep('trials'), ...
    arg_norep('targets'), ...
    arg({'plambda','Lambda','lambda'}, [], [0 1], 'Optional regularization parameter. Reasonable range: 0:0.1:1 - greater is stronger. Requires that the regularization mode is set to either "shrinkage" or "independence" (not necessary in "auto" mode).'), ...
    arg({'regularization','Regularizer','Regularization'}, 'auto', {'none','auto','shrinkage','independence'}, 'Type of regularization. Regularizes the robustness / flexibility of covariance estimates. Auto is analytical covariance shrinkage, shrinkage is shrinkage as selected via plambda, and independence is feature independence, also selected via plambda.'), ...
    arg({'weight_bias','WeightedBias'}, false, [], 'Account for class priors in bias. If you do have unequal probabilities for the different classes, this should be enabled.'), ...
    arg({'weight_cov','WeightedCov'}, false, [], 'Account for class priors in covariance. If you do have unequal probabilities for the different classes, it makes sense to enable this.'), ...
    arg({'robust','Robust'}, false, [], 'Use robust estimation. Uses geometric medians in place of means; can help if some trials are very noisy.'), ...
    arg({'votingScheme','VotingScheme'},'1vR',{'1v1','1vR'},'Voting scheme. If multi-class classification is used, this determine how binary classifiers are arranged to solve the multi-class problem. 1v1 gets slow for large numbers of classes (as all pairs are tested), but can be more accurate than 1vR.'));

% find the class labels
classes = unique(targets);
if length(classes) > 2
    % learn a voting arrangement of models...
    model = ml_trainvote(trials, targets, votingScheme, @ml_trainlda, @ml_predictlda, varargin{:},'weight_bias',true);    %#ok<*NODEF>
elseif length(classes) == 1
    error('BCILAB:only_one_class','Your training data set has no trials for one of your classes; you need at least two classes to train a classifier.\n\nThe most likely reasons are that one of your target markers does not occur in the data, or that all your trials of a particular class are concentrated in a single short segment of your data (10 or 20 percent). The latter would be a problem with the experiment design.');
else
    % pre-prune degenerate features
    retain = true(1,size(trials,2));
    for c = 1:length(classes)
        X = trials(targets==classes(c),:);
        n{c} = size(X,1);
        mu{c} = mean(X,1);
        v{c} = var(X,[],1);
        retain = retain & isfinite(mu{c}) & isfinite(v{c}) & (n{c}==1 | v{c} > eps);
    end
    % apply feature mask...
    trials = trials(:,retain);
    % estimate distribution of each class...
    for c = 1:length(classes)
        lams{c} = NaN;
        X = trials(targets==classes(c),:);
        n{c} = size(X,1);
        if robust
            mu{c} = geometric_median(X);
        else
            mu{c} = mean(X,1);
        end
        if n{c} == 1
            sig{c} = zeros(size(X,2));
        elseif strcmp(regularization,'auto')
            if robust
                % for lack of a better solution in this case we estimate lambda using a non-robust
                % Ledoit-Wolf estimator and then use it in a subsequent robust re-estimation
                [dummy,lams{c}] = cov_shrink(X); %#ok<ASGLU>
                sig{c} = cov_blockgeom(X,1);
                sig{c} = (1-lams{c})*sig{c} + lams{c}*eye(length(sig{c}))*abs(mean(diag(sig{c})));
            else
                [sig{c},lams{c}] = cov_shrink(X);
            end
        else
            if robust
                sig{c} = cov_blockgeom(X,1);
            else
                sig{c} = cov(X);
            end
            if ~isempty(plambda) && ~strcmp(regularization,'none')
                % plambda-dependent regularization
                if strcmp(regularization,'independence')
                    sig{c} = (1-plambda)*sig{c} + plambda * diag(diag(sig{c}));
                elseif strcmp(regularization,'shrinkage')
                    sig{c} = (1-plambda)*sig{c} + plambda*eye(length(sig{c}))*abs(mean(diag(sig{c})));
                else
                    error('unknown regularization mode');
                end
            end
        end
    end

    ns = quickif(weight_cov,n,{1 1});
    nb = quickif(weight_bias,n,{1 1});
    % do the math
    mu_both = (mu{1}*nb{2} + mu{2}*nb{1}) / (nb{1}+nb{2});
    sig_both = (sig{1}*ns{1} + sig{2}*ns{2}) / (ns{1}+ns{2});
    w = (mu{2} - mu{1}) / sig_both;
    w = w / (mu{2}*w' - mu_both*w');
    model = struct('w',{w}, 'b',{mu_both*w'}, 'classes',{classes},'featuremask',{retain},'lams',{lams});
    1;
end

repo: REST


text:
function gpcf = gpcf_noise(varargin)
%GPCF_NOISE  Create a independent noise covariance function
%
%  Description
%    GPCF = GPCF_NOISE('PARAM1',VALUE1,'PARAM2,VALUE2,...) creates
%    independent noise covariance function structure in which the
%    named parameters have the specified values. Any unspecified
%    parameters are set to default values.
%
%    GPCF = GPCF_NOISE(GPCF,'PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    modify a covariance function structure with the named
%    parameters altered with the specified values.
%
%    Parameters for independent noise covariance function [default]
%      noiseSigma2       - variance of the independent noise [0.1]
%      noiseSigma2_prior - prior for noiseSigma2 [prior_logunif]
%
%    Note! If the prior is 'prior_fixed' then the parameter in
%    question is considered fixed and it is not handled in
%    optimization, grid integration, MCMC etc.
%
%  See also
%    GP_SET, GPCF_*, PRIOR_*

% Copyright (c) 2007-2010 Jarno Vanhatalo
% Copyright (c) 2010 Aki Vehtari

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

  ip=inputParser;
  ip.FunctionName = 'GPCF_NOISE';
  ip.addOptional('gpcf', [], @isstruct);
  ip.addParamValue('noiseSigma2',0.1, @(x) isscalar(x) && x>0);
  ip.addParamValue('noiseSigma2_prior',prior_logunif, @(x) isstruct(x) || isempty(x));
  ip.parse(varargin{:});
  gpcf=ip.Results.gpcf;

  if isempty(gpcf)
    init=true;
    gpcf.type = 'gpcf_noise';
  else
    if ~isfield(gpcf,'type') && ~isequal(gpcf.type,'gpcf_noise')
      error('First argument does not seem to be a valid covariance function structure')
    end
    init=false;
  end

  % Initialize parameter
  if init || ~ismember('noiseSigma2',ip.UsingDefaults)
    gpcf.noiseSigma2=ip.Results.noiseSigma2;
  end

  % Initialize prior structure
  if init
    gpcf.p=[];
  end
  if init || ~ismember('noiseSigma2_prior',ip.UsingDefaults)
    gpcf.p.noiseSigma2=ip.Results.noiseSigma2_prior;
  end

  if init
    % Set the function handles to the subfunctions
    gpcf.fh.pak = @gpcf_noise_pak;
    gpcf.fh.unpak = @gpcf_noise_unpak;
    gpcf.fh.lp = @gpcf_noise_lp;
    gpcf.fh.lpg = @gpcf_noise_lpg;
    gpcf.fh.cfg = @gpcf_noise_cfg;
    gpcf.fh.ginput = @gpcf_noise_ginput;
    gpcf.fh.cov = @gpcf_noise_cov;
    gpcf.fh.trcov  = @gpcf_noise_trcov;
    gpcf.fh.trvar  = @gpcf_noise_trvar;
    gpcf.fh.recappend = @gpcf_noise_recappend;
  end

end

function [w, s] = gpcf_noise_pak(gpcf)
%GPCF_NOISE_PAK  Combine GP covariance function parameters into
%                one vector.
%
%  Description
%    W = GPCF_NOISE_PAK(GPCF) takes a covariance function data
%    structure GPCF and combines the covariance function
%    parameters and their hyperparameters into a single row
%    vector W. This is a mandatory subfunction used for example
%    in energy and gradient computations.
%
%       w = [ log(gpcf.noiseSigma2)
%             (hyperparameters of gpcf.magnSigma2)]'
%
%
%  See also
%    GPCF_NOISE_UNPAK

  w = []; s = {};
  if ~isempty(gpcf.p.noiseSigma2)
    w(1) = log(gpcf.noiseSigma2);
    s = [s 'log(noise.noiseSigma2)'];
    % Hyperparameters of noiseSigma2
    [wh sh] = gpcf.p.noiseSigma2.fh.pak(gpcf.p.noiseSigma2);
    w = [w wh];
    s = [s sh];
  end

end

function [gpcf, w] = gpcf_noise_unpak(gpcf, w)
%GPCF_NOISE_UNPAK  Sets the covariance function parameters
%                     into the structure
%
%  Description
%    [GPCF, W] = GPCF_NOISE_UNPAK(GPCF, W) takes a covariance
%    function data structure GPCF and a hyper-parameter vector W,
%    and returns a covariance function data structure identical
%    to the input, except that the covariance hyper-parameters
%    have been set to the values in W. Deletes the values set to
%    GPCF from W and returns the modified W. This is a mandatory
%    subfunction used for example in energy and gradient computations.
%
%    Assignment is inverse of
%       w = [ log(gpcf.noiseSigma2)
%             (hyperparameters of gpcf.magnSigma2)]'
%
%  See also
%    GPCF_NOISE_PAK

  if ~isempty(gpcf.p.noiseSigma2)
    gpcf.noiseSigma2 = exp(w(1));
    w = w(2:end);

    % Hyperparameters of lengthScale
    [p, w] = gpcf.p.noiseSigma2.fh.unpak(gpcf.p.noiseSigma2, w);
    gpcf.p.noiseSigma2 = p;
  end
end

function lp = gpcf_noise_lp(gpcf)
%GPCF_NOISE_LP  Evaluate the log prior of covariance function parameters
%
%  Description
%    LP = GPCF_NOISE_LP(GPCF) takes a covariance function
%    structure GPCF and returns log(p(th)), where th collects the
%    parameters. This is a mandatory subfunction used for example
%    in energy computations.
%
%  See also
%   GPCF_NOISE_PAK, GPCF_NOISE_UNPAK, GPCF_NOISE_G, GP_E

% Evaluate the prior contribution to the error. The parameters that
% are sampled are from space W = log(w) where w is all the
% "real" samples. On the other hand errors are evaluated in the
% W-space so we need take into account also the Jacobian of
% transformation W -> w = exp(W). See Gelman et.al., 2004,
% Bayesian data Analysis, second edition, p24.

  lp = 0;
  gpp=gpcf.p;

  if ~isempty(gpcf.p.noiseSigma2)
    % Evaluate the prior contribution to the error.
    lp = gpp.noiseSigma2.fh.lp(gpcf.noiseSigma2, gpp.noiseSigma2) +log(gpcf.noiseSigma2);
  end
end

function lpg = gpcf_noise_lpg(gpcf)
%GPCF_NOISE_LPG  Evaluate gradient of the log prior with respect
%                to the parameters.
%
%  Description
%    LPG = GPCF_NOISE_LPG(GPCF) takes a covariance function
%    structure GPCF and returns LPG = d log (p(th))/dth, where th
%    is the vector of parameters. This is a mandatory subfunction
%    used for example in gradient computations.
%
%  See also
%    GPCF_NOISE_PAK, GPCF_NOISE_UNPAK, GPCF_NOISE_LP, GP_G

  lpg = [];
  gpp=gpcf.p;

  if ~isempty(gpcf.p.noiseSigma2)
    lpgs = gpp.noiseSigma2.fh.lpg(gpcf.noiseSigma2, gpp.noiseSigma2);
    lpg = [lpg lpgs(1).*gpcf.noiseSigma2+1 lpgs(2:end)];
  end
end

function DKff = gpcf_noise_cfg(gpcf, x, x2, mask, i1)
%GPCF_NOISE_CFG  Evaluate gradient of covariance function
%                with respect to the parameters
%
%  Description
%    DKff = GPCF_NOISE_CFG(GPCF, X) takes a covariance function
%    data structure GPCF, a matrix X of input vectors and returns
%    DKff, the gradients of covariance matrix Kff = k(X,X) with
%    respect to th (cell array with matrix elements). This is a
%    mandatory subfunction used in gradient computations.
%
%    DKff = GPCF_NOISE_CFG(GPCF, X, X2) takes a covariance
%    function data structure GPCF, a matrix X of input vectors
%    and returns DKff, the gradients of covariance matrix Kff =
%    k(X,X2) with respect to th (cell array with matrix
%    elements). This subfunction is needed when using sparse
%    approximations (e.g. FIC).
%
%    DKff = GPCF_NOISE_CFG(GPCF, X, [], MASK) takes a covariance
%    function data structure GPCF, a matrix X of input vectors
%    and returns DKff, the diagonal of gradients of covariance
%    matrix Kff = k(X,X2) with respect to th (cell array with
%    matrix elements). This subfunction is needed when using
%    sparse approximations (e.g. FIC).
%
%  See also
%    GPCF_NOISE_PAK, GPCF_NOISE_UNPAK, GPCF_NOISE_E, GP_G

  DKff = {};

  if ~isempty(gpcf.p.noiseSigma2)
    gpp=gpcf.p;
    DKff{1}=gpcf.noiseSigma2;
  end
  if nargin==4
    % Use memory save option
    if i1==0
      % Return number of hyperparameters
      DKff=1;
      return
    end
    DKff=DKff{1};
  end

end

function DKff = gpcf_noise_ginput(gpcf, x, t, i1)
%GPCF_NOISE_GINPUT  Evaluate gradient of covariance function with
%                   respect to x
%
%  Description
%    DKff = GPCF_NOISE_GINPUT(GPCF, X) takes a covariance
%    function data structure GPCF, a matrix X of input vectors
%    and returns DKff, the gradients of covariance matrix Kff =
%    k(X,X) with respect to X (cell array with matrix elements).
%    This subfunction is needed when computing gradients with
%    respect to inducing inputs in sparse approximations.
%
%    DKff = GPCF_NOISE_GINPUT(GPCF, X, X2) takes a covariance
%    function data structure GPCF, a matrix X of input vectors
%    and returns DKff, the gradients of covariance matrix Kff =
%    k(X,X2) with respect to X (cell array with matrix elements).
%    This subfunction is needed when computing gradients with
%    respect to inducing inputs in sparse approximations.
%
%  See also
%    GPCF_NOISE_PAK, GPCF_NOISE_UNPAK, GPCF_NOISE_E, GP_G

end

function C = gpcf_noise_cov(gpcf, x1, x2)
% GP_NOISE_COV     Evaluate covariance matrix between two input vectors
%
%  Description
%    C = GP_NOISE_COV(GP, TX, X) takes in covariance function of
%    a Gaussian process GP and two matrixes TX and X that contain
%    input vectors to GP. Returns covariance matrix C. Every
%    element ij of C contains covariance between inputs i in TX
%    and j in X. This is a mandatory subfunction used for example in
%    prediction and energy computations.
%
%  See also
%    GPCF_NOISE_TRCOV, GPCF_NOISE_TRVAR, GP_COV, GP_TRCOV

  if isempty(x2)
    x2=x1;
  end
  [n1,m1]=size(x1);
  [n2,m2]=size(x2);

  if m1~=m2
    error('the number of columns of X1 and X2 has to be same')
  end

  C = sparse([],[],[],n1,n2,0);
end

function C = gpcf_noise_trcov(gpcf, x)
%GP_NOISE_TRCOV  Evaluate training covariance matrix of inputs
%
%  Description
%    C = GP_NOISE_TRCOV(GP, TX) takes in covariance function of a
%    Gaussian process GP and matrix TX that contains training
%    input vectors. Returns covariance matrix C. Every element ij
%    of C contains covariance between inputs i and j in TX. This is
%    a mandatory subfunction used for example in prediction and
%    energy computations.
%
%
%  See also
%    GPCF_NOISE_COV, GPCF_NOISE_TRVAR, GP_COV, GP_TRCOV

  [n, m] =size(x);
  n1=n+1;

  C = sparse([],[],[],n,n,0);
  C(1:n1:end)=C(1:n1:end)+gpcf.noiseSigma2;

end

function C = gpcf_noise_trvar(gpcf, x)
% GP_NOISE_TRVAR     Evaluate training variance vector
%
%    Description
%    C = GP_NOISE_TRVAR(GPCF, TX) takes in covariance function
%    of a Gaussian process GPCF and matrix TX that contains
%    training inputs. Returns variance vector C. Every
%    element i of C contains variance of input i in TX. This is
%    a mandatory subfunction used for example in prediction and
%    energy computations.
%
%
%    See also
%    GPCF_NOISE_COV, GP_COV, GP_TRCOV

  [n, m] =size(x);
  C=ones(n,1)*gpcf.noiseSigma2;

end

function reccf = gpcf_noise_recappend(reccf, ri, gpcf)
%RECAPPEND Record append
%
%  Description
%    RECCF = GPCF_NOISE_RECAPPEND(RECCF, RI, GPCF) takes a
%    covariance function record structure RECCF, record index RI
%    and covariance function structure GPCF with the current MCMC
%    samples of the hyperparameters. Returns RECCF which contains
%    all the old samples and the current samples from GPCF.
%    This subfunction is needed when using MCMC sampling (gp_mc).
%
%  See also
%    GP_MC and GP_MC -> RECAPPEND

  if nargin == 2
    % Initialize the record
    reccf.type = 'gpcf_noise';

    % Initialize parameters
    reccf.noiseSigma2 = [];

    % Set the function handles
    reccf.fh.pak = @gpcf_noise_pak;
    reccf.fh.unpak = @gpcf_noise_unpak;
    reccf.fh.e = @gpcf_noise_lp;
    reccf.fh.lpg = @gpcf_noise_lpg;
    reccf.fh.cfg = @gpcf_noise_cfg;
    reccf.fh.cov = @gpcf_noise_cov;
    reccf.fh.trcov  = @gpcf_noise_trcov;
    reccf.fh.trvar  = @gpcf_noise_trvar;
    %  gpcf.fh.sampling = @hmc2;
    reccf.sampling_opt = hmc2_opt;
    reccf.fh.recappend = @gpcf_noise_recappend;
    reccf.p=[];
    reccf.p.noiseSigma2=[];
    if ~isempty(ri.p.noiseSigma2)
      reccf.p.noiseSigma2 = ri.p.noiseSigma2;
    end
  else
    % Append to the record
    gpp = gpcf.p;

    % record noiseSigma2
    reccf.noiseSigma2(ri,:)=gpcf.noiseSigma2;
    if isfield(gpp,'noiseSigma2') && ~isempty(gpp.noiseSigma2)
      reccf.p.noiseSigma2 = gpp.noiseSigma2.fh.recappend(reccf.p.noiseSigma2, ri, gpcf.p.noiseSigma2);
    end
  end
end