text:
function [F,F_r,F_u] = odo3(F,u)

% ODO3 3D Odometry evolution.
%   F = ODO3(F,U) performs one step on the pose F of a vehicle, given
%   odometry increments U=[DX;DV] in robot frame.
%   - F is a frame structure (see FRAME).
%   - Position increment DX is given in robot frame F.
%   - Orientation increment DV is given as a Rotation Vector in robot frame F.
%
%   [F,F_r,F_u] = ODO3(F,U) gives the full Jacobians wrt state and odometry
%   inputs.
%
%   See also FRAME, V2Q, QPROD, QUATERNION.

%   Copyright 2005-2009 Joan Sola @ LAAS-CNRS.

dv = u(4:6);
dx = u(1:3);

if nargout == 1

    x = fromFrame(F,dx); % Position update

    q  = F.x(4:end);
    q2 = qProd(q,v2q(dv)); % quaternion update

    F.x = [x;q2]; % frame update

else  % Jacobians

    [x,X_r,X_dx] = fromFrame(F,dx); % Position update and jacobians

    q               = F.x(4:end);
    [dq,DQ_dv]      = v2q(dv);
    [q2,Q2_q,Q2_dq] = qProd(q,dq); % quaternion update
    Q2_dv           = Q2_dq*DQ_dv;

    F.x = [x;q2]; % frame update

    F_r  = [X_r;zeros(4,3) Q2_q];
    F_u  = [X_dx zeros(3,3);zeros(4,3) Q2_dv];
end

F = updateFrame(F);

return

%% Jacobians

syms x y z a b c d real
syms dx dy dz dp dq dr real
Fi.x=[x;y;z;a;b;c;d];
Dx = [dx;dy;dz];
Dv = [dp;dq;dr];
u  = [Dx;Dv];

Fi = updateFrame(Fi);

[F,F_r,F_u] = odo3(Fi,u);

F_rs = jacobian(F.x,Fi.x);
F_us = jacobian(F.x,[Dx;Dv]);

simplify(F_r-F_rs)
simplify(F_u-F_us)

% ========== 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:
%dump_pdfbimagemodel.m
%written by: Duncan Po
%Date: Sept 27, 2002
%Saves a pdfb model into a file in a readable format
%Usage: dump_pdfbimagemodel(model, filename)
%       model - model of image to be saved
%       filename - name of the file (including full path) to save model in

function dump_pdfbimagemodel(model, filename)

nl = model{1}.nlevels;

for m = 1:nl
    levndir(m) = length(model{1}.stdv{m});
end;

for i = 1:length(model)
    dump_pdfbmodel_to_file(model{i}, filename, levndir, nl);
end;

author: Sable

repo: mcbench-benchmarks


text:
function [metabolite_structure,IDsAdded] = parseBiggID4VMH(metabolite_structure,startSearch,endSearch,grebMoreIDs)
% the problem is that by chance Bigg and VMH could have the same ID but for
% different metabolites -- I do not do any additional checks right now
% which is dangerous (hence I do not greb more ID's by default)

if ~exist('grebMoreIDs','var')
    grebMoreIDs = 0;
end

% try to guess bigg ID from VMH ID
Mets = fieldnames(metabolite_structure);
if ~exist('startSearch','var')
    startSearch = 1;
end
if ~exist('endSearch','var')
    endSearch = length(Mets);
end

annotationSource = 'Bigg website (matching)';
annotationType = 'automatic';

a = 1;
IDsAdded = '';

mapping = {
    'inchiKey'     'https://identifiers.org/inchikey/'
    'hmdb'  'http://identifiers.org/hmdb'
    'metanetx'  'http://identifiers.org/metanetx.chemical'
    'keggId'    'http://identifiers.org/kegg.compound/'
    'biocyc'    'http://identifiers.org/biocyc/'
    'reactome'  'http://identifiers.org/reactome/' % for the moment I greb only the first entry
    'cheBIId'   'http://identifiers.org/chebi/'
    };
for i = startSearch : endSearch

    %get VMHId
    MetsID= metabolite_structure.(Mets{i}).VMHId;
    % try if a valid website exists
    try

        url = strcat('http://bigg.ucsd.edu/universal/metabolites/',MetsID);
        [syst,success] = urlread(url);
        if ~success %no entry found
            MetsID = regexprep(MetsID,'_','__');
            url  = strcat('http://bigg.ucsd.edu/universal/metabolites/',MetsID);
            [syst,success] = urlread(url);
        end

        if success
            metabolite_structure.(Mets{i}).biggId = MetsID;
            metabolite_structure.(Mets{i}).biggId_source = [annotationSource,':',annotationType,':',datestr(now)];
            if  grebMoreIDs == 1
                for k = 1 : size(mapping,1)
                    [metabolite_structure,idNew] = getData(metabolite_structure,syst,Mets{i}, mapping(k,:),IDsAdded);
                    if ~isempty(idNew) && isempty(find(isnan(idNew),1))
                        if isempty(metabolite_structure.(Mets{i}).(mapping{k,1})) || ~isempty(find(isnan(metabolite_structure.(Mets{i}).(mapping{k,1})),1))
                            metabolite_structure.(Mets{i}).(mapping{k,1}) = (idNew);
                            metabolite_structure.(Mets{i}).(strcat(mapping{k,1},'_source')) = [annotationSource,':',annotationType,':',datestr(now)];
                            IDsAdded{a,1} = Mets{i};
                            IDsAdded{a,2} = mapping{k,1};
                            IDsAdded{a,3} = char(idNew);
                            a = a+1;
                        end
                    end
                end
            end

        end
    catch
        continue;
    end
end

function [metabolite_structure,idNew] = getData(metabolite_structure,syst,met, map,IDsAdded)
a = size(IDsAdded,1)+1;
idNew = '';
try
    startvalue=strfind(syst,map{1,2});
    string = syst(startvalue(1):startvalue(1)+200);
    [tok,rem] = strtok(string,'>');
    [tok2,rem2] = strtok(rem,'<');
    idNew = regexprep(tok2,'>','');
    if contains(idNew,'META:')
        idNew = regexprep(idNew,'META:','');
    end
    if contains(idNew,'CHEBI:')
        idNew = regexprep(idNew,'CHEBI:','');
    end
end

author: opencobra

repo: cobratoolbox


text:
function overtopo()
    % provide controls that would go around a topographic map

    ZG=ZmapGlobal.Data; % used by get_zmap_globals
    report_this_filefun();

    to1 = figure_w_normalized_uicontrolunits( ...
        'Name','Topographic Map',...
        'NumberTitle','off', ...
        'backingstore','on',...
        'NextPlot','add', ...
        'Visible','on', ...
        'Position',position_in_current_monitor(ZG.map_len(1), ZG.map_len(2)));

    add_controls()

    function add_controls()
        labelPos = [0.9 0.93 0.10 0.05];
        uicontrol(...
            'Style','popup',...
            'Units','normalized',...
            'Position',labelPos,...
            'Value',1,...
            'String',cb_shader(),...
            'BackgroundColor',[0.7 0.7 0.7]',...
            'callback',@cb_shader);

        labelPos = [0.9 0.75 0.10 0.05];
        uicontrol(...
            'Style','popup',...
            'Units','normalized',...
            'Position',labelPos,...
            'Value',1,...
            'String',cb_volcanoes(),...
            'BackgroundColor',[0.7 0.7 0.7]',...
            'callback',@cb_volcanoes);

        uicontrol('Units','normal',...
            'Position',[.0 .0 .28 .04],'String','Plot using Mapping Toolbox ',...
            'callback', @(~,~)zmaptopo.overmaptb)

        labelPos = [0.9 0.65 0.10 0.05];
        uicontrol(...
            'Style','popup',...
            'Units','normalized',...
            'Position',labelPos,...
            'Value',1,...
            'String',cb_eqplot(),...
            'BackgroundColor',[0.7 0.7 0.7]',...
            'callback',@cb_eqplot);

        labelPos = [0.9 0.55 0.10 0.05];
        uicontrol(...
            'Style','popup',...
            'Units','normalized',...
            'Position',labelPos,...
            'Value',1,...
            'String',cb_faultplot(),...
            'BackgroundColor',[0.7 0.7 0.7]',...
            'callback',@cb_faultplot);

        labelPos = [0.9 0.45 0.10 0.05];
        uicontrol(...
            'Style','popup',...
            'Units','normalized',...
            'Position',labelPos,...
            'Value',1,...
            'String',cb_coastplot(),...
            'BackgroundColor',[0.7 0.7 0.7]',...
            'callback',@cb_coastplot);

        labelPos = [0.9 0.35 0.10 0.05];
        uicontrol(...
            'Style','popup',...
            'Units','normalized',...
            'Position',labelPos,...
            'Value',1,...
            'String',cb_mainevent(),...
            'BackgroundColor',[0.7 0.7 0.7]',...
            'callback',@cb_mainevent);

        labelPos = [0.9 0.25 0.10 0.05];
        uicontrol(...
            'Style','popup',...
            'Units','normalized',...
            'Position',labelPos,...
            'Value',1,...
            'String',cb_stations(),...
            'BackgroundColor',[0.7 0.7 0.7]',...
            'callback',@cb_stations);
    end

    %% callbacks
    function labelList = cb_shader(mysrc, ~)
        labelList={'flat', 'interp', 'faceted'};
        if nargin>0
            shading(labelList{mysrc.Value});
        end
    end

    function labelList = cb_volcanoes(mysrc,~)
        labelList= {'Volcanoes on','off'};
        if nargin == 0
            return
        end
        switch mysrc.Value
            case 1 % on
                plot(vo.Longitude, vo.Latitude, '^r', 'LineWidth', 1.5, 'MarkerSize',8,...
                    'MarkerFaceColor', 'w', 'MarkerEdgeColor', 'r', 'Tag', 'overtopo_volcanoes');
            case 2
                delete(findobj(gcf, 'Tag', 'overtopo_volcanoes'));
        end
    end

    function labelList = cb_eqplot(mysrc,~, catalog)
        labelList = {'EQ (dot)', 'EQ (o)', 'No EQ'};
        if nargin == 0
            return
        end

        switch mysrc.Value
            case 1
                plot(catalog.Longitude,catalog.Latitude, '.r', 'MarkerSize', 1, 'Tag', 'overtopo_eq') ;
            case 2
                plot(catalog.Longitude,catalog.Latitude, 'or', 'MarkerSize', 2, 'Tag', 'overtopo_eq') ;
            case 3
                delete(findobj(gcf,'Tag','overtopo_eq'));
        end
    end

    function labelList = cb_faultplot(mysrc, ~)
        labelList = {'Faults', 'No Faults'};
        if nargin == 0
            return
        end
        switch mysrc.Value
            case 1
                plot(faults(:,1), faults(:,2), 'm', 'Linewidth', 2, 'Tag', 'overtopo_faults');
            case 2
                delete(findobj(gcf,'Tag','overtopo_faults')) ;
        end
    end

    function labelList = cb_coastplot(mysrc,~)
        labelList = {'Coast','No Coast'};
        if nargin == 0
            return
        end
        switch mysrc.Value
            case 1
                plot(coastline(:,1),coastline(:,2),'w','Linewidth',2, 'Tag', 'overtopo_coast');
            case 2
                delete(findobj(gcf,'Tag','overtopo_coast')) ;
        end
    end

    function labelList = cb_mainevent(mysrc,~)
        labelList = {' Main','No Main'};
        if nargin == 0
            return
        end
        switch mysrc.Value
            case 1
                plot(ZG.maepi.Longitude, ZG.maepi.Latitude, 'hm',...
                    'LineWidth', 1.5, 'MarkerSize', 12,...
                    'MarkerFaceColor', 'y', 'MarkerEdgeColor', 'k', 'Tag', 'overtopo_mainevts')
                te1 = text(ZG.maepi.Longitude, ZG.maepi.Latitude, stri2);
                set(te1,'FontWeight', 'bold', 'Color', 'k', 'FontSize', 9, 'Clipping', 'on',...
                    'Tag', 'overtopo_mainevts');

            case 2
                delete(findobj(gcf, 'Tag','overtopo_mainevts'));
        end
    end

    function labelList = cb_stations(mysrc,~)
        labelList = {'Stations','No Stations'};
        if mysrc.Value == 1
            h1 = h1topo;
            plotstations ;
        end
    end

end

author: CelsoReyes


text:
function [cap res] = pll_synth_3rd_order(ipump, vco_sensitivity, fout, fcomp, bandwidth, pm)

% Synthesizes PLL loop components for a 3rd order system using the
% topology shown below.
%
% Input - ipump is the charge pump current in Amperes
%       - vco_sensitivity is the VCO sensitivity in Hertz/Volt
%       - fout is the output frequency in Hertz
%       - fcomp is the comparison frequency in Hertz
%       - bandwidth is the open loop bandwidth in Hertz
%       - pm is the phase margin in degrees
% Output - cap is the capacitors of the loop in Farads [C1, C2, C3, 0]
%        - res is the resistors of the  loop in Ohms [ R2, R3, 0]
%
%
% The method used here is derived from that presented in Dean Banerjee's
% Book "PLL Performance, Simulation, and Design" 4th Ed available at
% National Semiconductors site www.national.com.  Most of the work is
% derived from Chapter 22 pp.178-185.
%
% No closed form solution exists for 3rd
% order systems, therefore some simplification of the model has been made
% to derive the results.  Due to the simplifications the results will
% also not always be 100% accurate.  It is recommended that the
% results are verified using the pll_simulation.m script, which does not
% use a simplified model.
% Loop Topology
%
%        +  _____                                   ______
% fcomp -->|Phase|-----------------R3---------------| VCO |---->fout
%          |Det. |     |      |         |           |     |  |
%           -----      |      |         |            -----   |
%            ^ -      C1     R2        C3                    |
%            |         |      |         |                    |
%            |        GND    C2        GND                   |
%            |                |                              |
%            |               GND         _______             |
%            ----------------------------| 1/N |--------------
%                                        -------
%
%   Author: Ben Gilbert
%   Homepage: http://nicta.com.au/people/gilbertb
%   Email: ben.gilbert (wibble) nicta.com.au
%   (c) 2009 by National ICT Australia (NICTA)
%   %

%% Design Parameters
% The following values will give sensible results but not necessarily
% optimal
T31 = 0.6; % T3/T1 ratio of time constants
gamma = 1.0; % optimization factor p172 and pp220-227 of Banerjee

% Conversion of parameters to more convenient units
Kpd = ipump/2/pi; % phase detector gain
Kvco = vco_sensitivity*2*pi; % vco gain
omega = 2*pi* bandwidth; % open loop bandwidth in radians/sec

%% Synthesis
%% Find T1 using numerical methods
t1 = (sec(pm*pi/180)-tan(pm*pi/180))/(omega*(1+T31)); % approximate value for T1

% T1 is the root of this equation
num1 = @(t) atan(gamma./(omega.*t*(1+T31)))-atan(omega.*t)-atan(omega.*t*T31)- pm*pi/180;
T1 = fzero(num1,t1);% find root numerically

%% Find T2 and T3
T3 = T1*T31;
T2 = gamma/omega^2/(T1+T3);

%% Solving for Loop Components
A0 = Kpd * Kvco * fcomp/fout / omega^2 * sqrt((1+omega^2*T2^2)/(1+omega^2*T1^2)/(1+omega^2*T3^2));
A1 = A0*(T1+T3);
A2 = A0*T1*T3;

C1 = A2/T2^2*(1+sqrt(1+T2/A2*(T2*A0-A1)));
C3 = (-T2^2*C1^2+T2*A1*C1-A2*A0)/(T2^2*C1-A2);
C2 = A0-C1-C3;
R2 = T2/C2;
R3 = A2/C1/C3/T2;

cap = [C1 C2 C3 0];
res = [R2 R3 0];


text:
classdef transform
% Implements different versions of wavelet transform
% n=1024, J=10, L=6
% L - J = 10 - 6 = 4.
% s_10 = s_6 + d_6 + d_7 + d_8 + d_9
% j=J-1:-1:L = [9, 8, 7, 6]
% n = 256, J=8, L=0
% s_8 = s_0 + d_0 + d_1 + d_2 + d_3 + d_4 + d_5 + d_6 + d_7.
% s_J = s_L + sum(L <= j < J) d_j.

methods(Static)

    function w = forward_periodized_orthogonal(qmf, x, L)
        % Computes the forward wavelet transform of x
        %
        % Uses the periodized version of x
        % with an orthogonal wavelet basis
        % length of x must be dyadic.
        if nargin < 3
            % We perform full wavelet decomposition
            L = 0;
        end

        % Let's get the dyadic length of x and verify that
        % length of x is a power of 2.
        [n, J, consistent] = spx.wavelet.dyad_length(x);
        if ~consistent
            error('x must be of dyadic length');
        end
        if L >= J
            error('L must be smaller than dyadic index of x');
        end
        % We will work with row vectors in this function.
        col = false;
        if iscolumn(x)
            col = true;
            % Convert it to row vector
            x = x';
        end
        % Create the storage for wavelet coefficients.
        w = zeros(1, n);
        for j=J-1:-1:L
            % Start from the finest level and keep going down.
            % identify the hipass component of x and downsample it.
            c = spx.wavelet.hi_pass_down_sample(qmf, x);
            % Identify the locations where the hipass component will be stored.
            indices = spx.wavelet.dyad(j);
            w(indices) = c;
            % Replace x with its low pass downsampled version
            x = spx.wavelet.lo_pass_down_sample(qmf, x);
        end
        % Store the remaining contents of x in the beginning of array
        w(1:(2^L)) = x;
        if col
            % Convert the wavelet coefficients to a column vector
            w = w';
        end
    end

    function x = inverse_periodized_orthogonal(qmf, w, L)
        % Computes the inverse wavelet transform of x
        %
        % Uses the periodized version of x
        % with an orthogonal wavelet basis
        % length of x must be dyadic.
        if nargin < 3
            % We perform full wavelet composition
            L = 0;
        end

        % Let's get the dyadic length of w and verify that
        % length of w is a power of 2.
        [n, J, consistent] = spx.wavelet.dyad_length(w);
        if ~consistent
            error('w must be of dyadic length');
        end
        if L >= J
            error('L must be smaller than dyadic index of w');
        end
        % We will work with row vectors in this function.
        col = false;
        if iscolumn(w)
            col = true;
            % Convert it to row vector
            w = w';
        end
        % initialize x with its coerce approximation
        x = w(1:2^L);
        for j=L:J-1
            % Identify the locations where the hipass component is stored.
            indices = spx.wavelet.dyad(j);
            c = w(indices);
            % Compute the low pass portion of the next level of approximation
            x_low = spx.wavelet.up_sample_lo_pass(qmf, x);
            if iscolumn(x_low)
                x_low = x_low';
            end
            % Compute the high pass portion of the next level of approximation
            x_hi = spx.wavelet.up_sample_hi_pass(qmf, c);
            if iscolumn(x_hi)
                x_hi = x_hi';
            end
            % Compute the next level approximation of x
            x = x_low + x_hi;
        end
        if col
            % Convert the resulting vector to a column vector
            x = x';
        end
    end

end

end

author: indigits

repo: sparse-plex


text:
function data = read_artinis_oxy55(filename, header, begsample, endsample, chanindx)
% reads Artinix oxy55-files into FieldTrip format
%
% use as
%   header = read_artinis_oxy5(filename)
% or
%   event  = read_artinis_oxy5(filename, read_event)
% where read_event is a Boolean with the value true, or
%   data   = read_artinis_oxy55(filename, header, [begsample], [endsample], [chanindx])
% where begsample, endsample and chanindx are optional.
%
% The returned variables will be in FieldTrip style.
%
% See also FT_READ_HEADER, FT_READ_DATA, READ_ARTINIS_OXU5

% You are using the FieldTrip NIRS toolbox developed and maintained by
% Artinis Medical Systems (http://www.artinis.com). For more information
% on FieldTrip, see http://www.fieldtriptoolbox.org
%
% This work is licensed under a Creative Commons Attribution-ShareAlike 4.0
% International License. To view a copy of this license, visit
% http://creativecommons.org/licenses/by-sa/4.0/ or send a letter to
% Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.
%
% Creative Commons Attribution-ShareAlike 4.0 International License:
% -----------------------------------
% You are free to:
%
%     Share - copy and redistribute the material in any medium or format
%     Adapt - remix, transform, and build upon the material
%     for any purpose, even commercially.
%
%     The licensor cannot revoke these freedoms as long as you follow the
%     license terms.
%
% Under the following terms:
%
%     Attribution - You must give appropriate credit, provide a link to
%                    the license, and indicate if changes were made. You
%                    may do so in any reasonable manner, but not in any way
%                    that suggests the licensor endorses you or your use.
%
%     ShareAlike - If you remix, transform, or build upon the material,
%                   you must distribute your contributions under the same
%                   license as the original.
%
%     No additional restrictions - You may not apply legal terms or
%                                   technological measures that legally
%                                   restrict others from doing anything the
%                                   license permits.
%
% -----------------------------------
%
% This toolbox is not to be used for medical or clinical purposes.
%
% Copyright (c) 2021 by Artinis Medical Systems.
% Contact: askforinfo@artinis.com

if nargin > 2 && islogical(header)
  % header must be defined if begsample and endsample are used
  error('wrong type of header variable.')
end

if nargin == 1 || nargin == 2 && islogical(header) && ~header
  data = read_oxy5_header(filename);
  if isfield(data, 'opto')
    % ensure that the optode definition is according to the latest standards
    data.opto = ft_datatype_sens(data.opto);
  end
elseif nargin == 2  && islogical(header) && header
  data = read_oxy3_event(filename);
else % nargin > 1 && ~islogical(header)
  if nargin < 5
    chanindx = 1:header.nChans;
    if nargin < 4
      endsample = header.nSamples;
      if nargin < 3
        begsample = 1;
      end
    end
  end
  data = read_oxy5_data(filename, 0, header.nSamples); % sample subselection does not work yet

  if endsample > size(data, 2)
    warning('Cannot deliver all requested samples, nan''ing %d sample(s)', endsample-size(data, 2));
    data(:, end:endsample) = nan;
  end

  data = data(chanindx, begsample:endsample);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function writeGraph(nodes, edges, fileName)
%WRITEGRAPH Write a graph to an ascii file
%
%   writeGraph(NODES, EDGES, FILENAME)
%
%   Example
%     % create a basic graph and save it to a file
%     nodes = [10 10;20 10;10 20;20 20;27 15];
%     edges = [1 2;1 3;2 4;2 5;3 4;4 5];
%     writeGraph(nodes, edges, 'simpleGraph.txt');
%
%   See also
%     readGraph
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2014-01-20,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2014 INRA - Cepia Software Platform.

%% File creation and init

% extract "sizes" of the graph
nNodes = size(nodes, 1);
nEdges = size(edges, 1);
nDims  = size(nodes, 2);

% file opening
f = fopen(fileName, 'wt');
if f == -1
    error(['could not open file for writing: ' fileName]);
end

% write header
fprintf(f, '# graph\n');

%% write nodes info

fprintf(f, '# nodes\n');
fprintf(f, '%d %d\n', nNodes, nDims);
format = ['%g' repmat(' %g', 1, nDims-1) '\n'];
for iNode = 1:nNodes
    fprintf(f, format, nodes(iNode, :));
end

%% write edges info

fprintf(f, '# edges\n');
fprintf(f, '%d\n', nEdges);

format = '%d %d\n';
for iEdge = 1:nEdges
    fprintf(f, format, edges(iEdge, :));
end

% close file
fclose(f);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function patch_handle_array = draw_objsrc(obj_array, src_array, grid3d, withinterp, withpml)

chkarg(istypesizeof(obj_array, 'EMObject', [1 0]), ...
	'"obj_array" should be row vector with EMObject as elements.');
chkarg(istypesizeof(grid3d, 'Grid3d'), '"grid3d" should be instance of Grid3d.');
chkarg(istypesizeof(withinterp, 'logical'), '"withinterp" should be logical.');
chkarg(istypesizeof(withpml, 'logical'), '"withpml" should be logical.');

if withinterp
	lplot = grid3d.lplot(GT.prim, withinterp, withpml);
else
	lplot = grid3d.lvoxelbound(GT.prim, withpml);
end

for w = Axis.elems
	if length(lplot{w}) <= 3
		lplot{w} = linspace(lplot{w}(1), lplot{w}(end), 10);
	end
end

boxes_pml = EMObject.empty(0,0);
gray = [0.5 0.5 0.5];
if withpml
	for w = Axis.elems
		for s = Sign.elems
			if grid3d.Npml(w,s) > 0
				bound = grid3d.bound;
				bound(w, alter(s)) = grid3d.lpml(w,s);
				box = Box(bound);
				pml = Material('PML', gray, 1.0);
				box_pml = EMObject(box, pml);
				boxes_pml = [boxes_pml(1:end), box_pml];
			end
		end
	end
end

srcobj_array = EMObject.empty(0, length(src_array));
i = 0;
for src = src_array
	i = i+1;
	green = 'g';  % green
	srcmat = Material('Source', green, 1.0);
	srcobj_array(i) = EMObject(src.shape, srcmat);
end

patch_handle_array = [];
lplotobj = cell(1, Axis.count);  % indices
for obj = [obj_array, srcobj_array, boxes_pml]
	shape = obj.shape;
	color = obj.material.color;
% 	if ~isequal(color, 'none') && ~istypesizeof(shape, 'Domain')
	if ~isequal(color, 'none')
		for w = Axis.elems
			bn = shape.bound(w, Sign.n);
			bp = shape.bound(w, Sign.p);
			in = find(lplot{w} >= bn, 1, 'first');
			ip = find(lplot{w} <= bp, 1, 'last');

			% Find the indices for a box that is slighly larger than "shape".
			if in > 1
				in = in - 1;
			end
			if ip < length(lplot{w})
				ip = ip + 1;
			end
			lw = lplot{w}(in:ip);
			Nmin = 20;  % minimun number of sampling points
			if ~isempty(lw) && length(lw) < Nmin
				lw = linspace(lw(1), lw(end), Nmin);
			end
			lplotobj{w} = lw;
		end

		if ~isempty(lplotobj{Axis.x}) && ~isempty(lplotobj{Axis.y}) && ~isempty(lplotobj{Axis.z})  % shape in inside domain
			if ~istypesizeof(shape, 'ZeroVolShape')
				lsf = shape.lsf;
			else
				lsf = @(x,y,z) shape.lsf(x, y, z, true);
			end
			[X, Y, Z] = meshgrid(lplotobj{:});
			Level = lsf(X, Y, Z);
			hp = patch(isosurface(X, Y, Z, Level, -eps));  % -eps instead of 0 for ZeroVolumeShape
			isonormals(X, Y, Z, Level, hp)

			if istypesizeof(shape, 'Point')
				set(hp, 'FaceColor', color, 'EdgeColor', color, 'LineWidth', 3);
			elseif istypesizeof(shape, 'Point') || istypesizeof(shape, 'Plane')
				set(hp, 'FaceColor', color, 'EdgeColor', color, 'LineWidth', 1);
			else
				set(hp, 'FaceColor', color, 'EdgeColor', 'none');
			end
			if isequal(obj.material.name, 'PML') || isequal(obj.material.name, 'Source')
				alpha(hp, 0.2);
			end

			patch_handle_array = [patch_handle_array(1:end), hp];
		end
	end
end

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function S1 = my_min2(S1, sig, varargin)
% 2D radius filtering

sig = ceil(sig);

xs = repmat([-sig:sig], 2*sig+1, 1);
ys = xs';

rs = (xs.^2 + ys.^2).^.5;

xs = xs(rs<=sig);
ys = ys(rs<=sig);

dsnew = size(S1);
S1 = S1(:,:,:);

[Ly, Lx, NT] = size(S1);

Smax = S1;
for j = 1:numel(xs)
    yc = [1:Ly] + ys(j);
    ig = ~(yc<1 | yc>Ly);
    yc = yc(ig);

    xc = [1:Lx] + xs(j);
    ig = ~(xc<1 | xc>Lx);
    xc = xc(ig);

    Smax(yc - ys(j),xc-xs(j), :) = min(Smax(yc - ys(j), xc-xs(j),:), S1(yc, xc, :));
end

S1 = reshape(Smax, dsnew);

author: cortex-lab


text:
% [INPUT]
% ds = A structure representing the dataset.
% category = A string representing the type of validation to perform (optional, default='').
%
% [OUTPUT]
% ds = A structure containing the validated dataset.

function ds = validate_dataset(varargin)

    persistent categories;

    if (isempty(categories))
        categories = {
            'BubblesDetection' 'Comparison' 'Component' 'Connectedness' 'CrossEntropy' 'CrossQuantilogram' ...
            'CrossSectional' 'Default' 'Liquidity' 'RegimeSwitching' 'Spillover' 'TailDependence'
        };
    end

    persistent ip;

    if (isempty(ip))
        ip = inputParser();
        ip.addRequired('ds',@(x)validateattributes(x,{'struct'},{'nonempty'}));
        ip.addOptional('category','',@(x)any(validatestring(x,categories)));
    end

    ip.parse(varargin{:});

    ipr = ip.Results;
    ds = ipr.ds;
    category = ipr.category;

    nargoutchk(1,1);

    ds = validate_dataset_internal(ds,category,categories);

end

function ds = validate_dataset_internal(ds,category,categories)

    validate_field(ds,'File',{'char'},{'nonempty' 'size' [1 NaN]});
    validate_field(ds,'Version',{'char'},{'nonempty' 'size' [1 NaN]});
    creation_date = validate_field(ds,'CreationDate',{'double'},{'real' 'finite' '>=' datenum(2014,1,1) 'scalar'}); %#ok<DATNM>

    result = validate_field(ds,'Result',{'char'},{'optional' 'nonempty' 'size' [1 NaN]});
    result_date = validate_field(ds,'ResultDate',{'double'},{'optional' 'real' 'finite' '>=' datenum(2014,1,1) 'scalar'}); %#ok<DATNM>
    result_analysis = validate_field(ds,'ResultAnalysis',{'function_handle'},{'optional' 'scalar'});
    result_serial = validate_field(ds,'ResultSerial',{'char'},{'optional' 'nonempty' 'size' [1 NaN]});

    if (isempty(result))
        if (~isempty(result_date))
            error(['The dataset field ''ResultDate'' is invalid.' new_line() 'Expected value to be empty.']);
        end

        if (~isempty(result_analysis))
            error(['The dataset field ''ResultAnalysis'' is invalid.' new_line() 'Expected value to be empty.']);
        end

        if (~isempty(result_serial))
            error(['The dataset field ''ResultSerial'' is invalid.' new_line() 'Expected value to be empty.']);
        end
    else
        if (~ismember(result,categories))
            error(['The dataset field ''Result'' is invalid.' new_line() 'Expected value to be a string containing a valid result category.']);
        end

        if (isempty(result_date) || (result_date < creation_date))
            error(['The dataset field ''ResultDate'' is invalid.' new_line() 'Expected value to be a numeric date greater than or equal to ' datestr(creation_date,'dd/mm/yyyy') '.']); %#ok<DATST>
        end

        if (isempty(result_analysis))
            error(['The dataset field ''ResultAnalysis'' is invalid.' new_line() 'Expected value to be a function handle.']);
        end

        if (isempty(result_serial))
            error(['The dataset field ''ResultSerial'' is invalid.' new_line() 'Expected value to be a function handle.']);
        end
    end

    validate_field(ds,'TimeSeries',{'cellstr'},{'nonempty' 'size' [1 8]});

    n = validate_field(ds,'N',{'double'},{'real' 'finite' 'integer' '>=' 3 'scalar'});
    t = validate_field(ds,'T',{'double'},{'real' 'finite' 'integer' '>=' 252 'scalar'});

    validate_field(ds,'DatesNum',{'double'},{'real' 'finite' 'integer' '>' 0 'nonempty' 'size' [t 1]});
    validate_field(ds,'DatesStr',{'cellstr'},{'nonempty' 'size' [t 1]});
    validate_field(ds,'MonthlyTicks',{'logical'},{'scalar'});

    validate_field(ds,'IndexName',{'char'},{'nonempty' 'size' [1 NaN]});
    validate_field(ds,'FirmNames',{'cellstr'},{'nonempty' 'size' [1 n]});

    validate_field(ds,'Index',{'double'},{'real' 'finite' 'nonempty' 'size' [t 1]});
    validate_field(ds,'Returns',{'double'},{'real' 'nanfinite' 'nonempty' 'size' [t n]});

    validate_field(ds,'Prices',{'double'},{'optional' 'real' 'nanfinite' 'nonnegative' 'nonempty' 'size' [t n]});
    validate_field(ds,'Volumes',{'double'},{'optional' 'real' 'nanfinite' 'nonnegative' 'nonempty' 'size' [t n]});
    validate_field(ds,'Capitalizations',{'double'},{'optional' 'real' 'nanfinite' 'nonnegative' 'nonempty' 'size' [t n]});

    validate_field(ds,'RiskFreeRate',{'double'},{'optional' 'real' 'finite' 'nonempty' 'size' [t 1]});
    validate_field(ds,'CDS',{'double'},{'optional' 'real' 'nanfinite' 'nonnegative' 'nonempty' 'size' [t n]});

    validate_field(ds,'Assets',{'double'},{'optional' 'real' 'nanfinite' 'nonnegative' 'nonempty' 'size' [t n]});
    validate_field(ds,'Equity',{'double'},{'optional' 'real' 'nanfinite' 'nonempty' 'size' [t n]});
    validate_field(ds,'Liabilities',{'double'},{'optional' 'real' 'nanfinite' 'nonnegative' 'nonempty' 'size' [t n]});
    validate_field(ds,'SeparateAccounts',{'double'},{'optional' 'real' 'nanfinite' 'nonnegative' 'nonempty' 'size' [t n]});

    state_variables = validate_field(ds,'StateVariables',{'double'},{'optional' 'real' 'finite' 'nonempty' 'size' [t NaN]});
    validate_field(ds,'StateVariablesNames',{'cellstr'},{'optional' 'nonempty' 'size' [1 size(state_variables,2)]});

    groups = validate_field(ds,'Groups',{'double'},{'real' 'finite' 'integer' '>=' 0 'scalar'});

    if (groups == 0)
        validate_field(ds,'GroupDelimiters',{'double'},{'size' [0 0]});
        validate_field(ds,'GroupNames',{'double'},{'size' [0 0]});
        validate_field(ds,'GroupShortNames',{'double'},{'size' [0 0]});
    else
        validate_field(ds,'GroupDelimiters',{'double'},{'real' 'finite' 'integer' 'positive' 'increasing' 'nonempty' 'size' [(groups - 1) 1]});
        validate_field(ds,'GroupNames',{'cellstr'},{'nonempty' 'size' [groups 1]});
        validate_field(ds,'GroupShortNames',{'cellstr'},{'nonempty' 'size' [groups 1]});
    end

    crises = validate_field(ds,'Crises',{'double'},{'real' 'finite' 'integer' '>=' 0 'scalar'});
    crises_type = validate_field(ds,'CrisesType',{'char'},{'nonempty' 'size' [1 1]});

    if (~strcmp(crises_type,'E') && ~strcmp(crises_type,'R'))
        error(['The dataset field ''CrisesType'' is invalid.' new_line() 'Expected value to be either ''E'' or ''R''.']);
    end

    if (crises == 0)
        validate_field(ds,'CrisesDummy',{'double'},{'size' [0 0]});
        validate_field(ds,'CrisisDates',{'double'},{'size' [0 0]});
        validate_field(ds,'CrisisDummies',{'double'},{'size' [0 0]});
        validate_field(ds,'CrisisNames',{'double'},{'size' [0 0]});
    else
        validate_field(ds,'CrisesDummy',{'double'},{'real' 'finite' 'integer' '>=' 0 '<=' 1 'nonempty' 'size' [t 1]});

        if (strcmp(crises_type,'E'))
            validate_field(ds,'CrisisDates',{'double'},{'real' 'finite' 'integer' '>' 0 'nonempty' 'size' [crises 1]});
            validate_field(ds,'CrisisDummies',{'double'},{'size' [0 0]});
        else
            validate_field(ds,'CrisisDates',{'double'},{'real' 'finite' 'integer' '>' 0 'nonempty' 'size' [crises 2]});
            validate_field(ds,'CrisisDummies',{'double'},{'real' 'finite' 'integer' '>=' 0 '<=' 1 'nonempty' 'size' [t crises]});
        end

        validate_field(ds,'CrisisNames',{'cellstr'},{'nonempty' 'size' [crises 1]});
    end

    validate_field(ds,'Defaults',{'double'},{'real' 'nonempty' 'offset' 'size' [1 n]});
    validate_field(ds,'Insolvencies',{'double'},{'real' 'nonempty' 'offset' 'size' [1 n]});

    validate_field(ds,'SupportsComponent',{'logical'},{'scalar'});
    validate_field(ds,'SupportsConnectedness',{'logical'},{'scalar'});
    validate_field(ds,'SupportsCrossEntropy',{'logical'},{'scalar'});
    validate_field(ds,'SupportsCrossSectional',{'logical'},{'scalar'});
    validate_field(ds,'SupportsDefault',{'logical'},{'scalar'});
    validate_field(ds,'SupportsLiquidity',{'logical'},{'scalar'});
    validate_field(ds,'SupportsRegimeSwitching',{'logical'},{'scalar'});
    validate_field(ds,'SupportsSpillover',{'logical'},{'scalar'});
    validate_field(ds,'SupportsTailDependence',{'logical'},{'scalar'});
    validate_field(ds,'SupportsComparison',{'logical'},{'scalar'});

    if (~isempty(category))
        if (strcmp(category,'Comparison'))
            if (~ds.SupportsComparison)
                error('The dataset cannot be used to compare systemic risk measures.');
            end
        else
            supports = ['Supports' category];

            if (~ds.(supports))
                category_e = category;
                indices = find(isstrprop(category,'upper'));

                for i = 1:numel(indices)
                    index = indices(i);

                    if (index == 1)
                        category_e(i) = lower(category(index));
                    else
                        category_e = [category_e(1:index-1) '-' lower(category_e(index)) category_e(index+1:end)];
                        indices(i+1:end) = indices(i+1:end) + 1;
                    end
                end

                error(['The dataset cannot be used to calculate ''' category_e ''' measures.']);
            end
        end
    end

end

function value = validate_field(ds,field_name,field_type,field_validator)

    if (~isfield(ds,field_name))
        error(['The dataset does not contain the field ''' field_name '''.']);
    end

    value = ds.(field_name);
    value_iscellstr = (numel(field_type) == 1) && strcmp(field_type{1},'cellstr');
    value_isfinite = any(strcmp(field_validator,'nanfinite'));
    value_isoffset = any(strcmp(field_validator,'offset'));
    value_isoptional = strcmp(field_validator{1},'optional');

    if (value_isoptional)
        empty = false;

        try
            validateattributes(value,{'double'},{'size' [0 0]});
            empty = true;
        catch
        end

        if (empty)
            return;
        end

        field_validator = field_validator(2:end);
    end

    if (value_iscellstr)
        if (~iscellstr(value) || any(cellfun(@length,value) == 0) || any(cellfun(@(x)size(x,1),value) ~= 1)) %#ok<ISCLSTR>
            error(['The dataset field ''' field_name ''' is invalid.' new_line() 'Expected value to be a cell array of non-empty character vectors.']);
        end

        field_type{1} = 'cell';
    else
        if (value_isfinite)
            index = strcmp(field_validator,'nanfinite');
            field_validator(index) = [];
        end

        if (value_isoffset)
            index = strcmp(field_validator,'offset');
            field_validator(index) = [];
        end
    end

    try
        validateattributes(value,field_type,field_validator);
    catch e
        error(['The dataset field ''' field_name ''' is invalid.' new_line() strrep(e.message,'Expected input','Expected value')]);
    end

    if (value_isfinite)
        value_check = value;
        value_check(isnan(value_check)) = 0;

        try
            validateattributes(value_check,field_type,{'finite'});
        catch
            error(['The dataset field ''' field_name ''' is invalid.' new_line() 'Expected value to be finite.']);
        end
    end

    if (value_isoffset)
        value_check = value;
        value_check(isnan(value_check)) = 2;

        try
            validateattributes(value_check,field_type,{'>',1});
        catch
            error(['The dataset field ''' field_name ''' is invalid.' new_line() 'Expected value to contain all the finite values > 1.']);
        end
    end

end

author: TommasoBelluzzo

repo: SystemicRisk


text:
function answer = isOnParallelWorker()
    try
        answer = ~isempty(getCurrentTask());
    catch err
        if ~strcmp(err.identifier, 'MATLAB:UndefinedFunction')
            rethrow(err);
        end
        answer = false;
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [J, grad] = linearRegCostFunction(X, y, theta, lambda)
%LINEARREGCOSTFUNCTION Compute cost and gradient for regularized linear
%regression with multiple variables
%   [J, grad] = LINEARREGCOSTFUNCTION(X, y, theta, lambda) computes the
%   cost of using theta as the parameter for linear regression to fit the
%   data points in X and y. Returns the cost in J and the gradient in grad

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost and gradient of regularized linear
%               regression for a particular choice of theta.
%
%               You should set J to the cost and grad to the gradient.
%

J = ((1 / (2 * m)) * sum(power((X * theta - y), 2))) + ((lambda / (2 * m)) * sum(power(theta(2 : end), 2)));
G = (lambda / m) .* theta;
G(1) = 0;
grad = ((1/m) .* X' * (X*theta - y)) + G;

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

grad = grad(:);

end

author: UtkarshPathrabe


text:
function A_inter = my_interp(A, filter_coef)
% Interpolation function
% This function performs the interpolation of the input image,
% using the specified filter
%
% Input:    A: input image
%           filter_coef: filter coefficients
%           example:    filter_coef = [1 1]/2 (Bilinear Filter)
%                       filter_coef = [1 -5 20 20 -5 1]/32 (6-tap filter)
% Output:   A_inter:  interpolated image
%
% Athanasopoulos Dionysios
% Postgraduate Student
% Computer Engineering and Informatics Dept.
% University of Patras, Greece

if (length(size(A)) == 3)
% if the input image is an RGB image or an YUV image
% or any 3D matrix just interpolate the three components separetely
    for i=1:3
        A_inter(:,:,i) = my_interp(A(:,:,i),filter_coef);
    end

else

    [m,n] = size(A);
    A_ = []; A_inter = [];
    % columns interpolation
    A_col= filter2(filter_coef,A);
    for i=1:n
        A_ = [A_ A(:,i) A_col(:,i)];
    end
    A_(:,end) = [];
    % rows interpolation
    A_rows = filter2(filter_coef,A_')';
    for i=1:m
        A_inter = [A_inter; A_(i,:); A_rows(i,:)];
    end
    A_inter(end,:) = [];

end


text:
function y = conj( x )

%Disciplined convex/geometric programming information for CONJ:
%   CONJ(X) imposes no convexity restrictions on its arguments. However,
%   since CONJ(X)=X when X is real, it is only useful for complex
%   affine expressions.

y = cvx( x.size_, conj( x.basis_ ) );

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function C=discSinCosTrans(x,type)
%%DISCSINCOSTRANS Compute the 1D discrete sine or cosine transformation of
%              a given array. All 16 standard types of sine and cosine
%              transforms are available plus a unitary version of the type
%              II even cosine transform.
%
%INPUTS: x  An (N+1)X1 or 1X(N+1) real or complex vector whose discrete
%           sine or cosine transform is desired.
%     type  A string specifying the type of discrete sine or cosine
%           transform desired. If omitted or an empty matrix is passed, the
%           default is 'CIIe'. Possible values are:
%           'CIe' Perform a type I even discrete cosine transform. The
%                 item in element k of the result is
%                 x(0)+(-1)^k*X(N-1)+2*sum_{j=1}^{N-2}x(j)*cos(pi*j*k/(N-1))
%                 where all indexation is form 0.
%           'CIIe' Perform a type II even discrete cosine transform. The
%                 item in element k of the result is the sum from n=0 to
%                 N-1 of
%                 2*x(n)*cos(pi*k*(n+1/2)/N), where all indexation is from
%                 0.
%         'CIIeU' This is the same as type CIIe, except for scaling. The
%                 scaling is
%                 CU(1)=CII(1)/(2*sqrt(N));
%                 CU(2:end)=CII(2:end)*0.5*sqrt(2/N); where CIIe is the
%                 type 2 transform and CU is the unitary transform.
%         'CIIIe' Perform a type III even discrete cosine transform.
%          'CIVe' Perform a type IV even discrete cosine transform.
%           'CIo' Perform a type I odd discrete cosine transform.
%          'CIIo' Perform a type II odd discrete cosine transform.
%         'CIIIo' Perform a type III odd discrete cosine transform.
%          'CIVo' Perform a type IV odd discrete cosine transform.
%           'SIe' Perform a type I even discrete sine transform.
%          'SIIe' Perform a type II even discrete sine transform.
%         'SIIIe' Perform a type III even discrete sine transform.
%          'SIVe' Perform a type IV even discrete sine transform.
%         'SIIIe' Perform a type III even discrete sine transform.
%          'SIVe' Perform a type IV even discrete sine transform.
%           'SIo' Perform a type I odd discrete sine transform.
%          'SIIo' Perform a type II odd discrete sine transform.
%         'SIIIo' Perform a type III odd discrete sine transform.
%          'SIVo' Perform a type IV odd discrete sine transform.
%
%OUTPUTS: C The (N+1)X1 or 1X(N+1) real discrete cosine or sine
%           transformation of x.
%
%The algorithms are implemented using the generalized discrete Fourier
%transform (GDFT) based upon the relations of the sequences and the GDFT
%that are given in [1]. The definition of the unitary version of the type
%II transformation is taken from [2].
%
%The inverse of this function is invDiscSinCosTrans.
%
%REFERENCES:
%[1] S. Martucci, "Symmetric convolution and the discrete sine and cosine
%    transforms," IEEE Transactions on Signal Processing, vol. 42, no. 5,
%    pp. 1038-1051, May 1994.
%[2] J. Makhoul, "A fast cosine transform in one and two dimensions," IEEE
%    Transactions on Acoustics, Speech, and Signal Processing, vol. ASSP-
%    28, no. 1, pp. 27-34, Feb. 1980.
%
%December 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2)
    type='CIIe';
end

N=length(x);
if(size(x,1)==N)
    colVec=true;
else
    colVec=false;
end

x=x(:);
switch(type)
    case 'CIIeU'
        %Perform the CIIe transform
        C=genDFT([x;flip(x)],0,1/2);
        C=C(1:N);

        %Modify the weighting.
        C(1)=C(1)/(2*sqrt(N));
        C(2:end)=C(2:end)*0.5*sqrt(2/N);
    case 'CIe'
        C=fft([x;x((N-1):-1:2)]);
        C=C(1:N);
    case 'CIIe'
        C=genDFT([x;flip(x)],0,1/2);
        C=C(1:N);
    case 'CIIIe'
        C=genDFT([x;0;-flip(x(2:end))],1/2,0);
        C=C(1:N);
    case 'CIVe'
        C=genDFT([x;-flip(x)],1/2,1/2);
        C=C(1:N);
    case 'CIo'
        C=fft([x;flip(x(2:end))]);
        C=C(1:N);
    case 'CIIo'
        C=genDFT([x;flip(x(1:(end-1)))],0,1/2);
        C=C(1:N);
    case 'CIIIo'
        C=genDFT([x;-flip(x(2:end))],1/2,0);
        C=C(1:N);
    case 'CIVo'
        C=genDFT([x;0;-flip(x)],1/2,1/2);
        C=C(1:N);
    case 'SIe'
        C=1j*fft([0;x;0;-flip(x)]);
        C=C(2:(N+1));
    case 'SIIe'
        C=1j*genDFT([x;-flip(x)],0,1/2);
        C=C(2:(N+1));
    case 'SIIIe'
        C=1j*genDFT([0;x(1:end);flip(x(1:(end-1)))],1/2,0);
        C=C(1:N);
    case 'SIVe'
        C=1j*genDFT([x;flip(x)],1/2,1/2);
        C=C(1:N);
    case 'SIo'
        C=1j*fft([0;x;-flip(x)]);
        C=C(2:(N+1));
    case 'SIIo'
        C=1j*genDFT([x;0;-flip(x)],0,1/2);
        C=C(2:(N+1));
    case 'SIIIo'
        C=1j*genDFT([0;x;flip(x)],1/2,0);
        C=C(1:N);
    case 'SIVo'
        C=1j*genDFT([x;flip(x(1:(end-1)))],1/2,1/2);
        C=C(1:N);
    otherwise
        error('Discrete cosine transform type not suported.')
end

%Get rid of finite precision errors. We know that if the purely input is
%real or imaginary, then the output must be the same.
if(isreal(x))
	C=real(C);
elseif(isImag(x))
    C=1j*imag(C);
end

if(~colVec)
    C=C.';
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
%COLAMD_DEMO demo for colamd, column approx minimum degree ordering algorithm
%
% Example:
%   colamd_demo
%
% The following m-files and mexFunctions provide alternative sparse matrix
% ordering methods for MATLAB.  They are typically faster (sometimes much
% faster) and typically provide better orderings than their MATLAB counterparts:
%
%       colamd          a replacement for colmmd.
%
%                       Typical usage:  p = colamd (A) ;
%
%       symamd          a replacement for symmmd.  Based on colamd.
%
%                       Typical usage:  p = symamd (A) ;
%
% For a description of the methods used, see the colamd.c file.
%
% http://www.cise.ufl.edu/research/sparse/colamd/
%
% See also colamd, symamd

% Minor changes:  in MATLAB 7, symmmd and colmmd are flagged as "obsolete".
% This demo checks if they exist, so it should still work when they are removed.

%    Copyright 1998-2007, Timothy A. Davis, and Stefan Larimore
%    Developed in collaboration with J. Gilbert and E. Ng.

%-------------------------------------------------------------------------------
% Print the introduction, the help info, and compile the mexFunctions
%-------------------------------------------------------------------------------

fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, 'Colamd2/symamd2 demo.') ;
fprintf (1, '\n-----------------------------------------------------------\n') ;
help colamd_demo ;

fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, 'Colamd help information:') ;
fprintf (1, '\n-----------------------------------------------------------\n') ;
help colamd2 ;

fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, 'Symamd help information:') ;
fprintf (1, '\n-----------------------------------------------------------\n') ;
help symamd2 ;

%-------------------------------------------------------------------------------
% Solving Ax=b
%-------------------------------------------------------------------------------

n = 100 ;
fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, 'Solving Ax=b for a small %d-by-%d random matrix:', n, n) ;
fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, '\nNote: Random sparse matrices are AWFUL test cases.\n') ;
fprintf (1, 'They''re just easy to generate in a demo.\n') ;

% set up the system

rand ('state', 0) ;
randn ('state', 0) ;
spparms ('default') ;
A = sprandn (n, n, 5/n) + speye (n) ;
b = (1:n)' ;

fprintf (1, '\n\nSolving via lu (PAQ = LU), where Q is from colamd2:\n') ;
q = colamd2 (A) ;
I = speye (n) ;
Q = I (:, q) ;
[L,U,P] = lu (A*Q) ;
fl = luflops (L, U) ;
x = Q * (U \ (L \ (P * b))) ;
fprintf (1, '\nFlop count for [L,U,P] = lu (A*Q):          %d\n', fl) ;
fprintf (1, 'residual:                                     %e\n', norm (A*x-b));

try
    fprintf (1, '\n\nSolving via lu (PAQ = LU), where Q is from colmmd:\n') ;
    q = colmmd (A) ;
    I = speye (n) ;
    Q = I (:, q) ;
    [L,U,P] = lu (A*Q) ;
    fl = luflops (L, U) ;
    x = Q * (U \ (L \ (P * b))) ;
    fprintf (1, '\nFlop count for [L,U,P] = lu (A*Q):          %d\n', fl) ;
    fprintf (1, 'residual:                                     %e\n', ...
	norm (A*x-b)) ;
catch
    fprintf (1, 'colmmd is obsolete; test skipped\n') ;
end

fprintf (1, '\n\nSolving via lu (PA = LU), without regard for sparsity:\n') ;
[L,U,P] = lu (A) ;
fl = luflops (L, U) ;
x = U \ (L \ (P * b)) ;
fprintf (1, '\nFlop count for [L,U,P] = lu (A*Q):          %d\n', fl) ;
fprintf (1, 'residual:                                     %e\n', norm (A*x-b));

%-------------------------------------------------------------------------------
% Large demo for colamd2
%-------------------------------------------------------------------------------

fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, 'Large demo for colamd2 (symbolic analysis only):') ;
fprintf (1, '\n-----------------------------------------------------------\n') ;

rand ('state', 0) ;
randn ('state', 0) ;
spparms ('default') ;
n = 1000 ;
fprintf (1, 'Generating a random %d-by-%d sparse matrix.\n', n, n) ;
A = sprandn (n, n, 5/n) + speye (n) ;

fprintf (1, '\n\nUnordered matrix:\n') ;
lnz = symbfact (A, 'col') ;
fprintf (1, 'nz in Cholesky factors of A''A:            %d\n', sum (lnz)) ;
fprintf (1, 'flop count for Cholesky of A''A:           %d\n', sum (lnz.^2)) ;

tic ;
p = colamd2 (A) ;
t = toc ;
lnz = symbfact (A (:,p), 'col') ;
fprintf (1, '\n\nColamd run time:                          %f\n', t) ;
fprintf (1, 'colamd2 ordering quality: \n') ;
fprintf (1, 'nz in Cholesky factors of A(:,p)''A(:,p):  %d\n', sum (lnz)) ;
fprintf (1, 'flop count for Cholesky of A(:,p)''A(:,p): %d\n', sum (lnz.^2)) ;

try
    tic ;
    p = colmmd (A) ;
    t = toc ;
    lnz = symbfact (A (:,p), 'col') ;
    fprintf (1, '\n\nColmmd run time:                          %f\n', t) ;
    fprintf (1, 'colmmd ordering quality: \n') ;
    fprintf (1, 'nz in Cholesky factors of A(:,p)''A(:,p):  %d\n', sum (lnz)) ;
    fprintf (1, 'flop count for Cholesky of A(:,p)''A(:,p): %d\n', ...
	sum (lnz.^2)) ;
catch
    fprintf (1, 'colmmd is obsolete; test skipped\n') ;
end

%-------------------------------------------------------------------------------
% Large demo for symamd2
%-------------------------------------------------------------------------------

fprintf (1, '\n-----------------------------------------------------------\n') ;
fprintf (1, 'Large demo for symamd2 (symbolic analysis only):') ;
fprintf (1, '\n-----------------------------------------------------------\n') ;

fprintf (1, 'Generating a random symmetric %d-by-%d sparse matrix.\n', n, n) ;
A = A+A' ;

fprintf (1, '\n\nUnordered matrix:\n') ;
lnz = symbfact (A, 'sym') ;
fprintf (1, 'nz in Cholesky factors of A:       %d\n', sum (lnz)) ;
fprintf (1, 'flop count for Cholesky of A:      %d\n', sum (lnz.^2)) ;

tic ;
p = symamd2 (A) ;
t = toc ;
lnz = symbfact (A (p,p), 'sym') ;
fprintf (1, '\n\nSymamd run time:                   %f\n', t) ;
fprintf (1, 'symamd2 ordering quality: \n') ;
fprintf (1, 'nz in Cholesky factors of A(p,p):  %d\n', sum (lnz)) ;
fprintf (1, 'flop count for Cholesky of A(p,p): %d\n', sum (lnz.^2)) ;

try
    tic ;
    p = symmmd (A) ;
    t = toc ;
    lnz = symbfact (A (p,p), 'sym') ;
    fprintf (1, '\n\nSymmmd run time:                   %f\n', t) ;
    fprintf (1, 'symmmd ordering quality: \n') ;
    fprintf (1, 'nz in Cholesky factors of A(p,p):  %d\n', sum (lnz)) ;
    fprintf (1, 'flop count for Cholesky of A(p,p): %d\n', sum (lnz.^2)) ;
catch
    fprintf (1, 'symmmd is obsolete\n') ;
end


text:
function [f,the_regimes,endog,data,tex]=load_filters(m)

n=numel(m);

if n>1

    f=cell(size(m));
    the_regimes=f;
    endog=f;
    data=f;
    tex=f;

    for ii=1:n

        [f{ii},the_regimes{ii},endog{ii},data{ii},tex{ii}]=load_filters(m(ii));

    end

    return

end

if isa(m,'abstvar')

    [~,~,~,f]=filter(m);

    the_regimes=generic.describe_regimes(m.markov_chain_info);

    data=m.estim_.data;

    endog=m.endogenous;

    tex=[];

else

    f=filter(m);

    the_regimes=generic.describe_regimes(m.markov_chains);

    endog=m.observables.name;

    data=m.options.data;

    if isa(data,'ts')

        data=pages2struct(data);

    end

    tex=get(m,'tex');

end

end

author: jmaih


text:
function [net, sName, sName_l2norm] = addOneLoop_forMeanShiftGrouping(net, sName, loopIdx, GaussianBandwidth, randSampleRatio)

if ~exist('GaussianBandwidth', 'var')
    GaussianBandwidth = 0.1;
end

if ~exist('randSampleRatio', 'var')
    randSampleRatio = 0.2;
end

%%
pre_l2_norm_layer = sName;

lName = sprintf('loop%d_meanshift_S_is_XX', loopIdx);
net.addLayer(lName, ...
    meanshift_S_is_XX(), ... softmaxlog logistic
    {sName}, lName);
sName = lName;

lName = sprintf('loop%d_meanshift_G_is_Gaussian', loopIdx);
net.addLayer(lName, ...
    meanshift_G_is_Gaussian('delta', GaussianBandwidth), ...
    {sName}, lName);
G_layer = lName;
sName = lName;

lName = sprintf('loop%d_meanshift_d_is_sumG', loopIdx);
net.addLayer(lName, ...
    meanshift_d_is_sumG(), ...
    {sName}, lName);
sName = lName;

lName = sprintf('loop%d_meanshift_q_is_inv_d', loopIdx);
net.addLayer(lName, ...
    meanshift_q_is_inv_d(), ...
    {sName}, lName);
sName = lName;

lName = sprintf('loop%d_meanshift_P_is_G_diag_q', loopIdx);
net.addLayer(lName, ...
    meanshift_P_is_G_diag_q(), ...
    {G_layer, sName}, lName);
sName = lName;

lName = sprintf('loop%d_meanshift_Y_is_XP', loopIdx);
net.addLayer(lName, ...
    meanshift_Y_is_XP(), ...
    {pre_l2_norm_layer, sName}, lName);
sName = lName;

lName = sprintf('loop%d_meanshift_Y_l2norm', loopIdx);
net.addLayer(lName, L2normalization(), sName, lName) ;
sName = lName;
sName_l2norm = lName;

lName = sprintf('loop%d_meanshift_cosSim', loopIdx);
gt_name =  sprintf('gt_ins');
net.addLayer(lName, cosineSimilarity_randSample('randSampleRatio', randSampleRatio), {sName, gt_name}, lName) ;
sName = lName;

% % add regression loss
% obj_name = sprintf('loop%d_instSeg_reg', loopIdx);
% net.addLayer(obj_name, ...
%     InstanceSegRegLoss_randSample('loss', 'cosinesimilarityabsregloss', 'lastLayerName', sName), ... softmaxlog logistic
%     {sName, gt_name}, obj_name);
%
%
% % add max-margin loss
% obj_name = sprintf('loop%d_instSeg_MM', loopIdx);
% input_name = sName;
% net.addLayer(obj_name, ...
%     InstanceSegMMLoss_randSample('loss', 'cosinesimilaritymmloss', 'marginAlpha_', 0.1, 'adaptiveMM', false, 'lastLayerName', sName), ...
%     {input_name, gt_name}, obj_name)

author: aimerykong


text:
function [Population] = EnvironmentalSelection(Population,N)
% Environmental selection

%------------------------------- Copyright --------------------------------
% Copyright 2017-2018 Yiping Liu
% Please contact {yiping0liu@gmail.com} if you have any problem.
%--------------------------------------------------------------------------

    %% Non-dominated sorting
    N = min(N,length(Population));
    [FrontNo,MaxFNo] = NDSort(Population.objs,Population.cons,N);
    Next = FrontNo < MaxFNo;

    %% Calculate the crowding distance of each solution
    CrowdDis = CrowdingDistance(Population.objs,FrontNo);

    %% Select the solutions in the last front based on their crowding distances
    Last     = find(FrontNo==MaxFNo);
    [~,Rank] = sort(CrowdDis(Last),'descend');
    Next(Last(Rank(1:N-sum(Next)))) = true;

    %% Population for next generation
    Population = Population(Next);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


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

%%
clear; close all; clc;

%% Syntax
% |[F,V,regionIndSub]=hemiSphereRegionMesh(hemiSphereStruct);|

%% Description
% The |hemiSphereRegionMesh| function creates the faces (F), vertices (or
% nodes, V) and the region indices (regeionIndSub) for a hemi-sphere
% according to the input structure hemiSphereStruct. The latter defines the
% sphere radius, the number of refinement steps for the regions and the
% number of refinement steps for the mesh. For more information on the
% refinement see the |geoSphere| and |subTri| functions and associated demo
% files.
% A complete sphere is first represented as an icosahedron which is then
% refined (subtriangulated) hemiSphereStruct.nRefineRegions times (whereby
% for each iteration each triangle is subdevided into 4 triangles). This
% initial subdevision defines the element regions. The next refinement step
% defines the number of triangles for each region. The field
% hemiSphereStruct.nRefineMesh defines how many times each mesh region is
% iteratively subtriangulated.

%% Examples

%%
clear; close all; clc;

%%
% Plot settings
fontSize=25;
faceAlpha=1;
lineWidth=1;
markerSize=5;

%% Example:  Creating a hemisphere mesh using the |hemiSphereRegionMesh| function
% Defining hemi-sphere parameters
hemiSphereStruct.sphereRadius=1; %Sphere radius
hemiSphereStruct.nRefineRegions=1; %Number of refinement steps for regions
hemiSphereStruct.nRefineMesh=2; %Number of refinement steps for mesh

% Get hemi-sphere mesh
[F,V,regionInd]=hemiSphereRegionMesh(hemiSphereStruct);

%%
% Plotting results

%Creating a random color for the each mesh region
cmap=hsv(max(regionInd(:)));
cmap=cmap(randperm(size(cmap,1)),:); %scramble colors

hf=cFigure; hold on;
gtitle('Half dome showing regions with subtriangulated mesh',fontSize);
gpatch(F,V,regionInd);
colormap(cmap);
axisGeom(gca,fontSize);
camlight headlight;
drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
classdef P1Function < FeFunction

    properties (Access = public)
    end

    properties (Access = private)
        interpolation
    end

    properties (Access = private)
        meshCoarse
        meshFine
    end

    methods (Access = public)

        function obj = P1Function(cParams)
            obj.init(cParams);
            obj.createInterpolation();
        end

        function fxV = evaluate(obj, xV)
            obj.interpolation.computeShapeDeriv(xV);
            shapes = obj.interpolation.shape;
            nNode  = size(shapes,1);
            nGaus  = size(shapes,2);
            nF     = size(obj.fValues,2);
            nElem  = size(obj.mesh.connec,1);
            fxV = zeros(nF,nGaus,nElem);
            for iGaus = 1:nGaus
                for iNode = 1:nNode
                    node = obj.mesh.connec(:,iNode);
                    Ni = shapes(iNode,iGaus);
                    fi = obj.fValues(node,:);
                    f(:,1,:) = Ni*fi';
                    fxV(:,iGaus,:) = fxV(:,iGaus,:) + f;
                end
            end
        end

        function N = computeShapeFunctions(obj, quad)
%             obj.mesh.computeInverseJacobian(quad,obj.interpolation);
            xV = quad.posgp;
            obj.interpolation.computeShapeDeriv(xV);
            N = obj.interpolation.shape;
        end

        function dNdx  = computeCartesianDerivatives(obj,quad)
            switch obj.mesh.type
                case 'LINE'
                    invJ  = obj.mesh.computeInverseJacobian(quad,obj.interpolation);
                    nElem = obj.mesh.nelem;
                    nNode = obj.interpolation.nnode;
                    nDime = obj.mesh.ndim;
                    nGaus = quad.ngaus;
                    deriv  = obj.mesh.interpolation.deriv(1,:,:,:);
                    dShapes = deriv;
                    dN = zeros(nDime,nNode,nElem,nGaus);
                    for iGaus = 1:nGaus
                        invJg = invJ(:,1);
                        for iDim = 1:nDime
                            for iNode = 1:nNode
                                dShapeI(:,1) = squeeze(dShapes(1,iNode,iGaus,:));
                                dN(iDim,iNode,:,iGaus) = invJg.*dShapeI;
                            end
                        end
                    end
                    dNdx = dN;

                otherwise
                    nElem = size(obj.mesh.connec,1);
                    nNode = obj.interpolation.nnode;
                    nDime = obj.interpolation.ndime;
                    nGaus = quad.ngaus;
                    invJ  = obj.mesh.computeInverseJacobian(quad,obj.interpolation);
                    dShapeDx  = zeros(nDime,nNode,nElem,nGaus);
                    for igaus = 1:nGaus
                        dShapes = obj.interpolation.deriv(:,:,igaus);
                        for jDime = 1:nDime
                            invJ_JI   = invJ(:,jDime,:,igaus);
                            dShape_KJ = dShapes(jDime,:);
                            dSDx_KI   = bsxfun(@times, invJ_JI,dShape_KJ);
                            dShapeDx(:,:,:,igaus) = dShapeDx(:,:,:,igaus) + dSDx_KI;
                        end
                    end
                    dNdx = dShapeDx;
            end
        end

        function gradFun = computeGradient(obj, quad)
            dNdx = obj.computeCartesianDerivatives(quad);
            nDimf = obj.ndimf;
            nDims = size(dNdx, 1); % derivX, derivY (mesh-related?)
            nNode = size(dNdx, 2);
            nElem = size(dNdx, 3);
            nGaus = size(dNdx, 4);

            grad = zeros(nDims,nDimf, nElem, nGaus);
            for iGaus = 1:nGaus
                dNdx_g = dNdx(:,:,:,iGaus);
                for iDims = 1:nDims
                    for iNode = 1:nNode
                        dNdx_i = squeeze(dNdx_g(iDims, iNode,:));
                        nodes  = obj.mesh.connec(:,iNode);
                        f      = obj.fValues(nodes,:);
                        dfdx(1,:,:) = (dNdx_i.*f)';
                        grad(iDims,:,:,iGaus) = grad(iDims,:,:,iGaus) + dfdx;
                    end
                end
            end
            fVR = reshape(grad, [nDims*nDimf,nElem, nGaus]);
%             s.fValues    = permute(fVR, [1 3 2]);
            s.fValues    = fVR;
            s.mesh       = obj.mesh;
            s.quadrature = quad;
            gradFun = FGaussDiscontinuousFunction(s);
        end

        function symGradFun = computeSymmetricGradient(obj,quad)
            grad = obj.computeGradient(quad);
            nDimf = obj.ndimf;
            nDims = size(grad.fValues, 1)/nDimf;
            nGaus = size(grad.fValues, 2);
            nElem = size(grad.fValues, 3);

            gradReshp = reshape(grad.fValues, [nDims,nDimf,nGaus,nElem]);
            gradT = permute(gradReshp, [2 1 3 4]);
            symGrad = 0.5*(gradReshp + gradT);

            rshp = reshape(symGrad, [nDims*nDimf,nGaus,nElem]);
            s.fValues    = permute(rshp, [1 3 2]);
            s.quadrature = quad;
            s.mesh       = obj.mesh;
            symGradFun = FGaussDiscontinuousFunction(s);
        end

        function divF = computeDivergence(obj,q)
            dNdx = obj.computeCartesianDerivatives(q);
            fV = obj.fValues;
            nodes = obj.mesh.connec;
            nNode = obj.mesh.nnodeElem;
            nDim  = obj.mesh.ndim;
            divV = zeros(q.ngaus,obj.mesh.nelem);
            for igaus = 1:q.ngaus
                for kNode = 1:nNode
                    nodeK = nodes(:,kNode);
                    for rDim = 1:nDim
                        dNkr = squeeze(dNdx(rDim,kNode,:,igaus));
                        fkr = fV(nodeK,rDim);
                        int(1,:) = dNkr.*fkr;
                        divV(igaus,:) = divV(igaus,:) + int;
                    end
                end
            end
            s.quadrature = q;
            s.mesh       = obj.mesh;
            s.fValues(1,:,:) = divV;
            divF = FGaussDiscontinuousFunction(s);
        end

        function fdivF = computeFieldTimesDivergence(obj,q)
            fG  = obj.evaluate(q.posgp);
            dfG = obj.computeDivergence(q);
            fdivFG = bsxfun(@times,dfG.fValues,fG);
            s.quadrature = q;
            s.mesh       = obj.mesh;
            s.fValues    = fdivFG;
            fdivF = FGaussDiscontinuousFunction(s);
        end

        function fFine = refine(obj,m,mFine)
            fNodes  = obj.fValues;
            fEdges  = obj.computeFunctionInEdges(m, fNodes);
            fAll    = [fNodes;fEdges];
            s.mesh    = mFine;
            s.fValues = fAll;
            fFine = P1Function(s);
        end

        function dofConnec = computeDofConnectivity(obj)
            conne  = obj.mesh.connec;
            nDimf  = obj.ndimf;
            nNode  = size(conne, 2);
            nDofsE = nNode*nDimf;
            dofsElem  = zeros(nDofsE,size(conne,1));
            for iNode = 1:nNode
                for iUnkn = 1:nDimf
                    idofElem   = nDimf*(iNode - 1) + iUnkn;
                    globalNode = conne(:,iNode);
                    idofGlobal = nDimf*(globalNode - 1) + iUnkn;
                    dofsElem(idofElem,:) = idofGlobal;
                end
            end
            dofConnec = dofsElem;
        end

        function plot(obj) % 2D domains only

            switch obj.mesh.type
                case {'TRIANGLE','QUAD'}
                    x = obj.mesh.coord(:,1);
                    y = obj.mesh.coord(:,2);
                    figure()
                    for idim = 1:obj.ndimf
                        subplot(1,obj.ndimf,idim);
                        z = obj.fValues(:,idim);
                        a = trisurf(obj.mesh.connec,x,y,z);
                        view(0,90)
                        %             colorbar
                        shading interp
                        a.EdgeColor = [0 0 0];
                        title(['dim = ', num2str(idim)]);
                    end
                case 'LINE'
                    x = obj.mesh.coord(:,1);
                    y = obj.fValues;
                    figure()
                    plot(x,y)
            end
        end

        function plotArrowVector(obj)
            figure()
            a = obj.fValues;
            x = obj.mesh.coord(:,1);
            y = obj.mesh.coord(:,2);
            ax = a(:,1);
            ay = a(:,2);
            q = quiver(x,y,ax,ay);
            q.ShowArrowHead = 'off';
        end

        function print(obj, s)
%             s.mesh
            s.mesh = obj.mesh;
            s.fun = {obj};
            p = FunctionPrinter(s);
            p.print();
        end

        function [res, pformat] = getDataToPrint(obj)
            nNods = size(obj.fValues, 1);
            s.nDimf   = obj.ndimf;
            s.nData   = nNods;
            s.nGroup  = nNods;
            s.fValues = obj.getFormattedFValues();
            fps = FunctionPrintingSettings(s);
            [res, pformat] = fps.getDataToPrint();
        end

    end

    methods (Access = public, Static)

        function p1 = create(mesh, ndimf)
            s.fValues = zeros(mesh.nnodes, ndimf);
            s.mesh    = mesh;
            p1 = P1Function(s);
        end

        function fS = times(f1,f2)
            fS = f1.fValues.*f2.fValues;
            s.fValues = fS;
            s.mesh    = f1.mesh;
            fS = P1Function(s);
        end

        function fS = sum(f1,f2)
            fS = f1.fValues+f2.fValues;
            s.fValues = fS;
            s.mesh    = f1.mesh;
            fS = P1Function(s);
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.mesh    = cParams.mesh;
            obj.fValues = cParams.fValues;
            obj.ndimf   = size(cParams.fValues,2);
            obj.order   = 'LINEAR';
            obj.meshCoarse = cParams.mesh;
        end

        function createInterpolation(obj)
            m.type = obj.mesh.type;
            obj.interpolation = Interpolation.create(m,'LINEAR');
        end

        % Printing
        function fM = getFormattedFValues(obj)
            fM = obj.fValues;
        end

        function f = computeFunctionInEdges(obj,m,fNodes)
            s.edgeMesh = m.computeEdgeMesh();
            s.fNodes   = fNodes;
            eF         = EdgeFunctionInterpolator(s);
            f = eF.compute();
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function varargout = flipud(varargin) %#ok<STOUT>
%FLIPUD   Flip/reverse an UNBNDFUN object.
%   FLIPUD is not support for FUN objects on unbounded domains.

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

error('CHEBFUN:UNBNDFUN:flipud:noSupport', ...
    'FLIPUD does not support UNBNDFUN objects.');

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function c=ref_edgtii_1(f,g,a,M)
%REF_EDGTII_1   Reference Even Discrete Gabor transform type II by DGT
%   Usage  c=ref_edgt(f,g,a,M);
%
%   If a is even, then the input window must be odd-centered of length 2L.
%
%   If a is odd, then the input window must be even-centered of length 2L.
%

L=size(f,1);
W=size(f,2);

N=L/a;

clong=ref_gdgt([f;conj(flipud(f))],g,a,M,.5,0,floor(a/2));

c=clong(1:M*N,:);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function pred = knnpred2(Xtest,X,class,class_exp,K,dist_type,pret_type)

% prediction of new samples with calculated model
%
% pred = knnpred(Xtest,X,class,K,dist_type,pret_type)
%
% ------------ INPUT ---------------------------------------------------
% Xtest:        dataset to be predicted [n_test x p] n objects, p variables
% X:            training data matrix (n x p)
% class:        training class vector (n x 1)
% K:            number of neighbors
% dist_type:    'euclidean' Euclidean distance
%               'mahalanobis' Mahalanobis distance
%               'cityblock' City Block metric
%               'minkowski' Minkowski metric
%               'sm' Sokal-Michener
%               'jt' Jaccard Tanimoto
%               'gle' Gleason-Dice
%               'ct4' Consonni-Todeschini
%               'ac' Austin-Colwell
% pret_type:    'cent' cenering
%               'scal' variance scaling
%               'auto' for autoscaling (centering + variance scaling)
%               'rang' range scaling (0-1)
%               'fp'   fingerprints
%
% ------------ OUTPUT --------------------------------------------------
% pred is a structure conyaining
% class_pred    predicted class vector [n_test x 1]
% neighbors     list of k neighbors for each predicted sample [n_test x k]
%
% version 1.0 - september 2009
% Davide Ballabio
% Milano Chemometrics and QSAR Research Group
% www.disat.unimib.it/chm

% version 2.0 - February 2012
% Kamel Mansouri
% Milano Chemometrics and QSAR Research Group
% www.disat.unimib.it/chm

% data check
if length(class)~=size(X,1)
    disp('the class input should be for the training set')
    %class_tr=input('class tr');
    %class=evalin(WS,);
    %keyboard
end

[n,p] = size(Xtest);
[X_scal_train,param] = data_pretreatment(X,pret_type);
X_scal = test_pretreatment(Xtest,param);
Xd = [X_scal;X_scal_train];
% D = pdist(Xd,model.set.dist_type);
% D = squareform(D);
D = knn_calc_dist(X_scal_train,X_scal,dist_type,pret_type);
neighbors = zeros(n,K);
w=zeros(n,K);
class_exp(find(isnan(class_exp)))=class(find(isnan(class_exp)));
for i=1:n
    D_in = D(i,:);
    [d_tmp,n_tmp] = sort(D_in);
    neighbors(i,:) = n_tmp(1:K);
    d_neighbors = d_tmp(1:K);
    if d_neighbors(1)<1e-5 %&& d_neighbors(2)>d_neighbors(1)
        d_neighbors(1)=0;
    end
    if d_neighbors(1)==0 %&& isnan(class_exp(neighbors(i,1)))%&& d_neighbors(2)~=0
        class_calc(i) = knnclass2(class(neighbors(i,:)),d_neighbors,max(class),K);
        class_calc_weighted(i)=class(neighbors(i,1));
        w(i,1)=1;
    else
        %class(find(~isnan(class_exp)))=class_exp(find(~isnan(class_exp)));
        class_calc(i) = knnclass2(class_exp(neighbors(i,:)),d_neighbors,max(class_exp),K);
        [yc(i),class_calc_weighted(i),w(i,:)] = nnrcalcy2(class_exp(neighbors(i,:)),d_neighbors,K);
        class_calc_weighted(i)=round(class_calc_weighted(i));
    end

    dc(i,:)=d_neighbors;
end

pred.neighbors  = neighbors;
pred.class_pred = class_calc';
pred.class_pred_w = class_calc_weighted';
pred.w=w;
pred.D=D;
pred.dc=dc;

author: kmansouri

repo: OPERA

save_path: github-repos/MATLAB/kmansouri-OPERA


text:
function [fProbability_] = callback_LogLikelihoodA1A2Value(fAValues, caCatalogs, mControl, fBValue)
% function [fProbability] = callback_LogLikelihoodAValue(fAValue, caCatalogs, mControl)
% -------------------------------------------------------------------------------------
% Helper callback-function for calc_MaxLikelihoodA.m
%   Computes the negative log-likelihood sum of a given a- and fixed b-value for a
%   set of given catalogs
%
% Input parameters:
%   fAValue        a-value
%   caCatalogs      Cell array containing the set of catalogs
%   mControl        Controlmatrix containing informations about the single catalogs
%                   mControl(n,:) contains information about caCatalogs{n}
%                   Column 1: Starting time of catalog
%                   Column 2: Magnitude of completeness
%                   Column 3: Starting magnitude bin
%                   Column 4: Magnitude bin stepsize (must be 0.1)
%   fBValue         Fixed b-value
%
% Output parameters:
%   fProbability    Negative log-likelihood of the given a- and b-value for the set of given catalogs
%
% Danijel Schorlemmer
% July 5, 2002

% Init variable
global fProbability;
vProbabilities = [];

[nRow_, nColumn_] = size(mControl);
fTotalLength_ = mControl(nRow_,6)-mControl(1,1);

% Loop over all catalogs
for nCnt_ = 1:length(caCatalogs)
  % Extract catalog from cell array
  mCatalog_ = caCatalogs{nCnt_};
  % Determine maximum magnitude of catalog
  fMaxMag_ = max(mCatalog_(:,6));
  % Set up vector of available magnitude bins
  vCnt_ = (mControl(nCnt_,2):mControl(nCnt_,4):(fMaxMag_+mControl(nCnt_,4)))'; % Add one more magnitude bin for later use of diff()
  % Compute lengths of periods and ajust thea-value
  fTimeLength_ = mControl(nCnt_,6) - mControl(nCnt_,1);
  fTimeRatio_ = fTimeLength_/fTotalLength_;
  % Compute the cumulative FMD
  if mControl(nCnt_,5) == 1;  % this is activity rate 1
      fA_ = fAValues(1) + log10(fTimeRatio_);
      vNumber_ = 10.^(fA_ - (fBValue * vCnt_));
  elseif  mControl(nCnt_,5) == 2;  % this is activity rate 1
      fA_ = fAValues(2) + log10(fTimeRatio_);
      vNumber_ = 10.^(fA_ - (fBValue * vCnt_));
  end
  % Determine the number of events in each magnitude bin
  mPredictionFMD_ = -diff(vNumber_);
  % Create the FMD for the period of observation
  vObservedFMD_ = histogram(mCatalog_(:,6), mControl(nCnt_,2):mControl(nCnt_,4):fMaxMag_);
  % Calculate the likelihoods for both of the models
  vProb_ = calc_log10poisspdf(vObservedFMD_', mPredictionFMD_);
  % Return the values (multiply by -1 to return the lowest value for the highest probability
  vProbabilities = [vProbabilities; sum(vProb_)];
end
% Sum the probabilities for all given catalogs
fProbability_ = (-1) * sum(vProbabilities);

fProbability = fProbability_;

author: CelsoReyes


text:
function[]=makefigs_morsebox
%MAKEFIGS_MORSEBOX  Makes a sample figure for MORSEBOX.

ga1=(1/3:.1:11);
be1=(1/3:.1:10);

[ga,be]=meshgrid(ga1,be1);
[fm,fe,fi,cf] = morsefreq(ga,be);
a=morsebox(ga,be);

figure
contourf(ga1,be1,a,(.5:.01:.6))
hold on, contour(ga1,be1,a,(.500:.002:.51),'k:')
colormap gray,flipmap,
axis([1 10 1 10])
xtick(1:10),ytick(1:10)
ax=gca;
hc=colorbar;

contour(ga1,be1,(fm-fe)./(2*pi),[ 0 0],'k','linewidth',2)
contour(ga1,be1,(fm-fi)./(2*pi),[ 0 0],'k','linewidth',2)
contour(ga1,be1,cf./(2*pi),[ 0 0],'k','linewidth',2)
caxis([.5 .6])
vlines(3,'k--')
plot(ga1,12./ga1,'k')

title('Morse Wavelet Area and Transitions')
xlabel('Gamma Parameter')
ylabel('Beta Parameter')
plot([1+sqrt(-1)*0 10+sqrt(-1)*9/2],'k','linewidth',3)
plot([1+sqrt(-1)*0 10+sqrt(-1)*9/2],'w--','linewidth',2)

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function labeled_output = PTKWatershedMeyerFromStartingPointsMatlab(image, starting_labels)
    % PTKWatershedMeyerFromStartingPointsMatlab. Implementation of Meyer flooding algorithm.
    %
    %     This is a Matlab-only implementation of the mex function
    %     PTKWatershedMeyerFromStartingPoints. Use the mex function is possible
    %     as this is faster.
    %
    %     Inputs
    %     ------
    %
    %     image = 16-bit ingeter image (int16). The watershed regions grow according to the minima of these points
    %
    %     starting_labels - 8-bit integer (int8). Labels of starting points for the watershed
    %
    %     labeled_output - 8-bit integer (int8). Labels of the image assigned to
    %                       watershed regions. Watershed points are given the label -2
    %
    %     The watershed starts from the positive-valued labels in
    %     starting_labels and grows out into the zero-valued points, with the
    %     result returned in labeled_output. Regions starting from points with
    %     the same label can merge together. Negative labels are treated as
    %     fixed barriers. The do not grow and other regions cannot grow into
    %     them.
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    disp('PTKWatershedMeyerFromStartingPointsMatlab');
    warning('PTKWatershedMeyerFromStartingPointsMatlab:MatlabVersion', 'This is a slow Matlab implementation of the mex function PTKWatershedMeyerFromStartingPoints. For improved performance, use PTKWatershedMeyerFromStartingPoints instead.');

    % Check inputs
    if (nargin ~= 2)
        error('Two inputs are required: the image and a label matrix of the staring points.');
    end

    if (nargout > 1)
         error('PTKWatershedMeyerFromStartingPointsMatlab produces one output but you have requested more.');
    end

    % Get the input images
    intensity_matrix = image;
    starting_indices = starting_labels;

    if ~isequal(size(image), size(starting_labels));
        error('The two input matrices must be of the same dimensions.');
    end

    intensity_data = intensity_matrix;
    startingpoints_data = starting_indices;

    % Our set of points uses a custom comparison function so it is automatically sorted by image intensity
    % but is guaranteed uniqueness in the voxel indices
    points_to_do = [];
    intensity_store = [];

    image_size = size(image);
    size_i = image_size(1);
    size_j = image_size(2);
    size_k = image_size(3);
    number_of_points = size_i*size_j*size_k;

    % Linear index offsets to nearest neighbours
    offsets = [1, -1, size_i, -size_i, size_i*size_j, -size_i*size_j];

    iteration_number = 0;
    max_iterations = 1000000000;

    % Initialise the output data
    output_data = startingpoints_data;

    % Populate the initial set of points
    initial_points = find(startingpoints_data > 0);
    for index = 1 : length(initial_points)
        point_index = initial_points(index);

        neighbours = point_index + offsets;
        neighbours = neighbours((neighbours > 0) & (neighbours <= number_of_points) & (output_data(neighbours) == 0));
        points_to_do = [points_to_do, neighbours];
        intensity_store = [intensity_store, intensity_data(neighbours)];
    end

    total_points = sum(starting_labels(:) == 0);

    % Iterate over remaining points
    while ~isempty(points_to_do)

        if (mod(iteration_number, 10000) == 0)
            points_left = sum(output_data(:) == 0);
            percentage_done = 100*(total_points - points_left)/total_points;
            disp(['Percentage done:' num2str(percentage_done)]);
        end

        [~, min_index] = min(intensity_store);
        point_index = points_to_do(min_index);
        points_to_do(min_index) = [];
        intensity_store(min_index) = [];

        label = output_data(point_index);

        % The point may already have been set
        if (label == 0)
            % Check nearest neighbours of this point
            neighbours = point_index + offsets;
            neighbours = neighbours((neighbours > 0) & (neighbours <= number_of_points));

            % Find label values of neighbours
            neighbour_labels = output_data(neighbours);

            % Ignore unset, watershed and boundary points
            neighbour_labels = unique(neighbour_labels(neighbour_labels > 0));

            % No labeled neighbour - this case should not be possible and
            % indidates a program error
            if isempty(neighbour_labels)
                error('No neighbouring point found - this case should never occur.');
            end

            if (numel(neighbour_labels) > 1)
                % More than one type of neighbouring label - mark as watershed point
                output_data(point_index) = -2;
            else

                % One neighbouring label found - mark this point
                output_data(point_index) = neighbour_labels;

                % Add neighbours to points-to-do
                neighbours = neighbours(output_data(neighbours) == 0);
                points_to_do = [points_to_do, neighbours];
                intensity_store = [intensity_store, intensity_data(neighbours)];
            end
        end

        iteration_number = iteration_number + 1;
        if (iteration_number > max_iterations)
            error('Error: Max Iteration number exceeded');
        end
    end
    labeled_output = output_data;
    disp(' - Completed PTKWatershedMeyerFromStartingPointsMatlab');
end

author: tomdoel

repo: pulmonarytoolkit


text:
function model = VGG16_for_Fast_RCNN_VOC0712(model)
% VGG 16layers (only finetuned from conv3_1)

model.solver_def_file        = fullfile(pwd, 'models', 'fast_rcnn_prototxts', 'vgg_16layers_conv3_1', 'solver_30k60k.prototxt');
model.test_net_def_file      = fullfile(pwd, 'models', 'fast_rcnn_prototxts', 'vgg_16layers_conv3_1', 'test.prototxt');

model.net_file               = fullfile(pwd, 'models', 'pre_trained_models', 'vgg_16layers', 'vgg16.caffemodel');
model.mean_image             = fullfile(pwd, 'models', 'pre_trained_models', 'vgg_16layers', 'mean_image');

end

author: ShaoqingRen


text:
function th=thresh(X)
count=imhist(X);
[M, N]=size(X);

th=M/5;
end

author: HuangCongQing


text:
function [ value, seed ] = uniform_01_sample ( seed )

%*****************************************************************************80
%
%% UNIFORM_01_SAMPLE is a portable random number generator.
%
%  Formula:
%
%    SEED = SEED * (7**5) mod ( 2**31 - 1 )
%    UNIFORM_01_SAMPLE = SEED * / ( 2**31 - 1 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 August 2007
%
%  Parameters:
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real VALUE, a random value between 0 and 1.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
%  Local parameters:
%
%    IA = 7**5
%    IB = 2**15
%    IB16 = 2**16
%    IP = 2**31-1
%
  ia = 16807;
  ib15 = 32768;
  ib16 = 65536;
  ip = 2147483647;
%
%  Don't let SEED be 0 or IP
%
  if ( seed == 0 | seed == ip )
    seed = floor ( ip / 2 );
  end
%
%  Get the 15 high order bits of SEED.
%
  ixhi = floor ( seed / ib16 );
%
%  Get the 16 low bits of SEED and form the low product.
%
  loxa = ( seed - ixhi * ib16 ) * ia;
%
%  Get the 15 high order bits of the low product.
%
  leftlo = floor ( loxa / ib16 );
%
%  Form the 31 highest bits of the full product.
%
  iprhi = ixhi * ia + leftlo;
%
%  Get overflow past the 31st bit of full product.
%
  k = floor ( iprhi / ib15 );
%
%  Assemble all the parts and presubtract IP.  The parentheses are essential.
%
  seed = ( ( ( loxa - leftlo * ib16 ) - ip ) + ( iprhi - k * ib15 ) * ib16 ) + k;
%
%  Add IP back in if necessary.
%
  if ( seed < 0 )
    seed = seed + ip;
  end
%
%  Multiply by 1 / (2**31-1).
%
  value = seed * 4.656612875E-10;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function data = getSubSet(data, id_subset)
% SYNTAX
%   data = getSubSet(data, id_subset)
%
% DESCRIPTION
%   return the array with the given ids

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

    data = data(id_subset);
end

author: goGPS-Project


text:
function [out1,out2,out3] = ginput4(arg1)

[out1,out2,out3] = ginput(arg1);

return;

%GINPUT Graphical input from mouse.
%   [X,Y] = GINPUT(N) gets N points from the current axes and returns
%   the X- and Y-coordinates in length N vectors X and Y.  The cursor
%   can be positioned using a mouse (or by using the Arrow Keys on some
%   systems).  Data points are entered by pressing a mouse button
%   or any key on the keyboard except carriage return, which terminates
%   the input before N points are entered.
%
%   [X,Y] = GINPUT gathers an unlimited number of points until the
%   return key is pressed.
%
%   [X,Y,BUTTON] = GINPUT(N) returns a third result, BUTTON, that
%   contains a vector of integers specifying which mouse button was
%   used (1,2,3 from left) or ASCII numbers if a key on the keyboard
%   was used.
%
%   Examples:
%       [x,y] = ginput;
%
%       [x,y] = ginput(5);
%
%       [x, y, button] = ginput(1);
%
%   See also GTEXT, UIRESTORE, UISUSPEND, WAITFORBUTTONPRESS.

%   Copyright 1984-2006 The MathWorks, Inc.
%   $Revision: 5.32.4.9 $  $Date: 2006/12/20 07:19:10 $

P = NaN*ones(16,16);
P(1:15,1:15) = 2*ones(15,15);
P(2:14,2:14) = ones(13,13);
P(3:13,3:13) = NaN*ones(11,11);
P(6:10,6:10) = 2*ones(5,5);
P(7:9,7:9) = 1*ones(3,3);

out1 = []; out2 = []; out3 = []; y = [];
c = computer;
if ~strcmp(c(1:2),'PC')
   tp = get(0,'TerminalProtocol');
else
   tp = 'micro';
end

if ~strcmp(tp,'none') && ~strcmp(tp,'x') && ~strcmp(tp,'micro'),
   if nargout == 1,
      if nargin == 1,
         out1 = trmginput(arg1);
      else
         out1 = trmginput;
      end
   elseif nargout == 2 || nargout == 0,
      if nargin == 1,
         [out1,out2] = trmginput(arg1);
      else
         [out1,out2] = trmginput;
      end
      if  nargout == 0
         out1 = [ out1 out2 ];
      end
   elseif nargout == 3,
      if nargin == 1,
         [out1,out2,out3] = trmginput(arg1);
      else
         [out1,out2,out3] = trmginput;
      end
   end
else

   fig = gcf;
   figure(gcf);

   if nargin == 0
      how_many = -1;
      b = [];
   else
      how_many = arg1;
      b = [];
      if  ischar(how_many) ...
            || size(how_many,1) ~= 1 || size(how_many,2) ~= 1 ...
            || ~(fix(how_many) == how_many) ...
            || how_many < 0
         error('MATLAB:ginput:NeedPositiveInt', 'Requires a positive integer.')
      end
      if how_many == 0
         ptr_fig = 0;
         while(ptr_fig ~= fig)
            ptr_fig = get(0,'PointerWindow');
         end
         scrn_pt = get(0,'PointerLocation');
         loc = get(fig,'Position');
         pt = [scrn_pt(1) - loc(1), scrn_pt(2) - loc(2)];
         out1 = pt(1); y = pt(2);
      elseif how_many < 0
         error('MATLAB:ginput:InvalidArgument', 'Argument must be a positive integer.')
      end
   end

   % Suspend figure functions
   state = uisuspend(fig);

   toolbar = findobj(allchild(fig),'flat','Type','uitoolbar');
   if ~isempty(toolbar)
        ptButtons = [uigettool(toolbar,'Plottools.PlottoolsOff'), ...
                     uigettool(toolbar,'Plottools.PlottoolsOn')];
        ptState = get (ptButtons,'Enable');
        set (ptButtons,'Enable','off');
   end

   %set(fig,'pointer','fullcrosshair');
   set(fig,'Pointer','custom','PointerShapeCData',P,'PointerShapeHotSpot',[8,8]);

   fig_units = get(fig,'units');
   char = 0;

   % We need to pump the event queue on unix
   % before calling WAITFORBUTTONPRESS
   drawnow

   while how_many ~= 0
      % Use no-side effect WAITFORBUTTONPRESS
      waserr = 0;
      try
	keydown = wfbp;
      catch
	waserr = 1;
      end
      if(waserr == 1)
         if(ishandle(fig))
            set(fig,'units',fig_units);
	    uirestore(state);
            error('MATLAB:ginput:Interrupted', 'Interrupted');
         else
            error('MATLAB:ginput:FigureDeletionPause', 'Interrupted by figure deletion');
         end
      end

      ptr_fig = get(0,'CurrentFigure');
      if(ptr_fig == fig)
         if keydown
            char = get(fig, 'CurrentCharacter');
            button = abs(get(fig, 'CurrentCharacter'));
            scrn_pt = get(0, 'PointerLocation');
            set(fig,'units','pixels')
            loc = get(fig, 'Position');
            % We need to compensate for an off-by-one error:
            pt = [scrn_pt(1) - loc(1) + 1, scrn_pt(2) - loc(2) + 1];
            set(fig,'CurrentPoint',pt);
         else
            button = get(fig, 'SelectionType');
            if strcmp(button,'open')
               button = 1;
            elseif strcmp(button,'normal')
               button = 1;
            elseif strcmp(button,'extend')
               button = 2;
            elseif strcmp(button,'alt')
               button = 3;
            else
               error('MATLAB:ginput:InvalidSelection', 'Invalid mouse selection.')
            end
         end
         pt = get(gca, 'CurrentPoint');

         how_many = how_many - 1;

         if(char == 13) % & how_many ~= 0)
            % if the return key was pressed, char will == 13,
            % and that's our signal to break out of here whether
            % or not we have collected all the requested data
            % points.
            % If this was an early breakout, don't include
            % the <Return> key info in the return arrays.
            % We will no longer count it if it's the last input.
            break;
         end

         out1 = [out1;pt(1,1)];
         y = [y;pt(1,2)];
         b = [b;button];
      end
   end

   uirestore(state);
   if ~isempty(toolbar) && ~isempty(ptButtons)
        set (ptButtons(1),'Enable',ptState{1});
        set (ptButtons(2),'Enable',ptState{2});
   end
   set(fig,'units',fig_units);

   if nargout > 1
      out2 = y;
      if nargout > 2
         out3 = b;
      end
   else
      out1 = [out1 y];
   end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function key = wfbp
%WFBP   Replacement for WAITFORBUTTONPRESS that has no side effects.

fig = gcf;
current_char = [];

% Now wait for that buttonpress, and check for error conditions
waserr = 0;
try
  h=findall(fig,'type','uimenu','accel','C');   % Disabling ^C for edit menu so the only ^C is for
  set(h,'accel','');                            % interrupting the function.
  keydown = waitforbuttonpress;
  current_char = double(get(fig,'CurrentCharacter')); % Capturing the character.
  if~isempty(current_char) && (keydown == 1)           % If the character was generated by the
	  if(current_char == 3)                       % current keypress AND is ^C, set 'waserr'to 1
		  waserr = 1;                             % so that it errors out.
	  end
  end

  set(h,'accel','C');                                 % Set back the accelerator for edit menu.
catch
  waserr = 1;
end
drawnow;
if(waserr == 1)
   set(h,'accel','C');                                % Set back the accelerator if it errored out.
   error('MATLAB:ginput:Interrupted', 'Interrupted');
end

if nargout>0, key = keydown; end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

repo: ekfmonoslam