text:
% function intersect = dtiFiberIntersectMesh(fibers, triangles, vertices)
%
% Given a set of 'fibers' and a surface mesh (defined by
% triangles/vertices), will find every triangle in the surface mesh
% that is intersected by the fibers.
%
% ARGUMENTS:
%
% fibers = an Nx1 cell array of fibers, where each cell contains a
%          3xN real array of points (XYZ order) that specify a fiber
%          path. A fiber must contain more than 1 point (or be empty).
%
% triangles = a 3xN uint32 array of triangles, where each entry is
%             an index into the vertices array.
%
% vertices = a 3xN double array of vertices in YXZ order.
%
% NOTE! We assume that the mesh vertices are X-Y swapped relative to
% the fiber vertices. (Because they are in our data representations.)
%
% RETURNS:
%
% intersect = an Nx1 cell array, where N = the number of fibers
%             (ie. length(fibers)). Each cell contains a 5xN array
%             where N is the number of intersections. It has the
%             following structure: [triangleIndex  fiberIndex X Y Z]
%
% NOTES:
% The hard work is done by the RAPID collision detection library. This
% is only free for non-commercial use, so you'll need to get the library
% from its maintainer if you want to rebuild the mex file. See:
% http://www.cs.unc.edu/~geom/OBB/OBBT.html or search google for
% "Robust and Accurate Polygon Interference Detection".
%
% To compile on Linux:
%    mex -O -I. dtiFiberIntersectMesh.cxx libRAPID.a
%
% on Windows:
%    mex -O -I. dtiFiberIntersectMesh.cxx RAPID.lib
%
% To make Rapid.lib on Windows, there is a directory
%    ..\mrDiffusion\src\RAPID_VStudio_201\MSVC_Compile
% with a Visual Studio Solution file.
%
% HISTORY:
% 2004.07.28 Bob Dougherty: wrote it.

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function varargout = add_lights()
  % ADD_LIGHTS  Add 5 lights  based on the current camera
  %
  %ls = add_lights()
  %
  % Outputs:
  %   ls  #ls list of light handles
  %
lights = @() { ...
light('Color',1.0*[1 1 1],'Position',(campos-camtarget),'Style','local'), ...
light('Color',0.5*teal,'Position', (campos-camtarget)*axisangle2matrix([0 0 1],pi/2),'Style','local'), ...
light('Color',0.5*teal,'Position',-(campos-camtarget)*axisangle2matrix([0 0 1],pi/2),'Style','local'), ...
light('Color',0.5*pink,'Position', (campos-camtarget)*axisangle2matrix([1 0 0],pi*0.9),'Style','local'), ...
light('Color',0.5*pink,'Position', (campos-camtarget)*axisangle2matrix([1 0 0],-pi*0.9),'Style','local') ...
};
  ls = lights();
  if nargout >= 1
    varargout{1} = ls;
  end
end

author: alecjacobson

repo: gptoolbox


text:
close all;
clear all;
clc;

data_file_path = 'bin/mmv_phase_transition_noiseless_s_2.mat';
options.export = true;
options.export_dir = 'bin';
options.export_name = 'mmv_noiseless_s_2';
options.chosen_ks = [2, 4, 8, 16, 32, 64];
options.subtitle = 'Noiseless, S = 2';
spx.pursuit.PhaseTransitionAnalysis.print_results(data_file_path, ...
    'CoSaMP', options);

author: indigits

repo: sparse-plex


text:
function [imax, xfin, s2fin, ufin, Cxx, uout] = OLOO(A, y, Q)
% SYNTAX:
%   [imax, xfin, s2fin, ufin, Cxx, uout] = OLOO(A, y, Q)
%
% INPUT:
%   A: design matrix
%   y: observations vector
%   Q: cofactor matrix
%
% OUTPUT:
%   imax:  index of the rejected blocks
%   x_fin: estimated parameters without outlier
%   s2fin: a posteriori sigma without outlier
%   ufin:  estimated residuals without outlier
%   Cxx:   parameters covariance of the final solution
%   uout:  residual of outlier observation
%
% DESCRIPTION:
%   perform LS on blocks of correlated observations
%   identify one (block) outlier
%   reject it
%   re-estimate unknowns
%   according to the theory in "L. Biagi and S. Caldera. An efficient leave one block out approach to identify outliers.Journal of Applied Geodesy, Volume 7, Issue 1, pages 11..19, 2013"
%
%   this version is optimized to manage l.o.o. of 1 observation at time, no blocks!
%
% CREDITS:
%   1.0: Stefano Caldera, 22.05.2014
%   1.1: Stefano Caldera, Andrea Gatti 08.12.2016 ( speedup improvements )

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:       Stefano Caldera 22.05.2014
%  Contributors:     Stefano Caldera,
%                    Andrea Gatti 08.12.2016 ( speedup improvements )
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

% when FTABLE is undefined, redefine it
global FTABLE; if isempty(FTABLE); FTABLE = finv(0.9995,1,1:size(A,1))'; end
if (size(FTABLE,1) < size(A,1)); FTABLE = finv(0.9995,1,1:size(A,1))'; end
n_blocks = length(y);
[m, n] = size(A);           % m: number of observations, n: number of unknowns
uout = NaN;

if m - n > 0

    % compute the global solution
    if isdiag(Q) % Q is tipically diagonal -> let's use this information
        Qm = diag(Q);
        Qm = Qm ./ min(Qm);
        invQ = diag(1./Qm);
        At_invQ = A' * invQ;
        Ninv = inv(At_invQ * A);
        xcap = Ninv * At_invQ * y;
        um = y - A * xcap;
        s2cap = um' * invQ * um/(m-n);
    else
        Q = Q ./ (min(diag(Q)));
        At_invQ = A'/Q;
        Ninv = inv(At_invQ * A);
        xcap = Ninv * At_invQ * y;
        um = y - A * xcap;
        s2cap = um' / Q * um/(m-n);
        Qm = diag(Q);
    end
    % convert A in a sparse matrix, faster and lighter

    use_sparse_approach = false;

    if sum(A(1,:)==0) > size(A,2) / 0.577     % if A is sparse enough
        use_sparse_approach  = true;
        A = sparse(A);
    end
    if m - n > 1
        %% start outliers rejection
        Im = eye(n);
        Bm = Ninv * A';
        Cm = diag(A * Bm);
        Km = Qm - Cm;
        Kminv = Km.^-1;
        wm = Qm .* Kminv .* um;
        s2m = ((m-n) .* s2cap - um .* Kminv .* um) ./ (m - n - 1);

        % original loop
        if use_sparse_approach
            % modified loop to exploit the sparse property of the A matrix
            Qw = zeros(n_blocks,1);
            tmp = (A .* repmat(Kminv,1,n));
            for i = 1 : n_blocks
                idOk = (tmp(i,:) ~= 0);
                Qw(i) = Qm(i) + Bm(idOk,i)' * (Im(idOk,idOk) + tmp(i,idOk)' * Bm(idOk,i)')  * A(i,idOk)';
            end
        else
            Qw = zeros(n_blocks,1);
            for i = 1 : n_blocks
                Qw(i) = Qm(i) + Bm(:,i)' * (Im + A(i,:)' * Kminv(i) * Bm(:,i)') * A(i,:)';
                %Qw3=Bm(:,i)'*(Im+A(i,:)'*Kminv(i)*Bm(:,i)');
            end
        end

        %toc
        Qwinv = Qw.^-1;
        deg2 = m - n - 1;

        F = wm .* Qwinv .* wm ./ s2m;
        Flim = FTABLE(deg2);

        %% apply final solution
        % find maximum F(i)/Flim(i)
        [Fmax, imax] = max(abs(F ./ Flim));

        if (Fmax < 1)
            % no outlier
            imax = 0;
            xfin = xcap;
            s2fin = s2cap;
            Ninvfin = Ninv;
            ufin = um;
            Cxx = s2fin * Ninvfin;

        else
            % if the maximum ratio exceedes the threshold, the observation is eliminated from the solution
            uout = um(imax);
            xfin = xcap - Bm(:, imax) * Kminv(imax) * um(imax);
            yfin = y;
            yfin(imax) = [];

            Afin = A;
            Afin(imax,:) = [];

            s2fin = s2m(imax);

            Ninvfin = Ninv + Bm(:, imax) * Kminv(imax) * Bm(:, imax)';
            ufin = yfin - Afin * xfin;
            Cxx = s2fin * Ninvfin;

        end

    else
        imax = 0;
        xfin = xcap;
        s2fin = s2cap;
        Cxx = s2cap * Ninv;
        ufin = um;
    end

else
    % detection is not possibile
    imax = 0;
    xfin = [];
    Cxx = [];
    ufin = [];
    s2fin = NaN;
end

author: goGPS-Project


text:
classdef comp2periodz < ZmapHGridFunction
    % COMP2PERIODZ compares seismicity rates for two time periods
    % The differences are as z- and beta-values and as percent change.

    properties
        periodA_start   datetime % start time for period 1
        periodA_end     datetime % end time for period 1
        periodB_start   datetime % start time for period 2
        periodB_end     datetime % end time for period 2
        binsize duration            = ZmapGlobal.Data.bin_dur;
    end

    properties(Constant)
        PlotTag = 'comp2periodz';
        ReturnDetails = cell2table({ ... VariableNames, VariableDescriptions, VariableUnits
            ...
            ... % these are returned by the calculation function
            'z_value',              'z-value',          '';...           #1 'valueMap'
            'pct_change',           'percent change',   'pct';... #2  'per'
            'beta_value',           'Beta value map',   '';...     #3 'beta_map'
            'Number_of_Events_1',   'Number of events in first period', '';... #4
            'Number_of_Events_2',   'Number of events in second period', '';... #5
            ...
            }, 'VariableNames', {'Names','Descriptions','Units'})

        CalcFields={...
            'z_value',              'pct_change',   'beta_value',...
            'Number_of_Events_1',   'Number_of_Events_2'}

        ParameterableProperties = ["periodA_start" "periodA_end" "periodB_start" "periodB_end" "binsize"];

        References="";
    end

    methods
        function obj=comp2periodz(zap, varargin)
            % COMP2PERIODZ compares seismicity rates for two time periods.
            % The differences are as z- and beta-values and as percent change.
            %   Stefan Wiemer 1/95
            %   Rev. R.Z. 4/2001

            obj@ZmapHGridFunction(zap, 'z_value');
            report_this_filefun();

            t0b = min(obj.RawCatalog.Date);
            teb = max(obj.RawCatalog.Date);
            obj.periodA_start = t0b;
            obj.periodB_end = teb;
            obj.periodA_end = t0b + (teb-t0b)/2;
            obj.periodB_start = obj.periodA_end+minutes(0.01);

            obj.parseParameters(varargin);

            obj.StartProcess();

        end

        function InteractiveSetup(obj)

            % get two time periods, along with grid and event parameters
            zdlg=ZmapDialog();
            zdlg.AddHeader('Please define two time periods to compare');
            zdlg.AddEdit('periodA_start','start period 1',  obj.periodA_start,'start time for period 1');
            zdlg.AddEdit('periodA_end',  'end period 1',    obj.periodA_end,'end time for period 1');
            zdlg.AddEdit('periodB_start','start period 2',  obj.periodB_start,'start time for period 2');
            zdlg.AddEdit('periodB_end',  'end period 2',    obj.periodB_end,'end time for period 2');
            zdlg.AddDurationEdit('binsize','Bin Size',      obj.binsize,'number of days in each bin',@days);
            obj.AddDialogOption(zdlg,'EventSelector');
            zdlg.Create('Name', 'Please choose rate change estimation option','WriteToObj',obj,'OkFcn',@obj.doIt);
        end

        function results=Calculate(obj)

            assert(obj.periodA_start < obj.periodA_end,'Period 1 starts before it ends');
            assert(obj.periodB_start < obj.periodB_end,'Period 2 starts before it ends');

            %  make grid, calculate start- endtime etc.  ...

            lt =  (obj.RawCatalog.Date >= obj.periodA_start &  obj.RawCatalog.Date < obj.periodA_end) ...
                | (obj.RawCatalog.Date >= obj.periodB_start &  obj.RawCatalog.Date <= obj.periodB_end);
            obj.RawCatalog = obj.RawCatalog.subset(lt);

            interval1_bins = obj.periodA_start : obj.binsize : obj.periodA_end; % starts
            interval2_bins = obj.periodB_start : obj.binsize : obj.periodB_end; % starts
            interval1_edges = [interval1_bins, interval1_bins(end)+obj.binsize];
            interval2_edges = [interval2_bins, interval2_bins(end)+obj.binsize];

            obj.gridCalculations(@calculation_function);

            obj.Result.period1.dateRange=[obj.periodA_start obj.periodA_end];
            obj.Result.period2.dateRange=[obj.periodB_start obj.periodB_end];

            if nargout
                results=obj.Result.values;
            end
            % save data

            % plot the results
            % old and valueMap (initially ) is the z-value matrix

            %det =  'ast';
            %ZG.shading_style = 'interp';
            % View the b-value map: view_ratecomp.m
            %    which could create a topography overlay ala dramap_z.m
            %
            %

            function out=calculation_function(catalog)
                % calulate values at a single point
                % calculate distance from center point and sort wrt distance

                idx_back =  catalog.Date >= obj.periodA_start &  catalog.Date < obj.periodA_end ;
                [cumu1, ~] = histcounts(catalog.Date(idx_back),interval1_edges);

                idx_after =  catalog.Date >= obj.periodB_start &  catalog.Date <= obj.periodB_end ;
                [cumu2, ~] = histcounts(catalog.Date(idx_after),interval2_edges);

                mean1 = mean(cumu1);        % mean seismicity rate in first interval
                mean2 = mean(cumu2);        % mean seismicity rate in second interval
                sum1 = sum(cumu1);          % number of earthquakes in the first interval
                sum2 = sum(cumu2);          % number of earthquakes in the second interval
                var1 = cov(cumu1);          % variance of cumu1
                var2 = cov(cumu2);          % variance of cumu2
                % remark (db): cov and var calculate the same value when applied to a vector
                ncu1 = length(interval1_bins);         % number of bins in first interval
                ncu2 = length(interval2_bins);         % number of bins in second interval

                % compute the z value "as":
                as = (mean1 - mean2)/ sqrt(var1/ncu1 +var2/ncu2);

                % calculate the percentage
                per = -((mean1-mean2)./mean1)*100;

                % beta nach reasenberg & simpson 1992, time of second interval normalised by time of first interval
                bet = (sum2-sum1*ncu2/ncu1)/sqrt(sum1*(ncu2/ncu1));

                out = [as  per bet sum1 sum2];

            end
        end
        function ModifyGlobals(obj)
            obj.ZG.bvg = obj.Result.values;
        end
    end %methods

    methods(Static)
        function h = AddMenuItem(parent, zapFcn, varargin)
            % create a menu item
            label = 'Compare two periods (z, beta, probabilty)';
            h =  uimenu(parent, 'Label', label,...
                'MenuSelectedFcn', @(~,~)XYfun.comp2periodz(zapFcn()),...
                varargin{:});
        end
    end % static methods
end % classdef

author: CelsoReyes


text:
function [ a_lu, pivot, info ] = r8ge_fa ( n, a )

%*****************************************************************************80
%
%% R8GE_FA performs a LINPACK style PLU factorization of an R8GE matrix.
%
%  Discussion:
%
%    The R8GE storage format is used for a general M by N matrix.  A storage
%    space is made for each logical entry.  The two dimensional logical
%    array is mapped to a vector, in which storage is by columns.
%
%    This is a simplified version of the LINPACK routine DGEFA.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Dongarra, Bunch, Moler, Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, real A(N,N), the matrix to be factored.
%
%    Output, real A_LU(N,N), an upper triangular matrix and
%    the multipliers used to obtain it.  The factorization
%    can be written A = L * U, where L is a product of
%    permutation and unit lower triangular matrices and U
%    is upper triangular.
%
%    Output, integer PIVOT(N), a vector of pivot indices.
%
%    Output, integer INFO, singularity flag.
%    0, no singularity detected.
%    nonzero, the factorization failed on the INFO-th step.
%
  a_lu(1:n,1:n) = a(1:n,1:n);

  info = 0;

  for k = 1 : n-1
%
%  Find L, the index of the pivot row.
%
    l = k;
    for i = k+1 : n
      if ( abs ( a_lu(l,k) ) < abs ( a_lu(i,k) ) )
        l = i;
      end
    end

    pivot(k) = l;
%
%  If the pivot index is zero, the algorithm has failed.
%
    if ( a_lu(l,k) == 0.0 )
      info = k;
      fprintf ( 1, '\n' );
      fprintf ( 1, 'R8GE_FA - Fatal error!\n' );
      fprintf ( 1, '  Zero pivot on step %d\n', info );
      return;
    end
%
%  Interchange rows L and K if necessary.
%
    if ( l ~= k )
      t         = a_lu(l,k);
      a_lu(l,k) = a_lu(k,k);
      a_lu(k,k) = t;
    end
%
%  Normalize the values that lie below the pivot entry A(K,K).
%
    a_lu(k+1:n,k) = -a_lu(k+1:n,k) / a_lu(k,k);
%
%  Row elimination with column indexing.
%
    for j = k+1 : n

      if ( l ~= k )
        t         = a_lu(l,j);
        a_lu(l,j) = a_lu(k,j);
        a_lu(k,j) = t;
      end

      a_lu(k+1:n,j) = a_lu(k+1:n,j) + a_lu(k+1:n,k) * a_lu(k,j);

    end
  end

  pivot(n) = n;

  if ( a_lu(n,n) == 0.0 )
    info = n;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8GE_FA - Fatal error!\n' );
    fprintf ( 1, '  Zero pivot on step %d\n', info );
    return;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function Model = hlp_sift_emptymodel(varargin)
% create a new (empty) Model datastructure with default fields initialized
% varargin is an optional set of <name,value> pairs indicating fields and
% associated values to store into the set
%
% MODEL
%   .morder             Autoregressive model order
%   .AR                 Autoregressive coefficient matrices [A1 A2 ... Ap]
%                       for model order p. This is of dimension [num_vars x num_vars*p]
%   .PE                 Prediction errors (noise covariance matrices) [E0 E1 .. Ep]
%                       This is of dimension [num_vars x num_vars*(p +1)]
%   .RC                 Reflection (parcorr) coefficients [num_vars x num_vars*(p +1)]
%   .mu                 Process mean [num_vars x 1]
%   .th                 Statistics structure (if using ARFIT)
%   .winStartTimes      Timestamps of window start points (sec)
%   .winstep            Sliding window step size (sec)
%   .winlen             Sliding window length (sec)
%   .algorithm          Algorithm used (see hlp_getMVARalgorithms)
%   .modelclass         Class of algorithm ('mvar')
%   .timeelapsed        CPU time required for each time window (sec)
%   .normalize          Data normalizion method used
%   .modelapproach      Approach used ('Segmentation VAR', 'SSM',...)
%   .taperFcn           Taper function used
%
% See Also: hlp_sift_emptyconn(), hlp_sift_emptyset()

if nargin > 1 && mod(length(varargin),2)
    error('SIFT:hlp_sift_emptymodel', ...
        ['Additional arguments must be a list of <name,value> pairs. ' ...
         'At least one name is missing its corresponding value.']);
end

Model.AR = {};
Model.PE = {};
Model.RC = {};
Model.mu = {};
Model.th = {};
Model.winStartTimes = []; %EEG.CAT.times(g.winStartIdx)/1000;
Model.morder        = [];
Model.winstep       = [];
Model.winlen        = [];
Model.algorithm     = '';
Model.modelclass    = '';
Model.timeelapsed   = [];
Model.normalize     = [];
Model.modelapproach = '';
Model.taperFcn      = '';

% append additional arguments
if nargin > 1
    Model = hlp_varargin2struct(varargin,Model);
end

repo: REST


text:
classdef EventSelectionParameters
    %holds basic values for specifying samples, either out to a distance, or a specified number
    %
    %CONSTRUCTION OPTIONS
    %  obj = EventSelectionParameters()
    %  obj = EventSelectionParameters('NumClosestEvents', N)
    %  obj = EventSelectionParameters('AllEventsInRadius', R)
    %  obj = EventSelectionParameters('NumClosestEventsUpToRadius', N, R)
    %  obj = EventSelectionParameters( ... , 'DistanceUnits', unitname)
    %
    %USE CASE
    %
    %  select based on distances, where D is a vector of distances of length N, in unit U.
    %  The unit is flexible, and can be any value accepted by either checkangleunits or ValidateLengthUnits
    %
    %    MASK = obj.SelectionFromDistances(D, U)
    %
    %  MASK will be a logical vector of length N. It is TRUE where D matches the conditions contained
    %  in obj.
    %
    %see also EventSelectionParameters/EventSelectionParameters

    properties
        % Samples will be collected out to this radius. To ignore radius, set as inf (default)
        MaxSampleRadius         (1,1) double {mustBePositive}    = inf
        % These many samples will be collected. To collect all, set as inf (default)
        NumClosestEvents        (1,1) double {mustBePositive}    = inf
        % Units of distance for the sample radius. Defaults to kilometers.
        DistanceUnits           (1,1) string {mustBeNonempty}    = standardizeDistanceUnits('km')
    end

    properties(Dependent)
        UseNumClosestEvents
        UseEventsInRadius
    end

    properties(Hidden)
        PrevMaxSampleRadius     (1,1) double    = inf;
        PrevNumClosestEvents    (1,1) double    = inf;
    end

    properties(Dependent,Hidden)
        % for backwards compatibility
        RadiusKm
    end

    methods
        function obj = EventSelectionParameters(style, varargin)
            %EventSelectionParameters constructor
            %obj = EventSelectionParameters()
            %obj = EventSelectionParameters('NumClosestEvents', N) choose the N closest events
            %
            %obj = EventSelectionParameters('AllEventsInRadius', R) choose ALL events within radius R
            %
            %obj = EventSelectionParameters('NumClosestEventsUpToRadius', N, R) choose the N
            %closest events, limited to radius R.
            %
            % The radius units can be specified as an optional Name,Value pair.
            %
            %obj = EventSelectionParameters( ... , 'DistanceUnits', unitname) where unitname is
            %unit understood by the validateLengthUnits or checkangleunits functions.

            if ~exist('style','var')
                style = 'none';
            end
            p=inputParser();
            switch style
                case 'AllEventsInRadius'
                    p.addRequired('MaxSampleRadius');
                    p.addParameter('DistanceUnits',obj.DistanceUnits);
                case 'NumClosestEvents'
                    p.addRequired('NumClosestEvents');
                case 'NumClosestEventsUpToRadius'
                    p.addRequired('NumClosestEvents');
                    p.addRequired('MaxSampleRadius');
                    p.addParameter('DistanceUnits',obj.DistanceUnits,@isStandardDistanceUnit);
                case 'none'
            end
            p.parse(varargin{:});
            obj=copyfields(obj, p.Results);
            %{
            fn =fieldnames(p.Results);
            for i=1:numel(fn)
                obj.(fn{i}) = p.Results.(fn{i});
            end
            %}

        end

        function rk = get.RadiusKm(obj)
            if isinf(obj.MaxSampleRadius)
                rk = obj.MaxSampleRadius;
            elseif obj.DistanceUnits == "kilometer"
                rk = obj.MaxSampleRadius;
            elseif obj.DistanceUnits== "degrees"
                rk = deg2km(obj.MaxSampleRadius);
            elseif obj.DistanceUnits == "radians"
                rk = rad2km(obj.MaxSampleRadius);
            else
                rk = unitsratio('kilometer', obj.DistanceUnits) * obj.MaxSampleRadius;
            end
        end
        function obj = set.RadiusKm(obj,val)
            obj.MaxSampleRadius = val;
            obj.DistanceUnits = 'kilometer';
        end

        function tf = get.UseNumClosestEvents(obj)
            tf = obj.NumClosestEvents < inf;
        end

        function tf = get.UseEventsInRadius(obj)
            tf = obj.MaxSampleRadius < inf;
        end

        function obj = set.UseNumClosestEvents(obj, val)
            if val
                obj.NumClosestEvents = obj.PrevNumClosestEvents;
            else
                if ~isinf(obj.NumClosestEvents)
                    obj.PrevNumClosestEvents = obj.NumClosestEvents;
                end
                obj.NumClosestEvents = inf;
            end
        end

        function obj = set.UseEventsInRadius(obj, val)
            if val
                obj.MaxSampleRadius = obj.PrevMaxSampleRadius;
            else
                if ~isinf(obj.MaxSampleRadius)
                    obj.PrevMaxSampleRadius = obj.MaxSampleRadius;
                end
                obj.MaxSampleRadius = inf;
            end
        end

        function obj = set.DistanceUnits(obj, units)
            obj.DistanceUnits = standardizeDistanceUnits(units);
        end

        function obj = set.MaxSampleRadius(obj, value)
            prev = obj.MaxSampleRadius;
            obj.MaxSampleRadius = value;
            if ~isinf(prev)
                obj.PrevMaxSampleRadius = prev; %#ok<MCSUP>
            end
        end
        function obj = set.NumClosestEvents(obj,value)
            if value==floor(value) || isinf(value)
                prev = obj.NumClosestEvents;
                obj.NumClosestEvents = value;
                if ~isinf(prev)
                    obj.PrevNumClosestEvents = prev; %#ok<MCSUP>
                end
            else
                error('NumClosestEvents must be either a positive integer value, or be inf');
            end
        end

        function mask = SelectionFromDistances(obj, distances, distUnits)
            %select values, based on distances
            %
            %MASK = obj.SelectionFromDistances(DISTANCES, UNITS) returns a logical vector of
            %length N that is TRUE for each value of DISTANCES that matches the conditions
            %contained in this object. The UNITS specify the units of distances, (eg. 'km',
            %'radians', 'degrees', etc..) and accepts any value understood by either the
            %checkangleunits or validateLengthUnit functions.
            %
            %This object's maximum radius will be compared (with the appropriate conversion)
            %with the provided distance vector.
            %
            % Example of automatic unit conversion
            %  >> esp = EventSelectionParameters('AllEventsInRadius', 1, 'DistanceUnits', 'degrees');
            %  >> esp.SelectionFromDistances([100 110 120], 'km')  % one degree of arc is ~111km
            %     ans =
            %       1×3 logical array
            %        1   1   0
            %
            %  >> esp.DistanceUnits = 'miles';
            %  >> esp.SelectionFromDistances([2.0, 1.5 1.0],'km'); % comparing against 1 mile
            %     ans =
            %       1×3 logical array
            %        0   1   1
            %
            %
            %When specifying a number of events and several distances fall across the boundary, it
            %is undefined which values will be selected. For example...
            %
            %  >> D = [ 1 1 1 1 1 1 1 1 ];  % 8 events, all at the same distance
            %  >> esp = EventSelectionParameters('NumClosestEvents',5);
            %  >> esp.SelectionFromDistances(D, 'km')
            %
            %   ans =
            %     1×8 logical array
            %     1   1   1   1   1   0   0   0
            %
            %see also checkangleunits, validateLengthUnit

            %% deal with the maximum radius

            if isinf(obj.MaxSampleRadius)                                 % distances do not matter
                mask = true; % scaler because unused later
            elseif distUnits == obj.DistanceUnits                       % simple case: Units match
                mask = distances <= obj.MaxSampleRadius;
            else
                distUnits = standardizeDistanceUnits(distUnits);
                if distUnits == obj.DistUnits
                    mask = distances <= obj.MaxSampleRadius;
                else
                    mask = [];
                    % convert Selection Units into provided units, then compare
                    OtherIsAngle = ismember(distUnits,["radians" , "degrees"]);
                    ThisIsAngle  = ismember(obj.DistanceUnits,["radians" , "degrees"]);
                end
            end

            % remember that UNITSRATIO arguments are (TO , FROM)
            if ~isempty(mask)
                % mask has been calculated already

            elseif OtherIsAngle == ThisIsAngle                      % both are length or both are angle
                maxDist = unitsratio(distUnits, obj.DistanceUnits) * obj.MaxSampleRadius;
                mask = distances < maxDist;

            elseif distUnits == "radians"
                maxDist = unitsratio(obj.DistanceUnits,'kilometer') * km2rad(obj.MaxSampleRadius);
                mask = distances < maxDist;

            elseif distUnits == "degrees"
                maxDist = unitsratio(obj.DistanceUnits,'kilometer') * km2deg(obj.MaxSampleRadius);
                mask = distances < maxDist;

            elseif obj.DistanceUnits == "radians"
                maxDist = unitsratio(distUnits, 'kilometer') * rad2km(obj.MaxSampleRadius);
                mask = distances < maxDist;

            elseif obj.DistanceUnits == "degrees"
                maxDist = unitsratio(distUnits, 'kilometer') * deg2km(obj.MaxSampleRadius);
                mask = distances < maxDist;
            end

            %% deal with the N closest events
            if ~isinf(obj.NumClosestEvents)
                [~,I]=mink(distances, obj.NumClosestEvents);
                % I'm sure this logic could be simplified...
                II = false(size(distances));
                II(I)=true;
                mask = mask & II;
            end

        end
    end

    methods(Static)
        function obj = fromStruct(s)
            %Backwards Compatibility with EventSelectionChoice
            if ~isstruct(s)
                error('expected a struct');
            end
            obj = EventSelectionParameters();

            if isfield(s,'UseEventsInRadius') && s.UseEventsInRadius
                if isfield(s,'RadiusKm')
                    obj.MaxSampleRadius = s.RadiusKm;
                end
            end
            if isfield(s,'UseNumNearbyEvents') && s.UseNumNearbyEvents
                if isfield(s,'NumClosestEvents')
                    obj.NumClosestEvents = s.NumClosestEvents;
                end
                if isfield(s,'MaxRadiusKm')
                    obj.MaxSampleRadius = s.MaxRadiusKm;
                end
            end
        end
    end
end


text:
function [delayDopPlot,Doppler,delay]=delayDopplerPlotNBPulseDop(x,y,M1,M2,wDoppler,T0)
%%DELAYDOPPLERPLOTNBPULSEDOP This function produces the complex
%               delay-Doppler plot for a narrowband signal assuming that
%               the same waveform is repeated every single pulse or that
%               each pulse ends with a region of non-broadcasting (>100%
%               duty cycle) and the maximum delay is less than that
%               duration. Range migration is not taken into account. Such
%               plots can be used for detection.
%
%INPUTS: x The NPulseX1 complex baseband waveform that is repeated every
%          pulse repetition interval (PRI). If NPulse<Ns, where Ns is the
%          number of samples in a PRI, then  it is assumed that the rest of
%          the signal is filled with zeros.
%        y The NsXNB set of complex baseband samples received over each of
%          the NB pulse repetition intervals. Ns is the number of samples
%          per interval. It is assumed that no gaps exist between
%          intervals.
%    M1,M2 These are the scaling factors to increase the resolution in
%          delay (M1) and Doppler (M2). These are integers. The number of
%          bins used in the plot increases by a factor of these. If omitted
%          or empty matrices are passed, a value of 1 is used.
% wDoppler This is an optional NBX1 or 1XNB window vector to lower the
%          Doppler sidelobes. One could use a window from, for example, the
%          windowFunSym function. If this parameter is omitted or an empty
%          matrix is passed, then the default of no windowing (the
%          rectangular window) is used. Range sidelobes should be lowered
%          either by windowing the pulse before broadcast or by other means
%          such as mismatched filtering.
%       T0 The sampling period. This parameter is only required if the
%          outputs Doppler and delay are desired.
%
%OUTPUTS: delayDopPlot This is an (M1*Ns)X(M2*NB) complex range-Doppler
%                      map. The values have been scaled so that peaks are
%                      approximate representations of the correct complex
%                      amplitude (not counting the effects of any
%                      windowing). However, as demonstrated in the example
%                      below, range migration will cause the complex
%                      amplitude of the peak to be biased from the true
%                      value.
%              Doppler An (M2*NB)X1 vector holding the values of the
%                      Doppler shifts used for each row on delayDopPlot.
%                      This requires that T0 be passed. These are positive
%                      and negative. These values are inverse seconds and
%                      when multiplied by the propagation speed (e.g. c=
%                      Constants.speedOfLight) divided by the carrier
%                      frequency fc will provide the range rate. If seconds
%                      are used for time and meters for distance, when
%                      multiplied by c/fc one gets range rate in meters per
%                      second. Positive is going away from the sensor.
%                delay An (M1*Ns)X1 vector holding the values of the delays
%                      used for each columns on delayDopPlot. This requires
%                      that T0 be passed. These are all non-negative.
%                      Multiplied by the propagation speed c this gives the
%                      range.
%
%Range-Doppler plots are discussed in many sources, but many aspects of the
%approximations are often omitted, so a brief derivation is given here.
%
%Let x(t) be the complex baseband waveform sent, concatenating all pulses
%together. The received passband signal is modeled as
%x_p(t)=real(exp(1j*2*pi*fc)*x(t))
%where x(t) is the complex baseband signal. To come to baseband, one
%would take the real part of exp(1j*2*pi*fc) times the complex baseband
%signal. In practice, one would use quadrature demodulation to take the
%real measured signal and obtain the in-phase (real) and quadrature
%(imaginary) parts of the complex baseband signal.
%
%The delayed and Doppler-shifted passband signal is
%y_p(t)=real(A*exp(1j*2*pi*fc*(t-tau-a*t)*x(t-tau-a*t))
%where tau is the delay measured from the time the pulse begins
%broadcasting, a is a Doppler shift, and A is a complex amplitude. This
%assumes a target with a constant range rate (constant Doppler frequency).
%For targets accelerating in range, additional terms are needed. The
%complex baseband received signal is thus
%y(t)=A*exp(-1j*2*pi*fc*(tau+a*t))*x(t-tau-a*t)
%
%The range-Doppler map is just a matched filter. A matched filter
%integrates the product of the complex conjugate of the ideal signal times
%the received signal. Here, the matched filter while observing over a time
%T is:
%integral_0^T exp(1j*2*pi*fc*(tau+a*t))*conj(x(t-tau-a*t))*y(t) dt
%In the problem at hand, the signal is approximated as a set of NB blocks
%or pulses of duration TB. Thus, T=NB*T. To simplify the signal
%processing, we assume either that all blocks have the same waveform or
%that each block ends with a period of no broadcasting (<100% duty cycle)
%and the maximum delay does not exceed that period. Additionally, the
%approximation is made that the delays due to the Doppler shift are
%constant over each block. All together, this makes the use of circular
%convolutions (and hence FFTs in the signal processing) possible.
%The matched filter under such a model is
%integral_0^T sum_{i=0}^{N_B-1} exp(1j*2*pi*fc*(tau+a*i*TB))*conj(x_i(t-tau-a*i*TB))*y_i(t) dt
%where x_i and y_i are the signals from the ith block. It is assumed that
%outside of a time windows of length T_B each of those signals is zero.
%To further simplify the signal processing, we assume that the Doppler
%shift delay in the x term is not significant (related to a narrowband
%approximation). Thus, the signal processing model for the matched filter
%is
%exp(1j*2*pi*fc*tau)*sum_{i=0}^{N_B-1}exp(1j*2*pi*fc*a*i*TB)*integral_0^T conj(x_i(t-tau))*y_i(t) dt
%As the signal y_i is sampled, we will approximate the integral in time
%using a Riemann sum. We will also say that tauHat=tau*T0 (T0 is the sample
%period) where tauHat is an integer >=0 for the discrete delay. Also, we
%will throw out the phase constant of exp(1j*2*pi*fc*tau) in front. The
%result is thus:
%sum_{i=0}^{N_B-1}exp(1j*2*pi*fc*a*i*TB)*sum_{k=0}^{Ns-1} conj(x(k-tauHat))*y(i,k)
%Thus, the integral has turned into a discrete sum over the samples in each
%block (the arguments of x and y access the discrete elements). There is no
%block marking on x, because it is assumed to be the same for all blocks.
%Looking at the form of the sum, it can be shown (for example, see Chapter
%5.7 of [1] for useful identities) that it is equivalent to the circular
%convolution of flipud(conj(x)) with y. This can be evaluated for all
%discrete delays tau and block i in parallel using ffts and an ifft. This
%is what is implemented below via the circConv function. The outer sum in
%the above can be seen to be NB times the ifft taken across all values of i
%when we replace a with a=aHat/(fc*TB*NB) where aHat is a discrete Doppler
%shift from 0 to (NB-1). Thus, the ifft function, as used below, evaluated
%all of the Doppler shifts in parallel for the matched filter, scaled by
%1/NB.
%
%Thus, as described above, we are evaluating
%(1/NB)*sum_{i=0}^{N_B-1}exp(1j*2*pi*aHat*i/NB)*sum_{k=0}^{Ns-1} conj(x(k-tauHat,i))*y(i,k)
%for discrete values of aHat and tauHat. However, it would be nice if the
%output could be related to the amplitude A. Well, if we assume that tauHat
%and aHat perfectly match the original signal AND the original signal was
%generated using the various narrowband simplications mentioned, then we
%have to multiply the result by 1/(x'*x). This is what is done below.
%However, opne must realize that the various approximations means that the
%amplitude estimates can still be biased. This is demonstrated in the
%example below.
%
%Often, the sidelobes of the signal in Doppler might be high. Thus, the
%input wDoppler can be used to apply tapering across each bin before the
%final ifft of the algorithm.
%
%The final thing to note is that the above method describes only processing
%for positive Doppler values. Due to the aliasing of large values to
%negative values with circular convolutions and FFTs, we use fftshift to
%get positive and negative values in the unambiguous Doppler region.
%
%Note that TB*c gives the maximum unambiguous range, where c is the speed
%of propagation in the medium, for example c=Constants.speedOfLight.
%Note that c/(2*fc*TB) is the magnitude of the maximum unambiguous range
%rate. One can get those as
% range=c*delay;
% rangeRate=Doppler*(c/fc);
%
%EXAMPLE:
%Our signal is a up-chirp. Here, we create the time-delayed Doppler shifted
%signal at baseband. The target signal is generated in two ways: Once
%including all of the phase shifts associated with range migration, and
%once with the ideal signal processing model that is used to derive the
%range-Doppler map using Fourier transforms to make it fast. We choose the
%target location to be EXACTLY on a delay-Doppler bin. It will be seen
%that the complex amplitude obtained using the ideal model matches the
%complex amplitude used in generating the signal (no noise is added).
%However, the complex amplitude obtained when using the more realistic
%model for the signal, including all range migration, is biased. This
%highlights how the lack of migration compensation biases complex amplitude
%determination. The range-Doppler plot of the realistic signal is
%displayed.
% fc=1e9;%1GHz carrier frequency.
% B=2e6;%2Mhz bandwidth.
% %Baseband start and end frequencies.
% fStart=-B/2;
% fEnd=B/2;
% %Sampling rate is two times the Nyquist rate.
% T0=1/(2*2*fEnd);%Sampling period in seconds.
% T=2e-5;%Chirp duration in seconds.
%
% PRF=2000;%Pulse repetition frequency (Hertz)
% TB=1/PRF;%The pulse repetition period.
% %The number of samples per PRI. The above parameters were chosen so that
% %this is an integer. Fix just deals with finite precision errors.
% Ns=fix(TB/T0);
%
% %Generate the reference signal. This is an up-chirp.
% x=LFMChirp(T,fStart,fEnd,T0);
% x=x(:);
%
% %We will use 64 pulse repetition intervals.
% NB=64;
%
% %True target parameters
% c=Constants.speedOfLight;
% rTrue=512*c*T0;%Target placed on
% tau=rTrue/c;%The true delay (s).
%
% %The true range rate (m/s).
% rrTrue=6/NB*(c/(fc*TB));
% a=rrTrue/c;
%
% %Allocate space for the received signal. The first dimensions is "fast
% %time"; the second dimension if "slow time".
% yReal=zeros(Ns,NB);
% yIdeal=zeros(Ns,NB);
%
% %Create the received signal, properly delayed and Doppler shifted for
% %each PRI. The same waveform is used in each PRI.
% t=0:T0:((Ns-1)*T0);%Sample times
% for i=0:(NB-1)
%     %The subtraction of i*TB deals with the start time of this pulse if
%     %it is not time zero.
%     tCur=t-a*t-tau-i*TB;
%     yReal(:,i+1)=512*exp(-1j*2*pi*fc*(tau+a*t)).*LFMChirp(T,fStart,fEnd,tCur);
%
%     tCur=t-tau-i*TB;
%     yIdeal(:,i+1)=512*exp(-1j*2*pi*fc*(tau+a*i*TB)).*LFMChirp(T,fStart,fEnd,tCur);
%
%     t=t+TB;%Increment to the next time step.
% end
%
% M1=1;
% M2=1;
% %If desired, one could use a window in the Doppler domain. For example,
% %by using wDoppler=windowFunSym(NB,'Blackman',1); Here, we choose to use
% %no window (same as a rectangular window).
% wDoppler=[];
% valsReal=delayDopplerPlotNBPulseDop(x,yReal,M1,M2,wDoppler,T0);
% [valsIdeal,Doppler,delay]=delayDopplerPlotNBPulseDop(x,yIdeal,M1,M2,wDoppler,T0);
%
% %Note that the amplitude of the ideal signal matches the true amplitude
% %used, but that of the real signal (with range migration) is biased.
% [magReal,idxReal]=max(abs(valsReal(:)));
% [magIdeal,idxIdeal]=max(abs(valsIdeal(:)));
% amplitudeMagReal=magReal
% amplitudeMagIdeal=magIdeal
%
% range=c*delay;
% rangeRate=Doppler*(c/fc);
%
% %We will make sure that the maximum point on the plots matches the true
% %inputs, since the true inputs were made to align with range-Doppler cells.
% [idxR,idxRR]=ind2sub(size(valsReal),idxReal);
% rDiffReal=range(idxR)-rTrue
% rrDiffReal=rangeRate(idxRR)-rrTrue
%
% [idxR,idxRR]=ind2sub(size(valsIdeal),idxIdeal);
% rDiffIdeal=range(idxR)-rTrue
% rrDiffIdeal=rangeRate(idxRR)-rrTrue
%
% figure(1)
% clf
% imagesc([rangeRate(1),rangeRate(end)],[range(1), range(end)]/1e3,10*log10(abs(valsReal)));
% set(gca,'YDir','normal')
% caxis([-30 27])
% colormap(jet(256))
% h1=xlabel('Range Rate (m/s)');
% h2=ylabel('Range (km)');
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%
%REFERENCES:
%[1] S. K. Mitra, Digital Signal Processing: A Computer-Based Approach,
%    3rd ed. Boston: McGraw Hill, 2006.
%
%November 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(M1))
    M1=1;
end

if(nargin<4||isempty(M2))
   M2=1;
end

NB=size(y,2);
Ns=size(y,1);

%The extra circshift aligns the zero-Doppler point to the position assumed
%in the delay output.
rangePlots=circshift(circConv(flipud([conj(x);zeros(Ns-length(x),1)]),y,Ns),1);

if(nargin>4&&~isempty(wDoppler))
    %Perform windowing to lower Doppler sidelobes.
    rangePlots=bsxfun(@times,wDoppler(:).',rangePlots);
end

delayDopPlot=ifftshift(ifft(rangePlots,NB*M2,2),2);

%Remove the scaling at the matched points so that we can get the true
%complex amplitude back.
delayDopPlot=M1*M2*delayDopPlot/(x'*x);

if(nargout>1)
    %T0 must be provided for these outputs.

    delay=(0:(T0/M1):((T0/M1)*(M1*Ns-1))).';
    %The pulse repetition interval.
    TB=Ns*T0;

    numNeg=(NB*M2-1)-ceil((NB*M2-1)/2)+1;
    Doppler=([(-numNeg):1:-1,0:(floor((NB*M2-1)/2)-1)].'/(NB*M2))*(1/TB);
end

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
%
% FILE NAME:    curateRxnAssoc.m
%
% PURPOSE: Curate rxnAssoc array structure by adding the reaction bounds
%          information extracted from HMR2 and Recon3D
%

% Load input
load('rxnAssoc.mat');       % rxnAssoc
load('humanGEM.mat');       % humanGEM v0.2.0
load('Recon3DRaven.mat');   % Recon3D

% Add fields for reaction bounds information
[~, indHMR]=ismember(rxnAssoc.rxnHMRID,ihuman.rxns);
rxnAssoc.lbHMR=ihuman.lb(indHMR);
rxnAssoc.ubHMR=ihuman.ub(indHMR);

[~, indRecon3D]=ismember(rxnAssoc.rxnRecon3DID,Recon3DRaven.rxns);
rxnAssoc.lbRecon3D=Recon3DRaven.lb(indRecon3D);
rxnAssoc.ubRecon3D=Recon3DRaven.ub(indRecon3D);

% Save to modelIntegration subfolder
save('rxnAssoc.mat','rxnAssoc');

author: SysBioChalmers

repo: Human-GEM


text:
function [phw,stat]=model_phw(time,sat,type,opt,rs,rr,vs,phw0)

stat=1; ds=zeros(1,3); dr=zeros(1,3);
if opt==0,phw=phw0;stat=0;return;end

%satellite yaw attitude model
[exs,eys]=sat_yaw(time,sat,type,opt,rs,vs);

%unit vector satellite to receiver
r=rr-rs;
ek=r/norm(r);

%unit vectors of receiver antenna
[~,E]=xyz2blh(rr);
exr(1)= E(2); exr(2)= E(5); exr(3)= E(8); %x = north
eyr(1)=-E(1); eyr(2)=-E(4); eyr(3)=-E(7); %y = west

%phase windup effect
eks=cross(ek,eys);
ekr=cross(ek,eyr);
for i=1:3
    ds(i)=exs(i)-ek(i)*dot(ek,exs')-eks(i);
    dr(i)=exr(i)-ek(i)*dot(ek,exr')+ekr(i);
end
cosp=dot(ds,dr)/norm(ds)/norm(dr);
if cosp<-1
    cosp=-1;
elseif cosp>1
    cosp=1;
end

ph=acos(cosp)/2/pi;
drs=cross(ds,dr);
tmp=dot(ek,drs);
if tmp<0
    ph=-ph;
end
phw=ph+floor(phw0-ph+0.5);

return

repo: GINav


text:
function title = p10_title ( )

%*****************************************************************************80
%
%% P10_TITLE returns a title for problem 10.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 February 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, a title for the problem.
%
  title = 'The cosine combo';

  return
end


text:
function Z = minus(X,Y)
%MINUS Binary subtraction (-) for tensors.
%
%   MINUS(A,B) is called for the syntax 'A - B' when A or B is a tensor. A
%   and B must have the same size, unless one is a scalar. A scalar can be
%   subtracted from a tensor of any size.
%
%   See also TENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

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

Z = tenfun(@minus,X,Y);

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
% Simulating 2*2 MIMO-LDPC Base-Band Systems
% Copyright (C2010-2015) Yang XIAO, Beijing Jiaotong University, Aug.10, 2010, E-Mail: yxiao@bjtu.edu.cn.
% This program can simulate 2*2 MIMO-LDPC base-band systems.
% The MIMO system's design can refer my following book.
% [1] Yang Xiao, MIMO Multiple Antenna Wireless Communication Systems, Press of Posts and Telecommunications, Beijing, 2009.
% Different from the MIMO scheme of IEEE 802.16e, our MIMO-LDPC system has no space-time coding,
% while it achieved good BER performance.
%---------------------------
clear;
p=87;   %p is the prime number of rank of sub-matrices of parity check matrix;
k2=3;   % k is the row weight of parity check matrix;
j2=6;   % j is the column weight of parity check matrix;
M=k2*p; % M is the number of row of parity check matrix;
N=j2*p;  % N is the number of column of parity check matrix;

NT1=1;
E0=eye(p);
EZ=zeros(p);
a=3; b=5;
R=0.5; % coding rate
frame_num =200;
Npf=2*N*frame_num
dlta=1/Npf;

EZ=zeros(p);
E0=eye(p);
% The design of parity check matrix of QC LDPC without girth_4, see my book.
E11=E0;                   E12=circshift(E0,a*NT1);   E13= circshift(E0,a^2*NT1);   E14= circshift(E0,a^3*NT1);   E15= circshift(E0,a^4*NT1);    E16= circshift(E0,a^5*NT1);
E21=EZ;                   E22=E0;                    E23=circshift(E0,a^2*b*NT1);  E24=circshift(E0,a^3*b*NT1);  E25= circshift(E0,a^4*b*NT1);  E26= circshift(E0,a^5*b*NT1);
E31=circshift(E0,b^2*NT1);E32=EZ;                    E33=E0;                       E34=circshift(E0,b^2*a^3*NT1);E35= circshift(E0,b^2*a^4*NT1);E36= circshift(E0,b^2*a^5*NT1);

h1=[E11 E12 E13 E14 E15 E16;
    E21 E22 E23 E24 E25 E26;
    E31 E32 E33 E34 E35 E36];
 A(:,1:M)=h1(:,1:M);
 B(:,1:M)=h1(:,M+1:N);
%

d=mod(inv_GF2(A)*B,2);
E00=eye(M);
% The generator matrix of QC LDPC, see my book.
G=[d' E00];
H=sparse(h1);

N1=0;
N2=16;
% Flat fading MIMO Channel
h11=0.9;
h12=0.3;
h21=0.4;
h22=0.73;

SNRindB1=N1:1:N2;
for i=1:length(SNRindB1)
error_count1= 0;
error_count2= 0;
    for f_n = 1:frame_num
    SNR=(10^(SNRindB1(i)/10));
  	sigma = 1/sqrt(2*R*SNR);
    x1 = (sign(randn(1,size(G,1)))+1)/2; % random information bits
    x2 = (sign(randn(1,size(G,1)))+1)/2; % random information bits
 y1 = mod(x1*G,2);                   % LDPC Encoding for signal_1 from atenna 1 of the transmitter
 y2 = mod(x2*G,2);                   % LDPC Encoding for signal_2 from atenna 2 of the transmitter
    z1=2*y1-1;                       % BPSK modulation
    z2=2*y2-1;                         % BPSK modulation
    z11=h11*z1+h12*z2+sigma*randn(1,size(G,2));    % AWGN transmission
    z22=h21*z1+h22*z2+sigma*randn(1,size(G,2));
HC=[h11 h12;h21 h22];   % Channel Matrix
HC1=HC^(-1);
u=HC1*[z11;z22];        % ZF spatial decoding
u1=u(1,:);              % recieved signal_1 from atenna 1 of the receiver
u2=u(2,:);              % recieved signal_2 from atenna 1 of the receiver
    f11=1./(1+exp(-2*u1/sigma^2));        % likelihoods
    f01=1-f11;
    [z1_hat, success, k] = ldpc_decode(f01,f11,H);   %  LDPC decoding for signal_1
    x1_hat = z1_hat(size(G,2)+1-size(G,1):size(G,2));
    x1_hat1 = x1_hat';
    f12=1./(1+exp(-2*u2/sigma^2));        % likelihoods
    f02=1-f12;
    [z2_hat, success, k] = ldpc_decode(f02,f12,H);   %  LDPC decoding for signal_2
    x2_hat = z2_hat(size(G,2)+1-size(G,1):size(G,2));
    x2_hat1 = x2_hat';
    error_count2= sum(xor(x1,x1_hat1))+sum(xor(x2,x2_hat1));    % bit error count
    error_count1 = error_count1 + error_count2;
    end
    BER(i)= error_count1/Npf
    if BER(i)<.1*dlta
    BER(i)=.08*dlta;
    end
end
figure(1)
semilogy(SNRindB1,BER);
xlabel('SNR(dB)');
ylabel('BER');
title('2*2 MIMO-LDPC Simulation');
axis([N1 N2 dlta 1]);
%diary off

author: Sable

repo: mcbench-benchmarks


text:
function [C] = tensor_product(A,B)
  Au = tensor_unfold(A);
  Bu = tensor_unfold(B);

  Auc = circmat(Au,size(A));

  AucBu = Auc*Bu;

  C = tensor_fold(AucBu,size(A));
end

repo: lrslibrary


text:
%% housekeeping
clc
%% RISE the model
linear=~true;
if linear
    m=rise('targets_lin');
else
    m=rise('targets','steady_state_file','sstate_model');
end
%% get the parameters
version='a';
% 'a' % model with an endogenous inflation target.
% 'b' % model with an exogenous inflation target.
% 'd' % model with backward-looking price setting.
start_at_mode=false;

[p,priors]=create_parameters(version,linear,start_at_mode);

%% push the parameters
m=set(m,'parameters',p);

%% create data
[data]=create_data();

%% estimate the model
ms=estimate(m,'data',data,'estim_priors',priors);

%% Impulse responses
myirfs=irf(ms);

%% plot responses
myvars={'Y','PAI','R','X'};
locs=locate_variables(myvars,ms.endogenous.name);
myvtex=ms.endogenous.tex_name(locs);
sstate=get(ms,'sstate');
ssdev=false;

shocks={'EPS_A','EPS_E','EPS_Z','EPS_V','EPS_PAI'};%m.exogenous.name;
locs=locate_variables(shocks,ms.exogenous.name);
myshtex=ms.exogenous.tex_name(locs);
close all
figure('name','Impulse reponses');
iter=0;
ss=1;
for ishock=1:numel(shocks)
    shock=shocks{ishock};
    for ivar=1:numel(myvars)
        vname=myvars{ivar};
        if ssdev
            ss=sstate.(vname);
        end
        iter=iter+1;
        subplot(5,4,iter)
        plot('0:16',myirfs.(shock).(vname)/ss,'linewidth',2)
        title([myvtex{ivar},' to ',myshtex{ishock}])
    end
end

author: jmaih


text:
function my_move_net(net_info, destination)

if strcmp(destination, 'gpu')

    for l=1:numel(net_info.ref.layers)
      switch net_info.ref.layers{l}.type
        case 'conv'
          for f = {'filters', 'biases', 'filtersMomentum', 'biasesMomentum'}
            f = char(f) ;
            if isfield(net_info.ref.layers{l}, f)
              net_info.ref.layers{l}.(f) = move_to_gpu_values(net_info.ref.layers{l}.(f)) ;
            end
          end
        otherwise
          % nothing to do ?
      end
    end

    net_info.ref.net_on_gpu=true;

end

if strcmp(destination, 'cpu')

    for l=1:numel(net_info.ref.layers)
      switch net_info.ref.layers{l}.type
        case 'conv'
          for f = {'filters', 'biases', 'filtersMomentum', 'biasesMomentum'}
            f = char(f) ;
            if isfield(net_info.ref.layers{l}, f)
              net_info.ref.layers{l}.(f) = move_to_cpu_values(net_info.ref.layers{l}.(f)) ;
            end
          end
        otherwise
          % nothing to do ?
      end
    end

    net_info.ref.net_on_gpu=false;

end

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function mtrScoreNonConTrackPaths(subjectDir, fgFile, roi1File, roi2File, smoothParam, shapeParam, lengthParam, optionsFile, fgOutFile)

% mtrScoreNonConTrackPaths(subjectDir, fgFile, roi1File, roi2File, smoothParam, shapeParam, lengthParam, optionsFile, fgOutFile)
%
% 06/10/2007: AJS created it.
%
% This function takes STT fibers (fgFile), takes only those that pass
% through two desired ROIs (roi1File, roi2File), clips them, then writes
% the clipped fibers and their ConTrack scores to file (fgOutFile).
%
% It requires certain parameters that Tony will hopefully document here
% soon (smoothParam, lengthParam, optionsFile).

%% Setup the options file
% Prepare the roi filename by concatenating the ROI names
[pathstr, file1, ext, versn] = fileparts(roi1File); %#ok<NASGU>
[pathstr, file2, ext, versn] = fileparts(roi2File); %#ok<NASGU>
maskFile = fullfile(pathstr,['mask_' file1 '_' file2 '.nii.gz']);

mtrCreateConTrackOptionsFromROIs(subjectDir, roi1File, roi2File, optionsFile, optionsFile, 'junk', maskFile);

%% Update options file with user parameters
dt6 = load(fullfile(subjectDir,'dt6.mat'),'xformToAcPc');
mtr = mtrLoad(fullfile(subjectDir,'conTrack',optionsFile),dt6.xformToAcPc);
mtr = mtrSet(mtr, 'smooth_std', smoothParam);
mtr = mtrSet(mtr, 'abs_normal', lengthParam);
mtr = mtrSet(mtr, 'shape_params_vec', [shapeParam 0.15 100]); %#ok<NASGU>

%% Clip the pathways to the supplied ROIs
mtrClipFiberGroupToROIs( fullfile(subjectDir,'fibers',fgFile), ...
                         fullfile(subjectDir,'dt6.mat'), ...
                         fullfile(subjectDir,'ROIs',roi1File), ...
                         fullfile(subjectDir,'ROIs',roi2File), ...
                         fullfile(subjectDir,'bin','wmMask.nii.gz'), ...
                         fullfile(subjectDir,'fibers',fgOutFile) );

%% Export to format for the executable
[pathstr, file, ext, versn] = fileparts(fgOutFile); %#ok<NASGU>
fgOutFileDat = fullfile(pathstr,[file '.dat']);
mtrExportFiberGroupToMetrotrac( fullfile(subjectDir,'fibers',fgOutFileDat), ...
                                fullfile(subjectDir,'fibers',fgOutFile), ...
                                fullfile(subjectDir,'bin','fa.nii.gz') );

%% Check environment we are running in
if(ispc)
    executable = which('resample_pdb.exe');
elseif(strcmp(computer,'GLNXA64'))
    executable = which('resample_pdb.glxa64');
elseif (strcmp(computer,'GLNX86'))
    executable = which('resample_pdb.glx');
else
    error('Not compiled for %s.',computer);
end

%% Run the executable
workDir = fullfile(subjectDir,'conTrack');
args = sprintf(' -i %s -o tempScoreNonCTRK  -lpn %d %s', optionsFile, intmax, fullfile(subjectDir,'fibers',fgOutFileDat));
cmd = [executable args];
disp(cmd); disp('...')
if(ispc)
    mtrExecuteWinBash(cmd,workDir);
else
    system(cmd,'-echo');
end
disp('Done')

if( exist(fullfile(subjectDir,'fibers',fgOutFileDat),'file') )
    delete(fullfile(subjectDir,'fibers',fgOutFileDat));
end
tempPaths = dir(fullfile(subjectDir,'conTrack','tempScoreNonCTRK*.dat'));
if(size(tempPaths,1)>1)
    fileList = tempPaths(1).name;
    for ii = 2:size(tempPaths,1); fileList = [fileList ', ' tempPaths(ii).name]; end
    error(['Special temp name already used for intermediate path creation please remove the following files ' fileList]);
end
movefile(fullfile(subjectDir,'conTrack',tempPaths(1).name),fullfile(subjectDir,'fibers',fgOutFileDat),'f');
disp(['Final pathway file has been written to ' fullfile(subjectDir,'fibers',fgOutFileDat)]);

return;


text:
function S = eliminator(n, ROWS)
%
%ELIMINATOR
%
% Returns a matrix S that is called eliminator.
% S is a matrix such that S*A (assuming that this product
% exists) will not contain rows listed in row vector ROWS.
% Numbers of rows in ROWS can be unsorted.
%
% Example:
%     A = hilb(6);
%     S = eliminator(6, [4 6 2]);
%     B = S * A;
% Matrix B will contain rows 1, 3, 5 from the matrix A,
% matrix B will not contain rows 2, 4, 6 from matrix A.
%
% (C) 2000 Igor Podlubny, Blas Vinagre, Tomas Skovranek
%
% See:
% [1] I. Podlubny, A.Chechkin, T. Skovranek, YQ Chen,
%     B. M. Vinagre Jara, "Matrix approach to discrete
%     fractional calculus II: partial fractional differential
%     equations". http://arxiv.org/abs/0811.1355
% [2] R.G. Cooke, Infinite Matrices and Sequence Spaces,
%     MacMillan and Co., London, 1950. 347 pp.

S = eye(n);
r = sort(ROWS);
m = size(r,2);

for k = m:(-1):1
   if r(k)-1 == 0
        S = S((r(k)+1):(size(S,1)-k+1),:);
   elseif r(k) == size(S,1)
        S = S(1:(r(k)-1),:);
   else
        S = [S(1:(r(k)-1),:);  S((r(k)+1):(size(S,1)),:)];
   end
end


text:
function [ SimilarityMatrix ] = similarity_neighbor_grid_further( x, side, types, dist)
%similarity_neighbor_grid_further Summary of this function goes here

    % this assumes that the patch is laid out with first column, then second
    % column, ... final column (column major)

%     dist = 2;
    SimilarityMatrix = eye(side*side);

    % types - 1 - horizontal, 2 - vertical, 3 - diagonal (bl-tr), 4 -
    % diagonal (br - tl)
    for t=1:numel(types)

        if(types(t) == 1)

            % for horizontal we want to link both neighbours
            % (which are offset from the points by height)
            i = 1:(side*side-side*dist);
            % create the neighboring links for i
            SimilarityMatrix(sub2ind([side^2, side^2], i, i+side*dist)) = 1;
            SimilarityMatrix(sub2ind([side^2, side^2], i+side*dist, i)) = 1;

        end
        if(types(t) == 2)

            % for vertical we want to link both neighbours except at edge
            % cases which are mod(y_loc,side) = 0 as they are at the edges
            i = 1:side*side;
            i_to_rem =[];
            for s=1:dist
                i_to_rem = union(i_to_rem, i(mod(i+s-1, side) == 0));
            end
            i_both = setdiff(i, i_to_rem);
            % create the neighboring links for i
            SimilarityMatrix(sub2ind([side^2, side^2], i_both+dist, i_both)) = 1;
            SimilarityMatrix(sub2ind([side^2, side^2], i_both, i_both+dist)) = 1;

        end
        if(types(t) == 3)

            % for diagonal to top right, and bottom left don't use right most column
            i = 1:(side^2)-dist * side;
            i_to_rem = [];
            for s=1:dist
                i_to_rem = union(i_to_rem,i(mod(i-s, side) == 0));
            end
            i_both = setdiff(i, i_to_rem);
            % create the neighboring links for i
            SimilarityMatrix(sub2ind([side^2, side^2], i_both+dist*side-dist, i_both)) = 1;
            SimilarityMatrix(sub2ind([side^2, side^2], i_both, i_both+dist*side-dist)) = 1;
        end
        if(types(t) == 4)

            % for diagonal to top left, and bottom right don't use right most column
            i = 1:(side^2)-dist*side;
            i_to_rem = [];
            for s=1:dist
                i_to_rem = union(i_to_rem, i(mod(i+s-1, side) == 0));
            end
            i_both = setdiff(i, i_to_rem);
            % create the neighboring links for i
            SimilarityMatrix(sub2ind([side^2, side^2], i_both+dist*side+dist, i_both)) = 1;
            SimilarityMatrix(sub2ind([side^2, side^2], i_both, i_both+dist*side+ dist)) = 1;

        end

    end
    assert(isequal(SimilarityMatrix, SimilarityMatrix'));
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function dist = getDistPeopleHist(expidx)

p = rcnn_exp_params(expidx);

[annolistAll, imgidxs, rectidxs, rectIgnore, groupidxs]= getAnnolist(expidx);

distAll = [];
nPeople = [];
imgidxsAll = [];

for i = 1:length(imgidxs)
    fprintf('.');
    rect = annolistAll(imgidxs(i)).annorect;
    for g = 1:length(groupidxs{imgidxs(i)})
        ridxs_group = intersect(rectidxs{imgidxs(i)},groupidxs{imgidxs(i)}{g});
        if (~isempty(ridxs_group))
            dist = inf(length(ridxs_group),length(ridxs_group));
            for ridx1 = 1:length(ridxs_group)
                for ridx2 = 1:length(ridxs_group)
                    if ((ridx1 == ridx2) || ...
                            ~isfield(rect(ridxs_group(ridx1)), 'annopoints') || isempty(rect(ridxs_group(ridx1)).annopoints) || ...
                            ~isfield(rect(ridxs_group(ridx2)), 'annopoints') || isempty(rect(ridxs_group(ridx2)).annopoints))
                        continue;
                    end
                    dist(ridx1,ridx2) = util_get_min_dist(rect(ridxs_group(ridx1)),rect(ridxs_group(ridx2)));
                end
            end
            md = min(dist,[],2);
            distAll = [distAll; mean(md)];
            nPeople = [nPeople; length(ridxs_group)];
            imgidxsAll = [imgidxsAll; i];
        end
    end

    if (~mod(i, 100))
        fprintf(' %d/%d\n',i,length(imgidxs));
    end
end
fprintf(' done\n');

edges = 0:4:40;
[~,bin] = histc(distAll,edges);
dist = zeros(max(nPeople),length(edges));

for i = 1:length(nPeople)
    dist(nPeople(i),bin(i)+1) = dist(nPeople(i),bin(i)+1) + 1;
end

nPeopleHist = zeros(max(nPeople),1);
for i = 1:length(nPeopleHist)
    nPeopleHist(i) = sum(nPeople(nPeople == i));
end

figure(100); clf;

set(0,'DefaultAxesFontSize', 16)
set(0,'DefaultTextFontSize', 16)

bar3(dist);
set(gca,'XLim',[0 max(distAll)],'XTickLabel',edges);
ylabel('# people');
xlabel('avg min dist, px');
print(gcf, '-dpng', [p.plotsDir '/distPeople-expidx' num2str(expidx) '.png']);

figure(100); clf;
bar(nPeopleHist);
set(gca,'XTickLabel',1:length(nPeopleHist));
ylabel('# people');
text(0.5:length(nPeopleHist)-0.5,nPeopleHist+max(nPeopleHist)/50,num2str(nPeopleHist/sum(nPeopleHist)*100,'%1.1f'))
print(gcf, '-dpng', [p.plotsDir '/nPeople-expidx' num2str(expidx) '.png']);

    function minDist = util_get_min_dist(rect1,rect2)

        points1 = rect1.annopoints.point;
        points2 = rect2.annopoints.point;
        minDist = Inf;
        for idx1 = 1:length(points1)
            for idx2 = 1:length(points2)
                d = 1/rect1.scale*norm([points1(idx1).x points1(idx1).y] - [points2(idx2).x points2(idx2).y]);
                if (minDist > d)
                    minDist = d;
                end
            end
        end

    end

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function [] = fitKnownShapes(shapeModelOpt,statesDir,inferredShapesDir)

%% Setup optimization parameters
params = get_params();
relaxOpt = params.opt.relaxOpt;
lambda = params.opt.lambda;
maxiter = params.opt.testIters;
Sbar = shapeModelOpt.S;
V = shapeModelOpt.V;
N = shapeModelOpt.normals;
tri = shapeModelOpt.tri;
shapePriorInstance = shapeModelOpt.shapePriorInstance;

%% Read file names to fit basis shapes to
fnames = getFileNamesFromDirectory(statesDir,'types',{'.mat'});
fnames = removeFlipNames(fnames);

%% Fit basis shapes
pBar =  TimedProgressBar( length(fnames), 30, ...
    'Fitting Basis shapes: Remaining ', ', Completed ', 'Fitted Basis Shapes in: ' );

parfor i=1:length(fnames)
    optimizeTestShape(fullfile(statesDir,fnames{i}),lambda,Sbar,V,N,tri,...
        shapePriorInstance,maxiter,fullfile(inferredShapesDir,fnames{i}),relaxOpt,params);
    pBar.progress();
end
pBar.stop();

end

%% Helper function to actually do the optimization
function [] = optimizeTestShape(fname,lambda,Sbar,V,N,tri,shapePriorInstance,maxiter,saveFile,relaxOpt,params)
%% PARAMS
stepsize = params.opt.testStepSize;
maxDelta = params.opt.maxDelta(3);

stateFile = fname;
state = load(stateFile);state = state.state;
K = size(V,2);
alpha = zeros(K,1);
numpoints = size(Sbar,1);

if(exist(saveFile,'file'))
    return;
end
%% iterations to learn alpha
initRot = state.cameraRot;
for iter=1:maxiter
    stepScale = 1;
    stepTr = 0.1;
    stepRot = 0.002;
    rotPriorStep = 0.1;

    if(det(state.cameraRot)>1.0001 || det(state.cameraRot)<0.9999)
        q = dcm2quat(state.cameraRot);
        state.cameraRot = quat2dcm(q/norm(q)); %wish I didn't have to do this - but such is life !
    end

    S = Sbar + reshape(V*alpha(:,1),numpoints,3);

    % Visualization
    if(mod(iter,10)==-1)
         vertices = S;
         subplot(2,2,1);
         verticesP = (state.cameraRot*vertices')';
         plot3(verticesP(:,1),verticesP(:,2),verticesP(:,3),'r.');
         axis equal;axis vis3d;view(0,-90);
         %disp(state.cameraRot);
         xlabel('x');ylabel('y');zlabel('z');

         subplot(2,2,2);
         p2 = transform2d(vertices,state.cameraRot,state.cameraScale,state.translation);
         imshow(color_seg(state.mask,state.im)); hold on; axis equal;axis off;axis vis3d;
         plot(p2(:,1),p2(:,2),'r.');

         subplot(2,2,3);
         imagesc(state.im);axis equal;axis off;axis vis3d;
         if(iter > maxiter - 10)
             disp('done');
             pause();close all;
         else
             pause(0.02);
         end
         clf
    end

    %% Image projection of points
    p2d = transform2d(S,state.cameraRot,state.cameraScale,state.translation);
    p2d = round(p2d);

    %% Occlusion Reasoning
    [occGrad,~] = occlusionGradOptimal(S,state,p2d);
    occGrad = lambda(1)*occGrad;

    %% Silhouette Consistency

    if(mod(iter-1,20)==0 || sum(ismember({'scale','translation'},relaxOpt)))
        % Finding projected points
        tree = vl_kdtreebuild(p2d');
    end
    silGrad = lambda(2)*silhouetteGrad(S,state,params, tree,p2d); %we also keep track of how many constraints each point violates
    silGrad = silGrad*500; % because we divide by number of silhouette points

    %% Keypoint Gradient
    if(ismember('kps',relaxOpt) && ~isempty(state.kps))
        kpGrad = lambda(6)*keypointGrad(S,state,params);
    else
        kpGrad = zeros(size(S));
    end

    %% Priors
    gradPriorAlpha = -(lambda(4)*numpoints)*alpha;
    shapePriorGrad = lambda(5)*shapePriorInstance(S);

    %% translation and scale gradients
    gradS2d = occGrad + silGrad + kpGrad;
    gradP3d = state.cameraScale*state.cameraRot*gradS2d';
    gradP2d = gradP3d(1:2,:);
    gradTr = mean(gradP2d,2);
    Sbar2d = bsxfun(@minus,p2d,mean(p2d,1));Sbar2d = Sbar2d';
    gradScale = mean(gradP2d(:).*Sbar2d(:))/state.cameraScale;

    %% Rotation gradients (experimental)
    gradRot = gradP3d*S/numpoints;
    %R = R_0 * expm(delta) = R_0(I+delta)
    gradDelta = state.cameraRot*gradRot;
    gradSkewDelta = (gradDelta-gradDelta')/2;
    deltaPrior = real(logm(state.cameraRot'*initRot));
    deltaPrior(abs(deltaPrior)<eps) = 0;
    normTwist = norm(gradSkewDelta,'fro');

    %% Deformation weight gradients
    gradShapePriorAlpha(:,1) = V'*shapePriorGrad(:);
    gradSilAlpha = V'*silGrad(:);
    gradOcclusionAlpha = V'*occGrad(:);
    gradKeypointAlpha = V'*kpGrad(:);
    gradAlpha = gradOcclusionAlpha + gradPriorAlpha + gradSilAlpha + gradShapePriorAlpha + gradKeypointAlpha; % + gradICPAlpha

    %% Updating transformation parameters
    step = min(maxDelta/max(max(abs(gradAlpha))),stepsize);
    if(iter > maxiter/2)
        alpha = alpha+step*gradAlpha;
    end
    stepScale = min(stepScale,0.02*state.cameraScale/abs(gradScale));
    if(ismember('scale',relaxOpt))
        state.cameraScale = state.cameraScale + gradScale*stepScale;
    end
    stepTr = min(stepTr,2/max(abs(gradTr)));
    if(ismember('translation',relaxOpt))
        state.translation = state.translation + gradTr*stepTr;
    end
    if(ismember('rotation',relaxOpt))
        stepRot = min(stepRot,pi/(50*normTwist));
        %stepRot = min(stepRot,pi/(25*normTwist));
        delta = expm(stepRot*gradSkewDelta + deltaPrior*rotPriorStep);
        if(det(delta) ~= 1)
            %disp('wtf');
	        %keyboard;
            %delta(1,:) = delta(1,:)/norm(delta(1,:));
            %delta(2,:) = delta(2,:)/norm(delta(2,:));
            %delta(3,:) = delta(3,:)/norm(delta(3,:));
        end
        state.cameraRot = state.cameraRot*delta;
    end
end

S = Sbar + reshape(V*alpha,numpoints,3);

%% Saving
state.alpha = alpha;
state.normals = N;
state.tri = tri;
state.S = S;
save(saveFile,'state');

end

repo: CategoryShapes


text:
function theta = homogeneous_to_full_GLG( L, mu, sigma, alpha, beta, kappa )

% Compute full parameter set from parameters from INLA
% The homogeneus GLG model used in INLA only contains five parameters; this
% function computes the remaining parameters in the 'full' GLG model.
%
% Syntax:
%   theta = homogeneus_to_full_GLG( L, mu, sigma, alpha, beta, kappa )
%
% Input:
%   L     : The number of levels in the transform
%
%   mu    : Mean of root node
%
%   sigma : Standard dev of root node
%
%   alpha : Intercept of transitions
%
%   beta  : Slope of transitions
%
%   kappa : Standard dev of transitions
%
% The dimension of the parameters determines the number of dimensions.
% If a parameters is empty, a value is picked at random.
%
%
% Output:
%   theta : A matrix with the full parameter set

% ----------------------------------------------------------------------
%                                               Fill in empty parameters
% ----------------------------------------------------------------------

if nargin == 1
    [mu, sigma, alpha, beta, kappa] = deal( [] );
end

parameters = {mu, sigma, alpha, beta, kappa};
dims = cellfun(@numel, parameters);

missing_vals = cellfun(@isempty, parameters);

% The default number of dimensions is 3
if all(missing_vals)
    D = 3;
else
    D = dims( find(dims, 1) );
end

% Simulate missing values
for j = find(missing_vals)
    switch j
        case {1, 3}
            parameters{j} = -1 + 2*rand(D, 1);

        case {2, 4, 5}
            parameters{j} = rand(D, 1);
    end
end

[mu, sigma, alpha, beta, kappa] = deal( parameters{:} );

% ----------------------------------------------------------------------
%                                                     Compute full theta
% ----------------------------------------------------------------------

theta = zeros(L, 5, D);

theta(1, 1, :) = mu;
theta(1, 2, :) = sigma;

% All alpha's, beta's and kappa's are the same
theta(2:end, 3, :) = repmat( reshape(alpha, [1 1 D]), [L-1 1 1] );
theta(2:end, 4, :) = repmat( reshape(beta, [1 1 D]), [L-1 1 1] );
theta(2:end, 5, :) = repmat( reshape(kappa, [1 1 D]), [L-1 1 1] );

% Compute remaining mu's and sigma's
for l = 2:L
    for d = 1:D
        if beta(d) == 1
            theta(l,1,d) = alpha(d)*(l-1) + mu(d);

            sigma2 = kappa(d)^2*l + sigma(d)^2;
            theta(l,2,d) = sqrt( sigma2 );
        else
            theta(l,1,d) = alpha(d)*(beta(d)^(l-1)-1)/(beta(d)-1) + beta(d)^(l-1)*mu(d);

            sigma2 = kappa(d)^2*(beta(d)^(2*(l-1))-1)/(beta(d)^2-1) + beta(d)^(2*(l-1))*sigma(d)^2;
            theta(l,2,d) = sqrt( sigma2 );
        end
    end
end

end


text:
function A11 = synA11Qmax(A00, sizeA11, cmin)
%-----------------------------------------------------------------------------
%
% For each point of colour 11 this function assigns the maximum value at the
% neighbouring gridpoints of colour 00.
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: December 12, 2001.
% (c) 1998-2002 Stichting CWI, Amsterdam
%-----------------------------------------------------------------------------
[n00, m00]=size(A00);
if nargin == 3
  o=[0 0];
  if ~all(size(o) == size(sizeA11))
    error(' synA11Qmax - unexpected dimensions of sizeA11 ')
  else
    clear o;
    n11=sizeA11(1);
    m11=sizeA11(2);
  end
elseif nargin == 2
  n11=n00-1;
  m11=m00-1;
else
  error(' synA11Qmax - number of arguments should be either 2 or 3 ')
end
%[n11, m11]=size(A11);
if m11 == m00-1
  if n11 == n00-1
    S = max(stripL(A00), stripR(A00));
    A11=max(stripD(S), stripU(S));
  elseif n11 == n00
    S = max(stripL(A00), stripR(A00));
    A11=max(S, stripU(extD(S, cmin)));
  else
    disp([' size A00 = ' int2str(size(A00)) ' size A11 = ' int2str([n11 m11])]);
    error(' synA11Qmax - A00 and target A11 do not match ');
  end
elseif m11 == m00
  if n11 == n00-1
    S = max(stripL(extR(A00, cmin)), A00);
    A11=max(stripD(S), stripU(S));
  elseif n11 == n00
    S = max(stripL(extR(A00, cmin)), A00);
    A11=max(S, stripU(extD(S, cmin)));
  else
    disp([' size A00 = ' int2str(size(A00)) ' size A11 = ' int2str([n11 m11])]);
    error(' synA11Qmax - A00 and target A11 do not match ');
  end
else
  disp([' size A00 = ' int2str(size(A00)) ' size A11 = ' int2str([n11 m11])]);
  error(' synA11Qmax - A00 and target A11 do not match ');
end
%------------------------------------------------------------------------------


text:
function I = esvm_show_top_exemplar_dets(test_struct, test_set, ...
                                         models, index)
% function I = esvm_show_top_exemplar_dets(test_struct, test_set, ...
%                                          models, index)
% Given the results of a detector, show the top exemplar detectiosn
% for exemplar id=[index]. Needs test_set to be able to load images
% for icons
%
% Copyright (C) 2011-12 by Tomasz Malisiewicz
% All rights reserved.
%
% This file is part of the Exemplar-SVM library and is made
% available under the terms of the MIT license (see COPYING file).
% Project homepage: https://github.com/quantombone/exemplarsvm

bbs = cat(1,test_struct.unclipped_boxes{:});
[aa,bb] = sort(bbs(:,end),'descend');
bbs = bbs(bb,:);

m = models{index};
m.model.svbbs = bbs;

m.model.svbbs = m.model.svbbs(m.model.svbbs(:,6)==index);
try
  m.model = rmfield(m.model,'svxs');
catch
end
m.train_set = test_set;

I = esvm_show_det_stack(m,4,4);

author: quantombone

repo: exemplarsvm


text:
function artifact = trl2artifact(trl)

% TRL2ARTIFACT converts between two representations of events or trials.
%
% FieldTrip uses a number of representations for events that are conceptually very similar
%   event    = structure with type, value, sample, duration and offset
%   trl      = Nx3 numerical array with begsample, endsample, offset
%   trl      = table with 3 columns for begsample, endsample, offset
%   artifact = Nx2 numerical array with begsample, endsample
%   artifact = table with 2 columns for begsample, endsample
%   boolvec  = 1xNsamples boolean vector with a thresholded TTL/trigger sequence
%   boolvec  = MxNsamples boolean matrix with a thresholded TTL/trigger sequence
%
% If trl or artifact are represented as a MATLAB table, they can have additional
% columns. These additional columns have to be named and are not restricted to
% numerical values.
%
% See also ARTIFACT2BOOLVEC, ARTIFACT2EVENT, ARTIFACT2TRL, BOOLVEC2ARTIFACT, BOOLVEC2EVENT, BOOLVEC2TRL, EVENT2ARTIFACT, EVENT2BOOLVEC, EVENT2TRL, TRL2ARTIFACT, TRL2BOOLVEC, TRL2EVENT

% Copyright (C) 2009, Ingrid Nieuwenhuis
% Copyright (C) 2020, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% trl and artifact are similar, except for the offset column
if isnumeric(trl)
  % do not keep any additional columns
  artifact = trl(:,1:2);
elseif istable(trl)
  % remove the offset, the additional columns can remain
  artifact = trl(:, [1 2 4:end]);
elseif iscell(trl)
  % there are multiple types of trials/artifacts, use recursion to loop over them
  artifact = cell(size(trl));
  for i=1:numel(trl)
    artifact{i} = artifact2trl(trl{i});
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% Import data from spreadsheet
% Script for importing data from the following spreadsheet:
%
%    Workbook: /Users/gt/Dropbox/KSC_rocket_launches_2016_onwards_v3.xlsx
%    Worksheet: Sheet1
%
% To extend the code for use with different selected data or a different
% spreadsheet, generate a function instead of a script.

% Auto-generated by MATLAB on 2017/11/27 16:46:08

%% Import the data
[~, ~, raw] = xlsread('/Users/gt/Dropbox/KSC_rocket_launches_2016_onwards_v3.xlsx','Sheet1');
raw = raw(2:end,:);
raw(cellfun(@(x) ~isempty(x) && isnumeric(x) && isnan(x),raw)) = {''};
stringVectors = string(raw(:,[3,4,5,8,13]));
stringVectors(ismissing(stringVectors)) = '';
raw = raw(:,[1,2,6,7,10,11,12,13]);

%% Replace non-numeric cells with NaN
R = cellfun(@(x) ~isnumeric(x) && ~islogical(x),raw); % Find non-numeric cells
raw(R) = {NaN}; % Replace non-numeric cells

%% Create output variable
data = reshape([raw{:}],size(raw));

%% Create table
cols = table;

%% Allocate imported array to column variable names
cols.DATE = data(:,1);
cols.TIME = data(:,2);
cols.COMPANY = categorical(stringVectors(:,1));
cols.ROCKETTYPE = categorical(stringVectors(:,2));
cols.PAYLOAD = stringVectors(:,3);
cols.MASSt = data(:,3);
cols.SLC = data(:,4);
cols.ORBIT = categorical(stringVectors(:,4));
cols.RECORDED = data(:,5);
cols.SPACEX = data(:,6);
cols.ALL = data(:,7);
cols.OTHER = data(:,8);

%% Clear temporary variables
clearvars data raw stringVectors R;

%% plot
dnum=cols.DATE-1 + datenum(1899,12,31) + cols.TIME;
figure
plot(dnum,cols.ALL,'b^-',dnum,cols.SPACEX,'g^-',dnum,cols.RECORDED,'r^-',dnum,cols.OTHER,'r*-');
datetick('x','keeplimits')
ylabel('Cumulative number')
xlabel('Date')
legend('Launches - any', 'Launches - SpaceX', 'Launches recorded', 'Other events recorded')

author: geoscience-community-codes

repo: GISMO


text:
function x = invpowspec(y, sr, wintime, steptime, excit)
%x = invpowspec(y, sr, wintime, steptime, excit)
%
% Attempt to go back from specgram-like powerspec to audio waveform
% by scaling specgram of white noise
%
% default values:
% sr = 8000Hz
% wintime = 25ms (200 samps)
% steptime = 10ms (80 samps)
% which means use 256 point fft
% hamming window
%
% excit is input excitation; white noise is used if not specified

% for sr = 8000
%NFFT = 256;
%NOVERLAP = 120;
%SAMPRATE = 8000;
%WINDOW = hamming(200);

[nrow, ncol] = size(y);

if nargin < 2
  sr = 8000;
end
if nargin < 3
  wintime = 0.025;
end
if nargin < 4
  steptime = 0.010;
end
if nargin < 5
  r = [];
else
  r = excit;
end

winpts = round(wintime*sr);
steppts = round(steptime*sr);

NFFT = 2^(ceil(log(winpts)/log(2)));

if NFFT ~= 2*(nrow-1)
  disp('Inferred FFT size doesn''t match specgram');
end

NOVERLAP = winpts - steppts;
SAMPRATE = sr;

% Values coming out of rasta treat samples as integers,
% not range -1..1, hence scale up here to match (approx)
%y = abs(specgram(x*32768,NFFT,SAMPRATE,WINDOW,NOVERLAP)).^2;

xlen = winpts + steppts*(ncol - 1);

if length(r) == 0
  r = randn(xlen,1);
end
r = r(1:xlen);

R = specgram(r/32768/12, NFFT, SAMPRATE, winpts, NOVERLAP);
R = R .* sqrt(y);;
x = ispecgram(R, NFFT, SAMPRATE, winpts, NOVERLAP);

author: stephencwelch

repo: Perceptual-Coding-In-Python


text:
function [adfreq, n, ts, fn, ad] = plx_ad_v(filename, channel)
% plx_ad_v(filename, channel): read all a/d data for the specified channel from a .plx or .pl2 file
%              returns a/d values in millivolts
%
% [adfreq, n, ts, fn, ad] = plx_ad_v(filename, channel)
% [adfreq, n, ts, fn, ad] = plx_ad_v(filename, 0)
% [adfreq, n, ts, fn, ad] = plx_ad_v(filename, 'FP01')
%
% INPUT:
%   filename - if empty string, will use File Open dialog
%   channel - 0-based channel number or channel name
%
%           a/d data come in fragments. Each fragment has a timestamp
%           and a number of a/d data points. The timestamp corresponds to
%           the time of recording of the first a/d value in this fragment.
%           All the data values stored in the vector ad.
%
% OUTPUT:
%   adfreq - digitization frequency for this channel
%   n - total number of data points
%   ts - array of fragment timestamps (one timestamp per fragment, in seconds)
%   fn - number of data points in each fragment
%   ad - array of a/d values converted to millivolts

adfreq = 0;
n = 0;
ts = -1;
fn = -1;
ad = -1;

if nargin ~= 2
    error 'Expected 2 input arguments';
end

[ filename, isPl2 ] = internalPL2ResolveFilenamePlx( filename );
if isPl2 == 1
    % return data from .pl2 file
    pl2Channel = channel;
    % PL2 uses 1-based channel numbers
    if ischar(channel) == 0
        pl2Channel = channel+1;
    end
    pl2ad = PL2Ad(filename, pl2Channel);
    if numel(pl2ad.Values) > 0
        adfreq = pl2ad.ADFreq;
        n = numel(pl2ad.Values);
        ts = pl2ad.FragTs;
        fn = pl2ad.FragCounts;
        ad = pl2ad.Values;
    end
    return;
end

channelNumber = plx_ad_resolve_channel(filename, channel);
if channelNumber == -1
    fprintf('\n plx_ad_v: no header for the specified A/D channel.');
    fprintf('\n    use plx_ad_info(filename) to print the list of a/d channels\n');
    return
end

[adfreq, n, ts, fn, ad] = mexPlex(17, filename, channelNumber);

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function output = callfiltersdsp(model)

options = model.options.filtersd;

model = yalmip2nonlinearsolver(model);

if ~model.derivative_available
    disp('Derivate-free call to filterSD not yet implemented')
    error('Derivate-free call to filterSD not yet implemented')
end
if model.options.savedebug
    save filtersddebug model
end
showprogress('Calling filterSDSP',model.options.showprogress);

cu = [ repmat(0,length(model.bnonlinineq),1);
    repmat(0,length(model.bnonlineq),1);
    repmat(0,length(model.b),1);
    repmat(0,length(model.beq),1)];

cl = [ repmat(-inf,length(model.bnonlinineq),1);
    repmat(0,length(model.bnonlineq),1);
    repmat(-inf,length(model.b),1);
    repmat(0,length(model.beq),1)];

if isempty(cu)
    Flow = [];
    Fupp = [];
end

% These are needed to avoid recomputation due to ipopts double call to get
% f and df, and g and dg
global latest_x_f
global latest_x_g
global latest_df
global latest_f
global latest_G
global latest_g
global latest_xevaled
global latest_x_xevaled
latest_G = [];
latest_g = [];
latest_x_f = [];
latest_x_g = [];
latest_xevaled = [];
latest_x_xevaled = [];

funcs.objective = @(x)ipopt_callback_f(x,model);
funcs.gradient = @(x)ipopt_callback_df(x,model);
if ~isempty(cu)
    funcs.constraints = @(x)ipopt_callback_g(x,model);
    funcs.jacobian  = @(x)ipopt_callback_dg(x,model);
else
    funcs.constraints = [];
    funcs.jacobian = [];
end

lb  = model.lb(:)';
ub = model.ub(:)';

if ~isempty(cu)
   nljacstr = jacobiansparsityfromnonlinear(model);
else
   nljacstr = [];
end

if ~model.options.warmstart
    model.x0 = (lb+ub)/2;
    model.x0(isinf(ub)) = lb(isinf(ub))+1;
    model.x0(isinf(lb)) = ub(isinf(lb))-1;
    model.x0(isinf(model.x0)) = 0;
end

options.display = model.options.verbose;
solvertime = tic;
[xout,fval,exitflag,stats,lambda] = filtersdsp(funcs.objective, funcs.gradient, model.x0, lb, ub, funcs.constraints, funcs.jacobian,nljacstr, cl, cu, options);
solvertime = toc(solvertime);

% Duals currently not supported
lambda = [];

x = RecoverNonlinearSolverSolution(model,xout);

switch exitflag
    case {0}
        problem = 0;
    case {1}
        problem = 2;
    case {2,3}
        problem = 1;
    case {5}
        problem = 3;
    case {4,9}
        problem = 4;
    case 105
        problem = 16;
    otherwise
        problem = -1;
end

% Internal format for duals
D_struc = [];

% Save all data sent to solver?
if model.options.savesolverinput
    solverinput.model = model;
else
    solverinput = [];
end

% Save all data from the solver?
if model.options.savesolveroutput
    solveroutput.x = xout;
    solveroutput.exitflag = exitflag;
    solveroutput.stats = stats;
else
    solveroutput = [];
end

% Standard interface
output = createOutputStructure(x,D_struc,[],problem,'filterSD',solverinput,solveroutput,solvertime);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [scoresAll, labelsAll, nGTall] = assignGTMulticutMulti_WAF(annolistWAF_det,annolistWAF_gt,pidxsAll,parts,thresh_gt)

if (nargin < 5)
    thresh_gt = 0.5;
end

% scoresAll = cell(length(pidxsAll),1);
% labelsAll = cell(length(pidxsAll),1);
% nGTall = zeros(length(pidxsAll),length(annolist));
%
% for i = 1:length(pidxsAll)
%     scoresAll{i} = cell(length(annolist),1);
%     labelsAll{i} = cell(length(annolist),1);
% end

assert(length(annolistWAF_det) == length(annolistWAF_gt));
bVis = false;

for imgidx = 1:length(annolistWAF_gt)

    fprintf('.');

%     dist = inf(length(keypointsAll(imgidx).det),length(annolist(imgidx).annorect),length(pidxsAll));
%     score = nan(length(keypointsAll(imgidx).det),length(annolist(imgidx).annorect),length(pidxsAll));
%     hasDet = false(length(keypointsAll(imgidx).det),length(annolist(imgidx).annorect),length(pidxsAll));
%     hasGT = false(length(keypointsAll(imgidx).det),length(annolist(imgidx).annorect),length(pidxsAll));

    for j = 1:length(annolistWAF_det(imgidx).stickmen) % clusters

        for ridx = 1:length(annolistWAF_gt(imgidx).stickmen) % gt people

            rect = annolist(imgidx).annorect(ridx);
            refDist = util_get_head_size(rect);
            points_gt = rect.annopoints.point;

            for i = 1:length(pidxsAll)
                pidx = pidxsAll(i);
                % part is a joint
                assert(parts(pidx+1).pos(1) == parts(pidx+1).pos(2));
                jidx = parts(pidx+1).pos(1);

                det = keypointsAll(imgidx).det{j}{jidx+1};
                if (~isempty(det))
                    point_det = det(1:2);
                    score(j,ridx,i) = det(3);
                    hasDet(j,ridx,i) = true;
                end
                pp = util_get_annopoint_by_id(points_gt, jidx);
                if (~isempty(pp))
                    hasGT(j,ridx,i) = true;
                end
                if (hasDet(j,ridx,i) && hasGT(j,ridx,i))
                    dist(j,ridx,i) = norm([pp.x pp.y] - point_det)/refDist;
                end
            end
        end
    end

    match = dist <= thresh_gt;
    pck = sum(match,3)./sum(hasGT,3);
    [val,idx] = max(pck,[],2);
    for j = 1:length(idx)
        pck(j,setdiff(1:size(pck,2),idx(j))) = 0;
    end
    [val,clusToGT] = max(pck,[],1);
    clusToGT(val == 0) = 0;

    for j = 1:length(keypointsAll(imgidx).det)
        if (ismember(j,clusToGT)) % match to GT
            ridx = find(clusToGT == j);

            if (bVis)
                rect = annolist(imgidx).annorect(ridx);
                points_gt = rect.annopoints.point;
                gt = cell(16,1);
                for i = 1:length(pidxsAll)
                    pidx = pidxsAll(i);
                    % part is a joint
                    assert(parts(pidx+1).pos(1) == parts(pidx+1).pos(2));
                    jidx = parts(pidx+1).pos(1);
                    pp = util_get_annopoint_by_id(points_gt, jidx);
                    if (~isempty(pp))
                        gt{jidx+1} = [pp.x pp.y -1];
                    end
                end

                img = imread(annolist(imgidx).image.name);
                figure(100); clf; subplot(1,2,1); imagesc(img); axis equal; hold on;
                vis_multicut2(keypointsAll(imgidx).det{j},pidxsAll,parts);
                subplot(1,2,2); imagesc(img); axis equal; hold on;
                vis_multicut2(gt,pidxsAll,parts);
            end

            s = squeeze(score(j,ridx,:));
            m = squeeze(match(j,ridx,:));
            hd = squeeze(hasDet(j,ridx,:));

            idxs = find(hd);
            for i = 1:length(idxs)
                scoresAll{idxs(i)}{imgidx} = [scoresAll{idxs(i)}{imgidx};s(idxs(i))];
                labelsAll{idxs(i)}{imgidx} = [labelsAll{idxs(i)}{imgidx};m(idxs(i))];
            end

        else % no matching to GT
            s = squeeze(score(j,1,:));
            m = false(size(match,3),1);
            hd = squeeze(hasDet(j,ridx,:));
            idxs = find(hd);
            for i = 1:length(idxs)
                scoresAll{idxs(i)}{imgidx} = [scoresAll{idxs(i)}{imgidx};s(idxs(i))];
                labelsAll{idxs(i)}{imgidx} = [labelsAll{idxs(i)}{imgidx};m(idxs(i))];
            end
        end
    end

    for ridx = 1:length(annolist(imgidx).annorect)
        hg = squeeze(hasGT(1,ridx,:));
%         idxs = find(hg);
%         for i = 1:length(idxs)
%             nGTall(idxs(i),imgidx) = nGTall(idxs(i),imgidx) + hg(idxs(i));
%         end
        nGTall(:,imgidx) = nGTall(:,imgidx) + hg;
    end

    if (~mod(imgidx, 100))
        fprintf(' %d/%d\n',imgidx,length(keypointsAll));
    end
end
fprintf(' done\n');
end


text:
function res = Fspikernel(seqS,seqT,N,lam,mu,p)

%spikernel kernel function
%
%   This code implements the spikernel function, defined in:
%   Spikernels: Embedding Spiking Neurons in Inner-Product Spaces
%   Lavi Shpigelman, Yoram Singer, Rony Paz and Eilon Vaadia
%   Advances in Neural Information Processing Systems (NIPS) 15
%   MIT Press, Cambridge, MA, 2003.
%
%   res= Fspikernel (seqS,seqT,N,lam,mu,p)
%
%   returns res = K(seqS,seqT) = sum_{i=1}^N p^i K_i(seqS,seqT) where the kernel parameters are:
%
%   Inputs
%       seqS and seqT are the two sequences to be compared.
%           Each row of these sequences is one 'letter' in the sequence (must be of same length)
%           Their column lengths are the sequence lengths (can be of different length).
%       N = max subsequence lengths (the kernel compares subsequences of lengths 1 to N and returns a sum of kernels)
%       lam=\lambda parameter from article
%       mu=\mu parameter from article
%       p= the q prarameter in the article - the parameter that is used in the sum of kernels to weigh them differently
%
%    This implementation assumes that the function d(x,y) is the squared \ell_2 norm
%

% ------------- parameter check -------------------------
lens=size(seqS,1);
lent=size(seqT,1);
q=size(seqS,2);

if (q~=size(seqT,2))
    error('the dimesions of the sequence letters in seqS and seqT must be the same')
end
if (length(N)~=1)
    error('N should be a scalar value');
end
if (length(lam)~=1)
    error('lam should be a scalar value');
end
if (length(mu)~=1)
    error('mu should be a scalar value');
end
if (length(p)~=1)
    error('p should be a scalar value');
end
% --------------------------------------------------------

dmu = mu.^(0.5*( repmat( sum( seqS.^2, 2), 1, lent) + ...
              ( repmat( sum(seqT.^2, 2), 1, lens))' - 2*seqS*seqT'))';

res = fspike( N, p, lam, dmu);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore