function curios(sys,task,method,options)
%%% summary: curios - Clearing Up Reducibility of Input-Output Systems
%%% project: emgr - EMpirical GRamian Framework ( https://gramian.de )
%%% authors: Christian Himpe ( 0000-0003-2194-6754 )
%%% license: BSD-2-Clause (2019)

    global ODE;
    global STAGES;

    persistent EV;
    persistent ev;
    persistent SV;

%% State-Space Structure

    if(isobject(sys))

        sys = struct(sys);

        if(isfield(sys,'A')), [sys.a] = sys.A; end;
        if(isfield(sys,'B')), [sys.b] = sys.B; end;
        if(isfield(sys,'C')), [sys.c] = sys.C; end;
        if(isfield(sys,'D')), [sys.d] = sys.D; end;
        if(isfield(sys,'E')), [sys.e] = sys.E; end;

        sys.M = size(sys.b,2);
        sys.N = size(sys.a,2);
        sys.Q = size(sys.c,1);
        sys.dt = 1.0/sys.N;
        sys.Tf = 1.0;

        if(isfield(sys,'e') && not(isempty(sys.e)))
            [L,U,P] = lu(sys.e);
            sys.f = @(x,u,p,t) L\(U\(P*(sys.a*x + sys.b*u)));
            sys.F = @(x,u,p,t) L\(U\(P*(sys.a*x))) + sys.c'*u;
        else
            sys.f = @(x,u,p,t) sys.a*x + sys.b*u;
            sys.F = @(x,u,p,t) sys.a'*x + sys.c'*u;
        end

        if(not(isfield(sys,'d')) || isempty(sys.d) || sys.d==0)
            sys.g = @(x,u,p,t) sys.c*x;
        else
            sys.g = @(x,u,p,t) sys.c*x + sys.d*u;
        end
    end

%% Default Arguments

    in = ones(sys.M,1);

    if(not(isfield(sys,'g'))),  sys.g = 1; end
    if(not(isfield(sys,'p'))),  sys.p = 0; end
    if(not(isfield(sys,'q'))),  sys.q = 0; end
    if(not(isfield(sys,'ut'))), sys.ut = 'i'; end
    if(not(isfield(sys,'vt'))), sys.vt = @(t) in*((t<=sys.dt)./sys.dt); end
    if(not(isfield(sys,'us'))), sys.us = 0; end
    if(not(isfield(sys,'xs'))), sys.xs = zeros(sys.N,1); end
    if(not(isfield(sys,'um'))), sys.um = 1; end
    if(not(isfield(sys,'xm'))), sys.xm = 1; end

%% Argument Check

    task = lower(task);
    method = lower(method);

    if(nargin<4 || isempty(options)), options = {'none'}; end

%% Internal Argument setup

    sysdim = [sys.M,sys.N,sys.Q];
    tdisc = [sys.dt,sys.Tf];
    config = zeros(1,12);
    proj = '';
    dp = @mtimes;

%% Utility Library

    picked = @(name) any(strcmp(options,name));
    gtimes = @(m) m*m';
    rcumsum = @(v) flipud(cumsum(flipud(v(:))));

%% Input Library:

    if(picked('step')),   sys.ut = 's'; end
    if(picked('chirp')),  sys.ut = 'c'; end
    if(picked('random')), sys.ut = 'r'; end

%% Kernel Library:

    DEG = 2;

    % Diagonal-ony Pseudo-kernel
    o_diag = @(x,y) sum(x.*y',2);

    % Trace-only pseudo kernel
    o_trac = @(x,y) sum(sum(x.*y'));

    % Time-weighted kernel
    if(picked('tweighted')),  dp = @(x,y) bsxfun(@times,x,(sys.dt * [1:1:size(x,2)]).^DEG) * y; end

    % Polynomial kernel
    if(picked('polynomial')), dp = @(x,y) (x*y).^DEG + 1.0; end

    % Sigmoid kernel
    if(picked('sigmoid')),    dp = @(x,y) tanh(x*y) + 1.0; end

    % Gaussian kernel
    if(picked('gauss')),      dp = @(x,y) exp(-gtimes(x - y')); end

    % Second-order position kernel
    if(picked('position')),   dp = @(x,y) x(1:sys.N/2,:)*y(:,1:sys.N/2); proj = 'secondo'; end

    % Second-order velocity kernel
    if(picked('velocity')),   dp = @(x,y) x(sys.N/2+1:end,:)*y(:,sys.N/2+1:end); proj = 'secondo'; end

%% Solver Library

    if(picked('rk1e')),       STAGES = 1; end
    if(picked('rk45')),       ODE = @rk45e; end

%% Configuration Library

    if(picked('jacobi')),   config(6) = 1; end
    if(picked('scaled')),   config(6) = 2; end
    if(picked('nonsym')),   config(7) = 1; end
    if(picked('active')),   config(8) = 1; end
    if(picked('linpar')),   config(9) = 1; end
    if(picked('logpar')),   config(9) = 2; end
    if(picked('coarse')),   config(10) = 1; end

    if(picked('steady')),   config(1) = 1; end
    if(picked('final')),    config(1) = 2; end
    if(picked('mean')),     config(1) = 3; end
    if(picked('rms')),      config(1) = 4; end
    if(picked('midr')),     config(1) = 5; end
    if(picked('gmean')),     config(1) = 6; end

    if(picked('linear')),   config(2:3) = 1; end
    if(picked('log')),      config(2:3) = 2; end
    if(picked('geom')),     config(2:3) = 3; end
    if(picked('sparse')),   config(2:3) = 4; end

    % Check for subplot configuration
    sc = cellfun(@(c) isnumeric(c) && numel(c)==3,options);
    if(any(sc))
        subpl = options{sc};
        options(sc) = [];
        if(isempty(options)), options = {'none'}; end
    else
        subpl = [];
    end

    if(picked('hold')), holding = 1; else, holding = 0; end

%% Backend Setup

    [EMGR,mdf] = sel();

%% Print Welcome

    fprintf('\n');
    disp('======== curios - Clearing Up Reducibility of Input-Output Systems ========');
    fprintf(' backend: emgr %.1f%s\n',EMGR('version'),mdf)
    fprintf(' %s: %s \n',task,method);
    fprintf(' options: ');
    for k=1:numel(options)-1, fprintf([options{k},', ']); end
    fprintf([options{end},'\n']);
    fprintf(' system dims: %d input(s), %d states, %d output(s) \n',sys.M,sys.N,sys.Q); 

%% Main body

    tic;

    switch(task)

%% State Reduction

        case 'state-reduction'

            switch(method)

                case 'controllability-truncation'

                    WC = EMGR(sys.f,sys.g,sysdim,tdisc,'c',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm,dp);
                    [UX,DX,VX] = svd(WC);

                case 'observability-truncation'

                    WO = EMGR(sys.f,sys.g,sysdim,tdisc,'o',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm,dp);
                    [UX,DX,VX] = svd(WO);

                case 'linear-direct-truncation'

                    WX = EMGR(sys.f,sys.F,sysdim,tdisc,'y',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm,dp);
                    [UX,DX,VX] = balco(WX);

                case 'nonlinear-direct-truncation'

                    if(picked('partitioned'))
                        config(11) = ceil(sys.N/8);
                        K = ceil(sys.N/config(11));
                        for k=1:K
                            config(12) = k;
                            wx{k} = EMGR(sys.f,sys.g,sysdim,tdisc,'x',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm,dp);
                        end
                        WX = cell2mat(wx);
                    else
                        WX = EMGR(sys.f,sys.g,sysdim,tdisc,'x',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm,dp);
                    end
                    [UX,DX,VX] = balco(WX);

                case 'linear-balanced-truncation'

                    WC = EMGR(sys.f,1,sysdim,tdisc,'c',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm,dp);
                    WO = EMGR(sys.F,1,sysdim,tdisc,'c',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm,dp);
                    [UX,DX,VX] = balco(WC,WO);

                case 'nonlinear-balanced-truncation'

                    WC = EMGR(sys.f,sys.g,sysdim,tdisc,'c',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm,dp);
                    WO = EMGR(sys.f,sys.g,sysdim,tdisc,'o',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm,dp);
                    [UX,DX,VX] = balco(WC,WO);

                case 'linear-dominant-subspaces'
                    WX = EMGR(sys.f,sys.F,sysdim,tdisc,'y',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm,dp);
                    [UX,DX,VX] = svd(WX);
                    UX = UX*DX;
                    VX = VX*DX;
                    proj = 'dominant';

                case 'nonlinear-dominant-subspaces'
                    WX = EMGR(sys.f,sys.g,sysdim,tdisc,'x',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm,dp);
                    [UX,DX,VX] = svd(WX);
                    UX = UX*DX;
                    VX = VX*DX;
                    proj = 'dominant';

                otherwise
                    error('curios: unknown state-reduction method.');
            end

            if(picked('gains')), [UX,DX,VX] = gains(sys,UX,DX,VX); end

            fprintf(' offline time: %.2f s\n',toc);
            [ord,err,nam] = assess(sys,UX,VX,[],[],proj);
            plot_error_1d(ord,err,nam,method,subpl);

            if(not(picked('noscore'))), morscore(ord{1},err{3},sys.N); end;

%% Parameter Reduction

        case 'parameter-reduction'

            switch(method)

                case 'controllability-based'

                    WI = EMGR(sys.f,sys.g,sysdim,tdisc,'s',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm,dp);

                case 'observability-based'

                    WI = EMGR(sys.f,sys.g,sysdim,tdisc,'i',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm,dp);

                case 'minimality-based'

                    if(picked('partitioned'))
                        config(11) = ceil((sys.N+size(sys.p,1))/8);
                        K = ceil((sys.N+size(sys.p,1))/config(11));
                        for k=1:K
                            config(12) = k;
                            wi{1,k} = EMGR(sys.f,sys.g,sysdim,tdisc,'j',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm,dp);
                        end
                        WI{1} = cell2mat(wi(1:ceil(sys.N/config(11))));
                        WI{2} = cell2mat(wi(ceil(sys.N/config(11))+1:end));
                        if(config(10))
                            WI{2} = -0.5*WI{2}'*WI{2};
                        else
                            WI{2} = -0.5*WI{2}'*ainv(WI{1}+WI{1}')*WI{2};
                        end
                    else
                        WI = EMGR(sys.f,sys.g,sysdim,tdisc,'j',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm,dp);
                    end

                otherwise
                    error('curios: unknown parameter-reduction method.');
            end

            [UP,DP,VP] = svd(WI{2});
            fprintf(' offline time: %.2f s\n',toc);
            [ord,err,nam] = assess(sys,[],[],UP,UP,proj);
            plot_error_1d(ord,err,nam,method,subpl);

            if(not(picked('noscore'))), morscore(ord{2},err{3},size(sys.p,1)); end;

%% Combined State and Parameter Reduction

        case 'combined-reduction'

            switch(method)

                case 'controllability-based'

                    WI = EMGR(sys.f,sys.g,sysdim,tdisc,'s',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm,dp);
                    WO = EMGR(sys.F,sys.g,sysdim,tdisc,'c',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm,dp);
                    [UX,DX,VX] = balco(WI{1},WO);

                case 'observability-based'

                    WC = EMGR(sys.f,sys.g,sysdim,tdisc,'c',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm,dp);
                    WI = EMGR(sys.f,sys.g,sysdim,tdisc,'i',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm,dp);
                    [UX,DX,VX] = balco(WC,WI{1});

                case 'minimality-based'

                    if(picked('partitioned'))
                        config(11) = ceil((sys.N+size(sys.p,1))/8);
                        K = ceil((sys.N+size(sys.p,1))/config(11));
                        for k=1:K
                            config(12) = k;
                            wi{1,k} = EMGR(sys.f,sys.g,sysdim,tdisc,'j',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm,dp);
                        end
                        WI{1} = cell2mat(wi(1:ceil(sys.N/config(11))));
                        WI{2} = cell2mat(wi(ceil(sys.N/config(11))+1:end));
                        if(config(10))
                            WI{2} = -0.5*WI{2}'*WI{2};
                        else
                            WI{2} = -0.5*WI{2}'*ainv(WI{1}+WI{1}')*WI{2};
                        end
                    else
                        WI = EMGR(sys.f,sys.g,sysdim,tdisc,'j',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm,dp);
                    end
                    if(picked('dominant'))
                        [UX,DX,VX] = svd(WI{1});
                        UX = UX*DX;
                        VX = VX*DX;
                        proj = 'dominant';
                    else
                        [UX,DX,VX] = balco(WI{1});
                    end

                otherwise
                    error('curios: unknown combined-reduction method.');
            end

            [UP,DP,VP] = svd(WI{2});
            fprintf(' offline time: %.2f s\n',toc);
            [ord,err,nam] = assess(sys,UX,VX,UP,UP,proj);
            plot_error_2d(ord,err,nam,method,subpl);

%% Sensitivity Analysis

        case 'sensitivity-analysis'

            switch(method)

                case 'input-state-based'

                    WS = EMGR(sys.f,sys.g,sysdim,tdisc,'s',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm);

                case 'input-output-based'

                    config(10) = 1;
                    WS = EMGR(sys.f,sys.g,sysdim,tdisc,'s',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm);

                otherwise
                    error('curios: unknown sensitivity-analysis method.');
            end

            fprintf(' offline time: %.2f s\n',toc);
            plot_mag_1d(WS{2}./max(WS{2}),method,'sensitivity',holding,subpl);

%% Parameter Identification

        case 'parameter-identification'

            switch(method)

                case 'state-output-based'

                    WI = EMGR(sys.f,sys.g,sysdim,tdisc,'i',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm,dp);

                case 'input-output-based'

                    WI = EMGR(sys.f,sys.g,sysdim,tdisc,'j',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm,dp);

                otherwise
                    error('curios: unknown parameter-identification method.');
            end

            EO = sort(abs(eig(WI{2})),'descend');
            fprintf(' offline time: %.2f s\n',toc);
            plot_mag_1d(abs(EO)./max(EO),method,'identifiability',holding,subpl);

%% Decentralized Control

        case 'decentralized-control'

            PM = zeros(sys.M,sys.Q);
            for m = 1:sys.M
                u0 = sparse(m,1,1,sys.M,1);
                f = @(x,u,p,t) sys.f(x,u0*u,p,t);
                for q = 1:sys.Q
                    y0 = sparse(1,q,1,1,sys.Q);

                    switch(method)

                        case 'linear'
                        F = @(x,u,p,t) sys.F(x,y0'*u,p,t);
                        WXij = EMGR(f,F,[1,sys.N,1],tdisc,'y',sys.p,config,1,0,sys.xs,1,sys.xm,o_diag);

                        case 'nonlinear'
                        g = @(x,u,p,t) y0*sys.g(x,u0*u,p,t);
                        WXij = EMGR(f,g,[1,sys.N,1],tdisc,'x',sys.p,config,1,0,sys.xs,1,sys.xm,o_diag);

                        otherwise
                            error('curios: unknown decentralized-control method.');
                    end

                    PM(m,q) = sum(WXij.^2);
                    if(picked('coherence')), PM(m,q) = sum(WXij)^2/PM(m,q); end
                end
            end

            PM = PM - min(PM(:));
            PM = PM./max(PM(:));
            fprintf(' offline time: %.2f s\n',toc);
            plot_iomat(PM,sys.M,sys.Q,method,subpl);

%% Nonlinearity Quantification

        case 'nonlinearity-quantification'

            K = 20;
            nl = zeros(1,K);
            sc = linspace(0.1,2.0,K);

            switch(method)

                case 'input-based'

                    for k = 1:K
                        nl(k) = EMGR(sys.f,sys.g,sysdim,tdisc,'c',sys.p,config,sys.ut,sys.us,sys.xs,sc(k)*sys.um,sys.xm,o_trac);
                    end

                case 'state-based'

                    for k = 1:K
                        nl(k) = EMGR(sys.f,sys.g,sysdim,tdisc,'x',sys.p,config,sys.ut,sys.us,sys.xs,sc(k)*sys.um,sc(k)*sys.xm,o_trac);
                    end

                case 'output-based'

                    for k = 1:K
                        nl(k) = EMGR(sys.f,sys.g,sysdim,tdisc,'o',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sc(k)*sys.xm,o_trac);
                    end

                otherwise
                    error('curios: unknown nonlinearity-quantification method.');
            end

            nl = abs(nl);
            nl = nl - min(nl);
            nl = nl./max(nl);
            fprintf(' offline time: %.2f s\n',toc);
            plot_mag_1d(nl,method,'nonlinearity',holding,subpl);

%% State Index

        case 'state-index'

            switch(method)

                case 'controllability'

                    w = EMGR(sys.f,sys.g,sysdim,tdisc,'c',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm,o_diag);

                case 'observability'

                    w = EMGR(sys.f,sys.g,sysdim,tdisc,'o',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm,o_diag);

                case 'minimality'

                    w = EMGR(sys.f,sys.g,sysdim,tdisc,'x',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm,o_diag);

                otherwise
                    error('curios: unknown state-index method.');
            end

            fprintf(' offline time: %.2f s\n',toc);
            plot_mag_1d(abs(w)./max(abs(w)),method,task,holding,subpl);

%% System Indices

        case 'system-index'

            if(not(picked('cached')) || isempty(EV))
                WX = EMGR(sys.f,sys.g,sysdim,tdisc,'x',sys.p,config,sys.ut,sys.us,sys.xs,sys.um,sys.xm,dp);
                EV = eig(WX);
                [~,IX] = sort(abs(EV),'descend');
                EV = real(EV(IX));
                ev = EV;
                ev(abs(ev)<(eps*max(abs(ev)))) = 0;
                SV = svd(WX);
            end

            switch(method)

                case 'cauchy-index'

                    in = sum(sign(ev)) - cumsum(sign(ev));

                case 'system-entropy'

                    in = 1.0 - ([1:sys.N]'*log(2*pi*exp(1)) + cumsum(log(abs(EV))))./(sys.N*log(2*pi*exp(1)) + sum(log(abs(EV))));

                case 'system-gain'

                    in = 1.0 - cumsum(EV)./sum(EV);
                    in = in./max(in);

                case 'hinf-bound'

                    in = 2.0 * rcumsum([abs(EV(2:end));0])./sum(abs(EV));
                    in = in./max(in);

                case 'hankel-bound'

                    in = [abs(EV(2:end));0]./sum(abs(EV));
                    in = in./max(in);

                case 'energy-fraction'

                    in = 1.0 - cumsum(SV)./sum(SV);
                    in = in./max(in);

                case 'storage-efficiency'

                    in = cumprod(abs(EV)).^(1.0./[1:numel(EV)]');
                    in = in./max(in);

                case 'ellipsoid-volume'

                    in = sqrt(cumprod(abs(EV)));
                    in = max(1e-16,in./max(in));

                case 'nyquist-area' % = System-Norm

                    in = 1.0 - sqrt(cumsum(EV.^2) ./ sum(EV.^2));
                    in = in./max(in);

                case 'io-coherence'

                    in = 1.0 - abs((cumsum(EV).^2 * sum(EV.^2)) ./ (cumsum(EV.^2) * sum(EV).^2));
                    in = in./max(in);

                otherwise
                    error('curios: unknown system-index method.');
            end

            fprintf(' offline time: %.2f s\n',toc);
            plot_mag_1d(in,method,task,holding,subpl);

        otherwise

            error('Unknown Task!');
    end

%% Clean up

    ODE = [];
    STAGEC = [];

    fprintf('\n');
end

%% LOCAL FUNCTION: select emgr backend
function [EMGR,mdf] = sel()

    if(exist('OCTAVE_VERSION','builtin'))
        assert(exist('emgr_oct')==2,'emgr_oct not found! Get emgr at: https://gramian.de');
        EMGR = @emgr_oct;
        mdf = '-oct';
    elseif(verLessThan('matlab','9.1'))
        assert(exist('emgr_lgc')==2,'emgr_lgc not found! Get emgr at: https://gramian.de');
        EMGR = @emgr_lgc;
        mdf = '-lgc';
    else

        assert(exist('emgr')==2,'emgr not found! Get emgr at: https://gramian.de');
        EMGR = @emgr;
        mdf = '';
    end
end

%% LOCAL FUNCTION: ainv (approximate inverse)
function x = ainv(m)

    d = diag(m);
    k = find(abs(d)>sqrt(eps));
    d(k) = 1.0./d(k);
    x = bsxfun(@times,m,-d);
    x = bsxfun(@times,x,d');
    x(1:numel(d)+1:end) = d;
end

%% LOCAL FUNCTION: balco (Balancing Controllability and Observability)
function [TR,HSV,TL] = balco(WC,WO)

    switch(nargin)

        case 1
            [LO,EX,LC] = eig(WC,'vector');
            [ex,IN] = sort(sqrt(abs(EX)),'descend');
            LO = LO(:,IN)*diag(ex);
            LC = LC(:,IN)*diag(ex);
        case 2
            [LC,EC] = eig(WC,'vector');
            [ec,IN] = sort(sqrt(abs(EC)),'descend');
            [LO,EO] = eig(WO,'vector');
            [eo,IN] = sort(sqrt(abs(EO)),'descend');
            LC = LC(:,IN)*diag(ec);
            LO = LO(:,IN)*diag(eo);
    end

    [U,HSV,V] = svd(LC'*LO);
    HSV = diag(HSV);
    y = sqrt(1.0./(HSV + 1e-14));
    TR = LO*V*diag(y);
    TL = LC*U*diag(y);
end

%% LOCAL FUNCTION: gains (Balanced Gains Sorting)
function [TR,HSV,TL] = gains(sys,tr,hs,tl)

    B = zeros(sys.N,sys.M);
    for k = 1:sys.M
        B(:,k) = sys.f(sys.xs,(1:sys.M==k)',sys.p,0);
    end
    C = sys.g(1,sys.us,sys.p,0);
    HSV = abs(sum((tr'*B).*(C*tl)',2)).*hs;
    [TMP,IX] = sort(HSV,'descend');
    TL = tl(:,IX);
    TR = tr(:,IX);
end

%% LOCAL FUNCTION: assess (Assess Reduced Order Model)
function [orders,errors,names] = assess(sys,UX,VX,UP,VP,m)

    global ODE;

    P = size(sys.p,1);
    Q = size(sys.q,2);

    if(not(isempty(UX)) && not(isempty(UP)))
        if(sys.N>10), xskip = round(sys.N./10); else, xskip = 1; end
        if(P>10), pskip = round(P./10); else, pskip = 1; end
    else
        if(sys.N>100), xskip = round(sys.N./100); else, xskip = 1; end
        if(P>100), pskip = round(P./100); else, pskip = 1; end
    end

    xtest = [1,xskip:xskip:min(size(UX,2),sys.N)-1];
    ptest = [1,pskip:pskip:min(size(UP,2),P)-1];

    if(isempty(UX) || isempty(VX))
        UX = 1;
        VX = 1;
        xtest = 1;
    end

    if(isempty(UP) || isempty(VP))
        UP = 1;
        VP = 1;
        ptest = 1;
    end

    l0norm = @(y) sum(prod(abs(y),1).^(1.0/size(y,1)),2);
    l1norm = @(y) norm(y(:),1);
    l2norm = @(y) norm(y(:),2);
    l8norm = @(y) norm(y(:),Inf);

    l0 = zeros(numel(xtest),numel(ptest));
    l1 = zeros(numel(xtest),numel(ptest));
    l2 = zeros(numel(xtest),numel(ptest));
    l8 = zeros(numel(xtest),numel(ptest));

    for q = 1:Q

        Y = ODE(sys.f,sys.g,[sys.dt,sys.Tf],sys.xs,@(t) sys.us + sys.vt(t),sys.q(:,q));
        n0 = l0norm(Y);
        n1 = l1norm(Y);
        n2 = l2norm(Y);
        n8 = l8norm(Y);

        for n = xtest

            switch(m)

                case 'dominant'
                    [ux,~,~] = svd([UX(:,1:n),VX(:,1:n)],'econ');
                    ux = ux(:,1:n);
                    vx = ux';

                case 'secondo'
                    nn = numel(sys.xs)/2;
                    ux = [UX(:,1:n),zeros(nn,n);zeros(nn,n),UX(:,1:n)];
                    vx = [VX(:,1:n),zeros(nn,n);zeros(nn,n),VX(:,1:n)]';

                otherwise
                    ux = UX(:,1:n);
                    vx = VX(:,1:n)';
            end

            ix = find(xtest==n);

            for p = ptest

                up = UP(:,1:p);
                vp = VP(:,1:p)';

                ip = find(ptest==p);

                y = ODE(@(x,u,p,t) vx*sys.f(ux*x,u,up*p,t), ...
                        @(x,u,p,t)    sys.g(ux*x,u,up*p,t), ...
                        [sys.dt,sys.Tf], vx*sys.xs, @(t) sys.us + sys.vt(t), vp*sys.q(:,q));

                e0 = l0norm(Y - y) / n0;
                e1 = l1norm(Y - y) / n1;
                e2 = l2norm(Y - y) / n2;
                e8 = l8norm(Y - y) / n8;

                if(Q > 1)
                    e0 = e0 .* e0;
                    e1 = e1 .* e1;
                    e2 = e2 .* e2;
                    e8 = e8 .* e8;
                end

                l0(ix,ip) = l0(ix,ip) + e0;
                l1(ix,ip) = l1(ix,ip) + e1;
                l2(ix,ip) = l2(ix,ip) + e2;
                l8(ix,ip) = l8(ix,ip) + e8;
            end
        end
    end

    if(Q>1)
        l0 = sqrt(l0 ./ Q);
        l1 = sqrt(l1 ./ Q);
        l2 = sqrt(l2 ./ Q);
        l8 = sqrt(l8 ./ Q);
    end

    l0 = min(l0,1.0);
    l1 = min(l1,1.0);
    l2 = min(l2,1.0);
    l8 = min(l8,1.0);

    orders = {xtest, ptest};
    errors = {l0, l1, l2, l8};
    names = {'L_1 Error', 'L_2 Error', 'L_\infty Error', 'L_0 Error'};
end

%% LOCAL FUNCTION: morscore (Model Reduction Scoring)
function morscore(orders,errors,N)

    nx = orders./(N-1);
    ny = log10(errors)./16.0 + 1.0;
    if(nx(end)~=1), nx(end+1) = 1; ny(end+1) = ny(end); end;
    ms = 1.0 - trapz(nx(:),ny(:));
    fprintf(' MOR score: %.2f',ms);
    text(0.5,0.5,num2str(ms,'%.2f'),'Units','normalized');
end

%% LOCAL FUNCTION: rk45e (Embedded Runge Kutta 4th / 5th Solver)
function y = rk45e(f,g,t,x0,u,p)

    [S,x] = ode45(@(t,x) f(x,u(t),p,t),[0,t(2)],x0,'InitialStep',t(1)); % Compute State Trajectory

    K = numel(S);
    z = g(x(1,:)',u(S(1)),p,S(1));
    z(end,K) = 0;

    for k = 2:K
        tk = S(k);
        z(:,k) = g(x(k,:)',u(tk),p,tk); % Compute Output Trajectory
    end

    y = interp1(S,z',0:t(1):t(2))';
end

%% LOCAL FUNCTION: plot_error_1d (1D Error Plot)
function plot_error_1d(orders,errors,names,ident,subpl)

    if(orders{1}==1), dom = orders{2}; else, dom = orders{1}; end

    if(isempty(subpl))
        figure('Name',ident,'NumberTitle','off');
    else
        if(subpl(3)==1), figure; end
        subplot(subpl(1),subpl(2),subpl(3));
    end

    semilogy(dom,errors{2},'r','linewidth',2); hold on;
    semilogy(dom,errors{3},'g','linewidth',2);
    semilogy(dom,errors{4},'b','linewidth',2);
    semilogy(dom,errors{1},'k--','linewidth',2); hold off;
    xlim([dom(1),dom(end)]);
    ylim([1e-16,1]);
    pbaspect([2,1,1]);

    if(isempty(subpl))
        xlabel('Reduced Dimension');
        ylabel('Relative Error');
        legend(names{1},names{2},names{3},names{4},'location','northeast');
    else
        set(gca,'Ytick',[]);
        set(gca,'Xtick',[]);
        if(subpl(3)<=subpl(2)), title(ident); end
        set([gca; findall(gca,'Type','text')],'FontSize',4);
    end
    set(gca,'YGrid','on');
end

%% LOCAL FUNCTION: plot_error2d (2D Error Plot)
function plot_error_2d(orders,errors,names,ident,subpl)

    if(isempty(subpl))
        figure('Name',ident,'NumberTitle','off');
    else
        if(subpl(3)==1), figure; end
        subplot(subpl(1),subpl(2),subpl(3));
    end

    h = surf(orders{1},orders{2},min(1.0,errors{3}));
    set(gca,'ZScale','log');
    set(h,'CData',log10(get(h,'CData')));
    set(gca,'CLim',log10(get(gca,'ZLim')));
    view(135,15);

    if(isempty(subpl))
        zl = zlim();
        zlim([zl(1),1]);
        ylabel('State Dimension')
        xlabel('Parameter Dimension');
        zlabel('Relative Error');
    else
        zlim([1e-16,1]);
        set(gca,'Ytick',[]);
        set(gca,'Xtick',[]);
        if(subpl(3)<=subpl(2)), title(ident); end
        set([gca; findall(gca,'Type','text')],'FontSize',4);
    end
end

%% LOCAL FUNCTION: plot_iomat (Plot Participation Matrix)
function plot_iomat(PM,M,Q,ident,subpl)

    if(isempty(subpl))
        figure('Name',ident,'NumberTitle','off');
    else
        if(subpl(3)==1), figure; end
        subplot(subpl(1),subpl(2),subpl(3));
    end

    imagesc(PM);

    if(isempty(subpl))
        colorbar;
        set(gca,'XTick',1:1:M,'YTick',1:1:Q);
    else
        set(gca,'Ytick',[]);
        set(gca,'Xtick',[]);
        set(gca,'Ztick',[]);
        title(ident);
    end
end

%% LOCAL FUNCTION: plot_mag_1d (1D Magnitude Plot)
function plot_mag_1d(mag,name,ident,holding,subpl)

    persistent c;

    if(isempty(c) || (not(holding) && not(isempty(c))) )
        c = 1;
    end

    dom = numel(mag);

    if(holding)
        col = lines(11);
        hold on;
    else
        col = zeros(11,3);

        if(isempty(subpl))
            figure('Name',ident,'NumberTitle','off');
        else
            if(subpl(3)==1), figure; end
            subplot(subpl(1),subpl(2),subpl(3));
        end
    end

    if(strcmp(name,'cauchy-index') || strcmp(name,'system-entropy'))

        plot(1:dom,mag,'Color',col(c,:),'linewidth',2,'DisplayName',name);
    else
        mag(mag==0) = min(mag(mag>0));
        semilogy(1:dom,mag,'Color',col(c,:),'linewidth',2,'DisplayName',name);
    end

    if(holding)
        hold off;
        c = c + 1;
    end

    if(isempty(subpl))
        legend show
        set(legend,'location','northeast');
    else
        title(ident);
    end

    xlim([1,dom]);
    yl = ylim();
    ylim([min([1e-4;yl(1);mag(:)]),max([1;yl(2);mag(:)])]);
    pbaspect([2,1,1]);
    set(gca,'YGrid','on');
end