text:
function dx = f_Wrapper(u,x,p,parIdx)
% dx = f(u,x,p)
    if coder.target('MATLAB')
        % Specify your own f(u,x,p) function for normal excution
    else
        % Specify your own f(u,x,p) function for code generation
        coder.cinclude('iiwa14.h');
        q = x(1:7,1);
        qd = x(8:end,1);
        qdd = zeros(7,1);
        tau = u(1:7,1);
        coder.ceval('qdd_cal',  coder.ref(q),...
                                coder.ref(qd),...
                                coder.ref(qdd),...
                                coder.ref(tau),...
                                parIdx);
        dx  = [qd;qdd];
    end

end

author: deng-haoyang

repo: ParNMPC


text:
classdef NodesConnector < handle

    properties (Access = private)
        nodes
        coord
        div
        theta
    end

    properties (Access = public)
        connec
    end

    methods (Access = public)

        function obj = NodesConnector(cParams)
            obj.init(cParams);
        end

        function computeConnections(obj)
            obj.connec = delaunay(obj.coord);
            if obj.nodes.vert == 6
                obj.deleteExtraElementsCaseA();
            elseif obj.theta(2) ~= 90
                obj.deleteExtraElementsCaseB();
            end
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.nodes = cParams.nodes;
            obj.coord = cParams.coord;
            obj.theta = cParams.theta;
            obj.div = cParams.div;
        end

        function deleteExtraElementsCaseA(obj)
            cont = 1;
            rowsToDelete = [];
            for i = 1:size(obj.connec,1)
                found = zeros(1,3);
                elementNodes = obj.connec(i,:);
                for j = 1:size(obj.connec,2)
                    if elementNodes(j) < obj.nodes.bound+1
                        found(j) = found(j)+1;
                    end
                end
                if sum(found) == 3
                    rowsToDelete(cont,1) = i;
                    cont = cont+1;
                end
            end
            obj.connec(rowsToDelete,:) = [];
        end

        function deleteExtraElementsCaseB(obj)
            cont = 1;
            rowsToDelete = [];
            for i = 1:size(obj.connec,1)
                found = zeros(1,3);
                noDeletionA = zeros(1,3);
                noDeletionB = zeros(1,3);
                noDeletionC = zeros(1,3);
                noDeletionD = zeros(1,3);
                elementNodes = obj.connec(i,:);
                for j = 1:size(obj.connec,2)
                    if (elementNodes(j) < obj.nodes.bound+1)
                        found(j) = found(j)+1;
                        if (elementNodes(j) == 1) || (elementNodes(j) == obj.nodes.vert+1) || (elementNodes(j) == obj.nodes.bound)
                            noDeletionA(j) = noDeletionA(j)+1;
                        end
                        if (elementNodes(j) == 2) || (elementNodes(j) == obj.nodes.vert+obj.div(1)-1) || (elementNodes(j) == obj.nodes.vert+obj.div(1))
                            noDeletionB(j) = noDeletionB(j)+1;
                        end
                        if (elementNodes(j) == 3) || (elementNodes(j) == obj.nodes.vert+obj.div(1)+obj.div(2)-2) || (elementNodes(j) == obj.nodes.vert+obj.div(1)+obj.div(2)-1)
                            noDeletionC(j) = noDeletionC(j)+1;
                        end
                        if (elementNodes(j) == 4) || (elementNodes(j) == obj.nodes.vert+2*obj.div(1)+obj.div(2)-3) || (elementNodes(j) == obj.nodes.vert+2*obj.div(1)+obj.div(2)-2)
                            noDeletionD(j) = noDeletionD(j)+1;
                        end
                    end
                end
                if (sum(found) == 3) && (sum(noDeletionA) ~= 3) && (sum(noDeletionB) ~= 3) && (sum(noDeletionC) ~= 3) && (sum(noDeletionD) ~= 3)
                    rowsToDelete(cont,1) = i;
                    cont = cont+1;
                end
            end
            obj.connec(rowsToDelete,:) = [];
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function sparse_grid_test01 ( dim_min, dim_max, level_max_min, ...
  level_max_max )

%*****************************************************************************80
%
%% TEST01 tests SPARSE_GRID_OFN_SIZE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 December 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_MIN, the minimum spatial dimension to consider.
%
%    Input, integer DIM_MAX, the maximum spatial dimension to consider.
%
%    Input, integer LEVEL_MAX_MIN, the minimum value of LEVEL_MAX to consider.
%
%    Input, integer LEVEL_MAX_MAX, the maximum value of LEVEL_MAX to consider.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01\n' );
  fprintf ( 1, '  SPARSE_GRID_OFN_SIZE returns the number of\n' );
  fprintf ( 1, '  distinct points in a sparse grid made up of all \n' );
  fprintf ( 1, '  product grids formed from completely nested open \n' );
  fprintf ( 1, '  quadrature rules.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The sparse grid is the sum of all product grids\n' );
  fprintf ( 1, '  of order LEVEL, with\n' );
  fprintf ( 1, '    0 <= LEVEL <= LEVEL_MAX.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  LEVEL is the sum of the levels of the 1D rules,\n' );
  fprintf ( 1, '  the order of the 1D rule is 2^(LEVEL+1) - 1,\n' );
  fprintf ( 1, '  the region is [-1,1]^DIM_NUM.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  For this kind of rule, there is complete nesting,\n' );
  fprintf ( 1, '  that is, a sparse grid of a given level includes\n' );
  fprintf ( 1, '  ALL the points on grids of lower levels.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   DIM: ' );

  for dim_num = dim_min : dim_max
    fprintf ( 1, '  %10d', dim_num );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   LEVEL_MAX\n' );
  fprintf ( 1, '\n' );

  for level_max = level_max_min : level_max_max
    fprintf ( 1, '    %4d', level_max );
    for dim_num = dim_min : dim_max
      point_num = sparse_grid_ofn_size ( dim_num, level_max );
      fprintf ( 1, '  %10d', point_num );
    end
    fprintf ( 1, '\n' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ consumption ] = infer_consumption(result, evaluation_and_training_days, usage_duration, appliance, setup)

    %% infer consumption from average consumption of corresponding appliance
    num_measurements = size(evaluation_and_training_days{1},1) * 86400 / setup.granularity;
    consumption = zeros(1, num_measurements);

    %% Fridge, Freezer, TV, Stereo
    if strcmp(appliance, 'Fridge') == 1 ...
            || strcmp(appliance, 'Freezer') == 1 ...
            || strcmp(appliance, 'TV') == 1 ...
            || strcmp(appliance, 'Stereo') == 1

        warning('off');
        training_days = evaluation_and_training_days{2};
        [avg_runtime, plug_events] = get_average_runtime_from_plug_data(appliance, training_days, setup);
        average_power = get_average_power_from_plug_data(appliance, training_days, setup, plug_events);
        %% HERE

        times = result.events(:,1);
        on_events = result.events(:,3) > 0;
        times_on_events = times(on_events);
        times_off_events = times(~on_events);
        events_to_skip = [];
        for i = 1:length(times)
            if any(events_to_skip == i)
                continue;
            end
            % is on event
            % on event: check if off_event is within runtime
            % yes: on_event:off_event, increase by off_event
            % no: on_event:on_event+runtime, increase by runtime [ check boundaries ]
            if on_events(i) == 1
                switch_on = times(i);
                potential_off_events = times_off_events > switch_on & times_off_events < switch_on + avg_runtime;
                if any(potential_off_events)
                    off_idx = min(find(potential_off_events));
                    switch_off = times_off_events(off_idx);
                    consumption(switch_on:switch_off) = average_power;
                    % skip off event
                    idx_off_event = find(times == switch_off);
                    events_to_skip = [events_to_skip, i:idx_off_event];
                else
                    if switch_on + avg_runtime <= num_measurements
                        consumption(switch_on:switch_on+avg_runtime) = average_power;
                    else
                        consumption(switch_on:end) = average_power;
                    end
                end
            % off event: check if on_event is within runtime
            % yes: on_event:off_event;
            % [ not possible, this case is handled above ]
            % no: off_event-runtime:off_event [ check boundaries ]
            else
                switch_off = times(i);
                if switch_off - avg_runtime >= 1
                    consumption(switch_off-avg_runtime:switch_off) = average_power;
                else
                    consumption(1:switch_off) = average_power;
                end
            end
        end
        warning('on');

    elseif strcmp(appliance, 'Dishwasher') == 1
        for times = result.usage_times_start
            average_power = 1000; %1 kWh
            consumption(times:times+usage_duration-1) = average_power / usage_duration * 3600;
        end

    elseif strcmp(appliance, 'Water kettle') == 1
        average_power = 1800;
        times = result.events(:,1)
        on_events = result.events(:,3) > 0
        times_on_events = times(on_events);
        times_off_events = times(~on_events);

        for i = 1:length(times_on_events)
           switch_on = times_on_events(i);
           potential_off_events = times_off_events > switch_on & times_off_events < switch_on + usage_duration;
           if any(potential_off_events)
               switch_off = times_off_events(max(find(potential_off_events)));
               consumption(switch_on:switch_off) = average_power;
           end
        end

    else
        error('Determining consumption for appliance %s not implemented yet');
    end
end

author: beckel

repo: nilm-eval

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


text:
% VOLTYPE determines the type of volume conduction model
%
% Use as
%   [type] = voltype(vol)
% to get a string describing the type, or
%   [flag] = voltype(vol, desired)
% to get a boolean value.
%
% See also COMPUTE_LEADFIELD

% Copyright (C) 2007, Robert Oostenveld
%
function [type] = voltype(vol, desired)

if isfield(vol, 'type')
  type = vol.type;
elseif isfield(vol, 'r') && prod(size(vol.r))==1
  type = 'singlesphere';
elseif isfield(vol, 'r') && isfield(vol, 'o') && all(size(vol.r)==size(vol.o))
  type = 'multisphere';
elseif isfield(vol, 'r')
  type = 'concentric';
elseif isfield(vol, 'bnd')
  type = 'bem';
end

if nargin>1
  type = strcmp(type, desired);
end

author: PatternRecognition

repo: OpenBMI


text:
function mbasis = basis_matrix_overhauser_nul ( alpha )

%*****************************************************************************80
%
%% BASIS_MATRIX_OVERHAUSER_NUL sets up the nonuniform left Overhauser spline basis matrix.
%
%  Discussion:
%
%    This basis matrix assumes that the data points P1, P2, and
%    P3 are not uniformly spaced in T, and that P1 corresponds to T = 0,
%    and P2 to T = 1. (???)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ALPHA.
%    ALPHA = || P2 - P1 || / ( || P3 - P2 || + || P2 - P1 || )
%
%    Output, real MBASIS(3,3), the basis matrix.
%
  mbasis(1,1) =   1.0 / alpha;
  mbasis(1,2) = - 1.0 / ( alpha * ( 1.0 - alpha ) );
  mbasis(1,3) =   1.0 / ( 1.0 - alpha );

  mbasis(2,1) = - ( 1.0 + alpha ) / alpha;
  mbasis(2,2) =   1.0 / ( alpha * ( 1.0 - alpha ) );
  mbasis(2,3) = - alpha / ( 1.0 - alpha );

  mbasis(3,1) =   1.0;
  mbasis(3,2) =   0.0;
  mbasis(3,3) =   0.0;

  return
end


text:
function [C,logicConverged]=triSurfPermuteColor(varargin)

%% Parse input

switch nargin
    case 2
        F=varargin{1};
        V=varargin{2};
        nc=4;
        maxIter=1000;
    case 3
        F=varargin{1};
        V=varargin{2};
        nc=varargin{3};
        maxIter=1000;
    case 4
        F=varargin{1};
        V=varargin{2};
        nc=varargin{3};
        maxIter=varargin{4};
    otherwise
        error('Wrong number of input arguments');
end

C=randi(nc,size(F,1),1);

%%

TR = triangulation(F,V);
N = neighbors(TR);

%%
numFaces=size(F,1);

%% Random iterations

q=0;
logicConverged=0;
numPrev=numFaces;
while q<maxIter && logicConverged==0

    q=q+1;

    CN=C(N); %Colors for neighbours

    logicDoubleColor=CN==C(:,ones(size(CN,2),1));
    logicRep=any(logicDoubleColor,2);

    C(logicRep)=randi(nc,nnz(logicRep),1);

    %     set(hp,'CData',C);
    %     drawnow;
    numNow=nnz(logicRep);
    p=abs(numPrev-numNow)./numFaces;

    numPrev=nnz(logicRep);

    logicConverged=nnz(logicRep)==0;

    if q==maxIter
        break
    end
end

logicConverged=nnz(logicRep)==0;

%% Structured iterations

if ~logicConverged
    while logicConverged==0

        CN=C(N); %Colors for neighbours

        logicDoubleColor=CN==C(:,ones(size(CN,2),1));
        logicRep=any(logicDoubleColor,2);

        Ir=find(logicRep,1);

        c=C(Ir);
        cn=CN(Ir,:);

        cFound=0;
        s=1;
        while cFound
            if ~any(cn-s==0);
                C(Ir)=s;
                cFound=1;
            elseif s==nc
                cFound=1;
            end
            s=s+1;
        end
        %         set(hp,'CData',C);
        %         drawnow;

        if nnz(logicRep)==0
            logicConverged=1;
        else
            q=q+1;
        end

    end
end
logicConverged=nnz(logicRep)==0;

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%% Vehicle Articulated Linear
% Linear articulated bicycle model with 4 degrees of freedom.
%
%
%
%
%
%% Sintax
% |dx = _VehicleModel_.Model(t,states,tspan)|
%
% |dx = _VehicleModel_.MassMatrix(t,states,tspan)|
%
%% Arguments
% The following table describes the input arguments:
%
%

%
 <tt>t</tt>	 Time

%
 <tt>states</tt>	 Model state variables: [XT YT PSI PHI VT ALPHAT dPSI dPHI]

%
 <tt>tspan</tt>	 Time span

%

%
%% See Also
%
% <https://www.mathworks.com/matlabcentral/fileexchange/58683-vehicle-dynamics-lateral Home>
%

author: andresmendes

repo: Vehicle-Dynamics-Lateral


text:
function degree = grNodeOuterDegree(node, edges)
%GRNODEOUTERDEGREE Outer degree of a node in a graph
%
%   DEG = grNodeOuterDegree(NODE, EDGES);
%   Returns the outer degree of a node in the given edge list, i.e. the
%   number of edges emanating from it.
%   NODE is the index of the node, and EDGES is a liste of couples of
%   indices (origin and destination node).
%
%   Note: Also works when node is a vector of indices
%
%   See Also:
%   grNodeDegree, grNodeInnerDegree
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2006-01-17
% Copyright 2006 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas).
%

%   HISTORY
%   2008-08-07 pre-allocate memory, update doc

% allocate memory
N = size(node, 1);
degree = zeros(N, 1);

% compute outer degree of each vertex
for i=1:N
    degree(i) = sum(edges(:,1)==node(i));
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function y  = laprnd(m, n, mu, sigma)
%LAPRND generate i.i.d. laplacian random number drawn from laplacian distribution
%   with mean mu and standard deviation sigma.
%   mu      : mean
%   sigma   : standard deviation
%   [m, n]  : the dimension of y.
%   Default mu = 0, sigma = 1.
%   For more information, refer to
%   http://en.wikipedia.org./wiki/Laplace_distribution

%   Author  : Elvis Chen (bee33@sjtu.edu.cn)
%   Date    : 01/19/07

%Check inputs
if nargin < 2
    error('At least two inputs are required');
end

if nargin == 2
    mu = 0; sigma = 1;
end

if nargin == 3
    sigma = 1;
end

% Generate Laplacian noise
u = rand(m, n)-0.5;
b = sigma / sqrt(2);
% y = mu - b.* sign(u).* log(1- 2.* abs(u));
y = mu + sign(u) .* ( pi - b.* log( exp(pi./b) + (2-2.*exp(pi./b)) .* abs(u) ) );

author: yuxianghao


text:
function yml=mbmling(xml,plt)

%MBMLING_CELL converts any XML into MbML (Matlab Markup Language) compliant string
%
%Syntax: yml=mbmling(xml,plt)
%
%Description:
%   Converts any XML syntax into Matlab Markup Language (MbML)
%
% Jonas Almeida, almeidaj@musc.edu, 18 May 2002, MAT4NAT Tbox

if nargin<2;plt=0;end
if plt==1;disp('MbMLing progress report');disp('---------BEGIN----(9 steps)----------');end

% Remove non-content lines if they exist
if plt==1;disp('1. Removing non-content lines if they exist');end
xml=regexprep(xml,'<[?!].*?>','');
% Remove empty spaces between tags
if plt==1;disp('2. Removing empty spaces between tags');end
%xml=regexprep(xml,'>[ ]+?<','><');
xml=strrep(xml,char([13,10]),' '); %replace chariage returns by spaces
n=length(xml)+1;
while n>length(xml)
    n=length(xml);
    xml=strrep(xml,'> ','>');
end
% Replace symbols that may conflict with matlab variable naming by underscore characters
if plt==1;disp('3. Replacing symbols that may conflict with matlab variable naming by underscore characters');end
xml=regexprep(xml,'<[\/]{0,1}(\w*)[^\w\/> ]+(.*?)([ >])','<$1_$2$3','tokenize');
% Replace one-tag format by open / close tagging
if plt==1;disp('4. Replace one-tag format by open / close tagging');end
xml_=regexprep(xml,'<(\w+)([^>]*?)/>','<$1$2></$1>','tokenize');
% Turn attributes into contents
if plt==1;disp('5. Turn attributes into contents');end
xml=regexprep(xml_,'<([^>]+) +(\w+)="(.*?)" *>','<$1><$2>$3</$2>','tokenize');
while ~strcmp(xml,xml_)
    %disp('...')
    xml_=xml;xml=regexprep(xml_,'<([^>]+) +(\w+)="(.*?)" *>','<$1><$2>$3</$2>','tokenize');
end
% remove leftover spaces in tag names
if plt==1;disp('6. Removing leftover spaces in tag names');end
xml=regexprep(xml,'(<\w+) +>','$1>','tokenize');
% Tag untagged contents
if plt==1;disp('7. Tagging untagged contents');end
yml=regexprep(xml,'(</\w+>)([^<>]+)</(\w+)','$1<$3>$2</$3></$3','tokenize');
% Tag cell structures
if plt==1;disp('8. Tag cell structures');end
yml=regexprep(yml,'<(\w+)>','<$1 class="cell"><cell>','tokenize');
yml=regexprep(yml,'<(/\w+)>','</""""><$1>','tokenize');
yml=regexprep(yml,'<(/\w+)>(<\w+ )','<$1></""""><cell>$2','tokenize');
% Remove cell tag arround content values
if plt==1;disp('9. Remove cell tag arround content values');end
yml=regexprep(yml,'class="cell"><cell>([^<]*)</"""">','class="char">$1','tokenize');
yml=strrep(yml,'</"""">','</cell>');
if plt==1;disp('------------END----------------------');end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%MEANSHIFT  Finds an object on a back projection image
%
%     window = cv.meanShift(probImage, window)
%     [window,iter] = cv.meanShift(probImage, window)
%     [...] = cv.meanShift(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __probImage__ Back projection of the object histogram. See
%   cv.calcBackProject for details.
% * __window__ Initial search window `[x,y,w,h]`.
%
% ## Output
% * __window__ Converged CAMSHIFT window `[x,y,w,h]`.
% * __iter__ Number of iterations CAMSHIFT took to converge.
%
% ## Options
% * __Criteria__ Stop criteria for the iterative search algorithm. Accepts a
%   struct with 'type', 'maxCount', and 'epsilon' fields. Default
%   `struct('type','Count+EPS', 'maxCount',100, 'epsilon',1.0)`
%
% The function implements the iterative object search algorithm. It takes the
% input back projection of an object and the initial position. The mass center
% in window of the back projection image is computed and the search window
% center shifts to the mass center. The procedure is repeated until the
% specified number of iterations `Criteria.maxCount` is done or until the
% window center shifts by less than `Criteria.epsilon`. The algorithm is used
% inside cv.CamShift and, unlike cv.CamShift, the search window size or
% orientation do not change during the search. You can simply pass the output
% of cv.calcBackProject to this function. But better results can be obtained
% if you pre-filter the back projection and remove the noise. For example, you
% can do this by retrieving connected components with cv.findContours,
% throwing away contours with small area (cv.contourArea), and rendering the
% remaining contours with cv.drawContours.
%
% See also: cv.CamShift, cv.calcBackProject
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function m_demo(index)
% M_DEMO  Demonstration program showing various maps in M_Map package
%         Dig into this to look for examples of things you want to do.
%
%         M_DEMO runs all the demos.
%         M_DEMO(NUM) runs examples NUM (1<=NUM<=10), pausing between examples.
%
%         Some demos may require you to install GSHHS or TerrainBase datafiles
%         (see documentation)

% Rich Pawlowicz (rich@ocgy.ubc.ca) 7/May/1997
% (thanks to Art Newhall for putting these examples into an m-file,
% and the Chuck Denham for enhancing the interface)
%
% 27/July/98 - more examples.
% 17/Aug/98     "
% 15/Nov/98  - another example, better interface.
% 23/Dec/98  - another example.

%
% This software is provided "as is" without warranty of any kind. But
% it's mine, so you can't sell it.

global MAP_PROJECTION

N_EXAMPLES=15;

if nargin==0
 index=1:N_EXAMPLES;
end

for i=index

clf;
switch i

  case 1

    m_proj('ortho','lat',48','long',-123');
    m_coast('patch','r');
    m_grid('linestyle','-','xticklabels',[],'yticklabels',[],'ytick',[-80:40:80]);
    xlabel('Orthographic Projection','visible','on');

  case 2

    m_proj('lambert','long',[-160 -40],'lat',[30 80]);
    m_coast('patch',[1 .85 .7]);
    [CS,CH]=m_elev('contourf',[500:500:4000]);
 %   m_elev('pcolor');
    m_grid('box','fancy','tickdir','in');
    colormap(flipud(copper));
    xlabel('Conic Projection of North America with elevations','visible','on');
    m_contfbar([0 .3],.9,CS,CH);

  case 3

    m_proj('stereographic','lat',90,'long',30,'radius',25);
    m_elev('contour',[-3500:1000:-500],'linecolor','b');
    m_grid('xtick',12,'tickdir','out','ytick',[70 80],'linestyle','-');
    m_coast('patch',[.7 .7 .7],'edgecolor','r');
    xlabel('Polar Stereographic Projection with bathymetry','visible','on');

  case 4

    subplot(211);
    Slongs=[-100 0;-75 25;0  45; 25 145;45 100;145 295;100 295];
    Slats= [  8 80;-80  8; 8 80;-80   8; 8  80;-80   0;  0  80];
    for l=1:7
     m_proj('sinusoidal','long',Slongs(l,:),'lat',Slats(l,:));
   %  colormap(m_colmap('blues'));caxis([-6000 0]);
   %  m_elev('shadedrelief');
     m_grid('fontsize',6,'xticklabels',[],'xtick',[-180:30:360],...
            'ytick',[-80:20:80],'yticklabels',[],'linestyle','-','color',[.9 .9 .9]);
     m_coast('patch','g');
    end
    xlabel('Interrupted Sinusoidal Projection of World Oceans');
    % In order to see all the maps we must undo the axis limits set by m_grid calls:
    set(gca,'xlimmode','auto','ylimmode','auto');

    subplot(212);
    Slongs=[-100 43;-75 20; 20 145;43 100;145 295;100 295];
    Slats= [  0  90;-90  0;-90   0; 0  90;-90   0;  0  90];
    for l=1:6
     m_proj('mollweide','long',Slongs(l,:),'lat',Slats(l,:));
   %  colormap(m_colmap('blues'));caxis([-6000 0]);
   %  m_elev('shadedrelief');
     m_grid('fontsize',6,'xticklabels',[],'xtick',[-180:30:360],...
            'ytick',[-80:20:80],'yticklabels',[],'linestyle','-','color','k');

     m_coast('patch',[.6 .6 .6]);
    end
    xlabel('Interrupted Mollweide Projection of World Oceans');
    set(gca,'xlimmode','auto','ylimmode','auto');

  case 5

    %% Nice looking data
    [lon,lat]=meshgrid([-136:2:-114],[36:2:54]);
    u=sin(lat/6);
    v=sin(lon/6);

    m_proj('oblique','lat',[56 30],'lon',[-132 -120],'aspect',.8);

    subplot(121);
    m_coast('patch',[.9 .9 .9],'edgecolor','none');
    m_grid('tickdir','out','yaxislocation','right',...
	   'xaxislocation','top','xlabeldir','end','ticklength',.02);
    hold on;
    m_quiver(lon,lat,u,v);
    xlabel('Simulated surface winds');

    subplot(122);
    m_coast('patch',[.9 .9 .9],'edgecolor','none');
    m_grid('tickdir','out','yticklabels',[],...
	   'xticklabels',[],'linestyle','none','ticklength',.02);
    hold on;
    [cs,h]=m_contour(lon,lat,sqrt(u.*u+v.*v));
    clabel(cs,h,'fontsize',8);
    xlabel('Simulated something else');

  case 6

    % Plot a circular orbit
    lon=[-180:180];
    lat=atan(tan(60*pi/180)*cos((lon-30)*pi/180))*180/pi;

    m_proj('miller','lat',82);
    m_coast('color',[0 .6 0]);
    m_line(lon,lat,'linewidth',3,'color','r');
    m_grid('linestyle','none','box','fancy','tickdir','out');

  case 7

    m_proj('lambert','lon',[-10 20],'lat',[33 48]);
    if MAP_PROJECTION.IsOctave
       [CS,CH]=m_etopo2('contourf',[-5000:500:0 250:250:3000],'linecolor','none');
    else
       [CS,CH]=m_etopo2('contourf',[-5000:500:0 250:250:3000],'edgecolor','none');
    end
    m_grid('linestyle','none','tickdir','out','linewidth',3);

    colormap([ m_colmap('blues',80); m_colmap('gland',48)]);
    brighten(.5);

    ax=m_contfbar(1,[.5 .8],CS,CH);
    title(ax,{'Level/m',''}); % Move up by inserting a blank line

  case 8

    m_vec;
    colormap(jet);

  case 9

    % Example showing the default coastline and all of the GSHHS coastlines.

    axes('position',[.35 .6 .37 .37]);
    m_proj('albers equal-area','lat',[40 60],'long',[-90 -50],'rect','on');
    m_coast('patch',[0 1 0]);
    m_grid('linestyle','none','linewidth',2,'tickdir','out','xaxisloc','top','yaxisloc','right','fontsize',6);
    m_text(-69,51,'Standard coastline','color','r','fontweight','bold');
    m_ruler([.5 .9],.8,3,'fontsize',8);
    drawnow;

    axes('position',[.09 .5 .37 .37]);
    m_proj('albers equal-area','lat',[40 54],'long',[-80 -55],'rect','on');
    m_gshhs_c('patch',[.2 .8 .2]);
    m_grid('linestyle','none','linewidth',2,'tickdir','out','xaxisloc','top','fontsize',6);
    m_text(-80,52.5,'GSHHS\_C (crude)','color','m','fontweight','bold');
    m_ruler([.5 .9],.8,2,'fontsize',8);
    drawnow;

    axes('position',[.13 .2 .37 .37]);
    m_proj('albers equal-area','lat',[43 48],'long',[-67 -58],'rect','on');
    m_gshhs_l('patch',[.4 .6 .4]);
    m_grid('linestyle','none','linewidth',2,'tickdir','out','fontsize',6);
    m_text(-66.5,43.5,'GSHHS\_L (low)','color','m','fontweight','bold');
    m_ruler([.5 .9],.8,3,'fontsize',8);
    drawnow;

    axes('position',[.35 .05 .37 .37]);
    m_proj('albers equal-area','lat',[45.8 47.2],'long',[-64.5 -62],'rect','on');
    m_gshhs_i('patch',[.5 .6 .5]);
    m_grid('linestyle','none','linewidth',2,'tickdir','out','yaxisloc','right','fontsize',6);
    m_text(-64.4,45.9,'GSHHS\_I (intermediate)   ','color','m','fontweight','bold','horizontalalignment','right');
    m_ruler([.5 .8],.1,3,'fontsize',8);
    drawnow;

    axes('position',[.5 .1 .37 .37]);
    m_proj('albers equal-area','lat',[46.375 46.6],'long',[-64.2 -63.7],'rect','on');
    m_gshhs_h('patch',[.6 .7 .6]);
    m_grid('linestyle','none','linewidth',2,'tickdir','out','xaxisloc','top','yaxisloc','right','fontsize',6);
    m_text(-64.18,46.4,'GSHHS\_H (high)','color','m','fontweight','bold');
    m_ruler([.5 .8],.2,3,'fontsize',8);
    drawnow;

    axes('position',[.55 .35 .37 .37]);
    m_proj('albers equal-area','lat',[46.55 46.65],'long',[-63.97 -63.77],'rect','on');
    m_gshhs_f('patch',[.7 .9 .7]);
    m_grid('linestyle','none','linewidth',2,'tickdir','out','xaxisloc','top','yaxisloc','right','fontsize',6);
    m_text(-63.95,46.56,'GSHHS\_F (full)','color','m','fontweight','bold');
    m_ruler([.5 .8],.2,3,'fontsize',8);

  case 10

    % Example showing a trackline plot

    clf
    m_proj('UTM','long',[-72 -68],'lat',[40 44]);
    m_gshhs_i('color','k');
    m_grid('box','fancy','tickdir','in');
    m_ruler(1.2,[.5 .8]);

    % fake up a trackline
    lons=[-71:.1:-67];
    lats=60*cos((lons+115)*pi/180);
    dates=datenum(1997,10,23,15,1:41,zeros(1,41));

    m_track(lons,lats,dates,'ticks',0,'times',4,'dates',8,...
           'clip','off','color','r','orient','upright');

  case 11

    % example showing range rings

    clf
    m_proj('hammer','clong',170);
    m_grid('xtick',[],'ytick',[],'linestyle','-');
    m_coast('patch','g');
    m_line(100.5,13.5,'marker','s','color','r');
    m_range_ring(100.5,13.5,[1000:1000:15000],'color','b','linewidth',2);
    xlabel('1000km range rings from Bangkok');

  case 12

    % Example showing speckle

    bndry_lon=[-128.8 -128.8 -128.3 -128 -126.8 -126.6 -128.8];
    bndry_lat=[49      50.33  50.33  50   49.5   49     49];

    clf;
    m_proj('lambert','long',[-130 -121.5],'lat',[47 51]);
    m_gshhs_i('color','k');
    m_gshhs_i('speckle','color','k');
    m_line(bndry_lon,bndry_lat,'linewidth',2,'color','k');     % Area outline ...
    m_hatch(bndry_lon,bndry_lat,'single',30,5,'color','k'); % ...with hatching added.

    m_grid('linewidth',2,'linestyle','none');
    title({'Speckled Boundaries','for nice B&W presentation','(best in postscript format)'});
    m_text(-128,48,{'Pacific','Ocean'},'fontsize',18);

  case 13

    % Colouring the ocean blue

    clf
    m_proj('miller','lat',[-77 77]);
 %   set(gca,'color',[.9 .99 1]);
    m_coast('patch',[.7 1 .7],'edgecolor','none');
    m_grid('box','fancy','linestyle','-','gridcolor','w','backcolor',[.2 .65 1]);

    cities={'Cairo','Washington','Buenos Aires'};
    lons=[ 30+2/60  -77-2/60   -58-22/60];
    lats=[ 31+21/60  38+53/60  -34-45/60];

    for k=1:3
      [range,ln,lt]=m_lldist([-123-6/60 lons(k)],[49+13/60  lats(k)],40);
      m_line(ln,lt,'color','r','linewidth',2);
      m_text(ln(end),lt(end),sprintf('%s - %d km',cities{k},round(range)));
    end

    % set(gcf,'color','w');  % To defeat the tendency of print to turn white into black

    title('Great Circle Routes','fontsize',12,'fontweight','bold');

    case 14
        clf
        m_proj('lambert','long',[-130 -122],'lat',[48 52.5],'rect','on');
        if MAP_PROJECTION.IsOctave
           [CS,CH]=m_elev('contourf',[-3000:500:-500 -200 -100 -50 -1 2 20 50 100 250:250:2000],'linecolor','none');
        else
           [CS,CH]=m_etopo2('contourf',[-3000:500:-500 -200 -100 -50 -1 2 20 50 100 250:250:2000],'edgecolor','none');
        end

        m_grid('linewi',2,'tickdir','out','yaxisloc','right');

        if MAP_PROJECTION.IsOctave
            ax=m_contfbar(-.03,[.5 .8],CS,CH,'linecolor','none');
        else
            ax=m_contfbar(-.03,[.5 .8],CS,CH,'edgecolor','none');
        end
        title(ax,{'meters',''}); % Move up by inserting a blank line

        colormap([m_colmap('blues',96);m_colmap('gland',64)]);
        caxis([-3000 2000]);

    case 15
        clf;
        m_proj('azimuthal equal-area','radius',156,'lat',-46,'long',-95,'rot',30);

        ax1=subplot(2,2,1,'align');
        m_coast('patch','r');
        m_grid('xticklabel',[],'yticklabel',[],'linestyle','-','ytick',[-60:30:60]);

        ax2=subplot(2,2,2,'align');
        if MAP_PROJECTION.IsOctave
            m_elev('contourf',[-7000:1000:0 500:500:3000],'linecolor','none');
        else
             m_elev('contourf',[-7000:1000:0 500:500:3000],'edgecolor','none');
        end
        colormap(ax2,[m_colmap('blues',70);m_colmap('gland',30)]);
        caxis(ax2,[-7000 3000]);
        m_grid('xticklabel',[],'yticklabel',[],'linestyle','-','ytick',[-60:30:60]);

        ax3=subplot(2,2,3,'align');
        colormap(ax3,[m_colmap('blues',70);m_colmap('gland',30)]);
        caxis(ax3,[-7000 3000]);
        m_elev('image');
        m_grid('xticklabel',[],'yticklabel',[],'linestyle','-','ytick',[-60:30:60]);

        ax4=subplot(2,2,4,'align');
        colormap(ax4,[m_colmap('blues')]);
        caxis(ax4,[-8000 000]);
        m_elev('shadedrelief','gradient',.5);
        m_coast('patch',[.7 .7 .7],'edgecolor','none');
        m_grid('xticklabel',[],'yticklabel',[],'linestyle','-','ytick',[-60:30:60]);

        ha = axes('Position',[0 0 1 1],'Xlim',[0 1],'Ylim',[0  1],'Box','off','Visible','off','Units','normalized', 'clipping' , 'off');
        text(0.5, 0.98,'This projection shows all oceans connected to each other','horiz','center','fontsize',20);

end

 if i<length(index)
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   disp('  hit return to continue');
   pause
   disp('        ...drawing');
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 end
end

author: goGPS-Project


text:
function out=sparseconv(x,f)

% SPARSECONV Speedy convolution of sparse vectors
%
%    OUT = SPARSECONV(X, F) computes the convolution of the
%    (sparse) input vectors X and F.
%
%    Gert Cuypers 4/10/2001
%    Esat Sista KULeuven

[rij,kol]=size(x);
if kol>1
   x=x.';
end
[rij,kol]=size(x);
if kol>1
   error('x moet vector zijn')
end

[rij,kol]=size(f);
if kol>1
   f=f.';
end

[rij,kol]=size(f);
if kol>1
   error('f moet vector zijn');
end
xs=sparse(x);
[xrow, xcol, xval]=find(xs);
nietnulx=length(xrow);
fs=sparse(f);
[frow, fcol, fval]=find(fs);
nietnulf=length(frow);
xpos=zeros(nietnulf*nietnulx,1);
ypos=zeros(nietnulf*nietnulx,1);
elementen=zeros(nietnulf*nietnulx,1);
dummy=(frow)*ones(1,nietnulx);
ypos(:)=dummy;
xpos(:)=dummy+ones(nietnulf,1)*xrow.'-1;
elementen(:)=ones(nietnulf,1)*xval.';
out=full(sparse(xpos,ypos,elementen,length(x)+length(f)-1,length(f))*fs);

author: Sable

repo: mcbench-benchmarks


text:
function obj = createRaster(obj, varargin)
% CREATERASTER Create raster from point cloud.

% Input parsing ----------------------------------------------------------------

validRasterValue = {'average' 'linear' 'mask' 'plane'};

p = inputParser;
p.addParameter('CellSize'     , 1        , @(x) isnumeric(x) && x>0);
p.addParameter('RasterValue'  , 'average', @(x) any(strcmpi(x, validRasterValue)));
p.addParameter('SearchRadius' , 3        , @(x) isnumeric(x) && x>0); % for 'RasterValue' 'average' and 'plane'
p.parse(varargin{:});
p = p.Results;

% Start ------------------------------------------------------------------------

procHierarchy = {'POINTCLOUD' 'RASTER'};
msg('S', procHierarchy);
msg('I', procHierarchy, sprintf('Point cloud label = ''%s''', obj.label));
msg('I', procHierarchy, sprintf('IN: RasterValue = ''%s''', p.RasterValue));
msg('V', p.CellSize, 'IN: CellSize', 'Prec', 2);
if any(strcmpi(p.RasterValue, {'average' 'plane'}))
    msg('V', p.SearchRadius, 'IN: SearchRadius', 'Prec', 2);
end

% Create query points ----------------------------------------------------------

% Origin (lower left point) of points
lim.min = min(obj.X(obj.act,[1 2]));
% Round origin (raster maps have coincident cell centers if mod(100, p.CellSize) == 0)
lim.min = (floor(lim.min/100))*100;

% Find raster cells in which points lie (starts with 0/0 at lower left corner (o))
% +-----+-----+-----+
% + 1/0 | 1/1 | 1/2 |
% +-----+-----+-----+
% + 0/0 | 0/1 | 0/2 |
% o-----+-----+-----+
rowCol = [floor((obj.X(obj.act,1)-lim.min(1))/p.CellSize) ...
          floor((obj.X(obj.act,2)-lim.min(2))/p.CellSize)];

% Remove multiple points
rowCol = unique(rowCol, 'rows');

% Transformation of indices to coordinate system ('query points')
xq = lim.min(1)+p.CellSize/2+rowCol(:,1)*p.CellSize;
yq = lim.min(2)+p.CellSize/2+rowCol(:,2)*p.CellSize;

% Assign RasterValue -----------------------------------------------------------

% Moving average
if strcmpi(p.RasterValue, 'average')
    idxAct = find(obj.act);
    ns = createns(obj.X(idxAct,[1 2]));
    idxSearchRadius = ns.rangesearch([xq yq], p.SearchRadius);
    zq = cellfun(@(x) mean(obj.X(idxAct(x),3)), idxSearchRadius);
end

% Linear interpolation
if strcmpi(p.RasterValue, 'linear')
    F = scatteredInterpolant(obj.X(obj.act,1), obj.X(obj.act,2), obj.X(obj.act,3), 'linear', 'none');
    zq = F(xq, yq);
end

% Mask
if strcmpi(p.RasterValue, 'Mask')
    zq = 1;
end

% Moving planes
if strcmpi(p.RasterValue, 'plane')
    idxAct = find(obj.act);
    msg('S', {procHierarchy{:} 'BUILD 2D KD-TREE'});
    ns = createns(obj.X(idxAct,[1 2]));
    msg('E', {procHierarchy{:} 'BUILD 2D KD-TREE'});
    msg('S', {procHierarchy{:} 'NNSEARCH'});
    idxSearchRadius = ns.rangesearch([xq yq], p.SearchRadius);
    msg('E', {procHierarchy{:} 'NNSEARCH'});

    msg('S', {procHierarchy{:} 'PLANE INTERPOLATION'});
    for c = 1:numel(xq)

        x = obj.X(idxAct(idxSearchRadius{c}),1);
        y = obj.X(idxAct(idxSearchRadius{c}),2);
        z = obj.X(idxAct(idxSearchRadius{c}),3);

        xm = mean(x);
        ym = mean(y);
        zm = mean(z);

        C = cov([x-xm, y-ym, z-zm]);
        [V, ~] = eig(C);
        n = V(:,1);

        if numel(n) ~= 3
            zq(c) = NaN;
        else
            zq(c) = zm - (n(1)*(xq(c)-xm) + n(2)*(yq(c)-ym)) / n(3);
        end

    end
    msg('E', {procHierarchy{:} 'PLANE INTERPOLATION'});

end

% Create GeoTiff ---------------------------------------------------------------

% Size of data array
sizeA = max(rowCol)-min(rowCol)+1;

% World matrix
W = [p.CellSize           0 min(xq)
              0 -p.CellSize max(yq)];

% Create spatialref.MapRasterReference object for GeoTiff
R = maprasterref(W, [sizeA(2) sizeA(1)]); % cols -> x, rows -> y

% Initialize data array with noData values
A = NaN(sizeA);

% Linear indices
idxLin = sub2ind(size(A), rowCol(:,1)-min(rowCol(:,1))+1, ...
                          rowCol(:,2)-min(rowCol(:,2))+1);

% Set data values
A(idxLin) = zq;

% Transform A, so that element 1,1 is upper, left corner (min(xq), max(yq))
A = flipud(A');

% Raster reduction point
if strcmpi(p.RasterValue, 'mask')
    RedPoi = [obj.redPoi(1) obj.redPoi(2) 0]; % no height reduction
else
    RedPoi = obj.redPoi;
end

% Add grid object to pointCloud object
obj.raster = raster({A R}, ...
                    obj.label, ...
                    'RedPoi', RedPoi);

% End --------------------------------------------------------------------------

msg('E', procHierarchy);

% 4Debug -----------------------------------------------------------------------

% save2('n:\Dropbox\Matlab\Scripts\2013_05\testgrid.xyz', [xq yq zeros(numel(xq),1)], '%.3f %.3f %.3f\n');

end

author: pglira


text:
function W = FormingConvFilter(Datacube,ConvFilter,num_of_filters)
% For each pixel in the Datacube, extract a small
% cube: ConvFilter(Height*Width*Channel), vectorizing into a column vector
% For all pixels in the Datacube, performing PCA/WPCA to gather
% num_of_filters eigenvecs, then transforming each eigenvecs which corresponds to
% a column of W into a small cube as our 3D ConvFilter.

X = im2colstep(Datacube,[ConvFilter.PatchSize,ConvFilter.PatchSize,ConvFilter.Channel]); % Dominated by column to extract
mu = mean(X,2);
% mu = mean(X,1);
X = bsxfun(@minus, X, mu);
Rx = X*X'/size(X,2);
[E,D] = eig(Rx);
[~, ind] = sort(diag(D),'descend');
% W = E(:,ind(1:num_of_filters));  % principal eigenvectors
W = E(:,ind(1:3));
% WPCA
% D = diag(D);
% D = D(ind(1:num_of_filters));
% W = bsxfun(@rdivide,W,sqrt(D)');

% W = cell(num_of_filters,1);
% for i=1:size(V,2)
%     W{i} = reshape(V(:,i),[Height,Width,Channel]);
% end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [ x, seed ] = circular_normal_sample ( a, b, seed )

%*****************************************************************************80
%
%% CIRCULAR_NORMAL_SAMPLE samples the Circular Normal PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A(2), a parameter of the PDF, the mean value.
%
%    Input, real B, a parameter of the PDF, the standard deviation.
%
%    Input/output, integer SEED, a seed for the random
%    number generator.
%
%    Output, real X(2), a sample of the PDF.
%
  [ v1, seed ] = r8_uniform_01 ( seed );
  [ v2, seed ] = r8_uniform_01 ( seed );

  r = sqrt ( - 2.0 * log ( v1 ) );

  x(1) = a(1) + b * r * cos ( 2.0 * pi * v2 );
  x(2) = a(2) + b * r * sin ( 2.0 * pi * v2 );

  return
end


text:
function model = gpsimMapUpdateG(model)

% GPSIMMAPUPDATEG Update the nonlinear transformation of f.
% FORMAT
% DESC updates the fields of model associated with the non-linear
% transformation of f, namely g, g_grad and g_grad2.
% ARG model : the model to be updated.
% ARG model : the model with the updated g representation.
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% MODIFIED : Pei Gao, 2008
%
% SEEALSO : gpsimMapFunctionalExpandParam, gpsimMapCreate

% SHEFFIELDML

% Remove mean function value from m (if mean function present).
if isfield(model, 'meanFunction') & ~isempty(model.meanFunction)
  f = model.f + modelOut(model.meanFunction, model.mapt);
else
  f = model.f;
end

% make sure f is a column vector
f = reshape(f, length(f), 1);

if model.ngParam == 0

switch model.nonLinearity
   case 'linear'
    model.g=f;        %linear case
    model.g_grad=ones(size(model.g));
    model.g_grad2=zeros(size(model.g));
    model.g_grad3=zeros(size(model.g));
    model.isConcave = true; % Is the log likelihood concave?
   case 'exp'
    model.g=exp(f);   %positive TF concentrations
%   model.gParam
%   if ~isempty(model.gParam)
%   model.ngParam = length(model.gParam(1,:));
%   model.dg = [];
%   model.dg2 = [];
%   end
    model.g_grad=model.g;
    model.g_grad2=model.g;
    model.g_grad3=model.g;
    model.isConcave = true;
   case 'quadratic'
    model.g=f.*f;
    model.g_grad=2*f;
    model.g_grad2=2*ones(size(f));
    model.isConcave = false;
   case 'negLogLogit'
    model.g=log(1+exp(f));   %positive TF concentrations
    model.g_grad=sigmoid(f);
    model.g_grad2=model.g_grad.*(1-model.g_grad);
    model.isConcave = true;
   case 'sigmoid'
    model.g = sigmoid(f);
    model.g_grad = model.g.*(1 - model.g);
    model.g_grad2 = model.g.*(1-model.g) - 2*model.g.*(model.g.*(1-model.g));
    model.isConcave = false;
   case 'repression'
    model.g = gpsimModelFunctions(model,'g');
    model.g_grad = gpsimModelFunctions(model,'grad');
    model.g_grad2 = gpsimModelFunctions(model,'grad2');
    model.g_grad3 = gpsimModelFunctions(model,'grad3');;
    model.isConcave = true; %% ??
   case 'activation'
    model.g = gpsimModelFunctions(model,'g');
    model.g_grad = gpsimModelFunctions(model,'grad');
    model.g_grad2 = gpsimModelFunctions(model,'grad2');
    model.g_grad3 = gpsimModelFunctions(model,'grad3');
    model.isConcave = true; %% ??
    if ~isempty(model.gParam)
      model.ngParam = length(model.gParam(1,:));
      model.dg = gpsimModelFunctions(model,'paramGrad');
      model.dg2 = gpsimModelFunctions(model,'paramGrad2');
    end
   otherwise
    error('Invalid non-linearity.')
end

else
  model.g = gpsimModelFunctions(model,'g');
  model.g_grad = gpsimModelFunctions(model,'grad');
  model.g_grad2 = gpsimModelFunctions(model,'grad2');
  model.g_grad3 = gpsimModelFunctions(model,'grad3');
  model.isConcave = true;
  model.dg = gpsimModelFunctions(model,'paramGrad');
  model.dggrad = gpsimModelFunctions(model,'paramGgrad');
  model.dggrad2 = gpsimModelFunctions(model,'paramGgrad2');
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function p = prior_sqrtt(varargin)
%PRIOR_SQRTT  Student-t prior structure for the square root of the parameter
%
%  Description
%    P = PRIOR_SQRTT('PARAM1', VALUE1, 'PARAM2', VALUE2, ...)
%    creates for the square root of the parameter Student's
%    t-distribution prior structure in which the named parameters
%    have the specified values. Any unspecified parameters are set
%    to default values.
%
%    P = PRIOR_SQRTT(P, 'PARAM1', VALUE1, 'PARAM2', VALUE2, ...)
%    modify a prior structure with the named parameters altered
%    with the specified values.
%
%    Parameters for Student-t prior [default]
%      mu       - location [0]
%      s2       - scale [1]
%      nu       - degrees of freedom [4]
%      mu_prior - prior for mu [prior_fixed]
%      s2_prior - prior for s2 [prior_fixed]
%      nu_prior - prior for nu [prior_fixed]
%
%  See also
%    PRIOR_t, PRIOR_*

% Copyright (c) 2000-2001,2010 Aki Vehtari
% Copyright (c) 2009 Jarno Vanhatalo
% Copyright (c) 2010 Jaakko Riihim�ki

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  ip=inputParser;
  ip.FunctionName = 'PRIOR_SQRTT';
  ip.addOptional('p', [], @isstruct);
  ip.addParamValue('mu',0, @(x) isscalar(x));
  ip.addParamValue('mu_prior',[], @(x) isstruct(x) || isempty(x));
  ip.addParamValue('s2',1, @(x) isscalar(x) && x>0);
  ip.addParamValue('s2_prior',[], @(x) isstruct(x) || isempty(x));
  ip.addParamValue('nu',4, @(x) isscalar(x) && x>0);
  ip.addParamValue('nu_prior',[], @(x) isstruct(x) || isempty(x));
  ip.parse(varargin{:});
  p=ip.Results.p;

  if isempty(p)
    init=true;
    p.type = 'Sqrt-t';
  else
    if ~isfield(p,'type') && ~isequal(p.type,'Sqrt-t')
      error('First argument does not seem to be a valid prior structure')
    end
    init=false;
  end

  % Initialize parameters
  if init || ~ismember('mu',ip.UsingDefaults)
    p.mu = ip.Results.mu;
  end
  if init || ~ismember('s2',ip.UsingDefaults)
    p.s2 = ip.Results.s2;
  end
  if init || ~ismember('nu',ip.UsingDefaults)
    p.nu = ip.Results.nu;
  end
  % Initialize prior structure
  if init
    p.p=[];
  end
  if init || ~ismember('mu_prior',ip.UsingDefaults)
    p.p.mu=ip.Results.mu_prior;
  end
  if init || ~ismember('s2_prior',ip.UsingDefaults)
    p.p.s2=ip.Results.s2_prior;
  end
  if init || ~ismember('nu_prior',ip.UsingDefaults)
    p.p.nu=ip.Results.nu_prior;
  end

  if init
    % set functions
    p.fh.pak = @prior_sqrtt_pak;
    p.fh.unpak = @prior_sqrtt_unpak;
    p.fh.lp = @prior_sqrtt_lp;
    p.fh.lpg = @prior_sqrtt_lpg;
    p.fh.recappend = @prior_sqrtt_recappend;
  end

end

function [w, s] = prior_sqrtt_pak(p)

  w=[];
  s={};
  if ~isempty(p.p.mu)
    w = p.mu;
    s=[s; 'Sqrt-Student-t.mu'];
  end
  if ~isempty(p.p.s2)
    w = [w log(p.s2)];
    s=[s; 'log(Sqrt-Student-t.s2)'];
  end
  if ~isempty(p.p.nu)
    w = [w log(p.nu)];
    s=[s; 'log(Sqrt-Student-t.nu)'];
  end
end

function [p, w] = prior_sqrtt_unpak(p, w)

  if ~isempty(p.p.mu)
    i1=1;
    p.mu = w(i1);
    w = w(i1+1:end);
  end
  if ~isempty(p.p.s2)
    i1=1;
    p.s2 = exp(w(i1));
    w = w(i1+1:end);
  end
  if ~isempty(p.p.nu)
    i1=1;
    p.nu = exp(w(i1));
    w = w(i1+1:end);
  end
end

function lp = prior_sqrtt_lp(x, p)

  lJ = -log(2*sqrt(x));  % log(1/(2*x^(1/2))) log(|J|) of transformation
  xt = sqrt(x);          % transformation
  lp = sum(gammaln((p.nu+1)./2) - gammaln(p.nu./2) - 0.5*log(p.nu.*pi.*p.s2) - (p.nu+1)./2.*log(1+(xt-p.mu).^2./p.nu./p.s2) + lJ);

  if ~isempty(p.p.mu)
    lp = lp + p.p.mu.fh.lp(p.mu, p.p.mu);
  end
  if ~isempty(p.p.s2)
    lp = lp + p.p.s2.fh.lp(p.s2, p.p.s2) + log(p.s2);
  end
  if ~isempty(p.p.nu)
    lp = lp + p.p.nu.fh.lp(p.nu, p.p.nu) + log(p.nu);
  end
end

function lpg = prior_sqrtt_lpg(x, p)

  lJg = -1./(2*x);        % gradient of log(|J|) of transformation
  xt  = sqrt(x);          % transformation
  xtg = 1./(2*sqrt(x));   % derivative of transformation
  lpg = xtg.*(-(p.nu+1).*(xt-p.mu)./(p.nu.*p.s2 + (xt-p.mu).^2)) + lJg;

  if ~isempty(p.p.mu)
    lpgmu = sum((p.nu+1).*(xt-p.mu)./(p.nu.*p.s2 + (xt-p.mu).^2)) + p.p.mu.fh.lpg(p.mu, p.p.mu);
    lpg = [lpg lpgmu];
  end
  if ~isempty(p.p.s2)
    lpgs2 = (sum(-1./(2.*p.s2)+((p.nu + 1)*(p.mu - xt)^2)./(2*p.s2*((p.mu-xt)^2 + p.nu*p.s2))) + p.p.s2.fh.lpg(p.s2, p.p.s2)).*p.s2 + 1;
    lpg = [lpg lpgs2];
  end
  if ~isempty(p.p.nu)
    lpgnu = (0.5*sum(digamma1((p.nu+1)./2)-digamma1(p.nu./2)-1./p.nu-log(1+(xt-p.mu).^2./p.nu./p.s2)+(p.nu+1)./(1+(xt-p.mu).^2./p.nu./p.s2).*(xt-p.mu).^2./p.s2./p.nu.^2) + p.p.nu.fh.lpg(p.nu, p.p.nu)).*p.nu + 1;
    lpg = [lpg lpgnu];
  end
end

function rec = prior_sqrtt_recappend(rec, ri, p)
% The parameters are not sampled in any case.
  rec = rec;
  if ~isempty(p.p.mu)
    rec.mu(ri,:) = p.mu;
  end
  if ~isempty(p.p.s2)
    rec.s2(ri,:) = p.s2;
  end
  if ~isempty(p.p.nu)
    rec.nu(ri,:) = p.nu;
  end
end


text:
function y = goertzel_classic(x,indvec)
% GOERTZEL_CLASSIC(X,INDVEC) computes DFT of one-dimensional signal X at indices
% contained in INDVEC, using the traditional second-order Goertzel algorithm.
% The indices must be integer values from 0 to N-1, where N is the length of
% X. (Index 0 corresponds to the DC component.)
%
% The output is a column complex vector of length LENGTH(INDVEC) containing
% the desired DFT coefficients.
%
% See also: goertzel_general_shortened.

% (c) 2009-2012, Pavel Rajmic, Brno University of Technology, Czech Rep.

%% Check the input arguments
if nargin < 2
    error('Not enough input arguments')
end
if ~isvector(x) || isempty(x)
    error('X must be a nonempty vector')
end

if ~isvector(indvec) || isempty(indvec)
    error('INDVEC must be a nonempty vector')
end
if ~isreal(indvec)
    error('INDVEC must contain real numbers')
end
% if ~isinteger(indvec)
%     error('INDVEC must contain only integer values')
% end

lx = length(x);
x = reshape(x,lx,1); %forcing x to be column

%% Initialization
no_freq = length(indvec); %number of frequencies to compute
y = zeros(no_freq,1); %memory allocation for the output coefficients

%% Computation via second-order system
% loop over the particular frequencies
for cnt_freq = 1:no_freq

    %for a single frequency:
    %a/ precompute the constants
    pik_term = 2*pi*(indvec(cnt_freq))/(lx);
    cos_pik_term2 = cos(pik_term) * 2;
    cc = exp(-1i*pik_term); % complex constant
    %b/ state variables
    s0 = 0;
    s1 = 0;
    s2 = 0;
    %c/ 'main' loop
    for ind = 1:lx %number of loops is the same as the length of signal
        %new state
        s0 = x(ind) + cos_pik_term2*s1 - s2;  % (*)
        %shifting the state variables
        s2 = s1;
        s1 = s0;
    end
    %d/ final computations
    s0 = cos_pik_term2*s1 - s2; %correspond to one extra performing of (*), where we set x(N+1)=0
    y(cnt_freq) = s0 - s1*cc; %resultant complex DFT coefficient
end


text:
function vol = ft_headmodel_bemcp(geom, varargin)

% FT_HEADMODEL_BEMCP creates a volume conduction model of the head
% using the boundary element method (BEM) for EEG. This function
% takes as input the triangulated surfaces that describe the boundaries
% and returns as output a volume conduction model which can be used
% to compute leadfields.
%
% The implementation of this function is based on Christophe Phillips'
% MATLAB code, hence the name "bemcp".
%
% Use as
%   vol = ft_headmodel_bem_cp(geom, ...)
%
% See also FT_PREPARE_VOL_SENS, FT_COMPUTE_LEADFIELD

ft_hastoolbox('bemcp', 1);

% get the optional arguments
hdmfile         = keyval('hdmfile', varargin);
conductivity    = keyval('conductivity', varargin);

% start with an empty volume conductor
vol = [];

if ~isempty(hdmfile)
  hdm = ft_read_vol(hdmfile);
  % copy the boundary of the head model file into the volume conduction model
  vol.bnd = hdm.bnd;
  if isfield(hdm, 'cond')
    % also copy the conductivities
    vol.cond = hdm.cond;
  end
else
  % copy the boundaries from the geometry into the volume conduction model
  vol.bnd = geom.bnd;
end

% determine the number of compartments
numboundaries = length(vol.bnd);

if ~isfield(vol, 'cond')
  % assign the conductivity of each compartment
  vol.cond = conductivity;
end

% determine the nesting of the compartments
nesting = zeros(numboundaries);
for i=1:numboundaries
  for j=1:numboundaries
    if i~=j
      % determine for a single vertex on each surface if it is inside or outside the other surfaces
      curpos = vol.bnd(i).pnt(1,:); % any point on the boundary is ok
      curpnt = vol.bnd(j).pnt;
      curtri = vol.bnd(j).tri;
      nesting(i,j) = bounding_mesh(curpos, curpnt, curtri);
    end
  end
end

if sum(nesting(:))~=(numboundaries*(numboundaries-1)/2)
  error('the compartment nesting cannot be determined');
end

% for a three compartment model, the nesting matrix should look like
%    0 1 1     the first is nested inside the 2nd and 3rd, i.e. the inner skull
%    0 0 1     the second is nested inside the 3rd, i.e. the outer skull
%    0 0 0     the third is the most outside, i.e. the skin
[~, order] = sort(-sum(nesting,2));

fprintf('reordering the boundaries to: ');
fprintf('%d ', order);
fprintf('\n');

% update the order of the compartments
vol.bnd    = vol.bnd(order);
vol.cond   = vol.cond(order);
vol.skin_surface   = numboundaries;
vol.source = 1;

% do some sanity checks
if length(vol.bnd)~=3
  error('this only works for three surfaces');
end
if vol.skin_surface~=3
  error('the third surface should be the skin');
end
if vol.source~=1
  error('the first surface should be the inside of the skull');
end

% Build Triangle 4th point
vol = triangle4pt(vol);

% 2. BEM model estimation, only for the scalp surface

defl =[ 0 0 1/size(vol.bnd(vol.skin_surface).pnt,1)];
% ensure deflation for skin surface, i.e. average reference over skin

% NOTE:
% Calculation proceeds by estimating each submatrix C_ij and combine them.
% There are 2 options:
% - calculating the matrices once, as it takes some time, keep them in
%   memory and use them the 2-3 times they're needed.
% - calculating the matrices every time they're needed, i.e. 2-3 times
% The former option is faster but requires more memory space as up to *8*
% square matrices of size C_ij have to be kept in memory at once.
% The latter option requires less memory, but would take much more time to
% estimate.
% This faster but memory hungry solution is implemented here.

% Deal first with surface 1 and 2 (inner and outer skull
%--------------------------------

% NOTE:
% C11st/C22st/C33st are simply the matrix C11/C22/C33 minus the identity
% matrix, i.e. C11st = C11-eye(N)

weight = (vol.cond(1)-vol.cond(2))/((vol.cond(1)+vol.cond(2))*2*pi);
C11st  = bem_Cii_lin(vol.bnd(1).tri,vol.bnd(1).pnt, weight,defl(1),vol.bnd(1).pnt4);
weight = (vol.cond(1)-vol.cond(2))/((vol.cond(2)+vol.cond(3))*2*pi);
C21    = bem_Cij_lin(vol.bnd(2).pnt,vol.bnd(1).pnt,vol.bnd(1).tri, weight,defl(1));
tmp1   = C21/C11st;

weight = (vol.cond(2)-vol.cond(3))/((vol.cond(1)+vol.cond(2))*2*pi);
C12    = bem_Cij_lin(vol.bnd(1).pnt,vol.bnd(2).pnt,vol.bnd(2).tri, weight,defl(2));
weight = (vol.cond(2)-vol.cond(3))/((vol.cond(2)+vol.cond(3))*2*pi);
C22st  = bem_Cii_lin(vol.bnd(2).tri,vol.bnd(2).pnt, weight,defl(2),vol.bnd(2).pnt4);
tmp2   = C12/C22st;

% Try to spare some memory:
tmp10 = - tmp2 * C21 + C11st;
clear C21 C11st
tmp11 = - tmp1 * C12 + C22st;
clear C12 C22st

% Combine with the effect of surface 3 (scalp) on the first 2
%------------------------------------------------------------
weight = (vol.cond(1)-vol.cond(2))/(vol.cond(3)*2*pi);
C31    = bem_Cij_lin(vol.bnd(3).pnt,vol.bnd(1).pnt,vol.bnd(1).tri, weight,defl(1));
%   tmp4   = C31/(- tmp2 * C21 + C11st );
%   clear C31 C21 C11st
tmp4 = C31/tmp10;
clear C31 tmp10

weight = (vol.cond(2)-vol.cond(3))/(vol.cond(3)*2*pi);
C32    = bem_Cij_lin(vol.bnd(3).pnt,vol.bnd(2).pnt,vol.bnd(2).tri, weight,defl(2));
%   tmp3   = C32/(- tmp1 * C12 + C22st );
%   clear  C12 C22st C32
tmp3 = C32/tmp11;
clear C32 tmp11

tmp5 = tmp3*tmp1-tmp4;
tmp6 = tmp4*tmp2-tmp3;
clear tmp1 tmp2 tmp3 tmp4

% Finally include effect of surface 3 on the others
%--------------------------------------------------
% As the gama1 intermediate matrix is built as the sum of 3 matrices, I can
% spare some memory by building them one at a time, and summing directly
weight = vol.cond(3)/((vol.cond(1)+vol.cond(2))*2*pi);
Ci3    = bem_Cij_lin(vol.bnd(1).pnt,vol.bnd(3).pnt,vol.bnd(3).tri, weight,defl(3));
gama1  = - tmp5*Ci3; % gama1 = - tmp5*C13;

weight = vol.cond(3)/((vol.cond(2)+vol.cond(3))*2*pi);
Ci3    = bem_Cij_lin(vol.bnd(2).pnt,vol.bnd(3).pnt,vol.bnd(3).tri, weight,defl(3));
gama1  = gama1 - tmp6*Ci3; % gama1 = - tmp5*C13 - tmp6*C23;

weight = 1/(2*pi);
Ci3    = bem_Cii_lin(vol.bnd(3).tri,vol.bnd(3).pnt, weight,defl(3),vol.bnd(3).pnt4);
gama1  = gama1 - Ci3; % gama1 = - tmp5*C13 - tmp6*C23 - C33st;
clear Ci3

% Build system matrix
%--------------------
i_gama1 = inv(gama1);
vol.mat = [i_gama1*tmp5 i_gama1*tmp6 i_gama1];

% remember the type
vol.type = 'bemcp';


text:
function [Rob,Lmk,Trj,Frm,Fac] = addKeyFrame(Rob,Lmk,Trj,Frm,Fac,factorRob, factorType)

% ADDKEYFRAME Add key frame to trajectory.
%   [Rob,Lmk,Trj,Frm,Fac] = ADDKEYFRAME(Rob,Lmk,Trj,Frm,Fac,factorRob,factorType)
%   adds a new frame Frm at the head of trajectory Trj. This new frame has
%   the pose of Rob. If the trajectory is full, the oldest frame is
%   deleted, and the graph is updated with the removal of all the necessary
%   pointers, factors, and, eventually, also the landmarks that have become
%   orphelin of any factor. It also creates the necessary factors linking
%   the new frame to the graph, and reserves space in the Map for the error
%   states of the new frame.
%
%   See also PRINTGRAPH, CHECKGRAPHINTEGRITY, ADDFRMTOTRJ.

%   Copyright 2015 Joan Sola @ IRI-UPC-CSIC.

% Add frame to trajectory
[Lmk,Trj,Frm,Fac] = addFrmToTrj(...
    Lmk,...
    Trj,...
    Frm,...
    Fac);

% Update new frame with Rob info
[Rob, Frm(Trj.head)] = rob2frm(...
    Rob,...
    Frm(Trj.head));

% Create motion factor
fac = find([Fac.used] == false, 1, 'first');

switch factorType

    case 'absolute'
        [Frm(Trj.head), Fac(fac)] = makeAbsFactor(...
            Frm(Trj.head), ...
            Fac(fac), ...
            Rob);

    case 'motion'
        head_old = mod(Trj.head - 2, Trj.maxLength) + 1;
        [Frm(head_old), Frm(Trj.head), Fac(fac)] = makeMotionFactor(...
            Frm(head_old), ...
            Frm(Trj.head), ...
            Fac(fac), ...
            factorRob);

    otherwise
        error('??? Unknown or invalid factor type ''%s''.', facType)
end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
y               = load('data/purse.dat')';
model           = [ssm_poisson ssm_llm];
randn('state', [1105946959; 3715058465]);
[model logL]    = estimate(y, model, exp(-6), [], 'fmin', 'bfgs', 'disp', 'iter');
fprintf(1, 'Loglikelihood: %g\nEta variance: %g\n', logL, model.param);

[alpha irr]     = fastsmo(y, model);
figure('Name', 'Purse data w/ poisson analysis');
plot(y, 'r:', 'DisplayName', 'Purse data'), hold all, plot(alpha, 'b', 'DisplayName', 'Estimated signal'), hold off; legend('show');

repo: REST


text:
function m = read_genesis_image(fname)
% m = read_ge(fname)

%
% read_genesis_image.m
%
% Original Author: Bruce Fischl
% CVS Revision Info:
%    $Author: nicks $
%    $Date: 2011/03/02 00:04:12 $
%    $Revision: 1.3 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

fid = fopen(fname, 'r', 'b');
if fid < 0
	disp(sprintf('read_genesis_image: error opening file %s', fname));
	return;
end

fseek(fid, 4, 'bof');
image_offset = fread(fid, 1, 'int');
rows = fread(fid, 1, 'int');
cols = fread(fid, 1, 'int');

fseek(fid, image_offset, 'bof');
m = fread(fid, [rows cols], 'short');

fclose(fid);

% eof

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2014 National Renewable Energy Laboratory and National
% Technology & Engineering Solutions of Sandia, LLC (NTESS).
% Under the terms of Contract DE-NA0003525 with NTESS,
% the U.S. Government retains certain rights in this software.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
% http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% wecSimTest runs WEC-Sim continuous integration testing suite.
%
%   results = wecSimTest returns a matlab.unittest.TestResult object
%
%   results = wecSimTest(..., NAME1, VALUE1, NAME2, VALUE2, ...) returns
%   a matlab.unittest.TestResult object, depending on the values of the
%   optional parameter name/value pairs. See Parameters below.
%
%   Parameters
%   ----------
%   'bemioTest'       Run tests for BEMIO. Default is true.
%
%   'regressionTest'  Run regression tests. These tests check if values
%                     have changed from a previous run. Default is true.
%
%   'compilationTest' Run compilation tests. These tests do not check
%                     correctness of the results. Default is true.
%
%   'runFromSimTest'  Run tests for execution of WECSim from Simulink.
%                     Default is true.
%
%   'rotationTest'    Run rotation tests. Default is true.
%
%   Users should also run the appropriate applications tests when
%   creating a PR into the WEC-Sim repository.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function results = wecSimTest(options)
    arguments
        options.bemioTest = true
        options.regressionTest = true
        options.compilationTest = true
        options.runFromSimTest = true
        options.rotationTest = true
    end
    % Import MATLAB unittest
    import matlab.unittest.TestSuite
    import matlab.unittest.Test
    import matlab.unittest.TestRunner
    import matlab.unittest.plugins.DiagnosticsRecordingPlugin
    import matlab.unittest.plugins.CodeCoveragePlugin
    % Define the tests
    suites = Test.empty();
    if options.bemioTest
        suites = [suites TestSuite.fromFile('tests/bemioTest.m')];
    end
    if options.regressionTest
        suites = [suites TestSuite.fromFile('tests/regressionTest.m')];
    end
    if options.runFromSimTest
        suites = [suites TestSuite.fromFile('tests/runFromSimTest.m')];
    end
    if options.rotationTest
        suites = [suites TestSuite.fromFile('tests/rotationTest.m')];
    end
    % Create TestRunner
    runner = TestRunner.withTextOutput; % Contains TestRunProgressPlugin, DiagnosticsOutputPlugin
    runner.addPlugin(DiagnosticsRecordingPlugin);
    runner.addPlugin(CodeCoveragePlugin.forFolder('./source','IncludingSubfolders',true));
    % Run the tests
    results = runner.run(suites);
    results.table
end

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
function x = dge_sl ( n, a_lu, pivot, b, job )

%*****************************************************************************80
%
%% DGE_SL solves a system factored by DGE_FA.
%
%  Discussion:
%
%    The DGE 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.
%
%    DGE_SL is a simplified version of the LINPACK routine DGESL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, real A_LU(N,N), the LU factors from DGE_FA.
%
%    Input, integer PIVOT(N), the pivot vector from DGE_FA.
%
%    Input, real B(N), the right hand side vector.
%
%    Input, integer JOB, specifies the operation.
%    0, solve A * x = b.
%    nonzero, solve A' * x = b.
%
%    Output, real X(N), the solution vector.
%
  x(1:n) = b(1:n);
%
%  Solve A * x = b.
%
  if ( job == 0 )
%
%  Solve PL * Y = B.
%
    for k = 1 : n-1

      l = pivot(k);

      if ( l ~= k )
        t    = x(l);
        x(l) = x(k);
        x(k) = t;
      end

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

    end
%
%  Solve U * X = Y.
%
    for k = n : -1 : 1
      x(k) = x(k) / a_lu(k,k);
      x(1:k-1) = x(1:k-1) - a_lu(1:k-1,k)' * x(k);
    end
%
%  Solve A' * X = B.
%
  else
%
%  Solve U' * Y = B.
%
    for k = 1 : n
      x(k) = ( x(k) - x(1:k-1) * a_lu(1:k-1,k) ) / a_lu(k,k);
    end
%
%  Solve ( PL )' * X = Y.
%
    for k = n-1 : -1 : 1

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

      l = pivot(k);

      if ( l ~= k )
        t    = x(l);
        x(l) = x(k);
        x(k) = t;
      end

    end

  end

  return
end

author: johannesgerer


text:
function [ar,e,dc]=v_lpccovar(s,p,t,w)
%V_LPCCOVAR performs covariance LPC analysis [AR,E,DC]=(S,P,T)
%
%  Inputs:  S(NS)    is the input signal
%           P        is the order (default: 12)
%           T(NF,:)  specifies the frames size details: each row specifies one frame
%                    T can be a cell array if rows have unequal numbers of values
%                       T(:,1) gives the start of the analysis interval: must be >P
%                       T(:,2) gives the end of the anaylsis interval [default: t(:+1,1)-1]
%                       subsequent pairs can be used to specify multiple disjoint segments
%                    If T is omitted, T(1,1)=P+1, T(1,2)=NS;
%                    The elements of t need not be integers.
%           W(NS)    The error at each sample is weighted by W^2 (default: 1)
%
% Outputs:  AR(NF,P+1)  are the AR coefficients with AR(:,1) = 1
%           E(NF,4)     each row is [Er Es Pr Ps] and gives the energy ("E") and power ("P")
%                       in the input signal window ("s") and in the LPC residual "r".
%                       The 'gain' of the LPC filter is g=sqrt(Pr); x=filter(g,ar,randn(:,1)) will
%                       generate noise with approximately the same power spectrum as the input s.
%           DC          is the DC component of the signal S. If this output is included,
%                       the LPC equations are modified to include a DC offset.

% Notes:
%
% (1a) If no DC output is specified AR(j,:)*S(n-(0:P)) ~ 0 or, equivalently,
%      S(n) ~ -AR(j,2:P)*S(n-(1:P)) where T(j,1) <= n <= T(j,2).
% (1b) If a DC output is specified AR(j,:)*(S(n-(0:P))-DC) ~ 0 or, equivalently,
%      S(n) ~ DC - AR(j,2:P)*(S(n-(1:P))-DC) = DC*sum(AR,j,:)) - AR(j,2:P)*S(n-(1:P))
%      where T(j,1) <= n <= T(j,2).
%
% (2) For speech processing P should be at least 2*F*L/C where F is the sampling
%     frequency, L the vocal tract length and C the speed of sound. For a typical
%     male (l=17 cm) this gives f/1000.
%
% (3) Each analysis frame should contain at least 2P samples. If note (1) is followed
%     this implies at least 2 ms of speech signal per frame.
%
% (4) It can be advantageous to restrict the analysis regions to time intervals
%     when the glottis is closed (closed-phase analysis). This can be achieved by
%     setting the T input parameter appropriately. If the closed-phase is shorter than
%     2 ms then two or more successive closed-phases should be used by defining 4 or more
%     elements in the corresponding row of T.
%
% (5) A previous version of this routine allowed T() to have a single row which would
%     be replicated for the entire file length. This has been removed because it gave rise
%     to an ambiguity.

%  Bugs: should really detect a singular matrix and reduce the order accordingly

%	   Copyright (C) Mike Brookes 1995
%      Version: $Id: v_lpccovar.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
s = s(:); % make it a column vector
if nargin < 2 p=12; end;
if nargin < 3 t=[p+1 length(s)]; end;
wq = nargin>3;
[nf,ng]=size(t);
if iscell(t)
    t{nf+1}=length(s)+1;
else
    if rem(ng,2)
        t(:,end+1)=[t(2:nf,1)-1; length(s)];
    end
end
ar=zeros(nf,p+1);
ar(:,1)=1;
e=zeros(nf,4);
dc=zeros(nf,1);
d0=nargout >2;
rs=(1:p);
for jf=1:nf
    if iscell(t)
        tj=t{jf};
        if rem(length(tj),2)
            tj(end+1)=t{jf+1}(1)-1;
        end
    else
        tj=t(jf,:);
    end

    ta = ceil(tj(1));
    tb = floor(tj(2));
    cs = (ta:tb).';
    for js=3:2:length(tj)
        ta = ceil(tj(js));
        tb = floor(tj(js+1));
        cs = [cs; (ta:tb).'];
    end
    %disp(cs([logical(1); (cs(2:end-1)~=cs(1:end-2)+1)|(cs(2:end-1)~=cs(3:end)-1); logical(1)])');
    nc = length(cs);
    pp=min(p,nc-d0);
    dm=zeros(nc,pp);	% predefine shape
    dm(:) = s(cs(:,ones(1,pp))-rs(ones(nc,1),1:pp));
    if nargout>2
        if wq
            dm = [ones(nc,1) dm].*w(cs(:,ones(1,1+pp)));
            sc=(s(cs).*w(cs));
            aa = (dm\sc).';
        else
            dm = [ones(nc,1) dm];
            sc=s(cs);
            aa = (dm\sc).';
        end
        ar(jf,2:pp+1) = -aa(2:pp+1);
        e(jf,1)=sc.'*(sc - dm*aa.');
        e(jf,2)=sc.'*sc;
        e(jf,3:4)=e(jf,1:2)/nc;
        dc(jf)=aa(1)/sum(ar(jf,:));
    else
        if wq
            dm = dm.*w(cs(:,ones(1,pp)));
            sc=(s(cs).*w(cs));
            aa = (dm\sc).';
        else
            sc=s(cs);
            aa = (dm\sc).';
        end;
        ar(jf,2:pp+1) = -aa;
        if nargout~=1
            e(jf,1)=real(sc'*(sc - dm*aa.'));
            e(jf,2)=real(sc'*sc);
            e(jf,3:4)=e(jf,1:2)/nc;
        end
    end
end
if ~nargout
    v_lpcar2ff(repmat(sqrt(e(:,3).^(-1)),1,p+1).*ar,255);
    ylabel('Power (dB)');
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function x = spm_invGcdf(F,h,l,tol)
% Inverse Cumulative Distribution Function (CDF) of Gamma distribution
% FORMAT x = spm_invGcdf(p,h,l,tol)
%
% F   - CDF (lower tail p-value)
% h   - Shape parameter (h>0)
% l   - Scale parameter (l>0)
% x   - Gamma ordinates at which CDF F(x)=F
% tol - tolerance [default 10^-6]
%__________________________________________________________________________
%
% spm_invGcdf implements the inverse Cumulative Distribution Function
% for Gamma distributions.
%
% Definition:
%--------------------------------------------------------------------------
% The Gamma distribution has shape parameter h and scale l, and is
% defined for h>0 & l>0 and for x in [0,Inf) (See Evans et al., Ch18,
% but note that this reference uses the alternative parameterisation of
% the Gamma with scale parameter c=1/l)
%
% The Cumulative Distribution Function (CDF) F(x) is the probability
% that a realisation of a Gamma random variable X has value less than
% x. F(x)=Pr{X<x}:  (See Evans et al., Ch18)
%
% Variate relationships: (Evans et al., Ch8 & Ch18)
%--------------------------------------------------------------------------
% For natural (strictly +ve integer) shape h this is an Erlang distribution.
%
% The Standard Gamma distribution has a single parameter, the shape h.
% The scale taken as l=1.
%
% The Chi-squared distribution with v degrees of freedom is equivalent
% to the Gamma distribution with scale parameter 1/2 and shape parameter v/2.
%
% Algorithm:
%--------------------------------------------------------------------------
% Interval bisection to find the inverse CDF to within tol.
%
% References:
%--------------------------------------------------------------------------
% Evans M, Hastings N, Peacock B (1993)
%       "Statistical Distributions"
%        2nd Ed. Wiley, New York
%
% Abramowitz M, Stegun IA, (1964)
%       "Handbook of Mathematical Functions"
%        US Government Printing Office
%
% Press WH, Teukolsky SA, Vetterling AT, Flannery BP (1992)
%       "Numerical Recipes in C"
%        Cambridge
%__________________________________________________________________________
% Copyright (C) 1992-2011 Wellcome Trust Centre for Neuroimaging

% Andrew Holmes
% $Id: spm_invGcdf.m 4182 2011-02-01 12:29:09Z guillaume $

%-Parameters
%--------------------------------------------------------------------------
Dtol   = 10^-8;
maxIt = 10000;

%-Format arguments, note & check sizes
%--------------------------------------------------------------------------
if nargin<4, tol = Dtol; end
if nargin<3, error('Insufficient arguments'), end

ad = [ndims(F);ndims(h);ndims(l)];
rd = max(ad);
as = [[size(F),ones(1,rd-ad(1))];...
      [size(h),ones(1,rd-ad(2))];...
      [size(l),ones(1,rd-ad(3))]];
rs = max(as);
xa = prod(as,2)>1;
if sum(xa)>1 && any(any(diff(as(xa,:)),1))
    error('non-scalar args must match in size');
end

%-Computation - Undefined and special cases
%--------------------------------------------------------------------------
%-Initialise result to zeros
x = zeros(rs);

%-Only defined for F in [0,1] & strictly positive h & l.
% Return NaN if undefined.
md = ( F>=0  &  F<=1  &  h>0  &  l>0 );
if any(~md(:))
    x(~md) = NaN;
    warning('Returning NaN for out of range arguments');
end

%-Special cases: x=0 when F=0, x=Inf when F=1
x(md & F==1) = Inf;

%-Compute where defined & not special case
%--------------------------------------------------------------------------
Q  = find( md  &  F>0  &  F<1 );
if isempty(Q), return, end
if xa(1), FQ=F(Q); FQ=FQ(:); else FQ=F*ones(length(Q),1); end
if xa(2), hQ=h(Q); hQ=hQ(:); else hQ=h*ones(length(Q),1); end
if xa(3), lQ=l(Q); lQ=lQ(:); else lQ=l*ones(length(Q),1); end
%-?Q=?Q(:) stuff is to avoid discrepant orientations of vector arguments!

%-Compute initial interval estimates [0,mean] & Grow to encompass F(QF);
%--------------------------------------------------------------------------
a  = zeros(length(Q),1);
b  = hQ./(2*lQ);
QQ = 1:length(Q);
while ~isempty(QQ)
    b(QQ) = 2*b(QQ);
    QQ    = find(spm_Gcdf(b,hQ,lQ) < FQ);
end

%-Interval bisection
%--------------------------------------------------------------------------
fa = spm_Gcdf(a,hQ,lQ)-FQ;
fb = spm_Gcdf(b,hQ,lQ)-FQ;
i  = 0;
xQ = a+(b-a)/2;
QQ = 1:length(Q);

while ~isempty(QQ) &&  i<maxIt
    fxQQ        = spm_Gcdf(xQ(QQ),hQ(QQ),lQ(QQ))-FQ(QQ);
    mQQ         = fa(QQ).*fxQQ > 0;
    a(QQ(mQQ))  = xQ(QQ(mQQ));   fa(QQ(mQQ))  = fxQQ(mQQ);
    b(QQ(~mQQ)) = xQ(QQ(~mQQ));  fb(QQ(~mQQ)) = fxQQ(~mQQ);
    xQ(QQ)      = a(QQ) + (b(QQ)-a(QQ))/2;
    QQ          = QQ( (b(QQ)-a(QQ))>tol );
    i           = i+1;
end

if i==maxIt
    warning('convergence criteria not reached - maxIt reached');
end

x(Q) = xQ;

author: spm


text:
%CODEGENERATOR.FFINDREPLACE Find and replace all occurrences of string in a file.
%
%   CGen.ffindreplace(fName, oText, nText, varargin)
%   FNAME is the absolute or relative path to the file to replace the text in.
%   OTEXT is the text passage to replace.
%   NTEXT is the new text.
%
% Notes::
% The function opens and sweeps the text file FNAME. All occurrences of
% OTEXT are replaced by NTEXT.
%
% Authors::
%  Joern Malzahn, (joern.malzahn@tu-dortmund.de)
%
% See also CodeGenerator.finsertfront.

% Copyright (C) 1993-2012, by Peter I. Corke
% Copyright (C) 2012-2013, by Joern Malzahn
%
% This file is part of The Robotics Toolbox for Matlab (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
%
% The code generation module emerged during the work on a project funded by
% the German Research Foundation (DFG, BE1569/7-1). The authors gratefully
% acknowledge the financial support.

function [ ] = ffindreplace(CGen, fName, oText, nText, varargin)

fid = fopen(fName,'r'); % open the file

% determine number of lines in the file
nLines = getNLines(fid);

% read all lines to cell
oldLines = cell(1,nLines);
newLines = {1,nLines};
for iLines = 1:nLines
    oldLines{iLines} = fgets(fid);
end

% close the file again
fclose(fid);

% replace all occurrences of oText by nText in each line
for iLines = 1:nLines
    newLines{iLines}= strrep(oldLines{iLines}, oText, nText);
end

% rewrite the file
fid = fopen(fName,'w');
for iLines = 1:nLines
   fprintf(fid,'%s',newLines{iLines});
end
fclose(fid);

end

function [nLines] = getNLines(fid)
nLines = 0;
while (~feof(fid))      % go through each line until eof is reached
    fgets(fid);
    nLines = nLines +1; % increment line counter
end
frewind(fid);           % set file indicator back to the beginning of the file
end

author: petercorke

repo: robotics-toolbox-matlab


text:
classdef Enableable < handle
    % Mixin for component with Enable property

    % Copyright 2020-2021 The MathWorks Inc.

    %% Properties
    properties (AbortSet)

        % Enable of the component
        Enable (1,1) matlab.lang.OnOffSwitchState = 'on'

    end %properties

    %% Internal properties
    properties (AbortSet, Transient, NonCopyable, ...
            Access = {?wt.abstract.BaseWidget, ?wt.test.BaseWidgetTest} )

        % List of graphics controls to apply to
        EnableableComponents (:,1) matlab.graphics.Graphics

    end %properties

    %% Accessors
    methods

        function set.Enable(obj,value)
            obj.Enable = value;
            obj.updateEnableableComponents()
        end

        function set.EnableableComponents(obj,value)
            obj.EnableableComponents = value;
            obj.updateEnableableComponents()
        end

    end %methods

    %% Methods
    methods (Access = protected)

        function updateEnableableComponents(obj)

            wt.utility.fastSet(obj.EnableableComponents,"Enable",obj.Enable);

        end %function

    end %methods

end %classdef

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ params ] = rnn_stack2params( stack, eI, W_t, sum_tied )
%RNN_STACK2PARAMS converts stack structure of RNN weights to single vector
%   Takes a stack strcutre with stack{l}.W and stack{l}.b for each layer
%   Also takes single matrix of temporal weights W_t
%   The flag sum_tied will sum tied encoder and decoder weights.
%       This is useful for gradient aggregation
%   Verifies the stack structure conforms to their descriptions in eI
%   Namely checks eI.layerSizes, eI.inputDim, and eI.temporalLayer

%% assume no weight tieing if parameter unset
if ~isfield(eI, 'tieWeights')
    eI.tieWeights = 0;
end;
if ~exist('sum_tied','var')
    sum_tied = false;
end;
%% default short circuits to false
if ~isfield(eI, 'shortCircuit')
    eI.shortCircuit = 0;
end;
% check short circuit consistency
assert( ~xor(eI.shortCircuit, isfield(stack{end},'W_ss')));
%% check first layer dimensions
assert( size(stack{1}.W,1) == eI.layerSizes(1));
assert( size(stack{1}.W,2) == eI.inputDim);
assert( size(stack{1}.b,1) == eI.layerSizes(1));
%% stack first layer
params = [ stack{1}.W(:); stack{1}.b(:)];
%% check and stack all layers. no special treatment of output layer
for l = 2 : numel(eI.layerSizes)
    assert( size(stack{l}.W,1) == eI.layerSizes(l));
    assert( size(stack{l}.W,2) == eI.layerSizes(l-1));
    assert( size(stack{l}.b,1) == eI.layerSizes(l));
    if ~eI.tieWeights || (l <= numel(eI.layerSizes)/2 ...
            || l == numel(eI.layerSizes))
        % untied layer, save the weights
        if eI.tieWeights && sum_tied && l < numel(eI.layerSizes)
            % sum decoder weights if its a tied encoder layer
            lDec = numel(eI.layerSizes) - l + 1 ;
            params = [ params; reshape(stack{l}.W + stack{lDec}.W',[],1)];
        else
            params = [ params; stack{l}.W(:)];
        end;
    end;
    % always aggregate bias
    params = [ params; stack{l}.b(:)];
end
%% append temporal weight matrix
if ~isempty(W_t) || eI.temporalLayer
    assert(size(W_t,1) == eI.layerSizes(eI.temporalLayer));
    assert(size(W_t,2) == eI.layerSizes(eI.temporalLayer));
    params = [ params; W_t(:)];
end
%% append short circuit matrix
if eI.shortCircuit
    params = [params; stack{end}.W_ss(:)];
end;

author: amaas

repo: rnn-speech-denoising


text:
function fmdplot(obj, tabgrouptag)
    % control the FMD plot in the main window

    Tags.xs = 'fmd xs contextmenu';

    % the fmd object is actually stored in the axes into which it plots
    myTab = findOrCreateTab(obj.fig, obj.fig, tabgrouptag, 'FMD');
    ax=findobj(myTab,'Type','axes');
    if isempty(ax)
        ax=axes(myTab);
        ylabel(ax,'Cum # events > M(x)');
        xlabel(ax,'Magnitude');
        if ~isempty(obj.catalog)
            bdiffobj=bdiff2(obj.getCurrentCatalog,'ax',ax,'InteractiveMode',false);
            bAnalysisWin = AnalysisBvalues(ax,bdiffobj);
            ax.UserData = bAnalysisWin;
            % ax.UserData=bdiffobj; %stash this, but keep it with the ZMapMainWindow.
            uimenu(ax.UIContextMenu,'Label','Cut catalog at Mc','MenuSelectedFcn', @(s,v)crop_to_mc(s, v, bdiffobj) );
        end

    elseif isempty(obj.catalog)
        cla(ax);
        ylabel(ax,'Cum # events > M(x)');
        xlabel(ax,'Magnitude');
    else
        % axes exists, and the catalog is not empty

        bAnalysisWin = ax.UserData;
        if isempty(bAnalysisWin)
            bdiffobj=bdiff2(obj.getCurrentCatalog,'ax',ax,'InteractiveMode',false);
            bAnalysisWin = AnalysisBvalues(ax,bdiffobj);
            ax.UserData = bAnalysisWin;
            % ax.UserData=bdiffobj;
        else
            bAnalysisWin.bobj.RawCatalog = obj.catalog;
            bAnalysisWin.bobj.Calculate();
            bAnalysisWin.bobj.updatePlot();
        end

        if isempty(ax.UIContextMenu)
            c = uicontextmenu(obj.fig);
            ax.UIContextMenu = c;
        end
        addLegendToggleContextMenuItem(ax.UIContextMenu,'bottom','above');
    end

    analy_win = ax.UserData;

    %% plot & synchronize cross sections

    plotted_xs = get(ax.Children,'Tag');
    if ~isempty(plotted_xs)
        plotted_xs = plotted_xs(startsWith(plotted_xs, 'Xsection fmd'));
    end
    if ~isempty(plotted_xs)
        existing_xs = obj.XSectionTitles;
        plotted_xs2 = extractAfter(plotted_xs,'Xsection fmd ');

        toDelIdx = ~startsWith(plotted_xs2,existing_xs); % all tags
        todel = string(plotted_xs(toDelIdx));
        todel = todel(~endsWith(todel,'line') & ~endsWith(todel,'Mc'));

        % delete cross sections that shouldn't exist
        %todel = plotted_xs(~startsWith(plotted_xs,existing_xs));
        %todel = plotted_xs(~ismember(plotted_xs,existing_xs));
        analy_win.remove_series(todel);
    end

    %% if necessary, add context menu to figure
    cxs=findobj(obj.fig,'Tag',Tags.xs);
    if isempty(cxs)
        cxs=uicontextmenu(obj.fig,'tag',Tags.xs);
        uimenu(cxs,'Label','Open in new window','MenuSelectedFcn',@cb_xstimeplot);
    end

    obj.plot_xsections(@xsplotter, 'Xsection fmd');
    clear_empty_legend_entries(obj.fig);

    %%
    function h = xsplotter(xs, xscat)
        xsProps.LineWidth   = 1.5;
        xsProps.DisplayName = xs.Name;
        xsProps.Color       = xs.Color;
        xsProps.Marker = 'x';
        mytag               = ['Xsection fmd ' xs.Name];
        h = analy_win.add_series(xscat, mytag, xsProps);

    end

    function cb_xstimeplot(~,~)
        % CB_XSTIMEPLOT shows the TIMEPLOT for the currently selected cross-section
        myName = get(gco,'DisplayName');
        idx = strcmp(obj.XSectionTitles,myName);
        ZG=ZmapGlobal.Data;
        ZG.newt2=obj.catalog.subset(obj.CrossSections(idx).inside(obj.catalog));
        ZG.newt2.Name=sprintf('Events within %g km of %s',obj.CrossSections(idx).name);
    end

    function crop_to_mc(~,~, bdiffobj)
        % should this crop the raw catalog?
        zdlg = ZmapDialog();
        zdlg.AddEdit('mc', "Cut Magnitude [Mc:" + bdiffobj.Result.Mc_value + "]", bdiffobj.Result.Mc_value,...
            'Choose magnitude to cut the catalog');
        [res,okpressed] = zdlg.Create('Name', 'Choose Cut Magnitude');
        if okpressed
            obj.rawcatalog = obj.rawcatalog.subset(res.mc <= obj.rawcatalog.Magnitude);
            obj.CatalogManager.RawCatalog = obj.rawcatalog;
            obj.replot_all;
        end
    end
end

author: CelsoReyes


text:
function maskM = fastPolyFill(pointsM,optS)
%function maskM = fastPolyFill(pointsM,optS)
%fastPolyFill:  fills in images given polygons defining edges.
%Intention is to only fill a voxel if it's center is in the polygon.
%Faster than previous version by a significant factor (about x10).
%JOD; 17 June 05.
%Fixed bug, 5 July 05.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

imageSizeV = optS.ROIImageSize;

xOffset = optS.xCTOffset;
yOffset = optS.yCTOffset;

numRows = double(imageSizeV(1));
numCols = double(imageSizeV(2));

xInV = pointsM(:,1);
yInV = pointsM(:,2);

maskM = zeros(numRows,numCols);
%could also be a double.

%preallocate a table of voxel indices
%'next to' refers to voxel centers just to the right of where polygon
%edges cross a line connecting voxel centers.
nextToPtsM = sparse(numRows * 4, numCols * 4); %hard to imagine bigger (sparse is faster than full)
numNextToPtsV = zeros(numRows,1);

%convert to "row and col space", that is, a continuous space where
%s is a coord that runs from 1 to numCols (along x axis)
%t is a coord that runs from 1 to numRows (along -y axis)
sInV = (xInV - xOffset)/optS.ROIxVoxelWidth + (numCols + 1)/2;
tInV = (yInV - yOffset)/optS.ROIyVoxelWidth + (numRows + 1)/2;

% % APA: snap points to integer-grid using nearest neighbor interpolation and
% % delete duplicates
% [S,T] = meshgrid(1:numCols+1,1:numRows+1);
% indSnap = interp2(S,T,reshape(1:(numCols+1)*(numRows+1),numRows+1,numCols+1),sInV,tInV,'*nearest');
% STpts = [S(ind2sub([numCols+1 numRows+1],indSnap)) T(ind2sub([numCols+1 numRows+1],indSnap))];
% goodPtsV = [1;any(diff(STpts),2)];
% sInV = STpts(find(goodPtsV),1);
% tInV = STpts(find(goodPtsV),2);

%For each polygonal element, compute and add one to 'row starters,' centers of
%voxels (assumed at image coords) whose rows cross the polygonal elements
%and are the first element to the right of that polygonal element's row crossing.
%
%Use the parameterization col = n * deltaRow + c + rowURPt    (URPt = upper
%right point).
%
%Algorithm:
%1.  Determine n
%2.  Determine c
%3.  Determine vector of row crossings
%4.  Determine vector of column crossings
%5.  Add one to row starters.

%Loop over polygonal elements
%For polygonal elements:

shift_sV = [sInV(end); sInV(1:end-1)];
shift_tV = [tInV(end); tInV(1:end-1)];

sM = [sInV(:), shift_sV];
tM = [tInV(:), shift_tV];

%Loop over polygonal edges, put ones where voxel centers are inside
%polygon.
for i = 1 : length(sInV)

    if tM(i,1) ~= tM(i,2)   %skip horizontal lines

        tMax = max([tM(i,1),tM(i,2)]);
        tMin = min([tM(i,1),tM(i,2)]);

        %determine n (line parameterization: s = n * t + c )
        n = (sM(i,1) - sM(i,2))/(tM(i,1) - tM(i,2));
        %determine c, could vectorize these two
        c = sM(i,1) - n * tM(i,1);

        %get delta_tV
        tPtsV = ceil(tMin) : floor(tMax);
%         tPtsV = tMin+1 : tMax; % APA
        delta_tV = tPtsV - tMax;
        sPtsV = n * delta_tV + c + n * tMax; %these are s values at edge 'crossings'

        %derive s values 'next to the right' of the crossings
        sVoxelsV = ceil(sPtsV);

        %catalogue
        for j = 1: length(sVoxelsV)
            num = numNextToPtsV(tPtsV(j));
            numNextToPtsV(tPtsV(j)) = num + 1;
            nextToPtsM(tPtsV(j),num+1) = sVoxelsV(j);
        end

    end

end

%This should be faster than the oft-used cumsum trick:
for i = 1 : numRows
    num = numNextToPtsV(i);
    if num ~=0
        %get 'em
        nextToPtsV = nextToPtsM(i,1:num);
        %sort 'em
        sortPtsV = sort(full(nextToPtsV));
        %fill image
        for j = 1 : length(sortPtsV)/2
            indV = sortPtsV(2*j-1) : sortPtsV(2*j) - 1;
            maskM(i * ones(1,length(indV)),indV) = 1;
        end
    end
end

%lastly, correct an oversight in writing the code: rows need to be flipped:
maskM = flipud(maskM);  %fast operation.

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR