text:
function error = vennX( data, resolution )
%
% function error = vennX( data, resolution )
%
% vennX - draws an area proportional venn diagram
%
% Draws a venn diagram (either two or three set) using
% circles, where the area of each region is proportional
% to the input values.
%
% INPUT:
%      data - a vector of counts for each set partition
%
%          For a two circle diagram:
%            data is a three element vector of:
%              |A|
%              |A and B|
%              |B|
%
%          For a three circle diagram:
%            data is a seven element vector of:
%               |A|
%               |A and B|
%               |B|
%               |B and C|
%               |C|
%               |C and A|
%               |A and B and C|
%
%       resolution - A measure of accuracy on the image,
%            typical values are within 1/100 to 1/1000 of
%            the maximum partition count.  Note that smaller
%            resolutions take longer compute time.
%
% OUTPUT:
%     error - the difference in area of each partition
%             between the actual area and the input vector
%
% EXAMPLES:
%
%     vennX( [ 106 26 257 ], .05 )
%
%     vennX( [ 75 143 210 ], .1 )
%
%     vennX( [ 16 3 10 6 19 8 3 ], .05 )
%
%
% COMMENTS:
%
%     The implementation is trivial, for the two circle case, two circles
%       are drawn to scale and moved closer and closer together until the
%       overlap is 'near' to the desired intersection. For the three
%       circle case, it is repeated three times, once for each pair of
%       circles.  Hence the two circle case is almost exact, whereas the
%       three circle case has much more error since the area |A and B and C|
%       is derived.  This means that large variations from random, especially
%       close to zero, will have larger errors, for example
%
%           vennX( [ 20 10 20 10 20 10 0], .1 )
%
%       as opposed to
%
%           vennX( [ 20 10 20 10 20 10 10], .1 )
%
% ENHANCEMENTS
%
%     The implementation could be sped up tremendously using a MRA
%     (multi-resolutional analysis) type algorithm.  e.g. start with a
%     resolution of .5 and find the distance between the circles, then use
%     that as a seed for a resolution of .1, then .05, .01, etc.
%
%     The error vector could be used as a measure to 'perturb' the position
%     of the third circle as to minimize the error.  This could be done
%     with a simple gradient descent method.  This would help the
%     exceptions described above where the distribution deviates from
%     random.
%
%     When small mishapen areas are drawn, the text does not match up, e.g.
%        vennX( [ 15 143 210 ], .1 )
%
%
%  Original implementation and method by Jeremy Heil, for the Order of
%  the Red Monkey, and the Tengu
%
%  Oct. 2004
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   figure;

   if length( data ) == 3
      dist = venn2( data(1), data(2), data(3), resolution );
      error = plot_venn2( data(1), data(2), data(3), resolution, dist );
      error = data - error';
   elseif length( data ) == 7
       %get the pairwise distance of each circle center from each other
       dist_A_B = venn2( data(1)+data(6), data(2)+data(7), data(3)+data(4), resolution );
       dist_B_C = venn2( data(3)+data(2), data(4)+data(7), data(5)+data(6), resolution );
       dist_A_C = venn2( data(1)+data(2), data(6)+data(7), data(4)+data(5), resolution );

       error = plot_venn3( data(1), data(2), data(3), data(4), data(5), data(6), data(7), ...
           resolution, dist_A_B, dist_B_C, dist_A_C );
       error = data - error';
   else
       'vennX error, data vector must be of length 3 or 7'
   end

    %change the colormap so that the background is white
    k = colormap;
    k = [ 1 1 1; k ];
    colormap(k)
    axis off

function error = plot_venn3( a, b, c, d, e, f, g, resolution, dist_A_B, dist_B_C, dist_A_C )

    r1 = sqrt( (a+b+f+g)/pi );
    r2 = sqrt( (b+c+d+g)/pi );
    r3 = sqrt( (d+e+f+g)/pi );

    %
    % Using a little geometry, think of the three circle's centers
    % as vertecies of a triangle.
    %
    y = ( dist_A_C^2 - dist_B_C^2 + dist_A_B^2 ) / 2 / dist_A_B;

    size_x = max( r1 + dist_A_B + r2, 2*r3 );
    size_y = max( r1, r2 ) + sqrt( dist_A_C^2 - y^2 ) + r3;

    %find the circle centers
    center1_x = r1;
    center1_y = max( r1, r2 );
    center2_x = r1 + dist_A_B;
    center2_y = center1_y;
    center3_x = r1 + y;
    center3_y = center1_y + sqrt( dist_A_C^2 - y^2 );

    [X,Y] = meshgrid( 0:resolution:size_x, 0:resolution:size_y );

    %draw the circles
    img = zeros( size(Y,1), size(Y,2) );

    img = img + 2 .* ( (X - center1_x).^2 + (Y - center1_y).^2 < r1^2 );
    img = img + 4 .* ( (X - center2_x).^2 + (Y - center2_y).^2 < r2^2 );
    img = img + 6 .* ( (X - center3_x).^2 + (Y - center3_y).^2 < r3^2 );

    clf
    imagesc(img)
    hold on

    %add the numbers and compute the error for each partition
    error = [];

    tmp = and( and( ...
           (X - center1_x).^2 + (Y - center1_y).^2 < r1^2, ...
           (X - center2_x).^2 + (Y - center2_y).^2 > r2^2 ), ...
           (X - center3_x).^2 + (Y - center3_y).^2 > r3^2 );
    [ i,j ] = find( tmp > 0 );
    error = [ error; sum(sum(tmp)) * resolution^2 ];
    text_x = mean(j);
    text_y = mean(i);
    h = text( text_x, text_y, num2str( a ) );
    set( h, 'FontWeight', 'bold' )

    tmp = and( and( ...
           (X - center1_x).^2 + (Y - center1_y).^2 < r1^2, ...
           (X - center2_x).^2 + (Y - center2_y).^2 < r2^2 ), ...
           (X - center3_x).^2 + (Y - center3_y).^2 > r3^2 );
    [ i,j ] = find( tmp > 0 );
    error = [ error; sum(sum(tmp)) * resolution^2 ];
    text_x = mean(j);
    text_y = mean(i);
    h = text( text_x, text_y, num2str( b ) );
    set( h, 'FontWeight', 'bold' )

    tmp = and( and( ...
           (X - center1_x).^2 + (Y - center1_y).^2 > r1^2, ...
           (X - center2_x).^2 + (Y - center2_y).^2 < r2^2 ), ...
           (X - center3_x).^2 + (Y - center3_y).^2 > r3^2 );
    [ i,j ] = find( tmp > 0 );
    error = [ error; sum(sum(tmp)) * resolution^2 ];
    text_x = mean(j);
    text_y = mean(i);
    h = text( text_x, text_y, num2str( c ) );
    set( h, 'FontWeight', 'bold' )

    tmp = and( and( ...
           (X - center1_x).^2 + (Y - center1_y).^2 > r1^2, ...
           (X - center2_x).^2 + (Y - center2_y).^2 < r2^2 ), ...
           (X - center3_x).^2 + (Y - center3_y).^2 < r3^2 );
    [ i,j ] = find( tmp > 0 );
    error = [ error; sum(sum(tmp)) * resolution^2 ];
    text_x = mean(j);
    text_y = mean(i);
    h = text( text_x, text_y, num2str( d ) );
    set( h, 'FontWeight', 'bold' )

    tmp = and( and( ...
           (X - center1_x).^2 + (Y - center1_y).^2 > r1^2, ...
           (X - center2_x).^2 + (Y - center2_y).^2 > r2^2 ), ...
           (X - center3_x).^2 + (Y - center3_y).^2 < r3^2 );
    [ i,j ] = find( tmp > 0 );
    error = [ error; sum(sum(tmp)) * resolution^2 ];
    text_x = mean(j);
    text_y = mean(i);
    h = text( text_x, text_y, num2str( e ) );
    set( h, 'FontWeight', 'bold' )

    tmp = and( and( ...
           (X - center1_x).^2 + (Y - center1_y).^2 < r1^2, ...
           (X - center2_x).^2 + (Y - center2_y).^2 > r2^2 ), ...
           (X - center3_x).^2 + (Y - center3_y).^2 < r3^2 );
    [ i,j ] = find( tmp > 0 );
    error = [ error; sum(sum(tmp)) * resolution^2 ];
    text_x = mean(j);
    text_y = mean(i);
    h = text( text_x, text_y, num2str( f ) );
    set( h, 'FontWeight', 'bold' )

    tmp = and( and( ...
           (X - center1_x).^2 + (Y - center1_y).^2 < r1^2, ...
           (X - center2_x).^2 + (Y - center2_y).^2 < r2^2 ), ...
           (X - center3_x).^2 + (Y - center3_y).^2 < r3^2 );
    [ i,j ] = find( tmp > 0 );
    error = [ error; sum(sum(tmp)) * resolution^2 ];
    text_x = mean(j);
    text_y = mean(i);
    h = text( text_x, text_y, num2str( g ) );
    set( h, 'FontWeight', 'bold' )

function error = plot_venn2( a, b, c, resolution, dist )

    r1 = sqrt( (a+b)/pi );
    r2 = sqrt( (b+c)/pi );

    size_x = r1 + dist + r2;
    size_y = max( 2*r1, 2*r2 );

    center1_x = r1;
    center1_y = size_y/2;
    center2_x = r1 + dist;
    center2_y = size_y/2;;

    [X,Y] = meshgrid( 0:resolution:size_x, 0:resolution:size_y );

    %draw the two circles and the overlap region
    img = zeros( size(Y,1), size(Y,2) );

    img = img + 2 .* ( (X - center1_x).^2 + (Y - center1_y).^2 < r1^2 );
    img = img + 4 .* ( (X - center2_x).^2 + (Y - center2_y).^2 < r2^2 );

    imagesc(img)
    hold on

    %
    % We want to draw the numbers at the center of mass
    % do this by computing the average x and y coordinates of
    % the center of each partition piece.
    %
    % Compute the error for each partition as the difference
    % between the area we meant to draw and the actual area
    % that was drawn
    %
    error = [];
    tmp = and( ...
           (X - center1_x).^2 + (Y - center1_y).^2 < r1^2, ...
           (X - center2_x).^2 + (Y - center2_y).^2 > r2^2 );
    [ i,j ] = find( tmp > 0 );
    error = [ error; sum(sum(tmp)) * resolution^2 ];
    text_x = mean(j);
    text_y = mean(i);
    h = text( text_x, text_y, num2str( a ) );
    set( h, 'FontWeight', 'bold' )

    tmp = and( ...
           (X - center1_x).^2 + (Y - center1_y).^2 < r1^2, ...
           (X - center2_x).^2 + (Y - center2_y).^2 < r2^2 );
    [ i,j ] = find( tmp > 0 );
    error = [ error; sum(sum(tmp)) * resolution^2 ];
    text_x = mean(j);
    text_y = mean(i);
    h = text( text_x, text_y, num2str( b ) );
    set( h, 'FontWeight', 'bold' )

    tmp = and( ...
           (X - center1_x).^2 + (Y - center1_y).^2 > r1^2, ...
           (X - center2_x).^2 + (Y - center2_y).^2 < r2^2 );
    [ i,j ] = find( tmp > 0 );
    error = [ error; sum(sum(tmp)) * resolution^2 ];
    text_x = mean(j);
    text_y = mean(i);
    h = text( text_x, text_y, num2str( c ) );
    set( h, 'FontWeight', 'bold' )

function dist = venn2( a, b, c, resolution )
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %VENN2
    %  dist = venn2( a, b, c, resolution )
    %
    % Computes the distance between the centers of
    % the two venn diagram circles.
    %
    % a          - values of A
    % b          - value of A and B
    % c          - value of B
    % resolution - measure of error
    %
    % dist       - the distance between the two centers
    %
    %  Does this by plotting the two circles in an
    % image with the specified resolution and
    % moving the centers towards each other until
    % the area of intersection is nearest the value
    % of b
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    r1 = sqrt( (a+b)/pi );
    r2 = sqrt( (b+c)/pi );

    size_x = 2*r1+2*r2;
    size_y = max( 2*r1, 2*r2 );

    center1_x = r1;
    center1_y = size_y/2;
    center2_x = 2*r1 + r2;
    center2_y = size_y/2;;

    [X,Y] = meshgrid( 0:resolution:size_x, 0:resolution:size_y );

    for new_center = (2*r1 + r2):-resolution:r1

        img = zeros( size(Y,1), size(Y,2) );
        img = and( (X - center1_x).^2 + (Y - center1_y).^2 < r1^2, ...
            (X - new_center).^2 + (Y - center2_y).^2 < r2^2 );

        if sum(sum(img)) * resolution^2 > b
            break
        end
    end

    dist = new_center - center1_x;

author: opencobra

repo: cobratoolbox


text:
function sF = drho(sF)
% first derivative in direction rho

s = size(sF);
sF = reshape(sF, []);

fhat = zeros(size(sF.fhat));
for m = 0:sF.bandwidth
  fhat(m*(m+1)+(-m:m)+1, :) = 1i*(-m:m)'.*sF.fhat(m*(m+1)+(-m:m)+1, :);
end
sF = S2FunHarmonic(fhat);

sF = reshape(sF, s);

end

author: mtex-toolbox

repo: mtex

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


text:
% This example shows how to calculate and plot both the fundamental
% quasi-TE eigenmode and quasi-TM eigenmode of an example 3-layer
% ridge waveguide using the semivectorial eigenmode solver.

% 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

[x,y,xc,yc,nx,ny,eps,edges] = waveguidemesh([n1,n2,n3],[h1,h2,h3], ...
                                            rh,rw,side,dx,dy);

% First, consider the quasi-TE mode:

[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

% Next, consider the quasi-TM mode:

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

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

figure(2);
contourmode(x,y,Ey);
title('Ey (TM mode)'); xlabel('x'); ylabel('y');
for v = edges, line(v{:}); end

author: Sable

repo: mcbench-benchmarks


text:
function [ diff_median_mask]= diff_median( i01, window_side, thres)

i01_size= size( i01);
diff_median_mat= zeros( [ i01_size]);

for j= 1+ window_side: i01_size(1)- window_side                             % rows
    for k= 1+ window_side: i01_size(2)- window_side                         % columns

        i01_temp= i01( j- window_side: j+ window_side, k- window_side:...
            k+ window_side);
        center= i01_temp( window_side+ 1, window_side+ 1);
        median_val= median( i01_temp( :));
        diff_median_mat( j, k)= abs( median_val- center);

    end
end

diff_median_mask= diff_median_mat>= ( thres* max( diff_median_mat( :)));

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function [cx, cy, w, h] = get_axis_aligned_BB(region)
% GETAXISALIGNEDBB extracts an axis aligned bbox from the ground truth REGION with same area as the rotated one
nv = numel(region);
assert(nv==8 || nv==4);

if nv==8
    cx = mean(region(1:2:end));
    cy = mean(region(2:2:end));
    x1 = min(region(1:2:end));
    x2 = max(region(1:2:end));
    y1 = min(region(2:2:end));
    y2 = max(region(2:2:end));
    A1 = norm(region(1:2) - region(3:4)) * norm(region(3:4) - region(5:6));
    A2 = (x2 - x1) * (y2 - y1);
    s = sqrt(A1/A2);
    w = s * (x2 - x1) + 1;
    h = s * (y2 - y1) + 1;
else
    x = region(1);
    y = region(2);
    w = region(3);
    h = region(4);
    cx = x+w/2;
    cy = y+h/2;
end

author: bertinetto

repo: cfnet

save_path: github-repos/MATLAB/bertinetto-cfnet


text:
function dst_scalars=nsst_scalars(L,shear_f,lpfilt)
% Since the nonsubsampled descrite shearlet transform is not orthogonal
% this function computes the noise level scalars of the transform with
% assigned parameters.
%
% Inputs:
%
% L                      - size of image decomposition
%
% shear_f                - the cell array containing the shearing filters
%
%
% lpfilt                 - lpfilt is the filter to be used for the Laplacian
%                          Pyramid/ATrous decomposition using the codes
%                          written by Arthur Cunha
%
% Output:
%
% dst_scalars            - the cell array containing the scalars of
%                          estimated noise levels for a white Gaussian noise
%                          of standard deviation 1 transform coefficients
%                          using Monte Carlo method with one iteration.
%
% Code contributors: Glenn R. Easley, Demetrio Labate, and Wang-Q Lim.
% Copyright 2011 by Glenn R. Easley. All Rights Reserved.
%

noise=randn(L,L);
level=length(shear_f);

% LP decomposition
y_noise = atrousdec(noise,lpfilt,level);

dst_scalars=cell(1,level+1);
dst_noise=y_noise{1};
dst_scalars{1}=median(abs(dst_noise(:) - median(dst_noise(:))))/.6745;

for i=1:level,
    l=size(shear_f{i},3);
    for k=1:l,
        dst_noise=conv2p(shear_f{i}(:,:,k),y_noise{i+1});
        dst_scalars{i+1}(k)=median(abs(dst_noise(:) - median(dst_noise(:))))/.6745;
    end
end

author: thomaskuestner


text:
function [V,F] = subdivided_sphere(iters,varargin)
  % SUBDIVIDED_SPHERE Generate a sphere by iteratively subdividing a
  % icosahedron in-plane and normalizing vertex locations to lie on the sphere.
  %
  % [V,F] = subdivided_sphere(iters,varargin)
  %
  % Inputs:
  %   iters  number of subdivision iterations (0 produces 12-vertex
  %     icosahedron)
  %   Optional:
  %     'Base'  followed by base shape.
  %     'Radius'  followed by scalar radius (multiplied against final V) {1}
  %     'SubdivisionMethod' followed by either:
  %       'loop'
  %       'sqrt3'
  %       {'upsample'}
  % Outputs:
  %   V  #V by 3 list of mesh vertices
  %   F  #F by 3 list of face indices into V
  %
  % See also: upsample, loop
  %

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

  switch base
  case 'icosahedron'
    [V,F] = icosahedron();
  case 'octahedron'
    [V,F] = octahedron();
  otherwise
    % Compute the 12 vertices
    phi = (1+sqrt(5))/2;  % Golden ratio
    V = [0   1  phi
            0  -1  phi
            0   1 -phi
            0  -1 -phi
            1  phi  0
          -1  phi  0
            1 -phi  0
          -1 -phi  0
            phi 0   1
          -phi 0   1
            phi 0  -1
          -phi 0  -1];
    % Scale to required radius
    V = V/(sqrt(1+phi^2));
    % Define the adjacency matrix
    F = [1  2  9
          1  9  5
          1  5  6
          1  6  10
          1  10 2
          2  7  9
          9  7  11
          9  11 5
          5  11 3
          5  3  6
          6  3  12
          6  12 10
          10 12 8
          10 8  2
          2  8  7
          4  7  8
          4  8  12
          4  12 3
          4  3  11
          4  11 7];
  end

  V = normalizerow(V);
  for iter = 1:iters
    switch subdivision_method
    case 'upsample'
      [V,F] = upsample(V,F);
    case 'loop'
      [V,F] = loop(V,F);
    case 'sqrt3'
      [V,F] = sqrt3(V,F);
    end
    V = normalizerow(V);

  end
  V = V*radius;

end

author: alecjacobson

repo: gptoolbox


text:
function [vname]=readcpx(fname,width,lines,endian,skipbytes,precision)
%READCPX Read complex interleaved files
%   A = READCPX(FILE_NAME, WIDTH, LINES, ENDIAN, SKIPBYTES, PRECISION)
%       FILE_NAME and WIDTH are compulsory
%
%   Andy Hooper, Nov 2003

if nargin < 2
  error('syntax is A = readcpx(FILE_NAME, WIDTH, LINES, ENDIAN, SKIPBYTES, PRECISION)')
end

if nargin < 3
    lines = inf;
end

if nargin < 4
    endian = 'n';
end

if nargin < 5
    skipbytes = 0;
end

if nargin < 6
    precision = 'float';
end

fid=fopen(fname,'r',endian);

fseek(fid,skipbytes,-1);

vname=fread(fid,[width*2,lines],[precision,'=>single']).';

vname=(complex(vname(:,1:2:end),vname(:,2:2:end)));

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


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

%%
clear; close all; clc;

%% Syntax
% |[Es,Vs,Cs]=element2lattice(E,V,cPar);|

%% Description
% This function converts an element description (elements and vertices i.e
% nodes into a lattice structure. The lattice structure is returned as a
% hexahederal mesh.

%% Examples
%

%%
% Plot settings
cMap=gjet(4);
fontSize=15;

%% Example: Creating a lattice structure on a hexahedral element

%%
% Creating example geometry.
boxDim=[1 1 1];
boxEl=[1 1 1];
[meshStruct]=hexMeshBox(boxDim,boxEl);
E=meshStruct.E;
V=meshStruct.V;
F=meshStruct.F;
[indBoundary]=tesBoundary(F);

%%
% Create lattice structure
controlParameter.shrinkFactor=0.2; %Strut sides are formed by shrinking the input mesh faces by this factor
controlParameter.numDigitKeep=5; %used for merging nodes
controlParameter.meshType='quad'; %desired output mesh type
controlParameter.indBoundary=indBoundary; %indices of the boundary faces
controlParameter.latticeSide=1; %1=side 1 the edge lattice, 2=side 2 the dual lattice to the edge lattice
[Fs1,Vs1,Cs1]=element2lattice(E,V,controlParameter); %Get lattice structure

%%
% Visualizing input mesh and lattice structures

cFigure;
title('An edge lattice structure on a hexahedral element','fontSize',fontSize)
hold on;
hp1=gpatch(F,V,0.5*ones(1,3),'k',0.25,4);
hp2=gpatch(Fs1(Cs1==0,:),Vs1,Cs1(Cs1==0));
hp3=gpatch(Fs1(Cs1==1,:),Vs1,Cs1(Cs1==1));
legend([hp1,hp2,hp3],'Input mesh','Inner faces','Boundary faces')
colormap(cMap);
cLim=caxis;
axisGeom(gca,fontSize);
camlight headlight; lighting flat;

drawnow;

%% Example: Different lattice sides

%%
% Creating example geometry.
boxDim=[1 1 1];
boxEl=[1 1 1];
[meshStruct]=hexMeshBox(boxDim,boxEl);
E=meshStruct.E;
V=meshStruct.V;
F=meshStruct.F;
[indBoundary]=tesBoundary(F);

%%
% Compute other "lattice side"
controlParameter.shrinkFactor=0.2; %Strut sides are formed by shrinking the input mesh faces by this factor
controlParameter.numDigitKeep=5; %used for merging nodes
controlParameter.meshType='quad'; %desired output mesh type
controlParameter.indBoundary=indBoundary; %indices of the boundary faces
controlParameter.latticeSide=2; %1=side 1 the edge lattice, 2=side 2 the dual lattice to the edge lattice
[Fs2,Vs2,Cs2]=element2lattice(E,V,controlParameter); %Get lattice structure

%%
% Visualizing input mesh and lattice structures

cFigure;
hs=subplot(2,2,1);
title('The input mesh','fontSize',fontSize)
hold on;
gpatch(F,V,0.5*ones(1,3),'k',0.25,3);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;
ha=axis; axis off;

subplot(2,2,2);
title('The two complementary lattice structures','fontSize',fontSize)
hold on;
gpatch(F,V,0.5*ones(1,3),'k',0.25,3);
gpatch(Fs1,Vs1,Cs1);
gpatch(Fs2,Vs2,Cs2);
colormap(cMap);
cLim=caxis;
axisGeom(gca,fontSize);
camlight headlight; lighting flat;
axis(ha); axis off;

subplot(2,2,3);
title('Lattice side 1','fontSize',fontSize)
hold on;
gpatch(F,V,0.5*ones(1,3),'k',0.25,3);
gpatch(Fs1,Vs1,Cs1);
colormap(cMap);
caxis(cLim);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;
axis(ha); axis off;

subplot(2,2,4);
title('Lattice side 2','fontSize',fontSize)
hold on;
gpatch(F,V,0.5*ones(1,3),'k',0.25,3);
gpatch(Fs2,Vs2,Cs2);
colormap(cMap);
caxis(cLim);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;
axis(ha); axis off;

drawnow;

%% Example: Refining the lattice structure

%%
% Refining a lattice structure by first refining the input mesh

cFigure;
gtitle('Lattice structure refinement',fontSize);

for q=1:4
    if q>1
        [E,V]=subHex(E,V,1); %Refine input mesh
        [F]=element2patch(E); %Patch data for plotting
        [indBoundary]=tesBoundary(F);
        controlParameter.indBoundary=indBoundary; %indices of the boundary faces
    end

    [Fs1,Vs1,Cs1]=element2lattice(E,V,controlParameter); %Get lattice structure

    % Visualizing input mesh and lattice structures

    subplot(2,2,q);
    hold on; title(['Split iteration ',num2str(q)]);
    gpatch(Fs1,Vs1,Cs1);
    colormap(cMap);
    axisGeom(gca,fontSize);
    camlight headlight; lighting flat;

end
drawnow;

%% Example: Creating a lattice structure on tetrahedral elements

%%
% Creating example geometry.

[V,~]=platonic_solid(1,1); %A single tetrahedron
E=[2 3 4 1]; %The element description
[E,V]=subTet(E,V,1); %Refine the tetrahedron once
[F]=element2patch(E); %Patch data for plotting
[indBoundary]=tesBoundary(F); %Get boundary face indices

%%
% Create lattice structure
controlParameter.shrinkFactor=0.2; %Strut sides are formed by shrinking the input mesh faces by this factor
controlParameter.numDigitKeep=5; %used for merging nodes
controlParameter.meshType='quad'; %desired output mesh type
controlParameter.indBoundary=indBoundary; %indices of the boundary faces

controlParameter.latticeSide=1; %1=side 1 the edge lattice, 2=side 2 the dual lattice to the edge lattice
[Fs1,Vs1,Cs1]=element2lattice(E,V,controlParameter); %Get lattice structure

%%
% Visualizing input mesh and lattice structures

cFigure;
title('An edge lattice structure on a hexahedral element','fontSize',fontSize)
hold on;
hp1=gpatch(F,V,0.5*ones(1,3),'k',0.25,4);
hp2=gpatch(Fs1(Cs1==0,:),Vs1,Cs1(Cs1==0));
hp3=gpatch(Fs1(Cs1==1,:),Vs1,Cs1(Cs1==1));
legend([hp1,hp2,hp3],'Input mesh','Inner faces','Boundary faces')
colormap(cMap);
cLim=caxis;
axisGeom(gca,fontSize);
camlight headlight; lighting flat;

drawnow;

%% Example: Different lattice sides

%%
% Compute other "lattice side"
controlParameter.latticeSide=2; %1=side 1 the edge lattice, 2=side 2 the dual lattice to the edge lattice
[Fs2,Vs2,Cs2]=element2lattice(E,V,controlParameter); %Get lattice structure

%%
% Visualizing input mesh and lattice structures

cFigure;
hs=subplot(2,2,1);
title('The input mesh','fontSize',fontSize)
hold on;
gpatch(F,V,0.5*ones(1,3),'k',0.25,3);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;
ha=axis; axis off;

subplot(2,2,2);
title('The two complementary lattice structures','fontSize',fontSize)
hold on;
gpatch(F,V,0.5*ones(1,3),'k',0.25,3);
gpatch(Fs1,Vs1,Cs1);
gpatch(Fs2,Vs2,Cs2);
colormap(cMap);
cLim=caxis;
axisGeom(gca,fontSize);
camlight headlight; lighting flat;
axis(ha); axis off;

subplot(2,2,3);
title('Lattice side 1','fontSize',fontSize)
hold on;
gpatch(F,V,0.5*ones(1,3),'k',0.25,3);
gpatch(Fs1,Vs1,Cs1);
colormap(cMap);
caxis(cLim);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;
axis(ha); axis off;

subplot(2,2,4);
title('Lattice side 2','fontSize',fontSize)
hold on;
gpatch(F,V,0.5*ones(1,3),'k',0.25,3);
gpatch(Fs2,Vs2,Cs2);
colormap(cMap);
caxis(cLim);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;
axis(ha); axis off;

drawnow;

%% Example: Changing lattice strut thickness (and porosity)

%%
% The strut thickness of the lattice depends on the shrinkfactor.

% Create lattice structure
controlParameter.numDigitKeep=5; %used for merging nodes
controlParameter.meshType='quad'; %desired output mesh type
controlParameter.indBoundary=indBoundary; %indices of the boundary faces

shrinkFactorSet=linspace(0.1,0.5,4);

for latticeSide=1:2
    controlParameter.latticeSide=latticeSide;
    cFigure;
    gtitle('Lattice structure porosity control',fontSize);

    for q=1:4
        controlParameter.shrinkFactor=shrinkFactorSet(q); %Strut sides are formed by shrinking the input mesh faces by this facto

        [Fs1,Vs1,Cs1]=element2lattice(E,V,controlParameter); %Get lattice structure

        % Visualizing input mesh and lattice structures
        subplot(2,2,q);
        hold on; title(['Shrink factor ',num2str(pround(shrinkFactorSet(q),3))]);
        gpatch(Fs1,Vs1,Cs1);
        colormap(cMap);
        axisGeom(gca,fontSize);
        camlight headlight; lighting flat;
    end
    drawnow;
end

%% Example: Changing output surface mesh type

%%
% Different output mesh types are available i.e. quadrilateral and
% triangular faces and also hexahedral elements.

meshTypeSet={'quad','tri'};

boxDim=[1 1 1];
boxEl=[1 1 1];
[meshStruct]=hexMeshBox(boxDim,boxEl);
E=meshStruct.E;
V=meshStruct.V;
F=meshStruct.F;
[indBoundary]=tesBoundary(F);

clear controlParameter

% Create lattice structure
controlParameter.latticeSide=2;
controlParameter.numDigitKeep=5; %used for merging nodes
controlParameter.indBoundary=indBoundary; %indices of the boundary faces
controlParameter.shrinkFactor=0.3;

cFigure;
gtitle('Lattice structure mesh output type control',fontSize);

for q=1:numel(meshTypeSet)
    controlParameter.meshType=meshTypeSet{q}; %The current mesh type

    [Fs1,Vs1,Cs1]=element2lattice(E,V,controlParameter); %Get lattice structure

    % Visualizing input mesh and lattice structures
    subplot(1,numel(meshTypeSet),q);
    hold on; title(['Mesh type ',meshTypeSet{q}]);
    gpatch(Fs1,Vs1,Cs1,'k',0.8,2);
    colormap(gjet(250));
    axisGeom(gca,fontSize);
    camlight headlight; lighting flat;
end
drawnow;

%% Example: Exporting hexahedral elements instead of surface elements

%%
% The |element2patch| function can also be used to export hexahedral
% elements directly by setting the meshType parameter to 'hex'.
% Furthermore, it is possible to subdevide the ellongated hexahedral
% elements using a certain number of split iterations. Below is an example
% of hexahedral element output with increasing number of split iterations
% used for the ellongated hexahedral elements.

boxDim=[1 1 1];
boxEl=[1 1 1];
[meshStruct]=hexMeshBox(boxDim,boxEl);
E=meshStruct.E;
V=meshStruct.V;
F=meshStruct.F;
[indBoundary]=tesBoundary(F);

clear controlParameter

% Create lattice structure
controlParameter.latticeSide=2;
controlParameter.numDigitKeep=5; %used for merging nodes
controlParameter.indBoundary=indBoundary; %indices of the boundary faces
controlParameter.shrinkFactor=0.25;
controlParameter.meshType='hex';

hexSplitSet=[0 1 2];
cFigure;
gtitle('Lattice structure hex element output and element mesh refinement',fontSize);

c=1; %counter for plotting
for latticeSide=1:2
    controlParameter.latticeSide=latticeSide;
    for q=1:numel(hexSplitSet)
        controlParameter.hexSplit=hexSplitSet(q); %The current mesh type

        [Es1,Vs1,Cs1]=element2lattice(E,V,controlParameter); %Get lattice structure

        [Fs1,Cs1F]=element2patch(Es1,Cs1); %Patch data for plotting

        % Visualizing input mesh and lattice structures
        subplot(2,numel(hexSplitSet),c);
        hold on; title(['hexSplit=',num2str(hexSplitSet(q))]);
        gpatch(Fs1,Vs1,Cs1F,'k',0.8,2);
        colormap(gjet(250));
        axisGeom(gca,fontSize);
        camlight headlight; lighting flat;
        c=c+1;
    end
end
drawnow;

%%
%

[V,F]=platonic_solid(1,1);
E=[1 2 4 3];

%
% boxDim=[1 1 1];
% boxEl=[1 1 1];
% [meshStruct]=hexMeshBox(boxDim,boxEl);
% E=meshStruct.E;
% V=meshStruct.V;
% F=meshStruct.F;
[indBoundary]=tesBoundary(F);

clear controlParameter

% Create lattice structure
controlParameter.latticeSide=2;
controlParameter.numDigitKeep=5; %used for merging nodes
controlParameter.indBoundary=indBoundary; %indices of the boundary faces
controlParameter.shrinkFactor=0.25;
controlParameter.meshType='hex';

hexSplitSet=[0 1 2];
cFigure;
gtitle('Lattice structure hex element output and element mesh refinement',fontSize);

c=1; %counter for plotting
for latticeSide=1:2
    controlParameter.latticeSide=latticeSide;
    for q=1:numel(hexSplitSet)
        controlParameter.hexSplit=hexSplitSet(q); %The current mesh type

        [Es1,Vs1,Cs1]=element2lattice(E,V,controlParameter); %Get lattice structure

        [Fs1,Cs1F]=element2patch(Es1,Cs1); %Patch data for plotting

        % Visualizing input mesh and lattice structures
        subplot(2,numel(hexSplitSet),c);
        hold on; title(['hexSplit=',num2str(hexSplitSet(q))]);
        gpatch(Fs1,Vs1,Cs1F,'k',0.8,2);
        colormap(gjet(250));
        axisGeom(gca,fontSize);
        camlight headlight; lighting flat;
        c=c+1;
    end
end
drawnow;

%% Example: Creating lattice structure variations through input mesh conversion/subdevission

cPar.shrinkFactor=0.25;
cPar.numDigitKeep=5;
cPar.meshType='quad';

cParSmooth.Method='HC';
cParSmooth.n=10;

cFigure;
c=1;
for latticeSide=1:2
    for testCase=1:3

        switch testCase
            case 1
                boxEl=[1 1 1];
                [meshStruct]=hexMeshBox(boxDim,boxEl);
                E=meshStruct.E;
                V=meshStruct.V;
                F=meshStruct.F;
                Ft=F;
                Vt=V;
                [E,V]=subHex(E,V,1,1);
            case 2
                boxEl=[1 1 1];
                [meshStruct]=hexMeshBox(boxDim,boxEl);
                E=meshStruct.E;
                V=meshStruct.V;
                F=meshStruct.F;
                Ft=F;
                Vt=V;
                [E,V]=subHex(E,V,1,1);
                [E,V,~]=hex2tet(E,V,[],1);
            case 3
                boxEl=[1 1 1];
                [meshStruct]=hexMeshBox(boxDim,boxEl);
                E=meshStruct.E;
                V=meshStruct.V;
                F=meshStruct.F;
                Ft=F;
                Vt=V;
                [E,V,~]=hex2tet(E,V,[],1);
                [E,V]=tet2hex(E,V);
        end

        %Get boundary indices
        [F]=element2patch(E); %Patch data for plotting
        [indBoundary]=tesBoundary(F); %Boundary indices

        %Compute lattice
        cPar.latticeSide=latticeSide;
        cPar.indBoundary=indBoundary;
        [Fn,Vn,Cn]=element2lattice(E,V,cPar);

%         %Refine mesh
%         [Fn,Vn]=subQuad(Fn,Vn,1);
%         Cn=repmat(Cn,[4 1]); %Replicate color info
%
%         %Smoothen
%         indRigid=Fn(Cn==1,:);
%         indRigid=unique(indRigid(:)); %Indices for boundary elements to hold on to
%         cParSmooth.RigidConstraints=indRigid;
%         if cParSmooth.n>0
%             [Vn]=tesSmooth(Fn,Vn,[],cParSmooth); %Smoothen mesh
%         end

        %Visualize
        subplot(2,3,c); hold on;
        gpatch(Ft,Vt,0.5*ones(1,3),'k',0.25);
        gpatch(Fn,Vn,Cn,'k',1);
        colormap(cMap);
        axisGeom(gca,fontSize);
        camlight headlight;
        lighting flat;

        c=c+1;
    end
end
drawnow;

%% Example: Create lattice structures on arbitry input meshes
% Below is an example of a general tetrahedral mesh

%%

testCase=2;
switch testCase
    case 1
        [F,V,~]=geoSphere(2,1); % Building a geodesic dome surface model
    case 2
        [F,V]=stanford_bunny('g'); %Bunny
        V_mean=mean(V,1);
        V=V-V_mean(ones(size(V,1),1),:);
end

%%
% Using tetgen to create a tetrahedral mesh (see |HELP_runTetGen|)

stringOpt='-pq1.2AaY';

inputStruct.stringOpt=stringOpt;
inputStruct.Faces=fliplr(F);
inputStruct.Nodes=V;
inputStruct.holePoints=[];
inputStruct.faceBoundaryMarker=ones(size(F,1),1); %Face boundary markers
inputStruct.regionPoints=getInnerPoint(F,V); %region points
inputStruct.regionA=tetVolMeanEst(F,V);
inputStruct.minRegionMarker=2; %Minimum region marker

[meshOutput]=runTetGen(inputStruct); %Run tetGen

Fb=meshOutput.facesBoundary;
Cb=meshOutput.boundaryMarker;
V=meshOutput.nodes;
CE=meshOutput.elementMaterialID;
E=meshOutput.elements;
F=meshOutput.faces;
[indBoundary]=tesBoundary(F); %Boundary indices

%%
% Create lattice structure

% Define spatially varying (on nodes) shrink factor
s=V(:,1);
s=s-min(s(:));
s=s./max(s(:));
s=(s*0.6)+0.05;

clear controlParameter
controlParameter.shrinkFactor=s; %Strut sides are formed by shrinking the input mesh faces by this factor
controlParameter.numDigitKeep=5; %used for merging nodes
controlParameter.meshType='quad'; %desired output mesh type
controlParameter.indBoundary=indBoundary; %indices of the boundary faces
controlParameter.latticeSide=1; %1=side 1 the edge lattice, 2=side 2 the dual lattice to the edge lattice

[Fs1,Vs1,Cs1]=element2lattice(E,V,controlParameter); %Get lattice structure

%%
% PLOTTING MODEL

%Selecting half of the model to see interior
Y=V(:,2); YE=mean(Y(E),2);
logicCutView=YE>mean(Y);
[Fs,Cs]=element2patch(E(logicCutView,:),CE(logicCutView),'tet4');

cFigure;
hold on;
title('Cut view of tetrahedral mesh model','FontSize',fontSize);
gpatch(Fb,V,0.5*ones(1,3),'none',0.5);
gpatch(Fs,V,Cs,'k',1);
camlight headlight;
axisGeom(gca,fontSize);
axis off;
colormap(cMap);
drawnow;

%%
% Visualize lattice structure
cFigure;
hold on;
title('Lattice structure on arbitrary input mesh','FontSize',fontSize);
% gpatch(Fb,V,0.5*ones(1,3),'none',0.5);
gpatch(Fs1,Vs1,Cs1,'none',1);
camlight headlight; lighting flat;
axisGeom(gca,fontSize);
axis off;
colormap(cMap);
drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function pos = findhashsorted(x,searchfor)

if numel(x) < 1000
    % Short enough use built-in despite being sorted...
    pos = find(x == searchfor);
else
    % Long enough to warrant a bisection search, JITs fine
    if searchfor>x(end)
        pos = [];
        return
    elseif searchfor < x(1)
        pos = [];
        return
    elseif searchfor == x(1)
        pos = 1;
        return
    elseif searchfor == x(end)
        pos = numel(x);
        return
    end
    low=1;
    high=numel(x);
    pos=[];
    d=numel(x);
    while high-low > 1
        mid = floor((low+high)/2);
        val = x(mid);
        if val > searchfor
            high = mid;
        elseif val < searchfor
            low = mid;
        else
            pos = mid;
            break
        end
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function J = calcJ(camera, imuState_k, camStates_k)
% Jacobian of feature observations w.r.t. feature locations

    C_CI = quatToRotMat(camera.q_CI);
    C_IG = quatToRotMat(imuState_k.q_IG);

    J = zeros(6, 12 + 6*size(camStates_k,2));
    J(1:3,1:3) = C_CI;
    J(4:6,1:3) = crossMat(C_IG' * camera.p_C_I);
    J(4:6,10:12) = eye(3);

end

author: utiasSTARS

repo: msckf-swf-comparison


text:
function [param, names] = gpsimCandidateExtractParam(model)

% GPSIMCANDIDATEEXTRACTPARAM Extract the parameters of a GPSIM model.
% FORMAT
% DESC extracts the model parameters from a structure containing
% the information about a Gaussian process for single input motif
% modelling.
% ARG model : the model structure containing the information about
% the model.
% RETURN params : a vector of parameters from the model.
%
% SEEALSO : gpsimCreate, gpsimAddCandidate, gpsimCandidateExpandParam, modelExtractParam
%
% COPYRIGHT : Neil D. Lawrence, 2007

% SHEFFIELDML

if nargout>1
  [param, names] = kernExtractParam(model.candidate.kern);
  for i=1:length(model.candidate.mu);
    names{end+1}=['Basal transcription ' num2str(i)];
  end
else
  param = kernExtractParam(model.candidate.kern);
end
fhandle = str2func([model.candidate.bTransform 'Transform']);
param = [param fhandle(model.candidate.B, 'xtoa')];

if isfield(model, 'fix')
  for i = 1:length(model.candidate.fix)
    param(model.candidate.fix(i).index) = model.candidate.fix(i).value;
  end
end
param = real(param);

% Remove main kernel parameters.
for i = model.kern.nParams:-1:1
  if nargout>1
    names(i)=[];
  end
  param(i) = [];
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%RNA Decay
%
% @wholeCellModelID Process_RNADecay
% @name             RNA Decay
% @description
%   Biology
%   ===============
%   In presence of ribonucleases such as ribonuclease R (MG_104_MONOMER) RNAs
%   have relatively short half lives compared to that of other macromolecules
%   (eg. protein, DNA) and the M. genitalium cell cycle length. The relatively
%   short half lives of RNAs enables the small M. genitalium with its very small
%   pool of RNAs and particularly mRNAs to sample a broader range of
%   configurations of the RNA pool over a shorter period that would be possible
%   with longer half lifes. This helps the cell more finely tune the expression
%   of proteins, more efficiently execute cell-cycle dependent events, and
%   respond to the external environment. This enhanced fitness due to short RNA
%   half lifes comes at a large energetic cost however.
%
%   In addition to ribonucleases, aminoacylated RNAs require peptidyl tRNA
%   hydrolase (MG_083_MONOMER) to release their conjugated amino acids.
%
%   This process decays all species of RNA, and at all maturation states
%   including aminoacylated states.
%
%   Knowledge Base
%   ===============
%   The knowledge base contains experimentally measured half lifes of many RNA
%   species measured largely in E. coli and mapped to M. genitalium by homology.
%   These half lifes are refined, by simulation.fitConstants to make them
%   consistent with other experimental data used to fit the model. Prior to
%   fitting missing half lifes are imputed either as the average of that of all
%   measured RNA species.
%
%      Type   Avg Half Life (m)
%      ====   =================
%      mRNA   4.5 +/- 2.0
%      rRNA   150
%      sRNA   89
%      tRNA   45
%
%   Representation
%   ===============
%   The substrates, enzymes, and RNAs properties represent the counts of
%   metabolites, ribonuclease R and peptidyl tRNA hydrolase enzymes, and RNAs.
%   This process contains no intermediate representation of RNA degradation; RNA
%   degradation is treated as an all-or-nothing event that either proceeds to
%   complete with a time step or doesn't progress at all.
%
%   decayRates represents the decay rate of each RNA species in seconds.
%   decayRates is informed by experimentally measured RNA half lifes organized
%   in the knowledge base, and fit by simulation.fitConstants. decayReactions
%   represents the metabolites required to decay each RNA species, and the
%   metabolic byproducts of the decay of each RNA species. decayReactions is
%   computed by the knowledge RNA classes based on the sequence, processing, and
%   modifications of each RNA species.
%
%   Initialization
%   ===============
%   All RNAs are initialized to the mature state. This is accomplished by the
%   simulation class initializeState method.
%
%   Simulation
%   ===============
%   This process models RNA decay as an enyzme-dependent poisson process with
%   rate parameter:
%     lambda = RNAs .* decayRates * stepSizeSec
%
%   Algorithm
%   +++++++++++++++
%   1. Stochastically select RNAs to decay based on poission distribution with
%      lambda = RNAs .* decayRates * stepSizeSec
%   2. (Ignore limits to decay posed by availability of metabolite reactants
%      since the only reactant is water, and water is abundantly available)
%   3. Limit RNA decay by available enzyme activity
%      a. All RNAs require ribonuclease R to decay
%      b. Additionally, only decay aminoacylated tRNAs up to the limit of
%         available peptidyl tRNA hydrolase activity.
%   4. Update counts of RNAs
%   5. Update counts of metabolic byproducts of RNA decay
%
% Author: Markus Covert, mcovert@stanford.edu
% Author: Jayodita Sanghvi, jayodita@stanford.edu
% Author: Jonathan Karr, jkarr@stanford.edu
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 7/30/2010
classdef RNADecay < edu.stanford.covert.cell.sim.Process

    %property annotations
    properties (Constant)
        optionNames__              = {}; %names of option properties
        fixedConstantNames__       = {   %names of fixed constant properties
		    'peptidylTRNAHydrolaseSpecificRate';
            'ribonucleaseRFragmentLength';
            'decayReactions';
            };
        fittedConstantNames__      = {}; %names of fitted constant properties
        localStateNames__          = {   %names of simulation state properties redundant with timecourses in this or other processes or the simulation
            'RNAs'};
    end

    %IDs, names, and local indices
    properties
        stimuliWholeCellModelIDs = {};   %whole cell model IDs of stimuli

        substrateWholeCellModelIDs = {}; %whole cell model IDs of substrates
        substrateIndexs_hydrogen         %index within substrates of hydrogen
        substrateIndexs_water            %index within substrates of water
        substrateIndexs_methionine       %index within substrates of methionine
        substrateIndexs_fmethionine      %index within substrates of formylmethionine
        substrateIndexs_glutamate        %index within substrates of glutamate
        substrateIndexs_glutamine        %index within substrates of glutamine
        substrateIndexs_formate          %index within substrates of formate
        substrateIndexs_ammonia          %index within substrates of ammonia
        substrateIndexs_aminoAcids       %index within substrates of amino acids
        substrateIndexs_nmps             %index within substrates of NMPs

        enzymeWholeCellModelIDs = {      %enzyme whole cell model ids
            'MG_104_MONOMER';            %ribonuclease R
            'MG_083_MONOMER'};           %peptidyl-tRNA hydrolase
        enzymeIndexs_ribonucleaseR         = 1; %index within enzymes of ribonuclease R
        enzymeIndexs_peptidylTRNAHydrolase = 2; %index within enzymes of peptidyl-tRNA hydrolase

        matureTRNAIndexs                 %indices of mature tRNAs within RNAs
        matureTMRNAIndexs                %indices of mature tmRNAs within RNAs
    end

    %fixed biological constants
    properties
        peptidylTRNAHydrolaseSpecificRate   %0.700 [PUB_0026]
        ribonucleaseRFragmentLength         %5 [PUB_0039]
        decayReactions                      %adjacency matrix -- RNAs X (reactants and products of RNA decay)
    end

    %global state (copied locally for convenience)
    properties
        RNAs                                %counts of RNAs
    end

    %global state (referenced locally for convenience)
    properties
        transcripts                         %New Transcripts state class
    end

    %constructor
    methods
        function this = RNADecay(wholeCellModelID, name)
            this = this@edu.stanford.covert.cell.sim.Process(wholeCellModelID, name);
        end
    end

    %communication between process/simulation
    methods
        %set references to state objects
        function storeObjectReferences(this, simulation)
            this.storeObjectReferences@edu.stanford.covert.cell.sim.Process(simulation);

            this.transcripts = simulation.state('Transcript');
            this.states = [this.states; {this.transcripts}];
        end

        %initialize constants
        function initializeConstants(this, knowledgeBase, simulation, varargin)
            s = this.rna;
            g = this.gene;

            %include all metabolites involved in RNA decay
            decayReactions_nascentRNA       = knowledgeBase.nascentRNAs.decayReactions;
            decayReactions_processedRNA     = knowledgeBase.processedRNAs.decayReactions;
            decayReactions_intergenicRNA    = knowledgeBase.intergenicRNAs.decayReactions;
            decayReactions_matureRNA        = knowledgeBase.matureRNAs.decayReactions;
            decayReactions_aminoacylatedRNA = knowledgeBase.aminoacylatedRNAs.decayReactions;

            this.substrateWholeCellModelIDs = unique([simulation.state('Metabolite').wholeCellModelIDs(...
                any(decayReactions_nascentRNA,       1) | ...
                any(decayReactions_processedRNA,     1) | ...
                any(decayReactions_intergenicRNA,    1) | ...
                any(decayReactions_matureRNA,        1) | ...
                any(decayReactions_aminoacylatedRNA, 1));
                'H';'H2O';'NH3';'FOR';
                'ALA';'ARG';'ASN';'ASP';'CYS';'GLN';'GLU';'GLY';'HIS';'ILE';'LEU';'LYS';'MET';'PHE';'PRO';'SER';'THR';'TRP';'TYR';'VAL';'FMET']);

            %super class method
            this.initializeConstants@edu.stanford.covert.cell.sim.Process(...
                knowledgeBase, simulation, varargin{:});

            %substrate indices
            this.substrateIndexs_hydrogen    = this.substrateIndexs({'H'});
            this.substrateIndexs_water       = this.substrateIndexs({'H2O'});
            this.substrateIndexs_methionine  = this.substrateIndexs({'MET'});
            this.substrateIndexs_fmethionine = this.substrateIndexs({'FMET'});
            this.substrateIndexs_glutamate   = this.substrateIndexs({'GLU'});
            this.substrateIndexs_glutamine   = this.substrateIndexs({'GLN'});
            this.substrateIndexs_ammonia     = this.substrateIndexs({'NH3'});
            this.substrateIndexs_formate     = this.substrateIndexs({'FOR'});
            this.substrateIndexs_aminoAcids  = this.substrateIndexs({'ALA';'ARG';'ASN';'ASP';'CYS';'GLN';'GLU';'GLY';'HIS';'ILE';'LEU';'LYS';'MET';'PHE';'PRO';'SER';'THR';'TRP';'TYR';'VAL';'FMET'});
            this.substrateIndexs_nmps        = this.substrateIndexs({'AMP';'CMP';'GMP';'UMP'});

            this.matureTRNAIndexs = find(any(s.matureRNAGeneComposition(g.tRNAIndexs, :), 1))';
            this.matureTMRNAIndexs = s.getIndexs('MG_0004');

            this.decayReactions = zeros(numel(this.rna.wholeCellModelIDs), numel(this.substrateWholeCellModelIDs));
            this.decayReactions(s.nascentIndexs, :)       = decayReactions_nascentRNA(:,       this.substrateMetaboliteGlobalIndexs);
            this.decayReactions(s.processedIndexs, :)     = decayReactions_processedRNA(:,     this.substrateMetaboliteGlobalIndexs);
            this.decayReactions(s.intergenicIndexs, :)    = decayReactions_intergenicRNA(:,    this.substrateMetaboliteGlobalIndexs);
            this.decayReactions(s.matureIndexs, :)        = decayReactions_matureRNA(:,        this.substrateMetaboliteGlobalIndexs);
            this.decayReactions(s.boundIndexs, :)         = decayReactions_matureRNA(:,        this.substrateMetaboliteGlobalIndexs);
            this.decayReactions(s.misfoldedIndexs, :)     = decayReactions_matureRNA(:,        this.substrateMetaboliteGlobalIndexs);
            this.decayReactions(s.damagedIndexs, :)       = decayReactions_matureRNA(:,        this.substrateMetaboliteGlobalIndexs);
            this.decayReactions(s.aminoacylatedIndexs, :) = decayReactions_aminoacylatedRNA(:, this.substrateMetaboliteGlobalIndexs);
        end

        %retrieve state from simulation
        function copyFromState(this)
            this.copyFromState@edu.stanford.covert.cell.sim.Process();

            this.RNAs = this.rna.counts(:, this.compartment.cytosolIndexs, :);
        end

        %send state to simulation
        function copyToState(this)
            this.copyToState@edu.stanford.covert.cell.sim.Process();

            this.rna.counts(:, this.compartment.cytosolIndexs, :) = this.RNAs;
        end
    end

    %memory alloction for state
    methods
        function allocateMemoryForState(this, numTimePoints)
            this.allocateMemoryForState@edu.stanford.covert.cell.sim.Process(numTimePoints);

            this.RNAs = zeros(size(this.rna.counts, 1), 1, numTimePoints);
        end
    end

    %model
    methods
        %Calculate
        %- contribution to FBA objective
        %- minimum expression consistent with cell cycle length
        function [bmProd, byProd, minEnzExp, maxEnzExp] = calcResourceRequirements_LifeCycle(this, ~, states)
            import edu.stanford.covert.util.ComputationUtil;
            invMat = this.rna.intergenicRNAMatrix * ...
                ComputationUtil.invertCompositionMatrix(this.rna.nascentRNAMatureRNAComposition);

            %% substrate and byproducts
            %data
            matureRNAReactions = this.decayReactions(this.rna.aminoacylatedIndexs, :); %same as mature for (m,r)RNA and for non-aminoacylated sRNA
            intergenicRNAReactions = this.decayReactions(this.rna.intergenicIndexs, :);
            intergenicRNADecays = invMat * states.rnaProductions;

            %RNA decay
            bmProd = ...
                + max(0, -matureRNAReactions)'     * states.rnaDecays ...
                + max(0, -intergenicRNAReactions)' * intergenicRNADecays;
            byProd = ...
                + max(0,  matureRNAReactions)'     * states.rnaDecays ...
                + max(0,  intergenicRNAReactions)' * intergenicRNADecays;

            %% enzymes

            %data
            rnaDecays = states.rnaDecays0;
            aminoacylatedRNADecays = rnaDecays([this.matureTRNAIndexs; this.matureTMRNAIndexs]);
            intergenicRNADecays = invMat * states.rnaProductions0;

            %RNA decay
            minEnzExp = zeros(size(this.enzymeWholeCellModelIDs));
            minEnzExp(this.enzymeIndexs_ribonucleaseR) = ...
                2 * (sum(rnaDecays) + sum(intergenicRNADecays));
            minEnzExp(this.enzymeIndexs_peptidylTRNAHydrolase) = ...
                2 * sum(aminoacylatedRNADecays) / this.peptidylTRNAHydrolaseSpecificRate;
            maxEnzExp = Inf(size(this.enzymeWholeCellModelIDs));
        end

        %initialization: RNAs intialized to mature/aminoacylated state by
        %simulation intializeState method
        function initializeState(~)
        end

        %resource requirements
        function result = calcResourceRequirements_Current(this)
            if this.enzymes(this.enzymeIndexs_ribonucleaseR) == 0
                result = zeros(size(this.substrates));
                return;
            end

            result = max(0, -this.decayReactions' * (min(1, this.rna.decayRates) .* this.RNAs));
            if ~isempty(this.transcripts.abortedTranscripts)
                result(this.substrateIndexs_water) = ...
                    + result(this.substrateIndexs_water) ...
                    + sum(this.transcripts.abortedTranscripts(:, 2) - 1);
            end
        end

        %simulation
        function evolveState(this)
            % import classes
            import edu.stanford.covert.cell.kb.ssRNA;

            % numbers of enzymes
            ribonucleaseR = this.enzymes(this.enzymeIndexs_ribonucleaseR);
            peptidylTRNAHydrolase = this.randStream.stochasticRound(...
                this.enzymes(this.enzymeIndexs_peptidylTRNAHydrolase) ...
                * this.peptidylTRNAHydrolaseSpecificRate ...
                * this.stepSizeSec);

            %Ribonuclease R required for decay, terminate early if no ribonuclease R
            if ribonucleaseR == 0
                return;
            end

            %% decay all aborted transcripts
            abortedSeqs = this.transcripts.abortedSequences;
            abortedTfs = false(size(abortedSeqs));
            for i = 1:numel(abortedSeqs)
                substrateCost = ssRNA.computeDecayReaction(ssRNA.computeBaseCount(...
                    abortedSeqs{i}, numel(this.substrates), this.substrateIndexs_nmps), ...
                    numel(abortedSeqs{i}), 'linear', ...
                    this.substrateIndexs_water, this.substrateIndexs_hydrogen)';
                if any(this.substrates < -substrateCost)
                    break;
                end
                abortedTfs(i) = true;
                this.substrates = this.substrates + substrateCost;
            end
            this.transcripts.abortedTranscripts = this.transcripts.abortedTranscripts(~abortedTfs, :);

            %% Stochastically decay free RNA as poisson process
            decayingRNAs = min(this.randStream.random('poisson', ...
                this.RNAs .* min(1e6, this.rna.decayRates * this.stepSizeSec)), ...
                this.RNAs);
            if ~any(decayingRNAs)
                return;
            end

            %Require peptidyl tRNA hydrolase to decay aminoacylated tRNAs
            tmp = decayingRNAs(this.rna.aminoacylatedIndexs);
            tmp2 = zeros(size(tmp));
            while any(tmp)
                if peptidylTRNAHydrolase <= 0
                    break;
                end
                idx = this.randStream.randsample(numel(tmp), 1, true, tmp);
                tmp(idx) = tmp(idx) - 1;
                tmp2(idx) = tmp2(idx) + 1;
                peptidylTRNAHydrolase = peptidylTRNAHydrolase - 1;
            end
            decayingRNAs(this.rna.aminoacylatedIndexs) = tmp2;

            %require substrates (water) to decay RNAs
            tmp = decayingRNAs;
            tmp2 = zeros(size(tmp));
            water = this.substrates(this.substrateIndexs_water);
            waterReqs = max(0, -this.decayReactions(:, this.substrateIndexs_water));
            while any(tmp)
                idx = this.randStream.randsample(numel(tmp), 1, true, tmp);
                if water < waterReqs(idx);
                    break;
                end
                water = water - waterReqs(idx);
                tmp(idx) = tmp(idx) - 1;
                tmp2(idx) = tmp2(idx) + 1;
            end
            decayingRNAs = tmp2;

            %update numbers of RNAs
            this.RNAs = this.RNAs - decayingRNAs;

            %update metabolites
            %- water, hydrogen for hydrolysis
            %- nucleotide monophosphate salvage
            this.substrates = this.substrates + this.decayReactions' * decayingRNAs;

            this.rna.counts(:, this.compartment.cytosolIndexs) = this.RNAs;
            if any(any(this.rna.updateExternalState(-decayingRNAs, true)))
                throw(MException('RNADecay:error', 'All RNAs should have been degraded'));
            end
            this.RNAs = this.rna.counts(:, this.compartment.cytosolIndexs);
        end
    end

    %get methods of dependent local state
    methods
        function value = getDryWeight(this)
            if size(this.RNAs, 3) == 1
                value = this.getDryWeight@edu.stanford.covert.cell.sim.Process() + ...
                    this.rna.molecularWeights' * this.RNAs / edu.stanford.covert.util.ConstantUtil.nAvogadro;
            else
                value = this.getDryWeight@edu.stanford.covert.cell.sim.Process() + ...
                    permute(this.rna.molecularWeights' * permute(this.RNAs,[1 3 2]),[1 3 2]) / edu.stanford.covert.util.ConstantUtil.nAvogadro;
            end
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function lambda = tris_eigenvalues ( n, x, y, z )

%*****************************************************************************80
%
%% TRIS_EIGENVALUES returns the eigenvalues of the TRIS matrix.
%
%  Discussion:
%
%    The eigenvalues will be complex if X * Z < 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, real X, Y, Z, the scalars that define A.
%
%    Output, complex LAMBDA(N,1), the eigenvalues.
%
  lambda = zeros ( n, 1 );

  for i = 1 : n
    angle = i * pi / ( n + 1 );
    lambda(i,1) = y + 2.0 * sqrt ( x * z ) * cos ( angle );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%
% Work in progress.  Needs to be written to pay heed to sparsity.
%
function [x] = Back_Solve(A,b)
x = A\b;


text:
function AA = int_tim_seg(tt,L,we,type, we_histo)
%function AA = int_tim_seg(tt,L,we,type, we_histo)
%  to be used with function fast_mr.m in @fast_mr
%  tt is time vector
%  L is number of segments, i.e. number of interpolation points
%  we is fieldmap, not required for linear and hanning types.
%  type is type of interpolator as follows:
%        1   is for ideal min-max using field map
%        2   is for histogram approximation, using we_histo
%  we_histo: only used when type 2 is selected
%            has two columns:
%                    column 1:  bin centers for histogram
%                    column 2:  histogram values at those bin centers

mint = min(tt(:));
maxt = max(tt(:));
rangt = maxt-mint;
minwe = min(we(:));
maxwe = max(we(:));
ndat = length(tt);
N = numel(we);

tau = (rangt+eps)/(L);

AA = zeros(L+1,ndat);

if L==0
  AA = ones(1,ndat);
  return
end

tt = tt-mint;

if (type == 1) % Exact LS interpolator
    gg = exp(i*we(:)*tau)*ones(1,L);  %CHECK MINUS SIGN HERE
    lll = ones(size(gg(:,1)))*[1:L];
    gl = gg.^lll;
    G = [ones(size(gg(:,1))),gl];
    glsum = sum(gl);
    GTG = zeros(L+1,L+1);
    GTG = GTG+diag(N*ones(L+1,1),0);
    for kk = 1:L
       GTG = GTG+diag(glsum(kk)*ones(L+1-kk,1),kk);
       GTG = GTG+diag(conj(glsum(kk))*ones(L+1-kk,1),-kk);
    end
       if (rcond(GTG)>10*eps)
           iGTGGT = inv(GTG)*G';
       else
           iGTGGT = pinv(GTG)*G';
           sprintf('used pinv instead')
      end
 %  for yy = 1:ndat
           cc = exp(i*we(:)*(tt(:))');
          % cc = exp(i*we(:)*tt(yy));
          % AA(:,yy) = (iGTGGT*cc)'.';
           AA = (iGTGGT*cc)'.';
 %   end
end
if (type == 2) % Approx Histogram Interpolator using we_histo
     %find bin size to give an integer value for KK
     rangwe = max(we_histo(:,1))-min(we_histo(:,1));
     minwe = min(we_histo(:,1));
     num_bins = length(we_histo(:,1));
     KK = floor(2*pi/((rangwe/num_bins)*tau));
     % bin size
     dwn = 2*pi/(KK*tau);
     bin_centers = we_histo(:,1);
     N_ap = we_histo(:,2).';
     %keyboard
     ftwe_ap = fft(N_ap.*(exp(i*dwn*[0:(num_bins-1)]*tau*(L))),KK);
     ftwe_ap = exp(-i*(minwe+dwn/2)*tau*([0:(KK-1)]-(L))).*ftwe_ap;
     GTGap_ap = zeros(L+1,L+1);
     for kk = 1:(2*L+1)
        GTGap_ap = GTGap_ap+diag(ftwe_ap(kk)*ones(L+1-abs((L+1)-kk),1),-(L+1)+kk);
     end
       if (rcond(GTGap_ap.')>10*eps)
           iGTGap_ap = inv(GTGap_ap.');
       else
           iGTGap_ap = pinv(GTGap_ap.');
           sprintf('used pinv instead')
      end
     for yy = 1:ndat
        ftc_ap= fft(N_ap.*(exp(i*[0:(num_bins-1)]*dwn*tt(yy))),KK);
        ftc_ap = exp(i*(minwe+dwn/2)*(tt(yy)-tau*[0:KK-1])).*(ftc_ap);
        GTc_ap = ftc_ap(1:L+1).';
        AA(:,yy) = (iGTGap_ap*GTc_ap)'.';
     end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [textures] = getNGTDMtextures(NGTDM,countValid)
% -------------------------------------------------------------------------
% function [textures] = getNGTDMtextures(NGTDM,countValid)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes texture features from an input Neighborhood
% Gray-Tone Difference Matrix (NGTDM).
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Amadasun, M., & King, R. (1989). Textural Features Corresponding to
%     Textural Properties. IEEE Transactions on Systems Man and Cybernetics,
%     19(5), 1264–1274.
% -------------------------------------------------------------------------
% INPUTS:
% - NGTDM: Neighborhood Gray-Tone Difference Matrix.
% - countValid: Number of valid voxels used in the NGTDM computation.
%               Required for the computation of texture features in
%               'getNGTDMtextures.m'
%
% ** 'NGTDM' and 'countValid' should be outputs from 'getNGTDM.m' **
% -------------------------------------------------------------------------
% OUTPUTS:
% - textures: Struture specifying the values of different NGTDM texture
%             features as defined below.
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: January 2013
% - Revision: May 2015
% -------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package 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 package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

% PRELIMINARY
nTot = sum(countValid);
countValid = countValid./nTot; % Now representing the probability of gray-level occurences
NL = length(NGTDM);
Ng = sum(countValid~=0);
pValid = find(countValid>0);
nValid = length(pValid);

% COMPUTATION OF TEXTURES
% 1. Coarseness, Ref.[1]
textures.Coarseness = (((countValid')*NGTDM) + eps)^(-1);

% 2. Contrast, Ref.[1]
val = 0;
for i = 1:NL
    for j = 1:NL
        val = val + countValid(i)*countValid(j)*(i-j)^2;
    end
end
textures.Contrast = val*sum(NGTDM)/(Ng*(Ng-1)*nTot);

% 3. Busyness, Ref.[1]
denom = 0;
for i = 1:nValid
    for j = 1:nValid
        denom = denom + abs(pValid(i)*countValid(pValid(i))-pValid(j)*countValid(pValid(j)));
    end
end
textures.Busyness = ((countValid')*NGTDM)/denom;

% 4. Complexity, Ref.[1]
val = 0;
for i = 1:nValid
    for j = 1:nValid
        val = val + (abs(pValid(i)-pValid(j))/(nTot*(countValid(pValid(i)) + countValid(pValid(j)))))*(countValid(pValid(i))*NGTDM(pValid(i)) + countValid(pValid(j))*NGTDM(pValid(j)));
    end
end
textures.Complexity = val;

% 5. Strength, Ref.[1]
val = 0;
for i = 1:nValid
    for j = 1:nValid
        val = val + (countValid(pValid(i))+countValid(pValid(j)))*(pValid(i)-pValid(j))^2;
    end
end
textures.Strength = val/(eps+sum(NGTDM));

end

author: mvallieres

repo: radiomics


text:
% Astar
% A* navigation class
%
% A concrete subclass of the Navigation class that implements the A*
% navigation algorithm. Methods included are for the standard case,
% multiobjective optimization (MOO) -- i.e. optimizes over several
% objectives/criteria -- and the A*-PO algorithms for MOO that utilizes
% Pareto optimality.
%
% Methods:
% 	plan            Compute the cost map given a goal and map
% 	path            Compute a path to the goal
% 	visualize       Display the obstacle map (deprecated)
% 	plot            Display the obstacle map
% 	costmap_modify 	Modify the costmap
% 	costmap_get     Return the current costmap
% 	costmap_set     Set the current costmap
% 	display         Print the parameters in human readable form
% 	char            Convert to string
%
% Properties:
% TBD
%
% Example 1::
%        load map1              % load map
%        goal = [50;30];
%        start=[20;10];
%        as = Astar(map);       % create Navigation object
%        as.plan(goal,2,3,0);   % setup costmap for specified goal;
%                               % standard D* algorithm w/ 2 objectives
%                               % and 3 costmap layers
%        as.path(start);        % plan solution path start-to-goal, animate
%        P = as.path(start);    % plan solution path start-to-goal, return
%                               % path
% Example 2::
%        goal = [100;100];
%        start = [1;1];
%        as = Astar(0);          % create Navigation object with pseudo-
%                                % random occupancy grid
%        ds.addCost(terrain);    % terrain is a 100x100 matrix of
%                                % elevations [0,1]
% 	     ds.plan(goal,3,4,0);    % setup costmap for specified goal
%                                % (3 and 4 include the added terrain cost)
%        as.path(start);         % plan solution path start-goal, animate
%        P = as.path(start);     % plan solution path start-goal, return
%                                % path
%
% Notes
% - Obstacles are represented by Inf in the costmap.
%
% References
% - A Pareto Optimal D* Search Algorithm for Multiobjective Path Planning,
%   A. Lavin.
% - A Pareto Front-Based Multiobjective Path Planning Algorithm, A. Lavin.
% - Robotics, Vision & Control, Sec 5.2.2, Peter Corke, Springer, 2011.
%
% See Also Navigation, Dstar

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

% Implementation notes:
%
% X is an index into the array of states.
% State pointers are kept as matlab array index rather than row,col format.

classdef Astar < Navigation

    properties (SetAccess=private, GetAccess=private)

        % essential world info
        costmap   % cost layers (1st is map)
        G         % index of goal point
        N         % number of objectives
        L         % number of cost layers

        % info kept per cell (state):
        b         % backpointer (0 means not set)
        t         % tag: NEW/OPEN/CLOSED

        algorithm % A*, A*-MOO, or A*-PO
        tie
        openlist  % priority queue with states and their costs
        niter
        changed
        openlist_maxlen
        quiet     % specifies verbosity

        % tag state values
        NEW = 0;
        OPEN = 1;
        CLOSED = 2;
    end

    methods  % start of public methods

        function as = Astar(world, varargin)
            %Astar.Astar A* constructor
            %
            % AS = Astar(MAP, OPTIONS) is a A* navigation object, and MAP
            % is an occupancy grid, a representation of a planar world as
            % a matrix whose elements are 0 (free space) or 1 (occupied).
            % The occupancy grid is coverted to a costmap with a unit cost
            % for traversing a cell.
            %
            % Options::
            % 'world' = 0   will call for a pseudo-random occupancy grid
            % 'goal',G      Specify the goal point (2x1)
            % 'metric',M    Specify the distance metric as 'Euclidean'
            %               (default) or 'cityblock'
            % 'inflate',K   Inflate all obstacles by K cells
            % 'quiet'       Don't display the progress spinner
            %
            % Other options are supported by the Navigation superclass.
            %
            % See also Navigation.Navigation.

            % Invoke the superclass constructor
            as = as@Navigation(world, varargin{:});  % includes the occgrid

            % Options
            opt.quiet = false;
            opt = tb_optparse(opt, varargin);
            as.quiet = opt.quiet;

            as.occgrid2costmap(as.occgrid);

            % Initialize A* state variables
            as.reset();
            if ~isempty(as.goal)
                as.goal_change();
            end
            as.changed = false;
        end

        function reset(as)
            %Astar.reset Reset the planner
            %
            % AS.reset() resets the A* planner.  The next instantiation
            % of AS.plan() will perform a global replan.

            % Build the matrices required to hold the state of each cell
            as.b = zeros(size(as.occgrid), 'uint32');     % backpointers
            as.t = zeros(size(as.occgrid), 'uint8');      % tags, all NEW=0
            as.costmap(:,:,2) = zeros(size(as.occgrid));  % path cost g
            as.costmap(:,:,3) = zeros(size(as.occgrid));  % path cost h

            % Priority queue has col for each open state, one row for the
            % state location, and rows for each cost layer
            as.openlist = zeros(as.L+1,0);
            as.openlist_maxlen = -Inf;
        end

        function goal_change(as)
            %Astar.goal_change Changes the costlayers due to new goal
            %position
            if isempty(as.b)
                return;
            end
            goal = as.goal;

            % Keep goal in index rather than row,col format
            as.G = sub2ind(size(as.occgrid), goal(2), goal(1));
            as.INSERT(as.G, as.projectCost(as.G), 'goalset');
            as.costmap(goal(2),goal(1),2) = 0;

            % If new goal modifies costs for a layer, recalculate here
            as.calcHeuristic(as.occgrid, as.goal);
        end

        function s = char(as)
            %Astar.char Convert Navigation object to string
            %
            % AS.char() is a string representing the state of the Astar
            % object in human-readable form.
            %
            % See also Astar.display, Navigation.char.

            % Work is done by the superclass
            s = char@Navigation(as);
        end

        function plot(as, varargin)
            %Astar.plot Visualize navigation environment
            %
            % AS.plot() displays the occupancy grid and the goal distance
            % in a new figure.  The goal distance is shown by intensity
            % which increases with distance from the goal.  Obstacles are
            % overlaid and shown in red.
            %
            % AS.plot(P) as above but also overlays a path given by the set
            % of points P (Mx2).
            %
            % See also Navigation.plot.

            plot@Navigation(as, 'distance', as.costmap(:,:,3), varargin{:});
        end

        function n = next(as, current)
            % Invoked by Navigation.step
            % Backpropagate from goal to start
            % Return [col;row] of previous step
            if as.changed
                error('Cost map has changed, replan');
            end
            X = sub2ind(size(as.occgrid), current(2), current(1));
            % Set X as the backpointer of X
            X = as.b(X);
            if X == 0
                % Goal (no further backpointer)
                n = [];
            else
                [r,c] = ind2sub(size(as.occgrid), X);
                n = [c;r];
            end
        end

        function plan(as, goal, N, layers, algorithm)
            %Astar.plan Prep the grid for planning.
            %
            % AS.plan() updates AS with a costmap of distance to the
            % goal from every non-obstacle point in the map.  The goal is
            % as specified to the constructor.
            %
            % Inputs:
            %   goal: goal state coordinates
            %   N: number of optimization objectives; standard A* is 2
            %   (i.e. distance and heuristic)
            %   layers: number of cost layers in costmap
            %   algorithm: specify standard A*(0), A*-MOO (1), A*-PO (2)

            % Setup parameters
            if nargin < 3
                N = 2;
            end
            if nargin < 4
                layers = 3;
            end
            if nargin < 5
                algorithm = 0;
            end
            as.N = N;
            as.L = layers;
            as.algorithm = algorithm;
            as.openlist = zeros(as.L+1,0);

            % Initialize cost layers
            for a = 2:as.L
                as.costmap(:,:,a) = zeros(size(as.occgrid));
            end

            % Cost priority/tiebreaker: layer 2 (distance to node)
            as.tie = 2;

            % Set goal
            if nargin > 1
                as.goal = goal;  % invokes superclass method set.goal()
            end
            if isempty(as.goal)
                error('must specify a goal point');
            end

            % Populate heuristic cost layer
            as.calcHeuristic(as.occgrid, as.goal);
        end

        function P = path(as, start)
        %Astar.path Find a path between two points
        %
        % AS.path(START) finds and displays a path from START to GOAL
        % which is overlaid on the occupancy grid.
        %
        % P = AS.path(START) returns the path (2xM) from START to GOAL.
            if nargin < 1
                error('must specify start state');
            end

            % Invoke the superclass path function, which iterates on our
            % next method
%             start = [start; 1];  % specifies backpropogation for NAV.path()
%             temp = start;
%             start = as.goal;
%             as.goal = temp;

            if nargout == 0
                path@Navigation(as, start);
            else
                P = path@Navigation(as, start);
            end
        end

        % Handler invoked by Navigation.path() to start the navigation
        % process -- calculate the solution path.
        % Line comments Ln reference A* pseudocode in Lavin's "A Pareto
        % Front-Based Multiobjective Path Planning Algorithm" where n is
        % the line number.
        function navigate_init(as, start)
            as.openlist = zeros(as.L+1,0);  % openlist must be empty
            % Begin search with the start node
            start = sub2ind(size(as.occgrid), start(2), start(1));
            as.openlist(1,1) = start;
            as.t(start) = as.OPEN;

            % Plan the A* path
            as.niter = 0; flag = 0;
            while ~isempty(as.openlist)                              % L4
                % Normalize costs on the open list, choose expansion state
                queue = normc(as.openlist(2:size(as.openlist,1),:)');
                if as.algorithm == 2
                    % Get Pareto optimal point off the open list
                    front = as.openlist(:,paretofront(queue));
                    [~,col] = min(front(as.tie+1,:));
                    X = front(1,col);                                % L5
                else
                    [~,ind]=min(sum(queue,2));
                    X = as.openlist(1,ind);                          % L5
                end
                    as.DELETE(X);                                    % L6

                as.niter = as.niter + 1;
                if ~as.quiet && mod(as.niter, 20) == 0
                    as.spinner();
                end

                % Populate the openlist
                for Y=as.neighbors(X)                                % L7,8
                    if(Y==as.G)                                      % L9
                        as.b(Y) = X;
                        as.updateCosts(Y,X,as.N)
                        flag = 1;  % flag for goal
                        break;
                    end
                    if as.t(Y)==as.NEW && as.costmap(Y)~=Inf
                        as.b(Y) = X;
                        as.updateCosts(Y,X,as.N);
                        % Project node's costs into objective space:
                        objspace = as.projectCost(Y,X);
                        as.INSERT(Y, objspace, '');
                    end
                end
                if as.verbose
                    disp(' ')
                end
                if flag==1 % goal found
                    break;
                end
            end
            if ~as.quiet
                fprintf('\r');
            end
            as.changed = false;
        end

        function layer = cost_get(as, layer)
        %Astar.cost_get Get the specified cost layer
            layer = as.costmap(:,:,layer);
        end

        function c = heurstic_get(as)
        %Astar.heuristice_get Get the current heuristic map
        %
        % C = AS.heuristice_get() is the current heuristic layer. It is
        % computed in Astar.plan.
        %
        % See also Astar.plan.
            c = as.costmap(:,:,3);
        end

        function c = costmap_get(as)
        %Astar.costmap_get Get the current costmap
        %
        % C = AS.costmap_get() is the current costmap.
        % The value of each element represents the cost of traversing the
        % cell.  It is autogenerated by the class constructor from the
        % occupancy grid such that:
        % - free cell (occupancy 0) has a cost of 1
        % - occupied cell (occupancy >0) has a cost of Inf
        %
        % See also Astar.costmap_set, Astar.costmap_modify.
            c = as.costmap;
        end

        function costmap_set(as, costmap)
        %Astar.costmap_set Set the current costmap
        %
        % AS.costmap_set(C) sets the current costmap.
        % This method accepts the full costmap -- i.e. all layers.
        %
        % Notes:
        % - After the cost map is changed the path should be replanned by
        %   calling AS.plan().
        %
        % See also Astar.costmap_get, Astar.costmap_modify.
            [i,j,k] = size(costmap);
            if ~all([i,j] == size(as.occgrid))
                error('costmap must be same size as occupancy grid');
            end
            as.L = k;  % set the number of cost layers
            as.costmap = costmap;
            as.changed = true;
        end

        function costmap_modify(as, point, newcost)
        %Astar.costmap_modify Modify cost map
        %
        % AS.costmap_modify(P, NEW) modifies the cost map at P=[X,Y] to
        % have the value NEW.  If P (2xM) and NEW (1xM) then the cost of
        % the points defined by the columns of P are set to the corresponding
        % elements of NEW.
        %
        % Notes::
        % - After one or more point costs have been updated the path
        %   should be replanned by calling AS.plan().
        %
        % See also Astar.costmap_set, Astar.costmap_get.
            if (newcost < 0) || (1 < newcost)
                error('new cost value must be normlaized [0,1]')
            end

            [i,j,k] = size(as.costmap);
            if (point(1) < 0) || (point(1) > i)
                error('1st dimension of point is out of bounds')
            end
            if (point(2) < 0) || (point(2) > j)
                error('2nd dimension of point is out of bounds')
            end
            if (point(2) < 0) || (point(3) > k)
                error('3rd dimension of point is out of bounds')
            end

            as.costmap(point) = newcost;
        end

        function addCost(as, values)
        %Astar.addCost Add an additional cost layer
        %
        % AS.addCost(values) adds the matrix specified by values as a
        % cost layer.
        % Inputs
        %   values: normalized matrix the size of the environment
            [i,j,k] = size(as.costmap);

            if [i,j]~=size(as.occgrid)
                error('layer size does not match the environment')
            end
            if max(max(values))~=1 || min(min(values))~=0
                error('layer values are not normalized [0,1]')
            end

            as.costmap(:,:,k+1) = values;
        end

        function flag = backProp(as)
            flag = 1;
        end

    end  % end of public methods

    methods (Access=protected)  % start of private methods

        function occgrid2costmap(as, og, cost)
            if nargin < 3
                cost = 1;
            end
            og(og==1) = Inf;  % occupied cells -> infinite path cost
            og(og==0) = cost;  % unoccupied cells -> path cost
            as.costmap(:,:,1) = og;
        end

        function calcHeuristic(as, grid, goal)
            as.costmap(:,:,3) = zeros(size(grid));
            for ii=1:size(grid,1)
                for jj=1:size(grid,2)
                    as.costmap(ii,jj,3) = sqrt((ii-goal(1))^2+(jj-goal(2))^2);
                end
            end
        end

        function k_new = updateCosts(as, a, b, obj)
            % NOTE: Only for costs that accumulate (i.e. sum) over the
            % path, and for dynamic costs.
            % E.g. the heuristic parameter only needs updating when the
            % goal state changes; its values are stored for each cell.
            %
            % Location moving from state b to a.
            %
            % The costs are coded to be (1) distance, (2) heuristic, (3)
            % elevation, (4) solar deviation, and (5) risk. If deviating
            % from these costs (in this order) you MUST EDIT THIS METHOD.
            [i,j,~] = size(as.costmap);

            if nargout > 0
                k_new = as.costmap(i*j+b) + as.dc(b,a);
                return
            end
            if obj == 0
                % Return what the new priority cost would be (k_new)
                return
            end
            if obj > 1
                % Standard A* search
                as.costmap(i*j+a) = as.costmap(i*j+b) + as.dc(b,a);
                % (no heuristic update needed)
            end
            if obj > 2
                % W/ elevation costs
                % (no elevation update needed)
            end
            if obj > 3
                % W/ solar costs
                % Rotate the solar vector 1rad per 100 steps
                sV = [cos(as.niter/100);sin(as.niter/100)];
                as.costmap(4*i*j+a) = dot(sV,as.vc(b,a));
            end
            if obj > 4
                % W/ risk costs
                % (no risk update needed)
            end
        end

        function pt = projectCost(as, a, b)
            % Returns the projection of state a into objective space. If
            % specified, location is moving from b to a (case 3).
            [i,j,k] = size(as.costmap);
            pt(1) = as.costmap(a);
            switch nargin
                case 2
                    pt(2) = as.costmap(i*j+a);
                case 3
                    pt(2) = as.costmap(i*j+b) + as.dc(a,b);
                otherwise
                    return
            end
            for n=3:k
                pt(n) = as.costmap((n-1)*i*j+a);
            end
        end

        function INSERT(as, X, pt, where)
            % Add state X to the openlist with objective space values
            % specified by pt.

            if nargin>2
                as.message('insert (%s) %d = %f\n', where, X, pt);
            end

            i = find(as.openlist(1,:) == X);
            if length(i) > 1
                error('A*:INSERT: state in open list %d times', X);
            end

            [i,j,~] = size(as.costmap);
            if (as.t(X) == as.OPEN || as.CLOSED) && ...
               (pt(as.tie) > as.costmap((as.tie-1)*i*j+X))
                % L13/14: If a node with same position as successor is in
                % the OPEN/CLOSED list & has a lower f than successor,
                % then skip this successor.
            else
                % Add a new column to the open list for this node
                as.openlist = [as.openlist [X; pt(:)]];
            end

            % Keep track of the max length of the openlist
            if numcols(as.openlist) > as.openlist_maxlen
                as.openlist_maxlen = numcols(as.openlist);
            end

            % Tag state X as open
            as.t(X) = as.OPEN;
        end

        function DELETE(as, X)
            as.message('delete %d\n', X);
            i = find(as.openlist(1,:) == X);
            if length(i) ~= 1
                error('A*:DELETE: state %d does not exist', X);
            end

            % Remove the column, close the state
            as.openlist(:,i) = [];
            as.t(X) = as.CLOSED;
        end

        function cost = dc(as, X, Y)
            % Return the distance cost of moving from state X to state Y
            [r,c] = ind2sub(size(as.occgrid), [X; Y]);
            dist = sqrt(sum(diff([r c]).^2));
            dcost = (as.costmap(X) + as.costmap(Y))/2;

            cost = dist * dcost;
        end

        function vector = vc(as, X, Y)
            % Return the robot unit vector -- direction of moving from
            % state X to state Y
            [Xi,Xj] = ind2sub(size(as.occgrid),X);
            [Yi,Yj] = ind2sub(size(as.occgrid),Y);
            vector = [Yi-Xi;Yj-Xj];
            vector = vector/norm(vector);
        end

        function Y = neighbors(as, X)
            % Return indices of neighbor states (max 8) as a row vector
            dims = size(as.occgrid);
            [r,c] = ind2sub(dims, X);

            % Of 8-way neighbors, only use those w/in grid bounds
            Y = [r-1 r-1 r-1 r r  r+1 r+1 r+1; c-1 c c+1 c-1 c+1 c-1 c c+1];
            k = (min(Y)>0) & (Y(1,:)<=dims(1)) & (Y(2,:)<=dims(2));
            Y = Y(:,k);
            Y = sub2ind(dims, Y(1,:)', Y(2,:)')';
        end

    end  % end of private methods

end

author: petercorke

repo: robotics-toolbox-matlab


text:
function gT = simwhiteKernDiagGradX(kern, t)

% SIMWHITEKERNDIAGGRADX Gradient of SIM-WHITE kernel's diagonal w.r.t. t.
% FORMAT
% DESC computes the gradient of the diagonal of the SIM-White (Single Input
% Motif - White) kernel matrix with respect to the elements of the input
% column vector given in t.
% ARG kern : the kernel structure for which gradients are being computed.
% ARG t : the input data in the form of a design matrix.
% RETURN gT : the gradients of the diagonal with respect to each element
% of t. The returned matrix has the same dimensions as t.
%
% SEEALSO : simwhiteKernParamInit, kernDiagGradX, simwhitekernGradX
%
% COPYRIGHT : David Luengo, 2009

% KERN

if (kern.isStationary == true)
    gT = zeros(size(t));
else
    gT = kern.variance * (kern.sensitivity^2) * exp(-2*kern.decay*t);
end


text:
function test_mu(varargin)
%
% demonstration file for NMFLibrary.
%
% This file illustrates how to use this library.
%
% This file is part of NMFLibrary.
%
% Created by H.Kasai on OCt. 27, 2017
%

    if nargin < 1
        clc;
        clear;
        close all;
        rng('default')

        m = 500;
        n = 100;
        V = rand(m,n);
        rank = 20;
        options = [];
        options.verbose = 2;
        options.max_epoch = 100;
        health_check_mode = false;
    else
        V = varargin{1};
        rank = varargin{2};
        options = varargin{3};
        health_check_mode = true;
    end

    % initialize factor matrices
    [x_init, ~] = generate_init_factors(V, rank, []);
    options.x_init = x_init;

    %% MU variants
    % MU
    options.alg = 'mu';
    [w_nmf_mu, infos_nmf_mu] = fro_mu_nmf(V, rank, options);

    % MU mod
    options.alg = 'mu_mod';
    [w_nmf_mu_mod, infos_nmf_mu_mod] = fro_mu_nmf(V, rank, options);

    % Accelerated MU
    options.alg = 'mu_acc';
    [w_nmf_mu_acc, infos_nmf_mu_acc] = fro_mu_nmf(V, rank, options);

    if ~health_check_mode
        %% plot
        display_graph('iter','cost', {'MU', 'MU-Mod', 'MU-ACC'}, {w_nmf_mu, w_nmf_mu_mod, w_nmf_mu_acc},  {infos_nmf_mu, infos_nmf_mu_mod, infos_nmf_mu_acc});
        display_graph('time','cost', {'MU', 'MU-Mod', 'MU-ACC'}, {w_nmf_mu, w_nmf_mu_mod, w_nmf_mu_acc},  {infos_nmf_mu, infos_nmf_mu_mod, infos_nmf_mu_acc});
    end

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
% This function can read in all the information from the HTK model file
% with specifications of AURORA2 baseline system.

% AUTHOR: XIAO XIONG
% CREATED: 7 Jul, 2006
% LAST MODIFIED: 7 Jul, 2006
function model = readGMMs(file_name)

FID = fopen(file_name, 'r');
if FID < 1
    fprintf('Failed to open the model file\n');
    return;
end

tmp = textscan(FID,'%s%d\n',1);
Dim = tmp{2};
fsearch('~v "varFloor1"',FID); fgetl(FID);fgetl(FID);
tmp = textscan(FID,'%n',Dim);
model.varFloor = tmp{1};
fgetl(FID);

modelCnt = 1;
while(fsearch('~g',FID)==0)
    tmp = textscan(FID,'%s',1);
    model.gmm{modelCnt}.name = regexprep(tmp{1}{1},'"','');
    fgetl(FID);fgetl(FID);
    [model.gmm{modelCnt}.prior, model.gmm{modelCnt}.mean, model.gmm{modelCnt}.var] = read_GMM(FID,1);
    fgetl(FID);
    modelCnt = modelCnt+1;
end

fclose(FID);

author: singaxiong

repo: SignalGraph


text:
function [params, names] = kbrExtractParam(model,dim);

% KBREXTRACTPARAM Extract parameters from the KBR model structure.
% FORMAT
% DESC extracts parameters from the kernel based regression model
% structure into a vector of parameters for optimisation.
% ARG model : the model structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the model.
%
% DESC extracts parameters and parameter names from the kernel based regression model structure.
% ARG model : the model structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the model.
% RETURN names : cell array of strings containing names for each parameter.
%
%
%
% SEEALSO : KBRCREATE, KBREXPANDPARAM, MODELEXTRACTPARAM, SCG, CONJGRAD
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006, 2008
%
% MODIFICATIONS : Carl Henrik Ek, 2007

% MLTOOLS

if(nargin<2)
  params = [model.A(:)' model.bias];
else
  params = model.A(:,dim);
  params = [params(:)' model.bias(dim)];
end

if nargout > 1
  % Add names to parameters
  counter = 0;
  for i = 1:model.numData
    for j = 1:model.outputDim
      counter = counter + 1;
      names{counter} = ['A(' num2str(i) ', ' num2str(j) ')'];
    end
  end
  for j = 1:model.outputDim
    counter = counter + 1;
    names{counter} = ['bias(' num2str(j) ')'];
  end
end


text:
function [R, nt] = buildr(S, varargin)

%BUILDR Builds response matrix for input to functions ENTROPY and
% INFORMATION
%
%   ------
%   SYNTAX
%   ------
%
%       [R, nt] = build_R(S, R1, R2, ..., RL)
%
%   -----------
%   DESCRIPTION
%   -----------
%   Let's consider an experiment in which, during each trial, a stimulus is
%   presented out of NS available stimuli and L distinct neural responses
%   are recorded simultaneously. Within this framework the input to BUILDR
%   is as follow:
%   - S is the stimulus-array, i.e., S(i) stores the value of the stimulus
%     presented during the i-th trial.
%   - Rj, j=1,...,L, is the j-th response array, i.e., Rj(i) stores the
%     value of the j-th response recorded during the i-th experiment.
%
%   The function will return the response matrix R which can be input to
%   the routines ENTROPY and INFORMATION. The routine also outputs the
%   trials per stimulus array, NT, which can be used as one of the
%   parameters of the option structure for ENTROPY and INFORMATION.
%
%   -------
%   REMARKS
%   -------
%
%   - Although the one described above is the framework which was kept in
%     mind while creating the toolbox functions, it has to be noted that
%     this function (and also the other in the toolbox) can be easily
%     applied to several other situations.
%
%   - The goal of this function is to help the user get familiar with the
%     structure of the response-matrix R which is input to the ENTROPY and
%     INFORMATION functions and, in particular, with the fact that the
%     stimulus values are not provided to these functions. This is because,
%     for the sake of information computation the only important parameter
%     concerning the stimulus is the number of times each stimulus was
%     presented: the value of each stimulus can thus be mapped to an
%     integer value and made implicit as the index to a page of the
%     response matrix R. It needs to be noted, however, that BUILDR, having
%     to be as generic as possible, will also be relatively slow. Often,
%     the way responses are recorded or computed allows building a response
%     matrix much more quickly than the way done by BUILR. For
%     computationally intensive tasks it is thus suggested that custom
%     routines are used instead of this built-in tool.
%
%   See also ENTROPY, INFORMATION

%   Copyright (C) 2009 Cesare Magri
%   Version: 1.0.4

% -------
% LICENSE
% -------
% This software is distributed free under the condition that:
%
% 1. it shall not be incorporated in software that is subsequently sold;
%
% 2. the authorship of the software shall be acknowledged and the following
%    article shall be properly cited in any publication that uses results
%    generated by the software:
%
%      Magri C, Whittingstall K, Singh V, Logothetis NK, Panzeri S: A
%      toolbox for the fast information analysis of multiple-site LFP, EEG
%      and spike train recordings. BMC Neuroscience 2009 10(1):81;
%
% 3.  this notice shall remain in place in each source file.

if ~isvector(S)
    error('Stimulus array must be a 1-D array');
end

totNt = length(S);

L = length(varargin);

R1toL = zeros(L, totNt);
for k=1:L
    if ~isvector(varargin{k})
        msg = 'Response arrays must be 1-D.';
        error('buildr:respNot1D', msg);
    end

    if length(varargin{k}) ~= totNt
        msg = 'Each response-array must have the same length as the stimulus array';
        error('buildr:differentTotNt', msg);
    end

    R1toL(k,:) = varargin{k};
end

uniqueS = unique(S);
Ns = length(uniqueS);

% Dispalying informations:
disp('Building R and nt:');
disp(['- number of stimuli = ' num2str(Ns)]);
disp(['- number of responses = ' num2str(L)]);

nt = zeros(Ns,1);
tFlag = false(totNt, Ns);
for s=1:Ns
    tFlag(:,s) = S==uniqueS(s);
	nt(s) = sum(tFlag(:,s));
end

maxNt = max(nt);
disp(['- maximum numer of trials = ' num2str(maxNt)]);
disp(['- minimum numer of trials = ' num2str(min(nt))]);
R = zeros(L, maxNt, Ns);
for s=1:Ns
    if nt(s)>0
        R(:,1:nt(s), s) = R1toL(:, tFlag(:,s));
    else
        msg = 'One or more stimuli with no corresponding response.';
        error('buildr:noResponseStimulus', msg);
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip