text:
function test_num = c8mat_exp_test_num ( )

%*****************************************************************************80
%
%% C8MAT_EXP_TEST_NUM returns the number of matrix exponential tests.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer TEST_NUM, the number of tests.
%
  test_num = 3;

  return
end


text:
function [x,zo]=overlapadd(f,win,inc)
%OVERLAPADD join overlapping frames together X=(F,WIN,INC)
%
% Inputs:  F(NR,NW) contains the frames to be added together, one
%                   frame per row.
%          WIN(NW)  contains a window function to multiply each frame.
%                   WIN may be omitted to use a default rectangular window
%                   If processing the input in chunks, WIN should be replaced by
%                   ZI on the second and subsequent calls where ZI is the saved
%                   output state from the previous call.
%          INC      gives the time increment (in samples) between
%                   succesive frames [default = NW].
%
% Outputs: X(N,1) is the output signal. The number of output samples is N=NW+(NR-1)*INC.
%          ZO     Contains the saved state to allow a long signal
%                 to be processed in chunks. In this case X will contain only N=NR*INC
%                 output samples.
%
% Example of frame-based processing:
%          INC=20       													% set frame increment
%          NW=INC*2     													% oversample by a factor of 2 (4 is also often used)
%          S=cos((0:NW*7)*6*pi/NW);								% example input signal
%          W=sqrt(hamming(NW+1)); W(end)=[];      % sqrt hamming window of period NW
%          F=enframe(S,W,INC);               			% split into frames
%          ... process frames ...
%          X=overlapadd(F,W,INC);           			% reconstitute the time waveform (omit "X=" to plot waveform)

%	   Copyright (C) Mike Brookes 2009
%      Version: $Id: overlapadd.m,v 1.2 2009/06/08 16:21:49 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[nr,nf]=size(f);            % number of frames and frame length
if nargin<2
    win=nf;                 % default increment
end
if isstruct(win)
    w=win.w;
    if ~numel(w) && length(w)~=nf
        error('window length does not match frames size');
    end
    inc=win.inc;
    xx=win.xx;
else
    if nargin<3
        inc=nf;
    end
    if numel(win)==1 && win==fix(win) && nargin<3       % win has been omitted
        inc=win;
        w=[];
    else
        w=win(:).';
        if length(w)~=nf
            error('window length does not match frames size');
        end
        if all(w==1)
            w=[];
        end
    end
    xx=[];      % partial output from previous call is null
end
nb=ceil(nf/inc);        % number of overlap buffers
no=nf+(nr-1)*inc;       % buffer length
z=zeros(no,nb);                      % space for overlapped output speech
if numel(w)
    z(repmat(1:nf,nr,1)+repmat((0:nr-1)'*inc+rem((0:nr-1)',nb)*no,1,nf))=f.*repmat(w,nr,1);
else
    z(repmat(1:nf,nr,1)+repmat((0:nr-1)'*inc+rem((0:nr-1)',nb)*no,1,nf))=f;
end
x=sum(z,2);
if ~isempty(xx)
    x(1:length(xx))=x(1:length(xx))+xx;     % add on leftovers from previous call
end
if nargout>1            % check if we want to preserve the state
    mo=inc*nr;          % completed output samples
    if no<mo
        x(mo,1)=0;
        zo.xx=[];
    else
        zo.xx=x(mo+1:end);
        zo.w=w;
        zo.inc=inc;
        x=x(1:mo);
    end
elseif ~nargout
    if isempty(xx)
        k1=nf-inc;  % dubious samples at start
    else
        k1=0;
    end
    k2=nf-inc;      % dubious samples at end
    plot(1+(0:nr-1)*inc,x(1+(0:nr-1)*inc),'>r',nf+(0:nr-1)*inc,x(nf+(0:nr-1)*inc),'<r', ...
        1:k1+1,x(1:k1+1),':b',k1+1:no-k2,x(k1+1:end-k2),'-b',no-k2:no,x(no-k2:no),':b');
    xlabel('Sample Number');
    title(sprintf('%d frames of %d samples with %.0f%% overlap = %d samples',nr,nf,100*(1-inc/nf),no));
end

author: decouples


text:
function value = i4_gcd ( i, j )

%*****************************************************************************80
%
%% I4_GCD finds the greatest common divisor of I and J.
%
%  Discussion:
%
%    Only the absolute values of I and J are
%    considered, so that the result is always nonnegative.
%
%    If I or J is 0, I4_GCD is returned as max ( 1, abs ( I ), abs ( J ) ).
%
%    If I and J have no common factor, I4_GCD is returned as 1.
%
%    Otherwise, using the Euclidean algorithm, I4_GCD is the
%    largest common factor of I and J.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 July 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, J, two numbers whose greatest common divisor
%    is desired.
%
%    Output, integer VALUE, the greatest common divisor of I and J.
%
  value = 1;
%
%  Return immediately if either I or J is zero.
%
  if ( i == 0 )
    value = max ( 1, abs ( j ) );
    return
  elseif ( j == 0 )
    value = max ( 1, abs ( i ) );
    return
  end
%
%  Set IP to the larger of I and J, IQ to the smaller.
%  This way, we can alter IP and IQ as we go.
%
  ip = max ( abs ( i ), abs ( j ) );
  iq = min ( abs ( i ), abs ( j ) );
%
%  Carry out the Euclidean algorithm.
%
  while ( 1 )

    ir = mod ( ip, iq );

    if ( ir == 0 )
      break;
    end

    ip = iq;
    iq = ir;

  end

  value = iq;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function spm_reslice(P,flags)
% Rigid body reslicing of images
% FORMAT spm_reslice(P,flags)
%
% P      - matrix or cell array of filenames {one string per row}
%          All operations are performed relative to the first image.
%          ie. Coregistration is to the first image, and resampling
%          of images is into the space of the first image.
%
% flags  - a structure containing various options.  The fields are:
%
%         mask   - mask output images (true/false) [default: true]
%                  To avoid artifactual movement-related variance the
%                  realigned set of images can be internally masked, within
%                  the set (i.e. if any image has a zero value at a voxel
%                  than all images have zero values at that voxel). Zero
%                  values occur when regions 'outside' the image are moved
%                  'inside' the image during realignment.
%
%         mean   - write mean image (true/false) [default: true]
%                  The average of all the realigned scans is written to
%                  an image file with 'mean' prefix.
%
%         interp - the B-spline interpolation method [default: 1]
%                  Non-finite values result in Fourier interpolation. Note
%                  that Fourier interpolation only works for purely rigid
%                  body transformations. Voxel sizes must all be identical
%                  and isotropic.
%
%         which  - values of 0, 1 or 2 are allowed [default: 2]
%                  0   - don't create any resliced images.
%                        Useful if you only want a mean resliced image.
%                  1   - don't reslice the first image.
%                        The first image is not actually moved, so it may
%                        not be necessary to resample it.
%                  2   - reslice all the images.
%                  If which is a 2-element vector, flags.mean will be set
%                  to flags.which(2).
%
%         wrap   - three values of either 0 or 1, representing wrapping in
%                  each of the dimensions. For fMRI, [1 1 0] would be used.
%                  For PET, it would be [0 0 0]. [default: [0 0 0]]
%
%         prefix - prefix for resliced images [default: 'r']
%
%__________________________________________________________________________
%
% The spatially realigned images are written to the original subdirectory
% with the same (prefixed) filename. They are all aligned with the first.
%
% Inputs:
% A series of images conforming to SPM data format (see 'Data Format'). The
% relative displacement of the images is stored in their header.
%
% Outputs:
% The routine uses information in their headers and writes the realigned
% image files to the same subdirectory with a prefix.
%__________________________________________________________________________
% Copyright (C) 1999-2017 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_reslice.m 7141 2017-07-26 09:05:05Z guillaume $

%__________________________________________________________________________
%
% The headers of the images contain a 4x4 affine transformation matrix 'M',
% usually affected by the `realignment' and `coregistration' modules.
% What these matrices contain is a mapping from the voxel coordinates
% (x0,y0,z0) (where the first voxel is at coordinate (1,1,1)), to
% coordinates in millimeters (x1,y1,z1).
%
% x1 = M(1,1)*x0 + M(1,2)*y0 + M(1,3)*z0 + M(1,4)
% y1 = M(2,1)*x0 + M(2,2)*y0 + M(2,3)*z0 + M(2,4)
% z1 = M(3,1)*x0 + M(3,2)*y0 + M(3,3)*z0 + M(3,4)
%
% Assuming that image1 has a transformation matrix M1, and image2 has a
% transformation matrix M2, the mapping from image1 to image2 is: M2\M1
% (ie. from the coordinate system of image1 into millimeters, followed
% by a mapping from millimeters into the space of image2).
%
% Several spatial transformations (realignment, coregistration,
% normalisation) can be combined into a single operation (without the
% necessity of resampling the images several times).
%__________________________________________________________________________
%
% Refs:
%
% Friston KJ, Williams SR, Howard R Frackowiak RSJ and Turner R (1995)
% Movement-related effect in fMRI time-series.  Mag. Res. Med. 35:346-355
%
% W. F. Eddy, M. Fitzgerald and D. C. Noll (1996) Improved Image
% Registration by Using Fourier Interpolation. Mag. Res. Med. 36(6):923-931
%
% R. W. Cox and A. Jesmanowicz (1999)  Real-Time 3D Image Registration
% for Functional MRI. Mag. Res. Med. 42(6):1014-1018
%__________________________________________________________________________

SVNid = '$Rev: 7141 $';

%-Say hello
%--------------------------------------------------------------------------
SPMid = spm('FnBanner',mfilename,SVNid);

%-Parameters
%--------------------------------------------------------------------------
if ~nargin,      P = spm_select([1 Inf],'image'); end
if iscellstr(P), P = char(P);    end
if ischar(P),    P = spm_vol(P); end

def_flags        = spm_get_defaults('realign.write');
def_flags.prefix = 'r';
if nargin < 2
    flags = def_flags;
else
    fnms = fieldnames(def_flags);
    for i=1:length(fnms)
        if ~isfield(flags,fnms{i})
            flags.(fnms{i}) = def_flags.(fnms{i});
        end
    end
end

if numel(flags.which) == 2
    flags.mean  = flags.which(2);
    flags.which = flags.which(1);
elseif ~isfield(flags,'mean')
    flags.mean  = 1;
end

%-Reslice
%--------------------------------------------------------------------------
if isempty(P)
    warning('Nothing to do.');
else
    reslice_images(P,flags);
end

fprintf('%-40s: %30s\n','Completed',spm('time'))                        %-#

%==========================================================================
%-function reslice_images(P,flags)
%==========================================================================
function reslice_images(P,flags)
% Reslice images volume by volume
% FORMAT reslice_images(P,flags)
% See main function for a description of the input parameters

if ~isfinite(flags.interp) % Use Fourier method
    % Check for non-rigid transformations in the matrixes
    for i=1:numel(P)
        pp = P(1).mat\P(i).mat;
        if any(abs(svd(pp(1:3,1:3))-1)>1e-7)
            fprintf('\n  Zooms  or shears  appear to  be needed');
            fprintf('\n  (probably due to non-isotropic voxels).');
            fprintf('\n  These  can not yet be  done  using  the');
            fprintf('\n  Fourier reslicing method.  Switching to');
            fprintf('\n  7th degree B-spline interpolation instead.\n\n');
            flags.interp = 7;
            break
        end
    end
end

if flags.mask || flags.mean
    spm_progress_bar('Init',P(1).dim(3),'Computing available voxels','planes completed');
    x1    = repmat((1:P(1).dim(1))',1,P(1).dim(2));
    x2    = repmat( 1:P(1).dim(2)  ,P(1).dim(1),1);
    if flags.mean
        Count    = zeros(P(1).dim(1:3));
        Integral = zeros(P(1).dim(1:3));
    end
    if flags.mask, msk = cell(P(1).dim(3),1);  end
    for x3 = 1:P(1).dim(3)
        tmp = zeros(P(1).dim(1:2));
        for i = 1:numel(P)
            tmp = tmp + getmask(inv(P(1).mat\P(i).mat),x1,x2,x3,P(i).dim(1:3),flags.wrap);
        end
        if flags.mask, msk{x3} = find(tmp ~= numel(P)); end
        if flags.mean, Count(:,:,x3) = tmp; end
        spm_progress_bar('Set',x3);
    end
end

nread = numel(P);
if ~flags.mean
    if flags.which == 1, nread = nread - 1; end
    if flags.which == 0, nread = 0; end
end
spm_progress_bar('Init',nread,'Reslicing','volumes completed');

[x1,x2] = ndgrid(1:P(1).dim(1),1:P(1).dim(2));
nread   = 0;
d       = [flags.interp*[1 1 1]' flags.wrap(:)];

for i = 1:numel(P)

    if (i>1 && flags.which==1) || flags.which==2
        write_vol = 1;
    else
        write_vol = 0;
    end
    if write_vol || flags.mean
        read_vol = 1;
    else
        read_vol = 0;
    end

    if read_vol
        if ~isfinite(flags.interp)
            v = abs(kspace3d(spm_bsplinc(P(i),[0 0 0 ; 0 0 0]'),P(1).mat\P(i).mat));
            for x3 = 1:P(1).dim(3)
                if flags.mean
                    Integral(:,:,x3) = ...
                        Integral(:,:,x3) + ...
                        nan2zero(v(:,:,x3) .* ...
                        getmask(inv(P(1).mat\P(i).mat),x1,x2,x3,P(i).dim(1:3),flags.wrap));
                end
                if flags.mask
                    tmp = v(:,:,x3); tmp(msk{x3}) = NaN; v(:,:,x3) = tmp;
                end
            end
        else
            C = spm_bsplinc(P(i), d);
            v = zeros(P(1).dim);
            for x3 = 1:P(1).dim(3)
                [tmp,y1,y2,y3] = getmask(inv(P(1).mat\P(i).mat),x1,x2,x3,P(i).dim(1:3),flags.wrap);
                v(:,:,x3)      = spm_bsplins(C, y1,y2,y3, d);
                % v(~tmp)      = 0;

                if flags.mean
                    Integral(:,:,x3) = Integral(:,:,x3) + nan2zero(v(:,:,x3));
                end
                if flags.mask
                    tmp = v(:,:,x3); tmp(msk{x3}) = NaN; v(:,:,x3) = tmp;
                end
            end
        end
        if write_vol
            VO         = P(i);
            VO.fname   = spm_file(P(i).fname, 'prefix',flags.prefix);
            VO.dim     = P(1).dim(1:3);
            VO.dt      = P(i).dt;
            VO.pinfo   = P(i).pinfo;
            VO.mat     = P(1).mat;
            VO.descrip = 'spm - realigned';
            VO = spm_write_vol(VO,v);
        end

        nread = nread + 1;
    end
    spm_progress_bar('Set',nread);
end

if flags.mean
    % Write integral image (16 bit signed)
    %----------------------------------------------------------------------
    Integral    = Integral./Count;
    PO          = P(1);
    PO.fname    = spm_file(P(1).fname, 'prefix','mean');
    PO          = rmfield(PO,'pinfo');
    PO.pinfo    = [max(max(max(Integral)))/32767 0 0]';
    PO.n        = [1 1];
    PO.descrip  = 'spm - mean image';
    PO.dt       = [spm_type('int16') spm_platform('bigend')];
    spm_write_vol(PO,Integral);
end

spm_progress_bar('Clear');

%==========================================================================
%-function v = kspace3d(v,M)
%==========================================================================
function v = kspace3d(v,M)
% 3D rigid body transformation performed as shears in 1D Fourier space
% FORMAT v = kspace3d(v,M)
% v        - image stored as a 3D array
% M        - rigid body transformation matrix
%
% v        - transformed image
%
% References:
% R. W. Cox and A. Jesmanowicz (1999)
% Real-Time 3D Image Registration for Functional MRI
% Magnetic Resonance in Medicine 42(6):1014-1018
%
% W. F. Eddy, M. Fitzgerald and D. C. Noll (1996)
% Improved Image Registration by Using Fourier Interpolation
% Magnetic Resonance in Medicine 36(6):923-931

[S0,S1,S2,S3] = shear_decomp(M);

d  = [size(v) 1 1 1];
g = 2.^ceil(log2(d));
if any(g~=d)
    tmp = v;
    v   = zeros(g);
    v(1:d(1),1:d(2),1:d(3)) = tmp;
    clear tmp;
end

% XY-shear
tmp1 = -sqrt(-1)*2*pi*([0:((g(3)-1)/2) 0 (-g(3)/2+1):-1])/g(3);
for j=1:g(2)
    t        = reshape( exp((j*S3(3,2) + S3(3,1)*(1:g(1)) + S3(3,4)).'*tmp1) ,[g(1) 1 g(3)]);
    v(:,j,:) = real(ifft(fft(v(:,j,:),[],3).*t,[],3));
end

% XZ-shear
tmp1 = -sqrt(-1)*2*pi*([0:((g(2)-1)/2) 0 (-g(2)/2+1):-1])/g(2);
for k=1:g(3)
    t        = exp( (k*S2(2,3) + S2(2,1)*(1:g(1)) + S2(2,4)).'*tmp1);
    v(:,:,k) = real(ifft(fft(v(:,:,k),[],2).*t,[],2));
end

% YZ-shear
tmp1 = -sqrt(-1)*2*pi*([0:((g(1)-1)/2) 0 (-g(1)/2+1):-1])/g(1);
for k=1:g(3)
    t        = exp( tmp1.'*(k*S1(1,3) + S1(1,2)*(1:g(2)) + S1(1,4)));
    v(:,:,k) = real(ifft(fft(v(:,:,k),[],1).*t,[],1));
end

% XY-shear
tmp1 = -sqrt(-1)*2*pi*([0:((g(3)-1)/2) 0 (-g(3)/2+1):-1])/g(3);
for j=1:g(2)
    t        = reshape( exp( (j*S0(3,2) + S0(3,1)*(1:g(1)) + S0(3,4)).'*tmp1) ,[g(1) 1 g(3)]);
    v(:,j,:) = real(ifft(fft(v(:,j,:),[],3).*t,[],3));
end

if any(g~=d), v = v(1:d(1),1:d(2),1:d(3)); end

%==========================================================================
%-function [S0,S1,S2,S3] = shear_decomp(A)
%==========================================================================
function [S0,S1,S2,S3] = shear_decomp(A)
% Decompose rotation and translation matrix A into shears S0, S1, S2 and
% S3, such that A = S0*S1*S2*S3. The original procedure is documented in:
% R. W. Cox and A. Jesmanowicz (1999)
% Real-Time 3D Image Registration for Functional MRI
% Magnetic Resonance in Medicine 42(6):1014-1018

A0 = A(1:3,1:3);
if any(abs(svd(A0)-1)>1e-7), error('Cannot decompose matrix.'); end

t  = A0(2,3); if t==0, t=eps; end
a0 = pinv(A0([1 2],[2 3])')*[(A0(3,2)-(A0(2,2)-1)/t) (A0(3,3)-1)]';
S0 = [1 0 0; 0 1 0; a0(1) a0(2) 1];
A1 = S0\A0;  a1 = pinv(A1([2 3],[2 3])')*A1(1,[2 3])';  S1 = [1 a1(1) a1(2); 0 1 0; 0 0 1];
A2 = S1\A1;  a2 = pinv(A2([1 3],[1 3])')*A2(2,[1 3])';  S2 = [1 0 0; a2(1) 1 a2(2); 0 0 1];
A3 = S2\A2;  a3 = pinv(A3([1 2],[1 2])')*A3(3,[1 2])';  S3 = [1 0 0; 0 1 0; a3(1) a3(2) 1];

s3 = A(3,4)-a0(1)*A(1,4)-a0(2)*A(2,4);
s1 = A(1,4)-a1(1)*A(2,4);
s2 = A(2,4);
S0 = [[S0 [0  0 s3]'];[0 0 0 1]];
S1 = [[S1 [s1 0  0]'];[0 0 0 1]];
S2 = [[S2 [0 s2  0]'];[0 0 0 1]];
S3 = [[S3 [0  0  0]'];[0 0 0 1]];

%==========================================================================
%-function [Mask,y1,y2,y3] = getmask(M,x1,x2,x3,dim,wrp)
%==========================================================================
function [Mask,y1,y2,y3] = getmask(M,x1,x2,x3,dim,wrp)
tiny = 5e-2; % From spm_vol_utils.c
y1   = M(1,1)*x1+M(1,2)*x2+(M(1,3)*x3+M(1,4));
y2   = M(2,1)*x1+M(2,2)*x2+(M(2,3)*x3+M(2,4));
y3   = M(3,1)*x1+M(3,2)*x2+(M(3,3)*x3+M(3,4));
Mask = true(size(y1));
if ~wrp(1), Mask = Mask & (y1 >= (1-tiny) & y1 <= (dim(1)+tiny)); end
if ~wrp(2), Mask = Mask & (y2 >= (1-tiny) & y2 <= (dim(2)+tiny)); end
if ~wrp(3), Mask = Mask & (y3 >= (1-tiny) & y3 <= (dim(3)+tiny)); end

%==========================================================================
%-function vo = nan2zero(vi)
%==========================================================================
function vo = nan2zero(vi)
vo = vi;
vo(~isfinite(vo)) = 0;

author: spm


text:
function val = get(obj, varargin)
%GET  Returns the value of a field.
%   VAL = OBJ.GET(FIELD1, FIELD2, ...)
%
%   Example:
%       var_order = obj.get('var', 'order');
%
%   See also OPT_MODEL.

%   MP-Opt-Model
%   Copyright (c) 2008-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

val = obj;
for k = 1:length(varargin)
    if ischar(varargin{k})
        val = val.(varargin{k});
    else
        val = val(varargin{k});
    end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
classdef mme_reserve_zone < mp.mm_element

%   MATPOWER
%   Copyright (c) 2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%     properties
%     end

    methods
        function name = name(obj)
            name = 'reserve_zone';
        end

        function obj = add_constraints(obj, mm, nm, dm, mpopt)
            dme = obj.data_model_element(dm);
            rgen_dme = dm.elements.reserve_gen;
            igr = rgen_dme.gen; %% index of online gen for each online reserve gen

            Areq = dme.zones(:, igr);
            lreq = dme.req;

            mm.add_lin_constraint('Rreq', Areq, lreq, [], {'R'});
        end

        function obj = data_model_update(obj, mm, nm, dm, mpopt)
            dme = obj.data_model_element(dm);

            %% get reserve zone prices
            ll = mm.get_idx('lin');
            prc = mm.soln.lambda.mu_l(ll.i1.Rreq:ll.iN.Rreq) / dm.base_mva;

            %% update in the data model
            dme.tab.prc(dme.on) = prc;
        end
    end     %% methods
end         %% classdef


text:
%
% data = read_matrix(FID,M,N)
%
% This function read a two dimensional matrix from the file specified. The
% dimension of the matrix need to be specified. The one dimensional vector
% can be accepted as a row matrix or column matrix.
%
% Author: Xiao Xiong
% Created: 6 Feb, 2007
% Last modified: 6 Feb, 2007

function data = read_matrix(FID,M,N)
% M and N are the number of rows and columns respectively
data = zeros(M,N);
for i=1:M
    line = fgetl(FID);
    temp = textscan(line, '%n', N);
    data(i,:) = temp{1};
end

author: singaxiong

repo: SignalGraph


text:
function R2 = RodriguesConversion(R1)
% Rodrigues - converts a Rodrigues rotation vector to rotation matrix or vice versa.
%
% Usage:
%           R2 = Rodrigues(R1)
%
% Input:
%           R1 : Rodrigues rotation vector (3x1 or 1x3) or rotation matrix (3x3).
%
% Output:
%           R2 : rotation matrix (3x3) or Rodrigues rotation vector (3x1 or 1x3).
%
% This code follows the algorithm given by
% [1] R. Hartley and A. Zisserman "Multiple View Geometry in Computer Vision,"
%     Cambridge, pp.583-585, 2003.
%
% Kim, Daesik
% Intelligent Systems Research Center
% Sungkyunkwan Univ. (SKKU), South Korea
% E-mail  : daesik80@skku.edu
% Homepage: http://www.daesik80.com
%
% July 2008  - Original version.

[r,c] = size(R1);

%% Rodrigues Rotation Vector to Rotation Matrix
if ((r == 3) && (c == 1)) || ((r == 1) && (c == 3))
    wx = [  0   -R1(3)  R1(2);
        R1(3)   0   -R1(1);
        -R1(2)  R1(1)   0   ];

    R1_norm = sqrt(R1(1)^2 + R1(2)^2 + R1(3)^2);

    if (R1_norm < eps)
        R2 = eye(3);
    else
        R2 = eye(3) + sin(R1_norm)/R1_norm*wx + (1-cos(R1_norm))/R1_norm^2*wx^2;
    end

    %% Rotation Matrix to Rodrigues Rotation Vector
elseif (r == 3) && (c == 3)
    w_norm = acos((trace(R1)-1)/2);
    if (w_norm < eps)
        R2 = [0 0 0]';
    else
        R2 = 1/(2*sin(w_norm))*[R1(3,2)-R1(2,3);R1(1,3)-R1(3,1);R1(2,1)-R1(1,2)]*w_norm;
    end
end

author: tobycollins

repo: IPPE

save_path: github-repos/MATLAB/tobycollins-IPPE


text:
function beampattern(N,d,w,steer_angle)
% beampattern(N,d,w,steer_angle)
% where:
% N = NO. OF ELEMENTS
% d = INTER-ELEMENT SPACINGS B/W ELEMENTS
% w = WEIGHTINGS
% steer_angle = STEERING ANGLE FOR BEAM PATTERN
% By: Muhammad Fahad
% 2005-MS-E-EE-38
if nargin < 4, steer_angle = 0; end    % for no steering direction
if nargin < 3, w = 1/N*ones(1,N); end   % for uniform linear arrays
if nargin < 2, d = 1/2; end   % for standard linear array

n = (-(N-1)/2:(N-1)/2).';
theta = pi*(-1:0.001:1);
u = cos(theta);
vv = exp(j*2*pi*d*n*u);
theta_T = steer_angle/180*pi;
ut=cos(theta_T);
W = w.*exp(j*n*pi*cos(theta_T))';
B = W*vv;
B = 10*log10(abs(B).^2);
figure
polardb(theta,B,-40)
return

author: Sable

repo: mcbench-benchmarks


text:
function y = im2jpeg(x, quality)
%IM2JPEG Compresses an image using a JPEG approximation.
%   Y = IM2JPEG(X, QUALITY) compresses image X based on 8 x 8 DCT
%   transforms, coefficient quantization, and Huffman symbol
%   coding. Input QUALITY determines the amount of information that
%   is lost and compression achieved.  Y is an encoding structure
%   containing fields:
%
%      Y.size      Size of X
%      Y.numblocks Number of 8-by-8 encoded blocks
%      Y.quality   Quality factor (as percent)
%      Y.huffman   Huffman encoding structure, as returned by
%                  MAT2HUFF
%
%   See also JPEG2IM.

%   Copyright 2002-2006 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.5 $  $Date: 2006/07/15 20:44:34 $
%
%   Revised: 3/20/06 by R.Woods to correct an 'eob' coding problem, to
%   check the 'quality' input for <= 0, and to fix a warning when
%   struct y is created.

error(nargchk(1, 2, nargin));             % Check input arguments
if ndims(x) ~= 2 | ~isreal(x) | ~isnumeric(x) | ~isa(x, 'uint8')
   error('The input must be a UINT8 image.');
end
if nargin ==2 && quality <= 0
   error('Input parameter QUALITY must be greater than zero.');
end
if nargin < 2
   quality = 1;   % Default value for quality.
end

m = [16 11  10  16  24  40  51  61        % JPEG normalizing array
     12  12  14  19  26  58  60  55       % and zig-zag redordering
     14  13  16  24  40  57  69  56       % pattern.
     14  17  22  29  51  87  80  62
     18  22  37  56  68  109 103 77
     24  35  55  64  81  104 113 92
     49  64  78  87  103 121 120 101
     72  92  95  98  112 100 103 99] * quality;

order = [1 9  2  3  10 17 25 18 11 4  5  12 19 26 33  ...
        41 34 27 20 13 6  7  14 21 28 35 42 49 57 50  ...
        43 36 29 22 15 8  16 23 30 37 44 51 58 59 52  ...
        45 38 31 24 32 39 46 53 60 61 54 47 40 48 55  ...
        62 63 56 64];

[xm, xn] = size(x);                % Get input size.
x = double(x) - 128;               % Level shift input
t = dctmtx(8);                     % Compute 8 x 8 DCT matrix

% Compute DCTs of 8x8 blocks and quantize the coefficients.
y = blkproc(x, [8 8], 'P1 * x * P2', t, t');
y = blkproc(y, [8 8], 'round(x ./ P1)', m);

y = im2col(y, [8 8], 'distinct');  % Break 8x8 blocks into columns
xb = size(y, 2);                   % Get number of blocks
y = y(order, :);                   % Reorder column elements

eob = max(y(:)) + 1;               % Create end-of-block symbol
r = zeros(numel(y) + size(y, 2), 1);
count = 0;
for j = 1:xb                       % Process 1 block (col) at a time
   i = max(find(y(:, j)));         % Find last non-zero element
   if isempty(i)                   % No nonzero block values
      i = 0;
   end
   p = count + 1;
   q = p + i;
   r(p:q) = [y(1:i, j); eob];      % Truncate trailing 0's, add EOB,
   count = count + i + 1;          % and add to output vector
end

r((count + 1):end) = [];           % Delete unusued portion of r

y           = struct;
y.size      = uint16([xm xn]);
y.numblocks = uint16(xb);
y.quality   = uint16(quality * 100);
y.huffman   = mat2huff(r);

author: Ultrasty

repo: Digital-Image-Processing


text:
function SenFig = drawSenFig(SenFig, Sen, Raw, Obs, FigOpt)

% DRAWSENFIG  Redraw one sensor figure.
% 	DRAWSENFIG(SENFIG, SEN, RAW, OBS, FIGOPT) updates all the handles in
% 	the handles structure SENFIG to reflect the observations OBS taken by
% 	sensor SEN, together with the raw data RAW. SENFIG is one sensor figure
% 	structure created with CREATESENFIG.
%
%   See also CREATESENFIG, DRAWMAPFIG.

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

% Sensor type:
% ------------
switch Sen.type

    % camera pinhole
    % --------------
    case {'pinHole','pinHoleDepth'}

        % 1. Raw data visualisation
        switch Raw.type
            case 'simu'
                drawRawPnts(SenFig, Raw);
                drawRawLines(SenFig, Raw);
            case 'dump'
                set(SenFig.img,'cdata',Raw.data.img);
                drawRawPnts(SenFig, Raw);
                drawRawLines(SenFig, Raw);
            case 'img'
                set(SenFig.img,'cdata',Raw.data.img);
            otherwise
        end

        % 2. Process only visible landmarks:
        % a - first erase lmks that were drawn but are no longer visible
        vis   = [Obs(:).vis];
        drawn = SenFig.drawn;
        erase = drawn & ~vis;
        if any(erase)
            set(SenFig.measure(erase),   'vis', 'off');
            set(SenFig.mean(erase),      'vis', 'off');
            set(SenFig.ellipse(erase,:), 'vis', 'off');
            set(SenFig.label(erase),     'vis', 'off');
            SenFig.drawn(erase) = false;
        end

        % b - now draw only visible landmarks
        for lmk = find(vis)

            SenFig.drawn(lmk) = true;
            % Landmark type:
            % --------------
            switch Obs(lmk).ltype

                case {'eucPnt'}  % Euclidean point
                    colors = FigOpt.sensor.colors.defPnt;
                    drawObsPnt(SenFig, Obs(lmk), colors, FigOpt.sensor.showEllip);

                case {'idpPnt','hmgPnt','ahmPnt','fhmPnt'}  % IDP and HMG points
                    colors = FigOpt.sensor.colors.othPnt;
                    drawObsPnt(SenFig, Obs(lmk), colors, FigOpt.sensor.showEllip);

                case {'plkLin','aplLin','idpLin','hmgLin', 'ahmLin'}  % lines
                    colors = FigOpt.sensor.colors.othLin;
                    drawObsLin(SenFig, Obs(lmk), colors, FigOpt.sensor.showEllip);

                    % ADD HERE FOR NEW LANDMARK TYPE
                case {'newLandmark'}
                    % do something

                otherwise
                    % Print an error and exit
                    error('??? Unable to display landmark ''%s'' with sensor ''%s''.',Obs(lmk).ltype,Sen.type);

            end % and of the "switch" on sensor type

        end

    % --------------
    case {'omniCam'}

        % 1. Raw data visualisation
        if strcmp(Raw.type,'simu')
            drawRawPnts(SenFig, Raw);
        else
            %drawRawImage(SenFig, Raw);
            imshow(Raw.img,'InitialMagnification', 'fit', 'Parent', SenFig.axes);
        end

        % 2. Process only visible landmarks:
        % a - first erase lmks that were drawn but are no longer visible
        vis   = [Obs(:).vis];
        drawn = SenFig.drawn;
        erase = drawn & ~vis;
        if any(erase)
            set(SenFig.measure(erase),   'vis', 'off');
            set(SenFig.mean(erase),      'vis', 'off');
            set(SenFig.ellipse(erase,:), 'vis', 'off');
            set(SenFig.label(erase),     'vis', 'off');
            SenFig.drawn(erase) = false;
        end
        % b - now draw only visible landmarks
        for lmk = find(vis)

            SenFig.drawn(lmk) = true;

            % Landmark type:
            % --------------
            switch Obs(lmk).ltype

                case {'eucPnt'}  % Euclidean point
                    colors = FigOpt.sensor.colors.defPnt;
                    drawObsPnt(SenFig, Obs(lmk), colors, FigOpt.sensor.showEllip);

                case {'idpPnt','hmgPnt','ahmPnt','fhmPnt'}  % IDP and HMG points
                    colors = FigOpt.sensor.colors.othPnt;
                    drawObsPnt(SenFig, Obs(lmk), colors, FigOpt.sensor.showEllip);

                case {'plkLin','aplLin','idpLin','hmgLin', 'ahmLin'}  % Plucker line
                    colors = FigOpt.sensor.colors.othLin;
                    drawObsLin(SenFig, Obs(lmk), colors, FigOpt.sensor.showEllip);

                    % ADD HERE FOR NEW LANDMARK TYPE
                case {'newLandmark'}
                    % do something

                otherwise
                    % Print an error and exit
                    error('??? Unable to display landmark ''%s'' with sensor ''%s''.',Obs(lmk).ltype,Sen.type);

            end % and of the "switch" on sensor type

        end

        % ADD HERE FOR INITIALIZATION OF NEW SENSORS's FIGURES
        % case {'newSensor'}
        % do something

        % unknown
        % -------
    otherwise
        error('??? Unknown sensor type ''%s''.',Sen.type);
end

%     end

% end

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

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   SLAMTB is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function A = slaffinitymat(X, X2, nnparams, varargin)
%SLAFFINITYMAT Constructs an affinity matrix
%
% $ Syntax $
%   - A = slaffinitymat(X, [], nnparams, ...)
%   - A = slaffinitymat(X, X2, nnparams, ...)
%
% $ Arguments $
%   - X:        The sample matrix of the (source) nodes
%   - X2:       The sample matrix of the (target) nodes
%   - nnparams: The cell array of parameters for finding nearest neighbors
%               in the form of {method, ...}. Please refer to slfindnn
%               for details of specifying nnparams.
%   - A:        The constructed affinity matrix
%
% $ Description $
%   - A = slaffinitymat(X, X2, nnparams, ...) constructs an affinity
%     matrix to represent the pairwise affinity between neighboring
%     samples. The affinity between non-neighboring samples is set to
%     zero. You can indicate a self-affinity matrix (affinity among the
%     the set of samples) by setting X2 to [].
%     By default, the function uses heated kernel to compute the affinity
%     as explained below. In addition, you can use other formulas to
%     translate the distances to affinity by supplying your tfunctor in
%     the properties.
%     The following is a table of properties that you can specify:
%       \*
%       \t   The Properties of Affinity Matrix construction         \\
%       \h     name         &      description                      \\
%             'sparse'      & Whether to construct sparse matrix
%                             (default = true)                      \\
%             'kernel'      & The kernel to compute affinity         \\
%                             - 'heat': the heated kernel (default)
%                               it uses the following formula to translate
%                               the Euclidean distances into affinities:
%                                 a = exp(- d^2 / (2 *sigma^2))
%                               Here sigma^2 is set to
%                                 diffusion^2 * avg(d^2).
%                               you can set the value of diffusion in the
%                               properties.
%                             - 'simple':  the simple kernel
%                               just set the affinity of all neighboring
%                               samples to 1.
%             'diffusion'   & The diffusion distance relative to
%                             the average distance. (default = 1)    \\
%             'tfunctor'    & The function to transform the distance
%                             values to edge values. (default = [])  \\
%             'sym'         & whether to symmetrizes the graph
%                             (default = true)                       \\
%             'symmethod'   & The method used to symmetrization
%                             (default = [])                          \\
%             'excludeself' & Whether to exclude the edges connecting
%                             the same node. (default = false).
%       \*
%
% $ Arguments $
%   - The properties sym, symmethod and excludeself only take effect
%     when input X2 = [].
%
% $ Remarks $
%   - It wrapps slnngraph to provide a convenient way to compute
%     typical affinity matrix.
%
% $ History $
%   - Created by Dahua Lin, on Sep 12nd, 2006
%

%% parse input and prepare parameters

if nargin < 3
    raise_lackinput('slaffinitymat', 3);
end

opts.sparse = true;
opts.kernel = 'heat';
opts.diffusion = 1;
opts.tfunctor = [];
opts.sym = true;
opts.symmethod = [];
opts.excludeself = false;
opts = slparseprops(opts, varargin{:});

if isempty(X2)
    if opts.excludeself
        X2 = [];
    else
        X2 = X;
    end
else
    opts.sym = false;
end

if isempty(opts.tfunctor)
    switch opts.kernel
        case 'heat'
            tfunctor = {@internal_compaffinity, opts.diffusion};
        case 'simple'
            tfunctor = @(x) ones(size(x));
        otherwise
            error('sltoolbox:invalidarg', ...
                'Invalid kernel name: %s', opts.kernel);
    end
else
    tfunctor = opts.tfunctor;
end

%% main wrapper

A = slnngraph(X, X2, nnparams, ...
    'valtype', 'numeric', ...
    'sparse', opts.sparse, ...
    'tfunctor', tfunctor, ...
    'sym', opts.sym, ...
    'symmethod', opts.symmethod);

%% The internal function to compute affinity

function affvals = internal_compaffinity(dists, diffusion)

sqs = dists .* dists;
sqs = sqs(:);

avgsq = sum(sqs) / length(sqs);
s = (diffusion^2) * avgsq * 2;

affvals = exp(-sqs / s);

author: lmthang


text:
function b = isPointInEllipse(point, ellipse, varargin)
%ISPOINTINELLIPSE Check if a point is located inside a given ellipse.
%
%   B = isPointInEllipse(POINT, ELLIPSE)
%   Returns true if point is located inside the given ellipse.
%
%   B = isPointInEllipse(POINT, ELLIPSE, TOL)
%   Specifies the tolerance value
%
%   Example:
%   isPointInEllipse([1 0], [0 0 2 1 0])
%   ans =
%       1
%   isPointInEllipse([0 0], [0 0 2 1 0])
%   ans =
%       1
%   isPointInEllipse([1 1], [0 0 2 1 0])
%   ans =
%       0
%   isPointInEllipse([1 1], [0 0 2 1 30])
%   ans =
%       1
%
%   See also
%     ellipses2d, isPointInCircle
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2011-03-11
% Copyright 2011-2022 INRA - TPV URPOI - BIA IMASTE

% extract computation tolerance
tol = 1e-14;
if ~isempty(varargin)
    tol = varargin{1};
end

% compute ellipse to unit circle transform
rot = createRotation(-deg2rad(ellipse(5)));
sca = createScaling(1./ellipse(3:4));
trans = sca * rot;

% transform points to unit circle basis
pTrans = bsxfun(@minus, point, ellipse(:,1:2));
pTrans = transformPoint(pTrans, trans);

% test if distance to origin smaller than 1
b = sqrt(sum(power(pTrans, 2), 2)) - 1 <= tol;

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function plotFlybyXferTraj(hAxis, departBodyInfo, flybyBodyInfo, arrivalBodyInfo, xferOrbitIn, xferOrbitOut, gmuXfr)
%plotFlybyXferTraj Summary of this function goes here
%   Detailed explanation goes here
    axes(hAxis);
    cla(hAxis, 'reset');
    set(hAxis,'Color',[0 0 0]);
    set(hAxis,'XTick',[], 'XTickMode', 'manual');
    set(hAxis,'YTick',[], 'YTickMode', 'manual');
    set(hAxis,'ZTick',[], 'ZTickMode', 'manual');

    legendHs = [];
    legendStrs = {};

    hold on;
    h=plotBodyOrbit(departBodyInfo, 'r', gmuXfr);
    legendHs(end+1) = h;
    legendStrs{end+1} = [cap1stLetter(lower(departBodyInfo.name)), ' Orbit'];
    hold on;
    h=plotBodyOrbit(flybyBodyInfo, 'g', gmuXfr);
    legendHs(end+1) = h;
    legendStrs{end+1} = [cap1stLetter(lower(flybyBodyInfo.name)), ' Orbit'];
    hold on;
    h=plotBodyOrbit(arrivalBodyInfo, 'b', gmuXfr);
    legendHs(end+1) = h;
    legendStrs{end+1} = [cap1stLetter(lower(arrivalBodyInfo.name)), ' Orbit'];

    hold on;
    h=plotOrbit('w', xferOrbitIn(1), xferOrbitIn(2), xferOrbitIn(3), xferOrbitIn(4), xferOrbitIn(5), AngleZero2Pi(xferOrbitIn(6)), AngleZero2Pi(xferOrbitIn(7)), gmuXfr);
    legendHs(end+1) = h;
    legendStrs{end+1} = ['Phase 1 Xfer Orbit'];
    hold on;
    h=plotOrbit('w', xferOrbitOut(1), xferOrbitOut(2), xferOrbitOut(3), xferOrbitOut(4), xferOrbitOut(5), AngleZero2Pi(xferOrbitOut(6)), AngleZero2Pi(xferOrbitOut(7)), gmuXfr,[],[],[],'--');
    legendHs(end+1) = h;
    legendStrs{end+1} = ['Phase 2 Xfer Orbit'];

    hold on;
    plot3(0,0,0, 'yo', 'MarkerEdgeColor', 'y', 'MarkerFaceColor', 'y', 'MarkerSize', 10);

    [rVect, vVect] = getStateAtTime(departBodyInfo, 0, gmuXfr);
    [sma, ecc, inc, longAscNode, ArgPeri, TrueAnom] = getKeplerFromState(rVect,vVect,gmuXfr);
    departBodyOrbit = [sma, ecc, inc, longAscNode, ArgPeri, TrueAnom];

    [rVect, vVect] = getStateAtTime(flybyBodyInfo, 0, gmuXfr);
    [sma, ecc, inc, longAscNode, ArgPeri, TrueAnom] = getKeplerFromState(rVect,vVect,gmuXfr);
    flybyBodyOrbit = [sma, ecc, inc, longAscNode, ArgPeri, TrueAnom];

    [rVect, vVect] = getStateAtTime(arrivalBodyInfo, 0, gmuXfr);
    [sma, ecc, inc, longAscNode, ArgPeri, TrueAnom] = getKeplerFromState(rVect,vVect,gmuXfr);
    arriveBodyOrbit = [sma, ecc, inc, longAscNode, ArgPeri, TrueAnom];

    hold on;
    plotLineOfNodes(hAxis, departBodyOrbit, flybyBodyOrbit, xferOrbitIn, gmuXfr);
    hold on;
    plotLineOfNodes(hAxis, flybyBodyOrbit, arriveBodyOrbit, xferOrbitOut, gmuXfr);

    hold on;
    plotXfrOrbitApses(hAxis, xferOrbitIn, AngleZero2Pi(xferOrbitIn(6)), AngleZero2Pi(xferOrbitIn(7)), gmuXfr);
    hold on;
    plotXfrOrbitApses(hAxis, xferOrbitOut, AngleZero2Pi(xferOrbitOut(6)), AngleZero2Pi(xferOrbitOut(7)), gmuXfr);

    [rVect,~]=getStatefromKepler(xferOrbitIn(1), xferOrbitIn(2), xferOrbitIn(3), xferOrbitIn(4), xferOrbitIn(5), AngleZero2Pi(xferOrbitIn(6)), gmuXfr);
    hold on;
    plot3(rVect(1), rVect(2), rVect(3), 'ro', 'MarkerEdgeColor', 'r', 'MarkerFaceColor', 'r', 'MarkerSize', 10);

    [rVect,~]=getStatefromKepler(xferOrbitIn(1), xferOrbitIn(2), xferOrbitIn(3), xferOrbitIn(4), xferOrbitIn(5), AngleZero2Pi(xferOrbitIn(7)), gmuXfr);
    hold on;
    plot3(rVect(1), rVect(2), rVect(3), 'go', 'MarkerEdgeColor', 'g', 'MarkerFaceColor', 'g', 'MarkerSize', 10);

    [rVect,~]=getStatefromKepler(xferOrbitOut(1), xferOrbitOut(2), xferOrbitOut(3), xferOrbitOut(4), xferOrbitOut(5), AngleZero2Pi(xferOrbitOut(7)), gmuXfr);
    hold on;
    plot3(rVect(1), rVect(2), rVect(3), 'bo', 'MarkerEdgeColor', 'b', 'MarkerFaceColor', 'b', 'MarkerSize', 10);

    hold on;
    axis equal;
    view([0 89]);
    hold off;

    title('');
    legend(legendHs, legendStrs, 'Location', 'Best', 'EdgeColor', 'w', 'TextColor', 'w', 'LineWidth', 1.5);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% ------------------------------------------------------
% SwarmOps - Heuristic optimization for Matlab
% Copyright (C) 2003-2010 Magnus Erik Hvass Pedersen.
% Please see the file license.txt for license details.
% SwarmOps on the internet: http://www.Hvass-Labs.org/
% ------------------------------------------------------

% Enforce boundaries:
%     if x<lower then y=lower
%     elseif x>upper then y=upper
%     else y=x
% The implementation below works for arrays as well.
% Parameters:
%     x; position to be bounded.
%     lower; lower boundary.
%     upper; upper boundary;
% Returns:
%     y; bounded position.
function y = bound(x, lower, upper)
    y = min(upper, max(lower, x));
end

% ------------------------------------------------------


text:
%%% Note: run the 'DataGen/A_data_generation.m' to generate
%%% training data first.
clear;
rng('default')
addpath(genpath('./.'));
addpath(genpath('../.'));
addpath(genpath('../DataGen'));
addpath('matconvnet');
vl_setupnn;

%%%-------------------------------------------------------------------------
%%% Data generation
%%%-------------------------------------------------------------------------
A_data_generation;

%%%-------------------------------------------------------------------------
%%% Configuration
%%%-------------------------------------------------------------------------

opts.gpus             = [1 2]; %%% this code can only support multi-GPU!

opts.numSubBatches    = 1;
opts.bnormLearningRate= 0;

%%% solver
opts.solver           = 'Adam';
opts.numberImdb       = 1;

opts.derOutputs       = {'objective', 1} ;

ext               =  {'*.jpg','*.png','*.bmp'};

opts.gradientClipping = false; %%% set 'true' to prevent exploding gradients in the beginning.
opts.backPropDepth    = Inf;
%%%-------------------------------------------------------------------------
%%%   Initialize model and load data
%%%----------------------------------------------------

%%%  load data

global CurTask;
CurTask = 'Denoising'; %% 'Deblocking' and 'SISR'

opts.imdbDir          = '../DataGen/DN_PATCH192';
filepaths           =  [];
for i = 1 : length(ext)
    filepaths = [filepaths; dir(fullfile(opts.imdbDir, ext{i}))];
end
imdb.imdbPath = opts.imdbDir;
imdb.filepaths = filepaths;
imdb.images.set = ones(numel(filepaths),1);
fprintf('-----------------------------------------------------------\n');
fprintf('--------------------Training Number %d---------------------\n', numel(filepaths));
fprintf('-----------------------------------------------------------\n');
image = imread(fullfile(opts.imdbDir, filepaths(1).name));
imdb.patch_size = size(image, 1);

opts.learningRate     = [logspace(-3, -3, 15) logspace(-3.8, -4, 20) logspace(-4.5, -5, 10)];

%% for denosining, set simga to [15, 25, 50]; for SISR, set simga to [2, 3, 4], and for deblocking, set sigma to [10, 20, 30, 40]
opts.sigma            = 50;

opts.modelName        = ['MWCNN_GDSigma' num2str(opts.sigma)]; %%% model name
imdb.modelName        = opts.modelName;

opts.expDir      = fullfile('data', opts.modelName);
opts.batchSize        = 32*numel(opts.gpus);
net = net_wavelet_haart_24;

[~, ~] = cnn_train(net, imdb, ...
    'expDir', opts.expDir, ...
    'learningRate',opts.learningRate, ...
    'derOutputs',opts.derOutputs, ...
    'bnormLearningRate',opts.bnormLearningRate, ...
    'numSubBatches',opts.numSubBatches, ...
    'numberImdb',opts.numberImdb, ...
    'backPropDepth',opts.backPropDepth, ...
    'batchSize', opts.batchSize, ...
    'modelname', opts.modelName, ...
    'sigma', opts.sigma, ...
    'gpus',opts.gpus) ;

repo: MWCNN


text:
function pass = test_mean( pref )

% Grab some preferences
if ( nargin == 0 )
    pref = chebfunpref();
end
tol = pref.techPrefs.chebfuneps;

% Example 1:
f = ballfun(@(x,y,z)1);
I = mean(f,1);
g = spherefun(1);
pass(1) = ( norm(I-g) < tol );

% Example 2:
f = ballfun(@(x,y,z)2);
I = mean(f,2);
g = diskfun(2);
pass(2) = ( norm(I-g) < tol );

% Example 3:
f = ballfun(@(x,y,z)3);
I = mean(f,3);
g = diskfun(3);
pass(3) = ( norm(I-g) < tol );
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function test_ft_prepare_singleshell

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_headmodel_singleshell ft_prepare_headmodel ft_prepare_singleshell

% function to test ft_headmodel_singleshell. this function is called
% by ft_prepare_headmodel

% the function should work either on an input (segmented) mri, or it
% should have a description of the geometry in the input, or it should
% have a hdmfile (string) that specifies which file to read

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% get the data which is needed

% read in the segmented mri
datadir = dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/beamformer');
load(fullfile(datadir, 'segmentedmri.mat'));
mri = segmentedmri; clear segmentedmri;

% specify the file for the headshape
hdmfile  = dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/Subject01.shape');

% read in the headshape
shape = ft_read_headshape(hdmfile);

%%%%%%%%%%%%%%%%%%%%%
% do the computations

% with an MRI in the input
cfg         = [];
cfg.method  = 'singleshell';
vol1        = ft_prepare_headmodel(cfg, mri);

cfg         = [];
vol1b       = ft_prepare_singleshell(cfg, mri);

% the following needs to be done to be able to make the comparison
vol1 = rmfield(vol1, 'cfg');
vol1b = rmfield(vol1b,'cfg');
vol1bu=ft_convert_units(vol1b,vol1.unit);

% Reason for dashboard failure:
% e.g.  vol1.bnd.pnt(5,:) is different than vol1bu.bnd.pnt(5,:)
% e.g.  vol1.bnd.pnt(8,:) is different than vol1bu.bnd.pnt(8,:)
vol1.bnd    = rmfield(vol1.bnd, 'cfg');
vol1bu.bnd  = rmfield(vol1bu.bnd, 'cfg');

if ~isequal(vol1, vol1bu)
  error('ft_prepare_singleshell and ft_prepare_headmodel gave different outputs');
end

% with a filename in the input
cfg         = [];
cfg.method  = 'singleshell';
% cfg.headshape = hdmfile;
% vol2        = ft_prepare_headmodel(cfg);
vol2        = ft_prepare_headmodel(cfg, shape);

cfg         = [];
cfg.headshape = hdmfile;
vol2b       = ft_prepare_singleshell(cfg);

% the following needs to be done to be able to make the comparison
vol2  = rmfield(vol2, 'cfg');
vol2b = rmfield(vol2b,'cfg');

vol2.bnd  = rmfield(vol2.bnd, 'cfg');
vol2b.bnd = rmfield(vol2b.bnd,'cfg');

[ok, msg] = isalmostequal(vol2, vol2b, 'abstol', 1e-5, 'diffabs', 1);

if ~ok
  disp(msg);
  error('ft_prepare_singleshell and ft_prepare_headmodel gave different outputs');
end

% with a point cloud in the input
cfg         = [];
cfg.method  = 'singleshell';
cfg.geom    = shape;
vol3        = ft_prepare_headmodel(cfg);

cfg         = [];
cfg.headshape = shape;
vol3b       = ft_prepare_singleshell(cfg);

% the following needs to be done to be able to make the comparison
vol3  = rmfield(vol3, 'cfg');
vol3b = rmfield(vol3b,'cfg');

vol3.bnd  = rmfield(vol3.bnd, 'cfg');
vol3b.bnd = rmfield(vol3b.bnd,'cfg');

[ok, msg] = isalmostequal(vol3, vol3b, 'abstol', 1e-5, 'diffabs', 1);
if ~ok
  disp(msg);
  error('ft_prepare_singleshell and ft_prepare_headmodel gave different outputs');
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function a_cum = r8vec_cum0 ( n, a )

%*****************************************************************************80
%
%% R8VEC_CUM0 computes the cumulutive sums of an R8VEC.
%
%  Example:
%
%    Input:
%
%      A = (/ 1.0, 2.0, 3.0, 4.0 /)
%
%    Output:
%
%      A_CUM = (/ 0.0, 1.0, 3.0, 6.0, 10.0 /)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, real A(N), the vector to be summed.
%
%    Output, real A_CUM(1:N+1), the cumulative sums.
%
  a_cum(1) = 0.0;

  for i = 1 : n
    a_cum(i+1) = a_cum(i) + a(i);
  end

  return
end


text:
function truefalse = ch_is_format_code ( c )

%*****************************************************************************80
%
%% CH_IS_FORMAT_CODE returns TRUE if a character is a FORTRAN format code.
%
%  Discussion:
%
%    The format codes accepted here are not the only legal format
%    codes in FORTRAN90.  However, they are more than sufficient
%    for my needs!
%
%  Table:
%
%    A  Character
%    B  Binary digits
%    D  Real number, exponential representation
%    E  Real number, exponential representation
%    F  Real number, fixed point
%    G  General format
%    I  Integer
%    L  Logical variable
%    O  Octal digits
%    Z  Hexadecimal digits
%    *  Free format
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 November 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character C, the character to be analyzed.
%
%    Output, logical TRUEFALSE, is TRUE if C is a FORTRAN format code
%     and FALSE otherwise.
%
      if ( ch_eqi ( c, 'A' ) )
    truefalse = 1;
  elseif ( ch_eqi ( c, 'B' ) )
    truefalse = 1;
  elseif ( ch_eqi ( c, 'D' ) )
    truefalse = 1;
  elseif ( ch_eqi ( c, 'E' ) )
    truefalse = 1;
  elseif ( ch_eqi ( c, 'F' ) )
    truefalse = 1;
  elseif ( ch_eqi ( c, 'G' ) )
    truefalse = 1;
  elseif ( ch_eqi ( c, 'I' ) )
    truefalse = 1;
  elseif ( ch_eqi ( c, 'L' ) )
    truefalse = 1;
  elseif ( ch_eqi ( c, 'O' ) )
    truefalse = 1;
  elseif ( ch_eqi ( c, 'Z' ) )
    truefalse = 1;
  elseif ( c == '*' )
    truefalse = 1;
  else
    truefalse = 0;
  end

  return
end


text:
%{
    Copyright (C) Cambridge Electronic Design Limited 2014
    Author: James Thompson
    Web: www.ced.co.uk email: james@ced.co.uk, softhelp@ced.co.uk

    This file is part of CEDS64ML, a MATLAB interface to the SON64 library.

    CEDS64ML 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.

    CEDS64ML 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 CEDS64ML.  If not, see <http://www.gnu.org/licenses/>.
%}

function [ i64Div ] = CEDS64ChanDiv( fhand, iChan )
%CEDS64CHANDIV The sample interval in Ticks
%   [ i64Div ] = CEDS64ChanDiv( fhand, iChan )
%   Inputs
%   fhand - Integer file handle
%   iChan - Channel number
%   Outputs
%   i64Div - The sample interval in Ticks or a negative error code
if(nargin == 2)
    i64Div = calllib('ceds64int', 'S64ChanDivide', fhand, iChan);
else
    i64Div = -22;
end
end

author: brainstorm-tools


text:
function write_coveralls_json(obj,output_fn)
    monitor=obj.monitor;
    notify(monitor,sprintf('Writing JSON file to %s', output_fn));

    json=get_coverage_json(obj);
    fid=fopen(output_fn,'w');
    cleaner=onCleanup(@()fclose(fid));
    fprintf(fid,'%s',json);

    notify(monitor,sprintf('Completed writing JSON file to %s', output_fn));

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [Ep M] = spm_induced_optimise_parameters(PARAMS)
% Demo routine that optimises free parameters
%==========================================================================
%
% This exemplar routine illustrates how one can adjust or tune prior
% parameter expectations to produce desired spectral responses as specified
% by the complex eigenvalue spectrum - or a reduced form that considers a
% small number of complex values (roots).
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_induced_optimise_parameters.m 6856 2016-08-10 17:55:05Z karl $

% Parameters to optimise
%--------------------------------------------------------------------------
if ~nargin, PARAMS = {'G','T','L'}; end

% spectral specification
%==========================================================================
J      = [3 7];                % indices of hidden states producing outputs
Np     = length(J);

% Target spectrum - gamma, beta and alpha
%--------------------------------------------------------------------------
Hz     = 1:128;
s(:,1) = -[128 64 32 64]'   + 1j*2*pi*[4 12 48 64]';
s(:,2) = -[128 64 256 256]' + 1j*2*pi*[8 16 48 64]';

% Model specification
%==========================================================================
Nc    = 1;
Ns    = 1;
options.spatial  = 'LFP';
options.model    = 'TFM';
M.dipfit.model = options.model;
M.dipfit.type  = options.spatial;
M.dipfit.Nc    = Nc;
M.dipfit.Ns    = Ns;
M.J            = J;
M.Hz           = Hz;

% get priors
%--------------------------------------------------------------------------
[pE pC] = spm_dcm_neural_priors({0 0 0},{},1,options.model);
[pE pC] = spm_L_priors(M.dipfit,pE,pC);
[pE pC] = spm_ssr_priors(pE,pC);
[x,f]   = spm_dcm_x_neural(pE,options.model);

% suppress measurement noise
%--------------------------------------------------------------------------
pE.a(2) =      - 2;
pE.b    = pE.b - 16;
pE.c    = pE.c - 16;

% a target data
%--------------------------------------------------------------------------
Gu    = spm_csd_mtf_gu(pE,M.Hz);
for i = 1:Np
    csd(:,i,i) = 512*full(Gu.*sum(spm_s2csd(s(:,i),Hz),2));
end

% orders and model
%==========================================================================
nx      = length(spm_vec(x));
nu      = Ns;
u       = sparse(1,nu);

% fix priors if a subset of parameters are specified
%--------------------------------------------------------------------------
pV    = spm_vec(pC);
V     = pV - pV;
for i = 1:length(PARAMS)
    V(spm_fieldindices(pE,PARAMS{i})) = 1;
end
pC    = spm_unvec(V.*pV,pC);

% create LFP model
%--------------------------------------------------------------------------
M.IS   = 'spm_csd_mtf';
M.FS   = 'spm_diag_array';
M.g    = @(x,u,P,M) P.L*x(M.J(:));

M.f    = f;
M.x    = x;
M.n    = nx;
M.pE   = pE;
M.pC   = pC;
M.hE   = 8;
M.hC   = 1/128;
M.m    = nu;
M.l    = Np;

% solve for steady state
%--------------------------------------------------------------------------
M.x    = spm_dcm_neural_x(pE,M);
M.u    = u;
M.Nmax = 32;

% Optimisation: Target (Y)
%==========================================================================
[Ep Cp] = spm_nlsi_GN(M,[],{csd});

% Characterise contributions of parameters to spectral representation
%==========================================================================

% Show results with current (prior) parameters
%--------------------------------------------------------------------------
spm_figure('GetWin','Spectral responses'); clf

% and plot spectra
%--------------------------------------------------------------------------
subplot(2,2,1)
Gp    = spm_csd_mtf(pE,M);
Gq    = spm_csd_mtf(Ep,M);
plot(Hz,abs(spm_diag_array(Gp{1})),'r'), hold on
plot(Hz,abs(spm_diag_array(Gq{1})),'b'),  hold on
plot(Hz,abs(spm_diag_array(csd)),  '--'),  hold off

title({'Spectral responses';'Before (red) and after (blue)'},'FontSize',16)
xlabel('Frequency')
ylabel('Spectral density')
axis square

% Show results with optimised parameters
%--------------------------------------------------------------------------
subplot(2,2,2)
sp  = spm_ssm2s(pE,M);
sq  = spm_ssm2s(Ep,M);
gp  = spm_s2csd(sp,Hz);
gq  = spm_s2csd(sq,Hz);
plot(Hz,gp,'r',Hz,gq,'b')

title({'Eigenmodes';'Before (red) and after (blue)'},'FontSize',16)
xlabel('Frequency')
ylabel('Spectral density of modes')
axis square

% Show old and new priors
%==========================================================================

% change in parameters (and conditional confidence)
%--------------------------------------------------------------------------
E     = spm_vec(Ep) - spm_vec(pE);
C     = diag(Cp);

% eliminate an interesting parameters and sought on basis of contribution
%--------------------------------------------------------------------------
[c,j] = sort(C,'ascend');
j     = j(find(abs(E(j)) > exp(-8) & C(j) > 0));

subplot(4,1,3)
spm_plot_ci(E(j),C(j))
title('Posterior updates','FontSize',16)
xlabel('Free parameters')
ylabel('Value')
set(gca,'XTick',1:length(j))
set(gca,'XTickLabel',spm_fieldindices(pE,j))

subplot(4,1,4)
bar(-log(C(j,:)))
title('Contribution (log precision)','FontSize',16)
xlabel('Free parameters')
ylabel('log precision')
set(gca,'XTick',1:length(j))
set(gca,'XLim',[0 (length(j) + 1)])
set(gca,'XTickLabel',spm_fieldindices(pE,j))

return

% examine Jacobian
%==========================================================================
spm_figure('GetWin','Jacobian'); clf

% get transfer functions and Jacobian
%--------------------------------------------------------------------------
[S,K,s,w,t,dfdx] = spm_dcm_mtf(Ep,M,[]);

% and plot
%--------------------------------------------------------------------------
subplot(2,2,1)
imagesc(dfdx)
title('Jacobian','FontSize',16)
xlabel('hidden state')
ylabel('hidden state')
axis square


text:
%Simply returns the next point which is closest to the startpoint and as close as possible on the line
function [x, dist, dir, searchPoints] = getNextPoint(searchPoints, x0, dir0, X, Y, MAXDIST, MAXANGLE, MAXDIST2LINE, doshow)

    %REturn values in case of its not working
    x = -1;
    dir = [0;0];
    dist = -1;
    %Returns the index of points which are close to the line
    [dl, ix_l,angles] = getClosestPointToLine(searchPoints, x0, dir0, 32);

    %for i=1:length(dl)
    %    dl(i)
    %end
    %Only take those points that are less than MAXPIXEL away from the line
    ix = find(dl < MAXDIST2LINE);
    ix_l = ix_l(ix);

    %Returns the index of points which are close to the start point
    [ix_n, d] = getClosestNeighbours(searchPoints(:, ix_l), x0, 16); %x0 important

    %Make a new index with closest points that are close to the line!
    ix = ix_l(ix_n);
    if(numel(ix)>0)
        d = norm(searchPoints(2:3, ix(1))-x0); %x0 important
        if(d < MAXDIST)
            %Check that there ar epoints, otherwise return a zero value
            if(length(ix) > 0 )
                %Take the best match
                x = searchPoints(2:3,ix(1));
                %Recompute direction
                dir = x-x0(1:2);  %x0 faux
                dist = norm(dir);
                dir = dir/dist;
                %Get the angle
                theta = acos(dir'*dir0)/(norm(dir)*norm(dir0))/pi*180;

%                 if(theta<85 | theta>105)
                if(theta<MAXANGLE)
                    %Mark it as used and store the point correspondance
                    searchPoints(1,ix(1)) = 0;
                    searchPoints(4:5, ix(1)) = [X;Y];
                    if(doshow)
                        text(searchPoints(2,ix(1))+5,searchPoints(3,ix(1))+5,[num2str(searchPoints(4,ix(1))) ',' num2str(searchPoints(5,ix(1)))],'Color','m')
                    end
                else
                    x = -1;
                end
            end
        end
    end

%
% %Simply returns the next point which is closest to the startpoint and as close as possible on the line
% function [x, dist, dir, searchPoints] = getNextPoint(searchPoints, x0, dir0, X, Y, MAXDIST, MAXANGLE, MAXDIST2LINE, doshow)
%
%     %REturn values in case of its not working
%     x = -1;
%     dir = [0;0];
%     dist = -1;
%     %Returns the index of points which are close to the line
%     [dl, ix_l,angles] = getClosestPointToLine(searchPoints, x0, dir0, 8);
%     %Only take those points that are less than MAXPIXEL away from the line
%     ix = find(dl < MAXDIST2LINE);
%     ix_l = ix_l(ix);
%
%     %Returns the index of points which are close to the start point
%     [ix_n, d] = getClosestNeighbours(searchPoints(:, ix_l), x0, 16);
%
%     %Make a new index with closest points that are close to the line!
%     ix = ix_l(ix_n);
%     if(numel(ix)>0)
%         d = norm(searchPoints(2:3, ix(1))-x0);
%         if(d < MAXDIST)
%             %Check that there ar epoints, otherwise return a zero value
%             if(length(ix) > 0 )
%                 %Take the best match
%                 x = searchPoints(2:3,ix(1));
%                 %Recompute direction
%                 dir = x-x0(1:2);
%                 dist = norm(dir);
%                 dir = dir/dist;
%                 %Get the angle
%                 theta = acos(dir'*dir0)/(norm(dir)*norm(dir0))/pi*180;
%
% %                 if(theta<85 | theta>105)
%                 if(theta<MAXANGLE)
%                     %Mark it as used and store the point correspondance
%                     searchPoints(1,ix(1)) = 0;
%                     searchPoints(4:5, ix(1)) = [X;Y];
%                     if(doshow)
%                         text(searchPoints(2,ix(1))+5,searchPoints(3,ix(1))+5,[num2str(searchPoints(4,ix(1))) ',' num2str(searchPoints(5,ix(1)))],'Color','m')
%                     end
%                 else
%                     x = -1;
%                 end
%             end
%         end
%     end
%
%
%
%
%
%
%

author: christianwengert


text:
function [mu, s2, F, EffIdx, NffIdx] = GPR_PredictFun(hyp, data, xs, ys, tData, tDist)
if size(data, 2) ~= 2
    data = data';
end
if size(data, 2) ~= 2
    error('Input for GPR_PredictFun is wrong!');
end
if isrow(xs)
    xs = xs';
end
if isrow(ys)
    ys = ys';
end
x = data(:, 1);
y = data(:, 2);
[mu s2] = myGPRFun(hyp, x, y, xs);
F = [mu+tData*sqrt(s2); flipdim(mu-tData*sqrt(s2),1)];
if ~isempty(ys)
    Ratio = (mu - ys)./sqrt(s2) / tData;
    %%%%%%%%% update idx.
    EffIdx = find( abs(Ratio) <= 1 & sqrt(s2) <= tDist );
    NffIdx = find( ~ismember( 1:1:length(xs), EffIdx) );
    if isrow(NffIdx)
        NffIdx = NffIdx';
    end
else
    EffIdx = [];
    NffIdx = [];
end
end

author: DrGabor

repo: LiDAR

save_path: github-repos/MATLAB/DrGabor-LiDAR


text:
function r8plu_inverse_test ( )

%*****************************************************************************80
%
%% R8PLU_INVERSE_TEST tests R8PLU_INVERSE;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  n = 5;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8PLU_INVERSE_TEST\n' );
  fprintf ( 1, '  R8PLU_INVERSE determines the inverse of a matrix\n' );
  fprintf ( 1, '  from its compressed PLU factors.\n' );
  fprintf ( 1, '  Using initial random number seed = %d\n', seed );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );
%
%  Set the matrix.
%
  [ a, seed ] = r8mat_uniform_01 ( n, n, seed );

  r8mat_print ( n, n, a, '  The matrix A:' );
%
%  Factor the matrix.
%
  [ pivot, lu, info ] = r8mat_to_r8plu ( n, a );

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST0852 - Fatal error!\n' );
    fprintf ( 1, '  R8MAT_TO_R8PLU declares the matrix is singular!\n' );
    fprintf ( 1, '  The value of INFO is %d\n', info );
    return
  end
%
%  Compute the inverse.
%
  b = r8plu_inverse ( n, pivot, lu );

  r8mat_print ( n, n, b, '  The inverse matrix B:' );
%
%  Compute the product C = A * B.
%
  c(1:n,1:n) = a(1:n,1:n) * b(1:n,1:n);

  r8mat_print ( n, n, c, '  The product matrix C = A * B:' );

  return
end

author: johannesgerer


text:
classdef BSPTestsBase < matlab.unittest.TestCase
    properties(TestParameter)
        % Pull in board permutations
        configs = board_variants;
        ignored_builds = {'AnalogDevices.adrv9361z7035.ccbox_lvds.modem.plugin_board'};
        SynthesizeDesign = {false};
    end

    properties
        Count = 0;
        TotalTests = 0;
        Folder = pwd;
    end

    methods(TestClassSetup)
        function disableWarnings(~)
            warning('off','hdlcommon:hdlcommon:InterfaceNotAssigned');
        end
        function testCount(testCase)
            testCase.TotalTests = length(testCase.configs);
            CountS = 0;
            save('tc.mat','CountS');
        end
        function setRuntimeFolder(testCase)
            testCase.Folder = tempname(pwd);
        end
    end

    methods(TestClassTeardown)
        function enableWarnings(~)
            warning('on','hdlcommon:hdlcommon:InterfaceNotAssigned');
        end
        function collectLogs(~)
            if ~exist([pwd,'/logs'],'dir')
                mkdir('logs');
            end
            system('cp *.log logs/');
        end
    end

    methods(TestMethodSetup)
        function loadTestCount(testCase)
            l = load('tc.mat');
            CountS = l.CountS + 1;
            testCase.Count = CountS;
            save('tc.mat','CountS');
        end
    end

    methods

        function CollectLogs(testCase,cfgb)
            disp('Log collector called');
            rdn = strrep(cfgb.ReferenceDesignName,'/','_');
            rdn = strrep(rdn,'(','');
            rdn = strrep(rdn,')','');
            system(join(["find '",testCase.Folder,"' -name 'workflow_task_VivadoIPPackager.log' | xargs -I '{}' cp {} ."],''));
            if exist('workflow_task_VivadoIPPackager.log','file')
                disp('Found workflow_task_VivadoIPPackager... copying');
                movefile('workflow_task_VivadoIPPackager.log',[rdn,'_VivadoIPPackager_',cfgb.mode,'.log']);
            end
            system(join(["find '",testCase.Folder,"' -name 'workflow_task_CreateProject.log' | xargs -I '{}' cp {} ."],''));
            if exist('workflow_task_CreateProject.log','file')
                disp('Found workflow_task_CreateProject... copying');
                movefile('workflow_task_CreateProject.log',[rdn,'_CreateProject_',cfgb.mode,'.log']);
            end
            system(join(["find '",testCase.Folder,"' -name 'workflow_task_BuildFPGABitstream.log' | xargs -I '{}' cp {} ."],''));
            if exist('workflow_task_BuildFPGABitstream.log','file')
                disp('Found workflow_task_BuildFPGABitstream... copying');
                movefile('workflow_task_BuildFPGABitstream.log',[rdn,'_BuildFPGABitstream_',cfgb.mode,'.log']);
            end
        end

        function cfg = ADRV9361_Variants(~,s)

            variants = {...
                'ccbob_cmos','ccbob_lvds',...
                'ccbox_lvds','ccfmc_lvds'};
            cfg = {};
            s = strjoin(s(1:end-2),'.');
            h1 = str2func([s,'.common.plugin_board']);h1 = h1();

            for k = 1:length(variants)

                mode = 'rx';
                h2 = str2func([s,'.',variants{k},'.plugin_rd_rx']);h2 = h2();
                ReferenceDesignName = h2.ReferenceDesignName;
                vivado_version = h2.SupportedToolVersion{:};
                cfg1 = struct('Board',h1,...
                    'ReferenceDesignName',ReferenceDesignName,...
                    'vivado_version',vivado_version,'mode',mode);

                mode = 'tx';
                h2 = str2func([s,'.',variants{k},'.plugin_rd_tx']);h2 = h2();
                ReferenceDesignName = h2.ReferenceDesignName;
                vivado_version = h2.SupportedToolVersion{:};
                cfg2 = struct('Board',h1,...
                    'ReferenceDesignName',ReferenceDesignName,...
                    'vivado_version',vivado_version,'mode',mode);

                mode = 'rx_tx';
                h2 = str2func([s,'.',variants{k},'.plugin_rd_rxtx']);h2 = h2();
                ReferenceDesignName = h2.ReferenceDesignName;
                vivado_version = h2.SupportedToolVersion{:};
                cfg3 = struct('Board',h1,...
                    'ReferenceDesignName',ReferenceDesignName,...
                    'vivado_version',vivado_version,'mode',mode);
                cfg = [cfg(:)',{cfg1},{cfg2},{cfg3}];

            end

        end

        function cfg = extractConfigs(~,config)
            s = strsplit(config,'.');
            modes = strsplit(s{end},'_');
            mode = modes{end};
            h1 = str2func(config);h1 = h1();

            if strcmp(s{2},'adrv9361z7035') && ~isempty(strfind(s{2},'modem'))
                assert(0);
            elseif strcmp(s{2},'adrv9361z7035') || ...
                    strcmp(s{2},'adrv9364z7020')
                h = str2func([strjoin(s(1:2),'.'),'.common.plugin_board']);
            else
                h = str2func([strjoin(s(1:end-1),'.'),'.plugin_board']);
            end
            board = h();

            ReferenceDesignName = h1.ReferenceDesignName;
            vivado_version = h1.SupportedToolVersion{:};
            cfg = struct('Board',board,...
                'ReferenceDesignName',ReferenceDesignName,...
                'vivado_version',vivado_version,'mode',mode);
        end

        function setVivadoPath(~,vivado)
            if ispc
                pathname = ['C:\Xilinx\Vivado\',vivado,'\bin\vivado.bat'];
            elseif isunix
                pathname = ['/opt/Xilinx/Vivado/',vivado,'/bin/vivado'];
            end
            assert(exist(pathname,'file')>0,'Correct version of Vivado is unavailable or in a non-standard location');
            hdlsetuptoolpath('ToolName', 'Xilinx Vivado', ...
                'ToolPath', pathname);
            pause(4);
        end
    end

    methods(Test)
        function testMain(testCase, configs, SynthesizeDesign)
            % Filter out ignored configurations
            if ismember(configs,testCase.ignored_builds)
                assumeFail(testCase);
            end
            % Extract board configuration
            cfgb = testCase.extractConfigs(configs);
            %             for cfg = cfgs
            if exist(testCase.Folder,'dir')
                rmdir(testCase.Folder,'s');
                pause(1);
            end
            % Set up vivado
            testCase.setVivadoPath(cfgb.vivado_version);
            % Build
            disp(repmat('/',1,80));
            disp(['Building: ',cfgb.Board.BoardName,' | ',cfgb.mode,...
                ' (',num2str(testCase.Count),' of ',num2str(testCase.TotalTests),')']);
            res = build_design(cfgb.Board,cfgb.ReferenceDesignName,...
                cfgb.vivado_version,cfgb.mode,cfgb.Board.BoardName,...
                SynthesizeDesign,testCase.Folder);
            % Check
            if isfield(res,'message') || isa(res,'MException')
                disp(['Build error: ', cfgb.ReferenceDesignName]);
                disp(res);
                disp(res.message);
                disp(res.stack);
                testCase.CollectLogs(cfgb);
                verifyEmpty(testCase,res,res.message);
            end
        end
    end
end

author: analogdevicesinc


text:
function [w,u] = biharmonicP3(node,elem,pde,bdFlag,option)
%   [w,u] = biharmonicP1(node,elem,pde,bdFlag) produces the mixed cubic finite element
%   approximation of the biharmonic equation, where w = laplace u
%   See also biharmonicP1, biharmonicP2, biharmonicP3.
%   Created by Jie Zhou.
%   Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if nargin<5, option = []; end
tic;
%% Construct Data Structure
[elem2dof,elem2edge,edge,fixedDof,freeDof] = dofP3(elem);
N = size(node,1);  NT = size(elem,1); NE=size(edge,1); Ndof = N+2*NE+NT;
%% Compute geometric quantities and gradient of local basis
[Dlambda,area] = gradbasis(node,elem);

%% Assemble stiffness matrix
% Since Dphi_i*Dphi_j is quadratic, numerical quadrature rule is used here
if ~isfield(option,'quadorder')
    option.quadorder = 6;   % default order
end
[lambda, weight] = quadpts(option.quadorder);
nQuad = size(lambda,1);
ii = zeros(55*NT,1); jj = zeros(55*NT,1); sA = zeros(55*NT,nQuad);sB = zeros(55*NT,nQuad);
index = 0;
for i = 1:10
    for j = i:10
        ii(index+1:index+NT) = double(elem2dof(:,i));
        jj(index+1:index+NT) = double(elem2dof(:,j));
        index = index + NT;
    end
end

for p = 1:nQuad
        % Dphi at quadrature points
    Dphi(:,:,1) = (27/2*lambda(p,1)*lambda(p,1)-9*lambda(p,1)+1).*Dlambda(:,:,1);
    Dphi(:,:,2) = (27/2*lambda(p,2)*lambda(p,2)-9*lambda(p,2)+1).*Dlambda(:,:,2);
    Dphi(:,:,3) = (27/2*lambda(p,3)*lambda(p,3)-9*lambda(p,3)+1).*Dlambda(:,:,3);
    Dphi(:,:,4) = 9/2*((3*lambda(p,2)*lambda(p,2)-lambda(p,2)).*Dlambda(:,:,3)+...
                   lambda(p,3)*(6*lambda(p,2)-1).*Dlambda(:,:,2));
    Dphi(:,:,5) = 9/2*((3*lambda(p,3)*lambda(p,3)-lambda(p,3)).*Dlambda(:,:,2)+...
                   lambda(p,2)*(6*lambda(p,3)-1).*Dlambda(:,:,3));
    Dphi(:,:,6) = 9/2*((3*lambda(p,3)*lambda(p,3)-lambda(p,3)).*Dlambda(:,:,1)+...
                   lambda(p,1)*(6*lambda(p,3)-1).*Dlambda(:,:,3));
    Dphi(:,:,7) = 9/2*((3*lambda(p,1)*lambda(p,1)-lambda(p,1)).*Dlambda(:,:,3)+...
                   lambda(p,3)*(6*lambda(p,1)-1).*Dlambda(:,:,1));
    Dphi(:,:,8) = 9/2*((3*lambda(p,1)*lambda(p,1)-lambda(p,1)).*Dlambda(:,:,2)+...
                   lambda(p,2)*(6*lambda(p,1)-1).*Dlambda(:,:,1));
    Dphi(:,:,9)  = 9/2*((3*lambda(p,2)*lambda(p,2)-lambda(p,2)).*Dlambda(:,:,1)+...
                   lambda(p,1)*(6*lambda(p,2)-1).*Dlambda(:,:,2));
    Dphi(:,:,10) = 27*(lambda(p,1)*lambda(p,2)*Dlambda(:,:,3)+lambda(p,1)*lambda(p,3)*Dlambda(:,:,2)+...
                   lambda(p,3)*lambda(p,2)*Dlambda(:,:,1));
    phi(:,1) = 0.5*(3*lambda(:,1)-1).*(3*lambda(:,1)-2).*lambda(:,1);
    phi(:,2) = 0.5*(3*lambda(:,2)-1).*(3*lambda(:,2)-2).*lambda(:,2);
    phi(:,3) = 0.5*(3*lambda(:,3)-1).*(3*lambda(:,3)-2).*lambda(:,3);
    phi(:,4) = 9/2*lambda(:,3).*lambda(:,2).*(3*lambda(:,2)-1);
    phi(:,5) = 9/2*lambda(:,3).*lambda(:,2).*(3*lambda(:,3)-1);
    phi(:,6) = 9/2*lambda(:,1).*lambda(:,3).*(3*lambda(:,3)-1);
    phi(:,7) = 9/2*lambda(:,1).*lambda(:,3).*(3*lambda(:,1)-1);
    phi(:,8) = 9/2*lambda(:,1).*lambda(:,2).*(3*lambda(:,1)-1);
    phi(:,9) = 9/2*lambda(:,1).*lambda(:,2).*(3*lambda(:,2)-1);
    phi(:,10) = 27*lambda(:,1).*lambda(:,2).*lambda(:,3);
    index = 0;
    for i = 1:10
    for j = i:10
        Bij = 0;
        Aij = 0;
                  Bij = Bij + weight(p)*dot(Dphi(:,:,i),Dphi(:,:,j),2);
                  Aij = Aij + weight(p)*dot(phi(p,i),phi(p,j),2);
        Bij = Bij.*area;
        Aij = Aij.*area;
        ii(index+1:index+NT) = double(elem2dof(:,i));
        jj(index+1:index+NT) = double(elem2dof(:,j));
        sB(index+1:index+NT,p) = Bij;
        sA(index+1:index+NT,p) = Aij;
        index = index + NT;
    end
    end
end

sA = sum(sA,2);
sB = sum(sB,2);
clear Aij Bij
diagIdx = (ii == jj);   upperIdx = ~diagIdx;
B = sparse(ii(diagIdx),jj(diagIdx),sB(diagIdx),Ndof,Ndof);
A = sparse(ii(diagIdx),jj(diagIdx),sA(diagIdx),Ndof,Ndof);
% A = spdiags(accumarray(ii(diagIdx),sA(diagIdx),[Ndof 1]),0,Ndof,Ndof);
BU = sparse(ii(upperIdx),jj(upperIdx),sB(upperIdx),Ndof,Ndof);
AU = sparse(ii(upperIdx),jj(upperIdx),sA(upperIdx),Ndof,Ndof);
B = B + BU + BU';
A = A + AU + AU';

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % subfunction Dphi
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     function s = Dphi(p,i) % gradient of basis phi
%     switch i
%         case 1
%             s =(27/2*lambda(p,1)*lambda(p,1)-9*lambda(p,1)+1).*Dlambda(:,:,1);
%         case 2
%             s =(27/2*lambda(p,2)*lambda(p,2)-9*lambda(p,2)+1).*Dlambda(:,:,2);
%         case 3
%             s =(27/2*lambda(p,3)*lambda(p,3)-9*lambda(p,3)+1).*Dlambda(:,:,3);
%          case 4
%             s =9/2*((3*lambda(p,2)*lambda(p,2)-lambda(p,2)).*Dlambda(:,:,3)+...
%                 lambda(p,3)*(6*lambda(p,2)-1).*Dlambda(:,:,2));
%          case 5
%             s =9/2*((3*lambda(p,3)*lambda(p,3)-lambda(p,3)).*Dlambda(:,:,2)+...
%                 lambda(p,2)*(6*lambda(p,3)-1).*Dlambda(:,:,3));
%
%          case 6
%             s =9/2*((3*lambda(p,3)*lambda(p,3)-lambda(p,3)).*Dlambda(:,:,1)+...
%                 lambda(p,1)*(6*lambda(p,3)-1).*Dlambda(:,:,3));
%          case 7
%             s =9/2*((3*lambda(p,1)*lambda(p,1)-lambda(p,1)).*Dlambda(:,:,3)+...
%                 lambda(p,3)*(6*lambda(p,1)-1).*Dlambda(:,:,1));
%          case 8
%             s =9/2*((3*lambda(p,1)*lambda(p,1)-lambda(p,1)).*Dlambda(:,:,2)+...
%                 lambda(p,2)*(6*lambda(p,1)-1).*Dlambda(:,:,1));
%
%           case 9
%             s =9/2*((3*lambda(p,2)*lambda(p,2)-lambda(p,2)).*Dlambda(:,:,1)+...
%                 lambda(p,1)*(6*lambda(p,2)-1).*Dlambda(:,:,2));
%          case 10
%             s = 27*(lambda(p,1)*lambda(p,2)*Dlambda(:,:,3)+lambda(p,1)*lambda(p,3)*Dlambda(:,:,2)+...
%             lambda(p,3)*lambda(p,2)*Dlambda(:,:,1));
%
%     end
%     end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % subfunction phi
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     function s = phi(p,i) % gradient of basis phi
%     switch i
%         case 1
%             s = 0.5*(3*lambda(p,1)-1).*(3*lambda(p,1)-2).*lambda(p,1);
%         case 2
%             s = 0.5*(3*lambda(p,2)-1).*(3*lambda(p,2)-2).*lambda(p,2);
%         case 3
%             s = 0.5*(3*lambda(p,3)-1).*(3*lambda(p,3)-2).*lambda(p,3);
%         case 4
%             s = 9/2*lambda(p,3).*lambda(p,2).*(3*lambda(p,2)-1);
%         case 5
%             s = 9/2*lambda(p,3).*lambda(p,2).*(3*lambda(p,3)-1);
%         case 6
%             s = 9/2*lambda(p,1).*lambda(p,3).*(3*lambda(p,3)-1);
%         case 7
%             s = 9/2*lambda(p,1).*lambda(p,3).*(3*lambda(p,1)-1);
%         case 8
%             s = 9/2*lambda(p,1).*lambda(p,2).*(3*lambda(p,1)-1);
%         case 9
%             s = 9/2*lambda(p,1).*lambda(p,2).*(3*lambda(p,2)-1);
%         case 10
%             s = 27*lambda(p,1).*lambda(p,2).*lambda(p,3);
%     end
%     end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Assemble right hand side by high order quadrature rule
% To reduce the effect of the error introduced by the numerical quadrature,
% the load term is computed using the 3rd order qudrature rule.
b = zeros(Ndof,1);
u = zeros(Ndof,1);
w = zeros(Ndof,1);

if ~isfield(option,'fquadorder')
    option.fquadorder = 6;   % default order
end
if ~isfield(pde,'f') || (isreal(pde.f) && (pde.f==0))
    pde.f = [];
end
if ~isempty(pde.f)
    % quadrature points in the barycentric coordinate
    [lambda,w] = quadpts(option.fquadorder);
    nQuad = size(lambda,1);
%     phi(:,1) = lambda(:,1).*(2*lambda(:,1)-1);
%     phi(:,2) = lambda(:,2).*(2*lambda(:,2)-1);
%     phi(:,3) = lambda(:,3).*(2*lambda(:,3)-1);
%     phi(:,4) = 4*lambda(:,2).*lambda(:,3);
%     phi(:,5) = 4*lambda(:,3).*lambda(:,1);
%     phi(:,6) = 4*lambda(:,1).*lambda(:,2);
    bt = zeros(NT,10);
    for p = 1:nQuad
        % quadrature points in the x-y coordinate
        pxy = lambda(p,1)*node(elem(:,1),:) ...
            + lambda(p,2)*node(elem(:,2),:) ...
            + lambda(p,3)*node(elem(:,3),:);
        if isfield(pde,'f') && isnumeric(pde.f)
            fp = pde.f;        % piecewise constant
        else
            fp = pde.f(pxy);   % function handle
        end
        for j = 1:10
            bt(:,j) = bt(:,j) + w(p)*phi(p,j)*fp;
        end
    end
    bt = bt.*repmat(area,1,10);
     b = accumarray(elem2dof(:),bt(:),[Ndof 1]);
end

[b1,u] = getbdP3(b);
% b1 is the first part of the right sides, b is the second part of right sides

 function [b1,u] = getbdP3(b)
    %% Boundary conditions for Poisson equation: P3 quadratic FEM.
    %
    % The set up of boundary condition consists of two parts:
    %

    %
    %  Modify the right hand side b. The Neumann boundary integral is added
    %  to b. For Dirichlet boundary ndoes, b(fixedDof) is the evaluation of
    % pde.g_D.
    %
    % Special attentation should be given for the pure Neumann boundary
    % condition. To enforce the compatible condition, the vector b should have
    % mean value zero. To avoid a singular matrix, the 1st node is chosen as
    % fixedDof.
    %
    % The order of assigning Neumann and Dirichlet boundary condition is
    % important to get the right setting at the intersection nodes of Dirichlet
    % and Neumann boundary edges.

    u = zeros(Ndof,1);

    %% Part 1: Find boundary edges and modify the load b
    % Neumann boundary condition
        if ~isfield(option,'gNquadorder')
            option.gNquadorder = 6;
        end
        b1 = zeros(Ndof,1);
        [lambdagN,weightgN] = quadpts1(option.gNquadorder);
                idxN = (bdFlag(:) == 1);      % all Neumann edges in bdFlag
        Neumannidx = elem2edge(idxN ); % index of Neumann and Robin edges
        % since boundary integral is also needed for Robin edges
        Neumann   = edge(Neumannidx,:);

        nQuadgN = size(lambdagN,1);
        % quadratic bases (1---3---4--2)
        bdphi = zeros(nQuadgN,4);
        bdphi(:,1) = 0.5*(3*lambdagN(:,1)-1).*(3*lambdagN(:,1)-2).*lambdagN(:,1);
        bdphi(:,2) = 0.5*(3*lambdagN(:,2)-1).*(3*lambdagN(:,2)-2).*lambdagN(:,2);
        bdphi(:,3) = 9/2*lambdagN(:,1).*lambdagN(:,2).*(3*lambdagN(:,1)-1);
        bdphi(:,4) = 9/2*lambdagN(:,1).*lambdagN(:,2).*(3*lambdagN(:,2)-1);
        % length of edge

        el = sqrt(sum((node(Neumann(:,1),:) - node(Neumann(:,2),:)).^2,2));
        ge = zeros(size(Neumann,1),4);
        for pp = 1:nQuadgN
            ppxy = lambdagN(pp,1)*node(Neumann(:,1),:) ...
                 + lambdagN(pp,2)*node(Neumann(:,2),:);
            gNu = pde.g_N(ppxy);
            ge(:,1) = ge(:,1) + weightgN(pp)*gNu*bdphi(pp,1);
            ge(:,2) = ge(:,2) + weightgN(pp)*gNu*bdphi(pp,2);
            ge(:,3) = ge(:,3) + weightgN(pp)*gNu*bdphi(pp,3);
            ge(:,4) = ge(:,4) + weightgN(pp)*gNu*bdphi(pp,4);
        end
        ge = ge.*repmat(el,1,4);
        b1(1:N) =  accumarray(Neumann(:), [ge(:,1); ge(:,2)],[N,1]);
        b1(N+2*Neumannidx-1) = b1(N+2*Neumannidx-1) + ge(:,3);
        b1(N+2*Neumannidx)   = b1(N+2*Neumannidx) + ge(:,4);

    %% Part 2: Find Dirichlet boundary edges and compute the boundary value
    % Dirichlet boundary conditions

        isDirichlet(elem2edge(bdFlag(:)==1)) = true;
        % interpolation
        idx = (fixedDof > N);         % index of edge nodes
        u(fixedDof(~idx)) = pde.g_D(node(fixedDof(~idx),:)); % bd value at vertex dofs
        % for P3,  we should divide the points of edge into two parts.
        bdEdgeIdx = fixedDof(idx) - N;
        %  First parts, the points  * is in  1---*------2
        bdEdgeMid = node(edge(isDirichlet,1),:)+(node(edge(isDirichlet,2),:) ...
                  - node(edge(isDirichlet,1),:))/3;
        u(N + bdEdgeIdx(1:2:end)) = pde.g_D(bdEdgeMid);
      %  Second parts, the points * is in  1------*---2
        bdEdgeMid = node(edge(isDirichlet,1),:)+2*(node(edge(isDirichlet,2),:)...
                  - node(edge(isDirichlet,1),:))/3;
        u(N + bdEdgeIdx(2:2:end)) = pde.g_D(bdEdgeMid);
        % modify the right hand side
         b1 = b1 - B*u;
         %b1(fixedDof) = u(fixedDof);

    end % end of getbdP3

B(:,fixedDof)=[];
Nu = Ndof-size(fixedDof,1);

b(fixedDof) = [];
     bigA = [A, B; ...
             B', sparse(Nu,Nu)];
     bigF = [b1; -b];
% Solver
%bigU=PFGMRES(bigA, bigF,sparse(Ndof+Nu,1), Ndof, 10, 1e-6, [],[]);
tic;
bigU=bigA\bigF;
toc;
%   norm(bigA*bigU-bigF)
%   norm(bigA*[p;u(freeDof)]-bigF)
   w=bigU(1:Ndof);
   u(freeDof)=bigU(Ndof+1:Ndof+Nu);

end                 % end of function PoissonP2

repo: Aero-matlab


text:
function [debug] = convertERTSimLogs(data_file, parser_config, save_file, slow_on)

% Author: Alexander Wischnewski     Date: 21-12-2018
%
% Description:
%   converts the logs of a GRT target based debug file (e.g. compiled
%   binary simulation on windows or unix) into the common log format.
%
% Input:
%   data_file:              Log data file
%   parser_config:          Parser config
%   save_file:              Output data file
%   slow_on:                true/false value for parsing extended slow
%   dataset
% Output:
%   debug:                  Struct with all timeseries in there

    % check if file ending is given and modify load files accordingly
    if(endsWith(data_file, '.txt'))
        file_standard = data_file;
        file_slow = [data_file(1:end-4) '_slow.txt'];
    else
        file_standard = [data_file '.txt'];
        file_slow = [data_file '_slow.txt'];
    end
    % check if parser config ending is given and modify load files accordingly
    if(endsWith(parser_config, '.txt'))
        config_standard = parser_config;
        config_slow = [parser_config(1:end-4) '_slow.txt'];
    else
        config_standard = [parser_config '.txt'];
        config_slow = [parser_config '_slow.txt'];
    end
    % load and parse data from standard logging file
    try
        data = dlmread(file_standard, ';');
        % skip last data point as it might not be fully written!
        debug = convertArrayToTs(data(1:end-1, :), config_standard);
        disp('Successfully parsed standard logging file!');
    catch
        disp('Could not parse standard logging file!');
        debug = 0;
    end

    % load and parse data from slow logging file
    if slow_on
        data_slow = dlmread(file_slow, ';');
        % skip last data point as it might not be fully written!
        debug_slow = convertArrayToTs(data_slow(1:end-1, :), config_slow);
        debug_slow = convertSeperateTsToArrayTs(debug_slow);
        disp('Successfully parsed slow logging file!');
            % write to file if its wanted by the user
        if(nargin == 4)
          save(save_file, 'debug', 'debug_slow', '-v7.3');
        end
    else
        % write to file if its wanted by the user
        if(nargin == 4)
          save(save_file, 'debug','-v7.3');
        end
    end
end

author: TUMFTM


text:
function list = shellSort(list)

    N = numel(list);
    increment = round(N/2);

    while increment > 0  %loop until increment becomes 0

        for i = (increment+1:N)
            temp = list(i);
            j = i;
            while (j >= increment+1) && (list(j-increment) > temp)
                list(j) = list(j-increment);
                j = j - increment;
            end

            list(j) = temp;

        end %for

        if increment == 2 %This case causes shell sort to become insertion sort
            increment = 1;
        else
            increment = round(increment/2.2);
        end
    end %while
end %shellSort

author: TheAlgorithms

repo: MATLAB-Octave


text:
%% Example
% real time control of the KUKA iiwa 7 R 800
% the impedence control is turned on
% Moving first joint of the robot, using a sinisoidal function

% An example script, it is used to show how to use the different
% functions of the KUKA Sunrise matlab toolbox

% First start the server on the KUKA iiwa controller
% Then run this script using Matlab

% Copyright: Mohammad SAFEEA, 8th of Nov 2017

% Important: Be careful when runnning the script, be sure that no human, nor obstacles
% are around the robot

close all,clear all;clc;

warning('off')

ip='172.31.1.147'; % The IP of the controller
% start a connection with the server
global t_Kuka;
t_Kuka=net_establishConnection( ip );

if ~exist('t_Kuka','var') || isempty(t_Kuka) || strcmp(t_Kuka.Status,'closed')
  warning('Connection could not be establised, script aborted');
  return;
else

      %% Get position roientation of end effector

%      jPos={0,pi/6,0,-pi/2,0,pi/2-pi/6,0};

        jPos={0,0,0,-pi/2,0,pi/2,0};
      %jPos={0,0,0,0,0,0,0};
      setBlueOff(t_Kuka); % turn Off blue light

      relVel=0.15;
      movePTPJointSpace( t_Kuka , jPos, relVel); % move to initial configuration
     %% Pause for 3 seocnds
     pause(3);
        %% Start direct servo in joint space
        massOfTool=0.5; % the mass of the tool attached to flange in Kg
        cOMx=0; % X coordinate of the center of mass of the tool in (mm)
        cOMy=0; % Y coordinate of the center of mass of the tool in (mm)
        cOMz=40; % Z coordinate of the center of mass of the tool in (mm)
        cStiness=900; % cartizian stifness
        rStifness=80; % rotational stifness
        nStifness=50; % null space stifness

        % Start the realtime control with impedence
        realTime_startImpedanceJoints(t_Kuka,massOfTool,cOMx,cOMy,cOMz,...
        cStiness,rStifness,nStifness);

       w=0.6; % motion constants, frequency rad/sec
       A=0.2; % motion constants, amplitude of motion

       a=datevec(now);
       t0=a(6)+a(5)*60+a(4)*60*60; % calculate initial time

       dt=0;
     precission=1000;
     precission_flag=true;
     tstart=t0;
     counter=0;
     duration=1*60; %1 minutes
     time_stamps=zeros(1,1000*duration);
       while(dt<duration)
         %% ferform trajectory calculation here
          a=datevec(now);
          time=a(6)+a(5)*60+a(4)*60*60;
          dt=time-t0;

          temp=A*(1-cos(w*dt));
            for dacount=1:7
                jPosCommand{dacount}=jPos{dacount}+temp;
            end
          counter=counter+1;
          %% Send joint positions to robot
          %sendJointsPositions( t ,jPos);
          sendJointsPositionsf( t_Kuka ,jPosCommand);
          time_stamps(counter)=dt;

       end
       tend=time;
       rate=counter/(tend-tstart);

       %% Stop the realtime control with impedence motion
       realTime_stopImpedanceJoints( t_Kuka );

       fprintf('\nThe rate of joint nagles update per second is: \n');
       disp(rate);
       fprintf('\n')
       pause(2);
       %% turn off light
       setBlueOff(t_Kuka);

      %% turn off the server
       net_turnOffServer( t_Kuka );

       fclose(t_Kuka);

       %% save time stamps
       save('timingdata.mat','time_stamps');

end
 warning('on')

repo: KST-Kuka-Sunrise-Toolbox


text:
function [Q,err] = quadgr(fun,a,b,tol,trace,varargin)
%QUADGR Gauss-Legendre quadrature with Richardson extrapolation.
%   [Q,ERR] = QUADGR(FUN,A,B,TOL) approximates the integral of a function
%   FUN from A to B with an absolute error tolerance TOL. FUN is a function
%   handle and must accept vector arguments. TOL is 1e-6 by default. Q is
%   the integral approximation and ERR is an estimate of the absolute error.
%
%   QUADGR uses a 12-point Gauss-Legendre quadrature. The error estimate is
%   based on successive interval bisection. Richardson extrapolation
%   accelerates the convergence for some integrals, especially integrals
%   with endpoint singularities.
%
%   QUADGR(FUN,A,B,TOL,TRACE) with non-zero TRACE displays the
%   extrapolation table.
%
%   QUADGR can also be used as the quadrature in DBLQUAD and TRIPLEQUAD.
%
%   Examples:
%     Q = quadgr(@(x) log(x),0,1)
%     [Q,err] = quadgr(@(x) exp(x),0,9999i*pi)
%     [Q,err] = quadgr(@(x) sqrt(4-x.^2),0,2,1e-12)
%     [Q,err] = quadgr(@(x) x.^-0.75,0,1)
%     [Q,err] = quadgr(@(x) 1./sqrt(1-x.^2),-1,1)
%     [Q,err] = quadgr(@(x) exp(-x.^2),-inf,inf,1e-9) % sqrt(pi)
%     [Q,err] = quadgr(@(x) cos(x).*exp(-x),0,inf,1e-9)
%
%   See also QUAD, QUADGK, DBLQUAD, TRIPLEQUAD

%   Author: Jonas Lundgren <splinefit@gmail.com> 2009

%   2009-03-17  First published
%   2010-04-14  Adapted to DBLQUAD and TRIPLEQUAD (varargin added)

if nargin < 3, help quadgr, return, end
if nargin < 4 || isempty(tol), tol = 1.e-6; end;
if nargin < 5 || isempty(trace), trace = 0; end;

% Order limits (required if infinite limits)
if a == b
    Q = b - a;
    err = b - a;
    return
elseif a > b
    reverse = true;
    atmp = a;
    a = b;
    b = atmp;
else
    reverse = false;
end

% Infinite limits
if isinf(a) || isinf(b)
    % Check real limits
    if ~isreal(a) || ~isreal(b) || isnan(a) || isnan(b)
        error('quadgr:inflim','Infinite intervals must be real.')
    end
    % Change of variable
    if isfinite(a) && isinf(b)
        % a to inf
        fun1 = @(t,varargin) fun(a + t./(1-t), varargin{:})./(1-t).^2;
        [Q,err] = quadgr(fun1,0,1,tol,trace,varargin{:});
    elseif isinf(a) && isfinite(b)
        % -inf to b
        fun2 = @(t,varargin) fun(b + t./(1+t), varargin{:})./(1+t).^2;
        [Q,err] = quadgr(fun2,-1,0,tol,trace,varargin{:});
    else % -inf to inf
        fun1 = @(t,varargin) fun(t./(1-t), varargin{:})./(1-t).^2;
        fun2 = @(t,varargin) fun(t./(1+t), varargin{:})./(1+t).^2;
        [Q1,err1] = quadgr(fun1,0,1,tol/2,trace,varargin{:});
        [Q2,err2] = quadgr(fun2,-1,0,tol/2,trace,varargin{:});
        Q = Q1 + Q2;
        err = err1 + err2;
    end
    % Reverse direction
    if reverse
        Q = -Q;
    end
    return
end

% Gauss-Legendre quadrature (12-point)
xq = [0.12523340851146894; 0.36783149899818018; 0.58731795428661748; ...
      0.76990267419430469; 0.9041172563704748; 0.98156063424671924];
wq = [0.24914704581340288, 0.23349253653835478, 0.20316742672306584, ...
      0.16007832854334636, 0.10693932599531818, 0.047175336386511842];
xq = [xq; -xq];
wq = [wq, wq];
nq = length(xq);

% Initiate vectors
maxit = 17;                 % Max number of iterations
Q0 = zeros(maxit,1);       	% Quadrature
Q1 = zeros(maxit,1);       	% First Richardson extrapolation
Q2 = zeros(maxit,1);       	% Second Richardson extrapolation

% One interval
hh = (b - a)/2;             % Half interval length
x = (a + b)/2 + hh*xq;      % Nodes
% Quadrature
Q0(1) = hh*wq*fun(x,varargin{:});

% Successive bisection of intervals
for k = 2:maxit

    % Interval bisection
    hh = hh/2;
    x = [x + a; x + b]/2;
    % Quadrature
    Q0(k) = hh*wq*sum(reshape(fun(x,varargin{:}),nq,[]),2);

    % Richardson extrapolation
    if k >= 5
        Q1(k) = richardson(Q0,k);
        Q2(k) = richardson(Q1,k);
    elseif k >= 3
        Q1(k) = richardson(Q0,k);
    end

    % Estimate absolute error
    if k >= 6
        Qv = [Q0(k), Q1(k), Q2(k)];
        Qw = [Q0(k-1), Q1(k-1), Q2(k-1)];
    elseif k >= 4
        Qv = [Q0(k), Q1(k)];
        Qw = [Q0(k-1), Q1(k-1)];
    else
        Qv = Q0(k);
        Qw = Q0(k-1);
    end
    [err,j] = min(abs(Qv - Qw));
    Q = Qv(j);

    % Convergence
    if err < tol || ~isfinite(Q)
        break;
    end

end

% Convergence check
if ~isfinite(Q)
    warning('quadgr:infnan','Integral approximation is Infinite or NaN.')
elseif err > tol
    warning('quadgr:maxiter','Max number of iterations reached without convergence.')
end

% The error estimate should not be zero
err = err + 2*eps(Q);
% Reverse direction
if reverse
	Q = -Q;
end

% Display convergence
if trace
    disp(' ')
    disp([Q0(1:k) Q1(1:k) Q2(1:k)])
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function R = richardson(Q,k)
% Richardson extrapolation with parameter estimation
if Q(k) ~= Q(k-1)
    c = real((Q(k-1)-Q(k-2))/(Q(k)-Q(k-1))) - 1;
else
    c = 1;
end
% The lower bound 0.07 admits the singularity x.^-0.9
c = max(c,0.07);
R = Q(k) + (Q(k) - Q(k-1))/c;


text:
function x = ProxLP(y,p,tau)

if p==Inf
    p = 1e3;
end

if length(y(:))>1
    for i=1:length(y(:))
        x(i)=ProxLP(y(i),p,tau);
    end
    x = reshape(x,size(y));
    return
end

x = linspace(0,abs(y)*10,1000);
[~,k] = min( 1/2*(x-abs(y)).^2 + tau*abs(x).^p );
x = x(k)*sign(y);

end

author: optimaltransport

path: github-repos/MATLAB/optimaltransport-optimaltransport.github.io/optimaltransport.github.io-2fa6db6e6a48ab9bd6676088db00bd7c5c8b4203/code/gradflow-metric/ProxLP (Copie en conflit de Gabriel Peyré 2017-11-26).m