text:
function sphere_monte_carlo_test01 ( )

%*****************************************************************************80
%
%% SPHERE_MONTE_CARLO_TEST01 tests SPHERE01_SAMPLE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 January 2014
%
%  Author:
%
%    John Burkardt
%
  e_test(1:3,1:7) = [ ...
    0, 0, 0; ...
    2, 0, 0; ...
    0, 2, 0; ...
    0, 0, 2; ...
    4, 0, 0; ...
    2, 2, 0; ...
    0, 0, 4 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01\n' );
  fprintf ( 1, '  Use SPHERE01_SAMPLE to estimate integrals over \n' );
  fprintf ( 1, '  the surface of the unit sphere.\n' );

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '         N        1              X^2             Y^2' )
  fprintf ( 1, '             Z^2             X^4           X^2Y^2           Z^4\n' );
  fprintf ( 1, '\n' );

  n = 1;

  while ( n <= 65536 )

    [ x, seed ] = sphere01_sample ( n, seed );

    fprintf ( 1, '  %8d', n );

    for j = 1 : 7

      e(1:3,1) = e_test(1:3,j);

      value(1:n,1) = monomial_value ( 3, n, e, x );

      result = sphere01_area ( ) * sum ( value(1:n) ) / n;

    fprintf ( 1, '  %14f', result );

    end

    fprintf ( 1, '\n' );

    n = 2 * n;

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '     Exact' );
  for j = 1 : 7

    e(1:3,1) = e_test(1:3,j);

    result = sphere01_monomial_integral ( e );

    fprintf ( 1, '  %14f', result );

  end
  fprintf ( 1, '\n' );

  return
end


text:
function cortMag = getCortMagParams(cortMag)
%
%  cortmag = getCortMagParams(cortMag)
%
%  Set up the initial cortMag structure, asking the user for
% information about different parameters. Returns 0 if the user
% cancels, an updated cortMag struct oitherwise.
%
% HISTORY:
%   2002.03.05 RFD (bob@white.stanford.edu) wrote it, based
%   on 'CortMagUI' by Wandell.
%   2002.12.11 RFD: added flatDataFlag.

prompt = {...
    'Bin distance      (mm):', ...
    'Ring Exp        (scan):',...
    'Stim Radius      (deg):',...
    'Co Thresh      (0-1.0):', ...
    'Stim Foveal Ph [0,2pi]:', ...
    'Stim Periph Ph [0,2pi]:', ...
    'Exclude ROIs:', ...
    'Shift Template #:',...
    'Use data from flat view:'};

% Default values
def={'4','2','20','0.20','4.71','4.71','','1','0'};

% If CortMag exists, use the existing values as the defaults
%
if isfield(cortMag,'binDist'),        def{1} = num2str(cortMag.binDist); end
if isfield(cortMag,'expNumber'),      def{2} = num2str(cortMag.expNumber); end
if isfield(cortMag,'stimulusRadius'), def{3} = num2str(cortMag.stimulusRadius); end
if isfield(cortMag,'coThresh'),       def{4} = num2str(cortMag.coThresh); end
if isfield(cortMag,'fovealPhase'),    def{5} = num2str(cortMag.fovealPhase); end
if isfield(cortMag,'peripheralPhase'),def{6} = num2str(cortMag.peripheralPhase); end
if isfield(cortMag,'excludeROIs'),    def{7} = num2str(cortMag.excludeROIs); end
if isfield(cortMag,'templateRoiNum'), def{8} = num2str(cortMag.templateRoiNum); end
if isfield(cortMag,'flatDataFlag'),   def{9} = num2str(cortMag.flatDataFlag); end

answer = inputdlg(prompt, 'CortMag Parameters', 1, def, 'on');

if ~isempty(answer)
   cortMag.binDist         = str2num(answer{1});
   cortMag.expNumber       = str2num(answer{2});
   cortMag.stimulusRadius  = str2num(answer{3});
   cortMag.coThresh        = str2num(answer{4});
   cortMag.fovealPhase     = str2num(answer{5});
   cortMag.peripheralPhase = str2num(answer{6});
   cortMag.excludeROIs     = str2num(answer{7});
   cortMag.templateRoiNum  = str2num(answer{8});
   cortMag.flatDataFlag    = str2num(answer{9});
else
   cortMag = 0;
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef TestTimelapser
    %TestTimelapser

    methods (Static)
        function test_1
            ttypes = {'AsIs', 'Crop'};
            for i=1:numel(ttypes)
                obj = cv.Timelapser(ttypes{i});
                typename = obj.typeid();
                validateattributes(typename, {'char'}, {'row', 'nonempty'});

                %{
                %TODO
                obj.initialize(corners, sizes);
                obj.process(img, mask, tl);
                dst = obj.getDst();
                %}
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [logCr,logr]=gencorint(x,dim,tau,logr,p,w,svd,q)
%Syntax: [logCr,logr]=gencorint(x,dim,tau,logr,p,w,svd,q)
%________________________________________________________
%
% Calculates the generalized Correlation Integral (Cr) of a time
% series x.
%
% logCr is the the value of log(Cr).
% logr is log(range).
% x is the time series.
% dim is the embedding dimension.
% tau is the time delay.
% p is defines the norm.
% w is the Theiler's correction.
% svd is the number of singular values taken into account.
% q is the generalization index.
%
%
% References:
%
% Grassberger P, Procaccia I (1983): Characterization of strange
% attractors. Physical Review Letters 50(5):346-349
%
% Theiler J (1986): Spurious dimension from correlation algorithms
% applied tolimited time-series data. Physical Review A 34(3):2427-
% 2432
%
% Albano A M, Muench J, Schwartz C, Mees A I, Rapp P E, (1988):
% Singular-value decomposition and the Grassberger-Procaccia algorithm.
% Physical Review A38:3017-3026
%
%
% Alexandros Leontitsis
% Department of Education
% University of Ioannina
% 45110 - Dourouti
% Ioannina
% Greece
%
% University e-mail: me00743@cc.uoi.gr
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% June 15, 2001.

if nargin<1 | isempty(x)==1
   error('You should provide a time series.');
else
   % x must be a vector
   if min(size(x))>1
      error('Invalid time series.');
   end
   x=x(:);
   % n is the time series length
   n=length(x);
end

if nargin<2 | isempty(dim)==1
   dim=2;
else
   % dim must be either a scalar or a vector
   if min(size(dim))>1
      error('dim must be a scalar or a vector.');
   end
   % dim must be an integer
   dim=round(dim);
   % dim values must be above 0
   dim=dim(find(dim>0));
end

if nargin<3 | isempty(tau)==1
   tau=1;
else
   % tau must be either a scalar or a vector
   if min(size(tau))>1
      error('tau must be a scalar or a vector.');
   end
   % tau must be an integer
   tau=round(tau);
   % tau values must be above 0
   tau=tau(find(tau>0));
end

if nargin<4  | isempty(logr)==1
   r=(max(x)-min(x))/10*(1:20)'/20;
   logr=log10(r);
else
    % logr must be either a scalar or a vector
    if min(size(dim))>1
        error('logr must be a scalar or a vector.');
    end
    % if it is a scalar, it determines the maximum range
    if length(logr)==1
        div=logr;
        % div must be positive
        if div<=0
            error('logr must be a positive scalar or vector');
        end
        r=(max(x)-min(x))/div*(1:20)'/20;
        logr=log10(r);
    else
        logr=logr(:);
        r=10.^logr;
    end
end

if nargin<5 | isempty(p)==1
   p=inf;
else
   % p must be either a scalar or a vector
   if min(size(dim))>1
      error('p must be a scalar.');
   end
   % p values must be positive
   p=p(find(p>0));
end

if nargin<6  | isempty(w)==1
   w=1;
else
   % w must be either a scalar or a vector
   if min(size(w))>1
      error('w must be either a scalar or a vector.');
   end
   % w must be an integer
   w=round(w);
   % w must be positive
   w=w(find(w>0));
end

if nargin<7 | isempty(svd)==1
    svd=[];
else
    % svd must be a scalar or a vector
    if min(size(svd))>1
        error('svd must be either a scalar or a vector.');
    end
    % svd must be an integer
    svd=round(svd);
    % svd must be positive
    svd=svd(find(svd>0));
end

if nargin<8 | isempty(q)==1
    q=2;
else
    % q must be either a scalar or a vector
    if min(size(q))>2
        error('q must be either a scalar or a vector.');
    end
end

% Only one of dim, tau, p, w, or q should be vector
l=[length(dim),length(tau),length(p),length(w),length(svd),length(q)];
if length(find(l>1))>1
   error('Only one of dim, tau, p, w, svd, or q should be vector.');
end

m=max(l);
dim=ones(1,m).*dim;
tau=ones(1,m).*tau;
p=ones(1,m).*p;
w=ones(1,m).*w;
if isempty(svd)==0
   svd=ones(1,m).*svd;
end
q=ones(1,m).*q;

for i=1:m

   % Reconstruct the time-delay phase-space
   [Y,T]=phasespace(x,dim(i),tau(i));
   if isempty(svd)==0
      svd(i)=min(svd(i),dim(i));
      % SVD on X
      [u,s,v]=svd(Y,0);

      % Calculate the Principal Components
      pc=Y*v;

      % Reconstruct the first svd Principal Components
      Y=pc(:,1:svd(i))*v(:,1:svd(i))';
   end

   % Initialize the logCr
   Cr=zeros(size(r));

   if q(i)==2 % ...fast

       for i1=1:T-w(i)
           for i2=i1+w(i):T
               dist=norm(Y(i1,:)-Y(i2,:),p(i));
               s=find(r>dist);
               Cr(s)=Cr(s)+1;
           end
       end
       Cr=2*Cr/T/(T-1);
       logCr(:,i)=log10(Cr);

   else % slow...

       for i1=1:T
           c1=zeros(size(r));
           for i2=1:T
               if i1<i2-(w(i)-1) | i1>i2+(w(i)-1)
                   dist=norm(Y(i1,:)-Y(i2,:),p(i));
                   s=find(r>dist);
                   c1(s)=c1(s)+1;
               end
           end
           c1=c1/(T-1);
           if q(i)~=1
               Cr=Cr+c1.^(q(i)-1);
           else
               Cr=Cr+c1;
           end
       end
       if q(i)~=1
           Cr=Cr/T;
           logCr(:,i)=log10(Cr.^(1/(q(i)-1)));
       else
           Cr=log10(Cr)/T;
           logCr(:,i)=Cr;
       end

   end
end

author: Sable

repo: mcbench-benchmarks


text:
function CPerm=permuteMatrix(C,order)
%%PERMUTEMATRIX Given an n1Xn2X...nS matrix, rearrange the dimensions of C
%      to be in the order specified by order. All elements of order must be
%      unique, real, positive, integer values from 1 to S. This functions
%      in the same manner as Matlab's permute function. The permutation
%      performed by this function is not done in-place. This type of
%      permutation is a type of tensor matrix transpose.
%
%INPUTS: C An n1Xn2X...XnS matrix.
%    order The length S order vector. This should hold integer values from
%          1 to S with no repeats.
%
%OUTPUTS: CPerm C with its dimensions permuted according to order.
%
%The algorithm is essentially brute force. It goes through each element in
%C and find the tuple of the associated element in CPerm to which the
%assignment must be performed.
%
%EXAMPLE:
% C=100*rand(3,4,5);
% order=[3,2,1];
% CPerm=permuteMatrix(C,order);
%One will ntoe that size(CPerm) is note [5,4,3].
%
%March 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

S=ndims(C);
numOrderEls=length(order);

if(isempty(order)&&~isempty(C))
   error('order cannot be empty with a non-empty C.')
end

if(numOrderEls<S)
    error('The ordering vector must be >=ndims(C).')
end

if(any(order>numOrderEls)||any(order<1))
    error('')
end

if(isempty(C))
    CPerm=[];
    return;
end

nVals=size(C);

if(numOrderEls>S)
   nVals=[nVals,ones(1,numOrderEls-S)];
   S=numOrderEls;
end

totalEls=numel(C);

cumProdNew=zeros(S,1);
nValsNew=zeros(S,1);
invPerm=zeros(S,1);
newIdx=zeros(S,1);%Initialization
idx=zeros(S,1);%Initialization

iNew=order(1);
nValsNew(1)=nVals(iNew);

invPerm(iNew)=1;
cumProdNew(1)=1;

for i=2:S
    iNew=order(i);

    nValsNew(i)=nVals(iNew);
    invPerm(iNew)=i;

    cumProdNew(i)=cumProdNew(i-1)*nValsNew(i-1);
end

CPerm=zeros(nValsNew(:)');

newLinIdx=0;
linIdx=0;
while(1)
    %Assign the current tuple.
    CPerm(newLinIdx+1)=C(linIdx+1);

    %Move on to the next tuple.
    linIdx=linIdx+1;

    if(linIdx>=totalEls)
        break;
    end

    %The code below is similar to that in getNextTuple, because we are
    %updating the tuples for newIdx.
    curLevel=1;
    while(1)
        %Try incrementing the order at this level.
        idx(curLevel)=idx(curLevel)+1;
        if(idx(curLevel)<nVals(curLevel))
            iNew=invPerm(curLevel);
            newIdx(iNew)=newIdx(iNew)+1;
            newLinIdx=newLinIdx+cumProdNew(iNew);

            idx(1:(curLevel-1))=0;
            break;
        else
            %If the value is invalid, then just keep ascending and
            %adjust newLinIdx.
            iNew=invPerm(curLevel);

            newLinIdx=newLinIdx-newIdx(iNew)*cumProdNew(iNew);
            newIdx(iNew)=0;

            curLevel=curLevel+1;
            continue;
        end
    end
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
classdef testStressRotationInVoigtNotation < handle

    properties (Access = public)
        tol =  1e-14;
    end

    properties (Access = protected)
        direction
        stress
        stressVoigt
        rotatedStressByVoigt
        rotatedStress
    end

    properties (Access = private)
        angle
    end

    methods (Access = public)

        function error = computeError(obj)
            rotStre        = obj.rotatedStress.getValue();
            rotStreByVoigt = obj.rotatedStressByVoigt.getValue();
            error = norm(double(rotStre) - double(rotStreByVoigt));
        end

    end

    methods (Access = protected)

        function compute(obj)
            obj.createAngle()
            obj.createDirection()
            obj.createStress()
            obj.createRotatedStress()
            obj.createRotatedStressWithVoigtNotation()
        end

        function createStress(obj)
            obj.stress = Stress3DTensor;
            obj.stress.createRandomTensor();
            obj.stressVoigt = Tensor2VoigtConverter.convert(obj.stress);
        end

        function createRotatedStress(obj)
            rotS = obj.rotateStress();
            obj.rotatedStress = Tensor2VoigtConverter.convert(rotS);
        end

        function rotS = rotateStress(obj)
            a  = obj.angle;
            d  = obj.direction;
            rotS = Rotator.rotate(obj.stress,a,d);
        end

    end

    methods (Access = private)

        function createAngle(obj)
            obj.angle = 2*pi*rand(1);
        end

        function createRotatedStressWithVoigtNotation(obj)
            theta = obj.angle;
            dir   = obj.direction;
            stre = obj.stressVoigt;
            rotStre = Rotator.rotate(stre,theta,dir);
            obj.rotatedStressByVoigt = rotStre;
        end

    end

    methods (Abstract,Access = protected)
        createDirection(obj)
    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function unitCell = calcUnitCell(xy,varargin)
% compute the unit cell for an EBSD data set
%
% Input
%  xy - spatial coordinates
%
% Output
%  unitCell - coordinates of the unit cell
%
% Options
%
% GridType - [automatic, hexagonal, rectangular]

% nothing to do -> skip
if isempty(xy)
  unitCell = [];
  return;
end

if size(xy,2) == 3
  unitCell = [calcUnitCell(xy(:,[1 2]), varargin{:});...
    calcUnitCell(xy(:,[1 3]), varargin{:}); ...
    calcUnitCell(xy(:,[2 3]), varargin{:})];
  return
end

% first estimate of the grid resolution
area = (max(xy(:,1))-min(xy(:,1)))*(max(xy(:,2))-min(xy(:,2)));
dxy = sqrt(area / length(xy));

% compensate for single line EBSD
if dxy==0
  lx = mean(diff(xy(:,1))); ly = mean(diff(xy(:,2)));
  if lx==0, lx=ly; else; ly=lx; end
  dxy= (ly+ly)/2;
end

% reduce data set
if length(xy)>10000, xy = subSample(xy,10000); end

% remove dublicates from the coordinates
xy = uniquetol(xy,0.01/sqrt(size(xy,1)),'ByRows',true);

try
  % compute Voronoi decomposition
  [v, c] = voronoin(xy,{'Qz'});

  % compute the area of all Voronoi cells
  areaf = @(x,y) abs(0.5.*sum(x(1:end-1).*y(2:end)-x(2:end).*y(1:end-1)));
  areaf = cellfun(@(c1) areaf(v([c1 c1(1)],1),v([c1 c1(1)],2)),c);

  % the unit cell should be the Voronoi cell with the smalles area
  [~, ci] = min(areaf);

  % compute vertices of the unit cell
  unitCell = [v(c{ci},1) - xy(ci,1),v(c{ci},2) - xy(ci,2)];
  % sometimes it happens that we have one point doubled, remove those
  ignore = [false;sqrt(sum(diff(unitCell,1).^2,2)) < max(sqrt(sum(diff(unitCell,1).^2,2)))/5];
  unitCell(ignore,:) = [];

  if isRegularPoly(unitCell,varargin)
    return
  end

  % second estimate of the grid resolution
  dxy2 = max(unitCell) - min(unitCell);

  if 100*dxy2 > dxy, dxy = dxy2;end

catch %#ok<CTCH>
end

% get grid options
dxy = get_option(varargin,'GridResolution',dxy);
cellType = get_option(varargin,'GridType','rectangular');
cellRot = get_option(varargin,'GridRotation',0*degree);

% otherwise take a regular unit cell
switch lower(cellType)

  case 'rectangular'

    unitCell = regularPoly(4,dxy,cellRot);

  case 'hexagonal'

    unitCell = regularPoly(6,dxy,cellRot);

  case 'circle'

    unitCell = regularPoly(16,dxy,cellRot);

  otherwise

    error('MTEX:plotspatial:UnitCell','Unknown unit cell type!')
end

% a regular polygon with s vertices, diameter d, and rotation rot
function unitCell = regularPoly(s,d,rot)

c = exp(1i*((pi/s:pi/(s/2):2*pi)+rot))./sqrt((s/2));
unitCell = [real(c(:)),imag(c(:))].*d;

function isRegular = isRegularPoly(unitCell,varargin)

sideLength = sqrt(sum((unitCell).^2,2));
sides      = numel(sideLength);

uC = complex(unitCell(:,1),unitCell(:,2));
nC = uC([2:end 1]);

enclosingAngle = uC./nC;
enclosingAngle = complex(abs(real(enclosingAngle)),...
  abs(imag(enclosingAngle)));

isRegular = any(sides == [4 6]) && ... % norm(sideLength - mean(sideLength))*dxy < 1e-5 && ...
  norm(enclosingAngle - mean(enclosingAngle)) < 0.05*degree;

% find a quare subset of about N points
function xy = subSample(xy,N)

xminmax = [min(xy(:,1));max(xy(:,1))];
yminmax = [min(xy(:,2));max(xy(:,2))];

% shrink range until only N points are inside
while length(xy) > N

  if diff(xminmax) > diff(yminmax)
    xminmax = [3 1;1 3] * xminmax ./ 4;
  else
    yminmax = [3 1;1 3] * yminmax ./ 4;
  end

  xy = xy(xy(:,1)>xminmax(1) & xy(:,1)<xminmax(2) & ...
    xy(:,2)>yminmax(1) & xy(:,2)<yminmax(2),:);

end

author: mtex-toolbox

repo: mtex

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


text:
function [gX, kern] = pathKernGradSym(kern, X, X2)

% PATHKERNGRADSYM Gradient of the symetric PATH kernel's parameters.
% FORMAT
% DESC computes the gradient of functions with respect to the
% path
% kernel's parameters. As well as the kernel structure and the
% input positions, the user provides a matrix PARTIAL which gives
% the partial derivatives of the function with respect to the
% relevant elements of the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x : the input locations for which the gradients are being
% computed.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The argument takes
% the form of a square matrix of dimension  numData, where numData is
% the number of rows in X.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters. The ordering of the vector should match
% that provided by the function kernExtractParam.
% RETURN kern : the updated kernel structure
%
% FORMAT
% DESC computes the derivatives as above, but input locations are
% now provided in two matrices associated with rows and columns of
% the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x1 : the input locations associated with the rows of the
% kernel matrix.
% ARG x2 : the input locations associated with the columns of the
% kernel matrix.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The matrix should
% have the same number of rows as X1 and the same number of columns
% as X2 has rows.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters.
% RETURN kern : the updated kernel structure
%
% SEEALSO pathKernParamInit, kernGradient, pathKernDiagGradient
%
% COPYRIGHT : Andrea Baisero, Carl Henrik Ek, 2013

% SHEFFIELDML

maxl=max(cellfun(@(x)size(x,1),[X,X2]));
kern=pathKernUpdateWMat(kern,maxl);

num=length(X);
num2=length(X2);

gX=cell(num,num2);
for i=1:num
    for j=1:num2
        gX{i,j}=pathKernGradSymSeq(kern,X{i},X2{j});
    end
end

function gX = pathKernGradSymSeq(kern, x, x2)

l=size(x,1);
l2=size(x2,1);
dim=size(x,2);

t=kernGradX(kern.gkern,x,x2);
kw=kern.wmat(1:l,1:l2);

gX=zeros(size(x));
for i=1:l
    gX(i,:)=sum((kw(i,:)'*ones(1,dim)).*t(:,:,i));
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function test_failed = test_ambiguityfunction
Lr = [1, 19, 20];

test_failed = 0;

disp(' ===============  TEST_AMBIGUITYFUNCTION ==============');

for ii = 1: length(Lr)
  L = Lr(ii);
    for n = 1:4

    if (n==1)
    type1 = 'auto';
    type2 = 'real';
    f = tester_rand(L,1);
    g = f;
    elseif (n==2)
    type1 = 'auto';
    type2 = 'complex';
    f = tester_crand(L,1);
    g = f;
    elseif (n==3)
    type1 = 'cross';
    type2 = 'real';
    f = tester_rand(L,1);
    g = tester_rand(L,1);
    elseif (n==4)
    type1 = 'cross';
    type2 = 'complex';
    f = tester_crand(L,1);
    g = tester_crand(L,1);
    end

    r1 = ref_ambiguityfunction(f,g);
    r2 = ambiguityfunction(f,g);

    res = norm(r1-r2);

    [test_failed, fail] = ltfatdiditfail(res, test_failed);
    s = sprintf('DAF %3s %3s L:%3i %0.5g %s', type1, type2, L, res, fail);
    disp(s);
    end
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function Xcq = cqt(x, B, fs, fmin, fmax, varargin)
%CQT  Constant-Q/Variable-Q transform
%   Usage:  Xcq = cqt(x, B, fs, fmin, fmax, varargin)
%
%   Input parameters:
%         x         : input signal
%         B         : number of bins per octave
%         fs        : sampling frequency
%         fmin      : lowest frequency to be analyzed
%         fmax      : highest frequency to be analyzed
%         varargin  : Optional input pairs (see table below)
%
%   Output parameters:
%         Xcq       : Struct consisting of
%           .c           : CQT coefficients
%           .cDC         : transform coefficients for f = 0
%           .cNyq        : transform coefficients for fs/2
%           .g           : cell array of analysis filters
%           .shift       : center frequencies of analysis filters
%           .M           : bandwidth of analysis filters
%           .xlen        : length of input signal
%           .phasemode   : 'local'  -> zero-centered filtered used
%                        : 'global' -> mapping function used
%           .rast        : time-frequency plane sampling scheme (full,
%                          piecewise, none)
%           .fmin
%           .fmax
%           .B
%           .format      : eighter 'cell' or 'matrix' (only applies for
%                          piecewise rasterization)
%
%   Optional input arguments arguments can be supplied like this:
%
%       Xcq = cqt(x, B, fs, fmin, fmax, 'rasterize', 'piecewise')
%
%   The arguments must be character strings followed by an
%   argument:
%
%     'rasterize':  can be set to (default is 'full');
%           - 'none':      Hop sizes are distinct for each frequency
%                          channel. Transform coefficients will be
%                          presented in a cell array.
%           - 'full':      The hop sizes for all freqency channels are
%                          set to the smallest hop size in the representa-
%                          tion. Transform coefficients will be presented
%                          in matrix format.
%           - 'piecewise': Hop sizes will be rounded down to be a power-of-
%                          two integer multiple of the smallest hop size in
%                          the representation. Coefficients will be
%                          presented either in a sparse matrix or as cell
%                          arrays (see 'format' option)
%
%     'phasemode':  can be set to (default is 'global')
%           - 'local':     Zero-centered filtered used
%           - 'global':    Mapping function used (see reference)
%
%     'format':     applies only for piecewise rasterization
%           - 'sparse':   Coefficients will be presented in a sparse matrix
%           - 'cell':     Coefficients will be presented in a cell array
%
%     'gamma':      the bandwidth of each filter is given by
%                            Bk = 1/Q * fk + gamma,
%                   where fk is the filters center frequency, Q is fully
%                   determined by the number of bins per octave and gamma
%                   is a bandwidth offset. If gamma = 0 the obtained
%                   filterbank is constant-Q. Setting gamma > 0 time
%                   resolution towards lower frequencies can be improved
%                   compared to the constant-Q case (e.g. ERB proportional
%                   bandwidths). See reference for more information.
%     'normalize':  coefficient normalization
%          - 'sine':    Filters are scaled such that a sinusoid with
%                       amplitude A in time domain will exhibit the same
%                       amplitude in the time-frequency representation.
%          - 'impulse': Filters are scaled such that an impulse in time
%                       domain will exhibit a flat response in the
%                       time-frequency representation (in the frame that
%                       centers the impulse)
%          - 'none':      ...
%     'winfun':        defines the window function that is used for filter
%                   design. See winfuns for more information.
%
%   See also:  nsgtf_real, winfuns
%
%   References:
%     C. Sch�rkhuber, A. Klapuri, N. Holighaus, and M. D�rfler. A Matlab
%     Toolbox for Efficient Perfect Reconstruction log-f Time-Frequecy
%     Transforms.
%
%     G. A. Velasco, N. Holighaus, M. D�rfler, and T. Grill. Constructing an
%     invertible constant-Q transform with non-stationary Gabor frames.
%     Proceedings of DAFX11, Paris, 2011.
%
%     N. Holighaus, M. D�rfler, G. Velasco, and T. Grill. A framework for
%     invertible, real-time constant-q transforms. Audio, Speech, and
%     Language Processing, IEEE Transactions on, 21(4):775-785, April 2013.
%
%
%
% Copyright (C) 2013 Christian Sch�rkhuber.
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 3.0 Unported
% License. To view a copy of this license, visit
% http://creativecommons.org/licenses/by-nc-sa/3.0/
% or send a letter to
% Creative Commons, 444 Castro Street, Suite 900,
% Mountain View, California, 94041, USA.

% Authors: Christian Sch�rkhuber
% Date: 20.09.13

%% check input arguments

%defaults
rasterize = 'full'; %fully rasterized
phasemode = 'global';
outputFormat = 'sparse'; %only applies if rasterize == 'octave'
normalize = 'sine';
windowFct = 'hann';
gamma = 0;

if nargin >= 6
    Larg = length(varargin);
    for ii=1:2:Larg
       switch varargin{ii}
           case {'rasterize'}
               rasterize = varargin{ii+1};
           case {'phasemode'}
               phasemode = varargin{ii+1};
           case {'format'}
               outputFormat = varargin{ii+1};
           case {'gamma'}
               gamma = varargin{ii+1};
           case {'normalize'}
               normalize = varargin{ii+1};
           case {'win'}
               windowFct = varargin{ii+1};
       end
    end
end

%% window design
[g,shift,M] = nsgcqwin(fmin,fmax,B,fs, length(x), ...
        'winfun', windowFct, 'gamma', gamma, 'fractional', 0);
 fbas = fs*cumsum(shift(2:end))./ length(x);
 fbas = fbas(1:size(M,1)/2-1);

%% compute coefficients
bins = size(M,1)/2 - 1;
switch rasterize
    case 'full'
        M(2:bins+1) = M(bins+1);
        M(bins+3:end) = M(bins+1:-1:2);

    case 'piecewise'
        temp = M(bins+1);
        octs = ceil(log2(fmax/fmin));
        %make sure that the number of coefficients in the highest octave is
        %dividable by 2 at least octs-times
        temp = ceil(temp/2^octs)*2^octs;
        mtemp = temp./ M;
        mtemp = 2.^( ceil(log2(mtemp)) -1);
        mtemp = temp./ mtemp;
        mtemp(bins+2) = M(bins+2); %don't rasterize Nyquist bin
        mtemp(1) = M(1); %don't rasterize DC bin
        M = mtemp;

    otherwise
end

switch normalize
    case {'sine','Sine','SINE','sin'}
        normFacVec = 2*M(1:bins+2)./length(x);
    case {'impulse','Impulse', 'IMPULSE','imp'}
        normFacVec = 2*M(1:bins+2)/cellfun(@length,g);
    case {'none','None','NONE','no'}
        normFacVec = ones(bins+2,1);
    otherwise
        error('Unkown normalization method!');
end

normFacVec = [normFacVec; normFacVec(end-1:-1:2)];
g = arrayfun(@(k) (g{k}*normFacVec(k)),1:(2*bins+2),'UniformOutput',0).';

c = nsgtf_real(x,g,shift,M,phasemode);

switch rasterize
    case 'full'
        cDC = cell2mat(c(1)).';
        cNyq = cell2mat(c(bins+2)).';
        c = cell2mat(c(2:bins+1).');
    case 'piecewise'
        cDC = cell2mat(c(1));
        cNyq = cell2mat(c(bins+2));
        if strcmp(outputFormat,'sparse')
            c = cqtCell2Sparse(c,M).';
        else
            c = c(2:end-1);
        end

    otherwise
        cDC = cell2mat(c(1));
        cNyq = cell2mat(c(end));
        c = c(2:end-1);
end

%% output
Xcq = struct('c', {c.'}, 'g', {g}, 'shift', shift, 'M', {M}, ...
    'xlen', length(x), 'phasemode', phasemode, 'rast', rasterize, ...
    'fmin', fmin, 'fmax', fmax, 'B', B, 'cDC', cDC, 'cNyq', cNyq, ...
    'format', outputFormat, 'fbas', fbas);

author: nesl


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [test_kernel,train_kernel,optimal_alpha] = DSK_optimization(train_data,train_label,test_data,opt)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input parameters:
% train_data: column cells containing the SPD matrices for training
% train_label: one column vector containing the labels for the training data
% test_data: test data with the same format as the training data
% opt:  a structure containing parameter settings
%       elements:
%       theta -- a kernel parameter
%       obj_method -- which criterion to be used
%       original_alpha  set to 1 to use original Stein kernel or 0 to use DSK
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Output parameters:
% test_kernel: the adjusted test_kernel
% train_kernel: the adjusted train_kernel
% optimal_alpha: the optimized adjustment parameters alpha
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Jianjia Zhang, jz163@uowmail.edu.au Dec, 2014, all rights reserved
% For implementation details, please refer to:
% "Learning Discriminative Stein Kernel for SPD Matrices and Its Applications."
% arXiv preprint arXiv:1407.1974 (2014).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [test_kernel,train_kernel,optimal_alpha] = DSK_optimization_new(TrainSet,TestSet,opt)
nmode = size(TrainSet.X_cov,1); % dimension of the SPD matrices
train_decomp = Decomposite_eig_new(TrainSet); % eigen decomposition of the training/test data
test_decomp = Decomposite_eig_new(TestSet);
%%%%%%%%%%%%%%%%%%%%%%%%%
if(~opt.original_alpha)
    initial_alpha = 1*ones(1,nmode); % the initial alpha corresponding to the original Stein kernel
    LB = 0.01*initial_alpha;

    options = optimset('Algorithm','interior-point'); % run interior-point algorithm
    options.Display = 'iter';
    options.Display = 'off';
    options.MaxIter = 100;
    options.TolFun = 1e-5;
    %tic
    optimal_alpha = fmincon(@(alpha) objfun_ff_new(alpha,TrainSet.y,train_decomp,opt.lambda,initial_alpha,opt.obj_method,opt.theta),initial_alpha,[],[],[],[],LB,[],[],options);
    %toc
else
    optimal_alpha = ones(1,nmode);
end
[S_test] = EigComp2SD_power_new(train_decomp,test_decomp,optimal_alpha); % compute the Stein divergence with the obtained adjustment parameter optimal_alpha
[S_train] = EigComp2SD_power_new(train_decomp,train_decomp,optimal_alpha);
test_kernel = exp(-1*opt.theta*S_test); % compute the kernel
train_kernel = exp(-1*opt.theta*S_train);
end

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function dispma3()
    % compare two rates, no fit
    % selects 4 times to define begining and end of two segments in
    % cumulative number curve and calls bvalnofit
    %
    % turned into function by Celso G Reyes 2017

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

    t1p = min(ZG.newt2.Date);
    t4p = max(ZG.newt2.Date);
    t2p = t4p - (t4p-t1p)/2;
    t3p = t2p;

    show_dialog()

    function show_dialog()
        % Input times t1p t2p t3p and t4p by editing or use cursor if desired
        %

        dlg=ZmapDialog();
        dlg.AddEdit('t1p','Time 1 (T1):', t1p,'enter begin segment 1');
        dlg.AddEdit('t2p','Time 2 (T2):', t2p,'enter end segment 1');
        dlg.AddEdit('t3p','Time 3 (T3):', t3p,'enter begin segment 2');
        dlg.AddEdit('t4p','Time 4 (T4):', t4p,'enter end segment 2');
        dlg.AddCheckbox('usecursor','Use Cursor', false,[],'select with mouse');

        par2 = 0.1 * ZG.newt2.Count;

        [myans, okpressed] =dlg.Create('Name', 'Select two time segments');
        if ~okpressed
            return
        end
        if myans.usecursor
            mouse_select()
        else
            bvanofit(ZG.newt2, [myans.t1p,myans.t2p] , [myans.t3p, myans.t4p]);
        end
    end

    function mouse_select()
        figure(findobj('Tag','cum','-and','Type','Figure'));

        seti = uicontrol('Units','normal','Position',[.4 .01 .2 .05],'String','Select T1  ');

        pause(0.5)

        par2 = 0.1 * ZG.newt2.Count;
        par3 = 0.12 * ZG.newt2.Count;
        t1 = ginput(1);
        t1p = [  t1 ; t1(1) t1(2)-par2; t1(1)   t1(2)+par2 ];
        plot(t1p(:,1),t1p(:,2),'r')
        text( t1(1),t1(2)+par3,['t1: ', num2str(t1p(1))] )
        seti.String = 'select T2';

        pause(0.5)

        t2 = [];
        t2 = ginput(1);
        t2p = [  t2 ; t2(1) t2(2)-par2; t2(1)   t2(2)+par2 ];
        plot(t2p(:,1),t2p(:,2),'r')
        text( t2(1),t2(2)+par3,['t2: ', num2str(t2p(1))] )
        seti.String = 'select T3';

        pause(0.5)

        t3 = [];
        t3 = ginput(1);
        t3p = [  t3 ; t3(1) t3(2)-par2; t3(1)   t3(2)+par2 ];
        plot(t3p(:,1),t3p(:,2),'r')
        text( t3(1),t3(2)+par3,['t3: ', num2str(t3p(1))] )
        seti.String = 'select T4';

        pause(0.5)

        t4 = [];
        t4 = ginput(1);
        t4p = [  t4 ; t4(1) t4(2)-par2; t4(1)   t4(2)+par2 ];
        plot(t4p(:,1),t4p(:,2),'r')
        text( t4(1),t4(2)+par3,['t4: ', num2str(t4p(1))] )

        delete(seti)

        show_dialog()
        pause(0.1)
    end
end

author: CelsoReyes


text:
% IM = mkAngle(SIZE, PHASE, ORIGIN)
%
% Compute a matrix of dimension SIZE (a [Y X] 2-vector, or a scalar)
% containing samples of the polar angle (in radians, CW from the
% X-axis, ranging from -pi to pi), relative to angle PHASE (default =
% 0), about ORIGIN pixel (default = (size+1)/2).

% Eero Simoncelli, 6/96.

function [res] = mkAngle(sz, phase, origin)

sz = sz(:);
if (size(sz,1) == 1)
  sz = [sz,sz];
end

% -----------------------------------------------------------------
% OPTIONAL args:

if (exist('origin') ~= 1)
  origin = (sz+1)/2;
end

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

[xramp,yramp] = meshgrid( [1:sz(2)]-origin(2), [1:sz(1)]-origin(1) );

res = atan2(yramp,xramp);

if (exist('phase') == 1)
  res = mod(res+(pi-phase),2*pi)-pi;
end

repo: SelfExSR


text:
%These tests are using the Human-GEM model.

%%%%%%%%%%%%%%%%%%%%%%%%%%
% TC0901 - Test that all reactions in the minimal model to be
%          used by the MILP can carry flux > 0.1
%%%%%%%%%%%%%%%%%%%%%%%%%%

%In the MILP, the flux of reactions with negative rxn scores cannot exceed 100
%At the same time, the flux of positive reactions need to be >= 0.1 to be able
%to turn them on. Likewise, we force the flux on essential reactions to be >= 0.1.
%Therefore, it is important that all reactions can get a flux > 0.1, while all
%reactions are limited to a flux <= 100. We can test this in a similar way to the
%haveFlux function in RAVEN, so that code is copied here and modified.

%This code is partly copied from the function haveFlux
%The prepData takes a couple of hours to generate. If you have it somewhere,
%better to load it
%cd C:\Work\MatlabCode\components\human-GEM\Human-GEMftINIT\Human-GEM %replace this with the root of the Human-GEM repo.
%ihuman = importYaml('model/Human-GEM.yml');
%save('model/Human-GEM.mat', 'ihuman');
%load('model/Human-GEM.mat')
%prepData = prepHumanModelForftINIT(ihuman, false);
prepData2 = load('model/PrepData2.mat').prepData;

upperLimit = 100;
lowerLimit = 0.1;

model = prepData2.minModel;
model = prepData.minModel;

model.lb(model.lb==-inf)=-upperLimit;
model.ub(model.ub==inf)=upperLimit;
model.lb(model.lb==-1000)=-upperLimit;
model.ub(model.ub==1000)=upperLimit;

%check
unique(model.lb) %ok, only 0 and -100
unique(model.ub) %ok, only 100

%First make a loop where we optimize for all
onFwd = false(numel(model.rxns,1));
onRev = false(numel(model.rxns,1));

iter = 1;
while true
    disp([num2str(iter) ': ' num2str(sum(onFwd))])
    model.c=ones(numel(model.c),1);
    model.c(onFwd) = 0; %don't include the ones already verified

    sol=solveLP(model);
    if isempty(sol.x)
        disp('Failed');
        break;
    end

    currRes = sol.x > lowerLimit;
    if sum(currRes & ~onFwd) == 0
        break;%we didn't turn any more rxns on
    end

    onFwd = onFwd | currRes;
    iter = iter + 1;
end

sum(~onFwd)%zero, perfect

%Now the other direction (reversible):
onRev = false(numel(model.rxns,1));

iter = 1;
while true
    disp([num2str(iter) ': ' num2str(sum(onRev))])
    model.c=-ones(numel(model.c),1);
    model.c(~model.rev) = 0;
    model.c(onRev) = 0; %don't include the ones already verified

    sol=solveLP(model);
    if isempty(sol.x)
        disp('Failed');
        break;
    end

    currRes = sol.x < -lowerLimit;
    if sum(currRes & ~onRev) == 0
        break;%we didn't turn any more rxns on
    end

    onRev = onRev | currRes;
    iter = iter + 1;
end

sum(~onRev(model.rev == 1))%

onRev = false(numel(model.rxns,1));

indLeft = find(model.rev == 1 & ~onRev);

constructEquations(model, model.rxns(indLeft))%no scaling here...

%go trough the last one by one
for i = 1:numel(indLeft)
    model.c=zeros(numel(model.c),1);
    model.c(indLeft(i)) = -1;

    sol=solveLP(model);
    if isempty(sol.x)
        disp('Failed');
        break;
    end

    currRes = sol.x(indLeft(i)) < -lowerLimit;
    disp([num2str(indLeft(i)) ': ' num2str(currRes)])
    onRev(indLeft(i)) = currRes;
    onRev = onRev | currRes;
    iter = iter + 1;
end

sum(~onRev(model.rev == 1))%0, perfect. This means that we can use the limits 0.1 and 100 in the MILP!

author: SysBioChalmers

repo: RAVEN


text:
% Compute geodesic distance along the mesh
%[V,F] = readOBJ('data/spot_good.obj'); % this one will work
[V,F] = readOBJ('data/spot_bad.obj');  % but this one will not!

% NOTE: The result may not be the same on the "bad" mesh after you repair
% it, because this code computes the distance from vertex "1", but which
% vertex is vertex "1" may change. That is fine, as long as it outputs
% reasonable a looking distance.

% Run the algorithm
% (I promise nothing is wrong with this algorithm, the problem is the mesh)
dists = heat_geodesic(V,F,1);

% Plot the resulting function
t = tsurf(F,V, 'CData',dists);
shading interp;
axis equal;
axis off;
colormap(cbrewer('Blues', 500));
light('Position',[-1.5 1 1],'Style','local');
lights = camlight;
set(t, 'FaceLighting','gouraud', 'FaceColor','interp');
set(t, 'DiffuseStrength',0.5, 'SpecularStrength',0.2, 'AmbientStrength',0.3);
camproj('perspective');

author: odedstein

repo: sgi-introduction-course


text:
% EASY6  computes  a baseline vector. With given  code  and carrier
%            phase observations we estimate the ambiguities using integer
%	         least-squares (ILS) or Goad algorithms and estimate the baseline
%            components by a Kalman filter. The code does not handle
%	        outliers.

% RINEX version3.03

%Kai Borre 27-07-2002
%Copyright (c) by Kai Borre
%$Revision: 1.1 $  $Date: 2005/01/20  $
%Total revision 2.0, February 15, 2016

% The following code is general, but with the following limitations

% 1. We assume the observational columns are identical and ordered in the
% same sequence in master and rover files. The present observations just
% by chance contain identical satellites which are ordered in the same
% sequence. So sats1 and sats2 are the same. Alternatively, we need to
% establish a matching between sats1 and sats2 and likewise for Obs1
% and Obs2

% 2. There are no satellites with elevation angle less than 10 degrees.
% Minimum is 12 degrees. There are no rising or setting satellites during
% the observation period

% 3. The satellite having highest elevation angle at the start becomes the
% second highest about the middle of the observation file.

% 4. We do not apply tropospheric correction as it is very small in case
% master and rover stations only are about 1 meter apart

% 5. the actual baseline lenght was measured as 0.645 meter.

% Initial setting and computation of constants
v_light = 299792458;	     % vacuum speed of light m/s
f1 = 154*10.23E6;		     % L1 frequency Hz
f2 = 120*10.23E6;			 % L2 frequency Hz
lambda1 = v_light/f1;	     % wavelength on L1:  .19029367  m
lambda2 = v_light/f2;	     % wavelength on L2:  .244210213 m

% Read RINEX ephemerides file and convert to internal Matlab format
rinexe('log_24h.15n','eph.dat');
Eph = get_eph('eph.dat');

% Selection of observation type
%  ss = 'C1W L1W C2W L2W' %; % pseudorange and carrier phase
ss = 'C1C L1C C2X L2X'; % just to test other obsrervation types
fprintf('\n %s', ss)

% Open the master observation file
ofile1 = 'log_24h.15o';
fid1 = fopen(ofile1,'rt');

% Open the rover observation file
ofile2 = 'log_r.15o';
fid2 = fopen(ofile2,'rt');

coo = [];
ant_delta = [];
linjer = 0;
ij = [];
OBS1 = [];
OBS2 = [];
x = zeros(3,1);

% Gobbling the master file header, and collecting useful information
while 1
    linjer = linjer +1;
    line = fgetl(fid1);
    answer = strfind(line,'END OF HEADER');
    if  ~isempty(answer), break; end;
    if (line == -1), eof = 1; break; end;

    answer = strfind(line,'REC #');
    if ~isempty(answer)
        var = textscan(fid1,'%s %d14','Delimiter','\n');
        var1 = var{1}{1};
        cox = str2double(var1(1:14));
        coy = str2double(var1(15:28));
        coz = str2double(var1(29:42));
    end

    answer = strfind(line,'ANT #');
    if ~isempty(answer)
        var = textscan(fid1,'%s %d14','Delimiter','\n');
        var1 = var{1}{1};
        ax = str2double(var1(1:14));
        ay = str2double(var1(15:28));
        az = str2double(var1(29:42));
    end

    answer = strfind(line,'SYS / # / OBS TYPES');
    if ~isempty(answer)
        tline1 = strsplit(line);
        line = fgetl(fid1);
        tline2 = strsplit(line);
        tt1 = horzcat(tline1,tline2);
        for qq = 1:4
            if qq == 1, i = strcmp(tt1,ss(1:3)); end
            if qq == 2, i = strcmp(tt1,ss(5:7)) ; end
            if qq == 3, i = strcmp(tt1,ss(9:11)); end
            if qq == 4, i = strcmp(tt1,ss(13:15)); end
            ii = find(i == 1) ;
            ii = ii-2;
            if ii > 15,  ii = ii-7; end;
            ij = [ij ii];
        end
        % the cell array of strings tline1 includes two strings in the start
        % which describe the system and number of observation types.
        % Both tline1 and tline2 termintes with six additional strings.
        % An extra string appears at  the start of tline2; it originates
        % from concatenation of  the two lines. The indexing does not
        % change even if you empty the cells. They remain as empty
        % cells and keep a place
        obs_col = ij;
    end;
    answer = strfind(line,'INTERVAL');
    if ~isempty(answer)
        interval = strtok(line);
        int = str2double(interval);
    end;
    answer = strfind(line,'# OF SATELLITES');
    if ~isempty(answer)
        NumberSV = strtok(line);
        numSV = str2double(NumberSV) ;
    end
end % end reading master header

% the string arrays for the tline1 and tline2 contain an integer after the
% carrier phase observation. We must account for this by the following
% correctional table.

% the following code may be done more elegantly
for qq = 1:4
    if qq == 1, a0 = 6; b0 = 7; end
    if qq == 2, a0 = 6; b0 = 7; end
    if qq == 3, a0 = 14; b0 =15; end
    if qq == 4, a0 =14; b0 =15; end
    if strcmp(ss(a0:b0), '1C'), obs_col(1,qq) = obs_col(1,qq)+1;
    elseif  strcmp(ss(a0:b0), '1W'), obs_col(1,qq) = obs_col(1,qq)+2;
    elseif strcmp(ss(a0:b0),  '2X'), obs_col(1,qq) = obs_col(1,qq) +3;
    elseif strcmp(ss(a0:b0),  '2W'), obs_col(1,qq) = obs_col(1,qq) +4;
    else  strcmp(ss(a0:b0),  '5X'), obs_col(1,qq) = obs_col(1,qq) +5;
    end;
end;

Pos = [];
Omc = [];
base = [];
x_acc = [];
innov = [];

[time,post] = textscan(fid1,'%s %d8','Delimiter','\n');
tid = time{1}{1};
year = str2double(tid(3:6));
month = str2double(tid(8:9));
day = str2double(tid(11:12));
hour = str2double(tid(14:15));
minute = str2double(tid(17:18));
second = str2double(tid(20:29));
% static = str2double(tid(31:32));
NoSvs1 = str2double(tid(34:36));
dt1  = str2double(tid(38:56));
h = hour+minute/60+second/3600;
jd = julday(year, month, day, h);
[~, sec_of_week] = gps_time(jd);
time1 = sec_of_week; % sow1

Obs1 = zeros(NoSvs1,length(obs_col));
for i = 1:NoSvs1
    [obs1,~] = textscan(fid1,'%s %d8','Delimiter','\n');
    obs1x = obs1{1}{1} ;
    obs1y = strsplit(obs1x);
    sat1 = obs1y{1} ;
    sats1(i) = str2double(sat1(2:3));
    Obs1(i,:) = str2double(obs1y(obs_col));
end

% Gobbling the rover file header, and collecting useful information
% the whole story repated for the rover with index 2
% we start by gobbeling the header
for qq = 1:linjer +3 % 1 plus 2 from applying textscan twice
    aa = fgetl(fid2);
end

% Reading header of observation record in rover file
[time,pos] = textscan(fid2,'%s %d8','Delimiter','\n');
tid = time{1}{1};
year = str2double(tid(3:6));
month = str2double(tid(8:9));
day = str2double(tid(11:12));
hour = str2double(tid(14:15));
minute = str2double(tid(17:18));
second = str2double(tid(20:29));
static = str2double(tid(31:32));
NoSvs2 = str2double(tid(34:36));
dt2  = str2double(tid(38:56));
h = hour+minute/60+second/3600;
jd = julday(year, month, day, h);
[~, sec_of_week] = gps_time(jd);
time2 = sec_of_week; % sow2

Obs2 = zeros(NoSvs2,length(obs_col));
for i = 1:NoSvs2
    [obs2,pos] = textscan(fid2,'%s %d8','Delimiter','\n');
    obs2y = obs2{1}{1};
    obs2 = strsplit(obs2y);
    sat2 = obs2{1};
    sats2(i,:) = str2double(sat2(2:3));
    Obs2(i,:) = str2double(obs2(obs_col));
end

% Establish synchronous reading of fid1 and fid2
% The epoch interval is int seconds
if  time1 > time2
    for j = 1:round(time1-time2)
        [time,pos] = textscan(fid2,'%s %d8','Delimiter','\n');
        tid = time{1}{1};
        year = str2double(tid(3:6));
        month = str2double(tid(8:9));
        day = str2double(tid(11:12));
        hour = str2double(tid(14:15));
        minute = str2double(tid(17:18));
        second = str2double(tid(20:29));
        static = str2double(tid(31:32));
        NoSvs2 = str2double(tid(34:36));
        dt2  = str2double(tid(38:56));
        h = hour+minute/60+second/3600;
        jd = julday(year, month, day, h);
        [~, sec_of_week] = gps_time(jd);
        time2 = sec_of_week; % sow2

        Obs2 = zeros(NoSvs2,length(obs_col));
        for i = 1:NoSvs2
            [obs2,pos] = textscan(fid2,'%s %d8','Delimiter','\n');
            obs2y = obs2{1}{1};
            obs2 = strsplit(obs2y);
            sat2 = obs2{1};
            sats2(i,:) = str2double(sat2(2:3));
            Obs2(i,:) = str2double(obs2(obs_col));
        end
    end
end

if  time1 < time2
    for j = 1:round(time2-time1)
    [time,pos] = textscan(fid1,'%s %d8','Delimiter','\n');
    tid = time{1}{1};
    year = str2double(tid(3:6));
    month = str2double(tid(8:9));
    day = str2double(tid(11:12));
    hour = str2double(tid(14:15));
    minute = str2double(tid(17:18));
    second = str2double(tid(20:29));
    static = str2double(tid(31:32));
    NoSvs1 = str2double(tid(34:36));
    dt1  = str2double(tid(38:56));
    h = hour+minute/60+second/3600;
    jd = julday(year, month, day, h);
    [~, sec_of_week] = gps_time(jd);
    time1 = sec_of_week; % sow1
    Obs1 = zeros(NoSvs1,length(obs_col));
    for i = 1:NoSvs1
        [obs1,pos] = textscan(fid1,'%s %d8','Delimiter','\n');
        obs1x = obs1{1}{1};
        obs1y = strsplit(obs1x);
        sat1 = obs1y{1};
        sats1(i) = str2double(sat1(2:3));
        Obs1(i,:) = str2double(obs1y(obs_col));
    end
    end
end % synchronization

fprintf('\nFirst Common Epoch Time1 :  %8.0f', time1)
fprintf('\nFirst Common Epoch Time2 :  %8.0f\n', time2)

b_acc = [];
innov_acc = [];
bk_acc = [];
epoch = 0;

disp(' ')
disp(' ... Now wait  a few minutes for further results and a plot!')

while 1
    epoch = epoch +1;
    % Reading header of observation block at Master
    [time,post] = textscan(fid1,'%s %d8','Delimiter','\n');
    tid = time{1}{1};
    year = str2double(tid(3:6));
    month = str2double(tid(8:9));
    day = str2double(tid(11:12));
    hour = str2double(tid(14:15));
    minute = str2double(tid(17:18));
    second = str2double(tid(20:29));
    % static = str2double(tid(31:32));
    NoSvs1 = str2double(tid(34:36));
    dt1  = str2double(tid(38:56));
    h = hour+minute/60+second/3600;
    jd = julday(year, month, day, h);
    [~, sec_of_week] = gps_time(jd);
    time1 = sec_of_week; % sow1

    % Reading observation block at Master
    Obs1 = zeros(NoSvs1,length(obs_col));
    for i = 1:NoSvs1
        obs1 = textscan(fid1,'%s %d8','Delimiter','\n');
        obs1x = obs1{1}{1} ;
        obs1y = strsplit(obs1x);
        sat1 = obs1y{1} ;
        sats1(i) = str2double(sat1(2:3));
        Obs1(i,:) = str2double(obs1y(obs_col));
    end

    % Reading header of observation record in rover file
    [time,pos] = textscan(fid2,'%s %d8','Delimiter','\n');
    tid = time{1}{1};
    year = str2double(tid(3:6));
    month = str2double(tid(8:9));
    day = str2double(tid(11:12));
    hour = str2double(tid(14:15));
    minute = str2double(tid(17:18));
    second = str2double(tid(20:29));
    static = str2double(tid(31:32));
    NoSvs2 = str2double(tid(34:36));
    dt2  = str2double(tid(38:56));
    h = hour+minute/60+second/3600;
    jd = julday(year, month, day, h);
    [~, sec_of_week] = gps_time(jd);
    time2 = sec_of_week; % sow2

    Obs2 = zeros(NoSvs2,length(obs_col));
    for i = 1:NoSvs2
        [obs2,pos] = textscan(fid2,'%s %d8','Delimiter','\n');
        obs2y = obs2{1}{1};
        obs2 = strsplit(obs2y);
        sat2 = obs2{1};
        sats2(i,:) = str2double(sat2(2:3));
        Obs2(i,:) = str2double(obs2(obs_col));
    end

    % Next we test if all observed SVs have an ephemeris.
    % sats1 contains the SVs as read in the observation lines
    % The intersect command delivers Sats in sorted order!
    % Therefore we must be careful in the follwing manipulations
    Sats = intersect(sats1,Eph(1,:));

    % The command ismember does not change the sequence of entries in sats1
    lia = ismember(sats1,Sats);
    % A 0 (zero) in lia indicates that a SV has been deleted. We delete the
    % corresponding row in the observations
    sats1(lia==0) = [];
    Obs1(lia==0,:) = []; % modfied because of two columns
    NoSV1 = length(sats1);

    % Next we test if all observed sats have an ephemeris.
    % sats contains the SVs as read in the observation lines.
    % The intersect command delivers Sats in sorted order!
    % Therefore we must be careful in the following manipulations
    Sats = intersect(sats2,Eph(1,:));

    % The command ismember does not change the sequence of entries in sats
    lia = ismember(sats2,Sats);
    % A 0 (zero) in lia indicates that a SV has been deleted. We delete the
    % corresponding row in the observations
    sats2(lia==0) = [];
    Obs2(lia==0, :) = [];
    NoSV2 = length(sats2);
    % end reading first record in rover file

    % By chance sats1 and sats 2 are the same in the sample data. Otherwise
    % we need to establish  a matching between sats1 and sats2 and
    % likewise for Obs1 and Obs2

    % We compute elevation angles for all sats1 and sats2
    [X_i, el] = recpo_ls(Obs1(:,3), sats1, time1, Eph);
    [X_j0,~] = recpo_ls(Obs2(:,3), sats2, time1, Eph);
    if epoch == 1
        X_j = X_j0(1:3,1);
    end
    % if epoch == 1
    % forcing the choice of reference satellite to be the second highest
    % at the start
    %  sats1
    % refSV = 5
    %  sats1(refSV) = [];
    % end

    if epoch == 1
        % The SV with largest elevation is taken as reference SV
        [y,refSV] = max(el) ;
        old_refSV = refSV;
        fprintf('\n Reference satellite is chosen as  PRN %3.0f\n', sats1(refSV))
        if refSV ~= old_refSV
            fprintf('\n New reference satellite is chosen as  PRN  %3.0f\n\n', sats1(refSV))
            sats1 = union(sats1, old_refSV);
            % no  PRN 12 must be inserted at the entry from where it was deleted!!!
        end
    end
    %The refSV is PRN 12 has index 6 in sats1
    % Finding columns in Eph for each SV
    m = NoSV1;
    col_Eph = zeros(1,m);
    for t = 1:m
        col_Eph(t) = find_eph(Eph,sats1(t),time1+dt1);
    end

    m1 = m-1; % m1 is all sv.s minus refSV
    X_a = [];
    X = zeros(3+2*m1,1);

    % creating a row of indices for all Svs minus the refSV
    t = 1: NoSV1;
    t1 = setdiff(t,refSV);
    % Computation of covariance matrix Sigma for double differenced observations
    D = [ones(m1,1) -eye(m1) -ones(m1,1) eye(m1)];
    Sigma = D*D';

    % A block for estimation of AMBIGUITIES
    % we accumulate the normals for the first 10 epochs
    switch(epoch)

        case {1, 2, 3, 4, 5, 6, 7, 8, 9,10}

            N = zeros(3+2*m1,3+2*m1);	  % initialization of normals
            rs = zeros(3+2*m1,1);	            % initialization of right side

            % Computing rho for refSV
            [~,rhok_j,~] = get_rho(time1, Obs2(refSV,1), Eph(:,col_Eph(refSV)), X_j);
            [~,rhok_i,Xk_ECF] = get_rho(time1, Obs1(refSV,1), Eph(:,col_Eph(refSV)), X_i);

            tt = 0;
            ts = length(t1);
            b = zeros(ts,1);
            bk = zeros(ts,1);
            A1 = [];
            for t = t1
                tt = tt+1;
                [~,rhol_j, ~] = get_rho(time1,Obs2(t,1), Eph(:,col_Eph(t)), X_j);
                [~,rhol_i, Xl_ECF] = get_rho(time1,Obs1(t,1), Eph(:,col_Eph(t)), X_i);
                A0 = [(Xk_ECF(1)-X_j(1))/rhok_j - (Xl_ECF(1)-X_j(1))/rhol_j  ...
                    (Xk_ECF(2)-X_j(2))/rhok_j - (Xl_ECF(2)-X_j(2))/rhol_j ...
                    (Xk_ECF(3)-X_j(3))/rhok_j - (Xl_ECF(3)-X_j(3))/rhol_j];
                A1 = [A1; A0];
                Phi1 = (Obs1(refSV,2)-Obs1(t,2)-Obs2(refSV,2)+Obs2(t,2))*lambda1;
                Phi2 = (Obs1(refSV,4)-Obs1(t,4)-Obs2(refSV,4)+Obs2(t,4))*lambda2;
                b(tt,:) =                   Phi1-lambda1*X(3+tt,1);
                b(length(t1)+tt,:) = Phi2-lambda2*X(3+length(t1)+tt,1);
                bk(tt,:) =                    rhok_i-rhok_j-rhol_i+rhol_j;
                bk(length(t1)+tt,:) =  rhok_i-rhok_j-rhol_i+rhol_j;
            end;
            A_modi = eye(m1);	    	          % modified coefficient matrix
            A_modi(:,refSV) = -ones(m1,1);
            A_aug = [A1 lambda1*A_modi 0*eye(m1); A1 0*eye(m1) lambda2*A_modi];
            N = N+A_aug'*kron(eye(2),Sigma)*A_aug;
            rs = rs+A_aug'*kron(eye(2),Sigma)*(b-bk);

        case 11

            PP = pinv(N);
            % X contains the three baseline components and next the float ambiguities
            X = PP*rs;
            % The next step is to convert the real ambiguities in X (all entries below the
            % fourth component and to the end) into integers. We indicate two algorithms:

            % the Integer Least-Squares (ILS) algorithm
            PP = 0.5*(PP+PP');  %to make PP symmetric
            %[a, sqnorm] = ILS( X(4:4+2*m1-1,1), PP(4:4+2*m1-1,4:4+2*m1-1),2);
            % Correcting to baseline vector as consequence of changing float ambiguities
            % to fixed ones
            % X(1:3,1) = X(1:3,1)-PP(1:3,4:4+2*m1-1)*inv(PP(4:4+2*m1-1,4:4+2*m1-1))*...
            %     (X(4:4+2*m1-1,1)-a(:,1)); %select SECOND  first set of candidates
            % X(4:4+2*m1-1,1) = a(:,1);
            %a = a(:,2); % forces the second set of N values in use

            % the Goad algorithm
            for j = 1:tt
                K1 = round(X(3+tt)-X(3+m1+tt));
                K2 = round(60*X(3+tt)-77*X(3+m1+tt));
                trueN2 = round((60*K1-K2)/17);   % (15.10)
                trueN1 = round(trueN2+K1);         % (15.11)
                amb(j,1:2) = [trueN1 trueN2];
            end
            a =  [amb(:,1);  amb(:,2) ];

            fprintf('\n N1 for PRN %3.0f: %3.0f',[sats1(t1); a(1:m1,1)'])
            fprintf('\n')
            fprintf('\n N2 for PRN %3.0f: %3.0f',[sats1(t1); a(m1+1:2*m1,1)'])
            fprintf('\n')

            % Setting covariances for the Kalman filter; the state vector contains (x,y,z)
            P =  10^8*eye(3);	                       		 % covariances of state vector
            Q = 0.1*eye(3);			                          % covariances of system
            R = 0.1*kron(eye(2),inv(Sigma));	    % covariances of observations
            [phi_j,lambda_j,h_j] = togeod(63781137, 298.257223563, X_j(1), X_j(2), X_j(3));
            h_i = h_j;

        otherwise

            % estimate of baseline components
            tt = 0;  % a counter
            A = zeros(m1,3); %zeros(2*m1,3);
            % Computing rho for refsv
            [torr2,rhok_j,~] = get_rho(time2, Obs2(refSV,3), Eph(:,col_Eph(refSV)), X_j);
            [torr1,rhok_i,Xk_ECF] = get_rho(time1, Obs1(refSV,3), Eph(:,col_Eph(refSV)), X_i);
            for t = t1
                tt = tt+1;
                [torr4,rhol_j, ~] = get_rho(time2,Obs2(t,3), Eph(:,col_Eph(t)), X_j);
                [torr3,rhol_i,Xl_ECF] = get_rho(time1,Obs1(t,3), Eph(:,col_Eph(t)), X_i);
                A0 = [(Xk_ECF(1)-X_j(1))/rhok_j - (Xl_ECF(1)-X_j(1))/rhol_j  ...
                    (Xk_ECF(2)-X_j(2))/rhok_j - (Xl_ECF(2)-X_j(2))/rhol_j ...
                    (Xk_ECF(3)-X_j(3))/rhok_j - (Xl_ECF(3)-X_j(3))/rhol_j];
                A(tt,:) = A0;
                A(m1+tt,:) = A0;
                %time_corr = (torr1-torr2-torr3-torr4);
                % Composing the right side

                % Tropospheric correction using standard meteorological parameters
                %[ ~,el_ki, ~] = topocent(X_i(1:3),Xk_ECF-X_i(1:3));
                %[ ~,el_li,~] = topocent(X_i(1:3),Xl_ECF-X_i(1:3));
                %[ ~,el_kj, ~] = topocent(X_j(1:3),Xk_ECF-X_j(1:3));
                %[az,el_lj,d] = topocent(X_j(1:3),Xl_ECF-X_j(1:3));
                %%el_ki,    el_li,    el_kj,    el_lj
                %t_corr = tropo(sin(el_lj*pi/180),...
                %    h_j*1.e-3,1013,293,50,0,0,0)...
                %    -tropo(sin(el_li*pi/180),....
                %    h_i*1.e-3,1013,293,50,0,0,0)...
                %   -tropo(sin(el_kj*pi/180),...
                %   h_j*1.e-3,1013,293,50,0,0,0)...
                %   +tropo(sin(el_ki*pi/180),...
                %   h_i*1.e-3,1013,293,50,0,0,0);
                t_corr = 0;
                b(tt,:) = (Obs1(refSV,2)-Obs1(t,2)-Obs2(refSV,2)+Obs2(t,2)  ...
                    -a(tt,1))*lambda1-t_corr;
                b(m1+tt,:)  = (Obs1(refSV,4)-Obs1(t,4)-Obs2(refSV,4)+Obs2(t,4) ...
                    -a(m1+tt,1))*lambda2-t_corr;
                bk(tt,:) =         ( rhok_i-rhok_j-rhol_i+rhol_j);
                bk(m1+tt,:) =  (rhok_i-rhok_j-rhol_i+rhol_j);
            end; % t

            % Kalman filter, see Table 8.2 on page 223 in Borre & Strang (2012):
            % Algorithms for Global Positioning, Wellesley-Cambridge Press
            P = P+Q;
            K = P*A'*inv(A*P*A'+R);
            x = x+ K*(b-A*x);
            P = (eye(3)-K*A)*P;
            b_old = b;
            b_acc = [b_acc b];
            x_acc = [x_acc x];
            X_j = X_i(1:3,:)+x;
    end % switch
    if feof(fid1) == 1, break, end;
    if feof(fid2) == 1, break; end;
end % while

% Transformation of geocentric baseline coordinates into topocentric coordinates
e = zeros(1,epoch-11);
n = zeros(1,epoch-11);
u = zeros(1,epoch-11);
for i = 1:epoch-11
    [e(i),n(i),u(i)] = xyz2enu(phi_j,lambda_j,x_acc(1,i),x_acc(2,i),x_acc(3,i));
end

fprintf('\n\nBaseline Components\n')
fprintf('\nX: %8.3f m,  Y: %8.3f m,  Z: %8.3f m\n', ...
    x_acc(1,end), x_acc(2,end), x_acc(3,end))
fprintf('\nE: %8.3f m,  N: %8.3f m,  U: %8.3f m\n',mean(e,2),mean(n,2),mean(u,2))

figure(1);
% we subtract e-, n-, and u-values for epoch 15 to start close to zero!
plot([(e-e(15))' (n-n(15))' (u-u(15))' ]*1000,'linewidth', .1)
title('Baseline Components from Kalman Filter')
ylabel('Baseline Relative to Initial Epoch [mm]','fontsize',14)
xlabel('Epochs [1 s interval]','fontsize',14)
legend('East','North','Up')
set(gca,'fontsize',14)
legend

print -dpdf easy61

figure(2);
plot(x_acc'*1000)
title('Baseline components  {\it(X}, {\itY}, {\itZ})')

ylabel('Absolute components [mm]','fontsize',14)

print -dpdf easy62

%%%%%%%%%% end easy6.m  %%%%%%%

author: yandld


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

%%
clear; close all; clc;

%% Syntax
% |[varargout]=loftLinQuad2hex(Fq,Vq,Vq2,numSteps);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
switch dataName
    case 'circle'
        [data4, data] = curves_with_angles(300);

         [data_noise, noise, temp_max_index] = addnoise(data, 0.4);
        feats = data_noise;

    case  'square'
          [data]=square(200);
           [data_noise, noise, temp_max_index] = addnoise(data, 0.05);
           feats=data_noise;
    case 'hole'
        [data] = generate.Generate_Ground_Truth(dataName,1000);
        [data_noise, noise, temp_max_index] = generate.addnoise(data, 0.1);
        feats = data_noise;
    case 'fishbowl'
        [data] = generate.Generate_Ground_Truth(dataName,600);
        [data_noise, noise, temp_max_index] = generate.addnoise(data, 0.1);%0.1
        feats = data_noise;
        [feats]=double(outlier_generation(feats));
       case 'sphere_helix'
        [data]=generate.Generate_Ground_Truth(dataName);
        [data_noise, noise, temp_max_index] = generate.addnoise(data, 0.4);
        feats = data_noise;

     case 'box_plane'
        box_plane;
        feats=data;
    case 'Hemisphere_plane'
        [data,data_clean]=Hemisphere_plane;
        feats=data;

    case 'planes'
    [feats,data]=pctest5 ();

case 'cyclo'
        load('test_cyclo.mat');
        data=X_iso;
        data=data(:,1:6000);
       var_noise=0.2;
       [data_noise,noise,temp_max_index] = generate.addnoise(data, var_noise);
       feats=data_noise;

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [samecost, sameassignment] = testassignment
%TESTASSIGNMENT  Test and compare assignment algorithms.
%		[SAMECOST, SAMEASSIGNMENT] = TESTASSIGN randomly generates distance
%		matrices and solves the assignment problem using different algorithms.
%		Edit the header of this file to change the simulation parameters.
%
%		assignment.html  File Exchange  Donate via PayPal
%
%		Markus Buehren
%		Last modified 05.07.2011
%
%		See also ASSIGNMENTOPTIMAL, ASSIGNMENTSUBOPTIMAL1,
%		ASSIGNMENTSUBOPTIMAL2, ASSIGNMENTALLPOSSIBLE.

% simulation time in seconds
testTime = 20;

% maximum matrix dimensions
maxOrders = [15, 25];

% If infAllowed in set to false, only distance matrices without infinite
%	costs are used
infAllowed = true;

% If the product of the dimensions of the distance matrix is smaller than
% maxDimProduct, the assignment function computing all possible assignments
% is used as reference. Set this to inf to use always or to 0 to never use.
maxDimProduct = 30;

% set recurstion limit
recursionLimit = get(0,'RecursionLimit');
set(0, 'RecursionLimit', 1000);

% use profiler or not
useProfiler = true;

% start profiler
if useProfiler
	profile clear
	profile on
end

% initialize
startTime = clock;
nassign = 0;
h = waitbar(0, 'Please wait', 'Name', mfilename);

while 1
	for dim1 = 1:maxOrders(1)
		for dim2 = 1:maxOrders(2)

			if ~infAllowed

				% generate distMatrix without infinite elements
				distMatrix = rand(dim1,dim2);

			else

				if rand(1) < 0.5
					% generate distMatrix with some infinite elements
					distMatrix = rand(dim1,dim2);

					if rand(1) < 0.5
						% set some elements to inf
						distMatrix(rand(dim1,dim2) > rand(1)) = inf;
					end

				else

					% generate distMatrix with many infinite elements
					distMatrix = repmat(inf, dim1, dim2);

					if rand(1) < 0.5
						for row=1:dim1
							if rand(1) < 0.8
								% set one element per row to finite number
								distMatrix(row, 1+floor(dim2*rand(1))) = rand(1);
							end
							if rand(1) < 0.3
								% set another element per row to finite number
								distMatrix(row, 1+floor(dim2*rand(1))) = rand(1);
							end
						end
					else
						for col=1:dim2
							if rand(1) < 0.8
								% set one element per column to finite number
								distMatrix(1+floor(dim1*rand(1)), col) = rand(1);
							end
							if rand(1) < 0.3
								% set another element per column to finite number
								distMatrix(1+floor(dim1*rand(1)), col) = rand(1);
							end
						end
					end

				end
			end

			% quantize distMatrix
			if rand(1) < 0.3
				distMatrix = round(max(10, 1000*rand(1)^2)*distMatrix);
			end

			% transpose distMatrix
			if rand(1) < 0.5
				distMatrix = distMatrix';
			end

			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
			[assignmentCell{1}, costCell{1}] = assignmentoptimal    (distMatrix); %#ok
			[assignmentCell{2}, costCell{2}] = assignmentsuboptimal1(distMatrix); %#ok
			[assignmentCell{3}, costCell{3}] = assignmentsuboptimal2(distMatrix); %#ok
			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

			% if dimensions are moderate, compute all possible solutions
			if dim1 * dim2 < maxDimProduct
				[assignment, cost] = assignmentallpossible(distMatrix);
			else
				assignment = assignmentCell{1};
				cost       = costCell{1};
			end

			% count trial runs
			nassign = nassign + 1;

			if ~exist('samecost', 'var')
				M = length(costCell);
				samecost        = zeros(M,1);
				sameassignment = zeros(M,1);
			end

			% compute penalty for non-assignments
			finiteIndex = isfinite(distMatrix);
			penalty = max(max(distMatrix(finiteIndex))) * dim1 * dim2;
			if ~isempty(penalty)
				cost = cost + length(find(~assignment)) * penalty;
			end

			for m=1:M

				% penalize non-assignments
				if ~isempty(penalty)
					costCell{m} = costCell{m} + length(find(~assignmentCell{m})) * penalty; %#ok
				end

				% compare costs
				if costCell{m} <= cost + 100*eps
					samecost(m) = samecost(m) + 1;
				end

				% compare assignments
				if all(assignmentCell{m} == assignment)
					sameassignment(m) = sameassignment(m) + 1;
				end

			end
		end
		% set waitbar
		curTime = etime(clock, startTime);
		waitbar(curTime/testTime, h);
	end

	% stop after given time
	if curTime > testTime
		break
	end
end
close(h);

% scale counters
samecost       = samecost/nassign;
sameassignment = sameassignment/nassign;

if useProfiler
	profile off
	profile report
end

% reset recurstion limit
set(0, 'RecursionLimit', recursionLimit);

author: yuxng


text:
%kkdbmupdate 'print, delete, or update (ascii) DBM key value'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kdbmupdate.pane file
%
% Parameters:
% InputFile: dbm 'DBM file', required: 'DBM file'
% String: key 'Key Name', default: ' ': 'DBM key'
% Toggle: printval 'Print Value', default: 0: 'print value of key'
% Toggle: delete 'Delete Key', default: 0: 'delete key'
% String: update 'Update Value', default: ' ': 'update value'
%
% Example: kkdbmupdate(dbm, {'dbm','';'key',' ';'printval',0;'delete',0;'update',' '})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kdbmupdate - print, delete, or update (ascii) DBM key value
%
%  DESCRIPTION
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

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

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function triangulation_test15 ( )

%*****************************************************************************80
%
%% TEST15 tests TRIANGULATION_ORDER3_BOUNDARY_NODE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 August 2006
%
%  Author:
%
%    John Burkardt
%
  dim_num = 2;
  node_num = 36;
  triangle_num = 41;

  triangle_node(1:3,1:triangle_num) = [ ...
     1,  8,  7; ...
     1,  2,  8; ...
     2,  9,  8; ...
     2,  3,  9; ...
     3, 10,  9; ...
     3,  4, 10; ...
     4, 11, 10; ...
     4,  5, 11; ...
     5, 12, 11; ...
     5,  6, 12; ...
     7, 14, 13; ...
     7,  8, 14; ...
     8, 15, 14; ...
     8,  9, 15; ...
    11, 18, 17; ...
    11, 12, 18; ...
    13, 20, 19; ...
    13, 14, 20; ...
    14, 21, 20; ...
    14, 15, 21; ...
    15, 22, 21; ...
    15, 16, 22; ...
    16, 23, 22; ...
    16, 17, 23; ...
    17, 24, 23; ...
    17, 18, 24; ...
    19, 26, 25; ...
    19, 20, 26; ...
    21, 28, 27; ...
    21, 22, 28; ...
    25, 30, 29; ...
    25, 26, 30; ...
    26, 31, 30; ...
    27, 32, 31; ...
    27, 28, 32; ...
    29, 34, 33; ...
    29, 30, 34; ...
    30, 35, 34; ...
    30, 31, 35; ...
    31, 36, 35; ...
    31, 32, 36 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST15\n' );
  fprintf ( 1, '  TRIANGULATION_ORDER3_BOUNDARY_NODE determines which\n' );
  fprintf ( 1, '  nodes lie on the boundary of a triangulation.\n' );

  node_boundary = triangulation_order3_boundary_node ( node_num, ...
    triangle_num, triangle_node );

  lvec_print ( node_num, node_boundary, '    Node  BN?' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [kN] = dyn2kN(dyn)
% Convert force from dyne to kilonewtons.
% Chad A. Greene 2012
kN = dyn*1e-8;
%
% C'mon, help me with this conversion--I'm dyne over here!

author: Sable


text:
function [llh] = tapas_mixedlinear_llh(data, theta, ptheta)
%% Likelihood of the nodes of a linear multivariate model with diagonal
% precisions on each parameter.
%
% It is the likelihood of
%
% p(y|x, beta, pi) = N(y - x*b, 1/pi)
%
% Note that this is a multivariate regression problem, so beta is a matrix
% in which each column correspond to one parameter of the model.
%

%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

y = data.y;
u = data.u;

x = data.u.x;
theta = theta.y;

[np, nc] = size(y);

llh = zeros(np, nc);

ln2pi = log(2 * pi);

for j = 1:nc
    % Mean
    beta = theta{j}.mu;
    % Precision
    pe = theta{j}.pe;
    % Regressors

    lpe = log(pe);

    % number of subjects x number of parameters
    py = u.x * beta;
    [ins, inp] = size(py);
    ty = y{j};

    % Residuals
    r = (reshape(ty, inp, ins)' - py);

    % Prediction error
    e = sum(bsxfun(@times, pe, r .* r), 2);

    % This is a vector. For each subject needs to include the precision.
    llh(:, j) = sum(-0.5 * ln2pi +  0.5 * lpe) - 0.5 * e;
end

end

author: translationalneuromodeling

repo: tapas


text:
function [map_osm] = load_osm_xml(filename)
% load OpenStreetMap XML file contents as a MATLAB structure
%
% download OpenStreetMap XML Data file to be parsed from:
%   http://www.openstreetmap.org/
%
% dependency
%   xml2struct, File Exchange ID = 28518, (c) 2010 by Wouter Falkena
%   http://www.mathworks.com/matlabcentral/fileexchange/28518-xml2struct
%
% note
%   xml2struct renamed to xml2struct_fex28518 to avoid name conflicts
%   with matlab bioinformatics toolbox
%
% 2010.11.06 (c) Ioannis Filippidis, jfilippidis@gmail.com
%
% See also PARSE_OPENSTREETMAP, PARSE_OSM.

if ~ischar(filename)
    error('Filename should be a string.')
end

if ~strfind(filename, '.osm')
    warning('Filename does not have the extension .osm')
end

map_osm = xml2struct_fex28518(filename); % downloaded osm file

author: johnyf

repo: openstreetmap


text:
function [x, mu, sigma] = zscore(x)
    mu=mean(x);
    sigma=max(std(x),eps);
	x=bsxfun(@minus,x,mu);
	x=bsxfun(@rdivide,x,sigma);
end

author: rasmusbergpalm

repo: DeepLearnToolbox


text:
function net = append_blocks(net, method, objectivetype, errortype, imchs)
% APPEND_BLOCKS function adds necessary blocks to the CNN.
%
%   NET = append_blocks(NET, METHOD, OBJECTIVE, ERROR, IMCHANNELS)
%
%   METHOD is a string that can contain one of the pooling strings:
%     'mac'   - global max pooling layer
%     'spoc'  - global average pooling layer
%     'gem'   - global generalized mean pooling layer with one p for all conv filters
%     'gemmp' - global generalized mean pooling layer with multiple p for each conv filter
%	If METHOD contains string 'edgefilter', the network will be pre-pended with the EdgeFilter layer.
%
%   OBJECTIVE & ERROR are a cell array with 2 values:
%     {'contrastiveloss', M} - contrastive loss with margin M
%     {'tripletloss', M}     - triplet loss with margin M
%     {'batchmap'}           - mean average precision for all queries in a batch
%
%   IMCHANNELS is a number of image channels that will be used for training, ie. 1 or 3.
%     If IMCHANNELS == 1 channels of first conv layers are aggregated using average pooling.
%
% Authors: F. Radenovic, G. Tolias, O. Chum. 2017.

	import dagnn.*;

	if ~exist('imchs'), imchs = 3; end
	if imchs == 1 	% sum aggregate first conv layer, if 1 input channel
		net.params(1).value = sum(net.params(1).value,3);
		net.layers(1).block.size(3) = 1;
	end
	method = lower(method);

	% Edge filtering layer
	if strfind(method, 'edgefilter')
		inputs = {'input'};
		outputs = {'x0'};
		params = struct('name', {}, 'value', {}, 'learningRate', [], 'weightDecay', []);
		name = 'edgelayer';
		block = EdgeFilter();

		params(1).name  = 'edgelayer_w'; params(1).value = single(10);
		params(2).name  = 'edgelayer_p'; params(2).value = single(0.5);
		params(3).name  = 'edgelayer_s'; params(3).value = single(500);
		params(4).name  = 'edgelayer_t'; params(4).value = single(0.1);

		net.addLayer(name, block, inputs, outputs, {params.name});

		% edgelayer_w params.value
		i = 1; f = net.getParamIndex(params(i).name);
		net.params(f).value = params(i).value; net.params(f).learningRate = 0; net.params(f).weightDecay = 0;
		% edgelayer_p params.value
		i = 2; f = net.getParamIndex(params(i).name);
		net.params(f).value = params(i).value; net.params(f).learningRate = 1; net.params(f).weightDecay = 0;
		% edgelayer_s params.value
		i = 3; f = net.getParamIndex(params(i).name);
		net.params(f).value = params(i).value; net.params(f).learningRate = 0; net.params(f).weightDecay = 0;
		% edgelayer_t params.value
		i = 4; f = net.getParamIndex(params(i).name);
		net.params(f).value = params(i).value; net.params(f).learningRate = 0.1; net.params(f).weightDecay = 0;

		net.layers(1).inputs{1} = outputs{1};
	end

	l = 0;

	% Global pooling layer
	l = l + 1;
	inputs = {sprintf('xx%d',l-1)};
	outputs = {'pooldescriptor'};
	params = struct('name', {}, 'value', {}, 'learningRate', [], 'weightDecay', []);
	name = 'pooldescriptor' ;
	if strfind(method, 'mac')
		block = MAC();
	elseif strfind(method, 'spoc')
		block = SPOC();
	elseif strfind(method, 'gem')
		block = GeM();
		if strfind(method, 'gemmp'), sz = [1 1 D 1]; else, sz = [1 1 1 1]; end
		params(1).name = 'gem';
		params(1).value = ones(sz, 'single') * 3;
		params(1).learningRate = 10;
		params(1).weightDecay  = 0;
	else
		error(sprintf('Unknown method %s\n', method));
	end
	net.addLayer(name, block, inputs, outputs, {params.name});
	if numel(strfind(method, 'gem'))
		findex = net.getParamIndex(params(1).name) ;
		if ~isempty(params(1).value)
			net.params(findex).value = params(1).value ;
			net.params(findex).learningRate = params(1).learningRate;
			net.params(findex).weightDecay = params(1).weightDecay;
		end
	end

	% L2 normalization layer
	inputs = {'pooldescriptor'};
	outputs = {'l2descriptor'};
	params = struct('name', {}, 'value', {}, 'learningRate', [], 'weightDecay', []);
	name = 'l2descriptor';
	block = L2N();
	net.addLayer(name, block, inputs, outputs, {params.name});

    % Objective layer
	inputs{1} = 'l2descriptor';
	inputs{2} = 'label';
	outputs = {'objective'};
	params = struct('name', {}, 'value', {}, 'learningRate', [], 'weightDecay', []);
	name = 'objective';
	switch objectivetype{1}
		case 'contrastiveloss'
			block = ContrastiveLoss('margin', objectivetype{2});
		case 'tripletloss'
			block = TripletLoss('margin', objectivetype{2});
		otherwise
			ferror('Unknow loss type\n');
	end
	net.addLayer(name, block, inputs, outputs, {params.name});

    % Error layer
	if exist('errortype') & numel(errortype)
		inputs{1} = 'l2descriptor';
		inputs{2} = 'label';
		outputs = {'error'};
		params = struct('name', {}, 'value', {}, 'learningRate', [], 'weightDecay', []);
		name = 'error';
		switch errortype{1}
			case 'contrastiveloss'
				block = ContrastiveLoss('margin', objectivetype{2});
			case 'tripletloss'
				block = TripletLoss('margin', objectivetype{2});
			case 'batchmap'
				block = BatchMAP();
			otherwise
				ferror('Unknow error type\n');
		end
		net.addLayer(name, block, inputs, outputs, {params.name});
	end

author: filipradenovic

repo: cnnimageretrieval


text:
function [x,fs,txt] = sapisynth(t,m)
%SAPISYNTH  synthesize text string or matrix [X,FS,TXT]=(T,M)
%
%  Usage:         sapisynth('Hello world');          % Speak text
%                 sapisynth([1 2+3i; -1i 4],'j');    % speak a matrix using 'j' for sqrt(-1)
%          [x,fs]=sapisynth('Hello world','k11');    % save waveform at 11kHz
%                 sapisynth('Hello world','fc');     % use a female child voice if available
%
%  Inputs: t  is either a text string or a matrix
%          m  is a mode string containing one or more of the
%             following options (# denotes an integer):
%
%             'l'   x=a cell array containing a list of talkers
%             'l#'  specify talker # (in the range 1:nvoices)
%             'r#'  speaking rate -10(slow) to +10(fast) [0]
%             'k#'  target sample rate in kHz [22]
%             'o'   audio output [default if no output arguments]
%             'j'   use 'j' rather than 'i' for complex numbers
%             'm','f' 'c','t','a','s' = Male Female Child, Teen, Adult, Senior
%                       specify any combination in order of priority
%             'v'   autoscale volumne to a peak value of +-1
%             'v#'  set volume (0 to 100) [100]
%             'p#'  set pitch -10 to +10 [0]
%
% Outputs: x    is the output waveform unless the 'l' option is chosen in
%               which case x is a cell array with one row per available
%               voice containing {'Name' 'Gender' 'Age'} where
%               Gender={Male,Female} and Age={Unknown,Child,Teen,Adult,Senior}
%          fs   is the actual sample frequency
%         txt   gives the actual text sring sent to the synthesiser
%
% The input text string can contain embedded command which are described
% in full at http://msdn.microsoft.com/en-us/library/ms717077(v=vs.85).aspx
% and summarised here:
%
% '... <bookmark mark="xyz"/> ...'               insert a bookmark
% '... <context id="date_mdy"> 03/04/01 </context> ...' specify order of dates
% '... <emph> ... </emph> ...'                   emphasise
% '... <volume level="50"> ... </volume> ...'    change volume level to 50% of full
% '... <partofsp part="noun"> XXX </partofsp> ...'      specify part of speech of XXX: unkown, noun, verb modifier, function, interjection
% '... <pitch absmiddle="-5"> ... </pitch> ...'  change pitch to -5 [0 is default pitch]
% '... <pitch middle="5"> ... </pitch> ...'      add 5 onto the pitch
% '... <pron sym="h eh 1 l ow "/> ...'           insert phoneme string
% '... <rate absspeed="-5"> ... </rate> ...'     change rate to -5 [0 is default rate]
% '... <rate speed="5"> ... </rate> ...'         add 5 onto the rate
% '... <silence msec="500"/> ...'                insert 500 ms of silence
% '... <spell> ... </spell> ...'                 spell out the words
% '... <voice required="Gender=Female;Age!=Child"> ...' specify target voice attributes to be Female non-child
%                                                         Age={Child, Teen, Adult, Senior}, Gender={Male, Female}
%
% Acknowledgement: This function was originally based on tts.m written by Siyi Deng

% Bugs/Suggestions:
%  (1) Allow the speaking of structures and cells
%  (2) Allow a blocking call to sound output and/or a callback procedure and/or a status call
%  (3) Have pitch and/or volume change to emphasise the first entry in a matrix row.
%  (4) extract true frequency from output stream

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

% Check that we are on a PC

if ~ispc, error('only works on a PC'); end

% decode the options

if nargin<2
    m='';
end
opts=zeros(52,3); % [exists+number specified, value]
lmode=length(m);
i=1;
while i<=lmode
    if i<lmode  % read a following integer if it exists
        [v,nv,e,ni]=sscanf(m(i+1:end),'%d',1);
    else
        nv=0;
        ni=1;
    end
    k=1+double(lower(m(i)))-'a'+26*(m(i)<'a');
    if k>=1 && k<=52
        opts(k,1)=1+nv;
        if nv
            opts(k,2)=v;
        end
        opts(k,3)=i;  % save position in mode string
    end
    i=i+ni;
end

S=actxserver('SAPI.SpVoice');
V=invoke(S,'GetVoices');  % get a list of voices from the registry
nv=V.Count;
if isempty(vv) || size(vvi,1)~=nv
    vv=cell(nv,3);
    vvi=zeros(nv,6);
    ages={'Senior' 'Adult' 'Teen' 'Child'};
    for i=1:nv
        VI=V.Item(i-1);
        vv{i,1}=VI.GetDescription;
        vv{i,2}=VI.GetAttribute('Gender');
        vvi(i,1)=MatchesAttributes(VI,'Gender=Male');
        vvi(i,2)=MatchesAttributes(VI,'Gender=Female');
        vv{i,3}='Unknown';
        for j=1:length(ages)
            if MatchesAttributes(VI,['Age=' ages{j}])
                vv{i,3}=ages{j};
                vvi(i,2+j)=1;
                break
            end
        end
    end
    vvj=vvi;
    % in the matrix below, the rows and columns are in the order Senior,Adult,Teen,Child.
    % Thus the first row gives the cost of selecting a voice with the wrong age when 'Senior'
    % was requested by the user. A voice of unkown age always scores 0 so entries with negative
    % values are preferred over 'unknown' while those with positive values are not.
    % Diagonal elements of the matrix are ignored (hence set to 0) since correct matches are
    % handled earlier with higher priority.
    vvj(:,3:6)=vvj(:,3:6)*[0 0 1 2; 1 0 2 3; 1 0 0 -1; 1 0 -1 0]'; % fuzzy voice attribute matching
end

% deal with the voice selection options

optv=opts([13 6 19 1 20 3],[3 1 2]);
if opts(12)   % if 'l' option specified - we need to get the voices
    if opts(12)>1
        S.Voice = V.Item(mod(opts(12,2)-1,nv));
    else
        x=vv;
        return
    end
elseif any(optv(:,2))
    optv(:,3)=(1:6)';
    optv=sortrows(optv(optv(:,2)>0,:));  % sort in order of occurrence in mode string
    no=size(optv,1);
    optp=zeros(nv,2*no+1);
    optp(:,end)=(1:nv)'; % lowest priority condition is original rank
    optp(:,1:no)=-vvi(:,optv(:,3));
    optp(:,no+1:2*no)=vvj(:,optv(:,3));
    optp=sortrows(optp);
    S.Voice = V.Item(optp(1,end)-1);
end

% deal with the 'r' option

if opts(18)>1  % 'r' option is specified with a number
    S.Rate=min(max(opts(18,2),-10),10);
end

% deal with the 'v' option

if opts(22)>1  % 'r' option is specified with a number
    S.Volume=min(max(opts(22,2),0),100);
end

% deal with the 'k' option

ff=[11025 12000 16000 22050 24000 32000 44100 48000]; % valid frequencies
if opts(11)>1  % 'k' option is specified with a number
    [v,jf]=min(abs(ff/1000-opts(11,2)));
else
    jf=4;  % default is 16kHz
end
fs=ff(jf);

M=actxserver('SAPI.SpMemoryStream');
M.Format.Type = sprintf('SAFT%dkHz16BitMono',fix(fs/1000));
S.AudioOutputStream = M;
if ischar(t)
    txt=t;
else
    txt='';
    if numel(t)
        sgns={' minus ', '', ' plus '};
        sz=size(t);
        w=permute(t,[2 1 3:numel(sz)]);
        sz(1:2)=sz(1)+sz(2)-sz(1:2); % Permute the first two dimensions for reading
        szp=cumprod(sz);
        imch='i'+(opts(10)>0);
        vsep='';
        for i=1:numel(w)
            wr=real(w(i));
            wi=imag(w(i));
            switch((wr~=0)+2*(wi~=0))+4*(abs(wi)==1)
                case {0,1}
                    txt=[txt sprintf('%s%g',vsep,wr)];
                case 2
                    txt=[txt sprintf('%s%g%c,',vsep,wi,imch)];
                case 3
                    txt=[txt sprintf('%s%g%s%g%c,',vsep,wr,sgns{2+sign(wi)},abs(wi),imch)];
                case 6
                    if wi>0
                        txt=[txt vsep imch ','];
                    else
                        txt=[txt vsep 'minus ' imch ','];
                    end
                case 7
                    txt=[txt sprintf('%s%g%s%c,',vsep,wr,sgns{2+sign(wi)},imch)];
            end
            % could use a <silence msec="???"/> command here
            vsep=[repmat('; ',1,find([0 mod(i,szp)]==0,1,'last')-1) ' '];
        end
    end
end

% deal with the 'p' option

if opts(16)>1  % 'r' option is specified with a number
    txt=[sprintf('<pitch absmiddle="%d"> ',min(max(opts(16,2),-10),10)) txt];
end

invoke(S,'Speak',txt);
x = mod(32768+reshape(double(invoke(M,'GetData')),2,[])'*[1; 256],65536)/32768-1;
delete(M);      % delete output stream
delete(S);      % delete all interfaces

if opts(22)==1 % 'v' option with no argument
    x=x*(1/max(abs(x))); % autoscale
end
if opts(15)>0 || ~nargout % 'o' option for audio output
    sound(x,fs);
end

author: decouples


text:
function [acc_stl] = STL(X_src,y_src,X_tar,y_tar,dim)
% Stratified Transfer Learning for Cross-domain Activity Recognition.
% Jindong Wang, Yiqiang Chen, Lisha hu, Xiaohui Peng, and Philip S. Yu.
% 2018 IEEE International Conference on Pervasive Computing and Communications (PerCom).

%% Inputs:
%%%% X_src : Source domain feature matrix, n_sample * n_dim
%%%% y_src : Source domain label vector, n_sample * 1
%%%% X_tar : Target domain feature matrix, n_sample * n_dim
%%%% y_src : Target domain label vector, n_sample * 1
%%%% dim   : Dimension after STL, an integer
%% Output:
%%%% acc_stl: accuracy of the algorithm

    label_set = unique(y_src(:));
    num_tree = 30;
    %voting
    %%RF
    [predict_label_rf,accuracy_rf] = jdrf(X_src,y_src,X_tar,y_tar,30);
    %%SVM
    [predict_label_svm,accuracy_svm,~] = jdsvm(X_src,y_src,X_tar,y_tar,100);
    %%KNN
    [predict_label_knn,accuracy_knn] = jdknn(X_src,y_src,X_tar,y_tar,5);
    %%voting
    label_voting = [];
    for i = 1 : size(y_tar,1)
        if predict_label_rf(i) == predict_label_svm(i)
            label_voting(i) = predict_label_rf(i);
        elseif predict_label_rf(i) == predict_label_knn(i)
            label_voting(i) = predict_label_rf(i);
        elseif predict_label_svm(i) == predict_label_knn(i)
            label_voting(i) = predict_label_svm(i);
        else
            label_voting(i) = -1;
        end
    end
    label_voting = label_voting';
    %get residual
    r_index = label_voting(:,1) == -1;
    X_tar_residual = X_tar(r_index,:);
    y_tar_residual = y_tar(r_index,:);

    %get candidate
    c_index = label_voting(:,1) ~= -1;
    X_tar_candidate = X_tar(c_index,:);
    y_tar_candidate = y_tar(c_index,:);
    label_candidate = label_voting(c_index);
    conmat1 = confusionmat(label_candidate,y_tar_candidate);
    %in-class transfer
    %%%new candidates and source domain
    X_tar_candidate_new = [];
    y_tar_candidate_new = [];
    X_src_new = [];
    y_src_new = [];
    X_tar_candidate_old = [];
    y_tar_candidate_old = [];
    for i = 1 : size(label_set,1)
        class_index = label_set(i);
        %source data in that class i
        src_index_i = y_src(:) == class_index;
        X_src_i = X_src(src_index_i,:);
        y_src_i = y_src(src_index_i,:);
        %candidate data in that class i
        tar_index_i = label_candidate(:) == class_index;
        X_tar_candidate_i = X_tar_candidate(tar_index_i,:);
        y_tar_candidate_i = y_tar_candidate(tar_index_i,:);
        X_tar_candidate_old = [X_tar_candidate_old;X_tar_candidate_i];
        y_tar_candidate_old = [y_tar_candidate_old;y_tar_candidate_i];
        %TCA
        [X_src_new_i,X_tar_candidate_new_i] = MyTCA(X_src_i,X_tar_candidate_i,dim);
        %merge
        X_src_new = [X_src_new;X_src_new_i];
        y_src_new = [y_src_new;y_src_i];
        X_tar_candidate_new = [X_tar_candidate_new;X_tar_candidate_new_i];
        y_tar_candidate_new = [y_tar_candidate_new;y_tar_candidate_i];
    end
    %train model on new candidate
    [predicted_label_candidate,acc_candidate] = ...
        jdrf(X_src_new,y_src_new,X_tar_candidate_new,y_tar_candidate_new,num_tree);
    conmat2 = confusionmat(predicted_label_candidate,y_tar_candidate_new);
    %train model on old candidate and predict on residual
    [predicted_label_residual,acc_residual] = jdrf(X_tar_candidate_old,predicted_label_candidate,X_tar_residual,y_tar_residual,num_tree);
    %calculate accuracy
    acc_stl = mean([predicted_label_candidate;predicted_label_residual] == [y_tar_candidate_new;y_tar_residual]);

end

author: jindongwang

repo: activityrecognition