text:
%This Matlab script can be used to reproduce Figure 3.4 in the monograph:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.0 (Last edited: 2017-11-04)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Empty workspace and close figures
close all;
clear;

%Define the range of BS antennas
Mvalues = [1 10 100];

%Angular standard deviation in the local scattering model (in degrees)
ASD = 10;

%Nominal angle of desired UE
thetaDesired = pi/6;

%Range of nominal angles of the interfering UE
varphiInterfererDegrees = -180:1:180;
varphiInterfererRadians = varphiInterfererDegrees*(pi/180);

%Define the antenna spacing (in number of wavelengths)
antennaSpacing = 1/2; %Half wavelength distance

%Define the effective SNR in (3.13) for the desired UE
SNR1dB = 10;
SNR1 = 10.^(SNR1dB/10);

%Define the effective SNR in (3.13) for the interfering UE
SNR2dB = 0;
SNR2 = 10.^(SNR2dB/10);

%Preallocate matrices for storing the simulation results
correlationcoeff = zeros(length(varphiInterfererRadians),length(thetaDesired),length(Mvalues));

%Compute the spatial correlation matrix of the desired UE
R1 = functionRlocalscattering(max(Mvalues),thetaDesired,ASD,antennaSpacing);

%% Go through all angles of interfering UE
for n = 1:length(varphiInterfererRadians)

    %Output simulation progress
    disp([num2str(n) ' angles out of ' num2str(length(varphiInterfererRadians))]);

    %Compute the spatial correlation matrix of the interfering UE
    R2 = functionRlocalscattering(max(Mvalues),varphiInterfererRadians(n),ASD,antennaSpacing);

    %Go through all number of antennas
    for m = 1:length(Mvalues)

        %Extract correlation matrices of the specified dimension
        R1m = R1(1:Mvalues(m),1:Mvalues(m));
        R2m = R2(1:Mvalues(m),1:Mvalues(m));

        %Compute the denominator in (3.18)
        normalization = sqrt(SNR1*SNR2*abs(trace(R1m*((SNR1*R1m+SNR2*R2m+eye(Mvalues(m)))\R1m)))*abs(trace(R2m*((SNR1*R1m+SNR2*R2m+eye(Mvalues(m)))\R2m))));

        %Compute absolute value of antenna-averaged correlation coefficient in (3.18)
        correlationcoeff(n,m) = sqrt(SNR1*SNR2)*abs(trace(R1m*((SNR1*R1m+SNR2*R2m+eye(Mvalues(m)))\R2m)))/normalization;

    end

end

%% Plot the simulation results
figure;
hold on; box on;

plot(varphiInterfererDegrees,correlationcoeff(:,1),'k-','LineWidth',1);
plot(varphiInterfererDegrees,correlationcoeff(:,2),'r--','LineWidth',1);
plot(varphiInterfererDegrees,correlationcoeff(:,3),'b-.','LineWidth',1);

xlabel('Angle of interfering UE [degree]');
ylabel('Antenna-averaged correlation coefficient');
xlim([-180 180]);
ylim([0 1.1]);

legend('M=1','M=10','M=100','Location','Best');

author: emilbjornson

repo: massivemimobook


text:
% Demo of use of TVDIP - total variation denoising algorithm
clear all;

% Generate test signal: piecewise smooth with noise corruption
y = [1.5*ones(10000,1); -0.7*ones(2000,1); 1.3*ones(100,1); -0.3*ones(5000,1)];
N = length(y);
y = y + 0.5*randn(N,1);

% Find the value of lambda greater than which the TVD solution is just the
% mean
lmax = tvdiplmax(y);

% Perform TV denoising for lambda across a range of values up to a small
% fraction of the maximum found above
lratio = [1e-4 1e-3 1e-2 1e-1];
L = length(lratio);
[x, E, status] = tvdip(y,lmax*lratio,1,1e-3,100);

% Plots - display original test signal y, and TVD results
close all;
figure;
for l = 1:L
    subplot(L,1,l);
    hold on;
    plot(y,'-','Color',0.8*[1 1 1]);
    plot(x(:,l),'k-');
    axis tight;
    title(sprintf('\\lambda=%5.2e, \\lambda/\\lambda_{max}=%5.2e',lmax*lratio(l),lratio(l)));
end
xlabel('n');
legend({'Input y_n','TVD x_n'});

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function r = ge(p,q)
% MEAS/GE  Implement p >= q for meas.

% find the difference between the two
diff = meas();
diff = p - q;

% see if p >= q at the 1-sigma level
if (diff.value > diff.error) | (abs(diff.value) < diff.error)
    r = 1;
else
    r = 0;
end

author: Sable

repo: mcbench-benchmarks


text:
function [fTmax, mRes] = calc_ChTime(mCatalog, fTimeWindow)
% [fTmax, mRes] = calc_ChTime(mCatalog, fTimeWindow)
% -------------------------------------------------------
% Function to identify times of maximum changes of seismicity in time using
% the absolute difference of the non-cumulative FMD in successive time windows
%
% Incoming:
% mCatalog    : EQ catalog
% fTimeWindow : Length of moving time window in dec. years
%
% Outgoing variable:
% fTmax : Time of maximum change
% mRes : Matrix containing time and absoulte changes
%
% J. Woessner
% last update: 21.01.2004

mRes = [];

% Initialize
fMinTime = min(mCatalog(:,3));
fMaxTime = max(mCatalog(:,3));
fTimePeriod = fMaxTime-fMinTime;

% Minimum bin
fMinBin = roundn(min(mCatalog(:,6)),-1);
fMaxBin = roundn(max(mCatalog(:,6)),-1);

% Seismicity of catalog normalized to time period
nNumevents = max(length(mCatalog(:,1)));
fNum = nNumevents/fTimePeriod;

% Starttime
fTime = fMinTime;

while fTime < fMaxTime-2*fTimeWindow
    vSel = (mCatalog(:,3) >= fTime & mCatalog(:,3) < fTime+fTimeWindow);
    nNumper = max(length(mCatalog(vSel,1)));
    fNumper = nNumper/fTimeWindow;
    fChTime = fNumper/fNum; % Relative change of entire data

    % Change in FMD
    [vFMD1, vBin1] = hist(mCatalog(vSel,6),fMinBin:0.1:fMaxBin);
    vSel2 = (mCatalog(:,3) >= fTime+fTimeWindow & mCatalog(:,3) < fTime+2*fTimeWindow);
    [vFMD2, vBin2] = hist(mCatalog(vSel2,6),fMinBin:0.1:fMaxBin);
    fChFMD = max(cumsum(abs(vFMD2-vFMD1)));

    mRes = [mRes; fTime fNumper fChTime fChFMD];

    % Shift time by half window width
    fTime = fTime+0.1; %fTimeWindow/10;
end

% Maximum positive change
vSel = (mRes(:,4) == max(mRes(:,4)));
mMax = mRes(vSel,:);
fTmax = mRes(vSel,1);

figure
hPlot1=plot(mRes(:,1),mRes(:,4)/mean(mRes(:,4)),'Color',[0 0 0],'Linewidth',2)
% sTitle = ['Time of maximum change: ' num2str(fTmax)]
% hTit=title(sTitle)
xlabel('Time [dec. year]','FontSize',12,'Fontweight','bold')
ylabel('Normalized \Delta_{FMD}','FontSize',12,'Fontweight','bold')
set(gca,'FontSize',12,'Fontweight','bold','Linewidth',2)
%set(hTit,'FontSize',12,'Fontweight','bold')

author: CelsoReyes


text:
function [N,P,M,abs_M] = readIPI(filename)
  % READIPI read ipiSofts skeleton format
  %
  % [N,P,M,abs_M] = readIPI(filename)
  %
  % Inputs:
  %   filename  path to ipiSoft file
  % Outputs:
  %   N  #B cell of bone names
  %   P  #B list of parent indices (0 means root)
  %   M  4 by 4 by #B list of relative transform matrices
  %   abs_M 4 by 4 by #B list of absolute transform matrices
  %
  fp = fopen(filename,'r');
  line = eat_comments(fp,'#');
  fscanf(fp,'\n');
  while true
    ii = 1;
    % read id number
    line = line(ii:end);
    [id,count,e,ii] = sscanf(line,'%d ',1);
    if count ~= 1
      error('Bad format');
      break;
    end

    % read name
    line = line(ii:end);
    [name,count,e,ii] = sscanf(line,'%s ',1);
    if count ~= 1
      error('Bad format');
      break;
    end
    N{id} = name;

    % read parent
    line = line(ii:end);
    [p,count,e,ii] = sscanf(line,'%d ',1);
    if count ~= 1
      error('Bad format');
      break;
    end
    P(id) = p;

    % read parent
    line = line(ii:end);
    [m,count,e,ii] = sscanf(line,'%g ',16);
    if count ~= 16
      error('Bad format');
      break;
    end
    M(:,:,id) = reshape(m,4,4)';

    line = fgets(fp);
    if line == -1
      break;
    end
  end
  fclose(fp);

  next = P;
  abs_M = M;
  while true
    % loop over bones
    for ii = 1:numel(P)
      % if not root
      if next(ii) ~= 0
        abs_M(:,:,ii) = M(:,:,next(ii)) * abs_M(:,:,ii);
        next(ii) = P(next(ii));
      end
    end
    if ~any(next)
      break;
    end
  end

  assert(all(P<=numel(P)));
end

author: alecjacobson

repo: gptoolbox


text:
function varargout = drawCube(cube, varargin)
%DRAWCUBE Draw a 3D centered cube, eventually rotated.
%
%   drawCube(CUBE)
%   Displays a 3D cube on current axis. CUBE is given by:
%   [XC YC ZC SIDE THETA PHI PSI]
%   where (XC, YC, ZC) is the CUBE center, SIDE is the length of the cube
%   main sides, and THETA PHI PSI are angles representing the cube
%   orientation, in degrees. THETA is the colatitude of the cube, between 0
%   and 90 degrees, PHI is the longitude, and PSI is the rotation angle
%   around the axis of the normal.
%
%   CUBE can be axis aligned, in this case it should only contain center
%   and side information:
%   CUBE = [XC YC ZC SIDE]
%
%   The function drawCuboid is closely related, but uses a different angle
%   convention, and allows for different sizes along directions.
%
%   Example
%   % Draw a cube with small rotations
%     figure; hold on;
%     drawCube([10 20 30  50  10 20 30], 'FaceColor', 'g');
%     axis equal;
%     view(3);
%
%   See also
%   meshes3d, polyhedra, createCube, drawCuboid
%

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-06-29,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

% default values
theta = 0;
phi   = 0;
psi   = 0;

%% Parses the input
if nargin == 0
    % no input: assumes cube with default shape
    xc = 0;	yc = 0; zc = 0;
    a = 1;

else
    % one argument: parses elements
    xc  = cube(:,1);
    yc  = cube(:,2);
    zc  = cube(:,3);
    a   = cube(:,4);

    % parses orientation if present
    k   = pi / 180;
    if size(cube, 2) >= 6
        theta = cube(:,5) * k;
        phi   = cube(:,6) * k;
    end
    if size(cube, 2) >= 7
        psi   = cube(:,7) * k;
    end
end

%% Compute cube coordinates

% create unit centered cube
[v, f] = createCube;
v = bsxfun(@minus, v, mean(v, 1));

% convert unit basis to cube basis
sca     = createScaling3d(a);
rot1    = createRotationOz(psi);
rot2    = createRotationOy(theta);
rot3    = createRotationOz(phi);
tra     = createTranslation3d([xc yc zc]);

% concatenate transforms
trans   = tra * rot3 * rot2 * rot1 * sca;

% transform mesh vertices
[x, y, z] = transformPoint3d(v, trans);

%% Process output
if nargout == 0
    % no output: draw the cube
    drawMesh([x y z], f, varargin{:});

elseif nargout == 1
    % one output: draw the cube and return handle
    varargout{1} = drawMesh([x y z], f, varargin{:});

elseif nargout == 3
    % 3 outputs: return computed coordinates
    varargout{1} = x;
    varargout{2} = y;
    varargout{3} = z;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function n = p36_n ( )

%*****************************************************************************80
%
%% P36_N returns the number of variables for problem 36.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer N.  If N is positive, N represents the
%    only legal value for N for this problem.  If N is
%    negative, the absolute value of N is the least legal
%    value of N, but other values are allowable.
%
  n = 2;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [etoe,idx] = EleToEle(t)
%--if an element shares an edge then it is a neighbor
%--returns element-to-element connectivity in crs format.
%--triangle ie is connected to etoe(idx(ie):idx(ie+1)-1,2) triangles
% kjr 2017, fixed August 2017 for self adj triangles (i.e., disjoint
% elements).
t = sort(t,2);
edges = [t(:,[1 2]);t(:,[1 3]);t(:,[2 3])];
nt = size(t,1);
trinum = repmat((1:nt)',3,1);
% use sort by converting edges to bytes for speed.
[~,tags]=sort(edges*[2^31;1]);
edges=edges(tags,:);
%[edges,tags] = sortrows(edges);
trinum = trinum(tags);
k = find(all(diff(edges,1)==0,2));
etoe=trinum([k,k+1]);
[~,dmy1]=sort(etoe(:,1));[~,dmy2]=sort(etoe(:,2));
% use sort by converting edges to bytes for speed
temp=[etoe(dmy1,:); fliplr(etoe(dmy2,:)); (1:nt)',(1:nt)'];
[~,tags]=sort(temp*[2^31;1]);
etoe=temp(tags,:);
%etoe=sortrows([etoe(dmy1,:); fliplr(etoe(dmy2,:)); (1:nt)',(1:nt)']);%added self adjs.
idx = find(diff(etoe(:,1))==1); idx=[idx;length(etoe)]; idx = [1;idx+1];
end

author: CHLNDDEV


text:
%BLUR  Smooths an image using the normalized box filter
%
%     dst = cv.blur(src)
%     dst = cv.blur(src, 'OptionName',optionValue, ...)
%
% ## Input
% * __src__ input image; it can have any number of channels, which are
%   processed independently, but the depth should be `uint8`, `uint16`,
%   `int16`, `single`, or `double`.
%
% ## Output
% * __dst__ output image of the same size and type as `src`.
%
% ## Options
% * __KSize__ blurring kernel size. default [5,5]
% * __Anchor__ Anchor point `[x,y]`. The default value `[-1,-1]` means that
%   the anchor is at the kernel center.
% * __BorderType__ Border mode used to extrapolate pixels outside of the
%   image. See cv.copyMakeBorder. default 'Default'
%
% The function smooths an image using the kernel:
%
%     K = ones(KSize) / prod(KSize)
%
% See also: cv.boxFilter, cv.bilateralFilter, cv.GaussianBlur, cv.medianBlur,
%  imfilter, fspecial
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function tagS = MR_image_module_tags_subset
%"MR_image_module_tags_subset"
%   Return a subset of the tags specified in C.7.3.1 in PS3.3
%   of the 2006 DICOM specification required for distinguishing MR series'.
%
%   Tags are returned in a struct array with 3 fields:
%
%   Tag: String containing hex DICOM tag of a field.
%  Type: String describing type of field, with 5 options:
%         '1' Field must exist, data must exist and be valid.
%         '2' Field must exist, data can exist and be valid, or be NULL.
%         '3' Field is optional, if the field exists data can exist and be
%             valid or be NULL.
%         '1C' Field must exist under certain conditions and contain valid
%             data.
%         '2C' Field must exist under certain conditions and can contain
%             valid data or be NULL.
%Children: For sequences, a tagS with the same format as this struct array.
%
%AI 11/20/2020
%
%Usage:
%   tagS = MR_image_module_tags_subset
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

%Initialize the tagS structure.
tagS = struct('tag', {}, 'tagdec', {}, 'type', {}, 'children', {});

%Create an empty tagS template for sequence creation.
template = tagS;

%Add tags based on PS3.3 attribute lists.

%Modality
tagS(end+1) = struct('tag', '00080060', 'tagdec', 524384, 'type', '1', 'children', []);

% Acquisition Time
tagS(end+1) = struct('tag', '00080032', 'tagdec', 524338, 'type', '3', 'children', []);

% b-value for MR scans
tagS(end+1) = struct('tag', '00431039', 'tagdec', 4395065, 'type', '4', 'children', []); % GE
tagS(end+1) = struct('tag', '00189087', 'tagdec', 1609863, 'type', '4', 'children', []); % Philips
tagS(end+1) = struct('tag', '0019100C', 'tagdec', 1642508, 'type', '4', 'children', []); % SIEMENS

% Temporal position ID
tagS(end+1) = struct('tag', '00200100', 'tagdec', 2097408, 'type', '3', 'children', []);

%Trigger time
tagS(end+1) = struct('tag', '00181060', 'tagdec', 1577056, 'type', '2C', 'children', []);

%Number of slices (for GE data)
tagS(end+1) = struct('tag', '0021104F', 'tagdec', 2166863, 'type', '3', 'children', []);

%Instance number
tagS(end+1) = struct('tag', '00200013', 'tagdec', 2097171, 'type', '1', 'children', []);

%Manufacturer
tagS(end+1) = struct('tag', '00080070', 'tagdec', 524400, 'type', '2', 'children', []);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
classdef (Abstract) BASEsettings < matlab.mixin.CustomDisplay
    %BASESETTINGS Abstract class with the common settings to all VARtypes
    %
    % BASEsettings Properties:
    %    frequency       - data frequency
    %    startdate       - sample start date
    %    enddate         - sample end date
    %    varendo         - endogenous variables;
    %    varexo          - exogenous variables
    %    lags            - number of lags
    %    const           - inclusion of a constant
    %    excelFile       - Excel file used for the inputs
    %    results_path    - path where there results file is stored
    %    results_sub     - name of the results file
    %    results         - save the results in the excel file (true/false)
    %    plot            - plot the results (true/false)
    %    workspace       - save the workspace as a .mat file (true/false)
    %    IRF             - impulse response functions
    %    IRFperiods      - number of periods for IRF
    %    F               - unconditional forecasts
    %    FEVD            - forecast error variance decomposition
    %    HD              - historical decomposition
    %    HDall           - plot the entire decomposition
    %    CF              - activate conditional forecasts
    %    IRFt            - structural identification
    %    Feval           - forecast evaluation
    %    CFt             - type of conditional forecasts
    %    Fstartdate      - start date for forecasts
    %    Fenddate        - start date for forecasts
    %    Fendsmpl        - start forecasts after the final sample period
    %    hstep           - step ahead evaluation
    %    window_size     - window_size for iterative forecasting
    %    evaluation_size - evaluation_size as percent of window_size
    %    cband           - confidence level for VAR coefficients
    %    IRFband         - confidence level for impusle response functions
    %    Fband           - confidence level for forecasts
    %    FEVDband        - confidence level for forecast error variance decomposition
    %    HDband          - confidence level for historical decomposition

    properties (SetAccess = private)

        VARtype   bear.VARtype = bear.VARtype.empty      % VAR model selected (1=OLS VAR, 2=BVAR, 3=mean-adjusted BVAR, 4=panel Bayesian VAR, 5=Stochastic volatility BVAR, 6=Time varying)

    end

    properties
        %% App settings on constant panel
        frequency (1,1) double  = 2;                     % data frequency (1=yearly, 2= quarterly, 3=monthly, 4=weekly, 5=daily, 6=undated)
        startdate               = '1974q1';              % sample start date; must be a string consistent with the date formats of the toolbox
        enddate                 = '2014q4';              % sample end date; must be a string consistent with the date formats of the toolbox
        varendo                 = 'YER HICSA STN';       % endogenous variables; must be a single string, with variable names separated by a space
        varexo                  = '';                    % exogenous variables, if any; must be a single string, with variable names separated by a space
        lags      (1,1) double  = 4;                     % number of lags
        const     (1,1) logical = true;                  % inclusion of a constant (1=yes, 0=no)

        excelFile    (1,:) char = '';                    % Excel file used for the inputs
        results_path (1,:) char = '';                    % path where there results file is stored
        results_sub  (1,:) char = 'results';             % name of the results file
        results      (1,1) logical = true;               % save the results in the excel file (true/false)
        plot         (1,1) logical = true;               % plot the results (true/false)
        workspace    (1,1) logical = true;               % save the workspace as a .mat file (true/false)

        % Model options
        IRF        (1,1) logical = true;  % activate impulse response functions (1=yes, 0=no)
        IRFperiods (1,1) double  = 20;    % number of periods for impulse response functions
        % structural identification:
        % 1=none,
        % 2=Cholesky,
        % 3=triangular factorisation
        % 4=sign, zero, magnitude, relative magnitude, FEVD, correlation restrictions,
        % 5=IV identification,
        % 6=IV identification & sign, zero, magnitude, relative magnitude, FEVD, correlation restrictions)
        IRFt  bear.IRFtype = bear.IRFtype(2);
    end

    properties (Dependent)
        FEVD       % activate forecast error variance decomposition (1=yes, 0=no)
        HD         % activate historical decomposition (1=yes, 0=no)
    end

    properties

        HDall      (1,1) logical = 0;     % if we want to plot the entire decomposition, all contributions (includes deterministic part) (1=yes, 0=no)
        F          (1,1) logical = true;  % activate unconditional forecasts (1=yes, 0=no)

        CF         (1,1) logical = false; % activate conditional forecasts (1=yes, 0=no)
        % type of conditional forecasts
        % 1 = standard (all shocks),
        % 2 = standard (shock-specific)
        % 3 = tilting (median), 4=tilting (interval)
        CFt bear.CFtype = bear.CFtype(1);

        % start forecasts immediately after the final sample period (1=yes, 0=no)
        % has to be set to 1 if start date for forecasts is not in-sample
        Fendsmpl (1,1) logical = false;
        Fstartdate = '2014q1'; % start date for forecasts (has to be an in-sample date; otherwise, ignore and set Fendsmpl=1)
        Fenddate   = '2016q4'; % end date for forecasts
        Feval (1,1) logical = false; % activate forecast evaluation (1=yes, 0=no)

        hstep           (1,1) double = 1;    % step ahead evaluation
        window_size     (1,1) double = 0;    % window_size for iterative forecasting 0 if no iterative forecasting                                            <                                                                                -
        evaluation_size (1,1) double = 0.5;  % evaluation_size as percent of window_size
        cband           (1,1) double = 0.68; % confidence/credibility level for VAR coefficients
        IRFband         (1,1) double = 0.68; % confidence/credibility level for impusle response functions
        Fband           (1,1) double = 0.68; % confidence/credibility level for forecasts
        FEVDband        (1,1) double = 0.68; % confidence/credibility level for forecast error variance decomposition
        HDband          (1,1) double = 0.68; % confidence/credibility level for historical decomposition

    end

    properties (Access = protected)
        FEVDinternal (1,1) logical = true;
        HDinternal   (1,1) logical = true;
    end

    methods

        function obj = BASEsettings(VARtype, excelPath)

            obj.VARtype = VARtype;
            obj.excelFile = excelPath;
            obj.results_path = pwd();

        end

        function obj = set.IRFt(obj, value)
            % call another function to check the value as desired,
            % and possibly even update it using some computation
            obj = checkIRFt(obj, value);

            % set set the property using the validated value
            % (only place we do assignment to avoid infinite recursion)
            obj.IRFt = value;
        end

        function value = get.FEVD(obj)
            value = getFEVD(obj);
        end

        function value = get.HD(obj)
            value = getHD(obj);
        end

        function obj = set.FEVD(obj, value)
            obj.FEVDinternal = value;
        end

        function obj = set.HD(obj, value)
            obj.HDinternal = value;
        end

    end

    methods (Access = protected)

        function obj = parseBEARSettings( obj, varargin )

            obj = bear.utils.pvset(obj, varargin{:});

        end

        function obj = checkIRFt(obj, ~)
        end

        function displayScalarObject(obj)

            % Grab property lists
            props = properties(obj)';
            meta = ?bear.settings.BASEsettings;

            % Get base properties
            baseProps = {meta.PropertyList.Name};

            specificProps = setdiff(props, baseProps);
            specificProps = props(ismember(props, specificProps)); % To keep original order

            mainProps = {'VARtype', 'frequency', 'startdate', ...
                'enddate', 'varendo', 'varexo', 'lags', 'const', ...
                'excelFile', 'results_path', 'results_sub', 'results', 'plot', 'workspace'};

            applicationProps = setdiff(baseProps, mainProps);
            applicationProps = props(ismember(props, applicationProps)); % To keep original order

            % header
            header = matlab.mixin.CustomDisplay.getSimpleHeader(obj);
            disp(header);

            % Preferences
            fprintf('\n --------- PREFERENCES ----------\n\n');
            mainProps = matlab.mixin.util.PropertyGroup(mainProps);
            matlab.mixin.CustomDisplay.displayPropertyGroups(obj, mainProps);

            % Specifications
            fprintf('\n --------- SPECIFICATIONS ----------\n\n');
            specificProps = matlab.mixin.util.PropertyGroup(specificProps);
            matlab.mixin.CustomDisplay.displayPropertyGroups(obj, specificProps);

            % Applications
            fprintf('\n --------- APPLICATIONS ----------\n\n');
            applicationProps = matlab.mixin.util.PropertyGroup(applicationProps);
            matlab.mixin.CustomDisplay.displayPropertyGroups(obj, applicationProps);

        end

    end

    methods (Access = protected)

        function value = getFEVD(obj)
            value = obj.FEVDinternal;
        end

        function value = getHD(obj)
            value = obj.HDinternal;
        end

    end

end

author: european-central-bank

repo: BEAR-toolbox


text:
function Z = zstream(action,D)
% Compress/decompress stream of bytes using Deflate/Inflate
% FORMAT Z = zstream('C',D)
% D        - data stream to compress (converted to uint8 if needed)
% Z        - compressed data stream (uint8 vector)
% FORMAT D = zstream('D',Z)
% Z        - data stream to decompress (uint8 vector)
% D        - decompressed data stream (uint8 vector)
%__________________________________________________________________________
%
% This C-MEX file relies on:
% * miniz, by Rich Geldreich
%   http://code.google.com/p/miniz/
% Fallback Java implementation is adapted from:
% * dzip/dunzip, by Michael Kleder
%   http://www.mathworks.com/matlabcentral/fileexchange/8899
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: zstream.m 6417 2015-04-21 16:03:44Z guillaume $

if exist('OCTAVE_VERSION','builtin')
    error('zstream.c not compiled - see Makefile');
end

switch upper(action)
    case 'C'
        D = typecast(D(:),'uint8');
        f = java.io.ByteArrayOutputStream();
        g = java.util.zip.DeflaterOutputStream(f);
        g.write(D);
        g.close;
        Z = typecast(f.toByteArray,'uint8');
        f.close;

    case 'D'
        import com.mathworks.mlwidgets.io.InterruptibleStreamCopier
        a   = java.io.ByteArrayInputStream(D);
        b   = java.util.zip.InflaterInputStream(a);
        isc = InterruptibleStreamCopier.getInterruptibleStreamCopier;
        c   = java.io.ByteArrayOutputStream;
        isc.copyStream(b,c);
        Z   = c.toByteArray;

    otherwise
        error('Unknown action "%s".',action);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% to test whether scg inference engine can handl dynameic BN
% Make a linear dynamical system
%   X1 -> X2
%   |     |
%   v     v
%   Y1    Y2

intra = zeros(2);
intra(1,2) = 1;
inter = zeros(2);
inter(1,1) = 1;
n = 2;

X = 2; % size of hidden state
Y = 2; % size of observable state

ns = [X Y];
dnodes = [];
onodes = [2];
eclass1 = [1 2];
eclass2 = [3 2];
bnet = mk_dbn(intra, inter, ns, dnodes, eclass1, eclass2);

x0 = rand(X,1);
V0 = eye(X);
C0 = rand(Y,X);
R0 = eye(Y);
A0 = rand(X,X);
Q0 = eye(X);

bnet.CPD{1} = gaussian_CPD(bnet, 1, 'mean', x0, 'cov', V0);
%bnet.CPD{2} = gaussian_CPD(bnet, 2, 'mean', zeros(Y,1), 'cov', R0, 'weights', C0, 'full', 'untied', 'clamped_mean');
%bnet.CPD{3} = gaussian_CPD(bnet, 3, 'mean', zeros(X,1), 'cov', Q0, 'weights', A0, 'full', 'untied', 'clamped_mean');
bnet.CPD{2} = gaussian_CPD(bnet, 2, 'mean', zeros(Y,1), 'cov', R0, 'weights', C0);
bnet.CPD{3} = gaussian_CPD(bnet, 3, 'mean', zeros(X,1), 'cov', Q0, 'weights', A0);

T = 5; % fixed length sequences

clear engine;
%engine{1} = kalman_inf_engine(bnet, onodes);
engine{1} = scg_unrolled_dbn_inf_engine(bnet, T, onodes);
engine{2} = jtree_unrolled_dbn_inf_engine(bnet, T);

N = length(engine);

% inference

ev = sample_dbn(bnet, T);
evidence = cell(n,T);
evidence(onodes,:) = ev(onodes, :);

t = 2;
query = [1 3];
m = cell(1, N);
ll = zeros(1, N);

engine{1} = enter_evidence(engine{1}, evidence);
[engine{2}, ll(2)] = enter_evidence(engine{2}, evidence);
m{1} = marginal_nodes(engine{1}, query);
m{2} = marginal_nodes(engine{2}, query, t);

% compare all engines to engine{1}
for i=2:N
  assert(approxeq(m{1}.mu, m{i}.mu));
  assert(approxeq(m{1}.Sigma, m{i}.Sigma));
%  assert(approxeq(ll(1), ll(i)));
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function s = imMedian(img, varargin)
%IMMEDIAN Median value of a grayscale image, or of each color component
%
%   S = imMedian(IMG)
%   Computes the median value of pixels in image IMG. If image is grayscale
%   image, the result is a scalar. If image is a color image, the result is
%   1-by-3 row vector, each componenent corresponding to one color of the
%   image.
%
%   S = imMedian(IMG, MASK)
%   Computes the median value only in the area specified by MASK.
%
%   S = imMedian(..., 'color', COL)
%   Forces the function to consider the image as color (if COL is TRUE) or
%   as grascale (if COL is FALSE). This can be useful for vector image with
%   more than 3 color components.
%
%
%   Example
%   % apply to cameraman image
%   img = imread('cameraman.tif');
%   imMedian(img)
%   ans =
%       144
%
%   % apply to a RGB image
%   img = imread('peppers.png');
%   imMedian(img)
%   ans =
%       90    43    60
%
%   See also
%   imMin, imMax, imMean, imMode
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-07-30,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

%% Process input arguments

% detect if image is color
color = isColorImage(img);

% check if user specified 'color' option
if length(varargin)>1
    var = varargin{end-1};
    if ischar(var)
        if strcmpi(var, 'color')
            color = varargin{end};
            varargin(end-1:end) = [];
        end
    end
end

%% Process color image

if color
    % If image is color, process each band separately

    % compute image size and dimension (including color dimension)
    dim = size(img);
    nd = length(dim);

    % create idnexing structure
    inds = cell(1, nd);
    for i=1:nd
        inds{i} = 1:dim(i);
    end

    % iterate on colors
    nc = dim(3);
    s = zeros(1, nc);
    for i=1:nc
        % modify the indexing structure to work on the i-th component
        inds{3} = i;
        s(i) = imMedian(img(inds{:}), varargin{:});
    end

    return;
end

%% process grayscale image

if isempty(varargin)
    % compute sum over all image
    s = median(img(:));
else
    % use first argument as mask
    s = median(img(varargin{1}));
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function net = net_init_cifar_cnn(opts)
% CNN_MNIST_LENET Initialize a CNN similar for MNIST

rng('default');
rng(0) ;

f=1/100 ;
net.layers = {} ;
% Block 1
net.layers{end+1} = struct('type', 'conv', ...
                           'weights', {{f*randn(5,5,3,32, 'single'), zeros(1, 32, 'single') }}, ...
                           'stride', 1, ...
                           'pad', 2) ;
net.layers{end+1} = struct('type', 'relu') ;

net.layers{end+1} = struct('type', 'pool','pool', 3, 'stride', 2,'pad', [1,1,1,1]) ;
% Block 2

net.layers{end+1} = struct('type', 'conv', ...
                           'weights', {{f*randn(5,5,32,32, 'single'), zeros(1,32, 'single')}}, ...
                           'stride', 1, ...
                           'pad', 2) ;
net.layers{end+1} = struct('type', 'relu') ;

net.layers{end+1} = struct('type', 'pool','pool', 3, 'stride', 2,'pad', [1,1,1,1]) ;

% Block 3

net.layers{end+1} = struct('type', 'conv', ...
                           'weights', {{f*randn(5,5,32,64, 'single'), zeros(1,64, 'single')}}, ...
                           'stride', 1, ...
                           'pad', 2) ;
net.layers{end+1} = struct('type', 'relu') ;

net.layers{end+1} = struct('type', 'pool','pool', 3, 'stride', 2,'pad', [1,1,1,1]) ;

% Block 4
net.layers{end+1} = struct('type', 'conv', ...
                           'weights', {{f*randn(4,4,64,64, 'single'), zeros(1,64, 'single')}}, ...
                           'stride', 1, ...
                           'pad', 0) ;
net.layers{end+1} = struct('type', 'relu') ;

% Block 5

net.layers{end+1} = struct('type', 'conv', ...
                           'weights', {{f*randn(1,1,64,10, 'single'), zeros(1, 10, 'single')}}, ...
                           'stride', 1, ...
                           'pad', 0) ;

net.layers{end+1} = struct('type', 'softmaxloss') ;

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet


text:
function line = sphere_llq_grid_lines ( nlat, nlong, line_num )

%*****************************************************************************80
%
%% SPHERE_LLQ_GRID_LINES returns lines for an LLQ grid on a sphere.
%
%  Discussion:
%
%    A SPHERE LLQ grid imposes a grid of quadrilaterals on a sphere,
%    using latitude and longitude lines.
%
%    The point numbering system is the same used in SPHERE_LLQ_GRID_POINTS,
%    and that routine may be used to compute the coordinates of the points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 April 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NLAT, NLONG, the number of latitude and longitude
%    lines to draw.  The latitudes do not include the North and South
%    poles, which will be included automatically, so NLAT = 5, for instance,
%    will result in points along 7 lines of latitude.
%
%    Input, integer LINE_NUM, the number of grid lines.
%
%    Output, integer LINE(LINE_NUM,2), contains pairs of point indices for
%    line segments that make up the grid.
%
  l = 0;
  line = zeros ( line_num, 2 );
%
%  "Vertical" lines.
%
  for j = 0 : nlong - 1

    old = 1;
    new = j + 2;
    l = l + 1;
    line(l,1:2) = [ old, new ];

    for i = 1 : nlat - 1

      old = new;
      new = old + nlong;
      l = l + 1;
      line(l,1:2) = [ old, new ];

    end

    old = new;
    l = l + 1;
    line(l,1:2) = [ old, 1 + nlat * nlong + 1 ];

  end
%
%  "Horizontal" lines.
%
  for i = 1 : nlat

    new = 1 + ( i - 1 ) * nlong + 1;

    for j = 0 : nlong - 2
      old = new;
      new = old + 1;
      l = l + 1;
      line(l,1:2) = [ old, new ];
    end

    old = new;
    new = 1 + ( i - 1 ) * nlong + 1;
    l = l + 1;
    line(l,1:2) = [ old, new ];

  end

  return
end


text:
% Copyright (C) Daphne Koller, Stanford University, 2012

rand('seed', 1);

% Construct the toy network
[toy_network, toy_factors] = ConstructToyNetwork(1,0.2);
toy_evidence = zeros(1, length(toy_network.names));
%toy_clique_tree = CreateCliqueTree(toy_factors, []);
%toy_cluster_graph = CreateClusterGraph(toy_factors,[]);

% Exact Inference
ExactM = ComputeExactMarginalsBP(toy_factors, toy_evidence, 0);
figure, VisualizeToyImageMarginals(toy_network, ExactM,1,'exact');

% Comment this in to run Approximate Inference on the toy network
% Approximate Inference
% % ApproxM = ApproxInference(toy_cluster_graph, toy_factors, toy_evidence);
% figure, VisualizeToyImageMarginals(toy_network, ApproxM);

% MCMC Inference
transition_names = {'Gibbs', 'MHUniform', 'MHGibbs', 'MHSwendsenWang1', 'MHSwendsenWang2'};

for j = 1:length(transition_names)
    samples_list = {};

    num_chains_to_run = 1;
    for i = 1:num_chains_to_run
        % Random Initialization
        A0 = ceil(rand(1, length(toy_network.names)) .* toy_network.card);

        % Initialization to all ones
        % A0 = i * ones(1, length(toy_network.names));

        [M, all_samples] = ...
            MCMCInference(toy_network, toy_factors, toy_evidence, transition_names{j}, 0, 500, 1, A0);
        samples_list{i} = all_samples;
        figure, VisualizeToyImageMarginals(toy_network, M, i, transition_names{j});
    end

    vis_vars = [3];
    VisualizeMCMCMarginals(samples_list, vis_vars, toy_network.card(vis_vars), toy_factors, ...
      500, ExactM(vis_vars),transition_names{j});
    disp(['Displaying results for MCMC with transition ', transition_names{j}]);
    disp(['Hit enter to continue']);
    pause;
end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
classdef MimWindowLevelBone < MimGuiPlugin
    % MimWindowLevelBone. Gui Plugin for using a preset bone window/level
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimWindowLevelBone is a Gui Plugin for the MIM Toolkit.
    %     The gui will create a button for the user to run this plugin.
    %     Running this plugin will change the window and level of the viewing
    %     panel to standard soft tissue values.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Bone'
        SelectedText = 'Bone'
        ToolTip = 'Changes the window and level settings to standard bone values (Window 2000HU Level 300HU)'
        Category = 'Window/Level Presets'
        Visibility = 'Dataset'
        Mode = 'Toolbar'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 1
        Icon = 'wl_bone.png'
        Location = 22
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            gui_app.ImagePanel.Window = 2000;
            gui_app.ImagePanel.Level = 300;
        end

        function enabled = IsEnabled(gui_app)
            enabled = gui_app.IsDatasetLoaded && gui_app.IsCT && ~strcmp(gui_app.ImagePanel.Mode, MimModes.View3DMode);
        end

        function is_selected = IsSelected(gui_app)
            is_selected = gui_app.ImagePanel.Window == 2000 && gui_app.ImagePanel.Level == 300;
        end

    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [detJ,J] = dtiJacobian(inData)
%
% [detJ,J] = dtiJacobian(inData)
%
% Computes the Jacobian (J) and the determinant of the Jacobian (detJ)
% given an XxYxZx3 or Nx3 input array. The array is assumed to represent a
% vector field. The Jacobian is a useful metric for interpreting
% deformations. E.g., abs(detJ) represents the local volume change for a
% deformation field.
%
% HISTORY:
% 2008.06.25 RFD wrote it.

sz = size(inData);

if((numel(sz)~=4 && numel(sz)~=2) || sz(end)~=3)
    error('inData must be XxYxZx3 or Nx3.');
end

% Approximate the Jacobian of the input array using grad.

if(numel(sz)==2)
    inData = reshape(inData,[1 1 1 sz(1) sz(2)]);
end

dim = size(inData);
J = zeros(dim(1),dim(2),dim(3),3,3);

for(ii=1:3)
    [gradX,gradY,gradZ] = gradient(inData(:,:,:,ii),1);
    J(:,:,:,ii,1) = gradX;
    J(:,:,:,ii,2) = gradY; % Or Y,X,Z?
    J(:,:,:,ii,3) = gradZ;
end
% Compute the determinant of J
% det(J) for a 3x3 [a b c; d e f; g h i] is: (aei+bfg+cdh)-(gec+hfa+idb)
% a=1,1; b=1,2; c=1,3; d=2,1; e=2,2; f=2,3; g=3,1; h=3,2; i=3,3;
detJ =  J(:,:,:,1,1).*J(:,:,:,2,2).*J(:,:,:,3,3) ...
      + J(:,:,:,1,2).*J(:,:,:,2,3).*J(:,:,:,3,1) ...
      + J(:,:,:,1,3).*J(:,:,:,2,1).*J(:,:,:,3,2) ...
      - J(:,:,:,3,1).*J(:,:,:,2,2).*J(:,:,:,1,3) ...
      - J(:,:,:,3,2).*J(:,:,:,2,3).*J(:,:,:,1,1) ...
      - J(:,:,:,3,3).*J(:,:,:,2,1).*J(:,:,:,1,2);

if(numel(sz)==2)
    J = squeeze(J);
    detJ = squeeze(detJ);
end
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function exportAISegToDICOM(planC,origScanNum,outputDicomPath,fname,dcmExportOptS,allLabelNamesC)
% function exportAISegToDICOM(planC,outputDicomPath,fname,dcmExportOptS,allLabelNamesC)
%
% APA, 4/28/2022

indexS = planC{end};

%Get custom options for DICOM export
% Handle special case of assignig reference UID fromanother structure
% (e.g. exporting registered images from MIM assistant changes their frameOfreferenceUID)
structRefForC = {};
count = 0;
if exist('dcmExportOptS','var') && isstruct(dcmExportOptS)
    structNameC = {planC{indexS.structures}.structureName};
    for iDcmOpt = 1:length(dcmExportOptS)
        if isfield(dcmExportOptS(iDcmOpt),'rt_struct') && ...
                isfield(dcmExportOptS(iDcmOpt).rt_struct,'referencedFrameOfReference')
            % Update to handle multiple structures
            toStructureName = dcmExportOptS(iDcmOpt).rt_struct.referencedFrameOfReference.toStructureName;
            fromStructureName = dcmExportOptS(iDcmOpt).rt_struct.referencedFrameOfReference.fromStructureName;
            strIndex = getMatchingIndex(fromStructureName,structNameC,'exact');
            assocScanV = getStructureAssociatedScan(strIndex,planC);
            strIndex = strIndex(assocScanV==origScanNum);
            if isempty(strIndex)
                continue
            end
            count = count + 1;
            structRefFrameOfReferenceUID = planC{indexS.structures}(strIndex).referencedFrameOfReferenceUID;
            refSeriesInstanceUID = planC{indexS.structures}(strIndex).referencedSeriesUID;
            sopClassUidC = {planC{indexS.structures}(strIndex).contour.referencedSopClassUID};
            sopInstanceUidC = {planC{indexS.structures}(strIndex).contour.referencedSopInstanceUID};
            structRefForC{count,1} = toStructureName;
            structRefForC{count,2} = structRefFrameOfReferenceUID;
            structRefForC{count,3} = refSeriesInstanceUID;
            structRefForC{count,4} = sopClassUidC;
            structRefForC{count,5} = sopInstanceUidC;
        end
    end
end

%Retain only user-input structures (in allLabelNamesC)
numStr = length(planC{indexS.structures});
allStrC = {planC{indexS.structures}.structureName};
keepStrNumV = zeros(length(allLabelNamesC),1);
for i=1:length(allLabelNamesC)
    matchStrNum = getMatchingIndex(allLabelNamesC{i},allStrC,'EXACT');
    keepStrNumV(i) =  max(matchStrNum);
end

allStrV = 1:numStr;
allStrV(keepStrNumV)=[];

planC = deleteStructure(planC,allStrV);

%Export DICOM RT structs to outputDicomPath
planC = generate_DICOM_UID_Relationships(planC,structRefForC);
export_RS_IOD(planC,outputDicomPath,fname);


text:
function [ x, w ] = rule04 ( n )

%*****************************************************************************80
%
%% RULE04 returns the rule of degree 4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 July 2014
%
%  Author:
%
%    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Hong Xiao, Zydrunas Gimbutas,
%    A numerical algorithm for the construction of efficient quadrature
%    rules in two and higher dimensions,
%    Computers and Mathematics with Applications,
%    Volume 59, 2010, pages 663-676.
%
%  Parameters:
%
%    Input, integer N, the number of nodes.
%
%    Output, real X(3,N), the coordinates of the nodes.
%
%    Output, real W(N), the weights.
%
  xs = [ ...
    -.1459612280979987E-01,-.2937244896309993, ...
    0.3230878337158274,-.8819654179396583E-01, ...
    0.1270203663314710,-.4414307688091184, ...
    0.3939418589633429,-.6034469714614210E-01, ...
    -.8914059834601185E-01,0.6545213033603760, ...
    -.7307642259677857 ];
  ys = [ ...
    0.2426564579199708E-02,0.1764396506764613, ...
    -.1932070410956556,0.2317884270105980E-01, ...
    0.5451410677215219,-.3848225631590180, ...
    0.2068744670639530,-.4573739074927080, ...
    0.7268092659599459,-.3970356243870571, ...
    -.2669420088648982 ];
  zs = [ ...
    0.7093565780103633,0.1108860494941134, ...
    0.1695426396704650,-.2011819391325586, ...
    0.1309503990759315,0.9225429679162532E-01, ...
    -.3111560426198242,-.3302215329322376, ...
    -.3507931737363739,-.2970424833951137, ...
    -.3861037570241846 ];
  ws = [ ...
    0.1033787090646894,0.1070356256090164, ...
    0.1504792582740940,0.1877987156186583, ...
    0.7395274312521298E-01,0.7712199925411270E-01, ...
    0.6745419368701999E-01,0.5819413648173244E-01, ...
    0.5378646711152148E-01,0.5366953183949744E-01, ...
    0.3811216334909158E-01 ];

  x(1,1:n) = xs(1:n);
  x(2,1:n) = ys(1:n);
  x(3,1:n) = zs(1:n);
  w(1:n) = ws(1:n);

  return
end


text:
function boundaryBenchHtml(baseDir)
% function boundaryBenchHtml(baseDir)
%
% See also boundaryBench, boundaryBenchGraphs.
%
% David Martin <dmartin@eecs.berkeley.edu>
% May 2003

bsdsURL = 'http://www.cs.berkeley.edu/projects/vision/grouping/segbench/BSDS300';
presentations = {'gray','color'};
presNames = {'Grayscale','Color'};
iidsTest = imgList('test');

% Infer list of algorithms from directories present.
for k = 1:length(presentations),
  pres = presentations{k};
  dirlist = dir(fullfile(baseDir,pres));
  algs{k} = {};
  for i = 1:length(dirlist),
    alg = dirlist(i).name;
    if alg(1)=='.', continue; end
    if strcmp(alg,'human'), continue; end
    if ~isdir(fullfile(baseDir,pres,alg)), continue; end
    if length(dir(fullfile(baseDir,pres,alg,'name.txt')))~=1, continue; end
    if length(dir(fullfile(baseDir,pres,alg,'score.txt')))~=1, continue; end
    fname = fullfile(baseDir,pres,alg,'scores.txt');
    if length(dir(fname))~=1, continue; end
    tmp = dlmread(fname); % iid,thresh,r,p,f
    if size(tmp,1)~=numel(iidsTest), continue; end
    algs{k}{1+numel(algs{k})} = alg;
  end
end

% Read in all the scores.
for k = 1:length(presentations),
  pres = presentations{k};
  for i = 1:numel(algs{k}),
    alg = algs{k}{i};
    fprintf(2,'Processing directory %s/%s...\n',pres,alg);
    fname = fullfile(baseDir,pres,alg,'scores.txt');
    iscores{k}{i} = dlmread(fname); % iid,thresh,r,p,f
    fname = fullfile(baseDir,pres,alg,'score.txt');
    ascores{k}{i} = dlmread(fname); % thresh,r,p,f
    fname = fullfile(baseDir,pres,alg,'name.txt');
    fp = fopen(fname,'r');
    names{k}{i} = fgetl(fp);
    fclose(fp);
  end

  fname = fullfile(baseDir,pres,'human','scores.txt');
  if numel(dir(fname))>0,
    hiscores{k} = dlmread(fname); % iid,r,p,f
  else
    hiscores{k} = zeros(10000,4);
  end

  fname = fullfile(baseDir,pres,'human','score.txt');
  if numel(dir(fname))>0,
    hscore{k} = dlmread(fname); % r,p,f
  else
    hscore{k} = zeros(3,1);
  end
end

% Sort algorithms by overall F measure.
fprintf(2,'Sorting algorithms...\n');
for k = 1:length(presentations),
  pres = presentations{k};
  f = [];
  for i = 1:numel(algs{k}), f = [f ascores{k}{i}(4)]; end
  [unused,aperm{k}] = sort(f);
  aperm{k} = aperm{k}(end:-1:1);
end

% Create algorithms page.
fprintf(2,'Creating algorithms page...\n');
[u1,u2,u3] = mkdir(baseDir,'html');
mainPage = fullfile(baseDir,'html','algorithms.html');
fp = fopen(mainPage,'w');
title = 'Boundary Detection Benchmark: Algorithm Ranking';
fprintf(fp,'\n');
fprintf(fp,'\n');
fprintf(fp,'<title>%s</title>\n',title);
fprintf(fp,'\n');
fprintf(fp,'\n');
fprintf(fp,'
Go to Main Page.\n');
fprintf(fp,'
Go to Image Ranking.\n');
fprintf(fp,'
%s
\n',title);
fprintf(fp,'
Summary Tables
\n',presNames{k});
fprintf(fp,'

\n');
fprintf(fp,'
\n');
for k = 1:length(presentations),
  pres = presentations{k};
  fprintf(fp,'%s	\n',presNames{k});
end
fprintf(fp,'
\n');
fprintf(fp,'
\n');
for k = 1:length(presentations),
  fprintf(fp,'\n');
  pres = presentations{k};
  fprintf(fp,'

\n');
  fprintf(fp,'
Rank	Score	Algorithm
\n');
  fprintf(fp,'
\n');
  fprintf(fp,'0	\n');
  fprintf(fp,'%4.2f	\n',hscore{k}(3));
  fprintf(fp,'Humans	\n');
  fprintf(fp,'
\n');
  for i = 1:numel(algs{k}),
    alg = algs{k}{aperm{k}(i)};
    name = names{k}{aperm{k}(i)};
    score = ascores{k}{aperm{k}(i)}(4);
    fprintf(fp,'
\n');
    fprintf(fp,'%d	\n',i);
    fprintf(fp,'%4.2f	\n',score);
    fprintf(fp,'%s	\n',pres,alg,name);
    fprintf(fp,'
\n');
  end
  fprintf(fp,'
\n');
  fprintf(fp,'	\n');
end
fprintf(fp,'
\n');
fprintf(fp,'
\n');
fprintf(fp,'
Comparison Graphs
\n',presNames{k});
fprintf(fp,'

\n');
fprintf(fp,'
\n');
for k = 1:length(presentations),
  pres = presentations{k};
  fprintf(fp,'%s	\n',presNames{k});
end
fprintf(fp,'
\n');
fprintf(fp,'
\n');
for k = 1:length(presentations),
  fprintf(fp,'\n');
  fprintf(fp,'
\n');
  pres = presentations{k};
  for i = 1:numel(algs{k}),
    for j = i+1:numel(algs{k}),
      ii = aperm{k}(i);
      jj = aperm{k}(j);
      fprintf(fp,'
%s vs. %s\n',pres,i,j,names{k}{ii},names{k}{jj});
    end
  end
  fprintf(fp,'
\n');
  fprintf(fp,'	\n');
end
fprintf(fp,'
\n');
fprintf(fp,'
\n');
fprintf(fp,'
Detail Tables
\n',presNames{k});
for k = 1:length(presentations),
  pres = presentations{k};
  fprintf(fp,'
%s
\n',presNames{k});
  fprintf(fp,'

\n');
  fprintf(fp,'
Rank	Score	Algorithm	PR Curve
\n');
  fprintf(fp,'
\n');
  fprintf(fp,'0	\n');
  fprintf(fp,'%4.2f	\n',hscore{k}(3));
  fprintf(fp,'Humans	\n');
  fprintf(fp,'	\n',pres);
  fprintf(fp,'
\n');
  for i = 1:numel(algs{k}),
    alg = algs{k}{aperm{k}(i)};
    name = names{k}{aperm{k}(i)};
    score = ascores{k}{aperm{k}(i)}(4);
    fprintf(fp,'
\n');
    fprintf(fp,'%d	\n',i);
    fprintf(fp,'%4.2f	\n',score);
    fprintf(fp,'%s	\n',pres,alg,name);
    fprintf(fp,'	\n',pres,alg);
    fprintf(fp,'
\n');
  end
  fprintf(fp,'
\n');
end
fprintf(fp,'
Page generated on %s.\n',datestr(now));
fprintf(fp,'\n');
fprintf(fp,'\n');
fclose(fp);

% Create algorithm pages.
for k = 1:length(presentations),
  pres = presentations{k};
  for i = 1:numel(algs{k}),
    alg = algs{k}{i};
    fprintf(2,'Creating page for ''%s''...\n',alg);
    fname = fullfile(baseDir,pres,alg,'main.html');
    fp = fopen(fname,'w');
    title = sprintf('[%s] Boundary Detection Benchmark: Algorithm "%s"',presNames{k},names{k}{i});
    fprintf(fp,'\n');
    fprintf(fp,'\n');
    fprintf(fp,'<title>%s</title>\n',title);
    fprintf(fp,'\n');
    fprintf(fp,'\n');
    fprintf(fp,'
Back to Algorithm Ranking.\n');
    fprintf(fp,'
%s
\n',title);
    fname2 = fullfile(baseDir,pres,alg,'about.html');
    if length(dir(fname2))==1,
      fprintf(fp,'
\n');
      fp2 = fopen(fname2,'r');
      fwrite(fp,char(fread(fp2)'));
      fclose(fp2);
    end
    fprintf(fp,'
\n');
    fprintf(fp,'
Click on an image for additional details.\n');
    %fprintf(fp,'
Test Images
\n');
    fprintf(fp,'

\n');
    nper = 5;
    for j = 1:numel(iidsTest),
      iid = iidsTest(j);
      score = iscores{k}{i}(j,5);
      if iid~=iscores{k}{i}(j,1), error('bug'); end
      if mod(j,nper)==1, fprintf(fp,'
\n'); end
      fprintf(fp,'\n');
      fprintf(fp,'

\n');
      fprintf(fp,'#%d (%d) F=%4.2f
 ',j,iid,score);
      %fprintf(fp,'#%d (%d)
\n',j,iid);
      fprintf(fp,'\n',iid);
      fprintf(fp,'\n',bsdsURL,pres,iid);
      fprintf(fp,'\n');
      fprintf(fp,'

\n');
      fprintf(fp,'	\n');
      if mod(j,nper)==0 | j==numel(iidsTest), fprintf(fp,'
\n'); end
    end
    fprintf(fp,'
\n');
    fprintf(fp,'
Page generated on %s.\n',datestr(now));
    fprintf(fp,'\n');
    fprintf(fp,'\n');
    fclose(fp);
    % Create image pages for this algorithm.
    for j = 1:numel(iidsTest),
      iid = iidsTest(j);
      score = iscores{k}{i}(j,5);
      fname = fullfile(baseDir,pres,alg,sprintf('%d.html',iid));
      fp = fopen(fname,'w');
      title = sprintf(...
          '[%s] Boundary Detection Benchmark: Algorithm "%s" Image #%d (%d)',presNames{k},names{k}{i},j,iid);
      fprintf(fp,'\n');
      fprintf(fp,'\n');
      fprintf(fp,'<title>%s</title>\n',title);
      fprintf(fp,'\n');
      fprintf(fp,'\n');
      fprintf(fp,'
Back to Algorithm "%s" page.\n',names{k}{i});
      fprintf(fp,'
%s
\n',title);
      %fprintf(fp,'
F=%4.2f',score);
      fprintf(fp,'
\n',bsdsURL,pres,iid);
      fprintf(fp,'
\n');
      fprintf(fp,'
\n');
      fprintf(fp,'
Machine
	>\n',iid);
      fprintf(fp,'	\n',iid);
      fprintf(fp,'
\n');
      fprintf(fp,'
\n');
      fprintf(fp,'
Human
	>\n',iid);
      fprintf(fp,'	\n',iid);
      fprintf(fp,'
\n');
      fprintf(fp,'
\n');
      fprintf(fp,'
Page generated on %s.\n',datestr(now));
      fprintf(fp,'\n');
      fprintf(fp,'\n');
      fclose(fp);
    end
  end
end

% Sort the images by best F measure across all algorithms.
% Remember which algorithm had the best score for each image.
fprintf(2,'Sorting images...\n');
for k = 1:length(presentations),
  if numel(algs{k})==0, continue; end
  bestf{k} = zeros(numel(iidsTest),1);
  bestalg{k} = zeros(numel(iidsTest),1);
  for j = 1:numel(iidsTest),
    f = zeros(numel(algs{k}),1);
    for i = 1:numel(algs{k}), f(i) = iscores{k}{i}(j,5); end
    [bestf{k}(j),bestalg{k}(j)] = max(f);
  end
  [unused,iperm{k}] = sort(bestf{k});
  iperm{k} = iperm{k}(end:-1:1);
end

% Create images page.
fprintf(2,'Creating images page...\n');
mainPage = fullfile(baseDir,'html','images.html');
fp = fopen(mainPage,'w');
title = 'Boundary Detection Benchmark: Image Ranking';
fprintf(fp,'\n');
fprintf(fp,'\n');
fprintf(fp,'<title>%s</title>\n',title);
fprintf(fp,'\n');
fprintf(fp,'\n');
fprintf(fp,'
Go to Main Page.\n');
fprintf(fp,'
Go to Algorithm Ranking.\n');
fprintf(fp,'
%s
\n',title);
fprintf(fp,'
Click on an image for additional details.\n');
fprintf(fp,'

\n');
fprintf(fp,'
	\n');
for k = 1:length(presentations),
  fprintf(fp,'%s	\n',presNames{k});
end
fprintf(fp,'
\n');
fprintf(fp,'
Rank	\n');
for k = 1:length(presentations),
  fprintf(fp,'ID	Image	Best Algorithm [Score]	\n');
end
fprintf(fp,'
\n');

for j = 1:numel(iidsTest),
  fprintf(fp,'
\n');
  fprintf(fp,'%d	\n',j);
  for k = 1:length(presentations),
    if numel(algs{k})==0,
      fprintf(fp,'	\n');
      fprintf(fp,'	\n');
      fprintf(fp,'	\n');
      continue;
    end
    pres = presentations{k};
    iid = iidsTest(iperm{k}(j));
    f = bestf{k}(iperm{k}(j));
    maxf = hiscores{k}(iperm{k}(j),4);
    alg = algs{k}{bestalg{k}(iperm{k}(j))};
    name = names{k}{bestalg{k}(iperm{k}(j))};
    fprintf(fp,'#%d (%d)	\n',iperm{k}(j),iid);
    fprintf(fp,'\n');
    fprintf(fp,'\n',iid,pres);
    fprintf(fp,'\n',bsdsURL,pres,iid);
    fprintf(fp,'\n');
    fprintf(fp,'	\n');
    fprintf(fp,'%s
[ %4.2f / %4.2f ]	\n',pres,alg,name,f,maxf);
  end
  fprintf(fp,'
\n');
end
fprintf(fp,'
\n');
fprintf(fp,'
Page generated on %s.\n',datestr(now));
fprintf(fp,'\n');
fprintf(fp,'\n');
fclose(fp);

% Create image pages.
for j = 1:numel(iidsTest),
  iid = iidsTest(j);
  fprintf(2,'Creating pages for image #%d (%d)...\n',j,iid);
  for k = 1:length(presentations),
    if numel(algs{k})==0, continue; end
    pres = presentations{k};
    fname = fullfile(baseDir,'html',sprintf('%d-%s.html',iid,pres));
    fp = fopen(fname,'w');
    title = sprintf('[%s] Boundary Detection Benchmark: Image #%d (%d) Rank=%d',presNames{k},j,iid,find(iperm{k}==j));
    fprintf(fp,'\n');
    fprintf(fp,'\n');
    fprintf(fp,'<title>%s</title>\n',title);
    fprintf(fp,'\n');
    fprintf(fp,'\n');
    fprintf(fp,'
Back to Image Ranking.\n');
    fprintf(fp,'
Go to %s page.\n',iid,presentations{3-k},presNames{3-k});
    fprintf(fp,'
%s
\n',title);

    f = zeros(numel(algs{k}),1);
    for i = 1:numel(algs{k}), f(i) = iscores{k}{i}(j,5); end
    [unused,perm] = sort(f);
    perm = perm(end:-1:1);

    fprintf(fp,'

\n');
    fprintf(fp,'
Rank	Score	Algorithm	\n');
    fprintf(fp,'
\n');
    fprintf(fp,'0	\n');
    fprintf(fp,'%4.2f	\n',hiscores{k}(j,4));
    fprintf(fp,'Humans	\n');
    fprintf(fp,'
\n');
    for i = 1:numel(algs{k}),
      fprintf(fp,'
\n');
      fprintf(fp,'%d	\n',i);
      fprintf(fp,'%4.2f	\n',f(perm(i)));
      fprintf(fp,'%s	\n',names{k}{perm(i)});
      fprintf(fp,'
\n');
    end
    fprintf(fp,'
\n');

    fprintf(fp,'
\n');
    fprintf(fp,'\n',bsdsURL,pres,iid);
    fprintf(fp,'\n',pres,iid);

    fprintf(fp,'

\n');
    fprintf(fp,'
Rank Algorithm (Score)	Pb	Precision/Recall
\n');
    fprintf(fp,'
\n');
    fprintf(fp,'0
Humans
(%4.2f)	\n',hiscores{k}(j,4));
  fprintf(fp,'	\n',pres,iid);
  fprintf(fp,'	\n',pres,iid);
  fprintf(fp,'
\n');
    for i = 1:numel(algs{k}),
      fprintf(fp,'
\n');
      fprintf(fp,'%d
%s
(%4.2f)	\n',i,names{k}{perm(i)},f(perm(i)));
      fprintf(fp,'	\n',pres,algs{k}{perm(i)},iid);
      fprintf(fp,'	\n',pres,algs{k}{perm(i)},iid);
      fprintf(fp,'
\n');
    end
    fprintf(fp,'
\n');

    fprintf(fp,'
Page generated on %s.\n',datestr(now));
    fprintf(fp,'\n');
    fprintf(fp,'\n');
    fclose(fp);
  end
end

author: Cloud-CV

repo: object-proposals


text:
% APM Specify Options
%
% response = apm_option(server,app,name,value)
%
% Function apm_option sends an option specification
%   to the APM server with the following arguments:
%
%   server = address of server
%      app = application name
%     name = option name
%    value = option value
%
% A response is returned indicating whether the option
%   was successfully recorded
%
% Either global options (NLC.option_name) or parameter/variable
%   options (NAME.option_name) are able to be specified with this
%   function. A full list of available options with their default
%   values are available here:
%
% Global Options:
%    http://apmonitor.com/wiki/index.php/Main/DbsGlobal
%
% Variable Options:
%    http://apmonitor.com/wiki/index.php/Main/DbsVariable
%
function response = apm_option(server,app,name,value)
   aline = ['option ' deblank(char(name)) ' = ' num2str(value)];
   response = apm(server,app,aline);

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
function generateDoc
%generateDoc Publish the example scripts in the doc directory

%   Steven L. Eddins
%   Copyright 2008-2009 The MathWorks, Inc.

doc_dir = fullfile(fileparts(which('mtest')), '..', 'doc');
addpath(doc_dir);
cd(doc_dir)
mfiles = dir('*.m');
for k = 1:numel(mfiles)
    publish(mfiles(k).name);
    cd(doc_dir)
end

author: ricedsp


text:
function ress = times(a,bb)
% performs only interpolation nufft

% for n=1:size(bb,3)
% b = bb(:,:,n);
% if a.adjoint
% 	b = b(:).*a.w(:);
% 	res = nufft_adj(b, a.st)/sqrt(prod(a.imSize));;
% 	res = reshape(res, a.imSize(1), a.imSize(2));
% 	res = fft2c(res);
% else
% 	b = reshape(b,a.imSize(1),a.imSize(2));
% 	b = ifft2c(b);
% 	res = (nufft(b, a.st)/sqrt(prod(a.imSize))).*a.w(:);
% 	res = reshape(res,a.dataSize(1),a.dataSize(2));
% end
% ress(:,:,n) = res;
% end
%

ress = mtimes(a,bb);

author: thomaskuestner


text:
% PURPOSE  : 	Clears ERPset(s)
%
% FORMAT   :
%
% >> ALLERP = pop_deleterpset( ALLERP, index);
%
% EXAMPLE  :
%
% >> ALLERP = pop_deleterpset( ALLERP, [3 5]);
%
% INPUTS   :
%
% ALLERP    - Includes all ERPsets in workspace
% Index     - ERPset(s) that you want to clear from the workspace
%
% OUTPUTS  :
%
% - updated (output) ALLERP. Will include all ERPsets, minus that deleted
%
% *** This function is part of ERPLAB Toolbox ***
% Author: Javier Lopez-Calderon
% Center for Mind and Brain
% University of California, Davis,
% Davis, CA
% 2009

function [ALLERP, erpcom] = pop_deleterpset(ALLERP, varargin)
erpcom = '';
if nargin<1
        help pop_deleterpset
end
if nargin==1
        try
                CURRENTERP = evalin('base', 'CURRENTERP');
        catch
                CURRENTERP = 0;
        end
        if  CURRENTERP == 0
                msgboxText =  'ERPsets menu is already empty...';
                title      =  'ERPLAB: no erpset(s)';
                errorfound(msgboxText, title);
                return
        end

        prompt    = {'Erpset(s) to clear:'};
        dlg_title = 'Delete erpset(s)';
        num_lines = 1;
        def = {num2str(CURRENTERP)}; %01-13-2009

        %
        % open window
        %
        answer = inputvalue(prompt,dlg_title,num_lines,def);

        if isempty(answer)
                disp('User selected Cancel')
                return
        end

        erpindex =  str2num(answer{1});
        nerpset  = length(ALLERP);
        erpindex = unique_bc2(erpindex);

        if isempty(erpindex)
                msgboxText =  'Wrong erpset index(es)';
                title      =  'ERPLAB: unrecognizable erpset(s)';
                errorfound(msgboxText, title);
                return
        end
        if max(erpindex)>nerpset || min(erpindex)<1
                erpm     = findobj('tag', 'linerp');
                nerpmenu = length(erpm);
                if max(erpindex)<=nerpmenu && nerpmenu>=1 && max(erpindex)>=1
                        %...
                else
                        msgboxText = ['Wrong erpset index(es)\n'...
                                'Check your erpset menu or write length(ALLERP) at command window for comprobation'];
                        title        =  'ERPLAB: pop_deleterpset not existing erpset(s)';
                        errorfound(sprintf(msgboxText), title);
                        return
                end
        end

        %
        % Somersault
        %
        [ALLERP, erpcom] = pop_deleterpset(ALLERP, 'Erpsets', erpindex, 'Saveas', 'on','History', 'gui');
        return
end

%
% Parsing inputs
%
p = inputParser;
p.FunctionName  = mfilename;
p.CaseSensitive = false;
p.addRequired('ALLERP');
% option(s)
p.addParamValue('Erpsets', 1); % erpset index or input file
p.addParamValue('Warning', 'off', @ischar); % history from scripting
p.addParamValue('Saveas', 'off', @ischar); % 'on', 'off'
p.addParamValue('History', 'script', @ischar); % history from scripting
p.parse(ALLERP, varargin{:});

erpindex = p.Results.Erpsets;
if strcmpi(p.Results.Warning,'on')
        warnop = 1;
else
        warnop = 0;
end
if ismember_bc2({p.Results.Saveas}, {'on','yes'})
        issaveas  = 1;
else
        issaveas  = 0;
end
if strcmpi(p.Results.History,'implicit')
        shist = 3; % implicit
elseif strcmpi(p.Results.History,'script')
        shist = 2; % script
elseif strcmpi(p.Results.History,'gui')
        shist = 1; % gui
else
        shist = 0; % off
end
nerpset = length(ALLERP);
erpindex   = unique_bc2(erpindex);

if isempty(erpindex)
        msgboxText =  'Wrong erpset index(es)';
        error(['ERPLAB says: ' msgboxText])
end
if max(erpindex)>nerpset || min(erpindex)<1
        msgboxText = 'Wrong erpset index(es)';
        error(['ERPLAB says: ' msgboxText])
end
detect   = ~ismember_bc2(1:nerpset,erpindex);
newindex = find(detect);
if isempty(newindex)
        ALLERP = [];
else
        ALLERP = ALLERP(newindex);
end
if issaveas
        updatemenuerp(ALLERP, -1)
        assignin('base','ALLERP',ALLERP);  % save to workspace. Dec 5, 2012
end
% erpcom = sprintf('ALLERP = pop_deleterpset( ALLERP, [%s]);', num2str(erpindex));
%
% History
%
skipfields = {'ALLERP', 'History'};
fn     = fieldnames(p.Results);
erpcom = sprintf( '%s = pop_deleterpset( %s ', inputname(1), inputname(1) );
for q=1:length(fn)
        fn2com = fn{q};
        if ~ismember_bc2(fn2com, skipfields)
                fn2res = p.Results.(fn2com);
                if ~isempty(fn2res)
                        if ischar(fn2res)
                                if ~strcmpi(fn2res,'off')
                                        erpcom = sprintf( '%s, ''%s'', ''%s''', erpcom, fn2com, fn2res);
                                end
                        else
                                if iscell(fn2res)
                                        if ischar([fn2res{:}])
                                                fn2resstr = sprintf('''%s'' ', fn2res{:});
                                        else
                                                fn2resstr = vect2colon(cell2mat(fn2res), 'Sort','on');
                                        end
                                        fnformat = '{%s}';
                                else
                                        fn2resstr = vect2colon(fn2res, 'Sort','on');
                                        fnformat = '%s';
                                end
                                if strcmpi(fn2com,'Criterion')
                                        if p.Results.Criterion<100
                                                erpcom = sprintf( ['%s, ''%s'', ' fnformat], erpcom, fn2com, fn2resstr);
                                        end
                                else
                                        erpcom = sprintf( ['%s, ''%s'', ' fnformat], erpcom, fn2com, fn2resstr);
                                end
                        end
                end
        end
end
erpcom = sprintf( '%s );', erpcom);

% get history from script. ERP
switch shist
        case 1 % from GUI
                displayEquiComERP(erpcom);
        case 2 % from script
                for i=1:length(ALLERP)
                        ALLERP(i) = erphistory(ALLERP(i), [], erpcom, 1);
                end
        case 3
                % implicit

                %for i=1:length(ALLERP)
                %        ALLERP(i) = erphistory(ALLERP(i), [], erpcom, 1);
                %end
                %fprintf('%%Equivalent command:\n%s\n\n', erpcom);
        otherwise %off or none
                erpcom = '';
end

%
% Completion statement
%
msg2end
eeglab redraw
return

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
%{
 * Copyright (C) 2020-2030, The Regents of The University of Michigan.
 * All rights reserved.
 * This software was developed in the Biped Lab (https://www.biped.solutions/)
 * under the direction of Jessy Grizzle, grizzle@umich.edu. This software may
 * be available under alternative licensing terms; contact the address above.
 * 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 OWNER 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.
 * The views and conclusions contained in the software and documentation are those
 * of the authors and should not be interpreted as representing official policies,
 * either expressed or implied, of the Regents of The University of Michigan.
 *
 * AUTHOR: Bruce JK Huang (bjhuang[at]umich.edu)
 * WEBSITE: https://www.brucerobot.com/
%}

function edge = findEdgePointsInIdealFrame(H_TL, X_ref, tag_length)
    % X_ref: points that have been transformed back to the ideal frame
    % X_ref(1,:) : x
    % X_ref(2,:) : y
    % X_ref(3,:) : z
    % X_ref(4,:) : I
    % X_ref(5,:) : R
    % X_ref(6,:) : scan
    % X_ref(7,:) : ExpNmbr
    top_ring = max(X_ref(5, :));
    bottom_ring = min(X_ref(5, :));
    num_scan = max(X_ref(6, :)) - min(X_ref(6, :));
    pick_points = max(1, ceil(num_scan));
    % Ideal frame (top-left, bottom-left, bottom-right, top-right)
%     ideal_frame.x = zeros(1,8);
%     ideal_frame.y = [-tag_length/2, -tag_length/2, -tag_length/2, tag_length/2, tag_length/2, tag_length/2, tag_length/2, -tag_length/2];
%     ideal_frame.z = [tag_length/2, -tag_length/2, -tag_length/2, -tag_length/2, -tag_length/2, tag_length/2, tag_length/2, tag_length/2];

    % left, bottom, right, top (y-axis, -z-axis, -y-axis, z-axis)
    edge_list(1).line = [-tag_length/2 0; -tag_length/2 1];
    edge_list(2).line = [0 -tag_length/2; 1 -tag_length/2];
    edge_list(3).line = [tag_length/2 0; tag_length/2 1];
    edge_list(4).line = [0 tag_length/2; 1 tag_length/2];
    edge(4).points = [];

    % choose one point per scan per ring
%     for j = 1:num_scan
%         current_scan_points = X_ref(:, X_ref(6,:)==j);
%         for i = bottom_ring:top_ring
%             if size(current_scan_points(2:3, current_scan_points(5,:)==i), 2) > 1
%                 current_ring_points = current_scan_points(1:3, current_scan_points(5,:)==i);
%                 extreme_points = [current_ring_points(:, 1), current_ring_points(:, end)];
%                 extreme_point1 = current_ring_points(:, 1);
%                 extreme_point2 = current_ring_points(:, end);
%                 dist1 = zeros(1, 4);
%                 dist2 = zeros(1, 4);
%                 for k = 1:4
%                     dist1(1, k) = pointToLineDistance(extreme_point1(2:3, :)', edge_list(k).line(1, :), edge_list(k).line(2, :));
%                     dist2(1, k) = pointToLineDistance(extreme_point2(2:3, :)', edge_list(k).line(1, :), edge_list(k).line(2, :));
%                 end
%                 [~, edge1_i] = min(dist1);
%                 [~, edge2_i] = min(dist2);
%                 edge_point_3D_1 = inv(H_TL) * [extreme_point1; 1];
%                 edge_point_3D_2 = inv(H_TL) * [extreme_point2; 1];
% %                 selected_edge_points = [current_ring_points(1:3, indices(selected_indices)); ones(1, length(selected_indices))];
%                 edge(edge1_i).points = [edge(edge1_i).points, edge_point_3D_1];
%                 edge(edge2_i).points = [edge(edge2_i).points, edge_point_3D_2];
%             end
%         end
%     end

     % choose a few points using all scans
    for i = bottom_ring:top_ring
        current_ring_points = X_ref(:, X_ref(5,:)==i);
%         dist = zeros(length(current_ring_points), 4);
        for k = 1:4
            dist = pointToLineDistance(current_ring_points(2:3,:)', edge_list(k).line(1, :), edge_list(k).line(2, :));
            [value, indices] = mink(dist, pick_points);
            selected_indices = find (value<tag_length/8);
            selected_edge_points = inv(H_TL) * [current_ring_points(1:3, indices(selected_indices)); ones(1, length(selected_indices))];
%             selected_edge_points = [current_ring_points(1:3, indices(selected_indices)); ones(1, length(selected_indices))];
            edge(k).points = [edge(k).points, selected_edge_points];
        end
    end

%     edge_list = [tag_length/2, -tag_length/2, -tag_length/2, tag_length/2];
%     edge(4).points = [];
%     for i = bottom_ring:top_ring
%         current_ring_points = X_ref(:, X_ref(5,:)==i);
%         first_point = current_ring_points(:, 1);
%         edge_distance = [abs(first_point(2) - edge_list(1)), ...
%                          abs(first_point(3) - edge_list(2)), ...
%                          abs(first_point(2) - edge_list(3)), ...
%                          abs(first_point(3) - edge_list(4))];
%         [~, first_index] = min(edge_distance);
%         first_point = inv(H_TL) * [first_point(1:3); 1];
%         edge(first_index).points = [edge(first_index).points, first_point];
%         last_point = current_ring_points(:, end);
%         edge_distance = [abs(last_point(2) - edge_list(1)), ...
%                          abs(last_point(3) - edge_list(2)), ...
%                          abs(last_point(2) - edge_list(3)), ...
%                          abs(last_point(3) - edge_list(4))];
%         [~, last_index] = min(edge_distance);
%         last_point = inv(H_TL) * [last_point(1:3); 1];
%         edge(last_index).points = [edge(last_index).points, last_point];
%     end
%     fig_hangle = figure(99999);
%     cla(fig_hangle)
%     X_move_back = inv(H_TL) * [X_ref(1:3,:); ones(1, size(X_ref(1, :), 2))];
%     scatter3(X_move_back(1,:), X_move_back(2,:), X_move_back(3,:), 'k.')
% %     scatter3(X_ref(1,:), X_ref(2,:), X_ref(3,:), 'k.')
%     hold on
%     scatter3(edge(1).points(1,:), edge(1).points(2,:), edge(1).points(3,:), 'ro', 'filled')
%     scatter3(edge(2).points(1,:), edge(2).points(2,:), edge(2).points(3,:), 'go', 'filled')
%     scatter3(edge(3).points(1,:), edge(3).points(2,:), edge(3).points(3,:), 'bo', 'filled')
%     scatter3(edge(4).points(1,:), edge(4).points(2,:), edge(4).points(3,:), 'mo', 'filled')
%     title('edge parsing results')
%     xlabel("x")
%     ylabel("y")
%     zlabel("z")
%     axis equal
end

author: UMich-BipedLab


text:
function [C,group_sizes]=contour_group(varargin)

%function [C,group_sizes]=contour_group(X,Y,Z,M,Sx,Sy,Sz,c,interpMethod)

%% Parse input

switch nargin
    case 8
        X=varargin{1};
        Y=varargin{2};
        Z=varargin{3};
        M=varargin{4};
        Sx=varargin{5};
        Sy=varargin{6};
        Sz=varargin{7};
        c=varargin{8};
        interpMethod='linear';
    case 9
        X=varargin{1};
        Y=varargin{2};
        Z=varargin{3};
        M=varargin{4};
        Sx=varargin{5};
        Sy=varargin{6};
        Sz=varargin{7};
        c=varargin{8};
        interpMethod=varargin{9};
end
%% Create contours
% Open figure window with visibility off

ht=figure('Visible','off');
if ndims(M)==3

else %2D array
    X(:,:,2)=X; Y(:,:,2)=Y; Z(:,:,2)=Z+1; M(:,:,2)=M;
end

%Contourslice is used since it provides children in handles
h=contourslice(X,Y,Z,M,Sx,Sy,Sz,[c c],interpMethod); %This will plot in invisible window

%% Parse children
%

C=cell(1,numel(h));
for i=1:length(h)
    %Get coordinates
    x=get(h(i),'XData'); y=get(h(i),'YData'); z=get(h(i),'ZData');

    %Vertices matrix
    V=[x(:) y(:) z(:)];
    V=V(all(~isnan(V),2),:); %Removing NaN's

    %Store in cell array
    C{i}=V;
end

%Get group sizes
group_sizes = cell2mat(cellfun(@(x) size(x,1), C,'UniformOutput',0)');

%Remove contours of lenght 1
C=C(group_sizes>1);
group_sizes=group_sizes(group_sizes>1);

%Sort C array according to group size
[group_sizes,ind_sort]=sort(group_sizes);
C=C(ind_sort);

close(ht);

end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


function results = OTB_DEEP_settings(seq, res_path, bSaveImage, parameters)

close all

s_frames = seq.s_frames;

% Feature specific parameters
hog_params.cell_size = 4;
hog_params.compressed_dim = 10;

% grayscale_params.colorspace='gray';
% grayscale_params.cell_size = 1;
%
% cn_params.tablename = 'CNnorm';
% cn_params.useForGray = false;
% cn_params.cell_size = 4;
% cn_params.compressed_dim = 3;
%
% ic_params.tablename = 'intensityChannelNorm6';
% ic_params.useForColor = false;
% ic_params.cell_size = 4;
% ic_params.compressed_dim = 3;

cnn_params.nn_name = 'imagenet-vgg-m-2048.mat'; % Name of the network
cnn_params.output_layer = [3 14];               % Which layers to use
cnn_params.downsample_factor = [2 1];           % How much to downsample each output layer
cnn_params.compressed_dim = [16 64];            % Compressed dimensionality of each output layer
cnn_params.input_size_mode = 'adaptive';        % How to choose the sample size
cnn_params.input_size_scale = 1;                % Extra scale factor of the input samples to the network (1 is no scaling)

% Which features to include
params.t_features = {
    struct('getFeature',@get_cnn_layers, 'fparams',cnn_params),...
    struct('getFeature',@get_fhog,'fparams',hog_params),...
    ...struct('getFeature',@get_colorspace, 'fparams',grayscale_params),...
    ...struct('getFeature',@get_table_feature, 'fparams',cn_params),...
    ...struct('getFeature',@get_table_feature, 'fparams',ic_params),...
};

% Global feature parameters1s
params.t_global.normalize_power = 2;    % Lp normalization with this p
params.t_global.normalize_size = true;  % Also normalize with respect to the spatial size of the feature
params.t_global.normalize_dim = true;   % Also normalize with respect to the dimensionality of the feature

% Image sample parameters
params.search_area_shape = 'square';    % The shape of the samples
params.search_area_scale = 4.5;         % The scaling of the target size to get the search area
params.min_image_sample_size = 200^2;   % Minimum area of image samples
params.max_image_sample_size = 300^2;   % Maximum area of image samples

% Detection parameters
params.refinement_iterations = 1;       % Number of iterations used to refine the resulting position in a frame
params.newton_iterations = 5;           % The number of Newton iterations used for optimizing the detection score
params.clamp_position = false;          % Clamp the target position to be inside the image

% Learning parameters
params.output_sigma_factor = 1/12;		% Label function sigma
params.learning_rate = 0.011;	 	 	% Learning rate
params.nSamples = 50;                   % Maximum number of stored training samples
params.sample_replace_strategy = 'lowest_prior';    % Which sample to replace when the memory is full
params.lt_size = 0;                     % The size of the long-term memory (where all samples have equal weight)
params.train_gap = 5;                   % The number of intermediate frames with no training (0 corresponds to training every frame)
params.skip_after_frame = 1;            % After which frame number the sparse update scheme should start (1 is directly)
params.use_detection_sample = true;     % Use the sample that was extracted at the detection stage also for learning

% Factorized convolution parameters
params.use_projection_matrix = true;    % Use projection matrix, i.e. use the factorized convolution formulation
params.update_projection_matrix = true; % Whether the projection matrix should be optimized or not
params.proj_init_method = 'pca';        % Method for initializing the projection matrix
params.projection_reg = 2e-8;	 	 	% Regularization paremeter of the projection matrix

% Generative sample space model parameters
params.use_sample_merge = true;                 % Use the generative sample space model to merge samples
params.sample_update_criteria = 'Merge';        % Strategy for updating the samples
params.weight_update_criteria = 'WeightedAdd';  % Strategy for updating the distance matrix
params.neglect_higher_frequency = false;        % Neglect hiigher frequency components in the distance comparison for speed

% Conjugate Gradient parameters
params.CG_iter = 5;                     % The number of Conjugate Gradient iterations in each update after the first frame
params.init_CG_iter = 10*20;            % The total number of Conjugate Gradient iterations used in the first frame
params.init_GN_iter = 10;               % The number of Gauss-Newton iterations used in the first frame (only if the projection matrix is updated)
params.CG_use_FR = false;               % Use the Fletcher-Reeves (true) or Polak-Ribiere (false) formula in the Conjugate Gradient
params.CG_standard_alpha = true;        % Use the standard formula for computing the step length in Conjugate Gradient
params.CG_forgetting_rate = 80;	 	 	% Forgetting rate of the last conjugate direction
params.precond_data_param = 0.3;	 	% Weight of the data term in the preconditioner
params.precond_reg_param = 0.02;	 	% Weight of the regularization term in the preconditioner
params.precond_proj_param = 75;	 	 	% Weight of the projection matrix part in the preconditioner

% Regularization window parameters
params.use_reg_window = true;           % Use spatial regularization or not
params.reg_window_min = 1e-4;			% The minimum value of the regularization window
params.reg_window_edge = 10e-3;         % The impact of the spatial regularization
params.reg_window_power = 2;            % The degree of the polynomial to use (e.g. 2 is a quadratic window)
params.reg_sparsity_threshold = 0.05;   % A relative threshold of which DFT coefficients that should be set to zero

% Interpolation parameters
params.interpolation_method = 'bicubic';    % The kind of interpolation kernel
params.interpolation_bicubic_a = -0.75;     % The parameter for the bicubic interpolation kernel
params.interpolation_centering = true;      % Center the kernel at the feature sample
params.interpolation_windowing = false;     % Do additional windowing on the Fourier coefficients of the kernel

% Scale parameters for the translation model
% Only used if: params.use_scale_filter = false
params.number_of_scales = 5;            % Number of scales to run the detector
params.scale_step = 1.02;               % The scale factor

% Scale filter parameters
% Only used if: params.use_scale_filter = true
params.use_scale_filter = false;        % Use the fDSST scale filter or not (for speed)
% params.scale_sigma_factor = 1/16;       % Scale label function sigma
% params.scale_learning_rate = 0.025;		% Scale filter learning rate
% params.number_of_scales_filter = 17;    % Number of scales
% params.number_of_interp_scales = 33;    % Number of interpolated scales
% params.scale_model_factor = 1.0;        % Scaling of the scale model
% params.scale_step_filter = 1.02;        % The scale factor for the scale filter
% params.scale_model_max_area = 32*16;    % Maximume area for the scale sample patch
% params.scale_feature = 'HOG4';          % Features for the scale filter (only HOG4 supported)
% params.s_num_compressed_dim = 'MAX';    % Number of compressed feature dimensions in the scale filter
% params.lambda = 1e-2;					% Scale filter regularization
% params.do_poly_interp = true;           % Do 2nd order polynomial interpolation to obtain more accurate scale

% Other parameters
params.visualization = 0;               % Visualiza tracking and detection scores
params.debug = 0;                       % Do full debug visualization

% Initialize
params.init_sz = [seq.init_rect(1,4), seq.init_rect(1,3)];
params.init_pos = [seq.init_rect(1,2), seq.init_rect(1,1)] + (params.init_sz - 1)/2;
params.s_frames = s_frames;

% Run tracker
results = tracker(params);


text:
function plotCPs(CPs)
% DESCRIPTION
% Plot Contribution Plots
%
%    plotCPs(CPs)
%
% INPUT
%   CPs         Contribution Plots
%
% Created on 18th April 2019, by Kepeng Qiu.
%-------------------------------------------------------------%

%
figure
temp = mean(abs(CPs),1);
bar(temp/sum(temp,2));

% Axis settings
tgca = 16;  % font size
tfont = 'Helvetica'; % font type
% tfont = 'Arial'; % font type
% set(gca,'yscale','log')
set(gca,'FontSize',tgca,'FontName',tfont)

% legend settings
tlegend = tgca*0.9;
legend({'Contribution Plots'},'FontSize',tlegend , ...
    'FontWeight','normal','FontName',tfont)

% label settings
tlabel = tgca*1.1;
xlabel('Variable','FontSize',tlabel,'FontWeight','normal', ...
    'FontName',tfont,'Color','k')
ylabel('CPs','FontSize',tlabel,'FontWeight','normal', ...
    'FontName',tfont,'Color','k')

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function snow = DetectSnow(dim,band_green,band_nir,band_bt,ndsi)
%DETECSNOW Summary of this function goes here
%   Detailed explanation goes here

% %     snow=zeros(dim,'uint8'); % Snow mask
    % It takes every snow pixels including snow pixel under thin clouds or icy clouds
    snow=ndsi>0.15&band_nir>1100&band_green>1000;
    if ~isempty(band_bt)
        snow=snow&band_bt<1000;
    end
%     snow(ids_snow)=1;
end

author: GERSL

repo: Fmask

save_path: github-repos/MATLAB/GERSL-Fmask


text:
function x = line_grid ( n, a, b, c )

%*****************************************************************************80
%
%% LINE_GRID: grid points over the interior of a line segment in 1D.
%
%  Discussion:
%
%    In 1D, a grid is created using N points.
%
%    Over the interval [A,B], we have 5 choices for grid centering:
%      1: 0,   1/3, 2/3, 1
%      2: 1/5, 2/5, 3/5, 4/5
%      3: 0,   1/4, 2/4, 3/4
%      4: 1/4, 2/4, 3/4, 1
%      5: 1/8, 3/8, 5/8, 7/8
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 August 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points.
%
%    Input, real A, B, the endpoints for each dimension.
%
%    Input, integer C, the grid centering for each dimension.
%    1 <= C <= 5.
%
%    Output, real X(N), the points.
%

%
%  Create the 1D grids in each dimension.
%
  x = zeros ( n, 1 );

  for j = 1 : n

    if ( c == 1 )

      if ( n == 1 )
        x(j) = 0.5 * ( a + b );
      else
        x(j) = (   ( n - j     ) * a   ...
                 + (     j - 1 ) * b ) ...
                 / ( n     - 1 );
      end
    elseif ( c == 2 )
      x(j) = (   ( n - j + 1 ) * a   ...
               + (     j     ) * b ) ...
               / ( n     + 1 );
    elseif ( c == 3 )
      x(j) = (   ( n - j + 1 ) * a   ...
               + (     j - 1 ) * b ) ...
               / ( n         );
    elseif ( c == 4 )
      x(j) = (   ( n - j ) * a   ...
               + (     j ) * b ) ...
               / ( n     );
    elseif ( c == 5 )
      x(j) = (   ( 2 * n - 2 * j + 1 ) * a   ...
               + (         2 * j - 1 ) * b ) ...
               / ( 2 * n             );
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function lij = fem_basis_1d ( i, j, x )

%*****************************************************************************80
%
%% FEM_BASIS_1D evaluates an arbitrary 1D basis function.
%
%  Discussion:
%
%    Given the maximum degree D for the polynomial basis defined
%    on a reference interval, we have D + 1 monomials
%    of degree at most D.  In each barycentric coordinate, we define
%    D+1 points, so that 0 <= I, J <= D and I+J = D, with
%    (I,J) corresponding to
%    * the basis point X(I,J) = ( I/D );
%    * the basis monomial P(I,J)(X) = X^I.
%
%    For example, with D = 2, we have simply:
%
%      A---B---C
%
%    with
%
%       I J    X      P(I,J)(X)
%
%    A (0 2) ( 0.0 )  1
%    B (1 1) ( 0.5 )  x
%    C (2 0) ( 1.0 )  x^2
%
%    Now instead of the monomials P(I,J)(X), we want a set of
%    polynomials L(I,J)(X) which span the same space, but have
%    the Lagrange property, namely L(I,J) (X) is 1 if X is
%    equal to X(I,J), and 0 if X is equal to any other
%    of the basis points.
%
%    This is easily arranged.  Given an index (I,J), we compute
%    1) I factors of the form (X-0)   * (X-1/D)   * ... * (X-(I-1)/D);
%    2) J factors of the form ???
%
%    This results in the product of I+J linear factors, in other words,
%    a polynomial of degree D.  This polynomial is 0 at all basis points
%    except X(I,J).  If we divide this polynomial by its value at
%    the basis point, we arrive at the desired Lagrange polynomial
%    L(I,J)(X).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 October 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, J, the integer barycentric coordinates of the basis
%    function, 0 <= I, J.  The polynomial degree D = I + J
%
%    Input, real X, the evaluation point.
%
%    Output, real LIJ, the value of the basis function at X.
%
  d = i + j;
  lij = 1.0;
  cij = 1.0;
  for p = 0 : i - 1
    lij = lij * ( d * x - p );
    cij = cij * (     i - p );
  end
  for p = 0 : j - 1
    lij = lij * ( d * x - ( d - p ) );
    cij = cij * (     i - ( d - p ) );
  end

  lij = lij / cij;

  return
end


text:
classdef SingleTargetMeasurementSimulatorX < MeasurementSimulatorX
% SingleTargetMeasurementSimulatorX class
%
% Summary of SingleTargetMeasurementSimulatorX:
%   Class implementation of a single-target measurement simulator
%
% SingleTargetMeasurementSimulatorX Properties:
%   + Model - A StateSpaceModelX object defining the base models used by the simulator
%
% SingleTargetMeasurementSimulatorX Methods:
%   + SingleTargetMeasurementSimulatorX - Constructor method
%   + simulate() - Generate measurements from ground-truth data
%
% (+) denotes puplic properties/methods
%
% See also MeasurementModelX and ClutterModelX template classes

    methods
        function this = SingleTargetMeasurementSimulatorX(varargin)
        % SingleTargetMeasurementSimulatorX Construct a measurement simulator object
        %
        % Parameters
        % ----------
        % Model: StateSpaceModelX object
        %   A state-space model that should define the following models:
        %       + Measurement - Object handle to MeasurementModelX SubClass
        %       + Clutter - Object handle to ClutterModelX SubClass
        %       + Detection - Object handle to DetectionModelX SubClass

            % Call super-class
            this@MeasurementSimulatorX(varargin{:});

        end

        function MeasurementScans = simulate(this, Track, varargin)
        % simulate Simulate measurements given a set of ground truth data
        %
        % Parameters
        % ----------
        % Track: GroundTruthTrackX
        %   An object containing ground-truth information for a single-target.
        %   Detections will be generated on the basis of the Trajectory
        %   property for each track.
        %
        % Returns
        % -------
        % MeasurementScans: (1 x NumScans) MeasurementListX array
        %   An object array, whose elements contain the simulated measurement
        %   scan at each timestep.

            % Dummy Detection probability
            MeasurementScans = MeasurementListX.empty();

            % Initialise storage
            numTimesteps = numel(Track.Trajectory);

            for k = 1:numTimesteps

                targetState = Track.Trajectory(k);
                timestamp = targetState.Timestamp;

                % Compute number of tracks and clutter measurements
                detectionProbability = 1;
                if ~isempty(this.Model.Clutter)
                    detectionProbability = this.Model.Detection.pdf(targetState.Vector);
                end

                targetDetected = binornd(1,detectionProbability);
                numClutter = 0;
                if ~isempty(this.Model.Clutter)
                    numClutter = this.Model.Clutter.random('cardinality');
                end
                numMeasurements = targetDetected + numClutter;

                % Empty MeasurementX array
                measurements = MeasurementX.empty(0, numMeasurements);

                % Generate true measurement
                if targetDetected
                    measurementVector = this.Model.Measurement.feval(targetState.Vector,true);
                    measurements(end+1) = MeasurementX(measurementVector, timestamp, this.Model);
                end

                % Generate clutter measurements
                if(numClutter>0)
                    clutterVectors = this.Model.Clutter.random(numClutter);
                    for clutterVector = clutterVectors
                        measurements(end+1) = MeasurementX(clutterVector, timestamp, this.Model);
                    end
                end

                MeasurementScans(k) = MeasurementListX(measurements);
            end

        end
    end
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function display(N)
%DISPLAY   Print to command line for CHEBOP2.
% DISPLAY is called automatically when a statement that results in a
% CHEBOP2 output is not terminated with a semicolon.

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

loose = ~isequal(get(0, 'FormatSpacing'), 'compact');

if ( loose )
    fprintf('\n')
end

% Compact version.
if ( isempty(N) )
    fprintf('    empty chebop2\n\n')
    return
end

% Get information that we want to display.
dom = N.domain;
op = N.op;

% Display the information.
fprintf('chebop2 of rank %u: (linear bivariate operator)\n',rank(N))
fprintf('         domain               operator  \n');

fprintf(['[%3.2g,%3.2g] x [%3.2g,%3.2g]  ' func2str(op)  '\n'], dom(1),...
                                          dom(2), dom(3), dom(4));
printBoundaryConditions(N.lbc, 'left');
printBoundaryConditions(N.rbc, 'right');
printBoundaryConditions(N.ubc, 'top');
printBoundaryConditions(N.dbc, 'bottom');

if ( loose )
    fprintf('\n');
end

end

function printBoundaryConditions(BC, string)
% Display the boundary conditions.

    exists = ~isempty(BC);
    if ( exists )
        if ( isa(BC, 'double') )
            fprintf([string ' boundary = %u'], BC);
        elseif ( isa(BC, 'chebfun') )
            % Try to convert to a double:
            n = length(BC);
            if ( n == 1 )
                val = BC.values;
                printBoundaryConditions(val, string)
                return
            end
            fprintf([string ' boundary = function']);
        elseif ( isa(BC, 'function_handle') )
            fprintf([string ' boundary = ' func2str(BC)]);
        end
        fprintf('\n');
    end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun