text:
function write_nitf_header(obj)
%WRITE_NITF_HEADER Writes the NITF file header
%
% It's not expected that the typical user would call this function, rather
% it is called as part of a bigger NITF/SICD writer.
%
% Inputs:
%       obj:        This is the 'hidden' object-specific handle.  It is
%                   similar to 'this' in Java.
%
% References:
%       MIL-STD-2500C, Department of Defense Interface Standard
%                      National Imagery Transmission Format version 2.1
%                      01 May 2006
%
% Written by: Tom Krauss, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%
% Note, the hard-coded numbers in the writeN calls are NITF field widths
% (in characters).  If the spec changes or we choose to use something other
% than NITF 2.1 these values must change.
%

% Always writing to NITF version 2.1
fwriten(obj.FID, 'NITF02.10', 9);  % FHDR and FVER
% Complexity level depends on the "size" of the image.
% Size of image data.  Note, we're converting the integer row/col counts to
% double because the multiplication may overflow normal Matlab integer
% arithmetic.  We also need to get the pixel size which is
if obj.image_data_size < (50*(1024*1024)) % Less than 50MB
    complexity = 3;
elseif obj.image_data_size < (1*(1024*1024*1024)) % Less than 1GB
    complexity = 5;
elseif obj.image_data_size < (2*(1024*1024*1024)) % Less than 2GB
    complexity = 6;
else
    complexity = 7;
end
fwriten(obj.FID, uint8(complexity),  2);  % CLEVEL
fwriten(obj.FID, 'BF01',      4);  % STYPE
fwriten(obj.FID, 'Unknown',   10); % OSTAID (not supposed to be blank)

if isfield(obj.sicdmeta,'ImageCreation')&&isfield(obj.sicdmeta.ImageCreation,'DateTime')
    fdt = datestr(obj.sicdmeta.ImageCreation.DateTime,'yyyymmddHHMMSS'); % Creation time of original image
else
    fdt = datestr(now,'yyyymmddHHMMSS'); % Creation time of this NITF
end
fwriten(obj.FID, fdt,         14); % FDT

if isfield(obj.sicdmeta,'CollectionInfo')&&isfield(obj.sicdmeta.CollectionInfo,'CoreName')
    FTITLE = obj.sicdmeta.CollectionInfo.CoreName;
else
    FTITLE = 'Unknown';
end
if obj.is_complex
    FTITLE = ['SICD: ' FTITLE];
end
fwriten(obj.FID, FTITLE, 80); % FTITLE

obj.write_nitf_security_tags(); % Security tags go here
fwriten(obj.FID,   '00000', 5);  % FSCOP - File Copy Number
fwriten(obj.FID,   '00000', 5);  % FSCPYS - File Number of Copies
fwriten(obj.FID,   '0', 1); % ENCRYP - Encryption. (0=none)
fwrite(obj.FID, 0, 'uint8'); % FBKGC 3-byte binary write for
fwrite(obj.FID, 0, 'uint8'); % background color
fwrite(obj.FID, 0, 'uint8'); % (red, green, blue)
fwriten(obj.FID,   '', 24); % ONAME
fwriten(obj.FID,   '', 18); % OPHONE

% File and header lengths

fwriten(obj.FID, uint64(obj.total_file_length), 12);   % File length
fwriten(obj.FID, uint64(obj.NITF_header_length), 6);   % Header length - fixed size for NITF 2.1

% Image Segment description
fwriten(obj.FID, uint64(1), 3);    % Always going to write a single image segment.
fwriten(obj.FID, uint64(obj.IS_subheader_size), 6);  % Length of image subheader - fixed size (512) for NITF 2.1
fwriten(obj.FID, uint64(obj.image_data_size), 10);

% Graphic segments (none)
fwriten(obj.FID, uint64(0), 3);

% Reserved extension segments (none)
fwriten(obj.FID, uint64(0), 3);

% Text Segments (none)
fwriten(obj.FID, uint64(0), 3);

% Data Extension segments.  This is where the SICD XML goes.  We'll make
% the XML string here so we can get its length even though we won't write
% the string until later.
fwriten(obj.FID, uint64(1), 3);
fwriten(obj.FID, uint64(obj.DES_header_length), 4);
fwriten(obj.FID, uint64(length(obj.DES_data)), 9);   % Length of DES data

% Reserved Extension segments (none)
fwriten(obj.FID, uint64(0), 3);

% User defined headers (none)
fwriten(obj.FID, uint64(0), 5);

% Extended headers (none)
fwriten(obj.FID, uint64(0), 5);

end

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

author: ngageoint


text:
function grafica2000nudos
%%%%%%%%%%%%%%%%%%           SAP2000 v14           %%%%%%%%%%%%%%%%%%%%%
%programa registrado por MECAVIL ingenieria-ciencia
%bresler_lin@hotmail.com
%numera los nudos
%clf
%cla
graficaTRIDIMENSIONAL
[frini frfin elmprp ielmprp conecc inud nud]=frame3dsap2000longitud;
%%%%%%%%%%%%%%%%%%%%%%
%enumeracion frames%%%
%%%%%%%%%%%%%%%%%%%%%%
%M=(frfin(:,[2 3 4])+frini(:,[2 3 4]))./2;
%axis('auto')

%for i=1:1:ielmprp
% text(M(i,1),M(i,2),M(i,3),[{'\color{blue}'} [,int2str( conecc(i,1) ),]],'FontSize',7)
%end
%%%%%%%%%%%%%%%%%%%%%%
%enumeracion nudos%%%
%%%%%%%%%%%%%%%%%%%%%%
figure(1),hold on
for i=1:1:inud
 text(nud(i,2),nud(i,3),nud(i,4),[{'\color[rgb]{0 0 0}'} [,int2str(nud(i) ),]],'FontSize',8)
end
hold off

author: Sable

repo: mcbench-benchmarks


text:
function [ detection, Q, G ] = detector( face_image, model )
%DETECTOR Given normalized image frame face_image and model structure this
%function detects positions of significant points on human face and returns
%its positions.
%
% INPUT:
%   face_image      ...     normalized image frame (see function getNormalizedFrame)
%   model           ...     structure holding parameters of model (see
%                           exampleUsageDetector.m for details how to create it)
%
% OUTPUT:
%   detection       ...     matrix 2xM format [x_0 ... x_{M-1};
%                                              y_0 ... y_{M-1}]
%                           detected positions of significant points
%
% 18-08-10 Michal Uricar
% 08-08-11 Michal Uricar, argmax_mex modification

    model.data.Images = face_image(:);

%     Psi = model.fgetPsiMat(model.data, 1);
%     Psi = model.fgetPsiMat(model.data);

%     [ q, g ] = model.fgetQG(model.data.options, model.W, Psi, model.data.mapTable);
%     detection = model.fargmax(model.data.options, q, g);

    [Psi, PsiSparse] = model.fgetPsiMat(model.data, 1);

    if (nargout > 1)
        [ q, g ] = model.fgetQG(model.data.options, model.W, Psi, model.data.mapTable);
        [detection, Q, G] = model.fargmax(model.data.options, q, g);
    else
        %detection = model.fargmax(model.data.options, q, g);
        detection = model.argmax_mex(model.data.options, model.W, model.data.mapTable, PsiSparse);
    end

end

author: uricamic

repo: flandmark

save_path: github-repos/MATLAB/uricamic-flandmark


text:
function [param]=zn3tak(input)
% [param]=zn3pd(input)
% Takahashi's controller for processes of 3rd order.
% This function computes parameters of the controller (r0, q0, q1, q2, p1, p2).
% Output of the controller is calculated follows:
%
%                  r0                            q0 + q1*z^-1 + q2*z^-2
% U(z^-1) = ----------------------- * W(z^-1) - ------------------------ * Y(z^-1)
%            1 + p1*z^-1 + p2*z^-2                1 + p1*z^-1 + p2*z^-2
%
% where p1=-1, p2=0
%
% Transfer function of the controlled system is:
%
%               b1*z^-1 + b2*z^-2 + b3*z^-3
% Gs(z^-1) = ---------------------------------
%             1 + a1*z^-1 + a2*z^-2 + a3*z^-3
%
% Input: input ... input parameters
%                  input(1) ... a1
%                  input(2) ... b1
%                  input(3) ... a2
%                  input(4) ... b2
%                  input(5) ... a3
%                  input(6) ... b3
%                  input(7) ... sample time T0
% Output: param ... controller parameters
%                   param(1) ... r0
%                   param(2) ... q0
%                   param(3) ... q1
%                   param(4) ... q2 (0)
%                   param(5) ... p1 (0)
%                   param(6) ... p2 (0)

a1 = input(1);
b1 = input(2);
a2 = input(3);
b2 = input(4);
a3 = input(5);
b3 = input(6);
T0 = input(7);

% compute ultimate gain and frequency
[Kpu, Tu] =  ultim([b1 b2 b3],[a1 a2 a3],T0);

Kp = 0.6*Kpu*(1-T0/Tu);
Ti = Kp*Tu/(1.2*Kpu);
Td = 3*Kpu*Tu/(40*Kp);

r0 = Kp*T0/Ti;
q0 = Kp*(1+T0/Ti+Td/T0);
q1 = -Kp*(1+2*Td/T0);
q2 = Kp*(Td/T0);
p1 = -1;
p2 = 0;

param=[r0; q0; q1; q2; p1; p2];


text:
function x = equidistant3 ( n )

%*****************************************************************************80
%
%% EQUIDISTANT3 returns the Type 3 Equidistant points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2018
%
%  Author:
%
%    John Burkardt.
%
%  Parameters:
%
%    Input, integer N, the number of points.
%
%    Input, real X(N), the points.
%
  x(1:n) = ( 1 - n : 2 : n - 1 ) / n;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [param]=db3w(input)
% [param] = db3w(input)
% Dead beat controller for processes of 3rd order (weak version).
% This function computes parameters of the controller (r0, q0, q1, q2, p1, p2).
% Output of the controller is calculated follows:
%
%                  r0                            q0 + q1*z^-1 + q2*z^-2
% U(z^-1) = ----------------------- * W(z^-1) - ------------------------ * Y(z^-1)
%            1 + p1*z^-1 + p2*z^-2                1 + p1*z^-1 + p2*z^-2
%
% where q2=0 and p2=0
% The controller parameters are calculated for step change of W.
%
%% Transfer function of the controlled system is:
%
%               b1*z^-1 + b2*z^-2 + b3*z^-3
% Gs(z^-1) = ---------------------------------
%             1 + a1*z^-1 + a2*z^-2 + a3*z^-3
%
% Input: input ... input parameters
%                  input(1) ... a1
%                  input(2) ... b1
%                  input(3) ... a2
%                  input(4) ... b2
%                  input(5) ... a3
%                  input(6) ... b3
% Output: param ... controller parameters
%                   param(1) ... r0
%                   param(2) ... q0
%                   param(3) ... q1
%                   param(4) ... q2
%                   param(5) ... p1
%                   param(6) ... p2

a1=input(1);
b1=input(2);
a2=input(3);
b2=input(4);
a3=input(5);
b3=input(6);

r0 = 1/b1;
q0 = -a1/b1;
q1 = -a2/b1;
q2 = -a3/b1;
p1 = b2/b1;
p2 = b3/b1;

param=[r0; q0; q1; q2; p1; p2];


text:
%%% Example script showing how to perform a 3D Total-Variation filtering with proxTV

clear all
close all

% Load color image (3 dimensions: length, width and color)
X = imread('colors.png');

% Introduce noise
noiseLevel = 0.2;
N = double(imnoise(X,'gaussian',0,noiseLevel));

% Filter using 3D TV-L1: for 3D one needs to invoke prox_TVgen
lambda=100;
disp('Filtering image...');
tic;
F = prox_TVgen(N,      [lambda lambda],             [1 2],                    [1 1]);
%              Image | Penalty in each dimension |  Dimensions to penalize  | Norms to use
toc;

% Any dimension can be penalized under any norm. By also penalizing the color dimension under TV-L2 we get a "decolored" image
lambda2=5;
disp('Color filtering...');
tic;
F2 = prox_TVgen(N,      [lambda lambda lambda2],     [1 2 3],                  [1 1 2]);
%               Image | Penalty in each dimension |  Dimensions to penalize  | Norms to use
toc;

% Plot results
figure();
subplot(2,2,1);
imshow(X);
title('Original');
subplot(2,2,2);
imshow(uint8(N));
title('Noisy');
subplot(2,2,3);
imshow(uint8(F));
title('Filtered');
subplot(2,2,4);
imshow(uint8(F2));
title('Color filtered');

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function sgmga_product_weight_tests ( )

%****************************************************************************80
%
%% SGMGA_PRODUCT_WEIGHT_TESTS calls SGMGA_PRODUCT_WEIGHT_TEST.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2010
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../sandia_rules' );

  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SGMGA_PRODUCT_WEIGHT_TESTS\n' );
  fprintf ( 1, '  Call SGMGA_PRODUCT_WEIGHT_TEST with various arguments.\n' );

  dim_num = 2;
  order_1d = [ 3, 5 ]';
  order_nd = prod ( order_1d(1:dim_num) );
  rule = [ 1, 1 ]';
  np = [ 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  sgmga_product_weight_test ( dim_num, order_1d, order_nd, rule, np, p );

  dim_num = 2;
  order_1d = [ 3, 7 ]';
  order_nd = prod ( order_1d(1:dim_num) );
  rule = [ 1, 5 ]';
  np = [ 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  sgmga_product_weight_test ( dim_num, order_1d, order_nd, rule, np, p );

  dim_num = 2;
  order_1d = [ 3, 3 ]';
  order_nd = prod ( order_1d(1:dim_num) );
  rule = [ 3, 7 ]';
  np = [ 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  sgmga_product_weight_test ( dim_num, order_1d, order_nd, rule, np, p );

  dim_num = 2;
  order_1d = [ 5, 5 ]';
  order_nd = prod ( order_1d(1:dim_num) );
  rule = [ 1, 8 ]';
  np = [ 0, 1 ]';
  np_sum = sum ( np(1:dim_num) );
  p(1,1) = 1.5;
  sgmga_product_weight_test ( dim_num, order_1d, order_nd, rule, np, p );

  dim_num = 2;
  order_1d = [ 5, 5 ]';
  order_nd = prod ( order_1d(1:dim_num) );
  rule = [ 2, 9 ]';
  np = [ 0, 2 ]';
  np_sum = sum ( np(1:dim_num) );
  p(1,1) = 0.5;
  p(2,1) = 1.5;
  sgmga_product_weight_test ( dim_num, order_1d, order_nd, rule, np, p );

  dim_num = 2;
  order_1d = [ 7, 9 ]';
  order_nd = prod ( order_1d(1:dim_num) );
  rule = [ 6, 10 ]';
  np = [ 1, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p(1,1) = 2.0;
  sgmga_product_weight_test ( dim_num, order_1d, order_nd, rule, np, p );

  dim_num = 3;
  order_1d = [ 2, 3, 3 ]';
  order_nd = prod ( order_1d(1:dim_num) );
  rule = [ 1, 4, 5 ]';
  np = [ 0, 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  sgmga_product_weight_test ( dim_num, order_1d, order_nd, rule, np, p );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SGMGA_PRODUCT_WEIGHT_TESTS:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../sandia_rules' );

  return
end


text:
function [V,Q,QT,QF,QL] = voxel_surface(W,varargin)
  % VOXEL_SURFACE Compute the surface quad mesh of a 3d logical image.
  %
  % [V,Q,QT,QF,QL] = voxel_surface(W,varargin)
  %
  % Inputs:
  %   W  n by m by k binary matix
  %   Optional:
  %     'Centers' followed by n*m*k by 3 list of cell centers
  % Outputs:
  %   V  (n+1)*(m+1)*(k+1) by 3 list of cell corner positions
  %   Q  #Q by 4 list of quads indexing V
  %   QT  n*m*k by 4 list of "top" quad faces indexing V
  %   QF  n*m*k by 4 list of "front" quad faces indexing V
  %   QL  n*m*k by 4 list of "left" quad faces indexing V
  %
  % Examples:
  %   [V,Q] = voxel_surface(W);
  %   [V,IM] = remove_unreferenced(V,Q);
  %   Q = IM(Q);
  %   trisurf(Q,V(:,1),V(:,2),V(:,3));
  %   axis equal;
  %
  % See also: voxelize, voxel_grid

  BC = [];
  % default values
  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'Centers'}, ...
    {'BC'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  dim = ndims(W);

  switch dim
  case 2
    side = [size(W,2) size(W,1)];
    if isempty(BC)
      [X,Y] = meshgrid( ...
        linspace(0,1,side(1)), ...
        linspace(0,1,side(2)));
      % barycenters of cells
      BC = [X(:) Y(:)];
    end
    NV = min(BC);
    XV = max(BC);
    r = (XV-NV)./(side-1);
    [X,Y] = meshgrid( ...
      (NV(1)-0.5*r(1))+linspace(0,1,side(1)+1)*(XV(1)-NV(1)+r(1)), ...
      (NV(2)-0.5*r(2))+linspace(0,1,side(2)+1)*(XV(2)-NV(2)+r(2)));

    [Q,V] = surf2patch(X,Y,0*X);
    V = V(:,1:2);
    QT = [Q(:,[3 4]);Q(:,[1 2]);];
    QF = [Q(:,[2 3]); Q(:,[4 1])];
    if ~any(W(:))
      Q = [];
    else
      assert(false,'non trivial W not supported');
    end

  case 3
    side = [size(W,2) size(W,1) size(W,3)];
    if isempty(BC)
      [X,Y,Z] = meshgrid( ...
        linspace(0,1,side(1)), ...
        linspace(0,1,side(2)), ...
        linspace(0,1,side(3)));
      % barycenters of cells
      BC = [X(:) Y(:) Z(:)];
    end

    NV = min(BC);
    XV = max(BC);
    r = (XV-NV)./(side-1);
    [X,Y,Z] = meshgrid( ...
      (NV(1)-0.5*r(1))+linspace(0,1,side(1)+1)*(XV(1)-NV(1)+r(1)), ...
      (NV(2)-0.5*r(2))+linspace(0,1,side(2)+1)*(XV(2)-NV(2)+r(2)), ...
      (NV(3)-0.5*r(3))+linspace(0,1,side(3)+1)*(XV(3)-NV(3)+r(3)));
    % corners of cells
    V = [X(:) Y(:) Z(:)];

    [II,JJ,KK] = ind2sub([side(2) side(1) side(3)]+1,reshape(1:size(V,1),[side(2) side(1) side(3)]+1));

    QF = [ ...
        reshape(sub2ind([side(2) side(1) side(3)]+1,II(1:end-1,1:end-1,1:end),JJ(1:end-1,1:end-1,1:end),KK(1:end-1,1:end-1,1:end)),[],1) ...
        reshape(sub2ind([side(2) side(1) side(3)]+1,II(  2:end,1:end-1,1:end),JJ(  2:end,1:end-1,1:end),KK(  2:end,1:end-1,1:end)),[],1) ...
        reshape(sub2ind([side(2) side(1) side(3)]+1,II(  2:end,  2:end,1:end),JJ(  2:end,  2:end,1:end),KK(  2:end,  2:end,1:end)),[],1) ...
        reshape(sub2ind([side(2) side(1) side(3)]+1,II(1:end-1,  2:end,1:end),JJ(1:end-1,  2:end,1:end),KK(1:end-1,  2:end,1:end)),[],1) ...
        ];

    QL = fliplr([ ...
        reshape(sub2ind([side(2) side(1) side(3)]+1,II(1:end-1,1:end,1:end-1),JJ(1:end-1,1:end,1:end-1),KK(1:end-1,1:end,1:end-1)),[],1) ...
        reshape(sub2ind([side(2) side(1) side(3)]+1,II(  2:end,1:end,1:end-1),JJ(  2:end,1:end,1:end-1),KK(  2:end,1:end,1:end-1)),[],1) ...
        reshape(sub2ind([side(2) side(1) side(3)]+1,II(  2:end,1:end,  2:end),JJ(  2:end,1:end,  2:end),KK(  2:end,1:end,  2:end)),[],1) ...
        reshape(sub2ind([side(2) side(1) side(3)]+1,II(1:end-1,1:end,  2:end),JJ(1:end-1,1:end,  2:end),KK(1:end-1,1:end,  2:end)),[],1) ...
        ]);

    QT = [ ...
        reshape(sub2ind([side(2) side(1) side(3)]+1,II(1:end,1:end-1,1:end-1),JJ(1:end,1:end-1,1:end-1),KK(1:end,1:end-1,1:end-1)),[],1) ...
        reshape(sub2ind([side(2) side(1) side(3)]+1,II(1:end,  2:end,1:end-1),JJ(1:end,  2:end,1:end-1),KK(1:end,  2:end,1:end-1)),[],1) ...
        reshape(sub2ind([side(2) side(1) side(3)]+1,II(1:end,  2:end,  2:end),JJ(1:end,  2:end,  2:end),KK(1:end,  2:end,  2:end)),[],1) ...
        reshape(sub2ind([side(2) side(1) side(3)]+1,II(1:end,1:end-1,  2:end),JJ(1:end,1:end-1,  2:end),KK(1:end,1:end-1,  2:end)),[],1) ...
        ];

    if ~any(W(:))
      Q = [];
    else
      Wp = padarray(W,[1 1 1],0);
      Dy = diff(Wp,1,1);
      Dx = diff(Wp,1,2);
      Dz = diff(Wp,1,3);
      Q = [ ...
        QF(        Dz(2:end-1,2:end-1,1:end)>0.5,:); ...
        QL(        Dx(2:end-1,1:end,2:end-1)>0.5,:); ...
        QT(        Dy(1:end,2:end-1,2:end-1)>0.5,:); ...
        fliplr(QF( Dz(2:end-1,2:end-1,1:end)<-0.5,:)); ...
        fliplr(QL( Dx(2:end-1,1:end,2:end-1)<-0.5,:)); ...
        fliplr(QT( Dy(1:end,2:end-1,2:end-1)<-0.5,:)); ...
        ];
    end
  end
end

author: alecjacobson

repo: gptoolbox


text:
% mean functions to be use by Gaussian process functions. There are two
% different kinds of mean functions: simple and composite:
%
% simple mean functions:
%
%   meanZero      - zero mean function
%   meanOne       - one mean function
%   meanConst     - constant mean function
%   meanLinear    - linear mean function
%   meanPoly      - polynomial mean function
%   meanDiscrete  - precomputed mean for discrete data
%   meanGP        - predictive mean of another GP
%   meanGPexact   - predictive mean of a regression GP
%   meanNN        - nearest neighbor mean function
%
% composite covariance functions (see explanation at the bottom):
%
%   meanScale     - scaled version of a mean function
%   meanPow       - power of a mean function
%   meanProd      - products of mean functions
%   meanSum       - sums of mean functions
%   meanMask      - mask some dimensions of the data
%   meanPref      - difference mean for preference learning
%
% Naming convention: all mean functions are named "mean/mean*.m".
%
%
% 1) With no or only a single input argument:
%
%    s = meanNAME  or  s = meanNAME(hyp)
%
% The mean function returns a string s telling how many hyperparameters hyp it
% expects, using the convention that "D" is the dimension of the input space.
% For example, calling "meanLinear" returns the string 'D'.
%
% 2) With two input arguments:
%
%    m = meanNAME(hyp, x)
%
% The function computes and returns the mean vector where hyp are the
% hyperparameters and x is an n by D matrix of cases, where D is the dimension
% of the input space. The returned mean vector is of size n by 1.
%
% 3) With three input arguments:
%
%    dm = meanNAME(hyp, x, i)
%
% The function computes and returns the n by 1 vector of partial derivatives
% of the mean vector w.r.t. hyp(i) i.e. hyperparameter number i.
%
% See also doc/usageMean.m.
%
% Copyright (c) by Carl Edward Rasmussen and Hannes Nickisch, 2014-12-08.
%                                      File automatically generated using noweb.

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


function [] = create_gnuradio_zc(fft_size, symbol_number)
    zc = conj(create_zc(fft_size, symbol_number));

    fprintf('[');
    for idx=1:length(zc)
        fprintf('(%0.4f+%0.4fj)', real(zc(idx)), imag(zc(idx)));
        if (idx ~= length(zc))
            fprintf(',')
        end
    end
    fprintf(']\n');
end


text:
function option = femoption(option)
%% FEMOPTION provides default options for femPoisson
%
% option = femoption(option)
%
%  - options.elemType:  type of finite element
%  - options.maxIt
%  - options.maxN
%  - options.L0
%  - options.refType    mesh refinement type
%  - options.printlevel
%  - options.plotflag
%  - options.rateflag
%  - options.dispflag
%  - options.tol
%  - options.lumpflag

if ~isfield(option,'elemType')
    option.elemType = 'P1';
end

if ~isfield(option,'maxIt')
    option.maxIt = 4;
end

if ~isfield(option,'maxN')
    option.maxN = 2e5;
end

if ~isfield(option,'L0')
    option.L0 = 0;
end

if ~isfield(option,'refType')
    option.refType = 'red';
end

if ~isfield(option,'printlevel')
    option.printlevel = 1;
end

if ~isfield(option,'plotflag')
    option.plotflag = 1;
end

if ~isfield(option,'rateflag')
    option.rateflag = 1;
end

if ~isfield(option,'dispflag')
    option.dispflag = 1;
end

if ~isfield(option,'tol')
    option.tol = 1e-8;
end

if ~isfield(option,'lumpflag') % mass lumping
    option.lumpflag = 0;       % no mass lumping
end

repo: ifem


text:
function varargout = projPointOnPolygon(point, poly, varargin)
%PROJPOINTONPOLYGON  Compute position of a point projected on a polygon
%
%   POS = projPointOnPolygon(POINT, POLYGON)
%   Compute the position of the orthogonal projection of a point on a
%   polygon.
%   POINT is a 1-by-2 row vector containing point coordinates
%   POLYGON is a N-by-2 array containing coordinates of polygon vertices
%
%   When POINT is an array of points, returns a column vector with as many
%   rows as the number of points.
%
%   [POS, DIST] = projPointOnPolygon(...)
%   Also returns the distance between POINT and POLYGON. The distance is
%   negative if the point is located inside of the polygon.
%
%   Example
%     poly = [10 10; 20 10;20 20;10 20];
%     projPointOnPolygon([15 0], poly)
%     ans =
%         0.5000
%     projPointOnPolygon([0 16], poly)
%     ans =
%         3.4000
%
%   See also
%   points2d, polygons2d, polygonPoint, projPointOnPolyline
%   distancePointpolygon
%

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2009-04-30,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.

% eventually copy first point at the end to ensure closed polygon
if sum(poly(end, :) == poly(1,:)) ~= 2
    poly = [poly; poly(1,:)];
end

% compute position wrt outline
[pos, minDist] = projPointOnPolyline(point, poly);

% process output arguments
if nargout <= 1
    varargout{1} = pos;
elseif nargout == 2
    varargout{1} = pos;
    if inpolygon(point(:,1), point(:,2), poly(:,1), poly(:,2))
        minDist = -minDist;
    end
    varargout{2} = minDist;
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [params, names] = rbfwhiteKernExtractParam(kern)

% RBFWHITEKERNEXTRACTPARAM Extract parameters from the RBF-WHITE kernel
% structure.
% FORMAT
% DESC extracts parameters from the RBF-WHITE kernel structure into a vector
% of parameters for optimisation.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the kernel. If
% the field 'transforms' is not empty in the kernel structure, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
%
% FORMAT
% DESC extracts parameters and parameter names from the RBF-WHITE kernel
% structure.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the kernel. If the
% field 'transforms' is not empty in the kernel structure, the parameters
% will be transformed before optimisation (for example positive only
% parameters could be logged before being returned).
% RETURN names : cell array of strings containing names for each
% parameter.
%
% SEEALSO rbfwhiteKernParamInit, rbfwhiteKernExpandParam, kernExtractParam,
% scg, conjgrad
%
% COPYRIGHT : David Luengo, 2009
%
% KERN

params = [kern.inverseWidth kern.variance];
if nargout > 1
  names = {'inverse width', 'variance'};
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function bernoulli_poly_values_test ( )

%*****************************************************************************80
%
%% BERNOULLI_POLY_VALUES_TEST demonstrates the use of BERNOULLI_POLY_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BERNOULLI_POLY_VALUES_TEST\n' );
  fprintf ( 1, '  BERNOULLI_POLY_VALUES stores values of\n' );
  fprintf ( 1, '  the Bernoulli polynomials.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      N            X            FX\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, x, fx ] = bernoulli_poly_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %6d  %12f  %24.16f\n', n, x, fx );

  end

  return
end


text:
[x,y]=meshgrid([-3:0.2:3]);
z=(sin(x.*y)+eps)./(x.*y+eps);
surf(x,y,z)

author: Eurus-Holmes


text:
%% Pose Optimization Demo - Weak Perspective case
% Weak Perspective case of our Pose Optimization.
% For convenience, we assume that the heatmaps (keypoint localizations)
% are precomputed and provided in the folder demo.

clear
startup

datapath = 'demo/pascal3d-sample/';
annotfile = sprintf('%s/annot/valid.mat',datapath);
load(annotfile);

% determines shape model used for pose optimization
% if set to 1 then we use the particular cad model instance
% if set to 0 then we use the deformable shape model
cadSpecific = 1;

for ID = 1:length(annot.imgname)

    % input
    imgname = annot.imgname{ID};
    center = annot.center(ID,:);
    scale = annot.scale(ID);
    class = annot.class{ID};
    indices = annot.indices{ID};
    cadID = annot.cad_index(ID);

    cad = load(sprintf('cad/%s.mat',class));
    cad = cad.(class);
    cad = cad(cadID);

    if cadSpecific
        dict = getPascalTemplate(cad);
    else
        dict = load(sprintf('dict/pca-%s.mat',class));
    end
    % read heatmaps and detect maximum responses
    heatmap = h5read(sprintf('%s/exp/valid_%d.h5',datapath,ID),'/heatmaps');
    heatmap = permute(heatmap(:,:,indices(dict.kpt_id)),[2,1,3]);
    [W_hp,score] = findWmax(heatmap);

    % pose optimization - weak perspective
    output_wp = PoseFromKpts_WP(W_hp,dict,'weight',score,'verb',false,'lam',1);

    % visualization
    img = imread(sprintf('%s/images/%s.jpg',datapath,imgname));
    vis_wp(img,output_wp,heatmap,center,scale,cad,dict);
    pause
    close all

end

author: geopavlakos


text:
function ierror = part_sf_check ( n, npart, a )

%*****************************************************************************80
%
%% PART_SF_CHECK checks a standard form partition of an integer.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher, Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998,
%    ISBN: 0-8493-3988-X,
%    LC: QA164.K73.
%
%  Parameters:
%
%    Input, integer N, the integer to be partitioned.
%    N must be positive.
%
%    Input, integer NPART, the number of parts of the partition.
%    1 <= NPART <= N.
%
%    Input, integer A(NPART), contains the partition.
%    A(1) through A(NPART) contain the nonzero integers which
%    sum to N.  The entries must be in DESCENDING order.
%
%    Output, integer IERROR, error flag.
%    0, no error.
%    -1, N is illegal.
%    -2, NPART is illegal.
%    -3, the entries do not add up to N.
%    I, the I-th entry of A is illegal.
%
  ierror = 0;

  if ( n < 1 )
    ierror = -1;
    return
  end

  if ( npart < 1 || n < npart )
    ierror = -2;
    return
  end
%
%  Every entry must lie between 1 and N.
%
  for i = 1 : npart
    if ( a(i) < 1 || n < a(i) )
      ierror = i;
      return
    end
  end
%
%  The entries must be in descending order.
%
  for i = 2 : npart
    if ( a(i-1) < a(i) )
      ierror = i;
      return
    end
  end
%
%  The entries must add up to N.
%
  if ( sum ( a(1:npart) ) ~= n )
    ierror = -3;
  end

  return
end


text:
classdef BiGE < ALGORITHM
% <many> <real/integer/label/binary/permutation>
% Bi-goal evolution

%------------------------------- Reference --------------------------------
% M. Li, S. Yang, and X. Liu, Bi-goal evolution for many-objective
% optimization problems, Artificial Intelligence, 2015, 228: 45-65.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Generate random population
            Population = Problem.Initialization();

            %% Optimization
            while Algorithm.NotTerminated(Population)
                MatingPool = MatingSelection(Estimation(Population.objs,1/Problem.N^(1/Problem.M)));
                Offspring  = OperatorGA(Problem,Population(MatingPool));
                Population = EnvironmentalSelection([Population,Offspring],Problem.N);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function loglik = compute_log_lik(P, S_bar, V, E_z, E_zz, RO, Tr, sigma_sq)

[K, T] = size(E_z);
J = size(S_bar, 2);

M_t = zeros(2*J, K);

loglik = - 0.5*T * (2*J*log(sigma_sq));
for t = 1:T,
   R_t = RO{t};

   Sdef = S_bar;
   for kk = 1:K,
      Sdef = Sdef + E_z(kk,t)*V((kk-1)*3+[1:3],:);

      M_t(1:J, kk) = (R_t(1,:)*V((kk-1)*3+[1:3],:))';
      M_t(J+1:end, kk) = (R_t(2,:)*V((kk-1)*3+[1:3], :))';
   end;

   invSigmaSq_p = eye(2*J)./sigma_sq;

   f_bar_t = R_t(1:2,:)*S_bar;
   f_bar_t = [f_bar_t(1,:) f_bar_t(2,:)]';

   f_t = [P(t, :) P(t+T, :)]';
   t_vect_t = [Tr(t,1)*ones(J,1); Tr(t,2)*ones(J,1)];

   covZ_t = E_zz((t-1)*K+1:t*K,:) - E_z(:,t)*E_z(:,t)';
   loglik = loglik - 0.5*(((f_t-f_bar_t-t_vect_t)./sigma_sq)'*(f_t-f_bar_t-t_vect_t)) + (((f_t-f_bar_t-t_vect_t)'./sigma_sq)*M_t*E_z(:,t)) ...
      - 0.5*trace(((M_t./sigma_sq)'*M_t) * E_zz((t-1)*K+1:t*K,:)) - 0.5*log(det(covZ_t));
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function x = blend_ij_w_1d1 ( x, r, s, m1, m2 )

%*****************************************************************************80
%
%% BLEND_IJ_W_1D1 extends weighted indexed scalar data along edges into a table.
%
%  Diagram:
%
%    Instead of assuming that the data in the table is equally spaced,
%    the arrays R and S are supplied, which should behave as
%    "coordinates" for the data.
%
%            S(1)  S(2)  S(3)  S(4)  S(5)  S(6)  S(M2)
%
%    R(1)  ( X11,  X12,  X13,  X14,  X15,  X16,  X1M2  )
%    R(2)  ( X21,  ...,  ...,  ...,  ...,  ...,  X2M2  )
%    R(3)  ( X31,  ...,  ...,  ...,  ...,  ...,  X3M2  )
%    R(4)  ( X41,  ...,  ...,  ...,  ...,  ...,  X4M2  )
%    R(M1) ( XM11, XM12, XM13, XM14, XM15, XM16, XM1M2 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    William Gordon,
%    Blending-Function Methods of Bivariate and Multivariate Interpolation
%    and Approximation,
%    SIAM Journal on Numerical Analysis,
%    Volume 8, Number 1, March 1971, pages 158-177.
%
%    William Gordon and Charles Hall,
%    Transfinite Element Methods: Blending-Function Interpolation over
%    Arbitrary Curved Element Domains,
%    Numerische Mathematik,
%    Volume 21, Number 1, 1973, pages 109-129.
%
%    William Gordon and Charles Hall,
%    Construction of Curvilinear Coordinate Systems and Application to
%    Mesh Generation,
%    International Journal of Numerical Methods in Engineering,
%    Volume 7, 1973, pages 461-477.
%
%    Joe Thompson, Bharat Soni, Nigel Weatherill,
%    Handbook of Grid Generation,
%    CRC Press, 1999.
%
%  Parameters:
%
%    Input, real X(M1,M2).
%    On input, data is contained in the "edge entries" X(1,J), X(I,1),
%    X(M1,J) and X(I,M2), for I = 1 to M1, and J = 1 to M2.
%
%    Input, real R(M1), S(M2), are "coordinates" for the rows and
%    columns of the array.  The values in R, and the values in S, should
%    be strictly increasing or decreasing.
%
%    Input, integer M1, M2, the number of rows and columns in X.
%
%    Output, real X(M1,M2), all entries in X have been assigned a value.
%

%
%  Interpolate values in the interior.
%
  for i = 2 : m1 - 1

    rr = ( r(i) - r(1) ) / ( r(m1) - r(1) );

    for j = 2 : m2 - 1

      ss = ( s(j) - s(1) ) / ( s(m2) - s(1) );

      x(i,j) = blend_112 ( rr, ss, x(1,1), x(1,m2), x(m1,1), x(m1,m2), ...
        x(i,1), x(i,m2), x(1,j), x(m1,j) );

    end

  end

  return
end


text:
function runs = contiguous(A,varargin)
%   RUNS = CONTIGUOUS(A,NUM) returns the start and stop indices for contiguous
%   runs of the elements NUM within vector A.  A and NUM can be vectors of
%   integers or characters.  Output RUNS is a 2-column cell array where the ith
%   row of the first column contains the ith value from vector NUM and the ith
%   row of the second column contains a matrix of start and stop indices for runs
%   of the ith value from vector NUM.    These matrices have the following form:
%
%   [startRun1  stopRun1]
%   [startRun2  stopRun2]
%   [   ...        ...  ]
%   [startRunN  stopRunN]
%
%   Example:  Find the runs of '0' and '2' in vector A, where
%             A = [0 0 0 1 1 2 2 2 0 2 2 1 0 0];
%
%   runs = contiguous(A,[0 2])
%   runs =
%           [0]    [3x2 double]
%           [2]    [2x2 double]
%
%   The start/stop indices for the runs of '0' are given by runs{1,2}:
%
%           1     3
%           9     9
%          13    14
%
%   RUNS = CONTIGUOUS(A) with only one input returns the start and stop
%   indices for runs of all unique elements contained in A.
%
%   CONTIGUOUS is intended for use with vectors of integers or characters, and
%   is probably not appropriate for floating point values.  You decide.
%

if prod(size(A)) ~= length(A),
    error('A must be a vector.')
end

if isempty(varargin),
    num = unique(A);
else
    num = varargin{1};
    if prod(size(num)) ~= length(num),
        error('NUM must be a scalar or vector.')
    end
end

for numCount = 1:length(num),

    indexVect = find(A(:) == num(numCount));
    shiftVect = [indexVect(2:end);indexVect(end)];
    diffVect = shiftVect - indexVect;

    % The location of a non-one is the last element of the run:
    transitions = (find(diffVect ~= 1));

    runEnd = indexVect(transitions);
    runStart = [indexVect(1);indexVect(transitions(1:end-1)+1)];

    runs{numCount,1} = num(numCount);
    runs{numCount,2} = [runStart runEnd];

end


text:
function x = cs_dmsol (A,b)
%CS_DMSOL x=A\b using the coarse Dulmage-Mendelsohn decomposition.
%   x = cs_dmsol(A,b) computes x=A\b where A may be rectangular and/or
%   structurally rank deficient, and b is a full vector.
%
%   Example:
%       Prob = UFget ('HB/arc130') ; A = Prob.A ; b = rand (size (A,1),1) ;
%       x = cs_dmsol (A,b) ; norm (A*x-b)
%
%   See also CS_QRSOL, CS_LUSOL, CS_DMPERM, SPRANK, RANK.

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

[m n] = size (A) ;
[p q r s cc rr] = cs_dmperm (A) ;
C = A (p,q) ;
b = b (p) ;
x = zeros (n,1) ;
if (rr(3) <= m & cc(4) <= n)                                                %#ok
    x (cc(4):n) = cs_qrsol (C (rr(3):m, cc(4):n), b (rr(3):m)) ;
    b (1:rr(3)-1) = b (1:rr(3)-1) - C (1:rr(3)-1, cc(4):n) * x (cc(4):n) ;
end
if (rr(2) < rr (3) & cc(3) < cc(4))                                         %#ok
    x (cc(3):cc(4)-1) = ...
        cs_lusol (C (rr(2):rr(3)-1, cc(3):cc(4)-1), b (rr(2):rr(3)-1)) ;
    b (1:rr(2)-1) = ...
        b (1:rr(2)-1) - C (1:rr(2)-1, cc(3):cc(4)-1) * x (cc(3):cc(4)-1) ;
end
if (rr(2) > 1 & cc(3) > 1)                                                  %#ok
    x (1:cc(3)-1) = cs_qrsol (C (1:rr(2)-1, 1:cc(3)-1), b (1:rr(2)-1)) ;
end
x (q) = x ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% PDETIMEDEMO
% This program creates the left preconditionied time-dependent PDE example
% from Chapter 3
%
%
% The example has the form
%
% u_t = - L(u) + f, 0 < x,y < 1
% with homogeneous Dirichlet BC
%
% In terms of the program below, the equation is
% Lu = - (u_xx + u_yy) + 20*u*(u_x + u_y)
%
% Within the nonlinear map, pdetime.m, we apply fish2d.m as
% a preconditioner.
%
% The action of the Laplacian, partial wrt x, and partial wrt y
% are computed in the routines lapmf, dxmf, and dymf. Here mf
% stands for MATRIX FREE.
%
global rhsf uold dt;
m=63; n=m*m;
h=1/(m+1);
h2=(m+1)*(m+1);
h1=(m+1);
tolh=1/h2;
mres=3;
%
% set up the equation
%
x=1:m;
e=x';
x=x'*h;
%
% Form the steady-state solution: sol(x,y)= 10 x y (1-x) (1-y) exp(x^4.5)
%
zsoly=x.*(1-x);
zsolx=zsoly.*exp(x.^4.5);
solt=10*zsolx*zsoly';
sol=solt(:);
clear rhs; clear axt; clear at; clear ayt; clear solt;
%
% Fix the right side so that the steady-state solution is known
%
rhsf = lapmf(sol) + 20*sol.*(dxmf(sol)+dymf(sol));
%
% Set the initial data to zero.
%
u0=zeros(m*m,1);
uold=u0;
uinit=uold;
dt=.1;
nt=10;
%
%
% Use nsoli.m with GMRES as the solver.
% Converged result from previous time step is initial iterate.
%
for it=1:nt+1
    norm(uold-sol)
    tol=[1.d-1,1.d-1]*tolh;
    x=zeros(n,1); parms = [40,40,-.1,1];
    [unew, it_histg, ierr] = nsoli(uinit,'pdetime',tol,parms);
    uinit=uold;
    uold=unew;
end
norm(uold-sol)

repo: Aero-matlab


text:
%MDL_JACO Create model of Kinova Jaco manipulator
%
% MDL_JACO is a script that creates the workspace variable jaco which
% describes the kinematic characteristics of a Kinova Jaco manipulator
% using standard DH conventions.
%
% Also define the workspace vectors:
%   qz         zero joint angle configuration
%   qr         vertical 'READY' configuration
%
% Reference::
% - "DH Parameters of Jaco" Version 1.0.8, July 25, 2013.
%
% Notes::
% - SI units of metres are used.
% - Unlike most other mdl_xxx scripts this one is actually a function that
%   behaves like a script and writes to the global workspace.
%
% See also SerialLink, mdl_mico, mdl_puma560.

% MODEL: Kinova, Jaco, 6DOF, standard_DH

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

function r = mdl_jaco()

    deg = pi/180;

    % robot length values (metres)
    D1 = 0.2755;
    D2 = 0.4100;
    D3 = 0.2073;
    D4 = 0.0743;
    D5 = 0.0743;
    D6 = 0.1687;
    e2 = 0.0098;

    % alternate parameters
    aa = 30*deg;
    ca = cos(aa);
    sa = sin(aa);
    c2a = cos(2*aa);
    s2a = sin(2*aa);
    d4b = D3 + sa/s2a*D4;
    d5b = sa/s2a*D4 + sa/s2a*D5;
    d6b = sa/s2a*D5 + D6;

    % and build a serial link manipulator

    % offsets from the table on page 4, "Mico" angles are the passed joint
    % angles.  "DH Algo" are the result after adding the joint angle offset.

    robot = SerialLink([
        Revolute('alpha', pi/2,  'a', 0,  'd', D1,   'flip')
        Revolute('alpha', pi,    'a', D2, 'd', 0,    'offset', -pi/2)
        Revolute('alpha', pi/2,  'a', 0,  'd', -e2,  'offset', pi/2)
        Revolute('alpha', 2*aa,  'a', 0,  'd', -d4b)
        Revolute('alpha', 2*aa,  'a', 0,  'd', -d5b, 'offset', -pi)
        Revolute('alpha', pi,    'a', 0,  'd', -d6b, 'offset', 100*deg)
        ], ...
        'name', 'Jaco', 'manufacturer', 'Kinova');

    % place the variables into the global workspace
    if nargin == 1
        r = robot;
    elseif nargin == 0
        assignin('base', 'jaco', robot);
        assignin('base', 'qz', [0 0 0 0 0 0]); % zero angles
        assignin('base', 'qr', [270 180 180 0 0 0]*deg); % vertical pose as per Fig 2
    end
end

author: Allopart

repo: rbpf-gmapping


text:
function quad = fibonacci_lattice_q2 ( k, f )

%*****************************************************************************80
%
%% FIBONACCI_LATTICE_Q2 applies a Fibonacci lattice integration rule in 2D.
%
%  Discussion:
%
%    This is a modification of the algorithm in FIBONACCI_LATTICE_Q.
%    It uses a nonlinear transformation on the integrand, which makes
%    the lattice rule more suitable for nonperiodic integrands.
%
%    The transformation replaces the integration variable X by
%
%      PHI(X) = 3*X^3 *( 10 - 15 * X + 6 * X^2 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Ian Sloan, Stephen Joe,
%    Lattice Methods for Multiple Integration,
%    Oxford, 1994,
%    ISBN: 0198534728,
%    LC: QA311.S56
%
%  Parameters:
%
%    Input, integer K, the index of the Fibonacci number to be used.
%    K must be at least 3.
%
%    Input, external real F, the name of the user-supplied routine
%    which evaluates the function, of the form:
%    function f ( dim_num, x )
%    integer dim_num
%    real f
%    real x(dim_num)
%    f = ...
%
%    Output, real QUAD, the estimated integral.
%
  dim_num = 2;

  quad = 0.0;

  m = fibonacci ( k );

  z(1) = 1;
  z(2) = fibonacci ( k - 1 );

  for j = 0 : m - 1
    x(1:dim_num) = mod ( j * z(1:dim_num) / m, 1.0 );
    dphi = 1.0;
    for i = 1 : dim_num
      y(i) = ( 10.0 - 15.0 * x(i) + 6.0 * x(i)^2 ) * x(i)^3;
      dphi = dphi * 30.0 * ( 1.0 - x(i) )^2 * x(i)^2;
    end
    quad = quad + f ( dim_num, y ) * dphi;
  end

  quad = quad / m;

  return
end


text:
function [ readerobj ] = open_iceye_reader( filename )
%OPEN_ICEYE_READER Intiates a reader object for ICEYE SLC HDF5 file format.
%
% Written by: Tim Cox, NRL.  Based on the CSM/NISAR reader by Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%% Open HDF5
fid=H5F.open(filename, 'H5F_ACC_RDONLY', 'H5P_DEFAULT');
meta=meta2sicd_iceye(fid);

% Open HDF objects for reading real and imaginary pixels
dset_id(1) = H5D.open(fid,'s_i');
dspace_id(1) = H5D.get_space(dset_id(1));
[~, datasize{1}] = H5S.get_simple_extent_dims(dspace_id(1));
dset_id(2) = H5D.open(fid,'s_q');
dspace_id(2) = H5D.get_space(dset_id(2));
[~, datasize{2}] = H5S.get_simple_extent_dims(dspace_id(2));

if strcmpi(h5read(filename,'/look_side'),'left')
    symmetry=[0 1 1];
else
    symmetry=[0 0 1];
end

%% Define object methods
readerobj = chipfun2readerobj(@chipper_iceye, datasize{1}); %real/imag are the same size
readerobj.get_meta=@() meta;
readerobj.close=@() close_hdf;

    function out = chipper_iceye(varargin)
        [dim1range,dim2range,subsample]=check_chipper_args(datasize{1}, varargin{:});
        tempargs=reorient_chipper_args(symmetry, datasize{1}([2 1]), dim1range, dim2range, subsample);
        [dim1range,dim2range,subsample]=deal(tempargs{:});
        offset=[dim2range(1) dim1range(1)]-1;
        subsample([1 2]) = subsample([2 1]);
        slabsize=ceil(([diff(dim2range) diff(dim1range)]+1)./subsample);
        H5S.select_hyperslab(dspace_id(1),'H5S_SELECT_SET',...
            offset,subsample,slabsize,ones(size(offset)));
        memspacei_id = H5S.create_simple(length(slabsize),slabsize,slabsize);
        temp_out_real=H5D.read(dset_id(1),'H5ML_DEFAULT',memspacei_id,dspace_id(1),'H5P_DEFAULT');
        H5S.select_hyperslab(dspace_id(2),'H5S_SELECT_SET',...
            offset,subsample,slabsize,ones(size(offset)));
        memspaceq_id = H5S.create_simple(length(slabsize),slabsize,slabsize);
        temp_out_imag=H5D.read(dset_id(2),'H5ML_DEFAULT',memspaceq_id,dspace_id(2),'H5P_DEFAULT');

        H5S.close(memspacei_id);
        H5S.close(memspaceq_id);
        out=reorient_chipper_data(symmetry, complex(temp_out_real,temp_out_imag));
    end

    function close_hdf()
        H5S.close(dspace_id(1));
        H5D.close(dset_id(1));
        H5S.close(dspace_id(2));
        H5D.close(dset_id(2));
        H5F.close(fid);
    end

end

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


text:
%% SIMPLE FIXED WING AGENT (fixedWing.m) //////////////////////////////////
% This class is designed to contain the placeholder dynamics of a
% fixed-wing aircraft.

% Author: James A. Douthwaite

classdef fixedWing < agent
    properties

    end
    %% ///////////////////////// MAIN METHODS /////////////////////////////
    methods
        % Constructor
        function [this] = fixedWing(varargin)
            % Call the super class
            this@agent(varargin);                                           % Create the super class 'agent'
            % //////////////// Check for user overrides ///////////////////
            this = this.ApplyUserOverrides(varargin); % Recursive overrides
            % /////////////////////////////////////////////////////////////
        end
        % Setup - 3D DUBLINS CAR STATE [x;y;z;v;psi;theta]
        function [this] = setup(this,localXYZVelocity,localXYZrotations)
            % Infer the initial conditions from the scenario.
            p     = zeros(3,1);            % The position in the local frame
            speed = norm(localXYZVelocity);% The speed of the aircraft
            theta = localXYZrotations(2);  % The elevation of the aircraft
            psi   = 0;                     % The yaw of the aircraft
            % Define the initial state from the scenario
            this.localState = [p(1);p(2);p(3);speed;theta;psi];
            this = this.SetGLOBAL('priorState',this.localState);
        end
        % Main
        function [this] = main(this,ENV,varargin)
            % INPUTS:
            % varargin - Cell array of inputs
            % >dt      - The timestep
            % >objects - The detectable objects cell array of structures
            % OUTPUTS:
            % obj      - The updated project

            % //////////// CHECK FOR NEW INFORMATION UPDATE ///////////////
            % UPDATE THE AGENT WITH THE NEW ENVIRONMENTAL INFORMATION
            %[obj,obstacleSet,agentSet,waypointSet] = obj.getAgentUpdate(varargin{1});                   % IDEAL INFORMATION UPDATE

            u_k = [0.1;0.1;0.1]; % Acceleration, pitch rate, yaw rate
            x_k = this.localState;

            % UPDATE LOCAL STATE
            newState = this.UpdateLocalState(ENV,x_k,u_k);

            % UPDATE THE GLOBAL PROPERTIES OF THE AGENT
            [this] = this.GlobalUpdate(ENV.dt,newState);

            % /////////////// RECORD THE AGENT-SIDE DATA //////////////////
            this.DATA.inputNames = {'$v (m/s)$','$\dot{\theta} (rad/s)$','$\dot{\psi} (rad/s)$'};
            this.DATA.inputs(1:length(this.DATA.inputNames),ENV.currentStep) = newState(4:6);         % Record the control inputs
        end
    end
    methods
        % GET THE STATE UPDATE (USING ODE45)
        function [X] = UpdateLocalState(this,TIME,X0,u)
            % This function computes the state update for the current agent
            % using the ode45 function.

            % Is of the form: u = [v;a;psi_dot;theta_dot];
            X = X0;

            % DETERMINE THE INTEGRATION PERIOD
            if TIME.currentTime == TIME.timeVector(end)
                return
            else
                [~,Xset] = ode45(@(t,X) this.dynamics_nonholonomicDublinsCar(X,u),...
                    [0 TIME.dt],X0,...
                    odeset('RelTol',1e-2,'AbsTol',TIME.dt*1E-2));
                X = Xset(end,:)'; % Pass the state at the end of the sample period
            end
        end
        % DEFINE THE GLOBAL UPDATE PROCEDURE FOR A 3D DUBLINS CAR MODEL
        function [this] = GlobalUpdate(this,dt,eulerState)

            % Retrieve current properties
            p_k	= this.GetGLOBAL('position');
            v_k	= this.GetGLOBAL('velocity');
            q_k	= this.GetGLOBAL('quaternion');
            x_k	= this.GetGLOBAL('priorState');

            % USE THE 'RATE' STATES DIRECTLY
            v_k_plus = [eulerState(4);0;0];                         % Assumed forward in the local axes
            omega_k_plus   = [0;eulerState(5);eulerState(6)];             % Neglect roll only
            localAxisRates = (omega_k_plus - [0;x_k(5:6)])/dt;
            % MAP THE LOCAL RATES TO GLOBAL RATES AND INTEGRATE QUATERNION
            % Previous properties

            % PREVIOUS ROTATION-MATRIX
            %R_k = quat2rotm(quaternion_k');
            R_k = OMAS_geometry.quaternionToRotationMatrix(q_k);
            % THE GLOBAL AXIS RATES
            omega = R_k'*localAxisRates;
            % UPDATE THE QUATERNION POSE
            q_k_plus = OMAS_geometry.integrateQuaternion(q_k,omega,dt);
            % REDEFINE THE ROTATION MATRIX
            R_k_plus = quat2rotm(q_k_plus');
            % MAP THE LOCAL VELOCITY TO THE GLOBAL AXES
            v_k_plus = R_k_plus'*v_k_plus;
            p_k_plus = p_k + dt*v_k;
            % ///////////////// REASSIGN K+1 PARAMETERS ///////////////////
            [this] = this.GlobalUpdate_direct(...
                p_k_plus,...    % Global position at k plius
                v_k_plus,...    % Global velocity at k plus
                q_k_plus);      % Quaternion at k plus
        end
    end
    % STATIC FUNCTIONS
    methods (Static)
        % THE NON-LINEAR HOLONOMIC DUBLINS-CAR MODEL
        function [dX] = dynamics_nonholonomicDublinsCar(X,u)
            % The state:  [dx;dy;dz;dv;dpsi;dtheta]
            % The inputs: [a,theta_dot,psi_dot]
%             dX = zeros(6,1);  % GLOBAL SPACE
%             dX(1) = X(4)*cos(X(6))*cos(X(5));
%             dX(2) = X(4)*cos(X(6))*sin(X(5));
%             dX(3) = X(4)*sin(X(6));
%             dX(4) = a;          % Acceleration
%             dX(5) = psi_dot;    % Yaw rate
%             dX(6) = theta_dot;  % Pitch rate
            dX = zeros(6,1);
            dX(1) = X(4);
            dX(2) = 0;
            dX(3) = 0;
            dX(4) = u(1);  % Acceleration
            dX(5) = u(2);  % Pitch rate
            dX(6) = u(3);  % Yaw rate
        end
    end
end

repo: OpenMAS


text:
function [mu,isig] = spm_affine_priors(typ)
% Distribution of the priors used in affine registration
%
% The parameters for this distribution were derived empirically from 227
% scans, that were matched to the ICBM space.
%_______________________________________________________________________
% Copyright (C) 2003-2011 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_affine_priors.m 4155 2011-01-11 15:22:39Z guillaume $

%% Values can be derived by...
%sn = spm_select(Inf,'.*seg_inv_sn.mat$');
%X  = zeros(size(sn,1),12);
%for i=1:size(sn,1),
%    p  = load(deblank(sn(i,:)));
%    M  = p.VF(1).mat*p.Affine/p.VG(1).mat;
%    J  = M(1:3,1:3);
%    V  = sqrtm(J*J');
%    R  = V\J;
%    lV      =  logm(V);
%    lR      = -logm(R);
%    P       = zeros(12,1);
%    P(1:3)  = M(1:3,4);
%    P(4:6)  = lR([2 3 6]);
%    P(7:12) = lV([1 2 3 5 6 9]);
%    X(i,:)  = P';
%end;
%mu   = mean(X(:,7:12));
%XR   = X(:,7:12) - repmat(mu,[size(X,1),1]);
%isig = inv(XR'*XR/(size(X,1)-1))

switch deblank(lower(typ))

case 'mni', % For registering with MNI templates...
    mu   = [0.0667 0.0039 0.0008 0.0333 0.0071 0.1071]';
    isig = 1e4 * [
        0.0902   -0.0345   -0.0106   -0.0025   -0.0005   -0.0163
       -0.0345    0.7901    0.3883    0.0041   -0.0103   -0.0116
       -0.0106    0.3883    2.2599    0.0113    0.0396   -0.0060
       -0.0025    0.0041    0.0113    0.0925    0.0471   -0.0440
       -0.0005   -0.0103    0.0396    0.0471    0.2964   -0.0062
       -0.0163   -0.0116   -0.0060   -0.0440   -0.0062    0.1144];

case 'imni', % For registering with MNI templates...
    mu   = -[0.0667 0.0039 0.0008 0.0333 0.0071 0.1071]';
    isig = 1e4 * [
        0.0902   -0.0345   -0.0106   -0.0025   -0.0005   -0.0163
       -0.0345    0.7901    0.3883    0.0041   -0.0103   -0.0116
       -0.0106    0.3883    2.2599    0.0113    0.0396   -0.0060
       -0.0025    0.0041    0.0113    0.0925    0.0471   -0.0440
       -0.0005   -0.0103    0.0396    0.0471    0.2964   -0.0062
       -0.0163   -0.0116   -0.0060   -0.0440   -0.0062    0.1144];

case 'rigid', % Constrained to be almost rigid...
    mu   = zeros(6,1);
    isig = eye(6)*1e8; % spm_affreg used 1e9

case 'subj', % For inter-subject registration...
    mu   = zeros(6,1);
    isig = 1e3 * [
        0.8876    0.0784    0.0784   -0.1749    0.0784   -0.1749
        0.0784    5.3894    0.2655    0.0784    0.2655    0.0784
        0.0784    0.2655    5.3894    0.0784    0.2655    0.0784
       -0.1749    0.0784    0.0784    0.8876    0.0784   -0.1749
        0.0784    0.2655    0.2655    0.0784    5.3894    0.0784
       -0.1749    0.0784    0.0784   -0.1749    0.0784    0.8876];

case 'eastern', % For East Asian brains to MNI...
    mu   = [0.0719   -0.0040   -0.0032    0.1416    0.0601    0.2578]';
    isig = 1e4 * [
        0.0757    0.0220   -0.0224   -0.0049    0.0304   -0.0327
        0.0220    0.3125   -0.1555    0.0280   -0.0012   -0.0284
       -0.0224   -0.1555    1.9727    0.0196   -0.0019    0.0122
       -0.0049    0.0280    0.0196    0.0576   -0.0282   -0.0200
        0.0304   -0.0012   -0.0019   -0.0282    0.2128   -0.0275
       -0.0327   -0.0284    0.0122   -0.0200   -0.0275    0.0511];

case 'none', % No regularisation...
    mu   = zeros(6,1);
    isig = zeros(6);

otherwise
    error(['"' typ '" not recognised as type of regularisation.']);
end


text:
function [ Sim QTot ] = computeAnimSimilarity( anim, method, ...
  varargin )
% Compute similarities between frames of a video sequence
%
% This routine computes:
%  - the pair infimum of CSFM if method = -2
%  - the best 2-view reconstruction in MSFM if method = 2
%  - the best 3-view reconstruction in MSFM if method = 3
%
% USAGE
%  [ Sim QTot ] = computeAnimSimilarity( anim, method, Sim )
%
% INPUTS
%  anim     - anim object (see GENERATETOYANIMATION for details)
%  method   - [2] uses pairs of frames, 3, uses triplets
%  varargin   - list of paramaters in quotes alternating with their values
%    for triadic comparison
%       'Sim2' pairwise similarity computed using method = 2
%       'perc' percentage of triplets to consider
%             ( only for speed purposes )
%
% OUTPUTS
%  Sim      - [ nFrame x nFrame ] or [ nFrame x nFrame x nFrame ]similarity
%             matrix based on total reconstruction error
%  QTot     - [ 4 x nFrame x nFrame ] optimal quaternions for pairwise
%
% EXAMPLE
%
% See also GENERATETOYANIMATION, VIEWANIMSIMILARITY
%
% Vincent's Structure From Motion Toolbox      Version 3.1.1
% Copyright (C) 2008-2011 Vincent Rabaud.  [vrabaud-at-cs.ucsd.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the GPL [see external/gpl.txt]

if ~isa(anim,'Animation'); error('anim must be of class Animation'); end
nFrame=anim.nFrame;

if nargin<2; method=2; end

% Compute the similarities
switch abs( method )
  case 2
    anim = anim.centerW();

    [ Sim, QTot ] = computeCsfmInfimum(anim.W);

    if method==2 % Best 2-view reconstruction for MSFM
      Sim = Sim/2;
    end

  case 3
    %%%%%%%% using triplets of frames for MSFM
    [ Sim2 perc ] = getPrmDflt( varargin, { 'Sim2' [] 'perc' 1 } );
    Sim = cell(1,nFrame); for i=1:nFrame; Sim{i}=sparse(nFrame,nFrame); end
    isDone = cell(1,nFrame); for i=1:nFrame; isDone{i}=sparse(nFrame,nFrame); end
    minSpan = 5;

    ticId = ticStatus( 'Triadic Computed' );

    % compute the probability to choose another frame given a frame in the triplet
    % basically, we'll pick one frame, and the 2 other frames according to that
    % probability. This one simply is proportional to the pairwise reconstruction
    % error (3 frames are more likely to be together if, pairwise, they look alike)
    probS = Sim2;

    % remove everything too close to each other
    toep = logical(spdiags(ones(nFrame,2*minSpan+1,-minSpan:minSpan, nFrame,nFrame)));
    probS(toep) = 0;

    % compute the cumulative sum and normalize by one
    probS = cumsum(Sim2,2);
    probS = bsxfun(@rdivide,probS,probS(:,end));

    % create as many triplets as necessary
    nDone = 0;
    while nDone/nFrame^3*100 < perc
      % choose two other frames according to the distribution
      % in probS and far enough from each other
      sample = randSample(nFrame,3);
      i1 = sample(1); i2 = sample(2); i3 = sample(3);
      if isDone{i1}(i2,i3) || min( pdist( sample' ) )<=minSpan; continue; end

      [ disc disc err ] = computeSMFromW( anim.isProj, ...
        anim.W(:,:, [ i1 i2 i3]), 'K', anim.K, 'method',Inf);

      perm = [ i1 i2 i3; i1 i3 i2; i3 i1 i2; i3 i2 i1; i2 i1 i3; ...
        i2 i3 i1 ];
      for l=1:6
        Sim{perm(l,1)}(perm(l,2),perm(l,3)) = err(1);
        isDone{perm(l,1)}(perm(l,2),perm(l,3)) = 1;
        nDone = nDone + 1;
      end
      if rand()>0.99
        tocStatus(ticId,min(perc,100*nDone/nFrame^3) / perc);
      end
    end
end

author: vrabaud


text:
function S = mk_naive_struct(n,C)
%
% S = mk_naive_struct(Number_of_nodes, Class_node)
%
S = zeros(n);
S(C,setdiff(1:n,C)) = 1;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function slwritearray(A, filename)
%SLWRITEARRAY Writes an array to an array file
%
% $ Syntax $
%   - slwritearray(A, filename)
%
% $ Arguments $
%   - A:            The array to be written
%   - filename:     The filename of the array file
%
% $ Description $
%   - slwritearray(A, filename) writes an array A to the array file.
%
% $ History $
%   - Created by Dahua Lin, on Jul 26th, 2006
%

value_types = { ...
    'double', ...
    'single', ...
    'logical', ...
    'char', ...
    'int8', ...
    'uint8', ...
    'int16', ...
    'uint16', ...
    'int32', ...
    'uint32', ...
    'int64', ...
    'uint64'};

%% open file

fid = fopen(filename, 'w');
if fid <= 0
    error('sltoolbox:filefail', ...
        'Fail to open file %s', filename);
end

%% write header

% write tag
fwrite(fid, ['arr', 0], 'char');

% write value type and dimension number
[tf, typeidx] = ismember(class(A), value_types);
if ~tf
    error('Unknown type for A: %s', class(A));
end
d = ndims(A);
info = uint8([typeidx, d, 0, 0]);
fwrite(fid, info, 'uint8');

%% write size

siz = size(A);
fwrite(fid, uint32(siz), 'uint32');

%% write data

fwrite(fid, A, class(A));

%% close file

fclose(fid);

author: lmthang


text:
solinit=bvpinit(linspace(-1,1,20),@dropinit);
sol=bvp4c(@drop,@dropbc,solinit);
fill(sol.x,sol.y(1,:),[0.7,0.7,0.7])
axis([-1,1,0,1])
xlabel('x','FontSize',12)
ylabel('h','Rotation',0,'FontSize',12)


text:
%% PIVOT AND UNPIVOT: a user's guide
% The scope of these two functions is to organize a _flat_ dataset into a PIVOT table (and eventually undo it).
% This guide will illustrate how to use them.

%% CONCEPT: flat dataset
% A _flat_ dataset usually contains a column of values and "n" columns of features. So, each row
% hs its value with its features.
%
% *1. A _flat_ dataset with two features (Region and Year)*
load Demodata.mat
display([{'Year','Region','People'}; Demodata.ex1])
%% CONCEPT: pivot table
% A pivot table as intended by the T-SQL command _pivot_ or as implemented by spreadsheet softwares is a data
% visualization tool which groups a _flat_ dataset by two chosen dimension.
%
% *2. Format Demodata.ex1 into a _pivot_ table*
warning off warnPivot:funGroup      % we'll come back later on the warning
display(Pivot(Demodata.ex1))
%% PIVOT input constraints
% $$Pivot(In, Fun, noHeaders, Pad)$$
%
% Not every type of input is supported by PIVOT. Only a _flat_ dataset with max three columns is supported. The first
% column should be the grouping dimension which will appear as the column  header and the second column the dimension
% that will appear as the row header. The last column should contain the values that will be "intersected" and grouped.
%
% Here are listed all the constraints to the inputs (please mail me if something is left out):
%
    %%
    %
    % * IN should have 3 columns
    % * IN can be numeric or a cell array of scalars or strings if headers (no mixing)
    % * Nested cells are not allowed
    % * NaNs or empty cells in the first two columns of IN are not allowed
    % * If IN is a cell the column with the values must be a cell array of scalars
    % * FUN must be a function handle
    % * FUN applied to a vector must return a scalar

try Pivot([1,3;2,4]); catch; err = lasterror; display(err.message); end %#ok
try Pivot('try'); catch; err = lasterror; display(err.message); end %#ok
try Pivot({'1999',1,2;1999,1,2}); catch; err = lasterror; display(err.message); end %#ok
try Pivot({1,2,3;1,{2},3}); catch; err = lasterror; display(err.message); end %#ok
try Pivot({NaN,1,9;1,3,4}); catch; err = lasterror; display(err.message); end %#ok
try Pivot([1,NaN,9;1,3,4]); catch; err = lasterror; display(err.message); end %#ok
try Pivot({[],23,9;1,3,4}); catch; err = lasterror; display(err.message); end %#ok
try Pivot({23,[],9;1,3,4}); catch; err = lasterror; display(err.message); end %#ok
try Pivot({23,1,9;1,3,'9'}); catch; err = lasterror; display(err.message); end %#ok
try Pivot(Demodata.ex1(2:end,:),'sum'); catch; err = lasterror; display(err.message); end %#ok
try Pivot(Demodata.ex1(2:end,:),@diff); catch; err = lasterror; display(err.message); end %#ok
%% PIVOT examples
% How to use this function:
%
% *1. IN numeric or cellarray of scalars. No grouping involved. Padded by default with NaN*
disp(Pivot(Demodata.ex2))
disp(Pivot(num2cell(Demodata.ex2))) % OUT will be automatically converted to numeric
%%
% *2. Grouping involved using default FUN (@sum) or @min*
disp(Pivot(Demodata.ex3))
disp(Pivot(Demodata.ex3,@min))
%%
% *3. OUT without headers*
disp(Pivot(Demodata.ex3,[],true))
%%
% *4. Pad with zeros or with '#' and no headers (OUT will be cell).*
disp(Pivot(Demodata.ex3,[],false,0))
disp(Pivot(Demodata.ex3,[],1,'#'))
%%
% *5. The same syntaxes could be used with IN headers as cellarrays of strings*
disp(Pivot(Demodata.ex4))
disp(Pivot(Demodata.ex4,@max))
disp(Pivot(Demodata.ex4,@max,1))
disp(Pivot(Demodata.ex4,@max,1,'???'))
%%
% *6. Full output*
[Out,colHeader,rowHeader,Settings] = Pivot(Demodata.ex4,@max,1,'???');
disp([{'colHeader','rowHeader'};[colHeader,rowHeader]])
disp('    Settings')
disp(Settings)
%% PIVOT warning
% We set *warning off warnPivot:funGroup* before. This was meant to avoid the following message
disp(lastwarn)
    %%
    % being thrown each time PIVOT is applied to a dataset which has multiple values on the
    % same intersection between the two grouping dimensions. The warning is meant
    % to remember the user that some values are being grouped
%% UNPIVOT input constraints
% $$unPivot(In, dim, rmPad)$$
%
% Only _pivot_ tables are accepted as inputs. The first elements of IN should be NaN. The first row and the
% first column (except the IN(1) element) are the headers/features that will define the values in
% IN(2:end,2:end).
%
% Here are listed all the constraints to the inputs (please mail me if something is left out):
%
    %%
    %
    % * IN should be a _pivot_ table
    % * IN(1) must be NaN
    % * Nested cells are not allowed
    % * DIM can be 1 or 2 or empty
    % * RMPAD can be a char or a scalar
try unPivot([1,3]); catch; err = lasterror; display(err.message); end %#ok
try unPivot({'G1',1999;'G2',24}); catch; err = lasterror; display(err.message); end %#ok
try unPivot({NaN,1999;'G2',{24}}); catch; err = lasterror; display(err.message); end %#ok
try unPivot(Demodata.ex5,3); catch; err = lasterror; display(err.message); end %#ok
try unPivot(Demodata.ex5,'1'); catch; err = lasterror; display(err.message); end %#ok
try unPivot(Demodata.ex5,[],{23}); catch; err = lasterror; display(err.message); end %#ok
%% UNPIVOT examples
% How to use this function:
%
% *1. Simplest case*
    %%
    % IN numeric or cellarray of scalars. Removes by default the NaNs. Sorts by default
    % according the column headers (second column in the _unPivoted_ table)
A = unPivot(Demodata.ex5); disp(A)
B = unPivot(Demodata.ex5,[]);
C = unPivot(Demodata.ex5,1);
if isequal(A,B,C); fprintf('A, B and C are equivalent'); else error('something''s wrong'); end
%%
% *2. Sort by row header (first column in the _unPivoted_ table)*
disp(unPivot(Demodata.ex5,2))
%%
% *3. Remove padded value other than NaN*
disp(unPivot(Demodata.ex5,1,3174))
%%
% *4. Same syntaxes can be applied to a cell array. Empty cells removed by default.*
disp(unPivot(Demodata.ex6))
disp(unPivot(Demodata.ex6,2))
disp(unPivot(Demodata.ex6,[],20))
disp(unPivot(Demodata.ex6,2,'???'))
%% Performance issues
% Lets test the performances of PIVOT and UNPIVOT.
    %%
    % First we create our Test input (which isn't included in Demodata since its about 35 MB)
    Nums1 = ceil(rand(100,1)*100);
    Nums2 = num2str(ceil(rand(100,1)*100));
    Test = [cellstr(strcat('Year', num2str(Nums1+1900))),...
            cellstr(strcat('Group', Nums2)),...
            num2cell(Nums1)];
    Test = repmat(Test,1000,1);
    %%
    % Then using the profile...
    profile on
    Testpivoted = Pivot(Test,[],[],' ');
    Testback = unPivot(Testpivoted,1,' ');
    profile off
    stats = profile('info');
    stats = [{stats.FunctionTable.FunctionName}',{stats.FunctionTable.TotalTime}'];
    disp(stats([strmatch('Pivot',stats(:,1)); strmatch('unPivot',stats(:,1))],:))
    %%
    % While UNPIVOT is fast enough PIVOT is slowed down by cell2mat like behaviour.
    % Here is implemented cell2floatmod which is a compact version of the cell2float by Jos(10584)
    %%
    % Also notice that the Test that was originally created is different from the unpivoted
    % version of the same Test because of the grouping.
    fprintf('Are they equal: ''%d''', isequalwithequalnans(Test, Testback))
%% Other Pivot-like functions
% *Statistics Toolbox owners:*
%%
% pivoting could be accomplished on a dataset obj by the functions _stack_ and _unstack_.
% These functions use a similar engine approach to Pivot and unPivot and the performances
% lie on the same scale but the syntax is different.
%%
% *FEX submission 20963-reshape-a-matrix by Dimitri Shvorob (ID:17777):*
%%
% The functions are called _wide2tall_ and _tall2wide_ and are probably the foundations for the
% statistics TB functions because of the syntax and the names. Has drawbacks on performance.

author: Sable


text:
% High fiber diet defintion. For details, please see https://www.vmh.life/#nutrition
% Ines Thiele 2016-2019

Diet = {%'Reaction'	'Flux Value'
'EX_etoh(e)'	'22.79219595'
'EX_h2o(e)'	'181671.3368'
'EX_caro(e)'	'0.013131607'
'EX_retinol(e)'	'1.099662212'
'EX_thm(e)'	'7.235600966'
'EX_adpcbl(e)'	'2.56714E-06'
'EX_ribflv(e)'	'0.006894923'
'EX_nac(e)'	'0.479928584'
'EX_ncam(e)'	'0.216377372'
'EX_pnto_R(e)'	'0.054851141'
'EX_pydam(e)'	'0.00499406'
'EX_pydxn(e)'	'0.004994745'
'EX_pydx(e)'	'0.005054977'
'EX_btn(e)'	'0.000479218'
'EX_10fthf(e)'	'0.000254697'
'EX_5mthf(e)'	'0.000261905'
'EX_thf(e)'	'0.000270786'
'EX_ascb_L(e)'	'1.513851945'
'EX_vitd3(e)'	'1.53391E-05'
'EX_avite1(e)'	'0.068538616'
'EX_phyQ(e)'	'0.000936308'
'EX_ca2(e)'	'33.02011078'
'EX_cl(e)'	'231.4450117'
'EX_k(e)'	'116.4477739'
'EX_mg2(e)'	'31.15531784'
'EX_na1(e)'	'147.9118802'
'EX_pi(e)'	'22.33606621'
'EX_cu2(e)'	'0.04732792'
'EX_fe2(e)'	'0.424456084'
'EX_fe3(e)'	'0.424456084'
'EX_mn2(e)'	'0.174378247'
'EX_zn2(e)'	'0.244799633'
'EX_mnl(e)'	'2.691195893'
'EX_xylt(e)'	'0.276050998'
'EX_lcts(e)'	'26.57468014'
'EX_malt(e)'	'0.803397055'
'EX_sucr(e)'	'137.1179664'
'EX_fru(e)'	'157.3272926'
'EX_gal(e)'	'1.010236134'
'EX_cellul(e)'	'0.237688253'
'EX_ala_L(e)'	'59.06129761'
'EX_arg_L(e)'	'40.67170675'
'EX_asp_L(e)'	'80.1605212'
'EX_cys_L(e)'	'13.82679835'
'EX_glu_L(e)'	'130.8154595'
'EX_gly(e)'	'60.78249181'
'EX_urate(e)'	'6.319184488'
'EX_his_L(e)'	'17.21966348'
'EX_ile_L(e)'	'41.93225125'
'EX_leu_L(e)'	'63.31791856'
'EX_lys_L(e)'	'46.34676332'
'EX_met_L(e)'	'12.66995194'
'EX_phe_L(e)'	'29.69793406'
'EX_pro_L(e)'	'50.90466905'
'EX_ser_L(e)'	'52.11494434'
'EX_thr_L(e)'	'38.25894566'
'EX_trp_L(e)'	'6.192991854'
'EX_tyr_L(e)'	'21.2683752'
'EX_val_L(e)'	'50.97864807'
'EX_dca(e)'	'0.043443556'
'EX_ddca(e)'	'0.101851489'
'EX_ttdca(e)'	'0.787507202'
'EX_ttdcea(e)'	'0.031285071'
'EX_ptdca(e)'	'0.072103355'
'EX_hdca(e)'	'34.44714331'
'EX_hpdca(e)'	'0.202603374'
'EX_ocdca(e)'	'9.424142226'
'EX_ocdcea(e)'	'130.2124615'
'EX_lnlc(e)'	'66.47282487'
'EX_lnlnca(e)'	'22.76317966'
'EX_strdnc(e)'	'0'
'EX_arach(e)'	'0.667592228'
'EX_CE2510(e)'	'1.006715204'
'EX_CE4843(e)'	'0.038033671'
'EX_arachd(e)'	'0.35958741'
'EX_docosac(e)'	'0.394315728'
'EX_doco13ac(e)'	'0.049887472'
'EX_adrn(e)'	'0'
'EX_clpnd(e)'	'0.050167498'
'EX_crvnc(e)'	'0.384328462'
'EX_but(e)'	'0.07233298'
'EX_octa(e)'	'0.015502414'
'EX_chsterol(e)'	'1.053035858'
'EX_sbt_D(e)'	'6.5699521'
'EX_glc_D(e)'	'120.210995'
'EX_hdcea(e)'	'3.616984668'
'EX_lgnc(e)'	'0.166975493'
'EX_fol(e)'	'0.000272645'
'EX_strch1(e)'	'45.11396195'
'EX_i(e)'	'0.003660628'
'EX_starch1200(e)'	'0.417683121'

};

Diet = regexprep(Diet,'EX_','Diet_EX_');
Diet = regexprep(Diet,'\(e\)','\[d\]');

author: opencobra

repo: cobratoolbox


text:
function [ y, j, ierror ] = yj_check_republican ( y, j )

%*****************************************************************************80
%
%% YJ_CHECK_REPUBLICAN checks a Republican YJ date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, J, the YJ date.
%
%    Output, integer IERROR, is 0 if no serious error was found in
%    the date, and 1 otherwise.
%

%
%  Check the year.
%
  [ y, ierror ] = y_check_republican ( y );

  if ( ierror ~= 0 )
    return
  end
%
%  Make sure J is not too small or too big.
%
  [ y, j ] = j_borrow_republican ( y, j );

  [ y, j ] = j_carry_republican ( y, j );

  return
end

author: johannesgerer


text:
function [inpC,p] = regCorrContrast(inp,Limit,pinit);
% regCorrContrast - removes mean and normalizes by the standard deviation,
%                both estimated by minimizing a robust error measure
%                between the sampled histogram and a mixture of 2 gaussians.
%
% [inpC,p] = regCorrContrast(inp, <pinit>);
%
% INPUT:
% - inp: set of original inplanes
% - pinit: initial parameters of the two best fitting Gaussians (optional)
%          (pinit = [mu1 sigma1^2 mu2 sigma2^2 p1 p2])
%
% Oscar Nestares - 5/99
%

if nargin<2
   Limit = 2;
end

% building the histogram
II = find(~isnan(inp));
[h x] = regHistogram(double(inp(II)), 256);

% normalizing the histogram
h = h/(sum(h)*mean(diff(x)));

% initial parameters: if they are not specified, chosoe the first
% gaussian with the actual mean and variance, and the second gaussian
% around 1 and with 1/10 of the actual variance, both with weights of 0.5
if nargin<3
   pinit = [mean(inp(II)) var(inp(II)) 1 var(inp(II))/10 0.5 0.5];
end

% minimizing the robust measure of the error
%p=fmins('regErrGaussRob', pinit, [], [], x, h);
p = fminsearch('regErrGaussRob', double(pinit), [], x, h);

% selecting the mean closer to 1
if abs(p(1)-1)>abs(p(3)-1)
   mu = p(3); sigma2 = p(4);
else
   mu = p(1); sigma2 = p(2);
end

% renormalizing
inpC = (inp - mu)/sqrt(sigma2);

% saturating for low and high values
%Limit = 4;  %%% This is a reasonable value, 2*std (now std=1)
Low = inpC < -Limit;
High = inpC > Limit;
inpC = inpC .*((~Low) & (~High)) + (-Limit)*Low + Limit*High;

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function ezaxes(command)
    %
    %   GUI input for axes
    %
    %   typing 'ezaxes' from the command line opens an
    %   interactive dialogue box, allowing the user to specify
    %   grids, axes size, fonts, scaling etc.
    %   using pulldown menus.
    %
    %
    %   Richard G. Cobb    3/96
    %   rcobb@afit.af.mil
    %

    if nargin < 1
        %
        figold=gcf;
        haold=gca;
        newfig=figure_w_normalized_uicontrolunits('Color',[0.9 0.9 0.9],'NumberTitle','off','Name',...
            'EZ AXES',...
            'NextPlot','new',...
            'Units','Normalized',...
            'Position',[.1 .1  .4 .3]);
        sfactor=1;
        axis(axis);
        axis off
        line([.52 .52],[-.2 1.2],'Color',[0 1 0],'LineWidth',3);
        uicontrol('Parent',newfig,...
            'style','pushbutton',...
            'string','Close',...
            'units','normalized',...
            'position',[.1 .85 .15 .1*sfactor],...
             'Callback','delete(gcf)');

        uicontrol('Parent',newfig,...
            'style','pushbutton',...
            'string','Refresh',...
            'units','normalized',...
            'position',[.3 .85 .15 .1*sfactor],...
             'Callback','ezarrow(98)');

        ui(1)=figold; bi(1)=haold;
        % Box Button
        ui(2)=0;
        bi(2) = uicontrol('Style','pushbutton','Units','Normalized',...
            'Position',[.1 .7 .2 .1*sfactor],...
            'String','Boxed',...
            'Visible','on',...
            'Value',0,...
            'Callback','ezaxes(2)');
        %% grid button
        ui(3)=0;
        bi(3) = uicontrol('Style','pushbutton','Units','Normalized',...
            'Position',[.1 .55 .2 .1*sfactor],...
            'String','Grid',...
            'Visible','on',...
            'Value',1,...
            'Callback','ezaxes(3)');
        % axes menu
        T4=['Equal ';'Square';'Normal';'Image ';'On    ';'Off   '];
        ui(4) = uicontrol('Style','Popup','Units','Normalized',...
            'Position',[.1 .4 .2 .1*sfactor],...
            'String',T4,...
            'Visible','off',...
            'HorizontalAlignment', 'left',...
            'Callback','ezaxes(24)');
        bi(4) = uicontrol('Style','pushbutton','Units','Normalized',...
            'Position',[.1 .4 .2 .1*sfactor],...
            'String','Axes',...
            'Visible','on',...
            'Value',0,...
            'Callback','ezaxes(4)');

        % Tick direction button
        ui(5)=0;
        bi(5) = uicontrol('Style','pushbutton','Units','Normalized',...
            'Position',[.1 .25 .2 .1*sfactor],...
            'String','Tick Direction',...
            'Visible','on',...
            'Value',0,...
            'Callback','ezaxes(5)');

        % Ticklength edit
        ui(6)=uicontrol('Style','edit','Units','Normalized',...
            'Position',[.1 .1 .2 .1*sfactor],...
            'String',' ',...
            'Visible','off',...
            'Value',0,...
            'Callback','ezaxes(26)');
        bi(6) = uicontrol('Style','pushbutton','Units','Normalized',...
            'Position',[.1 .1 .2 .1*sfactor],...
            'String','Tick Length',...
            'Visible','on',...
            'Value',0,...
            'Callback','ezaxes(6)');

        % grid line style menu
        T7=[' - ';' --';' : ';' -.'];
        ui(7) = uicontrol('Style','Popup','Units','Normalized',...
            'Position',[.32 .7 .2 .1*sfactor],...
            'String',T7,...
            'Visible','off',...
            'HorizontalAlignment', 'left',...
            'Callback','ezaxes(27)');
        bi(7) = uicontrol('Style','pushbutton','Units','Normalized',...
            'Position',[.32 .7 .2 .1*sfactor],...
            'String','GridLines',...
            'Visible','on',...
            'Value',0,...
            'Callback','ezaxes(7)');

        % Font name menu
        T8=['Avant Garde           ';...
            'Bookman               ';...
            'Courier               ';...
            'Helvetica             ';...
            'Helvetica Narrow      ';...
            'New Century Schoolbook';...
            'Palatino              ';...
            'Symbol                ';...
            'Times                 ';...
            'Zapf Chancery         ';...
            'Zapf Dingbats         '];
        ui(8) = uicontrol('Style','Popup','Units','Normalized',...
            'Position',[.32 .55 .2 .1*sfactor],...
            'String',T8,...
            'Visible','off',...
            'HorizontalAlignment', 'left',...
            'Callback','ezaxes(28)');
        bi(8) = uicontrol('Style','pushbutton','Units','Normalized',...
            'Position',[.32 .55 .2 .1*sfactor],...
            'String','FontName',...
            'Visible','on',...
            'Value',0,...
            'Callback','ezaxes(8)');

        % font size menu
        T9=['10';'12';'16';'20';'24';'30'];
        ui(9) = uicontrol('Style','Popup','Units','Normalized',...
            'Position',[.32 .4 .2 .1*sfactor],...
            'String',T9,...
            'Visible','off',...
            'HorizontalAlignment', 'left',...
            'Callback','ezaxes(29)');
        bi(9) = uicontrol('Style','pushbutton','Units','Normalized',...
            'Position',[.32 .4 .2 .1*sfactor],...
            'String','FontSize',...
            'Visible','on',...
            'Value',0,...
            'Callback','ezaxes(9)');

        % font weight menu
        T10=['Normal ';...
            'Light  ';...
            'Demi   ';...
            'Bold   '];
        ui(10) = uicontrol('Style','Popup','Units','Normalized',...
            'Position',[.32 .25 .2 .1*sfactor],...
            'String',T10,...
            'Visible','off',...
            'HorizontalAlignment', 'left',...
            'Callback','ezaxes(30)');
        bi(10) = uicontrol('Style','pushbutton','Units','Normalized',...
            'Position',[.32 .25 .2 .1*sfactor],...
            'String','FontWeight',...
            'Visible','on',...
            'Value',0,...
            'Callback','ezaxes(10)');

        % font angle menu
        T11=['Normal ';'Italic ';'Oblique'];
        ui(11) = uicontrol('Style','Popup','Units','Normalized',...
            'Position',[.32 .1 .2 .1*sfactor],...
            'String',T11,...
            'Visible','off',...
            'HorizontalAlignment', 'left',...
            'Callback','ezaxes(31)');
        bi(11) = uicontrol('Style','pushbutton','Units','Normalized',...
            'Position',[.32 .1 .2 .1*sfactor],...
            'String','FontAngle',...
            'Visible','on',...
            'Value',0,...
            'Callback','ezaxes(11)');

        %% direction button
        ui(12)=0;
        bi(12) = uicontrol('Style','pushbutton','Units','Normalized',...
            'Position',[.8 .85 .05 .1*sfactor],...
            'String','X',...
            'Visible','on',...
            'Value',0,...
            'BackgroundColor',[ 0 0 0 ],...
            'ForegroundColor',[ 0 1 0 ],...
            'Callback','ezaxes(12)');
        %% flip button

        ui(13)=0;
        bi(13) = uicontrol('Style','pushbutton','Units','Normalized',...
            'Position',[.55 .7 .1 .1*sfactor],...
            'String','Flip',...
            'Visible','on',...
            'Value',0,...
            'Callback','ezaxes(13)');

        %% grid button
        ui(14)=0;
        bi(14) = uicontrol('Style','pushbutton','Units','Normalized',...
            'Position',[.7 .7 .1 .1*sfactor],...
            'String','Grid',...
            'Visible','on',...
            'Value',0,...
            'Callback','ezaxes(14)');

        %% scale button
        ui(15)=0;
        bi(15) = uicontrol('Style','pushbutton','Units','Normalized',...
            'Position',[.85 .7 .1 .1*sfactor],...
            'String','Scale',...
            'Visible','on',...
            'Value',0,...
            'Callback','ezaxes(15)');

        %% limits edit
        ui(16)=uicontrol('Style','edit','Units','Normalized',...
            'Position',[.55 .475 .44 .1*sfactor],...
            'String',vec2str(get(haold,'Xlim')),...
            'Visible','off',...
            'Value',0,...
            'Callback','ezaxes(36)');
        bi(16) = uicontrol('Style','pushbutton','Units','Normalized',...
            'Position',[.55 .475 .2 .1*sfactor],...
            'String','Limits',...
            'Visible','on',...
            'Value',0,...
            'Callback','ezaxes(16)');

        %% tick edits
        ui(17)=uicontrol('Style','edit','Units','Normalized',...
            'Position',[.55 .325 .44 .1*sfactor],...
            'String',vec2str(get(haold,'XTick')),...
            'Visible','off',...
            'Value',0,...
            'Callback','ezaxes(37)');
        bi(17) = uicontrol('Style','pushbutton','Units','Normalized',...
            'Position',[.55 .325 .2 .1*sfactor],...
            'String','Ticks',...
            'Visible','on',...
            'Value',0,...
            'Callback','ezaxes(17)');

        %% ticklabel edit
        ui(18)=uicontrol('Style','edit','Units','Normalized',...
            'Position',[.55 .175 .44 .1*sfactor],...
            'String',mat2str(get(haold,'XTickLabels')),...
            'Visible','off',...
            'Value',0,...
            'Callback','ezaxes(38)');
        bi(18) = uicontrol('Style','pushbutton','Units','Normalized',...
            'Position',[.55 .175 .2 .1*sfactor],...
            'String','Tick Labels',...
            'Visible','on',...
            'Value',0,...
            'Callback','ezaxes(18)');

        %%% text strings
        uicontrol(...
            'style','text',...
            'string','Direction:',...
            'units','normalized',...
            'HorizontalAlignment','Center',...
            'position',[.5 .85 .3 .1*sfactor],...
            'BackgroundColor',[ 0 0 0 ],...
            'ForegroundColor',[ 0 1 0 ]);

        set(newfig,'UserData',[ui' bi']);

    elseif command == 1

    elseif command == 2
        handles = get(gcf,'UserData');
        if strcmp(get(handles(1,2),'Box'),'off')
            set(handles(1,2),'Box','on')
        else
            set(handles(1,2),'Box','off')
        end
    elseif command == 3
        handles = get(gcf,'UserData');
        figure_w_normalized_uicontrolunits(handles(1,1))
        grid
    elseif command == 4
        handles = get(gcf,'UserData');
        set(handles(4,2),'Visible','off')
        set(handles(4,1),'Visible','on')
    elseif command == 24
        handles = get(gcf,'UserData');
        T4=['Equal ';'Square';'Normal';'Image ';'On    ';'Off   '];
        figure_w_normalized_uicontrolunits(handles(1,1))
        axis(deblank(lower(T4(get(handles(4,1),'Value'),:))))
        set(handles(4,2),'Visible','on')
        set(handles(4,1),'Visible','off')

    elseif command == 5
        handles = get(gcf,'UserData');
        if strcmp(get(handles(1,2),'TickDir'),'in')
            set(handles(1,2),'TickDir','out')
        else
            set(handles(1,2),'TickDir','in')
        end
    elseif command == 6
        handles = get(gcf,'UserData');
        set(handles(6,2),'Visible','off')
        set(handles(6,1),'String',vec2str(get(handles(1,2),'Ticklength')))
        set(handles(6,1),'Visible','on')
    elseif command == 26
        handles = get(gcf,'UserData');
        set(handles(1,2),'TickLength',eval(get(handles(6,1),'String')))
        set(handles(6,2),'Visible','on')
        set(handles(6,1),'Visible','off')
    elseif command == 7
        handles = get(gcf,'UserData');
        set(handles(7,2),'Visible','off')
        set(handles(7,1),'Visible','on')
    elseif command == 27
        handles = get(gcf,'UserData');
        T7=[' - ';' --';' : ';' -.'];
        set(handles(1,2),'GridLineStyle',T7(get(handles(7,1),'Value'),:))
        set(handles(7,2),'Visible','on')
        set(handles(7,1),'Visible','off')

    elseif command == 8
        handles = get(gcf,'UserData');
        set(handles(8,2),'Visible','off')
        set(handles(8,1),'Visible','on')
    elseif command == 28
        handles = get(gcf,'UserData');
        T8=['Avant Garde           ';...
            'Bookman               ';...
            'Courier               ';...
            'Helvetica             ';...
            'Helvetica Narrow      ';...
            'New Century Schoolbook';...
            'Palatino              ';...
            'Symbol                ';...
            'Times                 ';...
            'Zapf Chancery         ';...
            'Zapf Dingbats         '];
        set(handles(1,2),'FontName',deblank(T8(get(handles(8,1),'Value'),:)))
        set(handles(8,2),'Visible','on')
        set(handles(8,1),'Visible','off')

    elseif command == 9
        handles = get(gcf,'UserData');
        set(handles(9,2),'Visible','off')
        set(handles(9,1),'Visible','on')
    elseif command == 29
        T9=['10';'12';'16';'20';'24';'30'];
        handles = get(gcf,'UserData');
        set(handles(1,2),'FontSize',eval(T9(get(handles(9,1),'Value'),:)))
        set(handles(9,2),'Visible','on')
        set(handles(9,1),'Visible','off')

    elseif command == 10
        handles = get(gcf,'UserData');
        set(handles(10,2),'Visible','off')
        set(handles(10,1),'Visible','on')
    elseif command == 30
        T10=['Normal ';...
            'Light  ';...
            'Demi   ';...
            'Bold   '];
        handles = get(gcf,'UserData');
        set(handles(1,2),'FontWeight',deblank(T10(get(handles(10,1),'Value'),:)))
        set(handles(10,2),'Visible','on')
        set(handles(10,1),'Visible','off')

    elseif command == 11
        handles = get(gcf,'UserData');
        set(handles(11,2),'Visible','off')
        set(handles(11,1),'Visible','on')
    elseif command == 31
        T11=['Normal ';'Italic ';'Oblique'];
        handles = get(gcf,'UserData');
        set(handles(1,2),'FontAngle',deblank(T11(get(handles(11,1),'Value'),:)))
        set(handles(11,2),'Visible','on')
        set(handles(11,1),'Visible','off')

    elseif command == 12
        handles = get(gcf,'UserData');
        if strcmp(get(handles(12,2),'String'),'X')
            set(handles(12,2),'String','Y')
        elseif strcmp(get(handles(12,2),'String'),'Y')
            set(handles(12,2),'String','Z')
        else
            set(handles(12,2),'String','X')
        end

    elseif command == 13
        handles = get(gcf,'UserData');
        if strcmp(get(handles(12,2),'String'),'X')
            if strcmp(get(handles(1,2),'XDir'),'normal')
                set(handles(1,2),'XDir','reverse')
            else
                set(handles(1,2),'XDir','normal')
            end
        elseif strcmp(get(handles(12,2),'String'),'Y')
            if strcmp(get(handles(1,2),'YDir'),'normal')
                set(handles(1,2),'YDir','reverse')
            else
                set(handles(1,2),'YDir','normal')
            end
        else
            if strcmp(get(handles(1,2),'ZDir'),'normal')
                set(handles(1,2),'ZDir','reverse')
            else
                set(handles(1,2),'ZDir','normal')
            end
        end

    elseif command == 14
        handles = get(gcf,'UserData');
        if strcmp(get(handles(12,2),'String'),'X')
            if strcmp(get(handles(1,2),'XGrid'),'on')
                set(handles(1,2),'XGrid','off')
            else
                set(handles(1,2),'XGrid','on')
            end
        elseif strcmp(get(handles(12,2),'String'),'Y')
            if strcmp(get(handles(1,2),'YGrid'),'on')
                set(handles(1,2),'YGrid','off')
            else
                set(handles(1,2),'YGrid','on')
            end

        else
            if strcmp(get(handles(1,2),'ZGrid'),'on')
                set(handles(1,2),'ZGrid','off')
            else
                set(handles(1,2),'ZGrid','on')
            end

        end
    elseif command == 15
        handles = get(gcf,'UserData');
        if strcmp(get(handles(12,2),'String'),'X')
            if strcmp(get(handles(1,2),'XScale'),'linear')
                set(handles(1,2),'XScale','log')
            else
                set(handles(1,2),'XScale','linear')
            end
        elseif strcmp(get(handles(12,2),'String'),'Y')
            if strcmp(get(handles(1,2),'YScale'),'linear')
                set(handles(1,2),'YScale','log')
            else
                set(handles(1,2),'YScale','linear')
            end

        else
            if strcmp(get(handles(1,2),'ZScale'),'linear')
                set(handles(1,2),'ZScale','log')
            else
                set(handles(1,2),'ZScale','linear')
            end
        end
    elseif command == 16
        handles = get(gcf,'UserData');
        set(handles(16,2),'Visible','off')
        if strcmp(get(handles(12,2),'String'),'X')
            set(handles(16,1),'String',vec2str(get(handles(1,2),'Xlim')))
        elseif strcmp(get(handles(12,2),'String'),'Y')
            set(handles(16,1),'String',vec2str(get(handles(1,2),'Ylim')))
        else
            set(handles(16,1),'String',vec2str(get(handles(1,2),'Zlim')))
        end
        set(handles(16,1),'Visible','on')
    elseif command == 36
        handles = get(gcf,'UserData');
        if strcmp(get(handles(12,2),'String'),'X')
            set(handles(1,2),'Xlim',eval(get(handles(16,1),'String')))
        elseif strcmp(get(handles(12,2),'String'),'Y')
            set(handles(1,2),'Ylim',eval(get(handles(16,1),'String')))
        else
            set(handles(1,2),'Zlim',eval(get(handles(16,1),'String')))
        end

        set(handles(16,2),'Visible','on')
        set(handles(16,1),'Visible','off')

    elseif command == 17
        handles = get(gcf,'UserData');
        set(handles(17,2),'Visible','off')
        if strcmp(get(handles(12,2),'String'),'X')
            set(handles(17,1),'String',vec2str(get(handles(1,2),'XTick')))
        elseif strcmp(get(handles(12,2),'String'),'Y')
            set(handles(17,1),'String',vec2str(get(handles(1,2),'YTick')))
        else
            set(handles(17,1),'String',vec2str(get(handles(1,2),'ZTick')))
        end
        set(handles(17,1),'Visible','on')
    elseif command == 37
        handles = get(gcf,'UserData');
        if strcmp(get(handles(12,2),'String'),'X')
            set(handles(1,2),'XTick',eval(get(handles(17,1),'String')))
        elseif strcmp(get(handles(12,2),'String'),'Y')
            set(handles(1,2),'YTick',eval(get(handles(17,1),'String')))
        else
            set(handles(1,2),'ZTick',eval(get(handles(17,1),'String')))
        end

        set(handles(17,2),'Visible','on')
        set(handles(17,1),'Visible','off')

    elseif command == 18
        handles = get(gcf,'UserData');
        set(handles(18,2),'Visible','off')
        if strcmp(get(handles(12,2),'String'),'X')
            set(handles(18,1),'String',ms2strv(get(handles(1,2),'XTickLabels')))
        elseif strcmp(get(handles(12,2),'String'),'Y')
            set(handles(18,1),'String',ms2strv(get(handles(1,2),'YTickLabels')))
        else
            set(handles(18,1),'String',ms2strv(get(handles(1,2),'ZTickLabels')))
        end
        set(handles(18,1),'Visible','on')
    elseif command == 38
        handles = get(gcf,'UserData');
        if strcmp(get(handles(12,2),'String'),'X')
            set(handles(1,2),'XTickLabels',ms2strv(get(handles(18,1),'String')))
        elseif strcmp(get(handles(12,2),'String'),'Y')
            set(handles(1,2),'YTickLabels',ms2strv(get(handles(18,1),'String')))
        else
            set(handles(1,2),'ZTickLabels',ms2strv(get(handles(18,1),'String')))
        end

        set(handles(18,2),'Visible','on')
        set(handles(18,1),'Visible','off')

    elseif command == 19

    elseif command == 20
    elseif command == 98
        handles = get(gcf,'UserData');
        figure_w_normalized_uicontrolunits(handles(1,1))
        refresh
    elseif command == 99
    end

author: CelsoReyes


text:
function test_ft_determine_units

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_determine_units

if nargout
  % assume that this is called by RUNTESTS
  tests = functiontests(localfunctions);
else
  % assume that this is called from the command line
  func = localfunctions;
  for i=1:numel(func)
    fprintf('evaluating %s\n', func2str(func{i}));
    feval(func{i});
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function testElec(testCase)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

elec = [];
elec.elecpos = randn(32,3);
for i=1:32
  elec.elecpos(i,:) = 100*elec.elecpos(i,:) ./ norm(elec.elecpos(i,:));
end
elec.unit = 'mm';

elec_mm = ft_determine_units(rmfield(ft_convert_units(elec, 'mm'), 'unit'));
elec_cm = ft_determine_units(rmfield(ft_convert_units(elec, 'cm'), 'unit'));
elec_m  = ft_determine_units(rmfield(ft_convert_units(elec, 'm'), 'unit'));

assert(strcmp(elec_mm.unit, 'mm'));
assert(strcmp(elec_cm.unit, 'cm'));
assert(strcmp(elec_m.unit, 'm'));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function testSphere(testCase)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

headmodel = [];
headmodel.r = 120;
headmodel.o = [0 0 40];
headmodel.unit = 'mm';

headmodel_mm = ft_determine_units(rmfield(ft_convert_units(headmodel, 'mm'), 'unit'));
headmodel_cm = ft_determine_units(rmfield(ft_convert_units(headmodel, 'cm'), 'unit'));
headmodel_m  = ft_determine_units(rmfield(ft_convert_units(headmodel, 'm'), 'unit'));

assert(strcmp(headmodel_mm.unit, 'mm'));
assert(strcmp(headmodel_cm.unit, 'cm'));
assert(strcmp(headmodel_m.unit, 'm'));


text:
function Y = var(w,flag)
   %VAR Variance of a waveform's data
   %  Y = var(W), where W is an N-dimensional waveform object returns the
   %  variance of the data within each waveform, normalized by N-1
   %
   %  Y = var(W,1) returns the variance within each waveform,
   %  normalized by N.
   %
   %  Y = var(W,weights), where weights is a vector the same length as the
   %  number of data samples, applies the weights to the variance.
   %
   % See the help for MATLAB's built in VAR for more details.
   %
   % Note: NAN values are completely ignored.
   %
   %  See also WAVEFORM/MEAN, WAVEFORM/MEDIAN, WAVEFORM/STD, NANVAR

   % AUTHOR: Celso Reyes, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   Y = zeros(size(w));

   %if the statistics toolbox is installed, use the builtin nanvar function to
   %ignore NaN values during the variance calculation.
   if ~isempty(ver('stats'))
      if exist('flag','var')
         for n = 1:numel(Y)
            Y(n) = nanvar(w(n).data,flag);
         end
      else
         for n = 1:numel(Y)
            Y(n) = nanvar(w(n).data);
         end
      end
   else
      % the statistics toolbox is not installed, so any nan values will have
      % to be dealt with (ignored) manually.
      if exist('flag','var')
         for n = 1:numel(Y)
            d = w(n).data;
            %d = d(~isnan(d));
            Y(n) = var(d(~isnan(d)),flag);
         end
      else
         for n = 1:numel(Y)
            d = w(n).data;
            %d = d(~isnan(d));
            Y(n) = var(d(~isnan(d)));
         end
      end
   end
end

author: geoscience-community-codes

repo: GISMO


text:
function [CD0] = DetermineProfileDrag(airfoildata,geo,panel)
%Determine the wing profile drag from the airfoils' drag polars

airfoil_r = geo.rootindex;
airfoil_t = geo.tipindex;

%Root dimensions
chord_r = geo.c_r;

%Determine Reynolds number at root; using log10 because regression was
%done this way
Re_r = log10(geo.Re_r);

%Convert back to degrees because the polynomial regression of the XFOIL
%drag data is in terms of degrees
alpha_r = (geo.i_r + geo.alpha)*180/pi;
b = cell2mat(airfoildata{2}(airfoil_r)); %Load the coefficients
c = cell2mat(airfoildata{3}(airfoil_r));
d = cell2mat(airfoildata{4}(airfoil_r));
k = cell2mat(airfoildata{5}(airfoil_r));
cd0_r = b(1) + b(2)*Re_r + b(3)*Re_r^2;
cd1_r = c(1) + c(2)*Re_r + c(3)*Re_r^2;
cd2_r = d(1) + d(2)*Re_r + d(3)*Re_r^2;
alpha_stall_r = k(1) + k(2)*Re_r + k(3)*Re_r^2;

%Test to see if airfoil is past stall angle of attack
if alpha_r > alpha_stall_r
    cd2_r = 2*cd2_r;
    %Double the last coefficient to account for the large increase in drag
    %past stall
end

%Tip dimensions
chord_t = geo.c_r*geo.taper;

%Determine Reynolds number at tip, use log10 because performed regression
%using log10
Re_t = log10(geo.Re_t);

%Convert back to degrees because the polynomial regression of the XFOIL
%drag data is in terms of degrees
alpha_t = (geo.i_r + geo.twist + geo.alpha)*180/pi;
b = cell2mat(airfoildata{2}(airfoil_t)); %Load the coefficients
c = cell2mat(airfoildata{3}(airfoil_t));
d = cell2mat(airfoildata{4}(airfoil_t));
k = cell2mat(airfoildata{5}(airfoil_t));
cd0_t = b(1) + b(2)*Re_t + b(3)*Re_t^2;
cd1_t = c(1) + c(2)*Re_t + c(3)*Re_t^2;
cd2_t = d(1) + d(2)*Re_t + d(3)*Re_t^2;
alpha_stall_t = k(1) + k(2)*Re_t + k(3)*Re_t^2;

%Test to see if airfoil is past stall angle of attack
if alpha_t > alpha_stall_t
    cd2_t = 2*cd2_t;
    %Double the last coefficient to account for the large increase in drag
    %past stall
end

%eta represents the fraction from 0 to 1 at the center of the panel where 0
%corresponds to y = 0 and 1 corresponds to y = b
eta = (2*(1:geo.ns)-1)/(2*geo.ns);
chord = chord_r + eta*(chord_t-chord_r);
alpha = alpha_r + eta*(alpha_t-alpha_r);
cd0 = cd0_r + eta*(cd0_t-cd0_r);
cd1 = cd1_r + eta*(cd1_t-cd1_r);
cd2 = cd2_r + eta*(cd2_t-cd2_r);

Deltay = zeros(1,geo.ns);
for i = 1:geo.ns
    Deltay(i) = panel(i,1).BV2(2) - panel(i,1).BV1(2);
end

D_q = (cd0 + cd1.*alpha + cd2.*alpha.^2).*chord.*Deltay;

CD0 = 2*sum(D_q,2)/geo.S;

repo: mcbench-benchmarks