text:
%
report_this_filefun(mfilename('fullpath'));
%
%  This file display the original menu
%
%  Stefan Wiemer 12/94

%  Create new figure
% Find out of figure already exists
%
[existFlag,figNumber]=figure_exists('ZMAP 5.0 - Menu');
newmenuFlag=~existFlag;
my_dir = ' ';

% Set up the Seismicity Map window Enviroment
%
if newmenuFlag
    menu  = figure_w_normalized_uicontrolunits( ...
        'Name','ZMAP 5.0 - Menu',...
        'NumberTitle','off', ...
        'MenuBar','none', ...
        'NextPlot','new', ...
        'Color',color_fbg,...
        'Visible','off', ...
        'Position',[ wex wey-300 (ZmapGlobal.Data.map_len - [300 150]) ]);
    orient tall
    axis off

    uicontrol(gcf,...
        'Units','normalized',...
        'Callback','showweb(''new'');',...
        'String','Introduction',...
        'BackgroundColor' ,[0.8 0.8 0.8]',...
        'Position',[0.10 0.90 0.80 0.09]);

    uicontrol(gcf,...
        'Units','normalized',...
        'Callback','think;load_zmapfile()',...
        'String','Load *.mat  Datafile and go',...
        'BackgroundColor' ,[0.8 0.8 0.8]',...
        'Position',[0.10 0.80 0.80 0.09]);

    uicontrol(gcf,...
        'Units','normalized',...
        'Callback','think;setup',...
        'String','Create or modify *.mat Datafile',...
        'BackgroundColor' ,[0.8 0.8 0.8]',...
        'Position',[0.10 0.70 0.80 0.09]);

    uicontrol(gcf,...
        'Units','normalized',...
        'Callback','think;datinf',...
        'String','Current Dataset Info ',...
        'BackgroundColor' ,[0.8 0.8 0.8]',...
        'Position',[0.10 0.60 0.80 0.09]);

    % uicontrol(gcf,...
    %'Units','normalized',...
    %'Callback','think;loadcal;',...
    %'String','Load new NCEC catalog ',...
    %'BackgroundColor' ,[0.8 0.8 0.8]',...
    %'Position',[0.10 0.50 0.80 0.09]);

    uicontrol(gcf,...
        'Units','normalized',...
        'Callback','think;newccat=[];clusmenu',...
        'String','Load Cluster datafile ',...
        'BackgroundColor' ,[0.8 0.8 0.8]',...
        'Position',[0.10 0.50 0.80 0.09]);

    uicontrol(gcf,...
        'Units','normalized',...
        'Callback','web(''http://www.seismo.ethz.ch/staff/stefan/zmapnews/papers.html'');',...
        'String',' Papers written using ZMAP',...
        'BackgroundColor' ,[0.8 0.8 0.8]',...
        'Position',[0.10 0.40 0.80 0.09]);

    uicontrol(gcf,...
        'Units','normalized',...
        'Callback','think;dosl = ''newf'';slshow ',...
        'String','Sample Slide Show',...
        'BackgroundColor' ,[0.8 0.8 0.8]',...
        'Position',[0.10 0.30 0.80 0.09]);

    uicontrol(gcf,...
        'Units','normalized',...
        'Callback','think;startwebd',...
        'String','Start WEB Doumentation ',...
        'BackgroundColor' ,[0.8 0.8 0.8]',...
        'Position',[0.10 0.20 0.80 0.09]);

    uicontrol(gcf,...
        'Units','normalized',...
        'Callback','think;working_dir_in',...
        'String','Set working directory ',...
        'BackgroundColor' ,[0.8 0.8 0.8]',...
        'Position',[0.10 0.10 0.80 0.09]);

    uicontrol(gcf,...
        'Units','normalized',...
        'Callback','close(menu);zmap_message_center.set_info('' '','' '');done',...
        'String','Close this window ',...
        'BackgroundColor' ,[0.8 0.8 0.8]',...
        'Position',[0.15 0.00 0.70 0.09]);

end   % if fig exist

figure_w_normalized_uicontrolunits(menu)
set(gcf,'visible','on')


text:
function op = smooth_entropy()
%SMOOTH_ENTROPY The entropy function -sum( x_i log(x_i) )
op = @smooth_entropy_impl;

function [ v, g ] = smooth_entropy_impl( x )
if any( x < 0 ),
    v = -Inf;
    if nargout > 1,
        g = NaN * ones(size(x));
    end
else
    logx = log(max(x,realmin));
    v = - tfocs_dot( x, logx );
    if nargout > 1,
        g = - logx - 1;
    end
end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function res = sensors(this, type, newsens)
% Sets and gets sensor fields for EEG and MEG
% returns empty matrix if no sensors are defined.
% FORMAT res = sensors(this, type, newsens)
%   type - 'EEG' or 'MEG'
% _______________________________________________________________________
% Copyright (C) 2008-2012 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: sensors.m 6829 2016-07-07 10:16:46Z vladimir $

if nargin<2
    error('Sensor type (EEG or MEG) must be specified');
end

switch lower(type)
    case 'eeg'
        if nargin < 3
            if isfield(this.sensors, 'eeg')
                res = this.sensors.eeg;
            else
                res = [];
            end
        else
            this.sensors(1).eeg = newsens;
            res = check(this);
        end
    case 'meg'
        if nargin < 3
            if isfield(this.sensors, 'meg')
                res = this.sensors.meg;
            else
                res = [];
            end
        else
            this.sensors(1).meg = newsens;
            res = check(this);
        end
    case 'src'
        if nargin < 3
            if isfield(this.sensors, 'src')
                res = this.sensors.src;
            else
                res = [];
            end
        else
            this.sensors(1).src = newsens;
            res = check(this);
        end
    otherwise
        error('Unsupported sensor type');
end

if  nargin < 3 && ~isempty(res) && ismember(lower(type), {'eeg', 'meg'}) && this.montage.Mind > 0
    sens = res;
    montage = this.montage.M(this.montage.Mind);
    if  ~isempty(intersect(sens.label, montage.labelorg))
        sensmontage = montage;
        [sel1, sel2] = spm_match_str(sens.label, sensmontage.labelorg);
        sensmontage.labelorg = sensmontage.labelorg(sel2);
        sensmontage.tra = sensmontage.tra(:, sel2);
        selempty  = find(all(sensmontage.tra == 0, 2));
        sensmontage.tra(selempty, :) = [];
        sensmontage.labelnew(selempty) = [];

        if isfield(sensmontage, 'chantypeorg')
            sensmontage.chantypeorg = sensmontage.chantypeorg(sel2);
        end
        if isfield(sensmontage, 'chanunitorg')
            sensmontage.chanunitorg = sensmontage.chanunitorg(sel2);
        end
        if isfield(sensmontage, 'chantypenew')
            sensmontage.chantypenew(selempty) = [];
        end
        if isfield(sensmontage, 'chanunitnew')
            sensmontage.chanunitnew(selempty) = [];
        end

        chanunitorig = sens.chanunit(sel1);
        chantypeorig = sens.chantype(sel1);
        labelorg     = sens.label;

        keepunused   = 'no'; % not sure if this is good for all cases

        sens = ft_apply_montage(sens, sensmontage, 'keepunused', keepunused);

        if strcmpi(type, 'MEG')
            if isfield(sens, 'balance') && ~isequal(sens.balance.current, 'none')
                balance = ft_apply_montage(getfield(sens.balance, sens.balance.current), sensmontage, 'keepunused', keepunused);
            else
                balance = sensmontage;
            end

            sens.balance.custom = balance;
            sens.balance.current = 'custom';
        end

        % If all the original channels contributing to a new channel have
        % the same units, transfer them to the new channel. This might be
        % wrong if the montage itself changes the units by scaling the data.
        if isfield(sens, 'chanunit')
            chanunit = sens.chanunit;
        else
            chanunit = repmat({'unknown'}, numel(sens.label), 1);
        end

        if isfield(sens, 'chantype')
            chantype = sens.chantype;
        else
            chantype = repmat({'unknown'}, numel(sens.label), 1);
        end

        for j = 1:numel(sens.label)
            k = strmatch(sens.label{j}, sensmontage.labelnew, 'exact');
            if ~isempty(k)
                if isequal(chanunit{j}, 'unknown')
                    unit = unique(chanunitorig(~~abs(sensmontage.tra(k, :))));
                    if numel(unit)==1
                        chanunit(j) = unit;
                    elseif strcmpi(type, 'MEG')
                        chanunit{j} = 'T';
                    else
                        chanunit{j} = 'V';
                    end
                end

                if isequal(chantype{j}, 'unknown')
                    ctype = unique(chantypeorig(~~abs(sensmontage.tra(k, :))));
                    if numel(ctype)==1
                        chantype(j) = ctype;
                    elseif strcmpi(type, 'MEG')
                        chantype{j} = 'megmag';
                    else
                        chantype{j} = 'eeg';
                    end
                end
            else %channel was not in the montage, but just copied
                k = strmatch(sens.label{j}, labelorg, 'exact');
                if isequal(chanunit{j}, 'unknown')
                    chanunit(j) = chanunitorig(k);
                end
                if isequal(chantype{j}, 'unknown')
                    chantype(j) = chantypeorig(k);
                end
            end
        end

        sens.chanunit = chanunit;
        sens.chantype = chantype;
    end

    if strcmpi(type, 'MEG')
        res = ft_datatype_sens(sens, 'amplitude', 'T', 'distance', 'mm');
    else
        res = ft_datatype_sens(sens, 'amplitude', 'V', 'distance', 'mm');
    end
end

author: spm


text:
function femI = genP1IFEM3D(mesh,fem,btM,btP)
%% Usage: Generate Quadrature Information on Interface Elements
%         Each interface tetrahedron is cut into 4 to 9 small tetrahedra
%
% femI.t --- global index on small tetrahedra
% femI.bas --- nt-4-4 matrix stores basis functions on all small tetra
% femI.gx --- Gaussian x nodes
% femI.gy --- Gaussian y nodes
% femI.gz --- Gaussian z nodes
% femI.gw --- Gaussian weights
% femI.area --- Areas of all small tetra
% femI.quadT --- number of small tetras of each interface element
% femI.locID --- order of local index, to match with the basis function

% ***** The following info only used for partial penalty IFEM *****

% femI.bas1 --- nti-4-4 basis function on 1st piece of interface element
% femI.bas2 --- nti-4-4 basis function on 2nd piece of interface element
% femI.plusPC --- 1 or 2 denotes which piece corresponds to btP.

% Last Modified by Xu Zhang 08/02/2020

%%
ntI = -min(mesh.tLoc);
t = zeros(8*ntI,4); bas = zeros(8*ntI,4,4); A = zeros(8*ntI,1);
gx = zeros(8*ntI,4); gy = zeros(8*ntI,4); gz = zeros(8*ntI,4);
QelemID = zeros(8*ntI,1);
quadT = zeros(ntI,1); locIDvec = zeros(ntI,4);
bas1 = zeros(ntI,4,4); bas2 = zeros(ntI,4,4); plusPC = zeros(ntI,1);
intID = find(mesh.tLoc<0);
id = 0;
for i = 1:ntI
    tID = intID(i);
    t_e = mesh.t_e(tID,:);
    intpt = mesh.eIntP(-mesh.eLoc(t_e(mesh.eLoc(t_e)<0)),:);
    pLocK = mesh.pLoc(mesh.t(tID,:));
    vert0 = mesh.p(mesh.t(tID,:),:);
    if size(intpt,1) == 3 % Type 1 Interface Element (3 intersection pts)
        if pLocK(1) == -1 && pLocK(2) == 1 && pLocK(3) == 1 && pLocK(4) == 1
            locID = [1,2,3,4]; vert = vert0(locID,:); coef = [btM,btP]; plusPC(i) = 2;
        elseif pLocK(1) == 1 && pLocK(2) == -1 && pLocK(3) == 1 && pLocK(4) == 1
            locID = [2,3,4,1]; vert = vert0(locID,:); coef = [btM,btP]; plusPC(i) = 2;
        elseif pLocK(1) == 1 && pLocK(2) == 1 && pLocK(3) == -1 && pLocK(4) == 1
            locID = [3,4,1,2]; vert = vert0(locID,:); coef = [btM,btP]; plusPC(i) = 2;
        elseif pLocK(1) == 1 && pLocK(2) == 1 && pLocK(3) == 1 && pLocK(4) == -1
            locID = [4,1,2,3]; vert = vert0(locID,:); coef = [btM,btP]; plusPC(i) = 2;
        elseif pLocK(1) == 1 && pLocK(2) == -1 && pLocK(3) == -1 && pLocK(4) == -1
            locID = [1,2,3,4]; vert = vert0(locID,:); coef = [btP,btM]; plusPC(i) = 1;
        elseif pLocK(1) == -1 && pLocK(2) == 1 && pLocK(3) == -1 && pLocK(4) == -1
            locID = [2,3,4,1]; vert = vert0(locID,:); coef = [btP,btM]; plusPC(i) = 1;
        elseif pLocK(1) == -1 && pLocK(2) == -1 && pLocK(3) == 1 && pLocK(4) == -1
            locID = [3,4,1,2]; vert = vert0(locID,:); coef = [btP,btM]; plusPC(i) = 1;
        elseif pLocK(1) == -1 && pLocK(2) == -1 && pLocK(3) == -1 && pLocK(4) == 1
            locID = [4,1,2,3]; vert = vert0(locID,:); coef = [btP,btM]; plusPC(i) = 1;
        end
        p1 = [vert(1,:);intpt]; t1 = delaunay(p1);
        p2 = [vert(2:4,:);intpt];t2 = delaunay(p2);
    elseif size(intpt,1) == 4 % Type 2 Interface Element (4 intersection pts)
        if pLocK(1) == -1 && pLocK(2) == -1 && pLocK(3) == 1 && pLocK(4) == 1
            locID = [1,2,3,4]; vert = vert0(locID,:); coef = [btM,btP]; plusPC(i) = 2;
        elseif pLocK(1) == -1 && pLocK(2) == 1 && pLocK(3) == -1 && pLocK(4) == 1
            locID = [1,3,2,4]; vert = vert0(locID,:); coef = [btM,btP]; plusPC(i) = 2;
        elseif pLocK(1) == -1 && pLocK(2) == 1 && pLocK(3) == 1 && pLocK(4) == -1
            locID = [1,4,2,3]; vert = vert0(locID,:); coef = [btM,btP]; plusPC(i) = 2;
        elseif pLocK(1) == 1 && pLocK(2) == -1 && pLocK(3) == -1 && pLocK(4) == 1
            locID = [1,4,2,3]; vert = vert0(locID,:); coef = [btP,btM]; plusPC(i) = 1;
        elseif pLocK(1) == 1 && pLocK(2) == -1 && pLocK(3) == 1 && pLocK(4) == -1
            locID = [1,3,2,4]; vert = vert0(locID,:); coef = [btP,btM]; plusPC(i) = 1;
        elseif pLocK(1) == 1 && pLocK(2) == 1 && pLocK(3) == -1 && pLocK(4) == -1
            locID = [1,2,3,4]; vert = vert0(locID,:); coef = [btP,btM]; plusPC(i) = 1;
        end
        p1 = [vert(1:2,:);intpt]; t1 = delaunay(p1);
        p2 = [vert(3:4,:);intpt]; t2 = delaunay(p2);
    else
        stp = 1;
    end
    [BAS1,BAS2] = basIFE3DP1coef(vert,intpt,coef);
    X1 = p1(t1(:,1),:); X2 = p1(t1(:,2),:); X3 = p1(t1(:,3),:); X4 = p1(t1(:,4),:);
    [gx1, gy1, gz1] = gaussPtetra(X1,X2,X3,X4,fem.ng);
    A1 = tetraArea(X1,X2,X3,X4);

    X1 = p2(t2(:,1),:); X2 = p2(t2(:,2),:); X3 = p2(t2(:,3),:); X4 = p2(t2(:,4),:);
    [gx2, gy2, gz2] = gaussPtetra(X1,X2,X3,X4,fem.ng);
    A2 = tetraArea(X1,X2,X3,X4);

    nt1 = size(t1,1); nt2 = size(t2,1);
    bas(id+1:id+nt1,:,1) = repmat(BAS1(1,:),nt1,1);
    bas(id+1:id+nt1,:,2) = repmat(BAS1(2,:),nt1,1);
    bas(id+1:id+nt1,:,3) = repmat(BAS1(3,:),nt1,1);
    bas(id+1:id+nt1,:,4) = repmat(BAS1(4,:),nt1,1);
    bas(id+nt1+1:id+nt1+nt2,:,1) = repmat(BAS2(1,:),nt2,1);
    bas(id+nt1+1:id+nt1+nt2,:,2) = repmat(BAS2(2,:),nt2,1);
    bas(id+nt1+1:id+nt1+nt2,:,3) = repmat(BAS2(3,:),nt2,1);
    bas(id+nt1+1:id+nt1+nt2,:,4) = repmat(BAS2(4,:),nt2,1);
    gx(id+1:id+nt1+nt2,:) = [gx1;gx2];
    gy(id+1:id+nt1+nt2,:) = [gy1;gy2];
    gz(id+1:id+nt1+nt2,:) = [gz1;gz2];
    A(id+1:id+nt1+nt2,:) = [A1;A2];
    QelemID(id+1:id+nt1+nt2) = tID;

    temp = fem.t(tID,:);
    temp2 = temp(locID);
    t(id+1:id+nt1+nt2,:) = repmat(temp2,nt1+nt2,1);
    locIDvec(i,:) = locID;
    quadT(i) = nt1+nt2;

    bas1(i,:,1) = BAS1(1,:); bas2(i,:,1) = BAS2(1,:);
    bas1(i,:,2) = BAS1(2,:); bas2(i,:,2) = BAS2(2,:);
    bas1(i,:,3) = BAS1(3,:); bas2(i,:,3) = BAS2(3,:);
    bas1(i,:,4) = BAS1(4,:); bas2(i,:,4) = BAS2(4,:);

    id = id+nt1+nt2;
end
t(id+1:end,:) = []; bas(id+1:end,:,:) = []; A(id+1:end) = [];
gx(id+1:end,:) = []; gy(id+1:end,:) = []; gz(id+1:end,:) = [];
QelemID(id+1:end,:) = [];
femI = struct('t',t,'bas',bas,'gx',gx,'gy',gy,'gz',gz,'area',A,'gw',fem.gw,...
    'quadT',quadT,'locID',locIDvec,'bas1',bas1,'bas2',bas2,'plusPC',plusPC,...
    'QelemID',QelemID);

repo: ifem


text:
function [ x, c, m ] = dvandprg ( n, alpha, b, x, c, m )

%*****************************************************************************80
%
%% DVANDPRG solves a Vandermonde system A' * x = f progressively.
%
%  Discussion:
%
%    This function receives the solution to the system of equations A' * x = f
%    where A is a Vandermonde matrix for alpha(0) through alpha(n-1),
%    and new values alpha(n) and f(n).  It updates the solution.
%
%    To solve a system of Nbig equations, this function may be called repeatedly,
%    with N = 1, 2, ..., Nbig.  Each time, a solution to the current subsystem
%    is returned.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 April 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Ake Bjorck, Victor Pereyra,
%    Solution of Vandermonde Systems of Equations,
%    Mathematics of Computation,
%    Volume 24, Number 112, October 1970, pages 893-903.
%
%  Parameters:
%
%    Input, integer N, the new order of the matrix, which is 1 larger
%    than on the previous call.  For the first call, N must be 1.
%
%    Input, real ALPHA(N), the parameters that define the matrix.
%    The values should be distinct.  The value ALPHA(N) has just been
%    added to the system.
%
%    Input, real B(N), the right hand side of the linear system.
%
%    Input, real X(N-1), the previous solution of the linear system.
%
%    Input, real C(N-1), real M(N-1), factorization data from the previous call.
%
%    Output, real X(N), the updated solution to the linear system.
%
%    Output, real C(N), M(N), updated factorization data.
%
  c(n) = b(n);
  for j = n - 1 : -1 : 1
    c(j) = ( c(j+1) - c(j) ) / ( alpha(n) - alpha(j) );
  end

  if ( n == 1 )
    m(n) = 1.0;
  else
    m(n) = 0.0;
  end
  cn = c(1);
  x(n) = c(1);

  for j = n - 1 : -1 : 1
    m(j+1) = m(j+1) - alpha(n-1) * m(j);
    x(n-j) = x(n-j) + m(j+1) * cn;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [new_clu log] =renumberclu(clu,log)
% re-orders clu to start from index 1 and iterate through the number of
% clusters

% rewritten by David Tingley, 2017
new_clu = zeros(length(clu),1);

cluster_names = unique(clu);
for i=1:length(cluster_names)
    if cluster_names(i) ~= 0  % keep cluster 0 (noise) as 0
    new_clu(find(clu==cluster_names(i))) = i+2;
     log = [log sprintf('%d -> %d; reordering clusters\n',cluster_names(i),i+2)];
    end
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [err,contrfile] = ContrastTest (N_runs, polort, N_basis, Task, N_tasks, N_contr, El, N_base)
%
%   [err,] = ContrastTest ()
%
%Purpose:
%   Construct the contrast test matrix used in -glt option of 3dDeconvolve.
%
%
%Input Parameters:
%
%
%
%Output Parameters:
%   err : 0 No Problem
%       : 1  Problems
%
%
%
%Key Terms:
%
%More Info :
%
%
%
%
%     Author : Gang Chen
%     Date : Tue Jul 29 13:01:05 EDT 2003
%     SSCC/NIMH/ National Institutes of Health, Bethesda Maryland

%Define the function name for easy referencing
FuncName = 'ContrastTest';

%Debug Flag
DBG = 1;

%initailize return variables
err = 1;

%Create a matrix so that 3dDeconvolve uses for regressor test.

length = (polort+1)*N_runs+N_base;
for (iT = 1:1:N_tasks),
   for (j = 1:1:N_basis)    ,
	   length = length + (Task(iT).BasisOpt_struct(j).maxlag - Task(iT).BasisOpt_struct(j).minlag + 1);
	end
end

%mtrx = zeros(N_contr, (polort+1)*N_runs+N_basis*N_tasks);  %for each run, there are (polort+1) paramters in the baseline model

mtrx = zeros(N_contr, length);
fprintf(1, '\n\nYour contrast vectors for -glt options are:\n');

for (i = 1:1:N_contr),
	for (j = 1:1:El(i).cnt),
	   for (k = 1:1:N_tasks),
		   if El(i).order(j) == k,
			   for (m=1:1:N_basis),
				   for (p = 1:1:Task(k).BasisOpt_struct(m).maxlag - Task(k).BasisOpt_struct(m).minlag + 1)
%			         mtrx(i,m+(k-1)*N_basis+(polort+1)*N_runs) = El(i).sign(j);
                  mtrx(i,p+(m-1)*(Task(k).BasisOpt_struct(m).maxlag - Task(k).BasisOpt_struct(m).minlag + 1)...
					      +(k-1)*(Task(k).BasisOpt_struct(m).maxlag - Task(k).BasisOpt_struct(m).minlag + 1)*N_basis...
					      +(polort+1)*N_runs) = El(i).sign(j);    %Don't ask me how I worked this out. It is so tedious and nasty
							                                        %to housekeep the indices.
							                                        %Don't even remember how I figured it out now
					end
				end
	      end
		end
	end
	contrfile(i).name = sprintf('%s_contr.1D',El(i).labels);
   fid = fopen(contrfile(i).name, 'w');
	fprintf(1,'No. %i contrast vector = ', i);
	for (n = 1:1:length),
		fprintf(1, '%g	', mtrx(i,n));
		fprintf(fid, '%g ', mtrx(i,n));
	end
	fprintf(1, '\n');
	fclose(fid);
	fprintf(1, '\nThe vector for No. %i contrast is stored in file %s. \n\n', i, contrfile(i).name);
end

err = 0;
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% initlist.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [x0,f0,l,L,istar,ncall] = initlist(fcn,data,u,v)
% generates an initialization list with the aid of line searches
%
% Input:
% fcn = 'fun' 	name of function fun(data,x), x an n-vector
% data		data vector
% [u,v]       	box in which the optimization is carried out (u, v
%             	n-vectors)
% Output:
% x0		array with n rows and at least 3 columns; the ith row
%		contains the initialization list values for the ith
%		coordinate
% f0		corresponding function values
% l		x0(i,l(i)) is the ith coordinate of the initial point
% L
% istar		x0(i,istar(i)) is the ith coordinate of the final best
%		point
% ncall		number of function calls used in the program
%
% Uses the following functions/m-files:
% gls.m and its subprograms
%

function [x0,f0,l,L,istar,ncall] = initlist(fcn,data,u,v)
ncall = 0;
nloc = 5;
small = 0.1;
smaxls = 25;
n = length(u);
x = min(max(u,0),v);	% absolutely smallest point
f = feval(fcn,data,x);
ncall = ncall + 1;
for i = 1:n
  alist = 0;
  flist = f;
  p = zeros(n,1);
  p(i) = 1;
  [alist,flist,nfls] = gls(fcn,data,u,v,x,p,alist,flist,nloc,small,smaxls);
  ncall = ncall + nfls;
  [alist1,flist1] = lspost(alist,flist);
  if isempty(find(alist1==0))
    alist1 = [alist1 0];
    flist1 = [flist1 f];
  end
  if length(alist1) < 3
    if isempty(find(alist1==alist(length(alist))))
      alist1 = [alist1 alist(length(alist))];
      flist1 = [flist1 flist(length(alist))];
    end
    if length(alist1) < 3
      if isempty(find(alist1==alist(1)))
        alist1 = [alist1 alist(1)];
        flist1 = [flist1 flist(length(alist))];
      end
      if length(alist1) < 3
        k = round((1+length(alist))/2);
        alist1 = [alist1 alist(k)];
        flist1 = [flist1 flist(k)];
      end
    end
  end
  [alist,ind] = sort(alist1);
  flist = flist1(ind);
  l(i) = find(alist == 0);
  [f1,istar(i)] = min(flist);
  L(i) = length(alist);
  x0(i,1:L(i)) = alist + x(i);
  f0(1:L(i),i) = flist';
  x(i) = x0(i,istar(i));
  f = feval(fcn,data,x);
end

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


text:
function [g,geq,dg,dgeq,xevaled] = fmincon_con_liftlayer(x,model,xevaled)

if isempty(model.lift)
    [g,geq,dg,dgeq,xevaled] = fmincon_con(x,model);
else
    xlift = zeros(length(model.linearindicies),1);
    xlift(model.lift.linearIndex) = x;
    xlift(model.lift.liftedIndex) = model.lift.d + model.lift.T*x;

    % Call the computational kernel which works in the fully expanded
    % normalized format
    [g,geq,dg,dgeq,xevaled] = fmincon_con(xlift,model);

    [f,df,xevaledout] = fmincon_fun(xlift,model);
    %Now map gradient to exposed varaibles to fmincon
    if ~isempty(dg)
        dg = dg(model.lift.linearIndex,:) + model.lift.T'*dg(model.lift.liftedIndex,:);
    end
    if ~isempty(dgeq)
        dgeq = dgeq(model.lift.linearIndex,:) + model.lift.T'*dgeq(model.lift.liftedIndex,:);
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [f,g]=idwilt(c,g,Ls)
%IDWILT  Inverse discrete Wilson transform
%   Usage:  f=idwilt(c,g);
%           f=idwilt(c,g,Ls);
%
%   Input parameters:
%      c     : $2M \times N$ array of coefficients.
%      g     : Window function.
%      Ls    : Final length of function (optional)
%   Output parameters:
%      f     : Input data
%
%   `idwilt(c,g)` computes an inverse discrete Wilson transform with window *g*.
%   The number of channels is deduced from the size of the coefficient array *c*.
%
%   The window *g* may be a vector of numerical values, a text string or a
%   cell array. See the help of |wilwin| for more details.
%
%   `idwilt(f,g,Ls)` does the same, but cuts of zero-extend the final
%   result to length *Ls*.
%
%   `[f,g]=idwilt(...)` additionally outputs the window used in the
%   transform. This is usefull if the window was generated from a
%   description in a string or cell array.
%
%   See also:  dwilt, wilwin, dgt, wilorth
%
%   References: bofegrhl96-1 liva95

%   AUTHOR : Peter L. Søndergaard.
%   TESTING: TEST_DWILT
%   REFERENCE: OK

complainif_argnonotinrange(nargin,2,3,mfilename);

M=size(c,1)/2;
N=2*size(c,2);
W=size(c,3);

a=M;
L=M*N;

assert_L(L,0,L,a,2*M,'IDWILT');

[g,info]=wilwin(g,M,L,'IDWILT');

wasrow=0;
if (ndims(c)==2 && info.wasrow)
  wasrow=1;
end;

f=comp_idwilt(c,g);

% Check if Ls was specified.
if nargin==3
  f=postpad(f,Ls);
else
  Ls=L;
end;

f=comp_sigreshape_post(f,Ls,wasrow,[0; W]);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [unitPos, unitSize] = getFootPrintRez( rez )

%for testing standalone
% load('C:\Users\labadmin\Documents\emouse_drift\nodrift_sn_series\rand_pos_2X_KS2_release\rc03_rezFinal.mat')

%input rez from KS1 or KS2, return array of unit "sizes"
%sites included in unit defined as those with abs(1st component of U) >
% 0.1*max(abs(U(:,1))
% size = diagonal of rectangle that contains the center of the sites, i.e.
% sqrt( (maxY-minY)^2 + (maxX-minX)^2 );

U = rez.U;
[nChan, nClu, nPC] = size(U);
unitPos = zeros(nClu,1);
unitSize = zeros(nClu,1);

nDist = 30;     %neighbor distance, in um
incFrac = 0.1;    %included sites have projection onto first PC this fraction of max

for j = 1:nClu
    currU = U(:,j,1);
    [maxU, maxI] = max(currU);
    [minU, minI] = min(currU);
    maxP = maxU - minU;
    unitPos(j) = rez.yc(maxI);
    inclSite = find(abs(currU) > incFrac*maxP);
    inclDist=[];
    for k = 1:numel(inclSite)
        inclDist(k) = sqrt( (rez.yc(inclSite(k)) - rez.yc(maxI))^2 + ...
                       (rez.xc(inclSite(k)) - rez.xc(maxI))^2 );
    end
    [sortDist, ~] = sort(inclDist);

    nextNearest = 0;
    currI = 1;
    while ((nextNearest < nDist) && currI < numel(inclDist))
        currI = currI + 1;
        nextNearest = sortDist(currI)-sortDist(currI-1);
    end
    if( nextNearest > nDist )
        %then last site tested was not contiguous with other included sites
        currI = currI - 1;
    end
    %fprintf( 'cluIndex, numIncl, currI: %d, %d, %d\n', j, numel(inclSite),currI);
    %unit radius est = sortDist(currI)
    unitSize(j) = 2*sortDist(currI);

end

end

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function net = resnet52_new_hope()

%----------------------img cnn----------------------
netStruct = load('./data/imagenet-vgg-verydeep-19.mat') ;
net = dagnn.DagNN.loadobj(netStruct) ;
net.removeLayer('fc8');
net.removeLayer('prob');
net.removeLayer('relu7');

for i = 1:numel(net.params)
    if(mod(i,2)==0)
        net.params(i).learningRate= 0; %0.02;
    else net.params(i).learningRate= 0; %0.001;
    end
    net.params(i).weightDecay=0; %1;
end

net.params(1).learningRate = 0; %1e-5;

fc1Block = dagnn.Conv('size',[1 1 4096 4096],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc1',fc1Block,{'x40'},{'fc1'},{'fc1f'});
net.addLayer('fc1bn',dagnn.BatchNorm(),{'fc1'},{'fc1bn'},...
    {'fc1bn_w','fc1bn_b','fc1bn_m'});
net.addLayer('fc1x',dagnn.ReLU(),{'fc1bn'},{'fc1bnx'});
fc1Block = dagnn.Conv('size',[1 1 4096 2048],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc1_1',fc1Block,{'fc1bnx'},{'fc1_1'},{'fc1_1f'});
net.addLayer('fc1_1bn',dagnn.BatchNorm(),{'fc1_1'},{'fc1_1bn'},...
    {'fc1_1bn_w','fc1_1bn_b','fc1_1bn_m'});
net.addLayer('fc1_1x',dagnn.ReLU(),{'fc1_1bn'},{'fc1_1bnx'});
net.addLayer('dropout',dagnn.DropOut('rate',0.75),{'fc1_1bnx'},{'fc1_1bnxd'});

%----------------------char cnn----------------------
% input is 1*32*20074*16
fc2Block = dagnn.Conv('size',[1 1 20074 300],'hasBias',true,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc2',fc2Block,{'data2'},{'fc2'},{'fc2f','fc2b'});
net.addLayer('fc2bn',dagnn.BatchNorm(),{'fc2'},{'fc2bn'},...
    {'fc2bn_w','fc2bn_b','fc2bn_m'});
%net.addLayer('fc2x',dagnn.ReLU(),{'fc2bn'},{'fc2bnx'});
%net.addLayer('dropout_diction',dagnn.DropOut('rate',0.5),{'fc2bn'},{'fc2bnd'});
% 32*256
convBlock = dagnn.Conv('size',[1 1 300 128],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc2_1_1',convBlock,{'fc2bn'},{'fc2_1_1'},{'fc2_1_1f'});
net.addLayer('fc2_1_1bn',dagnn.BatchNorm(),{'fc2_1_1'},{'fc2_1_1bn'},...
    {'fc2_1_1bn_w','fc2_1_1bn_b','fc2_1_1bn_m'});
net.addLayer('fc2_1_1x',dagnn.ReLU(),{'fc2_1_1bn'},{'fc2_1_1bnx'});

convBlock = dagnn.Conv('size',[1 2 128 128],'hasBias',false,'stride',[1,1],'pad',[0,0,1,0]);
net.addLayer('fc2_1_2',convBlock,{'fc2_1_1bnx'},{'fc2_1_2'},{'fc2_1_2f'});
net.addLayer('fc2_1_2bn',dagnn.BatchNorm(),{'fc2_1_2'},{'fc2_1_2bn'},...
    {'fc2_1_2bn_w','fc2_1_2bn_b','fc2_1_2bn_m'});
net.addLayer('fc2_1_2x',dagnn.ReLU(),{'fc2_1_2bn'},{'fc2_1_2bnx'});

convBlock = dagnn.Conv('size',[1 1 128 256],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc2_1_3',convBlock,{'fc2_1_2bnx'},{'fc2_1_3'},{'fc2_1_3f'});
net.addLayer('fc2_1_3bn',dagnn.BatchNorm(),{'fc2_1_3'},{'fc2_1_3bn'},...
    {'fc2_1_3bn_w','fc2_1_3bn_b','fc2_1_3bn_m'});

convBlock = dagnn.Conv('size',[1 1 300 256],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc2_1b',convBlock,{'fc2bn'},{'fc2_1b'},{'fc2_1bf'});
net.addLayer('fc2_1bbn',dagnn.BatchNorm(),{'fc2_1b'},{'fc2_1bbn'},...
    {'fc2_1bbn_w','fc2_1bbn_b','fc2_1bbn_m'});

net.addLayer('fc2_1sum',dagnn.Sum(),{'fc2_1_3bn','fc2_1bbn'},...
    {'fc2_1sum'});
net.addLayer('fc2_1x',dagnn.ReLU(),{'fc2_1sum'},{'fc2_1sumx'});

% 32*256
for i  = 2:3
    convBlock = dagnn.Conv('size',[1 1 256 64],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
    net.addLayer( sprintf('fc2_%d_1',i),convBlock,{sprintf('fc2_%dsumx',i-1)},{sprintf('fc2_%d_1',i)}, ...
        {sprintf('fc2_%d_1f',i)});
    net.addLayer(sprintf('fc2_%d_1bn_1',i),dagnn.BatchNorm(),{sprintf('fc2_%d_1',i)},{sprintf('fc2_%d_1bn',i)},...
        {sprintf('fc2_%d_1bn_w',i),sprintf('fc2_%d_1bn_b',i),sprintf('fc2_%d_1bn_m',i)});
    net.addLayer(sprintf('fc2_%d_1x',i),dagnn.ReLU(),{sprintf('fc2_%d_1bn',i)},{sprintf('fc2_%d_1bnx',i)});

    convBlock = dagnn.Conv('size',[1 2 64 64],'hasBias',false,'stride',[1,1],'pad',[0,0,1,0]);
    net.addLayer( sprintf('fc2_%d_2',i),convBlock,{sprintf('fc2_%d_1bnx',i)},{sprintf('fc2_%d_2',i)}, ...
        {sprintf('fc2_%d_2f',i)});
    net.addLayer(sprintf('fc2_%d_2bn',i),dagnn.BatchNorm(),{sprintf('fc2_%d_2',i)},{sprintf('fc2_%d_2bn',i)},...
        {sprintf('fc2_%d_2bn_w',i),sprintf('fc2_%d_2bn_b',i),sprintf('fc2_%d_2bn_m',i)});
    net.addLayer(sprintf('fc2_%d_2x',i),dagnn.ReLU(),{sprintf('fc2_%d_2bn',i)},{sprintf('fc2_%d_2bnx',i)});

    convBlock = dagnn.Conv('size',[1 1 64 256],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
    net.addLayer( sprintf('fc2_%d_3',i),convBlock,{sprintf('fc2_%d_2bnx',i)},{sprintf('fc2_%d_3',i)}, ...
        {sprintf('fc2_%d_3f',i)});
    net.addLayer(sprintf('fc2_%d_3bn',i),dagnn.BatchNorm(),{sprintf('fc2_%d_3',i)},{sprintf('fc2_%d_3bn',i)},...
        {sprintf('fc2_%d_3bn_w',i),sprintf('fc2_%d_3bn_b',i),sprintf('fc2_%d_3bn_m',i)});

    net.addLayer(sprintf('fc2_%dsum',i),dagnn.Sum(),{sprintf('fc2_%dsumx',i-1),sprintf('fc2_%d_3bn',i)},...
        {sprintf('fc2_%dsum',i)});
    net.addLayer(sprintf('fc2_%dx',i),dagnn.ReLU(),{sprintf('fc2_%dsum',i)},{sprintf('fc2_%dsumx',i)});
end

%32*256
convBlock = dagnn.Conv('size',[1 1 256 512],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc2_4a_1',convBlock,{'fc2_3sumx'},{'fc2_4a_1'},{'fc2_4a_1f'});
net.addLayer('fc2_4a_1bn',dagnn.BatchNorm(),{'fc2_4a_1'},{'fc2_4a_1bn'},...
    {'fc2_4a_1bn_w','fc2_4a_1bn_b','fc2_4a_1bn_m'});
net.addLayer('fc2_4a_1x',dagnn.ReLU(),{'fc2_4a_1bn'},{'fc2_4a_1bnx'});

convBlock = dagnn.Conv('size',[1 2 512 512],'hasBias',false,'stride',[2,2],'pad',[0,0,1,0]);
net.addLayer('fc2_4a_2',convBlock,{'fc2_4a_1bnx'},{'fc2_4a_2'},{'fc2_4a_2f'});
net.addLayer('fc2_4a_2bn',dagnn.BatchNorm(),{'fc2_4a_2'},{'fc2_4a_2bn'},...
    {'fc2_4a_2bn_w','fc2_4a_2bn_b','fc2_4a_2bn_m'});
net.addLayer('fc2_4a_2x',dagnn.ReLU(),{'fc2_4a_2bn'},{'fc2_4a_2bnx'});

convBlock = dagnn.Conv('size',[1 1 512 512],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc2_4a_3',convBlock,{'fc2_4a_2bnx'},{'fc2_4a_3'},{'fc2_4a_3f'});
net.addLayer('fc2_4a_3bn',dagnn.BatchNorm(),{'fc2_4a_3'},{'fc2_4a_3bn'},...
    {'fc2_4a_3bn_w','fc2_4a_3bn_b','fc2_4a_3bn_m'});

convBlock = dagnn.Conv('size',[1 1 256 512],'hasBias',false,'stride',[2,2],'pad',[0,0,0,0]);
net.addLayer('fc2_4b',convBlock,{'fc2_3sumx'},{'fc2_4b'},{'fc2_4bf'});
net.addLayer('fc2_4bbn',dagnn.BatchNorm(),{'fc2_4b'},{'fc2_4bbn'},...
    {'fc2_4bbn_w','fc2_4bbn_b','fc2_4bbn_m'});

%16*512
net.addLayer('fc2_4sum',dagnn.Sum(),{'fc2_4a_3bn','fc2_4bbn'},...
    {'fc2_4sum'});
net.addLayer('fc2_4x',dagnn.ReLU(),{'fc2_4sum'},{'fc3_1sumx'});

for i  = 2:4
    convBlock = dagnn.Conv('size',[1 1 512 128],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
    net.addLayer( sprintf('fc3_%d_1',i),convBlock,{sprintf('fc3_%dsumx',i-1)},{sprintf('fc3_%d_1',i)}, ...
        {sprintf('fc3_%d_1f',i)});
    net.addLayer(sprintf('fc3_%d_1bn_1',i),dagnn.BatchNorm(),{sprintf('fc3_%d_1',i)},{sprintf('fc3_%d_1bn',i)},...
        {sprintf('fc3_%d_1bn_w',i),sprintf('fc3_%d_1bn_b',i),sprintf('fc3_%d_1bn_m',i)});
    net.addLayer(sprintf('fc3_%d_1x',i),dagnn.ReLU(),{sprintf('fc3_%d_1bn',i)},{sprintf('fc3_%d_1bnx',i)});

    convBlock = dagnn.Conv('size',[1 2 128 128],'hasBias',false,'stride',[1,1],'pad',[0,0,1,0]);
    net.addLayer( sprintf('fc3_%d_2',i),convBlock,{sprintf('fc3_%d_1bnx',i)},{sprintf('fc3_%d_2',i)}, ...
        {sprintf('fc3_%d_2f',i)});
    net.addLayer(sprintf('fc3_%d_2bn',i),dagnn.BatchNorm(),{sprintf('fc3_%d_2',i)},{sprintf('fc3_%d_2bn',i)},...
        {sprintf('fc3_%d_2bn_w',i),sprintf('fc3_%d_2bn_b',i),sprintf('fc3_%d_2bn_m',i)});
    net.addLayer(sprintf('fc3_%d_2x',i),dagnn.ReLU(),{sprintf('fc3_%d_2bn',i)},{sprintf('fc3_%d_2bnx',i)});

    convBlock = dagnn.Conv('size',[1 1 128 512],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
    net.addLayer( sprintf('fc3_%d_3',i),convBlock,{sprintf('fc3_%d_2bnx',i)},{sprintf('fc3_%d_3',i)}, ...
        {sprintf('fc3_%d_3f',i)});
    net.addLayer(sprintf('fc3_%d_3bn',i),dagnn.BatchNorm(),{sprintf('fc3_%d_3',i)},{sprintf('fc3_%d_3bn',i)},...
        {sprintf('fc3_%d_3bn_w',i),sprintf('fc3_%d_3bn_b',i),sprintf('fc3_%d_3bn_m',i)});

    net.addLayer(sprintf('fc3_%dsum',i),dagnn.Sum(),{sprintf('fc3_%dsumx',i-1),sprintf('fc3_%d_3bn',i)},...
        {sprintf('fc3_%dsum',i)});
    net.addLayer(sprintf('fc3_%dx',i),dagnn.ReLU(),{sprintf('fc3_%dsum',i)},{sprintf('fc3_%dsumx',i)});
end

convBlock = dagnn.Conv('size',[1 1 512 1024],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc3_5a_1',convBlock,{'fc3_4sumx'},{'fc3_5a_1'},{'fc3_5a_1f'});
net.addLayer('fc3_5a_1bn',dagnn.BatchNorm(),{'fc3_5a_1'},{'fc3_5a_1bn'},...
    {'fc3_5a_1bn_w','fc3_5a_1bn_b','fc3_5a_1bn_m'});
net.addLayer('fc3_5a_1x',dagnn.ReLU(),{'fc3_5a_1bn'},{'fc3_5a_1bnx'});

convBlock = dagnn.Conv('size',[1 2 1024 1024],'hasBias',false,'stride',[2,2],'pad',[0,0,1,0]);
net.addLayer('fc3_5a_2',convBlock,{'fc3_5a_1bnx'},{'fc3_5a_2'},{'fc3_5a_2f'});
net.addLayer('fc3_5a_2bn',dagnn.BatchNorm(),{'fc3_5a_2'},{'fc3_5a_2bn'},...
    {'fc3_5a_2bn_w','fc3_5a_2bn_b','fc3_5a_2bn_m'});
net.addLayer('fc3_5a_2x',dagnn.ReLU(),{'fc3_5a_2bn'},{'fc3_5a_2bnx'});

convBlock = dagnn.Conv('size',[1 1 1024 1024],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc3_5a_3',convBlock,{'fc3_5a_2bnx'},{'fc3_5a_3'},{'fc3_5a_3f'});
net.addLayer('fc3_5a_3bn',dagnn.BatchNorm(),{'fc3_5a_3'},{'fc3_5a_3bn'},...
    {'fc3_5a_3bn_w','fc3_5a_3bn_b','fc3_5a_3bn_m'});

convBlock = dagnn.Conv('size',[1 1 512 1024],'hasBias',false,'stride',[2,2],'pad',[0,0,0,0]);
net.addLayer('fc3_5b',convBlock,{'fc3_4sumx'},{'fc3_5b'},{'fc3_5bf'});
net.addLayer('fc3_5bbn',dagnn.BatchNorm(),{'fc3_5b'},{'fc3_5bbn'},...
    {'fc3_5bbn_w','fc3_5bbn_b','fc3_5bbn_m'});

%8*1024
net.addLayer('fc3_5sum',dagnn.Sum(),{'fc3_5a_3bn','fc3_5bbn'},...
    {'fc3_5sum'});
net.addLayer('fc3_5x',dagnn.ReLU(),{'fc3_5sum'},{'fc4_1sumx'});

for i  = 2:6
    convBlock = dagnn.Conv('size',[1 1 1024 256],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
    net.addLayer( sprintf('fc4_%d_1',i),convBlock,{sprintf('fc4_%dsumx',i-1)},{sprintf('fc4_%d_1',i)}, ...
        {sprintf('fc4_%d_1f',i)});
    net.addLayer(sprintf('fc4_%d_1bn_1',i),dagnn.BatchNorm(),{sprintf('fc4_%d_1',i)},{sprintf('fc4_%d_1bn',i)},...
        {sprintf('fc4_%d_1bn_w',i),sprintf('fc4_%d_1bn_b',i),sprintf('fc4_%d_1bn_m',i)});
    net.addLayer(sprintf('fc4_%d_1x',i),dagnn.ReLU(),{sprintf('fc4_%d_1bn',i)},{sprintf('fc4_%d_1bnx',i)});

    convBlock = dagnn.Conv('size',[1 2 256 256],'hasBias',false,'stride',[1,1],'pad',[0,0,1,0]);
    net.addLayer( sprintf('fc4_%d_2',i),convBlock,{sprintf('fc4_%d_1bnx',i)},{sprintf('fc4_%d_2',i)}, ...
        {sprintf('fc4_%d_2f',i)});
    net.addLayer(sprintf('fc4_%d_2bn',i),dagnn.BatchNorm(),{sprintf('fc4_%d_2',i)},{sprintf('fc4_%d_2bn',i)},...
        {sprintf('fc4_%d_2bn_w',i),sprintf('fc4_%d_2bn_b',i),sprintf('fc4_%d_2bn_m',i)});
    net.addLayer(sprintf('fc4_%d_2x',i),dagnn.ReLU(),{sprintf('fc4_%d_2bn',i)},{sprintf('fc4_%d_2bnx',i)});

    convBlock = dagnn.Conv('size',[1 1 256 1024],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
    net.addLayer( sprintf('fc4_%d_3',i),convBlock,{sprintf('fc4_%d_2bnx',i)},{sprintf('fc4_%d_3',i)}, ...
        {sprintf('fc4_%d_3f',i)});
    net.addLayer(sprintf('fc4_%d_3bn',i),dagnn.BatchNorm(),{sprintf('fc4_%d_3',i)},{sprintf('fc4_%d_3bn',i)},...
        {sprintf('fc4_%d_3bn_w',i),sprintf('fc4_%d_3bn_b',i),sprintf('fc4_%d_3bn_m',i)});

    net.addLayer(sprintf('fc4_%dsum',i),dagnn.Sum(),{sprintf('fc4_%dsumx',i-1),sprintf('fc4_%d_3bn',i)},...
        {sprintf('fc4_%dsum',i)});
    net.addLayer(sprintf('fc4_%dx',i),dagnn.ReLU(),{sprintf('fc4_%dsum',i)},{sprintf('fc4_%dsumx',i)});
end

convBlock = dagnn.Conv('size',[1 1 1024 2048],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc4_7a_1',convBlock,{'fc4_6sumx'},{'fc4_7a_1'},{'fc4_7a_1f'});
net.addLayer('fc4_7a_1bn',dagnn.BatchNorm(),{'fc4_7a_1'},{'fc4_7a_1bn'},...
    {'fc4_7a_1bn_w','fc4_7a_1bn_b','fc4_7a_1bn_m'});
net.addLayer('fc4_7a_1x',dagnn.ReLU(),{'fc4_7a_1bn'},{'fc4_7a_1bnx'});

convBlock = dagnn.Conv('size',[1 2 2048 2048],'hasBias',false,'stride',[1,1],'pad',[0,0,1,0]);
net.addLayer('fc4_7a_2',convBlock,{'fc4_7a_1bnx'},{'fc4_7a_2'},{'fc4_7a_2f'});
net.addLayer('fc4_7a_2bn',dagnn.BatchNorm(),{'fc4_7a_2'},{'fc4_7a_2bn'},...
    {'fc4_7a_2bn_w','fc4_7a_2bn_b','fc4_7a_2bn_m'});
net.addLayer('fc4_7a_2x',dagnn.ReLU(),{'fc4_7a_2bn'},{'fc4_7a_2bnx'});

convBlock = dagnn.Conv('size',[1 1 2048 2048],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc4_7a_3',convBlock,{'fc4_7a_2bnx'},{'fc4_7a_3'},{'fc4_7a_3f'});
net.addLayer('fc4_7a_3bn',dagnn.BatchNorm(),{'fc4_7a_3'},{'fc4_7a_3bn'},...
    {'fc4_7a_3bn_w','fc4_7a_3bn_b','fc4_7a_3bn_m'});

convBlock = dagnn.Conv('size',[1 1 1024 2048],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc4_7b',convBlock,{'fc4_6sumx'},{'fc4_7b'},{'fc4_7bf'});
net.addLayer('fc4_7bbn',dagnn.BatchNorm(),{'fc4_7b'},{'fc4_7bbn'},...
    {'fc4_7bbn_w','fc4_7bbn_b','fc4_7bbn_m'});

%8*2048
net.addLayer('fc4_7sum',dagnn.Sum(),{'fc4_7a_3bn','fc4_7bbn'},...
    {'fc4_7sum'});
net.addLayer('fc4_7x',dagnn.ReLU(),{'fc4_7sum'},{'fc5_1sumx'});

poolBlock = dagnn.Pooling('poolSize',[1 8]);
net.addLayer('fc5_1',poolBlock,{'fc5_1sumx'},{'fc5_1bnx'});
%net.addLayer('fc5_1bn',dagnn.BatchNorm(),{'fc5_1'},{'fc5_1bn'},...
 %   {'fc5_1bn_w','fc5_1bn_b','fc5_1bn_m'});
%net.addLayer('fc5_1x',dagnn.ReLU(),{'fc5_1bn'},{'fc5_1bnx'});

fc5_2Block = dagnn.Conv('size',[1 1 2048 2048],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc5_2',fc5_2Block,{'fc5_1bnx'},{'fc5_2'},{'fc5_2f'});
net.addLayer('fc5_2bn',dagnn.BatchNorm(),{'fc5_2'},{'fc5_2bn'},...
    {'fc5_2bn_w','fc5_2bn_b','fc5_2bn_m'});
net.addLayer('fc5_2x',dagnn.ReLU(),{'fc5_2bn'},{'fc5_2bnx'});
net.addLayer('dropout2',dagnn.DropOut('rate',0.75),{'fc5_2bnx'},{'fc5_2bnxd'});

%----------------------add share layer----------------------
%1

fc_imgBlock = dagnn.Conv('size',[1 1 2048 29783],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc_img',fc_imgBlock,{'fc1_1bnxd'},{'prediction_img'},{'fcsharef'});
net.addLayer('softmaxloss_img',dagnn.Loss('loss','softmaxlog'),{'prediction_img','label_img'},'objective_img');
net.addLayer('top1err_img', dagnn.Loss('loss', 'classerror'), ...
    {'prediction_img','label_img'}, 'top1err_img') ;
net.addLayer('top5err_img', dagnn.Loss('loss', 'topkerror', ...
    'opts', {'topK',5}), ...
    {'prediction_img','label_img'}, 'top5err_img') ;
%2
fc_txtBlock = dagnn.Conv('size',[1 1 2048 29783],'hasBias',false,'stride',[1,1],'pad',[0,0,0,0]);
net.addLayer('fc_txt',fc_txtBlock,{'fc5_2bnxd'},{'prediction_txt'},{'fcsharef'});
net.addLayer('softmaxloss_txt',dagnn.Loss('loss','softmaxlog'),{'prediction_txt','label_txt'},'objective_txt');
net.addLayer('top1err_txt', dagnn.Loss('loss', 'classerror'), ...
    {'prediction_txt','label_txt'}, 'top1err_txt') ;
net.addLayer('top5err_txt', dagnn.Loss('loss', 'topkerror', ...
    'opts', {'topK',5}), ...
    {'prediction_txt','label_txt'}, 'top5err_txt') ;
%}

%}
net.initParams();
%----------NOTICE--------------
first = net.getParamIndex('fc2f');

net.params(first).learningRate = 1e-3;  %w
net.params(first+1).learningRate = 1e-3;  %b

load('./flickr30k_dictionary.mat');
%m = mean(subset.features,2);
%subset.features = subset.features-repmat(m,1,20074);
net.params(first).value = reshape(single(subset.features'),1,1,20074,300);

%net.conserveMemory = false;
%net.move('gpu') ;
%net.eval({'data',gpuArray(single(rand(224,224,3,32))),'data2',gpuArray(single(rand(1,32,20074,32)))});
end

author: layumi

repo: Image-Text-Embedding


text:
function lambda = herndon_eigenvalues ( n )

%*****************************************************************************80
%
%% HERNDON_EIGENVALUES returns the eigenvalues of the HERNDON matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real LAMBDA(N,1), the eigenvalues.
%
  lambda = zeros ( n, 1 );

  p = 3.0 + sqrt ( ( 4 * n - 3 ) * ( n - 1 ) * 3 ) / ( n + 1 ) );

  lambda(1:n-2,1) = 1.0;
  lambda(n-1,1) = 6.0 / ( p * ( n + 1 ) );
  lambda(n,1) = p / ( n * ( 5 - 2 * n ) );

  return
end


text:
function r = acosh(a)
%ACOSH        Hessian (elementwise) inverse hyperbolic cosine
%

% written  04/04/04     S.M. Rump
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  K = prod(size(a.x));
  if K==1                   % scalar hessian

    r.x = acosh(a.x);
    f = sqr(a.x) - 1;
    fs = 1 / ( sqrt(a.x+1)*sqrt(a.x-1) );
    r.dx = a.dx * fs;
    r.hx = a.hx * fs - reshape( ( (0.5*a.x/f)*r.dx ) * a.dx.' , size(a.hx) );

  else                      % matrix hessian

    N = getappdata(0,'INTLAB_HESSIAN_NUMVAR');
    N2 = N^2;

    r.x = acosh(a.x);
    if issparse(a.hx)               % input sparse

      ax = full(a.x(:));
      f = sqr(ax) - 1;
      ax = 1 ./ ( sqrt(ax+1).*sqrt(ax-1) );
      sizeax = length(ax);
      [ia,ja,sa] = find(a.dx);
      % check for emptyness: cures Matlab bug
      % a=sparse([],[],[],2,1), [i,j,s]=find(a), s(i).*s(:)  yields error
      if isempty(ia)
        r.dx = sparse([],[],[],N,sizeax);
        r.hx = sparse([],[],[],N2,sizeax);
      else
        adx1 = ( -0.5*full(a.x(:)) ) ./ f;
        if isa(a.x,'intval')          % sparse intval
          rdx = times(ax(ja),sa(:),0);
          adx1 = times(adx1(ja),sa(:),0);
          if rdx.complex
            r.dx = intval( sparse(ia,ja,rdx.mid,N,sizeax) , sparse(ia,ja,rdx.rad,N,sizeax) , 'midrad' );
          else
            r.dx = intval( sparse(ia,ja,rdx.inf,N,sizeax) , sparse(ia,ja,rdx.sup,N,sizeax) , 'infsup' );
          end
          if adx1.complex
            adx1 = intval( sparse(ia,ja,adx1.mid,N,sizeax) , sparse(ia,ja,adx1.rad,N,sizeax) , 'midrad' );
          else
            adx1 = intval( sparse(ia,ja,adx1.inf,N,sizeax) , sparse(ia,ja,adx1.sup,N,sizeax) , 'infsup' );
          end
        else                          % sparse point
          r.dx = sparse(ia,ja,ax(ja).*sa(:),N,sizeax);
          adx1 = sparse(ia,ja,adx1(ja).*sa(:),N,sizeax);
        end
        r.hx = adx2rhx(N,sizeax,adx1,r.dx);
      end
      [ia,ja,sa] = find(a.hx);      % sparse point or intval
      % check for emptyness: cures Matlab bug
      % a=sparse([],[],[],2,1), [i,j,s]=find(a), s(i).*s(:)  yields error
      if ~isempty(ia)
        if isa(a.x,'intval')
          rhx = times(ax(ja),sa(:),0);
          if rhx.complex
            r.hx = r.hx + intval( sparse(ia,ja,rhx.mid,N2,sizeax) , sparse(ia,ja,rhx.rad,N2,sizeax) , 'midrad' );
          else
            r.hx = r.hx + intval( sparse(ia,ja,rhx.inf,N2,sizeax) , sparse(ia,ja,rhx.sup,N2,sizeax) , 'infsup' );
          end
        else
          r.hx = r.hx + sparse(ia,ja,ax(ja).*sa(:),N2,sizeax);
        end
      end

    else                            % input full

      r.x = acosh(a.x);
      ax = a.x(:).';
      f = sqr(ax) - 1;
      fs = 1 ./ ( sqrt(ax+1).*sqrt(ax-1) );
      fs = fs(ones(N*N,1),:);
      r.dx = a.dx .* fs(1:N,:);
      adx = repmat(0.5*ax./f,N,1) .* r.dx;
      r.hx = a.hx .* fs - adx(repmat(1:N,N,1),:) .* a.dx(repmat(1:N,1,N),:);

    end

  end

  r = class(r,'hessian');

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
% SYNTAX:
%    flag = flagShrink(flag, expand_size, <omit_borders = false>)
%
% DESCRIPTION:
%    shrink the flag array
%    if flag is a matrix this flagging expansion will work column by column
%
% INPUT:
%   flag          [n_obs x n_arrays]
%   expand_size   n_epochs with flags to activate at the border of a flagged interval
%   omit_borders  when true do not shrink flags at the beginning or end of the matrix

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:       Andrea Gatti
%  Contributors:     Andrea Gatti
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   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/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

function flag = flagShrink(flag, expand_size, omit_borders)
    if nargin < 3
        omit_borders = false;
    end
    for c = 1 : size(flag, 2)
        flag(:, c) = ~conv(single(~flag(:, c)), ones(2 * expand_size + 1, 1)', 'same') > 0;
    end
    if omit_borders
        % do nothing
    else
        flag([1:expand_size (end - expand_size + 1) : end], :) = 0;
    end
end

author: goGPS-Project


text:
function planC = chopStructure(structureNum,minZ,maxZ,planC)
%function planC = chopStructure(structureNum,minZ,maxZ,planC)
%
%This function creates a new function which lies within minZ and maxZ.
%Use-case: To chop Skin structure so that dose is computed only around the prostate
%region.
%
%Usage: planC = chopStructure(2,-0.2,12.7)
%
%APA, 04/28/2009
%
% 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/>.

global stateS

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

%Get associated scan number
scanNum = getStructureAssociatedScan(structureNum);

[xVals, yVals, zVals] = getScanXYZVals(planC{indexS.scan}(scanNum));

minSlice = max(find(zVals <= minZ));
maxSlice = min(find(zVals >= maxZ));

newStr = newCERRStructure(scanNum,planC);

numStructs = length(planC{indexS.structures});
toAdd = numStructs + 1;
for i=1:length(zVals)
    newStr.contour(i).segments.points = [];
end
newStr.contour(minSlice:maxSlice) = planC{indexS.structures}(structureNum).contour(minSlice:maxSlice);
newStr.structureName = [planC{indexS.structures}(structureNum).structureName, '_chopped'];
newStr.strUID = createUID('structure');
newStr.assocScanUID = planC{indexS.scan}(scanNum).scanUID;
newStr.visible = 1;
planC{indexS.structures} = dissimilarInsert(planC{indexS.structures}, newStr, toAdd);
planC = getRasterSegs(planC, toAdd);

planC = updateStructureMatrices(planC, toAdd);

stateS.structsChanged = 1;
sliceCallBack('refresh');

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% ------------------------------------------------------------------------
%      Analyze non-linear ROIs that were clicked as 'object points'
% ------------------------------------------------------------------------

%% ENTER PARAMETERS AND FILE LOCATION

% file location of object points
save_folder = 'C:\Drive\Histology\for tutorial\Richards\processed';

% directory of reference atlas files
annotation_volume_location = 'C:\Drive\Histology\annotation_volume_10um_by_index.npy'; % from the allen inst (see readme)
structure_tree_location = 'C:\Drive\Histology\structure_tree_safe_2017.csv'; % located in github repo
CCF_to_FP_location =  'C:\Drive\Histology\CCF_to_FP.csv'; % located in github repo
FP_table_location = 'C:\Drive\Histology\FP_table_Chon_2020.csv'; % located in github repo
chon_images_loc = 'C:\Drive\Histology\Suppl_File1_Labels'; % from chon et al (supplementary data 4, https://www.nature.com/articles/s41467-019-13057-w)

% name of the saved object points
object_save_name_suffix = '';

% either set to 'all' or a list of indices from the clicked objects in this file, e.g. [2,3]
objects_to_analyze = 'all';

% plane used to view when points were clicked ('coronal' -- most common, 'sagittal', 'transverse')
plane = 'coronal';

% brain figure black or white
black_brain = true;

%% LOAD THE REFERENCE ANNOTATIONS AND PROBE POINTS

% load the reference brain annotations
if ~exist('av','var') || ~exist('st','var')
    disp('loading reference atlas...')
    av = readNPY(annotation_volume_location);
    st = loadStructureTree(structure_tree_location);
end
if ~exist('CCFtoFPtable','var') || ~exist('FPtable','var')
    CCFtoFPtable = loadCCFtoFP(CCF_to_FP_location);
    FPtable = loadFPtable(FP_table_location);
end

% load object points
objectPoints = load(fullfile(save_folder, ['probe_points' object_save_name_suffix]));

% determine which objects to analyze
if strcmp(objects_to_analyze,'all')
    objects = 1:size(objectPoints.pointList.pointList,1);
else
    objects = objects_to_analyze;
end

%% BRING UP THE RELEVANT DATA FOR EACH PROBE POINTS, FOR FURTHER ANALYSIS

% initialize cell array containing info on each clicked point
if length(objects) > 1
    roi_annotation = cell(length(objects),1);
    roi_location = cell(length(objects),1);
end

% generate needed values
bregma = allenCCFbregma(); % bregma position in reference data space
atlas_resolution = 0.010; % mm

% plot brain grid
ProbeColors = [1 1 1; 1 .75 0;  .3 1 1; .4 .6 .2; 1 .35 .65; .7 .7 1; .65 .4 .25; .7 .95 .3; .7 0 0; .6 0 .7; 1 .6 0];
% order of colors: {'white','gold','turquoise','fern','bubble gum','overcast sky','rawhide', 'green apple','purple','orange','red'};
fwireframe = plotBrainGrid([], [], [], black_brain); hold on;
fwireframe.InvertHardcopy = 'off';

for object_num = objects

    selected_object = objects(object_num);

    % get the object points for the currently analyzed object
    if strcmp(plane,'coronal')
        curr_objectPoints = objectPoints.pointList.pointList{selected_object,1}(:, [3 2 1]);
    elseif strcmp(plane,'sagittal')
        curr_objectPoints = objectPoints.pointList.pointList{selected_object,1}(:, [1 2 3]);
    elseif strcmp(plane,'transverse')
        curr_objectPoints = objectPoints.pointList.pointList{selected_object,1}(:, [1 3 2]);
    end

    % plot points on the wire frame brain
    figure(fwireframe); hold on
    hp = plot3(curr_objectPoints(:,1), curr_objectPoints(:,3), curr_objectPoints(:,2), '.','linewidth',2, 'color',[ProbeColors(object_num,:) .2],'markers',10);

    % use the point's position in the atlas to get the AP, DV, and ML coordinates
    ap = -(curr_objectPoints(:,1)-bregma(1))*atlas_resolution;
    dv = (curr_objectPoints(:,2)-bregma(2))*atlas_resolution;
    ml = (curr_objectPoints(:,3)-bregma(3))*atlas_resolution;

    roi_location_curr = [ap dv ml];

    % initialize array of region annotations
    roi_annotation_curr = cell(size(curr_objectPoints,1),3);
    roi_annotation_curr_FP = cell(size(curr_objectPoints,1),3);

    % loop through every point to get ROI locations and region annotations
    for point = 1:size(curr_objectPoints,1)

        % find the annotation, name, and acronym of the current ROI pixel
        ann = av(curr_objectPoints(point,1),curr_objectPoints(point,2),curr_objectPoints(point,3));
        name = st.safe_name{ann};
        acr = st.acronym{ann};

        roi_annotation_curr{point,1} = ann;
        roi_annotation_curr{point,2} = name;
        roi_annotation_curr{point,3} = acr;

        % find the annotation, name, and acronym of the current ROI pixel
        [ann_FP, name_FP, acr_FP] = CCF_to_FP(curr_objectPoints(point,1), curr_objectPoints(point,2), curr_objectPoints(point,3), CCFtoFPtable, FPtable, chon_images_loc);

        roi_annotation_curr_FP{point,1} = ann_FP;
        roi_annotation_curr_FP{point,2} = name_FP;
        roi_annotation_curr_FP{point,3} = acr_FP;

    end

    % save results in cell array
    if length(objects) > 1
        roi_annotation{object_num} = roi_annotation_curr;
        roi_location{object_num} = roi_location_curr;
    else
        roi_annotation = roi_annotation_curr;
        roi_location = roi_location_curr;
    end

    % display results in a table
    disp(['Clicked points for object ' num2str(selected_object)])
    roi_table = table(roi_annotation_curr(:,2),roi_annotation_curr(:,3), roi_annotation_curr_FP(:,2),roi_annotation_curr_FP(:,3),...
                        roi_location_curr(:,1),roi_location_curr(:,2),roi_location_curr(:,3), roi_annotation_curr(:,1), roi_annotation_curr_FP(:,1),...
         'VariableNames', {'CCF_name', 'CCF_acronym', 'FP_name', 'FP_acronym', 'AP_location', 'DV_location', 'ML_location', 'avIndex', 'fpIndex'});
     disp(roi_table)

end

% now, use roi_location and roi_annotation for your further analyses

author: cortex-lab

repo: allenCCF

save_path: github-repos/MATLAB/cortex-lab-allenCCF


%% pathRRT
%%  - create a path from a start node to an end node
%%    using the RRT algorithm.
%%  - RRT = Rapidly-exploring Random Tree
%%
%% Last Modified - 6/8/2006 - R. Beard
%%               - 4/15/2010 - R. Beard

function path_out=planRRT(wpp_start, wpp_end, map)

    % standard length of path segments
    segmentLength = 100;

    % desired down position is down position of end node
    pd = wpp_end(3);
    chi = -9999;

    % specify start and end nodes from wpp_start and wpp_end
    start_node = [wpp_start(1), wpp_start(2), pd, chi, 0, 0, 0];
    end_node = [wpp_end(1), wpp_end(2), pd, chi, 0, 0, 0];
    % format:  [N, E, D, chi, cost, parent_idx, flag_connect_to_goal]

    % establish tree starting with the start node
    tree = start_node;

    % check to see if start_node connects directly to end_node
    if ( (norm(start_node(1:3)-end_node(1:3))<segmentLength )...
            &(collision(start_node,end_node,map)==0) )
        path = [start_node; end_node];
    else
        numPaths = 0;
        while numPaths<3,
            [tree,flag] = extendTree(tree,end_node,segmentLength,map,pd,chi);
            numPaths = numPaths + flag;
        end
    end

    % find path with minimum cost to end_node
    path = findMinimumPath(tree,end_node);
    path_out = smoothPath(path,map);
    plotmap(map,path,path_out,tree);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% generateRandomNode
%%   create a random node (initialize)
function node=generateRandomNode(map,pd,chi)

    % randomly pick configuration
    pn       = map.width*rand;
    pe       = map.width*rand;
    pd       = pd; % constant altitute paths
    cost     = 0;
    node     = [pn, pe, pd, chi, cost, 0, 0];
    % format:  [N, E, D, chi, cost, parent_idx, flag_connect_to_goal]

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% collision
%%   check to see if a node is in collsion with obstacles
function collision_flag = collision(start_node, end_node, map)

    collision_flag = 0;

    [X,Y,Z] = pointsAlongPath(start_node, end_node, 0.1);

    for i = 1:length(X),
        if Z(i) >= downAtNE(map, X(i), Y(i)),
            collision_flag = 1;
        end
    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% pointsAlongPath
%%   Find points along straight-line path separted by Del (to be used in
%%   collision detection)
function [X,Y,Z] = pointsAlongPath(start_node, end_node, Del)

    X = [start_node(1)];
    Y = [start_node(2)];
    Z = [start_node(3)];

    q = [end_node(1:3)-start_node(1:3)];
    L = norm(q);
    q = q/L;

    w = start_node(1:3);
    for i=2:floor(L/Del),
        w = w + Del*q;
        X = [X, w(1)];
        Y = [Y, w(2)];
        Z = [Z, w(3)];
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% downAtNE
%%   find the world down coordinate at a specified (n,e) location
function down = downAtNE(map, n, e)

      [d_n,idx_n] = min(abs(n - map.buildings_n));
      [d_e,idx_e] = min(abs(e - map.buildings_e));

      if (d_n<=map.BuildingWidth) && (d_e<=map.BuildingWidth),
          down = -map.heights(idx_e,idx_n);
      else
          down = 0;
      end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% extendTree
%%   extend tree by randomly selecting point and growing tree toward that
%%   point
function [new_tree,flag] = extendTree(tree,end_node,segmentLength,map,pd,chi)

  flag1 = 0;
  loop_count = 0;
  while flag1==0,
    % select a random point
    randomNode=generateRandomNode(map,pd,chi);

    % find leaf on node that is closest to randomPoint
    tmp = tree(:,1:3)-ones(size(tree,1),1)*randomNode(1:3);
    [dist,idx] = min(diag(tmp*tmp'));
    L = min(sqrt(dist), segmentLength);
    cost     = tree(idx,5) + L;
    tmp = randomNode(1:3)-tree(idx,1:3);
    new_point = tree(idx,1:3)+L*(tmp/norm(tmp));
    new_node = [new_point, chi, cost, idx, 0];
    loop_count = loop_count+1
    if collision(tree(idx,:), new_node, map)==0,
      new_tree = [tree; new_node];
      flag1=1;
    end
  end

  % check to see if new node connects directly to end_node
  if ( (norm(new_node(1:3)-end_node(1:3))<segmentLength )...
      &&(collision(new_node,end_node,map)==0) )
    flag = 1;
    new_tree(end,7)=1;  % mark node as connecting to end.
  else
    flag = 0;
  end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% findMinimumPath
%%   find the lowest cost path to the end node
function path = findMinimumPath(tree,end_node)

    % find nodes that connect to end_node
    connectingNodes = [];
    for i=1:size(tree,1),
        if tree(i,7)==1,
            connectingNodes = [connectingNodes; tree(i,:)];
        end
    end

    % find minimum cost last node
    [tmp,idx] = min(connectingNodes(:,5));

    % construct lowest cost path
    path = [connectingNodes(idx,:); end_node];
    parent_node = connectingNodes(idx,6);
    while parent_node>1,
        parent_node = tree(parent_node,6);
        path = [tree(parent_node,:); path];
    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% smoothPath
%%   smooth the waypoint path
function newPath = smoothPath(path,map)

    newPath = path(1,:); % add the start node
    ptr =2;  % pointer into the path
    while ptr <= size(path,1)-1,
        if collision(newPath(end,:), path(ptr+1,:), map)~=0, % if there is a collision
            newPath = [newPath; path(ptr,:)];  % add previous node
        end
        ptr=ptr+1;
    end
    newPath = [newPath; path(end,:)];

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% plotmap
%%   plot obstacles and path
function plotmap(map,path,smoothedPath,tree)

    % setup plot
    figure(3), clf
    axis([0,map.width,0,map.width,0,2*map.MaxHeight]);
    xlabel('E')
    ylabel('N')
    zlabel('h')
    hold on

    % plot buildings
    V = [];
    F = [];
    patchcolors = [];
    count = 0;
    for i=1:map.NumBlocks,
        for j=1:map.NumBlocks,
            [Vtemp,Ftemp,patchcolorstemp] = buildingVertFace(map.buildings_n(i),...
                map.buildings_e(j),map.BuildingWidth,map.heights(j,i));
            V = [V; Vtemp];
            Ftemp = Ftemp + count;
            F = [F; Ftemp];
            count = count + 8;
            patchcolors = [patchcolors;patchcolorstemp];
        end
    end

    patch('Vertices', V, 'Faces', F,...
                     'FaceVertexCData',patchcolors,...
                     'FaceColor','flat');

    % draw tree
    for i=2:size(tree,1),
        X = [tree(i,1), tree(tree(i,6),1)];
        Y = [tree(i,2), tree(tree(i,6),2)];
        Z = [tree(i,3), tree(tree(i,6),3)];
        plot3(Y,X,-Z,'g')
    end

    % draw path
    X = path(:,1);
    Y = path(:,2);
    Z = path(:,3);
    plot3(Y,X,-Z,'r','linewidth',2);

    % draw smooth path
    X = smoothedPath(:,1);
    Y = smoothedPath(:,2);
    Z = smoothedPath(:,3);
    plot3(Y,X,-Z,'k','linewidth',2);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% buildingVertFace(x,y,width,height)
%%   define patches for a building located at (x,y)
function [V,F,patchcolors] = buildingVertFace(n,e,width,height)

  % vertices of the building
  V = [...
        e+width/2, n+width/2, 0;...
        e+width/2, n-width/2, 0;...
        e-width/2, n-width/2, 0;...
        e-width/2, n+width/2, 0;...
        e+width/2, n+width/2, height;...
        e+width/2, n-width/2, height;...
        e-width/2, n-width/2, height;...
        e-width/2, n+width/2, height;...
        ];
  % define faces of fuselage
  F = [...
        1, 4, 8, 5;... % North Side
        1, 2, 6, 5;... % East Side
        2, 3, 7, 6;... % South Side
        3, 4, 8, 7;... % West Side
        5, 6, 7, 8;... % Top
        ];

  myred = [1, 0, 0];
  mygreen = [0, 1, 0];
  myblue = [0, 0, 1];
  myyellow = [1,1,0];
  mymagenta   = [0, 1, 1];

  patchcolors = [...
    mygreen;... % North
    mygreen;... % East
    mygreen;... % South
    mygreen;... % West
    myyellow;...  % Top
    ];

end


text:
% GEODOME   Generates geodesic sphere
%
% Usage:  [xyz, A, F] = geodome(frequency, radius)
%
% The sphere is generated from subdivisions of the faces of an icosahedron.
%
% Arguments:
%        frequency - The number of subdivisions of each edge of an
%                    icosahedron.  A frequency of 1 will give you an
%                    icosahedron. A high number will give you a more
%                    spherical shape. Defaults to 2.
%           radius - Radius of the sphere. Defaults to 1.
% Returns:
%              xyz - A nx3 matrix defining the vertex list, each row is the
%                    x,y,z coordinates of a vertex.
%                A - Adjacency matrix defining the connectivity of the
%                    vertices
%                F - A mx3 matrix defining the face list.  Each row of F
%                    specifies the 3 vertices that make up the face.
%
% Example:
%    [xyz, A, F] = geodome(3, 5);  % Generate a 3-frequency sphere with
%                                  % radius 5
%    gplot3d(A, xyz);              % Use adjacency matrix and vertex list to
%                                  % generate a wireframe plot.
%    drawfaces(xyz, F);            % Use face and vertex lists to generate
%                                  % surface patch plot.
%    axis vis3d, axis off, rotate3d on
%
% See also: ICOSAHEDRON, GPLOT3D, DRAWFACES

% Copyright (c) 2009 Peter Kovesi
% Centre for Exploration Targeting
% The University of Western Australia
% peter.kovesi at uwa edu au
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.

% May   2009
% April 2014  - Dome radius scaling fixed (thanks to Brad Keserich)

function [xyz, A, F] = geodome(frequency, radius)

    if ~exist('frequency','var')
        frequency = 2;
    end

    if ~exist('radius','var')
        radius = 1;
    end

    if frequency < 1
        error('Frequency must be an integer >= 1');
    end

    % Generate vertices of base icosahedron
    [icosXYZ, icosA, icosF] = icosahedron(1);

    % Compute number of vertices in geodesic sphere.  This is the 12 vertices
    % of the icosahedron + the extra vertices introduced on each of the 30
    % edges + the extra vertices in the interior of the of each of the 20
    % faces.  This latter is a sum of an arithmetic series [1 .. (frequency-2)]
    fm1 = frequency - 1;
    fm2 = frequency - 2;
    nVert = 12 + 30*fm1 + 20*fm2*fm1/2;

    xyz = zeros(nVert,3);

    if frequency == 1             % Just return the icosahedron
        xyz = icosXYZ;
        A = icosA;
        F = icosF;
    else                          % For all frequencies > 1
        xyz(1:12,:) = icosXYZ;    % Grab the vertices of the icosahedron
        offset = 13;

        % Find the nodes that connect edges.  Note we use the upper triangular
        % part of the adjacency matrix so that we only get each edge once.
        [i,j] = find(triu(icosA));

        % Generate extra vertices on every edge of the icoshedron.
        for n = 1:length(i)
            xyz(offset:(offset+fm2) ,:) = ...
                   divideEdge(icosXYZ(i(n),:), icosXYZ(j(n),:), frequency);
            offset = offset+fm1;
        end

        % Generate the extra vertices within each face of the icoshedron.
        for f = 1:length(icosF)

            % Re subdivide two of the edges of the face and get the vertices
            % (Yes, this is wasteful but it makes code logic easier)
            V1 = divideEdge(icosXYZ(icosF(f,1),:), icosXYZ(icosF(f,2),:), frequency);
            V2 = divideEdge(icosXYZ(icosF(f,1),:), icosXYZ(icosF(f,3),:), frequency);

            % Now divide the edges that join the new vertices along the
            % subdivided edges.
            for v = 2:fm1
               VF = divideEdge(V1(v,:), V2(v,:), v);
               xyz(offset:(offset+v-2),:) = VF;
               offset = offset+v-1;
            end
        end

    end

    xyz = xyz*radius;       % Scale vertices to required radius
    A = adjacency(xyz);
    F = faces(A);

%---------------------------------------------------------------------
% Function to divide an edge defined between vertices V1 and V2 into nSeg
% segments and return the coordinates of the vertices.
% This function simplistically divides the distance between V1 and V2

function vert = divideEdgeOld(V1, V2, nSeg)

    edge = V2 - V1;  % Vector along edge

    % Now add appropriate fractions of the edge length to the first node
    vert = zeros(nSeg-1, 3);
    for f = 1:(nSeg-1)
        vert(f,:) = V1 + edge * f/nSeg;
        vert(f,:) = vert(f,:)/norm(vert(f,:));   % Normalize to unit length
    end

%---------------------------------------------------------------------
% Function to divide an edge defined between vertices V1 and V2 into nSeg
% segments and return the coordinates of the vertices.
% This function divides the *angle* between V1 and V2
% rather than the distance.
function vert = divideEdge(V1, V2, nSeg)

    axis = cross(V1,V2);
    angle = atan(norm(axis)/dot(V1,V2));
    axis = axis/norm(axis);

    % Now add appropriate fractions of the edge length to the first node
    vert = zeros(nSeg-1, 3);
    for f = 1:(nSeg-1)
        Q = newquaternion(f*angle/nSeg, axis);
        vert(f,:) =  quaternionrotate(Q,V1);
        vert(f,:) = vert(f,:)/norm(vert(f,:));   % Normalize to unit length
    end

%-------------------------------------------------------------------------
% Function to build adjacency matrix for geodesic sphere by brute force

function A = adjacency(xyz)

    nVert = length(xyz);
    A = zeros(nVert);

    % Find distances between all pairs of vertices
    for n1 = 1:nVert-1
        A(n1,n1) = Inf;
        for n2 = n1+1:nVert
            A(n1,n2) = norm(xyz(n2,:) - xyz(n1,:));
        end
    end
    A(nVert,nVert) = Inf;

    A = A+A';   % Make A symmetric

    % Find min distance in first row.
    minD = min(A(1,:)');

    % Assume that no edge can be more than 1.5 times this minimum distance, use
    % this to decide connectivity of nodes.
    A = A < minD*1.5;

%-----------------------------------------------------------------------------
% Function to find the triplets of vertices that define the faces of the
% geodesic sphere

function  F = faces(A)

    % Strategy:  We are only after cycles of length 3 in graph defined by A.
    % For every node N0 (except the last two, which we will not need to visit)
    % - Get list of neighbours, N1
    % - For each neighbour N1i in N1 find its list of neighbours, N2.
    % - Any neighbour in N2 that is in N1 must form a cycle of length 3 with N0

    nVert = length(A);
    F = [];

    for N0 = 1:nVert-2
        N1 = find(A(N0,:));            % Neighbours of N0
        for N1i = N1
            N2 = find(A(N1i,:));
            cycle = intersect(N2,N1);  % Find the 2 nodes of N2 that are in N1
            F = [F; [N0 N1i cycle(1)]; [N0 N1i cycle(2)]];
        end
    end

    % Each face will be found multiple times, eliminate duplicates.
    F = sort(F,2);            % Sort rows of face list
    F = unique(F,'rows');     % ...then extract the unique rows.

author: alecjacobson

repo: gptoolbox


text:
function a = gt( x, y )

%Disciplined convex programming information for GT (>):
%   The right-hand side of a less-than constraint must be convex. The
%   left-hand side must be concave. Of course, real constant and affine
%   expressions are both convex and concave and can be used on either
%   side as well.
%
%Disciplined geometric programming information for GT (>):
%   The right-hand side of a less-than constraint must be log-convex---
%   including positive constants, monomials, posynomials, generalized
%   posynomials, and products thereof. The left-hand side must be log-
%   concave---including positive constants, monomials, reciprocals of
%   log-convex expressions, and products thereof.
%
%Note that CVX does not distinguish between strict greater-than (>) and
%greater-than-or-equal (>=) constraints; they are treated identically.
%Feasible interior-point solvers tend to return points which satisfy
%strict inequality, but not all solvers do.

cvx_warn_strict( '>' );
evalin( 'caller', 'cvx_verify' );
b = cvx_pushcnstr( x, y, '>' );
if nargout, a = b; end

% 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 out=acot(x)

out=atan(1./x);

author: opencobra

repo: cobratoolbox


text:
classdef CF7 < PROBLEM
% <multi> <real> <large/none> <constrained>
% Constrained benchmark MOP

%------------------------------- Reference --------------------------------
% Q. Zhang, A. Zhou, S. Zhao, P. N. Suganthan, W. Liu, and S. Tiwari,
% Multiobjective optimization test instances for the CEC 2009 special
% session and competition, School of CS & EE, University of Essex, Working
% Report CES-487, 2009.
%------------------------------- 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".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 10; end
            obj.lower    = [0,zeros(1,obj.D-1)-2];
            obj.upper    = [1,zeros(1,obj.D-1)+2];
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            D  = size(X,2);
            J1 = 3 : 2 : D;
            J2 = 2 : 2 : D;
            Y           = zeros(size(X));
            Y(:,J1)     = X(:,J1) - cos(6*pi*repmat(X(:,1),1,length(J1))+repmat(J1,size(X,1),1)*pi/D);
            Y(:,J2)     = X(:,J2) - sin(6*pi*repmat(X(:,1),1,length(J2))+repmat(J2,size(X,1),1)*pi/D);
            h           = 2*Y.^2 - cos(4*pi*Y) + 1;
            h(:,[2,4])  = Y(:,[2,4]).^2;
            PopObj(:,1) = X(:,1)        + sum(h(:,J1),2);
            PopObj(:,2) = (1-X(:,1)).^2 + sum(h(:,J2),2);
        end
        %% Calculate constraint violations
        function PopCon = CalCon(obj,X)
            PopCon(:,1) = -X(:,2) + sin(6*pi*X(:,1)+2*pi/size(X,2)) + sign(0.5*(1-X(:,1))-(1-X(:,1)).^2).*sqrt(abs(0.5*(1-X(:,1))-(1-X(:,1)).^2));
            PopCon(:,2) = -X(:,4) + sin(6*pi*X(:,1)+4*pi/size(X,2)) + sign(0.25*sqrt(1-X(:,1))-0.5*(1-X(:,1))).*sqrt(abs(0.25*sqrt(1-X(:,1))-0.5*(1-X(:,1))));
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = (1-R(:,1)).^2;
            temp1  = 0.5<R(:,1) & R(:,1)<=0.75;
            temp2  = 0.75<R(:,1);
            R(temp1,2) = 0.5*(1-R(temp1,1));
            R(temp2,2) = 0.25*sqrt(1-R(temp2,1));
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R = obj.GetOptimum(100);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [Population,FrontNo,CrowdDis] = EnvironmentalSelection(Population,N,Points,W,delta)
% The environmental selection of r-NSGA-II

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

    %% Non-r-dominated sorting
    [FrontNo,MaxFNo] = NrDSort(Population.objs,N,Points,W,delta);
    Next = FrontNo < MaxFNo;

    %% Calculate the crowding distance of each solution
    CrowdDis = CrowdingDistance(Population.objs,FrontNo);

    %% Select the solutions in the last front by their crowding distances
    Last     = find(FrontNo==MaxFNo);
    [~,Rank] = sort(CrowdDis(Last),'descend');
    Next(Last(Rank(1:N-sum(Next)))) = true;

    %% Population for next generation
    Population = Population(Next);
    FrontNo    = FrontNo(Next);
    CrowdDis   = CrowdDis(Next);
end


text:
function Jb = JacobianBody_sym(Blist, thetalist)
% *** CHAPTER 5: VELOCITY KINEMATICS AND STATICS ***
% Takes Blist: The joint screw axes in the end-effector frame when the
%              manipulator is at the home position, in the format of a
%              matrix with the screw axes as the columns,
%       thetalist: A list of joint coordinates.
% Returns the corresponding body Jacobian (6xn real numbers).
% Example Input:
%
% clear; clc;
% Blist = [[0; 0; 1;   0; 0.2; 0.2], ...
%        [1; 0; 0;   2;   0;   3], ...
%        [0; 1; 0;   0;   2;   1], ...
%        [1; 0; 0; 0.2; 0.3; 0.4]];
% thetalist = [0.2; 1.1; 0.1; 1.2];
% Jb = JacobianBody(Blist, thetalist)
%
% Output:
% Jb =
%   -0.0453    0.9950         0    1.0000
%    0.7436    0.0930    0.3624         0
%   -0.6671    0.0362   -0.9320         0
%    2.3259    1.6681    0.5641    0.2000
%   -1.4432    2.9456    1.4331    0.3000
%   -2.0664    1.8288   -1.5887    0.4000

Jb = sym(Blist);
T = eye(4);
for i = length(thetalist) - 1: -1: 1
    T = T * simplify(MatrixExp6_sym(VecTose3(Blist(:, i + 1)),-1*thetalist(i + 1)));
	Jb(:, i) = Adjoint(T) * Blist(:, i);
end
end

author: ShuoYangRobotics

repo: QuadrupedSim


text:
function f = spm_Dpdf(x,a)
% Probability Density Function (PDF) of Dirichlet distribution
% FORMAT f = spm_Dpdf(x,a)
%
% x - Dirichlet variate
% a - Dirichlet parameters (a>0)
% f - PDF of Dirichlet-distribution at point x
%__________________________________________________________________________
%
% spm_Dpdf implements the Probability Density Function for Dirichlet
% distribution.
%
% Definition:
%--------------------------------------------------------------------------
% See http://en.wikipedia.org/wiki/Dirichlet_distribution
%
% Algorithm:
%--------------------------------------------------------------------------
% Direct computation using logs and MATLAB's implementation of the log of
% the gamma function (gammaln).
%__________________________________________________________________________
% Copyright (C) 2008-2011 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_Dpdf.m 4182 2011-02-01 12:29:09Z guillaume $

%-Check enough arguments
%--------------------------------------------------------------------------
if nargin<2, error('Insufficient arguments'), end

%-Computation
%--------------------------------------------------------------------------
a = a(:);
x = x(:);

f = exp( gammaln(sum(a)) + sum((a-1).*log(x+eps)) - sum(gammaln(a)) );

author: spm


text:
function m=v_pesq2mos(p)
%V_PESQ2MOS convert PESQ speech quality scores to MOS m=(p)
%Inputs:    p  is a matrix of PESQ scores
%
%Outputs:   m  is a matrix, the same size as p, of MOS scores
%
% The PESQ measure is defined in [2]. The mapping function, defined in [3],
% converts raw PESQ scores (which lie in the range -0.5 to 4.5) onto the
% MOS-LQO (Mean Opinion Score - Listening Quality Objective [2]) scale in the
% range 1 to 5. The MOS scale is defined in [1] as
%           5=Excellent, 4=Good, 3=Fair, 2=Poor, 1=Bad.
%
% Refs: [1]	ITU-T. Methods for subjective determination of transmission quality.
%           Recommendation P.800, Aug. 1996.
%       [2]	ITU-T. Mean opinion score (MOS) terminology.
%           Recommendation P.800.1, July 2006.
%       [2]	ITU-T. Perceptual evaluation of speech quality (PESQ), an objective
%           method for end-to-end speech quality assessment of narrowband telephone
%           networks and speech codecs. Recommendation P.862, Feb. 2001.
%       [3]	ITU-T. Mapping function for transforming P.862 raw result scores to MOS-LQO.
%           Recommendation P.862.1, Nov. 2003.

%      Copyright (C) Mike Brookes 2012-2013
%      Version: $Id: v_pesq2mos.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
persistent a b c d
if isempty(a)
    a=0.999;
    b=4.999-a;
    c=-1.4945;
    d=4.6607;
end
if nargout>0
    m=a+b./(1+exp(c*p+d));
else
    if nargin<1 || isempty(p)
        pp=linspace(-0.5,4.5,100);
    else
        pp=p;
    end
    plot(pp,v_pesq2mos(pp));
    xlabel('PESQ (P.862)');
    ylabel('Mean Opimion Score (MOS)');
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function outputinchar = converttochar(inputinsym)

% This file is part of SOSTOOLS - Sum of Squares Toolbox ver 3.00.
%
% Copyright (C)2002, 2004, 2013  A. Papachristodoulou (1), J. Anderson (1),
%                                G. Valmorbida (1), S. Prajna (2),
%                                P. Seiler (3), P. A. Parrilo (4)
% (1) Department of Engineering Science, University of Oxford, Oxford, U.K.
% (2) Control and Dynamical Systems - California Institute of Technology,
%     Pasadena, CA 91125, USA.
% (3) Aerospace and Engineering Mechanics Department, University of
%     Minnesota, Minneapolis, MN 55455-0153, USA.
% (4) Laboratory for Information and Decision Systems, M.I.T.,
%     Massachusetts, MA 02139-4307
%
% Send bug reports and feedback to: sostools@cds.caltech.edu
%
% 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, see <http://www.gnu.org/licenses/>.

char1 = char(inputinsym);
verstring = char(version);
if verstring(1) =='7' || verstring(1) == '8'
    char2 = strrep(char1,'matrix([','');
else
    char2 = strrep(char1,'array([','');
end;
outputinchar = strrep(char2,'])','');
if outputinchar(1)~='['
    outputinchar = ['[',outputinchar,']'];
end;


text:
% DEMROBOTWIRELESS3 Wireless Robot data from University of Washington with dynamics and no back constraints.

% FGPLVM

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

dataSetName = 'robotWireless';
experimentNo = 3;

% load data
[Y, lbls] = lvmLoadData(dataSetName);

% Set up model
options = fgplvmOptions('ftc');

latentDim = 2;
d = size(Y, 2);

model = fgplvmCreate(latentDim, d, Y, options);

% Add dynamics model.
options = gpOptions('ftc');
options.kern = kernCreate(model.X, {'rbf', 'white'});
options.kern.comp{1}.inverseWidth = 0.2;
% This gives signal to noise of 0.1:1e-3 or 100:1.
options.kern.comp{1}.variance = 0.1^2;
options.kern.comp{2}.variance = 1e-3^2;
model = fgplvmAddDynamics(model, 'gp', options);

% Optimise the model.
iters = 1000;
display = 1;

model = fgplvmOptimise(model, display, iters);

% Save the results.
capName = dataSetName;;
capName(1) = upper(capName(1));
save(['dem' capName num2str(experimentNo) '.mat'], 'model');

% Load the results and display dynamically.
fgplvmResultsDynamic(dataSetName, experimentNo, 'vector')

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function strDate = str_date(s, dateFormat)
% STR_Date: Reformat date string to dd-MMM-yyyy.

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2018

% Parse inputs
if (nargin < 2) || isempty(dateFormat)
    dateFormat = [];
end
% Clean string
if ischar(s)
    s = strtrim(strrep(s, char(0), ''));
end
% Check various input formats
try
    if ~isempty(dateFormat) && strcmpi(dateFormat, 'posix')
        strDate = datestr(double(s) ./ 86400 + datenum(1970,1,1,0,0,0), 'dd-mmm-yyyy');
    elseif ~isempty(dateFormat)
        strDate = datestr(datenum(s, dateFormat), 'dd-mmm-yyyy');
    elseif isequal(find(s == '/'), [3 6]) && ((length(s) == 10) || (length(s) == 19))
        strDate = datestr(datenum(s, 'dd/mm/yyyy'), 'dd-mmm-yyyy');
    elseif isequal(find(s == '/'), [3 6]) && ((length(s) == 8) || (length(s) == 17))
        strDate = datestr(datenum(s, 'dd/mm/yy'), 'dd-mmm-yyyy');
    elseif isequal(find(s == '.'), [3 6]) && ((length(s) == 10) || (length(s) == 19))
        strDate = datestr(datenum(s, 'dd.mm.yyyy'), 'dd-mmm-yyyy');
    elseif isequal(find(s == '.'), [3 6]) && ((length(s) == 8) || (length(s) == 17))
        strDate = datestr(datenum(s, 'dd.mm.yy'), 'dd-mmm-yyyy');
    elseif isequal(find(s == '-'), [3 6]) && ((length(s) == 10) || (length(s) == 19))
        strDate = datestr(datenum(s, 'dd-mm-yyyy'), 'dd-mmm-yyyy');
    elseif isequal(find(s == '-'), [3 7]) && (length(s) == 11)
        strDate = datestr(datenum(s, 'dd-mmm-yyyy'), 'dd-mmm-yyyy');
    elseif isequal(find(s == '-'), [5 8]) && ((length(s) == 10) || (length(s) == 19))
        strDate = datestr(datenum(s, 'yyyy-mm-dd'), 'dd-mmm-yyyy');
    else
        strDate = [];
    end
catch
    strDate = [];
end

author: brainstorm-tools