text:
function GenCrossValidationData()
SPR = 0;

tag = sprintf('.SPR%d', SPR);

para.local.clean_wav_root = ChoosePath4OS({'D:\Data\NoiseData\Libri\LibriSpeech\dev-clean', '/media/xiaoxiong/OS/data1/G/REVERB_Challenge'});
para.local.clean_wav_ext = 'flac';
para.local.rir_wav_root = ChoosePath4OS({'D:\Data\NoiseData\ReverbEstimation\RIR_T60_1ch', '/media/xiaoxiong/OS/data1/G/REVERB_Challenge'});
para.local.rir_wav_ext = 'wav';
para.local.noise_wav_root = ChoosePath4OS({'D:\Data\NoiseData\musan\noise', '/media/xiaoxiong/OS/data1/G/REVERB_Challenge'});
para.local.noise_wav_ext = 'wav';
para.local.useFileName = 1;      % if set to 0, load all training data to memory. otherwise, only load file names.
para.local.parallel_wav_root = ChoosePath4OS({'D:\Data\NoiseData\Libri\LibriSpeech\dev-separation', '/media/xiaoxiong/OS/data1/G/REVERB_Challenge'});
para.local.parallel_wav_root = [para.local.parallel_wav_root tag];
para.singlePrecision = 1;

fs = 16000;
para.IO.DynamicMixture.fs = fs;
para.IO.DynamicMixture.nUtt4Iteration = 100;   % dynamically generate 10000 distorted sentences for each training iteration
para.IO.DynamicMixture.SPR_PDF = 'uniform';      % distribution of the power ratio between the two sources, defined as 10log10(P1/P2),
                                                 % where P1>P2 are the power of the two mixture signals
para.IO.DynamicMixture.SPR_para = [-1 1] * SPR;        % parameters of SPR PDF. If use uniform, it is the lowest and highest SPR allowed.
                                                 % if use normal, it is the mean and variance of SPR.
para.IO.DynamicMixture.addReverb = 0;            % whether to add reverberation to the mixed signal
para.IO.DynamicMixture.addNoise = 0;             % whether to add noise to the mixed signal
para.IO.DynamicMixture.SNR_PDF = 'uniform';      % distribution of SNR, choose [uniform|normal]. SNR is defined as the energy ratio P1/P2,
                                                 % where P1>P2 are the power of the two mixture signals
para.IO.DynamicMixture.SNR_para = [0 20];       % parameters of SNR PDF. If use uniform, it is the lowest and highest SNR allowed.
                                                 % if use normal, it is the mean and variance of SNR.

[base_data, para] = LoadWavRIRNoise_Libri(para, 1);
para.IO.DynamicMixture.inputFeature = para.IO.DynamicDistortion.inputFeature;
para.IO.DynamicMixture.fileReader = para.IO.DynamicDistortion.fileReader;

[data, para] = GenDynamicMixture(base_data, para);

distorted = data(1).data;
clean1 = data(2).data;
clean2 = data(3).data;

my_mkdir([para.local.parallel_wav_root '/mixed']);
my_mkdir([para.local.parallel_wav_root '/clean1']);
my_mkdir([para.local.parallel_wav_root '/clean2']);
for i=1:length(distorted)
    PrintProgress(i, length(distorted), 100);
    curr_distorted = double(distorted{i});
    curr_distorted = curr_distorted / max(abs(curr_distorted));
    distorted_filename = [para.local.parallel_wav_root '/mixed/' sprintf('utt%d.wav', i)];
    audiowrite(distorted_filename, curr_distorted, fs);

    curr_clean = double(clean1{i});
    curr_clean = curr_clean / max(abs(curr_clean));
    clean_filename = [para.local.parallel_wav_root '/clean1/' sprintf('utt%d.wav', i)];
    audiowrite(clean_filename, curr_clean, fs);

    curr_clean = double(clean2{i});
    curr_clean = curr_clean / max(abs(curr_clean));
    clean_filename = [para.local.parallel_wav_root '/clean2/' sprintf('utt%d.wav', i)];
    audiowrite(clean_filename, curr_clean, fs);
end

author: singaxiong

repo: SignalGraph


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is a demo execution of ECC image alignment algorithm
% It shows results in the case of a partial overlap and shows
% how the algortihm refines the results of feature-based matching
%
%
% 20/2/2013, Georgios Evangelidis, georgios.evangelidis@inria.fr
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

transform = 'homography';

NoI = 20; % number of iterations
NoL = 1;  % number of pyramid-levels

im_demo=imread('wallImage.png'); % ... or try your image

[A,B,C]=size(im_demo);

if C==3
    im_demo=rgb2gray(im_demo);
end
imu8 = im_demo;

im_demo=double(im_demo);

template_demo = imread('wallTemplate.png');
[A2, B2, C2] = size(template_demo);

if C2==3
    template_demo = rgb2gray(template_demo);
end
tempu8 = template_demo;
template_demo = double (template_demo);

%initialization by feature-based matching (use 1 iteration to see the
%initial registration)
    init = [0.7 -0.6 -34;
        0.7 0.8 -70;
        0.001 -0.001 1];

%call ECC
[results, final_warp, warped_image]=ecc(im_demo, template_demo, NoL, NoI, transform, init);
mask = spatial_interp(ones(A,B),final_warp, 'linear', transform, 1:B2, 1:A2);

figure; imshow(imu8);title('Input Image');
figure; imshow(tempu8); title('Template Image');
figure; imshow(uint8(warped_image)); title('Final Warped Image');

clear align1 align2;

% alignment based on initial warp
warped_init = spatial_interp(im_demo, init, 'linear', transform, 1:B2, 1:A2);
align1 = double(tempu8).*mask;
align1(:,:,2) = double(warped_init).*mask;
align1(:,:,3) = double(tempu8).*mask;
figure;imshow(uint8(align1)); title('Initial alignment');

% alignment based on final warp
align2 = align1;
align2(:,:,2) = double(warped_image).*mask;
figure;imshow(uint8(align2)); title('Final alignment');

author: Sable

repo: mcbench-benchmarks


text:
%% Import Script for ODF Data
%
% This script was automatically created by the import wizard. You should
% run the whoole script or parts of it in order to import your data. There
% is no problem in making any changes to this script.

%% Specify Crystal and Specimen Symmetries

% crystal symmetry
CS = {crystal symmetry};

% specimen symmetry
SS = {specimen symmetry};

% plotting convention
setMTEXpref('xAxisDirection',{xAxisDirection});
setMTEXpref('zAxisDirection',{zAxisDirection});

%% Specify File Names

% path to files
pname = {path to files};

% which files to be imported
fname = {file names};

%% Import the Data

% specify kernel
psi = {kernel name}('halfwidth',{halfwidth});

% load the ODF into the variable odf
odf = ODF.load(fname,CS,SS,{method},'kernel',psi,'resolution',{resolution},...
  'interface',{interface},{options});

%% Correct Data

rot = rotation.byEuler({phi1},{Phi},{phi2});
odf = rotate(odf,rot);

author: mtex-toolbox

repo: mtex

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


text:
%                                                            y = bwdpr1(Lden, b)
% BWDPR1  Solves "PROD_k L(pk,betak)' * y = b", where
%     L(p,beta) = eye(n) + tril(p*beta',-1).
%
% **********  INTERNAL FUNCTION OF SEDUMI **********
%
% See also sedumi, dpr1fact, fwdpr1

function y = bwdpr1(Lden, b) %#ok
%
% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA
%

disp('The SeDuMi binaries are not installed.')
disp('In Matlab, launch "install_sedumi" in the folder you put the SeDuMi files.')
disp('For more information see the file Install.txt.')
error(' ')

repo: REST


text:
function model = add_head_parts(model, X, num_parts, hi_res_size, low_res_size, a)
% Add subparts to the head part
%   model = add_head_parts(model, X, num_parts, hi_res_size, low_res_size, a)
%
% Return value
%   model           Updated person grammar model
%
% Arguments
%   model           Person grammar model
%   X               X is the person grammar's "head" symbol
%   num_parts       Number of subparts to add
%   hi_res_size     Size of high resolution subparts
%   low_res_size    Size of low resolution subparts
%   a               Subpart filter weight scale (used to reduce the influence
%                   of newly initialized parts that are added to the model)

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2011-2012 Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

% Terminal symbols for the two head filter subtypes
X_lf  = model.rules{X}(1).rhs(1);
X_rf  = model.rules{X}(2).rhs(1);
fid   = model.symbols(X_lf).filter;
% Add high resolution subparts
model = model_add_parts(model, X, 1, 2, fid, num_parts, hi_res_size, 1, a);

% Add two new rules: X -> X_lf | X_rf
% so that we can add low resolution parts to these
[model, rule] = model_add_struct_rule(model, X, X_lf, {[0 0 0]});
[model, rule] = model_add_struct_rule(model, X, X_rf, {[0 0 0]}, 'mirror_rule', rule);
model.blocks(rule.offset.blocklabel).learn = 0;
% Add low resolution subparts
model = model_add_parts(model, X, 3, 4, fid, num_parts, low_res_size, 0, a);
% Set the scale score parameters so that the low resolution parts have
% very low scores above the bottom octave of the pyramid
% (i.e., revent these rules from firing above the bottom octave)
model.blocks(model.rules{X}(3).loc.blocklabel).w = [0; -1000; -1000];
model.blocks(model.rules{X}(4).loc.blocklabel).w = [0; -1000; -1000];

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
function [ node_xyz, triangle_node ] = ply_to_tri_mesh ( ply_file_name )

%*****************************************************************************80
%
%% PLY_TO_TRI_MESH converts data from a PLY file to TRI_MESH data.
%
%  Discussion:
%
%    This routine reads a PLY file, searches for the data that defines
%    a polygonal surface, reduces polygons to triangles, and returns
%    a pair of arrays that define a TRI_MESH, that is a triangular
%    mesh of a surface in 3D.
%
%    The original version of this program had a mistake that meant it
%    did not properly triangulate files whose faces were not already triangular.
%    This has been corrected (JVB, 25 February 2007).
%
%    Glenn Ramsey pointed out and corrected a problem that occurred
%    with an uninitialized value of Type2, 27 August 2012.
%
%  Example:
%
%     [ node_xyz, triangle_node ] = ply_to_tri_mesh ( 'cow.ply' );
%
%     trisurf ( triangle_node', node_xyz(1,:), node_xyz(2,:), node_xyz(3,:) );
%
%     colormap ( gray );
%     axis equal;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 August 2012
%
%  Author:
%
%    Pascal Getreuer 2004
%    Modifications by John Burkardt
%
%  Parameters:
%
%    Input, string PLY_FILE_NAME, the name of the PLY file to be read.
%
%    Output, real NODE_XYZ(3,NODE_NUM), the coordinates of the nodes.
%
%    Output, integer TRIANGLE_NODE(3,TRIANGLE_NUM), the indices of nodes
%    that make up the triangles.
%
%  Local Parameters:
%
%    COMMENTS, any comments from the file.
%
%    ELEMENTCOUNT, the number of each type of element in file.
%
%    ELEMENTS, the element data.
%
%    LISTFLAG is 1 if there is a list type.
%
%    PROPERTYTYPES, the element property types.
%
%    SAMEFLAT is 1 if all types are the same.
%
%    SIZEOF, size in bytes of each type.
%

%
%  Open the input file in "read text" mode.
%
  [ fid, Msg ] = fopen ( ply_file_name, 'rt' );

  if ( fid == -1 )
    error ( Msg );
  end

  Buf = fscanf ( fid, '%s', 1 );

  if ( ~strcmp ( Buf, 'ply' ) )
    fclose ( fid );
    error('Not a PLY file.');
  end
%
%  Read the header.
%
  Position = ftell(fid);
  Format = '';
  NumComments = 0;
  Comments = {};
  NumElements = 0;
  NumProperties = 0;
  Elements = [];
  ElementCount = [];
  PropertyTypes = [];
  ElementNames = {};  % list of element names in the order they are stored in the file
  PropertyNames = [];  % structure of lists of property names

  while ( 1 )
%
%  Read a line from the file.
%
    Buf = fgetl ( fid );
    BufRem = Buf;
    Token = {};
    Count = 0;
%
%  Split the line into tokens.
%
    while ( ~isempty(BufRem) )

      [ tmp, BufRem ] = strtok(BufRem);
%
%  Count the tokens.
%
      if ( ~isempty ( tmp ) )
        Count = Count + 1;
        Token{Count} = tmp;
      end

    end
%
%  Parse the line.
%
    if ( Count )

      switch lower ( Token{1} )
%
%  Read the data format.
%
      case 'format'

        if ( 2 <= Count )

          Format = lower ( Token{2} );

          if ( Count == 3 & ~strcmp ( Token{3}, '1.0' ) )
            fclose ( fid );
            error('Only PLY format version 1.0 supported.');
          end
        end
%
%  Read a comment.
%
      case 'comment'

        NumComments = NumComments + 1;
        Comments{NumComments} = '';
        for i = 2 : Count
          Comments{NumComments} = [Comments{NumComments},Token{i},' '];
        end
%
%  Read an element name.
%
      case 'element'

        if ( 3 <= Count )

          if ( isfield(Elements,Token{2}) )
            fclose ( fid );
            error(['Duplicate element name, ''',Token{2},'''.']);
          end

          NumElements = NumElements + 1;
          NumProperties = 0;
          Elements = setfield(Elements,Token{2},[]);
          PropertyTypes = setfield(PropertyTypes,Token{2},[]);
          ElementNames{NumElements} = Token{2};
          PropertyNames = setfield(PropertyNames,Token{2},{});
          CurElement = Token{2};
          ElementCount(NumElements) = str2double(Token{3});

          if ( isnan(ElementCount(NumElements)) )
            fclose ( fid );
            error(['Bad element definition: ',Buf]);
          end

        else

          error(['Bad element definition: ',Buf]);

        end
%
%  Read an element property.
%
      case 'property'

        if ( ~isempty(CurElement) & Count >= 3 )

          NumProperties = NumProperties + 1;
          eval(['tmp=isfield(Elements.',CurElement,',Token{Count});'],...
            'fclose(fid);error([''Error reading property: '',Buf])');

          if ( tmp )
            error(['Duplicate property name, ''',CurElement,'.',Token{2},'''.']);
          end
%
%  Add property subfield to Elements.
%
          eval(['Elements.',CurElement,'.',Token{Count},'=[];'], ...
            'fclose(fid);error([''Error reading property: '',Buf])');
%
%  Add property subfield to PropertyTypes and save type.
%
          eval(['PropertyTypes.',CurElement,'.',Token{Count},'={Token{2:Count-1}};'], ...
            'fclose(fid);error([''Error reading property: '',Buf])');
%
%  Record property name order.
%
          eval(['PropertyNames.',CurElement,'{NumProperties}=Token{Count};'], ...
            'fclose(fid);error([''Error reading property: '',Buf])');

         else

           fclose ( fid );

           if ( isempty(CurElement) )
             error(['Property definition without element definition: ',Buf]);
           else
             error(['Bad property definition: ',Buf]);
           end

         end
%
%  End of header.
%
        case 'end_header'
          break;

      end
    end
  end
%
%  Set reading for specified data format.
%
  if ( isempty ( Format ) )
    warning('Data format unspecified, assuming ASCII.');
    Format = 'ascii';
  end

  switch Format

    case 'ascii'
      Format = 0;
    case 'binary_little_endian'
      Format = 1;
    case 'binary_big_endian'
      Format = 2;
    otherwise
      fclose ( fid );
      error(['Data format ''',Format,''' not supported.']);

  end
%
%  Read the rest of the file as ASCII data...
%
  if ( ~Format )
    Buf = fscanf ( fid, '%f' );
    BufOff = 1;
  else
%
%  ...or, close the file, and reopen in "read binary" mode.
%
    fclose ( fid );
%
%  Reopen the binary file as LITTLE_ENDIAN or BIG_ENDIAN.
%
    if ( Format == 1 )
      fid = fopen ( Path, 'r', 'ieee-le.l64' );
    else
      fid = fopen ( Path, 'r', 'ieee-be.l64' );
    end
%
%  Find the end of the header again.
%  Using ftell on the old handle doesn't give the correct position.
%
    BufSize = 8192;
    Buf = [ blanks(10), char(fread(fid,BufSize,'uchar')') ];
    i = [];
    tmp = -11;

    while ( isempty(i) )
%
%  Look for end_header + CR/LF
%
      i = findstr(Buf,['end_header',13,10]);
%
%  Look for end_header + LF
%
      i = [i,findstr(Buf,['end_header',10])];

      if ( isempty(i) )
        tmp = tmp + BufSize;
        Buf = [Buf(BufSize+1:BufSize+10),char(fread(fid,BufSize,'uchar')')];
      end

    end
%
%  Seek to just after the line feed
%
    fseek ( fid, i + tmp + 11 + (Buf(i + 10) == 13), -1 );

  end
%
%  Read element data.
%
%  PLY and MATLAB data types (for fread)
%
  PlyTypeNames = {'char','uchar','short','ushort','int','uint','float','double', ...
    'char8','uchar8','short16','ushort16','int32','uint32','float32','double64'};

  MatlabTypeNames = {'schar','uchar','int16','uint16','int32','uint32','single','double'};

  SizeOf = [1,1,2,2,4,4,4,8];

  for i = 1 : NumElements
%
%  get current element property information
%
    eval(['CurPropertyNames=PropertyNames.',ElementNames{i},';']);
    eval(['CurPropertyTypes=PropertyTypes.',ElementNames{i},';']);
    NumProperties = size(CurPropertyNames,2);
%
%  Read ASCII data.
%
    if ( ~Format )

      for j = 1 : NumProperties

        Token = getfield(CurPropertyTypes,CurPropertyNames{j});

        if ( strcmpi(Token{1},'list') )
          Type(j) = 1;
        else
          Type(j) = 0;
        end
%
%  Glenn Ramsey 20120827
%  Initialise Type2{} to prevent uninitialised value error.
%
        Type2{j} = '';

      end
%
%  Parse the buffer.
%
      if ( ~any(Type) )
%
%  No list types
%
        Data = reshape ( ...
          Buf(BufOff:BufOff+ElementCount(i)*NumProperties-1), ...
          NumProperties, ElementCount(i) )';

        BufOff = BufOff + ElementCount(i) * NumProperties;

      else

        ListData = cell(NumProperties,1);

        for k = 1 : NumProperties
          ListData{k} = cell(ElementCount(i),1);
        end
%
%  list type
%
        for j = 1 : ElementCount(i)
          for k = 1 : NumProperties

            if ( ~Type(k) )
              Data(j,k) = Buf(BufOff);
              BufOff = BufOff + 1;
            else
              tmp = Buf(BufOff);
              ListData{k}{j} = Buf(BufOff+(1:tmp))';
              BufOff = BufOff + tmp + 1;
            end

          end
        end

      end

%
%  Read binary data.
%
    else

%
%  Translate PLY data type names to MATLAB data type names
%
      ListFlag = 0;
      SameFlag = 1;

      for j = 1 : NumProperties

        Token = getfield(CurPropertyTypes,CurPropertyNames{j});
%
%  Non-list type.
%
        if ( ~strcmp(Token{1},'list' ) )

          tmp = rem(strmatch(Token{1},PlyTypeNames,'exact')-1,8)+1;

          if ( ~isempty(tmp) )

            TypeSize(j) = SizeOf(tmp);
            Type{j} = MatlabTypeNames{tmp};
            TypeSize2(j) = 0;
            Type2{j} = '';

            SameFlag = SameFlag & strcmp(Type{1},Type{j});

          else

            fclose(fid);
            error(['Unknown property data type, ''',Token{1},''', in ', ...
              ElementNames{i},'.',CurPropertyNames{j},'.']);

          end
%
%  List type.
%
        else

          if ( length(Token) == 3 )

            ListFlag = 1;
            SameFlag = 0;
            tmp = rem(strmatch(Token{2},PlyTypeNames,'exact')-1,8)+1;
            tmp2 = rem(strmatch(Token{3},PlyTypeNames,'exact')-1,8)+1;
            disp(tmp, tmp2)

            if ( ~isempty(tmp) & ~isempty(tmp2) )
              TypeSize(j) = SizeOf(tmp);
              Type{j} = MatlabTypeNames{tmp};
              TypeSize2(j) = SizeOf(tmp2);
              Type2{j} = MatlabTypeNames{tmp2};
              disp('set Type2{j}')
            else
              fclose(fid);
              error(['Unknown property data type, ''list ',Token{2},' ',Token{3},''', in ', ...
                ElementNames{i},'.',CurPropertyNames{j},'.']);
            end

          else

            fclose(fid);
            error(['Invalid list syntax in ',ElementNames{i},'.',CurPropertyNames{j},'.']);

          end

        end

      end
%
%  read file
%
      if ( ~ListFlag )
%
%  No list types, all the same type (fast)
%
        if ( SameFlag )

          Data = fread(fid,[NumProperties,ElementCount(i)],Type{1})';
%
%  No list types, mixed type.
%
        else

          Data = zeros(ElementCount(i),NumProperties);

          for j = 1 : ElementCount(i)
            for k = 1 : NumProperties
              Data(j,k) = fread(fid,1,Type{k});
            end
          end

        end

      else

        ListData = cell(NumProperties,1);

        for k = 1 : NumProperties
          ListData{k} = cell(ElementCount(i),1);
        end

        if ( NumProperties == 1 )

          BufSize = 512;
          SkipNum = 4;
          j = 0;
%
%  List type, one property (fast if lists are usually the same length)
%
          while ( j < ElementCount(i) )

            BufSize = min(ElementCount(i)-j,BufSize);
            Position = ftell(fid);
%
%  Read in BufSize count values, assuming all counts = SkipNum
%
            [Buf,BufSize] = fread(fid,BufSize,Type{1},SkipNum*TypeSize2(1));
%
%  Find first count that is not SkipNum
%
            Miss = find(Buf ~= SkipNum);
%
%  Seek back to after first count
%
            fseek(fid,Position + TypeSize(1),-1);

            if ( isempty(Miss) )
%
%  All counts are SkipNum.
%
              Buf = fread(fid,[SkipNum,BufSize],[int2str(SkipNum),'*',Type2{1}],TypeSize(1))';
%
%  Undo the last skip.
%
              fseek(fid,-TypeSize(1),0);

              for k = 1:BufSize
                ListData{1}{j+k} = Buf(k,:);
              end

              j = j + BufSize;
              BufSize = floor ( 1.5 * BufSize );

            else
%
%  Some counts are SkipNum.
%
              if ( 1 < Miss(1) )

                Buf2 = fread ( fid, [SkipNum,Miss(1)-1],[int2str(SkipNum),'*',Type2{1}],TypeSize(1))';

                for k = 1:Miss(1)-1
                  ListData{1}{j+k} = Buf2(k,:);
                end

                j = j + k;

              end
%
%  Read in the list with the missed count.
%
              SkipNum = Buf(Miss(1));
              j = j + 1;
              ListData{1}{j} = fread ( fid, [1,SkipNum],Type2{1} );
              BufSize = ceil ( 0.6 * BufSize );

            end
          end

        else
%
%  List type(s), multiple properties (slow)
%
          Data = zeros(ElementCount(i),NumProperties);

          for j = 1 : ElementCount(i)
            for k = 1 : NumProperties

              if ( isempty(Type2{k}) )
                Data(j,k) = fread(fid,1,Type{k});
              else
                tmp = fread(fid,1,Type{k});
                ListData{k}{j} = fread(fid,[1,tmp],Type2{k});
              end

            end
          end
        end
      end
    end
%
%  Put data into Elements structure
%
    for k = 1 : NumProperties

      if ( ( ~Format & ~Type(k) ) || (Format & isempty(Type2{k})) )
        eval(['Elements.',ElementNames{i},'.',CurPropertyNames{k},'=Data(:,k);']);
      else
        eval(['Elements.',ElementNames{i},'.',CurPropertyNames{k},'=ListData{k};']);
      end

    end

  end

  clear Data
  clear ListData;

  fclose ( fid );
%
%  Output the data as a triangular mesh pair.
%
%  Find vertex element field.
%
  Name = {'vertex','Vertex','point','Point','pts','Pts'};
  Names = [];

  for i = 1 : length(Name)

    if ( any ( strcmp ( ElementNames, Name{i} ) ) )
      Names = getfield(PropertyNames,Name{i});
      Name = Name{i};
      break;
    end

  end

  if ( any(strcmp(Names,'x')) & any(strcmp(Names,'y')) & any(strcmp(Names,'z')) )
    eval(['node_xyz=[Elements.',Name,'.x,Elements.',Name,'.y,Elements.',Name,'.z];']);
  else
    node_xyz = zeros(1,3);
  end
  node_xyz = node_xyz';

  triangle_node = [];

% Find face element field

  Name = {'face','Face','poly','Poly','tri','Tri'};
  Names = [];

  for i = 1 : length(Name)
    if ( any(strcmp(ElementNames,Name{i})) )
      Names = getfield(PropertyNames,Name{i});
      Name = Name{i};
      break;
    end
  end

  if ( ~isempty(Names) )
      % find vertex indices property subfield
    PropertyName = {'vertex_indices','vertex_indexes','vertex_index','indices','indexes'};

    for i = 1 : length(PropertyName)
      if ( any(strcmp(Names,PropertyName{i})) )
        PropertyName = PropertyName{i};
        break;
      end
    end
%
%  Convert face index list to triangular connectivity.
%
    if ( ~iscell(PropertyName) )

      eval(['FaceIndices=Elements.',Name,'.',PropertyName,';']);
      N = length(FaceIndices);
      triangle_node = zeros(3,N*2);
      Extra = 0;

      for k = 1 : N

        triangle_node(1:3,k) = FaceIndices{k}(1:3)';

        for j = 4 : length(FaceIndices{k})
          Extra = Extra + 1;
          triangle_node(1,N + Extra) = FaceIndices{k}(1);
          triangle_node(2,N + Extra) = FaceIndices{k}(j-1);
          triangle_node(3,N + Extra) = FaceIndices{k}(j);
        end

      end
%
%  Add 1 to each vertex value; PLY vertices are zero based.
%
      triangle_node = triangle_node(:,1:N+Extra) + 1;

    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% bsearch(x,var)
% Written by Aroh Barjatya
% Binary search for values specified in vector 'var' within data vector 'x'
% The data has to be pre-sorted in ascending or decending order
% There is no way to predict how the function will behave if there
% are multiple numbers with same value.
% returns the index values of the searched numbers

function index = bsearch(x,var)

xLen = length(x);
[xRow xCol] = size(x);
if x(1) > x(xLen)	% means x is in descending order
    if xRow==1
        x = fliplr(x);
    else
        x = flipud(x);
    end
    flipped = 1;
elseif x(1) < x(xLen)	% means x is in ascending order
    flipped = 0;
else
    'badly formatted data. Type ''help bsearch\'')';
    return;
end

for i = 1:length(var)
    low = 1;
    high = xLen;
    if var(i) <= x(low)
        index(i) = low;
        continue;
    elseif var(i) >= x(high)
        index(i) = high;
        continue;
    end
    flag = 0;
    while (low <= high)
        mid = round((low + high)/2);
        if (var(i) < x(mid))
            high = mid;
        elseif (var(i) > x(mid))
            low = mid;
        else
            index(i) = mid;
            flag = 1;
            break;
        end
        if (low == high - 1)
            break
        end
    end
    if (flag == 1)
        continue;
    end
    if (low == high)
        index(i) = low;
    elseif ((x(low) - var(i))^2 > (x(high) - var(i))^2)
        index(i) = high;
    else
        index(i) = low;
    end
end

if flipped
    index = xLen - index + 1;
end


text:
roi = [1 1  480 720];
%Region of interest
maxNumObj = 200; % Maximum number of objects to track
alarmCount = 75; % Max no of frames an object can remain stationary before alarm is raised
maxConsecutiveMiss = 7; %Min frames an object's centroid changes after which it is not tracked

% System object for reading the video with each frame of the type 'Single'
hVideoSrc = vision.VideoFileReader;
hVideoSrc.Filename = input('Enter the video to run:- ','s');
hVideoSrc.VideoOutputDataType = 'single';

% Offsets for drawing bounding boxes in original input video
PtsOffset = int32(repmat([roi(1), roi(2), 0, 0],[maxNumObj 1]));

% Converts RGB image to YCbCr
hColorConv = vision.ColorSpaceConverter('Conversion', 'RGB to YCbCr');

% Does background subtraction on 2 images
hAutothreshold = vision.Autothresholder('ThresholdScaleFactor', 1.0);

% Removes noise and really small blobs
hClosing = vision.MorphologicalClose('Neighborhood', strel('square',10));

% Finds the blobs in the segmented images, properties of blobs also
% specified
hBlob = vision.BlobAnalysis('MaximumCount', maxNumObj, 'ExcludeBorderBlobs', true);
hBlob.MinimumBlobArea = 100;
hBlob.MaximumBlobArea = 5000;

% Creating system objects for players with their locations
pos = [10 300 roi(3)+25 roi(4)+25];
hAbandonedObjects = vision.VideoPlayer('Name', 'Abandoned Objects', 'Position', pos);
pos(1) = 46+roi(3); % move the next viewer to the right
hAllObjects = vision.VideoPlayer('Name', 'All Objects', 'Position', pos);
pos = [80+2*roi(3) 300 roi(3)-roi(1)+25 roi(4)-roi(2)+25];
hThresholdDisplay = vision.VideoPlayer('Name', 'Threshold', 'Position', pos);
mPlayer1=vision.DeployableVideoPlayer('Location',[10,100]);
mPlayer2=vision.DeployableVideoPlayer('Location',[20,110]);
mPlayer3=vision.DeployableVideoPlayer('Location',[20,110]);
mPlayer4=vision.DeployableVideoPlayer('Location',[20,110]);
mPlayer5=vision.DeployableVideoPlayer('Location',[20,110]);

firsttime = true; % Initialisation for storing background
allBlobList=[]; % Used for tracking blobs, stores all information about all blobs
frame_no=0; % Current frame no

% Loop runs till all the video frames are completed
while ~isDone(hVideoSrc)
    frame_no=frame_no+1;
    Im = step(hVideoSrc); % Stores current frame in Im
    OutIm = Im(roi(2):end, roi(1):end, :); % Selects the region of interest from the original video
    YCbCr = step(hColorConv, OutIm); % Gives the YCbCr image of frame
    CbCr = complex(YCbCr(:,:,2), YCbCr(:,:,3)); % CbCr has the color components of YCbCr
    step(mPlayer2,YCbCr); % Dis[plays YCbCr image

    % Stores background
    if firsttime
        firsttime = false;
        BkgY = YCbCr(:,:,1);
        BkgCbCr = CbCr;
    end

    SegY = step(hAutothreshold, abs(YCbCr(:,:,1)-BkgY)); % Background subtraction on the luminosity part
    SegCbCr = abs(CbCr-BkgCbCr) > 0.05; % Background subtraction on chroma parts
    step(mPlayer3,SegY); % Luminosity subtracted image

    % Fill in small gaps in the detected objects and clubs the separated
    % image
    Segmented = step(hClosing, SegY | SegCbCr);
    step(mPlayer1,Segmented); % Foreground

    % Perform blob analysis
    [Area, Centroid, BBox] = step(hBlob, Segmented);
    %[x y hitCount latest_detected_frame blob_number starting_frame misscount]
    %[^ ^ ^         ^                      ^           ^              ^]
    for blob=1:size(Centroid,1)%Traverses through the list of Centroids of all the blobs
        %To map x and y to the largest multiple of 5 less than it
        %eg->(103,107)=(100,105)
        roundX=Centroid(blob,1)-mod(Centroid(blob,1),5);
        roundY=Centroid(blob,2)-mod(Centroid(blob,2),5);
        found=false;%to check if the given blob centoid already was there or not
        for x=1:size(allBlobList,1)
            if(allBlobList(x,1)==roundX && allBlobList(x,2)==roundY)%if the centroid was found
                allBlobList(x,3)=allBlobList(x,3)+1;%increasing the hit count
                found=true;
                allBlobList(x,4)=frame_no;%storing the latest detected frame no.
                allBlobList(x,5)=blob;%storing the blob number of the centroid of that frame
                allBlobList(x,7)=0;%since the blob as detected, miss count will be zero
            end
        end
        if(found==false)%else we need to add the centroid to our list of centroids
            allBlobList=[allBlobList;[roundX roundY 1 frame_no blob frame_no 0 ]];
        end
    end
    BlobCount = size(BBox,1);%to get the total no. o blobs detected in that frame

    BBoxOffset = BBox + int32(repmat([roi(1) roi(2) 0 0],[BlobCount 1]));
    Imr = insertShape(Im,'Rectangle',BBoxOffset,'Color','green');%inserting green rectangles for that detected blobs
    myImr=insertMarker(Imr,Centroid);%inserting marker at the centroid of each blob
    rowToRemove=[];%list of all the rows to be  removed since there miss count exceeded the maxMissCount
    for blob=1:size(allBlobList,1)
        %check for the blob to be abandoned
        if(allBlobList(blob,3)>alarmCount && allBlobList(blob,7)<=maxConsecutiveMiss && frame_no==allBlobList(blob,4) )
            %allBlobList(blob,7)=0;
            myImr=insertShape(myImr, 'FilledRectangle', BBox(allBlobList(blob,5), :), 'color','red', 'Opacity', 0.5);
            sound(randn(4096,1),8192);
        %increment the miss count
        else
            allBlobList(blob,7)=allBlobList(blob,7)+1;
            if(allBlobList(blob,7)>maxConsecutiveMiss)%if misscount>maxConsicutiveMiss delete it
                rowToRemove=[rowToRemove;blob];
            end
        end
    end
    %deleting the rows for misscount>maxConsicutiveMiss
    for row=1:size(rowToRemove,1)
        try
        allBlobList(row,:)=[];
        end
    end
    step(hAllObjects, myImr);
    % Display the segmented video
    SegBBox = PtsOffset;
end
release(hVideoSrc);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function createHdf5File(filename, datasetName, inputData, varargin)
	% Creates an HDF5 file at filename under given datasetName hierarchy and saves inputData.
	% Biafra Ahanonu
	% started: 2014.01.07
	% inputs
		%
	% outputs
		%

	ciapkg.hdf5.createHdf5File(filename, datasetName, inputData,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function prior = priorTest(priorType);

% PRIORTEST Run some tests on the specified prior.
% FORMAT
% DESC runs some checks on the specified prior (gradients etc.).
% ARG priorType : the prior type to check.
% RETURN prior : the prior structure that was created to test.
%
% SEEALSO : priorParamInit, priorExpandParam, priorExtractParam,
% priorLogProb, priorGradient
%
% COPYRIGHT : Neil D. Lawrence, 2003, 2004
%
% MODIFICATIONS : Antti Honkela, 2012, 2014

% PRIOR

prior.type = priorType;
prior = priorParamInit(prior);

% Set the parameters randomly.
params = priorExtractParam(prior);
params = randn(size(params))./sqrt(randn(size(params)).^2);
prior = priorExpandParam(prior, params);
if isfield(prior, 'isBounded') && prior.isBounded,
  prior.b = 1+9*rand(1);
  x = rand(1, 10) * (prior.b - prior.a) + prior.a;
else
  x = randn(1, 10);
end
epsilon = 1e-6;

if exist([priorType 'PriorGradientParams'])
  params = priorExtractParam(prior);
  origParams = params;
  for i = 1:length(params);
    params = origParams;
    params(i) = origParams(i) + epsilon;
    prior = priorExpandParam(prior, params);
    Lplus(i) = priorLogProb(prior, x);
    params(i) = origParams(i) - epsilon;
    prior = priorExpandParam(prior, params);
    Lminus(i) = priorLogProb(prior, x);
  end
  params = origParams;
  prior = priorExpandParam(prior, params);
  [void, names] = priorExtractParam(prior);
  gLDiff = .5*(Lplus - Lminus)/epsilon;
  g = priorGradient(prior, x);

  paramMaxDiff = max(max(abs(gLDiff-g)));
  if paramMaxDiff > 2*epsilon
    l = 0;
    for i = 1:length(names)
      if l < length(names{i})
        l = length(names{i});
      end
    end

    fprintf([char(repmat(32, 1, l)) '\tanalytic   diffs     delta\n']);
    for i = 1:length(names)
      spaceLen = l - length(names{i});
      space = char(repmat(32, 1, spaceLen));
      fprintf([space names{i} ':\t%4.6f\t%4.6f\t%4.6f\n'], ...
              g(i), gLDiff(i), gLDiff(i) - g(i));
    end
  end
end
if exist([priorType 'PriorGradient'])
  origX = x;
  for i = 1:length(x);
    x = origX;
    x(i) = origX(i) + epsilon;
    Lplus(i) = priorLogProb(prior, x);
    x(i) = origX(i) - epsilon;
    Lminus(i) = priorLogProb(prior, x);
  end
  x = origX;
  gLDiff = .5*(Lplus - Lminus)/epsilon;
  g = priorGradient(prior, x);

  paramMaxDiff = max(max(abs(gLDiff-g)));
%  if paramMaxDiff > 2*epsilon

    fprintf([ '\tanalytic   diffs     delta\n']);
    for i = 1:length(x)
      fprintf([num2str(i) ':\t%4.6f\t%4.6f\t%4.6f\n'], ...
              g(i), gLDiff(i), gLDiff(i) - g(i));
    end
%  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


% -------------------------------------------------------------------------
%   Description:
%       Demo script to calculate the BT scores
%       This script reproduces the results of Figure 5 in our paper.
%
%   Citation:
%       A Comparative Study for Single Image Blind Deblurring
%       Wei-Sheng Lai, Jia-Bin Huang, Zhe Hu, Narendra Ahuja, and Ming-Hsuan Yang
%       IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2016
%
%   Contact:
%       Wei-Sheng Lai
%       wlai24@ucmerced.edu
%       University of California, Merced
% -------------------------------------------------------------------------

%% input dataset and attribute
dataset = 'real';
% dataset = 'uniform';
% dataset = 'nonuniform';

attribute = 'all';
% attribute = 'manmade';
% attribute = 'natural';
% attribute = 'people';
% attribute = 'saturated';
% attribute = 'text';

%% Load list
list_filename = fullfile('list', 'method.txt');
method = load_list(list_filename);

attr_filename = fullfile('attributes', sprintf('%s_%s.txt', dataset, attribute));
fprintf('Load %s\n', attr_filename);
attr_list = dlmread(attr_filename);

num_img = length(attr_list);
num_method = 14;

%% load votes
vote_filename = fullfile('votes', sprintf('votes_%s_balance_%s.csv', dataset, attribute));
fprintf('Load %s\n', vote_filename);
M = csvread(vote_filename, 1, 0); % offset the first row to skip header

%% build winning matrix
score_matrix = zeros(num_img, num_method);

for i = 1:num_img

    row = find(M(:, 1) == attr_list(i));

    % convert M to winning matrix
    C = construct_winning_matrix(M(row, :), num_method);

    % compute BT scores
    score = BT_EM_exp(C);

    % avoid Nan and Inf when compute mean
    s = score;
    s = s(s == s); % remove NaN
    s = s(s ~= Inf); % remove Inf

    % normalize to zero mean
    score = score - mean(s);
    score_matrix(i, :) = score';

end

%% compute mean scores for each method
method_score = zeros(num_method, 1);
for i = 1:num_method
    s = score_matrix(:, i);
    s = s(s == s); % remove NaN
    s = s(s ~= Inf); % remove Inf
    method_score(i) = mean(s);
end

%% sort methods by mean scores
[~, order] = sort(method_score, 'descend');

%% plot cumulative frequency
[color, marker, line_style] = color_spec;

figure; hold on;
select_method = order(1:14);
legend_cmd = 'legend(';

x = -5:0.5:10;
for i = 1:length(select_method)
    m = select_method(i);
    s = score_matrix(:, m);
    s = s(s == s); % remove NaN
    s = s(s ~= Inf);  % remove Inf

    y = histc(s, x);
    z = cumsum(y)/sum(y);
    h = plot(x, z, 'LineWidth', 2);
    h.Color = color{m};
    h.LineStyle = line_style{m};

    legend_cmd = sprintf('%s method{%d}', legend_cmd, select_method(i));
    if( i < length(select_method) )
        legend_cmd = sprintf('%s, ', legend_cmd);
    end
end
legend_cmd = sprintf('%s);', legend_cmd);
l = eval(legend_cmd);
l.FontSize = 16;
l.Location = 'southeast';

if( strcmp(attribute, 'manmade') )
    title(sprintf('%s (man-made)', dataset));
else
    title(sprintf('%s (%s)', dataset, attribute));
end
xlabel('B-T Scores');
ylabel('Cumulative Frequency');

h = gca;
h.FontName = 'Times New Roman';
h.FontSize = 24;

hold off;

%% save results

% filename = sprintf('bt_cdf_%s_%s.pdf', dataset, attribute);
% saveas(h, filename);
% fprintf('Save %s\n', filename);


text:
function [connections]=findConnectionMatrix2(mesh)
% Finds the sparse connection matrix for a mesh
%
%  connections=findConnectionMatrix2(mesh)
%
% The mesh slots must be:
%     .vertices  (nPoints,3): 3D coords of all mesh vertices
%     .triangles (nFaces,3) : triplets of indices specifies one corner of a face triangle
%
% The returned connections is a sparse matrix (nPoints*nPoints)
%
% AUTHOR: Wade
% DATE : Last modified 020701

nVerts    = meshGet(mesh,'nVertices');
triangles = meshGet(mesh,'triangles')';

% faceIndexList defines vertices that are connected to each other:
% If two vertices are in the same row, they're in the same triangle and
% therefore are connected

% The values here are going to be used to create scalars
% (connection_locations) that index into the connections matrix.
% The multiplication by nVerts makes it possible to add together
% a row and column value (see below).
%
r1=(triangles(:,1))*nVerts;
r2=(triangles(:,2));
r3=(triangles(:,3));
r4=(triangles(:,2))*nVerts;

% Calculate the connection_locations
% These are connections between (1,2), (1,3), and (2,3)
%
connect_locations=[ (r1+r2+1) ;(r1+r3+1); (r4+r3+1) ];

[cly,clx]=ind2sub([nVerts,nVerts],connect_locations);
connections = sparse(cly,clx,1,nVerts,nVerts);

%Now make it symmetric because if a is connected to b then b is connected to a.
connections = connections+(connections');
connections = (connections~=0);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function f = cumsum(f, dim)
%CUMSUM   Indefinite integral of a SPHEREFUN.
%
%  F = CUMSUM(F) or F = CUMSUM(F, DIM) is not defined on a spherefun.
%
%   This is not allowed and returns an error.  This function exists so that the
%   error message is meaningful to a SPHEREFUN user.
%
% See also SPHEREFUN/CUMSUM2.

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

error('CHEBFUN:SPHEREFUN:CUMSUM:fail',...
      'The indefinite integral is not defined for a spherefun.');

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
classdef MeshMerger < handle

    properties (Access = private)
        nodesA
        nodesB
        allNodes
        mergedNodes
        mergedMesh
    end

    properties (Access = private)
        meshA
        meshB
        isMergedNodeA
        isMergedNodeB
    end

    methods (Access = public)

        function obj = MeshMerger(cParams)
            obj.init(cParams)
            obj.computeNodesA();
            obj.computeNodesB();
            obj.computeAllNodes();
        end

        function m = compute(obj)
            obj.computeMergedNodes();
            obj.computeMergedMesh();
            m = obj.computeCanonicalMergedMesh();
        end

        function rNodes = computeRemainingNodes(obj)
            aNodes  = obj.allNodes;
            remainingA = true(obj.meshA.nnodes,1);
            remainingB = ~obj.isMergedNodeB;
            remaining  = [remainingA,remainingB];
            rNodes   = aNodes(remaining);
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.meshA = cParams.meshA;
            obj.meshB = cParams.meshB;
            obj.isMergedNodeA = cParams.isMergedNodeA;
            obj.isMergedNodeB = cParams.isMergedNodeB;
        end

        function computeNodesA(obj)
            npnod = obj.meshA.nnodes;
            nodes = (1:npnod)';
            obj.nodesA = nodes;
        end

        function computeNodesB(obj)
            npnodA = obj.meshA.nnodes;
            npnodB = obj.meshB.nnodes;
            nodes = (1:npnodB)' + npnodA;
            obj.nodesB = nodes;
        end

        function computeAllNodes(obj)
            nA = obj.nodesA;
            nB = obj.nodesB;
            nodes = [nA;nB];
            obj.allNodes = nodes;
        end

        function computeMergedNodes(obj)
            nA   = obj.nodesA;
            nB   = obj.nodesB;
            isMa = obj.isMergedNodeA;
            isMb = obj.isMergedNodeB;
            mNodesA          = nA;
            mNodesB(isMb,1)  = nA(isMa);
            mNodesB(~isMb,1) = nB(~isMb);
            nodes = [mNodesA;mNodesB];
            obj.mergedNodes = nodes;
        end

       function computeMergedMesh(obj)
            s.connec = obj.computeMergedConnec();
            s.coord  = obj.computeMergedCoord();
            m = Mesh(s);
            obj.mergedMesh = m;
        end

        function mConnec = computeMergedConnec(obj)
            aConnec = obj.computeAllConnec();
            s.oldNodes = obj.allNodes;
            s.newNodes = obj.mergedNodes;
            c = ConnecRenumbering(s);
            mConnec = c.renumber(aConnec);
        end

        function connec = computeAllConnec(obj)
            npnodA  = obj.meshA.nnodes;
            connecA = obj.meshA.connec;
            connecB = obj.meshB.connec + npnodA;
            connec = [connecA;connecB];
        end

        function mCoord = computeMergedCoord(obj)
            coordA = obj.meshA.coord;
            coordB = obj.meshB.coord;
            mCoord = [coordA;coordB];
        end

        function m = computeCanonicalMergedMesh(obj)
            s.remainingNodes = obj.computeRemainingNodes();
            s.mesh           = obj.mergedMesh;
            c = CannonicalMeshComputer(s);
            m = c.compute();
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%% Creating one CUDA accelerated test radiography

%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% This file is part of the TIGRE Toolbox
%
% Copyright (c) 2015, University of Bath and
%                     CERN-European Organization for Nuclear Research
%                     All rights reserved.
%
% License:            Open Source under BSD.
%                     See the full license at
%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
%
% Contact:            tigre.toolbox@gmail.com
% Codes:              https://github.com/CERN/TIGRE/
% Coded by:           Stefanie Kaser
%--------------------------------------------------------------------------
%% pCT geometry definition

geo.dDetector = [0.25; 0.25]; % Pixel size of optimized pRad in mm.
geo.DSID = 245; % Distance between source and upstream detector.
geo.DSO = 300; % Distance between source and origin (center of rotation).
geo.DSD = 355; % Distance between source and downstream detector.
geo.hull = [0; 0; 0; 0]; % We won't use a convex hull here (all entries are
% set to 0).
geo.sDetector = [20; 20]; % Defines the size (in mm) of the optimized pRad.
geo.mode = 'parallel'; % or 'cone' in case of cone beam geometry

%% pCT initial beam energy
eIn = single(100);

%% Load pCT data set.
% The data set corresponds to one radiographic image (pRad) of an Aluminum
% stair profile as used in https://arxiv.org/abs/2106.12890. However, our
% test data set is based on Monte Carlo simulations (GATE:
% doi.org/10.1016/S0920-5632(03)90969-8). The data set contains the
% protons' upstream and downstream positions and directions as well as the
% single proton WEPLs (all in mm).
data = pCTdata();

%% Binning the data into an optimized pRad
% Finally, the single proton data are binned into an optimized pRad.
proj = pCTCubicSpline_mex(data.posIn, data.posOut, data.dirIn, ...
    data.dirOut, data.Wepl, eIn, geo);
% We are only creating one test pRad here. Creating pRads at multiple
% rotation angles would allow to use the collection as input for TIGRE's
% reconstruction algorithms (see d04_SimpleReconstruction.m).

%% Plot the result
imshow(proj, [0, 3], 'InitialMagnification', 500);

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
% This example shows how to use the stretchmesh function to expand the
% size of the computation window at the edges without increasing the
% number of gridpoints.  It uses the same example waveguide
% structure considered in 'basic-semivector.m'

% Refractive indices:
n1 = 3.34;          % Lower cladding
n2 = 3.44;          % Core
n3 = 1.00;          % Upper cladding (air)

% Vertical dimensions:
h1 = 2.0;           % Lower cladding
h2 = 1.3;           % Core thickness
h3 = 0.5;           % Upper cladding
rh = 1.1;           % Ridge height

% Horizontal dimensions:
rw = 1.0;           % Ridge half-width
side = 1.5;         % Space on side

% Grid size:
dx = 0.0125;        % grid size (horizontal)
dy = 0.0125;        % grid size (vertical)

lambda = 1.55;      % vacuum wavelength
nmodes = 1;         % number of modes to compute

fprintf (1,'generating index mesh...\n');
[x,y,xc,yc,nx,ny,eps,edges] = waveguidemesh([n1,n2,n3],[h1,h2,h3], ...
                                            rh,rw,side,dx,dy);

% Stretch the 40 gridpoints on the south and east edges of the
% computational window by factors of 4 and 3, respectively.

[x,y,xc,yc,dx,dy] = stretchmesh(x,y,[0,40,40,0],[1,4,3,1]);

[Ex,neff] = svmodes(lambda,n2,nmodes,dx,dy,eps,'000S','EX');

fprintf(1,'neff = %.6f\n',neff);

figure(1);
contourmode(x,y,Ex);
title('Ex (TE Mode)'); xlabel('x'); ylabel('y');
for v = edges, line(v{:}); end


text:
function a = softMaxSelection(Q, state, actions, T)

    n_actions = length(actions);
    P = zeros(1, n_actions);

    random = rand();

    % Boltzmann distribution
    P = exp(Q(state,:)./T)./sum(exp(Q(state,:)./T));

    if( find(isnan(P),1) ~= 0 )

       a = getBestAction(Q, state, actions);
    else

        if(random < P(1))

            a = 1;
        elseif(random >= P(1) && random < sum(P(1:2)) )

            a = 2;
        elseif(random >= sum(P(1:2)))

            a = 3;
        end
    end

end

author: kennydl

repo: Reinforcment-Learning-With-Q-Learning


text:
classdef mp_table
%MP_TABLE  Very basic TABLE class for Octave or older Matlab.

%   MATPOWER
%   Copyright (c) 2021, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

    properties (Access=protected)
        Properties = struct('VariableNames',  {{}}, ...
                            'VariableValues', {{}}, ...
                            'RowNames',       {{}}, ...
                            'DimensionNames', {{}}    );
    end     %% properties

    methods
        function obj = mp_table(varargin)
            args = varargin;

            if nargin
                %% extract named arguments
                [var_names, row_names, dim_names, args] = ...
                    extract_named_args(obj, args);

                %% set default variable names
                nv = length(args);          %% number of variables
                if length(var_names) < nv
                    for k = nv:-1:length(var_names)+1
                        var_names{k} = inputname(k);
                        if isempty(var_names{k})
                            var_names{k} = sprintf('Var%d', k);
                        end
                    end
                end

                %% set default dimension names
                if isempty(dim_names)
                    dim_names = {'Row', 'Variables'};
                end

                %% check dimensions
                if length(var_names) > nv
                    warning('mp_table: ignoring %d extra VariableNames provided', length(var_names)-nv);
                end
                nr = size(args{1}, 1);      %% number of rows
                for k = 2:nv
                    if size(args{k}, 1) ~= nr
                        error('mp_table: variable %d expected to have %d rows', k, nr);
                    end
                end

                %% set table properties
                obj.Properties.VariableNames = var_names;
                obj.Properties.DimensionNames = dim_names;
                obj.Properties.RowNames = row_names;
                obj.Properties.VariableValues = args;
                obj.Properties.Map = cell2struct(num2cell(1:nv), var_names, 2);
            end
        end

        function new_obj = copy(obj)
            new_obj = eval(class(obj));  %% create new object

            %% make copies of properties
            if have_feature('octave')
                s1 = warning('query', 'Octave:classdef-to-struct');
                warning('off', 'Octave:classdef-to-struct');
            end
            props = fieldnames(obj);
            if have_feature('octave')
                warning(s1.state, 'Octave:classdef-to-struct');
            end
            for k = 1:length(props)
                new_obj.(props{k}) = obj.(props{k});
            end

            %% make copies of values
            for k = 1:length(obj.Properties.VariableValues)
                if isobject(obj.Properties.VariableValues{k}) && ...
                        ismethod(obj.Properties.VariableValues{k}, 'copy')
                    new_obj.Properties.VariableValues{k} = ...
                        copy(obj.Properties.VariableValues{k});
                end
            end
        end

        %% not really useable until Octave implements istable()
        function TorF = istable(obj)
            TorF = true;
        end

        function varargout = size(obj, dim)
            varargout = cell(1, nargout);
            w = length(obj.Properties.VariableValues);
            if w
                h = size(obj.Properties.VariableValues{1}, 1);
            else
                h = 0;
            end
            if nargin == 2
                if dim == 1
                  varargout{1} = h;
                elseif dim == 2
                  varargout{1} = w;
                else
                  varargout{1} = 1;
                end
            elseif nargout == 0 || nargout == 1
                varargout{1} = [h, w];
            else
                varargout{1} = h;
                varargout{2} = w;
                [varargout{3:end}] = deal(1);
            end
        end

        function TorF = isempty(obj)
            TorF = prod(size(obj)) == 0;
        end

        function N = end(obj, k, n)
            N = size(obj, k);
        end

        function n = numArgumentsFromSubscript(obj, s, indexingContext)
            n = 1;
        end

        function n = numel(obj, varargin)
            %% needed to avoid error when doing T{r1:rN, c} = b;
            %% see https://github.com/apjanke/octave-tablicious/issues/80#issuecomment-855198281
            n = 1;
        end

        function b = subsref(obj, s)
            switch s(1).type
                case '.'
                    if strcmp(s(1).subs, 'Properties')
                        b = obj.Properties;
                    else
                        p = obj.Properties;
                        b = p.VariableValues{p.Map.(s(1).subs)};
                    end
                case '()'
                    r = s(1).subs{1};
                    if length(s(1).subs) > 1
                        c = s(1).subs{2};
                    else
                        c = ':';
                    end
                    var_names = obj.Properties.VariableNames(1,c);
                    %% Octave has issues with this ...
                    % var_vals = cellfun(@(x)x(r,:), ...
                    %     obj.Properties.VariableValues(c), 'UniformOutput', 0);
                    %% ... so we use explicit subsref instead ...
                    ss = struct('type', '()', 'subs', {{r, ':'}});
                    var_vals = cellfun(@(x)subsref(x, ss), ...
                        obj.Properties.VariableValues(c), 'UniformOutput', 0);
                    dim_names = obj.Properties.DimensionNames;
                    if isempty(obj.Properties.RowNames)
                        row_names = {};
                    else
                        row_names = obj.Properties.RowNames(r);
                    end
                    b = mp_table(var_vals{:}, ...
                        'VariableNames', var_names, 'RowNames', row_names, ...
                        'DimensionNames', dim_names);
                case '{}'
                    a = horzcat(obj.Properties.VariableValues{s(1).subs{2}});
                    b = a(s(1).subs{1}, :);
            end
            if length(s) > 1    %% recurse
                b = subsref(b, s(2:end));
            end
        end

        function obj = subsasgn(obj, s, b)
            R = length(s) > 1;      %% need to recurse
            switch s(1).type
                case '.'
                    if strcmp(s(1).subs, 'Properties')
                        if R
                            obj.Properties = ...
                                subsasgn(obj.Properties, s(2:end), b);
                        else
                            obj.Properties = b;
                        end
                    else
                        if isfield(obj.Properties.Map, s(1).subs)
                            idx = obj.Properties.Map.(s(1).subs);
                            if R
                                obj.Properties.VariableValues{idx} = ...
                                    subsasgn(obj.Properties.VariableValues{idx}, ...
                                        s(2:end), b);
                            else
                                obj.Properties.VariableValues{idx} = b;
                            end
                        else    %% new variable
                            if R
                                error('mp_table/subsagn: cannot index new variable ''%s''', s(1).subs);
                            end
                            %% check size compatibility
                            if size(b, 1) ~= size(obj, 1)
                                error('mp_table/subsagn: number of rows (%d) in new variable ''%s'' does not match height of table (%d)', size(b, 1), s(1).subs, size(obj, 1));
                            end
                            %% add variable
                            idx = size(obj, 2) + 1;
                            obj.Properties.VariableNames{idx} = s(1).subs;
                            obj.Properties.Map.(s(1).subs) = idx;
                            obj.Properties.VariableValues{idx} = b;
                        end
                    end
                case '()'
                    nv = size(b, 2);
                    idx = s(1).subs{2};
                    if ischar(idx)
                        if length(idx) == 1 && idx == ':'
                            idx = [1:nv];
                        else        %% for 'a:b' syntax as a char array
                            idx = eval(['[' idx ']'])
                        end
                    end
                    if nv
                        ss = struct('type', '{}', 'subs', {{':', 0}});
                        for k = 1:nv
                            ss(1).subs{2} = k;
                            obj.Properties.VariableValues{idx(k)}(s(1).subs{1}, :) = ...
                                subsref(b, ss);     %% = b{:, k};
                        end
                    else    %% delete variable
                        obj.Properties.VariableNames(idx) = [];
                        obj.Properties.VariableValues(idx) = [];
                        %% rebuild map
                        nv = size(obj, 2);
                        obj.Properties.Map = cell2struct(num2cell(1:nv), obj.Properties.VariableNames, 2);
                    end
                    if R
                        error('mp_table/subasgn: sub-indexing following {}-indexing not supported');
                    end
                case '{}'
                    idx = s(1).subs{2};
                    if ischar(idx)
                        if length(idx) == 1 && idx == ':'
                            idx = [1:size(obj, 2)];
                        else        %% for 'a:b' syntax as a char array
                            idx = eval(['[' idx ']'])
                        end
                    end
                    j = 1;  %% col index in b
                    for k = 1:length(idx)
                        %% get num cols for var k
                        nc = size(obj.Properties.VariableValues{idx(k)}, 2);
                        if nc > 1
                            obj.Properties.VariableValues{idx(k)}(s(1).subs{1}, :) = b(:, j:j+nc-1);
                        else
                            obj.Properties.VariableValues{idx(k)}(s(1).subs{1}, :) = b(:, j);
                        end
                        j = j+nc;
                    end
                    if R
                        error('mp_table/subasgn: sub-indexing following ()-indexing not supported');
                    end
            end
        end

        function obj = horzcat(obj, varargin)
            [nr, nv] = size(obj);
            nvn = length(obj.Properties.VariableNames); %% num variable names
            nrn = length(obj.Properties.RowNames);      %% num row names

            %% for each table to append
            for j = 1:length(varargin)
                Tj = varargin{j};       %% table being appended
                [nrj, nvj] = size(Tj);
                nrnj = length(Tj.Properties.RowNames);      %% num row names Tj

                %% check for consistent number of rows
                if nrj ~= nr
                    error('mp_table/horzcat: number of rows in arg %d (%d) does not match that in arg 1 (%d)', ...
                        j+1, nrj, nr);
                end

                %% check for row name mismatch
                if nrn          %% obj has row names
                    if nrnj && ~isequal(obj.Properties.RowNames, ...
                                         Tj.Properties.RowNames)
                        for k = 1:nrn
                            if ~strcmp(obj.Properties.RowNames{k}, ...
                                        Tj.Properties.RowNames{k})
                                error('mp_table/horzcat: name of row %d in arg %d (%s) does not match the one in arg 1 (%s)', ...
                                    k, j+1, Tj.Properties.RowNames{k}, ...
                                        obj.Properties.RowNames{k});
                            end
                        end
                    end
                elseif nrnj     %% Tj has row names, but obj does not (use them)
                    obj.Properties.RowNames = Tj.Properties.RowNames;
                    nrn = nrnj;
                end
                var_names = Tj.Properties.VariableNames;

                %% stack variable names and check for duplicates
                obj.Properties.VariableNames = ...
                    [obj.Properties.VariableNames, var_names];
                nvn = nvn + nvj;    %% update num variable names
                if length(unique(obj.Properties.VariableNames)) < nvn
                    [~, ia, ~] = unique(obj.Properties.VariableNames);
                    k = find(ismember(1:nvn, ia));
                    error('mp_table/horzcat: duplicate variable name ''%s''', ...
                        obj.Properties.VariableNames{k(1)});
                end

                %% horzcat var values
                obj.Properties.VariableValues = ...
                    [ obj.Properties.VariableValues, ...
                       Tj.Properties.VariableValues ];
                nv = length(obj.Properties.VariableValues);

                %% update Properties.Map
                obj.Properties.Map = cell2struct(num2cell(1:nv), ...
                    obj.Properties.VariableNames, 2);
            end
        end

        function obj = vertcat(obj, varargin)
            [nr, nv] = size(obj);
            nvn = length(obj.Properties.VariableNames); %% num variable names
            nrn = length(obj.Properties.RowNames);      %% num row names

            %% for each table to append
            for j = 1:length(varargin)
                Tj = varargin{j};       %% table being appended
                [nrj, nvj] = size(Tj);
                nvnj = length(Tj.Properties.VariableNames); %% num var names Tj

                %% check for consistent number of variables
                if nvj ~= nv
                    error('mp_table/vertcat: number of variables in arg %d (%d) does not match that in arg 1 (%d)', ...
                        j+1, nvj, nv);
                end

                %% check for variable name mismatch
                if ~isequal(obj.Properties.VariableNames, ...
                             Tj.Properties.VariableNames)
                    for k = 1:nvn
                        if ~strcmp(obj.Properties.VariableNames{k}, ...
                                    Tj.Properties.VariableNames{k})
                            error('mp_table/vertcat: name of variable %d in arg %d (%s) does not match the one in arg 1 (%s)', ...
                                k, j+1, Tj.Properties.VariableNames{k}, ...
                                    obj.Properties.VariableNames{k});
                        end
                    end
                end

                %% create missing row names if any exist
                row_names = Tj.Properties.RowNames;
                if nrn                      %% obj has row names, Tj does not
                    if length(row_names) == 0
                        %% create names for rows of Tj
                        row_names = cell(nrj, 1);
                        for kk = 1:nrj
                            row_names{kk} = sprintf('Row%d', kk+nrn);
                        end
                    end
                elseif length(row_names)    %% Tj has row names, obj does not
                    %% create names for rows of obj
                    obj.Properties.RowNames = cell(nr, 1);
                    for kk = 1:nr
                        obj.Properties.RowNames{kk} = sprintf('Row%d', kk);
                    end
                end
                nrn = length(obj.Properties.RowNames);

                %% stack row names and check for duplicates
                if nrn
                    obj.Properties.RowNames = ...
                        [obj.Properties.RowNames; row_names];
                    nrn = nrn + nrj;    %% update num row names
                    if length(unique(obj.Properties.RowNames)) < nrn
                        [~, ia, ~] = unique(obj.Properties.RowNames);
                        k = find(ismember(1:nrn, ia));
                        error('mp_table/horzcat: duplicate row name ''%s''', ...
                            obj.Properties.RowNames{k(1)});
                    end
                end

                %% for each variable
                for k = 1:nv
                    %% vertcat var values
                    obj.Properties.VariableValues{k} = ...
                        [ obj.Properties.VariableValues{k}; ...
                           Tj.Properties.VariableValues{k} ];
                end
            end
        end

        function [var_names, row_names, dim_names, args] = ...
                extract_named_args(obj, args)
            var_names = {};
            row_names = {};
            dim_names = {};
            for k = length(args)-1:-1:1
                if ischar(args{k})
                    switch args{k}
                        case 'VariableNames'
                            var_names = args{k+1};
                            args(k:k+1) = [];
                        case 'RowNames'
                            row_names = args{k+1};
                            args(k:k+1) = [];
                        case 'DimensionNames'
                            dim_names = args{k+1};
                            args(k:k+1) = [];
                    end
                end
            end
        end

        function display(obj)
            allrows = 1;

            [nr, nv] = size(obj);

            %% rows to display
            max_nr = 25;
            first_last = 10;
            if nr <= max_nr || allrows
                rr = [1:nr]';
            else
                rr = [1:first_last nr-first_last+1:nr]';
            end

            %% var names
            vn = obj.Properties.VariableNames;
            rn = obj.Properties.RowNames;

            %% col_widths
            if isempty(rn)
                cwr = 0;
            else
                cwr = max(cellfun(@length, rn));    %% col width for row names
            end
            for v = nv:-1:1
                cw(v) = length(vn{v});  %% actual column width
                cwv(v) = 0;             %% column width for displayed values
                vv = obj.Properties.VariableValues{v};
                if size(vv, 2) > 1 || isa(vv, 'mp_table') || isa(vv, 'table')
                    s = sprintf('1x%d %s' , size(vv, 2), class(vv));
                    len_s = length(s);
                    if len_s > cw(v);
                        cw(v) = len_s;
                    end
                    if len_s > cwv(v);
                        cwv(v) = len_s;
                    end
                else
                    if iscell(vv)
                        for r = 1:length(rr)
                            s = sprintf('''%s''', vv{rr(r)});
                            len_s = length(s);
                            if len_s > cw(v);
                                cw(v) = len_s;
                            end
                            if len_s > cwv(v);
                                cwv(v) = len_s;
                            end
                        end
                    else
                        for r = 1:length(rr)
                            if isempty(vv)
                                s = '[]';
                            else
                                s = sprintf('%g', vv(rr(r), 1));
                            end
                            len_s = length(s);
                            if len_s > cw(v);
                                cw(v) = len_s;
                            end
                            if len_s > cwv(v);
                                cwv(v) = len_s;
                            end
                        end
                    end
                end
                spc{v} = repmat(' ', 1, floor((cw(v)-cwv(v))/2));
            end

            fprintf('\n %dx%d %s\n\n', nr, nv, class(obj));

            %% display var names & lines
            fprintf('  %s', repmat(' ', 1, cwr));
            for v = 1:nv
                fmt = sprintf('  %%-%ds', cw(v));
                name = sprintf('%s', vn{v});
                spcs = repmat(' ', 1, floor((cw(v)-length(name))/2));
                fprintf(fmt, [spcs name]);
            end
            fprintf('\n');
            fprintf('  %s', repmat(' ', 1, cwr));
            for v = 1:nv
                s = repmat('-', 1, cw(v));
                fprintf('  %s', s);
            end
            fprintf('\n');

            %% display var values
            for r = 1:length(rr)
                %% row names
                if isempty(rn)
                    fprintf('  ');
                else
                    fmt = sprintf('  %%-%ds', cwr);
                    fprintf(fmt, rn{rr(r)});
                end

                %% variable values
                for v = 1:nv
                    vv = obj.Properties.VariableValues{v};
                    if size(vv, 2) > 1 || isa(vv, 'mp_table') || isa(vv, 'table')
                        val = sprintf('%s1x%d %s', spc{v}, size(vv, 2), class(vv));
                        fmt = sprintf('  %%-%ds', cw(v));
                    else
                        if iscell(vv)
                            val = sprintf('%s''%s''', spc{v}, vv{rr(r), 1});
                            fmt = sprintf('  %%-%ds', cw(v));
                        else
                            if isempty(vv)
                                val = sprintf('[]%s', spc{v});
                            else
                                val = sprintf('%g%s', vv(rr(r), 1), spc{v});
                            end
                            fmt = sprintf('  %%%ds', cw(v));
                        end
                    end
                    fprintf(fmt, val);
                end
                fprintf('\n');

                %% print continuation dots if not displaying all rows
                if nr > length(rr) && r == first_last
                    for k = 1:3
                        for v = 1:nv
                            spcs = repmat(' ', 1, floor(cwv(v)/2));
                            if iscell(vv)
                                dot = sprintf('%s%s.', spc{v}, spcs);
                                fmt = sprintf('  %%-%ds', cw(v));
                            else
                                dot = sprintf('.%s%s', spcs, spc{v});
                                fmt = sprintf('  %%%ds', cw(v));
                            end
                            fprintf(fmt, dot);
                        end
                        fprintf('\n');
                    end
                end
            end
            fprintf('\n');
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [eta,etaf]=v_finishat(frac,tol,fmt)
%V_FINISHAT print estimated finish time of a long computation (FRAC,TOL,FMT)
% Usage:  (1)  for i=1:many
%                  v_finishat((i-1)/many);  % initializes on first pass when i=1
%                  ... computation ...
%              end
%
%         (2)  for i=1:many
%                  v_finishat([i-1 many]);  % alternative argument format
%                  ... computation ...
%              end
%
%         (3)  v_finishat(0);               % explicit initialization before loop
%              for i=1:many
%                  ... computation ...
%                  v_finishat(i/many);      % calculate fraction completed
%              end
%
%         (4)  for i=1:NI
%                  for j=1:NJ
%                      for k=1:NK
%                          v_finishat([i NI; j NJ; k-1 NK]); % one row per nested loop
%                          ... computation ...
%                      end
%                  end
%              end
%
% Inputs: FRAC = fraction of total comutation that has been completed
%                Alternatively at start of inner loop: [i NI; j NJ; k-1 NK ...] where i, j, k are
%                loop indices and NI, NJ, NK their limits. Use k instead of k-1 if placed at
%                the end of the inner loop. As a special case, FRAC=0 initializes the routine.
%         TOL  = Tolerance in minutes. If the estimated time has changed by less
%                than this, then nothing will be printed. [default 10% of remaining time]
%         FMT  = Format string which should include %s for estimated finish time, %d for remaining minutes and %f for fraction complete
%
% Output: ETA  = string containing the expected finish time
%                specifying this will suppress printing message to std err (fid=2)
%         ETAF = expected finish time as a daynumber
%
% Example:       v_finishat(0);
%                for i=1:many
%                    long computation;
%                    v_finishat(i/many);
%                end

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

persistent oldt oldnw
if nargin<3
    fmt='Estimated finish at %s (%.2f done, %d min remaining)\n';

end
nf=size(frac,1);
if all(frac(:,1)<=[ones(nf-1,1); 0]) % initialize if fraction done is <=0
    oldt=0;
    eta='Unknown';
    tic;
else
    if size(frac,2)==2
        fp=cumprod(frac(:,2));
        frac=sum((frac(:,1)-1)./fp)+1/fp(end);
    end
    nw=now;                             % current time as serial number
    sectogo=(1/frac-1)*toc;      % seconds to go
    newt=nw+sectogo/86400;       % add estimated time in days
    if nargin<2 || ~numel(tol)
        tol=max(0.1*(newt-nw)*1440,1);
    end
    if ~exist('oldt','var') || oldt==0 || (abs(newt-oldt)>tol/1440 && (nw-oldnw)>10/86400) || (nw-oldnw)>10/1440 || nargout>0
        oldt=newt;
        if floor(oldt)==floor(nw)
            df='HH:MM';
        else
            df='HH:MM dd-mmm-yyyy';
        end
        eta=datestr(oldt,df);
        if ~nargout
            ix=find(fmt=='%',1);
            while ~isempty(ix)
                fprintf(2,fmt(1:ix-1));
                fmt=fmt(ix:end);
                ix=find(fmt>='a' & fmt<='z',1); % find letter
                switch fmt(ix)
                    case 's'
                        fprintf(2,fmt(1:ix),eta);
                    case 'd'
                        fprintf(2,fmt(1:ix),round(sectogo/60));
                    case 'f'
                        fprintf(2,fmt(1:ix),frac);
                end
                fmt=fmt(ix+1:end);
                ix=find(fmt=='%',1);
            end
            fprintf(2,fmt);
        end
        oldnw=nw;                           %
    end
end
etaf=oldt;

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function [ f, g, H ] = opt10_fgh ( x, flag )

%*****************************************************************************80
%
%% OPT10_FGH evaluates F, G and H for test case #10.
%
%  Discussion:
%
%    This example is known as the Wood function.
%
%    The optimizing value is
%
%      X* = (1,1,1,1)
%
%    for which
%
%      F(X*) = 0.
%
%  Modified:
%
%    09 January 2008
%
%  Author:
%
%    Jeff Borggaard,
%    Gene Cliff,
%    Virginia Tech.
%
%  Reference:
%
%    John Dennis, Robert Schnabel,
%    Numerical Methods for Unconstrained Optimization
%    and Nonlinear Equations,
%    SIAM, 1996,
%    ISBN13: 978-0-898713-64-0,
%    LC: QA402.5.D44.
%
%  Parameters:
%
%    Input, real X(4), the evaluation point.
%
%    Input, string FLAG, indicates what must be computed.
%    'f' means only the value of F is needed,
%    'g' means only the value of G is needed,
%    'all' means F, G and H (if appropriate) are needed.
%    It is acceptable to behave as though FLAG was 'all'
%    on every call.
%
%    Output, real F, the optimization function.
%
%    Output, real G(4,1), the gradient column vector.
%
%    Output, real H(4,4), the Hessian matrix.
%
  n = length ( x );

  if ( n ~= 4 )
    fprintf ( '\n' );
    fprintf ( 'OPT10_FGH - Fatal error!\n' );
    fprintf ( '  The input vector X should have length 4.\n'),
    fprintf ( '  Instead, it has length = %d.\n', n );
    keyboard
  end

  f = 100   * ( x(1)^2 - x(2)   )^2 ...
    +         ( 1 - x(1)        )^2 ...
    +  90   * ( x(3)^2 - x(4)   )^2 ...
    +         ( 1 - x(3)        )^2 ...
    +  10   * ( x(2) + x(4) - 2 )^2 ...
    +   0.1 * ( x(2) - x(4)     )^2;

  g = zeros(n,1);

  g(1,1) =   100   * 2 * ( x(1)^2 - x(2) ) * 2 * x(1)...
         -           2 * ( 1 - x(1) );

  g(2,1) = - 100   * 2 * ( x(1)^2 - x(2) )...
           +  10   * 2 * ( x(2) + x(4) - 2 )...
           +   0.1 * 2 * ( x(2) - x(4) );

  g(3,1) =    90   * 2 * ( x(3)^2 - x(4) ) * 2 * x(3)...
           -         2 * ( 1 - x(3) );

  g(4,1) = -  90   * 2 * ( x(3)^2 - x(4) )...
           +  10   * 2 * ( x(2) + x(4) - 2 )...
           -   0.1 * 2 * ( x(2) - x(4) );

  H = zeros(n,n);

  H(1,1) =   1200 * x(1)^2 - 400 * x(2) + 2;
  H(1,2) = - 400 * x(1);
  H(1,3) =     0;
  H(1,4) =     0;

  H(2,1) = - 400 * x(1);
  H(2,2) =   220.2;
  H(2,3) =     0;
  H(2,4) =    19.8;

  H(3,1) = H(1,3);
  H(3,2) = H(2,3);
  H(3,3) = 1080 * x(3)^2 - 360 * x(4) + 2;
  H(3,4) = - 360 * x(3);

  H(4,1) = H(1,4);
  H(4,2) = H(2,4);
  H(4,3) = H(3,4);
  H(4,4) = 200.2;


text:
% Concatenate the neighbouring frames vectors to higher dimensional vectors
% Input is Number of dimension x number of frames
function [y] = ExpandContext(x, context, window_type)
[dim nFr] = size(x);

if nargin<3
    window_type = 'null';
end

if context>1
    half_context = (context-1)/2;
    idx = [ones(1,half_context) 1:nFr ones(1,half_context)*nFr];
    x = x(:,idx);
    %x = [ repmat(x(:,1),1,half_context) x repmat(x(:,end), 1, half_context) ];
    if 0    % this implementation may be slow
        y = [];
        for i=1:context
            y = [y; x(:,i:end-context+i)];
        end
    else
        if strcmp(class(x), 'gpuArray')
            y = gpuArray.zeros(dim*context,nFr);
        else
            y = zeros(dim*context,nFr);    % allocate memory first
        end
        context = gather(context);
        for i=1:context
            y((i-1)*dim+1:i*dim,:) = x(:,i:end-context+i);
        end
    end
else
    y = x;
end

switch window_type
    case 'null'
        window = [];
    case 'Hamming'
        window = hamming(context);
        window = repmat(window, 1, dim);
        window = reshape(window', dim*context, 1);
        window = diag(window);
        y = window * y;
    otherwise
        fprintf('Unkown window type\n');
        return;
end


text:
function anat = makeFlatAnat(view)
%
% function anat = makeFlatAnat(view)
%
% djh, 8/1/99  Culled this stuff out of loadAnat.
% djh, 8/18/99 Modified to use curvature if available.

disp('Interpolating flat images...');

% Initialize to zeros
imSize = view.ui.imSize;
anat = NaN*ones([imSize,2]);

% Loop through hemispheres and make the images with myGriddata
for h=1:2

    % Load gLocs to get 2D coords and Z
    if h==1
        [gLocs2d,gLocs3d,curvature] = loadGLocs('left',view.leftPath);

    else
        [gLocs2d,gLocs3d,curvature] = loadGLocs('right',view.rightPath);

    end
    if ~isempty(gLocs2d)
        y = gLocs2d(1,:);
        x = gLocs2d(2,:);

        % If curvature is nonempty, then use it.  Otherwise use z coord

        if ~isempty(curvature)
            z = curvature;
        else
            z = gLocs3d(3,:);
        end

        % Remove NaN coords
        allFinite = isfinite(x) & isfinite(y);

        x = x(allFinite);
        y = y(allFinite);

        % Deal with new curvature coding: only l1 curvatures are now
        % encoded to make the flat maps look nicer. The following code
        % checks to see if we are building an anatomy from the old or the
        % new version of the flat files.
        if (length(z)<length(x))
            x=x(1:length(z));
            y=y(1:length(z));
        else
            z = z(allFinite);
        end

        % call griddata
        yi = [1:imSize(1)]';
        xi = [1:imSize(2)];
        warning off;

        if (checkML7)
            disp('Matlab 7 or higher detected..');
%             anat(:,:,h) = griddata(x,y,z,xi,yi,'linear',{'QJ'});
            anat(:,:,h) = griddata(x,y,z,xi,yi,'linear');
        else
            anat(:,:,h) = griddata(x,y,z,xi,yi,'linear');
        end

        warning backtrace;
    end
end


text:
% Generate MAR(2) and fit MAR model

disp('Generating data from known MAR(2) model');
d=2;
p=2;
T=100;
w=[0;0];

% Coeffs at lag 1
A1 = [ 0.4   1.2;   0.3   0.7 ];
% Coeffs at lag 2
A2 = [ 0.35 -0.3;  -0.4  -0.5 ];
A = [ A1 A2 ];

C = [ 1.00  0.50;   0.50  1.50 ];
lambda_true=inv(C);

%  Generate observations
x = spm_mar_gen (w, A, C, T);

logev=[];
for m=1:5,
    disp(sprintf('Fitting MAR model with %d components',m));
    mar=spm_mar(x,m);
    logev=[logev; mar.fm];
end
logev=logev-min(logev);

figure
subplot(2,1,1);
plot(x);
title('Bivariate time series from MAR(2) model');
subplot(2,1,2);
bar(logev);
xlabel('Number of time lags');
ylabel('Log Evidence');

% Specify prior - this is optional.
% spm_mar.m runs without the prior being set.
prior=spm_mar_prior(d,p,'global');
[mar,y,y_pred]=spm_mar(x,2,prior);

disp(' ');
disp('Estimates from fitting MAR(2) model');
disp('Lag 1');
disp('True coefficients');
disp(A1);
disp('Estimated coefficients');
disp(-mar.lag(1).a)

disp('Lag 2');
disp('True coefficients');
disp(A2);
disp('Estimated coefficients');
disp(-mar.lag(2).a)

disp('Noise covariance');
disp('True:');
disp(C);
disp('Estimated:');
disp(mar.noise_cov);

author: spm


text:
%% Housekeeping
clear
close all
clc
%%
tic
%% create data

run('tut01_data')

%% estimate models

run('tut02_estimation')

%% compare residuals

run('tut03_compare_residuals')

%% set identification

run('tut04_identification')

%% Compare structural shocks

run('tut04b_compare_structural_shocks')
%% simple choleski IRFs

run('tut05_compare_choleski_irfs')

%% more involved restrictions

run('tut06_compare_combos_irfs')

%% variance decomposition

run('tut07_variance_decomposition')

%% historical decomposition

run('tut08_historical_decomposition')

%% sample distribution of parameters

run('tut09_bootstrap')

%% distribution of variance decomposition

run('tut10_variance_decomposition_distribution')

%% distribution of historical decomposition

run('tut11_historical_decomposition_distribution')

%% distribution of IRFs

run('tut12_irf_distribution')

%% Bayesian estimation

run('tut13_bayestimation')

%% forecast

run('tut14_bayesforecast')

%% conditional forecast

run('tut15_conditionalforecast')

%% entropic tilting

%% clean up
toc

delete *.mat

author: jmaih


text:
classdef AvoidObstacles < simiam.controller.Controller

% Copyright (C) 2013, Georgia Tech Research Corporation
% see the LICENSE file included with this software

    properties

        % memory banks
        E_k
        e_k_1

        % gains
        Kp
        Ki
        Kd

        % plot support
        p

        % sensor geometry
        calibrated
        sensor_placement
    end

    properties (Constant)
        inputs = struct('v', 0);
        outputs = struct('v', 0, 'w', 0)
    end

    methods

        function obj = AvoidObstacles()
            obj = obj@simiam.controller.Controller('avoid_obstacles');
            obj.calibrated = false;

            obj.Kp = 5;
            obj.Ki = 0.01;
            obj.Kd = 0.01;

            obj.E_k = 0;
            obj.e_k_1 = 0;
        end

        function outputs = execute(obj, robot, state_estimate, inputs, dt)

            % Compute the placement of the sensors
            if(~obj.calibrated)
                obj.set_sensor_geometry(robot);
            end

            % Unpack state estimate
            [x, y, theta] = state_estimate.unpack();

            % Poll the current IR sensor values 1-5
            ir_distances = robot.get_ir_distances();

            % Interpret the IR sensor measurements geometrically
            ir_distances_wf = obj.apply_sensor_geometry(ir_distances, state_estimate);

            % 1. Compute the heading vector for obstacle avoidance

            sensor_gains = [1 1 0.5 1 1];
            u_i = (ir_distances_wf-repmat([x;y],1,5))*diag(sensor_gains);
            u_ao = sum(u_i,2);

            % 2. Compute the heading and error for the PID controller
            theta_ao = atan2(u_ao(2),u_ao(1));
            e_k = theta_ao-theta;
            e_k = atan2(sin(e_k),cos(e_k));

            e_P = e_k;
            e_I = obj.E_k + e_k*dt;
            e_D = (e_k-obj.e_k_1)/dt;

            % PID control on w
            v = inputs.v;
            w = obj.Kp*e_P + obj.Ki*e_I + obj.Kd*e_D;

            % Save errors for next time step
            obj.E_k = e_I;
            obj.e_k_1 = e_k;

            % plot
            obj.p.plot_2d_ref(dt, theta, theta_ao, 'g');

%             fprintf('(v,w) = (%0.4g,%0.4g)\n', v,w);

            outputs.v = v;
            outputs.w = w;
        end

        % Helper functions

        function ir_distances_wf = apply_sensor_geometry(obj, ir_distances, state_estimate)

            % 1. Apply the transformation to robot frame.

            ir_distances_rf = zeros(3,5);
            for i=1:5
                x_s = obj.sensor_placement(1,i);
                y_s = obj.sensor_placement(2,i);
                theta_s = obj.sensor_placement(3,i);

                R = obj.get_transformation_matrix(x_s,y_s,theta_s);
                ir_distances_rf(:,i) = R*[ir_distances(i); 0; 1];
            end

            % 2. Apply the transformation to world frame.

            [x,y,theta] = state_estimate.unpack();

            R = obj.get_transformation_matrix(x,y,theta);
            ir_distances_wf = R*ir_distances_rf;

            ir_distances_wf = ir_distances_wf(1:2,:);
        end

        function set_sensor_geometry(obj, robot)
            obj.sensor_placement = zeros(3,5);
            for i=1:5
                [x, y, theta] = robot.ir_array(i).location.unpack();
                obj.sensor_placement(:,i) = [x; y; theta];
            end
            obj.calibrated = true;
        end

        function R = get_transformation_matrix(obj, x, y, theta)
            R = [cos(theta) -sin(theta) x; sin(theta) cos(theta) y; 0 0 1];
        end

        function reset(obj)
            % Reset accumulated and previous error
            obj.E_k = 0;
            obj.e_k_1 = 0;
        end

    end

end

author: jdelacroix

repo: simiam

save_path: github-repos/MATLAB/jdelacroix-simiam


text: %% S2AxisFieldHarmonic


text:
function ldl_make
%LDL_MAKE compile LDL
%
% Example:
%   ldl_make       % compiles ldlsparse and ldlsymbol
%
% See also ldlsparse, ldlsymbol

% Copyright 2006-2007 by Timothy A. Davis, http://www.suitesparse.com

d = '' ;
if (~isempty (strfind (computer, '64')))
    d = '-largeArrayDims' ;
end
eval (sprintf ('mex -O %s -DLDL_LONG -I../../SuiteSparse_config -I../Include -output ldlsparse ../Source/ldl.c ldlmex.c', d)) ;
eval (sprintf ('mex -O %s -DLDL_LONG -I../../SuiteSparse_config -I../Include -output ldlsymbol ../Source/ldl.c ldlsymbolmex.c', d)) ;

eval (sprintf ('mex -O %s -DLDL_LONG -I../../SuiteSparse_config -I../Include -output ldlsparse_short ../Source/ldl.c ldlmex_short.c', d)) ;
eval (sprintf ('mex -O %s -DLDL_LONG -I../../SuiteSparse_config -I../Include -output ldlsymbol_extra ../Source/ldl.c ldlsymbol_extra.c', d)) ;
eval (sprintf ('mex -O %s -DLDL_LONG -I../../SuiteSparse_config -I../Include -output ldlnumeric ../Source/ldl.c ldlnumericmex.c', d)) ;
fprintf ('LDL successfully compiled.\n') ;

repo: gspbox


text:
function y = square( x )

%SQUARE    Square.
%   SQUARE(X) is an array of the same size as X, whose elements are the
%   squares of the elements of X.
%
%   Disciplined convex programming information:
%       If X is real, then SQUARE(X) is convex and nonmonotonic in X. If X
%       is complex, then SQUARE(X) is neither convex nor concave. Thus when
%       when use in CVX expressions, X must be real and affine.

narginchk(1,1);
y = x .* x;

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
%CONTENTS GRAPHS Simple Toolbox for manipulating Geometric Graphs.
% Version 0.6 04-Sep-2017 .
%
%   The aim of this package is to provide functions to easily create,
%   modify and display geometric graphs (geometric in a sense position
%   nodes are associated to geometric position in 2D or 3D).
%
%   Graph are represented by a structure with at least two arrays:
%   * NODES, which contains coordinates of each vertex
%   * EDGES, which contains indices of start and end vertex.
%   Some graph functions consider adjacency list, as a cell array where
%   each cell contains the indices of the neighbor vertices.
%
%   Others arrays may sometimes be used:
%   * FACES, which contains indices of vertices of each face (either a
%       double array, or a cell array)
%   * CELLS, which contains indices of faces of each cell.
%
%   An alternative representation is to use a structure, with fields:
%   * 'nodes'
%   * 'edges'
%   corresponding to the data described above.
%
%   Note that topological description of 2D graph is entirely contained in
%   EDGES array, and that NODES array is used only to display the graph.
%
%   Caution: this type of data structure is easy to create and to manage,
%   but may be very inefficient for some algorithms.
%
%   Graphs are usually considered as non-oriented in this package.
%
%
% Graph creation
%   delaunayGraph              - Graph associated to Delaunay triangulation of input points.
%   euclideanMST               - Build euclidean minimal spanning tree of a set of points.
%   prim_mst                   - Minimal spanning tree by Prim's algorithm.
%   knnGraph                   - Create the k-nearest neighbors graph of a set of points.
%   relativeNeighborhoodGraph  - Relative Neighborhood Graph of a set of points.
%   gabrielGraph               - Gabriel Graph of a set of points.
%
% Create graph from images
%   imageGraph                 - Create equivalent graph of a binary image.
%   imageBoundaryGraph         - Convert boundary of a 2D/3D binary image into a graph or mesh.
%
% Voronoi Graphs
%   voronoi2d                  - Compute a voronoi diagram as a graph structure.
%   boundedVoronoi2d           - Comptues a bounded voronoi diagram as a graph structure.
%   centroidalVoronoi2d        - Centroidal Voronoi tesselation within a polygon.
%   centroidalVoronoi2d_MC     - Centroidal Voronoi tesselation by Monte-Carlo.
%   boundedCentroidalVoronoi2d - Create a 2D Centroidal Voronoi Tesselation in a box.
%   cvtUpdate                  - Update germs of a CVT with given points.
%   cvtIterate                 - Update germs of a CVT using random points with given density.
%   meshEnergy                 - Computes the energy of a tesselation, as the sum of second area moments.
%
% Geodesic and shortest path operations
%   grShortestPath             - Find a shortest path between two nodes in the graph.
%   grPropagateDistance        - Propagates distances from a vertex to other vertices.
%   grVertexEccentricity       - Eccentricity of vertices in the graph.
%   graphDiameter              - Diameter of a graph.
%   graphPeripheralVertices    - Peripheral vertices of a graph.
%   graphCenter                - Center of a graph.
%   graphRadius                - Radius of a graph.
%   grFindGeodesicPath         - Find a geodesic path between two nodes in the graph.
%   grFindMaximalLengthPath    - Find a path that maximizes sum of edge weights.
%
% Graph processing (general applications)
%   adjacencyListToEdges       - Convert an adjacency list to an edge array.
%   pruneGraph                 - Remove all edges with a terminal vertex.
%   mergeGraphs                - Merge two graphs, by adding nodes, edges and faces lists.
%   grMergeNodes               - Merge two (or more) nodes in a graph.
%   grMergeMultipleNodes       - Simplify a graph by merging multiple nodes.
%   grMergeMultipleEdges       - Remove all edges sharing the same extremities.
%   grSimplifyBranches         - Replace branches of a graph by single edges.
%
% Filtering operations on Graph
%   grMean                     - Compute mean value from neighbour nodes.
%   grMedian                   - Compute median value from neighbour nodes.
%   grDilate                   - Morphological dilation on graph.
%   grErode                    - Morphological erosion on graph.
%   grClose                    - Morphological closing on graph.
%   grOpen                     - Morphological opening on graph.
%
% Operations for geometric graphs
%   grEdgeLengths              - Compute length of edges in a geometric graph.
%   grMergeNodeClusters        - Merge cluster of connected nodes in a graph.
%   grMergeNodesMedian         - Replace several nodes by their median coordinate.
%   clipGraph                  - Clip a graph with a rectangular area.
%   clipGraphPolygon           - Clip a graph with a polygon.
%   clipMesh2dPolygon          - Clip a planar mesh with a polygon.
%   addSquareFace              - Add a (square) face defined from its vertices to a graph.
%   grFaceToPolygon            - Compute the polygon corresponding to a graph face.
%   graph2Contours             - Convert a graph to a set of contour curves.
%
% Graph information
%   grNodeDegree               - Degree of a node in a (undirected) graph.
%   grNodeInnerDegree          - Inner degree of a node in a graph.
%   grNodeOuterDegree          - Outer degree of a node in a graph.
%   grAdjacentNodes            - Find list of nodes adjacent to a given node.
%   grAdjacentEdges            - Find list of edges adjacent to a given node.
%   grOppositeNode             - Return opposite node in an edge.
%   grLabel                    - Associate a label to each connected component of the graph.
%
% Graph management (low level operations)
%   grRemoveNode               - Remove a node in a graph.
%   grRemoveNodes              - Remove several nodes in a graph.
%   grRemoveEdge               - Remove an edge in a graph.
%   grRemoveEdges              - Remove several edges from a graph.
%
% Graph display
%   drawGraph                  - Draw a graph, given as a set of vertices and edges.
%   drawGraphEdges             - Draw edges of a graph.
%   fillGraphFaces             - Fill faces of a graph with specified color.
%   drawDigraph                - Draw a directed graph, given as a set of vertices and edges.
%   drawDirectedEdges          - Draw edges with arrow indicating direction.
%   drawEdgeLabels             - Draw values associated to graph edges.
%   drawNodeLabels             - Draw values associated to graph nodes.
%   drawSquareMesh             - Draw a 3D square mesh given as a graph.
%   patchGraph                 - Transform 3D graph (mesh) into a patch handle.
%
% Input/Output
%   readGraph                  - Read a graph from a text file.
%   writeGraph                 - Write a graph to an ascii file.

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2005-11-07
% Copyright 2005-2022 INRA - Cepia Software Platform

help(mfilename);

% Deprecated functions
%   grSimplifyBranches_old     - Replace branches of a graph by single edges.
%   grRemoveMultiplePoints     - Remove groups of close nodes in a graph.
%   boundaryGraph              - Get boundary of image as a graph.
%   gcontour2d                 - Creates contour graph of a 2D binary image.
%   gcontour3d                 - Create contour graph of a 3D binary image.

% Functions that requires further development
%   quiverToGraph              - Converts quiver data to quad mesh.

% Other functions

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom