text:
function[Translations]=decode_greedy(parameter,batch,filename)
    T=size(batch.source_sen_matrix,2);
    % total number of time-steps
    N=size(batch.source_sen_matrix,1);
    % number of sentences
    L=length(batch.source_smallBatch);
    Max_Length=1.5*max(batch.num_word_in_doc);

    result=Forward(batch,parameter,0);
    % obtain embeddings for inputs
    first_words=BeamStep(parameter,result.source_sen{parameter.layer_num,1},1);
    % decode first words
    beamHistory=oneMatrix([N,Max_Length]);
    % store in history
    beamHistory(:,1)=first_words(:);

    for ll=1:parameter.layer_num
        beamSenStates_h{ll}=result.source_sen{ll,1};
        beamWordStates_h{ll}=result.source_sen{ll,1};
        beamSenStates_c{ll}=result.c_t_source_sen{ll,T};
        beamWordStates_c{ll}=result.c_t_source_sen{ll,T};
    end
    % update current lstm states
    for position=1:Max_Length
        words=beamHistory(:,position);
        % get earliest predicted words
        for ll=1:parameter.layer_num
            if ll==1
                x_t=parameter.vect(:,words);
            else
                x_t=beamWordStates_h{ll-1};
            end
            h_t_1=beamWordStates_h{ll};
            c_t_1=beamWordStates_c{ll};
            [A,beamWordStates_h{ll},beamWordStates_c{ll}]=lstmUnit(parameter.Word_T{ll},parameter,x_t,h_t_1,c_t_1,ll,-1,0);
            % lstm unit calculation, get lsmt status for current time step at word level
        end
        end_sen_index=find(words==parameter.sen_stop);
        % if predicted words incode sentence-ending indicators
        if length(end_sen_index)~=0
            for ll=1:parameter.layer_num
                if ll==1
                    x_t=beamWordStates_h{ll}(:,end_sen_index);
                else
                    x_t=beamSenStates_h{ll-1}(:,end_sen_index);
                end
                h_t_1=beamSenStates_h{ll}(:,end_sen_index);
                c_t_1=beamSenStates_c{ll}(:,end_sen_index);
                [A,beamSenStates_h{ll}(:,end_sen_index),beamSenStates_c{ll}(:,end_sen_index)]=lstmUnit(parameter.Sen_T{ll},parameter,x_t,h_t_1,c_t_1,ll,-1,0);
                % update lsmt status for current time step at sentence level
                beamWordStates_h{ll}(:,end_sen_index)=beamSenStates_h{ll}(:,end_sen_index);
                beamWordStates_c{ll}(:,end_sen_index)=beamSenStates_c{ll}(:,end_sen_index);
            end
        end
        all_next_words=BeamStep(parameter,beamWordStates_h{parameter.layer_num},0);
        % predict next word
        beamHistory(:,position+1)=all_next_words;
    end
    for senId=1:N
        vector=beamHistory(senId,:);
        vector=vector(1:floor(1.5*batch.num_word_in_doc(senId)));
        % maximium length of 1.5 * length of orginal inputs
        stop_sign=find(vector==parameter.doc_stop);
        if length(stop_sign)==0
            vector=vector;
        else
            vector=vector(1:stop_sign-1);
        end
        M=find(vector==parameter.sen_stop);
        vector(M)=[];
        % remove sentence-ending token
        dlmwrite(filename,vector,'delimiter',' ','-append');
    end
end

function[select_words]=BeamStep(parameter,h_t,isFirst)
    if isFirst==1 scores=parameter.soft_W(1:end-2,:)*h_t;
    else scores=parameter.soft_W*h_t;
    end
    % could not predict sentence-ending token or document-ending token for first word
    mx = max(scores);
    scores = bsxfun(@minus, scores, mx);
    logP=bsxfun(@minus, scores, log(sum(exp(scores))));
    [sortedLogP, sortedWords]=sort(logP, 'descend');
    select_words=sortedWords(1, :);

    probs_=exp(scores);
    norms = sum(probs_, 1);
    probs=bsxfun(@rdivide, probs_, norms);
end

author: jiweil

repo: Hierarchical-Neural-Autoencoder


text:
%  Figure 5.45 Feedback Control of Dynamic Systems, 5e
%                   Franklin, Powell, Emami
%
clf
numG=[1 0];
denG=[1 1 4];
K1=0:.1:8;   % K = 2 at breakin point
K2=100;
Kt=[K1 K2];
axis('square')
r=rlocus(numG,denG,Kt);    % roots vs. Kt
plot(r,'--'),grid on;
axis([-6 2 -4 4])
hold on
p=roots(denG);             % find and plot OL poles
z=roots(numG);             % find and plot OL zeros
plot(z(1),0,'o')
plot(real(p(1)),imag(p(1)),'x',real(p(2)),imag(p(2)),'x')
title('Fig. 5.45  Root Locus vs. K with K_{T} = 1')
xlabel('Re(s)')
ylabel('Im(s)')

% now do RL vs. K

Kt1 = 1;
r1=rlocus(numG,denG,Kt1);    % root at Kt = 1
plot(r1,'*')
num = 1;
den = poly(r1);
K = [0 30]
r2=rlocus(num,den,K);    % root locus vs. K
plot(r2,'-')

hold off
axis('normal')

author: Sable

repo: mcbench-benchmarks


function price = Price_Spread_Option_BjerksundStensland_2D(K, S0_1, S0_2, T, r, rho, sigma_1, sigma_2, q_1, q_2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% About: calcuates Price of 2D Spread option, (S_1(T) - S_2(T) - K)+, using Bjerksund & Stensland approximation
%        Note: when K = 0, this agrees with Magrabes exact formula
% Author: Justin Lars Kirkby
%
% Reference: Bjerksund, P. and Stensland, G. (2006): "Closed form spread option valuation"
%
% -----------------
% Params
% -----------------
% S0_1    = initial asset price of first asset, e.g. S0_1 = 100
% S0_2    = initial asset price of second asset, e.g. S0_2 = 100
% T       = time to maturity in years (e.g. T=1)
% r       = interest rate
% rho     = instantaneous correlation between S_1 and S_2
% sigma_1 = volatility of first asset (annualized), e.g. sigma = 0.2
% sigma_2 = volatility of second asset (annualized), e.g. sigma = 0.2
% q_1     = dividend yield of first asset, e.g. q_1 = 0.02
% q_2     = dividend yield of second asset, e.g. q_2 = 0.02
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

F1 = S0_1*exp((r-q_1)*T);
F2 = S0_2*exp((r-q_2)*T);

a = F2 + K;
b = F2 / a;
rhosigs = rho*sigma_1*sigma_2;
sig = sqrt(sigma_1^2 - 2*rhosigs*b + b^2*sigma_2^2);
sigst = sig*sqrt(T);

d1 = (log(F1/a) + (0.5*sigma_1^2 - b*rhosigs + 0.5*b^2*sigma_2^2)*T) / sigst;
d2 =  (log(F1/a) + (-0.5*sigma_1^2 + rhosigs + (0.5*b^2 - b)*sigma_2^2)*T) / sigst;
d3 = (log(F1/a) + (-0.5*sigma_1^2 + 0.5*b^2*sigma_2^2)*T) / sigst;

price = exp(-r*T)*(F1*normcdf(d1) - F2*normcdf(d2) - K*normcdf(d3));

end


text:
function [Hes,R] = VAorthog(Z,n,varargin)  % Vand.+Arnoldi orthogonalization
%VAORTHOG  Vandermonde with Arnoldi orthogonalization.
%
% This code comes from Brubeck and Trefethen, "Lightning Stokes solver", arXiv 2020.
% For the mathematics, see Brubeck, Nakatsukasa, and Trefethen, "Vandermonde with
% Arnoldi", SIAM Review, 2021.
%
%  Input:   Z = column vector of sample points
%           n = degree of polynomial (>= 0)
%         Pol = cell array of vectors of poles (optional)
% Output: Hes = cell array of Hessenberg matrices (length 1+length(Pol))
%           R = matrix of basis vectors
M = length(Z); Pol = []; if nargin == 3, Pol = varargin{1}; end
% First orthogonalize the polynomial part
Q = ones(M,1); H = zeros(n+1,n);
for k = 1:n
   q = Z.*Q(:,k);
   for j = 1:k, H(j,k) = Q(:,j)'*q/M; q = q - H(j,k)*Q(:,j); end
   H(k+1,k) = norm(q)/sqrt(M); Q(:,k+1) = q/H(k+1,k);
end
Hes{1} = H; R = Q;
% Next orthogonalize the pole parts, if any
while ~isempty(Pol)
   pol = Pol{1}; Pol(1) = [];
   np = length(pol); H = zeros(np,np-1); Q = ones(M,1);
   for k = 1:np
      q = Q(:,k)./(Z-pol(k));
      for j = 1:k, H(j,k) = Q(:,j)'*q/M; q = q - H(j,k)*Q(:,j); end
      H(k+1,k) = norm(q)/sqrt(M); Q(:,k+1) = q/H(k+1,k);
   end
   Hes{length(Hes)+1} = H; R = [R Q(:,2:end)];
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function  F = FGG_3d(f,knots,N,accuracy,GridListx, GridListy, GridListz)
%Description:
%This code implements the 3D version of the "accelerated"
%Gaussian-gridding-based NUFFT described in Greengard and Lee [1]. The
%gridding approach is very similar to previous work by Nguyen and Liu [2];
%the only difference is the use of a different convolution kernel ([1]
%claims that the Gaussian kernel has computational advantages). Both
%algorithms allow the user to specify the numerical precision of the
%routine, but [1] provides a nice summary that would allow one to tabulate
%the appropriate variable values for each desired numerical precision. Code
%for [2] is also available upon request.
%
%This code performs NUFFTs that form rectangular cubes of size N=[Nx,Ny,Nz]
%(not counting frequency padding for image interpolation). The approximate
%DFT attains errors on the order of 1e-6 (for more or less accuracy, vary
%the optional "accuracy" parameter).
%
%Inputs:
%       f: frequency-domain data (a complex Mx1 vector) unwrapped from a
%           matrix knots: k-space locations at which the data were measured
%           (an Mx3 vector). This data must be in double format.
%       knots: the frequency locations of the data points. These locations
%           should be normalized to correspond to the grid boundaries
%           [-N/2, N/2 -1/N], N even. If the knots are not scaled
%           properly, they will be shifted and scaled into this normalized
%           form.
%       N = [Nx,Ny,Nz]: the 1x3 vector denoting the size of the spatial
%           grid in the image domain. This parameter will determine the
%           spatial extent of the image.
%       accuracy: (optional input parameter) a positive integer indicating
%           the desired number of digits of accuracy
%Optional Inputs (highly recommended):
%       GridListx:(column vector) The x-locations of the
%           frequency grid (not yet scaled by c or 2*pi) onto which the
%           data should be interpolated
%       GridListy:(column vector) The y-locations of the
%           frequency grid (not yet scaled by c or 2*pi) onto which the
%           data should be interpolated
%       GridListz:(column vector) The z-locations of the
%           frequency grid (not yet scaled by c or 2*pi) onto which the
%           data should be interpolated
%Outputs:
%       F: the 3D NUFFT (approximate DFT) of f, with dimension [Nx,Ny,Nz].
%
%
%Usage Notes:
%In order for this function to work, the C
%file "FGG_Convolution3D.c" must be compiled into a Matlab executable
%(cmex) with the following command in the command prompt:
%
%mex FGG_Convolution3D.c
%
%A note on the effect of M_sp on the algorithm's accuracy:
%[R,M_sp]=[2,3]  ==> 1e-3 accuracy
%[R,M_sp]=[2,6]  ==> single precision
%[R,M_sp]=[2,9]  ==> 1e-9 accuracy
%[R,M_sp]=[2,12] ==> double precision
%
%References:
%[1] L. Greengard and J. Lee, "Accelerating the Nonuniform Fast Fourier
%Transform," SIAM Review, Vol. 46, No. 3, pp. 443-454.
%[2] N. Nguyen and Q. H. Liu, "Nonuniform fast fourier transforms," SIAM J.
%Sci. Comput., 1999.
%
%Please note:
%This code is free to use, but we ask that you please reference the source,
%as this will encourage future funding for more free AFRL products. This
%code was developed through the AFOSR Lab Task "Moving-Target Radar Feature
%Extraction."
%Project Manager: Arje Nachman
%Principal Investigator: Matthew Ferrara
%Date: November 2008
%
%Code by (send correspondence to):
%Matthew Ferrara, Research Mathematician
%AFRL Sensors Directorate Innovative Algorithms Branch (AFRL/RYAT)
%Matthew.Ferrara@wpafb.af.mil

%Explanation of variables used
%bw             The bandwidth of the input data (fmax-fmin)
%D              A [Nx,Ny,Nz] deconvolution matrix formed from the 3D
%               kronecker product of E4_x, E4_y, and E4_z
%dgx,dgy,dgz    The separation (in the frequency domain) of the
%               user-defined k-space grid onto which the data should be
%               interpolated
%E1, E2, E3     factors of the Gaussian filter in the frequency domain
%               (the Gaussian is factored to eliminate redundant
%               exponential calculations)
%E4             The Gaussian deconvolution filter
%f              The Mx1 vector of nonuniformly-spaced frequency data given
%               as input to the NUFFT routine
%f_tau          The [R*Nx,R*Ny] matrix of uniformly-spaced fequency-domain
%               data values after "gridding"
%F_tau          The FFT of f_tau
%F              The approximate DFT of f (the deconvolved FFT of f_tau)
%fmean          The average knot values of the input data (1x3 vector)
%j              Index variable that indexes the convolution loop through
%               the data points (1 <= j <= M)
%kmin           The minimum knot values of the input data (1x3 vector)
%kmax           The maximum knot values of the input data (1x3 vector)
%knots          The Mx3 vector of frequency locations given as input to
%               the NUFFT.
%M              The number of (k-space) data points (the length of the
%               input data vector f)
%M_sp           Width of the frequency-domain box used in the approximate
%               interpolation of each data point onto the frequency grid
%               (M_sp=6 for single precision and M_sp=12 for double
%               precision)
%N              The image size of the NUFFT output (N=[Nx, Ny, Nz])
%Nx             The length of the image in the x dimension
%Ny             The length of the image in the y dimension
%Nz             The length of the image in the z dimension
%R              Oversampling ratio for gridding in the frequency (data)
%               domain
%S              The [Nz,Nx*Ny]-sized kronecker product of E4z, E4_x, and
%               E4_y
%scale          1x3 vector used to scale the input data locations into the
%               normalized form
%shift          1x3 vector used to shift the input data locations into the
%               normalized form
%tau            The 3x1 Gaussian kernel spreading factor
%End Explanation of variables

%% Step 1: Initialize constant variables:
M=length(f);%number of frequency-domain data points
if nargin<4, accuracy=6; end
if nargin<3, N=M; accuracy=6; end
if length(N)<2, N=N(1)*[1,1,1];accuracy=6; end
%The size parameters [Nx,Ny] determine the spatial extent of the image
Nx=N(1); Ny=N(2); Nz=N(3);
%R is the oversampling ratio(>1) for gridding in the frequency (data)
%domain. There are diminishing returns in accuracy after R=2 (M_sp has a
%more direct effect on accuracy).
R=2;
%M_sp is the length of the convolution kernel
M_sp=accuracy;%This gives roughly 6 digits of accuracy
%The variance, tau, of the Gaussian filter may be different in each
%dimension
%tau = M_sp./(N.^2);%I was initially using this value
tau = (pi*M_sp./(N.*N*R*(R-.5)));%Suggested value of tau by Greengard [1]
%The length of the oversampled grid
M_r = R*N;
%Scale knots (data locations) to [-N/2,N/2-1] (I don't initially use
%Greengard's [0,2*pi] convention, but instead assume users will most likely
%be thinking in terms of fs<=f<=fe)
kmin=min(knots);
kmax=max(knots);

if nargin <=4,%Simply choose the most convenient frequency-domain grid
    bw=kmax-kmin;
    scale=(N-1)./bw;
    shift=-N/2-kmin.*scale;
    knots=repmat(scale,[M,1]).*knots + repmat(shift,[M,1]);
else%we specify knot locations in terms of the user-defined grid (this
    %consequently specifies the image pixel locations)
    fmean=(kmin+kmax)/2;
    dgx=GridListx(2)-GridListx(1);
    dgy=GridListy(2)-GridListy(1);
    dgz=GridListz(2)-GridListz(1);
    %We choose to ensure that the data are perfectly centered on the grid:
    kminx=fmean(1)-(Nx/2)*dgx;
    kmaxx=kminx+(Nx-1)*dgx;
    kminy=fmean(2)-(Ny/2)*dgy;
    kmaxy=kminy+(Ny-1)*dgy;
    kminz=fmean(3)-(Nz/2)*dgz;
    kmaxz=kminz+(Nz-1)*dgz;
    kmin=[kminx,  kminy, kminz];
    kmax=[kmaxx,  kmaxy, kmaxz];
    %The BW that covers the whole k-space region when confined to the
    %specified grid spacing
    bw=[kmaxx-kminx,  kmaxy-kminy,  kmaxz-kminz];
    scale=(N-1)./bw;
    shift=-.5*[Nx,Ny,Nz]-[kminx,kminy,kminz].*scale;
    knots=repmat(scale,[M,1]).*knots + repmat(shift,[M,1]);
end
%Switch knot locations to [0,2*pi] convention (used by Greengard):
knots=mod(2*pi*knots./repmat(N,[M,1]),2*pi);%Makes NUFFT implementation
%more straightforward when notation is the same!

%Precompute E_3, the constant component of the (truncated) Gaussian:
E_3x(1,1:M_sp) = exp(-((pi*(1:M_sp)/M_r(1)).^2)/tau(1));
%don't waste (slow) exponential calculations
E_3x=[fliplr(E_3x(1:(M_sp-1))),1,E_3x];
E_3y(1,1:M_sp) = exp(-((pi*(1:M_sp)/M_r(2)).^2)/tau(2));
%don't waste (slow) exponential calculations
E_3y=[fliplr(E_3y(1:(M_sp-1))),1,E_3y];
E_3z(1,1:M_sp) = exp(-((pi*(1:M_sp)/M_r(3)).^2)/tau(3));
%don't waste (slow) exponential calculations
E_3z=[fliplr(E_3z(1:(M_sp-1))),1,E_3z];

%Precompute E_4 (for devonvolution after the FFT)
kx_vec = (-Nx/2):(Nx/2-1);
%The Hadamard Inverse of the Fourier Transform of the truncated Gaussian
E_4x(1:Nx,1)=sqrt(pi/tau(1))*exp(tau(1)*(kx_vec.^2));
ky_vec = (-Ny/2):(Ny/2-1);
%The Hadamard Inverse of the Fourier Transform of the truncated Gaussian
E_4y(1,1:Ny)=sqrt(pi/tau(2))*exp(tau(2)*(ky_vec.^2));%
kz_vec = (-Nz/2):(Nz/2-1);
%The Hadamard Inverse of the Fourier Transform of the truncated Gaussian
E_4z(1:Nz,1)=sqrt(pi/tau(3))*exp(tau(3)*(kz_vec.^2));%
%End initialization of constant variables

%% Step 2: Approximate convolution for each datum location, (x_j,y_j). This
% step is implemented in C and compiled into a Matlab-executable (cmex)
% file.
%Initialize convolved data matrix
f_tau=zeros(M_r(1),M_r(2),M_r(3));
f_taui=zeros(M_r(1)*M_r(2)*M_r(3),1);%Imaginary components of f_tau
f_taur=zeros(M_r(1)*M_r(2)*M_r(3),1);%Real components of f_tau
%Perform convolution onto finely-spaced grid

[f_taur,f_taui]=...
    FGG_Convolution3D(double(real(f(:))),double(imag(f(:))),...
    double(knots(:)),double(E_3x), double(E_3y), double(E_3z), ...
    [double(M_sp), double(tau(1)), double(tau(2)), double(tau(3)),...
 double(M_r(1)), double(M_r(2)), double(M_r(3))]);

f_tau = reshape(f_taur+sqrt(-1)*f_taui,[M_r(1), M_r(2), M_r(3)]);

%% Step 3: Perform FFT and deconvolve the result
%Perform FFT
F_tau=fftshift(fftn(ifftshift(f_tau)));
%Chop off excess pixels(the fine spacing in the frequency domain expanded
%the image in the transform domain)
F_tau(1:round(.5*(R-1)*Nx),:,:)=[];
F_tau(Nx+1:end,:,:)=[];
F_tau(:,1:round(.5*(R-1)*Ny),:)=[];
F_tau(:,Ny+1:end,:)=[];
F_tau(:,:,1:round(.5*(R-1)*Nz))=[];
F_tau(:,:,Nz+1:end)=[];

%Deconvolve the FFT by Hadamard multiplication with the Gaussian
%To form the deconvolution matrix S, we could do this with a loop:
%for i=1:Nx
%    for j=1:Ny
%        for k=1:Nz
%        S(i,j,k)=E_4x(i)*E_4y(j)*E_4z(k);
%        end
%    end
%end
%HOWEVER, Matlab is extremely slow with loops so we instead use the kron()
%function and reshape the result (alternatively, we could have written
%another mex function):
S=kron(E_4z,E_4x*E_4y);%Scalars for 3D kronecker product (matrix must be
                       %reshaped)
D=permute(reshape(S,[Nx,Nz,Ny]),[1,3,2]);
F=F_tau.*D/(M*R*R*R); %Note that the scaling factor is 1/M, which is
%apparently not 1/M_r as shown in eqn (9) in [1]. This makes sense because
%the scaling factor for the DFT we are attempting to approximate (in eqn
%(1) of [1]) is 1/M.

%% If desired, compare to DFT
% F_true=zeros(Nx,Ny,Nz);
% for m=(-Nx/2):(Nx/2 -1)
%     for n=(-Ny/2):(Ny/2 -1)
%           for o=(-Nz/2):(Nz/2 -1)
%         for k=1:M
%             F_true(m+Nx/2+1,n+Ny/2+1,o+Nz/2+1)=...
%               F_true(m+Nx/2+1,n+Ny/2+1,o+Nz/2+1)+...
%               f(k)*exp(sqrt(-1)*(m*knots(k,1) + n*knots(k,2) + ...
%               o*knots(k,3)));
%         end
%     end
% end
% F_true=F_true/M;
%
% figure
% isosurface(abs(F))
% title('F via NUFFT')
%
% figure
% isosurface(abs(F_true))
% title('F via DFT')
%
%error('done!')


text:
function  a = polygon_area(x,y)
% AREA  Area of a planar polygon.
%	AREA(X,Y) Calculates the area of a 2-dimensional
%	polygon formed by vertices with coordinate vectors
%	X and Y. The result is direction-sensitive: the
%	area is positive if the bounding contour is counter-
%	clockwise and negative if it is clockwise.
%
%	See also TRAPZ.

%  Copyright (c) 1995 by Kirill K. Pankratov,
%	kirill@plume.mit.edu.
%	04/20/94, 05/20/95

 % Make polygon closed .............
x = [x(:); x(1)];
y = [y(:); y(1)];

 % Calculate contour integral Int -y*dx  (same as Int x*dy).
lx = length(x);
a = -(x(2:lx)-x(1:lx-1))'*(y(1:lx-1)+y(2:lx))/2;
a = abs(a);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%CYLGRID Generate 3d cylindrical hexahedral grid.
%
%   [ GRID ] = CYLGRID( NS, NR, NZ, R, LZ, XP, AX ) Generates a hexahedral
%   grid for a cylindrical domain with NS+NR cells in the radial direction.
%   NS specifies the cell resolution of the inner square (default 4), and NR
%   ns number of cells in the radial direction of the outer layer (default 3).
%   NZ specity the numver of cells in the lengthwise direction and LZ the
%   corresponding length. The optional arguments R, XP, and AX specify the
%   radius, center coordinates of the cylinder (default R = 1 and XP = [0;0;0]),
%   and the axis of alignment (default 1 equals to the x-axis).
%
%   Examples:
%
%      1) A cylindrical grid with radius 1 and length 1:
%
%      grid = cylgrid();
%
%      2) Cylinder with radius 0.5 length 2 centered and extending
%         in the negative y-direction from the point [1 1 2].
%
%      grid = cylgrid( 3, 4, 10, 0.5, 2, [1;1;2], -2 );
%
%   See also BLOCKGRID, CIRCGRID, HOLEGRID, LINEGRID, RECTGRID, RINGGRID, SPHEREGRID

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function [] = baselineSIRFS(class,jobID)
%GETDEPTHIMAGE Summary of this function goes here
%   Detailed explanation goes here

globals
statesDir = jobDirs(class,jobID,'state');
sirfsstatedir = fullfile(cachedir,class,sprintf('baselineSIRFS%s',jobID));
mkdirOptional(sirfsstatedir);
fnames = getFileNamesFromDirectory(dmapDir,'types',{'.mat'});

fprintf('\n%%%%%%%%%%%% SIRFS %%%%%%%%%%%%\n');
%% Loading precomputed shape
for i=1:length(fnames)
%for i=1:length(fnames)
    sirfsdmapFile = fullfile(sirfsstatedir,fnames{i});
    if(exist(sirfsdmapFile,'file'))
        continue;
    end
    stateFile = fullfile(statesDir,fnames{i});
    state = load(stateFile);state=state.state;
    sirfsstate = getsirfsstate(state);
    savefunc(sirfsdmapFile,sirfsstate);
    %vis_PSIRFS(sirfsstate); pause
end
end

function savefunc(file,state)
    save(file,'state');
end

function sirfsstate = getsirfsstate(state)
    sigma = 3; % Controls the "bandwidth" of the input depth (the standard deviation of a Gaussian at which point the signal becomes reliable)
    mult = 5; % Controls the importance of the input depth (the multiplier on the loss)
    niters = 200;
    input_image = im2double(state.im);
    input_image(input_image<1/255) = 1/255;
    sirfsstate = SIRFS(input_image, (state.mask), [], ...
        ['params.DO_DISPLAY = 0; params.N_ITERS_OPTIMIZE = ' num2str(niters) ';']);
    sirfsstate.mask = state.mask;
    sirfsstate.im = input_image;
    sirfsstate.dmap = sirfsstate.height;
end

repo: CategoryShapes


text:
function tests = test_beylkin
  tests = functiontests(localfunctions);
end

function test_qmf(testCase)
    f = spx.wavelet.baylkin.quad_mirror_filter();
    verifyTrue(testCase, spx.norm.is_unit_norm_vec(f));
end

author: indigits

repo: sparse-plex


text:
function [ f1, f2, f3 ] = feature_all( img )

if(ndims(img)==3)
    img=rgb2gray(img);
end

im1=im2double(img);

% Scale factor is 3
h=fspecial('gaussian',3);
im_f=double(imfilter(im1,h));
im2 = im_f(2:2:end,2:2:end);

h=fspecial('gaussian',3);
im_f=double(imfilter(im2,h));
im3 = im_f(2:2:end,2:2:end);

t1=block_dct(im1);
t2=block_dct(im2);
t3=block_dct(im3);

f1=[t1 t2 t3];

f2=global_gsm(img);

col=im2col(im1,[5 5],'distinct');
t1=svd(col);
col=im2col(im2,[5 5],'distinct');
t2=svd(col);
col=im2col(im3,[5 5],'distinct');
t3=svd(col);
f3=[t1 t2 t3];

end

repo: sr-metric


text:
function F = triangles_from_strip(S)
  % TRIANGLES_FROM_STRIP Create a list of triangles from a stream of indices
  % along a strip.
  %
  % Inputs:
  %   S  #S list of indices
  % Outputs:
  %   F  #S-2 by 3 list of triangle indices
  %

  %F = zeros(numel(S)-2,3);
  %for s = 3:numel(S)
  %  if mod(s,2) == 0
  %    F(s-2,:) = fliplr([S(s-2) S(s-1) S(s)]);
  %  else
  %    F(s-2,:) = [S(s-2) S(s-1) S(s)];
  %  end
  %end

  S = S(:);
  % Maintain order and orientations
  F = [ ...
    S(1:2:end-2) S(2:2:end-1) S(3:2:end) ...
    S(3:2:end-1) S(2:2:end-2) S(4:2:end)];
  F = reshape(F',3,numel(S)-2)';

end

author: alecjacobson

repo: gptoolbox


text:
function [params, s_new] = denoise_pow3(params, s, state)
% Kurtosis based denoising function
%   [params, s_new] = denoise_pow3(params, s, state)
%     params  Function specific modifiable parameters
%     state   DSS algorithm state
%     s       Source signal estimate, matrix of row vector signals
%     s_new   Denoised signal estimate
%
%   Calculates third power of the signal. Equals kurtosis if source
%   signal has unit variance.

% Copyright (C) 2004, 2005 DSS MATLAB package team (dss@cis.hut.fi).
% Distributed by Laboratory of Computer and Information Science,
% Helsinki University of Technology. http://www.cis.hut.fi/projects/dss/.
% $Id$

if nargin<3 | ~isstruct(state)
    params.name = 'Kurtosis based denoising';
    params.description = '';
    params.param = {};
    params.param_value ={};
    params.param_type = {};
    params.param_desc = {};
    params.param_type = {};
    params.approach = {'defl','symm'};
    params.alpha = {};
    params.beta = {'beta_global','beta_pow3'};
    return;
end

s_new = s.^3;


text:
function [EEG g] = pre_selectComps(varargin)
%
% Remove some components from icaact matrix and dipfit structure.
%
% Inputs:
%
%   EEG:        EEG data structure
%
% Optional:     <'Name',Value> pairs
%
%     VerbosityLevel:   Verbosity level. 0 = no output, 1 = text, 2 = graphical
%                       Possible values: 0,1,2
%                       Default value  : 0
%                       Input Data Type: real number (double)
%
%     ComponentsToKeep: Select components to analyze
%                       This should be a cell array of strings containing the IDs of components you wish to keep
%                       Input Data Type: boolean
% Outputs:
%
%   EEG:        processed EEG structure
%   g:          argument specification structure
%
% See Also: pre_prepData()
%
% References:
%
% [1] Mullen T (2010) The Source Information Flow Toolbox (SIFT):
%   Theoretical Handbook and User Manual. Section 6.5.1
%   Available at: http://www.sccn.ucsd.edu/wiki/Sift
%
% Author: Tim Mullen 2010, SCCN/INC, UCSD.
% Email:  tim@sccn.ucsd.edu

% This function is part of the Source Information Flow Toolbox (SIFT)
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% extract the allowable component names
CAT = arg_extract(varargin,'CAT',1);
MyComponentNames = CAT.curComponentNames;

% setup the argument structure
g = arg_define([0 1],varargin, ...
    arg_norep('CAT',mandatory), ...
    arg({'verb','VerbosityLevel'},0,{int32(0) int32(1) int32(2)},'Verbosity level. 0 = no output, 1 = text, 2 = graphical'), ...
    arg({'ComponentsToKeep'},MyComponentNames,MyComponentNames,'Select components to analyze. This should be a cell array of strings containing the IDs of channels or components you wish to keep','type','logical') ...
    );

nbcomps = size(g.EEG.CAT.srcdata,1);
ComponentIndicesToKeep = find(ismember_bc(MyComponentNames,g.ComponentsToKeep));

if g.verb, fprintf('Selecting components...\n'); end

% select a subset of components from g.EEG
rmcomps = setdiff_bc(1:nbcomps,ComponentIndicesToKeep);

g.EEG.CAT.srcdata(rmcomps,:)     = [];
g.EEG.CAT.srcdata = reshape(g.EEG.CAT.srcdata, size(g.EEG.CAT.srcdata,1), g.EEG.CAT.pnts, g.EEG.CAT.trials);

g.EEG.CAT.curComps            = ComponentIndicesToKeep;
g.EEG.CAT.curComponentNames   = g.ComponentsToKeep;

EEG = g.EEG;

repo: REST


text:
function test_ipopt

  auxdata = {} ;

  options.lb = [ -Inf, -Inf ] ;  % Lower bound on the variables.
  options.ub = [  Inf, Inf ] ;  % Upper bound on the variables.

  % The constraint functions are bounded to zero
  options.cl = [ 0, 0, -Inf ]; %  constraints
  options.cu = [ Inf, Inf, 0];

  % Set up the auxiliary data.
  options.auxdata = auxdata ;

  % Set the IPOPT options.
  options.ipopt.jac_d_constant   = 'no';
  options.ipopt.hessian_constant = 'no';
  options.ipopt.mu_strategy      = 'adaptive';
  options.ipopt.max_iter         = 400;
  options.ipopt.tol              = 1e-10;
  options.ipopt.linear_solver = 'ma57';

  % The callback functions.
  funcs.objective         = @objective;
  funcs.constraints       = @constraints;
  funcs.gradient          = @gradient;
  funcs.jacobian          = @jacobian;
  funcs.jacobianstructure = @jacobianstructure;
  if false
    options.ipopt.derivative_test = 'first-order';
    funcs.hessian           = @hessian;
    funcs.hessianstructure  = @hessianstructure;
  else
    options.ipopt.hessian_approximation      = 'limited-memory';
    %options.ipopt.limited_memory_update_type = 'bfgs' ; % {bfgs}, sr1 = 6; % {6}
    %options.ipopt.limited_memory_update_type = 'sr1' ;
    options.ipopt.limited_memory_update_type = 'bfgs' ; % {bfgs}, sr1 = 6; % {6}
  end

  % Run IPOPT.
  x0 = [-2, 1] ;

  tic
  [x, info] = ipopt_auxdata(x0,funcs,options);
  elapsed = toc ;

  info;

  x

end

%%
% map the indices with the corresponding index in the spase matrix
function f = objective(x,auxdata)
  f = 100*(x(2)-x(1)^2)^2+(1-x(1))^2 ;
end

%%
% map the indices with the corresponding index in the spase matrix
function g = gradient(x,auxdata)
  g = [ 400*x(1)*(x(1)^2-x(2))+2*x(1)-2, ...
        200*(x(2)-x(1)^2) ] ;
end

function f = constraints(x,auxdata)
  f = zeros(3,1) ;
  f(1) = x(1)*x(2)-1 ; % = 0
  f(2) = x(1) + x(2)^2 ; % >= 0
  f(3) = x(1) ;
end

function jac = jacobian(x,auxdata)
  jac = sparse([ x(2), x(1)   ; ...
                 1,    2*x(2) ; ...
                 1,    0 ]) ;
end

function jac = jacobianstructure(auxdata)
  jac = sparse(ones(3,2)) ;
end

function H = hessian(x, sigma, lambda, auxdata)
  H = sigma * [ 1200*x(1)^2-400*x(2)+2, 0 ; ...
                 -400*x(1)               200] ;
  H = H + lambda(1) * [ 0 0 ; 1 0 ] ;
  H = H + lambda(2) * [ 0 0 ; 0 2 ] ;
  H = sparse(H) ;
end

function H = hessianstructure(auxdata)
  H = sparse([ 1 0 ; 1 1]) ;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
classdef RMMEDA_F6 < PROBLEM
% <multi> <real> <large/none>
% Benchmark MOP for testing RM-MEDA

%------------------------------- Reference --------------------------------
% Q. Zhang, A. Zhou, and Y. Jin, RM-MEDA: A regularity model-based
% multiobjective estimation of distribution algorithm, IEEE Transactions on
% Evolutionary Computation, 2008, 12(1): 41-63.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            g = 1 + 9*mean((X(:,2:end).^2-repmat(X(:,1),1,size(X,2)-1)).^2,2);
            PopObj(:,1) = sqrt(X(:,1));
            PopObj(:,2) = g.*(1-(PopObj(:,1)./g).^2);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = 1 - R(:,1).^2;
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R = obj.GetOptimum(100);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% produces synthetic Omori aftershock sequence using given epicentres
% Omori parameters can be chosen arbitraily
%
% Samuel Neukomm, 27.2.2004

[filename,pathname] = uigetfile('*.mat','Load earthquake sequence');
do = ['load ' pathname filename]; eval(do)

lon = a.Longitude; lat = a.Latitude; mag = a.Magnitude; dep = a.Depth;
[m_main, main] = max(a.Magnitude);
if size(a,2) == 9
    date_matlab = datenum(a.Date.Year,a.Date.Month,a.Date.Day,a.Date.Hour,a.Date.Minute,zeros(size(a,1),1));
else
    date_matlab = datenum(a.Date.Year,a.Date.Month, a.Date.Day,a.Date.Hour,a.Date.Minute,a(:,10));
end

date_main = date_matlab(main);
time_aftershock = date_matlab-date_main;

l = time_aftershock(:) > 0;
t_aftershock = time_aftershock(l);
eqcatalogue = a.subset(l);

ncum = (1:length(t_aftershock))';

prompt = {'p value:','c value:','k value:'};
def = {'0.9','0.05','1000'};
dlgTitle = 'Choose Omori parameters';
lineNo = 1;
answer = inputdlg(prompt,dlgTitle,lineNo,def);

pv = str2double(answer{1});
cv = str2double(answer{2});
kv = str2double(answer{3});

tasnew = (cv^(1-pv)-ncum*(pv-1)/kv).^(1/(1-pv))-cv;

date_mat_new = tasnew + date_main;
[yr,mon,day,hr,mn,sec] = datevec(date_mat_new);

datum = [yr mon day hr mn sec];
yr = decyear(datum);

if size(a,2) == 9
    a = [a(main,:) 0; eqcatalogue(:,1) eqcatalogue(:,2) yr mon day eqcatalogue(:,6) eqcatalogue(:,7) hr mn sec];
else
    a = [a(main,1:10); eqcatalogue(:,1) eqcatalogue(:,2) yr mon day eqcatalogue(:,6) eqcatalogue(:,7) hr mn sec];
end

[filename, pathname] = uiputfile('*.mat', 'Save synthetic catalog as');

try
    save(fullfile(pathname, filename),'a','pv','cv','kv');
catch
    disp('failed to save'); %complain
end

author: CelsoReyes


text:
%% Pole Figure Tutorial
%
%%
% This tutorial explains the basic concepts for ananyzing x-ray, synchrotron
% and neutron diffraction pole figure data.
%
%% Import pole figure diffraction data
% Click on <matlab:import_wizard('PoleFigure') Import pole figure data> to
% start the import wizard which is a GUI leading you through the import of
% pole figure data. After finishing the wizard you will end up with a
% script similar to the following one.

% This script was automatically created by the import wizard. You should
% run the whole script or parts of it in order to import your data. There
% is no problem in making any changes to this script.
%
% *Specify Crystal and Specimen Symmetries*

% crystal symmetry for this ZnCuTi data is hexagonal.  Here we define the crystallographic unit cell and how it relates to cartesian xyz axes.
CS = crystalSymmetry('6/mmm', [2.633 2.633 4.8], 'X||a*', 'Y||b', 'Z||c');

% specimen symmetry tells MTEX if a certain symmetry should be present in the plotted pole figures.  The command used here selects triclinic, the most flexible option.
SS = specimenSymmetry('1');

% plotting convention
setMTEXpref('xAxisDirection','north');
setMTEXpref('zAxisDirection','outOfPlane');

%%
% *Specify File Names*

% path to files downloaded with the MTEX package
pname = [mtexDataPath filesep 'PoleFigure' filesep 'ZnCuTi' filesep];

% which pole figure files are to be imported
fname = {...
  [pname 'ZnCuTi_Wal_50_5x5_PF_002_R.UXD'],...
  [pname 'ZnCuTi_Wal_50_5x5_PF_100_R.UXD'],...
  [pname 'ZnCuTi_Wal_50_5x5_PF_101_R.UXD'],...
  [pname 'ZnCuTi_Wal_50_5x5_PF_102_R.UXD'],...
  };

% defocusing correction to compensate for the equipment-dependent loss of intensity at certain angles.
fname_def = {...
  [pname 'ZnCuTi_defocusing_PF_002_R.UXD'],...
  [pname 'ZnCuTi_defocusing_PF_100_R.UXD'],...
  [pname 'ZnCuTi_defocusing_PF_101_R.UXD'],...
  [pname 'ZnCuTi_defocusing_PF_102_R.UXD'],...
  };

%%
% *Specify Miller Indices*

% These correspond to the files loaded, in order.
h = { ...
  Miller(0,0,2,CS),...
  Miller(1,0,0,CS),...
  Miller(1,0,1,CS),...
  Miller(1,0,2,CS),...
  };

%%
% *Import the Data*

% create a Pole Figure variable containing the data
pf = PoleFigure.load(fname,h,CS,SS,'interface','uxd');

% create a defocusing pole figure variable
pf_def = PoleFigure.load(fname_def,h,CS,SS,'interface','uxd');

% correct data by applying the defocusing compensation
pf = correct(pf,'def',pf_def);

%%
% After running the script the variable |pf| is created which contains all
% information about the pole figure data. You may plot the data using the
% command <PoleFigure.plot.html plot>

plot(pf)

%%
% By default pole figures are plotted as intensity-colored dots for every
% data point. There are many options to specify the way pole figures are
% plotted in MTEX. Have a look at the <PoleFigurePlot.html plotting
% section> for more information.
%
% After import make sure that the Miller indices are correctly assigned to
% the pole figures and that the alignment of the specimen coordinate
% system, i.e., X, Y, Z is correct. In case of outliers or misaligned data,
% you may want to correct your raw data. Have a look at the
% <PoleFigureCorrection.html correction section> for further information.
% MTEX offers several methods correcting pole figure data, e.g.
%
% * rotating pole figures
% * scaling pole figures
% * finding outliers
% * removing specific measurements
% * superposing pole figures
%
% As an example we set all negative intensities to zero

pf(pf.intensities<0) = 0;
plot(pf)

%% ODF Estimation
%
% Once your data is in good shape, i.e. defocusing correction has been
% done and few outliers are left you can reconstruct an ODF out of
% this data. This is done by the command <PoleFigure.calcODF.html
% calcODF>.

odf = calcODF(pf,'silent')

%%
% Note that reconstructing an ODF from pole figure data is a severely ill-
% posed problem, i.e., it does *not* provide a unique solution. A more
% through discussion on the ambiguity of ODF reconstruction from
% pole figure data can be found <PoleFigure2ODFAmbiguity.html here>. As a
% rule of thumb: the more pole figures you have and the more consistent your
% pole figure data the better your reconstructed ODF will be.
%
% To check how well your reconstructed ODF fits the measured pole figure
% data use

figure;plotPDF(odf,pf.h)

%%
% Compare the recalculated pole figures with the measured data.
% A quantitative measure for the fitting is the so called RP value. They
% can be computed for each imported pole figure with

calcError(odf,pf)

%%
% In the case of a bad fit, you may want to tweak the reconstruction
% algorithm. See <PoleFigure2ODF.html here> for more information.

%% Visualize the ODF
% Finally one can plot the resulting ODF

plot(odf)
mtexColorMap LaboTeX

author: mtex-toolbox

repo: mtex

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


text:
function Z = projectData(X, U, K)
%PROJECTDATA Computes the reduced data representation when projecting only
%on to the top k eigenvectors
%   Z = projectData(X, U, K) computes the projection of
%   the normalized inputs X into the reduced dimensional space spanned by
%   the first K columns of U. It returns the projected examples in Z.
%

% You need to return the following variables correctly.
Z = zeros(size(X, 1), K);

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the projection of the data using only the top K
%               eigenvectors in U (first K columns).
%               For the i-th example X(i,:), the projection on to the k-th
%               eigenvector is given as follows:
%                    x = X(i, :)';
%                    projection_k = x' * U(:, k);
%

subsetU = U(:, 1:K);
Z = X * subsetU;

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

end

author: vkosuri

repo: CourseraMachineLearning


text:
function pass = test_coeffs2vals(pref)
% test various coeffs2vals, vals2coeffs, vals2vals, and coeffs2coeffs codes.
if ( nargin == 0 )
    pref = chebfunpref();
end

tol = 1000*eps;

f = chebfun(@exp, pref);
N = length(f);
c_leg = legcoeffs(f);
c_cheb = chebcoeffs(f);
v_leg = f(legpts(N));
v_cheb1 = f(chebpts(N,1));
v_cheb2 = f(chebpts(N,2));

pass(1) = norm(c_leg - chebcoeffs2legcoeffs(c_cheb), inf) < tol;
pass(2) = norm(c_leg - legvals2legcoeffs(v_leg), inf) < tol;
pass(3) = norm(c_leg - chebvals2legcoeffs(v_cheb1, 1), inf) < tol;
pass(4) = norm(c_leg - chebvals2legcoeffs(v_cheb2, 2), inf) < tol;

pass(5) = norm(c_cheb - legcoeffs2chebcoeffs(c_leg), inf) < tol;
pass(6) = norm(c_cheb - legvals2chebcoeffs(v_leg), inf) < tol;
pass(7) = norm(c_cheb - chebvals2chebcoeffs(v_cheb1, 1), inf) < tol;
pass(8) = norm(c_cheb - chebvals2chebcoeffs(v_cheb2, 2), inf) < tol;

pass(9) = norm(v_leg - legcoeffs2legvals(c_leg), inf) < tol;
pass(10) = norm(v_leg - chebcoeffs2legvals(c_cheb), inf) < tol;
pass(11) = norm(v_leg - chebvals2legvals(v_cheb1, 1), inf) < tol;
pass(12) = norm(v_leg - chebvals2legvals(v_cheb2, 2), inf) < tol;

pass(13) = norm(v_cheb1 - legcoeffs2chebvals(c_leg, 1), inf) < tol;
pass(14) = norm(v_cheb1 - chebcoeffs2chebvals(c_cheb,1), inf) < tol;
pass(15) = norm(v_cheb1 - legvals2chebvals(v_leg, 1), inf) < tol;
pass(16) = norm(v_cheb1 - chebvals2chebvals(v_cheb2, 2, 1), inf) < tol;

pass(17) = norm(v_cheb2 - legcoeffs2chebvals(c_leg, 2), inf) < tol;
pass(18) = norm(v_cheb2 - chebcoeffs2chebvals(c_cheb,2), inf) < tol;
pass(19) = norm(v_cheb2 - legvals2chebvals(v_leg, 2), inf) < tol;
pass(20) = norm(v_cheb2 - chebvals2chebvals(v_cheb1, 1, 2), inf) < tol;

end


text:
%kkinsertseg 'Insert Selected Data Segments from Input 1 into Input 2'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kinsertseg.pane file
%
% Parameters:
% InputFile: i1 'Input 1', required: 'Source Input File'
% InputFile: i2 'Input 2', required: 'Destination Input File'
% OutputFile: o 'Output ', required: 'output file'
% Toggle: val 'Value', default: 0: 'If -val flag is selected, operate on Value segment'
% Toggle: map 'Map', default: 0: 'If -map flag is selected, operate on map segment'
% Toggle: mask 'Mask', default: 0: 'If -mask flag is selected, operate on Mask segment'
% Toggle: loc 'Location', default: 0: 'If -loc flag is selected, operate on Location segment'
% Toggle: time 'Time', default: 0: 'If -time flag is selected, operate on Time segment'
%
% Example: o = kkinsertseg({i1, i2}, {'i1','';'i2','';'o','';'val',0;'map',0;'mask',0;'loc',0;'time',0})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kinsertseg - Insert Selected Data Segments from Input 1 into Input 2
%
%  DESCRIPTION
% This program will insert the specified
% segment(s) for the first input object into
% the second input object producing a new
% or output object.  If -segment is
% "Insert All" all the segments and
% their associated attributes will be inserted.
% The -segment parameter can also specify
% multiple segments by typing a string
% containing the segment names separated
% by space(s):
%
%                "segment1 segment2 segment3"
%
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkinsertseg(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkinsertseg(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i1', '__input';'i2', '__input';'o', '__output';'val', 0;'map', 0;'mask', 0;'loc', 0;'time', 0};
maxval={0,0,0,0,0,0,0,0};
minval={0,0,0,0,0,0,0,0};
istoggle=[0,0,0,1,1,1,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','InputFile','OutputFile','Toggle','Toggle','Toggle','Toggle','Toggle'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'kinsertseg"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
% Test file for chebtech/cumsum.m

function pass = test_cumsum(pref)

% Get preferences.
if ( nargin < 1 )
    pref = chebtech.techPref();
end

% Generate a few random points to use as test values.
seedRNG(6178);
x = 2 * rand(100, 1) - 1;

for n = 1:2
    if ( n == 1 )
        testclass = chebtech1();
    else
        testclass = chebtech2();
    end

  %%
  % Spot-check antiderivatives for a couple of functions.  We verify that the
  % chebtech antiderivatives match the true ones up to a constant by checking
  % that the standard deviation of the difference between the two on a large
  % random grid is small. We also check that feval(cumsum(f), -1) == 0 each
  % time.

  f = testclass.make(@(x) exp(x) - 1, [],  pref);
  F = cumsum(f);
  F_ex = @(x) exp(x) - x;
  err = std(feval(F, x) - F_ex(x));
  tol = 20*vscale(F)*eps;
  pass(n, 1) = (err < tol) && (abs(feval(F, -1)) < tol);

  f = testclass.make(@(x) 1./(1 + x.^2), [], pref);
  F = cumsum(f);
  F_ex = @(x) atan(x);
  err = feval(F, x) - F_ex(x);
  tol = 10*vscale(F)*eps;
  pass(n, 2) = (std(err) < tol) && (abs(feval(F, -1)) < tol);

  f = testclass.make(@(x) cos(1e4*x), [], pref);
  F = cumsum(f);
  F_ex = @(x) sin(1e4*x)/1e4;
  err = feval(F, x) - F_ex(x);
  tol = 5e4*vscale(F)*eps;
  pass(n, 3) = (std(err) < tol) && (abs(feval(F, -1)) < tol);

  z = exp(2*pi*1i/6);
  f = testclass.make(@(t) sinh(t*z), [], pref);
  F = cumsum(f);
  F_ex = @(t) cosh(t*z)/z;
  err = feval(F, x) - F_ex(x);
  tol = 10*vscale(F)*eps;
  pass(n, 4) = (std(err) < tol) && (abs(feval(F, -1)) < tol);

  %%
  % Check that applying cumsum() and direct construction of the antiderivative
  % give the same results (up to a constant).

  f = testclass.make(@(x) sin(4*x).^2, [], pref);
  F = testclass.make(@(x) 0.5*x - 0.0625*sin(8*x), [], pref);
  G = cumsum(f);
  err = G - F;
  tol = 10*vscale(G)*eps;
  values = err.coeffs2vals(err.coeffs);
  pass(n, 5) = (std(values) < tol) && (abs(feval(G, -1)) < tol);

  %%
  % Check that diff(cumsum(f)) == f and that cumsum(diff(f)) == f up to a
  % constant.

  f = testclass.make(@(x) x.*(x - 1).*sin(x) + 1, [], pref);
  g = diff(cumsum(f));
  err = feval(f, x) - feval(g, x);
  tol = 10*vscale(g)*eps;
  pass(n, 6) = (norm(err, inf) < 100*tol);
  h = cumsum(diff(f));
  err = feval(f, x) - feval(h, x);
  tol = 10*vscale(h)*eps;
  pass(n, 7) = (std(err) < tol)  && (abs(feval(h, -1)) < tol);

  %%
  % Check operation for array-valued chebtech objects.

  f = testclass.make(@(x) [sin(x) x.^2 exp(1i*x)], [], pref);
  F_exact = testclass.make(@(x) [(-cos(x)) (x.^3/3) (exp(1i*x)/1i)], [], pref);
  F = cumsum(f);
  err = std(feval(F, x) - feval(F_exact, x));
  tol = 10*max(vscale(F)*eps);
  pass(n, 8) = (norm(err, inf) < tol)  && all(abs(feval(F, -1)) < tol);

end

end


text:
function [bnet, Qnodes, Fnodes, Onode] = mk_hhmm(varargin)
% MK_HHMM Make a Hierarchical HMM
% function [bnet, Qnodes, Fnodes, Onode] = mk_hhmm(...)
%
% e.g. 3-layer hierarchical HMM where level 1 only connects to level 2
% and the parents of the observed node are levels 2 and 3.
% (This DBN is the same as Fig 10 in my tech report.)
%
%   Q1 ---------->   Q1
%   |  \            ^ |
%   |   v          /  |
%   |    F2 ------/   |
%   |   ^ ^        \  |
%   |  /  |         \ |
%   | /   |          ||
%   v     |          vv
%   Q2----| --------> Q2
%  /| \   |          ^|
% / |  v  |         / |
% | |   F3 --------/  |
% | |   ^          \  |
% | v  /            v v
% | Q3 ----------->  Q3
% |  |
% \  |
%  v v
%   O
%
%
% Optional arguments in name/value format [default value in brackets]
%
% Qsizes      - sizes at each level [ none ]
% allQ       - 1 means level i connects to all Q levels below, 0 means just to i+1 [0]
% transprob  - transprob{d}(i,k,j) = P(Q(d,t)=j|Q(d,t-1)=i,Q(1:d-1,t)=k)  ['leftright']
% startprob  - startprob{d}(k,j) = P(Q(d,t)=j|Q(1:d-1,t)=k)  ['leftstart']
% termprob   - termprob{d}(k,j) = P(F(d,t)=2|Q(1:d-1,t)=k,Q(d,t)=j) for d>1 ['rightstop']
% selfprop   - prob of a self transition (termprob default = 1-selfprop) [0.8]
% Osize       - size of O node
% discrete_obs - 1 means O is tabular_CPD, 0 means gaussian_CPD [0]
% Oargs       - cell array of args to pass to the O CPD  [ {} ]
% Ops         - Q parents of O [Qnodes(end)]
% F1          - 1 means level 1 can finish (restart), else there is no F1->Q1 arc [0]
% clamp1     - 1 means we clamp the params of the Q nodes in slice 1 (Qt1params) [1]
%   Note: the Qt1params are startprob, which should be shared with other slices.
%   However, in the current implementation, the Qt1params will only be estimated
%   from the initial state of each sequence.
%
% For d=1, startprob{1}(1,j) is only used in the first slice and
% termprob{1} is ignored, since we assume the top level never resets.
% Also, transprob{1}(i,j) can be used instead of transprob{1}(i,1,j).
%
% leftstart means the model always starts in state 1.
% rightstop means the model can only finish in its last state (Qsize(d)).
% unif means each state is equally like to reach any other
% rnd means the transition/starting probs are random (drawn from rand)
%
% Q1:QD in slice 1 are of type tabular_CPD
% Q1:QD in slice 2 are of type hhmmQ_CPD.
% F(2:D-1) is of type hhmmF_CPD, FD is of type tabular_CPD.

args = varargin;
nargs = length(args);

% get sizes of nodes and topology
Qsizes = [];
Osize = [];
allQ = 0;
Ops = [];
F1 = 0;
for i=1:2:nargs
  switch args{i},
   case 'Qsizes', Qsizes = args{i+1};
   case 'Osize',  Osize = args{i+1};
   case 'allQ',   allQ = args{i+1};
   case 'Ops',    Ops = args{i+1};
   case 'F1',     F1 = args{i+1};
  end
end
if isempty(Qsizes), error('must specify Qsizes'); end
if Osize==0, error('must specify Osize'); end
D = length(Qsizes);
Qnodes = 1:D;

if isempty(Ops), Ops = Qnodes(end); end

[intra, inter, Qnodes, Fnodes, Onode] = mk_hhmm_topo(D, allQ, Ops, F1);
ss = length(intra);
names = {};

if F1
  Fnodes_ndx = Fnodes;
else
  Fnodes_ndx = [-1 Fnodes]; % Fnodes(1) is a dummy index
end

% set default params
discrete_obs = 0;
Oargs = {};
startprob = cell(1,D);
startprob{1} = 'unif';
for d=2:D
  startprob{d} = 'leftstart';
end
transprob = cell(1,D);
transprob{1} = 'unif';
for d=2:D
  transprob{d} = 'leftright';
end
termprob = cell(1,D);
for d=2:D
  termprob{d} = 'rightstop';
end
selfprob = 0.8;
clamp1 = 1;

for i=1:2:nargs
  switch args{i},
   case 'discrete_obs', discrete_obs = args{i+1};
   case 'Oargs',        Oargs = args{i+1};
   case 'startprob',    startprob = args{i+1};
   case 'transprob',    transprob = args{i+1};
   case 'termprob',     termprob = args{i+1};
   case 'selfprob',     selfprob = args{i+1};
   case 'clamp1',       clamp1 = args{i+1};
  end
end

ns = zeros(1,ss);
ns(Qnodes) = Qsizes;
ns(Onode) = Osize;
ns(Fnodes) = 2;

dnodes = [Qnodes Fnodes];
if discrete_obs
  dnodes = [dnodes Onode];
end
onodes = [Onode];

bnet = mk_dbn(intra, inter, ns, 'observed', onodes, 'discrete', dnodes, 'names', names);
eclass = bnet.equiv_class;

for d=1:D
  if d==1
    Qps = [];
  elseif allQ
    Qps = Qnodes(1:d-1);
  else
    Qps = Qnodes(d-1);
  end
  Qpsz = prod(ns(Qps));
  Qsz = ns(Qnodes(d));
  if isstr(startprob{d})
    switch startprob{d}
     case 'unif', startprob{d} = mk_stochastic(ones(Qpsz, Qsz));
     case 'rnd', startprob{d} = mk_stochastic(rand(Qpsz, Qsz));
     case 'leftstart', startprob{d} = zeros(Qpsz, Qsz); startprob{d}(:,1) = 1;
    end
  end
  if isstr(transprob{d})
    switch transprob{d}
     case 'unif', transprob{d} = mk_stochastic(ones(Qsz, Qpsz, Qsz));
     case 'rnd', transprob{d} = mk_stochastic(rand(Qsz, Qpsz, Qsz));
     case 'leftright',
      LR = mk_leftright_transmat(Qsz, selfprob);
      temp = repmat(reshape(LR, [1 Qsz Qsz]), [Qpsz 1 1]); % transprob(k,i,j)
      transprob{d} = permute(temp, [2 1 3]); % now transprob(i,k,j)
    end
  end
  if isstr(termprob{d})
    switch termprob{d}
     case 'unif', termprob{d} = mk_stochastic(ones(Qpsz, Qsz, 2));
     case 'rnd', termprob{d} = mk_stochastic(rand(Qpsz, Qsz, 2));
     case 'rightstop',
      %termprob(k,i,t) Might terminate if i=Qsz; will not terminate if i<Qsz
      stopprob = 1-selfprob;
      termprob{d} = zeros(Qpsz, Qsz, 2);
      termprob{d}(:,Qsz,2) = stopprob;
      termprob{d}(:,Qsz,1) = 1-stopprob;
      termprob{d}(:,1:(Qsz-1),1) = 1;
     otherwise, error(['unrecognized termprob ' termprob{d}])
    end
  elseif d>1 % passed in termprob{d}(k,j)
    temp = termprob{d};
    termprob{d} = zeros(Qpsz, Qsz, 2);
    termprob{d}(:,:,2) = temp;
    termprob{d}(:,:,1) = ones(Qpsz,Qsz) - temp;
  end
end

% SLICE 1

for d=1:D
  bnet.CPD{eclass(Qnodes(d),1)} = tabular_CPD(bnet, Qnodes(d), 'CPT', startprob{d}, 'adjustable', clamp1);
end

if F1
  d = 1;
  bnet.CPD{eclass(Fnodes_ndx(d),1)}  = hhmmF_CPD(bnet, Fnodes_ndx(d), Qnodes(d), Fnodes_ndx(d+1), ...
						 'termprob', termprob{d});
end
for d=2:D-1
  if allQ
    Qps = Qnodes(1:d-1);
  else
    Qps = Qnodes(d-1);
  end
  bnet.CPD{eclass(Fnodes_ndx(d),1)}  = hhmmF_CPD(bnet, Fnodes_ndx(d), Qnodes(d), Fnodes_ndx(d+1), ...
						 'Qps', Qps, 'termprob', termprob{d});
end
bnet.CPD{eclass(Fnodes_ndx(D),1)}  = tabular_CPD(bnet, Fnodes_ndx(D), 'CPT', termprob{D});

if discrete_obs
  bnet.CPD{eclass(Onode,1)} = tabular_CPD(bnet, Onode, Oargs{:});
else
  bnet.CPD{eclass(Onode,1)} = gaussian_CPD(bnet, Onode, Oargs{:});
end

% SLICE 2

%for d=1:D
%  bnet.CPD{eclass(Qnodes(d),2)} = hhmmQ_CPD(bnet, Qnodes(d)+ss, Qnodes, d, D, ...
%					    'startprob', startprob{d}, 'transprob', transprob{d}, ...
%					    'allQ', allQ);
%end

d = 1;
if F1
  bnet.CPD{eclass(Qnodes(d),2)} = hhmmQ_CPD(bnet, Qnodes(d)+ss, 'Fself', Fnodes_ndx(d), ...
					    'Fbelow', Fnodes_ndx(d+1), ...
					    'startprob', startprob{d}, 'transprob', transprob{d});
else
    bnet.CPD{eclass(Qnodes(d),2)} = hhmmQ_CPD(bnet, Qnodes(d)+ss, ...
					    'Fbelow', Fnodes_ndx(d+1), ...
					    'startprob', startprob{d}, 'transprob', transprob{d});
end
for d=2:D-1
  if allQ
    Qps = Qnodes(1:d-1);
  else
    Qps = Qnodes(d-1);
  end
  Qps = Qps + ss; % since all in slice 2
  bnet.CPD{eclass(Qnodes(d),2)} = hhmmQ_CPD(bnet, Qnodes(d)+ss, 'Fself', Fnodes_ndx(d), ...
					    'Fbelow', Fnodes_ndx(d+1), 'Qps', Qps, ...
					    'startprob', startprob{d}, 'transprob', transprob{d});
end
d = D;
if allQ
  Qps = Qnodes(1:d-1);
else
  Qps = Qnodes(d-1);
end
Qps = Qps + ss; % since all in slice 2
bnet.CPD{eclass(Qnodes(d),2)} = hhmmQ_CPD(bnet, Qnodes(d)+ss, 'Fself', Fnodes_ndx(d), ...
					  'Qps', Qps, ...
					  'startprob', startprob{d}, 'transprob', transprob{d});

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function x = p02_start ( i, x )

%*****************************************************************************80
%
%% P02_START returns a starting point for problem 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the index of the starting point.
%
%    Output, real X, the starting point.
%
  if ( i == 1 )
    x = 0.0;
  elseif ( i == 2 )
    x = 1.0;
  elseif ( i == 3 )
    x = -5.0;
  elseif ( i == 4 )
    x = 10.0;
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P02_START - Fatal error!\n' );
    fprintf ( 1, '  Illegal starting point index = %d\n', i );
    error ( 'P02_START - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [int, keepindex, repindex] = mesh_laplacian_interp(lap, index);

% MESH_LAPLACIAN_INTERP: Computes the zero Laplacian interpolation matrix
%
% Useage:   [int, keepindex, repindex] = mesh_laplacian_interp(lap, index)
%
% This function calculates an interpolation matrix that provides
% the coefficients for the calculation of potential values at all
% unknown vertices of a mesh, given known potential values at
% a subset of the mesh vertices (at 'index').  The interpolation
% solution is constrained by a minimal norm of the Laplacian
% of the mesh.  See the reference below for details.
%
% 'lap' is the laplacian matrix for the full mesh (see mesh_laplacian)
% 'int' is the matrix which interpolates from the points in 'index'
% to the full mesh.  'index' is a row vector of indices into a
% subset of the vertices used to calculate 'lap'.  This subset
% is where the electric potential is known and usually corresponds
% to the given electrode vertices, eg:
%
% index = dsearchn(scalpvert,elecvert)';
%
% If 'index' contains repeated indices, only the unique indices
% are useful.  The 'keepindex' array can be used to select these.
% The 'repindex' array is the repeated indices.
%
% Interpolations can be done using matrix 'int', eg:
%
% [int, keepindex, repindex] = mesh_laplacian_interp(lap,index);
% if isempty(repindex),
%   Vint = int * Vknown;
% else
%   Vint = int * Vknown(keepindex);
% end
%
% This implements interpolation method B (p. 336) of
% Oostendorp T, Oosterom A & Huiskamp G (1989),
% Interpolation on a triangulated 3D surface.
% Journal of Computational Physics, 80: 331-343.
%

% Licence:  GNU GPL, no implied or express warranties
% History:  (c) 04/2002 Robert Oostenveld
%           - agreed to release 'lapint' under GNU GPL
%           04/2002, Darren.Weber@flinders.edu.au
%           - introduced check for index replications and
%             adjusted calculations/output accordingly
%           - converted lap to sparse matrix and solution
%             of interpolation matrix with \ operator
%           - accepts sparse lap input and returns sparse int
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if size(lap,1)~=size(lap,2), error('MESH_LAPLACIAN_INTERP: lap matrix is not square'); end

if issparse(lap), lap = full(lap); end

tic

if size(index,1)~=1, index = index'; end

% Remove any replicate indices from 'index'
[KnownIndex, i, i] = union(index,index);
if ~isequal(length(KnownIndex),length(index)),
    fprintf('\a\n\nMESH_LAPLACIAN_INTERP: Warning:\n\nTrimming duplicate values from index, use keepindex!\n\n');
end
keepindex = sort(i);
repindex = setdiff(1:length(index),sort(i));

KnownIndex = index(keepindex); % unsort KnownIndex
clear index

k = length(KnownIndex);
n = length(lap);

fprintf('MESH_LAPLACIAN_INTERP: Calc Interpolation matrix for %5d to %5d vertices...',k,n);

% find 'unknown' indices of lap matrix
UnknownIndex = setdiff(1:n, KnownIndex);

% reshuffle rows & columns of lap matrix
lapi = [KnownIndex, UnknownIndex];
lap = lap(lapi, :); % rows
lap = lap(:, lapi); % columns

% Segregate known/unknown portions of lap
k = length(KnownIndex);
n = length(lap);

L11 = lap(1:k    ,1:k    );
L12 = lap(1:k    ,(k+1):n);
L21 = lap((k+1):n,1:k    );
L22 = lap((k+1):n,(k+1):n);

clear lap lapi; % tidy up some memory

% Convert to sparse for quicker computation
A = sparse([L12; L22]);
B = sparse([L11; L21]);

clear L11 L12 L21 L22;  % tidy up some memory

%int = -pinv(A) * B;  % cannot use pinv with sparse matrix
int = -A \ B;

% Convert result back to full matrix
int = full(int);

% append the interpolating piece to the identity matrix
% these take care of the known potentials
int = [eye(k); int];

% reshuffle the columns of the interpolating matrix
[tmp, order] = sort(KnownIndex);
int = int(:,order);

% reshuffle the rows of the interpolating matrix
[tmp, order] = sort([KnownIndex, UnknownIndex]);
int = int(order, :);

int = sparse(int);

t = toc;
fprintf('done (%6.2f sec).\n',t);

return


text:
function [node,elem,interface] = interfacemeshdoc(box,phi,h)
%% INTERFACEMESHDOC doc for interfacemesh
%
% [node,elem,interfaceData] = INTERFACEMESH(box,phi,h) generates an
% interface-fitted mesh [node,elem]. The box is a rectangle containing the
% interface and the interface is defined as the zero level set of function
% phi. The parameter h is the size of the backgroud uniform mesh.
%
% Information on the discrete interface is stored in the output
% interface which includes:
%   - interface.vSign: sign of phi at vertices;
%        1: outside;  -1: inside;  0: on the interface
%   - interface.tElem: traingles near the interface;
%   - interface.sElem: traingles away from the interface;
%   - interface.edge: edges approximating the interface;
%   - interface.node: vertices on the interface;
%
%  Examples
%   % circle
%     box = [ -1, 1, -1, 1];
%     h = 0.1;
%     phi = @(p) sum(p.^2, 2) - 0.5.^2;
%     [node,elem,interfaceData] = interfacemeshdoc(box,phi,h);
%
%    % heart
%     box = [ -1, 1, -1, 1];
%    [node,elem,interface] = interfacemeshdoc(box,@phiheart,0.05);
%    showmesh(node,elem);
%    findedge(node,interface.edge,'all','noindex','draw');
%
%
% Author: Huayi Wei <weihuayi@xtu.edu.cn>, and Long Chen.
%
%   Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Construct the initial structure mesh
[node,elem,T] = squarequadmesh(box,h);
edge = T.edge;
edge2elem = T.edge2elem;
clear T;
N = size(node,1);
NT = size(elem,1);
showmesh(node,elem);
pause;

%% Step 1: Find points near interface
% compute the phi value at each vertex
phiValue = phi(node);
phiValue(abs(phiValue) < eps*h) = 0;  % treat nearby nodes as on the interface
vSign = sign(phiValue);
findnode(node,vSign==0);

% Find the intersection points between edges and the interface
isCutEdge = (vSign(edge(:,1)).* vSign(edge(:,2))<0);
findedge(node,edge,isCutEdge,'noindex','draw');
A = node(edge(isCutEdge,1),:);
B = node(edge(isCutEdge,2),:);
cutNode = findintersectbisect(phi,A,B);
Ncut = size(cutNode, 1);
vSign(N+1:N+Ncut) = 0;
findnode(cutNode,'all','noindex','color','k','MarkerSize',12);
pause;

% find interface elem and nodes
isInterfaceElem = false(NT,1);
isInterfaceElem(edge2elem(isCutEdge,[1,2])) = true;
isInterfaceElem(sum(abs(vSign(elem)), 2) < 3) = true; % 2 vertices on interface
findelem(node,elem,isInterfaceElem,'noindex');

isInterfaceNode = false(N,1);
isInterfaceNode(elem(isInterfaceElem,:)) = true;
findnode(node,isInterfaceNode,'noindex');
pause;

% add centers of special elements (two vertices on the interface)
%  0 - -1    -1 - 0     0 -  0      1 - 0
%  1 -  0     0 - 1     1 - -1      1 - 0
isSpecialElem = (sum(abs(vSign(elem)),2) <= 2); % two or more vertices on interfaces
% isSpecialElem = (sum(vSign(elem),2) == 0) & (sum(abs(vSign(elem)),2) == 2);
specialElem = elem(isSpecialElem,:);
auxPoint = (node(specialElem(:,1),:) + node(specialElem(:, 3),:))/2.0;
Naux = size(auxPoint,1);
vSign(N+Ncut+1:N+Ncut+Naux) = sign(phi(auxPoint));
% find the first two cases
isDiagInterface = (vSign(specialElem(:,1)).*vSign(specialElem(:,3)) == -1) | ...
                  (vSign(specialElem(:,2)).*vSign(specialElem(:,4)) == -1);
vSign(N+Ncut+find(isDiagInterface)) = 0;
findnode(auxPoint,'all','noindex','color','m');
pause;

% add cut points and aux points
node = [node; cutNode; auxPoint];
nearInterfaceNode = [node(isInterfaceNode,:); cutNode; auxPoint];
interfaceNodeIdx = [find(isInterfaceNode); N+(1:Ncut)'; N+Ncut+(1:Naux)'];
% figure(2);
% showmesh(node,elem);
findnode(nearInterfaceNode,'all','noindex');
pause;

%% Step 2: generate a Delaunay triangulation of interface points
% construct the Delaunay triangulation of interfaceNode
% different versions of matlab using different delaunay triangulation
matlabversion = version();
if str2double(matlabversion(end-5:end-2)) <= 2013
    DT = DelaunayTri(nearInterfaceNode); %#ok<*DDELTRI>
    tElem = DT.Triangulation;
else
    DT = delaunayTriangulation(nearInterfaceNode);
    tElem = DT.ConnectivityList;
end
tElem = fixorder(nearInterfaceNode,tElem); % correct the orientation
% showmesh(node,tElem);

%% Step 3: Post-processing
% get rid of the unnecessary triangles
NI = sum(isInterfaceNode);  % number of pts of the background mesh near interface
haveNewPts = (sum(tElem > NI,2) > 0); % triangles containing new added vertices
% vSignNearInterface = zeros(size(nearInterfaceNode,1),1);
% vSignNearInterface(1:sum(isInterfaceNode)) = vSign(isInterfaceNode);
% haveInterfacePts = (sum(abs(vSignNearInterface(tElem)),2) < 3);
tElem = tElem(haveNewPts,:);
% tElem = tElem(haveInterfacePts,:);
tElem = interfaceNodeIdx(tElem);  % map interfaceNode index to node index
showmesh(node,tElem,'Facecolor','y');
findnode(node,vSign==0,'noindex');
pause;

% get the remainding quad elems
sElem = elem(~isInterfaceElem,:);
% merge into one triangulation
elem = [tElem; sElem(:,[2 3 1]); sElem(:,[4 1 3])];
T = auxstructure(tElem);
showmesh(node,elem); hold on;
showmesh(node,tElem,'Facecolor','y');
isInterfaceEdge = ((vSign(T.edge(:,1)) == 0) & vSign(T.edge(:,2)) == 0);
interfaceEdge = T.edge(isInterfaceEdge,:);
nearInterfaceNode = find(vSign == 0);
findedge(node,interfaceEdge,'all','noindex','draw');
pause;

%% Generate interfaceData
interface.vSign = vSign;
interface.tElem = tElem;
interface.sElem = sElem;
interface.edge = interfaceEdge;
interface.node = nearInterfaceNode;

repo: ifem


text:
function [a_opt,dt_opt] = refine(params,ref,tgt,varargin)
% function [anew,dtnew] = refine(params,ref,tgt,{cons,projmodel})
%
% Optimize initial estimates (from Hough accumulator) using non-linear
% minimization.
%
% For a variable alpha: params = [dt0,alpha0].
% For a fixed alpha: params = [dt0] and cons = [alpha].

% read in variable arguments
cons	= [];
	if (length(varargin)>0), cons = varargin{1}; end
projmodel	= 'affine';
	if (length(varargin)>1), projmodel = varargin{2}; end

% parameters+constraints must equal 2
if (size(params,2)+size(cons,2) ~= 2)
	error('Incorrect number of parameters/constraints');
end

% extract alpha and dt from params/cons
dt = params(:,1);
if isempty(cons),	a = params(:,2);
else,							a = cons(:,1);
end

% perform non-linear optimization of initial estimates using sub-frame cost
% function
cmax = inf;
for hyp	= 1:length(a)
	fprintf('  (a=%0.3f,dt=%0.3f) -> ',a(hyp),dt(hyp));
	if isempty(cons)
		[pmin,c]	=	fminsearch(@subframecost,[dt(hyp),a(hyp)],[],ref,tgt,[],projmodel);
		dtnew	= pmin(1);
		anew	= pmin(2);
	else
		[pmin,c]	=	fminsearch(@subframecost,[dt(hyp)],[],ref,tgt,[a(hyp)],projmodel);
		dtnew	= pmin(1);
		anew	= a(hyp);
	end
	fprintf('(a=%0.3f,dt=%0.3f)\n',anew,dtnew);

	if (c < cmax)
		cmax		= c;
		dt_opt	= dtnew;
		a_opt		= anew;
	end
end

repo: mcbench-benchmarks


text:
clear all;

f = @(p)[p-p^2+3 2*p]';

% sampling intervals for each parameter
domain = [-3 3];
% grid size: number of grid points for each parameter
gridsize = 19;

% sampling
D = sampling_vec(f, domain, gridsize);

% hosvd
[S U] = hosvd(D, [1 0]);
U = U{1};

[W V] = genhull(U,'box');
[Wc3 Vc3] = genhull(U);

close all
set(0, 'DefaultFigureWindowStyle', 'docked')
plothull(U);
plothull(W);
plothull(Wc3);

figure
hold on
plot(U(:,1), U(:,2))
plot(V(:,1), V(:,2), 'ko')
plot(Vc3(:,1), Vc3(:,2), 'ro')
plot(Vc3(:,1), Vc3(:,2), 'r-')
hold off


text:
%ISVALDFILE Test whether the argument is a valid datafile or dataset
%
% 	N = ISVALDFILE(A);
% 	N = ISVALDFILE(A,M);
% 	N = ISVALDFILE(A,M,C);
%
% INPUT
%   A     Input argument, to be tested on datafile or dataset
%   M     Minimum number of objects per class in A
%   C     Minimum number of classes in A
%
% OUTPUT
%   N  1/0 if A is / isn't a valid datafile or dataset
%
% DESCRIPTION
% The function ISVALDFILE tests if A is a datafile or dataset that has
% at least C classes. It is an extension of ISVALDSET and can be used it
% when datasets as well as datafiles are allowed.
% For datafiles(sets) with soft or targets labels it is tested whether A
% has at least M objects.
%
% SEE ALSO (PRTools Guide)
% ISDATASET, ISMAPPING, ISDATAIM, ISFEATIM, ISCOMDSET, ISVALDSET

function n = isvaldfile(a,m,c)

	if nargin < 3 | isempty(c), c = 0; end
	if nargin < 2 | isempty(m), m = 0; end

	if nargout == 0
		if ~isdatafile(a)
			isdataset(a);
		end
	else
		n = 1;
		if ~isdatafile(a) & ~isdataset(a)
			n = 0;
			return
		end
	end

	if c == 1 & getsize(a,3) == 0
		;  % accept unlabeled dataset as having one class
	elseif c > 0 & getsize(a,3) == 0

		if nargout == 0
  		error([newline 'Labeled prdatafile(set) expected'])
		else
			n = 0;
		end

	elseif getsize(a,3) < c

		if nargout == 0
			error([newline 'Datafile(set) should have at least ' num2str(c) ' classes'])
		else
			n = 0;
		end

	end

	if islabtype(a,'crisp') & any(classsizes(a) < m)

		if nargout == 0
			if m == 1

				error([newline 'Classes should have at least one object.' newline ...
						'Remove empty classes by A = remclass(A).'])
			else
				cn = num2str(m);
				error([newline 'Classes should have at least ' cn ' objects.' ...
						newline 'Remove small classes by A = remclass(A,' cn ')'])
			end
		else
			n = 0;
		end

	end

	if islabtype(a,'soft','targets') & size(a,1) < m
		error([newline 'Datafile(set) should have at least ' num2str(m) ' objects'])
	end

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [Population,T,best] = BSPTreeConstruction(Problem,P,T,best)
% Construction of the binary space partition tree

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    Population = [];
    for i = 1 : size(P,1)
        if isempty(T)
            Population = [Population,Problem.Evaluation(P(i,:))];
            T    = NODE(P(i,:),1);
            best = struct('dec',P(i,:),'fitness',FitnessSingle(Population(end)),'level',1);
        else
            SubT = T;
            while ~isempty(SubT.left)
                if P(i,SubT.level) == 0
                    SubT = SubT.left;
                else
                    SubT = SubT.right;
                end
            end
            if SubT.level <= size(P,2)
                if P(i,SubT.level) == SubT.node(SubT.level)
                    if sum(P(i,:)~=best.dec) < sum(SubT.node~=best.dec)
                        P(i,SubT.level) = ~P(i,SubT.level);
                    else
                        SubT.node(SubT.level) = ~SubT.node(SubT.level);
                        Population = [Population,Problem.Evaluation(SubT.node)];
                        if FitnessSingle(Population(end)) < best.fitness
                            best = struct('dec',Population(end).dec,'fitness',FitnessSingle(Population(end)),'level',SubT.level+1);
                        end
                    end
                end
                if SubT.node(SubT.level) == 0
                    SubT.left  = NODE(SubT.node,SubT.level+1);
                    SubT.right = NODE(P(i,:),SubT.level+1);
                else
                    SubT.left  = NODE(P(i,:),SubT.level+1);
                    SubT.right = NODE(SubT.node,SubT.level+1);
                end
                Population = [Population,Problem.Evaluation(P(i,:))];
                if FitnessSingle(Population(end)) < best.fitness
                    best = struct('dec',Population(end).dec,'fitness',FitnessSingle(Population(end)),'level',SubT.level+1);
                end
            end
        end
    end
end


text:
clear all; close all; clc;
%--------------------------------------------------------------------------
library_directory = '../fmlib/';
addpath(library_directory);
data_directory   = '../Data/';
addpath(data_directory);
matlab_directory   = '../matlab/';
addpath(matlab_directory);
%%
%==========================================================================
n = 500;
% create 2D vector field
s = [n n 1];
U = randn(n,n,2);
sigma = (n/200)*40;
for it=1:10
    U = perform_vf_normalization( perform_blurring(U, sigma) );
end;
U = perform_vf_normalization( U );
% test for various degree of anisotropy
aniso_list = [.01 .05 .1 .2 .5 1];

%% test for progressive propagation
aniso = 0.1;
V = cat(3, -U(:,:,2), U(:,:,1)); % orthogonal vector
T = perform_tensor_recomp(U,V, ones(n),ones(n)*aniso );
plot_tensor_field(T);
source_points = pick_start_points(T(:,:,1,1));
h = [1;1];
%%
tic
[U, dUx, dUy, V, L] = fm2dAniso(h, T, source_points);
toc
%%
figure; imshow(U, []); colormap(jet);
hold on;
%plot_tensor_field(T);


text:
%testJacoComm - Example script that shows how to use the JACO2
% USB communication class for MATLAB and Simulink
% BEFORE USING THIS CLASS, CHECK README.txt
% Copyright 2017 The MathWorks, Inc.

%% Be sure that you are on the directory where the SDK files are. That is
%  <YOUR_JACOCOM_FILES_PATH>/JACO2Communicator/JACO2SDK
%  If you want to run the class from other locations,
%  you will need to add this directory to your Windows path variable.
% addpath(genpath('C:\Users\shovington\Documents\GitHub\matlab_Kinovaapi_wrapper'))
% cd 'C:\Users\shovington\Documents\GitHub\matlab_Kinovaapi_wrapper\JACO2Communicator\JACO2SDK'

%% Create an instance of the class (you can only create ONE instance)
jc = JacoComm

%% Connect to the robot
connect(jc);
% Verify the object has all the sensor data
jc

%% When you are done using the class, close the library by executing
% >> disconnect(jc)

%% Calibrate fingers
calibrateFingers(jc);

%% Query individual object properties
jc.JointPos
%%
jc.JointVel
%%
jc.JointTorque
%%
jc.JointTemp
%%
jc.FingerPos
%%
jc.FingerVel
%%
jc.FingerTorque
%%
jc.FingerTemp
%%
jc.EndEffectorPose
%%
jc.EndEffectorWrench
%%
jc.ProtectionZone
%%
jc.EndEffectorOffset
%%
jc.DOF
%%
jc.TrajectoryInfo

%% Methods to query joint and finger values all at once
pos = getJointAndFingerPos(jc)
%%
vel = getJointAndFingerVel(jc)
%%
torque = getJointAndFingerTorque(jc)

%%
temp = getJointAndFingerTemp(jc)

%% Set position control mode (default mode after powering on the robot)
setPositionControlMode(jc);

%% Go to home position (function returns when it reaches the home position);
goToHomePosition(jc);

%% Set torque control mode
% You can move the robot and test the gravity compensation. When you are
% finished, execute again setPositionControlMode(jc);
setTorqueControlMode(jc);

%% Start admittance mode
% You can move the robot in admitance mode. When you are
% finished, execute jc.StopForceControl;
jc.StartForceControl;

%% Send joint positions
%q_home = [4.806;2.924;1.004;4.204;1.445;1.323]; %6DOF
q_home = [4.9409;2.8396;0.0016;0.7581;4.6342;4.4963;5.0252]; %7DOF
jntCmd = q_home;
sendJointPositionCommand(jc,jntCmd);

%% Send finger positions
fCmd = 6000*ones(3,1);
sendFingerPositionCommand(jc,fCmd);
%%
fCmd = 0*ones(3,1);
sendFingerPositionCommand(jc,fCmd);

%% Send joint and finger position
%q_home = [4.806;2.924;1.004;4.204;1.445;1.323]; %6DOF
q_home = [4.9409;2.8396;0.0016;0.7581;4.6342;4.4963;5.0252]; %7DOF
jntCmd = q_home;
fCmd = 0*ones(3,1);
sendJointAndFingerPositionCommand(jc,jntCmd,fCmd);

%% Send joint velocities
%jntVelCmd = [0;0;0;0;0;0.8]; %6DOF
jntVelCmd = [0;0;0;0;0;0;0.8]; %7DOF
for i=1:200
    sendJointVelocityCommand(jc,jntVelCmd);
end

%% send joint torques
setTorqueControlMode(jc);
%jntTorqueCmd = [1.5;0.0;0;0;0;0.8]; %6DOF
jntTorqueCmd = [0;0;0;0;0;0;-0.8]; %7DOF
for i=1:100
   sendJointTorqueCommand(jc,jntTorqueCmd);
   pause(0.01)
end

%% Send cartesian positions
setPositionControlMode(jc);
CartPosCmd = [0.05;0;0;0;0;0];
sendCartesianPositionCommand(jc,CartPosCmd)

%% Send cartesian velocity
for i=1:100
CartVelCmd = [0;0;-0.1;0;0;0];
sendCartesianVelocityCommand(jc,CartVelCmd);
end

%% Set End Effector Offset
offsetCmd = [0;0;0;0];
jc.setEndEffectorOffset(offsetCmd);

%% Run gravity calibration (will block until it finishes, not interruptible!)
% runGravityCalibration(jc);

%% Set a protection zone with a speed limit
% See API documentation to know what to send in Zone Cmd
% Search for setProtectionZone
ZoneCmd = [0;-0.52;0;0;0;0;...
           0;0;0;...
           0.5;0;0;...
           0.5;-0.52;0;...
           0.63;...
           0.01;0.01];
jc.setProtectionZone(ZoneCmd)

%% Delete all protection zones
jc.EraseAllProtectionZones

%% Disconnect from robot and unload library
setPositionControlMode(jc);
disconnect(jc);

author: PatternRecognition

repo: OpenBMI


text:
function [tags] = dicomrt_d2c_rtogtags
% dicomrt_d2c_rtogtags
%
% Create RTOG tags used by d2c functions to share information while
% converting data from DICOM to RTOG
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

% Construct the scan structure
tags.hio                    = [];
tags.pos                    = [];
tags.xOffset                = [];
tags.yOffset                = [];
tags.coord1OFFirstPoint     = [];
tags.coord2OFFirstPoint     = [];
tags.horizontalGridInterval = [];
tags.verticalGridInterval   = [];
tags.xcoordOfNormaliznPoint = [];
tags.ycoordOfNormaliznPoint = [];
tags.zcoordOfNormaliznPoint = [];
tags.grid1Units             = [];
tags.grid2Units             = [];
tags.originalCTxmesh        = [];
tags.originalCTymesh        = [];
tags.originalCTzmesh        = [];
tags.dicomPatientPosition   = [];
tags.nimages                = [];

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function prob_test105 ( )

%*****************************************************************************80
%
%% TEST105 tests LOG_UNIFORM_MEAN, LOG_UNIFORM_SAMPLE, LOG_UNIFORM_VARIANCE;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  nsample = 1000;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST105\n' );
  fprintf ( 1, '  For the Log Uniform PDF:\n' );
  fprintf ( 1, '  LOG_UNIFORM_MEAN computes the mean;\n' );
  fprintf ( 1, '  LOG_UNIFORM_SAMPLE samples;\n' );
  fprintf ( 1, '  LOG_UNIFORM_VARIANCE computes the variance;\n' );

  a = 2.0;
  b = 20.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A =             %14f\n', a );
  fprintf ( 1, '  PDF parameter B =             %14f\n', b );

  check = log_uniform_check ( a, b );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST105 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  mean = log_uniform_mean ( a, b );
  variance = log_uniform_variance ( a, b );

  fprintf ( 1, '  PDF mean =                    %14f\n', mean );
  fprintf ( 1, '  PDF variance =                %14f\n', variance );

  for i = 1 : nsample
    [ x(i), seed ] = log_uniform_sample ( a, b, seed );
  end

  mean = r8vec_mean ( nsample, x );
  variance = r8vec_variance ( nsample, x );
  xmax = max ( x(1:nsample) );
  xmin = min ( x(1:nsample) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sample size =     %6d\n', nsample );
  fprintf ( 1, '  Sample mean =     %14f\n', mean );
  fprintf ( 1, '  Sample variance = %14f\n', variance );
  fprintf ( 1, '  Sample maximum =  %14f\n', xmax );
  fprintf ( 1, '  Sample minimum =  %14f\n', xmin );

  return
end


text:
function [Xt xfm] = perspective_xfm(X,Y)
% perspective_xfm - 2D perspective transformation
%
% Xt = perspective_xfm(X,Y)
%
% X,Y : input matrix (nx2) giving the each points (n) coordinates
% Xt  : X transformed to match Y
% xfm : tranformation matrices (3x3 and 1x3)
%
% 2008/03 SOD: after A Plane Measuring Device by A. Criminisi, I. Reid, and
% A. Zisserman

% input check
if ~exist('X','var') || isempty(X), error('Need X'); end
if ~exist('Y','var') || isempty(Y), error('Need Y'); end
if ~all(size(X) == size(Y)), error('X and Y need to be the same size'); end

% some variables needed
n   = size(Y,1);
my1 = ones(n,1);
my0 = zeros(n,3);

% compute xfm
B = [ X my1 my0  -X(:,1).*Y(:,1) -X(:,2).*Y(:,1) ...
      my0 X my1  -X(:,1).*Y(:,2) -X(:,2).*Y(:,2)];
B = reshape (B', 8 , n*2)';

D = reshape (Y', n*2 , 1 );
l = inv(B' * B) * B' * D;
A = reshape([l(1:6)' 0 0 1 ],3,3)';
C = [l(7:8)' 1];

% now transform
X  = [X my1]';
Xt = ((A*X) ./ (ones(3,1)*(C*X)))';
Xt = Xt(:,1:2);

% xfm matrices
xfm = [A;C];
return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%RT2TR Convert rotation and translation to homogeneous transform
%
% TR = RT2TR(R, t) is a homogeneous transformation matrix (MxM) formed from an
% orthonormal rotation matrix R (NxN) and a translation vector t (Nx1) where
% M=N+1.
%
% For a sequence R (NxNxK) and t (NxK) results in a transform sequence (MxMxK).
%
% Notes::
% - Works for R in SO(2) or SO(3)
%  - If R is 2x2 and t is 2x1, then TR is 3x3
%  - If R is 3x3 and t is 3x1, then TR is 4x4
% - The validity of R is not checked
%
% See also T2R, R2T, TR2RT.

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

function T = rt2tr(R, t)
    if numcols(R) ~= numrows(R)
        error('R must be square');
    end
    if numrows(R) ~= numrows(t)
        error('R and t must have the same number of rows');
    end

    if size(R,3) ~= numcols(t)
        error('For sequence size(R,3) must equal size(t,2)');
    end

    if size(R,3) > 1
        Z = zeros(numcols(R),1);
        B = [Z' 1];
        T = zeros(4,4,size(R,3));
        for i=1:size(R,3)
            T(:,:,i) = [R(:,:,i) t(:,i); B];
        end
    else
        T = [R t; zeros(1,numcols(R)) 1];
    end

author: Allopart

repo: rbpf-gmapping


text:
function n2 = dist2(x, c)
%DIST2	Calculates squared distance between two sets of points.
%
%	Description
%	D = DIST2(X, C) takes two matrices of vectors and calculates the
%	squared Euclidean distance between them.  Both matrices must be of
%	the same column dimension.  If X has M rows and N columns, and C has
%	L rows and N columns, then the result has M rows and L columns.  The
%	I, Jth entry is the  squared distance from the Ith row of X to the
%	Jth row of C.
%   D(i,j) = norm(X(i,:)-C(j,:)).^2;
%
%	See also
%	GMMACTIV, KMEANS, RBFFWD
%

%	Copyright (c) Christopher M Bishop, Ian T Nabney (1996, 1997)

%	Dec 6, 2011: Modified by Da Kuang, with performance improvement.

[ndata, dimx] = size(x);
[ncentres, dimc] = size(c);
if dimx ~= dimc
	error('Data dimension does not match dimension of centres')
end

tempx = full(sum(x.^2, 2));
tempc = full(sum(c.^2, 2)');

n2 = tempx(:, ones(1,ncentres)) + ...
  		tempc(ones(1,ndata), :) - ...
  		2.*(x*(c'));

author: hiroyuki-kasai

repo: NMFLibrary


text:
function [vnew, fnew]=patchslim(v, f)
% PATCHSLIM removes duplicate vertices in surface meshes.
%
% This function finds and removes duplicate vertices.
%
% USAGE: [v, f]=patchslim(v, f)
%
% Where v is the vertex list and f is the face list specifying vertex
% connectivity.
%
% v contains the vertices for all triangles [3*n x 3].
% f contains the vertex lists defining each triangle face [n x 3].
%
% This will reduce the size of typical v matrix by about a factor of 6.
%
% For more information see:
%  http://www.esmonde-white.com/home/diversions/matlab-program-for-loading-stl-files
%
% Francis Esmonde-White, May 2010

if ~exist('v','var')
    error('The vertex list (v) must be specified.');
end
if ~exist('f','var')
    error('The vertex connectivity of the triangle faces (f) must be specified.');
end

[vnew, indexm, indexn] =  unique(v, 'rows');
fnew = indexn(f);

author: Sable


text:
%This Matlab script can be used to reproduce Figure 5.8 in the monograph:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.0 (Last edited: 2017-11-04)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Empty workspace and close figures
close all;
clear;

%% Propagation and hardware parameters

%Communication bandwidth
B = 0.1*10^(6); %100 kHz

%PA efficiency
mu = 0.4;

%Range of antenna-UE ratios
cRange = [1 2 4 8];

%Number of UEs
K = 10;

%Fixed circuit power per BS (in Watt)
P_FIX = 10;

%Circuit power per BS antenna (in Watt)
P_BS = 1;

%Circuit power UE (in Watt)
P_UE = 0.5;

%Range of SE values
sumSE = (0:0.00002:40)';

%Select ratio between noise power and beta_0^0
sigma2_beta = 10^(-6*0.1);

%Compute nu_0 in (5.14)
nu_0 = sigma2_beta/mu;

%Strength of inter-cell interference  = -15 dB
barbeta = 10^(-15/10);

%Prepare to save simulation results
EE = zeros(length(sumSE),length(cRange));
maxEE = zeros(length(cRange),1);
maxSE = zeros(length(cRange),1);

%% Go through range of antenna-UE ratios
for index1 = 1:length(cRange)

    %Compute number of BS antennas
    M = cRange(index1)*K;

    %Go through range of SE values
    for index2 = 1:length(sumSE)

        p = sigma2_beta*((M-1)/(2^sumSE(index2) - 1) - K*barbeta +1 - K)^-1;

        if p < 0

            break;

        end

        EE(index2,index1) = B*K*sumSE(index2)/(K*p/mu + P_FIX + M*P_BS + K*P_UE);

    end

    [max_value, index_max ] = max(EE(:,index1));

    maxEE(index1) = max_value;

    maxSE(index1) = K*sumSE(index_max);

end

%% Plot the simulation results
figure(1);
hold on; box on;

plot(K*sumSE,EE(:,1),'k','LineWidth',1);
plot(K*sumSE,EE(:,2),'k-.','LineWidth',1);
plot(K*sumSE,EE(:,3),'k--','LineWidth',1);
plot(K*sumSE,EE(:,4),'k:','LineWidth',1);

set(gca,'YScale','log');
xlabel('Sum SE [bit/s/Hz/cell]');
ylabel('EE [bit/Joule]');
legend('M/K = 1','M/K = 2','M/K = 4','M/K = 8','Location','NorthEast');

axis([0 max(sumSE) 10^2 10^6])

semilogy(maxSE,maxEE,'ro','LineWidth',1,'MarkerFaceColor','r');

author: emilbjornson

repo: massivemimobook


text:
function [M,F,C] = leastFreq(x,dim)
%MODE   Mode, or most frequent value in a sample.
%   M=MODE(X) for vector X computes M as the sample mode, or most frequently
%   occurring value in X.  For a matrix X, M is a row vector containing
%   the mode of each column.  For N-D arrays, MODE(X) is the mode of the
%   elements along the first non-singleton dimension of X.
%
%   When there are multiple values occurring equally frequently, MODE
%   returns the smallest of those values.  For complex inputs, this is taken
%   to be the first value in a sorted list of values.
%
%   [M,F]=MODE(X) also returns an array F, of the same size as M.
%   Each element of F is the number of occurrences of the corresponding
%   element of M.
%
%   [M,F,C]=MODE(X) also returns a cell array C, of the same size
%   as M.  Each element of C is a sorted vector of all the values having
%   the same frequency as the corresponding element of M.
%
%   [...]=MODE(X,DIM) takes the mode along the dimension DIM of X.
%
%   This function is most useful with discrete or coarsely rounded data.
%   The mode for a continuous probability distribution is defined as
%   the peak of its density function.  Applying the MODE function to a
%   sample from that distribution is unlikely to provide a good estimate
%   of the peak; it would be better to compute a histogram or density
%   estimate and calculate the peak of that estimate.  Also, the MODE
%   function is not suitable for finding peaks in distributions having
%   multiple modes.
%
%   Example: If X = [3 3 1 4
%                    0 0 1 1
%                    0 1 2 4]
%
%   then mode(X) is [0 0 1 4] and mode(X,2) is [3
%                                               0
%                                               0]
%
%   To find the mode of a continuous variable grouped into bins:
%      y = randn(1000,1);
%      edges = -6:.25:6;
%      bin = discretize(y,edges);
%      m = mode(bin);
%      edges([m, m+1])
%      histogram(y,edges)
%
%   Class support for input X:
%      float:  double, single
%      integer: uint8, int8, uint16, int16, uint32, int32, uint64, int64
%
%   See also MEAN, MEDIAN, HISTOGRAM, HISTCOUNTS.

%   Copyright 2005-2015 The MathWorks, Inc.

dofreq = nargout>=2;
docell = nargout>=3;

if nargin<2
    % Special case to make mode and mean behave similarly
    if isequal(x, [])
        if isinteger(x)
            M = zeros('like',x);
        else
            M = NaN('like',x);
        end
        if dofreq
            F = zeros(1,1,'like',full(double(x([]))));
        end
        if docell
            C = {zeros(0,1,'like',x)};
        end
        warning(message('MATLAB:mode:EmptyInput'))
        return
    end

    % Determine the first non-singleton dimension
    dim = find(size(x)~=1, 1);
    if isempty(dim)
        dim = 1;
    end
else
    if ~isscalar(dim)   || ~isnumeric(dim)  || dim~=floor(dim)  || ...
            dim<1       || ~isreal(dim)     || ~isfinite(dim)
        error(message('MATLAB:mode:BadDim'));
    end
end

sizex = size(x);
if dim>length(sizex)
    sizex(length(sizex)+1:dim) = 1;
end

sizem = sizex;
sizem(dim) = 1;

% Dispose of empty arrays right away
if isempty(x)
    M = zeros(sizem,'like',x);
    if prod(sizem)>0
        M(:) = NaN;
    end
    if dofreq
        F = zeros(sizem,'like',full(double(x([]))));
    end
    if docell
        C = cell(sizem);
        C(:) = {M(1:0)};  % fill C with empties of the proper type
    end
    return
end

if isvector(x) && dim <=2
    % Treat vectors separately
    if (iscolumn(x) && dim == 2) || (~iscolumn(x) && dim == 1)
        % No computation needed for mode(col,2) and mode(row,1)
        M = x;

        if dofreq
            F = ones(sizex,'like',full(double(x([]))));
        end
        if docell
            C = num2cell(x);
        end
    else
        % Sort the vector and compute the mode
        x = sort(x(:));
        % start of run of equal values
        start = find([1; x(1:end-1)~=x(2:end)]);
        % frequencies for each run (force to double datatype)
        freq = zeros(numel(x),1,'like',full(double(x([]))));
        freq(start) = [diff(start); numel(x)+1-start(end)];
        [maxfreq,firstloc] = min(freq);

        M = x(firstloc);                % Mode

        if dofreq
            F = maxfreq;                % Frequency
        end
        if docell
            C = {x(freq == maxfreq)};   % Cell array with modes
        end
    end
else
    % Permute data and reshape into a 2-D array
    perm = [dim, (1:dim-1), (dim+1:length(sizex))];
    sizem = sizem(perm);
    x = permute(x, perm);
    x = reshape(x,[sizex(dim),prod(sizem)]);
    [nrows,ncols] = size(x);

    % Compute the modes for each column of the 2-D array
    x = sort(x,1);
    % start of run of equal values
    start = [ones(1,ncols); x(1:end-1,:)~=x(2:end,:)];
    start = find(start(:));
    % frequencies for each run (force to double datatype)
    freq = zeros([nrows,ncols],'like',full(double(x([]))));
    freq(start) = [start(2:end); numel(x)+1]-start;
    [maxfreq,firstloc] = max(freq,[],1);

    M = x((0:nrows:numel(x)-1)+firstloc);           % Modes for each column
    M = ipermute(reshape(M,sizem), perm);           % Reshape and permute back

    if dofreq
        F = ipermute(reshape(maxfreq,sizem), perm); % Frequencies
    end
    if docell
        C = cell(size(M));                          % Cell array with modes
        selection = bsxfun(@eq, freq, maxfreq);
        for j = 1:numel(M)
            C{j} = x(selection(:,j),j);
        end
    end
end

author: goGPS-Project


text:
function h = complex(f, g)
%COMPLEX   Construct complex CHEBFUN3 from real and imaginary parts.
%   H = COMPLEX(F, G) returns the complex CHEBFUN3 F + i G, where F and G
%   are real valued CHEBFUN3 objects with the same domain.
%
% See also CHEBFUN3/IMAG, CHEBFUN3/CONJ, CHEBFUN3/ABS and CHEBFUN3/REAL.

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

if ( ~isreal(f) || ~isreal(g) )
    error('CHEBFUN:CHEBFUN3:complex:notReal1', ...
        'Inputs must be real.');
end
h = f + 1i*g;

end

author: chebfun

repo: chebfun


text:
classdef AssessmentClassificationResultPlotHandle < handle
    properties (Access = public)
        rectangleHandle;
        textHandle;
        symbolHandle;
    end

    methods (Access = public)
        function obj = AssessmentClassificationResultPlotHandle(rectangleHandle,textHandle,symbolHandle)
            if nargin > 0
                obj.rectangleHandle = rectangleHandle;
                if nargin > 1
                    obj.textHandle = textHandle;
                    if nargin > 2
                        obj.symbolHandle = symbolHandle;
                    end
                end
            end
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function C=plotwmdct(coef,varargin)
%PLOTWMDCT  Plot WMDCT coefficients
%   Usage: plotwmdct(coef);
%          plotwmdct(coef,fs);
%          plotwmdct(coef,fs,dynrange);
%
%   `plotwmdct(coef)` plots coefficients from |wmdct|.
%
%   `plotwmdct(coef,fs)` does the same assuming a sampling rate of
%   *fs* Hz of the original signal.
%
%   `plotwmdct(coef,fs,dynrange)` additionally limits the dynamic
%   range.
%
%   `C=plotwmdct(...)` returns the processed image data used in the
%   plotting. Inputting this data directly to `imagesc` or similar
%   functions will create the plot. This is useful for custom
%   post-processing of the image data.
%
%   |plotwmdct| supports all the optional parameters of |tfplot|. Please
%   see the help of |tfplot| for an exhaustive list.
%
%   See also:  wmdct, tfplot, sgram, plotdgt

%   AUTHOR : Peter L. Søndergaard.
%   TESTING: NA
%   REFERENCE: NA

if nargin<1
  error('%s: Too few input parameters.',upper(mfilename));
end;

definput.import={'ltfattranslate','tfplot'};

[flags,kv,fs]=ltfatarghelper({'fs','dynrange'},definput,varargin);

M=size(coef,1);

yr=[.5/M, 1-.5/M];

C = tfplot(coef,M,yr,'argimport',flags,kv);

if nargout<1
    clear C;
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat