text:
%kMulWave 'multiplies each element with exp(-ikx) '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros MulWave.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'First Input data object'
% InputFile: OrderFactors 'Order Factors', optional: 'a matrix which supplies float value order factors (positive or negative)'
% Integer: nx 'OrdersX', default: 2: 'number of orders along x (0 .. +- OrdersX)'
% Integer: ny 'OrdersY', default: 0: 'number of orders along y (0 .. +- OrdersY)'
% Toggle: auto 'Determine Grid automatically', default: 0: 'if selected, the grid values are ignored and the grid position, direction and shift is determined from the first order in X and Y'
% Integer: aorderx 'UseOrder X #', default: 1: 'which order to use for shift determination'
% Integer: aordery 'Y', default: 1: 'which y order to use for shift determination'
% Double: angle 'Angle', default: 0: 'angle, with which the pattern is oriented to the X-Axis'
% Double: periodx 'Grid Spacing X', default: 4: 'number of pixels from grid to grid element'
% Double: periody 'Grid Spacing Y', default: 4: 'number of pixels from grid to grid element'
% Toggle: autom 'Determine Gridshift automatically', default: 0: 'The grid zero position will be determined automatically'
% Double: mx 'Grid shift X [rad]', default: 0: 'exp(-i val) will be applied to first order peak in X'
% Double: my 'Y [rad]', default: 0: 'exp(-i val) will be applied to first order peak in Y'
% Double: of 'Order Enhancement Factor', default: 1: 'Will determine the enhancement of higher orders in the sum'
% Double: zw 'Zero Order Weight', default: 1: 'if 1.0, the zero order will not have and ignore radius, if ==0.0 the zero order will be ignored in the center like all other orders'
% Double: ir 'Ignore Zero Radius', default: 0: 'The reconstructed zero peak in each order is ignored by weigthing with a certain gaussian radius'
% Double: ar 'Apodize Radius', default: 0: 'relative max. frequency to apodize to by a cos window (linear edge), see Gabor 1946. If this is zero, no apodization will be performed and the orders will be plainly added'
% Double: ary 'Apo Radius Y', default: 0.9: 'If given an anisotropic apodization will be performed'
% Double: kxmax 'Order Kmax X', default: 0.5: 'This defines the border of each reconstructed peak, as a relative frequency'
% Double: kymax 'Order Kmax Y', default: 0.5: 'This defines the border of each reconstructed peak, as a relative frequency'
% OutputFile: o 'Output', required: 'Resulting output data object'
% OutputFile: sum 'Sum', required: 'The sum of all elements (reconstructed Fourier-components) will be stored here'
% OutputFile: pat 'Patterns', optional: 'If given, patterns will be generated for each class of distances and saved (each class along depth), moved patterns along elements.'
% Integer: PStepsX 'Pattern Steps X', default: 3: 'Number of steps along X for pattern movement'
% Integer: PStepsY 'Y', default: 3: 'Number of steps along Y for pattern movement'
%
% Example: [o, sum, pat] = kMulWave({i, OrderFactors}, {'i','';'OrderFactors','';'nx',2;'ny',0;'auto',0;'aorderx',1;'aordery',1;'angle',0;'periodx',4;'periody',4;'autom',0;'mx',0;'my',0;'of',1;'zw',1;'ir',0;'ar',0;'ary',0.9;'kxmax',0.5;'kymax',0.5;'o','';'sum','';'pat','';'PStepsX',3;'PStepsY',3})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% MulWave - multiplies each element with exp(-ikx)
%
%  DESCRIPTION
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1996-2003, Rainer Heintzmann,  All rights reserved.
%

function varargout = kMulWave(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,..] = kMulWave(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'OrderFactors', '__input';'nx', 2;'ny', 0;'auto', 0;'aorderx', 1;'aordery', 1;'angle', 0;'periodx', 4;'periody', 4;'autom', 0;'mx', 0;'my', 0;'of', 1;'zw', 1;'ir', 0;'ar', 0;'ary', 0.9;'kxmax', 0.5;'kymax', 0.5;'o', '__output';'sum', '__output';'pat', '__output';'PStepsX', 3;'PStepsY', 3};
maxval={0,1,1,1,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,1,1};
minval={0,1,1,1,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1};
istoggle=[0,1,0,0,1,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,1,0,0];
was_set=istoggle * 0;
paramtype={'InputFile','InputFile','Integer','Integer','Toggle','Integer','Integer','Double','Double','Double','Toggle','Double','Double','Double','Double','Double','Double','Double','Double','Double','OutputFile','OutputFile','OutputFile','Integer','Integer'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=2; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
principal = 1000;
alpha = 0.1;
debt = principal * (1 + alpha) ^ 2;
disp(debt);

author: anishLearnsToCode

repo: introduction-to-programming-with-matlab


text:
function kern = rbfhKernExpandParam(kern, params)

% RBFHKERNEXPANDPARAM Create kernel structure from RBFH kernel's parameters.
% FORMAT
% DESC returns a radial basis function HEAT kernel structure filled with
% the parameters in the given vector. This is used as a helper function to
% enable parameters to be optimised in, for example, the NETLAB
% optimisation functions.
% ARG kern : the kernel structure in which the parameters are to be
% placed.
% ARG param : vector of parameters which are to be placed in the
% kernel structure.
% RETURN kern : kernel structure with the given parameters in the
% relevant locations.
%
% SEEALSO : rbfhKernParamInit, rbfhKernExtractParam, kernExpandParam
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

kern.inverseWidthTime = params(1);
kern.inverseWidthSpace = params(2);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function val=perm(A,boolRowsSkip,boolColsSkip)
%%PERM Calculate the matrix permanent allowing for rows and columns to be
%      skipped if desired (operate on a submatrix). The permanent is
%      equivalent to calculating the determininant in the standard
%      summation manner taught in school, except all of the minus signs are
%      replaced with plus signs. Permanents play a role in combinatorics
%      and the ability to skip rows and columns can be useful when using
%      permanents to compute target-measurement association probabilities.
%
%INPUTS: A An mXn matrix. If m<=n, then the standard matrix permanent is
%         found. If m>n, then the permanent of A' is found to be consistent
%         with the permanents of A and A' being equal in square matrices.
%         Empty matrices have a permanent of one by definition.
% boolRowsSkip An optional list of rows in A that should be skipped when
%         computing the matrix permanent. This is an mX1 or 1Xm boolean
%         vector where a 1 means that row should be skipped. If omitted or
%         an empty matrix is passed, no rows are skipped.
% boolColsSkip An optional list of columns in A that should be skipped when
%         computing the matrix permanent. This is an nX1 or 1Xn boolean
%         vector where a 1 means that column should be skipped. If omitted
%         or an empty matrix is passed, no columns are skipped.
%
%OUTPUTS: val The matrix permanent of A, omitting any rows or columns as
%             necessary.
%
%Whereas polynomial-time algorithms exist for calculating determinants, it
%was proven in [1] that matrix permanents can not be computed in polynomial
%time unless P=NP. However, efficient non-polynomial time algorithms exist.
%This file implements the method in theorem 4.1 on page 26 in Chapter 2 of
%[2] when dealing with general rectangular matrices. When dealing with
%square matrices, an algorithm based on PERMAN from Chapter 23 of [3] is
%used, because it is more efficient and appears to be less susceptible to
%finite precision errors. When skipping rows or columns, only Ryser's
%algorithm is used.
%
%EXAMPLE 1:
%Suppose that one is given a boolean matrix where 1 indicates that a target
%gates with a measurement and a 0 indicates that it does not gate. To deal
%with the possibility of missed detections, one often adds dummy
%measurements, one per target, each gating with only one target. The total
%number of possible target-measurement (and missed detection) assignments
%equals the permanent of the matrix.
% %For 5 targets, 5 measurements (everything gates) and no missed
% %detections
% perm(ones(5,5))
% %One gets 120, there are 120 possible assignments (factorial(5)).
% %For 5 targets target where everything gates with the possibility of
% %missed detections, one has
% perm([ones(5,5),eye(5)])
% %Leading to 1546 possible assignments to measurements and missed
% %detections.
%
%EXAMPLE 2:
%This is an example with skip lists.
% numRow=12;
% numCol=4;
% A=magic(max(numRow,numCol));
% A=A(1:numRow,1:numCol);
% boolRowsSkip=false(numRow,1);
% boolColsSkip=false(numCol,1);
% boolRowsSkip(3)=true;
% boolColsSkip([1,4])=true;
% perm(A,boolRowsSkip,boolColsSkip)
% perm(A',boolColsSkip,boolRowsSkip)
% %The permanent values should both be 494938.
%
%REFERENCES:
%[1] L. G. Valiant, "The complexity of computing the permanent,"
%    Theoretical Computer Science, vol. 8, no. 2, pp. 189-201, 1979.
%[2] H. J. Ryser, Combinatorial Mathematics, ser. The Carus Mathematical
%    Monographs. The Mathematical Association of America, 1963, no. 14.
%[3] A. Nijenhuis and H. S. Wilf, Combinatorial Algorithms for Computers
%    and Calculators, 2nd ed. New York: Academic press, 1978.
%
%October 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    if(isempty(A))
        %Empty matrices have a permanent of 1 by definition.
        val=1;
        return;
    elseif(numel(A)==1)
        %Scalar values have a permanent equal to that value.
        val=A;
        return;
    end

    m=size(A,1);
    n=size(A,2);

    %Use a modified version of Ryser's algorithm if skip lists are
    %provided.
    if(nargin>1)
        mTotal=m;
        nTotal=n;
        if(nargin<3||isempty(boolColsSkip))
            boolColsSkip=false(nTotal,1);
        end

        if(isempty(boolRowsSkip))
            boolRowsSkip=false(mTotal,1);
        end

        numRowsSkipped=sum(boolRowsSkip);
        numColsSkipped=sum(boolColsSkip);

        m=mTotal-numRowsSkipped;
        n=nTotal-numColsSkipped;

        %Empty matrices have a permanent of 1 by definition.
        if(isempty(A)||m==0||n==0)
            val=1;
            return;
        end

        if(m>n)
            A=A';
            temp=m;
            m=n;
            n=temp;

            temp=mTotal;
            mTotal=nTotal;
            nTotal=temp;

            temp=boolRowsSkip;
            boolRowsSkip=boolColsSkip;
            boolColsSkip=temp;
        end

        %Set the mapping of indices of the rows in the submatrix to indices in
        %the full matrix.
        rows2Keep=zeros(m,1);%Allocate space
        cumSkip=0;
        curIdx=1;
        for curRow=1:mTotal
            if(boolRowsSkip(curRow)==true)
                cumSkip=cumSkip+1;
                continue;
            else
                rows2Keep(curIdx)=curIdx+cumSkip;
                curIdx=curIdx+1;
            end
        end

        %Set the mapping of indices of the columns in the submatrix to indices in
        %the full matrix.
        cols2Keep=zeros(n,1);%Allocate space
        cumSkip=0;
        curIdx=1;
        for curCol=1:nTotal
            if(boolColsSkip(curCol)==true)
                cumSkip=cumSkip+1;
                continue;
            else
                cols2Keep(curIdx)=curIdx+cumSkip;
                curIdx=curIdx+1;
            end
        end

        binomTerm=1;
        val=0;
        for x=0:(m-1)
            val=val+SigmaSSkip(A,n-m+x,rows2Keep,cols2Keep)*binomTerm;

            binomTerm=binomTerm*(1-m+n+x)/(1+x)*(-1);
        end

        return;
    end

    if(m~=n)%Use Ryser's algorithm
        if(m>n)
            A=A';
            temp=m;
            m=n;
            n=temp;
        end

        binomTerm=1;
        val=0;
        for x=0:(m-1)
            val=val+SigmaS(A,n-m+x)*binomTerm;

            binomTerm=binomTerm*(1-m+n+x)/(1+x)*(-1);
        end
    else%If the matrix is square, use the PERMAN algorithm.
        x=zeros(n,1);%Temporary storage space.
        p=0;

        for i=1:n
            sumVal=sum(A(i,:));
            x(i)=A(i,n)-sumVal/2;
        end

        sgn=-1;
        code=[];
        nCard=n-1;
        while(1)
            sgn=-sgn;
            prodVal=sgn;
            [code,nCard,isLast,j]=getNextGrayCode(code,nCard);

            if(nCard~=0)
                z=2*code(j)-1;
                x=x+z*A(:,j);
            end
            for i=1:n
                prodVal=prodVal*x(i);
            end
            p=p+prodVal;

            if(isLast)
                break;
            end
        end
        val=2*(2*mod(n,2)-1)*p;
    end
end

function val=S(A)
%This function gives us the product of the row sums of A.
    val=prod(sum(A,2));
end

function val=SigmaS(A,r)
    %This adds up all of the possible values of S(A) where r columns of A
    %have been replaced by zeros. We shall choose the
    %n-r columns of A that are NOT zero.
    n=size(A,2);
    combLen=n-r;
    curComb=1:combLen;
    val=0;
    while(~isempty(curComb))
        val=val+S(A(:,curComb));
        curComb=getNextCombo(curComb,n,1);
    end
end

function val=SigmaSSkip(A,r,rows2Keep,cols2Keep)
    %This adds up all of the possible values of S(A) where r columns of A
    %have been replaced by zeros. We shall choose the
    %n-r columns of A that are NOT zero.
    n=length(cols2Keep);
    combLen=n-r;
    curComb=0:(combLen-1);

    val=0;
    while(~isempty(curComb))
        val=val+S(A(rows2Keep,cols2Keep(curComb+1)));
        curComb=getNextCombo(curComb,n);
    end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function metadata = psdreadheader(fid)
%PSDREADHEADER  Read the metadata from an open PSD file.
%   METADATA = PSDREADHEADER(FID) fread metadata from an Adobe Photoshop
%   PSD file using the file handle FID, which must already be opened.

%   Author: Jeff Mather
%   Copyright 1984-2004 The MathWorks, Inc.

fids = fopen('all');

if (isempty(find(fid == fids)))
    error('psdreadheader:badFID', 'Bad file ID.');
end

% Read the PSD header.
metadata.FormatSignature = fread(fid, 4, 'uint8=>char');
if (~isequal(metadata.FormatSignature', '8BPS'))
    fclose(fid);
    error('impsdinfo:badSignature', 'Format signature mismatch (%s).', ...
          metadata.FormatSignature)
end

metadata.FormatVersion = fread(fid, 1, 'uint16=>uint16');
if (metadata.FormatVersion ~= 1)
    error('impsdinfo:badVersion', 'Bad PSD version number (%d).', ...
          metadata.FormatVersion)
end

fseek(fid, 6, 'cof');  % Reserved
metadata.NumSamples = fread(fid, 1, 'uint16');
metadata.Rows = fread(fid, 1, 'uint32');
metadata.Columns = fread(fid, 1, 'uint32');
metadata.BitsPerSample = fread(fid, 1, 'uint16=>uint16');
metadata.ColorMode = decodeColorMode(fread(fid, 1, 'uint16=>uint16'));

% Read color mode data, image resources, and layer/mask information.
blockLen = fread(fid, 1, 'uint32');
metadata.ColorModeData.Length = blockLen;
if (blockLen > 0)
    metadata.ColorModeData.Data = fread(fid, blockLen, 'uint8=>uint8');
else
    metadata.ColorModeData.Data = [];
end

blockLen = fread(fid, 1, 'uint32');
metadata.ImageResources.Length = blockLen;
if (blockLen > 0)
    metadata.ImageResources.Data = fread(fid, blockLen, 'uint8=>uint8');
else
    metadata.ImageResources.Data = [];
end

blockLen = fread(fid, 1, 'uint32');
metadata.LayersAndMasks.Length = blockLen;
if (blockLen > 0)
    metadata.LayersAndMasks.Data = fread(fid, blockLen, 'uint8=>uint8');
else
    metadata.LayersAndMasks.Data = [];
end

% Read the image data block.
metadata.Compression = fread(fid, 1, 'uint16');

function modeString = decodeColorMode(modeNum)
%DECODECOLORMODE  Convert a numeric color mode to a string.

switch (modeNum)
case 0
    modeString = 'Bitmap';
case 1
    modeString = 'Grayscale';
case 2
    modeString = 'Indexed color';
case 3
    modeString = 'RGB';
case 4
    modeString = 'CMYK';
case 7
    modeString = 'Multichannel color';
case 8
    modeString = 'Duotone';
case 9
    modeString = 'Lab';
otherwise
    modeString = modeNum;
end

author: Sable

repo: mcbench-benchmarks


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
% z-Transform properties

%linearity

syms n z

x1=n^2;
x2=2^n;
a1=3;
a2=4;

Le=a1*x1+a2*x2;
Left=ztrans(Le,z)

X1=ztrans(x1);
X2=ztrans(x2);
Right=a1*X1+a2*X2


text:
function [y,N,algo] = erf_(x,N,algo)
% error function for a real or complex array
%
% Note: The Fresnel integrals C(x) and S(x) can be computed for real x as
%   C(x)+1i*S(x) == ((1+1i)/2)*erf_((sqrt(pi)*(1-1i)/2)*x)
% (See demo example at the end of this file.)
%
% syntax:
%   y = erf_(x);
%   [y,N,algo] = erf_(x,N,algo);
%
% inputs:
%
%   x: numeric array, can be complex
%
% optional inputs:
%
%   N: positive integer - scalar or size-matched to x, or [] (optional; unspecified or
%   [] defaults to 1000)
%   maximum number of terms to use in series or continued-fraction expansion.  N can be
%   specified independently for each x element
%
%   algo: integers 1, 2, or 3 - scalar or size-matched to x, or [] (optional; default
%   is [])
%   Set algo = 1 to force use of Taylor series expansion (good for small arguments).
%   Set algo = 2 to force use of erf continued-fraction expansion (good for arguments
%   with large real part).
%   Set algo = 3 to force use of Dawson's-function continued-fraction expansion (good
%   for arguments with large imaginary part).
%   Leave algo unspecified or [] to make the choice automatically.
%   algo can be specified independently for each x element.
%
% outputs:
%
%   y: numeric array, size-matched to x
%   y = 2/sqrt(pi) * integral from 0 to x of exp(-t^2) dt.
%
% optional outputs:
%
%   N: integer array, size-matched to x
%   number of terms used in series or continued-fraction expansion for each x element
%
%   algo: integer array of values 1, 2, or 3; size-matched to x
%   This indicates which algorithm was used for each x. 1: Taylor series, 2: erf
%   continued fraction, 3: Dawson's-function continued fraction
%
%
% Author: Kenneth C. Johnson, kjinnovation@earthlink.net, kjinnovation.com
% Version: November 1, 2011
%
%
% BSD Copyright notice:
%
% Copyright 2011 by Kenneth C. Johnson (kjinnovation.com)
%
% Redistribution and use in source and binary forms, with or without modification, are
% permitted provided that this entire copyright notice is duplicated in all such copies.
%
% This software is provided "as is" and without any express or implied warranties,
% including, without limitation, the implied warranties of merchantibility and fitness
% for any particular purpose.
%

if nargin<3
    algo = [];
    if nargin<2
        N = [];
    end
end

size_ = size(x);

if isempty(N)
    N = 1000;
end
if isscalar(N)
    N = repmat(N,size_);
elseif ~isequal(size(N),size_)
    error('erf_:validation','Size mismatch between args N and x.')
end

if isempty(algo)
    % Algorithm selection for erf_(x):
    % Apply algo 1 (Taylor series) for |real(x)|<=1, |imag(x)|<=6.
    % Apply algo 3 (Dawson's-function continued fraction) for |real(x)|<=1, |imag(x)|>6.
    % Apply algo 2 (erf continued fraction) for |real(x)|>1.
    % See test code at the bottom of this file for checking algorithm continuity across
    % algo boundaries.
    algo = ones(size(x));
    algo(abs(imag(x))>6) = 3;
    algo(abs(real(x))>1) = 2;
else
    if isscalar(algo)
        algo = repmat(algo,size_);
    elseif ~isequal(size(algo),size_)
        error('erf_:validation','Size mismatch between args algo and x.')
    end
    if ~all(algo==1 | algo==2 | algo==3)
        error('erf_:validation','Invalid arg: algo values must be 1, 2, or 3.')
    end
end

x = x(:);
N = N(:);
algo = algo(:);
y = zeros(size(x));
sgn = ones(size(x));
sgn(real(x)<0) = -1;
x = sgn.*x; % Only apply algorithm with abs(angle(x))<=pi/2; use erf(x)==-erf(-x).

warn = false; % for possible non-convergence

i = find(algo==1);
if ~isempty(i)
    % Apply series expansion to x(i) (see http://mathworld.wolfram.com/Erf.html, Eq 6).
    n = 0; % summation index
    Ni = N(i);
    term = x(i);
    xx = term.*term;
    term = (2/sqrt(pi))*term; % n-th summation term
    psum = term; % partial sum up to n-th term
    prec = abs(psum)*eps; % estimated precision of partial sum
    while true
        n = n+1;
        term = (-(2*n-1)/(n*(2*n+1)))*term.*xx;
        psum = psum+term;
        test = (abs(term)<=prec);
        if any(test)
            y(i(test)) = psum(test);
            N(i(test)) = n;
            i(test) = [];
            if isempty(i)
                break
            end
            Ni(test) = [];
            term(test) = [];
            psum(test) = [];
            prec(test) = [];
            xx(test) = [];
        end
        test = (n==Ni);
        if any(test)
            warn = true;
            y(i(test)) = psum(test);
            N(i(test)) = n;
            i(test) = [];
            if isempty(i)
                break
            end
            Ni(test) = [];
            term(test) = [];
            psum(test) = [];
            prec(test) = [];
            xx(test) = [];
        end
        prec = max(prec,abs(psum)*eps);
    end
end

i = find(algo==2);
if ~isempty(i)
    % Apply erf continued fraction to x(i).  Use Eq 6.9.4 (second form) in Numerical
    % Recipes in C++, 2nd Ed. (2002).  Apply Lentz's method (Sec. 5.2) to top-level
    % denominator in Eq 6.9.4.
    n = 0;
    Ni = N(i);
    b = x(i);
    b = 2*b.*b+1;
    f = b;
    C = f;
    D = zeros(size(i));
    while true
        n = n+1;
        a = -(2*n-1)*(2*n);
        b = b+4;
        D = b+a*D;
        D(D==0) = eps^2;
        D = 1./D;
        C = b+a./C;
        C(C==0) = eps^2;
        Delta = C.*D;
        f = f.*Delta;
        test = (abs(Delta-1)<=eps);
        if any(test)
            y(i(test)) = f(test);
            N(i(test)) = n;
            i(test) = [];
            if isempty(i)
                break
            end
            Ni(test) = [];
            b(test) = [];
            f(test) = [];
            C(test) = [];
            D(test) = [];
        end
        test = (n==Ni);
        if any(test)
            warn = true;
            y(i(test)) = f(test);
            N(i(test)) = n;
            i(test) = [];
            if isempty(i)
                break
            end
            Ni(test) = [];
            b(test) = [];
            f(test) = [];
            C(test) = [];
            D(test) = [];
        end
    end
    i = find(algo==2);
    y(i) = 1-exp(-x(i).^2).*((2/sqrt(pi))*x(i)./y(i));
end

i = find(algo==3);
if ~isempty(i)
    % Apply Dawson's-function (F) continued fraction to x(i).  See
    % http://mathworld.wolfram.com/DawsonsIntegral.html, Eq 14.
    % erf(x) = i*erfi(-i*x) = (2i/sqrt(pi))*exp(-x^2)*F(-i*x)
    % Apply Lentz's method to top-level denominator in Eq 14.
    n = 0;
    Ni = N(i);
    xx = x(i);
    xx = -xx.*xx; % (-1i*x)^2
    b = 1+2*xx;
    f = b;
    C = f;
    D = zeros(size(i));
    while true
        n = n+1;
        a = -4*n*xx;
        b = 2+b;
        D = b+a.*D;
        D(D==0) = eps^2;
        D = 1./D;
        C = b+a./C;
        C(C==0) = eps^2;
        Delta = C.*D;
        f = f.*Delta;
        test = (abs(Delta-1)<=eps);
        if any(test)
            y(i(test)) = f(test);
            N(i(test)) = n;
            i(test) = [];
            if isempty(i)
                break
            end
            Ni(test) = [];
            xx(test) = [];
            b(test) = [];
            f(test) = [];
            C(test) = [];
            D(test) = [];
        end
        test = (n==Ni);
        if any(test)
            warn = true;
            y(i(test)) = f(test);
            N(i(test)) = n;
            i(test) = [];
            if isempty(i)
                break
            end
            Ni(test) = [];
            xx(test) = [];
            b(test) = [];
            f(test) = [];
            C(test) = [];
            D(test) = [];
        end
    end
    i = find(algo==3);
    % Current state: y(i)==(-1i*x(i))./F(-1i*x(i)).  Convert F(-1i*x) to erfi(x).
    y(i) = (2/sqrt(pi))*exp(-x(i).^2).*x(i)./y(i);
end

y = sgn.*y;
y = reshape(y,size_);
N = reshape(N,size_);
algo = reshape(algo,size_);

if warn
    warning('erf_:convergence','Possible non-convergence in erf_');
end

return

%--------------------------------------------------------------------------------------
% Demo example: Cornu Spiral (http://en.wikipedia.org/wiki/Cornu_spiral)
FresnelCS = @(x) ((1+1i)/2)*erf_((sqrt(pi)*(1-1i)/2)*x); % Fresnel integrals
CS = FresnelCS(-10:.001:10);
figure, plot(real(CS),imag(CS)), axis equal

%--------------------------------------------------------------------------------------
% Test code for checking algorithm discontinuity across algo boundaries:

% boundary between algo 1 and algo 2:
close all
x = (0:.02:3).'; y = [0 6]; z = repmat(x,size(y))+1i*repmat(y,size(x));
[e,N] = erf_(z,[],1); [e_,N_] = erf_(z,[],2); d = log10(abs(2*(e-e_)./(e+e_)));
figure, hold on, plot(x,N), plot([1 1],ylim,':k')
xlabel('x'), ylabel('N, algo 1'), legend('y = 0','y = 6')
figure, hold on, plot(x,N_), plot([1 1],ylim,':k')
xlabel('x'), ylabel('N, algo 2'), legend('y = 0','y = 6')
figure, hold on, plot(x,d), plot([1 1],ylim,':k')
xlabel('x'), ylabel('log10(algo 1-2 diff)'), legend('y = 0','y = 6')
title('erf\_(x+1i*y), algo 1 vs 2 (boundary at |x|==1, |y|<=6)')

% boundary between algo 3 and algo 2:
close all
x = (0:.02:3).'; y = [6 10]; z = repmat(x,size(y))+1i*repmat(y,size(x));
[e,N] = erf_(z,[],3); [e_,N_] = erf_(z,[],2); d = log10(abs(2*(e-e_)./(e+e_)));
figure, hold on, plot(x,N), plot([1 1],ylim,':k')
xlabel('x'), ylabel('N, algo 3'), legend('y = 6','y = 10')
figure, hold on, plot(x,N_), plot([1 1],ylim,':k')
xlabel('x'), ylabel('N, algo 2'), legend('y = 6','y = 10')
figure, hold on, plot(x,d), plot([1 1],ylim,':k')
xlabel('x'), ylabel('log10(algo 3-2 diff)'), legend('y = 6','y = 10')
title('erf\_(x+1i*y), algo 3 vs 2 (boundary at |x|==1, |y|>6)')

% boundary between algo 1 and algo 3:
close all
y = (0:.02:10).'; x = [0 1]; z = repmat(x,size(y))+1i*repmat(y,size(x));
[e,N] = erf_(z,[],1); [e_,N_] = erf_(z,[],3); d = log10(abs(2*(e-e_)./(e+e_)));
figure, hold on, plot(y,N), plot([6 6],ylim,':k')
xlabel('y'), ylabel('N, algo 1'), legend('x = 0','x = 1')
figure, hold on, plot(y,N_), plot([6 6],ylim,':k')
xlabel('y'), ylabel('N, algo 3'), legend('x = 0','x = 1')
figure, hold on, plot(y,d), plot([6 6],ylim,':k')
xlabel('y'), ylabel('log10(algo 1-3 diff)'), legend('x = 0','x = 1')
title('erf\_(x+1i*y), algo 1 vs 3 (boundary at |x|<=1, |y|==6)')


text:
function [N]=gold_n(wl)

% the function calculates the index of refraction of gold (complex value), in the range of wavelength
% between 188nm and 1942nm.
% the input may be a scalar or vector, but it have to be an integer and in the right range.
% the values were taken from the paper:
%       P. B. Johnson, and R. W. Christy, "Optical Constants of the Noble Metals," Physical Review B 6, 4370 (1972).
% (the paper is available here: http://prb.aps.org/abstract/PRB/v6/i12/p4370_1 )
% the values for wavelength calculated with linear interpolation.
%
% Moshe Lindner, August 2010 (C).

if any((max(wl)>1942) || (min(wl)<188))
    error('Wavelength is out of range! (should be in the range of 188nm - 1942nm)')
    return
elseif any(wl~=ceil(wl))
    error('Wavelength must be an integer!')
    return
end

% n = real part
n=[0.92 0.56 0.43 0.35 0.27 0.22 0.17 0.16 0.14 0.13 0.14 0.21 0.29 0.43 0.62 1.04 1.31 1.39 1.45 1.46 1.47 1.46 1.48 1.5 1.48 1.48 1.54 1.53 1.53 1.49 1.47 1.43 1.38 1.35 1.33 1.33 1.32 1.32 1.3 1.31 1.3 1.3 1.3 1.3 1.33 1.33 1.34 1.32 1.28]';
n=n(end:-1:1);
% k = imaginary part
k=[13.78 11.21 9.519 8.145 7.15 6.35 5.663 5.083 4.542 4.103 3.697 3.272 2.863 2.455 2.081 1.833 1.849 1.914 1.948 1.958 1.952 1.933 1.895 1.866 1.871 1.883 1.898 1.893 1.889 1.878 1.869 1.847 1.803 1.749 1.688 1.631 1.577 1.536 1.497 1.46 1.427 1.387 1.35 1.304 1.277 1.251 1.226 1.203 1.188]';
k=k(end:-1:1);
% ev = photon energy in electronvolts
ev=[0.64 0.77 0.89 1.02 1.14 1.26 1.39 1.51 1.64 1.76 1.88 2.01 2.13 2.26 2.38 2.50 2.63 2.75 2.88 3 3.12 3.25 3.37 3.5 3.62 3.74 3.87 3.99 4.12 4.24 4.36 4.49 4.61 4.74 4.86 4.98 5.11 5.23 5.36 5.48 5.6 5.73 5.85 5.98 6.1 6.22 6.35 6.47 6.6]';
ev=ev(end:-1:1);
h=6.62606896e-34; % [J*sec] (Plank const.)
c=3e8; %[m/sec]  (speed of light)
q=1.6e-19; % [C] (electron charge)
lambda=1e9*h*c./(ev.*q);% wavelength in nanometer

lam=188:1942;
n1=interp1(lambda,n,lam,'linear');
k1=interp1(lambda,k,lam,'linear');

for q=1:length(wl)
    N(q)=n1(lam==wl(q)) + i*k1(lam==wl(q));
end


text:
function [CQ] = calculate_CQ(Seconds, Throttle, RPM, L, gramsMeas, varargin)
% This function calculates the parameter for a zero intercept linear fit
% between RPM^2 and torque of the form:
%   Torque = (RPM^2)*CQ
% Output is given in the form of returned CQ value and plots.
%
% Argument definitions:
%   L = length of torque rig arm from center of motor to scale contact
%   point (MUST BE IN INCHES)
%   gramsMeas = user measured values direct from scale, MUST BE IN GRAMS!
%   These values must correspond to ALL unique throttle values, even those
%   outside of "LowCutoff" and "HighCutoff."
% (Note: values of gramsMeas corresponding to throttle setting outside
% "LowCutoff" and "HighCutoff" are not used in calculation, therefore if
% measurements are only available for the range of "LowCutoff" to
% "HighCutoff", the vector can be padded with zeros in the appropriate
% locations.
%
% See also calculate_CR_B, calculate_CT.

% Assign variable length input arguments
if ((nargin-5)==0) % if no input arguments are given for cutoff range:
LowCutoff = min(Throttle); % Use full range of input throttle data
HighCutoff = max(Throttle);
else
LowCutoff = varargin{1}; % Ignore data below this throttle value
HighCutoff = varargin{2}; % Ignore data above this throttle value
end

% Find unique Throttle Setting values (Not counting 0)
ThrottleU = unique(Throttle(Throttle~=0));

% Check length of inputs match
lengthGM = length(gramsMeas);
lengthTU = length(ThrottleU);
if lengthGM ~= lengthTU
    msgbox({['Length of measured grams values vector (' num2str(lengthGM) ')'];
        [' and number of unique Throttle values (' num2str(lengthTU) ')']; ...
       ' do not match! Check input data!'},'Argument Error','warn');
    return
end

% Find unique/filtered Throttle Setting values between low/high cutoff limits
ThrottleUF = ThrottleU(ThrottleU>=LowCutoff & ThrottleU<=HighCutoff);

% Find coresponding grams measurments
gramsMeasU = gramsMeas(ThrottleU>=LowCutoff & ThrottleU<=HighCutoff); % Unique scale values

% Find an average RPM value for each unique/filtered Throttle Setting
RPMAverages100 = zeros(length(ThrottleUF),1); % preallocate
for i = 1:length(ThrottleUF) % For all of the unique and filtered throttle settings:
    RPMUF = RPM(Throttle==ThrottleUF(i));  % get RPM values cooresponding to unique/filtered Throttle value
    if length(RPMUF)>=120 % check to make sure there are enough data points for filtering
        RPMfilt100 = RPMUF(100:end); % RPMUF minus the first 100 data points (to remove transient RPM fluctuation as steady state value is reached)
    else RPMfilt100 = RPMUF((fix(length(RPMUF)/2)):end); % If there's a limited number of points, just take the second half of the values as steady state
    end
    RPMAverages100(i,1) = mean(RPMfilt100); % Take the average RPM value from these filtered data points
end
RPMsq = RPMAverages100.^2;

g = 9.81; % Gravity (m/s^2)
newtonMetersTorque = g*L*gramsMeasU/(1000*39.3701); % Torque in N*m (Conversion: (g-->kg and in-->m))
CQ = RPMsq\newtonMetersTorque; % This is a technique for forcing a zero intercept using matrix division techniques
CQp = [CQ;0];

% Plot the linear relationship and display calculated CQ value
figure
plot(RPMsq, newtonMetersTorque, '*', RPMsq, polyval(CQp, RPMsq))
xlabel('RPM^2 (rev^2/min^2)','FontSize',12)
ylabel('Torque (N*m)','FontSize',12)
title({['Linear Fit of Torque vs. RPM^2'];...
    [' ']; ['CQ: ' num2str(CQ) ' (N*m/RPM^2)']},'FontSize',14)
xlims = get(gca,'xlim');
ylims = get(gca,'ylim');
grid on

% Plot Residuals for assesing general quality of fit
predictedTorque = polyval(CQp, RPMsq); % predict thrust for each throttle setting of interest
residuals = newtonMetersTorque - predictedTorque; % find absolute error between prediction and actual
figure
plot(predictedTorque,residuals,'r*',[min(predictedTorque) max(predictedTorque)],[0 0])
title('CQ Residuals Plot','FontSize',14)
xlabel('Predicted Thrust (N)','FontSize',12)
ylabel('Prediction Error (N)','FontSize',12)
grid on
axis tight

repo: Quad-Sim

path: github-repos/MATLAB/dch33-Quad-Sim/Quad-Sim-961bc69d4939c8d0661eeb9820de668994262f65/Quadcopter Dynamic Modeling and Simulation/Data Acquisition & Analysis/MATLAB Data Analysis/calculate_CQ.m


text:
x = uiuc_sample;
x = x(1:256, 1:256);

Wop = wavelet_factory_2d(size(x));
%%
[Sx, Ux] = scat(x, Wop);

%%
plot_meta(Sx);

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
%% cell fill_method_listbox(void)
%
function [list] = get_method_list()
  L = load_method_list();
  [m,n] = size(L{1,1});
  list = {};
  k = 1;
  for i = 1 : m
    method_id = strtrim(char(L{1,1}(i)));
    method_name = strtrim(char(L{1,2}(i)));

    if(method_id == '-')
      method = [method_id];
    else
      method = [method_id ' - ' method_name];
    end

    list{k,1} = method;
    k = k + 1;
  end
end

author: andrewssobral

repo: lrslibrary


text:
function [valid, message] = applyUserValidation(SBMLStructure, level, version, packages, pkgVersion)
%  [valid, message] = applyUserValidation(SBMLStructure, level, version, packages, pkgVersion)
%
% a script that can be customized by the user to add additional
%   validation to the isSBML_Model function
%
% Takes
%
% 1. SBMLstructure - a MATLAB_SBML structure
% 2. level - an integer representing the SBML level
% 3. version - an integer representing the SBML version
% 4. packages - a cell array of the package prefixes used
%                              e.g. {'fbc'}
% 5. pkgVersion - an array of integers representing the version of each
%                               package; note indexing should match the
%                               packages indexing e.g. [1]
%
% Returns
%
% 1. valid =
%   - 1, if the structure represents passes the tests imposed
%   - 0, otherwise
% 2. a message explaining any failure

%<!---------------------------------------------------------------------------
% This file is part of libSBML.  Please visit http://sbml.org for more
% information about SBML, and the latest version of libSBML.
%
% Copyright (C) 2013-2018 jointly by the following organizations:
%     1. California Institute of Technology, Pasadena, CA, USA
%     2. EMBL European Bioinformatics Institute (EMBL-EBI), Hinxton, UK
%     3. University of Heidelberg, Heidelberg, Germany
%
% Copyright (C) 2009-2013 jointly by the following organizations:
%     1. California Institute of Technology, Pasadena, CA, USA
%     2. EMBL European Bioinformatics Institute (EMBL-EBI), Hinxton, UK
%
% Copyright (C) 2006-2008 by the California Institute of Technology,
%     Pasadena, CA, USA
%
% Copyright (C) 2002-2005 jointly by the following organizations:
%     1. California Institute of Technology, Pasadena, CA, USA
%     2. Japan Science and Technology Agency, Japan
%
% This library is free software; you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as
% published by the Free Software Foundation.  A copy of the license
% agreement is provided in the file named "LICENSE.txt" included with
% this software distribution and also available online as
% http://sbml.org/software/libsbml/license.html
%----------------------------------------------------------------------- -->

valid = 1;
message = '';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% example check that all GeneProducts listed have an id value

% if isempty(packages{1})
%     present = 0;
% else
%     [present, index] = ismember('fbc', packages);
% end;
%
% if (present)
%     if (pkgVersion(index) == 2)
%         gp = SBMLStructure.fbc_geneProduct;
%         [a, num] = size(gp);
%         for i=1:num
%             if (isempty(gp(num).fbc_id))
%                 valid = 0;
%                 message = 'geneProduct is missing the id attribute';
%             end;
%         end;
%     end;
% end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: opencobra

repo: cobratoolbox


text:
function [ fea, out ] = ex_axistressstrain5( varargin )
%EX_AXISTRESSSTRAIN5 Axisymmetric vibration modes of a hollow cylinder.
%
%   [ FEA, OUT ] = EX_AXISTRESSSTRAIN5( VARARGIN ) Axisymmetric
%   vibration modes of a hollow cylinder (NAFEMS Free Vibration Benchmark 41).
%
%   Reference:
%
%   [1] F. Abassian, D.J. Dawswell, and N.C. Knowles, Free Vibration
%   Benchmarks, Volume 3, NAFEMS, Glasgow, 1987.
%
%   Accepts the following property/value pairs.
%
%       Input       Value/{Default}        Description
%       -----------------------------------------------------------------------------------
%       igrid       scalar 0/{2}           Cell type (>0=quadrilaterals, <0=triangles)
%       sfun        string {sflag1}        Shape function
%       iplot       scalar 0/{1}           Plot solution (=1)
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       fea         struct                 Problem definition struct
%       out         struct                 Output struct

% Copyright 2013-2022 Precise Simulation, Ltd.

cOptDef = { 'igrid',    2;
            'sfun',     'sflag1';
            'iplot',    1;
            'tol',      0.01;
            'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;

E   = 2e11;
nu  = 0.3;
rho = 8000;

% Geometry definition.
fea.sdim = {'r' 'z'};
gobj = gobj_rectangle( 1.8, 2.2, 0, 10, 'R1' );
fea.geom.objects = { gobj };

fea.grid = rectgrid( abs(opt.igrid)*2, abs(opt.igrid)*50, [ 1.8, 2.2; 0, 10 ]);
if( opt.igrid<0 )
  fea.grid = quad2tri( fea.grid );
end

% Equations and problem definition.
fea = addphys( fea, @axistressstrain );
fea.phys.css.eqn.coef{1,end} = { nu  };
fea.phys.css.eqn.coef{2,end} = { E   };
fea.phys.css.eqn.coef{3,end} = { rho };
fea.phys.css.sfun            = { opt.sfun, opt.sfun };

% Solve problem.
fea = parsephys( fea );
fea = parseprob( fea );

[fea.sol.u,fea.sol.l] = solveeig( fea, 'fid', fid );

% Postprocessing.
if( opt.iplot>0 )
  postplot( fea, 'surfexpr', 'sqrt((r*u)^2+w^2)', 'solnum', 2 )
end

out = [];
f = sqrt(max(0,fea.sol.l))/(2*pi);
f_ref = [ 0; 243.773387; 378.534723; 394.046384; 397.467494; 405.041753 ];
out.err  = norm(f_ref-f)/norm(f_ref);
out.pass = out.err < opt.tol;

if( nargout==0 )
  clear fea out
end

author: precise-simulation

repo: featool-multiphysics


text:
function a = milnes ( m, n, x )

%*****************************************************************************80
%
%% MILNES returns the MILNES matrix.
%
%  Formula:
%
%    If ( I <= J )
%      A(I,J) = 1
%    else
%      A(I,J) = X(J)
%
%  Example:
%
%    M = 5, N = 5, X = ( 4, 7, 3, 8 )
%
%    1 1 1 1 1
%    4 1 1 1 1
%    4 7 1 1 1
%    4 7 3 1 1
%    4 7 3 8 1
%
%    M = 3, N = 6, X = ( 5, 7 )
%
%    1 1 1 1 1
%    5 1 1 1 1
%    5 7 1 1 1
%
%    M = 5, N = 3, X = ( 5, 7, 8 )
%
%    1 1 1
%    5 1 1
%    5 7 1
%    5 7 8
%    5 7 8
%
%  Properties:
%
%    A is generally not symmetric: A' /= A.
%
%    det ( A ) = ( 1 - X(1) ) * ( 1 - X(2) ) * ... * ( 1 - X(N-1) ).
%
%    A is singular if and only if X(I) = 1 for any I.
%
%    The family of matrices is nested as a function of N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Robert Gregory, David Karney,
%    Example 3.14, Example 5.24,
%    A Collection of Matrices for Testing Computational Algorithms,
%    Wiley, 1969, page 52, page 105,
%    LC: QA263.G68.
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of A.
%
%    Input, real X(*), the lower column values.
%    If M <= N, then X should be dimensioned M-1.
%    If N < M, X should be dimensioned N.
%
%    Output, real A(M,N), the matrix.
%
  for i = 1 : m
    for j = 1 : n
      if ( i <= j )
        a(i,j) = 1.0;
      else
        a(i,j) = x(j);
      end
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function X = traj_opt3c(path, total_time, ts)
% 3rd order trajectory optimization (minimum acceleration)
% by solving 4m linear equations. Assumes constant speed
% @author       Yiren Lu
% @email        luyiren [at] seas [dot] upenn [dot] edu
%
% @input:       path                (m+1) by 3 planning trajectory
%               total_time          total time
% @output       X                   4mx3 solution vector, where
%                                   X(4*(k-1)+1:4*k,d) is the coefficients
%                                   of x_k(t), i.e.,
%                                   x_k(t)=X(4*(k-1)+1:4*k,d)'*[t^3;t^2,t,1]
%                                   where d \in [1,2,3] to indicate dimension xyz
   path0 = path;
   % generate the trajectory here, and decide the total_time
   [m,n] = size(path0); % n == 3
   % there we set m = m - 1 for convenience
   m = m-1;
   % there are now in total m+1 points in the path, which seperate the path into
   % m subpaths.

   % ts(k) = t_{k+1}, e.g. ts(1) = t_0,
   % time planning according to the segment length

   % ts(m+1) = t_m
   % In 3rd order trajectory optimization, for each subpath there are 4
   % parameters
   % x_k(t) = c_{k,3}*t^3 + c_{k,2}*t^2 + c_{k,1}*t^1 + c_{k,0}
   %    for k = 1..(m-1)

   % X contains the parameters
   X = zeros(4*m,n);
   A = zeros(4*m, 4*m, n);
   Y = zeros(4*m,n);

   for i = 3
       A(:,:,i) = eye(4*m)*eps;
       % constraint 1: x_k(t_k) = x_{k+1}(t_k) = p_k, where p_k is a
       % waypoint
       % x_k(t) = c_{k,3}*t^3 + c_{k,2}*t^2 + c_{k,1}*t^1 + c_{k,0}
       %    for k = 1..(m-1)
       % e.g.   x_1(t_1) = x_2(t_1) = p_1;
       % there are in total 2*(m-1) constraints
       idx = 1; % constraint counter
       for k = 1:(m-1)
           A(idx, 4*(k-1)+1:4*k, i) = [ts(k+1)^3, ts(k+1)^2, ts(k+1), 1];
           Y(idx,i) = path0(k+1,i);
           idx = idx + 1;
           A(idx, 4*(k)+1:4*(k+1), i) = [ts(k+1)^3, ts(k+1)^2, ts(k+1), 1];
           Y(idx,i) = path0(k+1,i);
           idx = idx + 1;
       end
       % constraint 2: \dot{x}_k(t_k) = \dot{x}_{k+1}(t_k)
       % \dot{x}_k(t) = 3*c_{k,3}*t^2 + 2*c_{k,2}*t + c_{k,1}
       % e.g.  \dot{x}_1(t_1) = \dot{x}_2(t_1)
       % there are in total m-1 constraints
       for k = 1:(m-1)
           A(idx, 4*(k-1)+1:4*k, i) = [3*ts(k+1)^2, 2*ts(k+1), 1, 0];
%            A(idx, 4*(k)+1:4*(k+1), i) = -[3*ts(k+1)^2, 2*ts(k+1), 1, 0];
           Y(idx,i) = 0;
           idx = idx + 1;
           A(idx, 4*(k)+1:4*(k+1), i) = [3*ts(k+1)^2, 2*ts(k+1), 1, 0];
           Y(idx,i) = 0;
           idx = idx + 1;
       end

%        % constraint 3: \ddot{x}_k(t_k) = \ddot{x}_{k+1}(t_k)
%        % \ddot{x}_k(t) =6*c_{k,3}*t + 2*c_{k,2}
%        % e.g. \ddot{x}_1(t_1) = \ddot{x}_2(t_1)
%        % there are in total m-1 constraints
%        for k = 1:(m-1)
%            A(idx, 4*(k-1)+1:4*k, i) = [6*ts(k+1), 2, 0, 0];
%            A(idx, 4*(k)+1:4*(k+1), i) = -[6*ts(k+1), 2, 0, 0];
%            Y(idx,i) = 0;
%            idx = idx + 1;
%        end

       % so far there are 2*(m-1) + (m-1) + (m-1) = 4m - 4 constraints
       % there are 4 left:
       %    x_1(t_0) = p_0
       %    x_T(t_T) = p_T
       %    \dot{x}_0(t_0) = 0
       %    \dot{x}_T(t_T) = 0
       k = 1;
       A(idx, 4*(k-1)+1:4*k, i) = [ts(k)^3, ts(k)^2, ts(k), 1];
       Y(idx,i) = path0(k,i);
       idx = idx + 1;
       A(idx, 4*(k-1)+1:4*k, i) = [3*ts(k)^2, 2*ts(k), 1, 0];
       Y(idx,i) = 0;
       idx = idx + 1;
       k = m;
       A(idx, 4*(k-1)+1:4*k, i) = [ts(k+1)^3, ts(k+1)^2, ts(k+1), 1];
       Y(idx,i) = path0(k+1,i);
       idx = idx + 1;
       A(idx, 4*(k-1)+1:4*k, i) = [3*ts(k+1)^2, 2*ts(k+1), 1, 0];
       Y(idx,i) = 0;
       idx = idx + 1;
%        assert(rank(A(:,:,i))==4*m);
       X(:,i) = A(:,:,i)\Y(:,i);
   end
%    for i = 1:n
%        A(:,:,i) = eye(4*m)*eps;
%        % constraint 1: x_k(t_k) = x_{k+1}(t_k) = p_k, where p_k is a
%        % waypoint
%        % x_k(t) = c_{k,3}*t^3 + c_{k,2}*t^2 + c_{k,1}*t^1 + c_{k,0}
%        %    for k = 1..(m-1)
%        % e.g.   x_1(t_1) = x_2(t_1) = p_1;
%        % there are in total 2*(m-1) constraints
%        idx = 1; % constraint counter
%        for k = 1:(m-1)
%            A(idx, 4*(k-1)+1:4*k, i) = [ts(k+1)^3, ts(k+1)^2, ts(k+1), 1];
%            Y(idx,i) = path0(k+1,i);
%            idx = idx + 1;
%            A(idx, 4*(k)+1:4*(k+1), i) = [ts(k+1)^3, ts(k+1)^2, ts(k+1), 1];
%            Y(idx,i) = path0(k+1,i);
%            idx = idx + 1;
%        end
%        % constraint 2: \dot{x}_k(t_k) = \dot{x}_{k+1}(t_k)
%        % \dot{x}_k(t) = 3*c_{k,3}*t^2 + 2*c_{k,2}*t + c_{k,1}
%        % e.g.  \dot{x}_1(t_1) = \dot{x}_2(t_1)
%        % there are in total m-1 constraints
%        for k = 1:(m-1)
%            A(idx, 4*(k-1)+1:4*k, i) = [3*ts(k+1)^2, 2*ts(k+1), 1, 0];
%            A(idx, 4*(k)+1:4*(k+1), i) = -[3*ts(k+1)^2, 2*ts(k+1), 1, 0];
%            Y(idx,i) = 0;
%            idx = idx + 1;
%        end
%
%        % constraint 3: \ddot{x}_k(t_k) = \ddot{x}_{k+1}(t_k)
%        % \ddot{x}_k(t) =6*c_{k,3}*t + 2*c_{k,2}
%        % e.g. \ddot{x}_1(t_1) = \ddot{x}_2(t_1)
%        % there are in total m-1 constraints
%        for k = 1:(m-1)
%            A(idx, 4*(k-1)+1:4*k, i) = [6*ts(k+1), 2, 0, 0];
%            A(idx, 4*(k)+1:4*(k+1), i) = -[6*ts(k+1), 2, 0, 0];
%            Y(idx,i) = 0;
%            idx = idx + 1;
%        end
%
%        % so far there are 2*(m-1) + (m-1) + (m-1) = 4m - 4 constraints
%        % there are 4 left:
%        %    x_1(t_0) = p_0
%        %    x_T(t_T) = p_T
%        %    \dot{x}_0(t_0) = 0
%        %    \dot{x}_T(t_T) = 0
%        k = 1;
%        A(idx, 4*(k-1)+1:4*k, i) = [ts(k)^3, ts(k)^2, ts(k), 1];
%        Y(idx,i) = path0(k,i);
%        idx = idx + 1;
%        A(idx, 4*(k-1)+1:4*k, i) = [3*ts(k)^2, 2*ts(k), 1, 0];
%        Y(idx,i) = 0;
%        idx = idx + 1;
%        k = m;
%        A(idx, 4*(k-1)+1:4*k, i) = [ts(k+1)^3, ts(k+1)^2, ts(k+1), 1];
%        Y(idx,i) = path0(k+1,i);
%        idx = idx + 1;
%        A(idx, 4*(k-1)+1:4*k, i) = [3*ts(k+1)^2, 2*ts(k+1), 1, 0];
%        Y(idx,i) = 0;
%        idx = idx + 1;
% %        assert(rank(A(:,:,i))==4*m);
%        X(:,i) = A(:,:,i)\Y(:,i);
%    end

end

author: yrlu

repo: quadrotor

save_path: github-repos/MATLAB/yrlu-quadrotor


text:
% 改进的欧拉法求解初值问题
clear;
f = @(x,y) (y - 2 * x / y);
h = 0.1;   % 步长
x = 0:h:2;  % 离散点
y(1) = 1;   % 初值
y_xn = (1 + 2 * x).^(1/2);   % 解析解

% 改进的欧拉法
for i = 2 : length(x)
    yp = y(i-1) + h * f(x(i-1),y(i-1));
    yc = y(i-1) + h * f(x(i),yp);
    y(i) = (yp + yc) / 2;
end
y
y_xn
plot(x,y,'o-',x,y_xn,'k-');
legend('改进的欧拉法','精确解');

repo: NumericalAnalysis


text:
function S = getStartFinish_StateVisits (vp,T,options)
% Decomposes a Viterbi path into a list of visits, each characterised by
% start time point and finish time point.
%
% S is a (N x K) cell, where N is the number of data time series
% (length of T) and K is the number of states; each element of S contains
% a matrix, which n-th row corresponds to the n-th visit to the
% corresponding state, with two columns (start and finish time points).
% The values of these matrices are in time points (not seconds), and are
% relative to the length of the state time courses, which will be a bit
% shorter if a MAR or a time-delay embedded observation model were used.
%
% T must be a vector with the length of each time series, and vp is the
% Viterbi path as returned by hmmmar or hmmdecode.
% The parameter 'options' must be the same than the one supplied to the
% hmmmar function for training.
%
% Author: Diego Vidaurre, OHBA, University of Oxford (2018)

if nargin<3, options = struct(); options.Fs = 1; options.downsample = 1; end
if ~isfield(options,'Fs'), options.Fs = 1; end
if ~isfield(options,'downsample'), options.downsample = options.Fs; end

if nargin < 2, T = length(vp); end

if iscell(T)
    if size(T,1) == 1, T = T'; end
    for i = 1:length(T)
        if size(T{i},1)==1, T{i} = T{i}'; end
    end
    T = cell2mat(T);
end
N = length(T);

r = 1;
if isfield(options,'downsample') && options.downsample>0
    r = (options.downsample/options.Fs);
end

if isfield(options,'order') && options.order > 0
    T = ceil(r * T);
    T = T - options.order;
elseif isfield(options,'embeddedlags') && length(options.embeddedlags) > 1
    d1 = -min(0,options.embeddedlags(1));
    d2 = max(0,options.embeddedlags(end));
    T = T - (d1+d2);
    T = ceil(r * T);
end

K = length(unique(vp));
S = cell(N,K);

for j = 1:N
    t0 = sum(T(1:j-1));
    ind = (1:T(j)) + t0;
    vp_j = vp(ind);
    for k = 1:K
       vp_jk = zeros(T(j),1);
       vp_jk(vp_j==k) = 1;
       vp_jk_diff = diff(vp_jk);
       Nvisits = sum(vp_jk_diff==1) + (vp_jk(1)==1);
       S{j,k} = zeros(Nvisits,2);
       if vp_jk(1)==1 % first visit is to this state
           S{j,k}(1,1) = 1;
           S{j,k}(2:end,1) = find(vp_jk_diff==1)+1;
       else
           S{j,k}(:,1) = find(vp_jk_diff==1)+1;
       end
       if vp_jk(end)==1 % last visit is to this state
           S{j,k}(1:end-1,2) = find(vp_jk_diff==(-1));
           S{j,k}(end,2) = T(j);
       else
           S{j,k}(:,2) = find(vp_jk_diff==(-1));
       end
    end
end

end

author: OHBA-analysis

repo: HMM-MAR


text:
function res = spp_test(spp_model, imdb, roidb, feat_cache, suffix, fast, evaluate)
% res = spp_test(spp_model, imdb, roidb, feat_cache, suffix, fast, evaluate)
%   Compute test results using the trained spp_model on the
%   image database specified by imdb. Results are saved
%   with an optional suffix.
%
% Adapted from spp code written by Ross Girshick
% AUTORIGHTS
% ---------------------------------------------------------
% Copyright (c) 2014, Shaoqing Ren
%
% This file is part of the SPP code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------
% Copyright (c) 2014, Ross Girshick
%
% This file is part of the R-CNN code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------

if nargin < 5
    suffix = '';
end

if nargin < 6
    fast = false;
end

if nargin < 7
    evaluate = true;
end

t_start = tic();

conf = spp_config('sub_dir', fullfile(spp_model.cache_name, imdb.name));
image_ids = imdb.image_ids;

% assume they are all the same
feat_opts = spp_model.training_opts;
feat_opts.feat_cache = feat_cache;
num_classes = length(spp_model.classes);
rois = roidb.rois;

if conf.use_gpu
    spp_model.cnn.layers = spp_layers_in_gpu(spp_model.cnn.layers);
end

if ~exist('suffix', 'var') || isempty(suffix)
  suffix = '';
else
  suffix = ['_' suffix];
end

try
  aboxes = cell(num_classes, 1);
  for i = 1:num_classes
    load([conf.cache_dir spp_model.classes{i} '_boxes_' imdb.name suffix]);
    aboxes{i} = boxes;
  end
catch
  aboxes = cell(num_classes, 1);
  box_inds = cell(num_classes, 1);
  for i = 1:num_classes
    aboxes{i} = cell(length(image_ids), 1);
    box_inds{i} = cell(length(image_ids), 1);
  end

  max_per_set = 5 * length(image_ids);
  max_per_image = 100;
  top_scores = cell(num_classes, 1);
  thresh = -1.5*ones(num_classes, 1);

  if ~isfield(spp_model, 'folds')
    folds{1} = 1:length(image_ids);
  else
    folds = spp_model.folds;
  end

  count = 0;
  for f = 1:length(folds)
    for i = folds{f}
      count = count + 1;
      fprintf('%s: test (%s) %d/%d ', procid(), imdb.name, count, length(image_ids));
      th = tic;
      d = rois(i);
      d.feat = spp_load_cached_poolX_features(spp_model.spp_pooler, feat_opts.feat_cache, ...
          imdb.name, image_ids{i}, d.boxes);
      if isempty(d.feat)
        continue;
      end
      d.feat = spp_poolX_to_fcX(d.feat, feat_opts.layer, spp_model, conf.use_gpu);
      d.feat = spp_scale_features(d.feat, feat_opts.feat_norm_mean);
      zs = bsxfun(@plus, spp_model.detectors(f).W * d.feat, spp_model.detectors(f).B)';

      for j = 1:num_classes
        boxes = d.boxes;
        scores = zs(:,j);
        I = find(~d.gt & scores > thresh(j));
        keep = nms(cat(2, single(boxes(I,:)), single(scores(I))), 0.3);
        I = I(keep);
        if ~isempty(I)
          [~, ord] = sort(scores(I), 'descend');
          ord = ord(1:min(length(ord), max_per_image));
          I = I(ord);
          boxes = boxes(I,:);
          scores = scores(I);
          aboxes{j}{i} = cat(2, single(boxes), single(scores));
          box_inds{j}{i} = I;
        else
          aboxes{j}{i} = zeros(0, 5, 'single');
          box_inds{j}{i} = [];
        end
      end

      if mod(count, 1000) == 0
          [aboxes{j}, box_inds{j}, thresh(j)] = ...
             keep_top_k(aboxes{j}, box_inds{j}, i, max_per_set, thresh(j));
      end
      fprintf(' time: %.3fs\n', toc(th));
      if mod(count, 1000) == 0
        disp(thresh);
      end

    end
  end

  for i = 1:num_classes

    top_scores{i} = sort(top_scores{i}, 'descend');
    if (length(top_scores{i}) > max_per_set)
        thresh(i) = top_scores{i}(max_per_set);
    end

    % go back through and prune out detections below the found threshold
    for j = 1:length(image_ids)
      if ~isempty(aboxes{i}{j})
        I = find(aboxes{i}{j}(:,end) < thresh(i));
        aboxes{i}{j}(I,:) = [];
        box_inds{i}{j}(I,:) = [];
      end
    end

    save_file = [conf.cache_dir spp_model.classes{i} '_boxes_' imdb.name suffix];
    boxes = aboxes{i};
    inds = box_inds{i};
    save(save_file, 'boxes', 'inds');
    clear boxes inds;
  end
end

fprintf('spp_test_spm in %f seconds.\n', toc(t_start));

% ------------------------------------------------------------------------
% Peform AP evaluation
% ------------------------------------------------------------------------
if ~evaluate
    res = [];
    return;
end

if isequal(imdb.eval_func, @imdb_eval_voc)
    if fast
        classes = spp_model.classes;
        cache_name = spp_model.cache_name;
        parfor model_ind = 1:num_classes
          cls = classes{model_ind};
          res(model_ind) = imdb.eval_func(cls, aboxes{model_ind}, imdb, cache_name, suffix, fast);
        end
    else
        for model_ind = 1:num_classes
          cls = spp_model.classes{model_ind};
          res(model_ind) = imdb.eval_func(cls, aboxes{model_ind}, imdb, spp_model.cache_name, suffix, fast);
        end
    end
else
% ilsvrc
    res = imdb.eval_func(aboxes, imdb, spp_model.cache_name, suffix, fast);
end

if ~isempty(res)
    fprintf('\n~~~~~~~~~~~~~~~~~~~~\n');
    fprintf('Results:\n');
    aps = [res(:).ap]' * 100;
    disp(aps);
    disp(mean(aps));
    fprintf('~~~~~~~~~~~~~~~~~~~~\n');
end

% ------------------------------------------------------------------------
function [boxes, box_inds, thresh] = keep_top_k(boxes, box_inds, end_at, top_k, thresh)
% ------------------------------------------------------------------------
% Keep top K
X = cat(1, boxes{1:end_at});
if isempty(X)
  return;
end
scores = sort(X(:,end), 'descend');
thresh = scores(min(length(scores), top_k));
for image_index = 1:end_at
  bbox = boxes{image_index};
  keep = find(bbox(:,end) >= thresh);
  boxes{image_index} = bbox(keep,:);
  box_inds{image_index} = box_inds{image_index}(keep);
end

author: ShaoqingRen


text:
function obj = loadDynamics(obj, file_path, mmat_names, mmat_ddx_names, vf_names, varargin)
    % load the symbolic expression of system dynamical equations from MX
    % binary files for fast loading
    %
    % Parameters:
    %  file_path: the path to the files @type char
    %  vf_names: the names of the drift vectors Fvec @type cellstr
    %  skip_load_vf: indicates whether to skip the loading of drift vectors
    %  which could takes a long time @type logical

    opts = struct(varargin{:});
    if isfield(opts, 'DelayCoriolisSet')
        delay_set = opts.DelayCoriolisSet;
    else
        delay_set = false;
    end

    x = obj.States.x;
    % load the mass matrix using the default name

    if ~isempty(mmat_names)

        for i=1:numel(mmat_names)
            mmat_name = mmat_names{i};
            mmat_ddx_name = mmat_ddx_names{i};
            M = SymFunction(mmat_name,[],{x});
            M = load(M,file_path);
            [nr,nc] = size(M);
            assert(nr==obj.numState && nc==obj.numState,...
                'The size of the mass matrix should be (%d x %d).',obj.numState,obj.numState);

            obj.Mmat{i} = M;

            if strcmp(obj.Type,'SecondOrder')
                ddx = obj.States.ddx;
                MmatDx = SymFunction(mmat_ddx_name,[],{x,ddx});
            else
                dx = obj.States.dx;
                MmatDx = SymFunction(mmat_ddx_name,[],{x,dx});
            end
            MmatDx = load(MmatDx,file_path);
            obj.MmatDx{i} = MmatDx;
        end
    else
        mmat_ddx_name = ['MmatDx_' obj.Name];
        obj.Mmat = [];
        if strcmp(obj.Type,'SecondOrder')
            ddx = obj.States.ddx;
            obj.MmatDx = SymFunction(mmat_ddx_name,-ddx,{x,ddx});
        else
            dx = obj.States.dx;
            obj.MmatDx = SymFunction(mmat_ddx_name,-dx,{x,dx});
        end
    end
    obj.MmatName_ = cellfun(@(f)f.Name, obj.Mmat,'UniformOutput',false);

    % load the drift vector
    if nargin > 4
        if ~iscell(vf_names), vf_names = {vf_names}; end
        vf = cell(size(vf_names));
        for i=1:numel(vf_names)

            if strcmp(obj.Type,'SecondOrder') % second order system
                sfun_vf = SymFunction(vf_names{i},[],{obj.States.x,obj.States.dx});
            else                         % first order system
                sfun_vf = SymFunction(vf_names{i},[],{obj.States.x});
            end
            if ~isempty(vf_names{i}) && ~delay_set
                sfun_vf = load(sfun_vf, file_path);
            end
            vf{i} = sfun_vf;

        end
        obj.Fvec = vf;
        obj.FvecName_ = cellfun(@(f)f.Name, vf,'UniformOutput',false);
    end

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


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

%Empty workspace and close figures
close all;
clear;

%Number of BS antennas
M = 100;

%Angular standard deviation in the local scattering model (in degrees)
ASDs = [10 30];

%Set the nominal angle of the desired UE
varphiDesired = pi/6;

%Set range of nominal angles of the interfering UE
varphiInterfererDegrees = -180:1:180;
varphiInterfererRadians = varphiInterfererDegrees*(pi/180);

%Define the antenna spacing (in number of wavelengths)
antennaSpacing = 1/2; %Half wavelength distance

%Preallocate matrix for storing the simulation results
variance = zeros(length(varphiInterfererRadians),length(ASDs));

%% Go through the range of ASDs
for n = 1:length(ASDs)

    %Output simulation progress
    disp([num2str(n) ' ASDs out of ' num2str(length(ASDs))]);

    %Compute spatial correlation matrix of the desired UE
    R1 = functionRlocalscattering(M,varphiDesired,ASDs(n),antennaSpacing);

    %Go through all angles of the interfering UE
    for r = 1:length(varphiInterfererRadians)

        %Compute spatial correlation matrix of the interfering UE
        R2 = functionRlocalscattering(M,varphiInterfererRadians(r),ASDs(n),antennaSpacing);

        %Compute variance of favorable propagation according to (2.19)
        variance(r,n) = real(trace(R1*R2)/(trace(R1)*trace(R2)));

    end

end

%% Plot the simulation results
figure;
hold on; box on;

plot(varphiInterfererDegrees,variance(:,1),'r--','LineWidth',1);
plot(varphiInterfererDegrees,variance(:,2),'b-.','LineWidth',1);
plot(varphiInterfererDegrees,1/M*ones(length(varphiInterfererDegrees),1),'k-','LineWidth',1);

xlabel('Angle of interfering UE [degree]');
ylabel('Variance in (2.19)');
xlim([-180 180]);

legend('Gaussian, ASD 10^o','Gaussian, ASD 30^o','Uncorrelated','Location','NorthWest');

author: emilbjornson

repo: massivemimobook


text:
classdef grassmannian < handle
    % Ways to generate Grassmannian frames.
    properties
    end

    methods(Static)

        function result = minimum_coherence(m, n)
            % Returns the minimum coherence of an m times n matrix
            num = n - m;
            den = m * (n - 1);
            result = sqrt(num / den);
        end

        function n = n_upper_bound(m)
            % for a given value of m, returns the upper bound on n.
            a = m * (m  + 1) / 2;
            n = m  * (m + 1) / 2 - 1;
            lower_limit = m;
            while true
                b = (n - m) * (n - m + 1) / 2;
                upper_limit = min (a, b);
                if n < upper_limit
                    break;
                end
                n = n- 1;
            end
        end
        function [ns, coherences] = min_coherence_max_n(ms)
            if ~isvector(ms)
                error('ms must be a vector');
            end
            mm = length(ms);
            ns = zeros(1, mm);
            coherences = zeros(1, mm);
            for i=1:mm
                m = ms(i);
                n = spx.dict.grassmannian.n_upper_bound(m);
                ns(i) = n;
                coherences(i) = spx.dict.grassmannian.minimum_coherence(m, n);
            end
        end

        function n = max_n_for_coherence(m, mu)
            % estimates a value of n for m x n matrix which achieves a
            % given coherence
            num = (mu^2 -1) * m;
            den = mu^2 * m - 1;
            n = floor(num / den);
            if n < 0
                error('A Grassmannian frame for this dimension will never have so high coherence.');
            end
        end

        function result = alternate_projections(dict, options)
            % Converts a given dictionary to Grassmannian frame via
            % alternate projections
            [n, d] = size(dict);
            % We assume that the dictionary contains unit norm atoms
            target_mu  = spx.dict.grassmannian.minimum_coherence(n, d);
            % Compute the Gram matrix
            gram_matrix = dict' * dict;
            iterations = 1000;
            entry_shrinkage_factor = 0.9;
            coherence_shrinkage_factor = 0.9;
            verbose = false;
            if nargin > 1 && isstruct(options)
                if isfield(options, 'iterations')
                    iterations = options.iterations;
                end
                if isfield(options, 'entry_shrinkage_factor')
                    entry_shrinkage_factor = options.entry_shrinkage_factor;
                end
                if isfield(options, 'coherence_shrinkage_factor')
                    coherence_shrinkage_factor = options.coherence_shrinkage_factor;
                end
                if isfield(options, 'verbose')
                    verbose = options.verbose;
                end
            end
            total_entries = d^2 - d; % excluding the main diagonal.
            % Identify the indices in the gram matrix which are off diagonal
            off_diagonal_indices = abs(gram_matrix(:) - 1)> 1e-6;
            coherence_array = zeros(iterations, 1);
            mean_coherence_array = zeros(iterations, 1);
            for k=1:iterations
                % Convert gram matrix to an array
                gram_array = gram_matrix(:);
                % Absolute Gram matrix
                abs_gram_array = abs(gram_array);
                % sort the inner products by their absolute values
                sorted_gram_array = sort(abs_gram_array);
                upper_threhold_index = round(entry_shrinkage_factor * total_entries);
                upper_threshold = sorted_gram_array(upper_threhold_index);
                % identify coherence values above the threshold
                above_indices = off_diagonal_indices & (abs_gram_array > upper_threshold);
                % map the boolean flags to positions.
                above_indices = find(above_indices);
                above_values = abs_gram_array(above_indices);
                current_mu = max(above_values);
                mean_mu = mean(above_values);
                coherence_array(k) = current_mu;
                mean_coherence_array(k) = mean_mu;
                if verbose
                    fprintf('%6d:, Coherence:: target: %12.8f, ', k, target_mu);
                    fprintf('current: %12.8f, threshold: %12.8f, above: %d, mean: %12.8f\n', ...
                        current_mu, upper_threshold, length(above_values), mean_mu);
                else
                    fprintf('.');
                    if mod(k, 100) == 0
                        fprintf('\n');
                    end
                end
                % Update off diagonal entries
                gram_matrix(above_indices)=gram_matrix(above_indices)*coherence_shrinkage_factor;

                % reduce the rank back to n
                [U,S,V]=svd(gram_matrix);
                % Ensure that all higher singular values are set to 0.
                S(n+1:end,n+1:end)=0;
                % Reconstruct the Gram matrix.
                gram_matrix=U*S*V';
                % Ensure that the diagonal elements of G continue to be 1.
                % Compute  d = diag(G)  d^{-1/2}.
                gram_diag_sqrt_inv = diag(1./sqrt(diag(gram_matrix)));
                % Compute G = d^{-1/2} G d^{-1/2}.
                gram_matrix = gram_diag_sqrt_inv*gram_matrix*gram_diag_sqrt_inv;
            end
            % final dictionary
            % perform SVD of the gram matrix.
            [U,S,V]=svd(gram_matrix);
            % Construct the dictionary
            dict=S(1:n,1:n).^0.5*U(:,1:n)';
            result.dictionary = dict;
            result.coherence_array = coherence_array;
            result.mean_coherence_array = mean_coherence_array;
            result.target_coherence = target_mu;
            result.iterations = iterations;
        end
    end
end

author: indigits

repo: sparse-plex


text:
function plot_MapsSeries4(sFilename)
%
%
%
load(sFilename)

%% plot results mean z values
cmin=-4;cmax=4;

for i=1:size(params.mPolZ,2)

figure_w_normalized_uicontrolunits('Name','Z-probability by overlap','Position',[100 25 400 400]);
pcolor(params.vX,params.vY,...
    reshape(calc_ProbColorbar2Value(1-params.mPolZ(:,i)),...
    size(params.vY,1),size(params.vX,1)));
xlabel('longitude');
ylabel('latitude');
title(sprintf('%6.1f-%6.1f vs. %6.1f-%6.1f',params.fTstart,...
    params.fT-params.mVar(i,1),params.fT-params.mVar(i,1),params.fT));
plot_ProbColorbar2(cmin, cmax);
% set(gca,'XAxisLocation','Top','XTick',0.5,'XTickLabel', 'P(z)')
shading interp;
hold on;plot(params.mCatalog(:,1),params.mCatalog(:,2),'k.','MarkerSize',0.1);
plot_WiemerWyss1994
xlim([-117.1 -115.6]);
ylim([33 35.2]);

end

end

function plot_WiemerWyss1994
mLatLon=[[-117.1 33];
    [-115.6 33];
    [-115.6 35.2];
    [-117.1 35.2];
    [-117.1 33]];
hold on;plot(mLatLon(:,1),mLatLon(:,2),'k--');
mLanders=[-116.4 34.3];
hold on;plot(mLanders(1),mLanders(2),'ko','MarkerSize',20)
end

author: CelsoReyes


text:
function asp = aspectRatio(grains,varargin)
% aspectratio = length / width
%
% the aspect ratio is the ratio between the two
% <grain2d.principalComponents.html,principal componentes> of a grain
%
% Input
%  g - @grain2d
%
% Output
%  asp   - aspect--ratio
%
% See also
% grain2d/principalcomponents

[~,a,b] = principalComponents(grains);
asp = abs(a ./ b);

% asp2 = diameter(grains).^2 ./ area(grains);

author: mtex-toolbox

repo: mtex

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


text:
function clenshaw_curtis_rule ( n, a, b, filename )

%*****************************************************************************80
%
%% CLENSHAW_CURTIS_RULE generates a Clenshaw Curtis rule.
%
%  Discussion:
%
%    This program computes a standard Clenshaw Curtis quadrature rule
%    and writes it to a file.
%
%    The user specifies:
%    * N, the number of points in the rule;
%    * A, the left endpoint;
%    * B, the right endpoint;
%    * FILENAME, the root name of the output files.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 February 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CLENSHAW_CURTIS_RULE\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Compute a Clenshaw Curtis rule for approximating\n' );
  fprintf ( 1, '    Integral ( A <= x <= B ) f(x) dx\n' );
  fprintf ( 1, '  of order N.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The user specifies N, A, B, and FILENAME.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  N is the number of points:\n' );
  fprintf ( 1, '  A is the left endpoint;\n' );
  fprintf ( 1, '  B is the right endpoint;\n' );
  fprintf ( 1, '  FILENAME is used to generate 3 files:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    filename_w.txt - the weight file\n' );
  fprintf ( 1, '    filename_x.txt - the abscissa file.\n' );
  fprintf ( 1, '    filename_r.txt - the region file.\n' );
%
%  Get N.
%
  if ( nargin < 1 )
    n = input ( '  Enter the rule order N:  ' );
  elseif ( ischar ( n ) )
    n = str2num ( n );
  end
%
%  Get A.
%
  if ( nargin < 2 )
    a = input ( '  Enter the left endpoint A:  ' );
  elseif ( ischar ( a ) )
    a = str2num ( a );
  end
%
%  Get B.
%
  if ( nargin < 3 )
    b = input ( '  Enter the right endpoint B:  ' );
  elseif ( ischar ( b ) )
    b = str2num ( b );
  end
%
%  Get FILENAME.
%
  if ( nargin < 4 )
    fprintf ( 1,  '\n' );
    fprintf ( 1,  '  FILENAME is the ''root name'' of the quadrature files).\n' );
    filename = input ( '  Enter FILENAME as a quoted string:  ' );
  end
%
%  Input summary.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  N = %d\n', n );
  fprintf ( 1, '  A = %f\n', a );
  fprintf ( 1, '  B = %f\n', b );
  fprintf ( 1, '  FILENAME = "%s".\n', filename );
%
%  Construct the rule and output it.
%
  r = [ a; b ];
  [ x, w ] = clenshaw_curtis_compute ( n );
%
%  Rescale the rule.
%
  [ x, w ] = rescale ( a, b, x, w );
%
%  Write the rule.
%
  rule_write ( n, filename, x, w, r );
%
%  Terminate.
%
  fprintf ( 1,  '\n' );
  fprintf ( 1,  'CLENSHAW_CURTIS_RULE:\n' );
  fprintf ( 1,  '  Normal end of execution.\n' );
  fprintf ( 1,  '\n' );
  timestamp ( );

  return
end
function [ x, w ] = clenshaw_curtis_compute ( n )

%*****************************************************************************80
%
%% CLENSHAW_CURTIS_COMPUTE computes a Clenshaw Curtis quadrature rule.
%
%  Discussion:
%
%    Our convention is that the abscissas are numbered from left to right.
%
%    The rule is defined on [-1,1].
%
%    The integral to approximate:
%
%      Integral ( -1 <= X <= 1 ) F(X) dX
%
%    The quadrature rule:
%
%      Sum ( 1 <= I <= N ) W(I) * F ( X(I) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the rule.
%    1 <= N.
%
%    Output, real X(N), the abscissas.
%
%    Output, real W(N), the weights.
%
  if ( n < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CLENSHAW_CURTIS_COMPUTE - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of N = %d\n', n );
    error ( 'CLENSHAW_CURTIS_COMPUTE - Fatal error!' );
  end

  w = zeros ( n, 1 );
  x = zeros ( n, 1 );

  if ( n == 1 )
    x(1) = 0.0;
    w(1) = 2.0;
    return
  end

  for i = 1 : n
    x(i) = cos ( ( n - i ) * pi / ( n - 1 ) );
  end

  x(1) = -1.0;
  if ( mod ( n, 2 ) == 1 )
    x((n+1)/2) = 0.0;
  end
  x(n) = +1.0;

  w(1:n) = 1.0;

  for i = 1 : n

    theta = ( i - 1 ) * pi / ( n - 1 );

    for j = 1 : ( n - 1 ) / 2

      if ( 2 * j == ( n - 1 ) )
        b = 1.0;
      else
        b = 2.0;
      end

      w(i) = w(i) - b * cos ( 2.0 * j * theta ) / ( 4 * j * j - 1 );

    end

  end

  w(1)     =       w(1)     / ( n - 1 );
  w(2:n-1) = 2.0 * w(2:n-1) / ( n - 1 );
  w(n)     =       w(n)     / ( n - 1 );

  return
end
function r8mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% R8MAT_WRITE writes an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real TABLE(M,N), the points.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'R8MAT_WRITE - Error!' );
  end
%
%  Write the data.
%
%  For smaller data files, and less precision, try:
%
%     fprintf ( output_unit, '  %14.6f', table(i,j) );
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %24.16f', table(i,j) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function [ x, w ] = rescale ( a, b, x, w )

%*****************************************************************************80
%
%% RESCALE rescales a Legendre quadrature rule from [-1,+1] to [A,B].
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 October 2009
%
%  Author:
%
%    John Burkardt.
%
%  Parameters:
%
%    Input, real A, B, the endpoints of the new interval.
%
%    Input, integer N, the order.
%
%    Input, real X(N), the abscissas for [-1,+1].
%
%    Input, real W(N), the weights for [-1,+1].
%
%    Output, real X(N), the abscissas for [A,B].
%
%    Output, real W(N), the weights for [A,B].
%
  x = 0.5 * ( ( x + 1.0 ) * b - ( x - 1.0 ) * a );

  w = 0.5 * ( b - a ) * w;

  return
end
function rule_write ( order, filename, x, w, r )

%*****************************************************************************80
%
%% RULE_WRITE writes a quadrature rule to a file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ORDER, the order of the rule.
%
%    Input, string FILENAME, specifies the output files.
%    write files 'filename_w.txt', 'filename_x.txt', 'filename_r.txt' defining
%    weights, abscissas, and region.
%
%    Input, real X(ORDER), the abscissas.
%
%    Input, real W(ORDER), the weights.
%
%    Input, real R(2), the region.
%
  filename_x = strcat ( filename, '_x.txt' );
  filename_w = strcat ( filename, '_w.txt' );
  filename_r = strcat ( filename, '_r.txt' );

  fprintf ( 1, '\n' );
  fprintf ( 1,'  Creating quadrature files.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  "Root" file name is   "%s".\n', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Weight file will be   "%s".\n', filename_w );
  fprintf ( 1, '  Abscissa file will be "%s".\n', filename_x );
  fprintf ( 1, '  Region file will be   "%s".\n', filename_r );

  r8mat_write ( filename_w, 1, order, w' );
  r8mat_write ( filename_x, 1, order, x' );
  r8mat_write ( filename_r, 1, 2,     r' );

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end


text:
function q = divine_overall_quality(r)

fullpath = mfilename('fullpath');
idx      = find(fullpath == filesep);
dfolder = fullpath(1:(idx(end)-1));

% Function to compute overall quality given feature vector 'r'

load(fullfile(dfolder,'data_live_trained.mat'))

%% Classification
atrain = repmat(a_class,[size(r,1) 1]);btrain = repmat(b_class,[size(r,1) 1]);
x_curr = atrain.*r+btrain;

% [pred_class acc p] = svmpredict(1,x_curr,model_class,'-b 1');
% use `git clone git://github.com/gregfreeman/libsvm.git -b new_matlab_interface`

[pred_class p] = svmpredict(x_curr,model_class,struct('output','probability'));

%% Regression
for i = 1:5
    atrain = repmat(a_reg(i,:),[size(r,1) 1]);btrain = repmat(b_reg(i,:),[size(r,1) 1]);
    x_curr = atrain.*r+btrain;
%    [q(i) reg_acc(i,:)] = svmpredict(1,x_curr,model_reg{i});
    [q(i)] = svmpredict(x_curr,model_reg{i});
end
%% Final Score
% clc
q = sum(p.*q);

author: dsoellinger


text:
function [x, infos] = robust_online_mu_nmf(V, rank, in_options)
% Robust online non-negative matrix factorization (ONMF) with outliers (RONMF) algorithm.
%
% Inputs:
%       matrix      V
%       rank        rank
%       in_options  options
% Output:
%       w           solution of w
%       infos       information
%
% References:
%       R. Zhao and Y. F. Tan,
%       "Online nonnegative matrix factorization with outliers,"
%       ICASSP, 2016.
%
%
% This file is part of NMFLibrary.
%
% Created by H.Sakai and H.Kasai on Feb. 12, 2017
%
% Change log:
%
%       May. 20, 2019 (Hiroyuki Kasai): Added initialization module.
%
%       Jul. 12, 2022 (Hiroyuki Kasai): Modified code structures.
%

    % set dimensions and samples
    [m, n] = size(V);

    % set local options
    local_options = [];
    local_options.lambda        = 1;
    local_options.x_init_robust = true;

    % check input options
    if ~exist('in_options', 'var') || isempty(in_options)
        in_options = struct();
    end
    % merge options
    options = mergeOptions(get_nmf_default_options(), local_options);
    options = mergeOptions(options, in_options);

    % initialize factors
    init_options = options;
    [init_factors, ~] = generate_init_factors(V, rank, init_options);
    Wt = init_factors.W;
    H = init_factors.H;
    R = init_factors.R;

    % initialize
    method_name = 'Robust-Online-MU';
    epoch = 0;
    grad_calc_count = 0;
    l = zeros(m, options.batch_size) + options.lambda;

    if options.verbose > 0
        fprintf('# %s: started ...\n', method_name);
    end

    % store initial info
    clear infos;
    [infos, f_val, optgap] = store_nmf_info(V, Wt, H, R, options, [], epoch, grad_calc_count, 0);

    if options.verbose > 1
        fprintf('%s: Epoch = 0000, cost = %.16e, optgap = %.4e\n', method_name, f_val, optgap);
    end

    % set start time
    start_time = tic();

    % main outer loop
    while true

        % check stop condition
        [stop_flag, reason, max_reached_flag] = check_stop_condition(epoch, infos, options);
        if stop_flag
            display_stop_reason(epoch, infos, options, method_name, reason, max_reached_flag);
            break;
        end

        % Reset sufficient statistic
        At = zeros(m, rank);
        Bt = zeros(rank, rank);
        Ct = zeros(m, rank);

        % main inner loop
        for t = 1 : options.batch_size : n - 1

            % Retrieve vt, ht and rt
            vt = V(:, t:t+options.batch_size-1);
            ht = H(:, t:t+options.batch_size-1);
            rt = R(:, t:t+options.batch_size-1);

            % update ht/rt
            ht = ht .* (Wt.' * vt) ./ (Wt.' * (Wt * ht + rt));
            ht = ht + (ht<eps) .* eps;
            rt = rt .* vt ./ (Wt * ht + rt + l);

            % update sufficient statistics
            At = At + vt *  ht';
            Bt = Bt + ht *  ht';
            Ct = Ct + rt *  ht';

            % Update Wt
            Wt = Wt .* At ./ (Wt * Bt + Ct);
            Wt = Wt + (Wt<eps) .* eps;

            % Update H
            H(:,t:t+options.batch_size-1) = ht;

            % Update R
            R(:,t:t+options.batch_size-1) = rt;

            grad_calc_count = grad_calc_count + m * options.batch_size;
        end

        % measure elapsed time
        elapsed_time = toc(start_time);

        % update epoch
        epoch = epoch + 1;

        % store info
        infos = store_nmf_info(V, Wt, H, R, options, infos, epoch, grad_calc_count, elapsed_time);

        % display info
        display_info(method, epoch, infos, options);

    end

    x.W = Wt;
    x.H = H;
    x.R = R;

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function [transform, coordsys] = ft_headcoordinates(fid1, fid2, fid3, fid4, coordsys)

% FT_HEADCOORDINATES returns the homogeneous coordinate transformation matrix
% that converts the specified fiducials in any coordinate system (e.g. MRI)
% into the rotated and translated headcoordinate system.
%
% Use as
%   [transform, coordsys] = ft_headcoordinates(fid1, fid2, fid3, coordsys)
% or
%   [transform, coordsys] = ft_headcoordinates(fid1, fid2, fid3, fid4, coordsys)
%
% Depending on the desired coordinate system, the order of the fiducials is
% interpreted as follows
%
%   fid1 = nas
%   fid2 = lpa
%   fid3 = rpa
%   fid4 = extra point (optional)
%
%   fid1 = ac
%   fid2 = pc
%   fid3 = midsagittal
%   fid4 = extra point (optional)
%
%   fid1 = pt1
%   fid2 = pt2
%   fid3 = pt3
%   fid4 = extra point (optional)
%
%   fid1 = bregma
%   fid2 = lambda
%   fid3 = midsagittal
%   fid4 = extra point (optional)
%
% The fourth argument fid4 is optional and can be specified as an an extra point
% which is assumed to have a positive Z-coordinate. It will be used to ensure correct
% orientation of the Z-axis (ctf, 4d, bti, eeglab, yokogawa, neuromag, itab) or
% X-axis (acpc, spm, mni, tal). The specification of this extra point may result in
% the handedness of the transformation to be changed, but ensures consistency with
% the handedness of the input coordinate system.
%
% The coordsys input argument is a string that determines how the location of the
% origin and the direction of the axis is to be defined relative to the fiducials:
%   according to CTF conventions:             coordsys = 'ctf'
%   according to 4D conventions:              coordsys = '4d' or 'bti'
%   according to EEGLAB conventions:          coordsys = 'eeglab'
%   according to NEUROMAG conventions:        coordsys = 'itab'
%   according to ITAB conventions:            coordsys = 'neuromag'
%   according to YOKOGAWA conventions:        coordsys = 'yokogawa'
%   according to ASA conventions:             coordsys = 'asa'
%   according to FTG conventions:             coordsys = 'ftg'
%   according to ACPC conventions:            coordsys = 'acpc'
%   according to SPM conventions:             coordsys = 'spm'
%   according to MNI conventions:             coordsys = 'mni'
%   according to Talairach conventions:       coordsys = 'tal'
%   according to PAXINOS conventions:         coordsys = 'paxinos'
% If the coordsys input argument is not specified, it will default to 'ctf'.
%
% The CTF, 4D, YOKOGAWA and EEGLAB coordinate systems are defined as follows:
%   the origin is exactly between lpa and rpa
%   the X-axis goes towards nas
%   the Y-axis goes approximately towards lpa, orthogonal to X and in the plane spanned by the fiducials
%   the Z-axis goes approximately towards the vertex, orthogonal to X and Y
%
% The TALAIRACH, SPM and ACPC coordinate systems are defined as:
%   the origin corresponds with the anterior commissure
%   the Y-axis is along the line from the posterior commissure to the anterior commissure
%   the Z-axis is towards the vertex, in between the hemispheres
%   the X-axis is orthogonal to the midsagittal-plane, positive to the right
%
% The NEUROMAG and ITAB coordinate systems are defined as follows:
%   the X-axis is from the origin towards the RPA point (exactly through)
%   the Y-axis is from the origin towards the nasion (exactly through)
%   the Z-axis is from the origin upwards orthogonal to the XY-plane
%   the origin is the intersection of the line through LPA and RPA and a line orthogonal to L passing through the nasion
%
% The ASA coordinate system is defined as follows:
%   the origin is at the orthogonal intersection of the line from rpa-lpa and the line trough nas
%   the X-axis goes towards nas
%   the Y-axis goes through rpa and lpa
%   the Z-axis goes approximately towards the vertex, orthogonal to X and Y
%
% The FTG coordinate system is defined as:
%   the origin corresponds with pt1
%   the x-axis is along the line from pt1 to pt2
%   the z-axis is orthogonal to the plane spanned by pt1, pt2 and pt3
%
% The PAXINOS coordinate system is defined as:
%   the origin is at bregma
%   the x-axis extends along the Medial-Lateral direction, with positive towards the right
%   the y-axis points from dorsal to ventral, i.e. from inferior to superior
%   the z-axis passes through bregma and lambda and points from cranial to caudal, i.e. from anterior to posterior
%
% See also FT_ELECTRODEREALIGN, FT_VOLUMEREALIGN, FT_INTERACTIVEREALIGN, FT_AFFINECOORDINATES, COORDSYS2LABEL

% Copyright (C) 2003-2022, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% figure out the input arguments
if nargin==3
  coordsys = 'ctf'; % default
  fid4     = [];
elseif nargin==4 && ischar(fid4)
  coordsys = fid4;
  fid4     = [];
elseif nargin==4 && isnumeric(fid4)
  coordsys = 'ctf'; % default
elseif nargin==5
  % do nothing
else
  ft_error('incorrect specification of input parameters');
end

if isnumeric(coordsys)
  % these are for backward compatibility, but should preferably not be used any more
  if coordsys==0
    coordsys = 'ctf';
  elseif coordsys==1
    coordsys = 'asa';
  elseif coordsys==2
    coordsys = 'ftg';
  else
    ft_error('if coordsys is numeric, it should assume one of the values 0/1/2');
  end
end

% ensure that they are row vectors
fid1 = fid1(:)';
fid2 = fid2(:)';
fid3 = fid3(:)';
fid4 = fid4(:)';

assert(numel(fid1)==3, 'incorrect specification of fiducial 1');
assert(numel(fid2)==3, 'incorrect specification of fiducial 2');
assert(numel(fid3)==3, 'incorrect specification of fiducial 3');
assert(isempty(fid4) || numel(fid4)==3, 'incorrect specification of fiducial 4');

% ensure that it is lower case
coordsys = lower(coordsys);

% compute the origin and direction of the coordinate axes in MRI coordinates
switch coordsys
  case {'ctf' '4d' 'bti' 'eeglab' 'yokogawa'}
    % rename the marker points for convenience
    nas = fid1; lpa = fid2; rpa = fid3; extrapoint = fid4; clear fid*
    origin = (lpa+rpa)/2;
    dirx = nas-origin;
    dirz = cross(dirx,lpa-rpa);
    diry = cross(dirz,dirx);
    dirx = dirx/norm(dirx);
    diry = diry/norm(diry);
    dirz = dirz/norm(dirz);
  case {'asa'}
    % rename the marker points for convenience
    nas = fid1; lpa = fid2; rpa = fid3; extrapoint = fid4; clear fid*
    dirz = cross(nas-rpa, lpa-rpa);
    diry = lpa-rpa;
    dirx = cross(diry,dirz);
    dirx = dirx/norm(dirx);
    diry = diry/norm(diry);
    dirz = dirz/norm(dirz);
    origin = rpa + dot(nas-rpa,diry)*diry;
  case {'neuromag' 'itab'}
    % rename the fiducials
    nas = fid1; lpa = fid2; rpa = fid3; extrapoint = fid4; clear fid*
    dirz = cross(rpa-lpa,nas-lpa);
    dirx = rpa-lpa;
    diry = cross(dirz,dirx);
    dirx = dirx/norm(dirx);
    diry = diry/norm(diry);
    dirz = dirz/norm(dirz);
    origin = lpa + dot(nas-lpa,dirx)*dirx;
  case 'ftg'
    % rename the marker points for convenience
    pt1 = fid1; pt2 = fid2; pt3 = fid3; extrapoint = fid4; clear fid*
    origin = pt1;
    dirx = pt2-origin;
    diry = pt3-origin;
    dirz = cross(dirx,diry);
    diry = cross(dirz,dirx);
    dirx = dirx/norm(dirx);
    diry = diry/norm(diry);
    dirz = dirz/norm(dirz);
  case {'acpc' 'spm' 'mni' 'tal'}
    % rename the marker points for convenience
    ac = fid1; pc = fid2; midsagittal = fid3; extrapoint = fid4; clear fid*
    origin = ac;
    diry   = ac-pc;
    dirz   = midsagittal-ac;
    dirx   = cross(diry,dirz);
    dirz   = cross(dirx,diry);
    dirx   = dirx/norm(dirx);
    diry   = diry/norm(diry);
    dirz   = dirz/norm(dirz);
  case 'paxinos'
    bregma = fid1; lambda = fid2; midsagittal = fid3; extrapoint = fid4; clear fid*
    origin = bregma;           % bregma is slightly above the brain
    dirz = lambda-bregma;      % lambda is slightly above the brain
    diry = bregma-midsagittal; % midsagittal should be somewhere in the middle of the brain
    dirx = cross(diry,dirz);   % compute the positive x-direction, i.e. towards the right
    diry = cross(dirz,dirx);   % update the y-direction
    dirx = dirx/norm(dirx);
    diry = diry/norm(diry);
    dirz = dirz/norm(dirz);
  otherwise
    ft_error('unrecognized headcoordinate system "%s"', coordsys);
end

% use the extra point to validate that it is a right-handed coordinate system
if ~isempty(extrapoint)
  dirq = extrapoint-origin;
  dirq = dirq/norm(dirq);
  if any(strcmp(coordsys, {'ctf' '4d' 'bti' 'eeglab' 'yokogawa' 'neuromag' 'itab'}))
    phi = dirq(:)'*dirz(:);
    if sign(phi)<0
      ft_warning('the input coordinate system seems left-handed, flipping z-axis to keep the transformation matrix consistent');
      dirz = -dirz;
    end
  elseif any(strcmp(coordsys, {'acpc' 'spm' 'mni' 'tal'}))
    phi = dirq(:)'*dirx(:);
    if sign(phi)<0
      ft_warning('the input coordinate system seems left-handed, flipping x-axis to keep the transformation matrix consistent');
      dirx = -dirx;
    end
  else
    ft_warning('the extra input coordinate is not used with coordsys "%s"', coordsys);
  end
end

% compute the rotation matrix
rot = eye(4);
rot(1:3,1:3) = inv(eye(3) / [dirx; diry; dirz]);
% compute the translation matrix
tra = eye(4);
tra(1:4,4)   = [-origin(:); 1];
% combine these to compute the full homogeneous transformation matrix
transform = rot * tra;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function DrawRoute(X, Y)
for i = 2: length(X)
    plot([X(i - 1), X(i)], [Y(i - 1), Y(i)], 'b-')
    hold on
    pause(0.1)
    frame = getframe(gcf);
    imind = frame2im(frame);
    [imind, cm] = rgb2ind(imind, 256);
    if i == 2
        imwrite(imind, cm, 'test.gif', 'gif', 'Loopcount', inf, 'DelayTime', 1e-4);
    else
        imwrite(imind, cm, 'test.gif', 'gif', 'WriteMode', 'append', 'DelayTime', 1e-4);
    end

end

author: LiYangSir

repo: Smart-Algorithm


text:
function ind=findn(arr);

%FINDN   Find indices of nonzero elements.
%   I = FINDN(X) returns the indices of the vector X that are
%   non-zero. For example, I = FINDN(A>100), returns the indices
%   of A where A is greater than 100. See RELOP.
%
%   This is the same as find but works for N-D matrices using
%   ind2sub function
%
%   It does not return the vectors as the third output arguement
%   as in FIND
%
%   The returned I has the indices (in actual dimensions)
%
%   x(:,:,1)            x(:,:,2)            x(:,:,3)
%       = [ 1 2 3           =[11 12 13        =[21 22 23
%           4 5 6             14 15 16          24 25 26
%           7 8 9]            17 18 19]         27 28 29]
%
%   I=find(x==25) will return 23
%   but findn(x==25) will return 2,2,3
%
%   Also see find, ind2sub

%   Loren Shure, Mathworks Inc. improved speed on previous version of findn
%   by Suresh Joel Mar 3, 2003

in=find(arr);
sz=size(arr);
if isempty(in), ind=[]; return; end;
[out{1:ndims(arr)}] = ind2sub(sz,in);
ind = cell2mat(out);

author: Sable

repo: mcbench-benchmarks