text:
function [y,err]=golaycodec(x,enc,ext)
% GOLAYCODEC - encode/decode a binary array using the Golay code with error correction.
%
% golaycodec encodes or decodes a binary array using the Golay code.
% Words of length 12 are encoded as codewords of length 23 (or 24 for the
% extended Golay code).  The Golay code is tolerant of up to 3 errors,
% i.e. for any word of length 23 there is a unique codeword that differs
% in at most 3 positions.
%
% The Golay code can be constructed as a polynomial code using one of the
% degree 11 factors of z^23+1 (mod 2).  Here we use
%   g = 1 + z^2 + z^4 + z^5 + z^6 + z^10 + z^11
% For more details see Wikipedia etc.
%
% Usage:
%   y = golaycodec(x)
%   [y,err] = golaycodec(x)
%   ... = golaycodec(x,enc)
%   ... = golaycodec(x,enc,ext)
%
% Inputs:
%   x can be an Mx12, Mx23 or Mx24 binary (logical) array or an Mx1 integer array.
%     If x is Mx1, the binary representation of each integer is used.
%   enc is a logical; if true x is to be encoded, otherwise x is to be decoded.
%     If not specified, enc is inferred from the size of x or set to true.
%   ext is a logical; if true use the extended Golay code (length 24),
%     otherwise use the non-extended Golay code (length 23).
%     If not specified, enc is inferred from the size of x or set to false.
%     NOTE: the algorithm has so far only been implemented for ext=false
%
% Outputs:
%   y is an array consisting of the coded or decoded rows of x
%     If x is Mx12 and ext=false then y is Mx23
%     If x is Mx12 and ext=true then y is Mx24
%     If x is Mx23 or Mx24 then y is Mx12
%     If x is Mx1 then y is Mx1
%   err (only available if decoding is done, e.g enc=false) is an array
%     with the errors that differ each word in x from its closest codeword)
%
% Example: encode some messages, add transmission errors and decode
%   n=10;
%   x=zeros(n,12);for i=1:n,x(i,ceil(12*rand(1,8)))=1;end; % random message
%   y=golaycodec(x);
%   err=zeros(n,23);for i=1:n,err(i,ceil(23*rand(1,3)))=1;end; % 3 random errors per message
%   y1=xor(y,err); % add transmission error
%   [x1,err1]=golaycodec(y1); % decode: should have x1==x and err1==err
%

% Author: Ben Petschel 18/3/2009
%
% Change history:
%  18/3/2009 - created
%  30/4/2009 - updated to implement extended Golay codes
%

encdef = true; % default value of "enc" if not specified for integer input
extdef = false; % default value of "ext" if not specified

% input checking: x
nx=size(x,2);
if any(~isfinite(x(:))) || any(~isreal(x(:)))
  error('golaycodec: input array must be finite real values');
end;
if nx==1
  if any(x~=round(x))
    error('golaycodec: column input must be integers');
  end;
  if any(x<0)
    % check later that integers are small enough
    error('golaycodec: column input must be non-negative integers');
  end;
  intout = true;

elseif nx>1
  if any((x(:)~=0)&(x(:)~=1))
    error('golaycodec: entries of row or matrix input must be 0 or 1');
  end;
  if ~any(nx==[12,23,24]),
    error('golaycodec: number of columns of input must be 1, 12, 23 or 24');
  end;
  intout = false; % binary output

else
  error('golaycodec: empty input');
end; % if nx==1 elseif ... else ... end;

% input checking: enc, ext must be booleans
if nargin>=2,
  if isnumeric(enc),
    enc = enc>0;
  end;
  if ~islogical(enc) || ~isscalar(enc)
    error('golaycodec: enc must be 1x1 logical value');
  end;
end;
if nargin>=3,
  if isnumeric(ext),
    ext = ext>0;
  end;
  if ~islogical(ext) || ~isscalar(ext)
    error('golaycodec: enc must be 1x1 logical value');
  end;
end;

% input checking: consistency of nx with enc and ext if specified
switch nx
  case 1
    if nargin == 1,
      enc = encdef; % set to default value (usually true)
    end;
    if nargin <= 2,
      ext = extdef; % set to default value (usually false)
    end;
  case 12
    if nargin == 1,
      enc = true; % encode 12-digit word
    elseif ~enc,
      error('golaycodec: cannot decode 12-digit words');
    end;
    if nargin <=2,
      ext = extdef; % set to default value (usually false)
    end;
  case 23
    if nargin == 1,
      enc = false; % decode 23-digit word
    elseif enc,
      error('golaycodec: cannot encode 23-digit words');
    end;
    if nargin <=2,
      ext = false; % not from the extended code
    elseif ext,
      error('golaycodec: 23 digit words are not in the extended code');
    end;
  case 24
    if nargin == 1,
      enc = false; % decode 24-digit word
    elseif enc,
      error('golaycodec: cannot encode 24-digit words');
    end;
    if nargin <=2,
      ext = true; % use the extended code
    elseif ~ext,
      error('golaycodec: 24-digit words are not in the non-extended code');
    end;
  otherwise
    error('golaycodec: input must have 1, 12, 23 or 24 columns');
end; % switch nx

% for integer inputs, check that the integers are not too large
if intout
  if enc
    if any(x(:)>=2^12),
      error('golaycodec: integers for coding must be less than 2^12');
    end;
  elseif any(x(:)>=2^23),
    error('golaycodec: integers for decoding must be less than 2^23');
  end;
end;

% output checking: nargout>1 only for decoding
if (nargout>2) || ((nargout>1) && enc)
  error('golaycodec: can only specify two outputs when decoding');
end;

if nx==1,
  % convert integers to the appropriate length
  if enc
    nx=12;
  elseif ext
    % decode extended code
    nx=24;
  else
    % decode non-extended code
    nx=23;
  end;
  x = dec2logi(x,nx);
end; % if nx==1

% retain generating polynomial, parity check and error parity table
persistent g h errtab

if isempty(g)
  % generating polynomial: g = 1+x^2+x^4+x^5+x^6+x^10+x^11
  g=zeros(1,12);g([0,2,4,5,6,10,11]+1)=1;
end;

if isempty(h)
  % parity check polynomial is the other factor of x^23-1 (mod 2)
  % h = (1+x)*(1+x+x^5+x^6+x^7+x^9+x^11)
  h=mod(deconv([1,zeros(1,22),1],g),2);
end;

if isempty(errtab),
  % set up the error table for decoding: determine all errors of up to 3 bits
  % the error will be determined by simple table lookup of the parity check
  % the table is a sparse array indexed by the integer representation of the parity bits
  errtab = geterrtab(h,3,23);
end;

% now encode or decode the array
if enc
  y = circprod(x,g,23);
  if ext
    % if using extended code, append a checksum bit (row sums mod 2)
    y = [y,mod(sum(y,2),2)];
  end;
else
  if ext
    % drop the checksum bit if using extended code
    c = x(:,24);
    x = x(:,1:23);
  end;
  % check parity
  p = circprod(x,h,23);
  % look up parity in error table
  err = dec2logi(errtab(logi2dec(p)+1),23);
  y = zeros(size(x,1),12);
  for i=1:size(x,1),
    y(i,:) = mod(deconv(x(i,:)-err(i,:),g),2);
  end;
  if ext
    % if using extended code, append the checksum bit error
    err = [err,mod(c-sum(y,2),2)];
  end;
end;

% convert back to integers, if necessary
if intout,
  y = logi2dec(y);
  if nargout == 2,
    err = logi2dec(err);
  end;
end;

end % main function golaycodec

function t=geterrtab(h,k,n)
% produces lookup table of parity checks of all errors with weight <=k on n digits
% y is a sparse matrix: index is 1 + (h*err converted to integer)
% and y(i) is (err) converted to binary
%
% h must be row vector of length <= n
% also 0 <= k <= n.

m=length(h);
if (2^n-1)>intmax,
  error('golaycodec:geterrtab: too many digits to convert to integer');
end;
if m>n,
  error('golaycodec:geterrtab: parity check polynomial is too long');
end;

t = sparse(1,1,0,2^n,1); % initialize with 0*h+1 -> 0

for i=1:k,
  ind = nchoosek(1:n,i); % all possible combinations of i error digits
  nck = size(ind,1);
  errbin = zeros(nck,n); % error in binary format
  errbin(sub2ind([nck,n],repmat((1:nck)',i,1),ind(:)))=1; % in row j place 1's in columns ind(j,:)
  errint = logi2dec(errbin); % error in integer format
  errh = circprod(errbin,h,n); % do parity check using h
  t(logi2dec(errh)+1)=errint;
end;

end % helper function geterrtab(...)

function z=circprod(x,h,n)
% does parity check of rows of x using polynomial h
m=size(x,1);

% do circular convolution by fourier transform (fft acts on columns)
% (could also use conv and wrap the tail elements, but conv isn't vectorized)
z = mod(round(ifft(fft(x',n).*repmat(fft(h',n),1,m))'),2);

end % helper function paritycheck

function y=logi2dec(x)
% converts each row of logical array from binary to a decimal number
n=size(x,2);
y=x*(2.^(n-1:-1:0))';

end % helper function logi2dec

function y=dec2logi(x,n)
% converts vector of integers to array of logicals (binary representation)
y=rem(floor(x(:)*2.^(1-n:0)),2); % see dec2bin.m (base matlab) for details

end % helper function dec2logi


text:
%% init java, run at begin, once
% necessary to  load java files
% clear java
% javaclasspath(pwd)
%% steps:
% - load maps
% - choose algorithm (DSL: D*Lite; FDS: Focussed D*, AS: A*)
% - do initial search
% -(optional) resolve path
% - update map
% - search in updated map
% - resolve path

%% load maps
mapsNames = {'a', 'e', 'b', 'c','map01'};
maps = {};
%scalling should be 1(original map) or 4(downscalled map for better performance)
scalling = 4;
for i=1:1: length(mapsNames)
        tmp = LoadMap(  strcat(mapsNames{i}, '.png'), scalling);
        start = tmp.start';
        goal = tmp.goal';
    maps{end+1} = tmp.map;
end

%% initial search
DSL = 0;
FDS = 1;
AS = 0;
bool = 0;
clc
close all
tic
    if DSL
        state = DSLInit(start, goal, maps{1}, scalling);
        state = DSLComputePath(state);
    elseif FDS
        state = FDSInit(start, goal, maps{1}, scalling);
        state = FDSComputePath(state);
    elseif AS
        state = ASInit(start, goal, maps{1}, scalling);
        state = ASComputePath(state);
    end
    bool = 1;
toc
gah = state;
%% update map, search in update map, run twice
state.kM = 50;
close all
tic
    if DSL
        if bool
            bool = 0;
            state = DSLUpdateMap(state, maps{2});
        else
            state = DSLComputePath(state);
        end
    elseif FDS
        if bool
            bool = 0;
            state = FDSUpdateMap(state, maps{2});
            state.kM = 50;
        else
            state = FDSComputePath(state);
        end
    end
toc
%% resolve path, insert map name here if you want to get image in original size when map was downscalled
    state.path = ResolvePath(state);
    resp = PlotPath(state, scalling, 'e');

%% show A-star graph
a = state.graph(:,:,1);
a(a(:,:)==-1) = 0; % unvisited
a(state.map(:,:)==0) = 0.2; % obstacles
a(a(:,:)==inf) = 0;
a(a(:,:)>=1) = 0.3;
imshow(a)
%% show D-star graph
% figure
% obstacles 0.2
% visited : 0.3
% unavailble, unvisited: 0
% in queue: 1
% in queue and g-val == inf: 0.7
a = state.graph(:,:,1);
a(a(:,:)~=inf) = 0.3;
a(a(:,:)==inf) = 0;
a(state.map(:,:)==0) = 0.2;
b = state.graph(:,:,3);
b( state.map(:,:)==0 ) = 0;
b( state.graph(:,:,2)==inf ) = 0;
b = b*0.7;
imshow(b+a, 'Border', 'tight')

author: LazyFalcon


text:
HEURISTICS = 0;

%%%%%%%%%%%%%%%%%%%%%%%%
%     Normalize        %
%%%%%%%%%%%%%%%%%%%%%%%%

Normalize = @(x) (x - min(x(:)))/(max(x(:)) - min(x(:)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Input images        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  Img = 'lena';
Img = 'clena';

  switch lower(Img)

    case{'lena'}
      Ig = double( imread('gray_imgs/lena_gray_512.png') ) / 255;

    case{'peppers'}
      Ig = double( imread('gray_imgs/peppers_gray.png') ) / 255;

    case{'goldhl.'}
      Ig = double( imread('gray_imgs/goldhill_gray.png') ) / 255;

    case{'clena'}
      Ig = double( imread('color_imgs/lena_color_512.png') ) / 255;

    case{'none'}
      disp(' ');
      disp('Select a test image (see code for an example)...');
      disp(' ');
      disp('Exiting test_adaptL1 code...');
      return;

    otherwise
      error('Not a valid image\n');

  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %       noisy images        %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      Ig_01L1 = imnoise(Ig, 'salt & pepper', 0.1);
      Ig_02L1 = imnoise(Ig, 'salt & pepper', 0.2);
      Ig_03L1 = imnoise(Ig, 'salt & pepper', 0.3);
      Ig_05L1 = imnoise(Ig, 'salt & pepper', 0.5);
      Ig_07L1 = imnoise(Ig, 'salt & pepper', 0.7);
      Ig_09L1 = imnoise(Ig, 'salt & pepper', 0.9);

      Ig_05L2 = imnoise(Ig,'gaussian', 0, ((5/255)^2) );
      Ig_10L2 = imnoise(Ig,'gaussian', 0, ((10/255)^2) );
      Ig_20L2 = imnoise(Ig,'gaussian', 0, ((20/255)^2) );
      Ig_30L2 = imnoise(Ig,'gaussian', 0, ((30/255)^2) );

      InCell{1} = imnoise(Ig_10L2, 'salt & pepper', 0.2);
      InCell{2} = imnoise(Ig_20L2, 'salt & pepper', 0.2);
      InCell{3} = imnoise(Ig_30L2, 'salt & pepper', 0.2);
      InCell{4} = imnoise(Ig_10L2, 'salt & pepper', 0.4);

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

      InCell{5} = imnoise(Ig_05L2, 'salt & pepper', 0.3);
      InCell{6} = imnoise(Ig_05L2, 'salt & pepper', 0.5);

      InCell{7} = imnoise(Ig_10L2, 'salt & pepper', 0.3);
      InCell{8} = imnoise(Ig_10L2, 'salt & pepper', 0.5);

%        snr_noisy12 = snr(Ig, In);

%      figure; imagesc( Normalize(In) );
%      axis image; axis off; colormap gray;
%      title(sprintf('Noisy - L1/L2. SNR: %4.1fdB.\n ', ...
%                    snr_noisy12));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% pre-processing

InputDims = size(Ig);

if(length(InputDims) == 2)
  lambdaS_l2 = [ 0.0185, 0.03, 0.045, 0.0185, 0.0075, 0.0075, 0.0185, 0.0185 ];
else
  lambdaS_l2 = 1.5*[ 0.0185, 0.03, 0.045, 0.0185, 0.0075, 0.0075, 0.0185, 0.0185 ];
end

for k=1:length(InCell)

In = InCell{k};

  tstart = tic;

S0 = adaptMedian(In);

if(length(InputDims) == 2)

  if(HEURISTICS)
    lambda_l1l2(:,:,1)  =  0.85*(S0 == 3) + ...
                    1.1*(S0 == 5) + 1.2*(S0 == 7) + 1.4*(S0 == 9);
  else
    lambda_l1l2(:,:,1)  = 1.5*(S0 ~= 0);
  end

  lambda_l1l2(:,:,2) = lambdaS_l2(k)*(S0==0);

  planes=1;
else

  if(HEURISTICS)

    for d=1:InputDims(3),
    lambda_l1l2(:,:,d)  =  1.05*(S0(:,:,d) == 3) + ...
                    1.3*(S0(:,:,d) == 5) + 1.4*(S0(:,:,d) == 7) + 1.6*(S0(:,:,d) == 9);
    end
  else

    for d=1:InputDims(3),
      lambda_l1l2(:,:,d)  =  1.75*(S0(:,:,d) ~= 0);
    end
  end

  for d=1:InputDims(3),
    lambda_l1l2(:,:,3+d) = lambdaS_l2(k)*(S0(:,:,d)==0);
  end

  planes=3;
end

% ---------- Setup for IRN_L1L2  ------------------

pars_irn_adapt = irntvInputPars('l1tv_l2tv');

pars_irn_adapt.adapt_epsR   = 1;
pars_irn_adapt.epsR_cutoff  = 0.01;
pars_irn_adapt.adapt_epsF   = 1;
pars_irn_adapt.epsF_cutoff  = 0.05;

pars_irn_adapt.pcgtol_ini = 1e-4;

pars_irn_adapt.loops      = 1;

pars_irn_adapt.U0      = In;   % necessary the for adapt case
%  pars_irn_adapt.U0      = lmu;   % necessary the for adapt case

  Ig_L1L2 = irntv(In, {}, lambda_l1l2, pars_irn_adapt);

loops = 6;
%
for m = 2:loops

   for d=1:InputDims(3)
      lambda_l1l2(:,:,3+d) = (0.975^(m-1))*lambda_l1l2(:,:,3+d);
   end
   Ipass = Ig_L1L2;
   dI = Ipass - In;
   lambda_l1l2(:,:,1:planes) = adaptLambda(dI,  lambda_l1l2(:,:,1:planes), 0.5, S0, 0.8);
   pars_irn_adapt.U0      = Ipass;
   Ig_L1L2 = irntv(In, {}, lambda_l1l2, pars_irn_adapt);

end
l1l2_time(k) = toc(tstart)
%
figure; imagesc( Normalize(Ig_L1L2) ); axis image; axis off; colormap gray;

psnr_l1l2(k) = psnr(Ig, Ig_L1L2)

end

psnr_bing = [32.84, 30.23, 27.33, 29.96];
psnr_xiao = [36.20, 33.93, 33.19, 31.51];

[ psnr_l1l2' [psnr_bing'; psnr_xiao']]

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function err = getL2error3RT0(node,elem,sigma,sigmah,markedElem)
%% GETL2ERROR3RT0  L2 norm of RT0 element in 3D.
%
%  err = getL2error3RT0(node,elem,sigma,sigmah,markedElem)
%
%  The input sigma can just be a function boundle. sigmah is the flux
%  through faces. markedElem is used to compute the error in certain region
%  only.
%
%  Note that the ascend ordering of elem is used.
%
% Example
%
%   [node,elem] = cubemesh([-1,1,-1,1,-1,1],1);
%   maxIt = 3;
%   pde = mixBCdata3;
%   err = zeros(maxIt,1);
%   h = zeros(maxIt,1);
%   for i = 1:maxIt
%      [node,elem] = uniformrefine3(node,elem);
%      uI = faceinterpolate3(pde.Du,node,elem);
%      err(i) = getL2error3RT0(node,elem,pde.Du,uI);
%      h(i) = 2^(-i);
%   end
%   figure;
%   showrateh(h,err,2,'-+','|| u - u_I ||');
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Construct Data Structure
elem = sortelem3(elem);
elem2face = dof3face(elem);
% [elem2face,dofSign] = dof3RT0(elem);
NT = size(elem,1);% Ndof = max(elem2dof(:)); %N = size(node,1);
[Dlambda,volume] = gradbasis3(node,elem);
locFace = [2 3 4; 1 3 4; 1 2 4; 1 2 3];

%% Compute square of the L2 error element-wise
[lambda,w] = quadpts3(3); % quadrature order is 3
nQuad = size(lambda,1);
err = zeros(NT,1);
for p = 1:nQuad
    % quadrature points in the x-y-z coordinate
    pxyz = lambda(p,1)*node(elem(:,1),:) ...
         + lambda(p,2)*node(elem(:,2),:) ...
         + lambda(p,3)*node(elem(:,3),:) ...
         + lambda(p,4)*node(elem(:,4),:);
    sigmap = sigma(pxyz);
    sigmahp = zeros(NT,3);
    for l = 1:4 % for each basis
        i = locFace(l,1); j = locFace(l,2); k = locFace(l,3);
        % phi_l = 2(lambda_i Dlambda_j x Dlambda_k +
        %           lambda_j Dlambda_k x Dlambda_i +
        %           lambda_k Dlambda_i x Dlambda_j)
        sigmahp = sigmahp + repmat(2*sigmah(elem2face(:,l)),1,3).*...
                   (lambda(p,i)*mycross(Dlambda(:,:,j),Dlambda(:,:,k)) + ...
                    lambda(p,j)*mycross(Dlambda(:,:,k),Dlambda(:,:,i)) + ...
                    lambda(p,k)*mycross(Dlambda(:,:,i),Dlambda(:,:,j)));
    end
    err = err + w(p)*sum((sigmap - sigmahp).^2,2);
end
err = err.*volume;
% modify the error
err(isnan(err)) = 0;
if (nargin == 5) && ~isempty(markedElem)
    err = err(markedElem); % L2 err on some marked region
end
err = sqrt(sum(err));

repo: ifem


text:
function s = y_to_s_common ( y )

%*****************************************************************************80
%
%% Y_TO_S_COMMON writes a Common year into a string.
%
%  Format:
%
%    YearNumber BCE
%    YearNumber CE
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 October 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year.
%
%    Output, string S, a representation of the year.
%
  if ( y < 0 )
    s = sprintf ( 'BCE %d', y );
  else
    s = sprintf ( 'CE %d', y );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


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

%%
clear; close all; clc;

%% Syntax
% |[A]=dihedralAngles(varargin);|

%% Description
% UNDOCUMENTED
%% Examples
%

%%
% Plot settings

fontSize=20;
faceAlpha1=0.8;

%% CREATING A MESHED BOX
boxDim=2*ones(1,3);
boxEl=3*ones(1,3);
[meshStruct]=hexMeshBox(boxDim,boxEl);
E=meshStruct.E;
V=meshStruct.V;

%%
d=eye(3,3);
d(1,2)=1;

V=V*d;

%%

[A,EE,AE]=dihedralAngles(E,V,'hex8');
A=180*(A./pi);
AE=180*(AE./pi);

A_max=max(A,[],2);
A_min=min(A,[],2);

[F,A_max_F]=element2patch(E,A_max);
[~,A_min_F]=element2patch(E,A_min);
%%

cFigure;
subplot(1,2,1); hold on;
title(['Max dihedral angle ',num2str(max(A_max_F))])
gpatch(F,V,A_max_F,'k',1,1);
axisGeom; camlight headlight;
colormap(gca,gjet(25)); colorbar;
clim([min(A(:)) max(A(:))]);

subplot(1,2,2); hold on;
title(['Min dihedral angle ',num2str(min(A_min_F))])
gpatch(F,V,A_min_F,'k',1,1);
axisGeom; camlight headlight;
colormap(gca,gjet(25)); colorbar;
clim([min(A(:)) max(A(:))]);
gdrawnow;

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function result = subaperturemem(im, varargin)
%SUBAPERTUREMEM subaperture processed imaging
%    result = subaperturemem(compleximage, 'PropertyName', PropertyValue, ...)
%
% Calculates the subaperture-processed image on complex data that is held in memory
%
%       Property name     Description
%       frames            number of frames (default = 7)
%       apfraction        fraction of aperture for each subaperture
%                            (default = .25)
%       method            'normal' (default), 'fullpixel', or 'minimal'
%       platformdir       platform direction, 'right' (default) or 'left'
%       dim               dimension over which to split subaperture
%                            (default = 1)
%       fill              fill factor (default = 1)
%
% Output is stored in a cell array, where each element of the array is on
% frame.
%
% Limitation: Assumes complex data and all frames can be held in
% memory at once.
%
% Written by: Tom Braun
%             Wade Schwartzkopf
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%% Parse and validate arguments if not already done
if isreal(im) || (~isnumeric(im)) || (ndims(im) > 2)
    error('Input image must be a single complex image');
end
if (nargin>1)&&isstruct(varargin{1})
    inargs=varargin{1};
else
    inargs=parsesubapertureinputs(varargin{:});
end

%% Compute outpute size
num_x=size(im,inargs.dim); % Length along processing direction
if strcmp(inargs.method, 'minimal')
    inargs.output_res = ceil(num_x ./ inargs.fill ./ inargs.frames);
elseif strcmp(inargs.method, 'fullpixel')
    inargs.output_res = num_x;
else
    inargs.output_res = ceil(inargs.apfraction*num_x);
end

%% Setup for the subaperture processing
if (inargs.dim==2), im = im.'; end;
IM = fft(im);
clear im;
cutoff = floor(size(IM,1) ./ inargs.fill ./ 2);
IM = [IM(end-cutoff+1:end, :); IM(1:cutoff, :)]; % FFTSHIFT + chop off zeropad
result = cell(inargs.frames,1);
num_sa = ceil(inargs.apfraction * 2 * cutoff);
step = floor((2 * cutoff - num_sa)./(inargs.frames - 1));
offset = floor((2 * cutoff - (step*(inargs.frames-1)+num_sa)) ./ 2);

%% Run the subaperture processing.  We do it a frame at a time instead of a line at a time -- either works
for f = 1:inargs.frames
    if inargs.reverse_frames
        frame_num = inargs.frames - f + 1;
    else
        frame_num = f;
    end
    result{frame_num} = ifft(IM(offset + step*(f-1) + (1:num_sa), :), inargs.output_res);
    if inargs.dim == 2
        result{frame_num}=result{frame_num}.';
    end
end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
classdef AnalyzingAmplificators < handle

    properties (Access = private)
       homog
       m1
       m2
    end

    properties (Access = private)
        fileName
        pNorm
        path
    end

    methods (Access = public)

        function obj = AnalyzingAmplificators()
            obj.init();
            obj.createM1M2();
            obj.createHomogenizedComputer();
            obj.computeAmplificators()
            obj.plot();
        end

    end

    methods (Access = private)

        function init(obj)
            obj.fileName = 'SuperEllipseQOptAnalytic';
            obj.pNorm    = 16;
            obj.path     = '/home/alex/Dropbox/GregoireMeeting19January/';
        end

        function createM1M2(obj)
            npoints = 100;
            obj.m1 = 0.5*ones(npoints,1);
            obj.m2 = linspace(0.8,0.97,npoints)';
        end

        function createHomogenizedComputer(obj)
           s.type = 'ByVademecum';
           s.fileName = obj.fileName;
           h = HomogenizedVarComputer.create(s);
           obj.homog = h;
        end

        function computeAmplificators(obj)
            x{1} = obj.m1;
            x{2} = obj.m2;
            x{3} = [ones(1,length(obj.m1));zeros(1,length(obj.m1))];
            obj.homog.computePtensor(x,obj.pNorm);
        end

        function plot(obj)
            obj.plot1111();
            obj.plot2222();
            obj.plot1212();
            obj.printFigure(f,h)
        end

        function plot1212(obj)
            comp = 165;
            name = 'P1212';
            obj.plotComponent(comp,name);
        end

        function plot1111(obj)
            comp = 1287;
            name = 'P1111';
            obj.plotComponent(comp,name);
        end

        function plot2222(obj)
            comp = 495;
            name = 'P2222';
            obj.plotComponent(comp,name);
        end

        function plotComponent(obj,comp,name)
            f = figure();
            Pp = obj.homog.Pp(comp,:);
            Ppn = (Pp).^(1/obj.pNorm);
            h{1} = semilogy(obj.m2,(Ppn),'+-');
            legend(name);
            obj.printFigure(f,h,name)
        end

        function printFigure(obj,f,h,name)
            outputName = [obj.path,name,'m2_05'];
            printer = plotPrinter(f,h);
            printer.print(outputName);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function simsetup_rot
motionobject = getMotionObject();
if isempty(motionobject)
   testmotion;
   motionobject = getMotionObject();
end
end

author: Sable

repo: mcbench-benchmarks


text:
t=0:1:17
delta=2
tau=1.05
n=3;
hrf=gamma1D(t,delta,tau,n); hrf(1)=0;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Universitat Politecnica de Catalunya BarcelonaTech (UPC) - Spain
%  University of California Berkeley (UCB) - USA
%
%  Jordi Pont-Tuset <jordi.pont@upc.edu>
%  Pablo Arbelaez <arbelaez@berkeley.edu>
%  June 2014
% ------------------------------------------------------------------------
% This file is part of the MCG package presented in:
%    Arbelaez P, Pont-Tuset J, Barron J, Marques F, Malik J,
%    "Multiscale Combinatorial Grouping,"
%    Computer Vision and Pattern Recognition (CVPR) 2014.
% Please consider citing the paper if you use this code.
% ------------------------------------------------------------------------
%
% Script to compute MCG or SCG candidates on a whole dataset
%
% ------------------------------------------------------------------------
function im2mcg_all(mode, database, gt_set)
if nargin==0
    mode = 'fast';
end
if nargin<2
    database = 'pascal2012';
end
if nargin<3
    gt_set = 'val2012';
end

% Create out folder
if strcmp(mode,'fast')
    res_dir = fullfile(root_dir,'datasets',database,'SCG');
elseif strcmp(mode,'accurate')
    res_dir = fullfile(root_dir,'datasets',database,'MCG');
else
    error('Unknown mode for MCG: Possibilities are ''fast'' or ''accurate''')
end
if ~exist(res_dir,'dir')
    mkdir(res_dir);
end

% Load which images to consider from the database (train, val, etc.)
im_ids = database_ids(database,gt_set);

% Sweep all images and process them in parallel
num_workers = 4;
if isempty(gcp('nocreate'))
    parpool('local', num_workers)
end
num_images = length(im_ids);
parfor im_id = 1:num_images
    % File to store the candidates
    res_file   = fullfile(res_dir,[im_ids{im_id} '.mat']);

    % Do not recompute if already computed
    if ~exist(res_file,'file')
        % Read the image
        image = get_image(database, im_ids{im_id});

        % Call the MCG code
        candidates = im2mcg(image,mode);

        % Store the masks results
        parsave(res_file,candidates);
    end
end
delete(gcp('nocreate'))
end

function parsave(res_file,candidates)
    scores = candidates.scores; %#ok<NASGU>
    bboxes = candidates.bboxes; %#ok<NASGU>
    superpixels = candidates.superpixels; %#ok<NASGU>
    labels = candidates.labels; %#ok<NASGU>
    save(res_file,'scores','bboxes','superpixels','labels');
end

author: jponttuset

repo: mcg

save_path: github-repos/MATLAB/jponttuset-mcg


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
% ##2
%==============================================================================
%
%   - data                 3d-brains, Omega=(0,20)x(0,10)x(0,20), level=3:6, m=[128,64,128]
%   - viewer               imgmontage
%   - interpolation        linearInter
%   - distance             SSD
%   - pre-registration     affine3D
%   - regularizer          mfElastic
%   - optimizer            Trust-Region
% ===============================================================================

setup3DbrainData;

% prepare the plot
FAIRplots('clear')
Dshow = @(T,R,omega,m) viewIP(abs(T-R),omega,m,'colormap',gray(256));
FAIRplots('set','Dshow',Dshow);

% finally:
% run the MultiLevel Non-Parametric Image Registration
% Trust Region (multigrid preconditioner for conjugate gradient solver)
% no pre-registration

NPIRpara = optPara('TrustRegion','preconditioner',@MGsolver)
[yc,wc,his] = MLIR(ML,'NPIRpara',NPIRpara,'parametric',0);

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


text:
function prior = truncatedGammaPriorParamInit(prior)

% TRUNCATEDGAMMAPRIORPARAMINIT Truncated gamma prior model's parameter initialisation.
% FORMAT
% DESC initialises the parameters of the truncated gamma prior with some
% default parameters.
% ARG prior : prior structure to be initialised.
% RETURN prior : prior structure with initial values in place.
%
% SEEALSO : priorCreate
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006
%
% COPYRIGHT : Antti Honkela, 2013

% PRIOR

prior.a = 1e-6;
prior.b = 1e-6;
prior.lbound = 0;
prior.ubound = 1;

prior.transforms.index = [1 2];
prior.transforms.type = optimiDefaultConstraint('positive');
prior.nParams = 2;
prior.isBounded = 1;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function x = NormL1_project(x,weights,tau)

if isreal(x)
   x = oneProjector(x,weights,tau);
else
   xa  = abs(x);
   idx = xa < eps;
   xc  = oneProjector(xa,weights,tau);
   xc  = xc ./ xa; xc(idx) = 0;
   x   = x .* xc;
end

author: mpf


text:
%CHAR	create string representation of quaternion object

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

%	Copright (C) Peter Corke 1999
function s = char(q)

	s = [num2str(q.s), ' <' num2str(q.v(1)) ', ' num2str(q.v(2)) ', '   num2str(q.v(3)) '>'];

author: petercorke

repo: robotics-toolbox-matlab


text:
function fem2d_scalar_display ( prefix )

%*****************************************************************************80
%
%% FEM2D_SCALAR_DISPLAY creates surface plots of 2D FEM scalar data.
%
%  Discussion:
%
%    This program assumes that you have computed the value of some scalar
%    quantity (such as pressure or temperature) at a set of nodes.
%
%    You may have determined an order 3 triangulation of these nodes,
%    but if you have not, the program will work that out internally.
%
%    This program can read that data, and display a color contour of the
%    solution data.
%
%  Usage:
%
%    fem2d_scalar_display ( 'prefix' )
%
%    where
%
%    * 'prefix'_nodes.txt contains the node coordinates;
%    * 'prefix'_elements.txt contains the element definitions
%      (this file is optional, and if missing, the elements will be generated
%      by the program);
%    * 'prefix'_values.txt contains the nodal values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 November 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string PREFIX, the common file prefix.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM2D_SCALAR_DISPLAY:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Surface plot of a scalar U(X,Y) on a triangulated region.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  This program expects three input files:\n' );
  fprintf ( 1, '  * a node file,      the node coordinates,\n' );
  fprintf ( 1, '  * an element file,  triples of nodes that form elements,\n' );
  fprintf ( 1, '  * a value file,     solution values.\n' );
%
%  The command line argument is the common filename prefix.
%
  if ( nargin < 1 )

    fprintf ( 1, '\n' );

    prefix = input ( ...
      'Please enter the filename prefix:' );

  end
%
%  Create the filenames.
%
  node_filename = strcat ( prefix, '_nodes.txt' );
  element_filename = strcat ( prefix, '_elements.txt' );
  value_filename = strcat ( prefix, '_values.txt' );
%
%  Read the node data.
%
  [ dim_num, node_num ] = r8mat_header_read ( node_filename );

% fprintf ( 1, '\n' );
% fprintf ( 1, '  Read the header of "%s".', node_filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension DIM_NUM = %d\n', dim_num );
  fprintf ( 1, '  Number of points NODE_NUM = %d\n', node_num );

  if ( dim_num ~= 2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FEM2D_SCALAR_DISPLAY - Fatal error!\n' );
    fprintf ( 1, '  Dataset must have spatial dimension 2.\n' );
    error ( 'FEM2D_SCALAR_DISPLAY - Fatal error!' );
  end

  node_xy = r8mat_data_read ( node_filename, dim_num, node_num );

% fprintf ( 1, '\n' );
% fprintf ( 1, '  Read the data in "%s".\n', node_filename );

% r8mat_transpose_print_some ( dim_num, node_num, node_xy, 1, 1, dim_num, 5, ...
%   '  First 5 nodes:' );
%
%  Read or create the element data.
%
  if ( file_exist ( element_filename ) )

    [ element_order, element_num ] = i4mat_header_read ( ...
      element_filename );

    if ( element_order ~= 3 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'FEM2D_SCALAR_DISPLAY - Fatal error!\n' );
      fprintf ( 1, '  Data is not for a 3-node triangulation.\n' );
      error ( 'FEM2D_SCALAR_DISPLAY - Fatal error!' );
    end

%   fprintf ( 1, '\n' );
%   fprintf ( 1, '  Read the header of "%s".\n', ...
%     element_filename );
%   fprintf ( 1, '\n' );
    fprintf ( 1, '  Element order = %d\n', element_order );
    fprintf ( 1, '  Number of elements ELEMENT_NUM  = %d\n', ...
      element_num );

    element_node = i4mat_data_read ( element_filename, ...
      element_order, element_num );

%   fprintf ( 1, '\n' );
%   fprintf ( 1, '  Read the data in "%s".\n', element_filename );

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Creating triangulation for data.\n' );
    element_node = delaunayn ( node_xy' );
    element_node = element_node';
    [ element_order, element_num ] = size ( element_node );
    i4mat_write ( element_filename, element_order, element_num, element_node );
    fprintf ( 1, '  Triangulation data written to "%s".\n', element_filename );
  end

% i4mat_transpose_print_some ( element_order, element_num, ...
%   element_node, 1, 1, element_order, 10, ...
%   '  First 10 elements:' );
%
%  Detect and correct 0-based indexing.
%
  element_node = mesh_base_one ( node_num, element_order, element_num, ...
    element_node );
%
%  Read the values.
%
  [ value_dim, value_num ] = r8mat_header_read ( value_filename );

% fprintf ( 1, '\n' );
% fprintf ( 1, '  Read the header of "%s".', value_filename );
% fprintf ( 1, '\n' );
% fprintf ( 1, '  Spatial dimension = %d\n', value_dim );
% fprintf ( 1, '  Number of values  = %d\n', value_num );

  if ( value_dim ~= 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FEM2D_SCALAR_DISPLAY - Fatal error!\n' );
    fprintf ( 1, '  VALUE data must be scalar.\n' );
    error ( 'FEM2D_SCALAR_DISPLAY - Fatal error!' );
  end

  value = r8mat_data_read ( value_filename, value_dim, value_num );
%
%  Call TRISURF to plot the data.
%
  trisurf ( element_node', node_xy(1,:)', node_xy(2,:)', value', ...
    'Edgecolor', 'None' )

  xlabel ( '<---X--->', 'Fontsize', 16 );
  ylabel ( '<---Y--->', 'Fontsize', 16 );
  zlabel ( '<---U(X,Y)--->', 'Fontsize', 16 );
  title ( prefix, 'Fontsize', 24 );
%
%  Save it as a PNG file.
%
  png_filename = strcat ( prefix, '.png' );
  print ( '-dpng', png_filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saving a PNG version of plot as "%s"\n', png_filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM2D_SCALAR_DISPLAY:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function column_num = file_column_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_COLUMN_COUNT counts the columns in the first line of a file.
%
%  Discussion:
%
%    The file is assumed to be a simple text file.
%
%    Most lines of the file are presumed to consist of COLUMN_NUM words,
%    separated by spaces.  There may also be some blank lines, and some
%    comment lines, which have a "#" in column 1.
%
%    The routine tries to find the first non-comment non-blank line and
%    counts the number of words in that line.
%
%    If all lines are blanks or comments, it goes back and tries to analyze
%    a comment line.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the file.
%
%    Output, integer COLUMN_NUM, the number of columns in the file.
%
  FALSE = 0;
  TRUE = 1;
%
%  Open the file.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_COLUMN_COUNT - Error!' );
  end
%
%  Read one line, but skip blank lines and comment lines.
%  Use FGETL so we drop the newline character!
%
  got_one = FALSE;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( s_len_trim ( line ) == 0 )

    elseif ( line(1) == '#' )

    else
      got_one = TRUE;
      break;
    end

  end

  fclose ( input_unit );

  if ( got_one == FALSE )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Warning!\n' );
    fprintf ( 1, '  The file does not seem to contain any data.\n' );
    column_num = -1;
    return;
  end

  column_num = s_word_count ( line );

  return
end
function value = file_exist ( file_name )

%*****************************************************************************80
%
%% FILE_EXIST reports whether a file exists.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character FILE_NAME, the name of the file.
%
%    Output, logical FILE_EXIST, is TRUE if the file exists.
%
  fid = fopen ( file_name );

  if ( fid == -1 )
    value = 0;
  else
    fclose ( fid );
    value = 1;
  end

  return
end
function row_num = file_row_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_ROW_COUNT counts the number of row records in a file.
%
%  Discussion:
%
%    Each input line is a "RECORD".
%
%    The records are divided into three groups:
%
%    * BLANK LINES (nothing but blanks)
%    * COMMENT LINES (begin with a '#')
%    * DATA RECORDS (anything else)
%
%    The value returned by the function is the number of data records.
%
%    By the way, if the MATLAB routine FGETS is used, instead of
%    FGETL, then the variable LINE will include line termination
%    characters, which means that a blank line would not actually
%    have zero characters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 December 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the input file.
%
%    Output, integer ROW_NUM, the number of rows found.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_ROW_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_ROW_COUNT - Error!' );
  end

  blank_num = 0;
  comment_num = 0;
  row_num = 0;

  record_num = 0;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    record_num = record_num + 1;
    record_length = s_len_trim ( line );

    if ( record_length <= 0 )
      blank_num = blank_num + 1;
    elseif ( line(1) == '#' )
      comment_num = comment_num + 1;
    else
      row_num = row_num + 1;
    end

  end

  fclose ( input_unit );

  return
end
function table = i4mat_data_read ( input_filename, m, n )

%*****************************************************************************80
%
%% I4MAT_DATA_READ reads data from an I4MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer M, N, the number of rows and columns in the data.
%
%    Output, integer TABLE(M,N), the point coordinates.
%
  table = zeros ( m, n );
%
%  Build up the format string for reading M real numbers.
%
  string = ' ';

  for i = 0 : m
    string = strcat ( string, ' %d' );
  end

  input_unit = fopen ( input_filename );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_DATA_READ - Error!\n' );
    fprintf ( 1, '  Could not open the input file.\n' );
    error ( 'I4MAT_DATA_READ - Error!' );
  end

  i = 0;

  while ( i < n )

    line = fgets ( input_unit );

    if ( line == -1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'I4MAT_DATA_READ - Error!\n' );
      fprintf ( 1, '  End of input while reading data.\n' );
      error ( 'I4MAT_DATA_READ - Error!' );
    end

    if ( line(1) == '#' )

    elseif ( s_len_trim ( line ) == 0 )

    else

      [ x, count ] = sscanf ( line, string );

      if ( count == m )
        i = i + 1;
        table(1:m,i) = x(1:m);
      end

    end

  end

  fclose ( input_unit );

  return
end
function [ m, n ] = i4mat_header_read ( input_filename )

%*****************************************************************************80
%
%% I4MAT_HEADER_READ reads the header from an I4MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Output, integer M, the spatial dimension.
%
%    Output, integer N, the number of points.
%
  m = file_column_count ( input_filename );

  if ( m <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data columns in\n' );
    fprintf ( 1, '  the file %s.\n', input_filename );
  end

  n = file_row_count ( input_filename );

  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data rows in\n' );
    fprintf ( 1, '  the file %s\n', input_filename );
  end

  return
end
function i4mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% I4MAT_TRANSPOSE_PRINT_SOME prints some of an I4MAT, transposed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 June 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, integer A(M,N), an M by N matrix to be printed.
%
%    Input, integer ILO, JLO, the first row and column to print.
%
%    Input, integer IHI, JHI, the last row and column to print.
%
%    Input, string TITLE, a title.
%
  incx = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );

  for i2lo = max ( ilo, 1 ) : incx : min ( ihi, m )

    i2hi = i2lo + incx - 1;
    i2hi = min ( i2hi, m );
    i2hi = min ( i2hi, ihi );

    inc = i2hi + 1 - i2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row: ' );
    for i = i2lo : i2hi
      fprintf ( 1, '%7d  ', i );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col\n' );
    fprintf ( 1, '\n' );

    j2lo = max ( jlo, 1 );
    j2hi = min ( jhi, n );

    for j = j2lo : j2hi

      fprintf ( 1, '%5d  ', j );
      for i2 = 1 : inc
        i = i2lo - 1 + i2;
        fprintf ( 1, '%7d  ', a(i,j) );
      end
      fprintf ( 1, '\n' );

    end

  end

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

%*****************************************************************************80
%
%% I4MAT_WRITE writes an I4MAT file.
%
%  Discussion:
%
%    An I4MAT is an array of I4's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 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, integer TABLE(M,N), the points.
%
%    Input, logical HEADER, is TRUE if the header is to be included.
%

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

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'I4MAT_WRITE - Error!' );
  end
%
%  Write the data.
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %12d', round ( table(i,j) ) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function element_node = mesh_base_one ( node_num, element_order, ...
  element_num, element_node )

%*****************************************************************************80
%
%% MESH_BASE_ONE ensures that the element definition is one-based.
%
%  Discussion:
%
%    The ELEMENT_NODE array contains nodes indices that form elements.
%    The convention for node indexing might start at 0 or at 1.
%    Since a MATLAB program will naturally assume a 1-based indexing, it is
%    necessary to check a given element definition and, if it is actually
%    0-based, to convert it.
%
%    This function attempts to detect 0-based node indexing and correct it.
%
%    Thanks to Feifei Xu for pointing out that I was subtracting 1 when I
%    should have been adding 1!  29 November 2012.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 November 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer ELEMENT_ORDER, the order of the elements.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input/output, integer ELEMENT_NODE(ELEMENT_ORDE,ELEMENT_NUM), the element
%    definitions.
%
  node_min = min ( min ( element_node(1:element_order,1:element_num) ) );
  node_max = max ( max ( element_node(1:element_order,1:element_num) ) );

  if ( node_min == 0 && node_max == node_num - 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MESH_BASE_ONE:\n' );
    fprintf ( 1, '  The element indexing appears to be 0-based!\n' );
    fprintf ( 1, '  This will be converted to 1-based.\n' );
    element_node(1:element_order,1:element_num) = ...
      element_node(1:element_order,1:element_num) + 1;
  elseif ( node_min == 1 && node_max == node_num )
%   fprintf ( 1, '\n' );
%   fprintf ( 1, 'MESH_BASE_ONE:\n' );
%   fprintf ( 1, '  The element indexing appears to be 1-based!\n' );
%   fprintf ( 1, '  No conversion is necessary.\n' );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MESH_BASE_ONE - Warning!\n' );
    fprintf ( 1, '  The element indexing is not of a recognized type.\n' );
    fprintf ( 1, '  NODE_MIN = %d\n', node_min );
    fprintf ( 1, '  NODE_MAX = %d\n', node_max );
    fprintf ( 1, '  NODE_NUM = %d\n', node_num );
  end

  return
end
function table = r8mat_data_read ( input_filename, m, n )

%*****************************************************************************80
%
%% R8MAT_DATA_READ reads data from an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer M, N, the number of rows and columns of data.
%
%    Output, real TABLE(M,N), the point coordinates.
%
  table = zeros ( m, n );
%
%  Build up the format string for reading M real numbers.
%
  string = ' ';

  for i = 0 : m
    string = strcat ( string, ' %f' );
  end

  input_unit = fopen ( input_filename );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_DATA_READ - Error!\n' );
    fprintf ( 1, '  Could not open the file.\n' );
    error ( 'R8MAT_DATA_READ - Error!' );
  end

  i = 0;

  while ( i < n )

    line = fgets ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( line(1) == '#' )

    elseif ( s_len_trim ( line ) == 0 )

    else

      [ x, count ] = sscanf ( line, string );

      if ( count == m )
        i = i + 1;
        table(1:m,i) = x(1:m);
      end

    end

  end

  fclose ( input_unit );

  return
end
function [ m, n ] = r8mat_header_read ( input_filename )

%*****************************************************************************80
%
%% R8MAT_HEADER_READ reads the header from an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Output, integer M, the spatial dimension.
%
%    Output, integer N, the number of points.
%
  m = file_column_count ( input_filename );

  if ( m <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data columns in\n' );
    fprintf ( 1, '  the file %s.\n', input_filename );
  end

  n = file_row_count ( input_filename );

  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data rows in\n' );
    fprintf ( 1, '  the file %s\n', input_filename );
  end

  return
end
function r8mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, real A(M,N), an M by N matrix to be printed.
%
%    Input, integer ILO, JLO, the first row and column to print.
%
%    Input, integer IHI, JHI, the last row and column to print.
%
%    Input, string TITLE, an optional title.
%
  incx = 5;

  if ( 0 < s_len_trim ( title ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, '%s\n', title );
  end

  for i2lo = max ( ilo, 1 ) : incx : min ( ihi, m )

    i2hi = i2lo + incx - 1;
    i2hi = min ( i2hi, m );
    i2hi = min ( i2hi, ihi );

    inc = i2hi + 1 - i2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row: ' );
    for i = i2lo : i2hi
      fprintf ( 1, '%7d       ', i );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col\n' );

    j2lo = max ( jlo, 1 );
    j2hi = min ( jhi, n );

    for j = j2lo : j2hi

      fprintf ( 1, '%5d ', j );
      for i2 = 1 : inc
        i = i2lo - 1 + i2;
        fprintf ( 1, '%12f', a(i,j) );
      end
      fprintf ( 1, '\n' );

    end

  end

  return
end
function len = s_len_trim ( s )

%*****************************************************************************80
%
%% S_LEN_TRIM returns the length of a character string to the last nonblank.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be measured.
%
%    Output, integer LEN, the length of the string up to the last nonblank.
%
  len = length ( s );

  while ( 0 < len )
    if ( s(len) ~= ' ' )
      return
    end
    len = len - 1;
  end

  return
end
function word_num = s_word_count ( s )

%*****************************************************************************80
%
%% S_WORD_COUNT counts the number of "words" in a string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be examined.
%
%    Output, integer WORD_NUM, the number of "words" in the string.
%    Words are presumed to be separated by one or more blanks.
%
  FALSE = 0;
  TRUE = 1;

  word_num = 0;
  s_length = length ( s );

  if ( s_length <= 0 )
    return;
  end

  blank = TRUE;

  for i = 1 : s_length

    if ( s(i) == ' ' )
      blank = TRUE;
    elseif ( blank == TRUE )
      word_num = word_num + 1;
      blank = FALSE;
    end

  end

  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:
classdef SwerlingIVSqLawD
%%SWERLINGIVSQLAWD This class implements various functions and detections
%                statistics related to the Swerling IV target model when
%                used with a square law detector. All methods are static
%                (the class holds no information and does not need to be
%                instantiated).
%
%Implemented methods are: mean, var, PDF, CDF, PD4Threshold, PD4PFA, rand
%
%Swerling models are given in [1] and in Chapter 11 of [2]. The Swerling IV
%model assumes that the observed power SNR of the target varies in terms of
%a central gamma distribution with shape parameter k=2 and scale parameter
%theta=avgSNR/2, that in a pulse train for detection the target SNR
%fluctuates from pulse to pulse, and that a square law detector is used to
%incoherently integrate multiple pulses The difference from a Swerling III
%model is the target amplitude fluctuations from pulse to pulse.

%The square law detector for N samples is
%y=sum_{i=1}^N r_i^2
%where r_i is the ith real amplitude. The model for the amplitude given a
%target amplitude or power is developed in Chapter 9 of [2]. The squared
%real amplitude is r^2_i=y_{I,i}^2+y_{Q,i}^2 where y_{I,i} and y_{Q,i} are
%the in-phase and quadrature components of the filter output. Define
%Rp=2*sampPowSNR, as in Equation 9.2-34. The sampled power SNR of the
%target in the Swerling IV model is different for each pulse. The sample
%power SNR is sampPowSNR=GammaD.rand(1,2,avgSNR/2) as mentioned in [1]. The
%value y=y_{I,i}+sqrt(-1)*y_{Q,i} conditioned on the sampled power SNR is
%modeled as being distributed circularly complex Gaussian with mean
%sqrt(Rp) and variance 2. (See, Equation 9.3-35a in [2]). The variance
%being 2 simply reflects having normalized the variance on the I and Q
%components each to 1.  As in Equation 10.4-13 of [2], the distribution of
%Y=y/2 conditioned on the target SNR values is noncentral chi-squared with
%a change of variables. The value y conditioned on the target SNR values is
%noncentral chi-squared with nu=2*N degrees of freedom and lambda=N*Rp as
%the noncentrality parameter. The expressions in [1] and [2] are in terms
%of Y and not y, so a transformation has to be performed to make things in
%terms of y.
%
%Under this model, for a given average power SNR value, one can generate a
%random sample from the distribution as
% sample=0;
% for curPulse=1:N
%     sampPowSNR=GammaD.rand(1,2,avgSNR/2);
%     Rp=2*sampPowSNR;%The peak power of the signal from the SNR, Equation
%                 %9.4-8 in [2].
%     A=sqrt(Rp);%Equation 9.4-8, the amplitude of the signal
%     sample=sample+sum(abs(ComplexGaussianD.rand(1,A,2)).^2);
% end
%
%REFERENCES:
%[1] P. Swerling, "Probability of detection for fluctuating targets," The
%    RAND Corporation, Santa Monica, CA, Tech. Rep. RM-1217, 1954.
%[2] J. V. Di Franco and W. L. Rubin, Radar Detection. SciTech Publishing
%    Inc., Rayliegh, NC: 2004.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

methods(Static)
function val=mean(avgSNR,N)
%%MEAN Obtain the mean of the distribution of the detection power
%      (normalized to a noise variance of 1) under a Swerling IV model in a
%      square-law detector.
%
%INPUTS:avgSNR The average power signal to noise ratio of the target.
%        N The number of pulses that are to be incoherently added for
%          detection (in a square-law detector). In a Swerling IV model, a
%          different realization of the target power is used across each
%          pulse and the noise corrupting the pulses varies. If this
%          parameter is omitted or an empty matrix is passed, N=1 is used.
%
%OUTPUTS: val The value of the mean for the given parameters.
%
%The mean is given in Equation 11.5-20 of [1] where equation 11.4-8 relates
%the parameterization used there to an average squared amplitude. To
%relate the average squared amplitude ABar^2 to avgSNR as used in [2], we
%use avgSNR=ABar^2/2. However, the first moment from that expression has to
%be multiplied by 2, because as seen in Equation 10.4-4 of [1], the
%definition is in terms of a scaled version of the square law detector.
%
%EXAMPLE:
%Here, we validate the mean by generating random samples and comparing the
%computed mean to the sample mean.
% avgSNR=2;
% N=4;
% numSamples=1e5;
% meanVal=SwerlingIVSqLawD.mean(avgSNR,N)
% meanSampVal=mean(SwerlingIVSqLawD.rand([numSamples,1],avgSNR,N))
%One will see that both mean values are about 24.
%
%REFERENCES:
%[1] J. V. Di Franco and W. L. Rubin, Radar Detection. SciTech Publishing
%    Inc., Rayliegh, NC: 2004.
%[2] P. Swerling, "Probability of detection for fluctuating targets," The
%    RAND Corporation, Santa Monica, CA, Tech. Rep. RM-1217, 1954.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(nargin<2||isempty(N))
        N=1;
    end

    val=2*N*(1+avgSNR);
end

function val=var(avgSNR,N)
%%VAR Obtain the variance of the distribution of the detection power
%     (normalized to a noise variance of 1) under a Swerling IV model in a
%     square-law detector.
%
%INPUTS:avgSNR The average power signal to noise ratio of the target.
%        N The number of pulses that are to be incoherently added for
%          detection (in a square-law detector). In a Swerling IV model, a
%          different realization of the target power is used across each
%          pulse and the noise corrupting the pulses varies. If this
%          parameter is omitted or an empty matrix is passed, N=1 is used.
%
%OUTPUTS: val The value of the variance for the given parameters.
%
%The variance is given in Equation 11.5-22b of [1] where equation 11.4-8
%relates the parameterization used there to an average squared amplitude.
%To relate the average squared amplitude ABar^2 to avgSNR as used in [2],
%we use avgSNR=ABar^2/2. However, the variance from that expression has to
%be multiplied by 4, because as seen in Equation 10.4-4 of [1], the
%definition in is in terms of a scaled version of the square law detector.'
%
%EXAMPLE:
%Here, we validate the variance by generating random samples and comparing
%the computed variance to the sample variance.
% avgSNR=2;
% N=4;
% numSamples=1e5;
% varVal=SwerlingIVSqLawD.var(avgSNR,N)
% varSampVal=var(SwerlingIVSqLawD.rand([numSamples,1],avgSNR,N))
%One will see that both variance values are about 112.
%
%REFERENCES:
%[1] J. V. Di Franco and W. L. Rubin, Radar Detection. SciTech Publishing
%    Inc., Rayliegh, NC: 2004.
%[2] P. Swerling, "Probability of detection for fluctuating targets," The
%    RAND Corporation, Santa Monica, CA, Tech. Rep. RM-1217, 1954.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(nargin<2||isempty(N))
        N=1;
    end

    val=4*N*(1+avgSNR.*(2+avgSNR/2));
end

function val=PDF(x,avgSNR,N)
%%PDF Evaluate the scalar probability density function (PDF) of the
%     distribution of the detection power (normalized to a noise variance
%     of 1) under a Swerling IV model in a square-law detector.
%
%INPUTS: x The point or points at which the PDF should be evaluated. v>=0
%          for nonzero PDF values.
%   avgSNR The average power signal to noise ratio of the target.
%        N The number of pulses that are to be incoherently added for
%          detection (in a square-law detector). In a Swerling IV model, a
%          different realization of the target power is used across each
%          pulse and the noise corrupting the pulses varies. If this
%          parameter is omitted or an empty matrix is passed, N=1 is used.
%
%OUTPUTS: val The values of the PDF evaluated at the given points.
%
%The PDF is given in Equation 11.5-14 where equation 11.4-8 relates the
%parameterization used there to an average squared amplitude. To relate the
%average squared amplitude ABar^2 to avgSNR as used in [2], we use
%avgSNR=ABar^2/2. However the PDF is derived is based on a scaled
%square-law detector, as in Equation 10.4-4 in [1]. Thus, in implementing
%the PDF, the scaling is removed.
%
%EXAMPLE:
%Here, we validate the PDF by generating random samples and comparing the
%PDF plot with a histogram of the random samples.
% avgSNR=2;
% N=4;
% numSamples=10e3;
%
% figure(1)
% clf
% histogram(SwerlingIVSqLawD.rand([numSamples,1],avgSNR,N),'Normalization','pdf')
% hold on
% numPoints=1000;
% x=linspace(0,80,numPoints);
% vals=SwerlingIVSqLawD.PDF(x,avgSNR,N);
% plot(x,vals,'linewidth',2)
% h1=xlabel('x');
% h2=ylabel('PDF(x)');
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%One will see that the histogram matches well with the plot.
%
%REFERENCES:
%[1] J. V. Di Franco and W. L. Rubin, Radar Detection. SciTech Publishing
%    Inc., Rayliegh, NC: 2004.
%[2] P. Swerling, "Probability of detection for fluctuating targets," The
%    RAND Corporation, Santa Monica, CA, Tech. Rep. RM-1217, 1954.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(nargin<3||isempty(N))
        N=1;
    end

    %Initial change of variable to undo the scaling from Equation 10.4-4 in
    %[1].
    x=x/2;

    xBar2=avgSNR/2;
    val=zeros(size(x));
    numEls=numel(x);

    for curEl=1:numEls
        v=x(curEl);

        if(v>0)
            val(curEl)=exp((N-1)*log(v)-v/(1+xBar2)-2*N*log(1+xBar2)-gammaln(N))*hypergeometric1F1(-N,N,-xBar2/(1+xBar2)*v);
        elseif(v==0&&N==1)
            val(curEl)=1/((1+xBar2)^(2*N)*factorial(N-1));
        end
    end

    %Adjust for the change of variable to undo the scaling from Equation
    %10.4-4 in [1].
    val=val/2;
end

function val=CDF(v,avgSNR,N)
%%CDF Evaluate the scalar cumulative distribution function (CDF) of the
%     distribution of the detection power (normalized to a noise variance
%     of 1) under a Swerling IV model in a square-law detector.
%
%INPUTS: v The point or points at which the CDF should be evaluated. v>0
%          for nonzero CDF values.
%   avgSNR The average power signal to noise ratio of the target.
%        N The number of pulses that are to be incoherently added for
%          detection (in a square-law detector). In a Swerling IV model, a
%          different realization of the target power is used across each
%          pulse and the noise corrupting the pulses varies. If this
%          parameter is omitted or an empty matrix is passed, N=1 is used.
%
%OUTPUTS: val The values of the CDF evaluated at the given points.
%
%The derivation of the detection probability in [1] is as 1-the value of
%the CDF. Thus, this function just evaluated
%1-SwerlingIVSqLawD.PD4Threshold(avgSNR,v,N).
%
%REFERENCES:
%[1] P. Swerling, "Probability of detection for fluctuating targets," The
%    RAND Corporation, Santa Monica, CA, Tech. Rep. RM-1217, 1954.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(nargin<3||isempty(N))
        N=1;
    end

    val=1-SwerlingIVSqLawD.PD4Threshold(avgSNR,v,N);
    val(v<=0)=0;
end

function PD=PD4Threshold(avgSNR,thresh,N)
%%PD4THRESHOLD Determine the detection probability of a Swerling IV target
%           given its signal to noise ratio and the value of the
%           detection threshold, assuming that a square-law detector is
%           used. The noise in each sample is assumed to be Gaussian with
%           variance 1 (normalized noise power).
%
%INPUTS: avgSNR A vector or matrix of average power signal to noise ratios
%               of the target at which one wishes to evaluate the detection
%               probability.
%        thresh The scalar normalized detection threshold to use. This is
%               the threshold to use if the noise variance is 1.
%             N The number of pulses that are to be incoherently added for
%               detection (in a square-law detector). In a Swerling IV
%               model, a different realization of the target power is used
%               across each pulse and the noise corrupting the pulses
%               varies. If this parameter is omitted or an empty matrix is
%               passed, N=1 is used.
%
%OUTPUTS: PD The detection probability of the target.
%
%This function implements the sum in Equations 11.5-19 in [1]. However, the
%expression used is derived based on a scaled square-law detector, as in
%Equation 10.4-4. Thus, the threshold in this function is scaled
%appropriately.
%
%EXAMPLE:
%Here, we validate the results by comparing the PD from this function to
%the PD computed using random samples.
% avgSNR=2;
% thresh=30;
% N=4;
% PD=SwerlingIVSqLawD.PD4Threshold(avgSNR,thresh,N)
% numSamples=1e5;
% PDSamp=mean(SwerlingIVSqLawD.rand([numSamples,1],avgSNR,N)>=thresh)
%One will see that both PD values are about 0.254.
%
%REFERENCES:
%[1] J. V. Di Franco and W. L. Rubin, Radar Detection. SciTech Publishing
%    Inc., Rayliegh, NC: 2004.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(nargin<3||isempty(N))
       N=1;
    end

    thresh=thresh/2;

    R=2*avgSNR;

    logR4=log(R/4);
    logCoeff=gammaln(N+1)-N*log(1+R/4);
    sumVal=0;
    for k=0:N
        curTerm=exp(logCoeff+k*logR4+log(gammainc(thresh/(1+R/4),k+N))-gammaln(k+1)-gammaln(N-k+1));

        sumVal=sumVal+curTerm;
    end

    PD=1-sumVal;
end

function PD=PD4PFA(avgSNR,PFA,N)
%%PD4PFA Determine the detection probability of a Swerling IV target given
%        its signal to noise ratio and the probability of false alarm
%        implied by a particular unspecified detection threshold, assuming
%        that a square law detector is used.
%
%INPUTS: avgSNR A vector or matrix of average power signal to noise ratios
%               of the target at which one wishes to evaluate the detection
%               probability.
%           PFA The probability of false alarm, 0<=PFA<1.
%             N The number of pulses that are to be incoherently added for
%               detection (in a square-law detector). In a Swerling IV
%               model, a different realization of the target power is used
%               across each pulse and the noise corrupting the pulses
%               varies. If this parameter is omitted or an empty matrix is
%               passed, N=1 is used.
%
%OUTPUTS: PD The detection probability of the target, 0<=PD<=1.
%
%This function just calls PFA2SquareLawThreshold to convert the false alarm
%rate into a normalized detection threshold after which it calls
%SwerlingIVSqLawD.PD4Threshold.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(nargin<3||isempty(N))
       N=1;
    end

    thresh=PFA2SquareLawThreshold(PFA,N);
    PD=SwerlingIVSqLawD.PD4Threshold(avgSNR,thresh,N);
end

function vals=rand(NDims,avgSNR,N)
%%RAND Generate random variables representing a target whose sampled signal
%      to average noise power ratio is generated according to a Swerling IV
%      model detected by a square-law detector.
%
%INPUTS: NDims If NDims is a scalar, then rand returns an NDimsXNDims
%          matrix of random variables. If NDims=[M,N1] is a two-element row
%          vector, then rand returns an MXN1 matrix of random variables.
%   avgSNR The average power signal to noise ratio of the target.
%        N The number of pulses that are to be incoherently added for
%          detection (in a square-law detector). In a Swerling IV model, a
%          different realization of the target power is used across each
%          pulse and the noise corrupting the pulses varies. If this
%          parameter is omitted or an empty matrix is passed, N=1 is used.
%
%OUTPUTS: vals A matrix whose dimensions are determined by N of the
%              generated random variables.
%
%Following the model of [1], a sample from the central gamma distribution
%is performed to obtain the complex signal to noise ratio of every single
%one of the N pulses. Then, the square law detector output conditioned on
%the signal to noise ratio is generated directly using NonFlucSqLawD.rand
%for each pulse and the values are added.
%
%REFERENCES:
%[1] P. Swerling, "Probability of detection for fluctuating targets," The
%    RAND Corporation, Santa Monica, CA, Tech. Rep. RM-1217, 1954.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(nargin<3||isempty(N))
        N=1;
    end

    if(isscalar(NDims))
        dims=[NDims, NDims];
    else
        dims=NDims;
    end

    vals=zeros(dims);

    numEls=numel(vals);
    for curEl=1:numEls
        for curPulse=1:N
            %Equation I.2 for  the input signal-to-noise power ratio. It is
            %different for each sample in the pulse train in the Swerling
            %IV model.

            curSNR=GammaD.rand(1,2,avgSNR/2);
            vals(curEl)=vals(curEl)+NonFlucSqLawD.rand(1,curSNR,1);
        end
    end
end
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [tscr] = triscr2(pp,tt)
%TRISCR2 calc. area-len. ratios for triangles in a 2-simplex
%triangulation in the two-dimensional plane.
%   [SCR2] = TRISCR2(VERT,TRIA) returns the area-len. ratios
%   where SCR2 is a T-by-1 vector, VERT is a V-by-2 array of
%   XY coordinates, and TRIA is a T-by-3 array of
%   vertex indexing, where each row defines a triangle, such
%   that VERT(TRIA(II,1),:), VERT(TRIA(II,2),:) and VERT(
%   TRIA(II,3),:) are the coordinates of the II-TH triangle.
%
%   See also TRIAREA, TRIANG2, TRIBAL2

%   Darren Engwirda : 2017 --
%   Email           : de2363@columbia.edu
%   Last updated    : 17/01/2017

%--------------------------- compute signed area-len. ratios
    scal = 4.0 * sqrt(3.0) / 3.0;

    area = triarea(pp,tt) ;             % also error checks!

    lrms = sum((pp(tt(:,2),:)        ...
              - pp(tt(:,1),:)).^2,2) ...
         + sum((pp(tt(:,3),:)        ...
              - pp(tt(:,2),:)).^2,2) ...
         + sum((pp(tt(:,3),:)        ...
              - pp(tt(:,1),:)).^2,2) ;

    lrms =(lrms / 3.0) .^ 1.00 ;

    tscr = scal * area ./ lrms ;

end

author: dengwirda


text:
function s = convert(tree,uid)
% XMLTREE/CONVERT Converter an XML tree in a structure
%
% tree      - XMLTree object
% uid       - uid of the root of the subtree, if provided.
%             Default is root
% s         - converted structure
%__________________________________________________________________________
%
% Convert an XMLTree into a structure, when possible.
% When several identical tags are present, a cell array is used.
% The root tag is not saved in the structure.
% If provided, only the structure corresponding to the subtree defined
% by the uid UID is returned.
%__________________________________________________________________________
% Copyright (C) 2002-2015  http://www.artefact.tk/

% Guillaume Flandin
% $Id: convert.m 6480 2015-06-13 01:08:30Z guillaume $

% Exemple:
% tree = 'field1<c>field2</c>field3';
% toto = convert(xmltree(tree));
% <=> toto = struct('b',{{'field1', 'field3'}},'c','field2')

%error(nargchk(1,2,nargin));

% Get the root uid of the output structure
if nargin == 1
    % Get the root uid of the XML tree
    root_uid = root(tree);
else
    % Uid provided by user
    root_uid = uid;
end

% Initialize the output structure
% struct([]) should be used but this only works with Matlab 6
% So we create a field that we delete at the end
%s = struct(get(tree,root_uid,'name'),''); % struct([])
s = struct('deletedummy','');

%s = sub_convert(tree,s,root_uid,{get(tree,root_uid,'name')});
s = sub_convert(tree,s,root_uid,{});

s = rmfield(s,'deletedummy');

%==========================================================================
function s = sub_convert(tree,s,uid,arg)
    type = get(tree,uid,'type');
    switch type
        case 'element'
            child = children(tree,uid);
            l = {};
            ll = {};
            for i=1:length(child)
                if isfield(tree,child(i),'name')
                    ll = { ll{:}, get(tree,child(i),'name') };
                end
            end
            for i=1:length(child)
                if isfield(tree,child(i),'name')
                    name = get(tree,child(i),'name');
                    nboccur = sum(ismember(l,name));
                    nboccur2 = sum(ismember(ll,name));
                    l = { l{:}, name };
                    if nboccur || (nboccur2>1)
                        arg2 = { arg{:}, name, {nboccur+1} };
                    else
                        arg2 = { arg{:}, name};
                    end
                else
                    arg2 = arg;
                end
                s = sub_convert(tree,s,child(i),arg2);
            end
            if isempty(child)
                s = sub_setfield(s,arg{:},'');
            end
            %- saving attributes : does not work with b
            %- but ok with <c>b</c>
%             attrb = attributes(tree,'get',uid);     %-
%             if ~isempty(attrb)                      %-
%                 arg2 = {arg{:} 'attributes'};       %-
%                 s = sub_setfield(s,arg2{:},attrb);  %-
%             end                                     %-
        case 'chardata'
            s = sub_setfield(s,arg{:},get(tree,uid,'value'));
            %- convert strings into their numerical equivalent when possible
            %- e.g. string '3.14159' becomes double scalar 3.14159
%             v = get(tree,uid,'value');              %-
%             cv = str2num(v);                        %-
%             if isempty(cv)                          %-
%                 s = sub_setfield(s,arg{:},v);       %-
%             else                                    %-
%                 s = sub_setfield(s,arg{:},cv);      %-
%             end                                     %-
        case 'cdata'
            s = sub_setfield(s,arg{:},get(tree,uid,'value'));
        case 'pi'
            % Processing instructions are evaluated if possible
            app = get(tree,uid,'target');
            switch app
                case {'matlab',''}
                    s = sub_setfield(s,arg{:},eval(get(tree,uid,'value')));
                case 'unix'
                    s = sub_setfield(s,arg{:},unix(get(tree,uid,'value')));
                case 'dos'
                    s = sub_setfield(s,arg{:},dos(get(tree,uid,'value')));
                case 'system'
                    s = sub_setfield(s,arg{:},system(get(tree,uid,'value')));
                otherwise
                    try
                        s = sub_setfield(s,arg{:},feval(app,get(tree,uid,'value')));
                    catch
                        warning('[XMLTree] Unknown target application');
                    end
            end
        case 'comment'
            % Comments are forgotten
        otherwise
            warning(sprintf('Type %s unknown : not saved',get(tree,uid,'type')));
    end

%==========================================================================
function s = sub_setfield(s,varargin)
% Same as setfield but using '{}' rather than '()'
%if (isempty(varargin) | length(varargin) < 2)
%    error('Not enough input arguments.');
%end

subs = varargin(1:end-1);
for i = 1:length(varargin)-1
    if (isa(varargin{i}, 'cell'))
        types{i} = '{}';
    elseif ischar(varargin{i})
        types{i} = '.';
        subs{i} = varargin{i}; %strrep(varargin{i},' ',''); % deblank field name
    else
        error('Inputs must be either cell arrays or strings.');
    end
end

% Perform assignment
try
   s = builtin('subsasgn', s, struct('type',types,'subs',subs), varargin{end});
catch
   error(lasterr)
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function region = p14_region ( )

%*****************************************************************************80
%
%% P14_REGION returns the name of the integration region for problem 14.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string REGION, the name of the integration region.
%
  region = 'box';

  return
end


text:
function failed_corrupt_matfiles

% MEM 8gb
% WALLTIME 10:00:00
% DEPENDENCY

files = findfiles(dccnpath('/home/common/matlab/fieldtrip/data/test'));
status = true(size(files));

for i=1:length(files)
  try
    tmp = load(files{i});
    clear tmp
    fprintf('passed reading %s\n', files{i});
  catch
    status(i) = false;
    warning('failed reading %s\n', files{i});
  end
end

if any(status==false)
  error('not all files could be read into MATLAB');
end

function files = findfiles(p)
d = dir(p);
d = d([d.bytes]<4*1e9); % skip files that are too large

files = {d(~[d.isdir]).name}';
files = files(~cellfun(@isempty, regexp(files, '\.mat$'))); % only select the mat files
files = files( cellfun(@isempty, regexp(files, '^\._')));   % exclude the OS X resource forks

for i=1:length(files)
  files{i} = fullfile(p, files{i});
end

dirs  = {d( [d.isdir]).name}';
dirs  = dirs(3:end); % the first two are . and ..

for i=1:length(dirs)
  files = cat(1, files, findfiles(fullfile(p, dirs{i})));
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author: Eugenio Alcala Baselga
% Date: 02/06/2018
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function dx = dynamic_complex_model(t, states, u)

% % %     % SEAT IRI Vehicle Parameters
% % %     m       = 1500 + 2*80 + 300; % mass, kg
% % %     Area    = 2.2;      % Superficie frontal, m^2
% % %     Cd      = 0.3;      % Cd aerodinamico
% % %     g       = 9.81;     % Gravity
% % %     ro      = 1.225;    % Densidad del aire en SI
% % %     a       = 4;
% % %     b       = 1.8;
% % %     I       = m*(a*a + b*b)/12;
% % %     FzF     = m/2 * 9.81;
% % %     FzR     = m/2 * 9.81;
% % %     muy     = 0.8;

    % CVC Vehicle Parameters
    m       = vehicle.m;
    I       = vehicle.I;
    a       = vehicle.a;
    b       = vehicle.b;
    ro      = vehicle.ro;
    Cd      = vehicle.Cd;
    Area    = vehicle.Area;
    g       = vehicle.g;
    muy     = vehicle.muy;
    FzF     = vehicle.FzF;
    FzR     = vehicle.FzR;

    mu_friction  = u(3);
    F_friction = mu_friction*m*g;  % Fuerzas que se oponen

    % States
    X       = states(1);
    Y       = states(2);
    PSI     = states(3);
    V       = states(4);        % Velocity [m/s]
    ALPHAT  = states(5);        % Slip angle [rad]
    OMEGA   = states(6);        % Yaw rate [rad/s]

    DELTA  = u(2);

    % Slip angles
    ALPHAF = atan2((V * sin(ALPHAT) + a * OMEGA), (V * cos(ALPHAT))) - DELTA;
    ALPHAR = atan2((V * sin(ALPHAT) - b * OMEGA), (V * cos(ALPHAT)));

    %%% Lateral forces (non-linear model)
    FyF = Fy_Tire_Model(ALPHAF, FzF, muy);
    FyR = Fy_Tire_Model(ALPHAR, FzR, muy);

%     Cf      = 15000;
%     FyF = Cf*(DELTA-ALPHAT-(a*OMEGA/V))
%     FyR = Cf*(-ALPHAT+(b*OMEGA/V))

    FxR = u(1);

    F_drag  = 0.5*ro*Cd*Area*V*V;

    % Equations of motion:
    dx(1,1) = V * cos(ALPHAT + PSI);    % X
    dx(2,1) = V * sin(ALPHAT + PSI);    % Y
    dx(3,1) = OMEGA;                     % dPSI
    dx(4,1) =  ( FxR * cos(ALPHAT) + FyF * sin(ALPHAT - DELTA) + FyR * sin(ALPHAT)...
                - F_drag - F_friction )/m ;                          %dV
    dx(5,1) = ( - FxR * sin(ALPHAT) +...
        FyF * cos(ALPHAT - DELTA) + FyR * cos(ALPHAT) - m * V * OMEGA) / (m * V);    %dALPHAT
    dx(6,1) = (FyF * a * cos(DELTA) - FyR * b) / I;          % Theta acceleration


text:
function [L,S,obj,err,iter] = rmsc(X,lambda,opts)

% Solve the Robust Multi-view Spectral Clustering (RMSC) problem by M-ADMM
%
% min_{L,S_i} ||L||_*+lambda*\sum_i ||S_i||_1,
% s.t. X_i=L+S_i, i=1,...,m, L>=0, L1=1.
% ---------------------------------------------
% Input:
%       X       -    d*n*m tensor
%       lambda  -    >0, parameter
%       opts    -    Structure value in Matlab. The fields are
%           opts.tol        -   termination tolerance
%           opts.max_iter   -   maximum number of iterations
%           opts.mu         -   stepsize for dual variable updating in ADMM
%           opts.max_mu     -   maximum stepsize
%           opts.rho        -   rho>=1, ratio used to increase mu
%           opts.DEBUG      -   0 or 1
%
% Output:
%       L       -    d*n matrix
%       S       -    d*n*m tensor
%       obj     -    objective function value
%       err     -    residual
%       iter    -    number of iterations
%
% version 1.0 - 19/06/2016
%
% Written by Canyi Lu (canyilu@gmail.com)
%

tol = 1e-8;
max_iter = 500;
rho = 1.1;
mu = 1e-4;
max_mu = 1e10;
DEBUG = 0;

if ~exist('opts', 'var')
    opts = [];
end
if isfield(opts, 'tol');         tol = opts.tol;              end
if isfield(opts, 'max_iter');    max_iter = opts.max_iter;    end
if isfield(opts, 'rho');         rho = opts.rho;              end
if isfield(opts, 'mu');          mu = opts.mu;                end
if isfield(opts, 'max_mu');      max_mu = opts.max_mu;        end
if isfield(opts, 'DEBUG');       DEBUG = opts.DEBUG;          end

[d,n,m] = size(X);
L = zeros(d,n);
S = zeros(d,n,m);
Z = L;
Y = S;
dY = S;
Y2 = L;
iter = 0;
for iter = 1 : max_iter
    Lk = L;
    Sk = S;
    Zk = Z;
    % first super block {Z,S_i}
    [Z,nuclearnormZ] = prox_nuclear(L+Y2/mu,1/mu);
    for i = 1 : m
        S(:,:,i) = prox_l1(-L+X(:,:,i)-Y(:,:,i)/mu,lambda/mu);
    end
    % second super block {L}
    temp = (sum(X-S-Y/mu,3)+Z-Y2/mu)/(m+1);
    L = project_simplex(temp);

    for i = 1 : m
        dY(:,:,i) = L+S(:,:,i)-X(:,:,i);
    end
    dY2 = L-Z;
    chgL = max(abs(Lk(:)-L(:)));
    chgZ = max(abs(Zk(:)-Z(:)));
    chgS = max(abs(Sk(:)-S(:)));
    chg = max([chgL chgS chgZ max(abs(dY(:))) max(abs(dY2(:)))]);
    if DEBUG
        if iter == 1 || mod(iter, 10) == 0
            obj = nuclearnormZ+lambda*norm(S(:),1);
            err = sqrt(norm(dY(:))^2+norm(dY2,'fro')^2);
            disp(['iter ' num2str(iter) ', mu=' num2str(mu) ...
                    ', obj=' num2str(obj) ', err=' num2str(err)]);
        end
    end

    if chg < tol
        break;
    end
    Y = Y + mu*dY;
    Y2 = Y2 + mu*dY2;
    mu = min(rho*mu,max_mu);
end
obj = nuclearnormZ+lambda*norm(S(:),1);
err = sqrt(norm(dY(:))^2+norm(dY2,'fro')^2);

author: canyilu

repo: LibADMM-toolbox


text:
function [ a, seed ] = runs_simulate ( m, n, seed )

%*****************************************************************************80
%
%% RUNS_SIMULATE simulates a case governed by the Runs PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the parameters of the PDF.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, integer A(M+N), a sequence of M 0's and N 1's chosen
%    uniformly at random.
%
%    Output, integer SEED, a seed for the random number generator.
%
  a(1:m) = 0;
  a(m+1:m+n) = 1;

  for i = 1 : m+n-1

    [ j, seed ] = i4_uniform ( i, m+n, seed );

    k    = a(i);
    a(i) = a(j);
    a(j) = k;

  end

  return
end


text:
function value = c8_le_l2 ( x, y )

%*****************************************************************************80
%
%% C8_LE_L2 := X <= Y for complex values, and the L2 norm.
%
%  Discussion:
%
%    The L2 norm can be defined here as:
%
%      C8_NORM2(X) = sqrt ( ( real (X) )**2 + ( imag (X) )**2 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, complex X, Y, the values to be compared.
%
%    Output, logical VALUE, is TRUE if X <= Y.
%
  if ( ( real ( x ) )^2 + ( imag ( x ) )^2 ...
    <= ( real ( y ) )^2 + ( imag ( y ) )^2 )
    value = 1;
  else
    value = 0;
  end

  return
end


text:
function r = showrateh(h,err,k,opt,str)
%% SHOWRATEH rate of an err sequence
%
%  r = SHOWRATEH(h,err) finds the number r such that err = N^r and plots the
%  err vs N in loglog scale.
%
%  r = SHOWRATEH(h,err,k) finds the number r such that err(k:end)=N(k:end)^r.
%
%  The function accepts standard plotting properting. For example, r =
%  showrate(N,err,[],'r') will plot the error curve in red.
%
% See also showrate2, showresult, showmesh, showsolution
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

N = 1./h;
if (nargin<=2)
    k = 1; opt = '-*';
end
if (nargin<=3)
    opt = '-*';
end
if ~exist('str','var')
    str = 'Error';
end
r = -showrate(N,err,k,opt);
h_legend = legend(str,['C_1h^{' num2str(r,2) '}'],'LOCATION','Best');
set(h_legend,'FontSize', 14);
xlabel('log(1/h)');
title(['Rate of convergence is Ch^{' num2str(r,2) '}'],'FontSize', 14);