text:
using Statistics

function Energysignal(data,Ts)
    return Energysignal_=sum(data.^2,dims=1).*Ts
end
function accRMS(data,step_max_perturb)
    return rms_=(sum(data.^2,dims=1)./(step_max_perturb+1)).^(1/2)
end
function accVDV(data,Ts)
    return vdv_=(sum(data.^4,dims=1).*Ts).^(1/4)
end
function crestFactor(data,accRMS)
    return crest_factor_=mapslices(maximum, data, dims=1)./accRMS
end
function convertPolar(data)
    # convert 2D cartesian to 2D polar coordinate
    data_=[sum(data.^2,dims=2).^(1/2) atand.(data[:,1]./data[:,2])]
end

function splitData(data,filename,plot_order)
    data_=[data[filename[i,1]][] for i in 1 : size(plot_order,1)]
    if size(filename,1)>size(plot_order,1)
        for k in 1:size(filename,1)÷size(plot_order,1)-1
            data_=hcat(data_,[data[filename[i,1]][] for i in 1+k*size(plot_order,1) : size(plot_order,1)+k*size(plot_order,1)])
        end
    end

    return data_
end

function splitData_2(data,filename,plot_order,k)
    data_=[data[filename[i,1]][k] for i in 1 : size(plot_order,1)]
    if size(filename,1)>size(plot_order,1)
        for k in 1:size(filename,1)÷size(plot_order,1)-1
            data_=hcat(data_,[data[filename[i,1]][k] for i in 1+k*size(plot_order,1) : size(plot_order,1)+k*size(plot_order,1)])
        end
    end

    return data_
end

function splitDataMTiXsens(data,filename,plot_order,n,k)
    data_=[data[filename[i,1]][n][1,k] for i in 1 : size(plot_order,1)]
    if size(filename,1)>size(plot_order,1)
        for k in 1:size(filename,1)÷size(plot_order,1)-1
            data_=hcat(data_,[data[filename[i,1]][n][1,k] for i in 1+k*size(plot_order,1) : size(plot_order,1)+k*size(plot_order,1)])
        end
    end

    return data_
end

function analysisPULSE(data,Ts,step_max_perturb)
    data_=data.-(sum(data[1:10,1],dims=1)/10)

    Energysignal_=Energysignal(data_,Ts)
    rms_=accRMS(data_,step_max_perturb)
    vdv_=accVDV(data_,Ts)
    crest_factor_=crestFactor(abs.(data_),rms_)
    std_=std(data_,dims=1)
    min_=minimum(data_,dims=1)
    max_=maximum(data_,dims=1)
    range_=max_.-min_

    return Energysignal_, rms_, vdv_, crest_factor_,std_,min_,max_,range_
end

function analysisPULSE_MTiXsens(AccFree_MTiXsens,
                                filename,
                                sensorXsensNames,
                                Ts_MTiXsens,
                                step_max_perturb,
                                f_coupure=15,
                                f_sampling=100)
    Energysignal_=Dict()
    rms_=Dict()
    vdv_=Dict()
    crest_factor_=Dict()

    Energysignal_polar_=Dict()
    rms_polar_=Dict()
    vdv_polar_=Dict()
    # crest_factor_polar_=Dict()
    for experiment=1:size(filename,1)
        Energysignal_[filename[experiment,1]]=Dict()
        rms_[filename[experiment,1]]=Dict()
        vdv_[filename[experiment,1]]=Dict()
        crest_factor_[filename[experiment,1]]=Dict()
        min_[filename[experiment,1]]=Dict()
        max_[filename[experiment,1]]=Dict()
        range_[filename[experiment,1]]=Dict()

        Energysignal_polar_[filename[experiment,1]]=Dict()
        rms_polar_[filename[experiment,1]]=Dict()
        vdv_polar_[filename[experiment,1]]=Dict()
        crest_factor_polar[filename[experiment,1]]=Dict()
        min_polar_[filename[experiment,1]]=Dict()
        max_polar_[filename[experiment,1]]=Dict()
        range_polar_[filename[experiment,1]]=Dict()

        # crest_factor_polar_[filename[experiment,1]]=Dict()
        for n in sensorXsensNames
            data_=convert(Matrix,AccFree_MTiXsens[filename[experiment,1]][n][:,2:4])
            # f_coupure=25
            # f_sampling=100
            designmethod = Butterworth(5)
            ff = digitalfilter(Lowpass(f_coupure/f_sampling),designmethod)
            data_=[filtfilt(ff,data_[:,1]) filtfilt(ff,data_[:,2]) filtfilt(ff,data_[:,3])]

            # acc_temp=convert(Matrix,AccFree_MTiXsens[filename[experiment,1]][n][:,2:4])
            acc_temp=[sqrt.(data_[:,1].^2+data_[:,2].^2) data_[:,3]]

            Energysignal_[filename[experiment,1]][n],
            rms_[filename[experiment,1]][n],
            vdv_[filename[experiment,1]][n],
            crest_factor_[filename[experiment,1]][n],
            min_[filename[experiment,1]][n],
            max_[filename[experiment,1]][n],
            range_[filename[experiment,1]][n]=analysisPULSE(acc_temp,Ts_MTiXsens,step_max_perturb)

            acc_temp_polar=convertPolar(acc_temp)

            Energysignal_polar_[filename[experiment,1]][n],
            rms_polar_[filename[experiment,1]][n],
            vdv_polar_[filename[experiment,1]][n],
            crest_factor_polar_[filename[experiment,1]][n],
            min_polar_[filename[experiment,1]][n],
            max_polar_[filename[experiment,1]][n],
            range_polar_[filename[experiment,1]][n]=analysisPULSE(acc_temp,Ts_MTiXsens,step_max_perturb)

            # Energysignal_polar_[filename[experiment,1]][n]=convertPolar(Energysignal_[filename[experiment,1]][n])
            # rms_polar_[filename[experiment,1]][n]=convertPolar(rms_[filename[experiment,1]][n])
            # vdv_polar_[filename[experiment,1]][n]=convertPolar(vdv_[filename[experiment,1]][n])
            # crest_factor_polar_[filename[experiment,1]][n]=
        end
    end
    return Energysignal_, rms_, vdv_, crest_factor_, Energysignal_polar_, rms_polar_, vdv_polar_
end

function analysisPULSE_All(data,filename,Ts,step_max_perturb,m)
    Energysignal=Dict()
    rms=Dict()
    vdv=Dict()
    crest_factor=Dict()
    std_=Dict()
    min_=Dict()
    max_=Dict()
    range_=Dict()
    for experiment=1:size(filename,1)
       Energysignal[filename[experiment,1]],
       rms[filename[experiment,1]],
       vdv[filename[experiment,1]],
       crest_factor[filename[experiment,1]],
       std_[filename[experiment,1]],
       min_[filename[experiment,1]],
       max_[filename[experiment,1]],
       range_[filename[experiment,1]]=analysisPULSE(data[filename[experiment,1]][:,m],Ts,step_max_perturb)
    end
    return Energysignal, rms, vdv, crest_factor,std_,min_,max_,range_
end

function analysisPULSE_IBS(data,filename,Ts,step_max_perturb,m)
    Energysignal=Dict()
    rms=Dict()
    vdv=Dict()
    crest_factor=Dict()
    std_=Dict()
    min_=Dict()
    max_=Dict()
    range_=Dict()
    for experiment=1:size(filename,1)
       data_=data[filename[experiment,1]][:,m]

       Energysignal[filename[experiment,1]],
       rms[filename[experiment,1]],
       vdv[filename[experiment,1]],
       crest_factor[filename[experiment,1]],
       std_[filename[experiment,1]],
       min_[filename[experiment,1]],
       max_[filename[experiment,1]],
       range_[filename[experiment,1]]=analysisPULSE(data_[data_.<0],Ts,step_max_perturb)
    end
    return Energysignal, rms, vdv, crest_factor,std_,min_,max_,range_
end

function addScatterPlotMTiXsens(pPlot,Energysignal,filename,plot_order,perturbIndices,n,k,perturbLabel,linecolor,perturbshape)
    if isempty(perturbIndices)
        data=splitDataMTiXsens(Energysignal,filename,plot_order,n,k)
        pPlot[n] = scatter!(
             plot_order,
            data,
            label=perturbLabel,
            markershape=perturbshape
        );
    else
        for ind in 1:size(perturbIndices,1)-1
            i1=perturbIndices[ind]
            i2=perturbIndices[ind+1]-1
            data=splitDataMTiXsens(Energysignal,filename[i1:i2],plot_order[i1:i2],n,k)
            pPlot[n] = scatter!(
                plot_order[i1:i2],
                data,
                label=perturbLabel[ind],
                markercolor=(if isempty(linecolor);:match;else;linecolor[ind];end),# markercolor=linecolor[ind],
                markershape=perturbshape[ind]
            );
        end
        ind=size(perturbIndices,1)
        i1=perturbIndices[ind]
        data=splitDataMTiXsens(Energysignal,filename[i1:end],plot_order[i1:end],n,k)
        pPlot[n] = scatter!(
            plot_order[i1:end],
            data,
            label=perturbLabel[ind],
            markercolor=(if isempty(linecolor);:match;else;linecolor[ind];end),# markercolor=linecolor[ind],
            markershape=perturbshape[ind]
        );
    end
end

function addScatterPlot_2(pPlot,Energysignal,filename,plot_order,perturbIndices,k,perturbLabel,linecolor,perturbshape)
    if isempty(perturbIndices)
        data=splitData_2(Energysignal,filename,plot_order,k)
        pPlot = scatter!(
             plot_order,
            data,
            label=perturbLabel,
            markershape=perturbshape
        );
    else
        for ind in 1:size(perturbIndices,1)-1
            i1=perturbIndices[ind]
            i2=perturbIndices[ind+1]-1
            data=splitData_2(Energysignal,filename[i1:i2],plot_order[i1:i2],k)
            pPlot = scatter!(
                plot_order[i1:i2],
                data,
                label=perturbLabel[ind],
                markercolor=(if isempty(linecolor);:match;else;linecolor[ind];end),# markercolor=linecolor[ind],
                markershape=perturbshape[ind]
            );
        end
        ind=size(perturbIndices,1)
        i1=perturbIndices[ind]
        data=splitData_2(Energysignal,filename[i1:end],plot_order[i1:end],k)
        pPlot = scatter!(
            plot_order[i1:end],
            data,
            label=perturbLabel[ind],
            markercolor=(if isempty(linecolor);:match;else;linecolor[ind];end),# markercolor=linecolor[ind],
            markershape=perturbshape[ind]
        );
    end
end

function addScatterPlot(pPlot,Energysignal,filename,plot_order,perturbIndices,perturbLabel,linecolor,perturbshape)
    if isempty(perturbIndices)
        data=splitData(Energysignal,filename,plot_order)
        pPlot = scatter!(
             plot_order,
            data,
            label=perturbLabel,
            markershape=perturbshape
        );
    else
        for ind in 1:size(perturbIndices,1)-1
            i1=perturbIndices[ind]
            i2=perturbIndices[ind+1]-1
            data=splitData(Energysignal,filename[i1:i2],plot_order[i1:i2])
            pPlot = scatter!(
                plot_order[i1:i2],
                data,
                label=perturbLabel[ind],
                markercolor=(if isempty(linecolor);:match;else;linecolor[ind];end),# markercolor=linecolor[ind],
                markershape=perturbshape[ind]
            );
        end
        ind=size(perturbIndices,1)
        i1=perturbIndices[ind]
        data=splitData(Energysignal,filename[i1:end],plot_order[i1:end])
        pPlot = scatter!(
            plot_order[i1:end],
            data,
            label=perturbLabel[ind],
            markercolor=(if isempty(linecolor);:match;else;linecolor[ind];end),# markercolor=linecolor[ind],
            markershape=perturbshape[ind]
        );
    end
end

function plotAnalysisPULSE_MTiXsens(
            Energysignal_MTiXsens,
            rms_MTiXsens,
            vdv_MTiXsens,
            crest_factor_MTiXsens,
            filenameMTiXsens,
            sensorXsensNames,
            axis_name=["x-y";"z"],
            perturbLabel=[],
            plot_order=[],
            perturbIndices=[],
            linecolor= [],
            perturbshape=[],
            addTitle="")

    # plot_order=[1 1 1 2 2 2 3 3 3]'
    # linecolor= [:blue :red];
    # perturbshape=[:circle :utriangle]
    # perturbLabel=["Planche 1" "Planche 2"]
    #each sensor
    for k in 1:size(axis_name,1)
        ####################################################################
        pMTiXsens_SignalEnergy=Dict()
        pMTiXsens_RMS=Dict()
        pMTiXsens_VDV=Dict()
        pMTiXsens_crest_factor=Dict()
        for n in sensorXsensNames
            pMTiXsens_SignalEnergy[n]=plot(
            title = string("MTiXsens signal energy in ",n," along ",axis_name[k],"-axis",addTitle),
            layout = (1, 1),
            # legend = false,
            legend=:best,
            #ylims = (0,30),
            xlabel = "vel [-]",
            ylabel = "Signal energy [m.s^-1]")

            addScatterPlotMTiXsens(pMTiXsens_SignalEnergy,Energysignal_MTiXsens,filenameMTiXsens,plot_order,perturbIndices,n,k,perturbLabel,linecolor,perturbshape)

            ####################################################################
            pMTiXsens_RMS[n]=plot(
            title = string("MTiXsens RMS in ",n," along ",axis_name[k],"-axis",addTitle),
            layout = (1, 1),
            # legend = false,
            legend=:best,
            #ylims = (0,30),
            xlabel = "vel [-]",
            ylabel = "RMS []")

            addScatterPlotMTiXsens(pMTiXsens_RMS,rms_MTiXsens,filenameMTiXsens,plot_order,perturbIndices,n,k,perturbLabel,linecolor,perturbshape)

            ####################################################################
            pMTiXsens_VDV[n]=plot(
            title = string("MTiXsens VDV in ",n," along ",axis_name[k],"-axis",addTitle),
            layout = (1, 1),
            # legend = false,
            legend=:best,
            #ylims = (0,30),
            xlabel = "vel [-]",
            ylabel = "VDV []")

            addScatterPlotMTiXsens(pMTiXsens_VDV,vdv_MTiXsens,filenameMTiXsens,plot_order,perturbIndices,n,k,perturbLabel,linecolor,perturbshape)

            ####################################################################
            pMTiXsens_crest_factor[n]=plot(
            title = string("MTiXsens Crest Factor in ",n," along ",axis_name[k],"-axis",addTitle),
            layout = (1, 1),
            # legend = false,
            legend=:best,
            #ylims = (0,30),
            xlabel = "vel [-]",
            ylabel = "Crest Factor []")

            addScatterPlotMTiXsens(pMTiXsens_crest_factor,crest_factor_MTiXsens,filenameMTiXsens,plot_order,perturbIndices,n,k,perturbLabel,linecolor,perturbshape)

            pMTiXsens_crest_factor[n]=plot!([0.8 3.2]',[9 9]',
                                                linestyle=:dot,
                                                linecolor=:black,
                                                label="threshold")
        end

        #savefig(plot_Acc_MTiXsens_norme[filenameMTiXsens[experiment,1]], string("save_figure/MTiXsens/Dossier_Xsens_", experiment, ".pdf"))
        for n in sensorXsensNames
            savefig(pMTiXsens_SignalEnergy[n],string("save_figure/figure_SE_RMS_VDV/SE_",n,"_",axis_name[k],"axis",".png"))
            savefig(pMTiXsens_RMS[n],string("save_figure/figure_SE_RMS_VDV/RMS_",n,"_",axis_name[k],"axis",".png"))
            savefig(pMTiXsens_VDV[n],string("save_figure/figure_SE_RMS_VDV/VDV_",n,"_",axis_name[k],"axis",".png"))

            savefig(pMTiXsens_crest_factor[n],string("save_figure/figure_SE_RMS_VDV/crest_factor_",n,"_",axis_name[k],"axis",".png"))

        end
    end

    # #tete-siege
    # for k in 1:size(axis_name,1)
    #     ####################################################################
    #     linecolor= [:blue :red];
    #
    #     pMTiXsens_SignalEnergy_tete_siege = scatter(
    #          plot_order,
    #         hcat([Energysignal_MTiXsens[filenameMTiXsens[i,1]][:tete][1,k] for i in 1 : 9],
    #                 [Energysignal_MTiXsens[filenameMTiXsens[i,1]][:tete][1,k] for i in 13 : 21],
    #                 [Energysignal_MTiXsens[filenameMTiXsens[i,1]][:siege][1,k] for i in 1 : 9],
    #                 [Energysignal_MTiXsens[filenameMTiXsens[i,1]][:siege][1,k] for i in 13 : 21]
    #             ),
    #         title = [string("MTiXsens signal energy in tete-siege along ",axis_name[k],"-axis")],
    #         layout = (1, 1),
    #         # legend = false,
    #         legend=:best,
    #         #ylims = (0,30),
    #         xlabel = "vel [-]",
    #         ylabel = "Signal energy [m.s^-1]",
    #         markercolor=linecolor,
    #         label=["Planche 1 tete" "Planche 2 tete" "Planche 1 siege" "Planche 2 siege"],
    #         markershape=[:circle :circle :rect :rect]
    #     );
    #
    #     pMTiXsens_RMS_tete_siege = scatter(
    #          plot_order,
    #         hcat([rms_MTiXsens[filenameMTiXsens[i,1]][:tete][1,k] for i in 1 : 9],
    #                 [rms_MTiXsens[filenameMTiXsens[i,1]][:tete][1,k] for i in 13 : 21],
    #                 [rms_MTiXsens[filenameMTiXsens[i,1]][:siege][1,k] for i in 1 : 9],
    #                 [rms_MTiXsens[filenameMTiXsens[i,1]][:siege][1,k] for i in 13 : 21]
    #             ),
    #         title = [string("MTiXsens RMS in tete-siege along ",axis_name[k],"-axis")],
    #         layout = (1, 1),
    #         # legend = false,
    #         legend=:best,
    #         #ylims = (0,7),
    #         xlabel = "vel [-]",
    #         ylabel = "RMS []",
    #         markercolor=linecolor,
    #         label=["Planche 1 tete" "Planche 2 tete" "Planche 1 siege" "Planche 2 siege"],
    #         markershape=[:circle :circle :rect :rect]
    #     );
    #
    #     pMTiXsens_VDV_tete_siege = scatter(
    #          plot_order,
    #         hcat([vdv_MTiXsens[filenameMTiXsens[i,1]][:tete][1,k] for i in 1 : 9],
    #                 [vdv_MTiXsens[filenameMTiXsens[i,1]][:tete][1,k] for i in 13 : 21],
    #                 [vdv_MTiXsens[filenameMTiXsens[i,1]][:siege][1,k] for i in 1 : 9],
    #                 [vdv_MTiXsens[filenameMTiXsens[i,1]][:siege][1,k] for i in 13 : 21]
    #             ),
    #         title = [string("MTiXsens VDV in tete-siege along ",axis_name[k],"-axis")],
    #         layout = (1, 1),
    #         # legend = false,
    #         legend=:best,
    #         #ylims = (0,8),
    #         xlabel = "vel [-]",
    #         ylabel = "VDV []",
    #         markercolor=linecolor,
    #         label=["Planche 1 tete" "Planche 2 tete" "Planche 1 siege" "Planche 2 siege"],
    #         markershape=[:circle :circle :rect :rect]
    #     );
    #
    #     savefig(pMTiXsens_SignalEnergy_tete_siege,string("save_figure/figure_SE_RMS_VDV/SE_tete-siege_",axis_name[k],"axis"))
    #     savefig(pMTiXsens_RMS_tete_siege,string("save_figure/figure_SE_RMS_VDV/RMS_tete-siege_",axis_name[k],"axis"))
    #     savefig(pMTiXsens_VDV_tete_siege,string("save_figure/figure_SE_RMS_VDV/VDV_tete-siege_",axis_name[k],"axis"))
    # end
    #
    # #coussin-siege
    # for k in 1:size(axis_name,1)
    #     ####################################################################
    #     linecolor= [:blue :red];
    #
    #     pMTiXsens_SignalEnergy_coussin_siege = scatter(
    #          plot_order,
    #         hcat([Energysignal_MTiXsens[filenameMTiXsens[i,1]][:coussin][1,k] for i in 1 : 9],
    #                 [Energysignal_MTiXsens[filenameMTiXsens[i,1]][:coussin][1,k] for i in 13 : 21],
    #                 [Energysignal_MTiXsens[filenameMTiXsens[i,1]][:siege][1,k] for i in 1 : 9],
    #                 [Energysignal_MTiXsens[filenameMTiXsens[i,1]][:siege][1,k] for i in 13 : 21]
    #             ),
    #         title = [string("MTiXsens signal energy in coussin-siege along ",axis_name[k],"-axis")],
    #         layout = (1, 1),
    #         # legend = false,
    #         legend=:best,
    #         #ylims = (0,30),
    #         xlabel = "vel [-]",
    #         ylabel = "Signal energy [m.s^-1]",
    #         markercolor=linecolor,
    #         label=["Planche 1 coussin" "Planche 2 coussin" "Planche 1 siege" "Planche 2 siege"],
    #         markershape=[:circle :circle :rect :rect]
    #     );
    #
    #     pMTiXsens_RMS_coussin_siege = scatter(
    #          plot_order,
    #         hcat([rms_MTiXsens[filenameMTiXsens[i,1]][:coussin][1,k] for i in 1 : 9],
    #                 [rms_MTiXsens[filenameMTiXsens[i,1]][:coussin][1,k] for i in 13 : 21],
    #                 [rms_MTiXsens[filenameMTiXsens[i,1]][:siege][1,k] for i in 1 : 9],
    #                 [rms_MTiXsens[filenameMTiXsens[i,1]][:siege][1,k] for i in 13 : 21]
    #             ),
    #         title = [string("MTiXsens RMS in coussin-siege along ",axis_name[k],"-axis")],
    #         layout = (1, 1),
    #         # legend = false,
    #         legend=:best,
    #         #ylims = (0,7),
    #         xlabel = "vel [-]",
    #         ylabel = "RMS []",
    #         markercolor=linecolor,
    #         label=["Planche 1 coussin" "Planche 2 coussin" "Planche 1 siege" "Planche 2 siege"],
    #         markershape=[:circle :circle :rect :rect]
    #     );
    #
    #     pMTiXsens_VDV_coussin_siege = scatter(
    #          plot_order,
    #         hcat([vdv_MTiXsens[filenameMTiXsens[i,1]][:coussin][1,k] for i in 1 : 9],
    #                 [vdv_MTiXsens[filenameMTiXsens[i,1]][:coussin][1,k] for i in 13 : 21],
    #                 [vdv_MTiXsens[filenameMTiXsens[i,1]][:siege][1,k] for i in 1 : 9],
    #                 [vdv_MTiXsens[filenameMTiXsens[i,1]][:siege][1,k] for i in 13 : 21]
    #             ),
    #         title = [string("MTiXsens VDV in coussin-siege along ",axis_name[k],"-axis")],
    #         layout = (1, 1),
    #         # legend = false,
    #         legend=:best,
    #         #ylims = (0,8),
    #         xlabel = "vel [-]",
    #         ylabel = "VDV []",
    #         markercolor=linecolor,
    #         label=["Planche 1 coussin" "Planche 2 coussin" "Planche 1 siege" "Planche 2 siege"],
    #         markershape=[:circle :circle :rect :rect]
    #     );
    #
    #     savefig(pMTiXsens_SignalEnergy_coussin_siege,string("save_figure/figure_SE_RMS_VDV/SE_coussin-siege_",axis_name[k],"axis"))
    #     savefig(pMTiXsens_RMS_coussin_siege,string("save_figure/figure_SE_RMS_VDV/RMS_coussin-siege_",axis_name[k],"axis"))
    #     savefig(pMTiXsens_VDV_coussin_siege,string("save_figure/figure_SE_RMS_VDV/VDV_coussin-siege_",axis_name[k],"axis"))
    # end
    #
    # #tete-dossier
    # for k in 1:size(axis_name,1)
    #     ####################################################################
    #     linecolor= [:blue :red];
    #
    #     pMTiXsens_SignalEnergy_tete_dossier = scatter(
    #          plot_order,
    #         hcat([Energysignal_MTiXsens[filenameMTiXsens[i,1]][:tete][1,k] for i in 1 : 9],
    #                 [Energysignal_MTiXsens[filenameMTiXsens[i,1]][:tete][1,k] for i in 13 : 21],
    #                 [Energysignal_MTiXsens[filenameMTiXsens[i,1]][:dossier][1,k] for i in 1 : 9],
    #                 [Energysignal_MTiXsens[filenameMTiXsens[i,1]][:dossier][1,k] for i in 13 : 21]
    #             ),
    #         title = [string("MTiXsens signal energy in tete-dossier along ",axis_name[k],"-axis")],
    #         layout = (1, 1),
    #         # legend = false,
    #         legend=:best,
    #         #ylims = (0,30),
    #         xlabel = "vel [-]",
    #         ylabel = "Signal energy [m.s^-1]",
    #         markercolor=linecolor,
    #         label=["Planche 1 tete" "Planche 2 tete" "Planche 1 dossier" "Planche 2 dossier"],
    #         markershape=[:circle :circle :rect :rect]
    #     );
    #
    #     pMTiXsens_RMS_tete_dossier = scatter(
    #          plot_order,
    #         hcat([rms_MTiXsens[filenameMTiXsens[i,1]][:tete][1,k] for i in 1 : 9],
    #                 [rms_MTiXsens[filenameMTiXsens[i,1]][:tete][1,k] for i in 13 : 21],
    #                 [rms_MTiXsens[filenameMTiXsens[i,1]][:dossier][1,k] for i in 1 : 9],
    #                 [rms_MTiXsens[filenameMTiXsens[i,1]][:dossier][1,k] for i in 13 : 21]
    #             ),
    #         title = [string("MTiXsens RMS in tete-dossier along ",axis_name[k],"-axis")],
    #         layout = (1, 1),
    #         # legend = false,
    #         legend=:best,
    #         #ylims = (0,7),
    #         xlabel = "vel [-]",
    #         ylabel = "RMS []",
    #         markercolor=linecolor,
    #         label=["Planche 1 tete" "Planche 2 tete" "Planche 1 dossier" "Planche 2 dossier"],
    #         markershape=[:circle :circle :rect :rect]
    #     );
    #
    #     pMTiXsens_VDV_tete_dossier = scatter(
    #          plot_order,
    #         hcat([vdv_MTiXsens[filenameMTiXsens[i,1]][:tete][1,k] for i in 1 : 9],
    #                 [vdv_MTiXsens[filenameMTiXsens[i,1]][:tete][1,k] for i in 13 : 21],
    #                 [vdv_MTiXsens[filenameMTiXsens[i,1]][:dossier][1,k] for i in 1 : 9],
    #                 [vdv_MTiXsens[filenameMTiXsens[i,1]][:dossier][1,k] for i in 13 : 21]
    #             ),
    #         title = [string("MTiXsens VDV in tete-dossier along ",axis_name[k],"-axis")],
    #         layout = (1, 1),
    #         # legend = false,
    #         legend=:best,
    #         #ylims = (0,8),
    #         xlabel = "vel [-]",
    #         ylabel = "VDV []",
    #         markercolor=linecolor,
    #         label=["Planche 1 tete" "Planche 2 tete" "Planche 1 dossier" "Planche 2 dossier"],
    #         markershape=[:circle :circle :rect :rect]
    #     );
    #
    #     savefig(pMTiXsens_SignalEnergy_tete_dossier,string("save_figure/figure_SE_RMS_VDV/SE_tete-dossier_",axis_name[k],"axis"))
    #     savefig(pMTiXsens_RMS_tete_dossier,string("save_figure/figure_SE_RMS_VDV/RMS_tete-dossier_",axis_name[k],"axis"))
    #     savefig(pMTiXsens_VDV_tete_dossier,string("save_figure/figure_SE_RMS_VDV/VDV_tete-dossier_",axis_name[k],"axis"))
    # end
    #
    # #STHT chest-siege
    # for k in 1:size(axis_name,1)
    #     ####################################################################
    #     linecolor= [:blue :red];
    #
    #     pMTiXsens_SignalEnergy_STHT_chest_siege = scatter(
    #          plot_order,
    #         hcat([Energysignal_MTiXsens[filenameMTiXsens[i,1]][:chest][1,k]/Energysignal_MTiXsens[filenameMTiXsens[i,1]][:siege][1,k] for i in 1 : 9],
    #                 [Energysignal_MTiXsens[filenameMTiXsens[i,1]][:chest][1,k]/Energysignal_MTiXsens[filenameMTiXsens[i,1]][:siege][1,k] for i in 13 : 21]
    #             ),
    #         title = [string("MTiXsens signal energy STHT chest-siege along ",axis_name[k],"-axis")],
    #         layout = (1, 1),
    #         # legend = false,
    #         legend=:best,
    #         #ylims = (0,2),
    #         xlabel = "vel [-]",
    #         ylabel = "Signal energy STHT [-]",
    #         markercolor=linecolor,
    #         label=perturbLabel,
    #         markershape
    #     );
    #
    #     pMTiXsens_VDV_STHT_chest_siege=plot!([0.8 3.2]',[1 1]',
    #                                         linestyle=:dot,
    #                                         linecolor=:black,
    #                                         label="threshold")
    #
    #     pMTiXsens_RMS_STHT_chest_siege = scatter(
    #          plot_order,
    #         hcat([rms_MTiXsens[filenameMTiXsens[i,1]][:chest][1,k]/rms_MTiXsens[filenameMTiXsens[i,1]][:siege][1,k] for i in 1 : 9],
    #                 [rms_MTiXsens[filenameMTiXsens[i,1]][:chest][1,k]/rms_MTiXsens[filenameMTiXsens[i,1]][:siege][1,k] for i in 13 : 21]
    #             ),
    #         title = [string("MTiXsens RMS STHT chest-siege along ",axis_name[k],"-axis")],
    #         layout = (1, 1),
    #         # legend = false,
    #         legend=:best,
    #         #ylims = (0,2),
    #         xlabel = "vel [-]",
    #         ylabel = "RMS STHT []",
    #         markercolor=linecolor,
    #         label=perturbLabel,
    #         markershape
    #     );
    #     pMTiXsens_RMS_STHT_chest_siege=plot!([0.8 3.2]',[1 1]',
    #                                         linestyle=:dot,
    #                                         linecolor=:black,
    #                                         label="threshold")
    #
    #
    #     pMTiXsens_VDV_STHT_chest_siege = scatter(
    #          plot_order,
    #         hcat([vdv_MTiXsens[filenameMTiXsens[i,1]][:chest][1,k]/vdv_MTiXsens[filenameMTiXsens[i,1]][:siege][1,k] for i in 1 : 9],
    #                 [vdv_MTiXsens[filenameMTiXsens[i,1]][:chest][1,k]/vdv_MTiXsens[filenameMTiXsens[i,1]][:siege][1,k] for i in 13 : 21],
    #             ),
    #         title = [string("MTiXsens VDV STHT chest-siege along ",axis_name[k],"-axis")],
    #         layout = (1, 1),
    #         # legend = false,
    #         legend=:best,
    #         #ylims = (0,2),
    #         xlabel = "vel [-]",
    #         ylabel = "VDV STHT []",
    #         markercolor=linecolor,
    #         label=perturbLabel,
    #         markershape
    #     );
    #     pMTiXsens_VDV_STHT_tete_siege=plot!([0.8 3.2]',[1 1]',
    #                                         linestyle=:dot,
    #                                         linecolor=:black,
    #                                         label="threshold")
    #
    #     savefig(pMTiXsens_SignalEnergy_STHT_chest_siege,string("save_figure/figure_SE_RMS_VDV/SE_STHT_chest-siege_",axis_name[k],"axis"))
    #     savefig(pMTiXsens_RMS_STHT_chest_siege,string("save_figure/figure_SE_RMS_VDV/RMS_STHT_chest-siege_",axis_name[k],"axis"))
    #     savefig(pMTiXsens_VDV_STHT_chest_siege,string("save_figure/figure_SE_RMS_VDV/VDV_STHT_chest-siege_",axis_name[k],"axis"))
    # end
    #
    # #STHT chest-dossier
    # for k in 1:size(axis_name,1)
    #     ####################################################################
    #     linecolor= [:blue :red];
    #
    #     pMTiXsens_SignalEnergy_STHT_chest_dossier = scatter(
    #          plot_order,
    #         hcat([Energysignal_MTiXsens[filenameMTiXsens[i,1]][:chest][1,k]/Energysignal_MTiXsens[filenameMTiXsens[i,1]][:dossier][1,k] for i in 1 : 9],
    #                 [Energysignal_MTiXsens[filenameMTiXsens[i,1]][:chest][1,k]/Energysignal_MTiXsens[filenameMTiXsens[i,1]][:dossier][1,k] for i in 13 : 21]
    #             ),
    #         title = [string("MTiXsens signal energy STHT chest-dossier along ",axis_name[k],"-axis")],
    #         layout = (1, 1),
    #         # legend = false,
    #         legend=:best,
    #         #ylims = (0,2),
    #         xlabel = "vel [-]",
    #         ylabel = "Signal energy STHT [-]",
    #         markercolor=linecolor,
    #         label=perturbLabel,
    #         markershape
    #     );
    #
    #     pMTiXsens_VDV_STHT_chest_dossier=plot!([0.8 3.2]',[1 1]',
    #                                         linestyle=:dot,
    #                                         linecolor=:black,
    #                                         label="threshold")
    #
    #     pMTiXsens_RMS_STHT_chest_dossier = scatter(
    #          plot_order,
    #         hcat([rms_MTiXsens[filenameMTiXsens[i,1]][:chest][1,k]/rms_MTiXsens[filenameMTiXsens[i,1]][:dossier][1,k] for i in 1 : 9],
    #                 [rms_MTiXsens[filenameMTiXsens[i,1]][:chest][1,k]/rms_MTiXsens[filenameMTiXsens[i,1]][:dossier][1,k] for i in 13 : 21]
    #             ),
    #         title = [string("MTiXsens RMS STHT chest-dossier along ",axis_name[k],"-axis")],
    #         layout = (1, 1),
    #         # legend = false,
    #         legend=:best,
    #         #ylims = (0,2),
    #         xlabel = "vel [-]",
    #         ylabel = "RMS STHT []",
    #         markercolor=linecolor,
    #         label=perturbLabel,
    #         markershape
    #     );
    #     pMTiXsens_RMS_STHT_chest_dossier=plot!([0.8 3.2]',[1 1]',
    #                                         linestyle=:dot,
    #                                         linecolor=:black,
    #                                         label="threshold")
    #
    #
    #     pMTiXsens_VDV_STHT_chest_dossier = scatter(
    #          plot_order,
    #         hcat([vdv_MTiXsens[filenameMTiXsens[i,1]][:chest][1,k]/vdv_MTiXsens[filenameMTiXsens[i,1]][:dossier][1,k] for i in 1 : 9],
    #                 [vdv_MTiXsens[filenameMTiXsens[i,1]][:chest][1,k]/vdv_MTiXsens[filenameMTiXsens[i,1]][:dossier][1,k] for i in 13 : 21],
    #             ),
    #         title = [string("MTiXsens VDV STHT chest-dossier along ",axis_name[k],"-axis")],
    #         layout = (1, 1),
    #         # legend = false,
    #         legend=:best,
    #         #ylims = (0,2),
    #         xlabel = "vel [-]",
    #         ylabel = "VDV STHT []",
    #         markercolor=linecolor,
    #         label=perturbLabel,
    #         markershape
    #     );
    #     pMTiXsens_VDV_STHT_chest_dossier=plot!([0.8 3.2]',[1 1]',
    #                                         linestyle=:dot,
    #                                         linecolor=:black,
    #                                         label="threshold")
    #
    #     savefig(pMTiXsens_SignalEnergy_STHT_chest_dossier,string("save_figure/figure_SE_RMS_VDV/SE_STHT_chest-dossier_",axis_name[k],"axis"))
    #     savefig(pMTiXsens_RMS_STHT_chest_dossier,string("save_figure/figure_SE_RMS_VDV/RMS_STHT_chest-dossier_",axis_name[k],"axis"))
    #     savefig(pMTiXsens_VDV_STHT_chest_dossier,string("save_figure/figure_SE_RMS_VDV/VDV_STHT_chest-dossier_",axis_name[k],"axis"))
    # end
end

function plotAnalysisPULSE_All(
            Energysignal,
            rms,
            vdv,
            crest_factor,
            filename,
            sensorName=[],
            perturbLabel=[],
            plot_order=[],
            perturbIndices=[],
            linecolor= [],
            perturbshape=[],
            addTitle="")

    #each sensor
    ####################################################################
    # plot_order=[1 1 1 2 2 2 3 3 3]',
    # linecolor= [:blue :red]
    # perturbshape=[:circle :utriangle]

    pSignalEnergy=plot(
    title = string(sensorName, " signal energy",addTitle),
    layout = (1, 1),
    # legend = false,
    legend=:best,
    #ylims = (0,30),
    xlabel = "vel [-]",
    ylabel = "SE []")

    addScatterPlot(pSignalEnergy,Energysignal,filename,plot_order,perturbIndices,perturbLabel,linecolor,perturbshape)

    ####################################################################
    pRMS=plot(
    title = string(sensorName, " RMS",addTitle),
    layout = (1, 1),
    # legend = false,
    legend=:best,
    #ylims = (0,30),
    xlabel = "vel [-]",
    ylabel = "RMS []")

    addScatterPlot(pRMS,rms,filename,plot_order,perturbIndices,perturbLabel,linecolor,perturbshape)

    ####################################################################
    pVDV=plot(
    title = string(sensorName, " VDV",addTitle),
    layout = (1, 1),
    # legend = false,
    legend=:best,
    #ylims = (0,30),
    xlabel = "vel [-]",
    ylabel = "VDV []")

    addScatterPlot(pVDV,vdv,filename,plot_order,perturbIndices,perturbLabel,linecolor,perturbshape)

    ####################################################################
    pcrest_factor=plot(
    title = string(sensorName, " Crest Factor",addTitle),
    layout = (1, 1),
    # legend = false,
    legend=:best,
    #ylims = (0,30),
    xlabel = "vel [-]",
    ylabel = "Crest Factor []")

    addScatterPlot(pcrest_factor,crest_factor,filename,plot_order,perturbIndices,perturbLabel,linecolor,perturbshape)

    pcrest_factor=plot!([0.8 3.2]',[9 9]',
                                        linestyle=:dot,
                                        linecolor=:black,
                                        label="threshold")

    savefig(pSignalEnergy,string("save_figure/figure_SE_RMS_VDV/SE_",sensorName,".png"))
    savefig(pRMS,string("save_figure/figure_SE_RMS_VDV/RMS_",sensorName,".png"))
    savefig(pVDV,string("save_figure/figure_SE_RMS_VDV/VDV_",sensorName,".png"))
    savefig(pcrest_factor,string("save_figure/figure_SE_RMS_VDV/crest_factor_",sensorName,".png"))
end

function plotAnalysisPULSE_All_2(
            std,
            min,
            max,
            range,
            filename,
            sensorName=[],
            axis_name=["x";"y"],
            perturbLabel=[],
            plot_order=[],
            perturbIndices=[],
            linecolor= [],
            perturbshape=[],
            addTitle="")

    #each sensor
    ####################################################################
    # plot_order=[1 1 1 2 2 2 3 3 3]',
    # linecolor= [:blue :red]
    # perturbshape=[:circle :utriangle]
    for k in 1:size(axis_name,1)
        pStd=Dict()
        pMTiXsens_RMS=Dict()
        pMTiXsens_VDV=Dict()
        pMTiXsens_crest_factor=Dict()

        pStd=plot(
        title = string(sensorName, " standard deviation of CoP"," along ",axis_name[k],"-axis",addTitle),
        layout = (1, 1),
        # legend = false,
        legend=:best,
        #ylims = (0,30),
        xlabel = "vel [-]",
        ylabel = "CoP [cm]")

        # addScatterPlot(pStd,std,filename,plot_order,perturbIndices,perturbLabel,linecolor,perturbshape)
        addScatterPlot_2(pStd,std,filename,plot_order,perturbIndices,k,perturbLabel,linecolor,perturbshape)
        ####################################################################
        pMin=plot(
        title = string(sensorName, " min of CoP"," along ",axis_name[k],"-axis",addTitle),
        layout = (1, 1),
        # legend = false,
        legend=:best,
        #ylims = (0,30),
        xlabel = "vel [-]",
        ylabel = "CoP [cm]")

        # addScatterPlot(pMin,min,filename,plot_order,perturbIndices,perturbLabel,linecolor,perturbshape)
        addScatterPlot_2(pMin,min,filename,plot_order,perturbIndices,k,perturbLabel,linecolor,perturbshape)

        ####################################################################
        pMax=plot(
        title = string(sensorName, " max of CoP"," along ",axis_name[k],"-axis",addTitle),
        layout = (1, 1),
        # legend = false,
        legend=:best,
        #ylims = (0,30),
        xlabel = "vel [-]",
        ylabel = "CoP [cm]")

        # addScatterPlot(pMax,max,filename,plot_order,perturbIndices,perturbLabel,linecolor,perturbshape)
        addScatterPlot_2(pMax,max,filename,plot_order,perturbIndices,k,perturbLabel,linecolor,perturbshape)

        ####################################################################
        pRange=plot(
        title = string(sensorName, " range of CoP"," along ",axis_name[k],"-axis",addTitle),
        layout = (1, 1),
        # legend = false,
        legend=:best,
        #ylims = (0,30),
        xlabel = "vel [-]",
        ylabel = "CoP [cm]")

        # addScatterPlot(pRange,range,filename,plot_order,perturbIndices,perturbLabel,linecolor,perturbshape)
        addScatterPlot_2(pRange,range,filename,plot_order,perturbIndices,k,perturbLabel,linecolor,perturbshape)

        ## #
        savefig(pStd,string("save_figure/figure_SE_RMS_VDV/STD_",sensorName,"_",axis_name[k],".png"))
        savefig(pMin,string("save_figure/figure_SE_RMS_VDV/MIN_",sensorName,"_",axis_name[k],".png"))
        savefig(pMax,string("save_figure/figure_SE_RMS_VDV/MAX_",sensorName,"_",axis_name[k],".png"))
        savefig(pRange,string("save_figure/figure_SE_RMS_VDV/RANGE_",sensorName,"_",axis_name[k],".png"))
    end

end

function splitFileName(filename,shiftFactor_perturb=0,shiftFactor_wheel=0)
    splitName=[]
    plot_order=[]
    perturbLabel=[]
    wheelType=[]
    for k in filename
            splitName_temp=split(k,'_')
            append!(splitName,[splitName_temp])
            append!(plot_order,parse(Int,splitName_temp[5][2]))
            append!(perturbLabel,[splitName_temp[4]])
            append!(wheelType,[splitName_temp[3]])
    end

    perturbLabel_unique=unique(perturbLabel)
    perturbIndices=[]
    for k in perturbLabel_unique
            append!(perturbIndices,findfirst(perturbLabel.==k))
    end
    for k in 1:size(perturbIndices,1)-1
            i1=perturbIndices[k]
            i2=perturbIndices[k+1]-1
            plot_order[i1:i2]=plot_order[i1:i2].+(shiftFactor_perturb*(k-1))
    end
    plot_order[perturbIndices[end]:end]=plot_order[perturbIndices[end]:end].+(shiftFactor_perturb*(size(perturbIndices,1)-1))

    wheelLabel_unique=unique(wheelType)
    wheelIndices=[]
    for k in wheelLabel_unique
            append!(wheelIndices,findfirst(wheelType.==k))
    end
    for k in 1:size(wheelIndices,1)-1
            i1=wheelIndices[k]
            i2=wheelIndices[k+1]-1
            plot_order[i1:i2]=plot_order[i1:i2].+(shiftFactor_wheel*(k-1))
    end
    plot_order[wheelIndices[end]:end]=plot_order[wheelIndices[end]:end].+(shiftFactor_wheel*(size(wheelIndices,1)-1))

    return splitName, plot_order, perturbLabel, perturbLabel_unique, perturbIndices, wheelType, wheelLabel_unique, wheelIndices
end


text:
#### Balance Law Interface

"""
    BalanceLaw

An abstract type representing a PDE balance law of the form

elements for balance laws of the form

```math
q_{,t} + Σ_{i=1,...d} F_{i,i} = s
```

Subtypes `L` should define the methods below
"""
abstract type BalanceLaw end # PDE part

"""
    vars_state_conservative(::L, FT)

a tuple of symbols containing the state variables
given a float type `FT`.
"""
function vars_state_conservative end

"""
    vars_state_auxiliary(::L, FT)

a tuple of symbols containing the auxiliary variables
given a float type `FT`.
"""
function vars_state_auxiliary end

"""
    vars_state_gradient(::L, FT)

a tuple of symbols containing the transformed variables
of which gradients are computed given a float type `FT`.
"""
function vars_state_gradient end

"""
    vars_state_gradient_flux(::L, FT)

a tuple of symbols containing the diffusive variables
given a float type `FT`.
"""
function vars_state_gradient_flux end

"""
    vars_gradient_laplacian(::L, FT)

a tuple of symbols containing the transformed variables
of which gradients of laplacian are computed, they must
be a subset of `vars_state_gradient`, given a float type `FT`.
"""
vars_gradient_laplacian(::BalanceLaw, FT) = @vars()

"""
    vars_hyperdiffusive(::L, FT)

a tuple of symbols containing the hyperdiffusive variables
given a float type `FT`.
"""
vars_hyperdiffusive(::BalanceLaw, FT) = @vars()

"""
    vars_integrals(::L, FT)

a tuple of symbols containing variables to be integrated
along a vertical stack, given a float type `FT`.
"""
vars_integrals(::BalanceLaw, FT) = @vars()

"""
    vars_reverse_integrals(::L, FT)

a tuple of symbols containing variables to be integrated
along a vertical stack, in reverse, given a float type `FT`.
"""
vars_reverse_integrals(::BalanceLaw, FT) = @vars()

"""
    init_state_conservative!(
      ::L,
      state_conservative::Vars,
      state_auxiliary::Vars,
      coords,
      args...)

Initialize the conservative state variables at ``t = 0``
"""
function init_state_conservative! end

"""
    init_state_auxiliary!(
      ::L,
      state_auxiliary::Vars,
      coords,
      args...)

Initialize the auxiliary state, at ``t = 0``
"""
function init_state_auxiliary! end

"""
    flux_first_order!(
        ::L,
        flux::Grad,
        state_conservative::Vars,
        state_auxiliary::Vars,
        t::Real
    )

Compute first-order flux terms in balance law equation
"""
function flux_first_order! end

"""
    flux_second_order!(
        ::L,
        flux::Grad,
        state_conservative::Vars,
        state_gradient_flux::Vars,
        hyperdiffusive::Vars,
        state_auxiliary::Vars,
        t::Real
    )

Compute second-order flux terms in balance law equation
"""
function flux_second_order! end

"""
    source!(
        ::L,
        source::Vars,
        state_conservative::Vars,
        diffusive::Vars,
        state_auxiliary::Vars,
        t::Real
    )

Compute non-conservative source terms in balance law equation
"""
function source! end

"""
    compute_gradient_argument!(
        ::L,
        transformstate::Vars,
        state_conservative::Vars,
        state_auxiliary::Vars,
        t::Real
    )

transformation of state variables to variables of which gradients are computed
"""
compute_gradient_argument!(::BalanceLaw, args...) = nothing

"""
    compute_gradient_flux!(
        ::L,
        state_gradient_flux::Vars,
        ∇transformstate::Grad,
        state_auxiliary::Vars,
        t::Real
    )

transformation of gradients to the diffusive variables
"""
compute_gradient_flux!(::BalanceLaw, args...) = nothing

"""
    transform_post_gradient_laplacian!(
        ::L,
        Qhypervisc_div::Vars,
        ∇Δtransformstate::Grad,
        state_auxiliary::Vars,
        t::Real
    )

transformation of laplacian gradients to the hyperdiffusive variables
"""
function transform_post_gradient_laplacian! end

"""
    wavespeed(
        ::L,
        n⁻,
        state_conservative::Vars,
        state_auxiliary::Vars,
        t::Real
    )

wavespeed
"""
function wavespeed end

"""
    boundary_state!(
        ::NumericalFluxGradient,
        ::L,
        state_conservative⁺::Vars,
        state_auxiliary⁺::Vars,
        normal⁻,
        state_conservative⁻::Vars,
        state_auxiliary⁻::Vars,
        bctype,
        t
    )
    boundary_state!(
        ::NumericalFluxFirstOrder,
        ::L,
        state_conservative⁺::Vars,
        state_auxiliary⁺::Vars,
        normal⁻,
        state_conservative⁻::Vars,
        state_auxiliary⁻::Vars,
        bctype,
        t
    )
    boundary_state!(
        ::NumericalFluxSecondOrder,
        ::L,
        state_conservative⁺::Vars,
        state_gradient_flux⁺::Vars,
        state_auxiliary⁺:
        Vars, normal⁻,
        state_conservative⁻::Vars,
        state_gradient_flux⁻::Vars,
        state_auxiliary⁻::Vars,
        bctype,
        t
    )

Apply boundary conditions for

 - `NumericalFluxGradient` numerical flux (internal method)
 - `NumericalFluxFirstOrder` first-order unknowns
 - `NumericalFluxSecondOrder` second-order unknowns

"""
function boundary_state! end

"""
    update_auxiliary_state!(
        dg::DGModel,
        m::HBModel,
        Q::MPIStateArray,
        t::Real,
        elems::UnitRange,
    )

Update the auxiliary state variables
"""
function update_auxiliary_state! end

"""
    update_auxiliary_state_gradient!

Update the auxiliary state gradient variables
"""
function update_auxiliary_state_gradient! end

"""
    integral_load_auxiliary_state!

Specify how to compute integrands. Can be functions of the conservative state and auxiliary variables.
"""
function integral_load_auxiliary_state! end

"""
    integral_set_auxiliary_state!

Specify which auxiliary variables are used to store the output of the integrals.
"""
function integral_set_auxiliary_state! end

"""
    reverse_integral_load_auxiliary_state!

Specify auxiliary variables need their integrals reversed.
"""
function reverse_integral_load_auxiliary_state! end

"""
    reverse_integral_set_auxiliary_state!

Specify which auxiliary variables are used to store the output of the reversed integrals.
"""
function reverse_integral_set_auxiliary_state! end

# Internal methods
number_state_conservative(m::BalanceLaw, FT) =
    varsize(vars_state_conservative(m, FT))
number_state_auxiliary(m::BalanceLaw, FT) = varsize(vars_state_auxiliary(m, FT))
number_state_gradient(m::BalanceLaw, FT) = varsize(vars_state_gradient(m, FT))
number_state_gradient_flux(m::BalanceLaw, FT) =
    varsize(vars_state_gradient_flux(m, FT))
num_gradient_laplacian(m::BalanceLaw, FT) =
    varsize(vars_gradient_laplacian(m, FT))
num_hyperdiffusive(m::BalanceLaw, FT) = varsize(vars_hyperdiffusive(m, FT))
num_integrals(m::BalanceLaw, FT) = varsize(vars_integrals(m, FT))
num_reverse_integrals(m::BalanceLaw, FT) =
    varsize(vars_reverse_integrals(m, FT))

ext: jl


text:
using Electron

app = Application();
w = Window(app);

# OLD 1

# function handle(message)
#     @info message["cmd"]
# end

# Electron.set_msg_handler(handle,w)

# run(w, "const {ipcRenderer} = require('electron')") //done for me already
# run(w, "ipcRenderer.send('msg-for-julia-process', {cmd: 'message for Julia'})") //done for me already

# OLD 2

# @async begin
#     while true
#         m = read(w)
#         println(m)
#     end
# end

# CURRENT
@async for msg in msgchannel(w)
    @info msg
end

run(w, "sendMessageToJulia({cmd: 'message for Julia'})")

read(w) #to get the message


text:
function contact_forces(u, p, t)
    @unpack footang, g, iff, ihat, irf, ish, ith, k1, k2, k3, k4, k5, k6, k7, k8, lff, lffo, lhat, lmtpxi, lrf, lrff, lrffo, lrfo, lsh, lsho, lth, ltho, ltoexi, mff, mhat, mrf, msh, mth, mtpb, mtpk, rmtpxi, rtoexi, u4, u5, u6, u7, u8, u9, z = p
    @inbounds q1, q2, q3, u1, u2, u3 = u

    # specified variables
    rmtp = _rmtp(t)
    ra = _ra(t)
    rk = _rk(t)
    rh = _rh(t)
    la = _la(t)
    lmtp = _lmtp(t)
    lh = _lh(t)
    lk = _lk(t)
    lkp = _lkp(t)
    lap = _lap(t)
    lhp = _lhp(t)
    lmtpp = _lmtpp(t)
    rhp = _rhp(t)
    rkp = _rkp(t)
    rap = _rap(t)
    rmtpp = _rmtpp(t)

    # grf
    pop2y = q2 - lff * sin(la + lh + lmtp - lk - q3)
    pop10y = q2 + lth * sin(lh - q3) + lsh * sin(rh - rk - q3) + lrf * sin(la + lh - lk - q3) - lth * sin(rh - q3) - lsh * sin(lh - lk - q3) - lrf * sin(ra + rh - rk - q3) - lff * sin(la + lh + lmtp - lk - q3)
    pop11y = q2 + lth * sin(lh - q3) + lsh * sin(rh - rk - q3) + lrf * sin(la + lh - lk - q3) + lff * sin(ra + rh + rmtp - rk - q3) - lth * sin(rh - q3) - lsh * sin(lh - lk - q3) - lrf * sin(ra + rh - rk - q3) - lff * sin(la + lh + lmtp - lk - q3)

    # left foot
    if q2 < 0.0
        dp1x = q1 - ltoexi
        lry1 = -k3 * q2 - k4 * abs(q2) * u2
        lrx1 = lry1 * (-k1 * dp1x - k2 * u1) # (-k1*(q1-pop1xi) - k2*u1)*ry1
    else
        lrx1 = 0.0
        lry1 = 0.0
    end

    if pop2y < 0.0
        pop2x = q1 + lff * cos(la + lh + lmtp - lk - q3)
        vop2x = u1 - lff * sin(la + lh + lmtp - lk - q3) * (lap + lhp + lmtpp - lkp - u3 - u5 - u7 - u9)
        vop2y = u2 - lff * cos(la + lh + lmtp - lk - q3) * (lap + lhp + lmtpp - lkp - u3 - u5 - u7 - u9)
        dp2x = pop2x - lmtpxi
        lry2 = -k7 * pop2y - k8 * abs(pop2y) * vop2y
        lrx2 = lry2 * (-k5 * dp2x - k6 * vop2x)
    else
        lrx2 = 0.0
        lry2 = 0.0
    end

    # right foot
    if pop11y < 0.0
        pop11x = q1 + lth * cos(rh - q3) + lsh * cos(lh - lk - q3) + lrf * cos(ra + rh - rk - q3) + lff * cos(la + lh + lmtp - lk - q3) - lth * cos(lh - q3) - lsh * cos(rh - rk - q3) - lrf * cos(la + lh - lk - q3) - lff * cos(ra + rh + rmtp - rk - q3)
        vop11x = u1 + lth * sin(lh - q3) * (lhp - u3 - u5) + lsh * sin(rh - rk - q3) * (rhp - rkp - u3 - u4 - u6) + lrf * sin(la + lh - lk - q3) * (lap + lhp + lmtpp - lkp - u3 - u5 - u7 - u9) + lff * sin(ra + rh + rmtp - rk - q3) * (rap + rhp + rmtpp - rkp - u3 - u4 - u6 - u8) - lth * sin(rh - q3) * (rhp - u3 - u4) - lsh * sin(lh - lk - q3) * (lhp - lkp - u3 - u5 - u7) - lrf * sin(ra + rh - rk - q3) * (rap + rhp - rkp - u3 - u4 - u6 - u8) - lff * sin(la + lh + lmtp - lk - q3) * (lap + lhp + lmtpp - lkp - u3 - u5 - u7 - u9)
        vop11y = u2 + lth * cos(lh - q3) * (lhp - u3 - u5) + lsh * cos(rh - rk - q3) * (rhp - rkp - u3 - u4 - u6) + lrf * cos(la + lh - lk - q3) * (lap + lhp + lmtpp - lkp - u3 - u5 - u7 - u9) + lff * cos(ra + rh + rmtp - rk - q3) * (rap + rhp + rmtpp - rkp - u3 - u4 - u6 - u8) - lth * cos(rh - q3) * (rhp - u3 - u4) - lsh * cos(lh - lk - q3) * (lhp - lkp - u3 - u5 - u7) - lrf * cos(ra + rh - rk - q3) * (rap + rhp - rkp - u3 - u4 - u6 - u8) - lff * cos(la + lh + lmtp - lk - q3) * (lap + lhp + lmtpp - lkp - u3 - u5 - u7 - u9)
        dp11x = pop11x - rtoexi
        rry1 = -k3 * pop11y - k4 * abs(pop11y) * vop11y
        rrx1 = rry1 * (-k1 * dp11x - k2 * vop11x) # (-k1*(q1-pop1xi) - k2*u1)*ry1
    else
        rrx1 = 0.0
        rry1 = 0.0
    end

    if pop10y < 0.0
        pop10x = q1 + lth * cos(rh - q3) + lsh * cos(lh - lk - q3) + lrf * cos(ra + rh - rk - q3) + lff * cos(la + lh + lmtp - lk - q3) - lth * cos(lh - q3) - lsh * cos(rh - rk - q3) - lrf * cos(la + lh - lk - q3)
        vop10x = u1 + lth * sin(lh - q3) * (lhp - u3 - u5) + lsh * sin(rh - rk - q3) * (rhp - rkp - u3 - u4 - u6) + lrf * sin(la + lh - lk - q3) * (lap + lhp + lmtpp - lkp - u3 - u5 - u7 - u9) - lth * sin(rh - q3) * (rhp - u3 - u4) - lsh * sin(lh - lk - q3) * (lhp - lkp - u3 - u5 - u7) - lrf * sin(ra + rh - rk - q3) * (rap + rhp - rkp - u3 - u4 - u6 - u8) - lff * sin(la + lh + lmtp - lk - q3) * (lap + lhp + lmtpp - lkp - u3 - u5 - u7 - u9)
        vop10y = u2 + lth * cos(lh - q3) * (lhp - u3 - u5) + lsh * cos(rh - rk - q3) * (rhp - rkp - u3 - u4 - u6) + lrf * cos(la + lh - lk - q3) * (lap + lhp + lmtpp - lkp - u3 - u5 - u7 - u9) - lth * cos(rh - q3) * (rhp - u3 - u4) - lsh * cos(lh - lk - q3) * (lhp - lkp - u3 - u5 - u7) - lrf * cos(ra + rh - rk - q3) * (rap + rhp - rkp - u3 - u4 - u6 - u8) - lff * cos(la + lh + lmtp - lk - q3) * (lap + lhp + lmtpp - lkp - u3 - u5 - u7 - u9)
        dp10x = pop10x - rmtpxi
        rry2 = -k7 * pop10y - k8 * abs(pop10y) * vop10y
        rrx2 = rry2 * (-k5 * dp10x - k6 * vop10x)
    else
        rrx2 = 0.0
        rry2 = 0.0
    end

    return lrx1, lry1, lrx2, lry2, rrx1, rry1, rrx2, rry2
end


text:
using CSV
using DataFrames
using LinearAlgebra
using RecipesBase
using StateSpaceModels
using Statistics
using Test

# Functions that are used in different tests
include("utils.jl")

# Core functionality
include("systems.jl")
include("prints.jl")

# Models
include("models/locallevel.jl")
include("models/locallineartrend.jl")
include("models/locallevelcycle.jl")
include("models/locallevelexplanatory.jl")
include("models/basicstructural.jl")
include("models/basicstructural_explanatory.jl")
include("models/basicstructural_multivariate.jl")
include("models/sarima.jl")
include("models/unobserved_components.jl")
include("models/linear_regression.jl")
include("models/exponential_smoothing.jl")
include("models/dar.jl")

# Visualization
include("visualization/forecast.jl")
include("visualization/components.jl")
include("visualization/backtest.jl")
include("visualization/diagnostics.jl")


text:
# Preamble ########################################################################################
using Pkg
Pkg.activate(Base.current_project())

cd("..") # Go to the parent directory
###################################################################################################
# Packages needed #################################################################################
using UnitRots
using DelimitedFiles
using ImageCore
using ImageIO
using ImageMagick
using FileIO
using PyPlot
using BenchmarkTools
###################################################################################################
# Angle of rotation ###############################################################################

const θ = π/4

###################################################################################################
# Retrieve of data from image #####################################################################

@info "Retrieving the image and extracting the two-dimensional field"
input = Int64.(channelview(Gray.(load("data/benchmark/delta31.png"))))

###################################################################################################
# Identification of the grayscale field dimensions ################################################

j = (size(input)[1] - 1)/2;

N = Integer(2*j + 1);

@info "The size of the field is $N"

###################################################################################################
# Rotation of the image ###########################################################################

@info "Performing rotation at θ = $θ"

output = real(Emodes(input, θ))

###################################################################################################
begin
	fig, (ax1, ax2) = subplots(1, 2, figsize = (20,10), constrained_layout = true)

	ax1.imshow(real(input), cmap = "gray")
    ax1.set_xticks([-0.5:1:N;])
    ax1.set_yticks([-0.5:1:N;])
    ax1.set_xticklabels([])
    ax1.set_yticklabels([])
    ax1.tick_params(axis = "both", length = 0)
    ax1.grid(which = "major", color = "gray", linestyle = "-", linewidth = 1.5, alpha = 0.5)

    ax2.imshow(nnorm(output), cmap = "gray")
    ax2.set_xticks([-0.5:1:N;])
    ax2.set_yticks([-0.5:1:N;])
    ax2.set_xticklabels([])
    ax2.set_yticklabels([])
    ax2.tick_params(axis = "both", length = 0)
    ax2.grid(which = "major", color = "gray", linestyle = "-", linewidth = 1.5, alpha = 0.5)

    savefig("figures/benchmark.png", dpi = 300)
end

max_stars_repo_name: rurz/UnitRots

max_issues_repo_name: rurz/UnitRots

max_forks_repo_name: rurz/UnitRots


text:
# Load the packages we need
using JLD2, Plots, DataFrames

# Load the parameters data frame
@load "params.jld2" params

# Create empty arrays that will hold the R2 values
R2_train = []
R2_val = []

# Let's load all the results we want, and put them into the arrays
for i = 1:nrow(params)
	experiment_id = params[i, :experiment_id]
	@load "results/results-$experiment_id.jld2" R2tr R2va
	push!(R2_train, R2tr)
	push!(R2_val, R2va)
end

# Add the R2 values to the dataframe
params[!, :R2_train] = R2_train
params[!, :R2_val] = R2_val

# Subset the dataframe to only some of the results
dataforplot = filter(row -> row[:n_estimators] == 160 && row[:max_features] == 2, params)

# Make the plot
plot(dataforplot[!, :max_depth], [dataforplot[!, :R2_train] dataforplot[!, :R2_val]],
        label = ["Training" "Validation"], xlabel="Maximum tree depth",
        ylabel = "R2")
Plots.savefig("plot.pdf")


text:
struct SSB <: abstractBarycenter end

naif_id(::Type{SSB}) = 0


text:
__precompile__(true)
module Reversi

export
    # types
    Board, State, Player, RandomPlayer, Minimax, StochasticMinimax,

    # board
    encode, legal_actions_bits, next_state_bits, score, decode, legal_actions,
    game_over,

    # state
    do_action, do_pass, opponent,

    # player
    select_action

const MARKERS = Dict(0 => "   ", 2 => " \u25cf ", 1 => " \u25cb ")

include("board.jl")
include("state.jl")
include("players.jl")
include("gameplay.jl")

include("edax.jl")
# package code goes here

end # module


text:
# calculates the mean from N stochastic A->B reactions at different rates
# this really tests different ways of constructing the jump problems
using DiffEqBase, DiffEqJump, Statistics
using Test

# using BenchmarkTools
# dobenchmark = true

doprint     = false
dotest      = true
Nrxs        = 16
Nsims       = 8000
tf          = .1
baserate    = .1
A0          = 100
exactmean   = (t,ratevec) -> A0 * exp(-sum(ratevec) * t)
SSAalgs     = [Direct(),RSSA()]

spec_to_dep_jumps = [collect(1:Nrxs),[]]
jump_to_dep_specs = [[1,2] for i=1:Nrxs]
namedpars = (vartojumps_map=spec_to_dep_jumps, jumptovars_map=jump_to_dep_specs)
rates = ones(Float64, Nrxs) * baserate;
cumsum!(rates, rates)
exactmeanval = exactmean(tf, rates)

function runSSAs(jump_prob)
    Asamp = zeros(Int,Nsims)
    for i in 1:Nsims
        sol = solve(jump_prob, SSAStepper())
        Asamp[i] = sol[1,end]
    end
    mean(Asamp)
end

# uses constant jumps as a tuple within a JumpSet
function A_to_B_tuple(N, method)
    # jump reactions
    jumpvec = []
    for i in 1:N
        ratefunc = (u,p,t) -> rates[i] * u[1]
        affect!  = function (integrator)
            integrator.u[1] -= 1
            integrator.u[2] += 1
        end
        push!(jumpvec, ConstantRateJump(ratefunc, affect!))
    end

    # convert jumpvec to tuple to send to JumpProblem...
    jumps     = ((jump for jump in jumpvec)...,)
    jset      = JumpSet((), jumps, nothing, nothing)
    prob      = DiscreteProblem([A0,0], (0.0,tf))
    jump_prob = JumpProblem(prob, method, jset; save_positions=(false,false), namedpars...)

    jump_prob
end

# uses constant jumps as a vector within a JumpSet
function A_to_B_vec(N, method)
    # jump reactions
    jumps = Vector{ConstantRateJump}()
    for i in 1:N
        ratefunc = (u,p,t) -> rates[i] * u[1]
        affect!  = function (integrator)
            integrator.u[1] -= 1
            integrator.u[2] += 1
        end
        push!(jumps, ConstantRateJump(ratefunc, affect!))
    end

    # convert jumpvec to tuple to send to JumpProblem...
    jset      = JumpSet((), jumps, nothing, nothing)
    prob      = DiscreteProblem([A0,0], (0.0,tf))
    jump_prob = JumpProblem(prob, method, jset; save_positions=(false,false), namedpars...)

    jump_prob
end

# uses a single mass action jump to represent all reactions
function A_to_B_ma(N, method)
    reactstoch = Vector{Vector{Pair{Int,Int}}}();
    netstoch   = Vector{Vector{Pair{Int,Int}}}();
    for i = 1:N
        push!(reactstoch,[1 => 1])
        push!(netstoch,[1 => -1, 2=>1])
    end

    majumps   = MassActionJump(rates, reactstoch, netstoch)
    jset      = JumpSet((), (), nothing, majumps)
    prob      = DiscreteProblem([A0,0], (0.0,tf))
    jump_prob = JumpProblem(prob, method, jset; save_positions=(false,false), namedpars...)

    jump_prob
end

# uses one mass action jump to represent half the reactions and a vector
# of constant jumps for the other half. Stores them in a JumpSet
function A_to_B_hybrid(N, method)
    # half reactions are treated as mass action and half as constant jumps
    switchidx = (N//2).num

    # mass action reactions
    reactstoch = Vector{Vector{Pair{Int,Int}}}();
    netstoch   = Vector{Vector{Pair{Int,Int}}}();
    for i in 1:switchidx
        push!(reactstoch,[1 => 1])
        push!(netstoch,[1 => -1, 2=>1])
    end

     # jump reactions
     jumps = Vector{ConstantRateJump}()
     for i in (switchidx+1):N
         ratefunc = (u,p,t) -> rates[i] * u[1]
         affect!  = function (integrator)
             integrator.u[1] -= 1
             integrator.u[2] += 1
         end
         push!(jumps, ConstantRateJump(ratefunc, affect!))
     end

    majumps   = MassActionJump(rates[1:switchidx] , reactstoch, netstoch)
    jset      = JumpSet((), jumps, nothing, majumps)
    prob      = DiscreteProblem([A0,0], (0.0,tf))
    jump_prob = JumpProblem(prob, method, jset; save_positions=(false,false), namedpars...)

    jump_prob
end

# uses a mass action jump to represent half the reactions and a vector
# of constant jumps for the other half. Passes them to JumpProblem as a splatted tuple
function A_to_B_hybrid_nojset(N, method)
    # half reactions are treated as mass action and half as constant jumps
    switchidx = (N//2).num

    # mass action reactions
    reactstoch = Vector{Vector{Pair{Int,Int}}}();
    netstoch   = Vector{Vector{Pair{Int,Int}}}();
    for i in 1:switchidx
        push!(reactstoch,[1 => 1])
        push!(netstoch,[1 => -1, 2=>1])
    end

     # jump reactions
     jumpvec = Vector{ConstantRateJump}()
     for i in (switchidx+1):N
         ratefunc = (u,p,t) -> rates[i] * u[1]
         affect!  = function (integrator)
             integrator.u[1] -= 1
             integrator.u[2] += 1
         end
         push!(jumpvec, ConstantRateJump(ratefunc, affect!))
     end
    constjumps = (jump for jump in jumpvec)
    majumps   = MassActionJump(rates[1:switchidx] , reactstoch, netstoch)
    jumps     = (constjumps...,majumps)
    prob      = DiscreteProblem([A0,0], (0.0,tf))
    jump_prob = JumpProblem(prob, method, jumps...; save_positions=(false,false), namedpars...)

    jump_prob
end

# uses a vector of mass action jumps of vectors to represent half the reactions and a vector
# of constant jumps for the other half. Passes them to JumpProblem as a JumpSet
function A_to_B_hybrid_vecs(N, method)
    # half reactions are treated as mass action and half as constant jumps
    switchidx = (N//2).num

    # mass action reactions
    majumps = Vector{MassActionJump}()
    for i in 1:switchidx
        push!(majumps, MassActionJump([rates[i]], [[1 => 1]], [[1 => -1, 2=>1]] ))
    end

     # jump reactions
     jumpvec = Vector{ConstantRateJump}()
     for i in (switchidx+1):N
         ratefunc = (u,p,t) -> rates[i] * u[1]
         affect!  = function (integrator)
             integrator.u[1] -= 1
             integrator.u[2] += 1
         end
         push!(jumpvec, ConstantRateJump(ratefunc, affect!))
     end
    jset      = JumpSet((), jumpvec, nothing, majumps)
    prob      = DiscreteProblem([A0,0], (0.0,tf))
    jump_prob = JumpProblem(prob, method, jset; save_positions=(false,false), namedpars...)

    jump_prob
end

# uses a vector of scalar mass action jumps to represent half the reactions and a vector
# of constant jumps for the other half. Passes them to JumpProblem as a JumpSet
function A_to_B_hybrid_vecs_scalars(N, method)
    # half reactions are treated as mass action and half as constant jumps
    switchidx = (N//2).num

    # mass action reactions
    majumps = Vector{MassActionJump}()
    for i in 1:switchidx
        push!(majumps, MassActionJump(rates[i], [1 => 1], [1 => -1, 2=>1] ))
    end

     # jump reactions
     jumpvec = Vector{ConstantRateJump}()
     for i in (switchidx+1):N
         ratefunc = (u,p,t) -> rates[i] * u[1]
         affect!  = function (integrator)
             integrator.u[1] -= 1
             integrator.u[2] += 1
         end
         push!(jumpvec, ConstantRateJump(ratefunc, affect!))
     end
    jset      = JumpSet((), jumpvec, nothing, majumps)
    prob      = DiscreteProblem([A0,0], (0.0,tf))
    jump_prob = JumpProblem(prob, method, jset; save_positions=(false,false), namedpars...)

    jump_prob
end

# uses a vector of scalar mass action jumps to represent half the reactions and a vector
# of constant jumps for the other half. Passes them to JumpProblem as a single splatted tuple.
function A_to_B_hybrid_tups_scalars(N, method)
    # half reactions are treated as mass action and half as constant jumps
    switchidx = (N//2).num

    # mass action reactions
    majumpsv = Vector{MassActionJump}()
    for i in 1:switchidx
        push!(majumpsv, MassActionJump(rates[i], [1 => 1], [1 => -1, 2=>1] ))
    end

    # jump reactions
    jumpvec = Vector{ConstantRateJump}()
    for i in (switchidx+1):N
        ratefunc = (u,p,t) -> rates[i] * u[1]
        affect!  = function (integrator)
            integrator.u[1] -= 1
            integrator.u[2] += 1
        end
        push!(jumpvec, ConstantRateJump(ratefunc, affect!))
    end

    jumps     = ((maj for maj in majumpsv)..., (jump for jump in jumpvec)...)
    prob      = DiscreteProblem([A0,0], (0.0,tf))
    jump_prob = JumpProblem(prob, method, jumps...; save_positions=(false,false), namedpars...)

    jump_prob
end

# uses a mass action jump to represent half the reactions and a tuple
# of constant jumps for the other half. Passes them to JumpProblem as a JumpSet.
function A_to_B_hybrid_tups(N, method)
    # half reactions are treated as mass action and half as constant jumps
    switchidx = (N//2).num

    # mass action reactions
    majumps = Vector{MassActionJump}()
    for i in 1:switchidx
        push!(majumps, MassActionJump([rates[i]], [[1 => 1]], [[1 => -1, 2=>1]] ))
    end

     # jump reactions
     jumpvec = Vector{ConstantRateJump}()
     for i in (switchidx+1):N
        ratefunc = (u,p,t) -> rates[i] * u[1]
        affect!  = function (integrator)
            integrator.u[1] -= 1
            integrator.u[2] += 1
        end
        push!(jumpvec, ConstantRateJump(ratefunc, affect!))
     end
    jumps    = ((jump for jump in jumpvec)...,)
    jset      = JumpSet((), jumps, nothing, majumps)
    prob      = DiscreteProblem([A0,0], (0.0,tf))
    jump_prob = JumpProblem(prob, method, jset; save_positions=(false,false), namedpars...)

    jump_prob
end

 jump_prob_gens = [A_to_B_tuple, A_to_B_vec, A_to_B_ma, A_to_B_hybrid, A_to_B_hybrid_nojset,
                   A_to_B_hybrid_vecs, A_to_B_hybrid_vecs_scalars, A_to_B_hybrid_tups,A_to_B_hybrid_tups_scalars]
#jump_prob_gens = [A_to_B_tuple, A_to_B_ma, A_to_B_hybrid, A_to_B_hybrid_vecs, A_to_B_hybrid_vecs_scalars,A_to_B_hybrid_tups_scalars]

for method in SSAalgs
    for jump_prob_gen in jump_prob_gens
        jump_prob = jump_prob_gen(Nrxs, method)
        meanval   = runSSAs(jump_prob)
        if doprint
            println("Method: ", method, ", Jump input types: ", jump_prob_gen,
                    ", sample mean = ", meanval, ", actual mean = ", exactmeanval)
        end
        @test abs(meanval - exactmeanval) < 1.

        # if dobenchmark
        #     @btime (runSSAs($jump_prob);)
        # end
    end
end

# for dependency graph methods just test with mass action jumps
SSAalgs        = [NRM(), SortingDirect(), DirectCR()]
jump_prob_gens = [A_to_B_ma]
for method in SSAalgs
    for jump_prob_gen in jump_prob_gens
        jump_prob = jump_prob_gen(Nrxs, method)
        meanval   = runSSAs(jump_prob)
        if doprint
            println("Method: ", method, ", Jump input types: ", jump_prob_gen,
                    ", sample mean = ", meanval, ", actual mean = ", exactmeanval)
        end
        @test abs(meanval - exactmeanval) < 1.

        # if dobenchmark
        #     @btime (runSSAs($jump_prob);)
        # end
    end
end


text:
@doc raw"""
    CompressibleEulerMulticomponentEquations2D(; gammas, gas_constants)

!!! warning "Experimental code"
    This system of equations is experimental and can change any time.

Multicomponent version of the compressible Euler equations
```math
\partial t
\begin{pmatrix}
\rho v_1 \\ \rho v_2 \\ E \\ \rho_1 \\ \rho_2 \\ \vdots \\ \rho_{n}
\end{pmatrix}
+
\partial x
\begin{pmatrix}
\rho v_1^2 + p \\ \rho v_1 v_2 \\ (E+p) v_1 \\ \rho_1 v_1 \\ \rho_2 v_1 \\ \vdots \\ \rho_{n} v_1
\end{pmatrix}
+
\partial y
\begin{pmatrix}
\rho v_1 v_2 \\ \rho v_2^2 + p \\ (E+p) v_2 \\ \rho_1 v_2 \\ \rho_2 v_2 \\ \vdots \\ \rho_{n} v_2
\end{pmatrix}
=
\begin{pmatrix}
0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \vdots \\ 0
\end{pmatrix}
```
for calorically perfect gas in two space dimensions.

In case of more than one component, the specific heat ratios `gammas` and the gas constants
`gas_constants` in [kJ/(kg*K)] should be passed as tuples, e.g., `gammas=(1.4, 1.667)`.

The remaining variables like the specific heats at constant volume 'cv' or the specific heats at
constant pressure 'cp' are then calculated considering a calorically perfect gas.
"""
struct CompressibleEulerMulticomponentEquations2D{NVARS, NCOMP, RealT<:Real} <: AbstractCompressibleEulerMulticomponentEquations{2, NVARS, NCOMP}
  gammas            ::SVector{NCOMP, RealT}
  gas_constants     ::SVector{NCOMP, RealT}
  cv                ::SVector{NCOMP, RealT}
  cp                ::SVector{NCOMP, RealT}

  function CompressibleEulerMulticomponentEquations2D{NVARS, NCOMP, RealT}(gammas       ::SVector{NCOMP, RealT},
                                                                           gas_constants::SVector{NCOMP, RealT}) where {NVARS, NCOMP, RealT<:Real}

    NCOMP >= 1 || throw(DimensionMismatch("`gammas` and `gas_constants` have to be filled with at least one value"))

    cv = gas_constants ./ (gammas .- 1)
    cp = gas_constants + gas_constants ./ (gammas .- 1)

    new(gammas, gas_constants,cv, cp)
  end
end

function CompressibleEulerMulticomponentEquations2D(; gammas, gas_constants)

  _gammas        = promote(gammas...)
  _gas_constants = promote(gas_constants...)
  RealT          = promote_type(eltype(_gammas), eltype(_gas_constants))

  NVARS = length(_gammas) + 3
  NCOMP = length(_gammas)

  __gammas        = SVector(map(RealT, _gammas))
  __gas_constants = SVector(map(RealT, _gas_constants))

  return CompressibleEulerMulticomponentEquations2D{NVARS, NCOMP, RealT}(__gammas, __gas_constants)
end

@inline Base.real(::CompressibleEulerMulticomponentEquations2D{NVARS, NCOMP, RealT}) where {NVARS, NCOMP, RealT} = RealT

function varnames(::typeof(cons2cons), equations::CompressibleEulerMulticomponentEquations2D)

  cons  = ("rho_v1", "rho_v2", "rho_e")
  rhos  = ntuple(n -> "rho" * string(n), Val(ncomponents(equations)))
  return (cons..., rhos...)
end

function varnames(::typeof(cons2prim), equations::CompressibleEulerMulticomponentEquations2D)

  prim  = ("v1", "v2", "p")
  rhos  = ntuple(n -> "rho" * string(n), Val(ncomponents(equations)))
  return (prim..., rhos...)
end

# Set initial conditions at physical location `x` for time `t`

"""
    initial_condition_convergence_test(x, t, equations::CompressibleEulerMulticomponentEquations2D)

A smooth initial condition used for convergence tests in combination with
[`source_terms_convergence_test`](@ref)
(and [`boundary_condition_convergence_test`](@ref) in non-periodic domains).
"""
function initial_condition_convergence_test(x, t, equations::CompressibleEulerMulticomponentEquations2D)
  c       = 2
  A       = 0.1
  L       = 2
  f       = 1/L
  omega   = 2 * pi * f
  ini     = c + A * sin(omega * (x[1] + x[2] - t))

  v1      = 1.0
  v2      = 1.0

  rho     = ini

  # Here we compute an arbitrary number of different rhos. (one rho is double the next rho while the sum of all rhos is 1)
  prim_rho  = SVector{ncomponents(equations), real(equations)}(2^(i-1) * (1-2)/(1-2^ncomponents(equations)) * rho for i in eachcomponent(equations))

  prim1 = rho * v1
  prim2 = rho * v2
  prim3 = rho^2

  prim_other = SVector{3, real(equations)}(prim1, prim2, prim3)

  return vcat(prim_other, prim_rho)
end

"""
    source_terms_convergence_test(u, x, t, equations::CompressibleEulerMulticomponentEquations2D)

Source terms used for convergence tests in combination with
[`initial_condition_convergence_test`](@ref)
(and [`boundary_condition_convergence_test`](@ref) in non-periodic domains).
"""
@inline function source_terms_convergence_test(u, x, t, equations::CompressibleEulerMulticomponentEquations2D)
  # Same settings as in `initial_condition`
  c       = 2
  A       = 0.1
  L       = 2
  f       = 1/L
  omega   = 2 * pi * f

  gamma  = totalgamma(u, equations)

  x1, x2  = x
  si, co  = sincos((x1 + x2 - t)*omega)
  tmp1    = co * A * omega
  tmp2    = si * A
  tmp3    = gamma - 1
  tmp4    = (2*c - 1)*tmp3
  tmp5    = (2*tmp2*gamma - 2*tmp2 + tmp4 + 1)*tmp1
  tmp6    = tmp2 + c

  # Here we compute an arbitrary number of different rhos. (one rho is double the next rho while the sum of all rhos is 1
  du_rho  = SVector{ncomponents(equations), real(equations)}(2^(i-1) * (1-2)/(1-2^ncomponents(equations)) * tmp1 for i in eachcomponent(equations))

  du1 = tmp5
  du2 = tmp5
  du3 = 2*((tmp6 - 1.0)*tmp3 + tmp6*gamma)*tmp1

  du_other  = SVector{3, real(equations)}(du1, du2, du3)

  return vcat(du_other, du_rho)
end

"""
    initial_condition_shock_bubble(x, t, equations::CompressibleEulerMulticomponentEquations2D{5, 2})

A shock-bubble testcase for multicomponent Euler equations
- Ayoub Gouasmi, Karthik Duraisamy, Scott Murman
  Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations
  [arXiv: 1904.00972](https://arxiv.org/abs/1904.00972)
"""
function initial_condition_shock_bubble(x, t, equations::CompressibleEulerMulticomponentEquations2D{5, 2})
  # bubble test case, see Gouasmi et al. https://arxiv.org/pdf/1904.00972
  # other reference: https://www.researchgate.net/profile/Pep_Mulet/publication/222675930_A_flux-split_algorithm_applied_to_conservative_models_for_multicomponent_compressible_flows/links/568da54508aeaa1481ae7af0.pdf
  # typical domain is rectangular, we change it to a square, as Trixi can only do squares
  @unpack gas_constants = equations

  # Positivity Preserving Parameter, can be set to zero if scheme is positivity preserving
  delta   = 0.03

  # Region I
  rho1_1  = delta
  rho2_1  = 1.225 * gas_constants[1]/gas_constants[2] - delta
  v1_1    = zero(delta)
  v2_1    = zero(delta)
  p_1     = 101325

  # Region II
  rho1_2  = 1.225-delta
  rho2_2  = delta
  v1_2    = zero(delta)
  v2_2    = zero(delta)
  p_2     = 101325

  # Region III
  rho1_3  = 1.6861 - delta
  rho2_3  = delta
  v1_3    = -113.5243
  v2_3    = zero(delta)
  p_3     = 159060

  # Set up Region I & II:
  inicenter = SVector(zero(delta), zero(delta))
  x_norm = x[1] - inicenter[1]
  y_norm = x[2] - inicenter[2]
  r = sqrt(x_norm^2 + y_norm^2)

  if (x[1] > 0.50)
    # Set up Region III
    rho1    = rho1_3
    rho2    = rho2_3
    v1      = v1_3
    v2      = v2_3
    p       = p_3
  elseif (r < 0.25)
    # Set up Region I
    rho1    = rho1_1
    rho2    = rho2_1
    v1      = v1_1
    v2      = v2_1
    p       = p_1
  else
    # Set up Region II
    rho1    = rho1_2
    rho2    = rho2_2
    v1      = v1_2
    v2      = v2_2
    p       = p_2
  end

  return prim2cons(SVector(v1, v2, p, rho1, rho2), equations)
end

"""
    initial_condition_shock_bubble(x, t, equations::CompressibleEulerMulticomponentEquations2D{6, 3})

Adaption of the shock-bubble testcase for multicomponent Euler equations to 3 components
- Ayoub Gouasmi, Karthik Duraisamy, Scott Murman
  Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations
  [arXiv: 1904.00972](https://arxiv.org/abs/1904.00972)
"""
function initial_condition_shock_bubble(x, t, equations::CompressibleEulerMulticomponentEquations2D{6, 3})
  # bubble test case, see Gouasmi et al. https://arxiv.org/pdf/1904.00972
  # other reference: https://www.researchgate.net/profile/Pep_Mulet/publication/222675930_A_flux-split_algorithm_applied_to_conservative_models_for_multicomponent_compressible_flows/links/568da54508aeaa1481ae7af0.pdf
  # adapted to 3 component testcase
  # typical domain is rectangular, we change it to a square, as Trixi can only do squares

  @unpack gas_constants = equations

  # Positivity Preserving Parameter, can be set to zero if scheme is positivity preserving
  delta   = 0.05

  # Region Ia
  rho1_1a   = delta
  rho2_1a   = 1.225 * gas_constants[1]/gas_constants[2] - delta
  rho3_1a   = delta
  v1_1a     = zero(delta)
  v2_1a     = zero(delta)
  p_1a      = 101325

  # Region Ib
  rho1_1b   = delta
  rho2_1b   = delta
  rho3_1b   = 1.225 * gas_constants[1]/gas_constants[3] - delta
  v1_1b     = zero(delta)
  v2_1b     = zero(delta)
  p_1b      = 101325

  # Region II
  rho1_2    = 1.225-delta
  rho2_2    = delta
  rho3_2    = delta
  v1_2      = zero(delta)
  v2_2      = zero(delta)
  p_2       = 101325

  # Region III
  rho1_3    = 1.6861 - delta
  rho2_3    = delta
  rho3_3    = delta
  v1_3      = -113.5243
  v2_3      = zero(delta)
  p_3       = 159060

  # Set up Region I & II:
  inicenter_a = SVector(zero(delta), 1.8)
  x_norm_a    = x[1] - inicenter_a[1]
  y_norm_a    = x[2] - inicenter_a[2]
  r_a         = sqrt(x_norm_a^2 + y_norm_a^2)

  inicenter_b = SVector(zero(delta), -1.8)
  x_norm_b    = x[1] - inicenter_b[1]
  y_norm_b    = x[2] - inicenter_b[2]
  r_b         = sqrt(x_norm_b^2 + y_norm_b^2)

  if (x[1] > 0.5)
    # Set up Region III
    rho1    = rho1_3
    rho2    = rho2_3
    rho3    = rho3_3
    v1      = v1_3
    v2      = v2_3
    p       = p_3
  elseif (r_a < 0.25)
    # Set up Region I
    rho1    = rho1_1a
    rho2    = rho2_1a
    rho3    = rho3_1a
    v1      = v1_1a
    v2      = v2_1a
    p       = p_1a
  elseif (r_b < 0.25)
    rho1    = rho1_1b
    rho2    = rho2_1b
    rho3    = rho3_1b
    v1      = v1_1b
    v2      = v2_1b
    p       = p_1b
  else
    # Set up Region II
    rho1    = rho1_2
    rho2    = rho2_2
    rho3    = rho3_2
    v1      = v1_2
    v2      = v2_2
    p       = p_2
  end

  return prim2cons(SVector(v1, v2, p, rho1, rho2, rho3), equations)
end

"""
    initial_condition_shock_bubble(x, t, equations::CompressibleEulerMulticomponentEquations2D{10, 7})

Adaption of the shock-bubble testcase for multicomponent Euler equations to 7 components
- Ayoub Gouasmi, Karthik Duraisamy, Scott Murman
  Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations
  [arXiv: 1904.00972](https://arxiv.org/abs/1904.00972)
"""
function initial_condition_shock_bubble(x, t, equations::CompressibleEulerMulticomponentEquations2D{10, 7})
  # bubble test case, see Gouasmi et al. https://arxiv.org/pdf/1904.00972
  # other reference: https://www.researchgate.net/profile/Pep_Mulet/publication/222675930_A_flux-split_algorithm_applied_to_conservative_models_for_multicomponent_compressible_flows/links/568da54508aeaa1481ae7af0.pdf
  # adapted to 7 component testcase
  # typical domain is rectangular, we change it to a square, as Trixi can only do squares

  @unpack gas_constants = equations

  # Positivity Preserving Parameter, can be set to zero if scheme is positivity preserving
  delta   = 0.05

  # Region Ia
  rho1_1a   = delta
  rho2_1a   = 1.274 * gas_constants[1]/gas_constants[2] - delta
  rho3_1a   = delta
  rho4_1a   = delta
  rho5_1a   = delta
  rho6_1a   = delta
  rho7_1a   = delta
  v1_1a     = zero(delta)
  v2_1a     = zero(delta)
  p_1a      = 101325

  # Region Ib
  rho1_1b   = delta
  rho2_1b   = delta
  rho3_1b   = 1.274 * gas_constants[1]/gas_constants[3] - delta
  rho4_1b   = delta
  rho5_1b   = delta
  rho6_1b   = delta
  rho7_1b   = delta
  v1_1b     = zero(delta)
  v2_1b     = zero(delta)
  p_1b      = 101325

  # Region Ic
  rho1_1c   = delta
  rho2_1c   = delta
  rho3_1c   = delta
  rho4_1c   = 1.274 * gas_constants[1]/gas_constants[4] - delta
  rho5_1c   = delta
  rho6_1c   = delta
  rho7_1c   = delta
  v1_1c     = zero(delta)
  v2_1c     = zero(delta)
  p_1c      = 101325

  # Region Id
  rho1_1d   = delta
  rho2_1d   = delta
  rho3_1d   = delta
  rho4_1d   = delta
  rho5_1d   = 1.274 * gas_constants[1]/gas_constants[5] - delta
  rho6_1d   = delta
  rho7_1d   = delta
  v1_1d     = zero(delta)
  v2_1d     = zero(delta)
  p_1d      = 101325

  # Region Ie
  rho1_1e   = delta
  rho2_1e   = delta
  rho3_1e   = delta
  rho4_1e   = delta
  rho5_1e   = delta
  rho6_1e   = 1.274 * gas_constants[1]/gas_constants[6] - delta
  rho7_1e   = delta
  v1_1e     = zero(delta)
  v2_1e     = zero(delta)
  p_1e      = 101325

  # Region If
  rho1_1f   = delta
  rho2_1f   = delta
  rho3_1f   = delta
  rho4_1f   = delta
  rho5_1f   = delta
  rho6_1f   = delta
  rho7_1f   = 1.274 * gas_constants[1]/gas_constants[7] - delta
  v1_1f     = zero(delta)
  v2_1f     = zero(delta)
  p_1f      = 101325

  # Region II
  rho1_2    = 1.225-delta
  rho2_2    = delta
  rho3_2    = delta
  rho4_2    = delta
  rho5_2    = delta
  rho6_2    = delta
  rho7_2    = delta
  v1_2      = zero(delta)
  v2_2      = zero(delta)
  p_2       = 101325

  # Region III
  rho1_3    = 1.6861 - delta
  rho2_3    = delta
  rho3_3    = delta
  rho4_3    = delta
  rho5_3    = delta
  rho6_3    = delta
  rho7_3    = delta
  v1_3      = -113.5243
  v2_3      = zero(delta)
  p_3       = 159060

  # Set up Region I & II:
  #inicenter_a = SVector(zero(delta), 1.8)
  inicenter_a = SVector(zero(delta), 2.55)
  x_norm_a    = x[1] - inicenter_a[1]
  y_norm_a    = x[2] - inicenter_a[2]
  r_a         = sqrt(x_norm_a^2 + y_norm_a^2)

  #inicenter_b = SVector(zero(delta), -1.8)
  inicenter_b = SVector(zero(delta), 1.60)
  x_norm_b    = x[1] - inicenter_b[1]
  y_norm_b    = x[2] - inicenter_b[2]
  r_b         = sqrt(x_norm_b^2 + y_norm_b^2)

  inicenter_c = SVector(zero(delta), 0.65)
  x_norm_c    = x[1] - inicenter_c[1]
  y_norm_c    = x[2] - inicenter_c[2]
  r_c         = sqrt(x_norm_c^2 + y_norm_c^2)

  inicenter_d = SVector(zero(delta), -0.65)
  x_norm_d    = x[1] - inicenter_d[1]
  y_norm_d    = x[2] - inicenter_d[2]
  r_d         = sqrt(x_norm_d^2 + y_norm_d^2)

  inicenter_e = SVector(zero(delta), -1.60)
  x_norm_e    = x[1] - inicenter_e[1]
  y_norm_e    = x[2] - inicenter_e[2]
  r_e         = sqrt(x_norm_e^2 + y_norm_e^2)

  inicenter_f = SVector(zero(delta), -2.55)
  x_norm_f    = x[1] - inicenter_f[1]
  y_norm_f    = x[2] - inicenter_f[2]
  r_f         = sqrt(x_norm_f^2 + y_norm_f^2)

  if (x[1] > 0.4)
    # Set up Region III
    rho1    = rho1_3
    rho2    = rho2_3
    rho3    = rho3_3
    rho4    = rho4_3
    rho5    = rho5_3
    rho6    = rho6_3
    rho7    = rho7_3
    v1      = v1_3
    v2      = v2_3
    p       = p_3
  elseif (r_a < 0.25)
    # Set up Region I
    rho1    = rho1_1a
    rho2    = rho2_1a
    rho3    = rho3_1a
    rho4    = rho4_1a
    rho5    = rho5_1a
    rho6    = rho6_1a
    rho7    = rho7_1a
    v1      = v1_1a
    v2      = v2_1a
    p       = p_1a
  elseif (r_b < 0.25)
    # Set up Region I
    rho1    = rho1_1b
    rho2    = rho2_1b
    rho3    = rho3_1b
    rho4    = rho4_1b
    rho5    = rho5_1b
    rho6    = rho6_1b
    rho7    = rho7_1b
    v1      = v1_1b
    v2      = v2_1b
    p       = p_1b
  elseif (r_c < 0.25)
    # Set up Region I
    rho1    = rho1_1c
    rho2    = rho2_1c
    rho3    = rho3_1c
    rho4    = rho4_1c
    rho5    = rho5_1c
    rho6    = rho6_1c
    rho7    = rho7_1c
    v1      = v1_1c
    v2      = v2_1c
    p       = p_1c
  elseif (r_d < 0.25)
    # Set up Region I
    rho1    = rho1_1d
    rho2    = rho2_1d
    rho3    = rho3_1d
    rho4    = rho4_1d
    rho5    = rho5_1d
    rho6    = rho6_1d
    rho7    = rho7_1d
    v1      = v1_1d
    v2      = v2_1d
    p       = p_1d
  elseif (r_e < 0.25)
    # Set up Region I
    rho1    = rho1_1e
    rho2    = rho2_1e
    rho3    = rho3_1e
    rho4    = rho4_1e
    rho5    = rho5_1e
    rho6    = rho6_1e
    rho7    = rho7_1e
    v1      = v1_1e
    v2      = v2_1e
    p       = p_1e
  elseif (r_f < 0.25)
    # Set up Region I
    rho1    = rho1_1f
    rho2    = rho2_1f
    rho3    = rho3_1f
    rho4    = rho4_1f
    rho5    = rho5_1f
    rho6    = rho6_1f
    rho7    = rho7_1f
    v1      = v1_1f
    v2      = v2_1f
    p       = p_1f
  else
    # Set up Region II
    rho1    = rho1_2
    rho2    = rho2_2
    rho3    = rho3_2
    rho4    = rho4_2
    rho5    = rho5_2
    rho6    = rho6_2
    rho7    = rho7_2
    v1      = v1_2
    v2      = v2_2
    p       = p_2
  end

  return prim2cons(SVector(v1, v2, p, rho1, rho2, rho3, rho4, rho5, rho6, rho7), equations)
end

"""
    initial_condition_weak_blast_wave(x, t, equations::CompressibleEulerMulticomponentEquations2D)

A for multicomponent adapted weak blast wave taken from
- Sebastian Hennemann, Gregor J. Gassner (2020)
  A provably entropy stable subcell shock capturing approach for high order split form DG
  [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044)
"""
function initial_condition_weak_blast_wave(x, t, equations::CompressibleEulerMulticomponentEquations2D)
  # From Hennemann & Gassner JCP paper 2020 (Sec. 6.3)
  # Set up polar coordinates
  inicenter         = SVector(0.0, 0.0)
  x_norm            = x[1] - inicenter[1]
  y_norm            = x[2] - inicenter[2]
  r                 = sqrt(x_norm^2 + y_norm^2)
  phi               = atan(y_norm, x_norm)
  sin_phi, cos_phi  = sincos(phi)

  prim_rho          = SVector{ncomponents(equations), real(equations)}(r > 0.5 ? 2^(i-1) * (1-2)/(1-2^ncomponents(equations))*1.0 : 2^(i-1) * (1-2)/(1-2^ncomponents(equations))*1.1691 for i in eachcomponent(equations))

  v1                = r > 0.5 ? 0.0 : 0.1882 * cos_phi
  v2                = r > 0.5 ? 0.0 : 0.1882 * sin_phi
  p                 = r > 0.5 ? 1.0 : 1.245

  prim_other         = SVector{3, real(equations)}(v1, v2, p)

  return prim2cons(vcat(prim_other, prim_rho),equations)
end

# Calculate 1D flux for a single point
@inline function flux(u, orientation::Integer, equations::CompressibleEulerMulticomponentEquations2D)
  rho_v1, rho_v2, rho_e  = u

  rho = density(u, equations)

  v1    = rho_v1/rho
  v2    = rho_v2/rho
  gamma = totalgamma(u, equations)
  p     = (gamma - 1) * (rho_e - 0.5 * rho * (v1^2 + v2^2))

  if orientation == 1
    f_rho = densities(u, v1, equations)
    f1  = rho_v1 * v1 + p
    f2  = rho_v2 * v1
    f3  = (rho_e + p) * v1
  else
    f_rho = densities(u, v2, equations)
    f1  = rho_v1 * v2
    f2  = rho_v2 * v2 + p
    f3  = (rho_e + p) * v2
  end

  f_other  = SVector{3, real(equations)}(f1, f2, f3)

  return vcat(f_other, f_rho)
end

"""
    flux_chandrashekar(u_ll, u_rr, orientation, equations::CompressibleEulerMulticomponentEquations2D)

Entropy conserving two-point flux by
- Ayoub Gouasmi, Karthik Duraisamy (2020)
  "Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations""
  arXiv:1904.00972v3 [math.NA] 4 Feb 2020
"""
@inline function flux_chandrashekar(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerMulticomponentEquations2D)
  # Unpack left and right state
  @unpack gammas, gas_constants, cv = equations
  rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll
  rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr
  rhok_mean   = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i+3], u_rr[i+3]) for i in eachcomponent(equations))
  rhok_avg    = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i+3] + u_rr[i+3]) for i in eachcomponent(equations))

  # Iterating over all partial densities
  rho_ll      = density(u_ll, equations)
  rho_rr      = density(u_rr, equations)

  # extract velocities
  v1_ll       = rho_v1_ll/rho_ll
  v2_ll       = rho_v2_ll/rho_ll
  v1_rr       = rho_v1_rr/rho_rr
  v2_rr       = rho_v2_rr/rho_rr
  v1_avg      = 0.5 * (v1_ll + v1_rr)
  v2_avg      = 0.5 * (v2_ll + v2_rr)
  v1_square   = 0.5 * (v1_ll^2 + v1_rr^2)
  v2_square   = 0.5 * (v2_ll^2 + v2_rr^2)
  v_sum       = v1_avg + v2_avg

  enth      = zero(v_sum)
  help1_ll  = zero(v1_ll)
  help1_rr  = zero(v1_rr)

  for i in eachcomponent(equations)
    enth      += rhok_avg[i] * gas_constants[i]
    help1_ll  += u_ll[i+3] * cv[i]
    help1_rr  += u_rr[i+3] * cv[i]
  end

  T_ll        = (rho_e_ll - 0.5 * rho_ll * (v1_ll^2 + v2_ll^2)) / help1_ll
  T_rr        = (rho_e_rr - 0.5 * rho_rr * (v1_rr^2 + v2_rr^2)) / help1_rr
  T           = 0.5 * (1.0/T_ll + 1.0/T_rr)
  T_log       = ln_mean(1.0/T_ll, 1.0/T_rr)

  # Calculate fluxes depending on orientation
  help1       = zero(T_ll)
  help2       = zero(T_rr)
  if orientation == 1
    f_rho       = SVector{ncomponents(equations), real(equations)}(rhok_mean[i]*v1_avg for i in eachcomponent(equations))
    for i in eachcomponent(equations)
      help1     += f_rho[i] * cv[i]
      help2     += f_rho[i]
    end
    f1 = (help2) * v1_avg + enth/T
    f2 = (help2) * v2_avg
    f3 = (help1)/T_log - 0.5 * (v1_square + v2_square) * (help2) + v1_avg * f1 + v2_avg * f2
  else
    f_rho       = SVector{ncomponents(equations), real(equations)}(rhok_mean[i]*v2_avg for i in eachcomponent(equations))
    for i in eachcomponent(equations)
      help1     += f_rho[i] * cv[i]
      help2     += f_rho[i]
    end
    f1 = (help2) * v1_avg
    f2 = (help2) * v2_avg + enth/T
    f3 = (help1)/T_log - 0.5 * (v1_square + v2_square) * (help2) + v1_avg * f1 + v2_avg * f2
  end
  f_other  = SVector{3, real(equations)}(f1, f2, f3)

  return vcat(f_other, f_rho)
end

# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation
@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerMulticomponentEquations2D)
  rho_v1_ll, rho_v2_ll, rho_e_ll = u_ll
  rho_v1_rr, rho_v2_rr, rho_e_rr = u_rr

  # Calculate primitive variables and speed of sound
  rho_ll   = density(u_ll, equations)
  rho_rr   = density(u_rr, equations)
  gamma_ll = totalgamma(u_ll, equations)
  gamma_rr = totalgamma(u_rr, equations)

  v1_ll = rho_v1_ll / rho_ll
  v2_ll = rho_v2_ll / rho_ll
  v_mag_ll = sqrt(v1_ll^2 + v2_ll^2)
  p_ll = (gamma_ll - 1) * (rho_e_ll - 1/2 * rho_ll * v_mag_ll^2)
  c_ll = sqrt(gamma_ll * p_ll / rho_ll)

  v1_rr = rho_v1_rr / rho_rr
  v2_rr = rho_v2_rr / rho_rr
  v_mag_rr = sqrt(v1_rr^2 + v2_rr^2)
  p_rr = (gamma_rr - 1) * (rho_e_rr - 1/2 * rho_rr * v_mag_rr^2)
  c_rr = sqrt(gamma_rr * p_rr / rho_rr)

  λ_max = max(v_mag_ll, v_mag_rr) + max(c_ll, c_rr)
end

@inline function max_abs_speeds(u, equations::CompressibleEulerMulticomponentEquations2D)
  rho_v1, rho_v2, rho_e = u

  rho   = density(u, equations)
  v1    = rho_v1 / rho
  v2    = rho_v2 / rho

  gamma = totalgamma(u, equations)
  p     = (gamma - 1) * (rho_e - 1/2 * rho * (v1^2 + v2^2))
  c     = sqrt(gamma * p / rho)

  return (abs(v1) + c, abs(v2) + c, )
end

# Convert conservative variables to primitive
@inline function cons2prim(u, equations::CompressibleEulerMulticomponentEquations2D)
  rho_v1, rho_v2, rho_e = u

  prim_rho = SVector{ncomponents(equations), real(equations)}(u[i+3] for i in eachcomponent(equations))

  rho   = density(u, equations)
  v1    = rho_v1 / rho
  v2    = rho_v2 / rho
  gamma = totalgamma(u, equations)
  p     = (gamma - 1) * (rho_e - 0.5 * rho * (v1^2 + v2^2))
  prim_other =  SVector{3, real(equations)}(v1, v2, p)

  return vcat(prim_other, prim_rho)
end

# Convert conservative variables to entropy
@inline function cons2entropy(u, equations::CompressibleEulerMulticomponentEquations2D)
  @unpack cv, gammas, gas_constants = equations
  rho_v1, rho_v2, rho_e = u

  rho       = density(u, equations)

  v1        = rho_v1 / rho
  v2        = rho_v2 / rho
  v_square  = v1^2 + v2^2
  gamma     = totalgamma(u, equations)
  p         = (gamma - 1) * (rho_e - 0.5 * rho * v_square)
  s         = log(p) - gamma*log(rho)
  rho_p     = rho / p

  # Multicomponent stuff
  help1 = zero(v1)

  for i in eachcomponent(equations)
    help1 += u[i+3] * cv[i]
  end

  T         = (rho_e - 0.5 * rho * v_square) / (help1)

  entrop_rho  = SVector{ncomponents(equations), real(equations)}( -1.0 * (cv[i] * log(T) - gas_constants[i] * log(u[i+3])) + gas_constants[i] + cv[i] - (v_square / (2*T)) for i in eachcomponent(equations))

  w1        = v1/T
  w2        = v2/T
  w3        = -1.0/T

  entrop_other = SVector{3, real(equations)}(w1, w2, w3)

  return vcat(entrop_other, entrop_rho)
end

# Convert primitive to conservative variables
@inline function prim2cons(prim, equations::CompressibleEulerMulticomponentEquations2D)
  @unpack cv, gammas = equations
  v1, v2, p = prim

  cons_rho = SVector{ncomponents(equations), real(equations)}(prim[i+3] for i in eachcomponent(equations))
  rho     = density(prim, equations)
  gamma   = totalgamma(prim, equations)

  rho_v1  = rho * v1
  rho_v2  = rho * v2
  rho_e   = p/(gamma-1) + 0.5 * (rho_v1 * v1 + rho_v2 * v2)

  cons_other = SVector{3, real(equations)}(rho_v1, rho_v2, rho_e)

  return vcat(cons_other, cons_rho)
end

"""
    totalgamma(u, equations::CompressibleEulerMulticomponentEquations2D)

Function that calculates the total gamma out of all partial gammas using the
partial density fractions as well as the partial specific heats at constant volume.
"""
@inline function totalgamma(u, equations::CompressibleEulerMulticomponentEquations2D)
  @unpack cv, gammas = equations

  help1 = zero(u[1])
  help2 = zero(u[1])

  for i in eachcomponent(equations)
    help1 += u[i+3] * cv[i] * gammas[i]
    help2 += u[i+3] * cv[i]
  end

  return help1/help2
end

@inline function density_pressure(u, equations::CompressibleEulerMulticomponentEquations2D)
  rho_v1, rho_v2, rho_e = u

  rho          = density(u, equations)
  gamma        = totalgamma(u, equations)
  rho_times_p  = (gamma - 1) * (rho * rho_e - 0.5 * (rho_v1^2 + rho_v2^2))

  return rho_times_p
end

@inline function density(u, equations::CompressibleEulerMulticomponentEquations2D)
  rho = zero(u[1])

  for i in eachcomponent(equations)
    rho += u[i+3]
  end

  return rho
 end

 @inline function densities(u, v, equations::CompressibleEulerMulticomponentEquations2D)

  return SVector{ncomponents(equations), real(equations)}(u[i+3]*v for i in eachcomponent(equations))
 end