function moretests(ut)
%%% summary: moretests
%%% project: emgr - EMpirical GRamian Framework ( https://gramian.de )
%%% authors: Christian Himpe ( 0000-0003-2194-6754 )
%%% license: BSD-2-Clause (2019)

    if( (nargin<1)  || isempty(ut) ), ut = 1; end

    if(exist('OCTAVE_VERSION','builtin'))
        EMGR = @emgr_oct;
        fprintf('emgr (version: %1.1f%s)\n',EMGR('version'),'-oct');
    elseif(verLessThan('matlab','9.1'))
        EMGR = @emgr_lgc;
        fprintf('emgr (version: %1.1f%s)\n',EMGR('version'),'-lgc');
    else
        EMGR = @emgr;
        fprintf('emgr (version: %1.1f%s)\n',EMGR('version'),'');
    end

%% SYSTEM SETUP
    M = 4;				% number of inputs
    N = M*M;				% number of states
    Q = M;				% number of outputs
    dt = 0.01;				% time step size
    Tf = 1.0;				% time horizon
    X0 = linspace(0,1,N)';		% initial state
    P = 0.5+0.5*cos(1:N)';		% parameter
    R = [0.5*ones(N,1),ones(N,1)];	% parameter range

    A = diag(N*(sin(1:N)-1));		% system matrix
    B = sin(toeplitz(1:N,1:M)) + 1.0;	% input matrix
    C = cos(toeplitz(1:M,1:N)) + 1.0;	% output matrix

    LIN = @(x,u,p,t) A*x + B*u + p;	% vector field
    ADJ = @(x,u,p,t) A'*x + C'*u;	% adjoint vector field
    OUT = @(x,u,p,t) C*x;		% output functional

%% Controllability Gramian

    fprintf('Empirical Controllability Gramian: ');
    W = {};

    % Baseline
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'c',P,[0,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % Centering
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'c',P,[1,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'c',P,[2,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'c',P,[3,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'c',P,[4,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'c',P,[5,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'c',P,[6,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % Input scales
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'c',P,[0,1,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'c',P,[0,2,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'c',P,[0,3,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#'); % fast drop
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'c',P,[0,4,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % Input Rotations
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'c',P,[0,0,0,1,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % Normalization
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'c',P,[0,0,0,0,0,1,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'c',P,[0,0,0,0,0,2,0,0,0,0,0,0],ut,0,X0)); fprintf('#');

    plotsingvals(W,[2,6,1],'Controllability Gramian',10); fprintf('\n');

%% Observability Gramian

    fprintf('Empirical Observability Gramian:   ');
    W = {};

    % Baseline
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'o',P,[0,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % Centering
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'o',P,[1,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'o',P,[2,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'o',P,[3,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'o',P,[4,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'o',P,[5,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'o',P,[6,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % State scales
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'o',P,[0,0,1,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'o',P,[0,0,2,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'o',P,[0,0,3,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#'); % fast drop
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'o',P,[0,0,4,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % State Rotations
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'o',P,[0,0,0,0,1,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % Normalization
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'o',P,[0,0,0,0,0,1,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'o',P,[0,0,0,0,0,2,0,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % Average Observability Gramian
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'o',P,[0,0,0,0,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % Extra Input
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'o',P,[0,0,0,0,0,0,0,1,0,0,0,0],ut,0,X0)); fprintf('#');

    plotsingvals(W,[2,6,2],'Observability Gramian',10); fprintf('\n');

%% Linear Cross Gramian

    fprintf('Empirical Linear Cross Gramian:    ');
    W = {};

    % Baseline
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[0,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % Centering
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[1,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[2,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[3,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[4,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[5,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[6,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % Scales
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[0,1,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[0,2,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[0,3,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#'); % fast drop
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[0,4,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % Rotations
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[0,0,0,1,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % Normalization
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[0,0,0,0,0,1,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[0,0,0,0,0,2,0,0,0,0,0,0],ut,0,X0)); fprintf('#');

    plotsingvals(W,[2,6,3],'Linear Cross Gramian',10);

    W = {};

    % Non-Symmetric Cross Gramian
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[0,0,0,0,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % Centering
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[1,0,0,0,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[2,0,0,0,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[3,0,0,0,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[4,0,0,0,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[5,0,0,0,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[6,0,0,0,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % Scales
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[0,1,0,0,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[0,2,0,0,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[0,3,0,0,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#'); % fast drop
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[0,4,0,0,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % Rotations
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[0,0,0,1,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#'); % deep drop

    % Normalization
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[0,0,0,0,0,1,1,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,ADJ,[M,N,Q],[dt,Tf],'y',P,[0,0,0,0,0,2,1,0,0,0,0,0],ut,0,X0)); fprintf('#');

    plotsingvals(W,[2,6,4],'Linear Cross Gramian (Non-Symmetric)',12); fprintf('\n');

%% Cross Gramian

    fprintf('Empirical Cross Gramian:           ');
    W = {};

    % Baseline
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % Centering
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[1,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[2,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[3,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[4,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[5,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[6,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % Input scales
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,1,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,2,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,3,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,4,0,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % State scales
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,0,1,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,0,2,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,0,3,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,0,4,0,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % Input Rotations
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,0,0,1,0,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % State Rotations
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,0,0,0,1,0,0,0,0,0,0,0],ut,0,X0)); fprintf('#'); % fast drop

    % Normalization
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,0,0,0,0,1,0,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,0,0,0,0,2,0,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % Extra Input
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,0,0,0,0,0,0,1,0,0,0,0],ut,0,X0)); fprintf('#');

    plotsingvals(W,[2,6,5],'Cross Gramian',17);
    W = {};

    % Baseline
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,0,0,0,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % Centering
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[1,0,0,0,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[2,0,0,0,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[3,0,0,0,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[4,0,0,0,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[5,0,0,0,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[6,0,0,0,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % Input scales
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,1,0,0,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,2,0,0,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,3,0,0,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,4,0,0,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % State scales
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,0,1,0,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,0,2,0,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,0,3,0,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,0,4,0,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % Input Rotations
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,0,0,1,0,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#');

    % State Rotations
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,0,0,0,1,0,1,0,0,0,0,0],ut,0,X0)); fprintf('#'); % fast drop

    % Normalization
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,0,0,0,0,1,1,0,0,0,0,0],ut,0,X0)); fprintf('#');
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,0,0,0,0,2,1,0,0,0,0,0],ut,0,X0)); fprintf('#'); 

    % Extra Input
    W{end+1} = svd(EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'x',P,[0,0,0,0,0,0,1,1,0,0,0,0],ut,0,X0)); fprintf('#');

    plotsingvals(W,[2,6,6],'Cross Gramian (Non-Symmetric)',17); fprintf('\n');

%% Sensitivity Gramian

    fprintf('Empirical Sensitivity Gramian:     ');
    W = {};

    % Baseline
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[0,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');

    % Centering
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[1,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[2,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[3,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[4,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[5,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[6,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');

    % Input scales
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[0,1,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[0,2,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[0,3,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[0,4,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');

    % Input Rotations
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[0,0,0,1,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');

    % Extra Input
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[0,0,0,0,0,0,0,1,0,0,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');

    % Parameter Scaling
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[0,0,0,0,0,0,0,0,1,0,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#'); % fast drop
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[0,0,0,0,0,0,0,0,2,0,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');

    plotsingvals(W,[2,6,9],'Sensitivity Gramian',14);
    W = {};

    % Baseline
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[0,0,0,0,0,0,0,0,0,1,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');

    % Centering
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[1,0,0,0,0,0,0,0,0,1,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[2,0,0,0,0,0,0,0,0,1,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[3,0,0,0,0,0,0,0,0,1,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[4,0,0,0,0,0,0,0,0,1,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[5,0,0,0,0,0,0,0,0,1,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[6,0,0,0,0,0,0,0,0,1,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');

    % Input scales
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[0,1,0,0,0,0,0,0,0,1,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[0,2,0,0,0,0,0,0,0,1,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[0,3,0,0,0,0,0,0,0,1,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[0,4,0,0,0,0,0,0,0,1,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');

    % State scales
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[0,0,1,0,0,0,0,0,0,1,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[0,0,2,0,0,0,0,0,0,1,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[0,0,3,0,0,0,0,0,0,1,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[0,0,4,0,0,0,0,0,0,1,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');

    % Input Rotations
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[0,0,0,1,0,0,0,0,0,1,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');

    % State Rotations
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[0,0,0,0,1,0,0,0,0,1,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');

    % Extra Input
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[0,0,0,0,0,0,0,1,0,1,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');

    % Parameter Scaling
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[0,0,0,0,0,0,0,0,1,1,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'s',R,[0,0,0,0,0,0,0,0,2,1,0,0],ut,0,X0); W{end+1} = sort(w{2},'descend'); fprintf('#'); % fast drop

    plotsingvals(W,[2,6,10],'Sensitivity Gramian (Input-Output)',20); fprintf('\n');

%% Identifiability Gramian

    fprintf('Empirical Identifiability Gramian: ');
    W = {};

    % Baseline
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'i',R,[0,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');

    % Centering
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'i',R,[1,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'i',R,[2,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'i',R,[3,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'i',R,[4,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'i',R,[5,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'i',R,[6,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');

    % State scales
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'i',R,[0,0,1,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'i',R,[0,0,2,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'i',R,[0,0,3,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'i',R,[0,0,4,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');

    % State Rotations
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'i',R,[0,0,0,0,1,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');

    % Averaged Identifiability
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'i',R,[0,0,0,0,0,0,1,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');

    % Extra Input
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'i',R,[0,0,0,0,0,0,0,1,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#'); % fast drop

    % Parameter Scaling
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'i',R,[0,0,0,0,0,0,0,0,1,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#'); % fast drop (r only)
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'i',R,[0,0,0,0,0,0,0,0,2,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#'); % fast drop (r only)

    % Schur Complement
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'i',R,[0,0,0,0,0,0,0,0,0,1,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');

    plotsingvals(W,[2,6,8],'Identifiability Gramian',13); fprintf('\n');

%% Joint Gramian

    fprintf('Empirical Joint Gramian:           ');
    W = {};

    % Baseline
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');

    % Centering
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[1,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[2,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[3,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[4,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[5,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[6,0,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');

    % Input scales
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,1,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,2,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,3,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,4,0,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');

    % State scales
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,0,1,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,0,2,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,0,3,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,0,4,0,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');

    % Input Rotations
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,0,0,1,0,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');

    % State Rotations
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,0,0,0,1,0,0,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');

    % Extra Input
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,0,0,0,0,0,0,1,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#'); % fast drop

    % Parameter Scaling
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,0,0,0,0,0,0,0,1,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,0,0,0,0,0,0,0,2,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');

    % Schur Complement
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,0,0,0,0,0,0,0,0,1,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');

    plotsingvals(W,[2,6,11],'Cross-Identifiability Gramian',18);
    W = {};

    % Baseline
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,0,0,0,0,0,1,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');

    % Centering
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[1,0,0,0,0,0,1,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[2,0,0,0,0,0,1,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[3,0,0,0,0,0,1,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[4,0,0,0,0,0,1,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[5,0,0,0,0,0,1,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[6,0,0,0,0,0,1,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');

    % Input scales
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,1,0,0,0,0,1,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,2,0,0,0,0,1,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,3,0,0,0,0,1,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,4,0,0,0,0,1,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');

    % State scales
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,0,1,0,0,0,1,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,0,2,0,0,0,1,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,0,3,0,0,0,1,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,0,4,0,0,0,1,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');

    % Input Rotations
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,0,0,1,0,0,1,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');

    % State Rotations
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,0,0,0,1,0,1,0,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');

    % Extra Input
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,0,0,0,0,0,1,1,0,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#'); % fast drop

    % Parameter Scaling
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,0,0,0,0,0,1,0,1,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,0,0,0,0,0,1,0,2,0,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');

    % Schur Complement
    w = EMGR(LIN,OUT,[M,N,Q],[dt,Tf],'j',R,[0,0,0,0,0,0,1,0,0,1,0,0],ut,0,X0); W{end+1} = svd(w{2}); fprintf('#');

    plotsingvals(W,[2,6,12],'Cross-Identifiability Gramian (Non-symmetric)',18); fprintf('\n\n');
end

function plotsingvals(W,s,l,b)

    if(s(3)==1), figure; end;
    subplot(s(1),s(2),s(3))
    c = winter(numel(W)-1);
    semilogy(W{2}./max(W{2}),'Linewidth',2,'color',c(1,:));
    hold on;
    for k=3:numel(W)

        semilogy(W{k}./max(W{k}),'Linewidth',2,'color',c(k-1,:));
    end
    semilogy(W{1}./max(W{1}),'Linewidth',2,'color',[1,0,0]);
    semilogy(W{b}./max(W{b}),'Linewidth',2,'color',[1,0,1]);
    hold off;
    xlim([1,numel(W{1})]);
    yl = ylim();
    ylim([min(yl(1),0.1),1]);
    xlabel(l)
    set([gca; findall(gca,'Type','text')],'FontSize',6);
end