text:
function [seg,glo]=snrseg(s,r,fs,m,tf)
%SNRSEG Measure segmental and global SNR [SEG,GLO]=(S,R,FS,M,TF)
%
%Usage: (1) seg=snrseg(s,r,fs);                  % s & r are noisy and clean signal
%       (2) seg=snrseg(s,r,fs,'wz');             % no VAD or inerpolation used ['Vq' is default]
%       (3) [seg,snr]=snrseg(s,r,fs,'Vq',0.03);  % 30 ms frames
%
% Inputs:    s  test signal
%            r  reference signal
%           fs  sample frequency (Hz)
%            m  mode [default = 'V']
%                 w = No VAD - use whole file
%                 v = use sohn VAD to discard silent portions
%                 V = use P.56-based VAD to discard silent portions [default]
%                 a = A-weight the signals
%                 b = weight signals by BS-468
%                 q = use quadratic interpolation to remove delays +- 1 sample
%                 z = do not do any alignment
%                 p = plot results
%           tf  frame increment [0.01]
%
% Outputs: seg = Segmental SNR in dB
%          glo = Global SNR in dB (typically 7 dB greater than SNR-seg)
%
% This function compares a noisy signal, S, with a clean reference, R, and
% computes the segemntal signal-to-noise ratio (SNR) in dB. The signals,
% which must be of the same length, are split into non-overlapping frames
% of length TF (default 10 ms) and the SNR of each frame in dB is calculated.
% The segmental SNR is the average of these values, i.e.
%         SEG = mean(10*log10(sum(Ri^2)/sum((Si-Ri)^2))
% where the mean is over frames and the sum runs over one particular frame.
% Two optional modifications can be made to this basic formula:
%
%    (a) Frames are excluded if there is no significant energy in the R
%        signal. The idea is to limit the calculation to frames in which
%        speech is active. By default, the voicebox function "activlev" is
%        used to detect the inactive frames (the 'V' mode option).
%
%    (b) In each frame independently, the reference signal is shifted by up
%        to +- 1 sample to find the alignment than minimizes the noise
%        component (S-R)^2. This shifting accounts for small misalignments
%        and/or sample frequency differences between the two signals. For
%        larger shifts, you can use the voicebox function "sigalign".
%        Accurate alignemnt is especially important at high SNR values.
%
% If no M argument is specified, both these modifications will be applied;
% this is equivalent to specifying M='Vq'.

% Bugs/suggestions
% (1) Optionally restrict the bandwidth to the smaller of the two
%     bandwidths either with an extra parameter or automatically determined

%      Copyright (C) Mike Brookes 2011
%      Version: $Id: snrseg.m,v 1.4 2011/07/08 08:12:34 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<4 || ~ischar(m)
    m='Vq';
end
if nargin<5 || ~numel(tf)
    tf=0.01; % default frame length is 10 ms
end
snmax=100;  % clipping limit for SNR

% filter the input signals if required

if any(m=='a')  % A-weighting
    [b,a]=stdspectrum(2,'z',fs);
    s=filter(b,a,s);
    r=filter(b,a,r);
elseif any(m=='b') %  BS-468 weighting
    [b,a]=stdspectrum(8,'z',fs);
    s=filter(b,a,s);
    r=filter(b,a,r);
end

mq=~any(m=='z');
nr=min(length(r), length(s));
kf=round(tf*fs); % length of frame in samples
ifr=kf+mq:kf:nr-mq; % ending sample of each frame
ifl=ifr(end);
nf=numel(ifr);
rf=sum(reshape(r(mq+1:ifl).^2,kf,nf),1);
ef=sum(reshape((s(mq+1:ifl)-r(mq+1:ifl)).^2,kf,nf),1);
if mq
    efm=sum(reshape((s(3:ifl+1)-r(2:ifl)).^2,kf,nf),1);
    efp=sum(reshape((s(1:ifl-1)-r(2:ifl)).^2,kf,nf),1);
    efa=0.5*(efp+efm)-ef;
    efb=0.5*(efp-efm);
    efmk=(abs(efb)<2*efa) & (efa>0); % mask for frames with a valid minimum
    if any(efmk)
        ef(efmk)=ef(efmk)-0.25*efb(efmk).^2./efa(efmk);
    end
    ef=min(min(ef,efm),efp);
end

em=ef==0; % mask for zero noise frames
rm=rf==0; % mask for zero reference frames
snf=10*log10((rf+rm)./(ef+em));
snf(rm)=-snmax;
snf(em)=snmax;

% select the frames to include

if any(m=='w')
    vf=ones(1,nf); % include all frames
elseif any(m=='v');
    vs=vadsohn(r,fs,'na');
    nvs=length(vs);
    [vss,vix]=sort([ifr'; vs(:,2)]);
    vjx=zeros(nvs+nf,5);
    vjx(vix,1)=(1:nvs+nf)'; % sorted position
    vjx(1:nf,2)=vjx(1:nf,1)-(1:nf)'; % prev VAD frame end (or 0 or nvs+1 if none)
    vjx(nf+1:end,2)=vjx(nf+1:end,1)-(1:nvs)'; % prev snr frame end (or 0 or nvs+1 if none)
    dvs=[vss(1)-mq; vss(2:end)-vss(1:end-1)];  % number of samples from previous frame boundary
    vjx(:,3)=dvs(vjx(:,1)); % number of samples from previous frame boundary
    vjx(1:nf,4)=vs(min(1+vjx(1:nf,2),nvs),3); % VAD result for samples between prev frame boundary and this one
    vjx(nf+1:end,4)=vs(:,3); % VAD result for samples between prev frame boundary and this one
    vjx(1:nf,5)=1:nf; % SNR frame to accumulte into
    vjx(vjx(nf+1:end,2)>=nf,3)=0;  % zap any VAD frame beyond the last snr fram
    vjx(nf+1:end,5)=min(vjx(nf+1:end,2)+1,nf); % SNR frame to accumulate into
    vf=full(sparse(1,vjx(:,5),vjx(:,3).*vjx(:,4),1,nf))>kf/2; % accumulate into SNR frames and compare with threshold
else  % default is 'V'
    [lev,af,fso,vad]=activlev(r,fs);    % do VAD on reference signal
    vf=sum(reshape(vad(mq+1:ifl),kf,nf),1)>kf/2; % find frames that are mostly active
end
seg=mean(snf(vf));
glo=10*log10(sum(rf(vf))/sum(ef(vf)));

if ~nargout || any (m=='p')
    subplot(311);
    plot((1:length(s))/fs,s);
    ylabel('Signal');
    title(sprintf('SNR = %.1f dB, SNR_{seg} = %.1f dB',glo,seg));
    axh(1)=gca;
    subplot(312);
    plot((1:length(r))/fs,r);
    ylabel('Reference');
    axh(2)=gca;
    subplot(313);
    snv=snf;
    snv(~vf)=NaN;
    snu=snf;
    snu(vf>0)=NaN;
    plot([1 nr]/fs,[glo seg; glo seg],':k',((1:nf)*kf+(1-kf)/2)/fs,snv,'-b',((1:nf)*kf+(1-kf)/2)/fs,snu,'-r');
    ylabel('Frame SNR');
    xlabel('Time (s)');
    axh(3)=gca;
    linkaxes(axh,'x');
end

author: decouples


text:
if (  isempty(which('units.mat')) ...
   || isempty(which('nn_prepare')) ...
   || isempty(which('genbyode')) ...
   || isempty(which('opentstool')) ...
   || isempty(which('nn_prepare')) )
	settspath
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function movex_text(h_txt,x_pos)
%MOVEX_TEXT moves text horizontally.

% Copyright 2004-2010 The MathWorks, Inc.

FmtSpec=getappdata(get(get(h_txt,'parent'),'parent'),'FmtSpec');
msg=sprintf(FmtSpec,x_pos);
pos=get(h_txt,'position');
pos(1)=x_pos;
set(h_txt,'Position',pos,'String',msg)

author: Sable

repo: mcbench-benchmarks


text:
function phi = basis_brick27 ( n, p )

%*****************************************************************************80
%
%% BASIS_BRICK20: BRICK20 basis functions at natural coordinates.
%
%  Discussion:
%
%        8----19---7
%       /|         /
%     20 |   26   /|
%     /          / |
%    5----17----6  |
%    |   |      |  |
%    |  16---24-|-15
%    |  /|      | /|
%    |25 |  27  |23|        t
%    |/         |/ |        |   s
%   13----22---14  |        |  /
%    |   |      |  |        | /
%    |   |      |  |        |/
%    |   4--11--|--3        0---------r
%    |  /       | /
%    | 12   21  |10
%    |/         |/
%    1----9-----2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real P(3,N), natural coordinates of evaluation points.
%
%    Output, real PHI(27,N), the basis function values.
%
  rm(1:n) = p(1,1:n) + 1.0;
  rz(1:n) = p(1,1:n);
  rp(1:n) = p(1,1:n) - 1.0;

  sm(1:n) = p(2,1:n) + 1.0;
  sz(1:n) = p(2,1:n);
  sp(1:n) = p(2,1:n) - 1.0;

  tm(1:n) = p(3,1:n) + 1.0;
  tz(1:n) = p(3,1:n);
  tp(1:n) = p(3,1:n) - 1.0;

  phi(1,1:n)  =         rz .* rp       .* sz .* sp       .* tz .* tp / 8.0;
  phi(2,1:n)  =   rm .* rz             .* sz .* sp       .* tz .* tp / 8.0;
  phi(3,1:n)  =   rm .* rz       .* sm .* sz             .* tz .* tp / 8.0;
  phi(4,1:n)  =         rz .* rp .* sm .* sz             .* tz .* tp / 8.0;
  phi(5,1:n)  =         rz .* rp       .* sz .* sp .* tm .* tz       / 8.0;
  phi(6,1:n)  =   rm .* rz             .* sz .* sp .* tm .* tz       / 8.0;
  phi(7,1:n)  =   rm .* rz       .* sm .* sz       .* tm .* tz       / 8.0;
  phi(8,1:n)  =         rz .* rp .* sm .* sz       .* tm .* tz       / 8.0;

  phi(9,1:n)  = - rm       .* rp       .* sz .* sp       .* tz .* tp / 4.0;
  phi(10,1:n) = - rm .* rz       .* sm       .* sp       .* tz .* tp / 4.0;
  phi(11,1:n) = - rm       .* rp .* sm .* sz             .* tz .* tp / 4.0;
  phi(12,1:n) = -       rz .* rp .* sm       .* sp       .* tz .* tp / 4.0;
  phi(13,1:n) = -       rz .* rp       .* sz .* sp .* tm       .* tp / 4.0;
  phi(14,1:n) = - rm .* rz             .* sz .* sp .* tm       .* tp / 4.0;
  phi(15,1:n) = - rm .* rz       .* sm .* sz       .* tm       .* tp / 4.0;
  phi(16,1:n) = -       rz .* rp .* sm .* sz       .* tm       .* tp / 4.0;
  phi(17,1:n) = - rm       .* rp       .* sz .* sp .* tm .* tz       / 4.0;
  phi(18,1:n) = - rm .* rz       .* sm       .* sp .* tm .* tz       / 4.0;
  phi(19,1:n) = - rm       .* rp .* sm .* sz       .* tm .* tz       / 4.0;
  phi(20,1:n) = -       rz .* rp .* sm       .* sp .* tm .* tz       / 4.0;

  phi(21,1:n) =   rm       .* rp .* sm       .* sp       .* tz .* tp / 2.0;
  phi(22,1:n) =   rm       .* rp       .* sz .* sp .* tm       .* tp / 2.0;
  phi(23,1:n) =   rm .* rz       .* sm       .* sp .* tm       .* tp / 2.0;
  phi(24,1:n) =   rm       .* rp .* sm .* sz       .* tm       .* tp / 2.0;
  phi(25,1:n) =         rz .* rp .* sm       .* sp .* tm       .* tp / 2.0;
  phi(26,1:n) =   rm       .* rp .* sm       .* sp .* tm .* tz       / 2.0;

  phi(27,1:n) = - rm       .* rp .* sm       .* sp .* tm       .* tp;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
load('.\faust_synthetic\distance_matrix\tr_reg_095.mat')
load('.\faust_synthetic\shapes\tr_reg_095.mat')
trisurf(TRIV,VERT(:,1),VERT(:,2),VERT(:,3),D(1,:)); colorbar; colormap flag; axis equal; shading interp; axis off
hold
scatter3(VERT(1,1),VERT(1,2),VERT(1,3),30,'g','filled');
title('Geodesic Distance from Source Point','FontSize',28)

author: OshriHalimi


text:
function one_result=seg_eva_gen_result_from_con_mat(one_con_mat, exclude_class_idxes)

        class_num=size(one_con_mat, 1);

        accuracy_classes=zeros(class_num, 1);
        inter_union_score_classes=zeros(class_num, 1);

        if nargin<2
            exclude_class_idxes=[];
        end

        valid_gt_class_sel=true(class_num, 1);

        global_pos_num=0;
        global_predict_num=0;

        for c_idx=1:class_num

            one_true_pos_num=one_con_mat(c_idx,c_idx);
            one_gt_pos_num=sum(one_con_mat(c_idx,:));
            one_predict_pos_num=sum(one_con_mat(:, c_idx));

            global_pos_num=global_pos_num+one_true_pos_num;
            global_predict_num=global_predict_num+one_predict_pos_num;

            if one_gt_pos_num>0
                one_accuracy=one_true_pos_num/(one_gt_pos_num+eps);
                one_inter_union_score=one_true_pos_num/...
                    (one_gt_pos_num+one_predict_pos_num-one_true_pos_num+eps);

                accuracy_classes(c_idx)=one_accuracy;
                inter_union_score_classes(c_idx)=one_inter_union_score;
            else
                valid_gt_class_sel(c_idx)=false;
            end

        end

        valid_class_sel=true(class_num, 1);
        valid_class_sel(exclude_class_idxes)=false;
        valid_class_sel=valid_class_sel & valid_gt_class_sel;

        accuracy_per_class=mean(accuracy_classes(valid_class_sel));
        inter_union_score_per_class=mean(inter_union_score_classes(valid_class_sel));
        accuracy_global=global_pos_num/global_predict_num;

        one_result.confusion_mat=one_con_mat;
        one_result.class_num=class_num;

        one_result.accuracy_global=accuracy_global;
        one_result.accuracy_per_class=accuracy_per_class;
        one_result.inter_union_score_per_class=inter_union_score_per_class;
        one_result.accuracy_classes=accuracy_classes;
        one_result.inter_union_score_classes=inter_union_score_classes;

        one_result.valid_class_num=nnz(valid_class_sel);
        one_result.valid_class_idxes=find(valid_class_sel);
        one_result.exclude_class_idxes=exclude_class_idxes;

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function fh = isolate_axes(ah, vis)
%ISOLATE_AXES Isolate the specified axes in a figure on their own
%
% Examples:
%   fh = isolate_axes(ah)
%   fh = isolate_axes(ah, vis)
%
% This function will create a new figure containing the axes/uipanels
% specified, and also their associated legends and colorbars. The objects
% specified must all be in the same figure, but they will generally only be
% a subset of the objects in the figure.
%
% IN:
%    ah - An array of axes and uipanel handles, which must come from the
%         same figure.
%    vis - A boolean indicating whether the new figure should be visible.
%          Default: false.
%
% OUT:
%    fh - The handle of the created figure.

% Copyright (C) Oliver Woodford 2011-2013

% Thank you to Rosella Blatt for reporting a bug to do with axes in GUIs
% 16/03/12: Moved copyfig to its own function. Thanks to Bob Fratantonio
%           for pointing out that the function is also used in export_fig.m
% 12/12/12: Add support for isolating uipanels. Thanks to michael for suggesting it
% 08/10/13: Bug fix to allchildren suggested by Will Grant (many thanks!)
% 05/12/13: Bug fix to axes having different units. Thanks to Remington Reid for reporting
% 21/04/15: Bug fix for exporting uipanels with legend/colorbar on HG1 (reported by Alvaro
%           on FEX page as a comment on 24-Apr-2014); standardized indentation & help section
% 22/04/15: Bug fix: legends and colorbars were not exported when exporting axes handle in HG2

    % Make sure we have an array of handles
    if ~all(ishandle(ah))
        error('ah must be an array of handles');
    end
    % Check that the handles are all for axes or uipanels, and are all in the same figure
    fh = ancestor(ah(1), 'figure');
    nAx = numel(ah);
    for a = 1:nAx
        if ~ismember(get(ah(a), 'Type'), {'axes', 'uipanel'})
            error('All handles must be axes or uipanel handles.');
        end
        if ~isequal(ancestor(ah(a), 'figure'), fh)
            error('Axes must all come from the same figure.');
        end
    end
    % Tag the objects so we can find them in the copy
    old_tag = get(ah, 'Tag');
    if nAx == 1
        old_tag = {old_tag};
    end
    set(ah, 'Tag', 'ObjectToCopy');
    % Create a new figure exactly the same as the old one
    fh = copyfig(fh); %copyobj(fh, 0);
    if nargin < 2 || ~vis
        set(fh, 'Visible', 'off');
    end
    % Reset the object tags
    for a = 1:nAx
        set(ah(a), 'Tag', old_tag{a});
    end
    % Find the objects to save
    ah = findall(fh, 'Tag', 'ObjectToCopy');
    if numel(ah) ~= nAx
        close(fh);
        error('Incorrect number of objects found.');
    end
    % Set the axes tags to what they should be
    for a = 1:nAx
        set(ah(a), 'Tag', old_tag{a});
    end
    % Keep any legends and colorbars which overlap the subplots
    % Note: in HG1 these are axes objects; in HG2 they are separate objects, therefore we
    %       don't test for the type, only the tag (hopefully nobody but Matlab uses them!)
    lh = findall(fh, 'Tag', 'legend', '-or', 'Tag', 'Colorbar');
    nLeg = numel(lh);
    if nLeg > 0
        set([ah(:); lh(:)], 'Units', 'normalized');
        try
            ax_pos = get(ah, 'OuterPosition'); % axes and figures have the OuterPosition property
        catch
            ax_pos = get(ah, 'Position'); % uipanels only have Position, not OuterPosition
        end
        if nAx > 1
            ax_pos = cell2mat(ax_pos(:));
        end
        ax_pos(:,3:4) = ax_pos(:,3:4) + ax_pos(:,1:2);
        try
            leg_pos = get(lh, 'OuterPosition');
        catch
            leg_pos = get(lh, 'Position');  % No OuterPosition in HG2, only in HG1
        end
        if nLeg > 1;
            leg_pos = cell2mat(leg_pos);
        end
        leg_pos(:,3:4) = leg_pos(:,3:4) + leg_pos(:,1:2);
        ax_pos = shiftdim(ax_pos, -1);
        % Overlap test
        M = bsxfun(@lt, leg_pos(:,1), ax_pos(:,:,3)) & ...
            bsxfun(@lt, leg_pos(:,2), ax_pos(:,:,4)) & ...
            bsxfun(@gt, leg_pos(:,3), ax_pos(:,:,1)) & ...
            bsxfun(@gt, leg_pos(:,4), ax_pos(:,:,2));
        ah = [ah; lh(any(M, 2))];
    end
    % Get all the objects in the figure
    axs = findall(fh);
    % Delete everything except for the input objects and associated items
    delete(axs(~ismember(axs, [ah; allchildren(ah); allancestors(ah)])));
end

function ah = allchildren(ah)
    ah = findall(ah);
    if iscell(ah)
        ah = cell2mat(ah);
    end
    ah = ah(:);
end

function ph = allancestors(ah)
    ph = [];
    for a = 1:numel(ah)
        h = get(ah(a), 'parent');
        while h ~= 0
            ph = [ph; h];
            h = get(h, 'parent');
        end
    end
end

author: optimaltransport


text:
function [pkg, filenames] = MavenDownload(varargin)
    %MAVENDOWNLOAD  Downloads JAR artifacts from Maven Central
    %
    %     pkg = MavenDownload(groupId, artifactId)
    %     pkg = MavenDownload(..., 'OptionName',optionValue, ...)
    %     [pkg, filenames] = MavenDownload(...)
    %
    % ## Input
    % * __groupId__ group Id (organization name), may be empty if there is no
    %   ambiguity.
    % * __artifactId__ artifact Id (module name).
    %
    % ## Options
    % * __Version__ (optional) package version. If not specified, the
    %   latest version is requested.
    % * __SkipDependencies__ Decide if dependencies should not be processed.
    %   Default false.
    % * __Download__ Decide if packages should be downloaded. Default true.
    % * __Verbose__ Verbosity messages. Default true.
    % * __OutDir__ Directory where to place downloaded JAR files. By default,
    %   the current directory is used.
    %
    % ## Output
    % * __pkg__ struct describing the package and its dependencies
    % * __filenames__ Full paths to downloaded JAR files.
    %
    % This is a poor man's implementation of a Java package manager.
    % It searches the [Maven Central Repository](http://search.maven.org/)
    % and downloads a package JAR file and its dependencies.
    %
    % The Maven REST API is documented [here](http://search.maven.org/#api).
    %
    % ## Example
    %
    %     MavenDownload('org.jsoup', 'jsoup', 'OutDir','jsoup')
    %     MavenDownload('org.commonjava.googlecode.markdown4j', 'markdown4j', 'OutDir','markdown4j')
    %     MavenDownload('org.markdownj', 'markdownj-core', 'OutDir','markdownj')
    %     MavenDownload('org.tautua.markdownpapers', 'markdownpapers-core', 'OutDir','markdownpapers')
    %     MavenDownload('org.pegdown', 'pegdown', 'OutDir','pegdown')
    %     MavenDownload('com.github.rjeschke', 'txtmark', 'OutDir','txtmark')
    %     MavenDownload('com.atlassian.commonmark', 'commonmark', 'OutDir','commonmark-java')
    %     MavenDownload('com.atlassian.commonmark', 'commonmark-ext-gfm-tables', 'OutDir','commonmark-java')
    %
    % See also: webread, websave, urlencode, jsondecode
    %

    % parse inputs
    nargoutchk(0,2);
    opts = parse_options(varargin{:});

    % find requested package
    pkg = struct();
    pkg.groupId = opts.groupId;
    pkg.artifactId = opts.artifactId;
    pkg.version = opts.version;

    % search for latest version if none specified
    if isempty(pkg.groupId) || isempty(pkg.version)
        pkg = package_search(pkg, opts);
    end

    % find dependencies
    if ~opts.SkipDependencies
        pkg = package_dependencies(pkg, opts);
    end

    % download JAR files (including all dependencies)
    if opts.Download
        filenames = package_download(pkg, opts);
    end
end

function opts = parse_options(varargin)
    %PARSE_OPTIONS  Help function to parse function inputs

    % helper function to validate true/false arguments
    isbool = @(x) isscalar(x) && (islogical(x) || isnumeric(x));

    p = inputParser();
    p.addRequired('groupId', @ischar);
    p.addRequired('artifactId', @(x) ischar(x) && ~isempty(x));
    p.addParameter('version', '', @ischar);
    p.addParameter('SkipDependencies', false, isbool);
    p.addParameter('Download', true, isbool);
    p.addParameter('Verbose', true, isbool);
    p.addParameter('OutDir', '', @ischar);
    p.parse(varargin{:});

    opts = p.Results;
    opts.SkipDependencies = logical(opts.SkipDependencies);
    opts.Download = logical(opts.Download);
    opts.Verbose = logical(opts.Verbose);
end

function pkg = package_search(pkg, opts)
    %PACKAGE_SEARCH  Search for package from Maven by artifact Id

    % search by artifact id
    if opts.Verbose, disp(['Searching... ' pkg.artifactId]); end
    url = 'http://search.maven.org/solrsearch/select';
    query = {'q',sprintf('a:"%s"',pkg.artifactId), 'rows',5, 'wt','json'};
    data = webread(url, query{:}, weboptions('ContentType','json'));

    % check number of results
    assert(~isempty(data), 'No data received');
    res = data.response;
    if res.numFound == 0
        error('mexopencv:err', ...
            'artifactId = "%s" was not found', pkg.artifactId);
    elseif res.numFound > 1 && isempty(pkg.groupId)
        %disp({res.docs.g})
        error('mexopencv:err', ...
            'Ambiguous artifactId = "%s", specify groupId', pkg.artifactId);
    end

    % match groupId if specified
    idx = 1;
    if ~isempty(pkg.groupId)
        g = {res.docs.g};
        if res.numFound > 1
            [~,idx] = ismember(validatestring(pkg.groupId, g), g);
        else
            assert(strcmp(pkg.groupId, g{1}), ...
                'groupId = "%s" does not match "%s"', pkg.groupId, g{1});
        end
    end

    % unique project name
    pkg.groupId = res.docs(idx).g;
    pkg.artifactId = res.docs(idx).a;
    if isempty(pkg.version)
        pkg.version = res.docs(idx).latestVersion;
    end
end

function pkg = package_dependencies(pkg, opts)
    %PACKAGE_DEPENDENCIES  Fill package dependencies

    %TODO: cache dependencies to avoid fetching same POM multiple times
    %TODO: cyclic dependencies

    % fetch POM XML file
    fpath = package_filepath(pkg, '.pom');
    if opts.Verbose, disp(['Fetching... ' fpath]); end
    url = 'http://search.maven.org/remotecontent';
    query = {'filepath',fpath};
    pom = webread(url, query{:}, weboptions('ContentType','xmldom'));

    % parse it for dependencies (and transitive dependencies)
    pkg.deps = parse_pom(pom);
    for i=1:numel(pkg.deps)
        pkg.deps(i) = package_dependencies(pkg.deps(i), opts);
    end
end

function filenames = package_download(pkg, opts)
    %PACKAGE_DOWNLOAD  Download JAR package and its dependencies

    % fetch JAR file
    [fpath, filename] = package_filepath(pkg, '.jar');
    if ~isempty(opts.OutDir)
        filename = fullfile(opts.OutDir, filename);
        if ~isdir(opts.OutDir)
            mkdir(opts.OutDir);
        end
    end
    if exist(filename, 'file') ~= 2
        if opts.Verbose, disp(['Downloading... ' fpath]); end
        url = 'http://search.maven.org/remotecontent';
        query = {'filepath',fpath};
        filename = websave(filename, url, query{:});
    else
        if opts.Verbose, disp(['Skipped ' fpath]); end
    end
    filenames = filename;

    % download dependencies
    if ~opts.SkipDependencies
        fnames = cell(size(pkg.deps));
        for i=1:numel(pkg.deps)
            fnames{i} = package_download(pkg.deps(i), opts);
        end

        % combine and flatten names
        filenames = [filenames fnames];
        while any(cellfun(@iscell, filenames))
            filenames = [filenames{:}];
        end
        filenames = unique(filenames, 'stable');
    end
end

function [fpath, filename] = package_filepath(pkg, packaging)
    %PACKAGE_FILEPATH  Construct path to fetch package from Maven Central

    % fpath = org/github/com/my-project/1.0.0/my-project-1.0.0.jar
    assert(~isempty(pkg.groupId));
    assert(~isempty(pkg.artifactId));
    assert(~isempty(pkg.version));
    packaging = validatestring(packaging, ...
        {'.pom', '.jar', '-sources.jar', '-javadoc.jar'});
    filename = sprintf('%s-%s%s', pkg.artifactId, pkg.version, packaging);
    fpath = sprintf('%s/%s/%s/%s', ...
        strrep(pkg.groupId,'.','/'), pkg.artifactId, pkg.version, filename);
end

function deps = parse_pom(pom)
    %PARSE_POM  Extract dependencies from POM document DOM
    %
    % See also: xmlread
    %

    % output structure-array of dependencies
    deps = struct('groupId',{}, 'artifactId',{}, 'version',{}, 'deps',{});
    k = 1;  % insertion index

    % debug
    %{
    disp(pom.saveXML([]));
    keyboard
    %}

    % get <dependencies> tag
    nodes = pom.getElementsByTagName('dependencies').item(0);
    if isempty(nodes)
        return;  % no dependencies
    end

    % check if contains a <parent> tag
    node = pom.getElementsByTagName('parent').item(0);
    if ~isempty(node)
        p_groupId = getElementTextByTagName(node, 'groupId');
        p_artifactId = getElementTextByTagName(node, 'artifactId');
        p_version = getElementTextByTagName(node, 'version');
    else
        p_groupId = '';
        p_artifactId = '';
        p_version = '';
    end

    % loop over child <dependency> tags
    for i=1:nodes.getLength()
        node = nodes.item(i-1);

        % skip anything other than <dependency> (text/comments/etc. nodes)
        if node.getNodeType() ~= 1 || ~strcmp(node.getNodeName(), 'dependency')
            continue;
        end

        % check dependency <scope> and ignore test dependencies
        scope = getElementTextByTagName(node, 'scope');
        if ~isempty(scope) && strcmp(scope, 'test')
            continue;
        end

        % get dependency <groupId>, <artifactId>, <version>
        deps(k).groupId = getElementTextByTagName(node, 'groupId');
        deps(k).artifactId = getElementTextByTagName(node, 'artifactId');
        deps(k).version = getElementTextByTagName(node, 'version');

        %HACK: fill empty fields with parent fields
        if isempty(deps(k).groupId)
            deps(k).groupId = p_groupId;
        end
        if isempty(deps(k).artifactId)
            deps(k).artifactId = p_artifactId;
        end
        if isempty(deps(k).version)
            deps(k).version = p_version;
        end

        % prepare for next entry
        k = k + 1;
    end
end

function txt = getElementTextByTagName(node, tag)
    %GETELEMENTTEXTBYTAGNAME  Helper function to get text content of tag node

    % equivalent to: txt = node.select('tag').first().text()
    nodeTAG = node.getElementsByTagName(tag).item(0);
    if ~isempty(nodeTAG)
        txt = char(nodeTAG.getFirstChild().getTextContent());
    else
        txt = '';
    end
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function X = real(X)
%REAL (overloaded)

X.basis = real(X.basis);
X = clean(X);
if isa(X,'sdpvar')
   X.conicinfo = [0 0];
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function g = p35_g ( n, x )

%*****************************************************************************80
%
%% P35_G evaluates the gradient for problem 35.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the values of the variables.
%
%    Output, real G(N), the gradient of the objective function.
%
  g = zeros ( n, 1 );

  factor1 = 0.0;
  df1dx1 = 0.0;
  for i = 1 : 5
    y = i;
    factor1 = factor1 + y * cos ( ( y + 1.0 ) * x(1) + y );
    df1dx1 = df1dx1 - y * ( y + 1.0 ) * sin ( ( y + 1.0 ) * x(1) + y );
  end

  factor2 = 0.0;
  df2dx2 = 0.0;
  for i = 1 : 5
    y = i;
    factor2 = factor2 + y * cos ( ( y + 1.0 ) * x(2) + y );
    df2dx2 = df2dx2 - y * ( y + 1.0 ) * sin ( ( y + 1.0 ) * x(2) + y );
  end

  g(1) = df1dx1 * factor2;
  g(2) = factor1 * df2dx2;

  return
end


text:
function [pol, res, zer] = prz(zj, fj, wj)
%PRZ   Computes poles, residues, and zeros of a rational function in
%      barycentric form.
%   [POL, RES, ZER] = PRZ(ZJ, FJ, WJ) returns vectors of poles POL,
%   residues RES, and zeros ZER of the rational function defined by
%   support points ZJ, function values FJ, and barycentric weights WJ.
%
% See also AAA, PRZTRIG, REVAL.

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

m = length(wj);

% Compute poles via generalized eigenvalue problem:
B = eye(m+1);
B(1,1) = 0;
E = [0 wj.'; ones(m, 1) diag(zj)];
pol = eig(E, B);
% Remove zeros of denominator at infinity:
pol = pol(~isinf(pol));

% Compute residues via formula for res of quotient of analytic functions:
N = @(t)(1./bsxfun(@minus,t,zj.')) * (fj.*wj);
Ddiff = @(t) -((1./bsxfun(@minus,t,zj.')).^2) * wj;
res = N(pol)./Ddiff(pol);

% Compute zeros via generalized eigenvalue problem:
E = [0 (wj.*fj).'; ones(m, 1) diag(zj)];
zer = eig(E, B);
% Remove zeros of numerator at infinity:
zer = zer(~isinf(zer));

end % End of PRZ().

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [ n_data, x, fx ] = sinh_values ( n_data )

%*****************************************************************************80
%
%% SINH_VALUES returns some values of the hyperbolic sine function.
%
%  Discussion:
%
%    In Mathematica, the function can be evaluated by:
%
%      Sinh[x]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964,
%    ISBN: 0-486-61272-4,
%    LC: QA47.A34.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Cambridge University Press, 1999,
%    ISBN: 0-521-64314-7,
%    LC: QA76.95.W65.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 18;

  fx_vec = [ ...
      -74.203210577788758977, ...
       -1.1752011936438014569, ...
        0.00000000000000000000, ...
        0.10016675001984402582, ...
        0.20133600254109398763, ...
        0.30452029344714261896, ...
        0.41075232580281550854, ...
        0.52109530549374736162, ...
        0.63665358214824127112, ...
        0.75858370183953350346, ...
        0.88810598218762300657, ...
        1.0265167257081752760, ...
        1.1752011936438014569, ...
        3.6268604078470187677, ...
       10.017874927409901899, ...
       27.289917197127752449, ...
       74.203210577788758977, ...
    11013.232874703393377 ];

  x_vec = [ ...
   -5.0, ...
   -1.0, ...
    0.0, ...
    0.1, ...
    0.2, ...
    0.3, ...
    0.4, ...
    0.5, ...
    0.6, ...
    0.7, ...
    0.8, ...
    0.9, ...
    1.0, ...
    2.0, ...
    3.0, ...
    4.0, ...
    5.0, ...
   10.0 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0.0;
    fx = 0.0;
  else
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function [ signal ] = msk_mod( bit_stream,frequency,Tb,Eb )
% %UNTITLED Summary of this function goes here
% %   Detailed explanation goes here

% clc;
% clear all;
% close all;
%
% % bitstream to be tx
% bit_stream=[1 0 0 1 0 0 1 1];
%
f=frequency;
%
% % carrier frequency
% f=30;
%
% % bit energy
%
% Eb=1;
% %
% % bit period
% Tb=0.001;

%time
t=0.001 : .001 : 1;

% time for p(t)

% t2=0.01:.01:(2*Tb);

nn=length(bit_stream);

amp=sqrt(2*Eb);

% make the data stream length even
if (mod(nn,2)~=0)
    bit_stream(1,nn+1)=0;
    nn=nn+1;
end

% length of bit stream
N=length(bit_stream);

m_i=zeros(1,N/2);% Even bit stream
m_q=zeros(1,N/2);% Odd bit stream

a=1;
% Generate the odd and even stream from the input bit stream
% in bipolar nrz form
for i=1:1:N/2
    if (bit_stream(1,a)==0)
        m_i(1,i)=-1;
    else
        m_i(1,i)=1;
    end
    a=a+1;
      if (bit_stream(1,a)==0)
        m_q(1,i)=-1;
    else
        m_q(1,i)=1;
    end
    a=a+1;
end

Smsk=zeros(N/2,1000);

% for loop=1:N/2
%     if (m_i(1,loop)==1)
%         fi_k=0;
%     else
%         fi_k=pi;
%     end
%
%     Smsk(loop,:)=amp*cos(2*pi*f*t-(m_i(1,loop)*m_q(1,loop)*pi*(t/(2*Tb)))+fi_k);
% end
%
% signal=column_to_row(Smsk);

%
for loop =1:N/2
    Smsk(loop,:)=(amp*(m_i(1,loop)*sin(2*pi*t/(4*Tb))).*cos(2*pi*f*t)) + (amp*(m_q(1,loop)*cos(2*pi*(t/(4*Tb)))).*sin(2*pi*f*t));
end
signal=column_to_row(Smsk);

end


text:
function add_legend_entry(plot_obj,lh,legend_string,append_name)
% Places a (curve handle, legend string) pair in a list for when
% display_legend is called.
% Inputs:
%   lh: A curve handle.
%   legend_string: The string to be displayed in the legend for the
%     curve corresponding to lh.
%   append_name: A boolean indicating whether to prepend the system
%     name (set in 'set_system') to the 'legend_string'.

if ~strcmp(legend_string,'')
    plot_obj.handles_vec = [plot_obj.handles_vec lh];
    if exist('append_name','var') && append_name
	plot_obj.legend_strings = {plot_obj.legend_strings{:},[plot_obj.sys_name ' ' legend_string]};
    else
	plot_obj.legend_strings = {plot_obj.legend_strings{:},legend_string};
    end
end

end

author: nesl


text:
% this function return the transformer's (input argument) loss
function FValue=SP_tlosscost1(InValue)

global KWTransformersTypes OpenCircuitLosses
Tag = find(KWTransformersTypes == InValue);
FValue = OpenCircuitLosses(1,Tag);


text:
function obj = eq_pnp_func_new(in1,in2,in3,in4,in5)
coef_f0_q_sym1 = in1(:,1);
coef_f0_q_sym2 = in1(:,2);
coef_f1_q_sym1 = in2(:,1);
coef_f0_q_sym3 = in1(:,3);
coef_f1_q_sym2 = in2(:,2);
coef_f2_q_sym1 = in3(:,1);
coef_f0_q_sym4 = in1(:,4);
coef_f1_q_sym3 = in2(:,3);
coef_f2_q_sym2 = in3(:,2);
coef_f3_q_sym1 = in4(:,1);
coef_f0_q_sym5 = in1(:,5);
coef_f1_q_sym4 = in2(:,4);
coef_f2_q_sym3 = in3(:,3);
coef_f3_q_sym2 = in4(:,2);
coef_f0_q_sym6 = in1(:,6);
coef_f1_q_sym5 = in2(:,5);
coef_f2_q_sym4 = in3(:,4);
coef_f3_q_sym3 = in4(:,3);
coef_f0_q_sym7 = in1(:,7);
coef_f1_q_sym6 = in2(:,6);
coef_f2_q_sym5 = in3(:,5);
coef_f3_q_sym4 = in4(:,4);
coef_f0_q_sym8 = in1(:,8);
coef_f1_q_sym7 = in2(:,7);
coef_f2_q_sym6 = in3(:,6);
coef_f3_q_sym5 = in4(:,5);
coef_f0_q_sym9 = in1(:,9);
coef_f1_q_sym8 = in2(:,8);
coef_f2_q_sym7 = in3(:,7);
coef_f3_q_sym6 = in4(:,6);
coef_f1_q_sym9 = in2(:,9);
coef_f2_q_sym8 = in3(:,8);
coef_f3_q_sym7 = in4(:,7);
coef_f2_q_sym9 = in3(:,9);
coef_f3_q_sym8 = in4(:,8);
coef_f3_q_sym9 = in4(:,9);
coef_f0_q_sym10 = in1(:,10);
coef_f0_q_sym11 = in1(:,11);
coef_f0_q_sym20 = in1(:,20);
coef_f1_q_sym10 = in2(:,10);
coef_f0_q_sym12 = in1(:,12);
coef_f0_q_sym21 = in1(:,21);
coef_f1_q_sym11 = in2(:,11);
coef_f1_q_sym20 = in2(:,20);
coef_f2_q_sym10 = in3(:,10);
coef_f0_q_sym13 = in1(:,13);
coef_f0_q_sym22 = in1(:,22);
coef_f1_q_sym12 = in2(:,12);
coef_f1_q_sym21 = in2(:,21);
coef_f2_q_sym11 = in3(:,11);
coef_f2_q_sym20 = in3(:,20);
coef_f3_q_sym10 = in4(:,10);
coef_f0_q_sym14 = in1(:,14);
coef_f0_q_sym23 = in1(:,23);
coef_f1_q_sym13 = in2(:,13);
coef_f1_q_sym22 = in2(:,22);
coef_f2_q_sym12 = in3(:,12);
coef_f2_q_sym21 = in3(:,21);
coef_f3_q_sym11 = in4(:,11);
coef_f3_q_sym20 = in4(:,20);
coef_f0_q_sym15 = in1(:,15);
coef_f0_q_sym24 = in1(:,24);
coef_f1_q_sym14 = in2(:,14);
coef_f1_q_sym23 = in2(:,23);
coef_f2_q_sym13 = in3(:,13);
coef_f2_q_sym22 = in3(:,22);
coef_f3_q_sym12 = in4(:,12);
coef_f3_q_sym21 = in4(:,21);
coef_f0_q_sym16 = in1(:,16);
coef_f1_q_sym15 = in2(:,15);
coef_f1_q_sym24 = in2(:,24);
coef_f2_q_sym14 = in3(:,14);
coef_f2_q_sym23 = in3(:,23);
coef_f3_q_sym13 = in4(:,13);
coef_f3_q_sym22 = in4(:,22);
coef_f0_q_sym17 = in1(:,17);
coef_f1_q_sym16 = in2(:,16);
coef_f2_q_sym15 = in3(:,15);
coef_f2_q_sym24 = in3(:,24);
coef_f3_q_sym14 = in4(:,14);
coef_f3_q_sym23 = in4(:,23);
coef_f0_q_sym18 = in1(:,18);
coef_f1_q_sym17 = in2(:,17);
coef_f2_q_sym16 = in3(:,16);
coef_f3_q_sym15 = in4(:,15);
coef_f3_q_sym24 = in4(:,24);
coef_f0_q_sym19 = in1(:,19);
coef_f1_q_sym18 = in2(:,18);
coef_f2_q_sym17 = in3(:,17);
coef_f3_q_sym16 = in4(:,16);
coef_f1_q_sym19 = in2(:,19);
coef_f2_q_sym18 = in3(:,18);
coef_f3_q_sym17 = in4(:,17);
coef_f2_q_sym19 = in3(:,19);
coef_f3_q_sym18 = in4(:,18);
coef_f3_q_sym19 = in4(:,19);
q0 = in5(1,:);
q1 = in5(2,:);
q2 = in5(3,:);
q3 = in5(4,:);
t92 = q0.^2;
t93 = q1.^2;
t94 = q2.^2;
t95 = q3.^2;
t96 = t92.^2;
obj = [-coef_f1_q_sym11.*t92+coef_f0_q_sym18.*t93-coef_f1_q_sym1.*t96+coef_f0_q_sym12.*t93.^2+coef_f0_q_sym11.*q0.*q1-coef_f1_q_sym18.*q0.*q1-coef_f1_q_sym22.*q0.*q2+coef_f0_q_sym22.*q1.*q2-coef_f1_q_sym24.*q0.*q3+coef_f0_q_sym24.*q1.*q3+coef_f0_q_sym2.*t92.*t93-coef_f1_q_sym5.*t92.*t93-coef_f1_q_sym8.*t92.*t94-coef_f1_q_sym10.*t92.*t95+coef_f0_q_sym15.*t93.*t94+coef_f0_q_sym17.*t93.*t95+coef_f0_q_sym1.*q0.*q1.*t92-coef_f1_q_sym2.*q0.*q1.*t92+coef_f0_q_sym5.*q0.*q1.*t93-coef_f1_q_sym12.*q0.*q1.*t93+coef_f0_q_sym8.*q0.*q1.*t94-coef_f1_q_sym15.*q0.*q1.*t94+coef_f0_q_sym10.*q0.*q1.*t95-coef_f1_q_sym17.*q0.*q1.*t95-coef_f1_q_sym3.*q0.*q2.*t92+coef_f0_q_sym6.*q0.*q2.*t93-coef_f1_q_sym13.*q0.*q2.*t93-coef_f1_q_sym19.*q0.*q2.*t94-coef_f1_q_sym21.*q0.*q2.*t95+coef_f0_q_sym3.*q1.*q2.*t92-coef_f1_q_sym4.*q0.*q3.*t92-coef_f1_q_sym6.*q1.*q2.*t92+coef_f0_q_sym7.*q0.*q3.*t93+coef_f0_q_sym13.*q1.*q2.*t93-coef_f1_q_sym14.*q0.*q3.*t93-coef_f1_q_sym20.*q0.*q3.*t94+coef_f0_q_sym19.*q1.*q2.*t94+coef_f0_q_sym21.*q1.*q2.*t95-coef_f1_q_sym23.*q0.*q3.*t95+coef_f0_q_sym4.*q1.*q3.*t92-coef_f1_q_sym7.*q1.*q3.*t92+coef_f0_q_sym14.*q1.*q3.*t93+coef_f0_q_sym20.*q1.*q3.*t94+coef_f0_q_sym23.*q1.*q3.*t95-coef_f1_q_sym9.*q2.*q3.*t92+coef_f0_q_sym16.*q2.*q3.*t93+coef_f0_q_sym9.*q0.*q1.*q2.*q3-coef_f1_q_sym16.*q0.*q1.*q2.*q3;-coef_f2_q_sym11.*t92+coef_f0_q_sym22.*t94-coef_f2_q_sym1.*t96+coef_f0_q_sym19.*t94.^2-coef_f2_q_sym18.*q0.*q1+coef_f0_q_sym11.*q0.*q2-coef_f2_q_sym22.*q0.*q2-coef_f2_q_sym24.*q0.*q3+coef_f0_q_sym18.*q1.*q2+coef_f0_q_sym24.*q2.*q3-coef_f2_q_sym5.*t92.*t93+coef_f0_q_sym3.*t92.*t94-coef_f2_q_sym8.*t92.*t94-coef_f2_q_sym10.*t92.*t95+coef_f0_q_sym13.*t93.*t94+coef_f0_q_sym21.*t94.*t95-coef_f2_q_sym2.*q0.*q1.*t92-coef_f2_q_sym12.*q0.*q1.*t93+coef_f0_q_sym6.*q0.*q1.*t94-coef_f2_q_sym15.*q0.*q1.*t94-coef_f2_q_sym17.*q0.*q1.*t95+coef_f0_q_sym1.*q0.*q2.*t92-coef_f2_q_sym3.*q0.*q2.*t92+coef_f0_q_sym5.*q0.*q2.*t93-coef_f2_q_sym13.*q0.*q2.*t93+coef_f0_q_sym8.*q0.*q2.*t94-coef_f2_q_sym19.*q0.*q2.*t94+coef_f0_q_sym10.*q0.*q2.*t95-coef_f2_q_sym21.*q0.*q2.*t95+coef_f0_q_sym2.*q1.*q2.*t92-coef_f2_q_sym4.*q0.*q3.*t92-coef_f2_q_sym6.*q1.*q2.*t92+coef_f0_q_sym12.*q1.*q2.*t93-coef_f2_q_sym14.*q0.*q3.*t93+coef_f0_q_sym9.*q0.*q3.*t94-coef_f2_q_sym20.*q0.*q3.*t94+coef_f0_q_sym15.*q1.*q2.*t94-coef_f2_q_sym23.*q0.*q3.*t95+coef_f0_q_sym17.*q1.*q2.*t95-coef_f2_q_sym7.*q1.*q3.*t92+coef_f0_q_sym16.*q1.*q3.*t94+coef_f0_q_sym4.*q2.*q3.*t92-coef_f2_q_sym9.*q2.*q3.*t92+coef_f0_q_sym14.*q2.*q3.*t93+coef_f0_q_sym20.*q2.*q3.*t94+coef_f0_q_sym23.*q2.*q3.*t95+coef_f0_q_sym7.*q0.*q1.*q2.*q3-coef_f2_q_sym16.*q0.*q1.*q2.*q3;-coef_f3_q_sym11.*t92+coef_f0_q_sym24.*t95-coef_f3_q_sym1.*t96+coef_f0_q_sym23.*t95.^2-coef_f3_q_sym18.*q0.*q1-coef_f3_q_sym22.*q0.*q2+coef_f0_q_sym11.*q0.*q3-coef_f3_q_sym24.*q0.*q3+coef_f0_q_sym18.*q1.*q3+coef_f0_q_sym22.*q2.*q3-coef_f3_q_sym5.*t92.*t93-coef_f3_q_sym8.*t92.*t94+coef_f0_q_sym4.*t92.*t95-coef_f3_q_sym10.*t92.*t95+coef_f0_q_sym14.*t93.*t95+coef_f0_q_sym20.*t94.*t95-coef_f3_q_sym2.*q0.*q1.*t92-coef_f3_q_sym12.*q0.*q1.*t93-coef_f3_q_sym15.*q0.*q1.*t94+coef_f0_q_sym7.*q0.*q1.*t95-coef_f3_q_sym17.*q0.*q1.*t95-coef_f3_q_sym3.*q0.*q2.*t92-coef_f3_q_sym13.*q0.*q2.*t93-coef_f3_q_sym19.*q0.*q2.*t94+coef_f0_q_sym9.*q0.*q2.*t95-coef_f3_q_sym21.*q0.*q2.*t95+coef_f0_q_sym1.*q0.*q3.*t92-coef_f3_q_sym4.*q0.*q3.*t92-coef_f3_q_sym6.*q1.*q2.*t92+coef_f0_q_sym5.*q0.*q3.*t93-coef_f3_q_sym14.*q0.*q3.*t93+coef_f0_q_sym8.*q0.*q3.*t94-coef_f3_q_sym20.*q0.*q3.*t94+coef_f0_q_sym10.*q0.*q3.*t95+coef_f0_q_sym16.*q1.*q2.*t95-coef_f3_q_sym23.*q0.*q3.*t95+coef_f0_q_sym2.*q1.*q3.*t92-coef_f3_q_sym7.*q1.*q3.*t92+coef_f0_q_sym12.*q1.*q3.*t93+coef_f0_q_sym15.*q1.*q3.*t94+coef_f0_q_sym17.*q1.*q3.*t95+coef_f0_q_sym3.*q2.*q3.*t92-coef_f3_q_sym9.*q2.*q3.*t92+coef_f0_q_sym13.*q2.*q3.*t93+coef_f0_q_sym19.*q2.*q3.*t94+coef_f0_q_sym21.*q2.*q3.*t95+coef_f0_q_sym6.*q0.*q1.*q2.*q3-coef_f3_q_sym16.*q0.*q1.*q2.*q3;t92+t93+t94+t95-1.0];


text:
%% HILBERT_FSU uses the fsuClusterMatlab command to run the HILBERT "computation".
%
%  Discussion:
%
%    The arguments to fsuClusterMatlab mean:
%
%    '~/matlab/hilbert_parfor' allows us to specify an output directory;
%
%    [] allows us to specify switches to the MOAB queue handler, such as more time.
%
%    'p" indicates this is a "parfor" job;
%
%    'w' indicates that the MATLAB session should wait for completion;
%
%    workers is the number of processors requested;
%
%    @hilbert_fun indicates the function we want to execute;
%
%    { n } contains input arguments to the function.
%
%    The function "hilbert_fun" must correspond to a MATLAB M-file "hilbert_fun.m"
%    and that file must be in MATLAB's path.
%
%    The easiest way to ensure that the file is in MATLAB's path is to
%    create a subdirectory called "matlab" immediately below your home directory,
%    and place the function file there (as well as any other files needed
%    by the job).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2012
%
%  Author:
%
%    John Burkardt
%
  m = 100;
  n = 100;
  workers = 4;

  tic

  fsuClusterMatlab ( '~/matlab/hilbert_parfor', [], 'p', 'w', workers, ...
    @hilbert_fun, { m, n } )

  toc


text:
function hello_test ( )

%*****************************************************************************80
%
%% HELLO_TEST tests the HELLO program.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2012
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HELLO_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Demonstrate the use of MATLAB''s MEX compiler.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Compile the C file "hello.c" with the command\n' );
  fprintf ( 1, '    mex hello.c\n' );

  mex hello.c

  fprintf ( 1, '  A compiled file called "hello.mexmaci64" has been created.\n' );
  fprintf ( 1, '  I''ll convince you of this by typing "ls".\n' );
  fprintf ( 1, '\n' );

  ls

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now, the command "hello ( )" will execute the\n' );
  fprintf ( 1, '  compiled file.\n' );

  hello ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now I will delete the "hello.mexmaci64" file.\n' );

  delete hello.mexmaci64
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HELLO_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function cdf = gompertz_cdf ( x, a, b )

%*****************************************************************************80
%
%% GOMPERTZ_CDF evaluates the Gompertz CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Johnson, Kotz, and Balakrishnan,
%    Continuous Univariate Distributions, Volume 2, second edition,
%    Wiley, 1994, pages 25-26.
%
%  Parameters:
%
%    Input, real X, the argument of the CDF.
%
%    Input, real A, B, the parameters of the PDF.
%    1 < A, 0 < B.
%
%    Output, real CDF, the value of the CDF.
%
  if ( x <= 0.0 )
    cdf = 0.0;
  else
    cdf = 1.0 - exp ( - b * ( a^x - 1.0 ) / log ( a ) );
  end

  return
end


text:
function [yNext,didConverge]=implicitWeakRungeKStep(y,t,a,B,deltaT,algorithm,aCur,BCur,useNewton,useGaussian,maxIter,RelTol,AbsTol)
%%IMPLICITWEAKRUNGEKSTEP Perform a single step of an implicit weak
%           stochastic Runge-Kutta method under It� calculus. This
%           integrates d-dimensional stochastic differential equation of
%           the form:
%           dy=a(y,t)*dt+B(y,t)*dW
%           where dW is the differential of an m-dimensional Wiener
%           process. As the stepsize used decreases, weak methods converge
%           such that integrals with the random process are a measure are
%           correct (for example, to determine moments). However, they do
%           not converge to the optimal path, unlike strong methods.
%
%INPUTS: y The dX1 initial value of the random process.
%        t The scalar initial time of the random process. If an empty
%          matrix is passed, t=0 is used.
%        a A function handle to the drift function. This is called as
%          a(y,t) and returns a dX1 vector. If one wishes to use Newton's
%          method for the implicit iteration, then the calling format is
%          [aVal,papy]=a(y,t), where papy is the dXd matrix of partial
%          derivatives of a with respect to the elements of y papy(:,i) is
%          the derivative with respect to the ith component of y.
%        B A function handle to the diffusion matrix function. This is
%          called as B(y,t) and returns a dXm matrix.
%   deltaT The time increment over which the step is taken.
% algorithm A parameter specifying the algorithm to use. Possible values
%          are:
%          0 Use the implicit order 2.0 weak scheme for scalar noise from
%            Equation 5.12 of Chapter 15.5 of [1], which is the same as
%            Equation 4.12 of Chapter 15.4 and is the implicit form of
%            Equation 1.1 in Chapter 15.1 of [1]. This requires that m=1.
%          1 Use the autonomous implicit order 2.0 weak scheme from
%            Equation 5.14 of Chapter 15.5 of [1], which is the same as
%            Equation 4.13 of Chapter 15.4 and is the implicit form of
%            Equation 1.3 in Chapter 15.1 of [1]. This requires that a and
%            B not depend on t.
%          2 Use the autonomous explicit order 2.0 weak scheme for additive
%            noise that comes as a special case of Equation 5.14 of Chapter
%            15.5 of [1]. This requires that a and B not depend on t and
%            that B does not depend on x.
% aCur, BCur Often one might already have the values a(x,t) and B(x,t). If
%          so, then they should be provided as the dX1 and dXm aCur and
%          BCur to avoid recalculation. If unavailable, these values can be
%          omitted or empty matrices passed.
% useNewton Indicates whether the implicit iteration should be performed
%          using Newton's method or fixed-point iteration. The default if
%          omitted or an empty matrix is passed is false (fixed-point
%          iteration). If Newton's method is used, then papy must be
%          returned by the function a.
% useGaussian Algorithms 0 and 1 have a choice of how the random component
%          is generated. If useGaussian=true, then Gaussian random
%          variables will be used. Otherwise, simpler random variables
%          having the same desired moment properties will be used. the
%          default if omitted or an empty matrix is passed is true.
%  maxIter The maximum number of iterations to perform. The default if
%          omitted or an empty matrix is passed is 2.
% RelTol, AbsTol The relative and absolute tolerances on the iterations
%          before declaring convergence. If these are set to 0 (the default
%          if omitted or empty matrices are passed), then the algorithm
%          will just iterate for the maximum number of iterations. the
%          tolerances apply to each element of x. Convergence is declared
%          if all(diff<=AbsTol)||all(diff<=RelTol*abs(yNext)).
%
%OUTPUTS: yNext The estimated value of the process after taking a step of
%               deltaT. This is a random value.
%   didConverge If RelTol and/or AbsTol are not zero and maxIter>0, then
%               this indicates whether the iterations converged to the
%               desired accuracy. Otherwise, this is just an empty matrix.
%
%EXAMPLE 1:
%This is an example of a nonlinear scalar problem with non-additive noise
%where an explicit solution is available as a basis of comparison. In
%Chapter 4.4 of [1], the stochastic differential equation and its solution
%are from Equation 4.40. We compare the performance of the explicit
%solution (taking the expected value using cubature integration) and the
%implicit solution.
% rng(1)%Make exact run repeatable.
% numMC=1e3;
% numSteps=1;
% algorithm=1;
% useNewton=true;
% deltaT=1.1;
% y0=0.5;
% aDrift=@(y,t)((1/3)*y^(1/3));
% BDiff=@(y,t)(y^(2/3));
% papy=@(y)1/(9*y^(2/3));
% aFun=@(y,t)dealRobust(aDrift(y,t),papy(y));
% explSim=@(W)(y0^(1/3)+(1/3)*W)^3;
%
% %Take the expected value of the explicit solution using quadrature
% %integration.
% [xi,w]=quadraturePoints1D(6);%2*6-1=11th order.
% numPts=length(w);
% xi=sqrt(deltaT)*xi;
% muCub=0;
% for k=1:numPts
%     muCub=muCub+w(k)*explSim(xi(:,k));
% end
%
% valsRKImp=zeros(1,numMC);
% valsRK=zeros(1,numMC);
% for curMC=1:numMC
%     y=y0;
%     yExp=y0;
%     t=0;
%     for curStep=1:numSteps
%         s=rng();%Record the state prior to generating the random variables.
%         y=implicitWeakRungeKStep(y,t,aFun,BDiff,deltaT/numSteps,algorithm,[],[],useNewton);
%
%         rng(s)%Drive the explicit step with the same random process as the
%               %implicit one.
%         yExp=weakRungeKStep(yExp,t,aDrift,BDiff,deltaT/numSteps,algorithm);
%         t=t+deltaT/numSteps;
%     end
%     valsRKImp(curMC)=y;
%     valsRK(curMC)=yExp;
% end
% muRKImp=mean(valsRKImp);
% muRK=mean(valsRK);
% abs((muRKImp-muCub)./muCub)%Relative mean error, implicit.
% abs((muRK-muCub)./muCub)%Relative mean error, explicit.
%The implicit error will be about 0.0086 and the explicit error will be
%about 0.0153. Thus, the algorithm improves the explicit method in this
%instance.
%
%EXAMPLE 2:
%Here, we compare the implicit strong Taylor method to the explicit method
%on a linear model. This is done with the same noise driving both
%processes.
% algorithm=2;
% useNewton=false;
% numMC=1e3;
% deltaT=1/3;
% numSteps=5;
% y0=[1/4;-12];
% A=[1.1,0.1;
%    -0.2,2.2];
% D=[1.5,-0.4;
%    0.1,1];
% d=size(D,1);
% [F,Q]=linDynMod2Disc(deltaT,A,D);
% mu=F*y0;
% P=Q;
% aDrift=@(x,t)(A*x);
% BDiff=@(x,t)(D);
%
% aFun=@(y,t)dealRobust(aDrift(y,t),A);
%
% valsRKImp=zeros(d,numMC);
% valsRK=zeros(d,numMC);
% for curMC=1:numMC
%     yExpl=y0;
%     y=y0;
%     t=0;
%     for curStep=1:numSteps
%         s=rng();%Record the state prior to generating the random variables.
%         y=implicitWeakRungeKStep(y,t,aFun,BDiff,deltaT/numSteps,algorithm,[],[],useNewton);
%
%         rng(s)%Drive the explicit step with the same random process as the
%               %implicit one.
%         yExpl=weakRungeKStep(yExpl,t,aDrift,BDiff,deltaT/numSteps,algorithm);
%         t=t+deltaT/numSteps;
%     end
%     valsRKImp(:,curMC)=y;
%     valsRK(:,curMC)=yExpl;
% end
% [muRKImp,PRKImp]=calcMixtureMoments(valsRKImp);
% [muRK,PRK]=calcMixtureMoments(valsRK);
% norm(muRKImp-mu)./norm(mu)%Relative mean error, implicit.
% norm(PRKImp-P,'fro')./norm(P,'fro')%Relative variance error,implicit.
% norm(muRK-mu)./norm(mu)%Relative mean error, explicit.
% norm(PRK-P,'fro')./norm(P,'fro')%Relative variance error,explicit.
%One will typically see that the relative errors of the implicit method are
%notably better than the explicit method.
%
%EXAMPLE 3:
%This is an example of what is considered a "stiff" problem in Section 12.2
%of [1] with non-additive noise. The problem has an explicit solution in
%terms of the Wiener process W. We use quadrature integration to obtain the
%mean of that solution, which is used as the true moments for comparison.
%We compare this implicit method with an explicit step. This is a 2X1
%process with scalar noise.
% rng(1)%Make the exact run repeatable.
% algorithm=0;
% useNewton=false;
% useGaussian=false;
% numMC=1e3;
% a=25;
% b=2;
% deltaT=1/10;
% numSteps=10;
% A=[-a,a;
%     a,-a];
% B=[b,0;
%    0,b];
% aDrift=@(y,t)A*y;
% BDiff=@(y,t)B*y;
% d=2;
% papy=A;
% aFun=@(y,t)dealRobust(aDrift(y,t),papy);
% %The explicit solution; equation 2.5 in Chapter 12.2 of [1].
% expSol=@(W,deltaT,y0)(expm((A-(1/2)*B^2)*deltaT+B*W)*y0);
% y0=[1;0.1];
%
% %Take the expected value of the explicit solution using quadrature
% %integration.
% [xi,w]=quadraturePoints1D(6);%2*6-1=11th order.
% numPts=length(w);
% xi=sqrt(deltaT)*xi;
% muCub=zeros(d,1);
% for k=1:numPts
%     muCub=muCub+w(k)*expSol(xi(:,k),deltaT,y0);
% end
%
% valsRK=zeros(d,numMC);
% valsRKExplicit=zeros(d,numMC);
% for curMC=1:numMC
%     y=y0;
%     yExp=y0;
%     t=0;
%     for curStep=1:numSteps
%         s=rng();%Record the state prior to generating the random variables.
%         y=implicitWeakRungeKStep(y,t,aFun,BDiff,deltaT/numSteps,algorithm,[],[],useNewton,useGaussian);
%
%         rng(s)%Drive the explicit step with the same random process as the
%               %implicit one.
%         yExp=weakRungeKStep(yExp,t,aDrift,BDiff,deltaT/numSteps,algorithm,[],[],useGaussian);
%         t=t+deltaT/numSteps;
%     end
%     valsRK(:,curMC)=y;
%     valsRKExplicit(:,curMC)=yExp;
% end
% muRK=mean(valsRK,2);
% muRKExp=mean(valsRKExplicit,2);
% norm((muRK-muCub)./norm(muCub))%Relative mean error, implicit.
% norm((muRKExp-muCub)./norm(muCub))%Relative mean error, explicit.
%One will get a relative implicit error of about 0.0119 and a relative
%explicit error of about 0.0121. Getting rid of the random seed, the
%implicit solution is typically more accurate than the explicit. Using the
%Gaussian noise option also improves performance.
%
%REFERENCES:
%[1] P. E. Kloeden and E. Platen, Numerical Solution of Stochastic
%    Differential Equations. Berlin: Springer, 1999.
%
%November 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<13||isempty(AbsTol))
    AbsTol=0;
end

if(nargin<12||isempty(RelTol))
    RelTol=0;
end

if(nargin<11||isempty(maxIter))
    maxIter=2;
end

if(nargin<10||isempty(useGaussian))
    useGaussian=true;
end

if(nargin<9||isempty(useNewton))
    useNewton=false;
end

if(nargin<8||isempty(BCur))
    BCur=B(y,t);
end

if(nargin<7||isempty(aCur))
    aCur=a(y,t);
end

if(isempty(t))
    t=0;
end

if(algorithm>=0&&algorithm<=2)
    [yNext,fixedTerm]=weakRungeKStep(y,t,a,B,deltaT,algorithm,aCur,BCur,useGaussian);

    tNext=t+deltaT;
    if(useNewton)
        [yNext,didConverge]=implicitNewtonIter(y,a,deltaT/2,yNext,tNext,fixedTerm,maxIter,RelTol,AbsTol);
    else
        [yNext,didConverge]=implicitFixedPointIter(y,a,deltaT/2,yNext,tNext,fixedTerm,maxIter,RelTol,AbsTol);
    end
else
    error('Unknown algorithm specified.')
end
end

function [yNext,didConverge]=implicitFixedPointIter(yOld,a,coeff,yNext,tNext,fixedTerm,maxIter,RelTol,AbsTol)
%IMPLICITFIXEDITER This iterates the function
%                  yNext=a(yNext,tNext)*coeff+fixedTerm
%                  Iterations continue until the RelTol and or AbsTol
%                  conditions are met or maxIter iterations have elapsed.
%
%December 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(RelTol==0&&AbsTol==0)
        %Fixed point iteration for the maximum number of steps.
        for curIter=1:maxIter
            yNext=a(yNext,tNext)*coeff+fixedTerm;
        end
        didConverge=[];
    else%Iterate until meeting the convergence bound.
        didConverge=false;
        for curIter=1:maxIter
            yNext=a(yNext,tNext)*coeff+fixedTerm;

            diff=abs(yOld-yNext);
            if(all(diff<=AbsTol)||all(diff<=RelTol*abs(yNext)))
                didConverge=true;
                break;
            end

            yOld=yNext;
        end
    end
end

function [yNext,didConverge]=implicitNewtonIter(yOld,a,coeff,yNext,tNext,fixedTerm,maxIter,RelTol,AbsTol)
%%IMPLICITNEWTONITER This function uses Newton's method to try to solve
%                    a(yNext,tNext)*coeff+fixedTerm-yNext=0
%                    for yNext.Iterations continue until the RelTol and or
%                    AbsTol conditions are met or maxIter iterations have
%                    elapsed.
%
%December 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.

    d=size(yNext,1);

    I=eye(d,d);
    if(RelTol==0&&AbsTol==0)
        %Newton's iteration for the maximum number of steps.
        for curIter=1:maxIter
            [aVal,papy]=a(yNext,tNext);
            F=yNext-fixedTerm-coeff*aVal;
            dF=I-coeff*papy;
            yNext=yNext-dF\F;
        end
        didConverge=[];
    else%Iterate until meeting the convergence bound.
        didConverge=false;
        for curIter=1:maxIter
            [aVal,papy]=a(yNext,tNext);
            F=aVal-fixedTerm-coeff*aVal;
            dF=I-coeff*papy;
            yNext=yNext-dF\F;

            diff=abs(yOld-yNext);
            if(all(diff<=AbsTol)||all(diff<=RelTol*abs(yNext)))
                didConverge=true;
                break;
            end

            yOld=yNext;
        end
    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:
function val = get(A, prop_name)
%
%  Get psfMatrix properties from specified object and return
%  the value.
%
%     val = get(A, prop_name);
%
%  Valid choices for prop_name are:
%    'psf', 'matdata', 'type', 'boundary', 'transpose'
%

%  J. Nagy  5/2/01

switch prop_name
case 'psf'
  val = A.psf;
case 'matdata'
  val = A.matdata;
case 'type'
  val = A.type;
case 'boundary'
  val = A.boundary;
case 'transpose'
  val = A.transpose;
otherwise
  error([prop_name, 'Is not valid psfMatrix property'])
end

repo: IRtools


text:
function net = updatemodel(tr_x, tr_y, Params, net)

    [N,V]=size(tr_x);
    batchsize=V;
    run=8000;
    index=round(1+rand(1,run*batchsize)*(N-1));
    for j=1:run
        x=tr_x(index((j-1)*batchsize+1:j*batchsize),:);
        y=tr_y(index((j-1)*batchsize+1:j*batchsize),:);
        net=trainNet(x, y, Params, net);
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function line = fitLine3d(points)
%FITLINE3D Fit a 3D line to a set of points.
%
%   LINE = fitLine3d(PTS)
%
%   Example
%     pts = randn(300, 3);
%     pts = transformPoint3d(pts, createScaling3d([6 4 2]));
%     pts = transformPoint3d(pts, createRotationOx(pi/6));
%     pts = transformPoint3d(pts, createRotationOy(pi/4));
%     pts = transformPoint3d(pts, createRotationOz(pi/3));
%     pts = transformPoint3d(pts, createTranslation3d([5 4 3]));
%     elli = equivalentEllipsoid(pts);
%     figure; drawPoint3d(pts); axis equal;
%     hold on; drawEllipsoid(elli, ...
%         'drawEllipses', true, 'EllipseColor', 'b', 'EllipseWidth', 3);
%     line = fitLine3d(pts);
%     drawLine3d(line, 'color', 'm', 'LineWidth', 4);
%
%   See also
%     lines3d, equivalentEllipsoid, fitPlane
%

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2012-11-11, using Matlab 7.9.0.529 (R2009b)
% Copyright 2012-2022 INRA - Cepia Software Platform

% number of points
n = size(points, 1);

% compute centroid
center = mean(points);

% compute the covariance matrix
covPts = cov(points)/n;

% perform a principal component analysis with 2 variables,
% to extract inertia axes
[U, S] = svd(covPts);

% sort axes from greater to lower
[dummy, ind] = sort(diag(S), 'descend'); %#ok<ASGLU>

% format U to ensure first axis points to positive x direction
U = U(ind, :);
if U(1,1) < 0
    U = -U;
    % keep matrix determinant positive
    U(:,3) = -U(:,3);
end

line = [center U(:,1)'];

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [w1rp, Tau] = compute_w1rp(Pulse)
%compute_w1rms Compute the equivalent power of a rectangular pulse of
%duration of the FWHM of the shaped pulse

Trf = Pulse.Trf;
omega2 = Pulse.omega2;
if moxunit_util_platform_is_octave
    int = quad(omega2, 0, Trf);
else
    int = integral(omega2, 0, Trf);
end

if strcmp(Pulse.shape,'hard')
    Tau = Trf;
else
    x = 0:Trf/1000:Trf;
    y = omega2(x);
    Tau = fwhm(x,y);
end

w1rp = sqrt( int / Tau );

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
classdef C10MOP7 < EvoXBenchProblem
% <multi> <real> <large/none> <expensive/none>

%------------------------------- Reference --------------------------------
% Z. Lu, R. Cheng, Y. Jin, K. C. Tan, and K. Deb, Neural architecture
% search as multiobjective optimization benchmarks: Problem formulation and
% performance assessment, IEEE Transactions on Evolutionary Computation,
% 2023.
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            config.name = 'nb201';
            config.args.fidelity = 200;
            config.args.objs = ['err&params&flops&edgegpu_latency&' ...
                'edgegpu_energy&eyeriss_latency&eyeriss_energy&' ...
                'eyeriss_arithmetic_intensity'];
            config.args.dataset = 'cifar10';
            config.args.normalized_objectives = true;
            obj.Setting@EvoXBenchProblem(config);
        end
    end
end


text:
function pass = test_coeffs_vals( pref )

% Grab some preferences
if ( nargin == 0 )
    pref = chebfunpref();
end
tol = 1e3*pref.techPrefs.chebfuneps;
j = 1;
%grab some coeffs
u = diskfun(@(x,y) exp(-cos(pi*(x+y))));
v = diskfun(@(x,y) x.*sin(x.*y));
f = [u;v];

% test coeffs2:
[x, y] = coeffs2(f);
pass(j) = norm(diskfun.coeffs2diskfun(x)-u, inf) < tol;
pass(j+1) = norm(diskfun.coeffs2diskfun(y)-v, inf) < tol;
j = j+2;

%parameters
[x, y] = coeffs2(f, 50, 60);
pass(j) = norm(x-coeffs2(u, 50, 60), inf) < tol;
pass(j+1) = norm(y-coeffs2(v, 50, 60), inf) < tol;
j = j+2;

% test coeffs2diskfunv:
f2 = diskfunv.coeffs2diskfunv(coeffs2(u), coeffs2(v));
pass(j) = norm(f - f2) < tol;
j = j+1;

%test coeffs2vals:
[x, y] = coeffs2(f2);
[u, v] = diskfunv.coeffs2vals(x,y);
pass(j) = norm(diskfun.coeffs2vals(x)-u, 'inf') < tol;
pass(j+1) = norm(diskfun.coeffs2vals(y)-v, 'inf')< tol;
j= j+2;

%testvals2coeffs:
[a,b] = diskfunv.vals2coeffs(u,v);
pass(j) = norm(x - a, 'inf');
pass(j+1) = norm(y - b, 'inf');

if (nargout > 0)
    pass = all(pass(:));
end
end


text:
function [T_request_Nm, OverallRatio] = control_ForceRequestConverter(F_x_request_N, ...
                        tire_radius_m, i_c, i_f, i_r, i_gearset, ...
                        gear_in)
%__________________________________________________________________________
%% Documentation
%
% Author:       Dominik St�rk (dominik.staerk@tum.de)
%
% Start Date:   21.12.2020
%
% Description:
% This function converts the requested overall longitudinal force
% in a torque request for the engine. The effects of gear-ratios and
% -efficiencies as well as differentials and acceleration-losses are taken
% into account.
%
% Inputs:
%   F_x_request_N   - requested overall longitudinal force [N]
%   tire_radius_m   - tire`s radius [m]
%   i_c [2x1]       - central differential ratio and efficiency [-, -]
%   i_f [2x1]       - front differential ratio and efficiency [-, -]
%   i_r [2x1]       - rear differential ratio and efficiency [-, -]
%   i_gearset [4x7] - gear information
%   gear_in         - currently selected gear
%
% Outputs:
%   T_request_Nm    - requested engine torque [Nm]

p_torque_rear = 1;

%% request engine torque
T_Nm = F_x_request_N * tire_radius_m;

if gear_in ~= 0
    T_request_Nm = T_Nm / (i_c(1)*i_c(2) * ...
      ((1-p_torque_rear)*i_f(1)*i_f(2) + p_torque_rear*i_r(1)*i_r(2)) * ...
      i_gearset(2, gear_in+1)*i_gearset(4, gear_in+1));
    OverallRatio = (i_c(1)*i_c(2) * ...
      ((1-p_torque_rear)*i_f(1)*i_f(2) + p_torque_rear*i_r(1)*i_r(2)) * ...
      i_gearset(2, gear_in+1)*i_gearset(4, gear_in+1))/tire_radius_m;
else
    % use first gear even if in neutral
    T_request_Nm = T_Nm / (i_c(1)*i_c(2) * ...
      ((1-p_torque_rear)*i_f(1)*i_f(2) + p_torque_rear*i_r(1)*i_r(2)) * ...
      i_gearset(2, gear_in+2)*i_gearset(4, gear_in+1));
    % no brake torque when in neutral
    OverallRatio = 0;
end

author: TUMFTM


text:
% MP_DG_PENALTY: apply the interior penalty method between patches. To be
%  used in multipatch geometries with RT and NDL spaces (div-preserving).
%
% USAGE:
%
%  A = mp_dg_penalty (space_v, msh, interfaces, visc, Cpen)
%
% INPUT:
%
%    space_v:    multipatch space, formed by several tensor product spaces plus the connectivity (see sp_multipatch).
%    msh:        multipatch mesh, consisting of several Cartesian meshes (see msh_multipatch)
%    interfaces: interface information (see mp_geo_load)
%    visc:       function handle to compute the viscosity. For now it is
%                     assumed to be the same for all patches
%    Cpen:       penalization constant
%
% OUTPUT:
%
%    A: computed matrix, to be added to the global matrix (see mp_solve_stokes_div_conforming)
%
%   For more details, see:
%    J.A. Evans, T.J.R. Hughes
%    Isogeometric divergence-conforming B-splines for the Darcy-Stokes-Brinkman equations
%    Math. Models Meth. Appl. Sci., 2012
%
% Copyright (C) 2014, 2015, 2020 Rafael Vazquez
%
%    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/>.

function A = mp_dg_penalty (space, msh, interfaces, visc, Cpen, varargin)

if (nargin ~= 5)
  error (['The function MP_DG_PENALTY has changed in version 3, to work with multipatch classes. ' ...
        'The old version, for a cell-array of spaces, can be called with MP_DG_PENALTY_OLD'])
end

rA = []; cA = []; vA = [];

ndim = msh.ndim;
rdim = msh.rdim;

for iref = 1:numel(interfaces)
  patch(1) = interfaces(iref).patch1;
  patch(2) = interfaces(iref).patch2;
  side(1) = interfaces(iref).side1;
  side(2) = interfaces(iref).side2;

  msh_side(1) = msh_eval_boundary_side (msh.msh_patch{patch(1)}, side(1));
  msh_side(2) = msh_eval_boundary_side (msh.msh_patch{patch(2)}, side(2));
  msh_side_int    = msh_boundary_side_from_interior (msh.msh_patch{patch(1)}, side(1));
  msh_side_int(2) = msh_boundary_side_from_interior (msh.msh_patch{patch(2)}, side(2));

  sp_aux = space.sp_patch{patch(1)}.constructor (msh_side_int(1));
  sp_bnd(1) = struct (sp_precompute (sp_aux, msh_side_int(1), 'value', true, 'gradient', true));
  sp_aux = space.sp_patch{patch(2)}.constructor (msh_side_int(2));
  sp_bnd(2) = struct (sp_precompute (sp_aux, msh_side_int(2), 'value', true, 'gradient', true));

  [sp_bnd(2), msh_side(2)] = reorder_elements_and_quad_points (sp_bnd(2), msh_side(2), interfaces(iref), ndim);

% I assume there are no dicontinuities in the viscosity
  for idim = 1:rdim
    x{idim} = reshape (msh_side(1).geo_map(idim,:,:), msh_side(1).nqn, msh_side(1).nel);
  end
  coeff_at_qnodes = visc (x{:});

  charlen = (sum (msh_side(1).charlen, 1).^(1/(ndim-1)));
  charlen = repmat (charlen, msh_side(1).nqn, 1);

  for ii = 1:2
    for jj = 1:2
      [rB, cB, vB] = op_gradu_v_otimes_n (sp_bnd(jj), sp_bnd(ii), msh_side(ii), coeff_at_qnodes/2);
      [rC, cC, vC] = op_u_otimes_n_v_otimes_n (sp_bnd(jj), sp_bnd(ii), msh_side(jj), msh_side(ii), ...
          coeff_at_qnodes ./ charlen);
      vB = space.dofs_ornt{patch(ii)}(rB)' .* vB .* space.dofs_ornt{patch(jj)}(cB)';
      vC = space.dofs_ornt{patch(ii)}(rC)' .* vC .* space.dofs_ornt{patch(jj)}(cC)' * Cpen;
      rB = space.gnum{patch(ii)}(rB); cB = space.gnum{patch(jj)}(cB);
      rC = space.gnum{patch(ii)}(rC); cC = space.gnum{patch(jj)}(cC);

      rA = [rA rB cB rC];
      cA = [cA cB rB cC];
      vA = [vA -vB -vB vC];
    end
  end
end

A = sparse (rA, cA, vA, space.ndof, space.ndof);
end

function [sp_bnd, msh_side] = reorder_elements_and_quad_points (sp_bnd, msh_side, interface, ndim)

  if (ndim == 2)
    if (interface.ornt == -1)
      elem = msh_side.nel:-1:1;
      qpoints = msh_side.nqn:-1:1;
    else
      elem = 1:msh_side.nel;
      qpoints = 1:msh_side.nqn;
    end
  elseif (ndim == 3)
    elem = reshape (1:msh_side.nel, msh_side.nel_dir);
    qpoints = reshape (1:msh_side.nqn, msh_side.nqn_dir);
    if (interface.flag == -1)
      elem = elem';
      qpoints = qpoints';
    end
    if (interface.ornt1 == -1)
      elem = flipud (elem);
      qpoints = flipud (qpoints);
    end
    if (interface.ornt2 == -1)
      elem = fliplr (elem);
      qpoints = fliplr (qpoints);
    end
    elem = elem(:)';
    qpoints = qpoints(:)';
  end

  msh_side.charlen      = msh_side.charlen(qpoints,elem);
  msh_side.normal       = msh_side.normal(:,qpoints,elem);
  msh_side.jacdet       = msh_side.jacdet(qpoints,elem);
  msh_side.geo_map      = msh_side.geo_map(:,qpoints,elem);
  msh_side.geo_map_jac  = msh_side.geo_map_jac(:,:,qpoints,elem);
  msh_side.quad_weights = msh_side.quad_weights(qpoints,elem);
  msh_side.quad_nodes   = msh_side.quad_nodes(:,qpoints,elem);

  sp_bnd.nsh                      = sp_bnd.nsh(elem);
  sp_bnd.connectivity             = sp_bnd.connectivity(:,elem);
  sp_bnd.shape_functions          = sp_bnd.shape_functions(:,qpoints,:,elem);
%   sp_bnd.shape_function_divs      = sp_bnd.shape_function_divs(qpoints,:,elem);
  sp_bnd.shape_function_gradients = sp_bnd.shape_function_gradients(:,:,qpoints,:,elem);

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [M_stimrange,stimrange] = GetStimRange(option,i_fish)
% Fish 1-5: 16-stim PT/Black/White
% Fish 6-7: PT/White
% Fish 8:
%     stimset(1).name = 'PT';
%     stimset(2).name = 'OMR';
%     stimset(3).name = 'Spontaneous';
%     stimset(4).name = 'PT-shock';
%     stimset(5).name = 'OMR-shock';
%     stimset(6).name = 'RB-shock';
% Fish 9-11
%     stimset(1).name = 'PT';
%     stimset(2).name = 'OMR';
%     stimset(3).name = 'Spontaneous';
%     stimset(4).name = 'Dot/prey';
%     stimset(5).name = 'Looming';
%     stimset(6).name = 'Red|Blue';
%     stimset(7).name = 'PT-shock';
%     stimset(8).name = 'OMR-shock';
%     stimset(9).name = 'RB-shock';
% Fish 12-18 except 16
%     stimset(1).name = 'PT';
%     stimset(2).name = 'OMR';
%     stimset(3).name = 'DF';
%     stimset(4).name = 'Spontaneous';
%     stimset(5).name = 'Looming';
% Fish 16
%     stimset(1).name = 'PT';
%     stimset(2).name = 'OMR';
%     [stimset(3).name = may be 'Spontaneous' but not credible, not used]

if ~exist('option','var'),
    % default range:
    M_stimrange = {1,1,1,1,1,1,1,... 1-7
        1:3,... 8
        1:5,1:5,1:5,... 9-11
        1:5,1:5,1:5,1:5,...12-15
        1:3,...16
        1:5,1:5}; % 17-18
elseif isempty(option), % same as above,
    % default range:
    M_stimrange = {1,1,1,1,1,1,1,... 1-7
        1:3,... 8
        1:5,1:5,1:5,... 9-11
        1:5,1:5,1:5,1:5,...12-15
        1:3,...16
        1:5,1:5}; % 17-18

elseif option == 'F', % full range included in data
    M_stimrange = {1,1,1,1,1,1,1,... 1-7
        1:5,... 8
        1:9,1:8,1:9,... 9-11
        1:5,1:5,1:5,1:5,...12-15
        1:3,...16
        1:5,1:5}; % 17-18

elseif option == '2', % no spont, only PT, OMR and DF
    M_stimrange = {1,1,1,1,1,1,1,... 1-7
        1:2,... 8
        1:2,1:2,1:2,... 9-11
        1:3,1:3,1:3,1:3,...12-15
        1:2,...16
        1:3,1:3}; % 17-18

% elseif option == '3', %
%     M_stimrange = {[],[],[],[],[],[],[],... 1-7
%         1:3,... 8
%         1:3,1:3,1:3,... 9-11
%         1:3,1:3,1:3,1:3,...12-15
%         1:3,...16
%         1:3,1:3}; % 17-18

elseif option == '5', % first 5, PT/Black/White
    M_stimrange = {1,1,1,1,1,[],[],... 1-7
        [],... 8
        [],[],[],... 9-11
        [],[],[],[],...12-15
        [],...16
        [],[]};  % 17-18

elseif option == 'M', % M = multistim
    M_stimrange = {[],[],[],[],[],[],[],... 1-7
        1:3,... 8
        1:5,1:5,1:5,... 9-11
        1:5,1:5,1:5,1:5,...12-15
        1:3,...16
        1:5,1:5}; % 17-18

elseif option == 'P', % PT = phototaxis
    M_stimrange = {[],[],[],[],[],1,1,... 1-7
        1,... 8
        1,1,1,... 9-11
        1,1,1,1,...12-15
        1,...16
        1,1};  % 17-18

elseif option == 'O', % OMR
    M_stimrange = {[],[],[],[],[],[],[],... 1-7
        2,... 8
        2,2,2,... 9-11
        2,2,2,2,...12-15
        2,...16
        2,2};  % 17-18
elseif option == 'S', % Spontaneuous
        M_stimrange = {[],[],[],[],[],[],[],... 1-7
        3,... 8
        3,3,3,... 9-11
        4,4,4,4,...12-15
        [],...16
        4,4};  % 17-18

elseif option == 'D', % DF, Dark flash
    M_stimrange = {1,1,1,1,1,[],[],... % [],[],[],[],[],[],[],... 1-7
        [],... 8
        [],[],[],... 9-11
        3,3,3,3,...12-15
        [],...16
        3,3};  % 17-18

elseif option == 'L', % Looming
        M_stimrange = {[],[],[],[],[],[],[],... 1-7
        [],... 8
        5,5,5,... 9-11
        5,5,5,5,...12-15
        [],...16
        5,5};  % 17-18

elseif option == 'Y', % Dot/Prey
        M_stimrange = {[],[],[],[],[],[],[],... 1-7
        [],... 8
        4,4,4,... 9-11
        [],[],[],[],...12-15
        [],...16
        [],[]};  % 17-18
end

if exist('i_fish','var'),
    stimrange = M_stimrange{i_fish};
else % dummy
    stimrange = 1;
end
end

author: xiuyechen

repo: FishExplorer


text:
function value = halfplane_contains_point_2d ( p1, p2, p )

%*****************************************************************************80
%
%% HALFPLANE_CONTAINS_POINT_2D reports if a half-plane contains a point in 2d.
%
%  Discussion:
%
%    The halfplane is assumed to be all the points "to the left" of the
%    line that passes from P1 through P2.  Thus, one way to
%    understand where the point P = (X,Y) is, is to compute the signed
%    area of the triangle ( P1, P2, P ).
%
%    If this area is
%      positive, the point is strictly inside the halfplane,
%      zero, the point is on the boundary of the halfplane,
%      negative, the point is strictly outside the halfplane.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P1(2,1), P2(2,1), two distinct points
%    on the line defining the half plane.
%
%    Input, real P(2,1), the point to be checked.
%
%    Output, logical VALUE, is TRUE if the halfplane
%    contains the point.
%
  area_signed = 0.5 * ...
    ( p1(1,1) * ( p2(2,1) - p(2,1)  ) ...
    + p2(1,1) * ( p(2,1)  - p1(2,1) ) ...
    + p(1,1)  * ( p1(2,1) - p2(2,1) ) );

  value = ( 0.0 <= area_signed );

  return
end


text:
function linpack_s_test03 ( )

%*****************************************************************************80
%
%% TEST03 tests SCHUD.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 June 2009
%
%  Author:
%
%    John Burkardt
%
  p = 20;
  ldr = p;
  nz = 1;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST03\n' );
  fprintf ( 1, '  For single precision, general storage,\n' );
  fprintf ( 1, '  SCHUD updates a Cholesky decomposition.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  In this example, we use SCHUD to solve a\n' );
  fprintf ( 1, '  least squares problem R * b = z.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The number of equations is P = %d\n', p );
%
%  Initialize.
%
  r(1:p,1:p) = 0.0;
  z(1:p,1:nz) = 0.0;
  for i = 1 : p
    x(i) = i;
  end
%
%  Use SCHUD to form R, Z and RHO by adding X and Y a row at a time.
%  X is a row of the least squares matrix and Y the right hand side.
%
  seed = 123456789;

  for i = 1 : p
    [ row, seed ] = r4mat_uniform_01 ( 1, p, seed );
    y(1) = row(1:p) * x(1:p)';
    rho(1) = 0.0;
    [ r, z, rho, c, s ] = schud ( r, ldr, p, row, z, p, nz, y, rho );
  end
%
%  Generate the least squares solution, b = inverse ( R ) * Z.
%
  for j = 1 : nz

    b(1:p,1) = z(1:p,j);
    job = 01;

    [ b, info ] = strsl ( r, ldr, p, b, job );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Solution vector # %d\n', j );
    fprintf ( 1, '  (Should be (1,2,3...,n))\n' );
    fprintf ( 1, '\n' );

    for i = 1 : p
      if ( i <= 5 | p-5 < i )
        fprintf ( 1, '  %6d  %14f\n', i, b(i,1) );
      end
      if ( i == 5 )
        fprintf ( 1, '  ......  ..............\n' );
      end
    end

  end

  return
end


text:
function cvglmnetPlot(CVglmnet)
% Plot MSE w. errorbars for a cross validation fit of glmnet
% Call: cvglmnetPlot(CVglmnet)
%   CVglmnet: output from cvglmnet(...) function
    errorbar(log(CVglmnet.glmnetOptions.lambda),CVglmnet.cvm,CVglmnet.stderr,'linewidth',2)
    hold on
    plot(log(CVglmnet.glmnetOptions.lambda),CVglmnet.cvm,'r-o','linewidth',2)
    [val idx2] = min(fliplr(CVglmnet.cvm));
    vals=flipud(log(CVglmnet.glmnetOptions.lambda));
    plot(log([CVglmnet.lambda_min CVglmnet.lambda_min]),minmax([CVglmnet.cvlo CVglmnet.cvup]),'g--',...
        'linewidth',3)
    if CVglmnet.lambda_min ~=CVglmnet.lambda_1se
        plot(log([CVglmnet.lambda_1se CVglmnet.lambda_1se]),minmax([CVglmnet.cvlo CVglmnet.cvup]),'k--',...
            'linewidth',3)
    end
    xlabel('log(\lambda)','fontsize',12,'fontweight','bold')
    ylabel('Mean squared error +/- std. err','fontsize',12,'fontweight','bold')
    title('Nonzero elements','fontsize',12,'fontweight','bold') %use as x-top-label.
    axis tight
    H1=gca;
    ax=axis;
    H2=axes('position',get(H1,'position'));
    set(H2,'color','none')
    axis(ax);
    idx=floor(linspace(1,length(CVglmnet.glmnetOptions.lambda),7));
    labl=flipud(CVglmnet.glmnet_object.df);
    set(H2,'XAxisLocation','top')
    set(H2,'xtick',vals(idx))
    set(H2,'xticklabel',labl(idx))
    set(H2,'YAxisLocation','right')
    set(H2,'yticklabel','')
    set([H1 H2],'box','off','fontsize',12,'fontweight','bold')
    text(vals(idx2),val,['min, nz = ' num2str(labl(idx2))],'fontweight','bold','fontsize',14);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function A = catstruct(varargin)
%
% taken from MATLAB File Exchange
% https://de.mathworks.com/matlabcentral/fileexchange/7842-catstruct
%
% Copyright (c) 2015, Jos van der Geest
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% * Redistributions of source code must retain the above copyright notice, this
%   list of conditions and the following disclaimer.
%
% * Redistributions in binary form must reproduce the above copyright notice,
%   this list of conditions and the following disclaimer in the documentation
%   and/or other materials provided with the distribution
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
% DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
% FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
% DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
% SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
% CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
% OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
% OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%
% CATSTRUCT   Concatenate or merge structures with different fieldnames
%   X = CATSTRUCT(S1,S2,S3,...) merges the structures S1, S2, S3 ...
%   into one new structure X. X contains all fields present in the various
%   structures. An example:
%
%     A.name = 'Me' ;
%     B.income = 99999 ;
%     X = catstruct(A,B)
%     % -> X.name = 'Me' ;
%     %    X.income = 99999 ;
%
%   If a fieldname is not unique among structures (i.e., a fieldname is
%   present in more than one structure), only the value from the last
%   structure with this field is used. In this case, the fields are
%   alphabetically sorted. A warning is issued as well. An axample:
%
%     S1.name = 'Me' ;
%     S2.age  = 20 ; S3.age  = 30 ; S4.age  = 40 ;
%     S5.honest = false ;
%     Y = catstruct(S1,S2,S3,S4,S5) % use value from S4
%
%   The inputs can be array of structures. All structures should have the
%   same size. An example:
%
%     C(1).bb = 1 ; C(2).bb = 2 ;
%     D(1).aa = 3 ; D(2).aa = 4 ;
%     CD = catstruct(C,D) % CD is a 1x2 structure array with fields bb and aa
%
%   The last input can be the string 'sorted'. In this case,
%   CATSTRUCT(S1,S2, ..., 'sorted') will sort the fieldnames alphabetically.
%   To sort the fieldnames of a structure A, you could use
%   CATSTRUCT(A,'sorted') but I recommend ORDERFIELDS for doing that.
%
%   When there is nothing to concatenate, the result will be an empty
%   struct (0x0 struct array with no fields).
%
%   NOTE: To concatenate similar arrays of structs, you can use simple
%   concatenation:
%     A = dir('*.mat') ; B = dir('*.m') ; C = [A ; B] ;

%   NOTE: This function relies on unique. Matlab changed the behavior of
%   its set functions since 2013a, so this might cause some backward
%   compatibility issues when dulpicated fieldnames are found.
%
%   See also CAT, STRUCT, FIELDNAMES, STRUCT2CELL, ORDERFIELDS

% version 4.1 (feb 2015), tested in R2014a
% (c) Jos van der Geest
% email: jos@jasen.nl

% History
% Created in 2005
% Revisions
%   2.0 (sep 2007) removed bug when dealing with fields containing cell
%                  arrays (Thanks to Rene Willemink)
%   2.1 (sep 2008) added warning and error identifiers
%   2.2 (oct 2008) fixed error when dealing with empty structs (thanks to
%                  Lars Barring)
%   3.0 (mar 2013) fixed problem when the inputs were array of structures
%                  (thanks to Tor Inge Birkenes).
%                  Rephrased the help section as well.
%   4.0 (dec 2013) fixed problem with unique due to version differences in
%                  ML. Unique(...,'last') is no longer the deafult.
%                  (thanks to Isabel P)
%   4.1 (feb 2015) fixed warning with narginchk

narginchk(1,Inf) ;
N = nargin ;

if ~isstruct(varargin{end}),
    if isequal(varargin{end},'sorted'),
        narginchk(2,Inf) ;
        sorted = 1 ;
        N = N-1 ;
    else
        error('catstruct:InvalidArgument','Last argument should be a structure, or the string "sorted".') ;
    end
else
    sorted = 0 ;
end

sz0 = [] ; % used to check that all inputs have the same size

% used to check for a few trivial cases
NonEmptyInputs = false(N,1) ;
NonEmptyInputsN = 0 ;

% used to collect the fieldnames and the inputs
FN = cell(N,1) ;
VAL = cell(N,1) ;

% parse the inputs
for ii=1:N,
    X = varargin{ii} ;
    if ~isstruct(X),
        error('catstruct:InvalidArgument',['Argument #' num2str(ii) ' is not a structure.']) ;
    end

    if ~isempty(X),
        % empty structs are ignored
        if ii > 1 && ~isempty(sz0)
            if ~isequal(size(X), sz0)
                error('catstruct:UnequalSizes','All structures should have the same size.') ;
            end
        else
            sz0 = size(X) ;
        end
        NonEmptyInputsN = NonEmptyInputsN + 1 ;
        NonEmptyInputs(ii) = true ;
        FN{ii} = fieldnames(X) ;
        VAL{ii} = struct2cell(X) ;
    end
end

if NonEmptyInputsN == 0
    % all structures were empty
    A = struct([]) ;
elseif NonEmptyInputsN == 1,
    % there was only one non-empty structure
    A = varargin{NonEmptyInputs} ;
    if sorted,
        A = orderfields(A) ;
    end
else
    % there is actually something to concatenate
    FN = cat(1,FN{:}) ;
    VAL = cat(1,VAL{:}) ;
    FN = squeeze(FN) ;
    VAL = squeeze(VAL) ;

    [UFN,ind] = unique(FN, 'last') ;
    % If this line errors, due to your matlab version not having UNIQUE
    % accept the 'last' input, use the following line instead
    % [UFN,ind] = unique(FN) ; % earlier ML versions, like 6.5

    if numel(UFN) ~= numel(FN),
        warning('catstruct:DuplicatesFound','Fieldnames are not unique between structures.') ;
        sorted = 1 ;
    end

    if sorted,
        VAL = VAL(ind,:) ;
        FN = FN(ind,:) ;
    end

    A = cell2struct(VAL, FN);
    A = reshape(A, sz0) ; % reshape into original format
end


text:
function value = r8poly_value ( n, a, x )

%*****************************************************************************80
%
%% R8POLY_VALUE evaluates an R8 polynomial.
%
%  Discussion:
%
%    For sanity's sake, the value of N indicates the NUMBER of
%    coefficients, or more precisely, the ORDER of the polynomial,
%    rather than the DEGREE of the polynomial.  The two quantities
%    differ by 1, but cause a great deal of confusion.
%
%    Given N and A, the form of the polynomial is:
%
%      p(x) = a(1) + a(2) * x + ... + a(n-1) * x^(n-2) + a(n) * x^(n-1)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the polynomial.
%
%    Input, real A(1:N), the coefficients of the polynomial.
%    A(1) is the constant term.
%
%    Input, real X, the point at which the polynomial is to be evaluated.
%
%    Output, real VALUE, the value of the polynomial at X.
%
  value = 0.0;
  for i = n : -1 : 1
    value = value * x + a(i);
  end

  return
end


text:
function [U, S] = pca(X)
%PCA Run principal component analysis on the dataset X
%   [U, S, X] = pca(X) computes eigenvectors of the covariance matrix of X
%   Returns the eigenvectors U, the eigenvalues (on diagonal) in S
%

% Useful values
[m, n] = size(X);

% You need to return the following variables correctly.
U = zeros(n);
S = zeros(n);

% ====================== YOUR CODE HERE ======================
% Instructions: You should first compute the covariance matrix. Then, you
%               should use the "svd" function to compute the eigenvectors
%               and eigenvalues of the covariance matrix.
%
% Note: When computing the covariance matrix, remember to divide by m (the
%       number of examples).
%

Sigma = ((1.0/m) .* (X' * X));
[U, S, V] = svd(Sigma);

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

end

author: UtkarshPathrabe


text:
% Example 6.2: Robust regression using the Huber penalty
% Section 6.1.2, Figure 6.5
% Boyd & Vandenberghe "Convex Optimization"
% Original by Lieven Vandenberghe
% Adapted for CVX by Joelle Skaf - 09/07/05
%
% Compares the solution of regular Least-squares:
%           minimize    sum(y_i - alpha - beta*t_i)^2
% to the solution of the following:
%           minimize    sum( phi_h (y_i - alpha - beta*t_i)^2 )
% where phi_h is the Huber penalty function, (t_i,y_i) are data points in a
% plane.

% Input data
randn('seed',1);
rand('seed',1);

m=40;  n=2;
xex = [5;1];
pts = -10+20*rand(m,1);
A = [ones(m,1) pts];
b = A*xex + .5*randn(m,1);
outliers = [-9.5; 9];  outvals = [20; -15];
A = [A; ones(length(outliers),1), outliers];
b = [b; outvals];
m = size(A,1);
pts = [pts;outliers];

% Least Squares
fprintf(1,'Computing the solution of the least-squares problem...');

xls =  A\b;

fprintf(1,'Done! \n');

% Huber
fprintf(1,'Computing the solution of the huber-penalized problem...');

cvx_begin quiet
    variable xhub(n)
    minimize(sum(huber(A*xhub-b)))
cvx_end

fprintf(1,'Done! \n');

% Plots
figure(1);  hold off
plot(pts,b,'o', [-11; 11], [1 -11; 1 11]*xhub, '-', ...
     [-11; 11], [1 -11; 1 11]*xls, '--');
axis([-11 11 -20 25])
title('Least-square fit vs robust least-squares fit (Huber-penalized)');
xlabel('x');
ylabel('y');
legend('Data points','Huber penalty','Regular LS','Location','Best');
%print -deps robustls.eps

author: yu-jiang

repo: radpbook

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


text:
function [ x, y ] = p00_dat ( prob, data_num )

%*****************************************************************************80
%
%% P00_DAT returns the data vector for any problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROB, the number of the desired test problem.
%
%    Input, integer DATA_NUM, the number of data points.
%
%    Output, real X(DATA_NUM), the abscissa data.
%
%    Output, real Y(DATA_NUM), the ordinate data.
%
  if ( prob == 1 )
    [ x, y ] = p01_dat ( data_num );
  elseif ( prob == 2 )
    [ x, y ] = p02_dat ( data_num );
  elseif ( prob == 3 )
    [ x, y ] = p03_dat ( data_num );
  elseif ( prob == 4 )
    [ x, y ] = p04_dat ( data_num );
  elseif ( prob == 5 )
    [ x, y ] = p05_dat ( data_num );
  elseif ( prob == 6 )
    [ x, y ] = p06_dat ( data_num );
  elseif ( prob == 7 )
    [ x, y ] = p07_dat ( data_num );
  elseif ( prob == 8 )
    [ x, y ] = p08_dat ( data_num );
  elseif ( prob == 9 )
    [ x, y ] = p09_dat ( data_num );
  elseif ( prob == 10 )
    [ x, y ] = p10_dat ( data_num );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_DAT - Fatal error!\n' );
    fprintf ( 1, '  Illegal problem number = %d\n', prob );
    error ( 'P00_DAT - Fatal error!' );
  end

  return
end

author: johannesgerer


text:
function [hpb] = W2hpb(W)
% Convert power from watts to boiler horsepower.
% Chad A. Greene 2012
hpb = W*0.000101942;

author: Sable


text:
function varargout=wqrs(varargin)
%
% wqrs(recordName,N,N0,signal,threshold,findJ,powerLineFrequency,resample)
%
%    Wrapper to WFDB WQRS:
%         http://www.physionet.org/physiotools/wag/wqrs-1.htm
%
%    Creates a WQRS annotation file  at the current MATLAB directory.
%    The annotation file will have the same name as the recorName file,
%    but followed with the *.wqrs suffix. Use RDANN to read the annoations
%    into MATLAB's workspace in order to read the sample QRS locations.
%
%    If recordName is the path to a record at PhysioNet's database, then
%    the annation files will be stored in a subdirectory with the same relative
%    path as recordName and under the current directory.
%
%
%    NOTE: In order to read the generated annotation file using RDANN, it is
%    necessary to have the WFDB record (*.hea and *.dat) files in the same
%    directory as the annotation file.
%
%
% CITING CREDIT: To credit this function, please cite the following paper in your work:
%
% Zong, W., G. B. Moody, and D. Jiang."A robust open-source algorithm to detect onset
% and duration of QRS complexes." Computers in Cardiology, 2003. IEEE, 2003.
%
%
% Required Parameters:
%
% recorName
%       String specifying the name of the record in the WFDB path or
%       in the current directory.
%
%
% Optional Parameters are:
%
% N
%       A 1x1 integer specifying the sample number at which to stop reading the
%       record file (default read all = N).
% N0
%       A 1x1 integer specifying the sample number at which to start reading the
%       annotion file (default 1 = begining of the record).
%
% signal
%       A 1x1 integer. Specify the signal to obtain the annotation (default
%       = 1, first signal).
%
% threshold
%       A 1x1 double. Specify the detection threshold in millivolts
%      (default = 100). Use higher values to reduce false detections, or lower values to
%       reduce the number of missed beats.
%
% findJ
%       A 1x1 boolean. Find and annotate J-points (QRS ends) as well as QRS onsets.
%       Default is false.
%
%
% powerLineFrequency
%       A 1x1 double. Specify the power line main frequency in Hz (default=
%       60 Hz), wqrs will apply a notch filter of the specified frequency
%       to the input signal before length-transforming it.
%
% resample
%       A 1x1 boolean. Resample the input at 120 Hz if the power line frequency
%       is 60 Hz, or at 150 Hz otherwise (default= false : do not resample).
%
%
%
% Source code written by Wei Zong and George B. Moody.
%
% MATLAB wrapper written by Ikaro Silva, 2013
% Last Modified: -
% Version 1.0
% Since 0.0.1
%
% %Example - Requires write permission to current directory
%wqrs('challenge/2013/set-a/a01');

%endOfHelp

javaWfdbExec=getWfdbClass('wqrs');

%Set default pararamter values
inputs={'recordName','N','N0','signal','threshold', ...
    'findJ','powerLineFrequency','resample'};
N=[];
N0=1;
signal=[]; %use application default
threshold=[];%use application default
findJ=[];
powerLineFrequency=[];
resample=[];
for n=1:nargin
    if(~isempty(varargin{n}))
        eval([inputs{n} '=varargin{n};']);
    end
end

N0=num2str(N0-1); %-1 is necessary because WFDB is 0 based indexed.
wfdb_argument={'-r',recordName,'-f',['s' N0]};

if(~isempty(N))
    wfdb_argument{end+1}='-t';
    %-1 is necessary because WFDB is 0 based indexed.
    wfdb_argument{end+1}=['s' num2str(N-1)];
end
if(~isempty(signal))
    wfdb_argument{end+1}='-s';
    %-1 is necessary because WFDB is 0 based indexed.
    wfdb_argument{end+1}=num2str(signal-1);
end
if(~isempty(threshold))
    wfdb_argument{end+1}='-m';
    %-1 is necessary because WFDB is 0 based indexed.
    wfdb_argument{end+1}=num2str(threshold-1);
end
if(~isempty(findJ) && findJ)
    wfdb_argument{end+1}='-j';
end

if(~isempty(powerLineFrequency))
    wfdb_argument{end+1}='-p';
    %-1 is necessary because WFDB is 0 based indexed.
    wfdb_argument{end+1}=num2str(powerLineFrequency);
end

if(~isempty(resample) && resample)
    wfdb_argument{end+1}='-R';
end

javaWfdbExec.execToStringList(wfdb_argument);

author: ikarosilva

repo: wfdb-app-toolbox


text:
function criterion = gbc_angle(q,CS,Dl,Dr,threshold,varargin)
%

d = angle(orientation(q(Dl),CS),orientation(q(Dr),CS));

if length(threshold) == 1

  criterion = d < threshold;

else

  criterion = 0.5 * ((d < threshold(1)) + (d < threshold(2)));

end

% now check whether the have a misorientation heigher or lower than a
% threshold
%m = inv(q(Dl)).*q(Dr);

%criterion = abs(dot(m,quaternion.id)) > cos(threshold/2);

%if any(~criterion)
%  qcs = quaternion(CS.properGroup);
%  criterion(~criterion) = max(abs(dot_outer(m(~criterion),qcs)),[],2) > cos(threshold/2);
%end

% o_Dl = orientation(q(Dl),CS,symmetry);
% o_Dr = orientation(q(Dr),CS,symmetry);
% criterion = dot(o_Dl,o_Dr) > cos(threshold/2);

author: mtex-toolbox

repo: mtex

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


text:
classdef TestGetOptimalNewCameraMatrix
    %TestGetOptimalNewCameraMatrix

    methods (Static)
        function test_1
            camMtx = [diag(rand(2,1)*100) rand(2,1)*50; 0 0 1];
            distCoeffs = ones(1,4)*1e-4;
            imgSz = [100 100];
            for alpha=0:0.25:1
                A = cv.getOptimalNewCameraMatrix(camMtx, distCoeffs, imgSz, ...
                    'Alpha',alpha);
                validateattributes(A, {class(camMtx)}, ...
                    {'2d', 'real', 'size',[3 3]});
            end
        end

        function test_2
            camMtx = single([diag(rand(2,1)*100) rand(2,1)*50; 0 0 1]);
            [A,validPixROI] = cv.getOptimalNewCameraMatrix(camMtx, ...
                zeros(1,4), [100 100], 'Alpha',0.5, ...
                'NewImageSize',[80 80], 'CenterPrincipalPoint',true);
            validateattributes(A, {class(camMtx)}, ...
                {'2d', 'real', 'size',[3 3]});
            validateattributes(validPixROI, {'numeric'}, ...
                {'vector', 'numel',4, 'integer'});
        end

        function test_error_argnum
            try
                cv.getOptimalNewCameraMatrix();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu


text:
function R = ProjectToSO3(mat)
% *** CHAPTER 3: RIGID-BODY MOTIONS ***
% Takes mat: A matrix near SO(3) to project to SO(3).
% Returns R representing the closest rotation matrix that is in SO(3).
% This function uses singular-value decomposition (see
% http://hades.mech.northwestern.edu/index.php/Modern_Robotics_Linear_Algebra_Review)
% and is only appropriate for matrices close to SO(3).
% Example Inputs:
%
% clear; clc;
% mat = [ 0.675, 0.150,  0.720;
%         0.370, 0.771, -0.511;
%        -0.630, 0.619,  0.472];
% R = ProjectToSO3(mat)
%
% Output:
% R =
%    0.6790    0.1489    0.7189
%    0.3732    0.7732   -0.5127
%   -0.6322    0.6164    0.4694

[U, S, V] = svd(mat);
R = U * V';
if det(R) < 0
    % In this case the result may be far from mat.
    R = [R(:, 1: 2); -R(:, 3)];
end
end

author: ShuoYangRobotics

repo: QuadrupedSim