text:
clear all;
figure

xvals=[];
startval = 1;
endval = 4;
for r=startval:.00025:endval
    r
    x = 0.5;
    for i=1:10000
        x = logistic(x,r);
    end
    xss = x;  % steady state x
    for i=1:1000
        x = logistic(x,r);
        xvals(1,length(xvals)+1) = r;
        xvals(2,length(xvals)) = x;
        if(abs(x-xss)<.001)
            break
        end
    end

end
figure(1)
plot(xvals(2,:),xvals(1,:),'.','LineWidth',.1,'MarkerSize',1.2,'Color',[0 0 0])
axis([0 1 1 endval])
axis ij
figure(2)
plot(xvals(2,:),xvals(1,:),'.','LineWidth',.1,'MarkerSize',1.2,'Color',[0 0 0])
axis([0 1 3.45 4])
axis ij

author: dynamicslab


text:
function F = diff(F, k, dim)
%DIFF   Derivative of a SEPARABLEAPPROX object.
%   DIFF(F) is the derivative of F along the y direction.
%
%   DIFF(F, N) is the Nth derivative of F in the y direction.
%
%   DIFF(F, N, DIM) is the Nth derivative of F along the dimension DIM.
%     DIM = 1 (default) is the derivative in the y-direction.
%     DIM = 2 is the derivative in the x-direction.
%
%   DIFF(F, [NX NY]) is the partial derivative of NX of F in the first variable,
%   and NY of F in the second derivative. For example, DIFF(F,[1 2]) is
%   d^3F/dxd^2y.
%
% See also GRADIENT, SUM, PROD.

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

% Empty check:
if ( isempty(F) )
    return
end

% Default to first derivative:
if ( nargin < 2 )
    k = 1;
end

% Default to partial derivative in y:
if ( nargin < 3 )
    dim = 1;
elseif ( numel(dim) ~= 1 )
    error('CHEBFUN:SEPARABLEAPPROX:diff:dim1', 'Dim should be either 1 or 2.');
end

% Diff the individual column and row slices:
if ( numel(k) == 2 && nargin < 3 )
   F.cols =  diff(F.cols, k(2));
   F.rows = diff(F.rows, k(1));
elseif ( dim == 1 )
    F.cols = diff(F.cols, k);
elseif ( dim == 2 )
    F.rows = diff(F.rows, k);
else
    error('CHEBFUN:SEPARABLEAPPROX:diff:dim2', ...
        'Can compute derivative in x or y only.');
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function h=drawEllipse(z,A,gammaVal,varargin)
%DRAWELLIPSE Draw 2D ellipses or 3D ellipsoids on the current plot, where a
%            point zp on the ellipse/ ellipsoid satisfies the equation
%            (zp-z)'*A*(zp-z)=gammaVal. If z is a 2D/3D location estimate
%            and inv(A) a Gaussian covariance matrix associated with the
%            estimate, then the ellipse is a probability region, where
%            gammaVal determines what amount of probability is in the
%            region. If omitted, gammaVal=16.2235 for 2D or
%            gammaVal=18.8049 for 3D, which corresponds approximately to
%            the 99.97% confidence region, is used (The value comes from
%            respectively ChiSquareD.invCDF(0.9997,2) and
%            ChiSquareD.invCDF(0.9997,3)).
%
%INPUTS: z A 2XN or 3XN vector corresponding to the centers of the N
%          ellipses or ellipsoids that should be drawn.
%        A A 2X2XN or 3X3XN set of N positive definite matrices that
%          specify the size and shape of the ellipse or ellipsoids, where
%          a point zp is on the ith ellipse/ ellipsoid if
%          (zp-z(:,i))'*A(:,:,i)*(zp-z(:,i))=gammaVal.
%          In three dimensions, A can have one zero eigenvalue, in which
%          case the fill3 command is used to draw an ellipse in 3D.
% gammaVal An optional parameter specifying the size of the ellipse/
%          ellipsoid. If omitted or an empty matrix is passed, then
%          gammaVal=16.2235 in 2D or gammaVal=18.8049 in 3D is used. These
%          are approximately the values for a 99.97% confidence region if A
%          are inverse covariance matrices of a Gaussian distribution.
%          gammaVal must be positive.
% varargin Sets of values that should be passed to the plot function to
%          format the ellipses or that will be passed to the surf function
%          to format the ellipsoids. For example, one could call the
%          function as drawEllipse(z,A,gammaVal,'--r','linewidth',2) to
%          plotellipses as thick, red lines. Note that Matlab will not
%          always properly render dashed lines due to the number of points
%          used to plot the shape. Also, if the ellipsoid is in 3D, but is
%          actually just a 2D ellipse, then if this parameter is omitted,
%          the fill3 command will be passed the option to color the
%          ellipse black. Otherwise, if this parameter is given, the user
%          must specify the ellipse color as it is not done here.
%
%OUTPUTS: h A NX1 cell array containing the plot objects for each of the
%           ellipses. This can be useful if, for example, one wishes to
%           change the transparency of object i to 50%, one can use the
%           alpha(h{i},0.5) command.
%
%When drawing 2D ellipses, the function getEllipsePoints is used to get the
%points on the ellipse, which are subsequently plotted.
%
%To find the points on an axis-aligned ellipsoid, we will use the ellipsoid
%function that is built into Matlab. Suppose that l=[x;y;z] and
%D=diag(a,b,c). Then the ellipsoid equation is x^2*a+y^2*b+z^2c=gammaVal.
%This means that the semi-axis lengths of this centered, axis-aligned
%ellipsoid are sqrt(gammaVal/a), sqrt(gammaVal/b), and sqrt(gammaVal/c).
%2501 vertices are used in the ellipsoids.
%
%July 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numEllipse=size(z,2);

%Save the value of hold on the plot so that it can be reset to its previous
%state after plotting all of the ellipses.
holdVal=ishold();

numDim=size(z,1);

if(nargout>0)
    h=cell(numEllipse,1);
end

switch(numDim)
    case 2
        if(nargin<3||isempty(gammaVal))
            gammaVal=16.2235;
        end

        numPoints=1000;
        for curEllip=1:numEllipse
            zp=getEllipsePoints(z(:,curEllip),A(:,:,curEllip),gammaVal,numPoints,false);

            %Make sure that all of the ellipses are printed.
            if(curEllip~=1)
                hold on
            end
            x=plot(zp(1,:),zp(2,:),varargin{:});

            if(nargout>0)
               h{curEllip}=x;
            end
        end
    case 3
        if(nargin<3||isempty(gammaVal))
            gammaVal=18.8049;
        end

        for curEllip=1:numEllipse
            if(rank(A(:,:, curEllip))<3)
                numPoints=500;
                %If the ellipsoid is actually just an ellipse, then we will
                %draw the 2D ellipse in 3D.

                %Perform an eigenvalue decomposition of A.
                [V,D]=eig(A(:,:,curEllip));
                %We will put the zero eigenvalue at the start, if it is
                %not already there.
                [~,idx]=sort(diag(D),'ascend');
                V=V(idx,idx);
                D=D(idx,idx);

                a=D(2,2);
                b=D(3,3);

                xBound=sqrt(gammaVal/a);
                x=linspace(-xBound,xBound,numPoints);
                %The real command deals with possible finite precision
                %issues.
                y=real(sqrt((gammaVal-x.^2*a)/b));
                %The centered, axis-aligned ellipse points.
                l=[zeros(1,2*numPoints);
                   x,fliplr(x);
                   y,-fliplr(y)];
                zp=bsxfun(@plus,V*l,z(idx,curEllip));
                idxInv=inversePermutation(idx);
                %Undo the reordering.
                zp=zp(idxInv);

                %Make sure that all of the ellipses are printed.
                if(curEllip~=1)
                    hold on
                end
                %Draw the ellipse as a filled shape.
                if(isempty(varargin))
                    x=fill3(zp(1,:),zp(2,:),zp(3,:),'k');
                else
                    x=fill3(zp(1,:),zp(2,:),zp(3,:),varargin{:});
                end
                if(nargout>0)
                    h{curEllip}=x;
                end
            else
                numPoints=50;
                %Perform an eigenvalue decomposition of A.
                [V,D]=eig(A(:,:,curEllip));
                zCur=z(:,curEllip);

                [xp,yp,zp]=ellipsoid(0,0,0,sqrt(gammaVal/D(1,1)),sqrt(gammaVal/D(2,2)),sqrt(gammaVal/D(3,3)),numPoints);

                %Rotate to the correct orientation
                xyzPoints=V*[xp(:)';yp(:)';zp(:)'];
                xp=reshape(zCur(1)+xyzPoints(1,:),numPoints+1,numPoints+1);
                yp=reshape(zCur(2)+xyzPoints(2,:),numPoints+1,numPoints+1);
                zp=reshape(zCur(3)+xyzPoints(3,:),numPoints+1,numPoints+1);

                %Make sure that all of the ellipses are printed.
                if(curEllip~=1)
                    hold on
                end

                x=surf(xp,yp,zp,varargin{:});
                if(nargout>0)
                    h{curEllip}=x;
                end
            end
        end
    otherwise
        error('Invalid Dimensionality')
end

%Restore the hold value to its original setting.
if(~holdVal)
    hold off
end

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function figureJ(varargin);
if nargin==0
    gcf=figure;
else
gcf=figure(varargin{1});
end;
set(gcf,'Color',[1 1 1]);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function predicted = badacostApply(X, classfr)
% Apply learned badacost weak learners ensemble
%
% USAGE
%  hs = badacostApply( X, model, [maxDepth], [minWeight], [nThreads] )
%
% INPUTS
%  X          - [FXN] N num vectors to classfy, F num feature vectors
%  model      - learned boosted tree classifier
%
% OUTPUTS
%  predicted  - [Nx1] predicted output labels
%
% EXAMPLE
%
% See also badacostTrain
%
% Author: Antonio Baldera, modified by Jose M. Buenaposada

n          = size(X,2);
margin_vec = zeros(classfr.num_classes, n);

for i=1:length(classfr.WEAK_LEARNERS)
  % z is a row vector with the labels
  z  = classfr.classify_weak_learner(classfr.WEAK_LEARNERS{i}, X);

  for j=1:n
    margin_vec(:,j) = margin_vec(:,j) + (classfr.WEIGHTS(i).*classfr.Y(:, z(j)));
  end
end;

% WARNING: Change to accomodate with theory (2016/11)
%[~, predicted] = min(classfr.Cprime' * margin_vec);
[~, predicted] = min(classfr.Cprime * margin_vec);
predicted = predicted(:);

end

author: MengyangPu

repo: EDTER

save_path: github-repos/MATLAB/MengyangPu-EDTER


text:
function moulton_set_test ( )

%*****************************************************************************80
%
%% MOULTON_SET_TEST tests MOULTON_SET.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 April 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MOULTON_SET_TEST\n' );
  fprintf ( 1, '  MOULTON_SET sets up Adams-Moulton quadrature;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Index       X             W\n' );
  fprintf ( 1, '\n' );

  for n = 1 : 10

    [ x, w ] = moulton_set ( n );

    fprintf ( 1, '\n' );

    for i = 1 : n
      fprintf ( 1, '  %2d  %12g  %12g\n', i, x(i), w(i) );
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Fig. 8.12   Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
%

clear all;
close all;

numGs=1;
denGs=[1 0 0];  % s^2

numDs=.81*[1 .2];
denDs=[1 2];

numC=conv(numGs,numDs);
denC=conv(denGs,denDs);

[numCL,denCL]=feedback(numC,denC,1,1);

tf=30;
t=0:.2:tf;
y=step(numCL,denCL,t);

axis([0 30 0 1.5])
plot(t,y),grid
hold on

% discrete design

T=1;
[numGz,denGz]=c2dm(numGs,denGs,T,'zoh');

numDz2=.389*[1 -.82];
denDz2=[1 -.135];

numz2=conv(numGz,numDz2);
denz2=conv(denGz,denDz2);

[numCLz2,denCLz2]=feedback(numz2,denz2,1,1);
N=tf/T+1;
td=0:1:tf;
yd2=dstep(numCLz2,denCLz2,N);

T2=0.5;
[numGz,denGz]=c2dm(numGs,denGs,T2,'zoh');

a=.2;
b=2;
Kc=0.81;
alpha=exp(-a*T2);
beta=exp(-b*T2);
Kd=Kc*a/b*(1-beta)/(1-alpha);
numDz3=Kd*[1 -alpha];
denDz3=[1 -beta];

numz3=conv(numGz,numDz3);
denz3=conv(denGz,denDz3);

[numCLz3,denCLz3]=feedback(numz3,denz3,1,1);
N3=61;
td2=0:T2:tf;
yd3=dstep(numCLz3,denCLz3,N3);
plot(td,yd2,'-',td,yd2,'*',td2,yd3,'o',td2,yd3,'-')
title('Fig. 8.12, Step responses of ghe continuous and digital implentations')
text(10,.6,'----------- Continuous design')
text(10,.4,'*----*----* Discrete equivalent design, T = 1 sec')
text(10,.2,'o----o----o Discrete equivalent design, T = 0.5 sec')
xlabel('Time (sec)')
ylabel('Plant output')
hold off

author: Sable

repo: mcbench-benchmarks


text:
%==========================================================================
% This is the testing code of SRMDNF for the widely-used <bicubic degradation>.
% For bicubic degradation, the basic setting is:
%   1. the blur kernel is delta kernel; "kernel = fspecial('gaussian',15,0.2)".
%   2. the downsampler is fixed to bicubic downsampler. You can either
%      train a new model with direct downsamper or use the estimated
%      kernel  k_b under direct downsampler.
%   3. there are three SRMDNF models, "SRMDNFx2.mat" for scale factor 2, "SRMDNFx3.mat"
%      for scale factor 3, and "SRMDNFx4.mat" for scale factor 4.
%==========================================================================
% The basic idea of SRMD is to learn a CNN to infer the MAP of general SISR, i.e.,
% solve x^ = arg min_x 1/2 ||(kx)\downarrow_s - y||^2 + lamda \Phi(x)
% via x^ = CNN(y,k;\Theta) or HR^ = CNN(LR,kernel;\Theta).
%
% There involves two important factors, i.e., blur kernel (k; kernel), in SRMDNF.
%
% For more information, please refer to the following paper.
%    @article{zhang2017learningsrmd,
%    title={Learning a Single Convolutional Super-Resolution Network for Multiple Degradations},
%    author={Kai, Zhang and Wangmeng, Zuo and Lei, Zhang},
%    year={2017},
%    }
%
% If you have any question, please feel free to contact with <Kai Zhang (cskaizhang@gmail.com)>.
%
% This code is for research purpose only.
%
% by Kai Zhang (Nov, 2017)
%==========================================================================
% clear; clc;
format compact;

addpath('utilities');
imageSets    = {'Set5','Set14','BSD100','Urban100'}; % testing dataset

%% select testing dataset, use GPU or not, ...
setTest      = imageSets([1,2,3,4]); % select the datasets for each tasks
showResult   = 0; % save restored images
pauseTime    = 0;
useGPU       = 1; % 1 or 0, true or false
method       = 'SRMDNF';
folderTest   = 'testsets';
folderResult = 'results';
if ~exist(folderResult,'file')
    mkdir(folderResult);
end

%% scale factor (2, 3, 4)
sf          = 4;

%% load model with scale factor sf
folderModel = 'models';
load(fullfile(folderModel,['SRMDNFx',int2str(sf),'.mat']));
%net.layers = net.layers(1:end-1);
net = vl_simplenn_tidy(net);
if useGPU
    net = vl_simplenn_move(net, 'gpu') ;
end

%% degradation parameter (kernel & noise level) setting
global degpar;
kernel = fspecial('gaussian',15, 0.2); % kernel is delta kernel. Note: the kernel size is fixed to 15X15.
degpar = single(net.meta.P*kernel(:));
tag    = ['_',method,'_x',num2str(sf),'_bicubic'];

for n_set = 1 : numel(setTest)

    %% search images
    setTestCur = cell2mat(setTest(n_set));
    disp('--------------------------------------------');
    disp(['    ----',setTestCur,'-----Super-Resolution-----']);
    disp('--------------------------------------------');
    folderTestCur = fullfile(folderTest,setTestCur);
    ext                 =  {'*.jpg','*.png','*.bmp'};
    filepaths           =  [];
    for i = 1 : length(ext)
        filepaths = cat(1,filepaths,dir(fullfile(folderTestCur, ext{i})));
    end

    %% prepare results
    eval(['PSNR_',setTestCur,'_x',num2str(sf),' = zeros(length(filepaths),1);']);
    eval(['SSIM_',setTestCur,'_x',num2str(sf),' = zeros(length(filepaths),1);']);
    folderResultCur = fullfile(folderResult, [setTestCur,tag]);
    if ~exist(folderResultCur,'file')
        mkdir(folderResultCur);
    end

    %% perform SISR
    for i = 1 : length(filepaths)

        HR  = imread(fullfile(folderTestCur,filepaths(i).name));
        C   = size(HR,3);
        if C == 1
            HR = cat(3,HR,HR,HR);
        end
        [~,imageName,ext] = fileparts(filepaths(i).name);
        HR  = modcrop(HR, sf);
        label_RGB = HR;
        %% bicubic degradation
        %blury_HR = imfilter(im2double(HR),double(kernel),'replicate'); % blur
        LR       = imresize(im2double(HR),1/sf,'bicubic'); % bicubic downsampling
        input    = im2single(LR);
        %input    = im2single(im2uint8(LR)); % another widely-used setting
        %tic
        if useGPU
            input = gpuArray(input);
        end
        res = vl_srmd(net, input,[],[],'conserveMemory',true,'mode','test','cudnn',true);
        %res = vl_srmd_concise(net, input); % a concise version of "vl_srmd".
        %res = vl_srmd_matlab(net, input); % you should also set "useGPU = 0;" and comment "net = vl_simplenn_tidy(net);"

        output_RGB = gather(res(end).x);

        %toc;
        if C == 1
            label  = mean(im2double(HR),3);
            output = mean(output_RGB,3);
        else
            label  = rgb2ycbcr(im2double(HR));
            output = rgb2ycbcr(double(output_RGB));
            label  = label(:,:,1);
            output = output(:,:,1);
        end

       %% calculate PSNR and SSIM
        [PSNR_Cur,SSIM_Cur] = Cal_PSNRSSIM(label*255,output*255,sf,sf); %%% single
        disp([setTestCur,'    ',int2str(i),'    ',num2str(PSNR_Cur,'%2.2f'),'dB','    ',filepaths(i).name]);
        eval(['PSNR_',setTestCur,'_x',num2str(sf),'(',num2str(i),') = PSNR_Cur;']);
        eval(['SSIM_',setTestCur,'_x',num2str(sf),'(',num2str(i),') = SSIM_Cur;']);
        if showResult
            imshow(cat(2,label_RGB,imresize(im2uint8(LR),sf),im2uint8(output_RGB)));
            drawnow;
            title(['SISR     ',filepaths(i).name,'    ',num2str(PSNR_Cur,'%2.2f'),'dB'],'FontSize',12)
            pause(pauseTime)
            imwrite(output_RGB,fullfile(folderResultCur,[imageName,'_x',int2str(sf),'_',int2str(PSNR_Cur*100),'.png']));% save results
        end

    end
    disp(['Average PSNR is ',num2str(mean(eval(['PSNR_',setTestCur,'_x',num2str(sf)])),'%2.2f'),'dB']);
    disp(['Average SSIM is ',num2str(mean(eval(['SSIM_',setTestCur,'_x',num2str(sf)])),'%2.4f')]);

    %% save PSNR and SSIM results
    save(fullfile(folderResultCur,['PSNR_',setTestCur,'_x',num2str(sf),'.mat']),['PSNR_',setTestCur,'_x',num2str(sf)]);
    save(fullfile(folderResultCur,['SSIM_',setTestCur,'_x',num2str(sf),'.mat']),['SSIM_',setTestCur,'_x',num2str(sf)]);

end

author: cszn

repo: SRMD

save_path: github-repos/MATLAB/cszn-SRMD


text:
function c = conv(c,varargin)

% c = CONV(c) convolves all traces with the final trace. By default the new
% correlation object has traces that are 2*N-1 samples long, where the
% original traces are of length N. Because all traces ...

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

% READ & CHECK ARGUMENTS
if (nargin>2)
    error('Wrong number of inputs');
end;

if ~isa(c,'correlation')
    error('First input must be a correlation object');
end

% GENERAL PARAMETERS
c = verify(c);
traces = get(c,'Traces');
keyTrace = traces;

X = fft(double(c.W));
Y = repmat( fft(double(c.W(keyTrace))) , 1 , traces );

save

Z = ifft(X.*conj(Y));
for n = 1:traces
    c.W(n) = set( c.W(n) , 'DATA' , Z(:,n) );
end

author: geoscience-community-codes

repo: GISMO


text:
function Z = ncvar_replace(X,Y,W)

% Call recurisvely for vector X
if length(X) > 1
    dimX = size(X);
    Z = [];
    X = X(:);
    for i = 1:length(X)
        Z = [Z;ncvar_replace(X(i),Y,W)];
    end
    Z = reshape(Z,dimX);
    return
end

% ...and recursively for vector Y
if length(Y) > 1
    Z = X;
    for i = 1:length(Y)
        Z = ncvar_replace(Z,Y(i),W(i));
    end
    return
end

nonCommutingTable         = yalmip('nonCommutingTable');
[monomtable,variabletype] = yalmip('monomtable');
if size(monomtable,1)>size(nonCommutingTable,1)
    nonCommutingTable((1+size(nonCommutingTable,1)):(size(monomtable,1)),1) = (1+size(nonCommutingTable,1)):(size(monomtable,1));
end

% Cast commutative variables as nc temporarily by adding them to the table
commuting = find(~any(nonCommutingTable,2));
nonCommutingTable(commuting,1) = commuting;

base = getbase(X);
vars = getvariables(X);
yvar = getvariables(Y);
Z = X;
newMonoms = [1];
for i = vars
    monTerms = nonCommutingTable(i,:);

    % Possibilities
    % [v(x) NC]
    % [NAN NC]
    if isa(Y,'sdpvar')
        if isnan(monTerms(1))
            % Purely non-commuting so no replacement can occur
            newMonoms = [newMonoms;recover(i)];
        else
            % call recursively with sdpvar replace
            if any(isnan(nonCommutingTable(monTerms(find(monTerms)))))
                [ct,vt] = coefficients(recovernc(i),Y);
            else
                [ct,vt] = coefficients(recover(i),Y);
            end
            newMonomial = ct;
            if degree(vt)>=1
                newMonomial = newMonomial*W^degree(vt);
            end
            newMonoms = [newMonoms;newMonomial];
        end
    else
        % Replacing a non-commuting with something else
        newMonomial = 1;
        for j = 2:max(find(monTerms))
            if monTerms(j)==yvar
                newMonomial = newMonomial*W;
            else
                newMonomial = newMonomial*recover(monTerms(j));
            end
        end
        newMonoms = [newMonoms;newMonomial];
    end
end
Z = reshape(base*newMonoms,size(X));

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
%
% rj(x, y, z, p, errtol)
%
% Inputs:
%
%   x       Input vector size 1xN.
%   y       Input vector size 1xN.
%   z       Input vector size 1xN.
%   p       Input vector size 1xN.
%   errtol  Error tolerance.
%
% Matlab function to compute Carlson's symmetric elliptic integral Rj.
% Implementation of Carlson's Duplication Algorithm 3 in "Computing
% Elliptic Integrals by Duplication," by B. C. Carlson, Numer. Math.
% 33, 1-16 (1979).
%
% Returns NaN's for any argument values outside input range.
%
% Algorithm is also from Carlson's ACM TOMS Algorithm 577.
%
% This code is a complete rewrite of the algorithm in vectorized form.
% It was not produced by running a FORTRAN to Matlab converter.
%
% The following text is copied from ACM TOMS Algorithm 577 FORTRAN code:
%
%   X AND Y ARE THE VARIABLES IN THE INTEGRAL RC(X,Y).
%
%   ERRTOL IS SET TO THE DESIRED ERROR TOLERANCE.
%   RELATIVE ERROR DUE TO TRUNCATION IS LESS THAN
%   16 * ERRTOL ** 6 / (1 - 2 * ERRTOL).
%
%   SAMPLE CHOICES:  ERRTOL     RELATIVE TRUNCATION
%                               ERROR LESS THAN
%                    1.D-3      4.D-18
%                    3.D-3      3.D-15
%                    1.D-2      4.D-12
%                    3.D-2      3.D-9
%                    1.D-1      4.D-6
%
% Note by TRH:
%
%   Absolute truncation error when the integrals are order 1 quantities
%   is closer to errtol, so be careful if you want high absolute precision.
%
% Thomas R. Hoffend Jr., Ph.D.
% 3M Company
% 3M Center Bldg. 236-GC-26
% St. Paul, MN 55144
% trhoffendjr@mmm.com
%

function f = rj(x, y, z, p, errtol)

% Argument limits as set by Carlson:
LoLim = 5.0 * realmin;
UpLim = 5.0 * realmax;

% Check input arguments for acceptability:
mask = (min([x; y; z]) >= 0) & ...
       (min([(x + y); (x + z); (y + z); p]) >= LoLim) & ...
       (max([x; y; z; p]) < UpLim);

% Define internally acceptable variable ranges for iterations:
Xi = x(mask);
Yi = y(mask);
Zi = z(mask);
Pi = p(mask);

% Carlson's duplication algorithm for Rj:
Xn = Xi;
Yn = Yi;
Zn = Zi;
Pn = Pi;
sigma = 0.0;
power4 = 1.0;
etolrc = 0.5 * errtol;
Mu = (Xn + Yn + Zn + Pn + Pn) * 0.2;
Xndev = (Mu - Xn) ./ Mu;
Yndev = (Mu - Yn) ./ Mu;
Zndev = (Mu - Zn) ./ Mu;
Pndev = (Mu - Pn) ./ Mu;
epslon = max( abs([Xndev Yndev Zndev Pndev]) );
while (epslon >= errtol)
    Xnroot = sqrt(Xn);
    Ynroot = sqrt(Yn);
    Znroot = sqrt(Zn);
    lambda = Xnroot .* (Ynroot + Znroot) + Ynroot .* Znroot;
    alpha = Pn .* (Xnroot + Ynroot + Znroot) + Xnroot .* Ynroot .* Znroot;
    alpha = alpha .* alpha;
    beta = Pn .* (Pn + lambda) .* (Pn + lambda);
    sigma = sigma + power4 .* rc(alpha, beta, etolrc);
    % Here we might need to shrink the size of the arrays if alpha, beta
    % out of range for Rc.
    mask = ~isnan(sigma);
    power4 = 0.25 * power4;
    Xn = 0.25 * (Xn(mask) + lambda(mask));
    Yn = 0.25 * (Yn(mask) + lambda(mask));
    Zn = 0.25 * (Zn(mask) + lambda(mask));
    Pn = 0.25 * (Pn(mask) + lambda(mask));
    Mu = (Xn + Yn + Zn + Pn + Pn) * 0.2;
    Xndev = (Mu - Xn) ./ Mu;
    Yndev = (Mu - Yn) ./ Mu;
    Zndev = (Mu - Zn) ./ Mu;
    Pndev = (Mu - Pn) ./ Mu;
    epslon = max( abs([Xndev Yndev Zndev Pndev]) );
end
C1 = 3.0 / 14.0;
C2 = 1.0 / 3.0;
C3 = 3.0 / 22.0;
C4 = 3.0 / 26.0;
EA = Xndev .* Yndev - Zndev .* Zndev;
EB = Xndev .* Yndev .* Zndev;
EC = Pndev .* Pndev;
E2 = EA - 3.0 * EC;
E3 = EB + 2.0 * Pndev .* (EA - EC);
S1 = 1.0 + E2 .* (-C1 + 0.75D0 * C3 * E2 - 1.5D0 * C4 * E3);
S2 = EB .* (0.5D0 * C2 + Pndev .* (-C3 -C3 + Pndev*C4));
S3 = Pndev .* EA .* (C2 - Pndev*C3) - C2*Pndev .* EC;
f(mask) = 3.0 * sigma + power4 * (S1 + S2 + S3) ./ (Mu .* sqrt(Mu));

% Return NaN's where input argument was out of range:
f(~mask) = NaN;


text:
%  INTERNAL FUNCTION: Computes shrinking and expansion objects for the manipulation of symmetric tensors
%
%  ::
%
%    [keep,expand,C,UC,B]=shrink_expand(n,k)
%    [keep,expand,C,UC,B]=shrink_expand(n,k,strategy)
%    [keep,expand,C,UC,B]=shrink_expand(n,k,strategy)
%
%  Args:
%
%     - **n** [scalar] : number of variables in the tensor
%     - **k** [scalar] : order of the tensor
%     - **strategy** ['up'|{'down'}] : alternative compression strategies
%
%       - 'up': i1<=i2<=i3....<=in
%       - 'down': i1>=i2>=i3....>=in
%
%  Returns:
%     :
%
%     - **keep** [cell array of k logical vectors] : n^k x 1 vector true for
%       the columns to be kept
%     - **expand** [cell array of k vectors] : n^k x 1 vector replicating the
%       compressed columns to form the grand tensor
%     - **C** [cell array of k matrices] : sparse compression matrix of size
%       n^k x g, where g=nchoosek(n+k-1,k) is the number of unique elements in
%       the tensor (matrix version of **keep**)
%     - **UC** [cell array of k matrices] : sparse expansion matrix of size
%       g x n^k, where g=nchoosek(n+k-1,k) is the number of unique elements in
%       the tensor (matrix version of **expand**)
%
%  Note:
%
%     - Essentially the results are given for orders from 1 to k
%     - useful for shrinking tensors of the form fvvv...v as used in higher-order
%       differentiation
%
%

author: jmaih


text:
function v=hypervolume(P,r,N)
% HYPERVOUME    Hypervolume indicator as a measure of Pareto front estimate.
%   V = HYPERVOLUME(P,R,N) returns an estimation of the hypervoulme (in
%   percentage) dominated by the approximated Pareto front set P (n by d)
%   and bounded by the reference point R (1 by d). The estimation is doen
%   through N (default is 1000) uniformly distributed random points within
%   the bounded hyper-cuboid.
%
%   V = HYPERVOLUMN(P,R,C) uses the test points specified in C (N by d).
%
% See also: paretofront, paretoGroup

% Version 1.0 by Yi Cao at Cranfield University on 20 April 2008

% Example
%{
% an random exmaple
F=(randn(100,3)+5).^2;
% upper bound of the data set
r=max(F);
% Approximation of Pareto set
P=paretofront(F);
% Hypervolume
v=hypervolume(F(P,:),r,100000);
%}

% Check input and output
error(nargchk(2,3,nargin));
error(nargoutchk(0,1,nargout));

P=P*diag(1./r);
[n,d]=size(P);
if nargin<3
    N=1000;
end
if ~isscalar(N)
    C=N;
    N=size(C,1);
else
    C=rand(N,d);
end

fDominated=false(N,1);
lB=min(P);
fcheck=all(bsxfun(@gt, C, lB),2);

for k=1:n
    if any(fcheck)
        f=all(bsxfun(@gt, C(fcheck,:), P(k,:)),2);
        fDominated(fcheck)=f;
        fcheck(fcheck)=~f;
    end
end

v=sum(fDominated)/N;


text:
% EXPeriment NAME
% returns the name of the current experiment
% this name indexes the CONFIGDATA
%
% $Id: expname.m,v 2.7 2005/05/23 16:23:35 svoboda Exp $

function name = expname()

% the name of a calibration fo BlueC must
% contain the following sub-string 'BigBlue'

% the name for the Hoenggeberg calibration must contain
% the sub-string Hoengg

% the name of the oscar setup must contain
% the string 'oscar'

name = 'strawlab_test';
%name = 'flydra';
%name = 'mamarama';
%name = 'humdra';
% name = '0801BlueCRZ';
% name = 'BlueCHoengg';
% name = 'TestData';
% name = 'OneSideH0105';
% name = 'OneSideH1416';
% name = 'BlueCRZ';
% name = 'ViRoom20030611';
% name = '2410ViRoom';
% name = 'Erlangen';
% name = 'oscar2c1p';
% name = 'oscardemo';
% name = 'ViRoom20030724';
% name = 'extern';

return

author: strawlab

repo: MultiCamSelfCal


text:
function spacing_matters ( )

%*****************************************************************************80
%
%  SPACING_MATTERS demonstrates a pernicious problem in defining vectors.
%
%  Discussion:
%
%    I don't like cramped coding styles that make it seem as though blank space
%    is wasted space.  So I prefer to expand function calls so that the name
%    of the function, the parentheses, and the arguments are all separated by
%    blank space.
%
%    I was doing this to some code that someone had given me when the code
%    suddenly stopped working.  The error message from MATLAB's code analyzer
%    was characteristically stern, cryptic and of no help whatsoever:
%
%      "Expression or statement is incorrect--possibly unbalanced (, {, or [."
%
%    It took me an hour to figure out what was happening, and I thought I
%    had a defective copy of MATLAB when I realized that I could make the
%    error come or go by adding or removing a blank space!
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 April 2011
%
%  Author:
%
%    John Burkardt
%

%
%  Question:
%
%    Which of the matrices A, B, C, D, E, and F is formed incorrectly,
%    according to MATLAB?
%

%
%  The following statement is illegal, according to MATLAB.
%
  a = [ ones ( 10, 1 ); zeros ( 10, 1 ) ];
%
%  The following, seemingly equivalent statement, is legal.
%
  b = [ ones(10,1); zeros(10,1) ];
%
%  This statement is legal.
%
  c = [ ones( 10, 1 ); zeros( 10, 1 ) ];
%
%  This statement is NOT legal!
%
  d = [ ones (10,1); zeros (10,1) ];
%
%  These statements, which can be regarded as equivalent to the definition
%  of A after "tokenizing" certain strings, are legal!
%
  e1 = ones ( 10, 1 );
  e2 = zeros ( 10, 1 );
  e = [ e1; e2 ];
%
%  And of course, these are legal.
%
  f1 = ones(10,1);
  f2 = zeros(10,1);
  f = [ f1; f2 ];
%
%  What is the problem?
%  MATLAB sees the words [ ones ( 10, 1 )... as an attempt to define
%  a vector out of a variable calles "ones", followed by an object ( 10, 1 )
%  that it can't interpret.  That's because MATLAB allows you to enter the
%  entries of a vector without separating commas.
%
%  Thus, we are allowed to say:
%
%    a = [ 1 2 3 ]
%
%  but this illusory advantage comes at the disastrous cost that things
%  that can SOMETIMES be separated by spaces CANNOT be separated by spaces
%  when occuring in a literal vector.  Strangely enough, MATLAB has no
%  problem correctly handling
%
%    g = [ 1 + 2  4 * 5 ]
%
  return
end


text:
function u=v_glotros(d,t,p)
%V_GLOTROS  Rosenberg glottal model U=(D,T,P)
% d is derivative of flow waveform
% t is in fractions of a cycle
% p has parameters
%	p(1)=closure time
%	p(2)=+ve/-ve slope ratio

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

if nargin < 2
  tt=(0:99)'/100;
else
  tt=mod(t,1);
end
u=zeros(size(tt));
de=[0.6 0.5]';
if nargin < 3
  p=de;
elseif length(p)<2
  p=[p(:); de(length(p)+1:2)];
end
pp=p(1)/(1+p(2));
ta=tt<pp;
tb=tt<p(1) & ~ta;
wa=pi/pp;
wb=0.5*pi/(p(1)-pp);
fb=wb*pp;
if d==0
  u(ta)=0.5*(1-cos(wa*tt(ta)));
  u(tb)=cos(wb*tt(tb)-fb);
elseif d==1
  u(ta)=0.5*wa*sin(wa*tt(ta));
  u(tb)=-wb*sin(wb*tt(tb)-fb);
elseif d==2
  u(ta)=0.5*wa^2.*cos(wa*tt(ta));
  u(tb)=-wb^2*cos(wb*tt(tb)-fb);
else
  error('Derivative must be 0,1 or 2');
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function match = birthday(groupsize)
% BIRTHDAY Simulates a single trial of the Birthday Paradox
%    MATCH = BIRTHDAY(GROUPSIZE) creates a randomly-selected birthday for
%    every member of a group of size GROUPSIZE and tests whether any of
%    the selected birthdays match.  MATCH is 1 if two or more members of
%    the group share the same birthday and 0 otherwise.
%
%    Example:
%    >> match = birthday(30)

% Match is zero until a birthday match is found
match = 0;

% Initialize list of taken birthdays
bdaylist = zeros(1, groupsize);

for person = 1:groupsize

    % Randomly select a birthdate for the individual (ignore leap years)
    birthdate = ceil(365*rand());

    % Check if someone else in the group shares the same birthday
    if (any(birthdate == bdaylist))

        % A match is found, return from the function
        match = 1;
        return;

    end

    % Add the birthdate to the list for the group
    bdaylist(person) = birthdate;

end

repo: Aero-matlab


text:
function display(X)
%DISPLAY   Display a BALLFUNV.
%   DISPLAY(F) outputs important information about the BALLFUNV F to the
%   command window, including its domain of definition and a summary of its
%   structure.
%
%   It is called automatically when the semicolon is not used at the end of a
%   statement that results in a BALLFUNV.

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

if ( isequal(get(0, 'FormatSpacing'), 'compact') )
	disp([inputname(1), ' =']);
	disp(X);
else
	disp(' ');
	disp([inputname(1), ' =']);
	disp(' ');
    disp(X);
end

end


text:
report_this_filefun(mfilename('fullpath'));

[file1,path1] = uigetfile([ '*.mat'],'Grid File (Referenz)');
load([path1 file1])
bvg1=bvg;

[file1,path1] = uigetfile([ '*.mat'],'Grid File');
load([path1 file1])
bvg2=bvg;

bvg = [bvg2(:,1)-bvg1(:,1) , bvg2(:,2)-bvg1(:,2) ,...
    bvg2(:,3)-bvg1(:,3) , bvg2(:,4)-bvg1(:,4) , ...
    bvg2(:,5)-bvg1(:,5) , bvg2(:,6)-bvg1(:,6) ,...
    bvg2(:,7)-bvg1(:,7) , bvg2(:,8)-bvg1(:,8) ,...
    bvg2(:,9)-bvg1(:,9) , bvg2(:,10)-bvg1(:,10) ,...
    bvg2(:,11)-bvg1(:,11) ];

normlap2=ones(length(tmpgri(:,1)),1)*nan;
normlap2(ll)= bvg(:,1);
re3=reshape(normlap2,length(yvect),length(xvect));

normlap2(ll)= bvg(:,5);
r=reshape(normlap2,length(yvect),length(xvect));

normlap2(ll)= bvg(:,6);
meg=reshape(normlap2,length(yvect),length(xvect));

normlap2(ll)= bvg(:,2);
old1=reshape(normlap2,length(yvect),length(xvect));

normlap2(ll)= bvg(:,7);
pro=reshape(normlap2,length(yvect),length(xvect));

normlap2(ll)= bvg(:,8);
avm=reshape(normlap2,length(yvect),length(xvect));

normlap2(ll)= bvg(:,9);
stanm=reshape(normlap2,length(yvect),length(xvect));

normlap2(ll)= bvg(:,10);
Prmap=reshape(normlap2,length(yvect),length(xvect));

old = re3;

% View the b-value map
view_bva_newwin

author: CelsoReyes


text:
function bpath_average ( seed )

%*****************************************************************************80
%
%% BPATH_AVERAGE: displays the average of 1000 Brownian paths.
%
%  Discussion:
%
%    This routine computes 1000 simulations of discretized Brownian
%    motion W(t) over the time interval [0,1] using 500 time steps.
%    The user specifies a random number seed.  Different values of
%    the seed will result in a different set of realizations of the path.
%
%    Actually, we are interested in a function u(W(t)):
%
%      u(W(t)) = exp ( t + W(t)/2 )
%
%    The routine plots 5 of the simulations, as well as the average
%    of all the simulations.
%
%    The plot of the average should be quite smooth.  Its expected
%    value is exp ( 9 * t / 8 ), and we compute the 'error', that is,
%    the difference between the averaged value and this expected
%    value.  This 'error' should decrease as the number of simulation
%    is increased.
%
%  Modified:
%
%    27 September 2006
%
%  Author:
%
%    Desmond Higham
%
%  Reference:
%
%    Desmond Higham,
%    An Algorithmic Introduction to Numerical Simulation of
%    Stochastic Differential Equations,
%    SIAM Review,
%    Volume 43, Number 3, September 2001, pages 525-546.
%
%  Parameters:
%
%    Input, integer SEED, a seed for the random number generator.
%

%
%  We will be computing M paths simultaneously.
%
  m = 1000;
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BPATH_AVERAGE\n' );
  fprintf ( 1, '  Average %d Brownian path simulations.\n', m );

  if ( nargin < 1 )
    seed = 123456789;
    fprintf ( 1, '  Using default value of SEED = %d\n', seed );
  end
%
%  Initialize the random number generator.
%
  rng ( seed )

  tmax = 1.0;
  n = 500;
  dt = tmax / n;
  t = [ dt: dt: 1 ];
%
%  Begin timing.
%
  tic
%
%  Set the increments.
%
  dw = sqrt ( dt ) * randn ( m, n );
%
%  Sum the increments.
%
  w0 = 0.0;
  w = cumsum ( dw, 2 );

  u = exp ( repmat ( t, [m 1] ) + 0.5 * w );

  umean = mean ( u );
%
%  End timing.
%
  toc
%
%  Plot the mean of the M paths.
%
  plot ( [0,t], [1,umean], 'b-', 'Linewidth', 2 );
  hold on
%
%  Overlay the plots of the first 5 paths.
%
  plot ( [0,t], [ones(5,1), u(1:5,:)], 'r--' );
  grid on
  hold off

  xlabel ( 't', 'FontSize', 16 )
  ylabel ( 'U(t)', 'FontSize', 16, 'Rotation', 0, 'HorizontalAlignment', 'right' )
  legend ( 'mean of 1000 paths', '5 individual paths', 2 )
  title ( '1000 paths computed by BPATH_AVERAGE' )

  filename = 'bpath_average.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Plot saved as "%s"\n', filename );
%
%  Sample error:
%
  averr = norm ( ( umean - exp ( 9 * t / 8 ) ), 'inf' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Maximum error in averaged data is %f\n', averr );

  return
end


text:
function value = f_05 ( x )

%*****************************************************************************80
%
%% F_05 evaluates (x+3)*(x-1)*(x-1).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 April 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the point at which F is to be evaluated.
%
%    Output, real VALUE, the value of the function at X.
%
  value = ( x + 3.0 ) .* ( x - 1.0 ) .* ( x - 1.0 );

  return
end


function [X,Y] = pairsToLabels(pairs, instance_matrix)
    imgs1 = [pairs.img1];
    imgs2 = [pairs.img2];

    keymat = [[imgs1.classId] [imgs2.classId]; [imgs1.id] [imgs2.id]];
    [ids, idsndx] = unique(keymat','rows');

    X = instance_matrix(:,idsndx);
    Y = ids(:,1)';
end


text:
function update_plots_nn_ps
global hnds
global r nn asz G ps dsm
global hpb ht hi

ax1=hnds.axes1;
ax2=hnds.axes2;

cla(ax1);
cla(ax2);

nn=str2num(get(hnds.nn,'string')); % number of cities
ps=str2num(get(hnds.ps,'string')); % population size

r=asz*rand(2,nn); % randomly distribute cities
% r(1,:) -x coordinaties of cities
% r(2,:) -y coordinaties of cities

dsm=zeros(nn,nn); % matrix of distancies
for n1=1:nn-1
    r1=r(:,n1);
    for n2=n1+1:nn
        r2=r(:,n2);
        dr=r1-r2;
        dr2=dr'*dr;
        drl=sqrt(dr2);
        dsm(n1,n2)=drl;
        dsm(n2,n1)=drl;
    end
end

% start from random closed pathes:
G=zeros(ps,nn); % genes, G(i,:) - gene of i-path, G(i,:) is row-vector with cities number in the path
for psc=1:ps
    G(psc,:)=randperm(nn);
end

% to plot best path:
hpb=plot(NaN,NaN,'r-','parent',ax1);
ht=title(ax1,' ');

set(ax1,'NextPlot','add');

% plot nodes numbers
for n=1:nn
    text(r(1,n),r(2,n),num2str(n),'color',[0.7 0.7 0.7],'parent',ax1);
end

plot(r(1,:),r(2,:),'k.','parent',ax1); % plot cities as black dots

axis(ax1,'equal');
set(ax1,'xlim',[-0.1*asz 1.1*asz]);
set(ax1,'ylim',[-0.1*asz 1.1*asz]);

hi=imagesc(G,'parent',ax2);
title(ax2,'color is city number');
colorbar('peer',ax2);
xlabel(ax2,'index in sequence of cities');
ylabel(ax2,'path number');


text:
function [M,K,C,P,phi]= ModalAnalysis(M,K,C,P,lambda)
% Modal Analysis of the System
%-------------------------------------------------------------------------
% Code written by :Siva Srinivas Kolukula                                 |
%                  Senior Research Fellow                                 |
%                  Structural Mechanics Laboratory                        |
%                  Indira Gandhi Center for Atomic Research               |
%                  INDIA                                                  |
% E-mail : allwayzitzme@gmail.com                                         |
%-------------------------------------------------------------------------
%--------------------------------------------------------------------------
% Purpose : To do the Modal Analysis of the system
%
% Synopsis : [M,K,C,P] = ModalAnalysis(M,K,C,P,lambda]
%
% Variable Description :
%           M - Mass Matrix
%           K - Stiffness Matrix
%           C - Damping Matrix
%           P - Force MAtrix
%           lambda - number of modes to be considered
%           phi - Reduced Normal modes of the system
%--------------------------------------------------------------------------

 % Modal Analysis
[V,D]=eig(K,M);
[W,k]=sort(diag(D));
V=V(:,k);
Factor=diag(V'*M*V);
Phi=V*inv(sqrt(diag(Factor)));
Omega=diag(sqrt(Phi'*K*Phi));

% selecting only first two Natural Frequencies and Mode shapes
w = diag(Omega(1:lambda)) ;
phi = Phi(:,1:lambda) ;

% Reducing the Degree's of Freedom of M, K and P
M = phi'*M*phi ;
C = phi'*C*phi ;
K = phi'*K*phi ;
P = phi'*P ;


text:
function [delay,weight] = driving_function_imp_wfs_ps(x0,nx0,xs,conf)
%DRIVING_FUNCTION_IMP_WFS_PS weights and delays for a point source in WFS
%
%   Usage: [delay,weight] = driving_function_imp_wfs_ps(x0,nx0,xs,conf)
%
%   Input parameters:
%       x0      - position  of secondary sources / m [nx3]
%       nx0     - direction of secondary sources [nx3]
%       xs      - position of point source [nx3]
%       conf    - configuration struct (see SFS_config)
%
%   Output parameters:
%       delay   - delay of the driving function / s
%       weight  - weight (amplitude) of the driving function
%
%   See also: sound_field_imp, sound_field_imp_wfs, driving_function_mono_wfs_ps
%
%   References:
%       Start (1997) - "Direct Sound Enhancement by Wave Field Synthesis",
%       TU Delft,
%       http://resolver.tudelft.nl/uuid:c80d5b58-67d3-4d84-9e73-390cd30bde0d
%
%       Wierstorf (2014) - "Perceptual Assessment of Sound Field Synthesis",
%       TU Berlin, https://doi.org/10.14279/depositonce-4310

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 4;
nargmax = 4;
narginchk(nargmin,nargmax);
isargmatrix(x0,nx0,xs);
isargstruct(conf);

%% ===== Configuration ==================================================
% Speed of sound
c = conf.c;
xref = conf.xref;
fs = conf.fs;
dimension = conf.dimension;
driving_functions = conf.driving_functions;

%% ===== Computation =====================================================
% Get the delay and weighting factors
if strcmp('2D',dimension) || strcmp('3D',dimension)

    % === 2- or 3-Dimensional ============================================

    switch driving_functions
    case 'default'
        % --- SFS Toolbox ------------------------------------------------
        % d using a point source as source model
        %
        %                   1  (x0-xs) nx0
        % d(x0,t) = h(t) * --- ----------- delta(t-|x0-xs|/c)
        %                  2pi  |x0-xs|^2
        %
        % https://sfs.rtfd.io/en/3.2/d_wfs/#equation-td-wfs-point
        %
        % r = |x0-xs|
        r = vector_norm(x0-xs,2);
        % Delay and amplitude weight
        delay = 1./c .* r;
        weight = 1./(2.*pi) .* vector_product(x0-xs,nx0,2) ./ r.^2;
        %
    case 'legacy'
        % --- Old SFS Toolbox default ------------------------------------
        % d using a point source as source model
        %
        %                   1   (x0-xs) nx0
        % d(x0,t) = h(t) * --- ------------- delta(t-|x0-xs|/c)
        %                  2pi |x0-xs|^(3/2)
        %
        % See Wierstorf (2014), eq.(2.63)
        %
        % r = |x0-xs|
        r = vector_norm(x0-xs,2);
        % Delay and amplitude weight
        delay = 1/c .* r;
        weight = 1/(2*pi) .* vector_product(x0-xs,nx0,2) ./ r.^(3/2);
        %
    otherwise
        error(['%s: %s, this type of driving function is not implemented', ...
            'for a point source.'],upper(mfilename),driving_functions);
    end

elseif strcmp('2.5D',dimension)

     % === 2.5-Dimensional ================================================

     % Reference point
     xref = repmat(xref,[size(x0,1) 1]);
     switch driving_functions
     case {'default', 'reference_point'}
         % Driving function with only one stationary phase approximation,
         % reference to one point in field
         %
         % r = |x0-xs|
         r = vector_norm(x0-xs,2);
         % 2.5D correction factor
         %         _____________________
         %        |      |xref-x0|
         % g0 = _ |---------------------
         %       \| |xref-x0| + |x0-xs|
         %
         % See Start (1997), eq. (3.11)
         %
         g0 = sqrt( vector_norm(xref-x0,2) ./ (vector_norm(x0-xref,2) + r) );
         %                                 ___
         %                                | 1    (x0-xs) nx0
         % d_2.5D(x0,t) = h_pre(t) * g0 _ |---  ------------- delta(t-|x0-xs|/c)
         %                               \|2pi  |x0-xs|^(3/2)
         %
         % https://sfs.rtfd.io/en/3.2/d_wfs/#equation-td-wfs-point-25d
         %
         % Delay and amplitude weight
         delay = 1./c .* r;
         weight = g0 ./ sqrt(2.*pi) .* vector_product(x0-xs,nx0,2) ./ r.^(3/2);
         %
     case 'reference_line'
         % Driving function with two stationary phase approximations,
         % reference to a line parallel to a LINEAR secondary source distribution
         %
         % Distance ref-line to linear ssd
         dref = vector_product(xref-x0,nx0,2);
         % Distance source and linear ssd
         ds = vector_product(xs-x0,nx0,2);
         %
         % 2.5D correction factor
         %        _______________________
         % g0 = \| d_ref / (d_ref - d_s)
         %
         % See Start (1997), eq. (3.16)
         %
         g0 = sqrt( dref ./ (dref - ds) );
         %                                 ___
         %                                | 1    (x0-xs) nx0
         % d_2.5D(x0,t) = h_pre(t) * g0 _ |---  ------------- delta(t-|x0-xs|/c)
         %                               \|2pi  |x0-xs|^(3/2)
         %
         % Inverse Fourier Transform of Start (1997), eq. (3.17)
         %
         % r = |x0-xs|
         r = vector_norm(x0-xs,2);
         % Delay and amplitude weight
         delay = 1./c .* r;
         weight = g0 ./ sqrt(2.*pi) .* vector_product(x0-xs,nx0,2) ./ r.^(3./2);
         %
     case 'legacy'
         % --- SFS Toolbox ------------------------------------------------
         % 2.5D correction factor
         %        ______________
         % g0 = \| 2pi |xref-x0|
         %
         g0 = sqrt(2*pi*vector_norm(xref-x0,2));
         %
         % d_2.5D using a point source as source model
         %
         %                        g0  (x0-xs) nx0
         % d_2.5D(x0,t) = h(t) * --- ------------- delta(t-|x0-xs|/c)
         %                       2pi |x0-xs|^(3/2)
         %
         % See Wierstorf (2014), eq.(2.64)
         %
         % r = |x0-xs|
         r = vector_norm(x0-xs,2);
         % Delay and amplitude weight
         delay = 1./c .* r;
         weight = g0 ./ (2.*pi) .* vector_product(x0-xs,nx0,2) ./ r.^(3./2);
         %
     otherwise
         error(['%s: %s, this type of driving function is not implemented', ...
           'for a 2.5D point source.'],upper(mfilename),driving_functions);
     end

else
    error('%s: the dimension %s is unknown.',upper(mfilename),dimension);
end

author: sfstoolbox

repo: sfs-matlab


text:
function [U,output] = cpd_nls(T,U0,options)
%CPD_NLS CPD by nonlinear least squares.
%   [U,output] = cpd_nls(T,U0) computes the factor matrices U{1}, ..., U{N}
%   belonging to a canonical polyadic decomposition of the N-th order
%   tensor T by minimizing 0.5*frob(T-cpdgen(U))^2. The algorithm is
%   initialized with the factor matrices U0{n}. The structure output
%   returns additional information:
%
%      output.Name  - The name of the selected algorithm.
%      output.<...> - The output of the selected algorithm.
%
%   cpd_nls(T,U0,options) may be used to set the following options:
%
%      options.Algorithm =   - The desired optimization method.
%      [@nls_gncgs| ...
%       {@nls_gndl}|@nls_lm]
%      options.LargeScale    - If true, the Gauss-Newton or Levenberg-
%      = sum(size(T))*R>1e2    Marquardt steps are computed using a
%                              preconditioned conjugate gradient algorithm.
%                              Otherwise, a direct solver is used.
%      options.M =           - The preconditioner to use when
%      [{'block-Jacobi'}|...   options.LargeScale is true.
%       'block-SSOR'|false]
%      options.PlaneSearch = - A function handle to the desired CPD plane
%      [{false},@cpd_eps]      search algorithm. Only applicable to
%                              algorithms with trust region globalization.
%      options.<...>         - Parameters passed to the selected method,
%                              e.g., options.TolFun, options.TolX and
%                              options.PlaneSearchOptions. See also help
%                              [options.Algorithm].
%
%   See also cpd_minf.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)
%
%   References:
%   [1] L. Sorber, M. Van Barel, L. De Lathauwer, "Optimization-based
%       algorithms for tensor decompositions: canonical polyadic
%       decomposition, decomposition in rank-(Lr,Lr,1) terms and a new
%       generalization," SIAM J. Opt., 2013.
%   [2] L. Sorber, M. Van Barel, L. De Lathauwer, "Unconstrained
%       optimization of real functions in complex variables," SIAM J. Opt.,
%       Vol. 22, No. 3, 2012, pp. 879-898.

% Format the tensor T.
T = fmt(T,true);
if isstruct(T), size_tens = T.size;
else size_tens = [size(T) ones(1,length(U0)-ndims(T))]; end

% Check the initial factor matrices U0.
N = length(U0);
R = size(U0{1},2);
if any(cellfun('size',U0,2) ~= R)
    error('cpd_nls:U0','size(U0{n},2) should be the same for all n.');
end

% Check the options structure.
isfunc = @(f)isa(f,'function_handle');
xsfunc = @(f)isfunc(f)&&exist(func2str(f),'file');
if nargin < 3, options = struct; end
if ~isfield(options,'Algorithm')
    funcs = {@nls_gndl,@nls_gncgs,@nls_lm};
    options.Algorithm = funcs{find(cellfun(xsfunc,funcs),1)};
end
if ~isfield(options,'CGMaxIter'), options.CGMaxIter = 10; end
if ~isfield(options,'Display'), options.Display = 0; end
if ~isfield(options,'JHasFullRank'), options.JHasFullRank = false; end
if ~isfield(options,'LargeScale')
    options.LargeScale = sum(size_tens)*R > 1e2;
end
if ~isfield(options,'M'), options.M = 'block-Jacobi'; end
if ~isfield(options,'TolLargeScale'), options.TolLargeScale = 0.02; end

% Adapt line/plane search if it is a CPD line/plane search.
if isfield(options,'LineSearch') && isfunc(options.LineSearch) && ...
   ~isempty(strfind(func2str(options.LineSearch),'cpd_'))
    linesearch = options.LineSearch;
    options.LineSearch = @ls;
end
if isfield(options,'PlaneSearch') && isfunc(options.PlaneSearch) && ...
   ~isempty(strfind(func2str(options.PlaneSearch),'cpd_'))
    planesearch = options.PlaneSearch;
    options.PlaneSearch = @ps;
end

% Call the optimization method.
cache.offset = cumsum([1 cellfun(@numel,U0(:).')]);
cache.T2 = frob(T,'squared');
if options.LargeScale, dF.JHJx = @JHJx; else dF.JHJ = @JHJ; end
dF.JHF = @grad;
switch options.M
    case 'block-SSOR', dF.M = @M_blockSSOR;
    case 'block-Jacobi', dF.M = @M_blockJacobi;
    case 'Jacobi', dF.M = @M_Jacobi;
    otherwise, if isa(options.M,'function_handle'), dF.M = options.M; end
end
state(U0,true);
[U,output] = options.Algorithm(@objfun,dF,U0(:).',options);
output.Name = func2str(options.Algorithm);

function state(z,firstrun)

    if nargin == 2 && firstrun
        % Store the fraction of known elements.
        if isstruct(T) && T.incomplete
            cache.scale = length(T.val)./prod(T.size);
        end
    end

    % Cache the factor matrices' Gramians.
    cache.UHU = zeros(N,R*R);
    for n = 1:N
        tmp = conj(z{n}'*z{n});
        cache.UHU(n,:) = tmp(:);
    end

    % Optionally cache the inverses of the Gramians for the preconditioner.
    % In a faster language, this should be the Cholesky factor instead.
    if ischar(options.M) || isa(options.M,'function_handle')
        cache.invW = cell(1,N);
        for n = 1:N
            tmp = cache.UHU([1:n-1 n+1:N],:);
            if N > 2, tmp = prod(tmp,1); end
            cache.invW{n} = inv(reshape(tmp,[R R]));
        end
    end

end

function fval = objfun(z)

    % CPD objective function.
    isincomplete = isstruct(T) && T.incomplete;
    issparse = isstruct(T) && T.sparse;
    if ~isincomplete || length(T.ind)/prod(T.size) > options.TolLargeScale
        fval = z{1}*kr(z(end:-1:2)).';
        if isincomplete, fval = fval(T.ind)-T.val;
        elseif issparse
            if ~isempty(T.ind), fval(T.ind) = fval(T.ind)-T.val; end
        else fval = fval-reshape(T,size(fval));
        end
    else
        fval = -T.val;
        for r = 1:R
            tmp = z{1}(T.sub{1},r);
            for n = 2:length(z), tmp = tmp.*z{n}(T.sub{n},r); end
            fval = fval+tmp;
        end
    end
    if isincomplete
        E = T; E.val = fval;
        if ~isempty(E.matrix)
            E.matrix = sparse(double(E.sub{1}), ...
                double(1+idivide(E.ind-1,int64(size(E.matrix,1)))), ...
                double(fval),size(E.matrix,1),size(E.matrix,2));
        end
        cache.residual = E;
    else
        cache.residual = reshape(fval,size_tens);
    end
    fval = 0.5*(fval(:)'*fval(:));

end

function grad = grad(z)

    % CPD scaled conjugate cogradient.
    state(z);
    E = cache.residual;
    offset = cache.offset;
    grad = nan(offset(end)-1,1);
    for n = 1:length(z)
        tmp = full(mtkrprod(E,z,n));
        grad(offset(n):offset(n+1)-1) = tmp(:);
    end

end

function JHJ = JHJ(z)

    % CPD Jacobian's Gramian.
    UHU = cache.UHU;
    JHJ = zeros(cache.offset(end)-1);
    for n = 1:N
        idxn = cache.offset(n):cache.offset(n+1)-1;
        Wn = reshape(prod(UHU([1:n-1 n+1:N],:),1),[R R]);
        JHJ(idxn,idxn) = kron(Wn,eye(size_tens(n)));
        for m = n+1:N
            idxm = cache.offset(m):cache.offset(m+1)-1;
            Wnm = reshape(prod(UHU([1:n-1 n+1:m-1 m+1:N],:),1),[R R]);
            JHJnm = bsxfun(@times,reshape(z{n},[size_tens(n) 1 1 R]), ...
                    reshape(conj(z{m}),[1 size_tens(m) R 1]));
            JHJnm = bsxfun(@times,JHJnm,reshape(Wnm,[1 1 R R]));
            JHJnm = permute(JHJnm,[1 3 2 4]);
            JHJnm = reshape(JHJnm,[size_tens(n)*R size_tens(m)*R]);
            JHJ(idxn,idxm) = JHJnm;
            JHJ(idxm,idxn) = JHJnm';
        end
    end

    % If incomplete, approximate the effect of missing entries.
    if isstruct(T) && T.incomplete
        JHJ = JHJ*cache.scale;
    end

end

function y = JHJx(z,x)

    % CPD fast Jacobian's Gramian vector product.
    % Ignores the fact that the tensor might be incomplete.
    offset = cache.offset;
    UHU = cache.UHU;
    XHU = zeros(R,R,N);
    y = nan(offset(end)-1,1);
    for n = 1:N
        Wn = UHU([1:n-1 n+1:N],:);
        if N > 2, Wn = prod(Wn,1); end
        tmp = reshape(x(offset(n):offset(n+1)-1),[],R);
        XHU(:,:,n) = conj(tmp'*z{n});
        y(offset(n):offset(n+1)-1) = tmp*reshape(Wn,[R R]);
    end
    for n = 1:N-1
        idxn = offset(n):offset(n+1)-1;
        Wn = zeros(R);
        for m = n+1:N
            idxm = offset(m):offset(m+1)-1;
            if N == 2
                Wn = Wn+XHU(:,:,m);
                JHJmnx = z{m}*XHU(:,:,n);
            else
                Wnm = UHU([1:n-1 n+1:m-1 m+1:N],:);
                if N > 3, Wnm = prod(Wnm,1); end
                Wnm = reshape(Wnm,[R R]);
                Wn = Wn+Wnm.*XHU(:,:,m);
                JHJmnx = z{m}*(Wnm.*XHU(:,:,n));
            end
            y(idxm) = y(idxm)+JHJmnx(:);
        end
        JHJnx = z{n}*Wn;
        y(idxn) = y(idxn)+JHJnx(:);
    end

    % If incomplete, approximate the effect of missing entries.
    if isstruct(T) && T.incomplete
        y = y*cache.scale;
    end

end

function x = M_blockJacobi(~,b)

    % Solve M*x = b, where M is a block-diagonal approximation for JHJ.
    % Equivalent to simultaneous ALS updates for each of the factors.
    x = nan(size(b));
    for n = 1:length(cache.offset)-1
        idx = cache.offset(n):cache.offset(n+1)-1;
        tmp = reshape(b(idx),[],size(cache.invW{1},1))*cache.invW{n};
        x(idx) = tmp(:);
    end

    % If incomplete, approximate the effect of missing entries.
    if isstruct(T) && T.incomplete
        x = x/cache.scale;
    end

end

function x = M_blockSSOR(z,b)

    % Solve Mx = b, where M is a block-Symmetric Successive Overrelaxation
    % preconditioner.
    % x = inv(D)*(U+D)*b
    B = cell(size(z));
    UHU = cache.UHU;
    BHU = nan(size(UHU));
    for n = 1:N
        B{n} = b(cache.offset(n):cache.offset(n+1)-1);
        BHU(:,:,n) = B{n}'*z{n};
    end
    X = B;
    for n = 1:N-1
        Wsum = zeros(R);
        for m = n+1:N
            Wnm = conj(reshape(prod(UHU([1:n-1 n+1:m-1 m+1:N],:),1),R,[]));
            Wsum = Wsum+Wnm.*conj(BHU(:,:,m));
        end
        Wn = conj(reshape(prod(UHU([1:n-1 n+1:N],:),1),[R R]));
        X{n} = X{n}+(z{n}*Wsum)/Wn;
    end
    % x = (L+D)*x
    B = X;
    for n = 1:N
        Wn = conj(reshape(prod(UHU([1:n-1 n+1:N],:),1),[R R]));
        BHU(:,:,n) = B{n}'*z{n};
        X{n} = B{n}*Wn;
    end
    for n = 2:N
        Wsum = zeros(R);
        for m = 1:n-1
            Wnm = conj(reshape(prod(UHU([1:n-1 n+1:m-1 m+1:N],:),1),R,[]));
            Wsum = Wsum+Wnm.*conj(BHU(:,:,m));
        end
        X{n} = X{n}+z{n}*Wsum;
    end
    x = cell2mat(cellfun(@(x)x(:),X(:),'UniformOutput',false));

    % If incomplete, approximate the effect of missing entries.
    if isstruct(T) && T.incomplete
        x = x/cache.scale;
    end

end

function [alpha,output] = ls(~,~,z,p,state,options)
    state.UHU = cache.UHU;
    state.T2 = cache.T2;
    [alpha,output] = linesearch(T,z,p,state,options);
end

function [alpha,output] = ps(~,~,z,p,q,state,options)
    state.UHU = cache.UHU;
    state.T2 = cache.T2;
    [alpha,output] = planesearch(T,z,p,q,state,options);
end

end

author: andrewssobral

repo: lrslibrary


text:
function[]=decode_beam(parameter,TestBatches,filename)
    disp('decode')
    parameter.beamSize=7;
    for batch_index=1:length(TestBatches)
        batch_index
        batch=TestBatches{batch_index};
        max_length=floor(batch.MaxLenSource*1);
        Word=batch.Word;
        N=size(Word,1);
        SourceLength=batch.SourceLength;
        [lstms,last_h_t,last_c_t]=Forward(batch,parameter,0);
        numElements =N*parameter.beamSize;
        translations=-zeroMatrix([N,max_length]);
        [first_scores,first_words]=BeamStep(parameter,last_h_t{parameter.layer_num},1);
        beamHistory=oneMatrix([numElements,max_length]);
        beamHistory(:,1)=first_words(:);
        beamScores=first_scores(:)';
        beamStates=cell(parameter.layer_num,1);
        for ll=1:parameter.layer_num
            beamStates{ll}.c_t=reshape(repmat(last_c_t{ll},parameter.beamSize,1),parameter.hidden,numElements);
            beamStates{ll}.h_t=reshape(repmat(last_h_t{ll},parameter.beamSize,1),parameter.hidden,numElements);
        end
        decoded=zeros(1,N);
        INDECE=oneMatrix([1,numElements]);
        for position=1:max_length
            words=beamHistory(:,position);
            for ll=1:parameter.layer_num
                if ll == 1
                    x_t=parameter.vect(:,words);
                else
                    x_t=beamStates{ll-1}.h_t;
                end
                h_t_1 = beamStates{ll}.h_t;
                c_t_1 = beamStates{ll}.c_t;
                [beamStates{ll}, h_t, c_t]=lstmUnit(parameter.W_T{ll},parameter,x_t,h_t_1, c_t_1, ll, -1,0);
                beamStates{ll}.h_t = h_t;
                beamStates{ll}.c_t = c_t;
            end
            [all_next_scores,all_next_words]=BeamStep(parameter,beamStates{parameter.layer_num}.h_t,0);
            all_next_scores=bsxfun(@plus,all_next_scores,beamScores);
            all_next_scores=reshape(all_next_scores,[parameter.beamSize*parameter.beamSize,N]);
            all_next_words=reshape(all_next_words,[parameter.beamSize*parameter.beamSize,N]);
            [all_next_scores,indices]=sort(all_next_scores,'descend');
            for sentId=1:N
                if decoded(sentId)==1 continue; end
                sorted_Indices = indices(:, sentId);
                sorted_next_words=all_next_words(sorted_Indices,sentId);
                end_index=find(sorted_next_words(1:parameter.beamSize)==parameter.stop,1);
                if sorted_next_words(1)==parameter.stop&&position>floor(SourceLength(sentId)/2)
                    decoded(sentId)=1;
                    end_index=sorted_Indices(1);
                    previous_index=floor((end_index-1)/parameter.beamSize)+1;
                    translations(sentId,1:position+1)=[beamHistory((sentId-1)*parameter.beamSize+previous_index,1:position),parameter.stop];
                    continue;
                end
                if position==max_length
                    end_index=sorted_Indices(1);
                    previous_index=floor((end_index-1)/parameter.beamSize)+1;
                    translations(sentId,1:position+1)=[beamHistory((sentId-1)*parameter.beamSize+previous_index,1:position),sorted_next_words(1)];
                    continue;
                end
                next_word_index=find(sorted_Indices~=parameter.stop,parameter.beamSize);
                beamScores(parameter.beamSize*(sentId-1)+1:parameter.beamSize*sentId)=all_next_scores(next_word_index,sentId);
                previous_index=(sentId-1)*parameter.beamSize+ floor((sorted_Indices(next_word_index)-1)/parameter.beamSize)+1;
                T=(sentId-1)*parameter.beamSize+1:(sentId-1)*parameter.beamSize+parameter.beamSize;
                beamHistory(T,:)=beamHistory(previous_index,:);
                beamHistory(T,position+1)=sorted_next_words(1:parameter.beamSize);
                INDECE(T)=previous_index;
            end
            if sum(decoded)==N
                break;
            end
            for ll=1:parameter.layer_num
                beamStates{ll}.c_t=beamStates{ll}.c_t(:,INDECE);
                beamStates{ll}.h_t=beamStates{ll}.h_t(:,INDECE);
            end
        end

        for senId=1:N
            vector=translations(senId,:);
            vector=vector(1:floor(1*batch.SourceLength(senId)));
            stop_sign=find(vector==parameter.stop);
            if length(stop_sign)==0
                dlmwrite(filename,vector,'delimiter',' ','-append');
            else
                dlmwrite(filename,vector(1:stop_sign-1),'delimiter',' ','-append');
            end
        end
    end
end

function[select_logP,select_words]=BeamStep(parameter,h_t,isFirst)
    if isFirst==1 scores=parameter.soft_W(2:end-1,:)*h_t;
    else scores=parameter.soft_W(2:end,:)*h_t;
    end
    mx = max(scores);
    scores = bsxfun(@minus, scores, mx);
    logP=bsxfun(@minus, scores, log(sum(exp(scores))));
    [sortedLogP, sortedWords]=sort(logP, 'descend');
    select_words=1+sortedWords(1:parameter.beamSize, :);
    select_logP=sortedLogP(1:parameter.beamSize, :);
    probs=exp(scores);
    norms = sum(probs, 1);
    probs=bsxfun(@rdivide, probs, norms);
end

author: jiweil

repo: Hierarchical-Neural-Autoencoder


text:
function d = xcorr1xr_orig(d)

   % This function performs cross correlations in sets by correlating one
   % trace against a row of traces at once. Generally it is faster than the
   % one at a time implimentation. All steps are included in this function.
   % That is, no calls to the Matlab built-in xcorr are used.

   % Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   % PREP NECESSARY TERMS
   [M,N] = size(d.w);
   pretrig = 86400*(d.trig-d.start);   % time between trace start and trigger
   l = (1./d.Fs)*[-M+1:M-1]';        % lag vector
   % next two lines are equivalent ways to determine normalization coefficients
   wcoeff = 1./sqrt(sum(d.w.*d.w));
   %for i = 1:size(d.w,2), wcoeff(i) = 1./norm(d.w(:,i)); end;
   d.corrmatrix = eye(length(d.trig),'single');
   d.lags = zeros(length(d.trig),'single');

   % GET FFT OF TRACES
   X = fft(d.w,2^nextpow2(2*M-1));
   Xc = conj(X);
   [MX,NX] = size(X);

   % TIME THE PROCESS
   starttime = now;
   tic;
   t0 = cputime;
   numcorr = ((length(d.trig))^2-length(d.trig))/2;
   count = 0;

   % LOOP THROUGH ROWS OF SIMILARITY MATRIX
   for n =1:N
      cols = n:N;
      % multiply fourier series and transform back to time domain
      CC = (X(:,n) * ones(1,length(cols))) .* Xc(:,cols);
      corr = ifft(CC);
      corr = [corr(end-M+2:end,:);corr(1:M,:)];

      % save only max correlation value
      [maxval,index] = max(corr);
      d.corrmatrix(n,cols)=maxval .* wcoeff(n) .* wcoeff(cols);   % normalized maximum correlation
      d.lags(n,cols)=l(index) + (pretrig(cols) - pretrig(n)); % lag in seconds

      if (n==100) | (n==1000)
         tclock = toc/86400;
         completed = sum(N-n:N-1)/((N*N-N)/2);
         fin = datestr(starttime + tclock/completed,16);
         if completed<0.5
            disp([num2str(100*completed,'%2.0f') '% completed. Estimated completion at ' fin '...']);
         end;
      end;
   end

   % FILL LOWER TRIANGULAR PART OF MATRICES
   d.corrmatrix = d.corrmatrix + d.corrmatrix' - eye(size(d.corrmatrix));
   d.lags = d.lags - d.lags';

   % DISPLAY RUN TIMES
   tclock = toc;
   tcpu = cputime-t0;
   if tclock>20
      disp(['Clock time to complete correlations: ' num2str(tclock,'%5.1f') ' s']);
      disp(['CPU time to complete correlations:   ' num2str(tcpu,'%5.1f') ' s']);
   end;
end


text:
 % Required for parsing the labels and numbers.
 files = dir('Data/*.jpg');
 for file = files
     [row,col] = size(file);
     for i=1:row
         % fprintf("File name %s", file(i).name);
         file_real_name = file(i).name;
         file_name = "Data/"+ file_real_name;
         label_name = "Labels/label_" + file_real_name + ".mat";
         label_name_num = "Labels/label_number_" + file_real_name + ".mat";
         img = imread(file_name{1});
         [labels, numlabels] = getSPLabels(img,200,1,150);
         save(label_name{1},'labels');
         save(label_name_num{1},'numlabels');
         labels = 0;
     end
 end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm

path: github-repos/MATLAB/zhangqianqianQQ-MachineVisionAlgorithm/MachineVisionAlgorithm-683338f6c3b1aab9fa2b80026915fe936aebf0ee/识别算法/Object Recognition based on super pixel/LabelParser.m


text:
close all;
clear all;
clc;

data_file_path = 'bin/ra_mmv_phase_transition_snr_30db_s_16.mat';
options.export = true;
options.export_dir = 'bin';
options.export_name = 'ra_mmv_snr_30_db_s_16';
options.chosen_ks = [2, 4, 8, 16, 32, 64];
options.subtitle = 'Rank Aware, SNR=30dB s=16';
spx.pursuit.PhaseTransitionAnalysis.print_results(data_file_path, ...
    'CoSaMP', options);

author: indigits

repo: sparse-plex


text:
% [INPUT]
% ds = A structure representing the dataset.

function analyze_dataset(varargin)

    persistent ip;

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

    ip.parse(varargin{:});

    ipr = ip.Results;
    ds = validate_dataset(ipr.ds);

    nargoutchk(0,0);

    analyze_dataset_internal(ds);

end

function analyze_dataset_internal(ds)

    safe_plot(@(id)plot_index(ds,id));

    safe_plot(@(id)plot_boxes(ds,'Returns',id));
    safe_plot(@(id)plot_sequence_returns(ds,id));

    if (~isempty(ds.Volumes))
        safe_plot(@(id)plot_boxes(ds,'Volumes',id));
        safe_plot(@(id)plot_sequence_other(ds,'Volumes',id));
    end

    if (~isempty(ds.Capitalizations))
        safe_plot(@(id)plot_boxes(ds,'Capitalizations',id));
        safe_plot(@(id)plot_sequence_other(ds,'Capitalizations',id));
    end

    if (~isempty(ds.CDS))
        safe_plot(@(id)plot_risk_free_rate(ds,id));
        safe_plot(@(id)plot_boxes(ds,'CDS',id));
        safe_plot(@(id)plot_sequence_other(ds,'CDS',id));
    end

    if (~isempty(ds.Assets) && ~isempty(ds.Equity))
        safe_plot(@(id)plot_boxes(ds,'Assets',id));
        safe_plot(@(id)plot_sequence_other(ds,'Assets',id));

        safe_plot(@(id)plot_boxes(ds,'Equity',id));
        safe_plot(@(id)plot_sequence_other(ds,'Equity',id));

        safe_plot(@(id)plot_boxes(ds,'Liabilities',id));
        safe_plot(@(id)plot_sequence_other(ds,'Liabilities',id));
    end

    if (ds.Crises > 0)
        safe_plot(@(id)plot_crises(ds,id));
    end

end

function plot_boxes(ds,target,id)

    n = ds.N;
    y = ds.(target);

    f = figure('Name',['Dataset > ' target ' (Box Plots)'],'Units','normalized','Position',[100 100 0.85 0.85],'Tag',id);

    boxplot(y,'Notch','on','Symbol','k.');
    set(findobj(f,'type','line','Tag','Median'),'Color','g');
    set(findobj(f,'-regexp','Tag','\w*Whisker'),'LineStyle','-');
    delete(findobj(f,'-regexp','Tag','\w*Outlier'));

    lower_av = findobj(f,'-regexp','Tag','Lower Adjacent Value');
    lower_av = cell2mat(get(lower_av,'YData'));
    y_low = min(lower_av(:));
    y_low = y_low - abs(y_low / 10);

    upper_av = findobj(f,'-regexp','Tag','Upper Adjacent Value');
    upper_av = cell2mat(get(upper_av,'YData'));
    y_high = max(upper_av(:));
    y_high = y_high + abs(y_high / 10);

    ax = gca();
    set(ax,'TickLength',[0 0]);
    set(ax,'XTick',1:n,'XTickLabels',ds.FirmNames,'XTickLabelRotation',45);
    set(ax,'YLim',[y_low y_high]);

    figure_title(f,target);

    maximize_figure(f);

end

function plot_crises(ds,id)

    f = figure('Name','Dataset > Crises','Units','normalized','Tag',id);

    if (strcmp(ds.CrisesType,'E'))
        cd = ds.CrisesDummy;
        cddn = ds.DatesNum(logical(cd));
        cddn_len = numel(cddn);

        plot(ds.DatesNum,nan(ds.T,1));

        hold on;
            for i = cddn_len:-1:1
                l = line(ones(2,1) .* cddn(i),[0 1],'Color',[1.000 0.400 0.400]);
                set(l,'Tag',num2str(i));
            end
        hold off;

        ax = gca();

        tooltips = ds.CrisisNames;
        tooltips_target = l;
    else
        cd = ds.CrisisDummies;
        k = size(cd,2);

        co = get(gca(),'ColorOrder');
        cor = ceil(k / size(co,1));
        co = repmat(co,cor,1);

        p = gobjects(k,1);

        hold on;
            for i = k:-1:1
                cddn = ds.DatesNum(logical(cd(:,i)));
                cddn_max = max(cddn);
                cddn_min = min(cddn);

                p(i) = patch('XData',[cddn_min cddn_max cddn_max cddn_min],'YData',[0 0 1 1],'EdgeAlpha',0.25,'FaceAlpha',0.40,'FaceColor',co(i,:));
                set(p(i),'Tag',num2str(i));
            end
        hold off;

        ax = gca();

        if (k <= 10)
            legend(ax,p,ds.CrisisNames,'Location','southoutside','Orientation','horizontal');

            tooltips = [];
            tooltips_target = [];
        else
            tooltips = ds.CrisisNames;
            tooltips_target = p;
        end
    end

    set(ax,'Box','on');
    set(ax,'XLim',[ds.DatesNum(1) ds.DatesNum(end)],'XTickLabelRotation',45);
    set(ax,'YLim',[0 1],'YTick',[]);

    if (ds.MonthlyTicks)
        date_ticks(ax,'x','mm/yyyy','KeepLimits','KeepTicks');
    else
        date_ticks(ax,'x','yyyy','KeepLimits');
    end

    if (strcmp(ds.CrisesType,'E'))
        figure_title(f,'Crises (Events)');
    else
        figure_title(f,'Crises (Ranges)');
    end

    maximize_figure(f);

    if (~isempty(tooltips))
        dcm = datacursormode(f);
        set(dcm,'Enable','on','SnapToDataVertex','off','UpdateFcn',@(targ,evtd)create_tooltip(targ,evtd,tooltips));

        if (verLessThan('MATLAB','9.7'))
            createDatatip(dcm,tooltips_target,[1 1]);
        else
            for i = 1:numel(tooltips_target)
                createDatatip(dcm,tooltips_target(i),[1 1]);
            end
        end
    end

    function tooltip = create_tooltip(~,evtd,tooltips)

        target = get(evtd,'Target');
        index = str2double(get(target,'Tag'));
        tooltip = tooltips{index};

    end

end

function plot_index(ds,id)

    y = ds.Index;

    yv = y(~isnan(y));
    y_obs = numel(yv);
    y_max = max(yv);
    y_min = min(yv);
    y_avg = mean(yv);
    y_med = median(yv);
    y_std = std(yv);
    y_ske = skewness(yv,0);
    y_kur = kurtosis(yv,0);

    r = yv - y_avg;
    [s_h,s_pval] = shapiro_test(r,0.05);
    [lbq_h,lbq_pval] = lbqtest(r.^2,'Alpha',0.05);
    [a_h,a_pval] = archtest(r,'Alpha',0.05);

    f = figure('Name','Dataset > Index','Units','normalized','Position',[100 100 0.85 0.85],'Tag',id);

    sub_1 = subplot(2,1,1);
    plot(sub_1,ds.DatesNum,ds.Index);
    set(sub_1,'XLim',[ds.DatesNum(1) ds.DatesNum(end)],'XTickLabelRotation',45);
    set(sub_1,'YLim',[(y_min - 0.01) (y_max + 0.01)]);
    set(sub_1,'XGrid','on','YGrid','on');
    t1 = title(sub_1,'Log Returns');
    set(t1,'Units','normalized');
    t1_position = get(t1,'Position');
    set(t1,'Position',[0.4783 t1_position(2) t1_position(3)]);

    if (ds.MonthlyTicks)
        date_ticks(sub_1,'x','mm/yyyy','KeepLimits','KeepTicks');
    else
        date_ticks(sub_1,'x','yyyy','KeepLimits');
    end

    sub_2 = subplot(2,1,2);
    hist = histogram(sub_2,ds.Index,50,'FaceColor',[0.749 0.862 0.933],'Normalization','pdf');
    edges = get(hist,'BinEdges');
    edges_max = max(edges);
    edges_min = min(edges);
    [hv,hp] = ksdensity(ds.Index);

    hold on;
        plot(sub_2,hp,hv,'-b','LineWidth',1.5);
    hold off;
    set(sub_2,'XLim',[(edges_min - (edges_min * 0.1)) (edges_max - (edges_max * 0.1))]);
    t2 = title(sub_2,'P&L Distribution');
    set(t2,'Units','normalized');
    t2_position = get(t2,'Position');
    set(t2,'Position',[0.4783 t2_position(2) t2_position(3)]);

    t = figure_title(f,['Index (' ds.IndexName ')']);
    t_position = get(t,'Position');
    set(t,'Position',[t_position(1) -0.0157 t_position(3)]);

    txt_obs = sprintf('Observations: %d',y_obs);
    txt_avg = sprintf('Mean: %.4f',y_avg);
    txt_med = sprintf('Median: %.4f',y_med);
    txt_std = sprintf('Standard Deviation: %.4f',y_std);
    txt_ske = sprintf('Skewness: %.4f',y_ske);
    txt_kur = sprintf('Kurtosis: %.4f',y_kur);

    if (s_h)
        txt_s = sprintf('Shapiro Test: T (%.4f)',s_pval);
    else
        txt_s = sprintf('Shapiro Test: F (%.4f)',s_pval);
    end

    if (lbq_h)
        txt_lbq = sprintf('LBQ Test: T (%.4f)',lbq_pval);
    else
        txt_lbq = sprintf('LBQ Test: F (%.4f)',lbq_pval);
    end

    if (a_h)
        txt_a = sprintf('ARCH Test: T (%.4f)',a_pval);
    else
        txt_a = sprintf('ARCH Test: F (%.4f)',a_pval);
    end

    txt = {txt_obs '' txt_avg txt_med txt_std txt_ske txt_kur '' txt_s txt_lbq txt_a};

    annotation('TextBox',(get(sub_2,'Position') + [0.01 -0.025 0 0]),'String',txt,'EdgeColor','none','FitBoxToText','on','FontSize',8);

    maximize_figure(f);

end

function plot_risk_free_rate(ds,id)

    rfr = ds.RiskFreeRate;
    y_limits_rfr = plot_limits(rfr,0.1);

    rfr_pc = [0; (((rfr(2:end) - rfr(1:end-1)) ./ rfr(1:end-1)) .* 100)];
    y_limits_rfr_pc = plot_limits(rfr_pc,0.1);

    f = figure('Name','Dataset > Risk-Free Rate','Units','normalized','Position',[100 100 0.85 0.85],'Tag',id);

    sub_1 = subplot(2,1,1);
    plot(sub_1,ds.DatesNum,smooth_data(rfr));
    set(sub_1,'XLim',[ds.DatesNum(1) ds.DatesNum(end)],'XTickLabelRotation',45);
    set(sub_1,'YLim',y_limits_rfr);
    set(sub_1,'XGrid','on','YGrid','on');
    t1 = title(sub_1,'Trend');
    set(t1,'Units','normalized');
    t1_position = get(t1,'Position');
    set(t1,'Position',[0.4783 t1_position(2) t1_position(3)]);

    sub_2 = subplot(2,1,2);
    plot(sub_2,ds.DatesNum,rfr_pc);
    set(sub_2,'XLim',[ds.DatesNum(1) ds.DatesNum(end)],'XTickLabelRotation',45);
    set(sub_2,'YLim',y_limits_rfr_pc);
    set(sub_2,'YTickLabels',arrayfun(@(x)sprintf('%.f%%',x),get(sub_2,'YTick'),'UniformOutput',false));
    set(sub_2,'XGrid','on','YGrid','on');
    t2 = title(sub_2,'Percent Change');
    set(t2,'Units','normalized');
    t2_position = get(t2,'Position');
    set(t2,'Position',[0.4783 t2_position(2) t2_position(3)]);

    if (ds.MonthlyTicks)
        date_ticks([sub_1 sub_2],'x','mm/yyyy','KeepLimits','KeepTicks');
    else
        date_ticks([sub_1 sub_2],'x','yyyy','KeepLimits');
    end

    figure_title(f,'Risk-Free Rate');

    maximize_figure(f);

end

function plot_sequence_returns(ds,id)

    n = ds.N;
    t = ds.T;
    dn = ds.DatesNum;
    mt = ds.MonthlyTicks;

    ts = ds.Returns;
    data = [repmat({dn},1,n); mat2cell(ts,t,ones(1,n))];

    plots_title = [repmat({'Log Returns'},1,n); repmat({'P&L Distribution'},1,n)];

    x_limits = [dn(1) dn(end)];

    core = struct();

    core.N = n;
    core.Data = data;
    core.Function = @(subs,data)plot_function(subs,data);

    core.OuterTitle = 'Dataset > Returns (Time Series)';
    core.InnerTitle = 'Returns (Time Series)';
    core.SequenceTitles = ds.FirmNames;

    core.PlotsAllocation = [2 1];
    core.PlotsSpan = {1 2};
    core.PlotsTitle = plots_title;

    core.XDates = {mt []};
    core.XGrid = {true false};
    core.XLabel = {[] []};
    core.XLimits = {x_limits []};
    core.XRotation = {45 []};
    core.XTick = {[] []};
    core.XTickLabels = {[] []};

    core.YGrid = {true false};
    core.YLabel = {[] []};
    core.YLimits = {[] []};
    core.YRotation = {[] []};
    core.YTick = {[] []};
    core.YTickLabels = {[] []};

    sequential_plot(core,id);

    function plot_function(subs,data)

        x = data{1};
        y = data{2};

        yv = y(~isnan(y));
        y_obs = numel(yv);
        y_max = max(yv);
        y_min = min(yv);
        y_avg = mean(yv);
        y_med = median(yv);
        y_std = std(yv);
        y_ske = skewness(yv,0);
        y_kur = kurtosis(yv,0);

        r = yv - y_avg;
        [s_h,s_pval] = shapiro_test(r,0.05);
        [lbq_h,lbq_pval] = lbqtest(r.^2,'Alpha',0.05);
        [a_h,a_pval] = archtest(r,'Alpha',0.05);

        d = find(isnan(y),1,'first');

        if (isempty(d))
            xd = [];
        else
            xd = x(d) - 1;
        end

        delete(findall(gcf,'type','annotation'));

		sub_1 = subs(1);
        plot(sub_1,x,y,'Color',[0.000 0.447 0.741]);
        set(sub_1,'YLim',[(y_min - 0.01) (y_max + 0.01)]);

        if (~isempty(xd))
            hold(sub_1,'on');
                plot(sub_1,[xd xd],get(sub_1,'YLim'),'Color',[1.000 0.400 0.400]);
            hold(sub_1,'off');
        end

		sub_2 = subs(2);
        hist = histogram(sub_2,y,50,'FaceColor',[0.749 0.862 0.933],'Normalization','pdf');
        edges = get(hist,'BinEdges');
        edges_max = max(edges);
        edges_min = min(edges);
        [dv,dp] = ksdensity(y);

        hold(sub_2,'on');
            plot(sub_2,dp,dv,'-b','LineWidth',1.5);
        hold(sub_2,'off');
        set(sub_2,'XLim',[(edges_min - (edges_min * 0.1)) (edges_max - (edges_max * 0.1))]);

        txt_obs = sprintf('Observations: %d',y_obs);
        txt_avg = sprintf('Mean: %.4f',y_avg);
        txt_med = sprintf('Median: %.4f',y_med);
        txt_std = sprintf('Standard Deviation: %.4f',y_std);
        txt_ske = sprintf('Skewness: %.4f',y_ske);
        txt_kur = sprintf('Kurtosis: %.4f',y_kur);

        if (s_h)
            txt_s = sprintf('Shapiro Test: T (%.4f)',s_pval);
        else
            txt_s = sprintf('Shapiro Test: F (%.4f)',s_pval);
        end

        if (lbq_h)
            txt_lbq = sprintf('LBQ Test: T (%.4f)',lbq_pval);
        else
            txt_lbq = sprintf('LBQ Test: F (%.4f)',lbq_pval);
        end

        if (a_h)
            txt_a = sprintf('ARCH Test: T (%.4f)',a_pval);
        else
            txt_a = sprintf('ARCH Test: F (%.4f)',a_pval);
        end

        txt = {txt_obs '' txt_avg txt_med txt_std txt_ske txt_kur '' txt_s txt_lbq txt_a};

        annotation('TextBox',(get(sub_2,'Position') + [0.01 -0.025 0 0]),'String',txt,'EdgeColor','none','FitBoxToText','on','FontSize',8);

    end

end

function plot_sequence_other(ds,target,id)

    n = ds.N;
    t = ds.T;
    dn = ds.DatesNum;
    mt = ds.MonthlyTicks;

    ts = smooth_data(ds.(target));
    data = [repmat({dn},1,n); mat2cell(ts,t,ones(1,n))];

    plots_title = repmat({' '},1,n);

    x_limits = [dn(1) dn(end)];

    core = struct();

    core.N = n;
    core.Data = data;
    core.Function = @(subs,data)plot_function(subs,data);

    core.OuterTitle = ['Dataset > ' target ' (Time Series)'];
    core.InnerTitle = [target ' (Time Series)'];
    core.SequenceTitles = ds.FirmNames;

    core.PlotsAllocation = [1 1];
    core.PlotsSpan = {1};
    core.PlotsTitle = plots_title;

    core.XDates = {mt};
    core.XGrid = {true};
    core.XLabel = {[]};
    core.XLimits = {x_limits};
    core.XRotation = {45};
    core.XTick = {[]};
    core.XTickLabels = {[]};

    core.YGrid = {true};
    core.YLabel = {[]};
    core.YLimits = {[]};
    core.YRotation = {[]};
    core.YTick = {[]};
    core.YTickLabels = {[]};

    sequential_plot(core,id);

    function plot_function(subs,data)

        x = data{1};
        y = data{2};

        d = find(isnan(y),1,'first');

        if (isempty(d))
            xd = [];
        else
            xd = x(d) - 1;
        end

		sub_1 = subs(1);
        plot(sub_1,x,y,'Color',[0.000 0.447 0.741]);

        if (~isempty(xd))
            hold(sub_1,'on');
                plot(sub_1,[xd xd],get(sub_1,'YLim'),'Color',[1.000 0.400 0.400]);
            hold(sub_1,'off');
        end

    end

end

author: TommasoBelluzzo

repo: SystemicRisk


text:
function [f]=comp_idwiltiii(coef,g)
%COMP_IDWILTIII  Compute Inverse discrete Wilson transform type III.
%
%   This is a computational routine. Do not call it
%   directly.

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

M=size(coef,1);
N=size(coef,2);
W=size(coef,3);
a=M;

L=N*M;

coef2=zeros(2*M,N,W,assert_classname(coef,g));

coef2(1:2:M,1:2:N,:)        = exp( i*pi/4)*coef(1:2:M,1:2:N,:);
coef2(2*M:-2:M+1,1:2:N,:)   = exp(-i*pi/4)*coef(1:2:M,1:2:N,:);

coef2(1:2:M,2:2:N,:)        = exp(-i*pi/4)*coef(1:2:M,2:2:N,:);
coef2(2*M:-2:M+1,2:2:N,:)   = exp( i*pi/4)*coef(1:2:M,2:2:N,:);

coef2(2:2:M,1:2:N,:)        = exp(-i*pi/4)*coef(2:2:M,1:2:N,:);
coef2(2*M-1:-2:M+1,1:2:N,:) = exp( i*pi/4)*coef(2:2:M,1:2:N,:);

coef2(2:2:M,2:2:N,:)        = exp( i*pi/4)*coef(2:2:M,2:2:N,:);
coef2(2*M-1:-2:M+1,2:2:N,:) = exp(-i*pi/4)*coef(2:2:M,2:2:N,:);

% Apply the generalized DGT and scale.
%f=comp_igdgt(coef2,g,a,2*M,L,0,.5,0,0)/sqrt(2);

f = comp_isepdgt(coef2,g,L,a,2*M,0);

halfmod=exp(pi*i*(0:L-1).'/(2*M))/sqrt(2);
f=f.*repmat(halfmod,1,W);

if isreal(coef) && isreal(g)
  f=real(f);
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function display(om, varargin)
%DISPLAY  Displays the object.
%   Called when semicolon is omitted at the command-line. Displays the details
%   of the variables, constraints, costs included in the model.
%
%   See also OPT_MODEL.

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

%% Due to a bug related to inheritance in constructors in
%% Octave 5.2 and earlier (https://savannah.gnu.org/bugs/?52614),
%% INIT_SET_TYPES() cannot be called directly in the
%% MP_IDX_MANAGER constructor, as desired.
%%
%% WORKAROUND:  Initialize MP_IDX_MANAGER fields here, if needed,
%%              after object construction, but before object use.
if isempty(om.var)          %% only if not already initialized
    om.init_set_types();
end

if nargin < 2
    more_set_types = {};
else
    more_set_types = varargin{1};
end

%% display details of each set type
set_types = {'var', 'nle', 'nli', 'lin', 'qdc', 'nlc', more_set_types{:}};
fprintf('\n');
for k = 1:length(set_types)
    om.display_set(set_types{k});
end

%% user data
fields = fieldnames(om.userdata);
if ~isempty(fields)
    fprintf('\nUSER DATA\n')
    fprintf('=========\n')
    fprintf('  name                               size       class\n');
    fprintf(' ------------------------------   -----------  --------------------\n');
    for k = 1:length(fields)
        f = om.userdata.(fields{k});
        [m, n] = size(f);
        fprintf('  %-31s %5dx%-5d   %s\n', fields{k}, m, n, class(f));
    end
else
    fprintf('USER DATA                   :  <none>\n');
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function n = triangle_right_lattice_point_num_2d ( a, b )

%*****************************************************************************80
%
%% TRIANGLE_RIGHT_LATTICE_POINT_NUM_2D: count lattice points.
%
%  Discussion:
%
%    The triangle is assumed to be a right triangle which, without loss
%    of generality, has the coordinates:
%
%    ( (0,0), (a,0), (0,b) )
%
%    The routine returns the number of integer lattice points that appear
%    inside the triangle or on its edges or vertices.
%
%    The formula for this function occurred to me (JVB) after some thought,
%    on 06 July 2009.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 July 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer A, B, define the vertices.
%
%    Output, integer N, the number of lattice points.
%
  n = ( ( a + 1 ) * ( b + 1 ) + i4_gcd ( a, b ) + 1 ) / 2;

  return
end

author: johannesgerer


text:
function [fc_2,Rc_2,Tc_2,H_2,distance,V_vert,V_hori,x_all_c,V_hori_pix,V_vert_pix,V_diag1_pix,V_diag2_pix]=Distor2Calib(k_dist,grid_pts_centered,n_sq_x,n_sq_y,Np,W,L,Xgrid_2,f_ini,N_iter,two_focal);

% Computes the calibration parameters knowing the
% distortion factor k_dist

% grid_pts_centered are the grid point coordinates after substraction of
% the optical center.

% can give an optional guess for the focal length f_ini (can set to [])
% can provide the number of iterations for the Iterative Vanishing Point Algorithm

% if the focal length is known perfectly, then, there is no need to iterate,
% and therefore, one can fix: N_iter = 0;

% California Institute of Technology
% (c) Jean-Yves Bouguet - October 7th, 1997

%keyboard;

if exist('two_focal'),
   if isempty(two_focal),
      two_focal=0;
   end;
else
   two_focal = 0;
end;

if exist('N_iter'),
   if ~isempty(N_iter),
      disp('Use number of iterations provided');
   else
      N_iter = 10;
   end;
else
   N_iter = 10;
end;

if exist('f_ini'),
   if ~isempty(f_ini),
      disp('Use focal provided');
      if length(f_ini)<2, f_ini=[f_ini;f_ini]; end;
      fc_2 = f_ini;
      x_all_c = [grid_pts_centered(1,:)/fc_2(1);grid_pts_centered(2,:)/fc_2(2)];
      x_all_c = comp_distortion(x_all_c,k_dist); % we can this time!!!
   else
     fc_2 = [1;1];
     x_all_c = grid_pts_centered;
   end;
else
   fc_2 = [1;1];
   x_all_c = grid_pts_centered;
end;

dX = W/n_sq_x;
dY = L/n_sq_y;

N_x = n_sq_x+1;
N_y = n_sq_y+1;

x_grid = zeros(N_x,N_y);
y_grid = zeros(N_x,N_y);

%%% Computation of the four vanishing points in pixels

   x_grid(:) = grid_pts_centered(1,:);
   y_grid(:) = grid_pts_centered(2,:);

   for k=1:n_sq_x+1,
      [U,S,V] = svd([x_grid(k,:);y_grid(k,:);ones(1,n_sq_y+1)]);
      vert(:,k) = U(:,3);
   end;

   for k=1:n_sq_y+1,
      [U,S,V] = svd([x_grid(:,k)';y_grid(:,k)';ones(1,n_sq_x+1)]);
      hori(:,k) = U(:,3);
   end;

   % 2 principle Vanishing points:
   [U,S,V] = svd(vert);
   V_vert = U(:,3);
   [U,S,V] = svd(hori);
   V_hori = U(:,3);

   % Square warping:

   vert_first = vert(:,1) - dot(V_vert,vert(:,1))/dot(V_vert,V_vert) * V_vert;
   vert_last = vert(:,n_sq_x+1) - dot(V_vert,vert(:,n_sq_x+1))/dot(V_vert,V_vert) * V_vert;

   hori_first = hori(:,1) - dot(V_hori,hori(:,1))/dot(V_hori,V_hori) * V_hori;
   hori_last = hori(:,n_sq_y+1) - dot(V_hori,hori(:,n_sq_y+1))/dot(V_hori,V_hori) * V_hori;

   x1 = cross(hori_first,vert_first);
   x2 = cross(hori_first,vert_last);
   x3 = cross(hori_last,vert_last);
   x4 = cross(hori_last,vert_first);

   x1 = x1/x1(3);
   x2 = x2/x2(3);
   x3 = x3/x3(3);
   x4 = x4/x4(3);

   [square] = Rectangle2Square([x1 x2 x3 x4],W,L);

   y1 = square(:,1);
   y2 = square(:,2);
   y3 = square(:,3);
   y4 = square(:,4);

   H2 = cross(V_vert,V_hori);

   V_diag1 = cross(cross(y1,y3),H2);
   V_diag2 = cross(cross(y2,y4),H2);

   V_diag1 = V_diag1 / norm(V_diag1);
   V_diag2 = V_diag2 / norm(V_diag2);

   V_hori_pix = V_hori;
   V_vert_pix = V_vert;
   V_diag1_pix = V_diag1;
   V_diag2_pix = V_diag2;

% end of computation of the vanishing points in pixels.

if two_focal, % only if we attempt to estimate two focals...
   % Use diagonal lines also to add two extra vanishing points (?)
   N_min = min(N_x,N_y);

   if N_min < 2,
      use_diag = 0;
      two_focal = 0;
      disp('Cannot estimate two focals (no existing diagonals)');
   else
      use_diag = 1;
      Delta_N = abs(N_x-N_y);
      N_extra = round((N_min - Delta_N - 1)/2);
      diag_list = -N_extra:Delta_N+N_extra;
      N_diag = length(diag_list);
      diag_1 = zeros(3,N_diag);
      diag_2 = zeros(3,N_diag);
   end;
else
   % Give up the use of the diagonals (so far)
   % it seems that the error is increased
   use_diag = 0;
end;

% The vertical lines: vert, Horizontal lines: hori
vert = zeros(3,n_sq_x+1);
hori = zeros(3,n_sq_y+1);

for counter_k = 1:N_iter, 	% the Iterative Vanishing Points Algorithm to
                                % estimate the focal length accurately

   x_grid(:) = x_all_c(1,:);
   y_grid(:) = x_all_c(2,:);

   for k=1:n_sq_x+1,
      [U,S,V] = svd([x_grid(k,:);y_grid(k,:);ones(1,n_sq_y+1)]);
      vert(:,k) = U(:,3);
   end;

   for k=1:n_sq_y+1,
      [U,S,V] = svd([x_grid(:,k)';y_grid(:,k)';ones(1,n_sq_x+1)]);
      hori(:,k) = U(:,3);
   end;

   % 2 principle Vanishing points:
   [U,S,V] = svd(vert);
   V_vert = U(:,3);
   [U,S,V] = svd(hori);
   V_hori = U(:,3);

   % Square warping:

   vert_first = vert(:,1) - dot(V_vert,vert(:,1))/dot(V_vert,V_vert) * V_vert;
   vert_last = vert(:,n_sq_x+1) - dot(V_vert,vert(:,n_sq_x+1))/dot(V_vert,V_vert) * V_vert;

   hori_first = hori(:,1) - dot(V_hori,hori(:,1))/dot(V_hori,V_hori) * V_hori;
   hori_last = hori(:,n_sq_y+1) - dot(V_hori,hori(:,n_sq_y+1))/dot(V_hori,V_hori) * V_hori;

   x1 = cross(hori_first,vert_first);
   x2 = cross(hori_first,vert_last);
   x3 = cross(hori_last,vert_last);
   x4 = cross(hori_last,vert_first);

   x1 = x1/x1(3);
   x2 = x2/x2(3);
   x3 = x3/x3(3);
   x4 = x4/x4(3);

   [square] = Rectangle2Square([x1 x2 x3 x4],W,L);

   y1 = square(:,1);
   y2 = square(:,2);
   y3 = square(:,3);
   y4 = square(:,4);

   H2 = cross(V_vert,V_hori);

   V_diag1 = cross(cross(y1,y3),H2);
   V_diag2 = cross(cross(y2,y4),H2);

   V_diag1 = V_diag1 / norm(V_diag1);
   V_diag2 = V_diag2 / norm(V_diag2);

   % Estimation of the focal length, and normalization:

   % Compute the ellipsis of (1/f^2) positions:
   % a * (1/fx)^2 + b * (1/fx)^2 = -c

   a1 = V_hori(1);
   b1 = V_hori(2);
   c1 = V_hori(3);

   a2 = V_vert(1);
   b2 = V_vert(2);
   c2 = V_vert(3);

   a3 = V_diag1(1);
   b3 = V_diag1(2);
   c3 = V_diag1(3);

   a4 = V_diag2(1);
   b4 = V_diag2(2);
   c4 = V_diag2(3);

   if two_focal,

      A = [a1*a2 b1*b2;a3*a4 b3*b4];
      b = -[c1*c2;c3*c4];

      f = sqrt(abs(1./(inv(A)*b)));

   else

      f = sqrt(abs(-(c1*c2*(a1*a2 + b1*b2) + c3*c4*(a3*a4 + b3*b4))/(c1^2*c2^2 + c3^2*c4^2)));

      f = [f;f];

   end;

   % REMARK:
   % if both a and b are small, the calibration is impossible.
   % if one of them is small, only the other focal length is observable
   % if none is small, both focals are observable

   fc_2 = fc_2 .* f;

   % DEBUG PART: fix focal to 500...
   %fc_2= [500;500]; disp('Line 293 to be earased in Distor2Calib.m');

   % end of focal compensation

   % normalize by the current focal:

   x_all = [grid_pts_centered(1,:)/fc_2(1);grid_pts_centered(2,:)/fc_2(2)];

   % Compensate by the distortion factor:

   x_all_c = comp_distortion(x_all,k_dist);

end;

% At that point, we hope that the distortion is gone...

x_grid(:) = x_all_c(1,:);
y_grid(:) = x_all_c(2,:);

for k=1:n_sq_x+1,
   [U,S,V] = svd([x_grid(k,:);y_grid(k,:);ones(1,n_sq_y+1)]);
   vert(:,k) = U(:,3);
end;

for k=1:n_sq_y+1,
   [U,S,V] = svd([x_grid(:,k)';y_grid(:,k)';ones(1,n_sq_x+1)]);
   hori(:,k) = U(:,3);
end;

% Vanishing points:
[U,S,V] = svd(vert);
V_vert = U(:,3);
[U,S,V] = svd(hori);
V_hori = U(:,3);

% Horizon:

H_2 = cross(V_vert,V_hori);

%   H_2 = cross(V_vert,V_hori);

% pick a plane in front of the camera (positive depth)
if H_2(3) < 0, H_2 = -H_2; end;

% Rotation matrix:

if V_hori(1) < 0, V_hori = -V_hori; end;

V_hori = V_hori/norm(V_hori);
H_2 = H_2/norm(H_2);

V_hori = V_hori - dot(V_hori,H_2)*H_2;

Rc_2 = [V_hori cross(H_2,V_hori) H_2];

Rc_2 = Rc_2 / det(Rc_2);

%omc_2 = rodrigues(Rc_2);

%Rc_2 = rodrigues(omc_2);

% Find the distance of the plane for translation vector:

xc_2 = [x_all_c;ones(1,Np)];

Zc_2 = 1./sum(xc_2 .* (Rc_2(:,3)*ones(1,Np)));

Xo_2 = [sum(xc_2 .* (Rc_2(:,1)*ones(1,Np))).*Zc_2 ; sum(xc_2 .* (Rc_2(:,2)*ones(1,Np))).*Zc_2];

XXo_2 = Xo_2 - mean(Xo_2')'*ones(1,Np);

distance_x = norm(Xgrid_2(1,:))/norm(XXo_2(1,:));
distance_y = norm(Xgrid_2(2,:))/norm(XXo_2(2,:));

distance = sum(sum(XXo_2(1:2,:).*Xgrid_2(1:2,:)))/sum(sum(XXo_2(1:2,:).^2));

alpha = abs(distance_x - distance_y)/distance;

if (alpha>0.1)&~two_focal,
   disp('Should use two focals in x and y...');
end;

% Deduce the translation vector:

Tc_2 = distance * H_2;

return;

   V_hori_pix/V_hori_pix(3)
   V_vert_pix/V_vert_pix(3)
   V_diag1_pix/V_diag1_pix(3)
   V_diag2_pix/V_diag2_pix(3)

repo: ekfmonoslam