text:
%{
Solves a simple quadratic problem
(with optimal solution known in closed-form)
with all the various solvers

the quadratic problem is:

min_x  c'x + x'Px/2, P > 0, hence the optimal solution is x = -inv(Q)*c

On a simple quadratic example like this, when there is no projection,
all the solver behave similarly. So this script is more useful for debugging
to make sure that all solvers are working

%}

%% Make a problem
N   = 100;
randn('state',2334);
P   = randn(N); P = P*P' + .1*eye(N); % ensure P is symmetric positive definite
c   = randn(N,1);

x_ref   = -P\c;

%% Solve with TFOCS
f   = smooth_quad( P, c );
affine    = [];
projector = [];
x0        = zeros(N,1);
opts      = [];
opts.printEvery     = 500;
opts.tol  = 1e-8;
opts.printStopCrit  = true;
opts.restart    = 350;  % use this since the problem is strongly convex
opts.errFcn     = @(f,x) norm( x-x_ref)/norm(x_ref);
opts.maxIts     = 750;

% The 'N83' and 'AT' methods can take advantage when you know
%   the exact value of the strong convexity parameter
%   (and the Lipschitz constant).  In this case, do not
%   use restart. But we won't test this right now,
%   since it is unusual that you have information
%   about the strong convexity constant.
% opts.restart    = Inf;
% opts.mu         = min(eig(P)); % strong convexity parameter. N83 can make use of this
% opts.Lexact     = max(eig(P));

solverList = { 'GRA', 'AT', 'LLM', 'N07', 'N83', 'TS' };
% solverList = { 'GRA', 'AT', 'N83', 'TS' };
% solverList = { 'LLM','N07'}; % these take two steps per iteration

figure(1); clf;

for k = 1:length( solverList )
    solver  = solverList{k};
    opts.alg    = solver;
    fprintf('\nSolver: %s\n\n', solver );
    [x,out,optsOut] = tfocs( f, affine, projector, x0, opts );
    semilogy( out.err );
    hold all
end
legend( solverList );

%% close all figures
close all

% 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:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function [Optionspreis] = BlackScholesPrice(S,K,r,T,sigma,type)
% calculates the black scholes price using matrices
    matS = repmat(S,1,length(T));   % S matrix
    matT = repmat(T',length(S),1);  % T matrix

    d1=(log(matS./K)+(r+sigma^2/2)*matT)./(sigma*sqrt(matT));
    d2=(log(matS./K)+(r-sigma^2/2)*matT)./(sigma*sqrt(matT));

    if type==0
        Optionspreis = K*exp(-r*matT).*normcdf(-d2)-matS.*normcdf(-d1);
    else
        Optionspreis = matS.*normcdf(d1)-K*exp(-r*matT).*normcdf(d2);
    end

end

author: Sable

repo: mcbench-benchmarks


text:
function EEG = highPassAndICA(EEG, varargin)
% Perform a high-pass filter and ICA on data that has been Prepped
%
% Parameters
%    EEG   EEGLAB EEG structure
%    varargin    name value pairs
%         detrendChannels
%         detrendCutoff
%         icatype
%         extended
try
    %% Convert name-value pair parameters to structure
    params = vargin2struct(varargin);

    if isfield(EEG.etc, 'noiseDetection') && ...
            isfield(EEG.etc.noiseDetection, 'detrend') && ...
            ~isfield(params, 'detrendChannels')
        params.detrendChannels = EEG.etc.noiseDetection.detrend.detrendChannels;
    end
    [EEG, detrend]  = removeTrend(EEG, params);
    EEG.etc.highPassAndICA.detrend = detrend;

    %% Perform ICA on the data
    icatype = 'runica';
    if isfield(params, 'icatype')
        icatype = params.icatype;
    end
    extended = 1;
    if isfield(params, 'extended')
        extended = params.extended;
    end
    if isfield(EEG.etc, 'noiseDetection') && ...
            isfield(EEG.etc.noiseDetection, 'reference')
        referenceChannels = ...
            EEG.etc.noiseDetection.reference.referenceChannels;
        interpolatedChannels = ...
            EEG.etc.noiseDetection.reference.interpolatedChannels.all;
        channelsLeft = setdiff(referenceChannels, interpolatedChannels);
        pcaDim = length(channelsLeft) - 1;
        fprintf('%s: pcaDim = %d\n', EEG.setname, pcaDim);
        EEG = pop_runica(EEG, 'icatype', icatype, 'extended', extended, ...
            'chanind', referenceChannels, 'pca', pcaDim);
        EEG.etc.highPassAndICA.ICA = ...
            ['extended infomax with pca Dim ' num2str(pcaDim)];
        EEG.etc.noiseDetection.reference = ...
                   cleanupReference(EEG.etc.noiseDetection.reference);
    else
        EEG = pop_runica(EEG, 'icatype', icatype, 'extended', extended);
        EEG.etc.highPassAndICA.ICA = [icatype '_extended_' extended];
    end

catch mex
    errorMessages.highPassAndICA = ['failed highPassAndICA: ' getReport(mex)];
    errorMessages.status = 'unprocessed';
    EEG.etc.highPassAndICA.errors = errorMessages;
    fprintf(2, '%s\n', errorMessages.highPassAndICA);
end

author: VisLab

repo: EEG-Clean-Tools


text:
%H2E Homogeneous to Euclidean
%
% E = H2E(H) is the Euclidean version (K-1xN) of the homogeneous
% points H (KxN) where each column represents one point in P^K.
%
% Reference::
% - Robotics, Vision & Control: Second Edition, P. Corke, Springer 2016; p604.
%
% See also E2H.

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

function e = h2e(h)

    if isvector(h)
        h = h(:);
    end
    e = h(1:end-1,:) ./ repmat(h(end,:), numrows(h)-1, 1);

author: petercorke

repo: spatialmath-matlab


text:
function x = ortega_eigen_right ( n, u, v, d )

%*****************************************************************************80
%
%% ORTEGA_EIGEN_RIGHT returns the right eigenvectors of the ORTEGA matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    James Ortega,
%    Generation of Test Matrices by Similarity Transformations,
%    Communications of the ACM,
%    Volume 7, 1964, pages 377-378.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    2 <= N.
%
%    Input, real U(N), V(N), vectors which define the matrix.
%    U'V must not equal -1.0.  If, in fact, U'V = 0, and U, V and D are
%    integers, then the matrix, inverse, eigenvalues, and eigenvectors
%    will be integers.
%
%    Input, real D(N), the desired eigenvalues.
%
%    Output, real X(N,N), the determinant.
%
  x = zeros(n,n);

  for j = 1 : n
    for i = 1 : n

      if ( i == j )
        x(i,j) = 1.0 + u(i) * v(j);
      else
        x(i,j) =       u(i) * v(j);
      end

    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%
% Usage:   Y =mexBayer(X,offset);
%
% Name: mexBayer
%
% Description: mexBayer applies a Bayer pattern to an image X.
%       There are four possible offsets.
%
% Inputs: X:  double m x n matrix
%         offset: scalar, 0,1,2 or 3
%
% Output: Y: double m x m matrix
%
% Author: Julien Mairal, 2009

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function setPlotModeROE(hObj,hEvt,hFig)
% Set ROE  plot mode
%
% AI 05/12/21

ud = guidata(hFig);
sel = get(hObj,'Value');
ud.plotMode = sel - 1;
if ud.plotMode==3
    set(ud.handle.modelsAxis(10),'String','Enter scale factor');
    set(ud.handle.modelsAxis(9),'Visible','On','Enable','On');
    set(ud.handle.modelsAxis(10),'Visible','On');
elseif ud.plotMode==4
    txt = sprintf('Enter\n \x0394nfrx');
    set(ud.handle.modelsAxis(10),'String',txt);
    set(ud.handle.modelsAxis(9),'Visible','On','Enable','On');
    set(ud.handle.modelsAxis(10),'Visible','On');
else
    set(ud.handle.modelsAxis(9),'Visible','Off');
    set(ud.handle.modelsAxis(10),'Visible','Off');
end

guidata(hFig,ud);

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function []=plot_dramapb(vResults)
% --------------------------------------------
% Plot db map on topography
%
% Incoming variables:
% vResults: Necessarily includes vResults.mPolygon, vResults.mValueGrid, vResults.fSpacingHorizontal
% tmap    : Topographic map data ETOPO5
% tmapleg : Topographic map data ETOPO5 legend
%
% Author: J. Woessner, woessner@seismo.ifg.ethz.ch
% last update: 27.02.03

% Get tmap from tbase.bin
[tmap,tmapleg] = tbase(3);

% Construct a map graticule mesh for surface object display
[lat,lon] = meshgrat(tmap,tmapleg);

% Get coordinates from calculated Mc values, set up matrices
vX=min(vResults.mPolygon(:,1)):vResults.fSpacingHorizontal:max(vResults.mPolygon(:,1));
vY=min(vResults.mPolygon(:,2)):vResults.fSpacingHorizontal:max(vResults.mPolygon(:,2));
[X , Y]  = meshgrid(vX,vY);

% % Set values of dMc to NaN
vSel = find(vResults.mValueGrid(:,10)> 0.5);
vResults.mValueGrid(vSel,10) = NaN;

% Reshape
mMcValues=reshape(vResults.mValueGrid(:,10),length(vY),length(vX));
% Interpolate Values
mMcGrid = interp2(X,Y,mMcValues,lon,lat);

% Create figure
exfig=figure_exists('db Worldmap',0);
if exfig==0
    ui=figure_w_normalized_uicontrolunits( 'Name','db Worldmap',...
        'NumberTitle','off');
    figure_w_normalized_uicontrolunits(ui)
end

clf
hold on;
axis off;
axesm('MapProjection','robinson')

% Magic coloring
mi = min(min(mMcGrid));
l =  isnan(mMcGrid);
mMcGrid(l) = mi-10;
ll = tmap < 0 & mMcGrid < 0;
mMcGrid(ll) = mMcGrid(ll)*0 + 10;
hMap=meshm(mMcGrid,tmapleg,size(tmap),tmap);

daspectm('m',20);
tightmap
view([0 90])
camlight; lighting phong
set(gca,'projection','perspective');

%j = jet;
j = hot(400000);
j = j(500:350000,:);
j = flipud(j);
j = [ [ 0.9 0.9 0.9 ] ; j; [ 0.5 0.5 0.5] ];

%caxis([ min(min(mMcValues)) max(max(mMcValues))]);
caxis([0 0.5])
colormap(j); brighten(0.1);

axis off; set(gcf,'color','k')

setm(gca,'ffacecolor','k')
setm(gca,'fedgecolor','w','flinewidth',1);

setm(gca,'mlabellocation',60)
setm(gca,'meridianlabel','on')
setm(gca,'plabellocation',25)
setm(gca,'parallellabel','on')
setm(gca,'Fontcolor','w','FontSize',10,'Fontweight','bold','Fontname','times','Labelunits','degrees')
setm(gca,'Labelformat','none')
setm(gca,'Grid','on')
setm(gca,'GLineStyle','-')

h5 = colorbar;
set(h5,'position',[0.82 0.35 0.01 0.3],'TickDir','out','Ycolor','w','Xcolor','w',...
    'Fontweight','bold','FontSize',10,'YTick',[0 0.1 0.2 0.3 0.4 0.5 0.6]);
set(gcf,'Inverthardcopy','off');

author: CelsoReyes


text:
% function [img]=jp2read(filename)
%
% reads a JPEG2000 compliant file and stores information
% in a matlab image

function img=jp2read(filename)

    import jj2000.disp.*;
    import jj2000.j2k.util.*;
    import jj2000.j2k.decoder.*;
    import java.awt.image.*;
    import java.awt.*;
    import java.io.*;

    defpl=ParameterList();
    params=Decoder.getAllParameters();

    for i=1:size(params,1)
        if ne(params(i,4),'')
            defpl.setProperty(params(i,1),params(i,4));
        end
    end

    defpl.setProperty('i',filename);
    defpl.setProperty('v','on');
    defpl.setProperty('debug','on');

    isp=ImgScrollPane();

    d=Decoder(ParameterList(defpl),isp);

    ml=StreamMsgLogger(FileOutputStream('out.log'), ...
        FileOutputStream('err.log'),78);

    FacilityManager.registerMsgLogger('',ml);

    d.run();

    image=isp.getImage();

    w=image.getWidth();
    h=image.getHeight();

    pg=PixelGrabber(image,0,0,w,h,false);
    pg.startGrabbing();
    a=pg.getPixels();

    img=java2im(a,h,w);
return


text:
function prob_test136 ( )

%*****************************************************************************80
%
%% TEST136 tests SECH_MEAN, SECH_SAMPLE, SECH_VARIANCE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  nsample = 1000;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST136\n' );
  fprintf ( 1, '  For the Hyperbolic Secant PDF:\n' );
  fprintf ( 1, '  SECH_MEAN computes the mean;\n' );
  fprintf ( 1, '  SECH_SAMPLE samples;\n' );
  fprintf ( 1, '  SECH_VARIANCE computes the variance.\n' );

  a = 3.0;
  b = 2.0;

  check = sech_check ( a, b );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST136 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  mean = sech_mean ( a, b );
  variance = sech_variance ( a, b );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A =             %14f\n', a );
  fprintf ( 1, '  PDF parameter B =             %14f\n', b );
  fprintf ( 1, '  PDF mean =                    %14f\n', mean );
  fprintf ( 1, '  PDF variance =                %14f\n', variance );

  for i = 1 : nsample
    [ x(i), seed ] = sech_sample ( a, b, seed );
  end

  mean = r8vec_mean ( nsample, x );
  variance = r8vec_variance ( nsample, x );
  xmax = max ( x(1:nsample) );
  xmin = min ( x(1:nsample) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sample size =     %6d\n', nsample );
  fprintf ( 1, '  Sample mean =     %14f\n', mean );
  fprintf ( 1, '  Sample variance = %14f\n', variance );
  fprintf ( 1, '  Sample maximum =  %14f\n', xmax );
  fprintf ( 1, '  Sample minimum =  %14f\n', xmin );

  return
end


text:
function v = norm(F)
%NORM   Frobenius norm of a CHEBFUN3V object.
%   V = NORM(F) returns the Frobenius norm of F, i.e.
%       V = sqrt(norm(F1).^2 + norm(F2).^2),
%   or
%       V = sqrt(norm(F1).^2 + norm(F2).^2 + norm(F3).^2) .
%   where F = [F1 F2 F3]^T.

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

% Empty check:
if ( isempty(F) )
    v = [];
    return
end

nF = F.nComponents;
v = 0;
for jj = 1:nF
    v = v + norm(F.components{jj})^2;
end
v = sqrt(v);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function W=calcKalmanGain(R,PzPred,otherInfo)
%%CALCKALMANGAIN Using the output of the Kalman filter measurement
%           prediction function KalmanMeasPred and the covariance
%           matrix of a measurement R, compute the gain of a Kalman filter.
%           This function can be useful for creating a gain based on some
%           type of "maximum" covariance matrix for a scene for purposes of
%           giving it to the calcMissedGateCov function to determine how
%           the covariance matrix of the target state prediction should be
%           increased for a missed-detection event when gating is
%           performed.
%
%INPUTS: R The zDimXzDim measurement covariance matrix.
%   PzPred The zDimXzDim covariance matrix of the measurement predicted by
%          the filter.
% otherInfo The structure returned by the KalmanMeasPred function that
%          includes various terms that can be reused.
%
%OUTPUTS: W The xDimXzDim gain matrix used in the filter.
%
%This function produces the gain in the same manner as is done in
%KalmanUpdateWithPred and KalmanUpdate. See the comments to those functions
%for more information.
%
%EXAMPLE:
%Here, we demonstrate that the gain used in a complete update is the same
%as that obtained using just the prediction via KalmanMeasPred and then
%calcKalmanGain for the gain without using the measurement z itself.
% xPred=[1e3;-2e3;100;200];
% PPred=[28,   3.5,    6,  8.5;
%      3.5,    23,  8.5,   11;
%        6,   8.5,   18, 13.5;
%      8.5,    11, 13.5,   13];
% z=1e3*[-5.498856156296510;
%        1.199241491470584];
% R=eye(2);
% H=[0, 4, 9, 8;
%    6, 3, 0, 6];
% [~,~,~,~,W]=KalmanUpdate(xPred,PPred,z,R,H);
% [~,PzPred,otherInfo]=KalmanMeasPred(xPred,PPred,H);
% W1=calcKalmanGain(R,PzPred,otherInfo);
% %One will see that the result below is true (1).
% all(W1(:)==W(:))
%
%June 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

Pxz=otherInfo.Pxz;

Pzz=PzPred+R;
W=Pxz/Pzz;

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
% inputdlg2() - inputdlg function clone with coloring and help for
%               eeglab().
%
% Usage:
%   >> Answer = inputdlg2(Prompt,Title,LineNo,DefAns,funcname);
%
% Inputs:
%   Same as inputdlg. Using the optional additionnal funcname parameter
%   the function will create a help button. The help message will be
%   displayed using the pophelp() function.
%
% Output:
%   Same as inputdlg
%
% Note: The advantage of this function is that the color of the window
%       can be changed and that it displays an help button. Edit
%       supergui to change window options. Also the parameter LineNo
%       can only be one.
%
% Author: Arnaud Delorme, CNL / Salk Institute, La Jolla, 11 August 2002
%
% See also: supergui(), inputgui()

% Copyright (C) Arnaud Delorme, CNL / Salk Institute, arno@salk.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, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [result] = inputdlg2(Prompt,Title,LineNo,DefAns,funcname);

if nargin < 4
   help inputdlg2;
   return;
end;
if nargin < 5
	funcname = '';
end;

if length(Prompt) ~= length(DefAns)
	error('inputdlg2: prompt and default answer cell array must have the smae size');
end;

geometry = {};
listgui = {};

% determine if vertical or horizontal
% -----------------------------------
geomvert = [];
for index = 1:length(Prompt)
	geomvert = [geomvert size(Prompt{index},1) 1];  % default is vertical geometry
end;
if all(geomvert == 1) & length(Prompt) > 1
	geomvert = []; % horizontal
end;

for index = 1:length(Prompt)
	if ~isempty(geomvert) % vertical
		geometry = { geometry{:} [ 1] [1 ]};
	else
		geometry = { geometry{:} [ 1 0.6 ]};
	end;
	listgui = { listgui{:} { 'Style', 'text', 'string', Prompt{index}}  ...
				{ 'Style', 'edit', 'string', DefAns{index} } };
end;

result = inputgui(geometry, listgui, ['pophelp(''' funcname ''');'], Title, [], 'normal', geomvert);

author: PatternRecognition

repo: OpenBMI


classdef yolov3Layer < nnet.layer.Layer
    % 该自定义类仍有很大不灵活性，表现在1、predict函数输出不能自由控制；2、类属性在predict函数里面无法改变；3、训练阶段forward函数不执行
    % 4、该类无法自定义其他普通成员函数。
    % 参考官方文档：Define Custom Deep Learning Layers
    properties
        mask       % 如当前yolov3层mask为[1，2，3]，matlab中从1开始，对应cfg文件中以0开始的mask
        anchors  % 所有的anchors,形如[width1,height1; weight2,height2; ...]，对应cfg文件中的anchors，大小为allnumAnchors*2
        num % total number of anchors,所有anchors的数量，对应cfg文件中的num
        ignore_thresh % keeps duplicated detections if IoU(detect, truth) >ignore_thresh, which will be fused during NMS (is used for training only)
        truth_thresh % adjusts duplicated detections if IoU(detect, truth) > truth_thresh, which will be fused during NMS (is used for training only)
        jitter % randomly crops and resizes images with changing aspect ratio from x(1 - 2*jitter) to x(1 + 2*jitter) (data augmentation parameter is used only from the last layer)
        random % randomly resizes network for each 10 iterations from 1/1.4 to 1.4(data augmentation parameter is used only from the last layer)
        classes   % 目标检测总共类别数量，对应cfg文件中的classes

        anchorsUse %当前yolov3层的anchors，n*2大小，[width, height]，宽高为相对输入原图大小,注意这个与matlab官网estimateAnchorBoxes的宽高顺序相反
        nAnchors % 当前yolov3层使用anchors的数量，一般为3
        numX          % 传入到该yolov3层特征图的宽，像素
        numY          % 传入到该yolov3层特征图的宽，像素
        imageSize     % 网络输入图像大小，[imageHeight,imageWidth]
        arc           % 方式,取'default'、'uCE'、'uBCE'中的一种
        stride       % 传入到该yolov3层特征下采样率，[w,h],如输入网络大小是416*416宽乘高，该特征层为13*13，则stride = 32
    end

    methods
        function layer = yolov3Layer(name, mask,allAnchors,nClasses,yoloIndex,imageSize,arc)
            % 输入参数：
            % name: 字符向量或字符串，层名字
            % mask: 见上面properties说明
            % allAnchors：同上面properties之anchors说明
            % nClasses： 1*1 标量，检测的类别数量
            % imageSize：[imageHeight,imageWidth],输入网络的图像大小
            % yoloIndex：1*1标量，从1开始的索引，输出yolo检测层的序号
            % arc:字符向量或字符串，计算损失的方式
            %
            % cuixingxing150@gmail.com
            % 2020.4.20
            %
            assert(size(allAnchors,2)==2,'allAnchors must have n*2 shape!');% nAchors*2,形如[width,height]

            layer.mask = mask;
            layer.anchors = allAnchors;
            layer.num = numel(allAnchors)/2;
            layer.classes = nClasses;
            layer.ignore_thresh = .7;
            layer.truth_thresh = 1;
            layer.jitter = .3;
            layer.random = 1;

            anchorsUse = allAnchors(mask,:);% na*2, 注意matlab中mask是从1开始的索引
            layer.Name = name;
            text = ['all number classes:', num2str(nClasses),...
                ',used anchor box:',mat2str(round(anchorsUse)),...
                ', yoloLayerID:',num2str(yoloIndex),...
                ', arc:',arc];
            layer.Description = text;
            layer.Type = 'yoloV3Layer';

            layer.numY = 1;
            layer.numX = 1;
            layer.stride =1;

            layer.anchorsUse = anchorsUse;
            layer.nAnchors = size(anchorsUse,1);
            layer.imageSize = imageSize;
            layer.arc = arc;
        end

         function Z = predict(layer, X) % training
            % matlab中上一层输入的特征X为H*W*C*N，即h*w*c*bs
            X = dlarray(X);
            layer.numY = size(X,1);
            layer.numX = size(X,2);
            bs = size(X,4);
            layer.stride = max(layer.imageSize)./max(layer.numX,layer.numY);

            Z = X;
            % 2020.4.23 实际写下面2行的话，待整个yolov3网络predict时候报错！！！故移到外面写
%             Z = X.reshape(layer.numY,layer.numX,layer.nAnchors,(5+layer.nClasses),bs);
%             Z = permute(Z,[5,3,1,2,4]);% 输出特征图，Z=bs*na*h*w*(5+nc),
         end

    end
end


text:
function [x,t] = ToolGenerateSineWave(fFreq, fLengthInS, fSampleRateInHz)

    [m n]   = size(fFreq);
    if (min(m,n)~=1)
        error('illegal frequency dimension')
    end
    if (m<n)
        fFreq   = fFreq';
    end

    t = linspace(0,fLengthInS-1/fSampleRateInHz,fSampleRateInHz*fLengthInS);

    x = sin(2*pi*fFreq*t);
end

author: alexanderlerch

repo: ACA-Code


text:
function pass = test_divide(pref)
% Test contour

if ( nargin == 0)
    pref = chebfunpref;
end
tol = 1000*pref.cheb3Prefs.chebfun3eps;

f = chebfun3(@(x,y,z) cos(x.*y.*z));
g = chebfun3(@(x,y,z) cos(x.*y.*z)./2);

pass(1) = norm( f./2 - g ) < tol;

pass(2) = norm( f/2 - g ) < tol;

pass(3) = norm( 2.\f - g ) < tol;

pass(4) = norm( 2\f - g ) < tol;

end


text:
function showFeatures(I, pts,  ptSyle, ptWidth)
%PLOTFEATURES Plot features using MATLAB
%
% INPUT:
%   - I(image): image corresponding to the points
%   - pts(1, N): array of structure conatining keypoints
%   - ptSyle: defines color and type of scatter plot; eg: '+r'
%   - ptWidth: defines thickness of scatter points

assert (size(pts, 2) > size(pts, 1), 'Input expected to array of shape (1, N)')

% get locations of points
loc_x = zeros(1, size(pts, 2));
loc_y = zeros(1, size(pts, 2));
class_pt = zeros(1, size(pts, 2));

index = 1;
for point = pts
    index = index + 1;
    loc_x(index) = point.location(1);
    loc_y(index) = point.location(2);
    class_pt(index) = point.class;
end

color = ['r', 'b', 'g', 'y'];

% plot image
imshow(I);
hold on;
for i = 1:4
    scatter(loc_y(class_pt == i), loc_x(class_pt == i), ptSyle, color(i), 'LineWidth', ptWidth);
end
legend('Blob mimimum', 'Blob maximum', 'Corner minimum', 'Corner maximum')

end

repo: Stereo-Odometry-SOFT


text:
function test_symm_clustering()
%
% demonstration file for NMFLibrary.
%
% This file illustrates how to use this library.
% This demonstrates Symm-ANLS algorithm and Symm-Newton algorithm.
%
% This file is part of NMFLibrary.
%
% Created by H.Kasai on Jun. 26, 2019

    clc;
    clear;
    close all;

    %% generate synthetic data of (mxn) matrix
    if 1
        input = importdata('../../data/ORL.mat');
        M = input.data;
        gnd = input.label;
    else
        input = importdata('../../data/COIL20.mat');
        M = (input.TrainSet.X)';
        gnd = (input.TrainSet.y)';
    end
    clear input;

    V = calcu_similarity_matrix(M);
    rank = length(unique(gnd));

    %W_init = 2 * full(sqrt(mean(mean(V)) / rank)) * rand(m, rank);
    %W_init = rand(m, rank);
    %options.x_init.W = W_init;
    %options.x_init.H = (options.x_init.W)';

    lambda = 0.66;

    %% initialize rank to be factorized
    options.verbose = 2;
    options.max_epoch = 100;
    options.calc_symmetry = true;
    options.calc_clustering_acc = true;
    options.clustering_gnd = gnd;
    options.clustering_classnum = rank;
    %options.clustering_eval_num = 10;
    %options.init_alg = 'symm_mean';

    %% perform factroization
    % Symm-ANLS
    options.alpha = lambda;
    [w_symm_anls, infos_symm_anls] = symm_anls(V, rank, options);
    % Symm-Newton
    [w_symm_newton, infos_symm_newton] = symm_newton(V, rank, options);
    % Symm-Hals
    options.lambda = lambda;
    [w_symm_halsacc, infos_symm_halsacc] = symm_halsacc(V, rank, options);

    %% plot
    display_graph('epoch','cost', {'Symm-ANLS', 'Symm-Newton', 'Symm-HALS'}, ...
        {w_symm_anls, w_symm_newton, w_symm_halsacc}, {infos_symm_anls, infos_symm_newton, infos_symm_halsacc});
    display_graph('time','cost', {'Symm-ANLS', 'Symm-Newton', 'Symm-HALS'}, ...
        {w_symm_anls, w_symm_newton, w_symm_halsacc}, {infos_symm_anls, infos_symm_newton, infos_symm_halsacc});

    %% symmetry
    display_graph('epoch','symmetry', {'Symm-ANLS', 'Symm-Newton', 'Symm-HALS'}, ...
        {w_symm_anls, w_symm_newton, w_symm_halsacc}, {infos_symm_anls, infos_symm_newton, infos_symm_halsacc});

    %% clustering
    display_graph('epoch','clustering_acc', {'Symm-ANLS', 'Symm-Newton', 'Symm-HALS'}, ...
        {w_symm_anls, w_symm_newton, w_symm_halsacc}, {infos_symm_anls, infos_symm_newton, infos_symm_halsacc});
    display_graph('epoch','clustering_nmi', {'Symm-ANLS', 'Symm-Newton', 'Symm-HALS'}, ...
        {w_symm_anls, w_symm_newton, w_symm_halsacc}, {infos_symm_anls, infos_symm_newton, infos_symm_halsacc});
    display_graph('epoch','clustering_purity', {'Symm-ANLS', 'Symm-Newton', 'Symm-HALS'}, ...
        {w_symm_anls, w_symm_newton, w_symm_halsacc}, {infos_symm_anls, infos_symm_newton, infos_symm_halsacc});

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function [VG,A,PPG]= colorgrad(f,T)
%COLORGRAD Computes the vector gradient of an RGB image.
%   [VG,VA,PPG] = COLORGRAD(F,T) computes the vector gradient, VG, and
%   corresponding angle array, VA, (in radians) of RGB image F. It also
%   computes PPG, the per-plane composite gradient obtained by summing
%   the 2-D gradients of the individual color planes. Input T is a
%   threshold in the range [0,1]. If it is included in the input,
%   outputs VG and PPG are logical matrices such that VG(x,y) = 0 for
%   values <= T and VG(x,y) = 1 otherwise. Similar comments apply to
%   PPG. If T is not included in the input argument, then both output
%   gradients are scaled to the double range [0,1].
%
%  All equations on which this function is based are from Chapter 7 of
%  DIPUM3E.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

if size(f,3) ~= 3
   error('Input image must be of size M-by-N-by-3');
end

if nargin == 2 && (T < 0 || T > 1)
   error('T must be in the range [0 1]')
end

% Compute the x and y derivatives of the three component images using
% Sobel operators.
sh = fspecial('sobel');
sv = sh';
Rx = imfilter(double(f(:,:,1)),sh,'replicate');
Ry = imfilter(double(f(:,:,1)),sv,'replicate');
Gx = imfilter(double(f(:,:,2)),sh,'replicate');
Gy = imfilter(double(f(:,:,2)),sv,'replicate');
Bx = imfilter(double(f(:,:,3)),sh,'replicate');
By = imfilter(double(f(:,:,3)),sv,'replicate');

% Compute the parameters of the vector gradient.
gxx = Rx.^2 + Gx.^2 + Bx.^2;
gyy = Ry.^2 + Gy.^2 + By.^2;
gxy = Rx.*Ry + Gx.*Gy + Bx.*By;
A = 0.5*(atan(2*gxy./(gxx - gyy + eps)));
G1 = 0.5*((gxx + gyy) + (gxx - gyy).*cos(2*A) + 2*gxy.*sin(2*A));

% Now repeat for angle + pi/2. Then select the maximum at each point.
A = A + pi/2;
G2 = 0.5*((gxx + gyy) + (gxx - gyy).*cos(2*A) + 2*gxy.*sin(2*A));
G1 = G1.^0.5;
G2 = G2.^0.5;

% Form VG by picking the maximum at each (x,y) and then scale
% to the range [0,1] using function mat2gray.
VG = mat2gray(max(G1,G2));

% Select the corresponding angles. Where G2 is greater than G1, pick
% values from array A + pi/2. Else, use values from A, which correspond
% to G1. The angles are in radians.
A(G2 > G1) = A(G2 > G1) + pi/2;

% Compute the per-plane gradients.
RG = hypot(Rx,Ry);
GG = hypot(Gx,Gy);
BG = hypot(Bx,By);
% Form the composite by adding the individual results and scale to
% [0,1].
PPG = mat2gray(RG + GG + BG);

% If T was provided, threshold the result.
if nargin == 2
   VG = VG > T;
   PPG = PPG > T;
end

author: dipum

repo: dipum-toolbox

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


text:
function [F,sE,sC] = VBA_spm_log_evidence(varargin)
% Return the log-evidence of a reduced model (under Laplace approximation)
% FORMAT [F,sE,sC] = spm_log_evidence(qE,qC,pE,pC,rE,rC)
% FORMAT [F,sE,sC] = spm_log_evidence(qE,qC,pE,pC,priorfun,varargin)
% FORMAT [F,sE,sC] = spm_log_evidence(qE,qC,pE,pC)
%
% qE,qC    - posterior expectation and covariance of full model
% pE,pC    - prior expectation and covariance of full model
% rE,rC    - prior expectation and covariance of reduced model
% or
% priorfun - inline function that returns prior moments
%            {rE rC} = priorfun(varargin{:})
%
% or (if omitted) rE = 0 and rC = 0;
%
% F        - reduced log-evidence: ln p(y|reduced model) - ln p(y|full model)
% [sE,sC]  - posterior expectation and covariance of reduced model
%
%--------------------------------------------------------------------------
% This routine assumes the reduced model is nested within a full model and
% that the posteriors (and priors) are Gaussian. Nested here means that the
% prior precision of the reduced model, minus the prior precision of the
% full model is positive definite. We additionally assume that the prior
% means are unchanged. The two input argument formats are for use with
% spm_argmax.
%__________________________________________________________________________
% Copyright (C) 2005 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_log_evidence.m 4281 2011-03-31 19:49:57Z karl $

% Compute reduced log-evidence
%==========================================================================

% check to see if priors are specified by a function
%--------------------------------------------------------------------------
qE = varargin{1};
qC = varargin{2};
pE = varargin{3};
pC = varargin{4};
try
    priors = varargin{5}(varargin{6:end});
    rE     = priors{1};
    rC     = priors{2};
catch
    try
        rE = varargin{5};
        rC = varargin{6};
    catch
        n  = size(qC,1);
        rE = sparse(n,1);
        rC = sparse(n,n);
    end
end

% reduced subspace
%--------------------------------------------------------------------------
qE  = VBA_spm_vec(qE);
pE  = VBA_spm_vec(pE);
rE  = VBA_spm_vec(rE);

if nargout < 2
    dE  = pE - rE;
    dC  = pC - rC;
    k   = find(dE | any(dC,2));
    if ~isempty(k)
        qE  = qE(k);
        pE  = pE(k);
        rE  = rE(k);
        qC  = qC(k,k);
        pC  = pC(k,k);
        rC  = rC(k,k);
    else

        % the reduced and full models are the same
        %------------------------------------------------------------------
        F   = 0;
        return
    end
end

% fix tolerance for matrix inversions
%--------------------------------------------------------------------------
TOL   = exp(-16);

% remove fixed parameters under full model
%--------------------------------------------------------------------------
i     = find(diag(pC));

% preliminaries
%--------------------------------------------------------------------------
qP    = VBA_spm_inv(qC(i,i),TOL);
pP    = VBA_spm_inv(pC(i,i),TOL);
rP    = VBA_spm_inv(rC(i,i),TOL);
sP    = qP + rP - pP;
sC    = VBA_spm_inv(sP,TOL);
sE    = qP*qE(i) + rP*rE(i) - pP*pE(i);

% log-evidence
%--------------------------------------------------------------------------
F     = VBA_spm_logdet(rP*qP*sC*pC(i,i)) ...
      - (qE(i)'*qP*qE(i) + rE(i)'*rP*rE(i) - pE(i)'*pP*pE(i) - sE'*sC*sE);
F     = F/2;

% restore full conditional density
%--------------------------------------------------------------------------
if nargout > 1
    rE(i)   = sC*sE;
    rC(i,i) = sC;
    sE      = VBA_spm_unvec(rE,varargin{1});
    sC      = rC;
end

author: MBB-team

repo: VBA-toolbox


text:
function [str,iter]=display(comp,tab,iter)

  global display_tree_array_indexing;
  if isempty(display_tree_array_indexing)
      display_tree_array_indexing=2;
  end;
  global display_tree_index_at_front;
  if isempty(display_tree_index_at_front)
      display_tree_index_at_front=1;
  end;

  if nargin==3
   if isempty(iter)
    iter{1}=[];
   else
      i=length(iter);
      iter{i+1}=iter{i};  %% <--- new member of tree
   end
   if display_tree_array_indexing==0
     if isempty(iter{length(iter)})
      index=[];
     else
      index=[num2str(length(iter)-1) ':'];
     end
     if display_tree_index_at_front
      str=[ index get_name(comp) '\n'];
     else
      str=[ get_name(comp) ' (' index(length(index))  ')\n'];
     end;
   else
     index=[];
     if ~isempty(iter{length(iter)})
      index= num2str(iter{length(iter)});

      a=findstr('  ',index);
      a=setdiff([1:length(index)],a); index=index(a);

      if display_tree_array_indexing==2
    	a=find(index==32);
        if isempty(a)
            a=0;
        end;
	    index=index(max(a)+1:length(index));
      else
        index=['[' index ']'];
      end
      index=[index ':'];
     end;
     if display_tree_index_at_front
       str=[ index  get_name(comp) '\n'];
     else
       str=[ get_name(comp)  ' ' index(1:length(index)-1)  '\n'];
     end;
   end

  else
   disp([get_name(comp)]);
  end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function check = sudoku_check ( s )

%*****************************************************************************80
%
%% SUDOKU_CHECK checks a partial or filled-in Sudoku puzzle.
%
%  Discussion:
%
%    The routine ensures that
%    1) each entry of S is either 0 (unfilled) or between 1 and 9.
%    2) each row contains no more than one occurrence of a digit.
%    3) each column contains no more than one occurrent of a digit.
%    4) each box contains no more than one occurrence of a digit.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer S(9,9) contains the Sudoku puzzle.  Unfilled
%    entries should be set to 0.
%    .
%    Output, integer CHECK, is
%    0, if no errors were detected.
%        10*I+J, if A(I,J) contains an illegal digit,
%    100+10*I+J, if A(I,J) violates the row condition.
%    200+10*I+J, if A(I,J) violates the column condition.
%    300+10*I+J, if A(I,J) violates the box condition.
%

%
%  Digit check.
%
  i = 0;
  for i3 = 1 : 3
    for i2 = 1 : 3
      i = i + 1;
      j = 0;
      for j3 = 1 : 3
        for j2 = 1 : 3
          j = j + 1;
          if ( s(i,j) < 0 | 9 < s(i,j) )
            check = 10 * i + j;
            return
          end
        end
      end
    end
  end
%
%  Row check.
%
  i = 0;
  for i3 = 1 : 3
    for i2 = 1 : 3
      i = i + 1;
      digit(1:9) = 0;
      j = 0;
      for j3 = 1 : 3
        for j2 = 1 : 3
          j = j + 1;
          k = s(i,j);
          if ( 1 <= k & k <= 9 )
            digit(k) = digit(k) + 1;
            if ( 1 < digit(k) )
              check = 100 + 10 * i + j;
              return
            end
          end
        end
      end
    end
  end
%
% Column check.
%
  j = 0;
  for j3 = 1 : 3
    for j2 = 1 : 3
      j = j + 1;
      digit(1:9) = 0;
      i = 0;
      for i3 = 1 : 3
        for i2 = 1 : 3
          i = i + 1;
          k = s(i,j);
          if ( 1 <= k & k <= 9 )
            digit(k) = digit(k) + 1;
            if ( 1 < digit(k) )
              check = 200 + 10 * i + j;
              return
            end
          end
        end
      end
    end
  end
%
%  Box check
%
  for i3 = 1 : 3
    for j3 = 1 : 3
      digit(1:9) = 0;
      i = ( i3 - 1 ) * 3;
      for i2 = 1 : 3
        j = ( j3 - 1 ) * 3;
        i = i + 1;
        for j2 = 1 : 3
          j = j + 1;
          k = s(i,j);
          if ( 1 <= k & k <= 9 )
            digit(k) = digit(k) + 1;
            if ( 1 < digit(k) )
              check = 300 + 10 * i + j;
              return
            end
          end
        end
      end
    end
  end
%
%  No errors discovered.
%
  check = 0;

  return
end


text:
function bytes = get_filesize(file_path)
%GET_FILESIZE Returns the size of the specified file in bytes.
% This is a wrapper for dir().
%
% Usage:
%   bytes = get_filesize(file_path)
%
% See also: dir

% Get file attributes
attributes = dir(file_path);
bytes = attributes.bytes;

if nargout == 0
    printf('%s: *%s*', get_filename(file_path), bytes2str(bytes))
    clear bytes
end

end

author: talmo

repo: leap

save_path: github-repos/MATLAB/talmo-leap


text:
% Test file for chebtech/clenshaw.m

function pass = test_clenshaw(varargin)

% Set a tolerance (pref.chebfuneps doesn't matter)
tol = 10*eps;

%%
% Test that a single coefficient is evaluated correctly:
% For a scalar evaluation:
c = sqrt(2);
v = chebtech.clenshaw(0, c);
pass(1) = c == v;

% For a vector evaluation:
x = [-.5 ; 1];
v = chebtech.clenshaw(x, c);
pass(2) = ( all(size(v) == [2, 1]) && all(c == v) );

% For a row vector evaluation with column coefficients:
c = [c, c, c];
v = chebtech.clenshaw(x, c);
pass(3) = ( all(size(v) == [2, 3]) && norm(repmat(c, 2, 1) - v) == 0);

%%
% Test that a vector coefficient is evaluated correctly:
% Some simple data :
c = (5:-1:1).';
x = [-.5 ; -.1 ; 1];

% Scalar coefficient
v = chebtech.clenshaw(x, c);
% Exact values:
vTrue = [3 ; 3.1728 ; 15];
pass(4) = norm(v - vTrue, inf) < tol;

% In vectorised form:
vTrue2 = [0 ; 3.6480 ; 15];
v = chebtech.clenshaw(x, [c, c(end:-1:1)]);
pass(5) = norm(v - [vTrue, vTrue2], inf) < tol;

end


text:
function value = zabs2 ( z )

%*****************************************************************************80
%
%% ZABS2 returns the L2 norm of a complex number.
%
%  Discussion:
%
%    The L2 norm of a complex number is the square root of the sum
%    of the squares of the real and imaginary components.
%
%    ZABS2 ( Z ) = sqrt ( ( real ( Z ) )**2 + ( imaginary ( Z ) )**2 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 May 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979.
%
%    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
%    Basic Linear Algebra Subprograms for FORTRAN usage,
%    ACM Transactions on Mathematical Software,
%    Volume 5, Number 3, pages 308-323, 1979.
%
%  Parameters:
%
%    Input, complex Z, the number whose norm is desired.
%
%    Output, real VALUE, the L2 norm of Z.
%
  value = sqrt ( ( real ( z ) ).^2 + ( imag ( z ) ).^2 );

  return
end


text:
function [K,MK,data] = hks(V,F,varargin)
  % HKS Heat Kernel Signature "A Concise and Provably Informative Multi-Scale
  % Signature Based on Heat Diffusion" [Sun et al. 2009]
  %
  % Inputs:
  %   V  #V by dim list of mesh vertex positions
  %   F  #F by ss list of element indices into V
  % Outputs:
  %   K  #V by #K list of HKS feature vectors for each vertex
  %   MK  #K list of mass-weighted spatial sums
  %   data
  %     lambda  #lambda list of eigen values
  %     phi  #V by #lambda list of eigen functions
  %     phi_sqr  #V by #lambda list of squared eigen functions
  %     exp_m_lambda  #lambda list of exp(- eigen values)
  %     M  #V by #V mass matrix
  %     L  #V by #V Laplacian matrix
  %
  % Example:
  %   [K,MK] = hks(V,F);
  %   % Compare all vertices' signatures to first using metric of [Sun et al.
  %   % 2009]
  %   tsurf(F,V,'CData',sqrt(sum(((K-K(1,:))./MK).^2,2)),fphong);

  tmin = [];
  tmax = [];
  data.t = [];
  data.exp_m_lambda = [];
  data.phi_sqr = [];
  data.M = [];
  k = 300;
  nt = 100;
  data = [];
  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'tmax','tmin','Data','NumModes','NumTimeSteps'}, ...
    {'tmax','tmin','data','k','nt'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  if isempty(data)
    data.L = cotmatrix(V,F);
    data.M = massmatrix(V,F);
    [data.phi,data.lambda] = eigs(-data.L,data.M,k,'sm');
    data.lambda = diag(data.lambda);
    data.phi_sqr = data.phi.^2;
    data.exp_m_lambda = exp(-data.lambda);
  end
  if isempty(tmin)
    tmin = 4*log(10)/max(data.lambda);
  end
  if isempty(tmax)
    tmax = 4*log(10)/min(data.lambda(data.lambda>min(data.lambda)));
  end

  data.t = exp(linspace(log(tmin),log(tmax),nt));

  K = (data.phi_sqr)*(data.exp_m_lambda.^data.t);
  MK = diag(data.M)'*K;

end

author: alecjacobson

repo: gptoolbox


text:
function varargout = ellipseAsPolygon(ellipse, N)
%ELLIPSEASPOLYGON Convert an ellipse into a series of points
%
%   Deprecated, use ellipseToPolygon instead.
%
%   P = ellipseAsPolygon(ELL, N);
%   converts ELL given as [x0 y0 a b] or [x0 y0 a b theta] into a polygon
%   with N edges. The result P is (N+1)-by-2 array containing coordinates
%   of the N+1 vertices of the polygon.
%   The resulting polygon is closed, i.e. the last point is the same as the
%   first one.
%
%   P = ellipseAsPolygon(ELL);
%   Use a default number of edges equal to 72. This result in one piont for
%   each 5 degrees.
%
%   [X Y] = ellipseAsPolygon(...);
%   Return the coordinates o fvertices in two separate arrays.
%
%   See also:
%   ellipses2d, circleAsPolygon, rectAsPolygon, drawEllipse
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 06/04/2005.
%

%   HISTORY
%   2011-03-30 use angles in degrees, add default value for N
%   2011-12-09 deprecate

warning('matGeom:deprecated', ...
    'function "ellipseAsCurve" is deprecated, use "ellipseToPolygon" instead');

% format output
if nargout <= 1
    varargout = {ellipseToPolygon(ellipse, N)};
else
    [x, y] = ellipseToPolygon(ellipse, N);
    varargout = {x, y};
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [gx,dgdx] = g_RFX(x,P,u,in)

gx = in.X*x;
dgdx = in.X';


text:
function [foldChange,standardError,solControl,solCondition] = eFlux(model,controlExpression,conditionExpression,varargin)
% Calculate the objective fold change according to the eFlux approach for
% expression integration as described in:
% Interpreting Expression Data with Metabolic Flux Models: Predicting Mycobacterium tuberculosis Mycolic Acid Production
% Colijn C, Brandes A, Zucker J, Lun DS, Weiner B, et al. (2009)
% PLOS Computational Biology 5(8): e1000489. https://doi.org/10.1371/journal.pcbi.1000489
%
% USAGE:
%    [foldChange,standardError] = eFlux(model,controlExpression,conditionExpression,varargin)
%
% INPUTS:
%    model:                  The COBRA model struct to use
%    controlExpression:      struct for the control expression with two fields required and one optional field:
%                              * .target       - the names of the target (rxns or genes)
%                              * .value        - the value for the target. Positive values for all constraint reactions, negative values for unconstraint reactions.
%                              * .preprocessed - Indicator whether the provided targets are genes (false), or reactions (true) Default: false
%    conditionExpression:    struct for the condition expression (fields are the same as controlExpression)
% OPTIONAL INPUT:
%    varargin:               parameters given as struct or parameter/value pairs.
%                             * testNoise - indicator whether to run multiple calculations with added noise to get a significance of the fold change. Requires either a noise function and a standard deviation of noise ('noiseFun' and 'noiseStd' respectively) or a controlData struct and a noise function.
%                             * noiseCount - number of noisy controls to create if noise is tested (default: 10)
%                             * noiseFun - The noise function to use, has to be a function handle taking 2 arguments (mean and std)
%                             * noiseStd - The standard deviation(s) to use. Either a single value (used for all values, or a vector with length equal to controlExpression.value).
%                             * controlData - a struct (like controlExpression which has a value matrix with multiple values per controlExpression to determine the noise distribution. If provided with testNoise == false, the values from this struct will be used to determine the noise.
%                             * minSum:           Switch for the processing of Genetic data. If false, ORs in the GPR will be treated as min. If true(default), ORs will be treated as addition.
%                             * softBounds:       Whether to use soft bounds for the infered constraints or to add flexibility variables (default: false).
%                             * weightFactor:     The weight factor for soft bounds (default: 1)
% OUTPUTS:
%    foldChange:             The fold change between the objective of the
%                            condition and the objective of the control
%                            expression
%    standardError:          The error if noise is being used.
%    solControl:             The solution of the given Control expression;
%    solCondition:           The solution of the given Condition expression;
%
% ..Author:    Thomas Pfau OCt 2018
%
% NOTE:
%    This si an implementation of the eFlux concept as presented in:
%    Interpreting Expression Data with Metabolic Flux Models: Predicting Mycobacterium tuberculosis Mycolic Acid Production
%    Colijn C, Brandes A, Zucker J, Lun DS, Weiner B, et al. (2009)
%    PLOS Computational Biology 5(8): e1000489. https://doi.org/10.1371/journal.pcbi.1000489
%    Please note, that this code does not perform any preprocessing expcept
%    for that described in the above paper after array normalization.

normFun = @(mean,std) normrnd(mean,std);

if ~isfield(controlExpression,'preprocessed')
    controlExpression.preprocessed = true;
end

parser = inputParser();
parser.addRequired('model',@(x) verifyModel(x,'simpleCheck',true));
parser.addRequired('controlExpression',@(x) verifyEFluxExpressionStruct(model,x));
parser.addRequired('conditionExpression',@(x) verifyEFluxExpressionStruct(model,x));
parser.addParameter('testNoise',false, @(x) islogical(x) || (isnumeric(x) && (x==1 || x==0)));
parser.addParameter('noiseCount',10, @isnumeric);
parser.addParameter('noiseFun',normFun, @(x) isa(x, 'function_handle'));
parser.addParameter('noiseStd',[], @isnumeric);
parser.addParameter('controlData',[], @(x) verifyEFluxExpressionStruct(model,x));
parser.addParameter('minSum',false,@(x) islogical(x) || (isnumeric(x) && (x==1 || x==0)));
parser.addParameter('softBounds',false,@(x) islogical(x) || (isnumeric(x) && (x==1 || x==0)));
parser.addParameter('weightFactor',1,@isnumeric);

parser.parse(model,controlExpression,conditionExpression,varargin{:});
testNoise = parser.Results.testNoise;
noiseCount = parser.Results.noiseCount;
noiseStd = parser.Results.noiseStd;
noiseFun = parser.Results.noiseFun;
controlData = parser.Results.controlData;
noisyControl = [];
% remove enforcing bounds (contradict eFlux concept)
if(any(model.lb > 0 | model.ub < 0))
    warning('Enforcing bounds for the following fluxes have been removed:\n%s', strjoin(model.rxns((model.lb > 0 | model.ub < 0)),'\n'));
    model.lb(model.lb > 0) = 0;
    model.ub(model.ub < 0) = 0;
end

if testNoise
    noisyControl = repmat(columnVector(controlExpression.value),1,noiseCount);
    if isempty(controlData) && isempty(noiseStd)
        error('To test noise, either a standard deviation, or appropriate controlData has to be provided')
    end
    if ~isempty(controlData)
        if ~isempty(noiseStd)
            error('To test noise, either a standard deviation, or appropriate controlData has to be provided but not both.')
        end
        % ok, we use the controlData and create noise based on its
        % standarddeviations.
        noiseStd = std(controlData.values');
    else
        if numel(noiseStd) == 1
            noiseStd = repmat(noiseStd,size(noisyControl,1),1);
        end
        controlData.target = controlExpression.target;
        controlData.preprocessed = controlExpression.preprocessed;
    end
    for i = 1:noiseCount
            % add the noise.
            noise = arrayfun(@(x) noiseFun(0,x),columnVector(noiseStd));
            noisyControl(:,i) = noisyControl(:,i) + columnVector(noise);
    end
else
    if ~isempty(controlData)
        noisyControl = controlData.values;
    end
end

% calculate the condition solution
conditionModel = applyEFluxConstraints(model,conditionExpression,varargin{:});
solCondition = optimizeCbModel(conditionModel);
conditionValue =  solCondition.f;

% calculate the default solution (for controlExpression
defaultControl = applyEFluxConstraints(model,controlExpression,varargin{:});
solControl = optimizeCbModel(defaultControl);
controlValue = zeros(1+size(noisyControl,2),1);
controlValue(1) = solControl.f;

for noisy = 1:size(noisyControl,2)
    % calculcate the noisy solutions
    if isfield(controlData,'preprocessed')
        exprStruct.target = controlData.target;
        exprStruct.value = noisyControl(:,noisy);
        exprStruct.preprocessed = controlData.preprocessed;
    else
        exprStruct.target = controlData.target;
        exprStruct.value = noisyControl(:,noisy);
    end
    defaultControl = applyEFluxConstraints(model,exprStruct,varargin{:});
    controlSol = optimizeCbModel(defaultControl);
    controlValue(1+noisy) = controlSol.f;
end
% set the outputs
foldChanges = conditionValue./controlValue;
foldChange = foldChanges(1);
standardError = std(foldChanges) / sqrt(numel(foldChanges));

end

author: opencobra

repo: cobratoolbox


text:
function [traj, infStates] = tapas_hgf_ar1(r, p, varargin)
% Calculates the trajectories of the agent's representations under the HGF
%
% This function can be called in two ways:
%
% (1) tapas_hgf_ar1(r, p)
%
%     where r is the structure generated by tapas_fitModel and p is the parameter vector in native space;
%
% (2) tapas_hgf_ar1(r, ptrans, 'trans')
%
%     where r is the structure generated by tapas_fitModel, ptrans is the parameter vector in
%     transformed space, and 'trans' is a flag indicating this.
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Transform paramaters back to their native space if needed
if ~isempty(varargin) && strcmp(varargin{1},'trans');
    p = tapas_hgf_ar1_transp(r, p);
end

% Number of levels
try
    l = r.c_prc.n_levels;
catch
    l = length(p)/6;

    if l ~= floor(l)
        error('tapas:hgf:UndetNumLevels', 'Cannot determine number of levels');
    end
end

% Unpack parameters
mu_0 = p(1:l);
sa_0 = p(l+1:2*l);
phi  = p(2*l+1:3*l);
m    = p(3*l+1:4*l);
ka   = p(4*l+1:5*l-1);
om   = p(5*l:6*l-2);
th   = exp(p(6*l-1));
al   = p(6*l);

% Add dummy "zeroth" trial
u = [0; r.u(:,1)];

% Number of trials (including prior)
n = length(u);

% Assume that if u has more than one column, the last contains t
try
    if r.c_prc.irregular_intervals
        if size(u,2) > 1
            t = [0; r.u(:,end)];
        else
            error('tapas:hgf:InputSingleColumn', 'Input matrix must contain more than one column if irregular_intervals is set to true.');
        end
    else
        t = ones(n,1);
    end
catch
    if size(u,2) > 1
        t = [0; r.u(:,end)];
    else
        t = ones(n,1);
    end
end

% Initialize updated quantities

% Representations
mu = NaN(n,l);
pi = NaN(n,l);

% Other quantities
muhat = NaN(n,l);
pihat = NaN(n,l);
v     = NaN(n,l);
w     = NaN(n,l-1);
da    = NaN(n,l);
dau   = NaN(n,1);

% Representation priors
% Note: first entries of the other quantities remain
% NaN because they are undefined and are thrown away
% at the end; their presence simply leads to consistent
% trial indices.
mu(1,:) = mu_0;
pi(1,:) = 1./sa_0;

% Representation update loop
% Pass through trials
for k = 2:1:n
    if not(ismember(k-1, r.ign))

        %%%%%%%%%%%%%%%%%%%%%%
        % Effect of input u(k)
        %%%%%%%%%%%%%%%%%%%%%%

        % 1st level
        % ~~~~~~~~~
        % Prediction
        muhat(k,1) = mu(k-1,1) +t(k) *phi(1) *(m(1) -mu(k-1,1));

        % Precision of prediction
        pihat(k,1) = 1/(1/pi(k-1,1) +t(k) *exp(ka(1) *mu(k-1,2) +om(1)));

        % Input prediction error
        dau(k) = u(k) -muhat(k,1);

        % Updates
        pi(k,1) = pihat(k,1) +1/al;
        mu(k,1) = muhat(k,1) +1/pihat(k,1) *1/(1/pihat(k,1) +al) *dau(k);

        % Volatility prediction error
        da(k,1) = (1/pi(k,1) +(mu(k,1) -muhat(k,1))^2) *pihat(k,1) -1;

        if l > 2
            % Pass through higher levels
            % ~~~~~~~~~~~~~~~~~~~~~~~~~~
            for j = 2:l-1
                % Prediction
                muhat(k,j) = mu(k-1,j) +t(k) *phi(j) *(m(j) -mu(k-1,j));

                % Precision of prediction
                pihat(k,j) = 1/(1/pi(k-1,j) +t(k) *exp(ka(j) *mu(k-1,j+1) +om(j)));

                % Weighting factor
                v(k,j-1) = t(k) *exp(ka(j-1) *mu(k-1,j) +om(j-1));
                w(k,j-1) = v(k,j-1) *pihat(k,j-1);

                % Updates
                pi(k,j) = pihat(k,j) +1/2 *ka(j-1)^2 *w(k,j-1) *(w(k,j-1) +(2 *w(k,j-1) -1) *da(k,j-1));

                if pi(k,j) <= 0
                    error('tapas:hgf:NegPostPrec', 'Negative posterior precision. Parameters are in a region where model assumptions are violated.');
                end

                mu(k,j) = muhat(k,j) +1/2 *1/pi(k,j) *ka(j-1) *w(k,j-1) *da(k,j-1);

                % Volatility prediction error
                da(k,j) = (1/pi(k,j) +(mu(k,j) -muhat(k,j))^2) *pihat(k,j) -1;
            end
        end

        % Last level
        % ~~~~~~~~~~
        % Prediction
        muhat(k,l) = mu(k-1,l) +t(k) *phi(l) *(m(l) -mu(k-1,l));

        % Precision of prediction
        pihat(k,l) = 1/(1/pi(k-1,l) +t(k) *th);

        % Weighting factor
        v(k,l)   = t(k) *th;
        v(k,l-1) = t(k) *exp(ka(l-1) *mu(k-1,l) +om(l-1));
        w(k,l-1) = v(k,l-1) *pihat(k,l-1);

        % Updates
        pi(k,l) = pihat(k,l) +1/2 *ka(l-1)^2 *w(k,l-1) *(w(k,l-1) +(2 *w(k,l-1) -1) *da(k,l-1));

        if pi(k,l) <= 0
            error('tapas:hgf:NegPostPrec', 'Negative posterior precision. Parameters are in a region where model assumptions are violated.');
        end

        mu(k,l) = muhat(k,l) +1/2 *1/pi(k,l) *ka(l-1) *w(k,l-1) *da(k,l-1);

        % Volatility prediction error
        da(k,l) = (1/pi(k,l) +(mu(k,l) -muhat(k,l))^2) *pihat(k,l) -1;
    else

        mu(k,:) = mu(k-1,:);
        pi(k,:) = pi(k-1,:);

        muhat(k,:) = muhat(k-1,:);
        pihat(k,:) = pihat(k-1,:);

        v(k,:)  = v(k-1,:);
        w(k,:)  = w(k-1,:);
        da(k,:) = da(k-1,:);

    end
end

% Remove representation priors
mu(1,:)  = [];
pi(1,:)  = [];

% Check validity of trajectories
if any(isnan(mu(:))) || any(isnan(pi(:)))
    error('tapas:hgf:VarApproxInvalid', 'Variational approximation invalid. Parameters are in a region where model assumptions are violated.');
else
    % Check for implausible jumps in trajectories
    dmu = diff(mu);
    dpi = diff(pi);
    rmdmu = repmat(sqrt(mean(dmu.^2)),length(dmu),1);
    rmdpi = repmat(sqrt(mean(dpi.^2)),length(dpi),1);

    jumpTol = 256;
    if any(abs(dmu(:)) > jumpTol*rmdmu(:)) || any(abs(dpi(:)) > jumpTol*rmdpi(:))
        error('tapas:hgf:VarApproxInvalid', 'Variational approximation invalid. Parameters are in a region where model assumptions are violated.');
    end
end

% Remove other dummy initial values
muhat(1,:) = [];
pihat(1,:) = [];
v(1,:)     = [];
w(1,:)     = [];
da(1,:)    = [];
dau(1)     = [];

% Create result data structure
traj = struct;

traj.mu     = mu;
traj.sa     = 1./pi;

traj.muhat  = muhat;
traj.sahat  = 1./pihat;

traj.v      = v;
traj.w      = w;
traj.da     = da;
traj.dau    = dau;

% Updates with respect to prediction
traj.ud = mu -muhat;

% Psi (precision weights on prediction errors)
psi        = NaN(n-1,l);
psi(:,1)   = 1./(al*pi(:,1));
psi(:,2:l) = pihat(:,1:l-1)./pi(:,2:l);
traj.psi   = psi;

% Epsilons (precision-weighted prediction errors)
epsi        = NaN(n-1,l);
epsi(:,1)   = psi(:,1) .*dau;
epsi(:,2:l) = psi(:,2:l) .*da(:,1:l-1);
traj.epsi   = epsi;

% Full learning rate (full weights on prediction errors)
wt        = NaN(n-1,l);
wt(:,1)   = psi(:,1);
wt(:,2:l) = 1/2 *(v(:,1:l-1) *diag(ka(1:l-1))) .*psi(:,2:l);
traj.wt   = wt;

% Create matrices for use by the observation model
infStates = NaN(n-1,l,4);
infStates(:,:,1) = traj.muhat;
infStates(:,:,2) = traj.sahat;
infStates(:,:,3) = traj.mu;
infStates(:,:,4) = traj.sa;

return;

author: translationalneuromodeling

repo: tapas


text:
function pde = DataDivCurlLinear(varargin)
%%
%           curl u = g    in \Omega,
%           div (eps* u) = f    in \Omega,
%           n \cdot (eps * u) = g_N  on \Gamma_0
% Reference: a linear vector field with modifiable curl and div

if isempty(varargin)
    alpha = 1; % strength of irrotational part (x,y,z)
    beta = 1; % strength of solenoidal part (z,x,y)
else
    alpha = varargin{1};
    beta = varargin{2};
end

pde.Eps = @(p) [1+p(:,1)*0 p(:,3)*0 p(:,3)*0; ...
    p(:,3)*0 1+p(:,2)*0 p(:,3)*0; ...
    0*p(:,1) 0*p(:,2) 1+p(:,3)*0];

pde.exactu = @(p) alpha*[p(:,1), p(:,2), p(:,3)] + ...
                   beta*[p(:,3), p(:,1), p(:,2)];
pde.f = @(p) 3*alpha*ones(size(p,1),1);
pde.g = @(p) beta*[ones(size(p,1),1),...
                  ones(size(p,1),1),...
                  ones(size(p,1),1)];

end

repo: ifem


text:
function [r, R, S] = imnoise3(M, N, C, A, B)
%IMNOISE3 Generates periodic noise.
%    [r, R, S] = IMNOISE3(M, N, C, A, B), generates a spatial
%    sinusoidal noise pattern, r, of size M-by-N, its Fourier
%    transform, R, and spectrum, S.  The remaining parameters are:
%
%    C is a K-by-2 matrix with K pairs of freq. domain coordinates (u,
%    v) that define the locations of impulses in the freq. domain. The
%    locations are with respect to the frequency rectangle center at
%    (floor(M/2) + 1, floor(N/2) + 1).  The impulse locations are spe-
%    cified as increments with respect to the center. For ex, if M =
%    N = 512, then the center is at (257, 257). To specify an impulse
%    at (280, 300) we specify the pair (23, 43); i.e., 257 + 23 = 280,
%    and 257 + 43 = 300. Only one pair of coordinates is required for
%    each impulse.  The conjugate pairs are generated automatically.
%
%    A is a 1-by-K vector that contains the amplitude of each of the
%    K impulse pairs. If A is not included in the argument, the
%    default used is A = ONES(1, K).  B is then automatically set to
%    its default values (see next paragraph).  The value specified
%    for A(j) is associated with the coordinates in C(j, 1:2).
%
%    B is a K-by-2 matrix containing the Bx and By phase components
%    for each impulse pair.  The default values for B are B(1:K, 1:2)
%    = 0.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.5 $  $Date: 2004/11/04 22:32:42 $

% Process input parameters.
[K, n] = size(C);
if nargin == 3
   A(1:K) = 1.0;
   B(1:K, 1:2) = 0;
elseif nargin == 4
   B(1:K, 1:2) = 0;
end

% Generate R.
R = zeros(M, N);
for j = 1:K
   u1 = M/2 + 1 + C(j, 1); v1 = N/2 + 1 + C(j, 2);
   R(u1, v1) = i * (A(j)/2) * exp(i*2*pi*C(j, 1) * B(j, 1)/M);
   % Complex conjugate.
   u2 = M/2 + 1 - C(j, 1); v2 = N/2 + 1 - C(j, 2);
   R(u2, v2) = -i * (A(j)/2) * exp(i*2*pi*C(j, 2) * B(j, 2)/N);
end

% Compute spectrum and spatial sinusoidal pattern.
S = abs(R);
r = real(ifft2(ifftshift(R)));

repo: weiyanmin


text:
function pts2 = grFaceToPolygon(varargin)
%GRFACETOPOLYGON Compute the polygon corresponding to a graph face.
%
%   PTS2 = grFaceToPolygon(NODES, EDGES, FACES, INDF)
%   PTS2 = grFaceToPolygon(NODES, FACES, INDF)
%   Where NODES, EDGES, and FACES are internal data of graph, and INDF is
%   the index of the face to extract. The result is the (ordered) set of
%   points composing the face.
%
%
%   PTS2 = grFaceToPolygon(GRAPH, INDF)
%   use structure representation for graph. The structure GRAPH must
%   contain data for fields 'nodes' and 'faces'.
%
%   If several indices face indices are specified, result is a cell array
%   of polygons.
%
%   The number of columns of PTS2 is the same as for NODES.
%

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2005-11-30
% Copyright 2005-2022 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas)

if length(varargin)==2
    % argument is a graph structure
    graph = varargin{1};
    nodes = graph.nodes;
    faces = graph.faces;
    indf  = varargin{2};

elseif length(varargin)==3
    % arguments are nodes, faces and indices
    nodes = varargin{1};
    faces = varargin{2};
    indf  = varargin{3};

elseif length(varargin)==4
    % arguments are nodes, edges, faces and indices, we forget edges
    nodes = varargin{1};
    faces = varargin{3};
    indf  = varargin{4};
end

if iscell(faces)
    % faces is a cell array
    if length(indf)==1
        face = faces{indf};
        pts2 = nodes(face, :);
    else
        pts2 = cell(length(indf), 1);
        for i=1:length(indf)
            face = faces{indf(i)};
            pts2{i} = nodes(face, :);
        end
    end
else
    % faces is an indices array: all faces have same number of vertices
    if length(indf)==1
        face = faces(indf, :);
        pts2 = nodes(face, :);
    else
        pts2 = cell(length(indf), 1);
        for i=1:length(indf)
            face = faces(indf(i), :);
            pts2{i} = nodes(face, :);
        end
    end
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
% remove all path with a given parent path

% varargin contains a list of path to exclude
function removepath(parentpath, varargin)

if isempty(parentpath), return; end;
folder = path;
if ispc, sep = ';'; else sep = ':'; end;
indSep = find(folder == sep);

indSep = [ 0 indSep length(folder)+1 ];
for iSep = 1:length(indSep)-1
    curPath = folder(indSep(iSep)+1:indSep(iSep+1)-1);
    if ~isempty(strfind(curPath, parentpath)) && ~any(strmatch(curPath, varargin, 'exact'))
        rmpath(curPath);
    end;
end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function sF = sqrt(sF, varargin)
% square root of a function
%
% Syntax
%   sF = sqrt(sF)
%   sF = sqrt(sF, 'bandwidth', bandwidth)
%
% Input
%  sF - @S2FunHarmonic
%
% Output
%  sF - @S2FunHarmonic
%
% Options
%  bandwidth - minimal degree of the spherical harmonic
%

sF = sF.quadrature(@(v) sqrt(sF.eval(v)),varargin{:});

end

author: mtex-toolbox

repo: mtex

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


text:
% The OxyGEN proyect by Profoty.xyz
% Together against the Covid-19
% V4.0 rev.17/03/2020
%
% Website/blog: oxygen.protofy.xyz
% Contact: oxygen@protofy.xyz

%GEOMETRICAL PARAMETERS

clear;

%Length from the bearing to the hinge
l0 = 27.5;

%Length of the vertical wall that supports the hinge
h0 = 6.5;

%Lenght of the vertical support of the bearing
hb = 4.0;

%Bearing radius
br = 1.1;

%Array of different max and min positions of the bearing
h1a = [9, 10, 11];
h2a = [20, 20, 20];

%Number of cycles per cam turn (1 to 3).
nc = 1;

%Minimum radius of the camshaft
rmin = 5.5;

%BREATHING CYCLE PARAMETERS

%Duration of the inhale cycle / duration of the whole cycle
lambda1 = 0.500;
lambda2 = 0.060;

%Soft transition between inhale and exhale cycles
dpsi21 = 0.12;
dpsi12 = 0.15;

%Adjust parameters of the inhale curve
ga1 = pi;
gb1 = 1.2;
ff1 = 20;
sp1 = 1;

%Adjust parametes of the exhale curve
ga2 = 9;
gb2 = .5;
ff2 = 20;

%GENERATION OF THE CURVES AND THE CAMSHAFT

%Generation of the geometry

l = sqrt(l0^2+hb^2);

for i = 1:numel(h1a);

    ymin(i) = h1a(i);
    ymax(i) = h2a(i);

    alphamin(i) = acos((h0 - h1a(i)) / l);
    alphamax(i) = acos((h0 - h2a(i)) / l);

    xmin(i) = l*sin(alphamin(i));
    xmax(i) = l*sin(alphamax(i));

    d(i) = sqrt((xmin(i)-xmax(i))^2 + (ymin(i)-ymax(i))^2);

    alphatan(i) = (alphamin(i) + alphamax(i)) / 2;

    xtan(i) = l*sin(alphatan(i));
    ytan(i) = h0 - l*cos(alphatan(i));

    xsup(i) = xtan(i) + (d(i)/2)*cos(alphatan(i));
    xinf(i) = xtan(i) - (d(i)/2)*cos(alphatan(i));

    ysup(i) = ytan(i) + (d(i)/2)*sin(alphatan(i));
    yinf(i) = ytan(i) - (d(i)/2)*sin(alphatan(i));

    xcam(i) = xtan(i) + (d(i)/2 + rmin + br)*cos(alphatan(i));
    ycam(i) = ytan(i) + (d(i)/2 + rmin + br)*sin(alphatan(i));

end;

%Time increment
dt = 0.01;

%time coordinate during the whole cycle
theta = 0:dt:2*pi;

%Generation of the soft transition between inhale and exhale curves
psi1 = [];
psi2 = [];

for i = 1:numel(theta);

    if theta(i) < (lambda2-dpsi21/2)*2*pi;
        psi1(i) = 0;

    elseif theta(i) < (lambda2+dpsi21/2)*2*pi;
        psi1(i) = 0.5 + 0.5*cos((theta(i)-(lambda2+dpsi21/2)*2*pi)/(2*dpsi21));

    else theta(i) < (lambda1-dpsi12/2)*2*pi;
        psi1(i) = 1;

    end;

    psi2(i) = 1 - psi1(i);

end;

%Generation of the complete breathing cycle rho(theta)

rho = [];
rho1 = [];
rho2 = [];
rho2next = [];

rhomin = 1000;
rhomax = 0;

for i = 1:numel(theta);

    %Inhale curve
    rho1(i) = ff1*normpdf(sp1*theta(i), ga1, gb1);
    rho1next(i) = ff1*normpdf(sp1*theta(i)+sp1*2*pi, ga1, gb1);

    %Exhale curve
    rho2(i) = ff2*gampdf(theta(i), ga2, gb2);
    rho2next(i) = ff2*gampdf(theta(i)+2*pi, ga2, gb2);

    rho(i) = rho1(i);

    %Capturing min and max in order to generate the normalized curve
    if rho(i) > rhomax
        rhomax = rho(i);
    end;

    if rho(i) < rhomin
        rhomin = rho(i);
    end;

end;

%Generation of a normalised curve and camshaft

rhonorm = [];
rhocam = [];

for i = 1:numel(theta);

    rhonorm(i) = (rho(i)-rhomin)/(rhomax-rhomin);

    for n = 1:numel(d)

        rhocam(n,i) = rmin + rhonorm(i)*d(n);

    end;

end;

%Generation of the first derivate of the camshaft geometry to analize and
%validate the design

drho = [];
drhonorm = [];
drhocam = [];

a = rmin;
b = 1/dt;

for i = 1 : numel(rho)-1;

    drho(i) = b*(rho(i+1)-rho(i));
    drhonorm(i) = b*(rhonorm(i+1)-rhonorm(i));
    drhocam(i) = b*(rhocam(i+1)-rhocam(i)) + a;

end;

drho(numel(rho)) = b*(rho(1)-rho(numel(rho)));
drhonorm(numel(rhonorm)) = b*(rhonorm(1)-rhonorm(numel(rhonorm)));
drhocam(numel(rhocam)) = b*(rhocam(1)-rhocam(numel(rhocam))) + a;

%X and Y coordinates during two cycles (for 2-cycle camshaft plot)

theta2 = [theta/2, (2*pi+theta)/2];

rhocam2 = [rhocam, rhocam];
drhocam2 = [drho, drho];

%X and Y coordinates during three cycles (for 3-cycle camshaft plot)

theta3 = [theta/3, (2*pi+theta)/3, (4*pi+theta)/3];

rhocam3 = [rhocam, rhocam, rhocam];
drhocam3 = [drho, drho, drho];

%GENERATION OF THE PLOTS

%Trying to print it out in real size (FAIL)
%set(gcf,'PaperUnits','centimeters');
%set(gcf,'PaperSize',[42 29.7]);

fpos = figure('Name', 'Dimensions', 'Units', 'centimeters', 'NumberTitle', 'off');
set(gcf,'PaperUnits','centimeters', 'PaperSize', [42/2 27.9/2]);
set(gcf, 'units', 'centimeters', 'position', [0, 0, 28, 20]);
fpos = gcf;
hold on

xlim([-5 35]);
ylim([0 30]);

for i = 1:numel(xmin)

    plot([0, xmin(i)], [h0, ymin(i)], '-x')
    plot([0, xmax(i)], [h0, ymax(i)], '-x')

    plot([0, xtan(i)], [h0, ytan(i)], '--xb')

    plot([xsup(i), xinf(i)], [ysup(i), yinf(i)], '--xr')
    scatter(xcam(i), ycam(i))

end;

hold off

%Plot of the breathing cycle interpolation by curves

fcycle = figure('Name', 'Breath Cycle curves', 'Units', 'centimeters', 'NumberTitle', 'off');

hold on
plot(theta, psi1, '-k')
plot(theta, psi2, '-k')
plot(theta, rho1, '-g')
plot(theta, rho1next, '-g')
plot(theta, rho2, '-g')
plot(theta, rho2next, '-g')
%plot(theta, drho, '-r')
plot(theta, rho, '-b')
hold off

set(gcf,'PaperUnits','centimeters', 'PaperSize', [42 27.9]);
set(gcf, 'units', 'centimeters', 'position', [0, 0, 27.9, 27.9]);
fcam3 = gcf;

%Plot of the normalized breathing cycle and first derivate

fnorm = figure('Name', 'Normalized Breath Cycle', 'Units', 'centimeters', 'NumberTitle', 'off');

hold on
%plot(theta, drhonorm, '-r')
plot(theta, rhonorm, '-b')
hold off

set(gcf,'PaperUnits','centimeters', 'PaperSize', [42/2 27.9/2]);
set(gcf, 'units', 'centimeters', 'position', [0, 0, 27.9, 27.9]);
fcam3 = gcf;

%Plot of a 1-cycle camshaft

for n = 1:numel(d)

    fcam1 = figure('Name', '1-Cycle Camshaft', 'Units', 'centimeters', 'NumberTitle', 'off');
    %hold on

    %polarplot(theta, drhocam, '-r')
    polar(theta, rhocam(n,:), '-b')

    %rlim([0 25]);
    hold off

    set(gcf,'PaperUnits','centimeters', 'PaperSize', [42 27.9]);
    set(gcf, 'units', 'centimeters', 'position', [0, 0, 27.9, 27.9]);
    fcam3 = gcf;

end

% %Plot of a 2-cycle camshaft
%
% fcam2 = figure('Name', '2-Cycle', 'Units', 'centimeters', 'NumberTitle', 'off');
%
% %hold on;
% %polar(theta, drhocam2, '-r')
% polar(theta2, rhocam2, '-b')
% %hold off;
%
% set(gcf,'PaperUnits','centimeters', 'PaperSize', [42 27.9]);
% set(gcf, 'units', 'centimeters', 'position', [0, 0, 27.9, 27.9]);
% fcam3 = gcf;
%
%
% %Plot of a 3-cycle camshaft
%
% fcam3 = figure('Name', '3-Cycle', 'Units', 'centimeters', 'NumberTitle', 'off');
%
% %hold on;
% %polar(theta,drhocam3)
% polar(theta3, rhocam3)
% %hold off;
%
% set(gcf,'PaperUnits','centimeters', 'PaperSize', [42 27.9]);
% set(gcf, 'units', 'centimeters', 'position', [0, 0, 27.9, 27.9]);
% fcam3 = gcf;

% The OxyGEN proyect by Profoty.xyz
% Together against the Covid-19
% V4.0 rev.17/03/2020
%
% Website/blog: oxygen.protofy.xyz
% Contact: oxygen@protofy.xyz

author: ProtofyTeam

repo: OxyGEN

save_path: github-repos/MATLAB/ProtofyTeam-OxyGEN


text:
function [status, results] = mrtrix_csdeconv(dwi_file, response_file, lmax, ...
                                   out_file, grad_file, mask_file, verbose)

%function [status, results] = mrtrix_csdeconv(dwi_file, response_file, lmax, ...
%           out_file, grad_file, [mask_file=entire volume], [verbose=true])
%
% Fit the constrained spherical deconvolution model to dwi data
%
% Parameters
% ----------
% dwi_file: The name of a dwi file in .mif format.
% response_file: The name of a .txt fiber response function file.
% lmax: The maximal harmonic order.
% out_file: The resulting .mif file.
% grad_file: a file containing gradient information in the mrtrix format.
% mask_file: a .mif file containing a mask. Default: the entire volume.
% verbose: Whether to display stdout to the command window.
%
% Returns
% -------
% status: whether (0) or not (1) the operation succeeded
% results: the results of the operation in stdout
%
% Notes
% -----
% http://www.brain.org.au/software/mrtrix/tractography/preprocess.html
%

if notDefined('verbose')
    verbose = true;
end

if notDefined('bkgrnd')
    bkgrnd = false;
end

if notDefined('mask_file')
cmd_str = sprintf('csdeconv %s %s -lmax %d %s -grad %s',...
    dwi_file, response_file, lmax, out_file, grad_file);

else
cmd_str = sprintf('csdeconv %s %s -lmax %d -mask %s %s -grad %s',...
    dwi_file, response_file, lmax, mask_file, out_file, grad_file);
end

[status,results] = mrtrix_cmd(cmd_str, bkgrnd, verbose);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
  function [out ii] = jf_pair_parse(cells, str, varargin)
%|function [out ii] = jf_pair_parse(cells, str, varargin)

arg.suffix = ' ';
arg = vararg_pair(arg, varargin);

str = [str arg.suffix];

tf = strncmp(str, cells, numel(str));
ii = find(tf);
if numel(ii) < 1
	fail('no match for "%s"', str)
end
if numel(ii) > 1
	fail('%d matches for "%s"', numel(ii), str)
end

out = cells{ii};
out = out((numel(str)+1):end);

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
classdef opf_model < opt_model
%OPF_MODEL  Constructor for OPF model class.
%   OM = OPF_MODEL(MPC)
%
%   This class implements the OPF model object used to encapsulate
%   a given OPF problem formulation. It allows for access to optimization
%   variables, constraints and costs in named blocks, keeping track of the
%   ordering and indexing of the blocks as variables, constraints and costs
%   are added to the problem.
%
%   This class is a subclass of OPT_MODEL that adds the 'mpc'
%   field for storing the MATPOWER case struct used to build the object
%   along with the get_mpc() method.
%
%   It also adds the 'cost' field and the following three methods for
%   implementing the legacy user-defined OPF costs:
%       add_legacy_cost
%       params_legacy_cost
%       eval_legacy_cost
%
%   The following is the structure of the data in the OPF model object.
%
%   om
%       <opt_model fields> - see OPT_MODEL for details
%       .cost       - data for legacy user-defined costs
%           .idx
%               .i1 - starting row index within full N matrix
%               .iN - ending row index within full N matrix
%               .N  - number of rows in this cost block in full N matrix
%           .N      - total number of rows in full N matrix
%           .NS     - number of cost blocks
%           .data   - data for each user-defined cost block
%               .N  - see help for ADD_LEGACY_COST for details
%               .H  -               "
%               .Cw -               "
%               .dd -               "
%               .rr -               "
%               .kk -               "
%               .mm -               "
%               .vs - cell array of variable sets that define xx for this
%                     cost block, where the N for this block multiplies xx
%           .order  - struct array of names/indices for cost blocks in the
%                     order they appear in the rows of the full N matrix
%               .name   - name of the block, e.g. R
%               .idx    - indices for name, {2,3} => R(2,3)
%       .mpc        - MATPOWER case struct used to create this model object
%           .baseMVA
%           .bus
%           .branch
%           .gen
%           .gencost
%           .A  (if present, must have l, u)
%           .l
%           .u
%           .N  (if present, must have fparm, H, Cw)
%           .fparm
%           .H
%           .Cw
%
%   See also OPT_MODEL.

%   MATPOWER
%   Copyright (c) 2008-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

    properties
        cost = [];
        mpc = struct();
    end     %% properties

    methods
        %% constructor
        function om = opf_model(mpc)
            args = {};
            have_mpc = 0;
            if nargin > 0
                if isa(mpc, 'opf_model')
                    args = { mpc };
                elseif isstruct(mpc)
                    have_mpc = 1;
                end
            end

            %% call parent constructor
            om@opt_model(args{:});

            if have_mpc
                om.mpc = mpc;
            end

            %% Due to a bug related to inheritance in constructors in
            %% Octave 5.2 and earlier (https://savannah.gnu.org/bugs/?52614),
            %% INIT_SET_TYPES() cannot be called directly in the
            %% MP_IDX_MANAGER constructor, as desired.
            %%
            %% WORKAROUND:  INIT_SET_TYPES() is called explicitly as needed
            %%              (if om.var is empty) in ADD_VAR(), DISPLAY() and
            %%              INIT_INDEXED_NAME(), after object construction,
            %%              but before object use.
        end

        function om = def_set_types(om)
            om.set_types = struct(...
                    'var',  'VARIABLES', ...
                    'lin',  'LINEAR CONSTRAINTS', ...
                    'nle',  'NONLIN EQ CONSTRAINTS', ...
                    'nli',  'NONLIN INEQ CONSTRAINTS', ...
                    'qdc',  'QUADRATIC COSTS', ...
                    'nlc',  'GEN NONLIN COSTS', ...
                    'cost', 'LEGACY COSTS'  ...
                );
        end

        function om = init_set_types(om)
            %% call parent to create base data structures for each type
            init_set_types@opt_model(om);

            %% finish initializing data structures for each type
            es = struct();  %% empty struct
            om.cost.data = struct( ...
                'N', es, ...
                'H', es, ...
                'Cw', es, ...
                'dd', es, ...
                'rh', es, ...
                'kk', es, ...
                'mm', es, ...
                'vs', es );
            om.cost.params = [];
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
addpath('../src/')
addpath('../src/utils/')

% fix random seed
rng(0);

% make your own discrete linear system with disturbance
A = [1 1; 0 1];
B = [0.5; 1];
Q = diag([1, 1]);
R = 0.1;

W_vertex = [0.15, 0.15; 0.15, -0.15; -0.15, -0.15; -0.15, 0.15]; % construct a convex set of disturbance (2dim here)
W = Polyhedron(W_vertex);

% construct disturbance Linear system
disturbance_system = DisturbanceLinearSystem(A, B, Q, R, W);

% constraints on state Xc and input Uc
Xc_vertex = [2, -2; 2 2; -10 2; -10 -2];
Uc_vertex = [1; -1];
Xc = Polyhedron(Xc_vertex);
Uc = Polyhedron(Uc_vertex);

% create a tube_mpc simulater
% if N_horizon is too small, the path will never reach inside the robust MPI-set X_mpi_robust in time step N_horizon, then the problem becomes infeasible.
N_horizon = 10;
w_min = [0; -0.10];
w_max = [0; 0.10];
mpc = TubeModelPredictiveControl(disturbance_system, Xc, Uc, N_horizon);

% The robust MPC guidances the path inside the robust MPI-set so that the path will reach the robust MPI-set in N_horizon.
x = [-7; -2];
savedir_name = './results/';
mkdir(savedir_name);

for i = 1:15
    disp(i)
    u_next = mpc.solve(x);
    x = disturbance_system.propagate(x, u_next); % additive disturbance is considered inside the method
    mpc.show_prediction();
    filename = strcat(savedir_name, 'tmpc_seq', number2string(i), '.png')
    saveas(gcf, char(filename)); % removing this line makes the code much faster
    clf;
end

author: HiroIshida

repo: robust-tube-mpc


text:
function [err, predictedDeg] = poirsonWandellCortMagErr(parm, cortMag)
%
% [err, predictedDeg] = poirsonWandellCortMagErr(parm, cortMag)
%
% AUTHOR: Wandell
% DATE:  11.17.00
% PURPOSE:
%    Fit both the distance scale factor and the foveal phase
% to a set of expanding ring data
%
%

dScale = parm(1);
dShift = parm(2);

% We restrict the range on the estimated fovealPhase to be
% a little past CORTMAG.fovealPhase.  This variable has the
% foveal phase of the stimulus.
% If the search parameter is smaller than the stimulus phase, well, that
% can't be right.  We belive in causality.
% If the parameter is more than 4-5 sec of hemodynamic delay,
% we also get unhappy and kick back a large error.
% We don't ask for the time, sigh.  So for now
% we assume the period is 36 sec and we demand that the phase be less than
% 2*pi/9;
% Finally, we need to make these measurements in complex phase representation
% to avoid wrapping problems near 0 and 2pi
%
phaseDifference = angle(exp(sqrt(-1)*(parm(3) - cortMag.fovealPhase)));
radPerSec = (2*pi)/36;
oneSec    = 1*radPerSec;
fiveSec   = 1.1*radPerSec;
if (phaseDifference < oneSec) | (phaseDifference > fiveSec)
    err = 10000000;
    if nargout == 2
        predictedDeg = NaN;
    end
    return;
else
    fovealPhase = parm(3);
end

% Use the basic function of distance to predict the visual
% field representation in degrees given the distance scale
% parameter, dScale.  This form keeps dist = 0 at 10 deg.
dist = cortMag.allCorticalDist - dShift;
predictedDeg = exp(dScale*dist + log(10));
% figure(testFig);
% plot(dist, predictedDeg,'o')

% Convert the predicted degrees into expected stimulus phase
% as a complex number.  This conversion assumes that the foveal
% phase is 0.
predictedRad = (2*pi)*(predictedDeg/cortMag.stimulusRadius);
predictedCx = exp(sqrt(-1)*predictedRad);
% figure(testFig);
% plot(dist, predictedRad,'o')
% plot(dist, angle(predictedCx),'o')

% Adjust the observed CX phases so that the foveal phase is 0.
% This is assumed by the function when we map from phase
% to degrees, later.
% WE SHOULD BE ABLE TO SET A RANGE OR EVEN FIX THE fovealPhase
% Figure out how to do this here!
%
observedCx = cortMag.allMeanPh/exp(sqrt(-1)*fovealPhase);

err = norm(predictedCx - observedCx);

if nargout == 2
   predictedRad = angle(predictedCx);
   l = find(predictedRad < 0);
   predictedRad(l) = predictedRad(l) + 2*pi;

   % Take into account the fact that we sometimes run the
   % phase map with a restricted range
   %
   phRange = phaseRange(CORTMAG);
   rad2deg = cortMag.stimulusRadius/phRange;
   predictedDeg = rad2deg*predictedRad;
   % figure(testFig);
   % plot(dist, predictedDeg,'o')
end

return;

% DEBUGGING AND NOTES


text:
%% 2's complement of binary no;

function twos_comp = twos_complement_of_binary(bin)
c=0;
temp=bin;
twos_comp=0;
while(temp>0 && rem(temp,10)==0)
    twos_comp=twos_comp*10 + rem(temp,10);
    temp=fix(temp/10);
    c=c+1;
end
if(temp>0)
    twos_comp=twos_comp*10 + rem(temp,10);
    temp=fix(temp/10);
    c=c+1;
end
while(temp>0)
    if(rem(temp,10)==1)
        twos_comp=twos_comp*10 + 0;
    else
        twos_comp=twos_comp*10 + 1;
    end
    temp=fix(temp/10);
    c=c+1;
end

temp=twos_comp;
twos_comp=0;
while(c>0)      %reversing the order;
    twos_comp=twos_comp*10 + rem(temp,10);
    temp=fix(temp/10);
    c=c-1;
end
end

author: TheAlgorithms

repo: MATLAB-Octave


text:
function for_euler
global sxyz sXYZ sa sN sav saa

set(sxyz,'Value',true);
set(sXYZ,'Value',true);
set(sa,'Value',true);
set(sN,'Value',true);
set(sav,'Value',false);
set(saa,'Value',false);

XYZ_visibility;
enter_eul_angs;

author: Sable

repo: mcbench-benchmarks


text:
function [vertices, landmarks_fitted, fit_params] = read_face_recon(reconstruction_dir)
    vertices = import_vertices_keypoints([reconstruction_dir, '/vertices.txt']);
    landmarks_fitted = import_vertices_keypoints([reconstruction_dir, '/landmarks_fitted.txt']);
    landmarks_fitted = [landmarks_fitted(59:end, :); landmarks_fitted(1, :); ...
        landmarks_fitted(51:58, :); landmarks_fitted(2:50, :)]; % reorder
    fit_params = import_params([reconstruction_dir, '/fit_params.txt']);
end

function vertices = import_vertices_keypoints(filename)
dataLines = [1, Inf];
opts = delimitedTextImportOptions("NumVariables", 3);
opts.DataLines = dataLines;
opts.Delimiter = ",";
opts.VariableNames = ["VarName1", "VarName2", "VarName3"];
opts.VariableTypes = ["double", "double", "double"];
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";
vertices = readtable(filename, opts);
vertices = table2array(vertices);
end

function params = import_params(filename, dataLines)
if nargin < 2
    dataLines = [1, Inf];
end
opts = delimitedTextImportOptions("NumVariables", 1);
opts.DataLines = dataLines;
opts.Delimiter = ",";
opts.VariableNames = "VarName1";
opts.VariableTypes = "double";
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";
params = readtable(filename, opts);
params = table2array(params);
end

author: zhixuany

repo: HUMBI

save_path: github-repos/MATLAB/zhixuany-HUMBI


text:
%FILENAMES Get filenames of datafile
%
%		[NAMES,DIRS,ROOT] = FILENAMES(A)
%
% INPUT
%   A         DATAFILE
%
% OUTPUT
%   NAMES     Names of the files in which the objects of A are stored.
%   DIRS      Names of the directories in which the objects of A are
%             stored.
%   ROOT      Rootdirectory of the datafile
%
% DESCRIPTION
% This routine facilitates the retrieval of the files where the objects of
% A are stored
%
% SEE ALSO (PRTools Guide)
% DATAFILES

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function a = shake();
global layout;
rand('twister',sum(100*clock));

a=1;
figure(1);
plot([0 1],[.25 .25],'k',[0 1],[.5 .5],'k',[0 1],[.75 .75],'k',[0 1],[1 1],'k'); hold on;
plot([.25 .25],[0 1],'k',[.5 .5],[0 1],'k',[.75 .75],[0 1],'k',[1 1],[0 1],'k');
axis([0 1 0 1],'equal','square');

Cube(1).letter = {'L','R','Y','T','T','E'};
Cube(2).letter = {'V','T','H','R','W','E'};
Cube(3).letter = {'E','G','H','W','N','E'};
Cube(4).letter = {'S','E','O','T','I','S'};
Cube(5).letter = {'A','N','A','E','E','G'};
Cube(6).letter = {'I','D','S','Y','T','T'};
Cube(7).letter = {'O','A','T','T','O','W'};
Cube(8).letter = {'M','T','O','I','C','U'};
Cube(9).letter = {'A','F','P','K','F','S'};
Cube(10).letter = {'X','L','D','E','R','I'};
Cube(11).letter = {'H','C','P','O','A','S'};
Cube(12).letter = {'E','N','S','I','E','U'};
Cube(13).letter = {'Y','L','D','E','V','R'};
Cube(14).letter = {'Z','N','R','N','H','L'};
Cube(15).letter = {'N','M','I','Qu','H','U'};
Cube(16).letter = {'O','B','B','A','O','J'};

rand(1,16);
[temp,order]=sort(rand(1,16));
Cube=Cube(order);

for i=1:4,
	for j=1:4,
		layout{i,j} = Cube(i+(j-1)*4).letter(ceil(rand*6));
		text(i/4-.15,j/4-.13,layout{i,j},'FontSize',20);
	end;
end;
set(gca,'XTick',[],'YTick',[]);


text:
% ------------------------------------------------------------------------
function stats = vid_image_stats(imdb_video, perc_training, base_path)
%VID_IMAGE_STATS
%	Compute basic colour stats for a random PERC_TRAINING of the dataset
%	Used for data augmentation during training.
% 	e.g. vid_image_stats(imdb_video, 0.1, /path/to/curated/ILSVRC15/')
% ------------------------------------------------------------------------
    % collect different stats for z and x crops (x contains more padding)
    z_sz = 127;
    x_sz = 255;
    samples_per_video = 16;
    imout_z = zeros(z_sz, z_sz, 3, samples_per_video, 'single');
    imout_x = zeros(x_sz, x_sz, 3, samples_per_video, 'single');
    crops_z_string = cell(1, samples_per_video);
    crops_x_string = cell(1, samples_per_video);
    n_video = numel(imdb_video.id);
    n_video_train = round(perc_training * n_video);
    avg_z = cell(1,n_video_train);
    rgbm1_z = cell(1,n_video_train);
    rgbm2_z = cell(1,n_video_train);
    n_z = samples_per_video*z_sz*z_sz;
    avg_x = cell(1,n_video_train);
    rgbm1_x = cell(1,n_video_train);
    rgbm2_x = cell(1,n_video_train);
    n_x = samples_per_video*x_sz*x_sz;
    for v=1:n_video_train
        n_obj = numel(imdb_video.objects{v});
        rand_objs = datasample(1:n_obj, samples_per_video);
        for o=1:samples_per_video
            crops_z_string{o} = [base_path strrep(imdb_video.objects{v}{rand_objs(o)}.frame_path, '.JPEG','') '.' num2str(imdb_video.objects{v}{rand_objs(o)}.track_id, '%02d') '.crop.z.jpg'];
            crops_x_string{o} = [base_path strrep(imdb_video.objects{v}{rand_objs(o)}.frame_path, '.JPEG','') '.' num2str(imdb_video.objects{v}{rand_objs(o)}.track_id, '%02d') '.crop.x.jpg'];
        end
        files = [crops_z_string crops_x_string];
        imgs = vl_imreadjpeg(files, 'numThreads', 12);
        crops_z = imgs(1:samples_per_video);
        crops_x = imgs(samples_per_video+1 : end);
        for o=1:samples_per_video
           imout_z(:,:,:,o) = crops_z{o};
           imout_x(:,:,:,o) = crops_x{o};
        end

        Z = reshape(permute(imout_z,[3 1 2 4]),3,[]);
        avg_z{end+1} = mean(imout_z, 4);
        rgbm1_z{end+1} = sum(Z,2)/n_z;
        rgbm2_z{end+1} = Z*Z'/n_z;

        X = reshape(permute(imout_x,[3 1 2 4]),3,[]);
        avg_x{end+1} = mean(imout_x, 4);
        rgbm1_x{end+1} = sum(X,2)/n_x;
        rgbm2_x{end+1} = X*X'/n_x;
        fprintf('Processed video %d/%d\n', v, n_video_train);
    end

    stats = struct();
    stats.z = struct();
    stats.x = struct();

    stats.z.averageImage = mean(cat(4,avg_z{:}),4);
    rgbm1_z = mean(cat(2,rgbm1_z{:}),2);
    stats.z.rgbm1 = rgbm1_z;
    rgbm2_z = mean(cat(3,rgbm2_z{:}),3);
    stats.z.rgbm2 = rgbm2_z;
    stats.z.rgbMean = rgbm1_z;
    stats.z.rgbCovariance = rgbm2_z - rgbm1_z*rgbm1_z';

    stats.x.averageImage = mean(cat(4,avg_x{:}),4);
    rgbm1_x = mean(cat(2,rgbm1_x{:}),2);
    stats.x.rgbm1 = rgbm1_x;
    rgbm2_x = mean(cat(3,rgbm2_x{:}),3);
    stats.z.rgbm2 = rgbm2_z;
    stats.x.rgbMean = rgbm1_x;
    stats.x.rgbCovariance = rgbm2_x - rgbm1_x*rgbm1_x';
end

author: bertinetto

repo: siamese-fc


text:
% Mediterranian diet defintion. For details, please see https://www.vmh.life/#nutrition
% Ines Thiele 2016-2019

Diet = {%'Reaction'	'Flux Value'
'EX_etoh(e)'	'355.9923939'
'EX_h2o(e)'	'220593.0188'
'EX_caro(e)'	'0.005979072'
'EX_retinol(e)'	'0.558558584'
'EX_thm(e)'	'8.328478195'
'EX_adpcbl(e)'	'1.22881E-05'
'EX_ribflv(e)'	'0.006111107'
'EX_nac(e)'	'0.445244325'
'EX_ncam(e)'	'0.211750949'
'EX_pnto_R(e)'	'0.085186526'
'EX_pydam(e)'	'0.004669003'
'EX_pydxn(e)'	'0.004669643'
'EX_pydx(e)'	'0.004725954'
'EX_btn(e)'	'0.00036173'
'EX_10fthf(e)'	'0.000264396'
'EX_5mthf(e)'	'0.000271879'
'EX_thf(e)'	'0.000281098'
'EX_ascb_L(e)'	'0.782052146'
'EX_vitd3(e)'	'3.18352E-05'
'EX_avite1(e)'	'0.054573176'
'EX_phyQ(e)'	'0.000301578'
'EX_ca2(e)'	'44.23174809'
'EX_cl(e)'	'543.2352692'
'EX_k(e)'	'101.3033815'
'EX_mg2(e)'	'25.58362477'
'EX_na1(e)'	'150.9506232'
'EX_pi(e)'	'23.46662249'
'EX_cu2(e)'	'0.04957039'
'EX_fe2(e)'	'0.200734175'
'EX_fe3(e)'	'0.200734175'
'EX_mn2(e)'	'0.136335394'
'EX_zn2(e)'	'0.261394922'
'EX_mnl(e)'	'0.779483982'
'EX_xylt(e)'	'0'
'EX_lcts(e)'	'35.94427638'
'EX_malt(e)'	'4.206879124'
'EX_sucr(e)'	'110.5182203'
'EX_fru(e)'	'125.4517893'
'EX_gal(e)'	'10.8426091'
'EX_cellul(e)'	'0.14529578'
'EX_ala_L(e)'	'61.65072082'
'EX_arg_L(e)'	'42.32821506'
'EX_asp_L(e)'	'80.87258611'
'EX_cys_L(e)'	'13.29163028'
'EX_glu_L(e)'	'168.7925032'
'EX_gly(e)'	'67.07363861'
'EX_urate(e)'	'5.780252608'
'EX_his_L(e)'	'18.55149638'
'EX_ile_L(e)'	'45.97634115'
'EX_leu_L(e)'	'73.5657289'
'EX_lys_L(e)'	'49.72502556'
'EX_met_L(e)'	'15.16359686'
'EX_phe_L(e)'	'34.89104305'
'EX_pro_L(e)'	'78.89334277'
'EX_ser_L(e)'	'59.78051737'
'EX_thr_L(e)'	'39.41362937'
'EX_trp_L(e)'	'6.728111908'
'EX_tyr_L(e)'	'23.68527804'
'EX_val_L(e)'	'54.75068355'
'EX_dca(e)'	'3.170503694'
'EX_ddca(e)'	'3.503490546'
'EX_ttdca(e)'	'12.65654482'
'EX_ttdcea(e)'	'1.377519391'
'EX_ptdca(e)'	'1.183522405'
'EX_hdca(e)'	'56.41911375'
'EX_hpdca(e)'	'0.886830484'
'EX_ocdca(e)'	'13.64408774'
'EX_ocdcea(e)'	'128.1607541'
'EX_lnlc(e)'	'73.59427421'
'EX_lnlnca(e)'	'7.87609621'
'EX_strdnc(e)'	'0.614547697'
'EX_arach(e)'	'1.876943078'
'EX_CE2510(e)'	'5.347026525'
'EX_CE4843(e)'	'0.130800704'
'EX_arachd(e)'	'0.268570142'
'EX_docosac(e)'	'2.557574808'
'EX_doco13ac(e)'	'4.673543672'
'EX_adrn(e)'	'0.025338442'
'EX_clpnd(e)'	'0.108954216'
'EX_crvnc(e)'	'3.860995608'
'EX_but(e)'	'8.165474895'
'EX_octa(e)'	'1.800374991'
'EX_chsterol(e)'	'0.687954461'
'EX_sbt_D(e)'	'8.441207695'
'EX_glc_D(e)'	'96.02127935'
'EX_hdcea(e)'	'5.649717068'
'EX_lgnc(e)'	'0.909993994'
'EX_fol(e)'	'0.000283028'
'EX_strch1(e)'	'41.02641436'
'EX_i(e)'	'0.003914204'
'EX_starch1200(e)'	'0.379838969'

    };
Diet = regexprep(Diet,'EX_','Diet_EX_');
Diet = regexprep(Diet,'\(e\)','\[d\]');

repo: cobratoolbox


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

%%
clear; close all; clc;

%% Syntax
% |[N]=numConnect(F,V);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function nav=uniqgeph(nav)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Copyright(c) 2020-2025, by Kai Chen, All rights reserved.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global gls;
geph0=nav.geph;
A=zeros(nav.ng,3);
geph1=repmat(gls.geph,nav.ng,1);

%sort geph
for i=1:nav.ng
    A(i,:)=[geph0(i).tof.time,geph0(i).toe.time,geph0(i).sat];
end
[~,index]=sortrows(A);
for i=1:nav.ng
    geph1(i)=geph0(index(i));
end
clear geph0;

%unique geph
j=1;
for i=2:nav.ng
    if geph1(i).sat~=geph1(j).sat||geph1(i).toe.time~=geph1(j).toe.time||geph1(i).svh~=geph1(j).svh
        j=j+1;
        geph1(j)=geph1(i);
    end
end

nav.ng=j;
geph2=repmat(gls.geph,nav.ng,1);
for i=1:nav.ng
    geph2(i)=geph1(i);
end
clear geph1;

nav=rmfield(nav,'geph');
nav.geph=geph2;

return

repo: GINav