text:
function  [coreness,kn] = kcoreness_centrality_bu(CIJ)
%KCORENESS_CENTRALITY_BU       K-coreness centrality
%
%   [coreness,kn] = kcoreness_centrality_bu(CIJ)
%
%   The k-core is the largest subgraph comprising nodes of degree at least
%   k. The coreness of a node is k if the node belongs to the k-core but
%   not to the (k+1)-core. This function computes the coreness of all nodes
%   for a given binary undirected connection matrix.
%
%   input:          CIJ,        connection/adjacency matrix (binary, undirected)
%
%   output:    coreness,        node coreness.
%                    kn,        size of k-core
%
%   References: e.g. Hagmann et al. (2008) PLoS Biology
%
%   Olaf Sporns, Indiana University, 2007/2008/2010/2012

N = size(CIJ,1);

% determine if the network is undirected - if not, compute coreness on the
% corresponding undirected network
CIJund = CIJ+CIJ';
if (any(CIJund(:)>1))
    CIJ = double(CIJund>0);
end;

coreness = zeros(1,N); kn = zeros(1,N);
for k=1:N
    [CIJkcore,kn(k)] = kcore_bu(CIJ,k);
    ss = sum(CIJkcore)>0;
    coreness(ss) = k;
end;

author: fieldtrip

repo: fieldtrip


text:
%
%
%
%      <>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> %
%      <                                                      > %
%      <                    DISCLAIMER:                       > %
%      <                                                      > %
%      <  THIS PROGRAM IS INTENDED FOR RESEARCH PURPOSES ONLY. > %
%      <  THIS PROGRAM IS IN NO WAY INTENDED FOR CLINICAL OR   > %
%      <                    OFFICIAL USE.                     > %
%      <                                                      > %
%      <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<> %
%
%
%

function [read,setup,sensorNames,sensorLocations,sensorOrientations] = readdata(folder,CHAN,TIME,TRIALS);

[setup,sensorIndex,sensorInfo] = readresfile(folder);

if nargin == 1
    CHAN = 'allchans';
    TIME = 'alltimes';
    TRIALS = 'alltrials';
end

[read,sensorNames,sensorLocations,sensorOrientations] = readmegfile(folder,setup,sensorIndex,sensorInfo,CHAN,TIME,TRIALS);

% Use function [read,setup,sensorNames,sensorLocations] = readdata(folder); to read in parts
% or all of data in a dataset. readdata.m will access readresfile.m to read
% in header, gain/offset, and channel information. It will access
% readmegfile.m to read in the data.
% EXAMPLE:
% [read,sensorNames] = readdata('/data/directory/datasetname.ds');
%
% Enter channel/time/trial parameters in on line 5 of readdata.m.
% FOR THE FUNCTION [read] = readmegfile(folder,setup,sensor,gain,'megsens','alltimes','alltrials');
%
% INPUTS---------------------------------------------------------------------
% folder:     In readdata.m, this should remain as 'folder' everywhere in the
%             program. It is in the command line of matlab where you indicate
%             the file and directory of the dataset you wish to view. It will
%             automatically be used in both readresfile.m and readmegfile.m.
%
% setup:      Enter as 'setup'; setup header information acquired by readresfile.m.
%             This input is needed for the program to work.
%
% sensorIndex:     Enter as 'sensorIndex'; sensor information acquired by readresfile.m,
%             such as megsens, othersens, refsens, and eegsens.
%
% sensorInfo:       Emter as 'gain'; gain/offset information acquired by readresfile.m,
%             including proper_gain, io_gain, and channel locations.
%
% CHAN:       ex: [30:35] - an interval of the desired channels to be read.
%             If CHAN = 'eegsens', only eeg channels/sensors will be read.
%             If CHAN = 'megsens', only meg channels/sensors will be read.
%             If CHAN = 'refsens', only reference channels/sensors will be read.
%             If CHAN = 'othersens' only the other channels/sensors will be
%             read.
%
% TIME:       ex. [0 5] - seconds 0 to 5: the desired time interval to be read.
%             If TIME = 'alltimes', the entire duration of the trial(s) will
%             be read (i.e. TIME = [1:setup.duration]).
%
% TRIALS:     If TRIALS = n, the nth trial will be read.
%             If TRIALS = [3,5,8] (for example), trials 3,5, and 8 will be
%             read and read{1} = data for trial three, read{2} = data for
%             trial 5, and read{3} = data for trial 8.
%             If TRIALS = [3:7] (for example), trials 3 through 7 will be
%             read.
%             If TRIALS = 'alltrials', the data for all of the trials will be
%             read (i.e. TRIALS = [1:setup.duration]).
%
% OUTPUTS--------------------------------------------------------------------
% read:   read contains all of the data. ex. type read{1} for it to display the first
%         set of data on the screen.
%
%----------------------------------------------------------------------------
%
% EXAMPLE:
% [read] = readmegfile(folder,setup,sensorIndex,sensorInfo,'megsens',[0 4],'alltrials');

author: PatternRecognition

repo: OpenBMI


text:
if ~exist('n_ima')|~exist('fc'),
   fprintf(1,'No calibration data available.\n');
   return;
end;

if ~exist('no_image'),
    no_image = 0;
end;

if ~exist('est_alpha'),
    est_alpha = 0;
end;

if ~exist('center_optim'),
    center_optim = 1;
end;

if ~exist('est_aspect_ratio'),
    est_aspect_ratio = 1;
end;

if ~exist('est_fc');
    est_fc = [1;1]; % Set to zero if you do not want to estimate the focal length (it may be useful! believe it or not!)
end;

if exist('est_dist'),
   if length(est_dist) == 4,
      est_dist = [est_dist ; 0];
   end;
end;

if exist('kc'),
   if length(kc) == 4;
      kc = [kc;0];
   end;
end;

if ~exist('fc_error'),
   fc_error = NaN*ones(2,1);
end;

if ~exist('kc_error'),
   kc_error = NaN*ones(5,1);
end;

if ~exist('cc_error'),
   cc_error = NaN*ones(2,1);
end;

if ~exist('alpha_c_error'),
   alpha_c_error = NaN;
end;

check_active_images;

if ~exist('solution_init'), solution_init = []; end;

for kk = 1:n_ima,
   if ~exist(['dX_' num2str(kk)]), eval(['dX_' num2str(kk) '= dX;']); end;
   if ~exist(['dY_' num2str(kk)]), eval(['dY_' num2str(kk) '= dY;']); end;
end;

if ~exist('solution'),
    solution = [];
end;

if ~exist('param_list'),
   param_list = solution;
end;

if ~exist('wintx'),
   wintx = [];
   winty = [];
end;

if ~exist('dX_default'),
   dX_default = [];
   dY_default = [];
end;

if ~exist('dX'),
    dX = [];
end;
if ~exist('dY'),
    dY = dX;
end;

if ~exist('Hcal'),
    Hcal = [];
end;
if ~exist('Wcal'),
    Wcal = [];
end;
if ~exist('map'),
    map = [];
end;

if ~exist('wintx_default')|~exist('winty_default'),
    wintx_default = max(round(nx/128),round(ny/96));
    winty_default = wintx_default;
end;

if ~exist('alpha_c'),
   alpha_c = 0;
end;

if ~exist('err_std'),
    err_std = [NaN;NaN];
end;

for kk = 1:n_ima,
   if ~exist(['y_' num2str(kk)]),
   	eval(['y_' num2str(kk) ' = NaN*ones(2,1);']);
   end;
   if ~exist(['n_sq_x_' num2str(kk)]),
   	eval(['n_sq_x_' num2str(kk) ' = NaN;']);
   	eval(['n_sq_y_' num2str(kk) ' = NaN;']);
   end;
   if ~exist(['wintx_' num2str(kk)]),
   	eval(['wintx_' num2str(kk) ' = NaN;']);
   	eval(['winty_' num2str(kk) ' = NaN;']);
   end;
   if ~exist(['Tc_' num2str(kk)]),
    eval(['Tc_' num2str(kk) ' = [NaN;NaN;NaN];']);
   end;
   if ~exist(['omc_' num2str(kk)]),
    eval(['omc_' num2str(kk) ' = [NaN;NaN;NaN];']);
    eval(['Rc_' num2str(kk) ' = NaN * ones(3,3);']);
   end;
   if ~exist(['omc_error_' num2str(kk)]),
      eval(['omc_error_' num2str(kk) ' = NaN*ones(3,1);']);
   end;
   if ~exist(['Tc_error_' num2str(kk)]),
      eval(['Tc_error_' num2str(kk) ' = NaN*ones(3,1);']);
   end;
   if ~exist(['Rc_' num2str(kk)]),
    eval(['Rc_' num2str(kk) ' = rodrigues(omc_' num2str(kk) ');']);
   end;
end;

if ~exist('small_calib_image'),
    small_calib_image = 0;
end;

if ~exist('check_cond'),
    check_cond = 1;
end;

if ~exist('MaxIter'),
    MaxIter = 30;
end;

save_name = 'Calib_Results';

if exist([ save_name '.mat'])==2,
    disp('WARNING: File Calib_Results.mat already exists');
    if exist('copyfile'),
        pfn = -1;
        cont = 1;
        while cont,
            pfn = pfn + 1;
            postfix = ['_old' num2str(pfn)];
            save_name = [ 'Calib_Results' postfix];
            cont = (exist([ save_name '.mat'])==2);
        end;
        copyfile('Calib_Results.mat',[save_name '.mat']);
        disp(['Copying the current Calib_Results.mat file to ' save_name '.mat']);
        if exist('Calib_Results.m')==2,
            copyfile('Calib_Results.m',[save_name '.m']);
            disp(['Copying the current Calib_Results.m file to ' save_name '.m']);
        end;
        cont_save = 1;
    else
        disp('The file Calib_Results.mat is about to be changed.');
        cont_save = input('Do you want to continue? ([]=no,other=yes) ','s');
        cont_save = ~isempty(cont_save);
    end;
else
    cont_save = 1;
end;

if cont_save,

    save_name = 'Calib_Results';

    if exist('calib_name'),

        fprintf(1,['\nSaving calibration results under ' save_name '.mat\n']);

        string_save = ['save ' save_name ' center_optim param_list active_images ind_active est_alpha est_dist est_aspect_ratio est_fc fc kc cc alpha_c fc_error kc_error cc_error alpha_c_error  err_std ex x y solution solution_init wintx winty n_ima type_numbering N_slots small_calib_image first_num image_numbers format_image calib_name Hcal Wcal nx ny map dX_default dY_default KK inv_KK dX dY wintx_default winty_default no_image check_cond MaxIter'];

        for kk = 1:n_ima,
            string_save = [string_save ' X_' num2str(kk) ' x_' num2str(kk) ' y_' num2str(kk) ' ex_' num2str(kk) ' omc_' num2str(kk) ' Rc_' num2str(kk) ' Tc_' num2str(kk) ' omc_error_' num2str(kk) ' Tc_error_' num2str(kk) ' H_' num2str(kk) ' n_sq_x_' num2str(kk) ' n_sq_y_' num2str(kk) ' wintx_' num2str(kk) ' winty_' num2str(kk) ' dX_' num2str(kk) ' dY_' num2str(kk)];
        end;

    else

        fprintf(1,['\nSaving calibration results under ' save_name '.mat (no image version)\n']);

        string_save = ['save ' save_name ' center_optim param_list active_images ind_active est_alpha est_dist est_aspect_ratio est_fc fc kc cc alpha_c fc_error kc_error cc_error alpha_c_error err_std ex x y solution solution_init wintx winty n_ima nx ny dX_default dY_default KK inv_KK dX dY wintx_default winty_default no_image check_cond MaxIter'];

        for kk = 1:n_ima,
            string_save = [string_save ' X_' num2str(kk) ' x_' num2str(kk) ' y_' num2str(kk) ' ex_' num2str(kk) ' omc_' num2str(kk) ' Rc_' num2str(kk) ' Tc_' num2str(kk) ' omc_error_' num2str(kk) ' Tc_error_' num2str(kk) ' H_' num2str(kk) ' n_sq_x_' num2str(kk) ' n_sq_y_' num2str(kk) ' wintx_' num2str(kk) ' winty_' num2str(kk) ' dX_' num2str(kk) ' dY_' num2str(kk)];
        end;

    end;

    %fprintf(1,'To load later click on Load\n');

    eval(string_save);

    saving_calib_ascii;

    fprintf(1,'done\n');

end;

repo: ekfmonoslam


text:
function J = redblue(m,whiteloc)
% Usage: J=redblue([m],[whiteloc])
% redblue: A colormap from red to white to blue. Optionally specify the
% colormap size in "m" and the normalized location of the white point in
% "whiteloc", with 0 being the minimum and 1 being the maximum.
%
%   Example:
%      load clown
%      imagesc(X);
%      colormap redblue

% Written 17 July 2009 by Douglas H. Kelley, dhk [at] dougandneely.com.

whitelocdefault=0.5;

if ~exist('m','var') || isempty(m)
   m = size(get(gcf,'colormap'),1);
end
if ~exist('whiteloc','var') || isempty(whiteloc)
    whiteloc=whitelocdefault;
end

if floor(whiteloc*m)==0
    J=[ ones(m,1) linspace(1,0,m)' linspace(1,0,m)' ];
elseif floor(whiteloc*m)==m
    J=[ linspace(0,1,m)' linspace(0,1,m)' ones(m,1)];
else
    J=[ [linspace(0,1,floor(whiteloc*m))'; ones(m-floor(whiteloc*m),1)]...
        [linspace(0,1,floor(whiteloc*m))'; linspace(1,0,m-floor(whiteloc*m))']...
        [ones(floor(whiteloc*m),1); linspace(1,0,m-floor(whiteloc*m))'] ];
end

author: Sable

repo: mcbench-benchmarks


text:
function vw = rmSelect(vw, loadModel, rmFile)
% rmSelect - select retinotopy model file and put in view struct
%
% vw = rmSelect([vw=current view], [loadModel=0], [rmFile=dialog]);
%
% The loadModel flag indicates whether the model file (a large file)
% should be loaded or not. By default it is 0: don't load it until needed
% (only the path to the selected rmFile file will be stored). If it's 1, it
% goes ahead and loads it.
%
% If the path to the retinotopy model file is provided as a third argument,
% will attempt to load it directly; otherwise, pops up a dialog.  You
% may also provide the string 'mostrecent' as the filename, in which case
% the code will look for the most recently-created model file, and
% select it (producing an error if one is not found).
%
% 2006/02 SOD: wrote it.
% ras 2006/10: added dialog.
% ras, 06/07: added 'mostrecent' flag to find the newest file
if ~exist('vw','var') || isempty(vw), vw = getCurView;  end;
if ~exist('loadModel','var') || isempty(loadModel), loadModel = true;   end;

% choose filename:
if ~exist('rmFile','var') || isempty(rmFile),
    rmFile = getPathStrDialog(dataDir(vw),...
        'Choose retinotopic model file name', ...
        '*.mat');
    drawnow;
elseif iscell(rmFile)
    rmFile = rmFile{1};
end

% if user just wants the newest file, check for it:
if ischar(rmFile) && ismember(lower(rmFile), {'newest' 'mostrecent'})
	pattern = fullfile( dataDir(vw), 'retModel-*.mat' );
	w = dir(pattern);
	if isempty(w)
		error('Most Recent File selected; no retModel-* files found.')
	end
	[dates, order] = sortrows( datevec({w.date}) ); % oldest -> newest
	rmFile = fullfile( dataDir(vw), w(order(end)).name );
end

% Locate rm file
if check4File(fullfile(dataDir(vw), rmFile))
    % first check in appropriate directory
    rmFile = fullfile(dataDir(vw), rmFile);
elseif check4File(fullfile(dataDir(vw), [rmFile '.mat']))
    % if not found, check whether we are missing the extension
    rmFile = fullfile(dataDir(vw), [rmFile '.mat']);
elseif exist(rmFile,'file')
    % if still not found, look on entire path
elseif exist([rmFile '.mat'], 'file')
    % if still not found, look on entire path with '.mat' added
    rmFile = [rmFile '.mat'];
else
    % We could not find rm file anywhere
    fprintf('[%s]:No file: %s\n',mfilename,rmFile);
    return;
end

% if the load model flag is 1, but the file's already selected, just load
% it and return:
if loadModel && checkfields(vw, 'rm' , 'retinotopyModelFile')
    if ~exist(rmFile, 'file')
        error('Model file %s not found.', rmFile);
    end;
    load(rmFile, 'model', 'params');
    vw = viewSet(vw, 'rmFile', rmFile);
    vw = viewSet(vw, 'rmModel', model);
    vw = viewSet(vw, 'rmParams', params);
    vw = viewSet(vw, 'rmModelNum', 1);
    return;
end;

% store rmFile filename:
vw = viewSet(vw,'rmFile',rmFile);

if loadModel==0
    % clear previous models but don't load them untill we need them:
    vw = viewSet(vw, 'rmModel', []);
else
    % go ahead and load
    load(rmFile, 'model', 'params');
    vw = viewSet(vw, 'rmModel', model);
    vw = viewSet(vw, 'rmParams', params);
    vw = viewSet(vw, 'rmModelNum', 1);
end;

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%% Main function to generate tests
function tests = givens_vs_svd()
tests = functiontests(localfunctions);
end

%% Test Functions
function testQRvsSVD(nullspace)
% compare 3 methods to compute null space basis vectors: svd,
% matlab null which is said to use svd, and qr decomposition
% qr decomposition may be done with givens or householder rotations

trials=1e1;
timers=zeros(3,1);
for shot=1:trials
    Hf=rand(10,3);

    tic
    [u,d,v]= svd(Hf');
    n1= v(:, 4:end);
    timers(1) = timers(1) + toc;
    assert(length(find(d>0))==3);
    %svd to get nullspace basis, https://inst.eecs.berkeley.edu/~ee127a/book/login/l_svd_mat_prop.html

    tic
    n2=null(Hf');
    timers(2) = timers(2) + toc;
    assert(max(abs(n1 - n2), [], 'all') < 1e-6)

    tic
    % compute nullspace with qr decomposition
    % http://stackoverflow.com/questions/2181418/computing-the-null-space-of-a-matrix-as-fast-as-possible
    [q, r]= qr(Hf);
    n3=q(:,4:end);
    timers(3) = timers(3) + toc;

    assert(max(abs(n1 - n3), [], 'all') < 1e-6)
end
disp('times consumed for computing nullspace by svd, matlab null, and qr');
disp(timers)
end

function testQRvsSVD2(projectToNullspace)
% To project a matrix or vector onto the nullspace of another matrix is
% useful in eliminating unnecessary variables. e.g., r= H_x*X+H_f*p +n
% where H_f is of size 2M*3. p is of size 3*1. To remove p, we need to
% project r and H_x onto the nullspace of H_f', Q_2.
% Note H_f =[Q_1, Q_2][R; 0]= Q[R; 0]
%
% To accomplish this, we can successively apply Givens rotations on
% [r, H_x, H_f] until it becomes [Q'*r, Q'*H_x, [R; 0]]. And remove its first
% 3 rows, we get [Q_2'*r, Q_2'*H_x, 0].
%
% This trick is used in Mourikis, Anastasios, and Stergios Roumeliotis.
% "A multi-state constraint Kalman filter for vision-aided inertial navigation."
% Robotics and Automation, 2007 IEEE International Conference on. IEEE, 2007.
% How to apply Givens rotations can be found in
% G. Golub and C. van Loan, Matrix computations. The Johns Hopkins
% University Press, London, 1996.

trials=1e2;
timers=zeros(2,1);
% For small matrices, svd outperforms qr +svd; for large thin tall
% matrices, svd is a little less efficient than qr +svd
m =10; %3000;
n =5; %1200;
for shot=1:trials
    assert(m>n);
    A= rand(m, n);
    tic
    [~, ~, v] = svd(A);
    x1= v(:, end);
    timers(1) = timers(1) + toc;

    tic
    [~, r] = qr(A);
    Th= r(1:n, 1:n);
    [~, ~, v2]= svd(Th);
    x2= v2(:, end);
    timers(2) = timers(2) + toc;

    assert(max(abs(x1 - x2), [], 'all') < 1e-6)
end

disp('times consumed for solving least squares by svd, and qr+svd');
disp(timers)
end


text:
function Population = subNSGAIII(Population,N,Z,Zmin)
% The environmental selection of NSGA-III

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Cheng He

    if isempty(Zmin)
        Zmin = ones(1,size(Z,2));
    end
    CC = min([length(Population),N]);
    %% Non-dominated sorting
    [FrontNo,MaxFNo] = NDSort(Population.objs,Population.cons,CC);
    Next = FrontNo < MaxFNo;

    %% Select the solutions in the last front
    Last   = find(FrontNo==MaxFNo);
    Choose = LastSelection(Population(Next).objs,Population(Last).objs,CC-sum(Next),Z,Zmin);
    Next(Last(Choose)) = true;
    % Population for next generation
    Population = Population(Next);
end

function Choose = LastSelection(PopObj1,PopObj2,K,Z,Zmin)
% Select part of the solutions in the last front

    PopObj = [PopObj1;PopObj2] - repmat(Zmin,size(PopObj1,1)+size(PopObj2,1),1);
    [N,M]  = size(PopObj);
    N1     = size(PopObj1,1);
    N2     = size(PopObj2,1);
    NZ     = size(Z,1);

    %% Normalization
    % Detect the extreme points
    Extreme = zeros(1,M);
    w       = zeros(M)+1e-6+eye(M);
    for i = 1 : M
        [~,Extreme(i)] = min(max(PopObj./repmat(w(i,:),N,1),[],2));
    end
    % Calculate the intercepts of the hyperplane constructed by the extreme
    % points and the axes
    Hyperplane = PopObj(Extreme,:)\ones(M,1);
    a = 1./Hyperplane;
    if any(isnan(a))
        a = max(PopObj,[],1)';
    end
    % Normalization
    PopObj = PopObj./repmat(a',N,1);

    %% Associate each solution with one reference point
    % Calculate the distance of each solution to each reference vector
    Cosine   = 1 - pdist2(PopObj,Z,'cosine');
    Distance = repmat(sqrt(sum(PopObj.^2,2)),1,NZ).*sqrt(1-Cosine.^2);
    % Associate each solution with its nearest reference point
    [d,pi] = min(Distance',[],1);

    %% Calculate the number of associated solutions except for the last front of each reference point
    rho = hist(pi(1:N1),1:NZ);

    %% Environmental selection
    Choose  = false(1,N2);
    Zchoose = true(1,NZ);
    % Select K solutions one by one
    while sum(Choose) < K
        % Select the least crowded reference point
        Temp = find(Zchoose);
        Jmin = find(rho(Temp)==min(rho(Temp)));
        j    = Temp(Jmin(randi(length(Jmin))));
        I    = find(Choose==0 & pi(N1+1:end)==j);
        % Then select one solution associated with this reference point
        if ~isempty(I)
            if rho(j) == 0
                [~,s] = min(d(N1+I));
            else
                s = randi(length(I));
            end
            Choose(I(s)) = true;
            rho(j) = rho(j) + 1;
        else
            Zchoose(j) = false;
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [rhsu] = AdvecRHS3D(u,time)

% function [rhsu] = AdvecRHS3D(u,time)
% Purpose  : Evaluate RHS flux in 3D advection

Globals3D;

% form field differences at faces
alpha=1; du = zeros(Nfp*Nfaces,K);
du(:) = 0.5*(u(vmapM)-u(vmapP)).*(nx(:) - alpha*abs(nx(:)));

% impose boundary condition at x=0
ubc  = exp(-1*( (Fx(mapB)-time).^2 + Fy(mapB).^2 + Fz(mapB).^2));
du(mapB) = 0.5*(u(vmapB)-ubc).*(nx(mapB)-alpha*abs(nx(mapB)));

% compute right hand sides of the semi-discrete PDE
rhsu = -(rx.*(Dr*u)+sx.*(Ds*u)+tx.*(Dt*u)) + LIFT*(Fscale.*(du));
return

repo: Aero-matlab


text:
function dist = getSlipDist(D)

%This function computes the amount that the stick slid in each direction

N = length(D.phase);
dist = zeros(N,1);

for i=1:N
   if strcmp(D.phase{i},'SLIDE_POS')
       dist(i) = range(D.raw(i).state.x);
   elseif strcmp(D.phase{i},'SLIDE_NEG')
       dist(i) = -range(D.raw(i).state.x);
   end
end

end

author: MatthewPeterKelly

repo: dscTutorials


text:
clear all
close all

disp('Compute Model Evidence for Logistic Regression models using');
disp('Thermodynamic Integration');
disp('Explore effect of number of chains');

% Logistic model
[M{1},U{1}] = mci_logistic_struct ('dct');
P=[1 10 10]';
[g,Y]=mci_logistic_gen(P,M{1},U{1});
figure;
plot(g,'k','LineWidth',2);
set(gca,'FontSize',16);
grid on
xlabel('Time,t');
ylabel('p(Y=1|t)');

% MCMC parameters
mcmc.J=64; % Number of temperatures/chains
mcmc.ntune=500;
mcmc.nsamp=250;
setinit=1;
if setinit
    % Start sampling from this point
    for j=1:mcmc.J,
        mcmc.init{j}=P;
    end
    mcmc.nscale=0;
else
    mcmc.nscale=500;
end
mcmc.remove_burn_in=1;

% No sharing of samples between chains
mcmc.gprob=0;

j=[4,8,16,32,64];
for i=1:length(j),
    mcmc.J=j(i);
    tic;
    [P,logev,D] = spm_mci_pop (mcmc,M,U,Y);
    els(i)=toc;
    L(i)=logev.ti;
    disp(sprintf('J=%d: TI Log Evidence = %1.2f',j(i),logev.ti));
end

figure
plot(j,L,'k','LineWidth',2);
grid on
set(gca,'FontSize',16);
xlabel('Number of Chains');
ylabel('Log Evidence');
title('Thermodynamic Integration');

% Annealed Importance Sampling
mcmc.inference='ais';
mcmc.anneal='geometric';
mcmc.prop='lmc';
mcmc.nprop=1;
mcmc.maxits=64;

j=[4,8,16,32,64,128,256,512,1024];
for i=1:length(j),
    mcmc.J=j(i);
    tic;
    post = spm_mci_ais (mcmc,M{1},U{1},Y);
    els(i)=toc;
    L(i)=post.logev;
    disp(sprintf('J=%d: AIS Log Evidence = %1.2f',j(i),L(i)));
end

figure
semilogx(j,L,'k','LineWidth',2);
grid on
set(gca,'FontSize',16);
xlabel('Number of Chains');
ylabel('Log Evidence');
title('Annealed Importance Sampling');

author: spm


text:
function [sys,x0,str,ts]=Nonlinear_PD(t,x,u,flag,Beta1,Beta2,Delta,a1,a2,h)

if flag==0

    sys = [0;1;1;2;0;1;1];
     x0 = [0];
    str = [ ];
     ts = [h 0];

elseif flag==2

    Bili=fal(u(1),a1,Delta);
    Weifen=fal(u(2),a2,Delta);
    sys=Beta1*Bili+Beta2*Weifen;

elseif flag==3
    sys=x;

elseif flag==4
    sys=sys+h;

else
    sys=[];
end
function f=fal(e,a,d)
    if abs(e)<d
        f=e*d^(a-1);
    else f=(abs(e))^a*sign(e);
    end

author: TianfaYao

repo: ADRC

save_path: github-repos/MATLAB/TianfaYao-ADRC


text:
function [L_b,lambda_b,h_b,v_eb_n,C_b_n] = ECEF_to_NED(r_eb_e,v_eb_e,C_b_e)
%ECEF_to_NED - Converts Cartesian  to curvilinear position, velocity
%resolving axes from ECEF to NED and attitude from ECEF- to NED-referenced
%
% Software for use with "Principles of GNSS, Inertial, and Multisensor
% Integrated Navigation Systems," Second Edition.
%
% This function created 2/4/2012 by Paul Groves
%
% Inputs:
%   r_eb_e        Cartesian position of body frame w.r.t. ECEF frame, resolved
%                 along ECEF-frame axes (m)
%   v_eb_e        velocity of body frame w.r.t. ECEF frame, resolved along
%                 ECEF-frame axes (m/s)
%   C_b_e         body-to-ECEF-frame coordinate transformation matrix
%
% Outputs:
%   L_b           latitude (rad)
%   lambda_b      longitude (rad)
%   h_b           height (m)
%   v_eb_n        velocity of body frame w.r.t. ECEF frame, resolved along
%                 north, east, and down (m/s)
%   C_b_n         body-to-NED coordinate transformation matrix

% Parameters
R_0 = 6378137; %WGS84 Equatorial radius in meters
% e = 0.0818191908425; %WGS84 eccentricity

ecc2=0.006694380022900;

% Copyright 2012, Paul Groves
% License: BSD; see license.txt for details

% Begins

% Convert position using Borkowski closed-form exact solution
% From (2.113)
lambda_b = atan2(r_eb_e(2),r_eb_e(1));

% From (C.29) and (C.30)
k1 = sqrt(1 - ecc2) * abs (r_eb_e(3));
k2 = ecc2 * R_0;
beta = sqrt(r_eb_e(1)^2 + r_eb_e(2)^2);
E = (k1 - k2) / beta;
F = (k1 + k2) / beta;

% From (C.31)
P = 4/3 * (E*F + 1);

% From (C.32)
Q = 2 * (E^2 - F^2);

% From (C.33)
D = P^3 + Q^2;

% From (C.34)
V = (sqrt(D) - Q)^(1/3) - (sqrt(D) + Q)^(1/3);

% From (C.35)
G = 0.5 * (sqrt(E^2 + V) + E);

% From (C.36)
T = sqrt(G^2 + (F - V * G) / (2 * G - E)) - G;

% From (C.37)
L_b = sign(r_eb_e(3)) * atan((1 - T^2) / (2 * T * sqrt (1 - ecc2)));

% From (C.38)
h_b = (beta - R_0 * T) * cos(L_b) +...
    (r_eb_e(3) - sign(r_eb_e(3)) * R_0 * sqrt(1 - ecc2)) * sin (L_b);

% Calculate ECEF to NED coordinate transformation matrix using (2.150)
cos_lat = cos(L_b);
sin_lat = sin(L_b);
cos_long = cos(lambda_b);
sin_long = sin(lambda_b);
C_e_n = [-sin_lat * cos_long, -sin_lat * sin_long,  cos_lat;...
                   -sin_long,            cos_long,        0;...
         -cos_lat * cos_long, -cos_lat * sin_long, -sin_lat];

% Transform velocity using (2.73)
v_eb_n = C_e_n * v_eb_e;

% Transform attitude using (2.15)
C_b_n = C_e_n * C_b_e;

% Ends

author: benzenemo

repo: TightlyCoupledINSGNSS


text:
function y = fix_filter(x)
% fix_filter (in case signal processing toolbox is not available).
% filters data x with an eliptic passband between [300 3000] Hz.

a = [1.0000 -2.3930  2.0859 -0.9413 0.2502];
b = [0.1966 -0.0167 -0.3598 -0.0167 0.1966];

x = x(:);
len = size(x,1);
b = b(:).';
a = a(:).';
nb = length(b);
na = length(a);
nfilt = max(nb,na);

nfact = 3*(nfilt-1);  % length of edge transients

rows = [1:nfilt-1  2:nfilt-1  1:nfilt-2];
cols = [ones(1,nfilt-1) 2:nfilt-1  2:nfilt-1];
data = [1+a(2) a(3:nfilt) ones(1,nfilt-2)  -ones(1,nfilt-2)];
sp = sparse(rows,cols,data);
zi = sp \ ( b(2:nfilt).' - a(2:nfilt).'*b(1) );

y = [2*x(1)-x((nfact+1):-1:2);x;2*x(len)-x((len-1):-1:len-nfact)];

if exist('FilterM','file')
    y = FilterM(b,a,y,[zi*y(1)]);
else
    y = filter(b,a,y,[zi*y(1)]);
end
y = y(length(y):-1:1);

%second filter, in the other way
if exist('FilterM','file')
    y = FilterM(b,a,y,[zi*y(1)]);
else
    y = filter(b,a,y,[zi*y(1)]);
end
y = y(length(y):-1:1);

y([1:nfact len+nfact+(1:nfact)]) = [];

y = y.';

author: csn-le


text:
function val = get(f, propName)
%GET    GET method for SPHEREFUN class.
%   P = GET(F, PROPNAME) returns the property P specified in the string
%   PROPNAME from the SPHEREFUN object F. Valid entries for the string
%   PROPNAME are:
%    'DOMAIN'
%    'COLS'
%    'ROWS'
%    'PIVOTVALUES'
%    'PIVOTLOCATIONS'
%    'BLOCKDIAG'
%    'IDXPLUS'
%    'IDXMINUS'
%    'NONZEROPOLES'

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

% Loop through an array of SPHEREFUN objects.
if ( numel(f) > 1 )
    val = cell(numel(f));
    for k = 1:numel(f)
        val{k} = get(f(k), propName);
    end
    return
end

% Get the properties.
switch ( propName )
    case 'domain'
        val = f.domain;
    case 'cols'
        val = f.cols;
    case 'rows'
        val = f.rows;
    case 'pivotValues'
        val = f.pivotValues;
    case 'pivotLocations'
        val = f.pivotLocations;
    case 'blockDiag'
        val = f.blockDiag;
    case 'idxPlus'
        val = f.idxPlus;
    case 'idxMinus'
        val = f.idxMinus;
    case 'nonZeroPoles'
        val = f.nonZeroPoles;
    otherwise
        error('CHEBFUN:SPHEREFUN:get:propName', ...
            [propName,' is not a valid SPHEREFUN property.'])
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function kern = ndsimKernExpandParamTransformSettings(kern, paramtransformsettings)

% NDSIMKERNEXPANDPARAMTRANSFORMSETTINGS Create kernel structure from SIM kernel's parameters' transform settings.
% FORMAT
% DESC returns a single input motif kernel structure filled with the
% parameters in the given vector. This is used as a helper function to
% enable parameters to be optimised in, for example, the NETLAB
% optimisation functions.
% ARG kern : the kernel structure in which the parameters are to be
% placed.
% ARG param : vector of parameters which are to be placed in the
% kernel structure.
% RETURN kern : kernel structure with the given parameters in the
% relevant locations.
%
% SEEALSO : simKernParamInit, simKernExtractParam, kernExpandParam
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% COPYRIGHT : Jaakko Peltonen, 2011

% KERN

%fprintf(1,'expanding SIM kernel param-transforms settings:\n');
%paramtransformsettings
%for k=1:length(paramtransformsettings),
%  paramtransformsettings{k}
%end;
%pause

if isfield(kern, 'options') ...
      && isfield(kern.options, 'isNegativeS') ...
      && kern.options.isNegativeS,
  if isfield(kern, 'gaussianInitial') && kern.gaussianInitial,
    % positive transformations for: inversewidth, initial variance
    for k=1:2,
      kern.transforms(k).transformsettings = paramtransformsettings{k};
    end;
  else
    % positive transformations for: inversewidth
    for k=1:1,
      kern.transforms(k).transformsettings = paramtransformsettings{k};
    end;
  end;

else
  if isfield(kern, 'gaussianInitial') && kern.gaussianInitial,
    % positive transformations for: inversewidth, kernel-variance,
    % initial variance
    for k=1:3,
      kern.transforms(k).transformsettings = paramtransformsettings{k};
    end;
  else
    % positive transformations for: inversewidth, kernel-variance
    for k=1:2,
      kern.transforms(k).transformsettings = paramtransformsettings{k};
    end;
  end;
end;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%% This script plots the loss and accuracy evolution during the
%%% training of a CNN by Caffe

% Path to the .txt file with the training output
%%% ObjDetectCNN
% file_path = '../Training_Results/output_training_ObjDetection_v1.txt';
% file_path = '../Training_Results/output_training_ObjDetection_finetunning_v2.2.txt';
% file_path = '../Training_Results/output_training_ObjDetection_finetunning_strict_v1.txt';
%%% MammoCNN
% file_path = '../../../MamoCNN/Training_Results/output_training_MammoCNN_v2.txt';
% file_path = '../../../MamoCNN/Training_Results/output_training_lr0-000000001.txt';
% file_path = '../../../MamoCNN/Training_Results/output_training_2-class_lr0-0000001.txt';
% file_path = '../../../MamoCNN/Training_Results/output_training_2-class_lr0-000001_3conv.txt';
% file_path = '../../../MamoCNN/Training_Results/output_training_2-class_lr0-0000001_2conv.txt';
% file_path = '../../../MamoCNN/Training_Results/output_training_5-class_finetunning_lr0-0000001.txt';
%%% FoodCNN
% file_path = '../../../FoodCNN/Training_Results/output_training_finetunning_v1.txt';
% file_path = '../../../FoodCNN/Training_Results/output_training_v1.txt';
% file_path = '../../../FoodCNN/Training_Results/output_training_finetunning_v2.txt';
% file_path = '../../../FoodCNN/Training_Results/output_training_finetunning_v2_1.txt';
%%% InformativeImagesDetectorCNN
% file_path = '../../../InformativeImagesDetector/Training_Results/output_training_finetunning_InformativeCNN_CV1_v1.txt';
% file_path = '../../../InformativeImagesDetector/Training_Results/output_training_finetunning_InformativeCNN_CV1_v2.txt';
% file_path = '../../../InformativeImagesDetector/Training_Results/output_training_finetunning_HybridPlaces_InformativeCNN_CV1_v1.txt';
% file_path = '../../../InformativeImagesDetector/Training_Results/output_training_finetunning_CNDS_InformativeCNN_CV1_v1.txt';
% file_path = '../../../InformativeImagesDetector/Training_Results/output_training_finetunning_CNDS_InformativeCNN_CV1_v2.txt';
file_path = '../../../InformativeImagesDetector/Training_Results/output_training_finetunning_CNDS_InformativeCNN_CV1_v3.txt';
% file_path = '../../../InformativeImagesDetector/Training_Results/output_training_finetunning_InformativeCNN_CV2_v1.txt';
% file_path = '../../../InformativeImagesDetector/Training_Results/output_training_finetunning_InformativeCNN_CV3_v1.txt';
% file_path = '../../../InformativeImagesDetector/Training_Results/output_training_finetunning_InformativeCNN_CV4_v1.txt';
%%% MNIST Test
% file_path = '../../../MNIST_tests/Training_Results/output_training_v1.txt';
% file_path = '../../../MNIST_tests/Training_Results/output_training_mod_v1.txt';
% file_path = '../../../MNIST_tests/Training_Results/output_training_mod_v2.txt';
% file_path = '../../../MNIST_tests/Training_Results/output_training_mod_v2_2-layers.txt';
% file_path = '../../../MNIST_tests/Training_Results/output_training_mod_v2_4-layers.txt';

test_accuracies = {'accuracy', 'accuracy_DSN_conv3'};
% test_accuracies = {'accuracy'};
% test_accuracies = {'prob'};

test_losses = {'loss', 'loss_DSN_conv3'};
% test_losses = {'loss'};

% Only pick 1 sample for each N
Nsubsample_loss = 5;
% Nsubsample_loss = 5;
Nsubsample_axis = 50;
% Nsubsample_axis = 50;
Nsubsample_accuracy = 1;

% Training loss, Test loss, Test accuracy and Max Accuracy colours
colours = {'k', 'winter', 'winter', 'r'};
lines_width = 2;

%% Get colours

%%% Test loss
test_loss_c = colormap(colours{2});
close(gcf);
test_loss_c = test_loss_c(round(linspace(1,size(test_loss_c,1)/2, length(test_losses))), :);

%%% Test accuracy
test_acc_c = colormap(colours{3});
close(gcf);
test_acc_c = test_acc_c(round(linspace(size(test_acc_c,1)-20,size(test_acc_c,1), length(test_accuracies))), :);

%% Read file
data = fileread(file_path);

%% Prepare figure
f = figure;
hold on;

%% Get training progress
loss = [];
loss_iter = [];

%%% Training loss
find_loss = regexp(data, 'Train net output #0: loss = ', 'split');
nSplits = length(find_loss);
for i = 2:nSplits
    this_loss = regexp(find_loss{i}, ' ', 'split');
    loss = [loss str2num(this_loss{1})];
end

%%% Training iteration
find_loss_iter = regexp(data, 'Iteration ', 'split');
nSplits = length(find_loss_iter);
for i = 2:nSplits
    this_loss_iter = regexp(find_loss_iter{i}, ', loss = ', 'split');
    if(length(this_loss_iter) > 1)
        loss_iter = [loss_iter str2num(this_loss_iter{1})];
    end
end

if(length(loss) == length(loss_iter)-1)
    loss_iter = loss_iter(1:end-1);
end

%% Get testing progress
accuracy = [];
loss_test = [];
accuracy_loss_iter = [];

%%% Testing accuracy
nTestAcc = length(test_accuracies);
for n = 1:nTestAcc
    find_accuracy = regexp(data, ['Test net output #.: ' test_accuracies{n} ' = '], 'split');
    nSplits = length(find_accuracy);
    for i = 2:nSplits
        this_accuracy = regexp(find_accuracy{i}, '\n', 'split');
        accuracy(n,i-1) = str2num(this_accuracy{1});
    end
end

%%% Testing loss
nTestLoss = length(test_losses);
for n = 1:nTestLoss
    find_loss_test = regexp(data, ['Test net output #.: ' test_losses{n} ' = '], 'split');
    nSplits = length(find_loss_test);
    for i = 2:nSplits
        this_loss_test = regexp(find_loss_test{i}, ' ', 'split');
        loss_test(n,i-1) = str2num(this_loss_test{1});
    end
end

%%% Testing iteration
find_accuracy_iter = regexp(data, ', Testing net (#0', 'split');
nSplits = length(find_accuracy_iter);
for i = 1:nSplits-1
    this_accuracy_iter = regexp(find_accuracy_iter{i}, 'Iteration ', 'split');
    accuracy_loss_iter = [accuracy_loss_iter str2num(this_accuracy_iter{end})];
end

%% Plot loss progress
data_plotted = {};

plot(loss_iter(1:Nsubsample_loss:end), loss(1:Nsubsample_loss:end), 'Color', colours{1}, 'LineWidth', lines_width);
data_plotted = {data_plotted{:}, 'Training loss'};

%% Plot accuracy progress
for n = 1:nTestLoss
    plot(accuracy_loss_iter(1:Nsubsample_accuracy:end), loss_test(n, 1:Nsubsample_accuracy:end), 'Color', test_loss_c(n,:), 'LineWidth', lines_width);
    data_plotted = {data_plotted{:}, ['Test ' strjoin(regexp(test_losses{n}, '_', 'split'), ' ')]};
end
for n = 1:nTestAcc
    plot(accuracy_loss_iter(1:Nsubsample_accuracy:end), accuracy(n, 1:Nsubsample_accuracy:end), 'Color', test_acc_c(n,:), 'LineWidth', lines_width);
    data_plotted = {data_plotted{:}, ['Test ' strjoin(regexp(test_accuracies{n}, '_', 'split'), ' ')]};
end

%% Plot max accuracy horizontal line
plot([0 loss_iter(end)], [1 1], 'Color', colours{4}, 'LineWidth', lines_width);
data_plotted = {data_plotted{:}, 'Max Accuracy'};

%% Plot max accuracy and min test loss
for n = 1:nTestAcc
    [val, pos] = max(accuracy(n,:));
    scatter(accuracy_loss_iter(pos), val, [], ...
                    'MarkerEdgeColor','k',...
                    'MarkerFaceColor',test_acc_c(n,:),...
                    'LineWidth',1.5);
    text(accuracy_loss_iter(pos)+0.02, val, sprintf('acc %0.3f\niter %d', val, accuracy_loss_iter(pos)), 'FontSize', 12);
end

for n = 1:nTestLoss
    [val, pos] = min(loss_test(n,:));
    scatter(accuracy_loss_iter(pos), val, [], ...
                    'MarkerEdgeColor','k',...
                    'MarkerFaceColor',test_loss_c(n,:),...
                    'LineWidth',1.5);
    text(accuracy_loss_iter(pos)+0.02, val, sprintf('loss %0.3f\niter %d', val, accuracy_loss_iter(pos)), 'FontSize', 12);
end

%% Plot general info
xticklabel_rotate(loss_iter(1:Nsubsample_axis:end), 45);
legend(data_plotted,3);
max_val_loss = max([max(loss) max(loss_test) 1]);
ylim([0 max_val_loss]);
set(gca, 'YTick', linspace(0, max_val_loss, 21));
xlabel('Iterations');
grid on;

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [t,X] = smdplot(t,X)
% See Article 2.6
% [t,X]=smdplot(t,X)
% This function plots the response and animates the
% motion of a damped linear harmonic oscillator
% characterized by the differential equation
% m*x''+c*x'+k*x=f1*cos(w*t)+f2*sin(w*t)
% with initial conditions x(0)=x0, x'(0)=v0.
% The animation depicts forced motion of a block
% attached to a wall by a spring. The block
% slides on a horizontal plane which provides
% viscous damping.

% Default data case used when no input arguments
% are given
if nargin==0
  m=1; c=.3; k=1; f1=1.5; f2=0; w=2; x0=0; v0=2;
  %tmax=30; nt=250; t=linspace(0,tmax,nt);
  tmax=25; nt=200; t=linspace(0,tmax,nt);
  X=smdsolve(m,c,k,f1,f2,w,x0,v0,t);
end

% Plot the displacement versus time
plot(t,X), xlabel('time'), ylabel('displacement')
title(...
'FORCED RESPONSE OF A DAMPED HARMONIC OSCILLATOR')
grid on, shg, pause(3)

% Add a block and a spring to the displacement
xmx=max(abs(X)); X=X/1.1/xmx;
xb=[0,0,1,1,0,0]/2; yb=[0,-1,-1,1,1,0]/2;

% Make an arrow tip
d=.08; h=.05;
xtip=[0,-d,-d,0]; ytip=[0,0,0,h,-h,0];

% Add a spring and a block to the response
[xs,ys]=spring; nm=length(X); ns=length(xs);
nb=length(xb); x=zeros(nm,ns+nb);y=[ys,yb];
for j=1:nm, x(j,:)=[-1+(1+X(j))*xs,X(j)+xb];end
r=[min(x(:)),max(x(:))+1,-2,2];
rx=r([1 1 2]); ry=[.5,-.5,-.5]; close;

% Plot the motion
for j=1:nm
   % Compute and scale the applied force
   f=f1*cos(w*t(j))+f2*sin(w*t(j));
   f=.3*f; fa=abs(f); sf=sign(f);
   xj=x(j,:); xmaxj=max(xj);
   if sf>0
      xforc=xmaxj+[0,fa,fa+xtip];
   else
      xforc=xmaxj+[fa,0,-xtip];
   end

   % Plot the spring, block, and force
   plot(xj,y,rx,ry,'k',xforc,ytip,'r')
   title('FORCED MOTION WITH DAMPING')
   axis(r), axis('off'), drawnow
   figure(gcf), pause(.05)
end

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

function [x,y] = spring(len,ht)
% This function generates a set of points
% defining a spring

if nargin==0, len=1; ht=.125; end
x=[0,.5,linspace(1,11,10),11.5,12];
y=[ones(1,5);-ones(1,5)];
y=[0;0;y(:);0;0]'; y=ht/2/max(y)*y;
x=len/max(x)*x;

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

function [x,v]=smdsolve(m,c,k,f1,f2,w,x0,v0,t)
% [x,v]=smdsolve(m,c,k,f1,f2,w,x0,v0,t)
% This function solves the differential equation
% m*x''(t)+c*x'(t)+k*x(t)=f1*cos(w*t)+f2*sin(w*t)
% with x(0)=x0 and x'(0)=v0
%
% m,c,k  - mass, damping and stiffness coefficients
% f1,f2  - magnitudes of cosine and sine terms in
%          the forcing function
% w      - frequency of the forcing function
% t      - vector of times to evaluate the solution
% x,v    - computed position and velocity vectors

ccrit=2*sqrt(m*k); wn=sqrt(k/m);

% If the system is undamped and resonance will
% occur, add a little damping
if c==0 & w==wn; c=ccrit/1e6; end;

% If damping is critical, modify the damping
% very slightly to avoid repeated roots
zeta=c/ccrit; if zeta==1, zeta=zeta+1e-6; end

% Forced response solution
a=(f1-i*f2)/(k-m*w^2+i*c*w);
X0=real(a); V0=real(i*w*a);
X=real(a*exp(i*w*t)); V=real(i*w*a*exp(i*w*t));

% Homogeneous solution
r=sqrt(zeta^2-1); s1=wn*(-zeta+r); s2=wn*(-zeta-r);
p=[1,1;s1,s2]\[x0-X0;v0-V0];

% Total solution satisfying the initial conditions
x=X+real(p(1)*exp(s1*t)+p(2)*exp(s2*t));
v=V+real(p(1)*s1*exp(s1*t)+p(2)*s2*exp(s2*t));


text:
function r = plus(a,b)
%PLUS         Taylor addition  a + b
%

% written  05/21/09     S.M. Rump
%

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

  if ~isa(a,'taylor')               % non-Taylor plus Taylor
    r = b;
    if isa(a,'intval')
      r.t = intval(r.t);
    end
    if isscalar(a)
      r.t(1,:) = r.t(1,:) + a;
    else
      sizeb = prod(b.size);
      if prod(sizeb)==1
        r.size = size(a);
        r.t = repmat(r.t,prod(r.size));
        r.t(1,:) = r.t(1,:) + a(:).';
      else
        if ~isequal(size(a),b.size)
          error('operands of different size')
        end
        r.t(1,:) = r.t(1,:) + a(:).';
      end
    end
  elseif ~isa(b,'taylor')           % Taylor plus non-Taylor
    r = a;
    if isa(b,'intval')
      r.t = intval(r.t);
    end
    if isscalar(b)
      r.t(1,:) = r.t(1,:) + b;
    else
      sizea = prod(a.size);
      if prod(sizea)==1
        r.size = size(b);
        r.t = repmat(r.t,prod(r.size));
        r.t(1,:) = r.t(1,:) + b(:).';
      else
        if ~isequal(size(b),a.size)
          error('operands of different size')
        end
        r.t(1,:) = r.t(1,:) + b(:).';
      end
    end
  else                              % Taylor plus Taylor
    r = a;
    if isa(b.t,'intval')
      r.t = intval(r.t);
    end
    sa = prod(a.size);
    sb = prod(b.size);
    if sa==1                        % a is scalar
      if sb~=1                      % b is not scalar
        r.size = b.size;
        a.t = repmat(a.t,1,sb);
      end
    else                            % a is not scalar
      if sb==1                      % b is scalar
        b.t = repmat(b.t,1,sa);
      else
        if ~isequal(a.size,b.size)
          error('operands of different size')
        end
      end
    end
    r.t = a.t + b.t;
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
%
% References:
%
% C. Lu. A Library of ADMM for Sparse and Low-rank Optimization. National University of Singapore, June 2016.
% https://github.com/canyilu/LibADMM.
% C. Lu, J. Feng, S. Yan, Z. Lin. A Unified Alternating Direction Method of Multipliers by Majorization
% Minimization. IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 40, pp. 527-541, 2018
%

addpath(genpath(cd))
clear

%% Examples for testing the sparse models
% For detailed description of the sparse models, please refer to the Manual.

%% generate toy data
d = 10;
na = 200;
nb = 100;

A = randn(d,na);
X = randn(na,nb);
B = A*X;
b = B(:,1);

opts.tol = 1e-6;
opts.max_iter = 1000;
opts.rho = 1.1;
opts.mu = 1e-4;
opts.max_mu = 1e10;
opts.DEBUG = 0;

%% l1
[X2,obj,err,iter] = l1(A,B,opts);
iter
obj
err
stem(X2(:,1))

%% group l1
g_num = 5;
g_len = round(na/g_num);
for i = 1 : g_num-1
    G{i} = (i-1)*g_len+1 : i*g_len;
end
G{g_num} = (g_num-1)*g_len+1:na;

[X2,obj,err,iter] = groupl1(A,B,G,opts);
iter
obj
err
stem(X2(:,1))

%% elastic net
lambda = 0.01;
[X2,obj,err,iter] = elasticnet(A,B,lambda,opts);
iter
obj
err
stem(X2(:,1))

%% fused Lasso
lambda = 0.01;
[x,obj,err,iter] = fusedl1(A,b,lambda,opts);
iter
obj
err
stem(x)

%% trace Lasso
[x,obj,err,iter] = tracelasso(A,b,opts);
iter
obj
err
stem(x)

%% k-support norm
k = 10;
[X,err,iter] = ksupport(A,B,k,opts);
iter
err
stem(X(:,1));

%% --------------------------------------------------------------

%% regularized l1
lambda = 0.01;
opts.loss = 'l1';
[X,E,obj,err,iter] = l1R(A,B,lambda,opts);
iter
obj
err
stem(X(:,1))

%% regularized group Lasso
g_num = 5;
g_len = round(na/g_num);

for i = 1 : g_num-1
    G{i} = (i-1)*g_len+1 : i*g_len;
end
G{g_num} = (g_num-1)*g_len+1:na;
lambda = 1;
opts.loss = 'l1';
[X,E,obj,err,iter] = groupl1R(A,B,G,lambda,opts);
iter
obj
err
stem(X(:,1))

%% regularized elastic net
lambda1 = 10;
lambda2 = 10;
opts.loss = 'l1';
[X,E,obj,err,iter] = elasticnetR(A,B,lambda1,lambda2,opts);
iter
obj
err
stem(X(:,1))
% stem(E(:,1))

%% regularized fused Lasso
lambda1 = 10;
lambda2 = 10;
opts.loss = 'l1';
[X,E,obj,err,iter] = fusedl1R(A,b,lambda1,lambda2,opts);
iter
obj
err
stem(X(:,1))
stem(E(:,1))

%% regularized trace Lasso
lambda = 0.1;
opts.loss = 'l1';
tic
[x,e,obj,err,iter] = tracelassoR(A,b,lambda,opts);
toc
iter
obj
err
stem(x)

%% regularized k-support norm
lambda = 0.1;
k = 10;
[X,E,err,iter] = ksupportR(A,B,lambda,k,opts);
iter
err
stem(X(:,1));

author: canyilu

repo: LibADMM-toolbox


text:
function [subs,vals] = find(t)
%FIND Find subscripts of nonzero elements in a sparse tensor.
%
%   [SUBS,VALS] = FIND(T) returns the subscripts and corresponding
%   values of the nonzero elements of T.
%
%   Note that unlike the standard MATLAB find function for an array,
%   find does not return linear indices. Instead, it returns an M x N
%   array where M is the number of nonzero values and N = ndims(T).
%   Thus, I(k,:) specifies the subscript of value V(k).
%
%   See also SPTENSOR, FIND.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

subs = t.subs;
vals = t.vals;


text:
function out = binavgcol(in, numinbin)
% Rebins image by averaging numinbin bins together along columns of the
% 2D data.
% out = binavgcol(in, numinbin)

di=0;
if strcmp(class(in),'dip_image')
    di=1;
    in=double(in);
end
tmp = in;

for ii=2:numinbin
    tmp = tmp + circshift(in,[-1*(ii-1) 0]);
end

nb = floor(size(in,1)/numinbin);
ix=numinbin*(1:nb)-(numinbin-1);
out = tmp(ix,:)/numinbin;
if di
    out=dip_image(out);
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function check_installation
% check whether all binaries are working properly

% check nfft
if ~check_nfft

  rmpath([mtex_path filesep 'extern' filesep 'nfft_openMP'])
  addpath([mtex_path filesep 'extern' filesep 'nfft'])

  hline()
  if check_nfft

    disp(' ');
    disp('  MTEX: Error running NFFT with openmp - I''m going to use it without openmp!');

  else

    disp(' ');
    disp('  MTEX: Error running NFFT!');
    disp('   ');
    disp('  The original error message was:');
    e = lasterror;
    disp(['  ' e.message]);

    disp(' ');
    disp(' I did not get NFFT working. This restricts the functionality of MTEX.')
    disp(' To overcome this restriction you may need to compile the NFFT your own.');
    if ismac
      disp(' Please have a look at https://mtex-toolbox.github.io/installation for possible workarounds')
    else
      disp(' Please have a look at https://github.com/mtex-toolbox/mtex/blob/develop/extern/nfft_openMP/readme.md');
    end

  end
  hline()
end

% check mex
check_mex;

end

% -------- check nfft

function e = check_nfft

e = 0;
try
  % check NFSFT
  f = @(x) x.theta;
  sF = S2FunHarmonic.quadrature(f,'bandwidth',12);
  assert(abs(f(xvector) - sF.eval(xvector))<1e-5);

  % check NFSOFT
  %ori = orientation.rand(100,crystalSymmetry);
  %odf = calcDensity(ori,'Fourier','bandwidth',12);
  %odf.eval(ori);
  e = 1;
end
end

% ----------- check mex files ---------------------------
function check_mex

% check for mex files
if fast_check_mex, return;end

hline('-');
disp('Checking mex files failed!');

disp('> Trying now to recompile mex files.');
opwd = pwd;
cd(fullfile(mtex_path,'mex'));
mex_install(mtex_path,false);
cd(opwd);

if fast_check_mex
  disp('> Hurray - everythink works!')
else
  hline()
  disp('  MTEX: Couldn''t get the mex files working!');
  disp('   ');
  disp('  The original error message was:');
  e = lasterror;
  disp(['  ' e.message]);

  disp(' ');
  disp(' Please have a look at https://mtex-toolbox.github.io/installation')
  disp(' or ask in the forum https://groups.google.com/forum/#!forum/mtexmail for help.');

  hline()
end
hline('-')
end

% ----------------------------------------------------------------------

function e = fast_check_mex

e = false;
mexpath = [mtex_path filesep 'mex'];

% check for existence of the files
mex = {'quaternion_*','S1Grid_*','S2Grid_*','SO3Grid_*','insidepoly_*','representationbased*'};
for k=1:numel(mex)
  cfile = dir(fullfile(mexpath, [mex{k} '.c']));
  for c=1:numel(cfile)
    [~,cname] = fileparts(cfile(c).name);
    mexfile = fullfile(mexpath, [cname,'.' mexext]);
    if ~exist(mexfile,'file')
       warning(['missing mex-file ' [cname,'.' mexext]])
      return
    end
  end
end

try
  S3G = equispacedSO3Grid(crystalSymmetry,specimenSymmetry,'points',100);
  dot_outer(S3G,quaternion.id,'epsilon',pi/4);

  SO3F = SO3FunHarmonic(rand(1000,1));
  SO3F.eval(rotation.rand(10));

  e = true;
catch
  e = false;
end

end

function hline(st)

if nargin < 1, st = '*'; end
fprintf('\n');
disp(repmat(st,1,80));

end

author: mtex-toolbox

repo: mtex

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


text:
function [theta] = normalEqn(X, y)
    %% NORMALEQN Computes the closed-form solution to linear regression
    %   NORMALEQN(X,y) computes the closed-form solution to linear
    %   regression using the normal equations.
    theta = pinv(X' * X) * X' * y;
end

author: worldveil

repo: coursera-ml


text:
function Header = spm_dicom_header(DicomFilename, DicomDictionary, Options)
% Read header information from a DICOM file
% FORMAT Header = spm_dicom_header(DicomFilename, DicomDictionary, Options)
% DicomFilename   - DICOM filename
% DicomDictionary - DICOM dictionary (see spm_dicom_headers)
% Options         - an (optional) structure containing fields
%                   abort      - if this is a function handle, it will
%                                be called with field name and value
%                                arguments.  If this function returns true,
%                                then reading the header will be aborted.
%                                [Default: false]
%                   all_fields - binary true/false, indicating what to do
%                                with fields that are not included in the
%                                DICOM dictionary.
%                                [Default: true]
%
% Header          - Contents of DICOM header
%
% Contents of headers are approximately explained in:
% http://medical.nema.org/standard.html
%
% This code may not work for all cases of DICOM data, as DICOM is an
% extremely complicated "standard".
%__________________________________________________________________________
% Copyright (C) 2002-2018 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_dicom_header.m 7374 2018-07-09 17:09:46Z guillaume $

if nargin < 3
    Options = struct('abort',false, 'all_fields',true);
end

Header = [];
DicomFilename = deblank(DicomFilename);
FID  = fopen(DicomFilename,'r','ieee-le');
if FID == -1
    warning('spm:dicom','%s: Cant open file.', DicomFilename);
    return;
end

fseek(FID,128,'bof');
dcm = char(fread(FID,4,'uint8')');
if ~strcmp(dcm,'DICM')
    % Try truncated DICOM file fomat
    fseek(FID,0,'bof');
    Tag.Group   = fread(FID,1,'ushort');
    Tag.Element = fread(FID,1,'ushort');
    if isempty(Tag.Group) || isempty(Tag.Element)
        fclose(FID);
        warning('spm:dicom','%s: Truncated file.', DicomFilename);
        return;
    end
    if isempty(find(DicomDictionary.group==Tag.Group & DicomDictionary.element==Tag.Element,1)) && ~(Tag.Group==8 && Tag.Element==0)
        % Entry not found in DICOM dictionary and not from a GE Twin+excite
        % that starts with with an 8/0 Tag that I can't find any
        % documentation for.
        fclose(FID);
        warning('spm:dicom','%s: Not a DICOM file.', DicomFilename);
        return;
    else
        fseek(FID,0,'bof');
    end
end

try
    Header = ReadDicom(FID, 'il', DicomDictionary, Options);
    if isempty(Header)
        fclose(FID);
        return;
    end
    Header.Filename = fopen(FID);
catch
    le = lasterror;
    fprintf('%s: Trouble reading DICOM file (%s), skipping.\n', fopen(FID), le.message);
end

fclose(FID);

%==========================================================================
% function [Header, BytesRead] = ReadDicom(FID, TransferSyntax, DicomDictionary, Options, NumBytes)
%==========================================================================
function [Header, BytesRead] = ReadDicom(FID, TransferSyntax, DicomDictionary, Options, NumBytes)
if nargin<5, NumBytes = 4294967295; end % FFFFFFFF
BytesRead = 0;
Header    = [];
while BytesRead < NumBytes
    Tag = ReadTag(FID, TransferSyntax, DicomDictionary, Options);
    if isempty(Tag), break; end

    %fprintf('(%.4X,%.4X) "%s" %d %d %s\n', Tag.Group, Tag.Element, Tag.VR, Tag.Length, Tag.ExtraBytes, Tag.Name);

    if Tag.Group==65534 && Tag.Element==57357 % FFFE,E00D ItemDelimitationItem
        break;
    end
    if Tag.Length>0
        % Handle particular fields as special cases, otherwise just use default treatment
        switch Tag.Name
            case {'GroupLength'}
                % Ignore it
                fseek(FID,Tag.Length,'cof');
            case {'PixelData'}
                Header.StartOfPixelData = ftell(FID);
                Header.SizeOfPixelData  = Tag.Length;
                Header.VROfPixelData    = Tag.VR;
                fseek(FID,Tag.Length,'cof');
            case {'CSAData'} % raw data
                Header.StartOfCSAData = ftell(FID);
                Header.SizeOfCSAData  = Tag.Length;
                fseek(FID,Tag.Length,'cof');
            case {'CSAImageHeaderInfo', 'CSASeriesHeaderInfo','CSANonImageHeaderInfoVA',...
                  'CSAMiscProtocolHeaderInfoVA','CSANonImageHeaderInfoVB','CSAMiscProtocolHeaderInfoVB'}
                Content  = DecodeCSA(FID,Tag.Length);
                Header.(Tag.Name) = Content;
            case {'TransferSyntaxUID'}
                Content = deblank(char(fread(FID,Tag.Length,'uint8')'));
                Header.(Tag.Name) = Content;
                switch Content
                    case {'1.2.840.10008.1.2'}      % Implicit VR Little Endian
                        TransferSyntax = 'il';
                    case {'1.2.840.10008.1.2.1'}    % Explicit VR Little Endian
                        TransferSyntax = 'el';
                    case {'1.2.840.10008.1.2.1.99'} % Deflated Explicit VR Little Endian
                        warning('spm:dicom','%s: Cant read Deflated Explicit VR Little Endian file.', fopen(FID));
                       %TransferSyntax = 'dl';
                        return;
                    case {'1.2.840.10008.1.2.2'}    % Explicit VR Big Endian
                        %warning('spm:dicom','%s: Cant read Explicit VR Big Endian file',fopen(FID));
                        TransferSyntax = 'eb'; % Unused
                    case {'1.2.840.10008.1.2.4.50','1.2.840.10008.1.2.4.51','1.2.840.10008.1.2.4.70',...
                          '1.2.840.10008.1.2.4.80','1.2.840.10008.1.2.4.90','1.2.840.10008.1.2.4.91'} % JPEG Explicit VR
                        TransferSyntax = 'el';
                        %warning('spm:dicom',['Cant read JPEG Encoded file "' fopen(FID) '".']);
                    otherwise
                       %TransferSyntax = 'el';
                        warning('spm:dicom','%s: Unknown Transfer Syntax UID (%s).',fopen(FID), Content);
                        return;
                end
            otherwise
                % Default treatment
                switch Tag.VR
                    case {'UN'}
                        % Unknown - read as char
                        Content = fread(FID,Tag.Length,'uint8')';
                    case {'AE', 'AS', 'CS', 'DA', 'DS', 'DT', 'IS', 'LO', 'LT',...
                            'PN', 'SH', 'ST', 'TM', 'UI', 'UT'}
                        % Character strings
                        Content = char(fread(FID,Tag.Length,'uint8')');

                        switch Tag.VR
                            case {'UI','ST'}
                                Content = deblank(Content);
                            case {'DS'}
                                try
                                    Content = textscan(Content,'%f','delimiter','\\')';
                                    Content = Content{1};
                                catch
                                    Content = textscan(Content,'%f','delimiter','/')';
                                    Content = Content{1};
                                end
                            case {'IS'}
                                Content = textscan(Content,'%d','delimiter','\\')';
                                Content = double(Content{1});
                            case {'DA'}
                                Content = strrep(Content,'.',' ');
                                Content = textscan(Content,'%4d%2d%2d');
                                [y,m,d] = deal(Content{:});
                                Content = datenum(double(y),double(m),double(d));
                            case {'TM'}
                                if any(Content==':')
                                    Content = textscan(Content,'%d:%d:%f');
                                    [h,m,s] = deal(Content{:});
                                    h       = double(h);
                                    m       = double(m);
                                else
                                    Content = textscan(Content,'%2d%2d%f');
                                    [h,m,s] = deal(Content{:});
                                    h       = double(h);
                                    m       = double(m);
                                end
                                if isempty(h), h = 0; end
                                if isempty(m), m = 0; end
                                if isempty(s), s = 0; end
                                Content = s+60*(m+60*h);
                            case {'LO'}
                                Content = strrep(Content,'+AF8-','_');
                            otherwise
                        end
                    case {'OB'}
                        % dont know if this should be signed or unsigned
                        Content = fread(FID,Tag.Length,'uint8')';
                    case {'US', 'AT', 'OW'}
                        Content = fread(FID,Tag.Length/2,'uint16')';
                    case {'SS'}
                        Content = fread(FID,Tag.Length/2,'int16')';
                    case {'UL'}
                        Content = fread(FID,Tag.Length/4,'uint32')';
                    case {'SL'}
                        Content = fread(FID,Tag.Length/4,'int32')';
                    case {'FL','OF'}
                        Content = fread(FID,Tag.Length/4,'float')';
                    case {'FD','OD'}
                        Content = fread(FID,Tag.Length/8,'double')';
                    case {'SQ'}
                        [Content, BytesJustRead] = ReadSQ(FID, TransferSyntax, DicomDictionary, Options, Tag.Length);
                        if BytesJustRead==-1 && isempty(Content)
                            Header    = [];
                            BytesRead = -1;
                            return;
                        end
                        Tag.Length = BytesJustRead;
                    otherwise
                        Content = fread(FID,Tag.Length,'uint8')';
                end
                if ~isempty(Tag.Name)
                    if isa(Options.abort, 'function_handle') && feval(Options.abort, Tag.Name, Content)
                        Header    = [];
                        BytesRead = -1;
                        return;
                    end
                    Header.(Tag.Name) = Content;
                end
        end
    end
    BytesRead = BytesRead + Tag.ExtraBytes + Tag.Length;
end

%==========================================================================
% function [Header, BytesRead] = ReadSQ(FID, TransferSyntax, DicomDictionary, Options, NumBytes)
%==========================================================================
function [Header, BytesRead] = ReadSQ(FID, TransferSyntax, DicomDictionary, Options, NumBytes)
Header = {};
n      = 0;
BytesRead = 0;
while BytesRead<NumBytes
    Tag.Group   = fread(FID,1,'ushort');
    Tag.Element = fread(FID,1,'ushort');
    Tag.Length  = fread(FID,1,'uint');
    if isempty(Tag.Length), return; end % End of file
    BytesRead   = BytesRead + 8;
    if (Tag.Group == 65534) && (Tag.Element == 57344)   % FFFE/E000 Item
        [Item, BytesJustRead] = ReadDicom(FID, TransferSyntax, DicomDictionary, Options, Tag.Length);
        if BytesJustRead==-1 && isempty(Item)
            Header    = [];
            BytesRead = -1;
            return;
        end

        BytesRead    = BytesRead + BytesJustRead;
        if ~isempty(Item)
            n         = n + 1;
            Header{n} = Item;
        end
    elseif (Tag.Group == 65279) && (Tag.Element == 224) % FEFF/00E0 Item (Byte-swapped)
        % Byte-swapped
        [Filename,Permission,MachineFormat]      = fopen(FID);
        TransferSyntaxSwapped = TransferSyntax;
        if TransferSyntax(2)=='b'
            TransferSyntaxSwapped(2) = 'l';
        else
            TransferSyntaxSwapped(2) = 'b';
        end
        [Item, BytesJustRead] = ReadDicom(FID, TransferSyntaxSwapped, DicomDictionary, Options, Tag.Length);
        if BytesJustRead==-1 && isempty(Item)
            BytesRead = -1;
            return;
        end

        BytesRead = BytesRead + BytesJustRead;
        n         = n + 1;
        Header{n} = Item;
        Position  = ftell(FID);
        fclose(FID);
        FID     = fopen(Filename, Permission, MachineFormat);
        fseek(FID, Position, 'bof');
    elseif (Tag.Group == 65534) && (Tag.Element == 57565) % FFFE/E0DD SequenceDelimitationItem
        break;
    elseif (Tag.Group == 65279) && (Tag.Element == 56800) % FEFF/DDE0 SequenceDelimitationItem (Byte-swapped)
        % Byte-swapped
        break;
    else
        warning('spm:dicom','%s: Tag (%.4X,%.4X) is unexpected in sequence.', fopen(FID), Tag.Group, Tag.Element);
    end
end

%==========================================================================
% function Tag = ReadTag(FID,TransferSyntax,DicomDictionary, Options)
%==========================================================================
function Tag = ReadTag(FID, TransferSyntax, DicomDictionary, Options)
[GroupElement, Bytesread] = fread(FID,2,'ushort');
if Bytesread < 2, Tag = []; return; end
Tag.Group   = GroupElement(1);
Tag.Element = GroupElement(2);
if Tag.Group == 2, TransferSyntax = 'el'; end
t           = find(DicomDictionary.group==Tag.Group & DicomDictionary.element==Tag.Element);
if ~isempty(t)
    Tag.Name = DicomDictionary.values(t).name;
    Tag.VR   = DicomDictionary.values(t).vr{1};
else
    if ~Options.all_fields
        % With a reduced DicomDictionary, this can speed things up considerably.
        Tag.Name = '';
    else
        if Tag.Element~=0
            % Odd Group numbers indicate private (manufacturer-specific) fields
            if rem(Tag.Group,2)
                % Shadow Group
                Tag.Name = sprintf('Private_%.4x_%.4x',Tag.Group,Tag.Element);
            else
                Tag.Name = sprintf('Tag_%.4x_%.4x',Tag.Group,Tag.Element);
            end
            Tag.VR   = 'UN';
        else
            Tag.Name = '';
            Tag.VR   = 'UN';
        end
    end
end

if TransferSyntax(2) == 'b'
    [Filename,Permission,MachineFormat] = fopen(FID);
    if strcmp(MachineFormat,'ieee-le') || strcmp(MachineFormat,'ieee-le.l64')
        Position = ftell(FID);
        fclose(FID);
        FID      = fopen(Filename,Permission,'ieee-be');
        fseek(FID,Position,'bof');
    end
end

Tag.ExtraBytes = 4; % Two shorts read so far

if TransferSyntax(1) =='e'
    Tag.VR         = char(fread(FID,2,'uint8')');
    Tag.ExtraBytes = Tag.ExtraBytes + 2;
    switch Tag.VR
        case {'OB','OW','OF','OD','SQ','UN','UT'}
            if ~strcmp(Tag.VR,'UN') || Tag.Group~=65534
                fread(FID,1,'ushort');
                Tag.ExtraBytes = Tag.ExtraBytes + 2;
            else
                warning('spm:dicom','%s: Possible problem with %s Tag (VR="%s").', fopen(FID), Tag.Name, Tag.VR);
            end
            Tag.Length     = double(fread(FID,1,'uint'));
            Tag.ExtraBytes = Tag.ExtraBytes + 4;

        case {'AE','AS','AT','CS','DA','DS','DT','FD','FL','IS','LO','LT','PN','SH','SL','SS','ST','TM','UI','UL','US'}
            Tag.Length     = double(fread(FID,1,'ushort'));
            Tag.ExtraBytes = Tag.ExtraBytes + 2;

        case char([0 0])
            if (Tag.Group == 65534) && (Tag.Element == 57357)    % ItemDeliminationItem
                % at least on GE, ItemDeliminationItem does not have a
                % VR, but 4 bytes zeroes as Length
                Tag.VR         = 'UN';
                fread(FID,1,'ushort'); % Should be zero
                Tag.Length     = 0;
                Tag.ExtraBytes = Tag.ExtraBytes + 2;

            elseif (Tag.Group == 65534) && (Tag.Element == 57565) % SequenceDelimitationItem
                Tag.VR         = 'UN';
                Tag.Length     = 0;
                fread(FID,1,'ushort'); % Should be zero
                Tag.ExtraBytes = Tag.ExtraBytes + 2;

            else
                warning('spm:dicom','%s: Don''t know how to handle VR of "\0\0" in %s.',fopen(FID),Tag.Name);
            end
        otherwise
            warning('spm:dicom','%s: Possible problem with %s Tag (VR="%s")\n', fopen(FID), Tag.Name, Tag.VR);
            fread(FID,1,'ushort');
            Tag.Length     = double(fread(FID,1,'uint'));
            Tag.ExtraBytes = Tag.ExtraBytes + 2+4;
    end
else
    Tag.Length     = double(fread(FID,1,'uint'));
    Tag.ExtraBytes = Tag.ExtraBytes + 4;
end

if isempty(Tag.VR) || isempty(Tag.Length)
    Tag = [];
    return;
end

if rem(Tag.Length,2)
    if Tag.Length==4294967295 % FFFFFFFF
        if strcmp(Tag.VR,'UN')
            warning('spm:dicom','%s: Unknown Value Representation of %s Tag, assuming ''SQ''.', fopen(FID), Tag.Name);
            Tag.VR = 'SQ';
        end
        return;
    else
        warning('spm:dicom','%s: Odd numbered Value Length in %s Tag (%X).', fopen(FID), Tag.Name, Tag.Length);
        Tag = [];
    end
end

%==========================================================================
% function t = DecodeCSA(FID,NumBytes)
%==========================================================================
function t = DecodeCSA(FID,NumBytes)
% Decode shadow information (0029,1010) and (0029,1020) from Siemens files
[Filename,Permission,MachineFormat] = fopen(FID);
Position = ftell(FID);
if strcmp(MachineFormat,'ieee-be') || strcmp(MachineFormat,'ieee-be.l64')
    fclose(FID);
    FID  = fopen(Filename,Permission,'ieee-le');
    fseek(FID,Position,'bof');
end

c   = fread(FID,4,'uint8');
fseek(FID,Position,'bof');

if all(c'==[83 86 49 48]) % "SV10"
    t = DecodeCSA2(FID,NumBytes);
else
    t = DecodeCSA1(FID,NumBytes);
end

if strcmp(MachineFormat,'ieee-be') || strcmp(MachineFormat,'ieee-be.l64')
    fclose(FID);
    FID  = fopen(Filename, Permission, MachineFormat);
end
fseek(FID, Position+NumBytes, 'bof');

%==========================================================================
% function t = DecodeCSA1(FID,NumBytes)
%==========================================================================
function t = DecodeCSA1(FID,NumBytes)
n   = fread(FID,1,'uint32');
if isempty(n) || n>1024 || n <= 0
    fseek(FID,NumBytes-4,'cof');
    t = struct('name','JUNK: Don''t know how to read this damned file format');
    return;
end
fread(FID,1,'uint32'); % Unused "M" or 77 for some reason
BytesRead = 2*4;
t(n)      = struct('name','', 'vm',[], 'vr','', 'syngodt',[], 'nitems',[], 'xx',[], 'item',struct('xx',{},'val',{}));
for i=1:n
    t(i).name = fread(FID,64,'uint8')';
    msk       = find(~t(i).name)-1;
    if ~isempty(msk)
        t(i).name = char(t(i).name(1:msk(1)));
    else
        t(i).name = char(t(i).name);
    end
    t(i).vm      = fread(FID,1,'int32')';
    t(i).vr      = fread(FID,4,'uint8')';
    t(i).vr      = char(t(i).vr(1:3));
    t(i).syngodt = fread(FID,1,'int32')';
    t(i).nitems  = fread(FID,1,'int32')';
    t(i).xx      = fread(FID,1,'int32')'; % 77 or 205
    BytesRead    = BytesRead + 64+4+4+4+4+4;
    if t(i).nitems > 0
        t(i).item(t(i).nitems) = struct('xx',[], 'val',[]);
    end
    for j=1:t(i).nitems
        % This bit is just wierd
        t(i).item(j).xx  = fread(FID,4,'int32')'; % [x x 77 x]
        BytesToRead      = t(i).item(j).xx(1)-t(1).nitems;
        if BytesToRead<0 || BytesToRead+Bytesread+4*4>NumBytes
            t(i).item(j).val = '';
            t(i).item        = t(i).item(1:j);
            BytesRead        = BytesRead + 4*4;
            break;
        end
        t(i).item(j).val = char(fread(FID,BytesToRead,'uint8')');
        fread(FID,4-rem(BytesToRead,4),'uint8');
        BytesRead        = BytesRead + 4*4+BytesToRead+(4-rem(BytesToRead,4));
    end
end

%==========================================================================
% function t = DecodeCSA2(FID, NumBytes)
%==========================================================================
function t = DecodeCSA2(FID, NumBytes)
fread(FID, 8, 'uint8'); % Unused
n        = fread(FID, 1, 'uint32');
if isempty(n) || n>1024 || n < 0
    fseek(FID, NumBytes-4, 'cof');
    t = struct('name','Don''t know how to read this damned file format');
    return;
end
fread(FID, 1, 'uint32'); % Unused "M" or 77 for some reason
Position = 16;
t(n)     = struct('name','', 'vm',[], 'vr','', 'syngodt',[], 'nitems',[], 'xx',[], 'item',struct('xx',{},'val',{}));
for i=1:n
    t(i).name = fread(FID, 64, 'uint8')';
    Position  = Position + 64;
    msk       = find(~t(i).name)-1;
    if ~isempty(msk)
        t(i).name = char(t(i).name(1:msk(1)));
    else
        t(i).name = char(t(i).name);
    end
    t(i).vm      = fread(FID, 1, 'int32')';
    t(i).vr      = fread(FID, 4, 'uint8')';
    t(i).vr      = char(t(i).vr(1:3));
    t(i).syngodt = fread(FID, 1, 'int32')';
    t(i).nitems  = fread(FID, 1, 'int32')';
    t(i).xx      = fread(FID, 1, 'int32')'; % 77 or 205
    Position     = Position + 20;
    if t(i).nitems > 0
        t(i).item(t(i).nitems) = struct('xx',[], 'val',[]);
    end
    for j=1:t(i).nitems
        t(i).item(j).xx  = fread(FID,4,'int32')'; % [x x 77 x]
        Position         = Position + 16;
        BytesToRead      = t(i).item(j).xx(2);
        if BytesToRead > NumBytes-Position
            BytesToRead      = NumBytes-Position;
            t(i).item(j).val = char(fread(FID, BytesToRead, 'uint8')');
            fread(FID,rem(4-rem(BytesToRead, 4), 4), 'uint8');
            t(i).item        = t(i).item(1:j);
            warning('spm:dicom','%s: Problem reading Siemens CSA field.', fopen(FID));
            return;
        end
        t(i).item(j).val = char(fread(FID, BytesToRead, 'uint8')');
        fread(FID,rem(4-rem(BytesToRead, 4), 4), 'uint8');
    end
end


text:
function [ r, w ] = oned_gauss ( rule )

%*****************************************************************************80
%
%% ONED_GAUSS sets up a Gauss quadrature rule in [-1,+1].
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 March 2011
%
%  Author:
%
%    Jeff Borggaard
%
%  Parameters:
%
%    Input, integer RULE, the number of Gauss points desired.
%
%    Output, real R(RULE), Gauss points located between (-1,1).
%
%    Output, real W(RULE), Gauss weights corresponding to r.
%
  r = zeros(rule,1);
  w = zeros(rule,1);

  if ( rule == 1 )
    r(1) = 0;
    w(1) = 2;
  elseif ( rule == 2 )
    r(1) =-1.0 / sqrt(3.0);
    r(2) =-r(1);
    w(1) = 1.0;
    w(2) = 1.0;
  elseif ( rule == 3 )
    r(1) =-sqrt(3.0/5.0);
    r(2) = 0.0;
    r(3) =-r(1);
    w(1) = 5.0 / 9.0;
    w(2) = 8.0 / 9.0;
    w(3) = w(1);
  elseif ( rule == 4 )
    r(1) =-sqrt((3.0+2.0*sqrt(6.0/5.0))/7.0);
    r(2) =-sqrt((3.0-2.0*sqrt(6.0/5.0))/7.0);
    r(3) =-r(2);
    r(4) =-r(1);
    w(1) = 0.5 - 1.0 / ( 6.0 * sqrt(6.0/5.0) );
    w(2) = 0.5 + 1.0 / ( 6.0 * sqrt(6.0/5.0) );
    w(3) = w(2);
    w(4) = w(1);
  elseif ( rule == 5 )
    r(1) =-sqrt(5.0+4.0*sqrt(5.0/14.0)) / 3.0;
    r(2) =-sqrt(5.0-4.0*sqrt(5.0/14.0)) / 3.0;
    r(3) = 0.0;
    r(4) =-r(2);
    r(5) =-r(1);
    w(1) = 161.0/450.0-13.0/(180.*sqrt(5.0/14.0));
    w(2) = 161.0/450.0+13.0/(180.*sqrt(5.0/14.0));
    w(3) = 128.0/225.0;
    w(4) = w(2);
    w(5) = w(1);
  elseif ( rule == 6 )
    r(1) = -0.2386191861;
    r(2) = -0.6612093865;
    r(3) = -0.9324695142;
    r(4) = - r(1);
    r(5) = - r(2);
    r(6) = - r(3);
    w(1) = .4679139346;
    w(2) = .3607615730;
    w(3) = .1713244924;
    w(4) = w(1);
    w(5) = w(2);
    w(6) = w(3);
  else
    error('Quadrature rule not supported')
    keyboard
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function Q = S_object(prediction,GT)
% S_object Computes the object similarity between foreground maps and ground
% truth(as proposed in "Structure-measure:A new way to evaluate foreground
% maps" [Deng-Ping Fan et. al - ICCV 2017])
% Usage:
%   Q = S_object(prediction,GT)
% Input:
%   prediction - Binary/Non binary foreground map with values in the range
%                [0 1]. Type: double.
%   GT - Binary ground truth. Type: logical.
% Output:
%   Q - The object similarity score
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% compute the similarity of the foreground in the object level
prediction_fg = prediction;
prediction_fg(~GT)=0;
O_FG = Object(prediction_fg,GT);

% compute the similarity of the background
prediction_bg = 1.0 - prediction;
prediction_bg(GT) = 0;
O_BG = Object(prediction_bg,~GT);

% combine the foreground measure and background measure together
u = mean2(GT);
Q = u * O_FG + (1 - u) * O_BG;

end

function score = Object(prediction,GT)

% check the input
if isempty(prediction)
    score = 0;
    return;
end
if isinteger(prediction)
    prediction = double(prediction);
end
if (~isa( prediction, 'double' ))
    error('prediction should be of type: double');
end
if ((max(prediction(:))>1) || min(prediction(:))<0)
    error('prediction should be in the range of [0 1]');
end
if(~islogical(GT))
    error('GT should be of type: logical');
end

% compute the mean of the foreground or background in prediction
x = mean2(prediction(GT));

% compute the standard deviations of the foreground or background in prediction
sigma_x = std(prediction(GT));

score = 2.0 * x./(x^2 + 1.0 + sigma_x + eps);
end

author: ArcherFMY


text:
classdef ConvergencePlotComputerMaterialDesignExperimentPoisson < handle

    properties (Access = public)

    end

    properties (Access = private)
        filesPath
        testNames
        linesData
        barData
        fieldsData
        outPutPlotPath
        alphas
    end

    properties (Access = private)

    end

    methods (Access = public)

        function obj = ConvergencePlotComputerMaterialDesignExperimentPoisson()
            obj.init();
            obj.loadFieldsData();
            obj.plotData();
        end

    end

    methods (Access = private)

        function init(obj)
            obj.filesPath = '/media/alex/MyPassport/MaterialDesign/CStar/';
            obj.outPutPlotPath = '/home/alex/Dropbox/MaterialDesign/CC/Poisson/';
            %obj.linesData = [1:2,5:10];
            %obj.barData = [3:4];
            obj.linesData = [1:6,9:15];
            obj.barData = [7:8];
            obj.testNames = {'NegPoissonNoPerimeter5x6/'};
        end

        function loadFieldsData(obj)
            for iCase = 1:numel(obj.testNames)
                testPath = fullfile(obj.filesPath,obj.testNames{iCase},'/');
                s.testPath  = testPath;
                s.linesData = obj.linesData;
                s.barData   = obj.barData;
                m = MonitoringDataLoader(s);
                obj.fieldsData{iCase} = m.obtainData();
            end
        end

        function plotData(obj)
            p{1} = obj.plotCost();
            p{2} = obj.plotVolum();
            f = figure(1);
            legend({'$||\bf{C}(\rho) - {C}^*||_2$','$\textrm{Vol}(\rho)$'},'Interpreter','latex','Location','Best');
            p = plotPrinter(f,p);
            p.print(fullfile(obj.outPutPlotPath,'CostAndVolume'))
        end

        function p = plotCost(obj)
            for iCase = 1:numel(obj.testNames)
                fieldData = obj.fieldsData{iCase};
                fieldToPlot = 'C - C not scaled';
                [x,y] = obj.obtainField(fieldToPlot,fieldData);
                p = semilogy(x,y);
                hold on
            end
        end

        function p = plotVolum(obj)
            figure(1)
            hold on
            yyaxis right
            for iCase = 1:numel(obj.testNames)
                fieldData = obj.fieldsData{iCase};
                fieldToPlot = 'Volum';
                [x,y] = obj.obtainField(fieldToPlot,fieldData);
                p = plot(x,y);
                ylim([0,1])
                hold on
            end
        end

    end

    methods (Access = private, Static)

        function [xV,yV] = obtainField(fieldName,fieldData)
            for iField = 1:numel(fieldData)
                title = fieldData{iField}.title;
                if strcmp(fieldName,title)
                    xV = fieldData{iField}.xValue;
                    yV = fieldData{iField}.yValue;
                end
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [gIX,B] = SortGroupIXbyScore(H,gIX,numU,descend) %#ok<INUSD>
% for ranking functions % new gIX is sorted based on H, size(H)=[numU,1];
if exist('descend','var'),
    [B,I] = sort(H,'descend');
else
    [B,I] = sort(H);
end

if ~exist('numU','var'),
    numU = length(unique(gIX));
end

gIX_last = gIX;
for i = 1:numU,
    gIX(gIX_last==I(i)) = i;
end
end

author: xiuyechen

repo: FishExplorer


text:
filename='Gripping_triangle_coarse';
ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'nonadjoint_compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
constraint_case = 'INEQUALITY';
optimizer = 'MMA'; incrementFactor = 1;
filterType = 'P1';

nsteps = 1;
Vfrac_final = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;


text:
function [Cn, PNR] = correlation_pnr_parallel(obj, frame_range)
%% compute correlation image and pnr imaging in parallel
%% input:
%   frame_range: 1 X 2 vector indicating the starting and ending frames

%% Output:
%   Cn:     correlation image
%   PNR:    peak to noise ratio
%% Author: Pengcheng Zhou, Columbia University, 2017
%% email: zhoupc1988@gmail.com

%% process parameters

% map data
mat_data = obj.P.mat_data;

% dimension of data
dims = mat_data.dims;
d1 = dims(1);
d2 = dims(2);
T = dims(3);
obj.options.d1 = d1;
obj.options.d2 = d2;

% parameters for patching information
patch_pos = mat_data.patch_pos;
block_pos = mat_data.block_pos;

% number of patches
[nr_patch, nc_patch] = size(patch_pos);

% downsampling
ssub = obj.options.ssub;
tsub = obj.options.tsub;

% frames to be loaded for initialization
if ~exist('frame_range', 'var')
    frame_range = obj.frame_range;
end
if isempty(frame_range)
    frame_range = [1, T];
else
    frame_range(frame_range<1) = 1;
    frame_range(frame_range>T) = T;
end
T = diff(frame_range) + 1;
obj.frame_range = frame_range;

% parameters for detrending
if isfield(obj.options, 'nk') % number of knots for creating spline basis
    nk = obj.options.nk;
else
    nk = 1;
end
detrend_method = obj.options.detrend_method;

% parameter for avoiding using boundaries pixels as seed pixels
options = obj.options;

%% compute correlation image and pnr image in parallel
Cn = zeros(d1, d2);
PNR = zeros(d1, d2);

results = cell(nr_patch, nc_patch);

for m=1:(nr_patch*nr_patch)
    fprintf('|');
end
fprintf('\n');
parfor mpatch=1:(nr_patch*nc_patch)
    % get the indices corresponding to the selected patch
    tmp_patch = patch_pos{mpatch};
    tmp_block = block_pos{mpatch};

    % boundaries pixels to be avoided for detecting seed pixels
    tmp_options = options;

    % patch dimension
    tmp_d1 = diff(tmp_block(1:2))+1;
    tmp_d2 = diff(tmp_block(3:4))+1;
    tmp_options.gSig = tmp_options.gSig/ssub;
    tmp_options.gSiz = ceil(tmp_options.gSiz/ssub);

    % load the patch data
    Ypatch = get_patch_data(mat_data, tmp_patch, frame_range, true);
    if (ssub~=1)||(tsub~=1)
        Ypatch = dsData(Ypatch, tmp_options);
    end
    [tmp_options.d1, tmp_options.d2, T] = size(Ypatch);
    Ypatch = double(reshape(Ypatch, [], T));
    if nk>1
        Ypatch_dt = detrend_data(Ypatch, nk, detrend_method); % detrend data
        [tmp_Cn, tmp_PNR] = correlation_image_endoscope(Ypatch_dt, tmp_options);
    else
        [tmp_Cn, tmp_PNR] = correlation_image_endoscope(Ypatch, tmp_options);
    end
    if (ssub~=1)
        tmp_Cn = imresize(tmp_Cn, [tmp_d1, tmp_d2]);
        tmp_PNR = imresize(tmp_PNR, [tmp_d1, tmp_d2]);
    end
    % put everthing into one struct variable
    results{mpatch} = struct('Cn', tmp_Cn, 'PNR', tmp_PNR);
    fprintf('.');
end
fprintf('\n');

%% collect results
for mpatch=1:(nr_patch*nc_patch)
    % get the indices corresponding to the selected patch
    [mr, mc] = ind2sub([nr_patch, nc_patch], mpatch);
    tmp_patch = patch_pos{mr, mc};
    tmp_block = block_pos{mr, mc};
    r0 = tmp_block(1);
    r1 = tmp_block(2);
    c0 = tmp_block(3);
    c1 = tmp_block(4);
    ind_patch = true(r1-r0+1, c1-c0+1);
    ind_patch((tmp_patch(1):tmp_patch(2))-r0+1, (tmp_patch(3):tmp_patch(4))-c0+1) = false;

    % unpack results
    tmp_results = results{mpatch};
    %     eval(sprintf('tmp_results=results_patch_%d;', mpatch));
    tmp_Cn = tmp_results.Cn;
    tmp_PNR = tmp_results.PNR;

    Cn(r0:r1, c0:c1) = max(Cn(r0:r1, c0:c1), tmp_Cn.*(1-ind_patch));
    PNR(r0:r1, c0:c1) = max(PNR(r0:r1, c0:c1), tmp_PNR.*(1-ind_patch));
end

author: zhoupc


text:
% VAR - variance of memory mapped underlying array
%
% Author: Arnaud Delorme, SCCN, INC, UCSD, Nov. 2008

% Copyright (C) 2008 Arnaud Delorme, SCCN, INC, UCSD
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function sumval = var(obj,flag,dim)

    if nargin < 2
        flag = 0;
    end
    if nargin < 3
        dim = 1;
    end

    meanvalsq = mean(obj,dim).^2;

    sumval = 0;
    s1  = size(obj);
    ss.type = '()';
    ss.subs(1:length(s1)) = { ':' };
    for index = 1:s1(dim)
        ss.subs{dim} = index;
        tmpdata = subsref(obj, ss);
        sumval = sumval + tmpdata.*tmpdata - meanvalsq;
    end
    if isempty(flag) || flag == 0
         sumval = sumval/(size(obj,dim)-1);
    else sumval = sumval/size(obj,dim);
    end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [um] = au2um(au)
% Convert length from astronomical units to microns.
% Chad A. Greene 2012
um = au*1.49597870691e+17;

author: Sable

repo: mcbench-benchmarks


text:
%HOUGHCIRCLES  Finds circles in a grayscale image using the Hough transform
%
%     circles = cv.HoughCircles(image)
%     circles = cv.HoughCircles(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __image__ 8-bit, single-channel, grayscale input image.
%
% ## Output
% * __circles__ Output vector of found circles. A cell-array of 3-element
%   floating-point vectors `{[x, y, radius], ...}`.
%
% ## Options
% * __Method__ Detection method. Currently, the only implemented method is
%   'Gradient' (default). One of:
%   * __Standard__ classical or standard Hough transform. Every line is
%     represented by two floating-point numbers `(rho,theta)`, where `rho` is
%     a distance between `(0,0)` point and the line, and `theta` is the angle
%     between x-axis and the normal to the line. Thus, the matrix must be (the
%     created sequence will be) of `single` type with 2-channels.
%   * __Probabilistic__ probabilistic Hough transform (more efficient in case
%     if the picture contains a few long linear segments). It returns line
%     segments rather than the whole line. Each segment is represented by
%     starting and ending points, and the matrix must be (the created sequence
%     will be) of the `int32` type with 4-channels.
%   * __MultiScale__ multi-scale variant of the classical Hough transform. The
%     lines are encoded the same way as 'Standard'.
%   * __Gradient__ basically 21HT, described in [Yuen90].
% * __DP__ Inverse ratio of the accumulator resolution to the image resolution.
%   For example, if `DP=1`, the accumulator has the same resolution as the
%   input image. If `DP=2`, the accumulator has half as big width and height.
%   default 1.
% * __MinDist__ Minimum distance between the centers of the detected circles.
%   If the parameter is too small, multiple neighbor circles may be falsely
%   detected in addition to a true one. If it is too large, some circles may
%   be missed. default `size(image,1)/8`.
% * __Param1__ First method-specific parameter. In case of 'Gradient', it is
%   the higher threshold of the two passed to the cv.Canny edge detector
%   (the lower one is twice smaller). default 100.
% * __Param2__ Second method-specific parameter. In case of 'Gradient', it is
%   the accumulator threshold for the circle centers at the detection stage.
%   The smaller it is, the more false circles may be detected. Circles,
%   corresponding to the larger accumulator values, will be returned first.
%   default 100.
% * __MinRadius__ Minimum circle radius. default 0.
% * __MaxRadius__ Maximum circle radius. If <= 0, uses the maximum image
%   dimension. If < 0, returns centers without finding the radius. default 0.
%
% The function finds circles in a grayscale image using a modification of the
% Hough transform.
%
% ### Note
% Usually the function detects the centers of circles well. However, it may
% fail to find correct radii. You can assist to the function by specifying the
% radius range (`MinRadius` and `MaxRadius`) if you know it. Or, you may set
% `MaxRadius` to a negative number to return centers only without radius
% search, and find the correct radius using an additional procedure.
%
% ## References
% [Yuen90]:
% > HK Yuen, John Princen, John Illingworth, and Josef Kittler.
% > "Comparative study of hough transform methods for circle finding".
% > Image and Vision Computing, 8(1):71-77, 1990.
%
% See also: cv.fitEllipse, cv.minEnclosingCircle, imfindcircles
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function cos_power_int_values_test ( )

%*****************************************************************************80
%
%% COS_POWER_INT_VALUES_TEST demonstrates the use of COS_POWER_INT_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 January 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'COS_POWER_INT_VALUES_TEST:\n' );
  fprintf ( 1, '  COS_POWER_INT_VALUES returns values of \n' );
  fprintf ( 1, '  the N-th power of the cosine function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      A      B      N      FX\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, a, b, n, fx ] = cos_power_int_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %8f  %8f  %6d  %24.16f\n', a, b, n, fx );

  end

  return
end


text:
function [events] = ma_updateEventFcnWithNewBodyInfo(events, newBodyInfo, celBodyData)
%ma_updateEventFcnWithNewBodyInfo Summary of this function goes here
%   Detailed explanation goes here

    for(i=1:length(events)) %#ok<*NO4LP>
        event = events{i};

        if(~isempty(strfind(func2str(event), 'ma_SoITransitionEvents')))
            event = @(T,Y) ma_SoITransitionEvents(T,Y, newBodyInfo, celBodyData);

        elseif(~isempty(strfind(func2str(event), 'ma_TrueAnomalyEvent')))
            [truTarget,~,refBody,~] = event('getInfo',[]);
            event = @(T,Y) ma_TrueAnomalyEvent(T,Y, truTarget, newBodyInfo, refBody);

        elseif(~isempty(strfind(func2str(event), 'ma_NodeEvents')))
            [node,~,refBody,~] = event('getInfo',[]);
            event = @(T,Y) ma_NodeEvents(T,Y, node, newBodyInfo, refBody);

        else
            error(['Unknown event function handle: ', func2str(event)]);
        end

        events{i} = event;
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function triangle_wandzura_rule_test02 ( )

%*****************************************************************************80
%
%% TRIANGLE_WANDZURA_RULE_TEST02 tests WANDZURA_RULE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 December 2006
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGLE_WANDZURA_RULE_TEST02\n' );
  fprintf ( 1, '  WANDZURA_RULE returns the points and weights\n' );
  fprintf ( 1, '  of a Wandzura rule for the triangle.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  In this test, we simply check that the weights\n' );
  fprintf ( 1, '  sum to 1.\n' );

  rule_num = wandzura_rule_num ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of available rules = %d\n', rule_num );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      Rule    Sum of weights\n' );
  fprintf ( 1, '\n' );

  for rule = 1 : rule_num

    order_num = wandzura_order_num ( rule );

    [ xy, w ] = wandzura_rule ( rule, order_num );

    w_sum = sum ( w(1:order_num) );

    fprintf ( 1, '  %8d  %14f\n', rule, w_sum );

  end

  return
end


text:
% this code is revised based on ILSVRC 2013 (http://www.image-net.org/challenges/LSVRC/2013/)

function  top_recall  = top_recall_Relationship(Nre, tuple_confs_cell, tuple_labels_cell, sub_bboxes_cell, obj_bboxes_cell)

load('gt.mat','gt_tuple_label','gt_obj_bboxes','gt_sub_bboxes');

%num_imgs = length(gt_tuple_label);
num_imgs = 1000;
for i=1:num_imgs
    [tuple_confs_cell{i}, ind] = sort(tuple_confs_cell{i},'descend');
    if length(ind) >= Nre
        tuple_confs_cell{i} = tuple_confs_cell{i}(1:Nre);
        tuple_labels_cell{i} = tuple_labels_cell{i}(ind(1:Nre),:);
        obj_bboxes_cell{i} = obj_bboxes_cell{i}(ind(1:Nre),:);
        sub_bboxes_cell{i} = sub_bboxes_cell{i}(ind(1:Nre),:);
    else
        tuple_labels_cell{i} = tuple_labels_cell{i}(ind,:);
        obj_bboxes_cell{i} = obj_bboxes_cell{i}(ind,:);
        sub_bboxes_cell{i} = sub_bboxes_cell{i}(ind,:);
    end
end

num_pos_tuple = 0;
for ii = 1 : num_imgs
    num_pos_tuple = num_pos_tuple + size(gt_tuple_label{ii},1);
end

tp_cell = cell(1,num_imgs);
fp_cell = cell(1,num_imgs);

gt_thr = 0.5;
iterc = 0;

% iterate over images
for i=1:num_imgs

    gt_tupLabel = gt_tuple_label{i};
    gt_objBox = gt_obj_bboxes{i};
    gt_subBox = gt_sub_bboxes{i};

    num_gt_tuple = size(gt_tupLabel,1);
    gt_detected = zeros(1,num_gt_tuple);

    labels = tuple_labels_cell{i};
    boxObj = obj_bboxes_cell{i};
    boxSub = sub_bboxes_cell{i};

    num_obj = size(labels,1);
    tp = zeros(1,num_obj);
    fp = zeros(1,num_obj);
    for j=1:num_obj

        bbO = boxObj(j,:);
        bbS = boxSub(j,:);
        ovmax = -inf;
        kmax = -1;

        for k=1:num_gt_tuple
            if norm(labels(j,:) - gt_tupLabel(k,:),2) ~= 0
                continue;
            end
            if gt_detected(k) > 0
                continue;
            end

            bbgtO = gt_objBox(k,:);
            bbgtS = gt_subBox(k,:);

            biO=[max(bbO(1),bbgtO(1)) ; max(bbO(2),bbgtO(2)) ; min(bbO(3),bbgtO(3)) ; min(bbO(4),bbgtO(4))];
            iwO=biO(3)-biO(1)+1;
            ihO=biO(4)-biO(2)+1;

            biS=[max(bbS(1),bbgtS(1)) ; max(bbS(2),bbgtS(2)) ; min(bbS(3),bbgtS(3)) ; min(bbS(4),bbgtS(4))];
            iwS=biS(3)-biS(1)+1;
            ihS=biS(4)-biS(2)+1;

            if iwO>0 & ihO>0 &  iwS>0 & ihS>0
                % compute overlap as area of intersection / area of union
                uaO=(bbO(3)-bbO(1)+1)*(bbO(4)-bbO(2)+1)+...
                   (bbgtO(3)-bbgtO(1)+1)*(bbgtO(4)-bbgtO(2)+1)-...
                   iwO*ihO;
                ovO =iwO*ihO/uaO;

                uaS=(bbS(3)-bbS(1)+1)*(bbS(4)-bbS(2)+1)+...
                   (bbgtS(3)-bbgtS(1)+1)*(bbgtS(4)-bbgtS(2)+1)-...
                   iwS*ihS;
                ovS =iwS*ihS/uaS;
                ov = min([ovO,ovS]);

                % makes sure that this object is detected according
                % to its individual threshold
                if ov >= gt_thr && ov > ovmax
                    ovmax=ov;
                    kmax=k;
                end
            end
        end

        if kmax > 0
            tp(j) = 1;
            gt_detected(kmax) = 1;

        else
            fp(j) = 1;
        end
    end

    % put back into global vector
    tp_cell{i} = tp;
    fp_cell{i} = fp;

end
t = tic;
tp_all = [];
fp_all = [];
confs = [];
for ii = 1 : num_imgs
tp_all = [tp_all; tp_cell{ii}(:) ];
fp_all = [fp_all; fp_cell{ii}(:) ];
confs = [confs; tuple_confs_cell{ii}(:)];
end

[confs, ind] = sort(confs,'descend');
tp_all = tp_all(ind);
fp_all = fp_all(ind);

tp = cumsum(tp_all );
fp = cumsum(fp_all );
recall =(tp/num_pos_tuple);
top_recall = recall(end);

end

author: Prof-Lu-Cewu

repo: Visual-Relationship-Detection