text:
classdef StaticSaliencyFineGrained < handle
    %STATICSALIENCYFINEGRAINED  The Fine Grained Saliency approach for Static Saliency
    %
    % The Fine Grained Saliency approach from [FGS].
    %
    % ## Saliency API
    %
    % Many computer vision applications may benefit from understanding where
    % humans focus given a scene. Other than cognitively understanding the way
    % human perceive images and scenes, finding salient regions and objects in
    % the images helps various tasks such as speeding up object detection,
    % object recognition, object tracking and content-aware image editing.
    %
    % About the saliency, there is a rich literature but the development is
    % very fragmented. The principal purpose of this API is to give a unique
    % interface, a unique framework for use and plug sever saliency
    % algorithms, also with very different nature and methodology, but they
    % share the same purpose, organizing algorithms into three main
    % categories:
    %
    % * Static Saliency
    % * Motion Saliency
    % * Objectness
    %
    % Saliency UML diagram:
    %
    % ![image](https://docs.opencv.org/3.3.1/saliency.png)
    %
    % To see how API works, try tracker demo: `computeSaliency_demo.m`.
    %
    % ## Static Saliency algorithms
    %
    % Algorithms belonging to this category, exploit different image features
    % that allow to detect salient objects in a non dynamic scenarios.
    %
    % ### Fine Grained Saliency
    %
    % This method calculates saliency based on center-surround differences.
    % High resolution saliency maps are generated in real time by using
    % integral images.
    %
    % ## References
    % [FGS]:
    % > Sebastian Montabone and Alvaro Soto. "Human detection using a mobile
    % > platform and novel features derived from a visual saliency mechanism".
    % > In Image and Vision Computing, Vol. 28 Issue 3, pages 391-402.
    % > Elsevier, 2010.
    %
    % See also: cv.StaticSaliencySpectralResidual
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    methods
        function this = StaticSaliencyFineGrained()
            %STATICSALIENCYFINEGRAINED  Constructor, creates a specialized saliency algorithm of this type
            %
            %     obj = cv.StaticSaliencyFineGrained()
            %
            % See also: cv.StaticSaliencyFineGrained.computeSaliency
            %
            this.id = StaticSaliencyFineGrained_(0, 'new');
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.StaticSaliencyFineGrained
            %
            if isempty(this.id), return; end
            StaticSaliencyFineGrained_(this.id, 'delete');
        end
    end

    %% Algorithm
    methods
        function clear(this)
            %CLEAR  Clears the algorithm state
            %
            %     obj.clear()
            %
            % See also: cv.StaticSaliencyFineGrained.empty,
            %  cv.StaticSaliencyFineGrained.load
            %
            StaticSaliencyFineGrained_(this.id, 'clear');
        end

        function b = empty(this)
            %EMPTY  Checks if detector object is empty
            %
            %     b = obj.empty()
            %
            % ## Output
            % * __b__ Returns true if the detector object is empty (e.g in the
            %   very beginning or after unsuccessful read).
            %
            % See also: cv.StaticSaliencyFineGrained.clear,
            %  cv.StaticSaliencyFineGrained.load
            %
            b = StaticSaliencyFineGrained_(this.id, 'empty');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm parameters to a file
            %
            %     obj.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in the specified
            % XML or YAML file.
            %
            % See also: cv.StaticSaliencyFineGrained.load
            %
            StaticSaliencyFineGrained_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     obj.load(fname)
            %     obj.load(str, 'FromString',true)
            %     obj.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from the specified XML or
            % YAML file (either from disk or serialized string). The previous
            % algorithm state is discarded.
            %
            % See also: cv.StaticSaliencyFineGrained.save
            %
            StaticSaliencyFineGrained_(this.id, 'load', fname_or_str, varargin{:});
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = obj.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.StaticSaliencyFineGrained.save,
            %  cv.StaticSaliencyFineGrained.load
            %
            name = StaticSaliencyFineGrained_(this.id, 'getDefaultName');
        end
    end

    %% Saliency
    methods
        function saliencyMap = computeSaliency(this, img)
            %COMPUTESALIENCY  Compute the saliency
            %
            %     saliencyMap = obj.computeSaliency(img)
            %
            % ## Input
            % * __img__ The input image, 8-bit 1 or 3-channel (internally
            %   converted to grayscale).
            %
            % ## Output
            % * __saliencyMap__ The computed saliency map, `uint8` matrix.
            %
            % See also: cv.StaticSaliencyFineGrained.computeBinaryMap
            %
            saliencyMap = StaticSaliencyFineGrained_(this.id, 'computeSaliency', img);
        end
    end

    %% StaticSaliency
    methods
        function binaryMap = computeBinaryMap(this, saliencyMap)
            %COMPUTEBINARYMAP  This function perform a binary map of given saliency map
            %
            %     binaryMap = obj.computeBinaryMap(saliencyMap)
            %
            % ## Input
            % * __saliencyMap__ the saliency map obtained through one of the
            %   specialized algorithms, `single` matrix.
            %
            % ## Output
            % * __binaryMap__ the binary map, `uint8` matrix with either 0 or
            %   255 values.
            %
            % This is obtained in this way:
            % In a first step, to improve the definition of interest areas and
            % facilitate identification of targets, a segmentation by
            % clustering is performed, using *K-means algorithm*. Then, to
            % gain a binary representation of clustered saliency map, since
            % values of the map can vary according to the characteristics of
            % frame under analysis, it is not convenient to use a fixed
            % threshold. So, *Otsu's algorithm* is used, which assumes that
            % the image to be thresholded contains two classes of pixels or
            % bi-modal histograms (e.g. foreground and back-ground pixels);
            % later on, the algorithm calculates the optimal threshold
            % separating those two classes, so that their intra-class variance
            % is minimal.
            %
            % See also: cv.StaticSaliencyFineGrained.computeSaliency
            %
            if isa(saliencyMap, 'uint8')
                % saliency map expected to be floating point in [0,1] range
                saliencyMap = single(saliencyMap) / 255;
            end
            binaryMap = StaticSaliencyFineGrained_(this.id, 'computeBinaryMap', saliencyMap);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


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

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

% Last Modified by GUIDE v2.5 03-Feb-2006 11:59:36

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

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

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

% Choose default menu params
set(handles.modelMenu,    'String',varargin{1});
set(handles.parameterMenu,'String',varargin{2});
set(handles.fieldMenu,    'String',varargin{3});

% default command line output
data.model      = get(handles.modelMenu(1),'Value');
data.parameter  = get(handles.parameterMenu(1),'Value');
data.field      = get(handles.fieldMenu(1),'Value');

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

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes rmLoadInterface wait for user response (see UIRESUME)
uiwait(handles.figure1);
return;

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

% quit
delete(handles.figure1);
return;

% --- Executes on selection change in modelMenu.
function modelMenu_Callback(hObject, eventdata, handles)
%contents           =  get(hObject,'String');
handles.data.model = get(hObject,'Value');
guidata(hObject,handles);
return;
function modelMenu_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
return;

% --- Executes on selection change in parameterMenu.
function parameterMenu_Callback(hObject, eventdata, handles)
%contents           =  get(hObject,'String');
%handles.data.parameter = contents{get(hObject,'Value')};
handles.data.parameter = get(hObject,'Value');
% here we change the (default) field the parameter goes to
contents  =  get(hObject,'String');
switch contents{get(hObject,'Value')},
 case {'variance explained' 'max of all possible log10(p)s' 'coherence'},
  % this one usually goes to the coherence field
  handles.data.field = find(strcmp(get(handles.fieldMenu,'String'),'co'));
 case {'polar-angle'},
  % this one usually goes to the phase field
  handles.data.field = find(strcmp(get(handles.fieldMenu,'String'),'ph'));
 otherwise,
  % rest will typically go to the map field
  handles.data.field = find(strcmp(get(handles.fieldMenu,'String'),'map'));
end;
% update in gui
set(handles.fieldMenu,'value',handles.data.field);
guidata(hObject,handles);
return;
function parameterMenu_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
return;

% --- Executes on selection change in fieldMenu.
function fieldMenu_Callback(hObject, eventdata, handles)
%contents           =  get(hObject,'String');
%handles.data.field = contents{get(hObject,'Value')};
handles.data.field = get(hObject,'Value');
guidata(hObject,handles);
return;
function fieldMenu_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
return;

% --- Executes on button press in cancelbutton.
function cancelbutton_Callback(hObject, eventdata, handles)
% clear output
handles.data = [];
guidata(hObject,handles);
% now continue and finish
uiresume(handles.figure1);
return;

% --- Executes on button press in gobutton.
function gobutton_Callback(hObject, eventdata, handles)
% now continue and finish
uiresume(handles.figure1);
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [v,s] = CircularVariance(angles,dim)

%CircularVariance - Estimate circular variance and standard deviation.
%
%  USAGE
%
%    [v,s] = CircularVariance(angles,dim)
%
%    angles         angles in radians
%    dim            optional dimension along which the mean should be computed
%
%  OUTPUT
%
%    v              variance
%    s              standard deviation
%
%  SEE
%
%    See also circ_mean, CircularConfidenceIntervals, Concentration,
%    ConcentrationTest.

% Copyright (C) 2004-2011 by Michaël Zugaro
%
% 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.

if nargin < 2,
	dim = 1;
end

isradians(angles);

angles = exp(i*angles);
r_bar = abs(mean(angles,dim)); % is this a bug? should it be circ_mean?
v = 1 - r_bar;
s = sqrt(-2*log(r_bar));

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [node,face,elem]=meshanellip(c0,rr,tsize,maxvol)
%
% [node,face,elem]=meshanellip(c0,rr,opt)
%
% create the surface and tetrahedral mesh of an ellipsoid
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%   c0:  center coordinates (x0,y0,z0) of the ellipsoid
%   rr:  radii of an ellipsoid,
%        if rr is a scalar, this is a sphere with radius rr
%        if rr is a 1x3 or 3x1 vector, it specifies the ellipsoid radii [a,b,c]
%        if rr is a 1x5 or 5x1 vector, it specifies [a,b,c,theta,phi]
%           where theta and phi are the rotation angles along z and x
%           axes, respectively. Rotation is applied before translation.
%   tsize: maximum surface triangle size on the sphere
%   maxvol: maximu volume of the tetrahedral elements
%
% output:
%   node: node coordinates, 3 columns for x, y and z respectively
%   face: integer array with dimensions of NB x 3, each row represents
%         a surface mesh face element
%   elem: integer array with dimensions of NE x 4, each row represents
%         a tetrahedron; if ignored, only produces the surface
%
% example:
%   [node,face,elem]=meshanellip([10,10,-10],[30,20,10,pi/4,pi/4],0.5,0.4);
%   plotmesh(node,elem,'x>10');axis equal;
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

if(nargin<3)
    error('you must at least provide c0, rr and tsize, see help for details');
end

rr=rr(:)';
if(length(rr)==1)
    rr=[rr,rr,rr];
elseif(length(rr)==3)
    % do nothing
elseif(length(rr)~=5)
    error('the number of elements for rr is not correct. see help for details');
end

rmax=min(rr(1:3));
if(nargout==3)
    if(nargin==3)
	maxvol=tsize*tsize*tsize;
    end
    [node,face,elem]=meshunitsphere(tsize/rmax,maxvol/(rmax*rmax*rmax));
else
    [node,face]=meshunitsphere(tsize/rmax);
end

node=node*diag(rr(1:3));
if(length(rr)==5)
   theta=rr(4);
   phi=rr(5);
   Rz=[cos(theta) sin(theta) 0;-sin(theta) cos(theta) 0;0 0 1];
   Rx=[1 0 0;0 cos(phi) sin(phi);0 -sin(phi) cos(phi)];
   node=(Rz*Rx*(node'))';
end
node=node+repmat(c0(:)',size(node,1),1);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function hActive = getActiveSegment
%Get the handle to the 'active' segment
%LM:  1 April 02, JOD.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

%Get handles to all the existing plotted segments
hSegsV = findobj('tag','newSegment');

hActive = [];

%Get active segment data:
if ~isempty(hSegsV)
  i = 1;
  done = 0;
  while ~done & i <= length(hSegsV)
    segS = get(hSegsV(i),'userdata');
    if segS.active == 1
      hActive = hSegsV(i);
      xV = segS.GUI_xV;
      yV = segS.GUI_yV;
      hPlot = segS.hPlot;
      done = 1;
    end
    i = i + 1;
  end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function wathen_test04 ( )

%*****************************************************************************80
%
%% WATHEN_TEST04 assembles, factor and solve using WATHEN_DAVIS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 June 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'WATHEN_TEST04\n' );
  fprintf ( 1, '  Assemble, factor and solve a Wathen system\n' );
  fprintf ( 1, '  defined by WATHEN_DAVIS.\n' );
  fprintf ( 1, '\n' );

  nx = 4;
  ny = 4;
  fprintf ( 1, '  Elements in X direction NX = %d\n', nx );
  fprintf ( 1, '  Elements in Y direction NY = %d\n', ny );
  fprintf ( 1, '  Number of elements = %d\n', nx * ny );
%
%  Compute the number of unknowns.
%
  n = wathen_order ( nx, ny );
  fprintf ( 1, '  Number of nodes N = %d\n', n );
%
%  Set up a random solution X1.
%
  seed = 123456789;
  [ x1, seed ] = r8vec_uniform_01 ( n, seed );
%
%  Compute the matrix.
%
  seed = 123456789;
  [ a, seed ] = wathen_davis ( nx, ny, n, seed );
%
%  Compute the corresponding right hand side B.
%
  b = a * x1;
%
%  Solve the linear system.
%
  x2 = a \ b;
%
%  Compute the maximum solution error.
%
  e = max ( abs ( x1 - x2 ) );
  fprintf ( 1, '  Maximum solution error is %g\n', e );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ucmThresh, meanNSP] = computeSPThreshold(imSet, numSP)

	imList = getImageSet(imSet);
	high = 1;
	low = 0.005;
	thresh = 0.005;

	while(high > low+thresh)
		mid = (low+high)/2;
		[meanNSP, nSP] = helper(imList, mid);
		if(meanNSP < numSP)
			high = mid;
		else
			low = mid;
		end
		fprintf('Between %0.3f, and %0.3f\n', low, high);
	end

	ucmThresh = mid;
end

function [meanNSP, nSP] = helper(imList, ucmThresh)
	parfor i = 1:length(imList),
		[~, ~, nSP(i), ~] = getSuperpixels(imList{i}, ucmThresh);
	end
	meanNSP = mean(nSP);
end

author: s-gupta

repo: rgbd

save_path: github-repos/MATLAB/s-gupta-rgbd


text:
xm    = 5;
Px    = 9;
kappa = 2;

[Xi W]      = SigmaPoints(xm, Px, kappa);
[xAvg xCov] = UT(Xi, W)

author: philbooks

repo: Kalman-Filter-for-Beginners


text:
%TEST8 test cs_cholsol, cs_lusol
%
% Example:
%   test8
% See also: testall

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

index = UFget ;
[ignore f] = sort (max (index.nrows, index.ncols)) ;
f = f (1:100) ;

% f = f(1)

for i = f
    Prob = UFget (i) ;
    disp (Prob) ;
    A = Prob.A ;
    [m n] = size (A) ;
    if (~isreal (A) | m ~= n)                                               %#ok
        continue
    end

    spd = 0 ;
    if (m == n)
        if (nnz (A-A') == 0)
            try
                p = amd (A) ;
            catch
                p = symamd (A) ;
            end
            [R,p] = chol (A (p,p)) ;
            spd = (p == 0) ;
        end
    end

    if (spd)
        C = A ;
    else
        C = A*A' + n*speye (n) ;
        try
            p = amd (C) ;
        catch
            p = symamd (C) ;
        end
        try
            R = chol (C (p,p)) ;
        catch
            continue
        end
    end

    b = rand (n,1) ;

    x1 = C\b ;
    x2 = cs_cholsol (C,b) ;
    r1 = norm (C*x1-b,1) / norm (C,1) ;
    r2 = norm (C*x2-b,1) / norm (C,1) ;
    err = abs (r1-r2) ;
    fprintf ('err %g\n', err) ;
    if (err > 1e-10)
        error ('!') ;
    end

    x2 = cs_lusol (C,b, 1, 0.001) ;
    r2 = norm (C*x2-b,1) / norm (C,1) ;
    err = abs (r1-r2) ;
    fprintf ('err %g (lu with amd(A+A'')\n', err) ;
    if (err > 1e-10)
        error ('!') ;
    end

    if (m ~= n)
        continue ;
    end

    x1 = A\b ;
    r1 = norm (A*x1-b,1) / norm (A,1) ;
    if (r1 < 1e-6)
        x2 = cs_lusol (A,b) ;
        r2 = norm (A*x2-b,1) / norm (A,1) ;
        fprintf ('lu resid %g %g\n', r1, r2) ;
    end
end


text:
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
%function ambifunt
%AMBIFUNT Unit test for the function ambifunb.

%	O. Lemoine - December 1995.
%	F. Auger - February 1996.

N=128;

% Ambiguity function of a pulse
sig=((1:N)'==40);
amb=ambifunb(sig);
[ik,jk]=find(amb~=0.0);
if any(jk~=N/2)|any(ik'-(1:N)),
 error('ambifunb test 1 failed');
end;

% Ambiguity function of a sine wave
sig=fmconst(N,0.2);
amb=ambifunb(sig);
[a b]=max(amb);
if any(b-odd(N/2)*ones(1,N-1)),
 error('ambifunb test 2 failed');
end;

% Energy
sig=noisecg(N);
amb=ambifunb(sig);
if abs(abs(amb(odd(N/2),N/2))-norm(sig)^2)>sqrt(eps),
 error('ambifunb test 3 failed');
end;

% Link with the Wigner-Ville distribution
sig=fmlin(N);
amb=ambifunb(sig);
amb=amb([(N+rem(N,2))/2+1:N 1:(N+rem(N,2))/2],:);
ambi=ifft(amb).';
tdr=zeros(N); 		% Time-delay representation
tdr(1:N/2,:)=ambi(N/2:N-1,:);
tdr(N:-1:N/2+2,:)=ambi(N/2-1:-1:1,:);
wvd1=real(fft(tdr));
wvd2=tfrwv(sig);
errors=max(max(abs(wvd1-wvd2)));
if errors>sqrt(eps),
 error('ambifunb test 4 failed');
end;

N=111;

% Ambiguity function of a pulse
sig=((1:N)'==40);
amb=ambifunb(sig);
[ik,jk]=find(amb~=0.0);
if any(jk~=(N+1)/2)|any(ik'-(1:N)),
 error('ambifunb test 5 failed');
end;

% Ambiguity function of a sine wave
sig=fmconst(N,0.2);
amb=ambifunb(sig);
[a b]=max(amb(:,2:N-1));
if any(b-((N+1)/2)*ones(1,N-2)),
 error('ambifunb test 6 failed');
end;

% Energy
sig=noisecg(N);
amb=ambifunb(sig);
if abs(abs(amb((N+1)/2,(N+1)/2))-norm(sig)^2)>sqrt(eps),
 error('ambifunb test 7 failed');
end;

% Link with the Wigner-Ville distribution
sig=fmlin(N);
amb=ambifunb(sig);
amb=amb([(N+1)/2:N 1:(N-1)/2],:);
ambi=ifft(amb).';
tdr=zeros(N); 		% Time-delay representation
tdr(1:(N+1)/2,:)      = ambi((N+1)/2:N,:);
tdr(N:-1:(N+1)/2+1,:) = ambi((N-1)/2:-1:1,:);
wvd1=real(fft(tdr));
wvd2=tfrwv(sig,1:N,N);
errors=max(max(abs(wvd1-wvd2)));
if errors>sqrt(eps),
 error('ambifunb test 8 failed');
end;

author: HeLiangHIT


text:
function coords = removeCoords(coords1,coords2)
%
% coords = removeCoords(coords1,coords2,dims)
%
% Removes coords1 from coords2, used for example to remove
% coordinates from an ROI.
%
% coords, coords1, and coords2: 3xN arrays of (y,x,z) coordinates
% dims is size of volume
%
% djh, 7/98
% djh, 2/2001, dumped coords2Indices & replaced with union(coords1',coords2','rows')
% MMS&AW 1/2005 changed the usage of the setdiff (line 16 and 17) to avoid resorting of the coord order

if ~isempty(coords1)
    [waste,indx]= setdiff(coords2',coords1','rows');
    coords=coords2(:,sort(indx));
else
    coords=coords2;
end
return


text:
function w = ymdf_to_weekday_english2 ( y, m, d, f )

%*****************************************************************************80
%
%% YMDF_TO_WEEKDAY_ENGLISH2 returns the weekday of an English YMDF date.
%
%  Discussion:
%
%    Lewis Carroll's method is used.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Gary Meisters,
%    Lewis Carroll's Day-of-the-Week Algorithm,
%    Math Horizons, November 2002, pages 24-25.
%
%    Lewis Carroll (Charles Dodgson),
%    To Find the Day of the Week for Any Given Date,
%    Nature, 31 March 1887.
%
%    John Conway,
%    Tomorrow is the Day After Doomsday,
%    Eureka, 1973.
%
%  Parameters:
%
%    Input, integer Y, M, D, real F, the YMDF date.
%
%    Output, integer W, is the week day number of the date, with
%    1 for Sunday, through 7 for Saturday.
%
  m_table = [ 0, 3, 3, 6, 1, 4, 6, 2, 5, 0, 3, 12 ];
%
%  Dates of 2 September 1752 and earlier fall under one scheme.
%
  y1 = 1752;
  m1 = 9;
  d1 = 3;
  f1 = 0.0;

  cmp1 = ymdf_compare ( y, m, d, f, y1, m1, d1, f1 );
%
%  Dates of 14 September 1752 and later fall under another.
%
  y1 = 1752;
  m1 = 9;
  d1 = 13;
  f1 = 0.0;

  cmp2 = ymdf_compare ( y, m, d, f, y1, m1, d1, f1 );

  if ( cmp1 == '<' )

    cval = 18 - floor ( y / 100 );

  elseif ( cmp2 == '>' )

    cval = 2 * ( 3 - mod ( floor ( y / 100 ), 4 ) );

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'YMDF_TO_WEEKDAY_ENGLISH2 - Fatal error!\n' );
    fprintf ( 1, '  The given date is illegal.\n' );
    error ( 'YMDF_TO_ENGLISH2 - Fatal error!' );

  end

  cc = floor ( y / 100 );
  yy = y - cc * 100;
  a = floor ( yy / 12 );
  b = yy - 12 * a;
  yval = a + b + floor ( b / 4 );

  mval = m_table(m);

  dval = d;

  if ( ( m == 1 || m == 2 ) && year_is_leap_english ( y ) )
    mval = mval + 6;
  end

  w = i4_modp ( cval + yval + mval + dval, 7 ) + 1;

  return
end


text:
function disc = extractBlock(disc, j, k)
%EXTRACTBLOCK   Extract the j-k block from a discretization.
%   DISCJK = EXTRACTBLOCK(DISC, J, K) extracts information relating to the
%   J-K block of the dscretization DISC and returns a new discretization
%   DISCJK. DISCJK.dimension will be modified by the amount specified in
%   DISC.dimAdjust(j,k).

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

% Extract the j-k block:
if ( iscell(disc.source.blocks) )
    disc.source.blocks = disc.source.blocks{j,k};
end

% Extract the dimension adjustment for this block and adjust the dimension:
if ( numel(disc.dimAdjust) > 1 )
    disc.dimAdjust = disc.dimAdjust(j,k);
end
disc.dimension = disc.dimension + disc.dimAdjust;

% Set the dimension adjustment to zero (as dimension was adjusted above)
disc.dimAdjust = 0;

% Extract the projection Order:
if ( numel(disc.projOrder) > 1 )
    disc.projOrder = disc.projOrder(j);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [params, names] = fileKernExtractParam(kern)

% FILEKERNEXTRACTPARAM Extract parameters from the FILE kernel structure.
% FORMAT
% DESC Extract parameters from the stored file kernel structure
% into a vector of parameters for optimisation.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the kernel. If
% the field 'transforms' is not empty in the kernel matrix, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
%
% FORMAT
% DESC Extract parameters and parameter names from the stored file
% kernel structure.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the kernel. If
% the field 'transforms' is not empty in the kernel matrix, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
% RETURN names : cell array of strings containing parameter names.
%
% SEEALSO fileKernParamInit, fileKernExpandParam, kernExtractParam, scg, conjgrad
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006
%
% KERN

params = [kern.variance];
if nargout > 1
  names ={'variance'};
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [X,info] = spqr_solve (A,B,opts)                                   %#ok
%SPQR_SOLVE solves a linear system or least squares problem via QR factorization.
% A is a sparse m-by-n matrix A and B is a sparse or full m-by-k matrix.
% If m == n, x = spqr_solve(A,B) solves A*X=B.  If m < n, a basic solution is
% found for the underdetermined system A*X=B.  If m > n, the least squares
% solution is found.  An optional third input argument specifies non-default
% options (see "help spqr" for details).  Only opts.tol, opts.ordering, and
% opts.solution are used; the others are implicitly set to opts.econ = 0,
% opts.Q = 'Householder', and opts.permutation = 'vector'.
%
% opts.solution:  'basic' (default) or 'min2norm'.  To obtain a minimum 2-norm
% solution to an undetermined system (m < n), use 'min2norm'.  For m >= n,
% these two options find the same solution.
%
% An optional second output provides statistics on the solution.
%
% Example:
%   x = spqr_solve (A,B)
%   [x,info] = spqr_solve (A,B)
%   x = spqr_solve (A,B,opts)
%   x = spqr_solve (A,B, struct ('solution','min2norm')) ;
%
% See also SPQR, SPQR_QMULT, QR, MLDIVIDE.

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

type spqr_solve
error ('spqr_solve mexFunction not found') ;


text:
function cvx_pop( p, do_warn )
evalin( 'caller', sprintf( 'cvx_pop(%d,%d,%d)', p.index_, p.id_, do_warn ) );

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

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


text:
% Function to perform Principle Component Analysis over a set of training
% vectors passed as a concatenated matrix.
%
% Usage:- [V,D,M] = pca(X,n)
%         [V,D] = pca(X,aM,n)
%
% where:-
%        <input>
%        X = concatenated set of column vectors
%        aM = assume that the mean is aM
%        n = number of principal components to extract (optional)
%        <output>
%        V = ensemble of column eigen-vectors
%        D = vector of eigen-values
%        M = mean of X (optional)
%
%
% by Simon Lucey 2006

function [v,d,varargout] = pca(X,varargin)

% Check arguments
neigs = min(size(X)) - 1; % Set the number of eigen-vectors
aM = mean(X,2); % Otherwise set to data mean
if nargin > 1
  % Check if second argument is a vector
  if size(varargin{1},1) == size(X,1)
    aM = varargin{1}; % Set to specified mean
    if nargin == 3
      neigs = varargin{2}; % Set the number of eigen-vectors
    end
  else
    neigs = varargin{1}; % Set the number of eigen-vectors
  end
end

N = size(X,2); % Get number of samples in training set
D = size(X,1); % Get dimension of training set

% Remove the mean from the data
nX = X - repmat(aM,[1,size(X,2)]);

if N > D
  [v,d] = pca_NgtD(nX,neigs);
else
  [v,d] = pca_NltD(nX,neigs);
end

% Set outputs
if nargout == 3
  varargout(1) = {aM};
end

% Function to do PCA if N is greater than (gt) D
function [v,d] = pca_NgtD(nX,n)

% Get the number of samples
N = size(nX,2);

% Calculate correlation matrix
R = (nX*nX')/N;

% Now calculate the eigen vectors of this matrix
[v,d] = eig(R);
d = diag(d); % Get main diagonal of eigen values
[dummy,i] = sort(d); % Sort values in ascending order
i = flipud(i); % Change to be in descending order

% Check if n is less than one
if n < 1
    n = num_neigs(d(i),n);
end

v = v(:,i(1:n)); % Set the eigen-vectors
d = d(i(1:n)); % Set the eigen-values

% Function to do PCA if N is less than (lt) D
function [v,d] = pca_NltD(nX,n)

N = size(nX,2);
D = size(nX,1);

% Now lets get the autocorrelation matrix for mxm not nxn matrix
S = (nX'*nX)/D;

% Now calculate eigenvectors and eigen values
[v,d] = eig(S);

% Set transpose of vectors
v = nX*v; % New eigen vectors as per pg 40 of Fukunaga
v = v*inv(sqrt(D*d)); % Re-weight the vectors

d = D*diag(d)/N; % Get main diagonal of eigen values
[dummy,i] = sort(d); % Sort values in ascending order
i = flipud(i); % Change to be in descending order

% Check if n is less than one
if n < 1
    n = num_neigs(d(i),n);
end

v = v(:,i(1:n)); % Set the eigen-vectors
d = d(i(1:n)); % Set the eigen-values

function n = num_neigs(d,p)

k = cumsum(d)./sum(d); % Get the cumulative sum
j = find(k >= p); % Find the index
n = j(1);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
classdef ParadigmCSP < ParadigmDataflowSimplified
    % Common Spatial Pattern(s) algorithm.
    %
    % The CSP paradigm is based on the design of the Berlin Brain-Computer Interface (BBCI) [1], more
    % comprehensively described in [2], which is mainly controlled by (sensori-)motor imagery. The
    % features exploited by this paradigm in its original form are Event-Related Synchronization and
    % Desynchronization [3] localized in the (sensori-)motor cortex, but the paradigm is not restricted
    % to these applications. CSP was originally introduced in [5] and first applied to EEG in [6].
    %
    % Due to its simplicity, speed and relative robustness, CSP is the bread-and-butter paradigm for
    % oscillatory processes, and if nothing else, can be used to get a quick estimate of whether the
    % data contains information of interest or not. Like para_bandpower, CSP uses log-variance features
    % over a single non-adapted frequency range (which may have multiple peaks), and neither temporal
    % structure (variations) in the signal is captured, nor are interactions between frequency bands.
    % The major strength of the paradigm is its adaptive spatial filter, which is computed using the CSP
    % algorithm.
    %
    % The paradigm is implemented as a standard sequence of signal (pre-)processing (spatial/spectral
    % filtering), feature extraction, and machine learning. The first preprocessing step is frequency
    % filtering, followed by an adaptively learned spatial filter (which is the defining propery of the
    % paradigm), followed by log-variance feature extraction and finally a (usually simple) machine
    % learning step applied to the log-variance features. The spatial filtering projects the channels of
    % the original signal down to a small set of (usually 4-6) surrogate channels, where the (linear)
    % mapping is optimized such that the variance in these channels is maximally informative w.r.t. to
    % the prediction task. The CSP filters can be obtained from the per-class signal covariance matrices
    % by solving a generalized eigenvalue problem (of the form [V,D]=eig(Cov1,Cov1+Cov2)). CSP can also
    % be applied to independent components to rate their importance or for better artifact robustness. A
    % wide range of classifiers can be used with CSP features, the most commonly used one being LDA.
    % There exists a large corpus of CSP variants and extensions, mostly to give better control over
    % spectral filtering, including multiband CSP (para_multiband_csp), Spectrally Weighted CSP
    % (para_speccsp), Invariant CSP, Common Spatio-Spectral Patterns (CSSP), Common Sparse Spectral
    % Spatial Pattern (CSSSP), Regularized CSP, and several others. A more advanced (but also
    % computationally more costly) paradigm than CSP is the Dual-Augmented Lagrange Paradigm
    % (para_dal/para_dal_hf). The length of the data epoch and the choice of a frequency band
    % (defaulting to motor imagery time scales and frequency ranges) are the parameters that are most
    % commonly tuned to the task, both of which can also be found via a small parameter search.
    %
    % Some application areas include detection of major brain rhythm modulations (e.g. alpha, beta), for
    % example related to relaxation/stress, aspects of workload, sensori-motor imagery, visual
    % processing vs. idling and other idle-rhythm-related questions, or emotion recognition. See also
    % [4].
    %
    % Examples:
    %   After an approach has been defined as in one of the following examples, a predictive model can be obtained
    %   (given a calibration data set and a specification of target markers) using bci_train:
    %   [loss,model,stats] = bci_train('Data',io_loadset('calibration_rec.eeg'),'Approach',myapproach','TargetMarkers',{'mymarker1','mymarker2'});
    %
    %   % define a basic CSP approach, using the defaults (7-30 Hz bandpass, shrinkage LDA classifier, 3 pattern pairs)
    %   myapproach = 'CSP';
    %
    %   % use an FIR filter restricted to the alpha band
    %   myapproach = {'CSP' 'SignalProcessing',{'FIRFilter',[6 8 14 15]}};
    %
    %   % use an IIR filter instead of the default FIR
    %   myapproach = {'CSP' 'SignalProcessing',{'FIRFilter','off', 'IIRFilter',[5 7 14 17]}};
    %
    %   % also restrict the model to a stationary subspace
    %   myapproach = {'CSP' 'SignalProcessing',{'FIRFilter','off', 'IIRFilter',[5 7 14 17], 'StationarySubspace',{'StationaryDim',-0.1,'Operation','keep_stationary'}}};
    %
    %   % use a sharp FFT band-pass filter
    %   myapproach = {'CSP' 'SignalProcessing',{'FIRFilter','off','SpectralSelection',[7 15]}};
    %
    %   % use a simple logistic regression classifier (variational Bayes) instead of the LDA
    %   myapproach = {'CSP' 'Prediction',{'MachineLearning',{'Learner','logreg'}}};
    %
    %   % use a simple logistic regression classifier (sparse variational Bayes)
    %   myapproach = {'CSP' 'Prediction',{'MachineLearning',{'Learner',{'logreg','variant','vb-ard'}}}};
    %
    %   % use the sparse logistic regression classifier but applied to a larger set of patterns
    %   myapproach = {'CSP' 'Prediction',{'MachineLearning',{'Learner',{'logreg','variant','vb-ard'}},'FeatureExtraction',{'PatternPairs',6}}};
    %
    %   % using quadratic discriminant analysis
    %   myapproach = {'CSP' 'Prediction',{'MachineLearning',{'Learner','qda'}}};
    %
    %   % using Gaussian mixture models (variational Bayesian Dirichlet process prior)
    %   myapproach = {'CSP' 'Prediction',{'MachineLearning',{'Learner','gmm'}}};
    %
    %   % using relevance vector machines (here with a fixed kernel scale for speed)
    %   myapproach = {'CSP' 'Prediction',{'MachineLearning',{'Learner',{'rvm','kernel','rbf','gamma',1}}}};
    %
    %   % use l1-regularized logreg (which involves a parameter search over the regularization parameter); takes about 2 minutes
    %   myapproach = {'CSP' 'Prediction',{'MachineLearning',{'Learner',{'logreg','variant','l1','lambda',search(2.^(-6:1.5:10))}},'FeatureExtraction',{'PatternPairs',6}}};
    %
    %   % use support vector machines (using the SVMlight package); note: also requires a parameter search (takes about 2 minutes)
    %   myapproach = {'CSP' 'Prediction',{'MachineLearning',{'Learner',{'svmlight','cost',search(2.^(-5:2.5:15)),'gamma',1}}}};
    %
    %   % optimize the size of the stationary subspace to retain; note: cross-validation takes place on continuous data (takes approx 2.5 minutes)
    %   myapproach = {'CSP' 'SignalProcessing',{'FIRFilter','off', 'IIRFilter',[5 7 14 17], 'StationarySubspace',{'StationaryDim',search(-0.5:0.1:-0.1),'Operation','keep_stationary'}}};
    %
    %   % optimize the location of the frequency band manually (note: 2-dimensional parameter space -- takes approx. 7 minutes)
    %   myapproach = {'CSP' 'SignalProcessing',{'FIRFilter','off', 'SpectralSelection',[search(7:10) search(14:2:30)]}};
    %
    % References:
    %  [1] Blankertz, B., Dornhege, G., Krauledat, M., M�ller, K., and Curio,G.
    %      "The non-invasive Berlin Brain-Computer interface: Fast acquisition of effective performance in untrained subjects."
    %      NeuroImage 37, 2 (Aug. 2007), 539?550.
    %  [2] Benjamin Blankertz, Ryota Tomioka, Steven Lemm, Motoaki Kawanabe, and Klaus-Robert Mueller.
    %      "Optimizing spatial filters for robust EEG single-trial analysis."
    %      IEEE Signal Process Mag, 25(1):41-56, January 2008
    %  [3] Pfurtscheller, G., and da Silva, L. "Event-related EEG/MEG synchronizaion and desynchronization: basic principles."
    %      Clin Neurophysiol 110 (1999), 1842-1857.
    %  [4] Buzsaki, G., "Rhythms of the brain"
    %      Oxford University Press US, 2006
    %  [5] Fukunaga K., "Introduction to Statistical Pattern Recognition"
    %      Academic Press, Computer Science and Scientific Computing Series, 1990
    %  [6] Ramoser, H., Gerking, M., and Pfurtscheller, G. "Optimal spatial filtering of single trial EEG during imagined hand movement."
    %      IEEE Trans. Rehab. Eng 8 (2000), 446, 441.
    %
    % Name:
    %   Common Spatial Patterns
    %
    %                            Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
    %                            2010-04-29

    methods

        function defaults = preprocessing_defaults(self)
            % define the default pre-processing parameters of this paradigm
            defaults = {'FIRFilter',{'Frequencies',[6 8 28 32],'Type','minimum-phase'}, 'EpochExtraction',[0.5 3.5], 'Resampling',100};
        end

        function model = feature_adapt(self,varargin)
            % adapt a feature representation using the CSP algorithm
            arg_define(varargin, ...
                arg_norep('signal'), ...
                arg({'patterns','PatternPairs'},3,uint32([1 1 64 10000]),'Number of CSP patterns (times two).','cat','Feature Extraction'),...
                arg({'dologtransform','LogTransform','logtransform'},true,[],'Perform log-transform. This is almost always the right thing to do.'),...
                arg({'shrinkage_cov','ShrinkageCovariance','ShrinkageCov'},false,[],'Shrinkage covariance estimator. Whether to use shrinkage to estimate the covariance matrices.'));

            if signal.nbchan < patterns
                error('CSP requires at least as many channels as you request output patterns. Please reduce the number of pattern pairs.'); end
            for k=1:2
                trials{k} = exp_eval(set_picktrials(signal,'rank',k));
                if shrinkage_cov
                    covar{k} = hlp_diskcache('featuremodels',@cov_shrink,reshape(trials{k}.data,size(trials{k}.data,1),[])');
                else
                    covar{k} = cov(reshape(trials{k}.data,size(trials{k}.data,1),[])');
                end
                covar{k}(~isfinite(covar{k})) = 0;
            end
            [V,D] = eig(covar{1},covar{1}+covar{2}); %#ok<NASGU>
            model.filters = V(:,[1:patterns end-patterns+1:end]);
            P = inv(V);
            model.patterns = P([1:patterns end-patterns+1:end],:);
            model.cov = cov(signal.data(:,:)');
            model.chanlocs = signal.chanlocs;
            model.logtransform = dologtransform;
        end

        function features = feature_extract(self,signal,featuremodel)
            % extract log-variance features according to CSP
            features = zeros(size(signal.data,3),size(featuremodel.filters,2));
            for t=1:size(signal.data,3)
                features(t,:) = var(signal.data(:,:,t)'*featuremodel.filters); end
            if featuremodel.logtransform
                features = log(features); end
        end

        function visualize_model(self,varargin) %#ok<*INUSD>
            args = arg_define([0 3],varargin, ...
                arg_norep({'myparent','Parent'},[],[],'Parent figure.'), ...
                arg_norep({'featuremodel','FeatureModel'},[],[],'Feature model. This is the part of the model that describes the feature extraction.'), ...
                arg_norep({'predictivemodel','PredictiveModel'},[],[],'Predictive model. This is the part of the model that describes the predictive mapping.'), ...
                arg({'patterns','PlotPatterns'},true,[],'Plot patterns instead of filters. Whether to plot spatial patterns (forward projections) rather than spatial filters.'), ...
                arg({'paper','PaperFigure'},false,[],'Use paper-style font sizes. Whether to generate a plot with font sizes etc. adjusted for paper.'), ...
                arg_nogui({'nosedir_override','NoseDirectionOverride'},'',{'','+X','+Y','-X','-Y'},'Override nose direction.'));
            arg_toworkspace(args);

            % determine nose direction for EEGLAB graphics
            try
                nosedir = args.fmodel.signal.info.chaninfo.nosedir;
            catch
                disp_once('Nose direction for plotting not store in model; assuming +X');
                nosedir = '+X';
            end
            if ~isempty(nosedir_override)
                nosedir = nosedir_override; end
            % number of pairs, and index of pattern per subplot
            np = size(featuremodel.patterns,1)/2;
            idx = [1:np 2*np:-1:np+1];
            % for each CSP pattern...
            for p=1:np*2
                subplot(2,np,p,'Parent',myparent);
                if args.patterns
                    plotdata = featuremodel.patterns(idx(p),:);
                else
                    plotdata = featuremodel.filters(:,idx(p));
                end
                topoplot(plotdata,featuremodel.chanlocs,'nosedir',nosedir);
                t = title(['CSP Pattern ' num2str(idx(p))]);
                if args.paper
                    set(t,'FontUnits','normalized');
                    set(t,'FontSize',0.1);
                end
            end
        end

        function layout = dialog_layout_defaults(self)
            % define the default configuration dialog layout
            layout = {'SignalProcessing.Resampling.SamplingRate', 'SignalProcessing.FIRFilter.Frequencies', ...
                'SignalProcessing.FIRFilter.Type', 'SignalProcessing.EpochExtraction', '', ...
                'Prediction.FeatureExtraction.PatternPairs', '', 'Prediction.MachineLearning.Learner'};
        end

        function tf = needs_voting(self)
            % standard CSP requires voting to handle more than 2 classes
            tf = true;
        end

    end
end

repo: REST


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

%% Syntax
% |[h]=cFigure(figStruct);|

%% Description
% Creates a custom figure using the input structure figStruct. The cFigure
% function provides easy control of background color, the color
% definitions, the figure window size (e.g. near maximal), and enables
% figure property cloning. It also allows users to create hidden figures
% which can be made visible for instance using the mfv command.

%% Examples

%%
clear; close all; clc;

%%
% Plot settings
fontSize=15;
faceAlpha=0.75;
edgeColor=0.3*ones(1,3);
edgeWidth=1.5;

%% Basic use of |cFigure| and specification of custom figure properties

%Modifying standard figure properties
figStruct.Name='GIBBON'; %Figure name
figStruct.Color='k'; %Figure background color

%Custom figure properties
figStruct.ColorDef='black'; %Setting colordefinitions to black
figStruct.ScreenOffset=0; %Setting spacing of figure with respect to screen edges

[H(1)]=cFigure(figStruct);

%Plot something
surf(peaks(25)); axis equal; axis tight;

%% Creating hidden figures which are made visible using |mfv|
% The below example demonstrates how to create 1 or more figures that
% remain hidden until the |mfv| command is used. This way the figures do
% not "steal focus". This can be an issue sometimes when multiple figures
% are used over time (and with some operational systems).

%Create property structure
figStructHid.Visible='off';

%Open a hidden figure
[H(1)]=cFigure(figStructHid);

%Plot something in it
hs=surf(peaks(25));
axis equal; axis tight; colormap jet;
drawnow;

%Open another hidden figure
[H(2)]=cFigure(figStructHid);

%Plot something in it
hs=surf(peaks(25));
axis equal; axis tight; colormap gray;
drawnow;

pause(2); %Add pause to demonstrate lack of visible figures

%%
% At this point all figures are hidden so have not stolen focus. The
% figures can be made visible using the |mfv| command:

mfv(H); %Now make all in H visible

%% Calling hidden figures using |scf|

%Create property structure
figStructHid.Visible='off';
figStruct.ScreenOffset=300;

%Open a hidden figure without plotting anything in it
[H(1)]=cFigure(figStructHid);

%Open another hidden figure
[H(2)]=cFigure(figStructHid);

%Plot something in it
hs=surf(peaks(25));
axis equal; axis tight; colormap autumn;
drawnow;

%Calling first figure again without making it visible
scf(H(1));

%Now plot something in it
hs=surf(peaks(25)); axis equal; axis tight; colormap winter;
drawnow;

%Now show the figures
mfv; %If this is used instead of mfv(H) then all hidden figures are shown

%% Using the preloaded view control widget |vcw|
% See also: |vcw|;
%%
% Location of the View Control Widget (|vcw|) activation button in the menubar:
%
% <<vcw_icon_menubar.png>>

%%
% Press v key to activate widget, press i to show help information
%
% <<vcw_help_show.png>>

%% Using the preloaded Export Figure Widget |efw|
% See also: |efw|;
%%
% Location of the Export Figure Widget (|efw|) button in the menubar:
%
% <<efw_icon_menubar.png>>

%%
%
% <<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 [f,J,Q] = spm_fx_mmc(x,u,P,M)
% state equations for a neural mass model of motor cortex
% Bhatt et al. 2016 Neuroimage
%
% FORMAT [f,J,D] = spm_fx_mmc(x,u,P,M)
% FORMAT [f,J]   = spm_fx_mmc(x,u,P,M)
% FORMAT [f]     = spm_fx_mmc(x,u,P,M)
% x      - state vector
%   x(:,1) - voltage     (middle pyramidal cells)
%   x(:,2) - conductance (middle pyramdidal cells)
%   x(:,3) - voltage     (superficial pyramidal cells)
%   x(:,4) - conductance (superficial pyramidal cells)
%   x(:,5) - current     (inhibitory interneurons)
%   x(:,6) - conductance (inhibitory interneurons)
%   x(:,7) - voltage     (deep pyramidal cells)
%   x(:,8) - conductance (deep pyramidal cells)
%
% f        - dx(t)/dt  = f(x(t))
% J        - df(t)/dx(t)
% D        - delay operator dx(t)/dt = f(x(t - d))
%                                    = D(d)*f(x(t))
%
% Prior fixed parameter scaling [Defaults]
%
% E  = (forward, backward, lateral) extrinsic rates
% G  = intrinsic rates
% D  = propagation delays (intrinsic, extrinsic)
% T  = synaptic time constants
% S  = slope of sigmoid activation function
%
%__________________________________________________________________________
% Copyright (C) 2016 Wellcome Trust Centre for Neuroimaging

% Bernadette van Wijk
% $Id: spm_fx_mmc.m 7185 2017-10-11 10:10:04Z spm $

% get dimensions and configure state variables
%--------------------------------------------------------------------------
x  = spm_unvec(x,M.x);            % neuronal states
n  = size(x,1);                   % number of sources

% [default] fixed parameters
%--------------------------------------------------------------------------
G  = [2 4 2 2 2 2 2 2 2 2 4 2 2 2]*200;         % intrinsic connections
T  = [3 2 12 18];                               % synaptic time constants [mp sp ii dp]

if isfield(M,'MMC_G'); G = M.MMC_G; end
if isfield(M,'MMC_T'); T = M.MMC_T; end

% [specified] fixed parameters
%--------------------------------------------------------------------------
if isfield(M,'pF')
    try, E = M.pF.E; end
    try, G = M.pF.G; end
    try, T = M.pF.T; end
    try, R = M.pF.R; end
end

% Extrinsic connections
%--------------------------------------------------------------------------
% mp = middle pyramidal
% sp = superficial pyramidal
% dp = deep pyramidal
% ii = inhibitory interneurons
%--------------------------------------------------------------------------
if n > 1
    E    = [1 1/2 1 1/2]*200;     % extrinsic (forward and backward)
    A{1} = exp(P.A{1})*E(1);      % forward  connections (sp -> mp)
    A{2} = exp(P.A{2})*E(2);      % forward  connections (sp -> sp)
    A{3} = exp(P.A{3})*E(3);      % backward connections (dp -> dp)

else
    A    = {0,0,0,0};
end

% detect and reduce the strength of reciprocal (lateral) connections
%--------------------------------------------------------------------------
for i = 1:length(A)
    L    = (A{i} > exp(-8)) & (A{i}' > exp(-8));
    A{i} = A{i}./(1 + 4*L);
end

% input connections
%--------------------------------------------------------------------------
C    = exp(P.C);

% pre-synaptic inputs: s(V)
%--------------------------------------------------------------------------
R    = 2/3;                      % slope of sigmoid activation function
B    = 0;                        % bias or background (sigmoid)
R    = R.*exp(P.S);              % gain of activation function
F    = 1./(1 + exp(-R*x + B));   % firing rate
S    = F - 1/(1 + exp(B));       % deviation from baseline firing

% input
%==========================================================================
if isfield(M,'u')

    % endogenous input
    %----------------------------------------------------------------------
    U = u(:)*1536;

else
    % exogenous input
    %----------------------------------------------------------------------
    U = C*u(:)*32;

end

% time constants and intrinsic connections
%==========================================================================
T    = ones(n,1)*T/1000;
G    = ones(n,1)*G;

% extrinsic connections
%--------------------------------------------------------------------------
% forward  (i)   2  sp -> mp (+ve)
% forward  (ii)  1  sp -> dp (+ve)
% backward (i)   2  dp -> sp (-ve)
% backward (ii)  1  dp -> ii (-ve)
%--------------------------------------------------------------------------
% free parameters on time constants and intrinsic connections
%--------------------------------------------------------------------------
% G(:,1)  mp -> mp (-ve self)  4
% G(:,2)  mp -> sp (+ve rec )  4
% G(:,3)  ii -> mp (-ve rec )  4
% G(:,4)  ii -> ii (-ve self)  4
% G(:,5)  mp -> ii (+ve rec )  4
% G(:,6)  dp -> ii (+ve rec )  2
% G(:,7)  sp -> sp (-ve self)  4
% G(:,8)  sp -> mp (+ve rec )  4
% G(:,9)  ii -> dp (-ve rec )  2
% G(:,10) dp -> dp (-ve self)  1
% G(:,11) sp -> dp (+ve rec)  2
% G(:,12) ii -> sp (-ve rec)  4
% G(:,13) sp -> ii (+ve rec)  4
% G(:,14) dp -> sp (+ve rec)  2
%--------------------------------------------------------------------------
% Neuronal states (deviations from baseline firing)
%--------------------------------------------------------------------------
%   S(:,1) - voltage     (middle pyramidal cells)
%   S(:,2) - conductance (middle pyramidal cells)
%   S(:,3) - voltage     (superficial pyramidal cells)
%   S(:,4) - conductance (superficial pyramidal cells)
%   S(:,5) - current     (inhibitory interneurons)
%   S(:,6) - conductance (inhibitory interneurons)
%   S(:,7) - voltage     (deep pyramidal cells)
%   S(:,8) - conductance (deep pyramidal cells)
%--------------------------------------------------------------------------
j     = [1 2 3 4];
for i = 1:size(P.T,2)
    T(:,j(i)) = T(:,j(i)).*exp(P.T(:,i));
end

j = 1:14;
for i = 1:size(P.G,2)
    G(:,j(i)) = G(:,j(i)).*exp(P.G(:,i));
end

% Modulatory effects of dp depolarisation on intrinsic connection j(1)
%--------------------------------------------------------------------------
if isfield(P,'M')
    G(:,j(1)) = G(:,j(1)).*exp(-P.M*32*S(:,7));
end

% Motion of states: f(x)
%--------------------------------------------------------------------------

% Conductance
%==========================================================================

% Middle layer (middle pyramidal cells): Hidden causes
%--------------------------------------------------------------------------
u      =   A{1}*S(:,3)+ U;
u      = - G(:,1).*S(:,1) - G(:,3).*S(:,5) + G(:,8).*S(:,3) + u;
f(:,2) =  (u - 2*x(:,2) - x(:,1)./T(:,1))./T(:,1);

% Supra-granular layer (superficial pyramidal cells): Hidden causes - error
%--------------------------------------------------------------------------
u      =  A{2}*S(:,3);% + U;
u      =  - G(:,7).*S(:,3) + G(:,2).*S(:,1) - G(:,12).*S(:,5) + G(:,14).*S(:,7) + u;
f(:,4) =  (u - 2*x(:,4) - x(:,3)./T(:,2))./T(:,2);

% Supra-granular layer (inhibitory interneurons): Hidden states - error
%--------------------------------------------------------------------------
u      =  0; %- A{4}*S(:,7);
u      =  - G(:,4).*S(:,5) + G(:,5).*S(:,1) + G(:,6).*S(:,7) + G(:,13).*S(:,3) + u;
f(:,6) =  (u - 2*x(:,6) - x(:,5)./T(:,3))./T(:,3);

% Infra-granular layer (deep pyramidal cells): Hidden states
%--------------------------------------------------------------------------
u      =   A{3}*S(:,7);% + U;
u      = - G(:,10).*S(:,7) - G(:,9).*S(:,5) + G(:,11).*S(:,3) + u;
f(:,8) =  (u - 2*x(:,8) - x(:,7)./T(:,4))./T(:,4);

% Voltage
%==========================================================================
f(:,1) = x(:,2);
f(:,3) = x(:,4);
f(:,5) = x(:,6);
f(:,7) = x(:,8);
f      = spm_vec(f);

if nargout < 2; return, end

% Jacobian
%==========================================================================
if isfield(M,'x'), x = spm_vec(M.x); else,  x = sparse(M.n,1); end
if isfield(M,'u'), u = spm_vec(M.u); else,  u = sparse(M.m,1); end
J  = spm_diff(M.f,x,u,P,M,1);

if nargout < 3; return, end

% delays
%==========================================================================
% Delay differential equations can be integrated efficiently (but
% approximately) by absorbing the delay operator into the Jacobian
%
%    dx(t)/dt     = f(x(t - d))
%                 = Q(d)f(x(t))
%
%    J(d)         = Q(d)df/dx
%--------------------------------------------------------------------------
% Implement: dx(t)/dt = f(x(t - d)) = inv(1 + D.*dfdx)*f(x(t))
%                     = Q*f = Q*J*x(t)
%--------------------------------------------------------------------------
Q  = spm_dcm_delay(P,M,J);

return

% notes and alpha function (kernels)
%==========================================================================
% x   = t*exp(k*t)
% x'  = exp(k*t) + k*t*exp(k*t)
%     = exp(k*t) + k*x
% x'' = 2*k*exp(k*t) + k^2*t*exp(k*t)
%     = 2*k*(x' - k*x) + k^2*x
%     = 2*k*x' - k^2*x

author: spm


text:
function qwgw_test ( )

%*****************************************************************************80
%
%% QWGW_TEST tests the QWGW library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'QWGW_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the QWGW library.\n' );

  qwgw_test01 ( );
  qwgw_test02 ( );
  qwgw_test03 ( );
  qwgw_test04 ( );
  qwgw_test05 ( );
  qwgw_test06 ( );
  qwgw_test07 ( );
  qwgw_test08 ( );
  qwgw_test09 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'QWGW_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function r = gt(p,q)
% this function compares the greater than condition of x, y, and z values of the structures that I use in
% the FVtool.
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Copyright (c) 2012-2016 Ali Akbar Eftekhari
% See the license file

if (isa(p, 'CellVector')&&isa(q, 'CellVector'))
    r=p;
    r.xvalue = p.xvalue>q.xvalue;
    r.yvalue = p.yvalue>q.yvalue;
    r.zvalue = p.zvalue>q.zvalue;
elseif isa(p, 'CellVector')
    r=p;
    r.xvalue = p.xvalue>q;
    r.yvalue = p.yvalue>q;
    r.zvalue = p.zvalue>q;
else
    r=q;
    r.xvalue = p>q.xvalue;
    r.yvalue = p>q.yvalue;
    r.zvalue = p>q.zvalue;
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function a = rutis5_inverse ( )

%*****************************************************************************80
%
%% RUTIS5_INVERSE returns the inverse of the RUTIS5 matrix.
%
%  Example:
%
%     105  167 -304  255
%     167  266 -484  406
%    -304 -484  881 -739
%     255  406 -739  620
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    John Todd,
%    Basic Numerical Mathematics, Volume 2: Numerical Algebra,
%    Academic Press, 1977.
%
%  Parameters:
%
%    Output, real A(4,4), the matrix.
%

%
%  Note that the matrix entries are listed by row.
%
  a(1:4,1:4) = [ ...
    105.0,  167.0, -304.0,  255.0; ...
    167.0,  266.0, -484.0,  406.0; ...
   -304.0, -484.0,  881.0, -739.0; ...
    255.0,  406.0, -739.0,  620.0  ];

  return
end


text:
function inside = circle_sector_contains_point_2d ( r, center, theta1, theta2, p )

%*****************************************************************************80
%
%% CIRCLE_SECTOR_CONTAINS_POINT_2D : is a point inside a circular sector?
%
%  Discussion:
%
%    A circular sector is formed by a circular arc, and the two straight line
%    segments that join its ends to the center of the circle.
%
%    A circular sector is defined by
%
%      ( X - CENTER(1) )**2 + ( Y - CENTER(2) )**2 = R**2
%
%    and
%
%      Theta1 <= Theta <= Theta2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the circle.
%
%    Input, real CENTER(2), the center of the circle.
%
%    Input, real THETA1, THETA2, the angles defining the arc,
%    in radians.  Normally, THETA1 < THETA2.
%
%    Input, real P(2), the point to be checked.
%
%    Output, logical INSIDE, is TRUE if the point is inside or on the
%    circular sector, FALSE otherwise.
%
  inside = 0;
%
%  Is the point inside the (full) circle?
%
  if ( ( p(1) - center(1) ) * ( p(1) - center(1) ) ...
     + ( p(2) - center(2) ) * ( p(2) - center(2) ) <= r * r )
%
%  Is the point's angle within the arc's range?
%  Try to force the angles to lie between 0 and 2 * PI.
%
    theta = atan4 ( p(2) - center(2), p(1) - center(1) );

    if ( r8_modp ( theta  - theta1,  2.0 * pi ) <= ...
         r8_modp ( theta2 - theta1,  2.0 * pi ) )

      inside = 1;

    end

  end

  return
end


text:
function [Index,MAX] = GpuGroup(numberOfGroups,xPrime,numberOfVariables)

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

    varsPerGroup = floor(numberOfVariables/numberOfGroups);
    if varsPerGroup ==1
        Index = gpuArray.linspace(1,numberOfVariables,numberOfVariables);
        MAX   = numberOfVariables;
    else
        B      = gpuArray.ones(1,varsPerGroup*numberOfGroups);
        remain = gpuArray.ones(1,(numberOfVariables-varsPerGroup*numberOfGroups))*(numberOfGroups+1);
        R      = reshape(B,varsPerGroup,numberOfGroups);
        k      = gpuArray.linspace(1,numberOfGroups,numberOfGroups);
        index  = R.*repmat(k,varsPerGroup,1);
        index  = reshape(index,1,varsPerGroup*numberOfGroups);
        INDEX  = [index remain];
        [~,I]  = sort(xPrime);
        Index(I) = INDEX;
        if(mod(numberOfVariables,numberOfGroups)==0)
            MAX = numberOfGroups;
        else
            MAX = numberOfGroups + 1;
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function export_h5(ebsd,fname,varargin)
% export EBSD data to a ascii file
%
% Input
%  ebsd - @EBSD
%  fname - filename
%
% Options
%  BUNGE   - Bunge convention (default)
%  ABG     - Matthies convention (alpha beta gamma)
%  DEGREE  - output in degree (default)
%  RADIANS - output in radians

if nargin > 2
  root = varargin{1};
else
  root = ['/' inputname(1)];
end

[p,f] = fileparts(fname);
fname = fullfile(p,[f '.h5']);

CSList = ebsd.CSList;
phaseMap = ebsd.phaseMap;

for k=1:numel(CSList)
  CS = CSList{k};

  cm = phaseMap(k);

  name = [root '/Phase/' num2str(cm) ];

  h5create(fname,name,1)
  h5write(fname,name,cm);

  if isa(CS,'symmetry')

    h5writeatt(fname,name,'Name',CS.pointGroup);
    h5writeatt(fname,name,'Mineral',CS.mineral);
    h5writeatt(fname,name,'Color',CS.color);
    h5writeatt(fname,name,'Laue',CS.LaueName);

    ax = norm(CS.axes);
    h5writeatt(fname,name,'a',ax(1));
    h5writeatt(fname,name,'b',ax(2));
    h5writeatt(fname,name,'c',ax(3));

    h5writeatt(fname,name,'alpha',CS.alpha);
    h5writeatt(fname,name,'beta',CS.beta);
    h5writeatt(fname,name,'gamma',CS.gamma);

    ali = CS.alignment;
    if numel(ali) >0
      ali(1:end-1) = strcat(ali(1:end-1),', ');
      ali = regexprep([ali{:}],',',', ');
      h5writeatt(fname,name,'Alignment',ali);
    end

  else

    h5writeatt(fname,name,'Name',CS);

  end

end

n = length(ebsd);

if n > 2^14
  chnk = 2^14;
else
  chnk = n;
end
opts =  {'Deflate',get_option(varargin,{'compression','Deflate'},9)};

[phi1,Phi,phi2] = Euler(ebsd.rotations,'ZXZ');

h5create(fname,[root '/Orientations'],[3 n],'Datatype','single','ChunkSize',[3 chnk],opts{:});
h5write(fname,[root '/Orientations'],double([phi1 Phi phi2]'));

h5writeatt(fname,[root '/Orientations'],'Parameterization','Euler');
h5writeatt(fname,[root '/Orientations'],'Convention','ZXZ');
h5writeatt(fname,[root '/Orientations'],'Unit','Radian');
h5writeatt(fname,[root '/Orientations'],'Mapping','Active');

h5create(fname,[root '/PhaseIndex'],[1 n],'Datatype','uint8','ChunkSize',[1 chnk],opts{:});
h5write(fname, [root '/PhaseIndex'],ebsd.phase(:).');

fn = fields(ebsd.prop);

if all(isfield(ebsd.prop,{'x','y','z'}))
  coords = [ebsd.prop.x,ebsd.prop.y,ebsd.prop.z];
elseif all(isfield(ebsd.prop,{'x','y'}))
  coords = [ebsd.prop.x,ebsd.prop.y];
end

if ~isempty(coords)
  h5create(fname,[root '/SpatialCoordinates'],[size(coords,2) n],'Datatype','single','ChunkSize',[size(coords,2) chnk],opts{:});
  h5write(fname,[root '/SpatialCoordinates'],double(coords'));

  u = ebsd.unitCell;
  if ~isempty(u)
    h5writeatt(fname,[root '/SpatialCoordinates'],'Unit','unknown')
    h5writeatt(fname,[root '/SpatialCoordinates'],'UnitCell',u)
  end
end

for k=1:numel(fn)

  p = ebsd.prop.(fn{k});

  if any(strcmp(fn{k},{'x','y','z'}))

    %     h5create(fname,['/EBSD/SpatialCoordinates/' fn{k}],n,opts{:});
    %     h5write(fname, ['/EBSD/SpatialCoordinates/' fn{k}],p);

  else
    if ~any(mod(p,1) ~= 0) % its an int

      if min(p) >= 0
        dtype = num2str(2^(sum(max(p) > 2.^(2.^(3:6))-1)+3),'uint%d');
      else
        dtype = num2str(2^(sum(max(p) > 2.^(2.^(3:6))-1)+3),'uint%d');
      end

    else
      dtype = 'single';
    end

    dim = min(size(p));
    h5create(fname,[root '/Properties/' fn{k}],[dim n],'Datatype',dtype,'ChunkSize',[dim chnk],opts{:});
    h5write(fname, [root '/Properties/' fn{k}],reshape(p',dim,[]));

  end

end

h5writeatt(fname,root,'DataType','EBSD');
h5writeatt(fname,root,'CreationTime',datestr(clock));

author: mtex-toolbox

repo: mtex

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


text:
%demo_SimpleFreeFieldHRSOS - Demonstrates the usage of the SimpleFreeFieldHRSOS conventions.

% #Author: Michael Mihocic
% #Author: Michael Mihocic: header documentation updated (28.10.2021)
%
% SOFA Toolbox - demo script
% Copyright (C) Acoustics Research Institute - Austrian Academy of Sciences
% Licensed under the EUPL, Version 1.2 or – as soon they will be approved by the European Commission - subsequent versions of the EUPL (the "License")
% You may not use this work except in compliance with the License.
% You may obtain a copy of the License at: https://joinup.ec.europa.eu/software/page/eupl
% Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing  permissions and limitations under the License.

%% Get an empy conventions structure
clear;
conventions='SimpleFreeFieldHRSOS';
disp(['Creating SOFA file with ' conventions 'conventions...']);
Obj = SOFAgetConventions(conventions);

%% Fill random data...
Obj.Data.SOS=rand(4800,2,6);

%% Update dimensions
Obj=SOFAupdateDimensions(Obj);

%% Fill with attributes
Obj.GLOBAL_ListenerShortName = 'dummy';
Obj.GLOBAL_History = 'created with a demo script';
Obj.GLOBAL_DatabaseName = 'none';
Obj.GLOBAL_ApplicationName = 'Demo of the SOFA Toolbox';
Obj.GLOBAL_ApplicationVersion = SOFAgetVersion('API');
Obj.GLOBAL_Organization = 'Acoustics Research Institute';
Obj.GLOBAL_AuthorContact = 'michael.mihocic@oeaw.ac.at';

%% save the SOFA file
SOFAfn=fullfile(SOFAdbPath,'sofatoolbox_test',[conventions '.sofa']);
disp(['Saving:  ' SOFAfn]);
Obj=SOFAsave(SOFAfn, Obj);

author: sofacoustics

repo: SOFAtoolbox


text:
function [x_g,w_g,phi,p_x,p_y,p_z] = threed_shape(x_local,r,s,t,w)
%-----------------------------------------------------------------------
%  threed_shape.m - computes test functions and derivatives on an
%                   element given element coordinates and Gauss points.
%
%  Copyright (c) 2002, Jeff Borggaard, Virginia Tech
%  Version: 1.0a
%
%  Usage:    [x_g,w_g,phi,p_x,p_y,p_z] = threed_shape(x_local,r,s,t,w)
%
%  Variables:     x_local
%                        Coordinates of the element nodes
%                 (r,s,t)
%                        Coordinates of Gauss points in unit tetrahedron
%                 w
%                        Gauss weights associated with (r,s,t)
%
%                 x_g
%                        Coordinates of Gauss points in the element
%                 w_g
%                        Gauss weights scaled by the element Jacobian
%                 phi
%                        Value of element shape functions at x_g
%                 p_x
%                 p_y
%                 p_z
%                        First spatial derivatives of phi
%                 p_xx, etc.
%                        Second spatial derivatives of phi
%                        (only applicable with 2nd order or higher)
%                        (currently not implemented)
%-----------------------------------------------------------------------
  x      = x_local;
  [n,t1] = size(x);  % t1 had better be 3
  rule   = length(r);

  if (n == 4)
    % Transform coordinates for linear elements
    c0 = x(1,:);
    c1 =-x(1,:) + x(2,:);
    c2 =-x(1,:)          + x(3,:);
    c3 =-x(1,:)                   + x(4,:);

    x_g(:,1) = c0(1) + c1(1)*r + c2(1)*s + c3(1)*t;
    xr       = c1(1);
    xs       = c2(1);
    xt       = c3(1);

    x_g(:,2) = c0(2) + c1(2)*r + c2(2)*s + c3(2)*t;
    yr       = c1(2);
    ys       = c2(2);
    yt       = c3(2);

    x_g(:,3) = c0(3) + c1(3)*r + c2(3)*s + c3(3)*t;
    zr       = c1(3);
    zs       = c2(3);
    zt       = c3(3);

    % Compute the Jacobian of the (r,s,t) -> (x,y,z) transformation
    jac = ( xr*ys*zt+xs*yt*zr+xt*yr*zs -xr*yt*zs-xs*yr*zt-xt*ys*zr );

    w_g = jac*w;

    % Invert derivatives of the mapping
    rx = ( ys*zt-yt*zs )/jac;
    ry = ( xt*zs-xs*zt )/jac;
    rz = ( xs*yt-xt*ys )/jac;

    sx = ( yt*zr-yr*zt )/jac;
    sy = ( xr*zt-xt*zr )/jac;
    sz = ( xt*yr-xr*yt )/jac;

    tx = ( yr*zs-ys*zr )/jac;
    ty = ( xs*zr-xr*zs )/jac;
    tz = ( xr*ys-xs*yr )/jac;

    % Compute shape function and derivatives at Gauss points
    phi = zeros(rule,n);
    phi(:,1) = 1 - r  - s  - t ;
    phi(:,2) =     r           ;
    phi(:,3) =          s      ;
    phi(:,4) =               t ;

    p_x = zeros(rule,n);
    p_x(:,1) =   - rx - sx - tx;
    p_x(:,2) =     rx          ;
    p_x(:,3) =          sx     ;
    p_x(:,4) =               tx;

    p_y = zeros(rule,n);
    p_y(:,1) =   - ry - sy - ty;
    p_y(:,2) =     ry          ;
    p_y(:,3) =          sy     ;
    p_y(:,4) =               ty;

    p_z = zeros(rule,n);
    p_z(:,1) =   - rz - sz - tz;
    p_z(:,2) =     rz          ;
    p_z(:,3) =          sz     ;
    p_z(:,4) =               tz;

  elseif (n == 10)
    % The following assumes straight sided elements
    c0 = x(1,:);
    c1 =-x(1,:) + x(2,:);
    c2 =-x(1,:)          + x(3,:);
    c3 =-x(1,:)                   + x(4,:);

    x_g(:,1) = c0(1) + c1(1)*r + c2(1)*s + c3(1)*t;
    xr       = c1(1);
    xs       = c2(1);
    xt       = c3(1);

    x_g(:,2) = c0(2) + c1(2)*r + c2(2)*s + c3(2)*t;
    yr       = c1(2);
    ys       = c2(2);
    yt       = c3(2);

    x_g(:,3) = c0(3) + c1(3)*r + c2(3)*s + c3(3)*t;
    zr       = c1(3);
    zs       = c2(3);
    zt       = c3(3);

    % Compute the Jacobian of the (r,s,t) -> (x,y,z) transformation
    jac = ( xr.*ys.*zt + xs.*yt.*zr + xt.*yr.*zs ...
           -xr.*yt.*zs - xs.*yr.*zt - xt.*ys.*zr );

    w_g = jac.*w;

    % Invert derivatives of the mapping
    rx = ( ys.*zt-yt.*zs )./jac;
    ry = ( xt.*zs-xs.*zt )./jac;
    rz = ( xs.*yt-xt.*ys )./jac;

    sx = ( yt.*zr-yr.*zt )./jac;
    sy = ( xr.*zt-xt.*zr )./jac;
    sz = ( xt.*yr-xr.*yt )./jac;

    tx = ( yr.*zs-ys.*zr )./jac;
    ty = ( xs.*zr-xr.*zs )./jac;
    tz = ( xr.*ys-xs.*yr )./jac;

    % Compute shape function and derivatives at Gauss points
    u  = 1 - r  - s  - t ;
    ux =   - rx - sx - tx;
    uy =   - ry - sy - ty;
    uz =   - rz - sz - tz;

    phi = zeros(rule,n);
    phi(:,1)  = u - 2*r.*u - 2*s.*u - 2*t.*u;
    phi(:,2)  = r - 2*r.*u - 2*r.*s - 2*r.*t;
    phi(:,3)  = s - 2*r.*s - 2*s.*u - 2*s.*t;
    phi(:,4)  = t - 2*r.*t - 2*s.*t - 2*t.*u;
    phi(:,5)  = 4*r.*u;
    phi(:,6)  = 4*r.*s;
    phi(:,7)  = 4*s.*u;
    phi(:,8)  = 4*r.*t;
    phi(:,9)  = 4*s.*t;
    phi(:,10) = 4*t.*u;

    p_x = zeros(rule,n);
    p_x(:,1)  = ux - 2*rx.*u - 2*r.*ux - 2*sx.*u ...
                   - 2*s.*ux - 2*tx.*u - 2*t.*ux    ;
    p_x(:,2)  = rx - 2*rx.*u - 2*r.*ux - 2*rx.*s ...
                   - 2*r.*sx - 2*rx.*t - 2*r.*tx    ;
    p_x(:,3)  = sx - 2*rx.*s - 2*r.*sx - 2*sx.*u ...
                   - 2*s.*ux - 2*sx.*t - 2*s.*tx    ;
    p_x(:,4)  = tx - 2*rx.*t - 2*r.*tx - 2*sx.*t ...
                   - 2*s.*tx - 2*tx.*u - 2*t.*ux    ;
    p_x(:,5)  = 4*rx.*u + 4*r.*ux;
    p_x(:,6)  = 4*rx.*s + 4*r.*sx;
    p_x(:,7)  = 4*sx.*u + 4*s.*ux;
    p_x(:,8)  = 4*rx.*t + 4*r.*tx;
    p_x(:,9)  = 4*sx.*t + 4*s.*tx;
    p_x(:,10) = 4*tx.*u + 4*t.*ux;

    p_y = zeros(rule,n);
    p_y(:,1)  = uy - 2*ry.*u - 2*r.*uy - 2*sy.*u ...
                   - 2*s.*uy - 2*ty.*u - 2*t.*uy    ;
    p_y(:,2)  = ry - 2*ry.*u - 2*r.*uy - 2*ry.*s ...
                   - 2*r.*sy - 2*ry.*t - 2*r.*ty    ;
    p_y(:,3)  = sy - 2*ry.*s - 2*r.*sy - 2*sy.*u ...
                   - 2*s.*uy - 2*sy.*t - 2*s.*ty    ;
    p_y(:,4)  = ty - 2*ry.*t - 2*r.*ty - 2*sy.*t ...
                   - 2*s.*ty - 2*ty.*u - 2*t.*uy    ;
    p_y(:,5)  = 4*ry.*u + 4*r.*uy;
    p_y(:,6)  = 4*ry.*s + 4*r.*sy;
    p_y(:,7)  = 4*sy.*u + 4*s.*uy;
    p_y(:,8)  = 4*ry.*t + 4*r.*ty;
    p_y(:,9)  = 4*sy.*t + 4*s.*ty;
    p_y(:,10) = 4*ty.*u + 4*t.*uy;

    p_z = zeros(rule,n);
    p_z(:,1)  = uz - 2*rz.*u - 2*r.*uz - 2*sz.*u ...
                   - 2*s.*uz - 2*tz.*u - 2*t.*uz    ;
    p_z(:,2)  = rz - 2*rz.*u - 2*r.*uz - 2*rz.*s ...
                   - 2*r.*sz - 2*rz.*t - 2*r.*tz    ;
    p_z(:,3)  = sz - 2*rz.*s - 2*r.*sz - 2*sz.*u ...
                   - 2*s.*uz - 2*sz.*t - 2*s.*tz    ;
    p_z(:,4)  = tz - 2*rz.*t - 2*r.*tz - 2*sz.*t ...
                   - 2*s.*tz - 2*tz.*u - 2*t.*uz    ;
    p_z(:,5)  = 4*rz.*u + 4*r.*uz;
    p_z(:,6)  = 4*rz.*s + 4*r.*sz;
    p_z(:,7)  = 4*sz.*u + 4*s.*uz;
    p_z(:,8)  = 4*rz.*t + 4*r.*tz;
    p_z(:,9)  = 4*sz.*t + 4*s.*tz;
    p_z(:,10) = 4*tz.*u + 4*t.*uz;

  end


text:
function [GAmp,GTime]=ADCCartesian(p)

global VCtl;
global VObj;

tMiddle=p.tMiddle;
Duplicates=max(1,p.Duplicates);
DupSpacing=max(0,p.DupSpacing);

% acqusition lobs
GxAmp=(1/VCtl.FOVFreq)/((VObj.Gyro/(2*pi))*(1/VCtl.BandWidth));
tHalf=1/(2*(VObj.Gyro/(2*pi))*GxAmp*VCtl.RFreq);
[GAmp1,GTime1]=StdTrap(tMiddle+VCtl.TEAnchorTime-tHalf-VCtl.MinUpdRate, ...
                       tMiddle+VCtl.TEAnchorTime+tHalf+VCtl.MinUpdRate, ...
                       tMiddle+VCtl.TEAnchorTime-tHalf,               ...
                       tMiddle+VCtl.TEAnchorTime+tHalf,               ...
                       1,2,VCtl.ResFreq,2);

GAmp=[GAmp1];
GTime=[GTime1];
[GTime,m,n]=unique(GTime);
GAmp=GAmp(m);

% Create Duplicates
if Duplicates~=1 & DupSpacing ~=0
    GAmp=repmat(GAmp,[1 Duplicates]);
    TimeOffset = repmat(0:DupSpacing:(Duplicates-1)*DupSpacing,[length(GTime) 1]);
    GTime=repmat(GTime,[1 Duplicates]) + (TimeOffset(:))';
end

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function centroids = kMeansInitCentroids(X, K)
%KMEANSINITCENTROIDS This function initializes K centroids that are to be
%used in K-Means on the dataset X
%   centroids = KMEANSINITCENTROIDS(X, K) returns K initial centroids to be
%   used with the K-Means on the dataset X
%

% You should return this values correctly
centroids = zeros(K, size(X, 2));

% ====================== YOUR CODE HERE ======================
% Instructions: You should set centroids to randomly chosen examples from
%               the dataset X
%
% Randomly reorder the indices of examples
randidx = randperm(size(X, 1));
% Take the first K examples as centroids
centroids = X(randidx(1:K), :);

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

end

author: zzlyw

repo: machine-learning-exercises


text:
function  Intan2PLX(parameterSets,tet,varargin)

% function converts an Intan file to a plx file that can be read by Offline
% Sorter etc..http://www.plexon.com/software-downloads
%
% calls read_intan_data available at http://www.intantech.com/downloads.html
% requires spike detection software package from http://urut.ch/new/serendipity/index.php?/pages/osort.html
% calls scripts for writing PLX: checkPLXData, writePLXFile.
% readPLXHeaders is not called but is a useful reference
%
% INPUTS
%
% parameterSets:  %1 -> wavelet, 2->power
% tet: true -> groups electrodes by 4 , false ->treats electrode individually
%filename:
% ex: Intan2PLX(2,true,filename)

%define parameters of spike sort

params=[];
params.bandPass=[300 3000];
params.nrNoiseTraces=0;
params.prewhiten=0;
params.samplingFreq=40000; %intan sampling frequency
params.limit=100;

%set up filter for LFP
n = 4;
Wn =params.bandPass/(params.samplingFreq/2);
[b,a] = butter(n,Wn);
Hd=[];
Hd{1}=b;
Hd{2}=a;

switch( parameterSets )
    case 1
        %wavelet detection
        params.extractionThreshold = 0;
        params.alignMethod=2; %1 pos, 2 neg, 3 mixed
        params.peakAlignMethod=4; %1 findPeak, 2 none, 3 power, 4 mteo
        params.detectionMethod=5;
        dp3=[];
        dp3.scalesRange=[0.2 1.0];
        dp3.waveletName='bior1.5';
        params.detectionParams=dp3;
    case 2
        %power detection
        params.extractionThreshold = 6;
        params.alignMethod=2; %1 pos, 2 neg, 3 mixed
        params.peakAlignMethod=1; %1 findPeak, 2 none, 3 power, 4 mteo
        params.detectionMethod=1; %1 power, 5 WDM
        dp1=[];
        dp1.kernelSize=18;
        params.detectionParams=dp1;
    otherwise
        error('invalid parameter set nr');
end

if  length(varargin)>0
    fileList{1} = varargin{1};
else
    dirName=uigetdir;
    fileList = getAllIntanFiles(dirName);
    answer = inputdlg('Enter basefile');
    fileList = fileList(cellfun(@(a) any(regexp(a,answer{1})),fileList));

    %sort files by date
    for i=1:length(fileList)

    file = dir(fileList{i});
    filets(i,:)=datevec(file.date);
    end

    [a,b]=sortrows(filets);
    fileList=fileList(b);
end

%determine plexon file destination
[filename, pathname] = uiputfile('*.plx', 'Select a plx file');
fname = strcat(pathname, filename(1:end-4));

for ii = 1:length(fileList)
    %read intan data
    [t,amps,temp,aux] = read_intan_data(fileList{ii});
    data=zeros(length(temp),16); %hardcoded to 16 channels
    data(:,amps)=temp; %fill in dead channels as zeros
    data = data/1000;  %convert from uV to mV (plx standard)

    %%load spikes%%

    %initialize data output
    plexon.ts = nan(1e6,1);
    plexon.waves= nan(1e6,32); %approaching max memory limit
    plexon.units =nan(1e6,1);
    plexon.chans  =nan(1e6,1);

    %find spikes

    num_spks=0;

    %read 4 at a time if trodal data
    if tet
        skip = 4;
    else
        skip = 1;
    end

    for i=1:skip:size(data,2)

        %extract spikes from filtered LFP

        if tet
            spikeTimestamps_e{4}=[];
            filteredSignal_e{4}=[];
            for k=1:4
                [~, filteredSignal_e{k}, ~,~, spikeTimestamps_e{k}, ~, ~, ~] = extractSpikes(data(:,i+k-1), Hd, params );

            end

            %enslave signal from other channels to spikes on any electrode
            [spikeWaveforms,spikeTimestamps] = getTetrodeData(filteredSignal_e,spikeTimestamps_e,params.samplingFreq);

        else

            %not tested, likely will not work
            [~, filteredSignal, ~,spikeWaveforms, spikeTimestamps, ~, ~, ~] = extractSpikes(data(:,i), Hd, params );
        end

        %store vector array

        nspk=length(spikeTimestamps)*4;

        plexon.ts(num_spks+1:num_spks+nspk) = repmat(spikeTimestamps',4,1); %repeat ts for each unit
        plexon.waves(num_spks+1:num_spks+nspk,:) = spikeWaveforms; %all waveforms
        plexon.units(num_spks+1:num_spks+nspk) = zeros(nspk,1); %no unit assignments
        plexon.chans(num_spks+1:num_spks+nspk) = upSample(i:i+skip-1,length(spikeTimestamps)); %channel assignments

        num_spks=num_spks+nspk;

    end

    %remove excess variable space
    plexon.ts = plexon.ts(~isnan(plexon.ts));
    plexon.waves = plexon.waves(~any(isnan(plexon.waves),2),:);
    plexon.units = plexon.units(~isnan(plexon.units));
    plexon.chans = plexon.chans(~isnan(plexon.chans));

    %generate default headings
    headers = checkPLXData(plexon);

    %change headings specific to Intan
    headers.spikePreAmpGain=1; %signal already corrected for system gains
    headers.bitsperslowsample=16;
    headers.bitsperspikesample=16; %resolution of spike (y-axis)
    headers.ADFrequency=params.samplingFreq;

    headers.spikeMaxMagnitudeMV=max(abs(plexon.waves(:)));
    headers.slowMaxMagnitudeMV=max(abs(plexon.waves(:)));
    headers.waveformfreq=params.samplingFreq;
    plexon.waves = (plexon.waves * (.5*2^headers.bitsperslowsample))/ headers.spikeMaxMagnitudeMV; %convert to 16 bit
    plexon.waves(end-1,:)=plexon.waves(end,:);
    % write file

    writePLXFile([fname '_f' num2str(ii) '.plx'], plexon,headers)
    disp(['Saved in ' fname])
end

end

function [spikeWaveforms_tet,spikeTimestamps_tet] = getTetrodeData(filteredSignal,spikeTimestamps,fs)

filteredSignal=cell2mat(filteredSignal);
t_ind=1:size(filteredSignal,1);
spikes=false(size(t_ind));

%combine all spike times from a tetrode and sort
t=cell2mat(spikeTimestamps);
t=t(t>16 & t<size(filteredSignal,1)-17);
t = sort(t);

%find blocks with spikes on multiple channels or within 32 samples
if any(t)
epoch=repmat(-15:16,length(t),1)+repmat(t',1,32);
spikes(epoch(:))=true;
onSet=find(diff(spikes)==1);
offSet=find(diff(spikes)==-1);

%select biggest spike event in that block
for i=1:length(onSet)

    inT=t(t>onSet(i) & t<offSet(i));
    inT=unique(inT);
    inBlock = ismember(1:size(filteredSignal,1),inT);
    maxBlock = filteredSignal(inBlock,:);
    [r,~] = find(maxBlock==min(maxBlock(:)));
    ditch=inT(inT~=inT(r));

    t(ismember(t,ditch))=nan; %only keep biggest spike

end
t(isnan(t))=[];
t=unique(t);
spikeTimestamps_tet=t/fs; %convert to seconds
t=repmat(-15:16,length(t),1)+repmat(t',1,32);

temp1=filteredSignal(:,1); %e1
temp1=temp1(t);
temp2=filteredSignal(:,2); %e2
temp2=temp2(t);
temp3=filteredSignal(:,3); %e3
temp3=temp3(t);
temp4=filteredSignal(:,4); %e4
temp4=temp4(t);

spikeWaveforms_tet=[temp1;temp2;temp3;temp4];

else
    spikeWaveforms_tet=zeros(0,16);
    spikeTimestamps_tet=[];
end

end

function X= upSample (input, bins)
bins=round(bins);
input=input(:);
X=nan(length(input)*bins,1);

for i = 0:bins-1
    X(1+i:bins:end)=input;
end

end

author: buzsakilab

repo: buzcode


text:
function [strshocks_record]=mastrshocks(beta_gibbs,delta_gibbs,D_record,Y,X,Z,n,k1,k3,It,Bu)

% first create the call storing the results
strshocks_record=cell(n,1);

% then loop over iterations
for ii=1:It-Bu

% recover the VAR coefficients, reshaped for convenience
B=reshape(beta_gibbs(:,ii),k1,n);
Delta=reshape(delta_gibbs(:,ii),k3,n);

% obtain the residuals from (XXX)
EPS=Y-X*B-Z*Delta;

% then recover the structural marix D
D=reshape(D_record(:,ii),n,n);

% obtain the structural disturbances from (XXX)
ETA=D\EPS';

   % save in struct_shocks_record
   for jj=1:n
   strshocks_record{jj,1}(ii,:)=ETA(jj,:);
   end

end

author: european-central-bank

repo: BEAR-toolbox


text:
function f_hat = calcFourier(SO3F,varargin)
% compute harmonic coefficients of SO3Fun
%
% Syntax
%
%  f_hat = calcFourier(SO3F)
%
%  f_hat = calcFourier(SO3F,'bandwidth',L)
%
% Input
%  SO3F - @SO3FunRBF
%     L - maximum harmonic degree / bandwidth
%
% Output
%  f_hat - harmonic/Fouier/Wigner-D coefficients
%

L = get_option(varargin,'bandwidth',min(SO3F.bandwidth,getMTEXpref('maxSO3Bandwidth')));

cs = SO3F.CS.properGroup;
ss = SO3F.SS.properGroup;

% the weights
c = SO3F.weights;

% the center orientations
ori = SO3F.center;

% for a few center symmetrize before computing c_hat
symCenter = 10*length(SO3F.center) * numSym(cs) * numSym(ss) < max(L^3,100);
if symCenter
  ori = symmetrise(SO3F.center,'proper');
  c = repmat(c(:).',size(ori,1),1) / numSym(cs) / numSym(ss);
end

SO3FH = SO3FunHarmonic.quadrature(ori,c,varargin{:},'bandwidth',L);

SO3FH = conv(SO3FH,SO3F.psi);

if ~symCenter, SO3FH = symmetrise(SO3FH); end

% add constant portion
SO3FH = SO3FH + SO3F.c0;

f_hat = SO3FH.fhat;


text:
function proposals=calclpoForIm(input, lpoconfig )
    if ~isfield(lpoconfig,'path')
        fprintf('Path to LPO does not exist. Please make sure you give a proper full path\n');
        return;
    else
        datapath=[lpoconfig.path, '/data/'];
        modelpath=[lpoconfig.path, '/models/'];
    end
    if(isstr(input))
        im = im2uint8(imread(input));
    else
        im = im2uint8(input); % Just to make input consistent
    end
    if(size(im, 3) == 1)
        im=repmat(im,[1,1,3]);
    end
    model=[modelpath, 'lpo_VOC_0.1.dat'];
    detector='MultiScaleStructuredForest';
    max_iou=0.9;
    if isfield(lpoconfig.params,'model')
        model=[modelpath, lpoconfig.params.model];
    end
    if isfield(lpoconfig.params,'detector')
        detector=lpoconfig.params.detector;
    end
    if isfield(lpoconfig.params,'max_iou')
        max_iou=lpoconfig.params.max_iou;
    end
    % Set a boundary detector by calling (before creating an OverSegmentation!):
    lpo_mex( 'setDetector', [detector '("'  datapath 'sf.dat")'] );

    p = LpoProposal( 'model', model, 'box_nms', 'max_iou', max_iou );
    os = LpoOverSegmentation( im );
    [segs,props] = p.propose( os );
    boxes = masksToBoxes( segs, props );
    if(isfield(lpoconfig.opts,'numProposals'))
        numProposals=lpoconfig.opts.numProposals;
        if(size(boxes,1)>=numProposals)
            boxes=boxes(1:numProposals,:);
        else
            fprintf('Only %d proposals were generated for the input image\n',size(boxes,1));
        end
    end
    proposals.boxes=boxes;
end

author: Cloud-CV

repo: object-proposals


text:
%% Setup
% folder setup
% saveFigFlag = 1;

outputDir = GetOutputDataDir;
% saveDir = [];
% saveDir{1} = fullfile(outputDir,'fig2_0422','PT_LR_stimrangePT');
% setDir(saveDir{1}); % make folder if doesn't exist
% saveDir{2} = fullfile(outputDir,'fig2_0422','LonRon_LR_stimrangePT');
% setDir(saveDir{2}); % make folder if doesn't exist
% saveDir{3} = fullfile(outputDir,'fig2_0422','motor_LR_stimrangePT');
% setDir(saveDir{3}); % make folder if doesn't exist

% params
reg_name = 'Fw_fc_reg0.5';
isStimNotMotor = 0; % 1

reg_range = 2; % 4
reg_thres = 0.5;
n_reg = 1;
% M_stimmotorflag = [1,1,0]; % 1 for stim and 0 for motor
% M_reg_range = {[3,2],[6,5],[4,4],[1,3]}; % [2,3] for PT, [8,9] for OMR for fish8; left/right pairs
% M_reg_thres = {0.5,0.5,0.5};
% n_reg = length(M_reg_thres);

%% Init load
hfig = figure;
InitializeAppData(hfig);
ResetDisplayParams(hfig);

%% Load fish
range = GetFishRange;%1:18;%
IM = cell(n_reg,max(range));
IM_full = cell(n_reg,max(range));
%%
i_set = 1;
for i_fish = range

    setappdata(hfig,'isMotorseed',0);

    ClusterIDs = GetClusterIDs('all');
    stimrange = 1;
    [cIX_load,gIX_load,M,stim,behavior,M_0] = LoadSingleFishDefault(i_fish,hfig,ClusterIDs,stimrange);

    %% regression
    setappdata(hfig,'isStimAvr',0);
    [M_0,M,behavior,stim] = UpdateTimeIndex(hfig);

    fishset = getappdata(hfig,'fishset');
    if isStimNotMotor
        [~,names,regressors] = GetStimRegressor(stim,fishset,i_fish);
    else
        [~,names,regressors] = GetMotorRegressor(behavior,i_fish);
    end

    %% regression
    Reg = regressors(reg_range,:);
    Corr = corr(Reg',M_0');

    %% choose top 2% of all cells
    if 0
        % keep best regression only
        [Corr_rows,corr_max] = ConvertCorrToBestRegRows(Corr);

        % top 2 %
        nCells_total = size(M_0,1);
        prct_const = 2; % default using 2
        [CIX,RegThres] = ChooseTopPercentOfFish(nCells_total,prct_const,Corr_rows);

        cIX = CIX{1};
        % get map color
        reg_thres = 0.25;
        gIX = MapXto1Dcolormap(corr_max(cIX),[reg_thres,1],64);

    else
        %% threshold by reg coeff
        cIX = (find(Corr>reg_thres))';
        wIX = Corr(cIX); % weight, i.e. corr.coeff

        gIX = MapXto1Dcolormap(wIX,[reg_thres,1],64);
    end
    %% set up custom colormap: hot (heatmap)
%     clrbins = 0.5:0.05:1;
%     clrmap0 = hot(round(1.5*length(clrbins)));
%     clrmap = clrmap0(end-length(clrbins):end,:);
%
%     gIX = interp1(clrbins,1:length(clrbins),wIX,'nearest');
%
    %% single color gradient (red)
    clr1 = [1,0,0];
    clr1_ = [0.5,0.4,0.4];
%     clr1_ = [0.7,0.5,0.5];
    numC = 64;
    clrmap = Make1DColormap([clr1_;clr1],numC);

    %% right plot
    I = LoadCurrentFishForAnatPlot(hfig,cIX,gIX);
    I.clrmap = clrmap;
    [h,im_full,im] = DrawCellsOnAnat(I);

    %% save figure
    close(h);
    IM{i_set,i_fish} = im;
    IM_full{i_set,i_fish} = im_full;
    %         savefolder = fullfile(saveDir{i_set},['regthres' num2str(reg_thres)]);
    %         setDir(savefolder);
    %         figName = ['Fish' num2str(i_fish)];
    % %         SaveFigureHelper(saveFigFlag, savefolder, figName);

end

%% save as tiff stack
i_reg = 1;
range_im = range;%1:18;
tiffdir = fullfile(outputDir,[reg_name '_allfish.tiff']);
% tiffdir = fullfile(outputDir,'White_1reg_allfish.tiff');

% display each plane and save as tif
h = figure;
for i_plane = range_im
    im = IM_full{i_reg,i_plane};
    image(im);axis equal; axis off
    drawnow;
    % save tiff
    if (i_plane == 1)
        imwrite(im, tiffdir, 'compression','none','writemode','overwrite')
    else
        imwrite(im, tiffdir, 'compression','none','writemode','append')
    end
    %     pause(0.2)
end
close(h)

%% [for later] plot from tiff stack
isPlotfromtiffstack = 0;
if isPlotfromtiffstack
    IM_full = cell(1,18);
    for i = 1:18
        im = double(imread(tiffdir,i))./255;
        IM_full{i} = im(317:1236,1:621,:);
    end
end

%% save avrage plot
i_reg = 1;
range_im = [1:3,5:18];
cellarray = IM_full(i_reg,range_im);

k_scale = 1; % this changes for every reg pair...
k_contrast = 1;

[h_anat,im_avr] = AverageAnatPlot(cellarray,k_contrast,k_scale);
imwrite(im_avr, fullfile(outputDir,[reg_name '_avr.tiff']), 'compression','none','writemode','overwrite');

%%
% anat_YX = getappdata(hfig,'anat_yx_norm');
% z_range_ventral = 40;
% z_range = z_range_ventral:345; % max 138*k_zres
% y_range = 81:1000;%81:990;% 81:1104; % max 1406
% dimv_yx3 = size(anat_YX(y_range,:,:));
% dimv_zx3 = [length(z_range),size(anat_YX,2),3];
% im_full = IM_full{1};
% im = IM_full{1}(dimv_zx3(1)+11:end,1:dimv_yx3(2),:);

%% Crop Anat to show only XY view or only XY/YZ view
IM_crop = CropFullAnat(IM_full,2);

cellarray = IM_crop(1,1:4);
[h_anat,im_avr] = AverageAnatPlot(cellarray,k_contrast,k_scale);

author: xiuyechen

repo: FishExplorer


text:
function mk=armijo(xk,dk )
    beta=0.5;
    sigma=0.2;
    m=0;
    mmax=20;
    while (m<=mmax)
        if(fun(xk+beta^m*dk)<=fun(xk)+sigma*beta^m*gfun(xk)'*dk)
            mk=m;
            break;
        end
    m=m+1;
    end
    alpha=beta^mk
    newxk=xk+alpha*dk
    fk=fun(xk)
    newfk=fun(newxk)

author: Hannah-Zhou