text:
function calculateResults(LDDEP::LDDEProblem, method::DiscretizationMethod, DiscretizationLength::Real; n_steps::Int64 = nStepOfLength(DiscretizationLength, method.Δt), calculate_additive::Bool = false, im::ItoIsometryMethod = Trapezoidal(20, method))
    result = Result(LDDEP, method, DiscretizationLength, n_steps = n_steps, calculate_additive = calculate_additive, im = im)
    calculateDetResults!(result)

    calculateStResults!(result)
    return (result)
end

function calculateStResults!(rst::AbstractResult)
    registerStResult!(rst,
    [rst.method(mx, rst) for mx in [rst.LDDEP.αs...,rst.LDDEP.βs...]])
    if rst.calculate_additive
        calculateStAdditiveResults!(rst)
    end
end
function calculateStAdditiveResults!(rst::AbstractResult)
    registerStAdditiveResults!(rst,
        rst.method.(rst.LDDEP.σs, Ref(rst)))
end
function registerStResult!(rst::AbstractResult, stsubmxs::Vector{<:Vector{<:stSubMX}})
    rst.stsubMXs .= stsubmxs;
end

function registerStAdditiveResults!(rst::AbstractResult, subvs::Vector{<:Vector{<:stSubV}})
    rst.stsubVs .= subvs;
end

function stDiscreteMappingSteps(rst::AbstractResult{d}) where d
    ns = calculate_noise_mxelems(rst)
    Iss = [[Vector{Int64}(undef, n) for n in ns] for x in 1:rst.n_steps]
    Jss = [[Vector{Int64}(undef, n) for n in ns] for x in 1:rst.n_steps]
    Vss = [[Vector{eltype(rst.stsubMXs[1][1].MXfun[1])}(undef, n) for n in ns] for x in 1:rst.n_steps]
    n0 = [fill(1, length(ns)) for i in 1:rst.n_steps]
    n1 = [zero(ns) for i in 1:rst.n_steps]
    for stsmx in rst.stsubMXs
        w = stsmx[1].nID
        for t in 1:rst.n_steps
            for i in eachindex(stsmx[t].ranges, stsmx[t].MXfun)
                n1[t][w] = n0[t][w] + d^2 - 1
                IJs = Iterators.product(stsmx[t].ranges[i]...) |> collect |> vec
                Iss[t][w][n0[t][w]:n1[t][w]] .= getindex.(IJs, 1)
                Jss[t][w][n0[t][w]:n1[t][w]] .= getindex.(IJs, 2)
                Vss[t][w][n0[t][w]:n1[t][w]] .= stsmx[t].MXfun[i][:]
                n0[t][w] = n1[t][w] + 1
            end
        end
    end
    if rst.calculate_additive
        ns = calculate_noise_velems(rst)
        IssV = [[Vector{Int64}(undef, n) for n in ns] for x in 1:rst.n_steps]
        VssV = [[Vector{eltype(rst.stsubVs[1][1].Vfun)}(undef, n) for n in ns] for x in 1:rst.n_steps]
        n0 = [fill(1, length(ns)) for i in 1:rst.n_steps]
        n1 = [zero(ns) for i in 1:rst.n_steps]
        for stsv in rst.stsubVs
            w = stsv[1].nID
            for t in 1:rst.n_steps
                n1[t][w] = n0[t][w] + d - 1
                IssV[t][w][n0[t][w]:n1[t][w]] .= Base.OneTo(d)
                VssV[t][w][n0[t][w]:n1[t][w]] .= stsv[t].Vfun
                n0[t][w] = n1[t][w] + 1
            end
        end
    else
        IssV = fill(fill(Vector{Int64}(undef, 0), length(ns)), rst.n_steps)
        VssV = fill(fill(Vector{Vector{Float64}}(undef, 0), length(ns)), rst.n_steps)
    end
    ([[sparse(Iss[t][w], Jss[t][w], Vss[t][w], rst.n, rst.n) for w in 1:rst.LDDEP.w] for t in 1:rst.n_steps], [[sparsevec(IssV[t][w], VssV[t][w], rst.n) for w in 1:rst.LDDEP.w] for t in 1:rst.n_steps]) # stMXs, stVs
end

function stDiscreteMapping(rst::AbstractResult)
    stDiscreteMapping(DiscreteMappingSteps(rst)...,
        stDiscreteMappingSteps(rst)...)
end

DiscreteMapping_M1(rst::AbstractResult) = DiscreteMapping_M1(DiscreteMappingSteps(rst)...)
DiscreteMapping_M1(rst::AbstractResult, idxs::AbstractArray{<:Integer}) = DiscreteMapping_M1(DiscreteMappingSteps(rst)..., idxs)
function DiscreteMapping_M1(LDDEP::LDDEProblem, method::DiscretizationMethod, DiscretizationLength::Real; args...)
    DiscreteMapping_M1(calculateResults(LDDEP, method, DiscretizationLength; args...))
end
function DiscreteMapping_M1(LDDEP::LDDEProblem, method::DiscretizationMethod, DiscretizationLength::Real, idxs; args...)
    DiscreteMapping_M1(calculateResults(LDDEP, method, DiscretizationLength; args...), idxs)
end

function DiscreteMapping_M2(stdm::stDiscreteMapping, rst::AbstractResult)
    M2_MXs = M2_Mapping_from_Sparse.(stdm.detMXs, Ref(rst)) .+
        [sum(M2_Mapping_from_Sparse(stMX, rst) for stMX in stMXsₜ) for stMXsₜ in stdm.stMXs] # H [time]
    M1toM2_MXs = fill(spzeros(eltype(eltype(M2_MXs)), size(M2_MXs, 1), rst.n), rst.n_steps)
    M2_Vs = fill(spzeros(eltype(eltype(M2_MXs)), size(M2_MXs, 1)), rst.n_steps)
    if rst.calculate_additive
        M1toM2_MXs .= M1toM2_Mapping_Generator_from_Sparse.(stdm.detMXs, stdm.detVs, Ref(rst)) .+
            [sum(M1toM2_Mapping_Generator_from_Sparse(stMX, stdm.stVs[t][i], rst) for (i, stMX) in enumerate(stMXsₜ)) for (t, stMXsₜ) in enumerate(stdm.stMXs)]
        M2_Vs .= M2_Additive_from_Sparse.(stdm.detVs, Ref(rst)) +
            [sum(M2_Additive_from_Sparse(stV, rst) for stV in stVₜ) for stVₜ in stdm.stVs]
    end
    DiscreteMapping_M2(stdm.ts, stdm.detMXs, stdm.detVs, M2_MXs, M1toM2_MXs, M2_Vs)
end

function DiscreteMapping_M2(stdm::stDiscreteMapping, idxs::AbstractArray{<:Integer}, rst::AbstractResult)
    M2_MXs = M2_Mapping_from_Sparse.(stdm.detMXs, Ref(idxs), Ref(rst))
    for t in eachindex(M2_MXs)
        M2_MXs[t] += sum(M2_Mapping_from_Sparse(stMX, idxs, rst) for stMX in stdm.stMXs[t])  # H [time]
    end
    M1toM2_MXs = fill(spzeros(eltype(eltype(M2_MXs)), size(M2_MXs, 1), length(idxs)), rst.n_steps)
    M2_Vs = fill(spzeros(eltype(eltype(M2_MXs)), size(M2_MXs, 1)), rst.n_steps)
    if rst.calculate_additive
        M1toM2_MXs .= M1toM2_Mapping_Generator_from_Sparse.(stdm.detMXs, stdm.detVs, Ref(idxs), Ref(rst)) .+
            [sum(M1toM2_Mapping_Generator_from_Sparse(stMX, stdm.stVs[t][i], idxs, rst) for (i, stMX) in enumerate(stMXsₜ)) for (t, stMXsₜ) in enumerate(stdm.stMXs)]
        M2_Vs .= M2_Additive_from_Sparse.(stdm.detVs, Ref(idxs), Ref(rst)) .+
            [sum(M2_Additive_from_Sparse(stV, idxs, rst) for stV in stVₜ) for stVₜ in stdm.stVs]
    end
    DiscreteMapping_M2(stdm.ts, getindex.(stdm.detMXs, Ref(idxs), Ref(idxs)), getindex.(stdm.detVs, Ref(idxs)), M2_MXs, M1toM2_MXs, M2_Vs)
end

DiscreteMapping_M2(rst::AbstractResult) = DiscreteMapping_M2(stDiscreteMapping(rst), rst)
DiscreteMapping_M2(rst::AbstractResult, idxs::AbstractArray{<:Integer}) = DiscreteMapping_M2(stDiscreteMapping(rst), idxs, rst)
function DiscreteMapping_M2(LDDEP::LDDEProblem, method::DiscretizationMethod, DiscretizationLength::Real; args...)
    DiscreteMapping_M2(calculateResults(LDDEP, method, DiscretizationLength; args...))
end
function DiscreteMapping_M2(LDDEP::LDDEProblem, method::DiscretizationMethod, DiscretizationLength::Real, idxs; args...)
    DiscreteMapping_M2(calculateResults(LDDEP, method, DiscretizationLength; args...), idxs)
end

# Some attributes of the mappings
function fixPointOfMapping(dm::DiscreteMapping_M1)
    (I - prodl(dm.M1_MXs)) \ Vector(reduce_additive(dm.M1_MXs, dm.M1_Vs))
end
function fixPointOfMapping(dm::DiscreteMapping_M1, idxs::AbstractVector{<:Integer})
    (I - prodl(dm.M1_MXs,idxs)) \ Vector(reduce_additive(dm.M1_MXs, dm.M1_Vs, idxs))
end

function fixPointOfMapping(dm::DiscreteMapping_M2)
    v1st = (I - prodl(dm.M1_MXs)) \ Vector(reduce_additive(dm.M1_MXs, dm.M1_Vs))
    (I - prodl(dm.M2_MXs)) \ Vector(reduce_additive(dm.M2_MXs, dm.M1toM2_MXs, dm.M1_MXs, dm.M2_Vs, dm.M1_Vs, v1st))
end

function spectralRadiusOfMapping(dm::DiscreteMapping_M1; args...)
    abs(eigs(prodl(dm.M1_MXs); args...)[1][1])
end
function spectralRadiusOfMapping(dm::DiscreteMapping_M1, idxs::AbstractVector{<:Integer}; args...)
    abs(eigs(prodl(dm.M1_MXs, idxs); args...)[1][1])
end

function spectralRadiusOfMapping(dm::DiscreteMapping_M2; args...)
    abs(eigs(prodl(dm.M2_MXs); args...)[1][1])
end


text:
using Statistics
using Dates
using PyPlot, Colors
using DataFrames

include("startup.jl")
include("manualMassLibrary.jl")

# *****************************    DEFINE PLOT TIME    *********************************************************************
selectedTime = true
toi = [Dates.DateTime(2018,5,10), Dates.DateTime(2018,5,22)]	# have to restart julia for new toi #
interpolateTraces = true
avMinutes = 30
intpolT = collect(Dates.datetime2unix(toi[1]):(60*avMinutes):Dates.datetime2unix(toi[2])) 	# 5min averages
intpolDT = Dates.unix2datetime.(intpolT)

# ********************************    LOAD VARS    ***************************************************************************
plotSMPS = true
plotNO3cims = true
plotACSM = true
plotMeteo = true
plotSRR = true
plotBC = false
#plotBLinfluence = false

if !(@isdefined createTotalPTR3data)
	println("load variables:")
	include("Bolivia-DefineVariables.jl")
else println("variables already defined!")
end

# ***************************** Define different masking times ***************************************************************

PacificFT = true

if PacificFT == true
	maskFT = (
			(DateTime(2018,2,19,02,25) .< intpolDT .< DateTime(2018,2,19,10,01))
			.+ (DateTime(2018,4,5,00,25) .< intpolDT .< DateTime(2018,4,5,09,01))
			.+ (DateTime(2018,5,11,07,20) .< intpolDT .< DateTime(2018,5,11,09,10))
			.+ (DateTime(2018,5,15,03,55) .< intpolDT .< DateTime(2018,5,15,08,31))
			.+ (DateTime(2018,5,18,04,25) .< intpolDT .< DateTime(2018,5,18,08,31))
			.+ (DateTime(2018,5,20,00,25) .< intpolDT .< DateTime(2018,5,20,08,01))
			)
	maskFT = findall(x -> x > 0, maskFT)
	maskBLnight = (
			(DateTime(2018,2,19,20,25) .< intpolDT .< DateTime(2018,2,20,05,01))
			.+ (DateTime(2018,4,6,00,25) .< intpolDT .< DateTime(2018,4,6,07,01))
			.+ (DateTime(2018,5,18,22,55) .< intpolDT .< DateTime(2018,5,19,03,01))
			.+ (DateTime(2018,5,15,23,55) .< intpolDT .< DateTime(2018,5,16,02,31))
			.+ (DateTime(2018,5,17,01,55) .< intpolDT .< DateTime(2018,5,17,08,01))
			.+ (DateTime(2018,5,22,00,55) .< intpolDT .< DateTime(2018,5,22,07,01))
			)
	maskBLnight = findall(x -> x > 0, maskBLnight)
	maskBLday = (
			(DateTime(2018,2,19,12,00) .< intpolDT .< DateTime(2018,2,19,20,01))
			.+ (DateTime(2018,4,5,09,25) .< intpolDT .< DateTime(2018,4,5,15,01))
			.+ (DateTime(2018,5,11,10,55) .< intpolDT .< DateTime(2018,5,11,16,01))
			.+ (DateTime(2018,5,12,11,55) .< intpolDT .< DateTime(2018,5,12,16,01))
			.+ (DateTime(2018,5,14,11,55) .< intpolDT .< DateTime(2018,5,14,16,01))
			.+ (DateTime(2018,5,15,11,55) .< intpolDT .< DateTime(2018,5,15,16,31))
			.+ (DateTime(2018,5,16,10,55) .< intpolDT .< DateTime(2018,5,16,16,31))
			.+ (DateTime(2018,5,20,11,55) .< intpolDT .< DateTime(2018,5,20,18,31))
			)
	maskBLday = findall(x -> x > 0, maskBLday)
elseif AmazonFT == true
maskFT = (
		(DateTime(2018,2,17,10,16) .< intpolDT .< DateTime(2018,2,17,09,31))
		)
maskFT = findall(x -> x > 0, maskFT)
maskBLnight = (
		(DateTime(2018,5,18,22,55) .< intpolDT .< DateTime(2018,5,19,03,01))
		.+ (DateTime(2018,5,15,23,55) .< intpolDT .< DateTime(2018,5,16,02,31))
		.+ (DateTime(2018,5,17,01,55) .< intpolDT .< DateTime(2018,5,17,08,01))
		.+ (DateTime(2018,5,22,00,55) .< intpolDT .< DateTime(2018,5,22,07,01))
		)
maskBLnight = findall(x -> x > 0, maskBLnight)
maskBLday = (
		(DateTime(2018,5,11,10,55) .< intpolDT .< DateTime(2018,5,11,16,01))
		.+ (DateTime(2018,5,12,11,55) .< intpolDT .< DateTime(2018,5,12,16,01))
		.+ (DateTime(2018,5,14,11,55) .< intpolDT .< DateTime(2018,5,14,16,01))
		.+ (DateTime(2018,5,16,10,55) .< intpolDT .< DateTime(2018,5,16,16,31))
		)
maskBLday = findall(x -> x > 0, maskBLday)
elseif MSAduringBL == true
	maskFT = ()
	maskFT = findall(x -> x > 0, maskFT)
	maskBLnight = (
			(DateTime(2018,5,18,22,55) .< intpolDT .< DateTime(2018,5,19,03,01))
			.+ (DateTime(2018,5,15,23,55) .< intpolDT .< DateTime(2018,5,16,02,31))
			.+ (DateTime(2018,5,17,01,55) .< intpolDT .< DateTime(2018,5,17,08,01))
			.+ (DateTime(2018,5,22,00,55) .< intpolDT .< DateTime(2018,5,22,07,01))
			)
	maskBLnight = findall(x -> x > 0, maskBLnight)
	maskBLday = (
			(DateTime(2018,2,17,10,55) .< intpolDT .< DateTime(2018,2,17,17,01))
			.+ (DateTime(2018,2,20,11,55) .< intpolDT .< DateTime(2018,2,20,16,01))
			.+ (DateTime(2018,3,15,11,55) .< intpolDT .< DateTime(2018,3,17,16,01))
			.+ (DateTime(2018,3,20,10,55) .< intpolDT .< DateTime(2018,3,20,16,31))
			.+ (DateTime(2018,3,26,11,55) .< intpolDT .< DateTime(2018,3,16,16,01))
			.+ (DateTime(2018,3,27,10,55) .< intpolDT .< DateTime(2018,3,27,16,31))
			)
end

# **************************** Select Variables and make histograms *********************************************************

# ****************************** make METEO figure *****************************************************************************************************************
fig = plt.figure() # new figure
WDstd_bp = [WDstdv[maskFT], WDstdv[maskBLnight], WDstdv[maskBLday]]
ax = fig.add_subplot(321)
ax.boxplot(WDstd_bp)
ax.set_ylabel("stdev(WD) [°]")
ax.set_xticklabels(["clean", "RL_n", "BL_day"])

WD_bp = [WD[maskFT], WD[maskBLnight], WD[maskBLday]]
ax = fig.add_subplot(322)
ax.boxplot(WD_bp)
ax.set_ylabel("WD [°]")
ax.set_xticklabels(["clean", "RL_n", "BL_day"])

WS_bp = [WS[maskFT], WS[maskBLnight], WS[maskBLday]]
ax = fig.add_subplot(323)
ax.boxplot(WS_bp)
ax.set_ylabel("WS [m/s]")
ax.set_xticklabels(["clean", "RL_n", "BL_day"])

Temp_bp = [TEMPstation[maskFT], TEMPstation[maskBLnight], TEMPstation[maskBLday]]
ax = fig.add_subplot(324)
ax.boxplot(Temp_bp)
ax.set_ylabel("Temp [°C]")
ax.set_xticklabels(["clean", "RL_n", "BL_day"])

RH_bp = [RHstation[maskFT], RHstation[maskBLnight], RHstation[maskBLday]]
ax = fig.add_subplot(325)
ax.boxplot(RH_bp)
ax.set_ylabel("RH [‰]")
ax.set_xticklabels(["clean", "RL_n", "BL_day"])

PRESSURE_bp = [PRESSUREstation[maskFT], PRESSUREstation[maskBLnight], PRESSUREstation[maskBLday]]
ax = fig.add_subplot(326)
ax.boxplot(PRESSURE_bp)
ax.set_ylabel("pressure [mbar]")
ax.set_xticklabels(["clean", "RL_n", "BL_day"])

fig.tight_layout(pad=1.5)
PyPlot.savefig("/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/FreeTrop_BL/meteo.png")

# ****************************** make particle figure *****************************************************************************************************************
fig = plt.figure() # new figure
ACSMsum = (ACSM_org.+ACSM_sulf.+ACSM_nitr.+ACSM_ammon.+ACSM_chlor)
ACSMsum_bp = [ACSMsum[maskFT], ACSMsum[maskBLnight], ACSMsum[maskBLday]]
ax = fig.add_subplot(221)
ax.boxplot(ACSMsum_bp)
ax.legend(["ACSMsum"])
ax.set_ylabel("particle mass [ug/m³]")
ax.set_yscale("log")
ax.set_xticklabels(["clean", "RL_n", "BL_day"])

ax = fig.add_subplot(222)
ACSM_FT = transpose(mean(ACSM_all[maskFT,:]./ACSMsum[maskFT], dims = 1))[:,1]
ACSM_BLn = transpose(mean(ACSM_all[maskBLnight,:]./ACSMsum[maskBLnight], dims = 1))[:,1]
ACSM_BLd = transpose(mean(ACSM_all[maskBLday,:]./ACSMsum[maskBLday], dims = 1))[:,1]
barWidth = 0.2
r2 = 1:length(ACSM_FT)
r1 = r2.-barWidth
r3 = r2.+barWidth
ax.bar(r1, ACSM_FT, width=barWidth, label = "clean")
ax.bar(r2, ACSM_BLn, width=barWidth, label = "RL night")
ax.bar(r3, ACSM_BLd, width=barWidth, label = "BL day")
ax.legend()
ax.set_ylabel("share [%]")
ax.set_xticks(r2)
ax.set_xticklabels(["Sulf", "Org", "Nitr", "Amm", "Cl"])

SMPS_totconc_bp = [SMPS_totconc[maskFT], SMPS_totconc[maskBLnight], SMPS_totconc[maskBLday]]
ax = fig.add_subplot(223)
ax.boxplot(SMPS_totconc_bp)
ax.legend(["SMPS_totconc"])
ax.set_ylabel("particle Nr.")
ax.set_yscale("log")
ax.set_xticklabels(["clean", "RL_n", "BL_day"])

ax = fig.add_subplot(224)
SMPS_binconc_FT = mean(SMPS_binconc[maskFT,:], dims=1)
SMPS_binconc_BLday = mean(SMPS_binconc[maskBLday,:], dims=1)
SMPS_binconc_BLnight = mean(SMPS_binconc[maskBLnight,:], dims=1)
ax.plot(SMPS_binsize, transpose(SMPS_binconc_FT), label = "clean")
ax.plot(SMPS_binsize, transpose(SMPS_binconc_BLnight), label = "BL night")
ax.plot(SMPS_binsize, transpose(SMPS_binconc_BLday), label = "BL day")
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_ylabel("dN/dlog(Dp)")
ax.legend()
fig.tight_layout(pad=1.5)
PyPlot.savefig("/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/FreeTrop_BL/particles.png")

# ****************************** make NitrateCIMS figure *****************************************************************************************************************
fig = plt.figure()  # new figure
Nitrophenol_bp = [Nitrophenol[maskFT], Nitrophenol[maskBLnight], Nitrophenol[maskBLday]]
ax = fig.add_subplot(221)
ax.boxplot(Nitrophenol_bp)
ax.legend(["Nitrophenol"])
ax.set_ylabel("signal")
ax.set_yscale("log")
ax.set_xticklabels(["clean", "RL_n", "BL_day"])

SO5trace_bp = [SO5trace[maskFT], SO5trace[maskBLnight], SO5trace[maskBLday]]
ax = fig.add_subplot(222)
ax.boxplot(SO5trace_bp)
ax.legend(["SO5"])
ax.set_ylabel("signal")
ax.set_xticklabels(["clean", "RL_n", "BL_day"])

MSAtrace_bp = [MSAtrace[maskFT], MSAtrace[maskBLnight], MSAtrace[maskBLday]]
ax = fig.add_subplot(223)
ax.boxplot(MSAtrace_bp)
ax.set_ylabel("signal")
ax.legend(["MSA"])
ax.set_xticklabels(["clean", "RL_n", "BL_day"])

SAtrace_bp = [SAtrace[maskFT], SAtrace[maskBLnight][.!(isnan.(SAtrace[maskBLnight]))], SAtrace[maskBLday]]
ax = fig.add_subplot(224)
ax.boxplot(SAtrace_bp)
ax.set_ylabel("signal")
ax.legend(["SA"])
ax.set_xticklabels(["clean", "RL_n", "BL_day"])

fig.tight_layout(pad=1.5)

PyPlot.savefig("/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/FreeTrop_BL/nitrateCIMSdata.png")

# ****************************** make SRR figure *****************************************************************************************************************
fig = plt.figure()  # new figure
SRRage_bp = [SRRage[maskFT], SRRage[maskBLnight], SRRage[maskBLday]]
ax = fig.add_subplot(221)
ax.boxplot(SRRage_bp)
ax.legend(["average age [h]"], loc = 2)
ax.set_xticklabels(["clean", "RL_n", "BL_day"])

SRRclockdir_bp = [SRRclockdir[maskFT], SRRclockdir[maskBLnight], SRRclockdir[maskBLday]]
ax = fig.add_subplot(222)
ax.boxplot(SRRclockdir_bp)
ax.legend(["clock dir. Wind"], loc = 2)
ax.set_xticklabels(["clean", "RL_n", "BL_day"])

#SRRheight_bp = [SRRheight[maskFT], SRRheight[maskBLnight], SRRheight[maskBLday]]
#ax = fig.add_subplot(223)
#ax.boxplot(SRRheight_bp)
#ax.legend(["av. hasl [m]"], loc = 2)
#ax.set_xticklabels(["clean", "RL_n", "BL_day"])

SRRconcAll_bp =  [SRRconcAll[maskFT], SRRconcAll[maskBLnight], SRRconcAll[maskBLday]]
ax = fig.add_subplot(224)
ax.boxplot(SRRconcAll_bp)
ax.legend(["total 'particles' in model"], loc = 2)
ax.set_xticklabels(["clean", "RL_n", "BL_day"])

SRR_bl_bp =  [SRR_bl[maskFT], SRR_bl[maskBLnight], SRR_bl[maskBLday]]
ax = fig.add_subplot(223)
ax.boxplot(SRR_bl_bp)
ax.legend(["influence from surface layer"], loc = 2)
ax.set_xticklabels(["clean", "RL_n", "BL_day"])

fig.tight_layout(pad=1.5)

PyPlot.savefig("/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/FreeTrop_BL/SRRs_overview.png")

# ****************************** make single SRR figures ********************************************************************************************************
# SRR plot 6 regions
figure()
SRRs_06_FT = transpose(mean(SRRs_06[maskFT,:], dims = 1))[:,1]
SRRs_06_BLn = transpose(mean(SRRs_06[maskBLnight,:], dims = 1))[:,1]
SRRs_06_BLd = transpose(mean(SRRs_06[maskBLday,:], dims = 1))[:,1]
barWidth = 0.2
r2 = 1:length(SRRs_06_FT)
r1 = r2.-barWidth
r3 = r2.+barWidth
PyPlot.bar(r1, SRRs_06_FT, width=barWidth, label = "clean")
PyPlot.bar(r2, SRRs_06_BLn, width=barWidth, label = "RL night")
PyPlot.bar(r3, SRRs_06_BLd, width=barWidth, label = "BL day")
legend()
ylabel("share [%]")
xticks(r2, SRR06names)
PyPlot.savefig("/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/FreeTrop_BL/SRRs_06.png")

# SRR plot all 18 regions
figure(figsize = (15,6))
SRRs_FT = transpose(mean(SRRs[maskFT,:], dims = 1))[:,1]
SRRs_BLn = transpose(mean(SRRs[maskBLnight,:], dims = 1))[:,1]
SRRs_BLd = transpose(mean(SRRs[maskBLday,:], dims = 1))[:,1]
barWidth = 0.2
r2 = 1:length(SRRs_FT)
r1 = r2.-barWidth
r3 = r2.+barWidth
PyPlot.bar(r1, SRRs_FT, width=barWidth, label = "clean")
PyPlot.bar(r2, SRRs_BLn, width=barWidth, label = "RL night")
PyPlot.bar(r3, SRRs_BLd, width=barWidth, label = "BL day")
legend()
ylabel("share [%]")
xticks(r2, SRRheader)
PyPlot.savefig("/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/FreeTrop_BL/SRRs_18.png")

# SRR plot divide by LR, MR, SM, SR
#=
figure()
SRRs_FT_range = [sum(SRRs_FT[findall(x -> occursin("LR", x), SRRheader)]), sum(SRRs_FT[findall(x -> occursin("MR", x), SRRheader)]), sum(SRRs_FT[findall(x -> occursin("SM", x), SRRheader)]), sum(SRRs_FT[findall(x -> occursin("SR", x), SRRheader)])]
SRRs_BLn_range = [sum(SRRs_BLn[findall(x -> occursin("LR", x), SRRheader)]), sum(SRRs_BLn[findall(x -> occursin("MR", x), SRRheader)]), sum(SRRs_BLn[findall(x -> occursin("SM", x), SRRheader)]), sum(SRRs_BLn[findall(x -> occursin("SR", x), SRRheader)])]
SRRs_BLd_range = [sum(SRRs_BLd[findall(x -> occursin("LR", x), SRRheader)]), sum(SRRs_BLd[findall(x -> occursin("MR", x), SRRheader)]), sum(SRRs_BLd[findall(x -> occursin("SM", x), SRRheader)]), sum(SRRs_BLd[findall(x -> occursin("SR", x), SRRheader)])]
barWidth = 0.2
r2 = 1:length(SRRs_FT_range)
r1 = r2.-barWidth
r3 = r2.+barWidth
PyPlot.bar(r1, SRRs_FT_range, width=barWidth, label = "clean")
PyPlot.bar(r2, SRRs_BLn_range, width=barWidth, label = "RL night")
PyPlot.bar(r3, SRRs_BLd_range, width=barWidth, label = "BL day")
ylabel("share [%]")
legend()
xticks(r2, ["LR", "MR", "SM", "SR"])
PyPlot.savefig("/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/FreeTrop_BL/SRRs_range.png")
=#
# ****************************** make PTR3 figures ********************************************************************************************************
# PTR3 data:

if !(@isdefined ptr3Traces)
	println("  -> load PTR3 Data")
	ptr3Peaks, ptr3Times, ptr3Traces = createTotalPTR3data()
	ptr3Masses = round.(ptr3Peaks[:,1]; digits = 3)
end

removeBG_mintype = true

if removeBG_mintype == true
	background = DatasetFunctions.nanmin(ptr3Traces, 1)
	background[background .< 0] .= 0
	ptr3Traces = ptr3Traces .- background
end

function ptr3Boxplot(fig, ax = subplot(111), loglin = "log" ; Cnrs=0:1:1, Onrs=0:1:1, Hnrs=0:2:10, Nnrs=0:1, Snrs=0:1, legString="")
	massesToPlot = round.([
			H[1]
			MasslistFunctions.createMassList2(; C=Cnrs, O=Onrs, H=Hnrs, N=Nnrs, S=Snrs, nHplus=1)[1]
			# always exclude C10H12O, C10H13O, C9H11ON, C6H12O4, C10H12O2, when looking at 12.5. - 15.5. !!!
			# always exclude C4H8 and C8H18O2 and C8H20O2 - Butanol !!!
			#MasslistFunctions.createMassList2(; C=1, O=1:5, H=1:10, N=0, S=0, nHplus=1, allowRadicals=true)[1]
			#MasslistFunctions.createMassList2(; C=9, O=1, H=10:8:18, N=0, S=0, nHplus=1, allowRadicals=true)[1]
			#MasslistFunctions.createMassList2(; C=7, O=1, H=12:2:14, N=0, S=0, nHplus=1, allowRadicals=true)[1]
			]; digits = 3)
	wantedIndices = findall(x -> x in massesToPlot, ptr3Masses)
	println(" exclude butanol influenced peaks")
	ptr3Traces[:, wantedIndices[findall(x -> (x in ["C4H8", "C4H10O", "C8H20O2", "C9H12", "C10H12O", "C10H13O", "C9H11ON", "C6H12O4", "C10H12O2", "C3H6O"]), ptr3Peaks[wantedIndices,2])]] .= NaN
	println("  -> interpolate chosen PTR3 data")
	ptr3Traces2plot = InterpolationFunctions.interpolateMatrix_1D(intpolT, ptr3Times, ptr3Traces[:,wantedIndices])
	ptr3sum = DatasetFunctions.nansum(ptr3Traces2plot,2)

	ptr3_bp = [ptr3sum[maskFT][.!(isnan.(ptr3sum[maskFT]))], ptr3sum[maskBLnight][.!(isnan.(ptr3sum[maskBLnight]))], ptr3sum[maskBLday][.!(isnan.(ptr3sum[maskBLday]))]]

	ax.boxplot(ptr3_bp)
	if (legString == "")
		if Nnrs != 0 && Snrs != 0
			legString = ["C$(Cnrs)H$(Hnrs)O$(Onrs)N$(Nnrs)S$(Snrs)"]
		elseif Nnrs != 0
			legString = ["C$(Cnrs)H$(Hnrs)O$(Onrs)N$(Nnrs)"]
		elseif Snrs != 0
			legString = ["C$(Cnrs)H$(Hnrs)O$(Onrs)S$(Snrs)"]
		else
			legString = ["C$(Cnrs)H$(Hnrs)O$(Onrs)"]
		end
	end
	ax.legend(legString, loc = 2)
	ax.set_xticklabels(["clean", "RL_n", "BL_day"])
	saveString = replace(legString[1], ":" => "-")
	ax.set_ylabel("signal [ppt]")
	ax.set_yscale(loglin)
end

function ptr3BoxplotNoO(fig, loglin = "log" ; Cnrs=0:1:1, Onrs=0:1:1, Hnrs=0:2:10, Nnrs=0:1, Snrs=0:1)
	massesToPlot = round.([
			H[1]
			MasslistFunctions.createMassList2(; C=Cnrs, O=Onrs, H=Hnrs, N=Nnrs, S=Snrs, nHplus=1)[1]
			# always exclude C10H12O, C10H13O, C9H11ON, C6H12O4, C10H12O2, when looking at 12.5. - 15.5. !!!
			# always exclude C4H8 and C8H20O2 - Butanol !!!
			#MasslistFunctions.createMassList2(; C=1, O=1:5, H=1:10, N=0, S=0, nHplus=1, allowRadicals=true)[1]
			#MasslistFunctions.createMassList2(; C=9, O=1, H=10:8:18, N=0, S=0, nHplus=1, allowRadicals=true)[1]
			#MasslistFunctions.createMassList2(; C=7, O=1, H=12:2:14, N=0, S=0, nHplus=1, allowRadicals=true)[1]
			]; digits = 3)
	wantedIndices = findall(x -> x in massesToPlot, ptr3MassesNoO)
	println(" exclude butanol influenced peaks")
	ptr3Traces[:, wantedIndices[findall(x -> (x in ["C4H8", "C4H10O", "C8H20O2", "C9H12", "C10H12O", "C10H13O", "C9H11ON", "C6H12O4", "C10H12O2", "C3H6O"]), ptr3Peaks[wantedIndices,2])]] .= NaN
	println("  -> interpolate chosen PTR3 data")
	ptr3Traces2plot = InterpolationFunctions.interpolateMatrix_1D(intpolT, ptr3TimesNoO, ptr3TracesNoO[:,wantedIndices])
	ptr3sum = DatasetFunctions.nansum(ptr3Traces2plot,2)

	ptr3_bp = [ptr3sum[maskFT][.!(isnan.(ptr3sum[maskFT]))], ptr3sum[maskBLnight][.!(isnan.(ptr3sum[maskBLnight]))], ptr3sum[maskBLday][.!(isnan.(ptr3sum[maskBLday]))]]

	ax = subplot(111)
	ax.boxplot(ptr3_bp)
	if Nnrs != 0 && Snrs != 0
		legString = ["C$(Cnrs)H$(Hnrs)O$(Onrs)N$(Nnrs)S$(Snrs)"]
	elseif Nnrs != 0
		legString = ["C$(Cnrs)H$(Hnrs)O$(Onrs)N$(Nnrs)"]
	elseif Snrs != 0
		legString = ["C$(Cnrs)H$(Hnrs)O$(Onrs)S$(Snrs)"]
	else
		legString = ["C$(Cnrs)H$(Hnrs)O$(Onrs)"]
	end
	ax.legend(legString, loc = 2)
	ax.set_xticklabels(["clean", "RL_n", "BL_day"])
	saveString = replace(legString[1], ":" => "-")
	ax.set_ylabel("signal [ppt]")
	ax.set_yscale(loglin)
end

function CHO_barplot(peaks, traces, timemask; Omin = 1)
	figure()
	ax = subplot(111)
	rename!(peaks, Symbol.(["mass", "name", "C", "C13", "H", "Hplus", "N", "O", "O18", "S"]))
	signal = DatasetFunctions.nanmean(traces[timemask,:],1)
	#signal[signal.<0] .= 0
	signal[isnan.(signal)] .= 0
	signal[findall(x -> (x in ["C4H8", "C4H10O", "C8H20O2", "C9H12", "C10H12O", "C10H13O", "C9H11ON", "C6H12O4", "C10H12O2", "C3H6O"]), peaks[:,2])] .= 0
	peaks.Signal = transpose(signal)[:,1]
	peaksCHO = sort!(peaks[peaks.N .== 0, :],[:O])
	peaksCHON = sort!(peaks[peaks.N .> 0, :],[:O])
	peaksCHO[peaksCHO.O .== 0, :]
	CnrSignalArray = zeros(29,20)
	#cols = range(colorant"red", stop=colorant"green", length=17)
	cols = distinguishable_colors(20, [RGB(1,1,1)])[2:end]
	pcols = map(col -> (red(col), green(col), blue(col)), cols)
	#for c in 1:29
	#	CnrSignalArray[c,1] = sum((peaksCHO[peaksCHO.O .== 0, :]).Signal[(peaksCHO[peaksCHO.O .== 0, :]).C .== c])
	#end
	#PyPlot.bar(collect(1:29), CnrSignalArray[:,1], color = pcols[1], label = "CxHz", align="center")
	for o in collect(Omin:maximum(peaksCHO.O))
		for c in 1:29
			CnrSignalArray[c, o] = sum((peaksCHO[peaksCHO.O .== o, :]).Signal[(peaksCHO[peaksCHO.O .== o, :]).C .== c])
		end
		PyPlot.bar(collect(1:29), CnrSignalArray[:,o], color = pcols[o], bottom=sum(CnrSignalArray[:,1:o-1], dims = 2)[:,1], align="center", label = "CxHyO$(o)")
	end
	legend()
end

function CHON_barplot(peaks, traces, timemask; Omin = 1)
	figure()
	ax = subplot(111)
	rename!(peaks, Symbol.(["mass", "name", "C", "C13", "H", "Hplus", "N", "O", "O18", "S"]))
	signal = DatasetFunctions.nanmean(traces[timemask,:],1)
	#signal[signal.<0] .= 0
	signal[isnan.(signal)] .= 0
	signal[findall(x -> (x in ["C4H8", "C4H10O", "C8H20O2", "C9H12", "C10H12O", "C10H13O", "C9H11ON", "C6H12O4", "C10H12O2", "C3H6O"]), peaks[:,2])] .= 0
	peaks.Signal = transpose(signal)[:,1]
	peaksCHON = sort!(peaks[peaks.N .> 0, :],[:O])
	peaksCHON[peaksCHON.O .== 0, :]
	CnrSignalArray = zeros(29,20)
	#cols = range(colorant"red", stop=colorant"green", length=17)
	cols = distinguishable_colors(20, [RGB(1,1,1)])[2:end]
	pcols = map(col -> (red(col), green(col), blue(col)), cols)
	for o in collect(Omin:maximum(peaksCHON.O))
		for c in 1:29
			CnrSignalArray[c, o] = sum((peaksCHON[peaksCHON.O .== o, :]).Signal[(peaksCHON[peaksCHON.O .== o, :]).C .== c]) # includes NH4+
		end
		PyPlot.bar(collect(1:29), CnrSignalArray[:,o], color = pcols[o], bottom=sum(CnrSignalArray[:,1:o-1], dims = 2)[:,1], align="center", label = "CxHyO$(o)N")
	end
	legend()
end

#=
CHO_barplot(DataFrame(ptr3Peaks), ptr3Traces, maskFT; Omin = 1)
title("free troposphere")
CHO_barplot(DataFrame(ptr3Peaks), ptr3Traces, maskBLnight; Omin = 1)
title("residual layer night")
CHO_barplot(DataFrame(ptr3Peaks), ptr3Traces, maskBLday; Omin = 1)
title("boundary layer day")

CHON_barplot(DataFrame(ptr3Peaks), ptr3Traces, maskFT; Omin = 3)
title("free troposphere")
CHON_barplot(DataFrame(ptr3Peaks), ptr3Traces, maskBLnight; Omin = 3)
title("residual layer night")
CHON_barplot(DataFrame(ptr3Peaks), ptr3Traces, maskBLday; Omin = 3)
title("boundary layer day")

saveString = "FT_noN"
fig = figure(saveString, figsize = (12,16))
ptr3Boxplot(fig, subplot(231), "linear"; Cnrs=2, Hnrs=6, Onrs=0, Nnrs=0, Snrs=1)
ptr3Boxplot(fig, subplot(232), "linear"; Cnrs=2, Hnrs=6, Onrs=1, Nnrs=0, Snrs=1)
ptr3Boxplot(fig, subplot(233), "linear"; Cnrs=2, Hnrs=6, Onrs=2, Nnrs=0, Snrs=1)
#ptr3Boxplot(fig, subplot(233), "linear"; Cnrs=2, Hnrs=2:8, Onrs=1:10, Nnrs=0, Snrs=0)
ptr3Boxplot(fig, subplot(234), "linear"; Cnrs=3, Hnrs=6, Onrs=2, Nnrs=0, Snrs=0)
ptr3Boxplot(fig, subplot(235), "linear"; Cnrs=4, Hnrs=8, Onrs=2, Nnrs=0, Snrs=0)
ptr3Boxplot(fig, subplot(236), "linear"; Cnrs=4, Hnrs=10, Onrs=2, Nnrs=0, Snrs=0)
fig.tight_layout(pad=1.5)
fig.savefig(string("/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/FreeTrop_BL/PTR3_" , saveString , ".png"))

saveString = "BL_noN"
fig = figure(saveString, figsize = (12,16))
ptr3Boxplot(fig, subplot(231), "linear"; Cnrs=5, Hnrs=6:8, Onrs=1:10, Nnrs=0, Snrs=0)
ptr3Boxplot(fig, subplot(232), "linear"; Cnrs=6, Hnrs=10, Onrs=1:10, Nnrs=0, Snrs=0)
ptr3Boxplot(fig, subplot(233), "linear"; Cnrs=7, Hnrs=12, Onrs=1:10, Nnrs=0, Snrs=0)
ptr3Boxplot(fig, subplot(234), "linear"; Cnrs=8, Hnrs=12, Onrs=1:10, Nnrs=0, Snrs=0)
ptr3Boxplot(fig, subplot(235), "linear"; Cnrs=9, Hnrs=14, Onrs=1:10, Nnrs=0, Snrs=0)
ptr3Boxplot(fig, subplot(236), "linear"; Cnrs=18, Hnrs=12:40, Onrs=3:10, Nnrs=0, Snrs=0)
fig.tight_layout(pad=1.5)
savefig(string("/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/FreeTrop_BL/PTR3_" , saveString , ".png"))

saveString = "FTvsBL_N"
fig = figure(saveString, figsize = (12,16))
ptr3Boxplot(fig, subplot(231), "linear"; Cnrs=2, Hnrs=5, Onrs=2, Nnrs=1, Snrs=0)
ptr3Boxplot(fig, subplot(232), "linear"; Cnrs=4, Hnrs=5, Onrs=1:10, Nnrs=1, Snrs=0)
ptr3Boxplot(fig, subplot(233), "linear"; Cnrs=4, Hnrs=7, Onrs=1:10, Nnrs=1, Snrs=0)
ptr3Boxplot(fig, subplot(234), "linear"; Cnrs=5, Hnrs=7, Onrs=1:10, Nnrs=1, Snrs=0)
ptr3Boxplot(fig, subplot(235), "linear"; Cnrs=7, Hnrs=11, Onrs=1:10, Nnrs=1, Snrs=0)
ptr3Boxplot(fig, subplot(236), "linear"; Cnrs=9, Hnrs=13, Onrs=1:10, Nnrs=1, Snrs=0)
fig.tight_layout(pad=1.5)
savefig(string("/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/FreeTrop_BL/PTR3_" , saveString , ".png"))

saveString = "BL_noN_all_log"
fig = figure(saveString, figsize = (12,16))
ptr3Boxplot(fig, subplot(111), "log"; Cnrs=5:18, Hnrs=6:50, Onrs=1:10, Nnrs=0, Snrs=0)
savefig(string("/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/FreeTrop_BL/PTR3_" , saveString , ".png"))

saveString = "BL_N_all_log"
fig = figure(saveString, figsize = (12,16))
ptr3Boxplot(fig, subplot(111), "log"; Cnrs=5:18, Hnrs=6:50, Onrs=1:10, Nnrs=1, Snrs=0)
savefig(string("/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/FreeTrop_BL/PTR3_" , saveString , ".png"))

# ### ###

PyPlot.figure()
PyPlot.bar(CnrsArray, totalMatrixCHO[:,1], color = pcols[1], label = "CxHyO1", align="center")
for o in collect(2:14)
	PyPlot.bar(CnrsArray, totalMatrixCHO[:,o], color = pcols[o], bottom=sum(totalMatrixCHO[:,1:o-1],2), align="center", label = "CxHyO$(o)")
end
axis("tight")
title("Contribution of diff CHOs, with O>$(Omin) during night 10-11.05.18")
grid("on")
PyPlot.xlabel("# C")
PyPlot.ylabel("concentration [ppt]")
PyPlot.legend()

PyPlot.figure()
PyPlot.bar(CnrsArray, totalMatrixCHON[:,1], color = pcols[1], label = "CxHyO1N", align="center")
for o in collect(2:14)
	PyPlot.bar(CnrsArray, totalMatrixCHON[:,o], color = pcols[o], bottom=sum(totalMatrixCHON[:,1:o-1],2), align="center", label = "CxHyO$(o)N")
end
axis("tight")
title("Contribution of diff CHONs, with O>$(Omin) during night 10-11.05.18")
grid("on")
PyPlot.xlabel("# C")
PyPlot.ylabel("concentration [ppt]")
PyPlot.legend()
=#

# *********** prepare nice plot for BAMS paper ********************
function BAMSfigure()
	saveString = "BAMS_marineFigure"
	N=3
	fig = figure(saveString, figsize = (12,16))
	ax1 = subplot(4,1,1)
	wantedIndices = findall(x -> x in ["C2H6S", "C2H5S", "C2H4OS", "CH4SO2", "C2H6OS"], ptr3Peaks[:,2])
	ptr3Traces2plot = InterpolationFunctions.interpolateMatrix_1D(intpolT, ptr3Times, ptr3Traces[:,wantedIndices])
	semilogy(intpolDT, ptr3Traces2plot./DatasetFunctions.nanmean(ptr3Traces2plot,1))
	semilogy(intpolDT, MSAtrace./DatasetFunctions.nanmean(MSAtrace), "--")
	semilogy(intpolDT, SAtrace./DatasetFunctions.nanmean(SAtrace), "--")
	legend(vcat(ptr3Peaks[wantedIndices,2], ["MSA", "SA"]))
# sulfur compounds row
	MSAtrace_bp = [MSAtrace[maskFT], MSAtrace[maskBLnight], MSAtrace[maskBLday]]
	ax = fig.add_subplot(4,3,4)
	ax.boxplot(MSAtrace_bp)
	ax.set_ylabel("signal")
	ax.legend(["MSA"])
	ax.set_yscale("log")
	ax.set_xticklabels(["clean", "RL_n", "BL_day"])

	SAtrace_bp = [SAtrace[maskFT], SAtrace[maskBLnight][.!(isnan.(SAtrace[maskBLnight]))], SAtrace[maskBLday]]
	ax = fig.add_subplot(4,3,5)
	ax.boxplot(SAtrace_bp)
	ax.set_ylabel("signal")
	ax.legend(["SA"])
	ax.set_yscale("log")
	ax.set_xticklabels(["clean", "RL_n", "BL_day"])

	ptr3Boxplot(fig, subplot(4,3,6), "log"; Cnrs=1:2, Hnrs=4:6, Onrs=0:1, Nnrs=0:1, Snrs=1, legString = ["\$ \\sum { C_{1,2}H_{4-6}O_{0,1}S_{1} }\$"])

# BL compounds row
	Nitrophenol_bp = [Nitrophenol[maskFT], Nitrophenol[maskBLnight], Nitrophenol[maskBLday]]
	ax = fig.add_subplot(4,3,7)
	ax.boxplot(Nitrophenol_bp)
	ax.legend(["Nitrophenol"])
	ax.set_ylabel("signal")
	ax.set_yscale("log")
	ax.set_xticklabels(["clean", "RL_n", "BL_day"])

	ptr3Boxplot(fig, subplot(4,3,8), "log"; Cnrs=5:18, Hnrs=6:50, Onrs=1:20, Nnrs=0:1, Snrs=0, legString = ["\$\\sum {C_{>4}H_{>6}O_{x}N_{0,1} }\$"])

	AH_bp = [AHstation[maskFT], AHstation[maskBLnight], AHstation[maskBLday]]
	ax = fig.add_subplot(4,3,9)
	ax.boxplot(AH_bp)
	ax.set_ylabel("AH [‰]")
	ax.set_xticklabels(["clean", "RL_n", "BL_day"])

# particles row
	ax = fig.add_subplot(4,3,10)
	SMPS_binconc_FT = mean(SMPS_binconc[maskFT,:], dims=1)
	SMPS_binconc_BLday = mean(SMPS_binconc[maskBLday,:], dims=1)
	SMPS_binconc_BLnight = mean(SMPS_binconc[maskBLnight,:], dims=1)
	ax.plot(SMPS_binsize, transpose(SMPS_binconc_FT), label = "clean")
	ax.plot(SMPS_binsize, transpose(SMPS_binconc_BLnight), label = "BL night")
	ax.plot(SMPS_binsize, transpose(SMPS_binconc_BLday), label = "BL day")
	ax.set_xscale("log")
	ax.set_yscale("log")
	ax.set_ylabel("dN/dlog(Dp)")
	ax.legend()
	fig.tight_layout(pad=1.5)

	ACSMsum = (ACSM_org.+ACSM_sulf.+ACSM_nitr.+ACSM_ammon.+ACSM_chlor)
	ACSMsum_bp = [ACSMsum[maskFT], ACSMsum[maskBLnight], ACSMsum[maskBLday]]
	ax = fig.add_subplot(4,3,11)
	ax.boxplot(ACSMsum_bp)
	ax.legend(["ACSMsum"])
	ax.set_ylabel("particle mass [ug/m³]")
	ax.set_yscale("log")
	ax.set_xticklabels(["clean", "RL_n", "BL_day"])

	ax = fig.add_subplot(4,3,12)
	ACSM_FT = transpose(mean(ACSM_all[maskFT,:]./ACSMsum[maskFT], dims = 1))[:,1]
	ACSM_BLn = transpose(mean(ACSM_all[maskBLnight,:]./ACSMsum[maskBLnight], dims = 1))[:,1]
	ACSM_BLd = transpose(mean(ACSM_all[maskBLday,:]./ACSMsum[maskBLday], dims = 1))[:,1]
	barWidth = 0.2
	r2 = 1:length(ACSM_FT)
	r1 = r2.-barWidth
	r3 = r2.+barWidth
	ax.bar(r1, ACSM_FT, width=barWidth, label = "clean")
	ax.bar(r2, ACSM_BLn, width=barWidth, label = "RL night")
	ax.bar(r3, ACSM_BLd, width=barWidth, label = "BL day")
	ax.legend()
	ax.set_ylabel("share [%]")
	ax.set_xticks(r2)
	ax.set_xticklabels(["Sulf", "Org", "Nitr", "Amm", "Cl"])

	savefig(string("/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/FreeTrop_BL/" , saveString , ".png"))
	savefig(string("/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/FreeTrop_BL/" , saveString , ".eps"))

end

#=
n=2
println("Plot Nr. $(n)")
SRRplot = subplot(N,1,n, sharex = ax1, label="SRR")
cmap = get_cmap("tab20")
colors = cmap(0:6)
SRRplot.patch.set_alpha(0)
SRRplot.semilogy(intpolDT, sum(SRRs_06[:,findall(x -> x in ["07_PW", "08_PW"],SRR06names)], dims = 2)./sum(SRRs_06,dims = 2))
SRRplot.semilogy(intpolDT, SRR_z0)
=#
#SRRplot.legend(SRR06names[findall(x -> x in ["07_PW", "08_PW"],SRR06names)]])

ax1 = subplot(211)
ax1.margins(0.05)           # Default margin is 0.05, value 0 means fit
ax1.plot(t1, f(t1))

ax2 = subplot(223)
ax2.margins(2, 2)           # Values >0.0 zoom out
ax2.plot(t1, f(t1))
ax2.set_title("Zoomed out")

ax3 = subplot(224)
ax3.margins(x=0, y=-0.25)   # Values in (-0.5, 0.0) zooms in to center
ax3.plot(t1, f(t1))
ax3.set_title("Zoomed in")

ext: jl


text:
function find_closest_n(
    dist_func::Function, collection;
    n = 10, ordering = Base.Order.Reverse,
    should_break::Nullable{Function} = nothing
)
    result = dist_func(first(collection))
    T = typeof(result)
    heap = BoundedBinaryHeap{T}(n, ordering)
    push!(heap, result)
    is_first = true
    for x in collection
        if is_first
            is_first = false
            continue
        end
        result = dist_func(x)
        isnothing(result) && continue
        push!(heap, result)
        if !isnothing(should_break) && should_break(result)
            break
        end
    end
    results = T[]
    while !isempty(heap)
        push!(results, pop!(heap))
    end
    @view results[end:-1:1]
end
find_closest(dist_func::Function, collection; kwargs...) =
    find_closest_n(dist_func, collection; n = 1, kwargs...)[1]

"""
    data_search_n(collection, x; n = 10, dist = Levenshtein())

Search for `n` closest matches in `collection` for `x`.

Returns a list of (distance, index)
"""
data_search_n(
    collection, x;
    n = 10,
    dist = StringDistances.Hamming()
) = find_closest_n(eachindex(collection), n = n) do i
    dist(x, collection[i]), i
end

data_search(collection, x; kwargs...) = data_search_n(collection, x; n = 1, kwargs...)[1][2]


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

"""A Cluster is a repeating pattern of indices defined in a lattice called the element lattice. The cluster of elements has its own lattice, which by definition is different from the element lattice unless each cluster consists of a single element. A LatticeCluster subtype must implement these methods:

clusterbasis(a::AbstractLatticeCluster) returns the AbstractBasis object for the cluster. This is not generally the same as the basis for the underlying lattice.

elementbasis(a::AbstractLatticeCluster) returns the AbstractBasis object the underlying lattice of the cluster.

clusterelements(a::AbstractLatticeCluster) returns the lattice indices, represented in the underlying lattice basis, for each of the elements in a unit cluster cell.

clustersize(a::AbstractLatticeCluster) returns the number of elements in a unit cluster cell.

For example this defines a Cluster of three hexagonal elements:

```
julia> function hex3cluster()
    clusterelts = SVector((0,0),(-1,0),(-1,1))
    eltlattice = HexBasis1()
    clusterbasis = LatticeBasis(( -1,2),(2,-1))
    return LatticeCluster(clusterbasis,eltlattice,clusterelts)
end

julia> lattice = hex3cluster()
LatticeCluster{3, 2, Float64, LatticeBasis{2, Int64}, HexBasis1{2, Float64}}(LatticeBasis{2, Int64}([-1 2; 2 -1]), HexBasis1{2, Float64}(), [(0, 0), (-1, 0), (-1, 1)])

julia> lattice[1,1]  #returns the locations of the 3 elements in the cluster at the cluster offset of 1,1
2×3 SMatrix{2, 3, Float64, 6} with indices SOneTo(2)×SOneTo(3):
 1.0  -0.5        1.0
 1.0   0.133975  -0.732051

```
"""
abstract type AbstractLatticeCluster{N1,N} end

"""returns the euclidean distance between lattice clusters (norm of the longest lattice basis vector). The companion function `latticediameter` returns the distance in lattice space which does not take into account the size the element basis underlying the cluster."""
euclideandiameter(a::AbstractLatticeCluster) = euclideandiameter(elementbasis(a))*latticediameter(a)
export euclideandiameter

"""Basic lattic cluster type. N1 is the number of tiles in the cluster, N is the dimension."""
struct LatticeCluster{N1, N, T<:Real, B1<:AbstractBasis{N,Int},B2<:AbstractBasis{N,T}} <: AbstractLatticeCluster{N1,N}
    clusterbasis::B1 #this basis defines the offsets of the clusters
    elementbasis::B2 #this basis defines the underlying lattice

    clusterelements::SVector{N1,NTuple{N,Int64}} #vector containing the lattice indices of the elements in the cluster. Each column is one lattice coordinate. These indices are assumed to be offsets from the origin of the lattice.

    LatticeCluster(clusterbasis::B1,eltlattice::B2,clusterelements::SVector{N1,NTuple{N,Int64}}) where{N1,N,T<:Real,B1<:AbstractBasis{N,Int64},B2<:AbstractBasis{N,T}} = new{N1,N,T,B1,B2}(clusterbasis,eltlattice,clusterelements)
end
export LatticeCluster

LatticeCluster(clusterbasis::B1,eltlattice::B2,clusterelements::SMatrix{N,N1,Int64}) where{N1,N,T<:Real,B1<:AbstractBasis{N,Int64},B2<:AbstractBasis{N,T}} = LatticeCluster(clusterbasis,eltlattice,SVector{N1}(reinterpret(reshape,NTuple{N,Int64},clusterelements)))

clusterelements(a::LatticeCluster) = a.clusterelements
export clusterelements
elementbasis(a::LatticeCluster) = a.elementbasis
export elementbasis
clustersize(a::LatticeCluster) = length(a.clusterelements)
clusterbasis(a::LatticeCluster) = a.clusterbasis
export clusterbasis

""" Lattice clusters have integer basis vectors. These represent the lattice in unit steps of the underlying element basis, not the Euclidean distance between lattice cluster centers. For example, the lattice cluster basis vectors for a hex3 lattice are (-1,2),(2,-1), where the units of the basis vectors represent steps in the underlying element basis. To get a Euclidean distance measure between cluster centers you need to multiply by the size of the element basis diameter. In the case of lenslet layout you need to multiply the cluster diameter by the lenslet diameter."""
latticediameter(a::Repeat.AbstractLatticeCluster) =   euclideandiameter(Repeat.basismatrix(Repeat.clusterbasis(a)))
export latticediameter

"""returns the positions of every element in a cluster given the cluster indices"""
function Base.getindex(A::LatticeCluster{N1,N,T,B1,B2}, indices::Vararg{Int, N}) where{N1,N,T,B1<:AbstractBasis{N,Int},B2<:AbstractBasis{N,T}}
    clusteroffset = A.clusterbasis[indices...]
    temp = MMatrix{N,N1,T}(undef)
    for i in 1:N1
        temp[:,i] = A.elementbasis[A.clusterelements[i]...] + clusteroffset
    end
    return SMatrix{N,N1,T}(temp) #positions of the cluster elements, not the lattice coordinates.

end

"""Can access any cluster in a cluster lattice so the range of indices is unlimited"""
function Base.size(::LatticeCluster{N1,N}) where{N1,N}
    return ntuple((i)->Base.IsInfinite(),N)
end

"""returns the integer coordinates of the tile at tileindex in the cluster at for a cluster with location cᵢ,cⱼ"""
function tilecoordinates(cluster::LatticeCluster{N1,N},cᵢ,cⱼ,tileindex) where{N1,N}
    @assert tileindex <= N1
    clustercenter = basismatrix(clusterbasis(cluster))*SVector(cᵢ,cⱼ)
    tileoffset = cluster.clusterelements[tileindex]
    return clustercenter .+ tileoffset
end
export tilecoordinates

Base.setindex!(A::AbstractBasis{N}, v, I::Vararg{Int, N}) where{N} = nothing #can't set lattice points. Might want to throw an exception instead.

""" returns the lattice indices of the elements in the cluster. These are generally not the positions of the elements"""
function clustercoordinates(a::LatticeCluster{N1,N},indices::Vararg{Int,N}) where{N1,N}
    clusteroffset = a.clusterbasis[indices...]
    temp = MMatrix{N,N1,Int}(undef)
    for i in 1:N1
        temp[:,i] = SVector{N,Int}(a.clusterelements[i]) + clusteroffset
    end
    return SMatrix{N,N1,Int}(temp) #the lattice coordinates in the underlying element basis of the cluster elements
end
export clustercoordinates

"""Given the (i,j) coordinates of a tile defined in the the underlying lattice basis of elements of the cluster compute the coordinates (cᵢ,cⱼ) of the cluster containing the tile, and the tile number of the tile in that cluster"""
function cluster_coordinates_from_tile_coordinates(cluster::S,i::Int,j::Int) where{N1,N, S<:AbstractLatticeCluster{N1,N}}
    found = false
    bmatrix = Rational.(basismatrix(clusterbasis(cluster)))

    local clustercoords::SVector{}
    tileindex = 0 #initialize to illegal value
    for coords in eachcol(clustercoordinates(cluster, 0,0)) #don't know which tile in the cluster the i,j coords come from so try all of them until find one that yields an integer value for the cluster coordinate.
        tileindex += 1
        possiblecenter = SVector{N,Rational}((i,j)) .- Rational.(coords)
        clustercoords = bmatrix \  possiblecenter
        if all(1 .== denominator.(clustercoords)) #If coordinates are integer this is the correct cluster value. If coordinates are not integer keep trying.
            found = true    #every tile position i,j corresponds to some cluster (cᵢ,cⱼ) so will always find an answer
            break
        end
    end
    @assert found == true #should always find a cluster corresponding to any i,j. If not then something is seriously wrong.

    return Int64.(clustercoords),tileindex  #return coordinates of the cluster and the index of the tile within that cluster
end
export cluster_coordinates_from_tile_coordinates

cluster_coordinates_from_tile_coordinates(cluster::S,coords::NTuple{2,Int64}) where{N1,N, S<:AbstractLatticeCluster{N1,N}} = cluster_coordinates_from_tile_coordinates(cluster,coords...)

"""
May want to have many properties associated with the elements in a cluster, which is why properties is represented as a DataFrame. The DataFrame in the properties field should have as many rows as there are elements in a cluster. At a minimum it must have a :Color and a :Name column.

Example: a lattice cluster of three hexagonal elements, labeled R,G,B with colors red,green,blue.

```
function hex3RGB()
    clusterelements = SVector((0,0),(-1,0),(-1,1)) #lattice indices of cluster elements represent in the element lattice basis
    colors = [colorant"red",colorant"green",colorant"blue"]
    names = ["R","G","B"]
    eltlattice = HexBasis1()
    clusterbasis = LatticeBasis(( -1,2),(2,-1)) #lattice indices, in the element lattice basis, representing the repeating pattern of the cluster
    lattice = LatticeCluster(clusterbasis,eltlattice,clusterelements)
    properties =  DataFrame(Color = colors, Name = names)
    return ClusterWithProperties(lattice,properties)
end
```
"""
struct ClusterWithProperties{N1,N,T} <: AbstractLatticeCluster{N1,N}
    cluster::LatticeCluster{N1,N,T}
    properties::DataFrame
end
export ClusterWithProperties

Base.getindex(A::ClusterWithProperties{N1,N,T}, indices::Vararg{Int,N}) where{N1,N,T} = cluster(A)[indices...]
Base.setindex!(A::ClusterWithProperties, v, I::Vararg{Int, N}) where{N} = nothing #can't set lattice points. Might want to throw an exception instead.

properties(a::ClusterWithProperties) = a.properties
export properties
cluster(a::ClusterWithProperties) = a.cluster
export cluster
clustercoordinates(a::ClusterWithProperties,indices...) = clustercoordinates(cluster(a),indices...)
elementbasis(a::ClusterWithProperties) = elementbasis(cluster(a))
export elementbasis
clustersize(a::ClusterWithProperties) = clustersize(a.cluster)
export clustersize
clusterbasis(a::ClusterWithProperties) = clusterbasis(a.cluster)
export clusterbasis


text:
using JuLIP, LinearAlgebra, DataDeps
using JuLIP.Potentials, JuLIP.Testing

# ##  EAM

# register(DataDep(
#     "JuLIP_testpots",
#     "A few EAM potentials for testing",
#     "https://www.dropbox.com/s/leub1c9ft1mm9fg/JuLIP_data.zip?dl=1",
#     post_fetch_method = file -> run(`unzip $file`)
#     ))

# test_pots = joinpath(datadep"JuLIP_testpots", "JuLIP_data") * "/"

# @info("Loading some interatomic potentials . .")
# calc = JuLIP.Potentials.EAM(test_pots * "pfe.plt",
#                               test_pots * "ffe.plt",
#                               test_pots * "F_fe.plt")

# # EAM Potential
# at = set_pbc!( bulk(:Fe, cubic = true), false ) * 2

## Si  SW

calc = JuLIP.Potentials.StillingerWeber()
at = set_pbc!( bulk(:Si, cubic = true), false ) * 2

##

@info("With threads")
JuLIP.usethreads!(true)
fdtest(calc, at, verbose=true)

@info("Serial")
JuLIP.usethreads!(false)
fdtest(calc, at, verbose=true)
##

JuLIP.usethreads!(true)
Et = energy(calc, at)
Ft = forces(calc, at)
Vt = virial(calc, at)

JuLIP.usethreads!(false)
Es = energy(calc, at)
Fs = forces(calc, at)
Vs = virial(calc, at)

@info("Compare threaded vs serial values")

@show abs(Et - Es)
@show norm(Vt - Vs)
@show maximum(norm.(Ft - Fs))


text:
using Pkg

include(joinpath(@__DIR__, "routines.jl"))

# We cannot be sure that Julia have already setup it's startup script
include(joinpath(@__DIR__, "..", "userdata", ".julia", "config", "juliawinconfig.jl"))

activate_binary("julia")
activate_binary("juliawin-prompt")
activate_binary("7z")

add_startup_script()

try
   using ArgParse
catch e
   Pkg.add("ArgParse")
   using ArgParse
end

s = ArgParseSettings(description = "This is the main commandline interface to Juliawin.")
@add_arg_table s begin
    "--install"
        help = "Choose specific package to install"
    "--install-dialog", "-d"
        help = "Enter the installation guide."
        action = :store_true
end

parsed_args = parse_args(ARGS, s)

function input(prompt::String="")::String
    print(prompt)
    return chomp(readline())
end

function ask_yn(message)
    answer = nothing
    while true
        answer = lowercase(input(message))
        if answer == "y" || answer == "n"
            return answer
        end
    end
end

function doflagset(app, dict)
    key = "juliawin-install$(app)"

    if key in keys(ENV)
        if lowercase(ENV[key]) == "1"
            return dict[app] = true
        elseif lowercase(ENV[key]) == "0"
            return dict[app] = false
        end
    end

    if (ans = lowercase(ask_yn("Install $app [Y/N]? "))) == "y"
        dict[app] = true
    elseif ans == "n"
        dict[app] = false
    end
end

if parsed_args["install"] !== nothing
    println("This is not implemented yet")

elseif parsed_args["install-dialog"]

    println("Note: For a good posix shell experience in in Julia, you will need a MinGW installation. "*
            "If you already have Git installed, you already have a proper MinGW installation, and can skip installation. "*
            "If you are unsure, go ahead and mark MinGW for installation.")
    println()

    dict = Dict()
    for app in ["MinGW", "VSCode", "Juno", "Pluto", "PyCall", "Jupyter" ]
        doflagset(app, dict)
    end

    if !isfile("$juliawinpackages/curl/bin/curl.exe")
        install_curl()
    end

    if dict["MinGW"]
        install_git()
    end

    if dict["VSCode"]
        install_vscode()
    end

    if dict["Juno"]
        install_atom()
        install_juno()
    end

    if dict["Pluto"]
        install_pluto()
    end

    if dict["PyCall"]
        @eval Pkg.add("PyCall")
        activate_binary("python")
    end

    if dict["Jupyter"]
       install_jupyter()
    end
end

max_stars_repo_name: heetbeet/julia-win

max_issues_repo_name: heetbeet/julia-win

max_forks_repo_name: heetbeet/julia-win


text:
module DataUtils

include("countmap.jl")
include("functionutils.jl")

end # module DataUtils


text:
length(ARGS)>=4 || error("Usage: julia zn11nnet.jl datafile lr drop1 drop2 hidden1 hidden2...")
using HDF5,JLD,KUnet
KUnet.gpu(true)
argn = 0
isdefined(:xtrn) || (@date @load ARGS[argn+=1])

lr = eval(parse(ARGS[argn+=1]))
drop1 = eval(parse(ARGS[argn+=1]))
drop2 = eval(parse(ARGS[argn+=1]))

net = Layer[]
while argn < length(ARGS)
    d = isempty(net) ? drop1 : drop2
    d > 0 && (net = push!(net, Drop(d)))
    h = eval(parse(ARGS[argn+=1]))
    net = append!(net, [Mmul(h), Bias(), Relu()])
end
d = isempty(net) ? drop1 : drop2
d > 0 && (net = push!(net, Drop(d)))
h = size(ytrn, 1)
net = append!(net, [Mmul(h), Bias(), XentLoss()])
setparam!(net, adagrad=1e-8, lr=lr)
display(net)
println("")

for epoch=1:20
    @date train(net, xtrn, ytrn)
    #@date zdev = predict(net,xdev)
    @date ztst = predict(net,xtst)
    @show (epoch, accuracy(ytst,ztst)) # accuracy(ydev,zdev)
end

# for h=7:16; h0=2^h
#     for g=3:8; g0=2f0^-g
#         @show (h0, g0)
#         net = [Mmul(h0), Bias(), Relu(), Mmul(size(ytrn,1)), Bias(), XentLoss()]
#         setparam!(net, adagrad=1e-8, lr=g0)
#         @date train(net, xtrn, ytrn)
#         @date zdev = predict(net,xdev)
#         @date ztst = predict(net,xtst)
#         @show (h0, g0, accuracy(ydev,zdev), accuracy(ytst,ztst))
#     end
# end

# cnet = cpucopy(net)
# @time train(net, xtrn, ytrn; iters=1000)
# net = gpucopy(cnet)
# gc()
# @time @profile train(net, xtrn, ytrn; iters=1000)

# for epoch=1:10
#     @show epoch
#     @date train(net, xtrn, ytrn; iters=1000)
#     @date zdev = predict(net,xdev)
#     @date ztst = predict(net,xtst)
#     @show (epoch, accuracy(ydev,zdev), accuracy(ytst,ztst))
# end

# Use this to optimize gamma:
# isdefined(:xtrn) || (@date @load "zn11cpv.jld")
# niters=2000
# nbatch=100
# nc = size(ytrn,1)
# for g0=0.10:0.01:0.15
# #for g=2:.5:4; g0 = 2f0^-g
#     println("")
#     @show g0
#     net = Layer[KPerceptron(nc, KUnet.kgauss, [g0])]
#     @date train(net, xtrn, ytrn; iters=niters, batch=nbatch)
#     @show size(net[1].s)
#     @date y = predict(net,xdev)
#     @show accuracy(ydev,y)
# end


text:
using Plots

"""
    randomwalks1d(m,n,stepfunction)

Simulate `m` random walks of length `n`, and plot the positions against time.

Each random step is determined by `stepfunction`.

# Example
```julia-repl
julia> randomwalks1d(5, 100, m -> rand([-1,1],m))
```
"""
function randomwalks1d(m,n,stepfunction)
	walks = zeros(m,n+1)
	for step ∈ 1:n
		walks[:,step+1] = walks[:,step] + stepfunction(m)
	end

	return plot(
		0:n,
		[walks[walk,:] for walk ∈ 1:m],
		legend = false
	)
end

"""
    randomwalk2d(n)

Simulate a two dimensional simple symmetric random walk of length `n`, and plot the path.

# Example
```julia-repl
julia> randomwalk2d(100)
```
"""
function randomwalk2d(n)
	walk = fill((0,0),(n+1))
	for step ∈ 1:n
		walk[step+1] = walk[step] .+ rand([(-1,0),(1,0),(0,-1),(0,1)])
	end

	return plot(walk, legend = false, aspect_ratio = :equal)
end

"""
    randomwalk3d(n)

Simulate a three dimensional simple symmetric random walk of length `n`, and plot the path.

# Example
```julia-repl
julia> randomwalk3d(100)
```
"""
function randomwalk3d(n)
	walk = fill((0,0,0),(n+1))
	for step ∈ 1:n
		walk[step+1] = walk[step] .+
			rand([(-1,0,0),(1,0,0),(0,-1,0),(0,1,0),(0,0,-1),(0,0,1)])
	end

	return plot(walk, legend = false, aspect_ratio = :equal)
end

"""
    Yₙdistribution(samplesize, n)

Simulate `samplesize` one dimensional simple symmetric random walks of length `n`, and plot a histogram of their endpoints compared to a standard normal distribution.

# Example
```julia-repl
julia> Yₙdistribution(10000, 10000)
```
"""
function Yₙdistribution(samplesize, n)
	Yₙs = [sum(rand([-1,1],n)) for _ ∈ 1:samplesize]
	distribution = stephist(
        Yₙs ./ √n,
        bins = -3:0.2:3,
        norm = :pdf,
        xticks = -3:1:3,
        showaxis = :x,
		legend = false,
        yticks = false
    )
	return plot!(
        distribution,
        x -> exp(-(x^2)/2)/√(2π),
        c = :black
    )
end

"""
    randomdirection(d)

Pick a random direction in one of the positive or negative coordinate axes in `d` dimensions.

# Examples
```julia-repl
julia> randomdirection(3)
3-element Vector{Float64}:
  0.0
  0.0
 -1.0

julia> randomdirection(4)
4-element Vector{Float64}:
 0.0
 1.0
 0.0
 0.0
```
"""
function randomdirection(d)
	x = zeros(d)
	x[rand(1:d)] = rand([1,-1])
	return x
end

"""
    returnprobability(samplesize, nₘₐₓ, d)

Simulate `samplesize` `d` dimensional simple symmetric random walks for up to `nₘₐₓ` steps.

Find when (if ever) they return to the origin, and return a vector representing this.

# Example
```julia-repl
julia> returnprobability(1000, 100, 2)
51-element Vector{Float64}:
 0.0
 0.234
 [...]
 0.585
 0.585
```
"""
function returnprobability(samplesize, nₘₐₓ, d)
    returnfrequencies = zeros(nₘₐₓ ÷ 2)
    for _ ∈ 1:samplesize
        y = zeros(d)
        for halfn ∈ 1:(nₘₐₓ ÷ 2)
            y += randomdirection(d) + randomdirection(d)
            all(y .== 0) && (returnfrequencies[halfn:end] .+= 1; break)
        end
    end
    return vcat(0, returnfrequencies ./ samplesize)
end

"""
    plotreturnprobabilities(samplesize, nₘₐₓ, ds)

Plot a graph of `n` against the approximate proportion of random walks of dimension `d` which return to the origin in at most `n` steps, for `n` up to `nₘₐₓ` and `d` in `ds`.

# Example
```julia-repl
plotreturnprobabilities(1000, 10000, 1:4)
```
"""
function plotreturnprobabilities(samplesize, nₘₐₓ, ds)
	returnprobs = [returnprobability(samplesize, nₘₐₓ, d) for d ∈ ds]
    plot(
        0:2:nₘₐₓ,
        returnprobs,
		label = hcat(
			["d = $d, p ≈ $(returnprobs[i][end])" for (d,i) ∈ enumerate(ds)]...),
        ylims = (0,1)
    )
end

"""
    solvequadratic(a,b,c)

Solve the quadratic `ax² + bx + c = 0`, outputting the two solutions as a tuple.

# Examples
```julia-repl
julia> solvequadratic(1,2,1)
(-1.0 + 0.0im, -1.0 - 0.0im)

julia> solvequadratic(4,5,6)
(-0.625 + 1.0532687216470449im, -0.625 - 1.0532687216470449im)
```
"""
solvequadratic(a,b,c) = Tuple([(-b + ε * √complex(b^2 - 4a*c))/2a for ε ∈ [1,-1] ])

"""
    laplace◯(x, samplesize, r, ϕ, σ)

Approximate the value of `f(x)` for the harmonic function `f` given for `abs(z) == r` by `f(z) = ϕ(z)`, using `samplesize` Gaussian random walks with steps of standard deviation `σ`.

# Example
```julia-repl
julia> laplace◯(0.2+0.3im, 1000, 1, imag, 0.1)
0.28950478317730216
```
"""
function laplace◯(x, samplesize, r, ϕ, σ)
    ϕtotal = 0
    for _ ∈ 1:samplesize
        z = x
        while abs(z) < r
            z += σ*complex(randn(2)...)
        end
        # Shrink z towards x
        w = solvequadratic(conj(x) - conj(z), x*conj(z) - conj(x)*z, r^2 * (z - x))[2]
        ϕtotal += ϕ(w)
    end
    return ϕtotal / samplesize
end

"""
    plotlaplace◯(gridsize, samplesize, r, ϕ, σ)

Plot the results of `laplace◯` as evaluated on a square grid with `gridsize` graduations, restricted to the circle of radius `r` centred at the origin.

# Example
```julia-repl
julia> plotlaplace◯(20, 1000, 1, imag, 0.1)
```
"""
function plotlaplace◯(gridsize, samplesize, r, ϕ, σ)
    graduations = range(-r, r, length = gridsize)
    solutionplot = wireframe(
        graduations,
        graduations,
        (x,y) -> hypot(x,y) ≤ r ? laplace◯(x+y*im, samplesize, r, ϕ, σ) : NaN,
        c = :black,
        legend = false
    )

    θs = 0:0.01π:2.01π
    z(θ) = r*exp(im*θ)
    return plot!(solutionplot, r.*cos.(θs), r.*sin.(θs), ϕ.(z.(θs)), c = :black)
end

"""
    rectanglefunction(xₘᵢₙ, xₘₐₓ, yₘᵢₙ, yₘₐₓ, ϕxₘᵢₙ, ϕxₘₐₓ, ϕyₘᵢₙ, ϕyₘₐₓ)

Create a function defined on the border of the rectangle in the complex plane bounded by `xₘᵢₙ`, `xₘₐₓ`, `yₘᵢₙ`, `yₘₐₓ`.

`ϕxₘᵢₙ`, `ϕxₘₐₓ`, `ϕyₘᵢₙ`, `ϕyₘₐₓ` give the value of the function along the corresponding edges. At the corners, `ϕxₘᵢₙ` and `ϕxₘₐₓ` take precedence if there is discontinuity.

# Example
```julia-repl
julia> f = rectanglefunction(-1,1,-1,1,cos,cos,∘(-,cos),∘(-,cos));

julia> f(1+0.8im)
0.6967067093471654

julia> f(-0.4+1im)
-0.9210609940028851
```
"""
function rectanglefunction(xₘᵢₙ, xₘₐₓ, yₘᵢₙ, yₘₐₓ, ϕxₘᵢₙ, ϕxₘₐₓ, ϕyₘᵢₙ, ϕyₘₐₓ)
    return z -> real(z) == xₘᵢₙ ? ϕxₘᵢₙ(imag(z)) :
                real(z) == xₘₐₓ ? ϕxₘₐₓ(imag(z)) :
                imag(z) == yₘᵢₙ ? ϕyₘᵢₙ(real(z)) :
                imag(z) == yₘₐₓ ? ϕyₘₐₓ(real(z)) :
                NaN
end

"""
    laplace□(x, samplesize, boundary, ϕ, σ)

Approximate the value of `f(x)` for the harmonic function `f` given on the rectangle defined by `boundary` by `f(z) = ϕ(z)`, using `samplesize` Gaussian random walks with steps of standard deviation `σ`.

# Example
```julia-repl
julia> laplace□(0.1+0.1im, 1000, (-1,1,-1,1), abs, 0.1)
1.0822962717434006
```
"""
function laplace□(x, samplesize, boundary, ϕ, σ)
    xₘᵢₙ, xₘₐₓ, yₘᵢₙ, yₘₐₓ = boundary
    (real(x) ∈ (xₘᵢₙ,xₘₐₓ) || imag(x) ∈ (yₘᵢₙ,yₘₐₓ)) && return ϕ(x)
    ϕtotal = 0
    for _ ∈ 1:samplesize
        z = x
        while (xₘᵢₙ < real(z) < xₘₐₓ) && (yₘᵢₙ < imag(z) < yₘₐₓ)
            z += σ * complex(randn(2)...)
        end
        # Shrink z towards x
        tvals = [
            (xₘᵢₙ - real(x))/(real(z) - real(x)), (xₘₐₓ - real(x))/(real(z) - real(x)),
            (yₘᵢₙ - imag(x))/(imag(z) - imag(x)), (yₘₐₓ - imag(x))/(imag(z) - imag(x))
        ]
        t = minimum(tvals[tvals .≥ 0])
        w = round(x + t*(z-x), digits = 3)
        ϕtotal += ϕ(w)
    end
    return ϕtotal / samplesize
end

"""
    plotlaplace□(gridsize, samplesize, boundary, ϕ, σ)

Plot the results of `laplace□` as evaluated on a rectangular grid with `gridsize` graduations.

# Example
```julia-repl
julia> plotlaplace□(20, 1000, (-1,1,-1,1), abs, 0.1)
```
"""
function plotlaplace□(gridsize, samplesize, boundary, ϕ, σ)
    return wireframe(
        range(boundary[1], boundary[2], length = gridsize),
        range(boundary[3], boundary[4], length = gridsize),
        (x,y) -> laplace□(x+y*im, samplesize, boundary, ϕ, σ),
        c = :black,
        legend = false
    )
end


text:
export orbitdiagram

"""
    orbitdiagram(ds::DiscreteDynamicalSystem, i, p_index, pvalues; kwargs...)
Compute the orbit diagram (also called bifurcation diagram) of the given system,
saving the `i` variable(s) for parameter values `pvalues`. The `p_index` specifies
which parameter of the equations of motion is to be changed.

`i` can be `Int` or `AbstractVector{Int}`.
If `i` is `Int`, returns a vector of vectors. Else
it returns vectors of vectors of vectors.
Each entry are the points at each parameter value.

## Keyword Arguments
* `Ttr::Int = 1000` : Transient steps;
  each orbit is evolved for `Ttr` first before saving output.
* `n::Int = 100` : Amount of points to save for each initial condition.
* `dt = 1` : Stepping time. Changing this will give you the orbit diagram of
  the `dt` order map.
* `u0 = get_state(ds)` : Initial condition. Besides a vector you can also give
  a vector of vectors such that `length(u0) == length(pvalues)`. Then each parameter
  has a different initial condition.

See also [`poincaresos`](@ref) and [`produce_orbitdiagram`](@ref).
"""
function orbitdiagram(ds::DDS{IIP, S, D}, idxs, p_index, pvalues;
    n::Int = 100, Ttr::Int = 1000, u0 = get_state(ds), dt = 1) where {IIP, S, D}

    p0 = ds.p[p_index]
    if D == 1
        idxs != 1 &&
        error("You have a 1D system and yet you gave `i=$i`. What's up with that!?")
    end

    typeof(u0) <: Vector{<:AbstractVector} && @assert length(u0)==length(p)

    i = typeof(idxs) <: Int ? idxs : SVector{length(idxs), Int}(idxs...)

    output = _initialize_output(ds.u0, i, n, length(pvalues))
    integ = integrator(ds)
    _fill_orbitdiagram!(output, integ, i, pvalues, p_index, n, Ttr, u0, dt)
    ds.p[p_index] = p0
    return output
end

function _initialize_output(u::S, i::Int, n, l) where {S}
    output = [zeros(eltype(S), n) for k in 1:l]
end
function _initialize_output(u::S, i::SVector, n, l) where {S}
    s = u[i]
    output = [Vector{typeof(s)}(undef, n) for k in 1:l]
end

function _fill_orbitdiagram!(output, integ, i, pvalues, p_index,
    n, Ttr, u0, dt)

    for (j, p) in enumerate(pvalues)

        integ.p[p_index] = p

        if typeof(u0) <: Vector{<:AbstractVector}
            st = u0[j]
        else
            st = u0
        end

        reinit!(integ, st)
        step!(integ, Ttr)

        for k in 1:n
            step!(integ, dt)
            @inbounds output[j][k] = integ.u[i]
        end
    end
    return output
end


text:
# using Pkg
# Pkg.activate("c:/git/CuCountmap")
using CuCountMap
using SortingAlgorithms
using Test

@testset "CuCountMap.jl - countmap" begin
    v = rand(UInt16, 1_000_000);
    @time a = cucountmap(v);
    @time aa = countmap(v);
    @test a == aa
end

# @testset "CuCountMap.jl - sort" begin
#     v = rand(UInt64, 1_000_000);
#     @time sorted_v = gpuradixsort(v);
#     @time sorted_v_cpu = sort(v, alg=RadixSort);
#     @test sorted_v == sorted_v_cpu
# end


text:
"""
    Rayleigh(σ)

The *Rayleigh distribution* with scale `σ` has probability density function

```math
f(x; \\sigma) = \\frac{x}{\\sigma^2} e^{-\\frac{x^2}{2 \\sigma^2}}, \\quad x > 0
```

It is related to the [`Normal`](@ref) distribution via the property that if ``X, Y \\sim \\operatorname{Normal}(0,\\sigma)``, independently, then
``\\sqrt{X^2 + Y^2} \\sim \\operatorname{Rayleigh}(\\sigma)``.

```julia
Rayleigh()       # Rayleigh distribution with unit scale, i.e. Rayleigh(1)
Rayleigh(σ)      # Rayleigh distribution with scale σ

params(d)        # Get the parameters, i.e. (σ,)
scale(d)         # Get the scale parameter, i.e. σ
```

External links

* [Rayleigh distribution on Wikipedia](http://en.wikipedia.org/wiki/Rayleigh_distribution)

"""
struct Rayleigh{T<:Real} <: ContinuousUnivariateDistribution
    σ::T
    Rayleigh{T}(σ::T) where {T<:Real} = new{T}(σ)
end

function Rayleigh(σ::T; check_args=true) where {T <: Real}
    check_args && @check_args(Rayleigh, σ > zero(σ))
    return Rayleigh{T}(σ)
end

Rayleigh(σ::Integer) = Rayleigh(float(σ))
Rayleigh() = Rayleigh(1.0, check_args=false)

@distr_support Rayleigh 0.0 Inf

#### Conversions

convert(::Type{Rayleigh{T}}, σ::S) where {T <: Real, S <: Real} = Rayleigh(T(σ))
convert(::Type{Rayleigh{T}}, d::Rayleigh{S}) where {T <: Real, S <: Real} = Rayleigh(T(d.σ), check_args=false)

#### Parameters

scale(d::Rayleigh) = d.σ
params(d::Rayleigh) = (d.σ,)
partype(::Rayleigh{T}) where {T<:Real} = T

#### Statistics

mean(d::Rayleigh) = sqrthalfπ * d.σ
median(d::Rayleigh{T}) where {T<:Real} = sqrt2 * sqrt(T(logtwo)) * d.σ # sqrt(log(4))
mode(d::Rayleigh) = d.σ

var(d::Rayleigh{T}) where {T<:Real} = (2 - T(π)/2) * d.σ^2
std(d::Rayleigh{T}) where {T<:Real} = sqrt(2 - T(π)/2) * d.σ

skewness(d::Rayleigh{T}) where {T<:Real} = 2 * sqrtπ * (T(π) - 3)/(4 - T(π))^(3/2)
kurtosis(d::Rayleigh{T}) where {T<:Real} = -(6*T(π)^2 - 24*T(π) +16)/(4 - T(π))^2

entropy(d::Rayleigh{T}) where {T<:Real} = 1 - T(logtwo)/2 + T(MathConstants.γ)/2 + log(d.σ)

#### Evaluation

function pdf(d::Rayleigh{T}, x::Real) where T<:Real
    σ2 = d.σ^2
    x > 0 ? (x / σ2) * exp(- (x^2) / (2σ2)) : zero(T)
end

function logpdf(d::Rayleigh{T}, x::Real) where T<:Real
    σ2 = d.σ^2
    x > 0 ? log(x / σ2) - (x^2) / (2σ2) : -T(Inf)
end

function logccdf(d::Rayleigh, x::Real)
    z = - x^2 / (2 * d.σ^2)
    return x > 0 ? z : isnan(x) ? oftype(z, NaN) : zero(x)
end
ccdf(d::Rayleigh, x::Real) = exp(logccdf(d, x))

cdf(d::Rayleigh, x::Real) = -expm1(logccdf(d, x))
logcdf(d::Rayleigh, x::Real) = log1mexp(logccdf(d, x))

quantile(d::Rayleigh, p::Real) = sqrt(-2d.σ^2 * log1p(-p))

#### Sampling

rand(rng::AbstractRNG, d::Rayleigh) = d.σ * sqrt(2 * randexp(rng))

#### Fitting

function fit_mle(::Type{<:Rayleigh}, x::AbstractArray{T}) where {T<:Real}
    # Compute MLE (and unbiasd estimator) of σ^2
    s2 = zero(T)
    for xi in x
        s2 += xi^2
    end

    s2 /= (2*length(x))
    return Rayleigh(sqrt(s2))
end


text:
macro newtontrace()
    quote
        if tracing
            dt = Dict()
            if o.extended_trace
                dt["x"] = copy(x)
                dt["g(x)"] = copy(gr)
                dt["h(x)"] = copy(H)
            end
            grnorm = vecnorm(gr, Inf)
            update!(tr,
                    iteration,
                    f_x,
                    grnorm,
                    dt,
                    o.store_trace,
                    o.show_trace,
                    o.show_every,
                    o.callback)
        end
    end
end

immutable Newton <: Optimizer
    linesearch!::Function
end

Newton(; linesearch!::Function = hz_linesearch!) =
  Newton(linesearch!)

function optimize{T}(d::TwiceDifferentiableFunction,
                     initial_x::Vector{T},
                     mo::Newton,
                     o::OptimizationOptions)
    # Print header if show_trace is set
    print_header(o)

    # Maintain current state in x and previous state in x_previous
    x, x_previous = copy(initial_x), copy(initial_x)

    # Count the total number of iterations
    iteration = 0

    # Track calls to function and gradient
    f_calls, g_calls = 0, 0

    # Count number of parameters
    n = length(x)

    # Maintain current gradient in gr
    gr = Array(T, n)

    # The current search direction
    # TODO: Try to avoid re-allocating s
    s = Array(T, n)

    # Buffers for use in line search
    x_ls, g_ls = Array(T, n), Array(T, n)

    # Store f(x) in f_x
    f_x_previous, f_x = NaN, d.fg!(x, gr)
    f_calls, g_calls = f_calls + 1, g_calls + 1

    # Store h(x) in H
    H = Array(T, n, n)
    d.h!(x, H)

    # Keep track of step-sizes
    alpha = alphainit(one(T), x, gr, f_x)

    # TODO: How should this flag be set?
    mayterminate = false

    # Maintain a cache for line search results
    lsr = LineSearchResults(T)

    # Trace the history of states visited
    tr = OptimizationTrace{typeof(mo)}()
    tracing = o.store_trace || o.show_trace || o.extended_trace || o.callback != nothing
    @newtontrace

    # Assess multiple types of convergence
    x_converged, f_converged = false, false
    g_converged = vecnorm(gr, Inf) < o.g_tol

    # Iterate until convergence
    converged = g_converged

    while !converged && iteration < o.iterations
        # Increment the number of steps we've had to perform
        iteration += 1

        # Search direction is always the negative gradient divided by
        # a matrix encoding the absolute values of the curvatures
        # represented by H. It deviates from the usual "add a scaled
        # identity matrix" version of the modified Newton method. More
        # information can be found in the discussion at issue #153.
        F, Hd = ldltfact!(Positive, H)
        s[:] = -(F\gr)

        # Refresh the line search cache
        dphi0 = vecdot(gr, s)
        clear!(lsr)
        push!(lsr, zero(T), f_x, dphi0)

        # Determine the distance of movement along the search line
        alpha, f_update, g_update =
          mo.linesearch!(d, x, s, x_ls, g_ls, lsr, alpha, mayterminate)
        f_calls, g_calls = f_calls + f_update, g_calls + g_update

        # Maintain a record of previous position
        copy!(x_previous, x)

        # Update current position # x = x + alpha * s
        LinAlg.axpy!(alpha, s, x)

        # Update the function value and gradient
        f_x_previous, f_x = f_x, d.fg!(x, gr)
        f_calls, g_calls = f_calls + 1, g_calls + 1

        # Update the Hessian
        d.h!(x, H)

        x_converged,
        f_converged,
        g_converged,
        converged = assess_convergence(x,
                                       x_previous,
                                       f_x,
                                       f_x_previous,
                                       gr,
                                       o.x_tol,
                                       o.f_tol,
                                       o.g_tol)

        @newtontrace
    end

    return MultivariateOptimizationResults("Newton's Method",
                                           initial_x,
                                           x,
                                           Float64(f_x),
                                           iteration,
                                           iteration == o.iterations,
                                           x_converged,
                                           o.x_tol,
                                           f_converged,
                                           o.f_tol,
                                           g_converged,
                                           o.g_tol,
                                           tr,
                                           f_calls,
                                           g_calls)
end


text:
###############################################################
# SymTensorValue Type
###############################################################

"""
Type representing a symmetric second-order tensor
"""
struct SymTensorValue{D,T,L} <: MultiValue{Tuple{D,D},T,2,L}
    data::NTuple{L,T}
    function SymTensorValue{D,T}(data::NTuple{L,T}) where {D,T,L}
        @assert L == D*(D+1)/2
        new{D,T,L}(data)
    end
end

###############################################################
# Constructors (SymTensorValue)
###############################################################

# Empty SymTensorValue constructor

SymTensorValue()                   = SymTensorValue{0,Int}(NTuple{0,Int}())
SymTensorValue{0}() where {T}      = SymTensorValue{0,Int}(NTuple{0,Int}())
SymTensorValue{0,T}() where {T}    = SymTensorValue{0,T}(NTuple{0,T}())
SymTensorValue(data::NTuple{0})    = SymTensorValue{0,Int}(data)
SymTensorValue{0}(data::NTuple{0}) = SymTensorValue{0,Int}(data)

# SymTensorValue single NTuple argument constructor

@generated function SymTensorValue(data::NTuple{L,T}) where {L,T}
  msg = "Invalid number of scalar arguments in SymTensorValue constructor"
  V = (sqrt(1+8*L)-1)/2
  @assert floor(Int,V) == ceil(Int,V) msg
  D = Int(V)
  quote
    SymTensorValue{$D,T}(data)
  end
end
SymTensorValue{D}(data::NTuple{L,T}) where {D,L,T}           = SymTensorValue{D,T}(data)
SymTensorValue{D,T1}(data::NTuple{L,T2}) where {D,L,T1,T2}   = SymTensorValue{D,T1}(NTuple{L,T1}(data))
SymTensorValue{D,T1,L}(data::NTuple{L,T2}) where {D,L,T1,T2} = SymTensorValue{D,T1}(NTuple{L,T1}(data))

# SymTensorValue single Tuple argument constructor

SymTensorValue(data::Tuple)                      = SymTensorValue(promote(data...))
SymTensorValue{D}(data::Tuple) where {D}         = SymTensorValue{D}(promote(data...))
SymTensorValue{D,T1}(data::Tuple) where {D,T1}   = SymTensorValue{D,T1}(NTuple{length(data),T1}(data))
SymTensorValue{D,T1,L}(data::Tuple) where {D,T1,L}   = SymTensorValue{D,T1}(NTuple{L,T1}(data))

# SymTensorValue Vararg constructor

SymTensorValue(data...)                    = SymTensorValue(data)
SymTensorValue{D}(data...)    where {D}    = SymTensorValue{D}(data)
SymTensorValue{D,T1}(data...) where {D,T1} = SymTensorValue{D,T1}(data)
SymTensorValue{D,T1,L}(data...) where {D,T1,L} = SymTensorValue{D,T1}(data)

# SymTensorValue single AbstractMatrix argument constructor

#From Square Matrices
@generated function _flatten_upper_triangle(data::AbstractArray,::Val{D}) where D
  str = ""
  for i in 1:D
    for j in i:D
      str *= "data[i,j], "
    end
  end
  Meta.parse("($str)")
end

SymTensorValue(data::AbstractMatrix{T}) where {T} = ((D1,D2)=size(data); SymTensorValue{D1}(data))
SymTensorValue{D}(data::AbstractMatrix{T}) where {D,T} = SymTensorValue{D,T}(_flatten_upper_triangle(data,Val{D}()))
SymTensorValue{D,T1}(data::AbstractMatrix{T2}) where {D,T1,T2} = SymTensorValue{D,T1}(_flatten_upper_triangle(data,Val{D}()))
SymTensorValue{D,T1,L}(data::AbstractMatrix{T2}) where {D,T1,T2,L} = SymTensorValue{D,T1,L}(_flatten_upper_triangle(data,Val{D}()))

###############################################################
# Conversions (SymTensorValue)
###############################################################

@generated function _SymTensorValue_to_array(arg::SymTensorValue{D,T,L}) where {D,T,L}
  str = ""
  for j in 1:D
    for i in 1:D
      p = _2d_sym_tensor_linear_index(D,i,j)
      str *= "arg.data[$p], "
    end
  end
  Meta.parse("SMatrix{D,D,T}(($str))")
end

# Direct conversion
convert(::Type{<:SymTensorValue{D,T}}, arg::AbstractArray) where {D,T} = SymTensorValue{D,T}(arg)
convert(::Type{<:SymTensorValue{D,T}}, arg::Tuple) where {D,T} = SymTensorValue{D,T}(arg)

# Inverse conversion
convert(::Type{<:MMatrix{D,D,T}}, arg::SymTensorValue) where {D,T} = MMatrix{D,D,T}(_SymTensorValue_to_array(arg))
convert(::Type{<:SMatrix{D,D,T}}, arg::SymTensorValue) where {D,T} = _SymTensorValue_to_array(arg)
convert(::Type{<:NTuple{L,T}}, arg::SymTensorValue) where {L,T} = NTuple{L,T}(Tuple(arg))

# Internal conversion
convert(::Type{<:SymTensorValue{D,T}}, arg::SymTensorValue{D}) where {D,T} = SymTensorValue{D,T}(Tuple(arg))
convert(::Type{<:SymTensorValue{D,T}}, arg::SymTensorValue{D,T}) where {D,T} = arg

###############################################################
# Other constructors and conversions (SymTensorValue)
###############################################################

@generated function zero(::Type{<:SymTensorValue{D,T}}) where {D,T}
  L=Int(D*(D+1)/2)
  quote
    SymTensorValue{D,T}(tfill(zero(T),Val{$L}()))
  end
end

zero(::Type{<:SymTensorValue{D,T,L}}) where {D,T,L} = SymTensorValue{D,T}(tfill(zero(T),Val{L}()))
zero(::SymTensorValue{D,T,L}) where {D,T,L} = zero(SymTensorValue{D,T,L})

@generated function one(::Type{<:SymTensorValue{D,T}}) where {D,T}
  str = join(["$i==$j ? one(T) : zero(T), " for i in 1:D for j in i:D])
  Meta.parse("SymTensorValue{D,T}(($str))")
end
one(::SymTensorValue{D,T}) where {D,T} = one(SymTensorValue{D,T})

mutable(::Type{<:SymTensorValue{D,T}}) where {D,T} = MMatrix{D,D,T}
mutable(::SymTensorValue{D,T}) where {D,T} = mutable(SymTensorValue{D,T})

change_eltype(::Type{SymTensorValue{D,T1,L}},::Type{T2}) where {D,T1,T2,L} = SymTensorValue{D,T2,L}
change_eltype(::SymTensorValue{D,T1,L},::Type{T2}) where {D,T1,T2,L} = change_eltype(SymTensorValue{D,T1,L},T2)

get_array(arg::SymTensorValue{D,T,L}) where {D,T,L} = convert(SMatrix{D,D,T}, arg)

###############################################################
# Introspection (SymTensorValue)
###############################################################

eltype(::Type{<:SymTensorValue{D,T}}) where {D,T} = T
eltype(::SymTensorValue{D,T}) where {D,T} = eltype(SymTensorValue{D,T})

size(::Type{<:SymTensorValue{D}}) where {D} = (D,D)
size(::SymTensorValue{D}) where {D} = size(SymTensorValue{D})

length(::Type{<:SymTensorValue{D}}) where {D} = D*D
length(::SymTensorValue{D}) where {D} = length(SymTensorValue{D})

num_components(::Type{<:SymTensorValue{D}}) where {D} = length(SymTensorValue{D})
num_components(::SymTensorValue{D}) where {D} = num_components(SymTensorValue{D})


text:
using MCMCChains, Test
using DataFrames

@testset "DataFrame constructor tests" begin
    val = rand(1000, 8, 4)
    chn = Chains(val,
                ["a", "b", "c", "d", "e", "f", "g", "h"],
                Dict(:internals => ["c", "d", "e", "f", "g", "h"]))

    df = DataFrame(chn, showall=true)
    @test size(df) == (4000, 8)

    df1 = DataFrame(chn, [:parameters])
    @test size(df1) == (4000, 2)

    df2 = DataFrame(chn, [:internals, :parameters])
    @test size(df2) == (4000, 8)

    df3 = DataFrame(chn[:a])
    @test size(df3) == (4000, 1)

    df4 = DataFrame(chn[:b])
    @test size(df4) == (4000, 1)

    df5 = DataFrame(chn, [:parameters], append_chains=false)
    @test size(df5) == (4, )
    @test size(df5[1]) == (1000, 2)

    df6 = DataFrame(chn, [:parameters, :internals], append_chains=false)
    @test size(df6) == (4, )
    @test size(df6[1]) == (1000, 8)

    df7 = DataFrame(chn, [:parameters, :internals],
        remove_missing_union=false)
    @test size(df7) == (4000, 8)

    df8 = DataFrame(chn, [:parameters, :internals],
        remove_missing_union=false, append_chains=false)
    @test size(df8) == (4, )
    @test size(df8[1]) == (1000, 8)
end
@testset "ChainDataFrame converter" begin
    val = rand(1000, 5, 4)
    chain = Chains(val)

    cdfs = describe(chain)

    dfs = DataFrame(cdfs)
    df1 = DataFrame(cdfs[1])
    df2 = DataFrame(cdfs[2])

    @test dfs[1] == df1
    @test dfs[2] == df2
end


text:
module Transliterate

export transliterate

include("replacements.jl")
"All languages defined in (`replacements`)[@ref]."
const languages = first.(replacements)

"""
    transliterate(str; languages=languages, custom_replacements=())

Converts non-ASCII characters into ASCII using [transliteration](https://en.wikipedia.org/wiki/Transliteration).

# Examples
```julia-repl
julia> transliterate("Déjà Vu!")
"Deja Vu!"

julia> transliterate("Привет")
"Privet"

julia> transliterate("ث س و"; languages=["ar"])
"th s w"

julia> transliterate("≠ ∉"; custom_replacements=("≠" => "not equal", "∉" => "not in"))
"not equal not in"
```
"""
function transliterate(str; languages=languages, custom_replacements=())
    str = foldl(replace, custom_replacements, init = str)

    if "la" ∉ languages
        str = foldl(replace, replacements[1].second, init = str)
    end
    for language ∈ languages
        language_id = findfirst(x -> x.first == language, replacements)
        d = replacements[language_id].second
        str = foldl(replace, d, init = str)
    end

    return str
end

end # module


text:
# module fiberDataProcLib
# export  CalDatainFile
using Statistics
using Printf
using Dates

function getRelPos(CumSumDataobj::spaceCumSumData ,vPos::Int )
    return vPos- CumSumDataobj.startID +1
end
function getAbsPos(CumSumDataobj::spaceCumSumData ,vPos::Int )
    return vPos- CumSumDataobj.startID +1

end
function getSpaceLoopRng(spaceCumSumDataObj::spaceCumSumData)
    # println("RNG:  $(spaceCumSumDataObj.startID) : $(spaceCumSumDataObj.endID) " )
    return spaceCumSumDataObj.startID : spaceCumSumDataObj.endID
end

function InitObj(spaceCumSumDataObj::spaceCumSumData , AlrmProcStartID1::Int,AlrmProcEndtID1::Int , startID1::Int ,endID1::Int,maxTimeCnt1::Int= 100 )
    spaceCumSumDataObj.startID= startID1
    spaceCumSumDataObj.endID =endID1
    SegdataNum = endID1-startID1+1

    spaceCumSumDataObj.AlrmProcStartID  = AlrmProcStartID1
    spaceCumSumDataObj.AlrmProcEndtID = AlrmProcEndtID1

    spaceCumSumDataObj.maxTimeCnt = maxTimeCnt1

    spaceCumSumDataObj.defaultSpaceCnt =10
    # Cache Data Number ：
    minutes = 10.0;
    spaceCumSumDataObj.defaultTimeCnt = floor(UInt ,minutes*60.0/0.42) #

    spaceCumSumDataObj.m1 = Array{UInt32}(undef,spaceCumSumDataObj.defaultTimeCnt,SegdataNum)
    spaceCumSumDataObj.m1 .= 0
    spaceCumSumDataObj.rectSum = Array{UInt64}(undef,SegdataNum)
    spaceCumSumDataObj.curRow = 0
    spaceCumSumDataObj.dt= ""
    spaceCumSumDataObj.isInit = false

end

########################
#计算当前的行号
###########################
function getCurRowCnt(  temobj::spaceCumSumData)
    idx = temobj.curRow +1
    dataCnt =  length( @view temobj.m1[:,1] )
    temobj.curRow = idx %dataCnt
    # @printf( "====%d === %d==============\n" , idx ,dataCnt  )
    return idx
 end

###########  计算做累加计算需要的 矩阵 的时间累加时  的  下标，
###因为 矩阵的数据结构为 循环填写有效数据，所以一共有两对   时间下标，每对包含起点、终点
####################
function getTimeRng(tcMax,curtc,vTcNum ) # time is Row
    tcNum = Int(vTcNum )
        if curtc < 1 ;curtc = tcMax ;end
        if curtc > tcMax ;curtc = tcMax ;end
        if tcNum > tcMax; tcNum =tcMax;end

        tctop1=0;tcbottom1 =0
        tctop2=0;tcbottom2 = 0

        tctop1= curtc
        tcbottom1  = tctop1-  tcNum + 1
        # @info ":::::::-::::::: ",tcbottom1 , tctop1,  tcNum , typeof(tcNum)
        if tcbottom1  <  1 ;
            tcbottom1 =1;
            cnt =  ( tctop1- tcbottom1)+1
            if  cnt < tcNum;
                tctop2= tcMax
                tcbottom2 = tctop2- (tcNum-cnt-1)
            end
            # @printf( "%d,%d,%d----,%d,%d,%d---%d\n",  tctop2,tcbottom2,(tc-cnt),  tctop1,tcbottom1, cnt  ,tc )
        end

      return (  (  tctop1,tcbottom1), (tctop2,tcbottom2,) )
end

#############################
# 根据光缆空间位置的 中心位置 vspacePos ,以及计算累加的总空间点数vsc,计算出  累加计算的两个下标位置。
# 当 vsc为偶数时，调整为单数。
##############

function getSpaceRng( cumsumObj::spaceCumSumData, vspacePos,vsc)
        vsideMax =cumsumObj.endID
        vsideMin = cumsumObj.startID

        if vspacePos <  vsideMin || vspacePos  >  vsideMax ;
            error("vspacePos 越界  ")
        end

        if vsc % 2 ==0
            vsc+=1
        end
        left = vsideMin ;
        right = vsideMax;
        halfsc = trunc(Int,vsc/2)

        left = vspacePos-halfsc-1 ;
        right =vspacePos+halfsc;
        if left < vsideMin ;
            left=0;
        end
        if  right >vsideMax ;
            right =vsideMax;
        end
       return  (left ,right)
end

function getRelativeSpaceRng( cumsumObj::spaceCumSumData, vposIndex,vsc)
    vsideMax =cumsumObj.endID
    vsideMin = cumsumObj.startID

    if vposIndex <  0 || vposIndex  >  vsideMax ;
        error("vspacePos 越界  ")
    end

    if vsc % 2 ==0
        vsc+=1
    end
    left = 0 ;
    right = vsideMax;
    halfsc = trunc(Int,vsc/2)

    left = vposIndex-halfsc-1 ;
    right =vposIndex+halfsc;
    if left < 0 ;
        left=0;
    end
    if  right >vsideMax  ;
        right =vsideMax;
    end
   return (left ,right)
end

function calRectAvg( CumSumDataobj::spaceCumSumData, vPos ,vsc=1 ,vTcNum=20 )
    temobj = CumSumDataobj
    rectsumObj =temobj.rectSum

    vsideMax  = length( rectsumObj )  #？

    timeRngMax1 = length(　temobj.m1[:,1])　 #？

    curtcIdx = temobj.curRow
    # (left,right )=getSpaceRng( vPos, vsideMax,vsc)
    (absleft1,absright1 )=getSpaceRng( temobj,vPos,vsc)
    # @info "-??-，vTcNum,curtcIdx ,timeRngMax1", vTcNum,curtcIdx ,timeRngMax1
    a = getTimeRng(timeRngMax1,curtcIdx,vTcNum )
    # @info "getTimeRng " , a
    # println(  "time rng : $a " )
    right = getRelPos(CumSumDataobj ,Int(absright1 )  )
    if absleft1 != 0;
        left = getRelPos(CumSumDataobj ,absleft1 )
    else
        left=0
    end

    numbercnt =0
    rows=0;
    cols=0
    if   left == 0 ;
        cols =  right
    else
        cols = right-left
    end

    sumv =0
    for ( top1,bottom1) in a;  ## time range loop
        # println( "absleft1,absright1,left,right,bottom1,top1::$absleft1,$absright1,$left,$right, $bottom1, $top1 ")
            # println( "VPos,absleft1,absright1,left,right,bottom1,top1::$vPos  :$absleft1,$absright1,$left,$right, $bottom1, $top1 ")
            if top1 == 0;
                continue;
            end

            rows += top1 - bottom1+1
            # println(  " ###:",size( temobj.m1) ," " ,right," " ,bottom1," ",top1)
            v1 = @view temobj.m1[ bottom1:top1,right  ]
            sum1= sum( v1)
            # println( sum1 )

            sum0=0 ;
            if left != 0;
                v0 =@view temobj.m1[bottom1:top1,left  ]
                sum0 = sum(  v0)
            end

            sumv += sum1-sum0

    end   # for ( top1,bottom1) in a;

        numbercnt =rows*cols
        if vTcNum*vsc != numbercnt
            # error(  "Logic Error vTcNum != numbercnt " )
            # println( "Logic Error vTcNum*vsc != numbercnt :$vTcNum * $vsc  != $numbercnt" )
        end
        avg = round( UInt,sumv /numbercnt)
        return avg
        # curtcIdx = getRelPos(CumSumDataobj ,vPos )
        # rectsumObj[ curtcIdx ]  = avg

end

######################
# vsc  计算均值时的 空间点颗粒 的 个数
# vTcNum  计算均值时的 时间点的 个数
# vPos  当前计算点的空间 位置，空间索引点 ，对应本结构下的 rectSum的下标位置,不是光缆的绝对参考 颗粒号
###############
function calRectSum( CumSumDataobj::spaceCumSumData, vPos ,vsc ,vTcNum )
    temobj = CumSumDataobj
    rectsumObj = temobj.rectSum
    avg = calRectAvg( CumSumDataobj, vPos ,vsc ,vTcNum)
    curtcIdx = getRelPos(CumSumDataobj ,vPos )
    # println("$curtcIdx:$avg")
    rectsumObj[ curtcIdx ]  = avg
end

function AlarmCal( CumSumDataobj::spaceCumSumData,threshold= 4000)
    println("Nothing to do  in AlarmCal")
    return
end

########### 计算空间距离上的累加数据，并保存到对应的数组中。
function cumSumProc(spaceCumSumDataobj::spaceCumSumData ,featuredata::Array{UInt16,1} )
    temobj = spaceCumSumDataobj
    idx = getCurRowCnt(  temobj)
    va = @view temobj.m1[idx , : ]  #
    pieceData = @view featuredata[ temobj.startID : temobj.endID ]
    # println( "length(pieceData) ,length(va ) ，length(featuredata）: $(length(pieceData)) ,$(length(va )),$(length(featuredata) )" )
    va .=   pieceData
    va .= cumsum(  va )

end

function setDTValue(cumObj::spaceCumSumData,dt  )
    cumObj.dt =  dt
end

function CalDynmicTh(CumSumDataobj::spaceCumSumData, vPos)
    segID = getAbsPos(CumSumDataobj ,vPos)
    CalTimePeriod,calCFactor ,minTh,maxTh = getSegDynParaValue( glbSegAlarmParasDictObj,segID )
    #  println("----",calCFactor )
    vsc=1
    avg = calRectAvg(CumSumDataobj, vPos ,vsc,CalTimePeriod  )
    if calCFactor != typemax( Int )
        dynTh1 = floor(Int,avg * (calCFactor/100) )
        if dynTh1 < minTh
            dynTh1 = minTh
        elseif dynTh1 > maxTh
            dynTh1 = maxTh
        end
        curAlarmsinSEG[segID].DynTh    = dynTh1
    else
        curAlarmsinSEG[segID].DynTh    = typemax( Int )
    end
    # @info " SefID ,Dyn , AVg  CalTimePeriod,calCFactor ,minTh,maxTh .is :  ",segID,curAlarmsinSEG[segID].DynTh ,Int(avg),CalTimePeriod,calCFactor ,minTh,maxTh
end

function AlarmRTCalAna( CumSumDataobj::spaceCumSumData)
    temobj = CumSumDataobj
    rectsumObj =temobj.rectSum
    viewa  = @view rectsumObj[ 1 : end]
    e = length( viewa )
     @threads for k = 1 : e
        temobj.startID > 1 ? segID = k - 1 + temobj.startID : segID = k
        if Int( viewa[ k ] ) > curAlarmsinSEG[segID].DynTh
            # @info "AlarmRTCalAna :",segID , Int(viewa[ k ] ), curAlarmsinSEG[segID].DynTh
            AlarmDynStatus[segID ] = true
        else
            AlarmDynStatus[segID ] = false
        end
    end
    # showResult( "$(temobj.dt[1:8]):$(temobj.dt[9:14])" ,AlarmDynStatus, AlarmStaticStatus,CumSumDataobj )
end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
workspace()
using HttpServer
using WebSockets
using JSON
#using JLD
#using Jukai
#using Jukai.Tokenization
#using Jukai.Tagging
#using Merlin

wsh = WebSocketHandler() do req, client
    println("Client: $(client.id) is connected.")
    #write(client, conf)
    while true
        #println("Request from $(client.id) recieved.")
        msg = bytestring(read(client))
        length(msg) > 1000 && continue
        #chars = convert(Vector{Char}, msg)
        #doc = Tokenization.decode(tokenizer, chars)
        #POSTagging.decode!(postagger, doc)
        #doc = map(x -> map(to_dict, x), doc)
        #res = JSON.json(doc)
        #println(res)
        #write(client, res)
    end
end

onepage = readall(joinpath(dirname(@__FILE__), "index.html"))
httph = HttpHandler() do req::Request, res::Response
    Response(onepage)
end
server = Server(httph, wsh)
run(server, 3000)

#=
const sample_tokens = [
Token(1, 6, 1, "Pierre", 1, "NNP", -2, 2),
Token(2, "Vinken", 2, "NNP", -1, 8),
Token(3, ",", 3, ",", -2, 2),
Token(4, "61", 4, "CD", -4, 5),
Token(5, "years", 5, "NNS", -3, 6),
Token(6, "old", 6, "JJ", -2, 2),
Token(7, ",", 7, ",", -2, 2),
Token(8, "will", 8, "MD", 0, 0),
Token(9, "join", 9, "VB", -1, 8),
Token(10, "the", 10, "DT", -3, 11),
Token(11, "board", 11, "NN", -2, 9),
Token(12, "Nov.", 12, "NNP", -2, 9),
Token(13, "29", 13, "CD", -3, 12),
Token(14, ".", 13, ".", -1, 8)]
=#

max_stars_repo_name: hshindo/JukaiWriter

max_issues_repo_name: hshindo/JukaiWriter

max_forks_repo_name: hshindo/JukaiWriter


text:
using Klimakoffer

NT = 48

mesh = Mesh()

co2_concentration = 315.0

model = Model(mesh, NT; co2_concentration = co2_concentration,compute_albedo=true, compute_sea_ice_extent = true, sea_ice_extent_year = 1979)

discretization = Discretization(mesh, model, NT)

GlobTemp = compute_equilibrium!(discretization, update_heat_capacity = true ,update_solar_forcing = true)


text:
using MaxwellSALT
using StaticArrays
using Test

include("gainobj.jl")


text:
using GemmKernels
using Test

macro test_if(label, expr)
    return quote
        if isempty(ARGS) || $(label) in ARGS
            $(esc(expr))
        else
            nothing
        end
    end
end

@testset "GemmKernels.jl" begin
    include("tiling.jl")
    include("matmul.jl")
    include("blas.jl")
end


text:
fname = :ARCMA97_abs
c = Combi(hessian_dense,PDataMA97,solve_modelARCDiagAbs,preprocessMA97,decreaseFact,Tparam())
include("Template.jl")


text:
# TODO: Check how other classification models handle the target type. Should I
#       simply enforce it to be an integer? Furthermore, a may consider enforcing
#       targets to be in the range [1, nclasses]. These could be label encoded
#       by the user or the model could detect and perform the encoding. If I did
#       this, `discriminators` could simply be a list, with the discriminator's
#       indices being the target they are associated to. `predict_response` return
#       would also be simplified to a vector or matrix
mutable struct MultiDiscriminatorClassifier{C, N <: Discriminator} <: AbstractModel
    n::Int
    seed::Union{Int,Nothing}
    partitioner::Union{RandomPartitioner,Nothing}
    discriminators::Dict{C,N}

    function MultiDiscriminatorClassifier{C,N}(n::Int; seed::Union{Nothing,Int}=nothing) where {C,N<:Discriminator}
        n < 1 && throw(DomainError(n, "tuple size `n` may not be less then 1"))
        !isnothing(seed) && seed < 0 && throw(DomainError(seed, "`seed` must be non-negative"))

        new{C,N}(n, seed, nothing, Dict{C,N}())
    end

    function MultiDiscriminatorClassifier{C,N}(width::Int, n::Int; seed::Union{Nothing,Int}=nothing) where {C,N<:Discriminator}
        partitioner = RandomPartitioner(width, n; seed)

        new{C,N}(n, seed, partitioner, Dict{C,N}())
    end
end

# Default discriminator type is Discriminator
MultiDiscriminatorClassifier{C}(args...; kargs...) where C =
    MultiDiscriminatorClassifier{C,Discriminator}(args...; kargs...)

# Default target type is 'Int'
MultiDiscriminatorClassifier(args...; kargs...) = MultiDiscriminatorClassifier{Int}(args...; kargs...)

function train!(model::MultiDiscriminatorClassifier{C,N}, X::AbstractVector{Bool}, y::C) where {C,N<:Discriminator}
    if isnothing(model.partitioner)
        model.partitioner = RandomPartitioner(length(X), model.n; seed=model.seed)
    end

    train!(get!(model.discriminators, y) do
        N(model.partitioner)
    end, X)
end

function train!(model::MultiDiscriminatorClassifier{C}, X::AbstractMatrix{Bool}, y::AbstractVector{C}) where {C}
    size(X, 1) != size(y, 1) && throw(
        DimensionMismatch("Number of observations (rows in X) must match the number of targets (elements in y)"))

    # I could create all necessary discriminators upfront, which could lead to a gain in performance. On the other hand,
    # the current design fits quite nicely with the repetitive training that is common in RL.

    for i in eachindex(y)
        train!(model, X[i, :], y[i])
    end
end

# TODO: Add predict and predict_response umbrella methods where the bleaching threshold is a keyword argument for
#       convenience.
predict(model::MultiDiscriminatorClassifier{C}, X; b=:none) where {C} = predict(model, X, b)
predict_response(model::MultiDiscriminatorClassifier{C}, X; b=:none) where {C} = predict_response(model, X, b)

# TODO: Check if model was trained
function predict(model::MultiDiscriminatorClassifier{C}, X::AbstractVector{Bool}, b::Int=0) where {C}
    largest_response = -1
    best_category = first(keys(model.discriminators))

    for (category, discriminator) in model.discriminators
        response = predict(discriminator, X; b)

        if response > largest_response
            largest_response = response
            best_category = category
        end
    end

    return best_category
end

function predict(model::MultiDiscriminatorClassifier{C}, X::AbstractMatrix{Bool}, b::Int) where {C}
    return C[predict(model, row; b) for row in eachrow(X)]
end

function predict_response(model::MultiDiscriminatorClassifier{C}, X::AbstractVecOrMat{Bool}, b::Int) where {C}
    Dict(k => predict(d, X; b) for (k,d) in model.discriminators)
end

function response_tie(responses)
    # Single-value vector cannot have a tie
    length(responses) == 1 && return false

    p = sortperm(responses; rev=true)

    # Check if the two largest responses are the same
    return responses[p[1]] == responses[p[2]]
end

const bleaching_policies = [:none, :linear]

function predict_response(model::MultiDiscriminatorClassifier{C}, X::AbstractVecOrMat{Bool}, b::Symbol) where {C}
    if b == :none
        return predict_response(model, X, 0)
    elseif b == :linear
        return predict_response_linear_bleaching(model, X)
    else
        throw(DomainError(b, "Unknown bleaching policy"))
    end
end

# Prediction with linear search for the bleaching threshold
# TODO: This could simply be called predict_response with a b=:linear argument
function predict_response_linear_bleaching(model::MultiDiscriminatorClassifier{C}, X::AbstractVector{Bool}) where {C}
    b = zero(Int)
    responses = predict_response(model, X, b)

    while true
        response_values = collect(values(responses))

        if all(isequal(0), response_values) || !response_tie(response_values)
            return responses
        end

        b += 1
        responses = predict_response(model, X, b)
    end
end

function predict_response_linear_bleaching(model::MultiDiscriminatorClassifier{C}, X::AbstractMatrix{Bool}) where {C}
    response = Dict(k => Int[] for k in keys(model.discriminators))

    mergewith!(push!, response, collect(predict_response_linear_bleaching(model, x) for x in eachrow(X))...)
end

function predict(model::MultiDiscriminatorClassifier{C}, X::AbstractVecOrMat{Bool}, b::Symbol) where {C}
    if b == :none
        return predict(model, X, 0)
    elseif b == :linear
        return predict_bleached(model, X)
    else
        throw(DomainError(b, "Unknown bleaching policy"))
    end
end

function predict_bleached(model::MultiDiscriminatorClassifier{C}, X::AbstractVector{Bool}) where {C}
    responses = predict_response_linear_bleaching(model, X)

    if all(isequal(0), values(responses))
        return rand(keys(responses))
    else
        return argmax(responses)
    end
end

# This is identical to other matrix-predicts. Maybe this could already be provided by the model interface
function predict_bleached(model::MultiDiscriminatorClassifier{C,BleachingDiscriminator}, X::AbstractMatrix{Bool}) where {C}
    return C[predict_bleached(model, row) for row in eachrow(X)]
end