text:
% Compiles all mex routines that are part of toolbox.
%
% If you get warnings about "omp.h", your compiler may lack default OpenMp
% support (for info about OpenMP see http://en.wikipedia.org/wiki/OpenMP).
% To compile without OpenMP alter code below to 'useOmp(:)=0'; and re-run
% toolboxCompile. Note that this will disable parallelization and make some
% routines (in particular training certain classifier) much slower.
%
% USAGE
%  toolboxCompile
%
% INPUTS
%
% OUTPUTS
%
% EXAMPLE
%
% See also
%
% Piotr's Computer Vision Matlab Toolbox      Version 3.30
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

% compile options including openmp support for C++ files
opts = {'-output'};
if(exist('OCTAVE_VERSION','builtin')), opts={'-o'}; end
if( ispc ), optsOmp={'OPTIMFLAGS="$OPTIMFLAGS','/openmp"'}; else
  optsOmp={'CXXFLAGS="\$CXXFLAGS','-fopenmp"'};
  optsOmp=[optsOmp,'LDFLAGS="\$LDFLAGS','-fopenmp"'];
end
optsOmp=[optsOmp '-DUSEOMP'];

% list of files (missing /private/ part of directory)
fs={'channels/convConst.cpp', 'channels/gradientMex.cpp',...
  'channels/imPadMex.cpp', 'channels/imResampleMex.cpp',...
  'channels/rgbConvertMex.cpp', 'classify/binaryTreeTrain1.cpp', ...
  'classify/fernsInds1.c', 'classify/forestFindThr.cpp',...
  'classify/forestInds.cpp', 'classify/meanShift1.c',...
  'detector/acfDetect1.cpp', 'images/assignToBins1.c',...
  'images/histc2c.c', 'images/imtransform2_c.c', ...
  'images/nlfiltersep_max.c', 'images/nlfiltersep_sum.c', ...
  'videos/ktComputeW_c.c', 'videos/ktHistcRgb_c.c', ...
  'videos/opticalFlowHsMex.cpp' };
n=length(fs); useOmp=zeros(1,n); if(~ismac), useOmp([6 9])=1; end

% compile every funciton in turn (special case for dijkstra)
disp('Compiling Piotr''s Toolbox.......................');
rd=fileparts(mfilename('fullpath')); rd=rd(1:end-9); tic;
errmsg=' -> COMPILE FAILURE: ''%s'' %s\n';
for i=1:n
  try %#ok<ALIGN>
    [d,f1,e]=fileparts(fs{i}); f=[rd '/' d '/private/' f1];
    if(useOmp(i)), optsi=[optsOmp opts]; else optsi=opts; end
    fprintf(' -> %s\n',[f e]); mex([f e],optsi{:},[f '.' mexext]);
  catch err, fprintf(errmsg,[f1 e],err.message); end
end
try %#ok<ALIGN>
  d=[rd '/matlab/private/']; fprintf(' -> %s\n',[d 'dijkstra1.cpp']);
  mex([d 'fibheap.cpp'], [d 'dijkstra1.cpp'], '-largeArrayDims', ...
    opts{:}, [d 'dijkstra1.' mexext]);
catch err, fprintf(errmsg,[f1 e],err.message); end
disp('..................................Done Compiling'); toc;

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function colorbarDataIndex = findColorbarData(obj, colorbarIndex, colorbarData)

    if nargin == 2

        %locate index of data associated with colorbar
        colorbarDataIndex = find( ...
            arrayfun( @(x)eq(x.AssociatedAxis, ...
                            obj.State.Colorbar(colorbarIndex).AssociatedAxis), ...
                            obj.State.Plot ...
                    ), ...
            1 );

        %if no matching data index found
        if isempty(colorbarDataIndex)
            colorbarDataIndex = max(min(colorbarIndex,obj.State.Figure.NumPlots),1);
        end

    elseif nargin == 3
        c = 1; a = 1;
        allAxesIndex = zeros(length(colorbarData.Parent.Children), 1);

        for n = 1:length(colorbarData.Parent.Children)
            if strcmp(colorbarData.Parent.Children(n).Type, 'colorbar')
                allColorbarIndex(c) = n;
                c = c + 1;
            elseif strcmp(colorbarData.Parent.Children(n).Type, 'axes')
                allAxesIndex(n) = a;
                a = a + 1;
            end
        end

        colorbarAxisIndex = allColorbarIndex(colorbarIndex) + 1;
        colorbarAxisIndex = allAxesIndex(colorbarAxisIndex);
        colorbarDataIndex = obj.State.Figure.NumAxes - colorbarAxisIndex + 1;

        colorbarDataIndex = colorbarDataIndex;
    end

end

author: plotly


text:
function x = log_series_cdf_inv ( cdf, a )

%*****************************************************************************80
%
%% LOG_SERIES_CDF_INV inverts the Logarithmic Series CDF.
%
%  Discussion:
%
%    Simple summation is used.  The only protection against an
%    infinite loop caused by roundoff is that X cannot be larger
%    than 1000.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 December 1999
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real CDF, the value of the CDF.
%
%    Input, real A, the parameter of the PDF.
%    0.0 < A < 1.0.
%
%    Output, real X, the argument of the CDF for which
%    CDF(X-1) <= CDF <= CDF(X).
%
  xmax = 1000;

  if ( cdf < 0.0 | 1.0 < cdf )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LOG_SERIES_CDF_INV - Fatal error!\n' );
    fprintf ( 1, '  CDF < 0 or 1 < CDF.\n' );
    error ( 'LOG_SERIES_CDF_INV - Fatal error!' );
  end

  cdf2 = 0.0;
  x = 1;

  while ( cdf2 < cdf & x < xmax )

    if ( x == 1 )
      pdf = - a / log ( 1.0 - a );
    else
      pdf = ( x - 1 ) * a * pdf / x;
    end

    cdf2 = cdf2 + pdf;

    x = x + 1;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function colorSegmentationSol ( config_name )

%
%  Create a replicated array of colors
%
  colors = {'red', 'yellow', 'green', 'magenta'};
%
%  Read the image into a replicated array
%
  image = imread('fabric.png');
%
%  Set up the matlabpool
%
if matlabpool('size') > 0; matlabpool('close'); end
%
% there are four colors, each lab filters one
%
if nargin < 1 || ~strcmp(config_name, 'ithaca');
    matlabpool('open', 'local', 4);
else
    matlabpool('open', 'ithaca_2009b', 4, ...
          'FileDependencies', {'colorSegmentationSol.m', 'colorFilter.m'});
end

%
%  Each lab filters a different color,
%  and writes the filtered result to a JPEG file.
%
  spmd

    filteredImage = colorFilter(image, colors{labindex});

    imwrite(filteredImage, ['fab' num2str(labindex) '.jpg']);

  end

  return
end


text:
function dt = CurvedCNSdt2D(Q, gamma, mu)

% function dt = CurvedCNSdt2D(Q, gamma, mu)
% Purpose: compute stable time step size for compressible Navier-Stokes solver

Globals2D;

% extract field variables
rho = Q(:,:,1); rhou = Q(:,:,2); rhov = Q(:,:,3); Ener = Q(:,:,4);

% evaluate fields at surface nodes
rho = rho(vmapM); rhou = rhou(vmapM); rhov = rhov(vmapM); Ener = Ener(vmapM);

% compute primitive variables
u = rhou./rho; v = rhov./rho;
p = (gamma-1.0)*(Ener - rho.*(u.^2+v.^2)/2);
c = sqrt(abs(gamma*p./rho));

h = 2./Fscale(:);
lam = sqrt ( u(:).^2 + v(:).^2 ) + c(:);
dt = 0.5*min(1./( (N+1)^2*lam./h + (N+1)^4*mu./(h.^2)))
return

repo: Aero-matlab


text:
function s = cvx_where

%CVX_WHERE    Returns the location of the CVX system.
%   CVX_WHERE returns a string containing the base directory of the CVX
%   modeling framework. Within that directory are some useful
%   subdirectories and files:
%       functions/    new functions
%       examples/     sample cvx models
%       LICENSE.txt   copyright information
%   The proper operation of this function assumes that it has not been
%   moved from its default position within the cvx distribution.

try
    s = dbstack('-completenames');
catch
    s = dbstack;
end
s = s(1);
if isfield( s, 'file' ),
    s = s.file;
else
    s = s.name;
end
if ispc,
    fs = '\';
else
    fs = '/';
end
temp = strfind( s, fs );
s( temp(end-1) + 1 : end ) = [];

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [ a_lu, info ] = r8pbu_fa ( n, mu, a )

%*****************************************************************************80
%
%% R8PBU_FA factors a R8PBU matrix.
%
%  Discussion:
%
%    The R8PBU storage format is for a symmetric positive definite band matrix.
%
%    To save storage, only the diagonal and upper triangle of A is stored,
%    in a compact diagonal format that preserves columns.
%
%    The diagonal is stored in row MU+1 of the array.
%    The first superdiagonal in row MU, columns 2 through N.
%    The second superdiagonal in row MU-1, columns 3 through N.
%    The MU-th superdiagonal in row 1, columns MU+1 through N.
%
%    The matrix A must be a positive definite symmetric band matrix.
%
%    Once factored, linear systems A*x=b involving the matrix can be solved
%    by calling R8PBU_SL.  No pivoting is performed.  Pivoting is not necessary
%    for positive definite symmetric matrices.  If the matrix is not positive
%    definite, the algorithm may behave correctly, but it is also possible
%    that an illegal divide by zero will occur.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 October 1998
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Bunch, Moler, Stewart,
%    LINPACK User's Guide,
%    SIAM, Philadelphia, 1979.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, integer MU, the number of superdiagonals of the matrix.
%    MU must be at least 0, and no more than N-1.
%
%    Input, real A(MU+1,N), the N by N matrix, stored in LINPACK
%    positive definite symmetric band matrix storage.
%
%    Output, real A_LU(MU+1,N), information describing a factored form
%    of the matrix, that can be used to solve linear systems
%    A*x=b, using R8PBU_SL.
%
%    Output, integer INFO, singularity flag.
%    0, the matrix is nonsingular.
%    nonzero, the matrix is singular.
%
  info = 0;
  a_lu(1:mu+1,1:n) = a(1:mu+1,1:n);

  for j = 1 : n

    ik = mu + 1;
    jk = max ( j - mu, 1 );
    mm = max ( mu + 2 - j, 1 );

    s = 0.0;

    for k = mm : mu

      a_lu(k,j) = ( a_lu(k,j) - a_lu(ik:ik+k-mm-1,jk)' * a_lu(mm:k-1,j) )...
        / a_lu(mu+1,jk);

      s = s + a_lu(k,j) * a_lu(k,j);

      ik = ik - 1;
      jk = jk + 1;

    end

    s = a_lu(mu+1,j) - s;

    if ( s <= 0.0 )
      info = j;
      fprintf ( 1, '\n' );
      fprintf ( 1, 'R8PBU_FA - Fatal error!\n' );
      fprintf ( 1, '  Nonpositive pivot on step %d\n', info );
      return;
    end

    a_lu(mu+1,j) = sqrt ( s );

  end

  return
end


text:
function [SE_MR,SE_RZF,SE_MMMSE,SE_ZF,SE_SMMSE] = functionComputeSE_DL_hardening(H,Hhat,C,R,tau_c,tau_p,nbrOfRealizations,M,K,L,p,rho)
%Compute DL SE for different transmit precoding schemes using Theorem 4.6.
%
%INPUT:
%H                 = M x nbrOfRealizations x K x L x L matrix with the
%                    exact channel realizations
%Hhat              = M x nbrOfRealizations x K x L x L matrix with the
%                    channel estimates
%C                 = M x M x K x L x L matrix with estimation error
%                    correlation matrices
%R                 = M x M x K x L x L matrix with spatial correlation
%                    matrices
%tau_c             = Length of coherence block
%tau_p             = Length of pilot sequences
%nbrOfRealizations = Number of channel realizations
%M                 = Number of antennas per BS
%K                 = Number of UEs per cell
%L                 = Number of BSs and cells
%p                 = Uplink transmit power per UE (same for everyone)
%rho               = Downlink transmit power per UE (same for everyone)
%
%OUTPUT:
%SE_MR    = K x L x length(tau_c) matrix where element (k,l,:) is the
%           downlink SE of UE k in cell l achieved with MR precoding for
%           different lengths of the coherence block
%SE_RZF   = Same as SE_MR but with RZF precoding
%SE_MMMSE = Same as SE_MR but with M-MMSE precoding
%SE_ZF    = Same as SE_MR but with ZF precoding
%SE_SMMSE = Same as SE_MR but with S-MMSE precoding
%
%
%This Matlab function was developed to generate simulation results to:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.01 (Last edited: 2020-05-15)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Store identity matrices of different sizes
eyeK = eye(K);
eyeM = eye(M);

if nargout > 2

    %Compute sum of all estimation error correlation matrices at every BS
    C_totM = reshape(p*sum(sum(C,3),4),[M M L]);

end

%Compute sum of intra-cell estimation error correlation matrices at every BS
if nargout>4

    CR_totS = zeros(M,M,L);

    for j = 1:L
        CR_totS(:,:,j) = p*(sum(C(:,:,:,j,j),3)+sum(sum(R(:,:,:,[1:j-1 j+1:end],j),3),4));
    end

end

%Compute the prelog factor assuming only downlink transmission
prelogFactor = (tau_c-tau_p)./(tau_c);

%Prepare to store simulation results for signal gains
signal_MR = zeros(K,L);

if nargout > 1
    signal_RZF = zeros(K,L);
end

if nargout > 2
    signal_MMMSE = zeros(K,L);
end

if nargout > 3
    signal_ZF = zeros(K,L);
end

if nargout > 4
    signal_SMMSE = zeros(K,L);
end

%Prepare to store simulation results for sum interference powers
interf_MR = zeros(K,L);

if nargout > 1
    interf_RZF = zeros(K,L);
end

if nargout > 2
    interf_MMMSE = zeros(K,L);
end

if nargout > 3
    interf_ZF = zeros(K,L);
end

if nargout > 4
    interf_SMMSE = zeros(K,L);
end

%% Go through all channel realizations
for n = 1:nbrOfRealizations

    %Go through all cells
    for j = 1:L

        %Extract channel realizations from all UEs to BS j
        Hallj = reshape(H(:,n,:,:,j),[M K*L]);

        %Extract channel realizations from all UEs to BS j
        Hhatallj = reshape(Hhat(:,n,:,:,j),[M K*L]);

        %Compute MR combining in (4.11)
        V_MR = Hhatallj(:,K*(j-1)+1:K*j);

        if nargout > 1 %Compute RZF combining in (4.9)
            V_RZF = (p*V_MR)/(p*(V_MR'*V_MR)+eyeK);
        end

        if nargout > 2 %Compute M-MMSE combining in (4.7)
            V_MMMSE = (p*(Hhatallj*Hhatallj')+C_totM(:,:,j)+eyeM)\(p*V_MR);
        end

        if nargout > 3 %Compute ZF combining in (4.10), with the small regularization term 1e-12 for numerical stability
            V_ZF = V_MR/(V_MR'*V_MR+1e-12*eyeK);
        end

        if nargout > 4 %Compute S-MMSE combining in (4.8)
            V_SMMSE = (p*(V_MR*V_MR')+CR_totS(:,:,j)+eyeM)\(p*V_MR);
        end

        %Go through all UEs in cell j
        for k = 1:K

            if norm(V_MR(:,k))>0

                %%MR precoding in (4.37) computed based on combining vector
                w = V_MR(:,k)/norm(V_MR(:,k));

                %Compute realizations of the terms inside the expectations
                %of the signal and interference terms of (4.26)
                signal_MR(k,j) = signal_MR(k,j) + (w'*H(:,n,k,j,j))/nbrOfRealizations;
                interf_MR = interf_MR + rho*reshape(abs(w'*Hallj).^2,[K L])/nbrOfRealizations;

                %%RZF precoding
                if nargout > 1

                    %RZF precoding in (4.37) computed based on combining vector
                    w = V_RZF(:,k)/norm(V_RZF(:,k));

                    %Compute realizations of the terms inside the expectations
                    %of the signal and interference terms of (4.26)
                    signal_RZF(k,j) = signal_RZF(k,j) + (w'*H(:,n,k,j,j))/nbrOfRealizations;
                    interf_RZF = interf_RZF + rho*reshape(abs(w'*Hallj).^2,[K L])/nbrOfRealizations;

                end

                %%M-MMSE precoding
                if nargout > 2

                    %M-MMSE precoding in (4.37) computed based on combining vector
                    w = V_MMMSE(:,k)/norm(V_MMMSE(:,k));

                    %Compute realizations of the terms inside the expectations
                    %of the signal and interference terms of (4.26)
                    signal_MMMSE(k,j) = signal_MMMSE(k,j) + (w'*H(:,n,k,j,j))/nbrOfRealizations;
                    interf_MMMSE = interf_MMMSE + rho*reshape(abs(w'*Hallj).^2,[K L])/nbrOfRealizations;

                end

                %%ZF precoding
                if nargout > 3

                    %ZF precoding in (4.37) computed based on combining vector
                    w = V_ZF(:,k)/norm(V_ZF(:,k));

                    %Compute realizations of the terms inside the expectations
                    %of the signal and interference terms of (4.26)
                    signal_ZF(k,j) = signal_ZF(k,j) + (w'*H(:,n,k,j,j))/nbrOfRealizations;
                    interf_ZF = interf_ZF + rho*reshape(abs(w'*Hallj).^2,[K L])/nbrOfRealizations;

                end

                %%S-MMSE precoding
                if nargout > 4

                    %S-MMSE precoding in (4.37) computed based on combining vector
                    w = V_SMMSE(:,k)/norm(V_SMMSE(:,k));

                    %Compute realizations of the terms inside the expectations
                    %of the signal and interference terms of (4.26)
                    signal_SMMSE(k,j) = signal_SMMSE(k,j) + (w'*H(:,n,k,j,j))/nbrOfRealizations;
                    interf_SMMSE = interf_SMMSE + rho*reshape(abs(w'*Hallj).^2,[K L])/nbrOfRealizations;

                end

            end

        end

    end

end

%% Prepare to compute the SEs with different pilot reuse factors
SE_MR = zeros(K,L,length(tau_c));

if nargout > 1
    SE_RZF = zeros(K,L,length(tau_c));
end

if nargout > 2
    SE_MMMSE = zeros(K,L,length(tau_c));
end

if nargout > 3
    SE_ZF = zeros(K,L,length(tau_c));
end

if nargout > 4
    SE_SMMSE = zeros(K,L,length(tau_c));
end

%Go through all lengths of the pilot reuse factors
for n = 1:length(tau_c)

    %Compute SEs according to Theorem 4.6
    SE_MR(:,:,n) = prelogFactor(n)*real(log2(1+(rho*abs(signal_MR).^2) ./ (interf_MR - rho*abs(signal_MR).^2 + 1)));

    if nargout > 1
        SE_RZF(:,:,n) = prelogFactor(n)*real(log2(1+(rho*abs(signal_RZF).^2) ./ (interf_RZF - rho*abs(signal_RZF).^2 + 1)));
    end

    if nargout > 2
        SE_MMMSE(:,:,n) = prelogFactor(n)*real(log2(1+(rho*abs(signal_MMMSE).^2) ./ (interf_MMMSE - rho*abs(signal_MMMSE).^2 +1)));
    end

    if nargout > 3
        SE_ZF(:,:,n) = prelogFactor(n)*real(log2(1+(rho*abs(signal_ZF).^2) ./ (interf_ZF - rho*abs(signal_ZF).^2 + 1)));
    end

    if nargout > 4
        SE_SMMSE(:,:,n) = prelogFactor(n)*real(log2(1+(rho*abs(signal_SMMSE).^2) ./ (interf_SMMSE - rho*abs(signal_SMMSE).^2 + 1)));
    end

end

author: emilbjornson

repo: massivemimobook


text:
function [su_hdr,ex_hdr,se_hdr,im_hdr,pix_hdr,im_offset] = GE_readHeader(IFileName)
%
% [su_hdr,ex_hdr,se_hdr,im_hdr,pix_hdr,im_offset] = GE_readHeader(IFileName)
% Reads the header info from a GE lx2 or 5.X file
%
%
% Souheil J. Inati
% Dartmouth College
% May 2000
% souheil.inati@dartmouth.edu
%

%%%% Open the IFile %%%%%
[fid,message] = fopen(IFileName,'r','b');         %note: 'b' = Big-Endian format
if fid == -1
  error(message)
end

% Check for the magic number at the first word
magic = fread(fid,1,'int32');
if magic == 1229801286
  % This is a 5.X format image
  byte_align = 0;
  pix_hdr_offset = 0;
  su_hdr_size = 114;
  ex_hdr_size = 1024;
  se_hdr_size = 1020;
  im_hdr_size = 1022;
else
  % Magic no. not at the 1st word, try at the LX2 position
  fseek(fid,3228,-1);
  magic = fread(fid,1,'int32');
  if magic == 1229801286
    % This is an LX2 format image
    byte_align = 1;
    pix_hdr_offset = 3228;
    su_hdr_size = 116;
    ex_hdr_size = 1040;
    se_hdr_size = 1028;
    im_hdr_size = 1044;
  else
    msg = 'This is not a 5.X or LX2 format image.  No Magic Number.';
    error(msg);
  end
end

% Load the pixel header
fseek(fid,pix_hdr_offset,-1);
pix_hdr = GE_readHeaderPixel(fid, byte_align);

% Compute the offsets
su_hdr_offset = pix_hdr.img_p_dbHdr;
ex_hdr_offset = su_hdr_offset + su_hdr_size;
se_hdr_offset = ex_hdr_offset + ex_hdr_size;
im_hdr_offset = se_hdr_offset + se_hdr_size;
im_offset = pix_hdr_offset + pix_hdr.img_hdr_length;

% Check for epirecon images
if (pix_hdr.img_l_dbHdr == 0 & byte_align==0)
  error('This is a epirecon image. No header')
end

% Load the suite header
fseek(fid,su_hdr_offset,-1);
su_hdr = GE_readHeaderSuite(fid, byte_align);

% Load the exam header
fseek(fid,ex_hdr_offset,-1);
ex_hdr = GE_readHeaderExam(fid, byte_align);

% Load the series header
fseek(fid,se_hdr_offset,-1);
se_hdr = GE_readHeaderSeries(fid, byte_align);

% Load the image header
fseek(fid,im_hdr_offset,-1);
im_hdr = GE_readHeaderImage(fid, byte_align);

% Close the file
fclose(fid);

return

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [calibration_matrix,corr_vector,fraction_vector,dates] = aps_spectrometer_PWV_meris_modis
% aps_spectrometer_PWV_meris_modis(batchfile)
% Scipt to load meris and modis PWV data and compare each SAR data for the
% cloud free pixels.
%
% Bekaert David -- University of Leeds
%
% modifications:
% DB    08/2014     Initial codings
% DB    03/2015     Add bootstrapping to get the uncertancy of the linear
%                   mapping between PWV of MERIS and MODIS.

save_dir= 'aps_m_M';
n_boot_runs = 1200;

curdir = pwd;

fontsize = 15;

% loading the date information
meris_datapath = getparm_aps('meris_datapath',1);
modis_datapath = getparm_aps('modis_datapath',1);
stamps_processed = getparm_aps('stamps_processed',1);
if strcmp(stamps_processed,'y')
   ps = load(getparm_aps('ll_matfile',1));
   dates = ps.day;
   fprintf('Stamps processed structure \n')
else
    ifgday_matfile = getparm_aps('ifgday_matfile',1);
    dates = load(ifgday_matfile);
    dates = dates.ifgday;
    dates = reshape(dates,[],1);
    dates = unique(dates);
end

%start loop here to calculate atmos correction for each date
ndates = length(dates);
fprintf('Starting the computation for each SAR date \n')
calibration_matrix = NaN([ndates 2]);            % this is the linear correction matrix [a b] between MERIS = a*MODIS+b (PWV);
corr_vector = NaN([ndates 1]);
fraction_vector =NaN([ndates 1]);
std_vector=NaN([ndates 2]);
calibration_matrix_boot = NaN([n_boot_runs 2 ndates]);

if exist([curdir filesep 'figures'],'dir')~=7
    mkdir([curdir filesep 'figures']);
end

for n = 1:ndates
    date_str = datestr(dates(n,1),'yyyymmdd');

    meris_filename_PWV_nointerp = [meris_datapath filesep date_str filesep date_str '_ZPWV_nointerp.xyz'];
    modis_filename_PWV_nointerp = [modis_datapath filesep date_str filesep date_str '_ZPWV_nointerp.xyz'];

    % checking if there is actual meris data for this date, if not just
    % leave NaN's in the matrix.
    if exist(meris_filename_PWV_nointerp,'file') ==2 && exist(modis_filename_PWV_nointerp,'file') ==2

        [xyz_meris_PWV,temp] = load_meris_SAR(meris_filename_PWV_nointerp);
        [xyz_modis_PWV,temp] = load_meris_SAR(modis_filename_PWV_nointerp);
        clear temp

        % keep only those points in common
        lonlat_modis = xyz_modis_PWV(:,1:2);
        lonlat_meris = xyz_meris_PWV(:,1:2);

        % number of original grid points, this also includes water!
        n_original = size(lonlat_meris,1);

        [lonlat, ix_modis, ix_meris ]= intersect(lonlat_modis,lonlat_meris,'rows','stable');

        PWV_modis=xyz_modis_PWV(ix_modis,3);
        PWV_meris=xyz_meris_PWV(ix_meris,3);

        ix_nan = sum(isnan([PWV_modis PWV_meris]),2)>0;
        lonlat(ix_nan,:)=[];
        PWV_modis(ix_nan,:)=[];
        PWV_meris(ix_nan,:)=[];

        if length(PWV_meris)>100
            % number of gridpoitns left
            n_used = size(lonlat,1);
            fraction_vector(n,1) = n_used./n_original;

            % getting the starts
            A = [PWV_modis  ones(size(PWV_meris)) ];
            calibration_matrix(n,:) = lscov(A,PWV_meris);
            corr_vector(n,1) = corr(PWV_modis,PWV_meris);

            % getting the uncertaincy through bootstrapping - Can just use
            % the powerlaw bootstrap function as its not technique dependent
            [std_vector(n,:),calibration_matrix_boot(:,:,n)]= aps_powerlaw_bootstrap(A,PWV_meris,n_boot_runs);

            limits = [floor(min([PWV_modis;PWV_meris]))   ceil(max([PWV_modis;PWV_meris]))];
            h1 = figure('name',['MODIS/MERIS PWV comparison ' date_str]);
            plot(limits, limits,'k--')
            hold on
            plot(limits,calibration_matrix(n,1).*limits+calibration_matrix(n,2),'r-','linewidth',2)
            hold on
            plot(PWV_modis,PWV_meris,'k.')
            xlim(limits)
            ylim(limits)
            hold on
            plot(limits,calibration_matrix(n,1).*limits+calibration_matrix(n,2),'r-','linewidth',2)
            hold on
            plot(limits, limits,'k--')
            xlabel('MODIS PWV','fontsize',fontsize)
            ylabel('MERIS PWV','fontsize',fontsize)
            title({[date_str ':'],['PWV_{MERIS} = ' num2str(round(calibration_matrix(n,1)*100)/100) ' [+-' num2str(round(std_vector(n,1)*100)/100)  '] *PWV_{MODIS} + ' num2str(round(calibration_matrix(n,2)*100)/100) ' [+-' num2str(round(std_vector(n,2)*100)/100) ']' ]},'fontsize',fontsize)
            legend('1-to-1 relation','Linear fit',0)
            set(gca,'fontsize',fontsize)
            print(h1,'-depsc',[curdir filesep 'figures' filesep 'PWV_MERIS_MODIS_' date_str '.eps'])
            print(h1,'-dpng',[curdir filesep 'figures' filesep 'PWV_MERIS_MODIS_' date_str '.png'])
            close(h1)
        end
        fprintf([num2str(n) ' completed out of ' num2str(ndates) '\n'])
    else
        fprintf([date_str ': no comparison\n'])
    end

end

save([modis_datapath filesep 'MODIS_calibration.mat'],'calibration_matrix','corr_vector','fraction_vector','dates','std_vector','calibration_matrix_boot');

author: dbekaert

repo: TRAIN

save_path: github-repos/MATLAB/dbekaert-TRAIN


text:
function events = in_events_oebin(sFile, EventFile)
% IN_EVENTS_OEBIN: Import events from a Open Ephys flat binary event file (timestamps.npy)
%
% USAGE:  events = in_events_oebin(sFile, EventFile)

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2020

% Read time stamps
evtTime = reshape(readNPY(EventFile), 1, []);
if isempty(evtTime)
    events = [];
    return
end

% Find other files, indicating the type of event
[evtDir, fName] = bst_fileparts(EventFile);
TextFile = bst_fullfile(evtDir, 'text.npy');
% ChanFile = bst_fullfile(evtDir, 'channels.npy');
ChanStateFile = bst_fullfile(evtDir, 'channel_states.npy');

% Get event labels
if file_exist(TextFile)
    disp('EOBIN> ERROR: Text events are not supported yet...');
    evtGroupLabel = {'TEXT'};
    evtGroupInd = {1:length(evtTime)};
    evtChan = [];

%     evtLabels = readNPY(TextFile);
%     % Create event groups
%     evtGroupLabel = unique(evtLabels);
%     evtGroupInd = cell(1, length(evtGroupLabel));
%     for iUnique = 1:length(evtGroupLabel)
%         evtGroupInd{iUnique} = find(strcmpi(evtLabels, evtGroupLabel{iUnique}));
%     end
%     % Read event channel
%     if file_exist(ChanFile)
%         evtChan = readNPY(ChanFile);
%     else
%         evtChan = [];
%     end
elseif file_exist(ChanStateFile)
    evtVal = readNPY(ChanStateFile);
    % Create event groups
    uniqueVal = unique(evtVal);
    evtGroupLabel = cell(1, length(uniqueVal));
    evtGroupInd = cell(1, length(uniqueVal));
    for iUnique = 1:length(uniqueVal)
        evtGroupLabel{iUnique} = num2str(uniqueVal(iUnique));
        evtGroupInd{iUnique} = find(evtVal == uniqueVal(iUnique));
    end
    % Do not get read channels
    evtChan = [];
else
    evtGroupLabel = {'Unknown'};
    evtGroupInd = {1:length(evtTime)};
    evtChan = [];
end

% Initialize list of events
events = repmat(db_template('event'), [1, length(evtGroupLabel)]);
% Get occurrences for each event
for iEvt = 1:length(evtGroupLabel)
    events(iEvt).label      = evtGroupLabel{iEvt};
    events(iEvt).times      = double(evtTime(evtGroupInd{iEvt})) ./ sFile.prop.sfreq;
    events(iEvt).epochs     = ones(1, length(events(iEvt).times));  % Epoch: set as 1 for all the occurrences
    events(iEvt).reactTimes = [];
    events(iEvt).select     = 1;
    events(iEvt).notes      = [];
    events(iEvt).channels   = [];
end

author: brainstorm-tools


text:
%% Transient Behavior and Location of Characteristic Roots in the Complex Plane
% Copyright 2011 The MathWorks, Inc.
%% Initial Value Problem
%
% Consider the system described by the following second-order,
% ordinary differential equation and associated initial conditions:
%
% $$\frac{d^2x}{dt^2} + 2\zeta\omega_n\frac{dx}{dt} + \omega_n^2x = 0; x(0)=x0,\frac{dx}{dt}(0)=0 $$
%
% We would like to examine the behavior of the free response of the
% system as the damping ratio changes.  This will be done first by solving
% the initial value problem for various damping ratios and plotting the
% results.  Then, we will explore changes in the damping ratio affect the
% characteristic roots.  (Recall that the characteristic roots dictate the
% free-response behavior.)

%% Set physical and derived parameters.
%
% To enable numerical solution of the differential equation, we will set
% specific values for the mass and (linear) spring rate.  We will then use
% these to derive the (circular) natural frequency.
%
param.m = 1;                            % Mass [kg]
param.k = 1;                            % Spring rate [N/m]
param.wn = sqrt(param.k/param.m);       % (Circular) natural frequency [rad/s]

%% Set initial conditions.
%
% To solve our 2nd order differential equation, we need two initial
% conditions: one on position, the other on velocity.  Moreover, recall
% from the accompanying PowerPoint slides that we have represented our
% 2nd order linear ODE into a system of 1st order ODE's amenable to
% solution using MATLAB. Having defined position and velocity as our
% system states, we must define an initial state vector.

x0 = 1;              % Position [m]
x0_dot = 0;          % Velocity [m/s]
z0 = [x0 ; x0_dot];  % Initial condition on state vector

%% Set simulation parameters.
%
% We will define the duration of our simulation to be approximately 4
% periods.  To accommodate changes in physical parameters, we will
% parameterize the duration by the natural frequency.  We then define the
% associated time step and vector of simulation time.  Note that the time
% step dt serves to define the times at which output will be provided, not
% the solver time step.

tend = 4*round((2*pi)/param.wn);    % Simulate for approximately 4 periods [s]
dt = round((2*pi)/param.wn)/1000;   % Time step [s]
tsim = 0:dt:tend;                   % Vector of simulation times [s]

%% Undamped system.
%
% We first consider an undamped system, characterized by a damping ratio of
% zero.  Physically, this represents a system with zero dissipation of
% energy.  We will plot the free response to confirm the expected behavior.
%
% To solve the system of 1st order ODE's, we will use the baseline ode45
% solver.  The solver requires a vector of state derivatives; in this
% case, we call upon a function to define them.

param.zeta = 0;     % Damping ratio

[t,z] = ode45(@(t,z) odesys( t,z,param ),tsim,z0);
x(:,1) = z(:,1);    % Save solution for future plotting.

plot(t,x(:,1),'linewidth',2);grid on;hold all

    ylim([-1 1])
    xlabel('Time [s]');
    ylabel('Position [m]')
    title('Free Response: Undamped 2nd Order System');

%% Consider several damping ratios.
%
% Now we will consider four distinct values of the damping ratio
% corresponding to four distinct types of system behavior: undamped,
% underdamped, critically damped, and overdamped.  These terms refer to the
% behavior of the system as it returns to its equilibrium state.  We will
% plot and compare the free responses for these cases.

zetavals = [0 1/sqrt(2) 1 1.5]; % Damping ratios to consider.

case1 = 'Undamped: zeta = 0' ;
case2 = 'Underdamped: zeta = 0.707' ;
case3 = 'Critically Damped: zeta = 1' ;
case4 = 'Overdamped: zeta = 1.5' ;

for i = 2: length(zetavals)

    param.zeta = zetavals(i); % Set current damping ratio.
    [t,z] = ode45(@(t,z) odesys( t,z,param ),tsim,z0); % Compute free response.
    x(:,i) = z(:,1); % Save solution for future plotting.

end;

figure      % Plot results

    plot(t,x,'linewidth',2);grid on
    xlabel('Time [s]');
    ylabel('Position [m]')
    title('Free Response of 2nd Order System for Various Damping Ratios');
    legend(case1,case2,case3,case4,0)

%% Track characteristic roots as damping ratio increases.
%
% Recalling that the free response of a system is governed by its
% characteristic roots, we will examine the changes in these roots as the
% damping ratio changes over a more finely spaced set of values.  Since the
% roots may have real and imaginary parts, we will plot them in the complex
% plane.  For each value of the damping ratio, we will compute the
% characteristic roots as the eigenvalues of the associated state matrix.

zeta_range = 0:0.01:1.5;          % Damping ratio: undamped to overdamped.
S = zeros(2,length(zeta_range));  % Initialize characteristic roots array.

for i = 1:length(zeta_range)

    S(:,i) = eig([0 1; -param.wn^2 -2*zeta_range(i)*param.wn]); % Characteristic roots.

end;

figure      % Plot characteristic roots in complex plane.

    plot(S,'bx','linewidth',2,'markersize',5); grid on;
    set(gca,'YAxisLocation','right')
    xlabel('Real Axis','Position',[-1.5 -0.1 1])
    ylabel('Imaginary Axis')
    title('Loci of Characteristic Roots in Complex Plane for Various Damping Ratios')

%% Show the "odesys" function.
%
% Here we record the function used to define the state matrix
% corresponding to our system of ordinary differential equations.

type odesys.m

author: Sable

repo: mcbench-benchmarks


text:
%% mff_valid.m
%  Matlab File
%  author Colin Davey
%  date 12/3/2013
%  Copyright 2013 EGI. All rights reserved.
%  Support routine for MFF Matlab code. Not intended to be called directly.
%
%  Tests whether a file is valid, and throws an exception if not. Gives an
%  informational warning if filename doesn't end in '.mff'.
%%
function mff_valid(filePath)
valid = false;
filePathExist = exist(filePath);
% Check that filePath exists.
if filePathExist == 0
    theException = MException('EGI_MFF:MFF_NO_EXIST', 'MFF does not exist.');
% Check that filePath isn't a regular file.
elseif filePathExist == 2
    theException = MException('EGI_MFF:MFF_NOT_DIR', 'MFF is not a folder or package.');
% Check that filePath is a directory.
elseif filePathExist ~= 7
    theException = MException('EGI_MFF:MFF_NOT_DIR_OR_FILE', 'MFF is not a folder, package or file.');
else
    % Check that filePath/info.xml exists and is a file.
    if exist([filePath filesep 'info.xml']) ~= 2
        theException = MException('EGI_MFF:INVALID_NO_INFO', 'MFF is not valid. There is no info.xml file.');
    % Check that filePath/info1.xml exists and is a file.
    elseif exist([filePath filesep 'info1.xml']) ~= 2
        theException = MException('EGI_MFF:INVALID_NO_INFO1', 'MFF is not valid. There is no info1.xml file.');
    % Check that filePath/signal1.bin exists and is a file.
    elseif exist([filePath filesep 'signal1.bin']) ~= 2
        theException = MException('EGI_MFF:INVALID_NO_SIGNAL1', 'MFF is not valid. There is no signal1.bin file.');
    % Check the version.
    else
        infoObj = mff_getObject(com.egi.services.mff.api.MFFResourceType.kMFF_RT_Info, 'info.xml', filePath);
        ver = infoObj.getMFFVersion;
        if ver ~= 3;
            theException = MException('EGI_MFF:WRONG_VER', 'MFF is not version 3. Please convert using EGI''s MFF File Converter. Contact supportteam@egi.com form more information.');
        else
            valid = true;
        end
    end
end
if ~valid
    throw(theException);
end
% Delete DS_Store file that gets generated when moving to a PC.
fullDS_StorePath = [filePath filesep '._.DSStore'];
if exist(fullDS_StorePath, 'file') == 2
    delete(fullDS_StorePath);
end

% Give a warning if name doesn't end in '.mff'.
mff_warning = '*** Filename does not end in ''.mff''. This can cause a problem for some versions of NetStation. ***';
if size(filePath,2) < size('.mff',2) + 1
%     warning('EGI_MFF:EXT_WARNING', mff_warning);
    fprintf('%s\n', mff_warning);
elseif ~(strcmp(lower(filePath(end-3:end)), '.mff'))
%     warning('EGI_MFF:EXT_WARNING', mff_warning);
    fprintf('%s\n', mff_warning);
end


text:
function [s,error] = openDEeph( s, denum, ineph, dontSave )
% open a binary JPL planetary ephemeris file and sets initial values
% short int ephem_open (char *ephem_name,
%                       double *jd_begin, double *jd_end,
%                       short int *de_number)
%
% ------------------------------------------------------------------------
%
%    PURPOSE:
%       This function opens a JPL planetary ephemeris file and
%       sets initial values.  This function must be called
%       prior to calls to the other JPL ephemeris functions.
%
%    REFERENCES:
%       Standish, E.M. and Newhall, X X (1988). "The JPL Export
%          Planetary Ephemeris"; JPL document dated 17 June 1988.
%
%    INPUT
%    ARGUMENTS:
%       *ephem_name (char)
%          Name of the direct-access ephemeris file.
%
%    OUTPUT
%    ARGUMENTS:
%       *jd_begin (double)
%          Beginning Julian date of the ephemeris file.
%       *jd_end (double)
%          Ending Julian date of the ephemeris file.
%       *de_number (short int)
%          DE number of the ephemeris file opened.
%
%    RETURNED
%    VALUE:
%       (short int)
%           0   ...file exists and is opened correctly.
%           1   ...file does not exist/not found.
%           2-10...error reading from file header.
%           11  ...unable to set record length; ephemeris (DE number)
%                  not in look-up table.
%
%    GLOBALS
%    USED:
%       SS                eph_manager.h
%       JPLAU             eph_manager.h
%       PC                eph_manager.h
%       VC                eph_manager.h
%       TWOT              eph_manager.h
%       EM_RATIO          eph_manager.h
%       BUFFER            eph_manager.h
%       IPT               eph_manager.h
%       LPT               eph_manager.h
%       NRL               eph_manager.h
%       KM                eph_manager.h
%       NP                eph_manager.h
%       NV                eph_manager.h
%       RECORD_LENGTH     eph_manager.h
%       EPHFILE           eph_manager.h
%
%    FUNCTIONS
%    CALLED:
%       fclose            stdio.h
%       free              stdlib.h
%       fopen             stdio.h
%       fread             stdio.h
%       calloc            stdlib.h
%
%    VER./DATE/
%    PROGRAMMER:
%       V1.0/06-90/JAB (USNO/NA)
%       V1.1/06-92/JAB (USNO/AA): Restructure and add initializations.
%       V1.2/07-98/WTH (USNO/AA): Modified to open files for different
%                                 ephemeris types. (200,403,404,405,406)
%       V1.3/11-07/WKP (USNO/AA): Updated prolog.
%       V1.4/09-10/WKP (USNO/AA): Changed ncon and denum variables and
%                                 sizeof ipt array to type 'int' for
%                                 64-bit system compatibility.
%       V1.5/09-10/WTH (USNO/AA): Added support for DE421, default case
%                                 for switch, close file on error.
%       V1.6/10-10/WKP (USNO/AA): Renamed function to lowercase to
%                                 comply with coding standards.
%
%    NOTES:
%       KM...flag defining physical units of the output states.
%          = 1, km and km/sec
%          = 0, AU and AU/day
%       Default value is 0 (KM determines time unit for nutations.
%                           Angle unit is always radians.)
%
% ------------------------------------------------------------------------
  persistent m;
  if isempty(m)
    m = containers.Map('KeyType','char','ValueType','any');
  end
  if nargin < 4
    dontSave = false;
  end
  s.de_number = 0;
  s.jd_begin = 0.0;
  s.jd_end = 0.0;
  if isnumeric(denum)
    denum = sprintf('%d',denum);
  end
  error = 0;
  if isempty(ineph)
    error = 1;
    return;
  elseif exist(ineph, 'file') ~= 2
    inasc = ineph;
    inasc(1:3) = 'asc';
    Ephem.asc2eph(denum,inasc,ineph);
  end
  if exist(ineph, 'file') == 2
    % found locally
    if m.isKey(ineph)
      st = m(ineph);
    else
      EPHFILE = fopen (ineph,'rb');
      if EPHFILE ~= -1
        for i=1:3
          [st.ttl{i},count] = fread (EPHFILE, Ephem.TTLsize, '*char');
          if count ~= Ephem.TTLsize
            error = 1;
            fprintf(2,'openDEeph: Error TTL retrieved %d from "%s", expected %d\n',...
              count,ineph,Ephem.TTLsize);
          end
          st.ttl{i} = strcat(transpose(st.ttl{i}));
        end
        % pos 252
        for i=1:Ephem.OLDCONMAX
          [st.CNAM{i},count] = fread (EPHFILE, Ephem.CNAMsize, '*char');
          if count ~= Ephem.CNAMsize
            error = 1;
            fprintf(2,'openDEeph: Error CNAM retrieved %d from "%s", expected %d\n',...
              count,ineph,Ephem.CNAMsize);
          end
          st.CNAM{i} = strcat(transpose(st.CNAM{i}));
        end
        % pos 2652
        [st.SS,count] = fread (EPHFILE, 3, 'double');
        if count ~= 3
          error = 1;
          fprintf(2,'openDEeph: Error SS retrieved %d from "%s", expected %d\n',...
            count,ineph,3);
        end
        [st.NCON,count] = fread (EPHFILE, 1, 'int32');
        if count ~= 1
          error = 1;
          fprintf(2,'openDEeph: Error NCON retrieved %d from "%s", expected %d\n',...
            count,ineph,1);
        end
        [st.JPLAU,count] = fread (EPHFILE, 1, 'double');
        if count ~= 1
          error = 1;
          fprintf(2,'openDEeph: Error JPLAU retrieved %d from "%s", expected %d\n',...
            count,ineph,1);
        end
        [st.EM_RATIO,count] = fread (EPHFILE, 1, 'double');
        if count ~= 1
          error = 1;
          fprintf(2,'openDEeph: Error EM_RATIO retrieved %d from "%s", expected %d\n',...
            count,ineph,1);
        end
        % pos 2652+8*3+4+8+8
        for i=1:12
          for j=1:3
            [st.IPT(j,i),count] = fread (EPHFILE, 1, 'int32');
            if count ~= 1
              error = 1;
              fprintf(2,'openDEeph: Error IPT retrieved %d from "%s", expected %d\n',...
                count,ineph,1);
            end
          end
        end
        % pos 2652+8*3+4+8+8+12*3*4
        [st.de_number,count] = fread (EPHFILE, 1, 'int32');
        st.de_number = sprintf('%d',st.de_number);
        if count ~= 1
          error = 1;
          fprintf(2,'openDEeph: Error de_number retrieved %d from "%s", expected %d\n',...
            count,ineph,1);
        end
        % pos 2652+8*3+4+8+8+12*3*4+4
        if strcmp(st.de_number,denum(1:min(length(st.de_number),length(denum)))) == false
          if s.output > 0
            fprintf(s.output,'openDEeph: Warning DE# retrieved %s from "%s", expected %s\n',...
              st.de_number,ineph,denum);
          end
        end
        [LPT,count] = fread (EPHFILE, 3, 'int32');
        if count ~= 3
          error = 1;
          fprintf(2,'openDEeph: Error LPT retrieved %d from "%s", expected %d\n',...
            count,ineph,3);
        end
        % pos 2652+8*3+4+8+8+12*3*4+4+3*4
        for i=Ephem.OLDCONMAX+1:st.NCON
          [st.CNAM{i},count] = fread (EPHFILE, Ephem.CNAMsize, '*char');
          if count ~= Ephem.CNAMsize
            error = 1;
            fprintf(2,'openDEeph: Error CNAM retrieved %d from "%s", expected %d\n',...
              count,ineph,Ephem.CNAMsize);
          end
          st.CNAM{i} = strcat(transpose(st.CNAM{i}));
        end
        [RPT,count] = fread (EPHFILE, 3, 'int32');
        if count ~= 3
          error = 1;
          fprintf(2,'openDEeph: Error RPT retrieved %d from "%s", expected %d\n',...
            count,ineph,3);
        end
        [TPT,count] = fread (EPHFILE, 3, 'int32');
        if count ~= 3
          error = 1;
          fprintf(2,'openDEeph: Error TPT retrieved %d from "%s", expected %d\n',...
            count,ineph,3);
        end
        for i=1:3
          st.IPT(i,13) = LPT(i);
          st.IPT(i,14) = RPT(i);
          st.IPT(i,15) = TPT(i);
        end
        st.NCOEFF = Ephem.numCoeff( st.IPT );
        st.RECORD_LENGTH = 8*st.NCOEFF;
        status = fseek(EPHFILE,(2-1)*st.RECORD_LENGTH,'bof');
        if status == -1
          fprintf(2,'openDEeph seek error at %d\n',2);
        end
        [st.CVAL,count] = fread (EPHFILE, st.NCON, 'double');
        if count ~= st.NCON
          error = 1;
          fprintf(2,'openDEeph: Error CVAL retrieved %d from "%s", expected %d\n',...
            count,ineph,st.NCON);
        end
        for i=1:st.NCON
          cnm = deblank(st.CNAM{i});
          indx = strfind(cnm,'.');
          for j=1:length(indx)
            cnm(indx(j)) = '_';
          end
          if ~isempty(cnm)
            st.(cnm) = st.CVAL(i);
          else
            st.('blank') = st.CVAL(i);
          end
        end
        st.masses(Ephem.Mercury) = st.GM1;
        st.masses(Ephem.Venus) = st.GM2;
        st.masses(Ephem.Mars) = st.GM4;
        st.masses(Ephem.Jupiter) = st.GM5;
        st.masses(Ephem.Saturn) = st.GM6;
        st.masses(Ephem.Uranus) = st.GM7;
        st.masses(Ephem.Neptune) = st.GM8;
        st.masses(Ephem.Pluto) = st.GM9;
        st.masses(Ephem.Sun) = st.GMS;
        st.masses(Ephem.EarthMoonBarycenter) = st.GMB;
        st.masses(Ephem.Moon) = st.GMB/(1+s.EM_RATIO);
        st.masses(Ephem.Earth) = st.GMB*s.EM_RATIO/(1+s.EM_RATIO);
        mm = 0;
        for i=Ephem.Mercury:Ephem.Sun
          mm = mm + st.masses(i);
        end
        st.masses(Ephem.SolarSystemBarycenter) = mm;

        pos = 1;
        status = fseek(EPHFILE,(3-1)*st.RECORD_LENGTH,'bof');
        if status == -1
          if s.output > 0
            fprintf(s.output,'openDEeph: Bad fseek of "%s"\n',ineph);
          end
          fclose (EPHFILE);
          error = 9;
          return;
        end
        st.scan = zeros(1,0);
        for NROUT = 3:281474976710655
          [sc,count] = fread (EPHFILE,st.NCOEFF,'double');
          if isempty(sc) || count == 0 % expected end
            break;
          end
          if st.NCOEFF ~= count
            fprintf(2,'openDEeph: Error COEFF retrieved %d from "%s", expected %d\n',...
              count,ineph,st.NCOEFF);
            %break;
          else
            for K=1:st.NCOEFF
              st.scan(pos) = sc(K);
              pos = pos + 1;
            end
          end
        end
        fclose (EPHFILE);
        if isempty(st.scan)
          if s.output > 0
            fprintf(s.output,'openDEeph: Bad textscan of "%s"\n',ineph);
          end
          error = 9;
          return;
        end
        if dontSave == false
          m(ineph) = st;
        end
      else
        if s.output > 0
          fprintf(s.output,'openDEeph: Could not open "%s"\n',ineph);
        end
        error = 10;
        return;
      end
    end
  else
    if s.output > 0
      fprintf(2,'openDEeph: Bad file access for "%s"\n',ineph);
    end
    error = 10;
    return;
  end

  s.scan = st.scan;
  s.SS = st.SS;
  s.de_number = denum;
  s.RECORD_LENGTH = st.RECORD_LENGTH;
  s.IPT = st.IPT;
  s.JPLAU = st.JPLAU;
  s.EM_RATIO = st.EM_RATIO;
  for i=Ephem.Mercury:Ephem.EarthMoonBarycenter
    s.masses(i) = st.masses(i);
  end

  s.jd_begin = s.SS(1);
  s.jd_end = s.SS(2);
  s.jd_inc = s.SS(3);
  s.numRecs =  1 + (s.jd_end-s.jd_begin)/s.jd_inc;
  if mod(s.numRecs,1.0) ~= 0
    if s.output > 0
      fprintf(s.output,'openDEeph: DE%s Bad jd_inc (%g) numRecs (%g) should be integer\n',...
        denum,s.indexInc,s.numRecs);
    end
    s.jd_begin = 0.0;
    s.jd_end = 0.0;
    s.jd_inc = 0.0;
    error = 12;
    return;
  end
  s.indexInc = s.RECORD_LENGTH / 8;
  s.indexOffset = -1;
end

author: goGPS-Project


text:
function [X_norm, mu, sigma] = featureNormalize(X)
%FEATURENORMALIZE Normalizes the features in X
%   FEATURENORMALIZE(X) returns a normalized version of X where
%   the mean value of each feature is 0 and the standard deviation
%   is 1. This is often a good preprocessing step to do when
%   working with learning algorithms.

% You need to set these values correctly
X_norm = X;
mu = zeros(1, size(X, 2));
sigma = zeros(1, size(X, 2));

% ====================== YOUR CODE HERE ======================
% Instructions: First, for each feature dimension, compute the mean
%               of the feature and subtract it from the dataset,
%               storing the mean value in mu. Next, compute the
%               standard deviation of each feature and divide
%               each feature by it's standard deviation, storing
%               the standard deviation in sigma.
%
%               Note that X is a matrix where each column is a
%               feature and each row is an example. You need
%               to perform the normalization separately for
%               each feature.
%
% Hint: You might find the 'mean' and 'std' functions useful.
%

mu=mean(X)
sigma=std(X)

for i = 1:size(X,2),
    X_norm(:,i) = (X(:,i) - mu(i)) / sigma(i);
end

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

end

author: xjwhhh

repo: AndrewNgMachineLearning


text:
classdef CoreCompiledFileInfo < handle
    % CoreCompiledFileInfo. Structure for caching information about compiled files.
    %
    %     CoreCompiledFileInfo is used by CoreMex and its related files to manage
    %     dependencies on mex files. You typically create a map of
    %     CoreMexInfo objects, each of which defines a file to be
    %     compiled, and pass it into CoreMex. You can also use CoreCudaInfo
    %     for cuda files.
    %
    %     The basic rule you must follow is to increment the version number
    %     of your CoreMexInfo object whenever you change the corresponding
    %     mex file. This will force recompilation of any older versions.
    %
    %
    %     Licence
    %     -------
    %     Part of CoreMat. https://github.com/tomdoel/coremat
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        Name
        Extension
        Path
        CompilerOptions
        OtherCompilerFiles
        CurrentVersion

        LastSuccessfulCompiledVersion
        LastSuccessfulCompileDatenum
        LastAttemptedCompiledVersion
        LastAttemptedCompileDatenum
        LastSuccessfulCompiler
        LastAttemptedCompiler
        LastCompileFailed

        NeedsRecompile
        StatusID
    end

    methods
        function obj = CoreCompiledFileInfo(current_version, filename, extension, file_path, compiler_options, other_files)
            if nargin > 0
                obj.Name = filename;
                obj.Extension = extension;
                obj.Path = file_path;
                obj.CompilerOptions = compiler_options;
                obj.OtherCompilerFiles = other_files;
                obj.CurrentVersion = current_version;
            end
        end
    end

end

author: tomdoel

repo: pulmonarytoolkit


text:
function cnn_compute_image_sizes( expidx, image_set )

p = exp_params(expidx);

% load annolist
if strcmp(image_set, 'test')
    load(p.testGT)
else
    load(p.trainGT)
end

num_images = length(annolist);

sizes = zeros(num_images, 2);
for i = 1:num_images
    tic_toc_print('imdb: %d/%d\n', i, num_images);
    info = imfinfo(annolist(i).image.name);
    sizes(i, :) = [info.Height info.Width];
end

max(sizes, [], 1)

data_dir = fullfile(p.expDir, p.shortName, 'data');
mkdir_if_missing(data_dir);
save(fullfile(data_dir, ['sizes_' image_set]), 'sizes');

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function discon_meas = localized_discontinuity_measure(topographies, channel_locations, n_ics)
% localized_discontinuity_measure calculates a measure of how discontinuous
% an IC scalp map is. A high value indicates that the IC might represent an
% artifact. A very similar measure was used for the ADJUST toolbox,
% http://www.unicog.org/pm/pmwiki.php/MEG/RemovingArtifactsWithADJUST
% Our weighting scheme for electrodes surrounding the electrode for which
% the measure is being calculated is a little different.
%
% Input:
% topographies: Matrix of topographies (scalp maps) of ICs. Each row
% corresponds to an IC, and each column to an electrode, so it is the
% transpose of the field icawinv in EEGLab data structures.
%
% channel_locations: Channel locations of the electrodes whose activations
% are given in topographies.
%
% n_ics: The number of ICs represented in topographies.
%
% Output:
% discon_meas: A vector of discontinuity measures for each IC. For each IC,
% the discontinuity measure is found by first calculating how peaked the
% spatial map is for each electrode. The discontinuity measure is then the
% highest value found.

% This code is based heavily on the code from the ADJUST toolbox. It was
% modified by Laura Froelich 24/4/2013.
%
% Copyright (C) 2009 Andrea Mognon and Marco Buiatti,
% Center for Mind/Brain Sciences, University of Trento, Italy
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

xpos=[channel_locations.X];ypos=[channel_locations.Y];zpos=[channel_locations.Z];
pos=[xpos',ypos',zpos'];
discon_meas=zeros(1,n_ics);

for ic=1:n_ics
    aux=zeros(1,length(channel_locations));
    for el=1:length(channel_locations)-1

        P=pos(el,:); %position of current electrode
        d=pos-repmat(P,length(channel_locations),1);
        dist=sqrt(sum((d.*d),2));

        [y,I]=sort(dist, 'ascend'); % sort electrodes in ascending order of distance, so closest are first in list
        repchas=I(2:end);
        weightdist=exp(-y(2:end)); % weights computed wrt distance
        weightdist=weightdist/sum(weightdist); % ensure that weights sum to one

        aux(el)=abs(topographies(ic,el)-sum(weightdist.*topographies(ic,repchas)')); % Since the sum
        % of weights is one, the quantity
        % sum(weightdist.*topographies(ic,repchas)'), is a weighted
        % average. Hence this calculation finds the difference between the
        % activation in the current electrode and the average activation
        % surrounding it.
    end

    discon_meas(ic)=max(aux);
end

repo: REST


text:
function normA = norm(A, n)
%NORM   Norm of a CHEBMATRIX object.
%   NORM(A) computes the Frobenius norm of the CHEBMATRIX object A, defined as
%   the sum of the squares of the 2-norms of each of the blocks.
%
%   NORM(A, 2) or NORM(A, 'fro') is the same as above.
%
%   NORM(A, INF) return the infinity norm of the CHEBMATRIX A, defined as the
%   maximum infinity norm of each of the blocks.
%
%   NORM(A, -INF) returns the negative infinity norm of the CHEBMATRIX A,
%   defined as the minimum of the negative infinity norm of each of the blocks.
%
% See also CHEBMATRIX, CHEBFUN/NORM.

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

% [TODO]: Add support for norms of operators.

% Empty CHEBMATRIX has norm 0.
if ( isempty(A) )
    normA = 0;
    return
end

if ( nargin == 1 )
    n = 'fro'; 	% Frobenius norm is the default.
end

% The norm of a CHEBMATRIX with inf x inf block(s) that are neither
% CHEBFUN2 nor CHEBFUN3 is not supported:
s = cellfun(@(b) min(size(b)), A.blocks);
t = cellfun(@(b) isa(b, 'chebfun2') || isa(b, 'chebfun3'), A.blocks);
if ( ~all(isfinite(s(:))) )
    if ( ~any(t) )
        error('CHEBFUN:CHEBMATRIX:norm:notSupported', ...
            'Norm of a chebmatrix with operator block(s) is not supported.')
    end
end

% Deal with different cases.
switch n

    case {'fro', 2}
        normA = cellfun(@(v) norm(v, 2), A.blocks);
        normA = sqrt(sum(normA(:).^2));

    case {inf, 'inf'}
        normA = cellfun(@(v) norm(v, inf), A.blocks);
        normA = max(normA(:));

    case {-inf, '-inf'}
        normA = cellfun(@(v) norm(v, -inf), A.blocks);
        normA = min(normA(:));

    otherwise
        if ( ~ischar(n) )
            n = num2str(n);
        end
        error('CHEBFUN:CHEBMATRIX:norm:unknown', ...
            'unsupported norm type ''%s''', n)

end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [F,TimeVector] = in_fread_mne(sFile, ChannelMat, iEpoch, SamplesBounds, iChannels)
% IN_READ_MNE:  Read a block of recordings from a MNE-Python object
%
% USAGE:  [F,TimeVector] = in_fread_mne(sFile, ChannelMat, iEpoch, SamplesBounds, iChannels)
%         [F,TimeVector] = in_fread_mne(sFile, ChannelMat, iEpoch, SamplesBounds)            : Read all the channels

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2019

% Get reference to python object
pyObj = sFile.filename;
% List of channels
if (nargin < 5) || isempty(iChannels)
    picks = 'all';
else
    picks = {ChannelMat.Channel(iChannels).Name};
end
% Times
FileBounds = [bst_py2mat(pyObj.first_samp), bst_py2mat(pyObj.last_samp)];
if (nargin < 4) || isempty(SamplesBounds)
    SamplesBounds = FileBounds;
end

% Raw data
% if ~py.isinstance(pyObj, py.sys.modules{'mne.io'}.BaseRaw) || strcmpi(class(pyObj), 'py.mne.io.fiff.raw.Raw')
if ismethod(pyObj, 'get_data')
    % Call reading function
    res = pyObj.get_data(pyargs('picks', picks, 'start', int32(SamplesBounds(1) - FileBounds(1)), 'stop', int32(SamplesBounds(2) - FileBounds(1) + 1), 'return_times', true));
    F = bst_py2mat(res{1});
    TimeVector = bst_py2mat(res{2}) + sFile.prop.times(1);

% Epoched data
else
    error('todo');
%     % Use data already read
%     if isfield(sFile.header, 'epochData') && ~isempty(sFile.header.epochData)
%         F = permute(sFile.header.epochData(iEpoch,:,:), [2,3,1]);
%         TimeVector = linspace(sFile.epochs(iEpoch).times(1), sFile.epochs(iEpoch).times(2), size(F,2));
%     % Read data from file
%     else
%         [F, TimeVector] = fif_read_evoked(sFile, sfid, iEpoch);
%     end
%     % Specific selection of channels
%     if ~isempty(iChannels)
%         F = F(iChannels, :);
%     end
%     % Specific time selection
%     if ~isempty(SamplesBounds)
%         iTime = SamplesBounds - round(sFile.epochs(iEpoch).times(1) .* sFile.prop.sfreq) + 1;
%         F = F(:, iTime(1):iTime(2));
%         TimeVector = TimeVector(iTime(1):iTime(2));
%     end
end


text:
function geometry_test0846 ()

%*****************************************************************************80
%
%% TEST0846 tests POLYLOOP_POINTS_ND.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2009
%
%  Author:
%
%    John Burkardt
%
  dim_num = 2;
  nk = 4;
  nt = 12;

  pk = [ ...
    0.0, 2.0; ...
    0.0, 0.0; ...
    1.0, 0.0; ...
    1.0, 2.0]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0846\n' );
  fprintf ( 1, '  POLYLOOP_POINTS_ND computes points on a polyloop.\n' );

  r8mat_transpose_print ( dim_num, nk, pk, '  The defining points:' );

  pt = polyloop_points_nd ( dim_num, nk, pk, nt );

  r8mat_transpose_print ( dim_num, nt, pt, '  The computed points:' );

  return
end


text:
function map = sc_update_uvMap(map, data, uvPixInd)

% SC_UPDATE_UVMAP: Update the map from the data and the indices

[imgH, imgW, nCh] = size(map);

offset = uint32((0:nCh-1)*imgH*imgW);

uvPixInd = bsxfun(@plus, uvPixInd, offset);
map(uvPixInd) = data;

% for ch = 1: nCh
%     map(uvPixInd + (ch-1)*imgH*imgW) = data(:,ch);
% end
end

repo: StructCompletion


text:
function [out] = saturation_8(p1,p2,S,Smax,In)
%saturation_8

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Saturation excess flow from a store with different degrees
%               of saturation (min-max linear variant)
% Constraints:  -
% @(Inputs):    p1   - minimum fraction contributing area [-]
%               p2   - maximum fraction contributing area [-]
%               S    - current storage [mm]
%               Smax - maximum contributing storage [mm]
%               In   - incoming flux [mm/d]

out = (p1+(p2-p1)*S/Smax)*In;

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function dObj = pathObjective(u)
% dObj = pathObjective(u)
%
% Computes the integrand of the objective function, in this case,
% minimizing the integral of acceleration squared
%
% INPUTS:
%   u = [u1;u2];
%
% OUTPUTS:
%   dOBj = ddx.^2
%

% u1 = u(1,:);  % Unused
u2 = u(2,:);

dObj = u2.^2;

end

author: MatthewPeterKelly

repo: OptimTraj


text:
function newfp = slchangefilepart(fp, varargin)
%SLCHANGEFILEPART Changes some parts of the file path
%
% $ Syntax $
%   - newfp = slchangefilepart(fp, partname1, part1, ...)
%
% $ Description $
%   - newfp = slchangefilepart(fp, partname1, part1, ...) changes the
%     specified part of a path to a new value to form a new path.
%     Please refer to slfilepart for part names
%
% $ Remarks $
%   - If you specify the name, then you should be specify title and ext.
%
% $ History $
%   - Created by Dahua Lin, on Aug 12nd, 2006
%

%% Main body

if isempty(varargin)
    newfp = fp;
else

    opts.parent = '';
    opts.name = '';
    opts.title = '';
    opts.ext = '';
    opts = slparseprops(opts, varargin{:});

    if ~isempty(opts.ext)
        if opts.ext(1) ~= '.'
            error('sltoolbox:invalidarg', ...
                'The extension string should start with a dot . ');
        end
    end

    [p.parent, p.title, p.ext] = fileparts(fp);

    if isempty(opts.name)
        p = updatefields(p, opts, {'parent', 'title', 'ext'});
        newfp = fullfile(p.parent, [p.title, p.ext]);
    else
        if ~isempty(opts.title) || ~isempty(opts.ext)
            error('sltoolbox:invalidarg', ...
                'When name is specified, title and ext should not be');
        end
        p = updatefields(p, opts, {'parent', 'name'});
        newfp = fullfile(p.parent, p.name);
    end

end

%% Auxiliary functions

function S = updatefields(S, newS, fns)

nf = length(fns);
for i = 1 : nf
    f = fns{i};
    if ~isempty(newS.(f))
        S.(f) = newS.(f);
    end
end

author: lmthang


text:
function v = cvx_value( x )
global cvx___
try
    tv = subsref( cvx___.problems( x.problem_ ).duals, x.name_ );
catch
    v = NaN;
    return
end
v = value_( tv, cvx___.y );
function y = value_( x, y )
switch class( x ),
    case 'struct', y = cell2struct( value_( struct2cell( x ), y ), fieldnames( x ), 1 );
    case 'cell',   y = cellfun( @(z)value_( z, y ), x, 'UniformOutput', false );
    case 'cvx',    y = value( x, y );
    otherwise,     y = x;
end

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

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function parInfo = sweepExample()
% Copyright 2010 The MathWorks, Inc.

%% Initialize Problem
m     =         5;  % mass
bVals =  .1:.1 :5;  % damping values
kVals = 1.5:.05:5;  % stiffness values
[kGrid, bGrid] = meshgrid(bVals, kVals);
peakVals = nan(size(kGrid));

%% Parameter Sweep
% disp('Computing ...');drawnow;

parInfo = Par(numel(kGrid));

parfor idx = 1:numel(kGrid)
   Par.tic;
  % Solve ODE
  [T,Y] = ode45(@(t,y) odesystem(t, y, m, bGrid(idx), kGrid(idx)), ...
    [0, 25], ...  % simulate for 25 seconds
    [0, 1]);      % initial conditions

  % Determine peak value
  peakVals(idx) = max(Y(:,1));
  parInfo(idx) = Par.toc;
end

stop(parInfo);

function dy = odesystem(t, y, m, b, k)
% 2nd-order ODE
%
%   m*X'' + b*X' + k*X = 0
%
% --> system of 1st-order ODEs
%
%   y  = X'
%   y' = -1/m * (k*y + b*y')

dy(1) = y(2);
dy(2) = -1/m * (k * y(1) + b * y(2));

dy = dy(:); % convert to column vector

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function obj = Mrate_method(obj)
% traditional SVD algorithm for rate maximization

global  H Ns V_ropt W_ropt;
t1 = clock;
[U,~,V] = svd(H);
V_ropt = V(:,1:Ns);
%power constraint
V_ropt = V_ropt / norm(V_ropt,'fro');
W_ropt = U(:,1:Ns);

t2 = clock;
runtime = etime(t2,t1);
obj.runtime = obj.runtime + runtime;
obj.V_B = V_ropt;
obj.W_B = W_ropt;
obj = get_metric(obj);

author: Zzhaoxingyu

repo: hybrid-beamforming-for-three-scenes


text:
function a = in( x, y )

narginchk(2,2);
b = newcnstr( evalin( 'caller', 'cvx_problem', '[]' ), x, y, '==' );
if nargout, a = b; end

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

author: cvxr


text:
function xn=linEqSolveFirstnRows(A,b,n)
%%LINEQSOLVEFIRSTNROWS This function solves the problem A*x=b for the first
%                      n rows of x. The matrix A CAN be singular as long as
%                      a unique solution exists for the first n rows of x.
%                      This function is useful for finding the position
%                      estimate in an information filter state even before
%                      estimates of all of the other target state
%                      components have become observable In such an
%                      instance A is the inverse covariance matrix, x is
%                      the target state and b is the information state.
%
%INPUTS: A The NXN matrix A in the equation A*x=b, where x is unknown. The
%          matrix can be singular, but the first n rows of x should be
%          observable.
%        b The NX1 column vector b in the equation A*x=b.
%        n The number of rows of x, starting from the first row and going
%          down, for which one wishes to solve in the equation A*b. n must
%          be less than or equal to the number of rows in A.
%
%OUTPUTS: xn The first n rows of the column vector of x solved from A*x=b.
%            If any of the components of xn are not finite, then the
%            problem was not completely observable. Because of how the
%            problem is solved, if the kth component is not finite, then
%            all subsequent components will not be finite, regardless of
%            whether they are observable.
%
%The linear equation is solved by performing a modified qr decomposition on
%the matrix A such that A=Q*R, where Q is an rotation matrix and R is a
%LOWER triangular matrix. A standard QR decomposition produces an upper
%triangular matrix. By flipping the rows and columns of a and then
%performing the inverse operations on Q and R, one can get a decomposition
%where R is a lower-triangular marix. One can then write R*x=Q'*b, since
%the transpose of a rotation matrix is its inverse. The first n rows of x
%can then be solved using forward substitution.
%
%The QR decomposition is in Matlab's built-in function qr. It is also
%discussed in Chapter 5.2 of [1].
%
%REFERENCES:
%[1] G. E. Golub and C. F. van Loan, Matrix Computations, 4rd ed.
%    Baltimore, MD: Johns Hopkins University Press, 2013.
%
%June 2014 David F.Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%Perform a lower-triangular qr decomposition.
[Q,R]=qr(rot90(A,2));
Q=rot90(Q,2);
R=rot90(R,2);
%Now, Q*R=A and R is lower-triangular.

b=Q'*b;

%Perform forward substitution to solve for the first n components of x.
xn=zeros(n,1);
xn(1)=b(1)/R(1,1);
for curRow=2:n
    xn(curRow)=(b(curRow)-sum(xn(1:(curRow-1))'.*R(curRow,1:(curRow-1))))/R(curRow,curRow);
end

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function base = stat_getBaselineDistrib(PConn,baseline,winCenterTimes,noavg)
% Return the distribution of an estimator averaged over a baseline and
% expanded out to the same dimensions as matrices in PConn connectivity
% object
%
% Inputs:
%
%       PConn:     connectivity distribution returned from stat_bootstrap()
%       baseline:  [min max] baseline interval relative to event (t=0 sec)
%
% Optional
%
%       winCenterTimes: vector of times corresponding to window centers
%                      (same units as baseline)
%       noavg:     if true, the baseline distribution is returned
% Outputs:
%
%       base:      the baseline mean distribution (expanded out to same dims as
%                  PConn)
%
% See Also: stat_bootstrap()
%
% References:
%
% [1] Mullen T (2010) The Source Information Flow Toolbox (SIFT):
%   Theoretical Handbook and User Manual. Chapter 6.8
%   Available at: http://www.sccn.ucsd.edu/wiki/SIFT
%
% Author: Tim Mullen, 2011, SCCN/INC, UCSD.
% Email:  tim@sccn.ucsd.edu

% This function is part of the Source Information Flow Toolbox (SIFT)
%
% 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, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

if nargin<4
    noavg = false;
end

if isstruct(PConn)
    % connectivity object
    connmethods = hlp_getConnMethodNames(PConn);

    for cnd=1:length(PConn)

        baseidx = getindex(PConn(cnd).erWinCenterTimes,baseline);

        for m=1:length(connmethods)
            sz = size(PConn(cnd).(connmethods{m}));
            if noavg
                base.(connmethods{m}) = PConn(cnd).(connmethods{m})(:,:,:,baseidx(1):baseidx(2),:);
            else
                basetmp = (mean(PConn(cnd).(connmethods{m})(:,:,:,baseidx(1):baseidx(2),:),4));  % collapse baseline over timepoints
                base.(connmethods{m}) = repmat(basetmp,[1 1 1 sz(4) 1]);  % expands out to appropriate size
            end
        end

    end

else
    if nargin<3
        error('SIFT:stat_getBaselineDistrib','You must provide winCenterTimes');
    end

    % single matrix
    baseidx = getindex(winCenterTimes,baseline);
    sz = size(PConn);
    if noavg
        base = PConn(:,:,:,baseidx(1):baseidx(2),:);
    else
        basetmp = (mean(PConn(:,:,:,baseidx(1):baseidx(2),:),4));   % collapse baseline over timepoints
        base = repmat(basetmp,[1 1 1 sz(4) 1]);  % expands out to appropriate size
    end
end


text:
% test for line integral convolution

n = 180;
path(path, 'toolbox/');

rep = 'results/lic/';
if not(exist(rep))
    mkdir(rep);
end

options.bound = 'per';
%% generate a random irrotational vector field
sigma = 30;
v = perform_blurring(randn(n,n,2), sigma, options);
v = perform_vf_normalization(v);

options.histogram = 'gaussian';
options.histogram = 'linear';
options.verb = 1;
% size of the features
options.spot_size = 1.3;

%% original image
name = 'mandrill-color';
name = 'flowers';
name = 'rand';
options.dt = 0.5;
w_list = [4 6 8 10 12 14];
if strcmp(name, 'rand')
    sigma = 1.2;
    M0 = perform_blurring(randn(n), sigma, options);
    M0 = perform_histogram_equalization(M0, options.histogram);
else
    options.dt = 1.5;
    w_list = w_list * 2;
    options.histogram = [];
    M0 = load_image(name,n);
    M0 = rescale( crop(M0,n) );
end

%% iterated lic
M = M0;
w = 12;
for i=1:4
    options.M0 = M;
    M = perform_lic(v, w, options);
end

%% lic for increasing times
close all;
clf;
options.M0 = M0;
for i=1:min(length(w_list),6)
    M = perform_lic(v, w_list(i), options);
    imageplot(M, '', 2,3, i);
end

saveas(gcf, [rep name '-lic-results.png'], 'png');

return;

%% display with overlapping vector field
if size(M,3)==1
    figure; clf;
    sub = 4;
    plot_vf(v(1:sub:end, 1:sub:end,:), M);
    colormap gray(256);
    saveas(gcf, [rep, 'lic-vector-field.png'], 'png');
end

author: gpeyre

repo: matlab-toolboxes


text:
function [w, infos] = sd(problem, in_options)
% Full steepest descent gradient algorithm.
%
% Inputs:
%       problem     function (cost/grad/hess)
%       in_options  options
% Output:
%       w           solution of w
%       infos       information
%
% This file is part of GDLibrary and SGDLibrary.
%
% Created by H.Kasai on Feb. 15, 2016
% Modified by H.Kasai on Mar. 23, 2018
% Modified by H.Kasai on Oct. 20, 2020

    % set dimensions and samples
    d = problem.dim;
    n = problem.samples;

    % set local options
    local_options = [];
    local_options.algorithm = 'SD';
    local_options.sub_mode = 'STANDARD';

    % merge options
    options = mergeOptions(get_default_options(d), local_options);
    options = mergeOptions(options, in_options);

    % initialise
    iter = 0;
    grad_calc_count = 0;
    w = options.w_init;
    w_old = w;
    prev_step = options.step_init;

    if ~isfield(options, 'S')
        if strcmp(options.step_alg, 'exact')
            options.S = eye(d);
        end
    else
        %
    end

    % initialize by BB step-size
    if strcmp(options.step_init_alg, 'bb_init')
        options.step_init = bb_init(problem, w);
    end

    % store first infos
    clear infos;
    [infos, f_val, optgap, grad, gnorm] = store_infos(problem, w, options, [], iter, grad_calc_count, 0);
    grad_old = grad;

    % display info
    if options.verbose
        if ~problem.prox_flag
            fprintf('SD: Iter = %03d, cost = %.24e, gnorm = %.4e, optgap = %.4e\n', iter, f_val, gnorm, optgap);
        else
            fprintf('PG: Iter = %03d, cost = %.24e, gnorm = %.4e, optgap = %.4e\n', iter, f_val, gnorm, optgap);
        end
    end

    % set start time
    start_time = tic();

    % main loop
    while (optgap > options.tol_optgap) && (gnorm > options.tol_gnorm) && (iter < options.max_epoch)

        options.iter = iter;
        [step, ~] = options.linesearchfun(options.step_alg, problem, w, w_old, grad, grad_old, prev_step, options);

        prev_step = step;
        w_old = w;
        if strcmp(options.sub_mode, 'SCALING')
            % diagonal scaling
            if isempty(options.S)
                h = problem.full_hess(w);
                options.S = diag(1./diag(h));
            end

            % update w
            w = w - step * options.S * grad;
        else
            % update w
            w = w - step * grad;
        end

        % proximal operator
        if problem.prox_flag
            w = problem.prox(w, step);
        end

        % store gradient
        grad_old = grad;

        % measure elapsed time
        elapsed_time = toc(start_time);

        % count gradient evaluations
        grad_calc_count = grad_calc_count + n;

        % update iter
        iter = iter + 1;

        % store infos
        [infos, f_val, optgap, grad, gnorm] = store_infos(problem, w, options, infos, iter, grad_calc_count, elapsed_time);

        % display infos
        if options.verbose
            if ~problem.prox_flag
                fprintf('SD: Iter = %03d, cost = %.24e, gnorm = %.4e, optgap = %.4e\n', iter, f_val, gnorm, optgap);
            else
                fprintf('PG: Iter = %03d, cost = %.24e, gnorm = %.4e, optgap = %.4e\n', iter, f_val, gnorm, optgap);
            end
        end
    end

    if gnorm < options.tol_gnorm
        fprintf('Gradient norm tolerance reached: tol_gnorm = %g\n', options.tol_gnorm);
    elseif optgap < options.tol_optgap
        fprintf('Optimality gap tolerance reached: tol_optgap = %g\n', options.tol_optgap);
    elseif iter == options.max_epoch
        fprintf('Max iter reached: max_epoch = %g\n', options.max_epoch);
    end

end

author: hiroyuki-kasai

repo: SGDLibrary


function currentScaleFactor = estimate_scale( im_gray, pos, currentScaleFactor)

global para

% extract the test sample feature map for the scale filter
xs = get_scale_sample(im_gray, pos, para.app_sz, para.scaleFactors*currentScaleFactor, para.scale_window, para.scale_model_sz);

% calculate the correlation response of the scale filter
xsf = fft(xs,[],2);
scale_response = real(ifft(sum(para.sf_num .* xsf, 1) ./ (para.sf_den + para.lambda)));

% find the maximum scale response
recovered_scale = find(scale_response == max(scale_response(:)), 1);

currentScaleFactor = currentScaleFactor*para.scaleFactors(recovered_scale);
if currentScaleFactor < para.min_scale_factor
    currentScaleFactor = para.min_scale_factor;
elseif currentScaleFactor > para.max_scale_factor
    currentScaleFactor = para.max_scale_factor;
end

% update the scale model
%===========================
% extract the training sample feature map for the scale filter
xs = get_scale_sample(im_gray, pos, para.app_sz, currentScaleFactor * para.scaleFactors, para.scale_window, para.scale_model_sz);

% calculate the scale filter update
xsf = fft(xs,[],2);
new_sf_num = bsxfun(@times, para.ysf, conj(xsf));
new_sf_den = sum(xsf .* conj(xsf), 1);

para.sf_den = (1 - para.interp_factor) * para.sf_den + para.interp_factor * new_sf_den;
para.sf_num = (1 - para.interp_factor) * para.sf_num + para.interp_factor * new_sf_num;

end


text:
function str = showdetail(item)

% function str = showdetail(item)
% Display details for a cfg_choice and all of its options.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: showdetail.m 4867 2012-08-30 13:04:51Z volkmar $

rev = '$Rev: 4867 $'; %#ok

str = showdetail(item.cfg_item);
str{end+1} = 'Class  : cfg_mchoice';
str{end+1} = ['Set one or more of the options listed as a struct with a ' ...
              'fields:'];
% Display detailed help for each cfg_choice value item
tags = tagnames(item, true);
for k = 1:numel(tags)
    str{end+1} = sprintf('.%s', tags{k});
end;

author: spm


text:
function [theta, J_history] = gradientDescentMulti(X, y, theta, alpha, num_iters)
%GRADIENTDESCENTMULTI Performs gradient descent to learn theta
%   theta = GRADIENTDESCENTMULTI(x, y, theta, alpha, num_iters) updates theta by
%   taking num_iters gradient steps with learning rate alpha

% Initialize some useful values
m = length(y); % number of training examples
J_history = zeros(num_iters, 1);

for iter = 1:num_iters

    % ====================== YOUR CODE HERE ======================
    % Instructions: Perform a single gradient step on the parameter vector
    %               theta.
    %
    % Hint: While debugging, it can be useful to print out the values
    %       of the cost function (computeCostMulti) and gradient here.
    %

    S = (1 / m) * (X' * (X * theta - y));
    theta = theta - alpha .* S;

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

    % Save the cost J in every iteration
    J_history(iter) = computeCostMulti(X, y, theta);

end

end

author: Ayatans

repo: Machine-Learning-homework


text:
% This is like skf_data_assoc_gmux, except the objects don't move.
% We are uncertain of their initial positions, and get more and more observations
% over time. The goal is to test deterministic links (0 covariance).
% This is like robot1, except the robot doesn't move and is always at [0 0],
% so the relative location is simply L(s).

nobj = 2;
N = nobj+2;
Xs = 1:nobj;
S = nobj+1;
Y = nobj+2;

intra = zeros(N,N);
inter = zeros(N,N);
intra([Xs S], Y) =1;
for i=1:nobj
  inter(Xs(i), Xs(i))=1;
end

Xsz = 2; % state space = (x y)
Ysz = 2;
ns = zeros(1,N);
ns(Xs) = Xsz;
ns(Y) = Ysz;
ns(S) = nobj;

bnet = mk_dbn(intra, inter, ns, 'discrete', S, 'observed', [S Y]);

% For each object, we have
% X(t+1) = F X(t) + noise(Q)
% Y(t) = H X(t) + noise(R)
F = eye(2);
H = eye(2);
Q = 0*eye(Xsz); % no noise in dynamics
R = eye(Ysz);

init_state{1} = [10 10]';
init_state{2} = [10 -10]';
init_cov = eye(2);

% Uncertain of initial state (position)
for i=1:nobj
  bnet.CPD{Xs(i)} = gaussian_CPD(bnet, Xs(i), 'mean', init_state{i}, 'cov', init_cov);
end
bnet.CPD{S} = root_CPD(bnet, S); % always observed
bnet.CPD{Y} = gmux_CPD(bnet, Y, 'cov', repmat(R, [1 1 nobj]), 'weights', repmat(H, [1 1 nobj]));
% slice 2
eclass = bnet.equiv_class;
for i=1:nobj
  bnet.CPD{eclass(Xs(i), 2)} = gaussian_CPD(bnet, Xs(i)+N, 'mean', zeros(Xsz,1), 'cov', Q, 'weights', F);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create LDS params

% X(t) = A X(t-1) + B U(t) + noise(Q)

% [L11]  = [1  ]  * [L1]       +   [Q ]
% [L2]     [  1]    [L2]           [ Q]

% Y(t)|S(t)=s  = C(s) X(t) + noise(R)
% Yt|St=1 = [1 0]  * [L1]  + R
%                    [L2]

nlandmarks = nobj;

% Create indices into block structure
bs = 2*ones(1, nobj); % sizes of blocks in state space
for i=1:nlandmarks
  landmark_block(:,i) = block(i, bs)';
end
Xsz = 2*(nlandmarks); % 2 values for each landmark plus robot
Ysz = 2; % observe relative location

% create block-diagonal trans matrix for each switch
A = zeros(Xsz, Xsz);
for i=1:nlandmarks
  bi = landmark_block(:,i);
  A(bi, bi) = eye(2);
end
A = repmat(A, [1 1 nlandmarks]); % same for all switch values

% create block-diagonal system cov
Qbig = zeros(Xsz, Xsz);
Qbig = repmat(Qbig, [1 1 nlandmarks]);

% create observation matrix for each value of the switch node
% C(:,:,i) = (0 ... I ...) where the I is in the i'th posn.
C = zeros(Ysz, Xsz, nlandmarks);
for i=1:nlandmarks
  C(:, landmark_block(:,i), i) = eye(2);
end

% create observation cov for each value of the switch node
Rbig = repmat(R, [1 1 nlandmarks]);

% initial conditions
init_x = [init_state{1}; init_state{2}];
init_V = zeros(Xsz, Xsz);
for i=1:nlandmarks
  bi = landmark_block(:,i);
  init_V(bi,bi) = init_cov;
end

%%%%%%%%%%%%%%%%
% Observe objects at random
T = 10;
evidence = cell(N, T);
data_assoc = sample_discrete(normalise(ones(1,nobj)), 1, T);
evidence(S,:) = num2cell(data_assoc);
evidence = sample_dbn(bnet, 'evidence', evidence);

% Inference
ev = cell(N,T);
ev(bnet.observed,:) = evidence(bnet.observed, :);
y = cell2num(evidence(Y,:));

engine = pearl_unrolled_dbn_inf_engine(bnet);
engine = enter_evidence(engine, ev);

loopy_est_pos = zeros(2, nlandmarks);
loopy_est_pos_cov = zeros(2, 2, nlandmarks);
for i=1:nobj
  m = marginal_nodes(engine, Xs(i), T);
  loopy_est_pos(:,i) = m.mu;
  loopy_est_pos_cov(:,:,i) = m.Sigma;
end

[xsmooth, Vsmooth] = kalman_smoother(y, A, C, Qbig, Rbig, init_x, init_V, 'model', data_assoc);

kf_est_pos = zeros(2, nlandmarks);
kf_est_pos_cov = zeros(2, 2, nlandmarks);
for i=1:nlandmarks
  bi = landmark_block(:,i);
  kf_est_pos(:,i) = xsmooth(bi, T);
  kf_est_pos_cov(:,:,i) = Vsmooth(bi, bi, T);
end

kf_est_pos
loopy_est_pos

kf_est_pos_time = zeros(2, nlandmarks, T);
for t=1:T
  for i=1:nlandmarks
    bi = landmark_block(:,i);
    kf_est_pos_time(:,i,t) = xsmooth(bi, t);
  end
end
kf_est_pos_time % same for all t since smoothed

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function lambda = chow_eigenvalues ( alpha, beta, n )

%*****************************************************************************80
%
%% CHOW_EIGENVALUES returns the eigenvalues of the CHOW matrix.
%
%  Example:
%
%    ALPHA = 2, BETA = 3, N = 5
%
%    9.49395943
%    6.10991621
%    3.0
%    3.0
%    3.0
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ALPHA, the ALPHA value.  A typical value is 1.0.
%
%    Input, real BETA, the BETA value.  A typical value is 0.0.
%
%    Input, integer N, the order of A.
%
%    Output, real LAMBDA(N,1), the eigenvalues of A.
%
  lambda = zeros ( n, 1 );

  k = n - round ( ( n + 1 ) / 2 );

  for i = 1 : k
    angle = i * pi / ( n + 2 );
    lambda(i,1) = beta + 4.0 * alpha * ( cos ( angle ) )^2;
  end

  lambda(k+1:n,1) = beta;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [massImbalance, imBalancedMass, imBalancedCharge, imBalancedRxnBool, elements, missingFormulaeBool, balancedMetBool] = checkMassChargeBalance(model, printLevel, modelName)
% Tests for a list of reactions if these reactions are
% mass-balanced by adding all elements on left hand side and comparing them
% with the sums of elements on the right hand side of the reaction.
% A reaction is considered elementally imbalanced if any of the molecular
% species involved is missing a chemical formula.
%
% USAGE:
%
%    [massImbalance, imBalancedMass, imBalancedCharge, imBalancedRxnBool, elements, missingFormulaeBool, balancedMetBool] = checkMassChargeBalance(model, printLevel, modelName)
%
% INPUT:
%    model:                  COBRA model structure:
%
%                              * .S - `m` x `n` stoichiometric matrix
%                              * .metForumlas - `m` x 1 cell array of metabolite formulas
%                              * .metCharges - `m` x 1 double array of charges
%                              * model.SIntRxnBool - Boolean of reactions heuristically though to be mass balanced. (optional)
%
% OPTIONAL INPUTS:
%    printLevel:             {-1, (0), 1} where:
%
%                            -1 = print out diagnostics on problem reactions to a file,
%                            0 = silent,
%                            1 = print elements as they are checked (display progress),
%                            2 = also print out diagnostics on problem reactions to screen
%    modelName:               name of the file
%
% OUTPUTS:
%    massImbalance:          `nRxn` x `nElement` matrix with mass imbalance
%                            for each element checked. 0 if balanced.
%    imBalancedMass:         `nRxn` x 1 cell with mass imbalance
%                            e.g. -3 H means three hydrogens disappear in the reaction.
%    imBalancedCharge:       `nRxn` x 1 vector with charge imbalance, empty if no imbalanced reactions
%    imBalancedRxnBool:      boolean vector indicating imbalanced reactions (including exchange reactions!)
%    elements:               `nElement` x 1 cell array of element abbreviations checked
%    missingFormulaeBool:    `nMet` x 1 boolean vector indicating metabolites without formulae
%    balancedMetBool:        boolean vector indicating metabolites exclusively involved in balanced reactions
%
%
% OPTIONAL OUTPUT FILES:
% modelName_mass_imbalanced_reactions.txt  provides a human readable summary
% for each mass imbalanced reaction.
%
% For each reaction it contains:
% j             the index of the column of the stoichiometric matrix corresponding to the reaction
% rxns{j}       the abbreviation of the reaction
% imbalance     the number of each element that is unbalanced
% equation      the reaction equation
%
% For each metabolite involved in a reaction it contains
% i               the index of the row of the stoichiometric matrix corresponding to the metabolite
% mets{i}         the abbreviation of the reaction
% S(i,j)          the stoichiometric coefficient
% metFormulas{i}  the chemical formula of the metabolite
%
% Example
% #Rxn;rxnAbbr;imbalance;equation
% 32;2AMACSULT;3 O, 1 S;2amac[c] + nadph[c] + paps[c]  -> nadp[c] + Lcyst[c] + pap[c]
%    58	       2amac[c]	-1	C3H5NO2
%    60	       nadph[c]	-1	C21H26N7O17P3
%    61	        nadp[c]	1	C21H25N7O17P3
%    62	        paps[c]	-1	C10H11N5O10P2
%    63	       Lcyst[c]	1	C3H6NO5S
%    64	         pap[c]	1	C10H11N5O10P2
%
% .. Authors:
%       - Ines Thiele 12/09
%       - IT, 06/10, Corrected some bugs and improved speed.
%       - RF, 09/09/10, Support for very large models and printing to file.
%       - RF, 18/12/14, Default is now to check balancing of all reactions.

[nMet, nRxn]=size(model.S);

if ~exist('printLevel', 'var')
    printLevel=0;
end

if ~isfield(model,'SIntRxnBool')
    model.SIntRxnBool=true(nRxn,1);%assume all reactions are supposed to be internal if no other info provided
end

if ~exist('modelName','var')
    modelName='';
end
if ~isfield(model,'rxns')
    for i=1:nRxn
        model.rxns{i,1}=['rxn' int2str(i)];
    end
end
% List of elements
elements = {'H','C', 'O', 'P', 'S', 'N', 'Mg','X','Fe','Zn','Co','Ca','Y','I','Na','Cl','K','R','FULLR'};

E=sparse(nMet, length(elements));
massImbalance=zeros(nRxn, length(elements));
missingFormulaeBool=cellfun(@isempty, model.metFormulas);
for j = 1 : length(elements)
    if j==1
        [dE, E_el, missingFormulaeBool]=checkBalance(model, elements{j}, printLevel,[],missingFormulaeBool);
        massImbalance(:, j)=dE;
        E(:, j)=E_el;
        if printLevel>0
            fprintf('%s\n', ['Checked element ' elements{j}]);
        end
    else
        %no need to print out for each element which metabolites have no
        %formula
        [massImbalance(:, j), E(:, j)]=checkBalance(model, elements{j}, 0);
        if printLevel>0
            fprintf('%s\n', ['Checking element ' elements{j}]);
        end
    end
end

%ignore mass imbalance of exchange reactions if the internal reactions have
%been identified at the input
if any(~model.SIntRxnBool)
    massImbalance(~model.SIntRxnBool,:)=1;
end

E=full(E);

% A reaction is considered elementally imbalanced if any of the molecular
% species involved is missing a chemical formula.
imBalancedRxnBool=any(massImbalance, 2) | any(model.S(missingFormulaeBool, :),1)';

imBalancedMass=cell(nRxn, 1);
for i = 1 : nRxn
    imBalancedMass{i, 1}='';
    if imBalancedRxnBool(i)
        for j = 1 : length(elements)
            if massImbalance(i, j) ~= 0
                if ~strcmp(imBalancedMass{i, 1}, '')
                    imBalancedMass{i, 1} = [imBalancedMass{i, 1} ', ' int2str(massImbalance(i, j)) ' ' elements{j}];
                else
                    imBalancedMass{i, 1} = [int2str(massImbalance(i, j)) ' ' elements{j}];
                end
            end
        end
        if strfind(imBalancedMass{i, 1}, 'NaN')
            imBalancedMass{i, 1}='NaN';
        end
    end
    if mod(i, 1000)==0
        fprintf('%n\t%s\n', i, ['reactions checked for ' elements{j} ' balance']);
    end
end
if printLevel==-1
    firstMissing=0;
    for p=1:nRxn
        %only print out for reactions supposed to be mass balanced
        if model.SIntRxnBool(p) && ~strcmp(imBalancedMass{p,1},'')
            %at the moment, ignore reactions with a metabolite that have
            %no formula
            if ~strcmp(imBalancedMass{p, 1}, 'NaN')
                if ~firstMissing
                    fid=fopen([modelName 'mass_imbalanced_reactions.txt'],'w');
                    fprintf(fid, '%s;%s;%s;%s\n', 'j', 'rxns{j}', 'imbalance', 'equation');
                    fprintf(fid, '%s;%s;%s;%s\n\n', 'i', 'mets{i}', 'S(i,j)', 'metFormulas{i}');
                    fprintf('%s\n',['There are mass imbalanced reactions, see ' modelName 'mass_imbalanced_reactions.txt'])
                    firstMissing=1;
                end
                equation=printRxnFormula(model, model.rxns(p), 0);
                fprintf(fid, '%s;%s;%s;%s\n', int2str(p), model.rxns{p}, imBalancedMass{p, 1}, equation{1});
                for m=1:size(model.S, 1)
                    if model.S(m, p) ~= 0
                        fprintf(fid,'%5u\t%15s\t%g\t%s\n',m,model.mets{m},full(model.S(m,p)),model.metFormulas{m});
                    end
                end
            end
        end
    end
    if firstMissing
        fclose(fid);
    end
end
if printLevel==2
    for p=1:nRxn
        %only print out for reactions supposed to be mass balanced
        if model.SIntRxnBool(p) && ~strcmp(imBalancedMass{p,1},'')
            %at the moment, ignore reactions with a metabolite that have
            %no formula
            if ~strcmp(imBalancedMass{p, 1}, 'NaN')
                equation=printRxnFormula(model, model.rxns(p), 0);
                fprintf('%6s\t%30s\t%10s\t%s\n', int2str(p), model.rxns{p}, imBalancedMass{p, 1}, equation{1});
                % for m=1:size(model.S, 1)
                %     if model.S(m, p) ~= 0
                %         fprintf('%5u\t%15s\t%g\t%s\n',m,model.mets{m},full(model.S(m,p)),model.metFormulas{m});
                %     end
                % end

            end
        end
    end
end

%
if nnz(strcmp('', imBalancedMass))==nRxn
    imBalancedMass=[];
end

% Check for charge balance (initialize with NaN, if the fields are not set
% this will make it clear.
imBalancedCharge = NaN * ones(nRxn, 1);
firstMissing=0;
if isfield(model, 'metCharges')
    for m=1:nMet
        if isnan(model.metCharges(m)) && ~isempty(model.metFormulas{m})
            if printLevel==2
                fprintf('%s\t%s\n', int2str(m), [model.mets{m} ' has no charge but has formula.'])
                if ~firstMissing
                    warning('model structure must contain model.metCharges field for each metabolite');
                end
                firstMissing=1;
            end
            if printLevel==-1
                if ~firstMissing
                    fid=fopen([modelName 'metabolites_without_charge.txt'],'w');
                end
                firstMissing=1;
                fprintf(fid, '%s\t%s\n', int2str(m), model.mets{m});
            end
        end
    end

    imBalancedCharge = model.S' * double(model.metCharges); % Matlab does not support this operation on two int values - one needs to be converted to double. The smaller matrix is selected.
end

if printLevel==-1 && isfield(model,'SIntRxnBool')
    firstMissing=0;
    if any(imBalancedCharge)
        for q=1:nRxn
            if model.SIntRxnBool(q) && strcmp(imBalancedMass{q, 1}, '') && imBalancedCharge(q) ~= 0
                if ~firstMissing
                    fid=fopen([modelName 'mass_balanced_charge_imbalanced_reactions.txt'],'w');
                    fprintf('%s\n',['There are mass balanced, but charge imbalanced reactions, see ' modelName 'charge_imbalanced_reactions.txt'])
                    firstMissing=1;
                end
                equation=printRxnFormula(model, model.rxns(q), 0);
                fprintf(fid, '%s\t%s\t%s\n', int2str(q), model.rxns{q}, equation{1});
                % for m=1:size(model.S, 1)
                %     if model.S(m, q) ~= 0
                %         fprintf(fid,'%5u\t%15s\t%g\t%g\t%s\n',m,model.mets{m},full(model.S(m,q)),model.metCharges(m),model.metFormulas{m});
                %     end
                % end
            end
        end
        if firstMissing
            fclose(fid);
        end
    end
end

if printLevel==2 && isfield(model,'SIntRxnBool')
    if any(imBalancedCharge)
        fprintf('%s\n', 'Mass balanced, but charged imbalanced reactions:')
        for q=1:nRxn
            if model.SIntRxnBool(q) && strcmp(imBalancedMass{p, 1}, '') && imBalancedCharge(q) ~= 0
                equation=printRxnFormula(model, model.rxns(q), 0);
                fprintf('%s\t%s\t%s\n', int2str(q), model.rxns{q}, equation{1});
                for m=1:size(model.S, 1)
                    if model.S(m, q) ~= 0
                        fprintf('%5u\t%15s\t%g\t%g\t%s\n',m,model.mets{m},full(model.S(m,q)),model.metCharges(m),model.metFormulas{m});
                    end
                end
            end
        end
    end
end

%If the field is set we should assume, that all values are defined.
if isfield(model, 'metCharges')
    imBalancedRxnBool = imBalancedRxnBool | imBalancedCharge ~= 0;
end

if isfield(model,'SIntRxnBool')
    imBalancedRxnBool = imBalancedRxnBool | ~model.SIntRxnBool;
end

%nonzero rows corresponding to completely mass balanced reactions
balancedMetBool = getCorrespondingRows(model.S,true(nMet,1),~imBalancedRxnBool,'exclusive');
massImbalance = sparse(massImbalance);

author: opencobra

repo: cobratoolbox