text:
### Edge and Gradient related tests ###

using Images, Base.Test, Color

EPS = 1e-14

## Checkerboard array, used to test image gradients

let
    white{T}(::Type{T}) = one(T)
    black{T}(::Type{T}) = zero(T)
    white{T<:Unsigned}(::Type{T}) = typemax(T)
    black{T<:Unsigned}(::Type{T}) = typemin(T)

    global checkerboard
    function checkerboard{T}(::Type{T}, sq_width::Integer, count::Integer)
        wh = fill(white(T), (sq_width,sq_width))
        bk = fill(black(T), (sq_width,sq_width))
        bw = [wh bk; bk wh]
        vert = repmat(bw, (count>>1), 1)
        isodd(count) && (vert = vcat(vert, [wh bk]))
        cb = repmat(vert, 1, (count>>1))
        isodd(count) && (cb = hcat(cb, vert[:,1:sq_width]))
        cb
    end

    checkerboard(sq_width::Integer, count::Integer) = checkerboard(Uint8, sq_width, count)
end

SZ=5

cb_array    = checkerboard(SZ,3)
cb_image_xy = grayim(cb_array)
cb_image_yx = grayim(cb_array)
cb_image_yx["spatialorder"] = ["y","x"]

cb_image_gray = convert(Image{Gray}, cb_image_xy)
cb_image_rgb = convert(Image{RGB}, cb_image_xy)
cb_image_rgb2 = convert(Image{RGB{Float64}}, cb_image_xy)

for method in ["sobel", "prewitt", "ando3", "ando4", "ando5", "ando4_sep", "ando5_sep"]
    ## Checkerboard array

    (agx, agy) = imgradients(cb_array, method)
    amag = magnitude(agx, agy)
    agphase = phase(agx, agy)
    @test (amag, agphase) == magnitude_phase(agx, agy)

    @test agx[1,SZ]   < 0.0   # white to black transition
    @test agx[1,2*SZ] > 0.0   # black to white transition
    @test agy[SZ,1]   < 0.0   # white to black transition
    @test agy[2*SZ,1] > 0.0   # black to white transition

    # Test direction of increasing gradient
    @test cos(agphase[1,SZ])   - (-1.0) < EPS   # increasing left  (=  pi   radians)
    @test cos(agphase[1,2*SZ]) -   1.0  < EPS   # increasing right (=   0   radians)
    @test sin(agphase[SZ,1])   -   1.0  < EPS   # increasing up    (=  pi/2 radians)
    @test sin(agphase[2*SZ,1]) - (-1.0) < EPS   # increasing down  (= -pi/2 radians)

    # Test that orientation is perpendicular to gradient
    aorient = orientation(agx, agy)
    @test all((cos(agphase).*cos(aorient) .+ sin(agphase).*sin(aorient) .< EPS) |
                ((agphase .== 0.0) & (aorient .== 0.0)))  # this part is where both are
                                                          # zero because there is no gradient

    ## Checkerboard Image with row major order

    (gx, gy) = imgradients(cb_image_xy, method)
    mag = magnitude(gx, gy)
    gphase = phase(gx, gy)
    @test (mag, gphase) == magnitude_phase(gx, gy)

    @test gx[SZ,1]   < 0.0   # white to black transition
    @test gx[2*SZ,1] > 0.0   # black to white transition
    @test gy[1,SZ]   < 0.0   # white to black transition
    @test gy[1,2*SZ] > 0.0   # black to white transition

    @test cos(gphase[SZ,1])   - (-1.0) < EPS   # increasing left  (=  pi   radians)
    @test cos(gphase[2*SZ,1]) -   1.0  < EPS   # increasing right (=   0   radians)
    @test sin(gphase[1,SZ])   -   1.0  < EPS   # increasing up    (=  pi/2 radians)
    @test sin(gphase[1,2*SZ]) - (-1.0) < EPS   # increasing down  (= -pi/2 radians)

    # Test that orientation is perpendicular to gradient
    orient = orientation(gx, gy)
    @test all((cos(gphase).*cos(orient) .+ sin(gphase).*sin(orient) .< EPS) |
                ((gphase .== 0.0) & (orient .== 0.0)))  # this part is where both are
                                                        # zero because there is no gradient

    ## Checkerboard Image with column-major order

    (gx, gy) = imgradients(cb_image_yx, method)
    mag = magnitude(gx, gy)
    gphase = phase(gx, gy)
    @test (mag, gphase) == magnitude_phase(gx, gy)

    @test gx[1,SZ]   < 0.0   # white to black transition
    @test gx[1,2*SZ] > 0.0   # black to white transition
    @test gy[SZ,1]   < 0.0   # white to black transition
    @test gy[2*SZ,1] > 0.0   # black to white transition

    # Test direction of increasing gradient
    @test cos(gphase[1,SZ])   - (-1.0) < EPS   # increasing left  (=  pi   radians)
    @test cos(gphase[1,2*SZ]) -   1.0  < EPS   # increasing right (=   0   radians)
    @test sin(gphase[SZ,1])   -   1.0  < EPS   # increasing up    (=  pi/2 radians)
    @test sin(gphase[2*SZ,1]) - (-1.0) < EPS   # increasing down  (= -pi/2 radians)

    # Test that orientation is perpendicular to gradient
    orient = orientation(gx, gy)
    @test all((cos(gphase).*cos(orient) .+ sin(gphase).*sin(orient) .< EPS) |
                ((gphase .== 0.0) & (orient .== 0.0)))  # this part is where both are
                                                        # zero because there is no gradient

    ## Checkerboard Image with Gray pixels

    (gxg, gyg) = imgradients(cb_image_gray, method)
    magg = magnitude(gxg, gyg)
    gphaseg = phase(gxg, gyg)
    @test (magg, gphaseg) == magnitude_phase(gxg, gyg)

    @test gxg[SZ,1]   < 0.0   # white to black transition
    @test gxg[2*SZ,1] > 0.0   # black to white transition
    @test gyg[1,SZ]   < 0.0   # white to black transition
    @test gyg[1,2*SZ] > 0.0   # black to white transition

    @test cos(gphaseg[SZ,1])   - (-1.0) < EPS   # increasing left  (=  pi   radians)
    @test cos(gphaseg[2*SZ,1]) -   1.0  < EPS   # increasing right (=   0   radians)
    @test sin(gphaseg[1,SZ])   -   1.0  < EPS   # increasing up    (=  pi/2 radians)
    @test sin(gphaseg[1,2*SZ]) - (-1.0) < EPS   # increasing down  (= -pi/2 radians)

    # Test that orientation is perpendicular to gradient
    orientg = orientation(gxg, gyg)
    @test all((cos(gphaseg).*cos(orientg) .+ sin(gphaseg).*sin(orientg) .< EPS) |
                ((gphaseg .== 0.0) & (orientg .== 0.0)))  # this part is where both are
                                                        # zero because there is no gradient

    ## Checkerboard Image with RBG pixels

    (gx_rgb, gy_rgb) = imgradients(cb_image_rgb, method)
    mag_rgb = magnitude(gx_rgb, gy_rgb)
    gphase_rgb = phase(gx_rgb, gy_rgb)
    @test (mag_rgb, gphase_rgb) == magnitude_phase(gx_rgb, gy_rgb)

    @test gx_rgb[1,SZ,1]   < 0.0   # white to black transition
    @test gx_rgb[1,2*SZ,1] > 0.0   # black to white transition
    @test gy_rgb[1,1,SZ]   < 0.0   # white to black transition
    @test gy_rgb[1,1,2*SZ] > 0.0   # black to white transition

    @test cos(gphase_rgb[1,SZ,1])   - (-1.0) < EPS   # increasing left  (=  pi   radians)
    @test cos(gphase_rgb[1,2*SZ,1]) -   1.0  < EPS   # increasing right (=   0   radians)
    @test sin(gphase_rgb[1,1,SZ])   -   1.0  < EPS   # increasing up    (=  pi/2 radians)
    @test sin(gphase_rgb[1,1,2*SZ]) - (-1.0) < EPS   # increasing down  (= -pi/2 radians)

    # Test that orientation is perpendicular to gradient
    orient_rgb = orientation(gx_rgb, gy_rgb)
    @test all((cos(gphase_rgb).*cos(orient_rgb) .+ sin(gphase_rgb).*sin(orient_rgb) .< EPS) |
                ((gphase_rgb .== 0.0) & (orient_rgb .== 0.0)))  # this part is where both are
                                                                # zero because there is no gradient

    ## Checkerboard Image with RBG{Float64} pixels

    (gx_rgb, gy_rgb) = imgradients(cb_image_rgb2, method)
    mag_rgb = magnitude(gx_rgb, gy_rgb)
    gphase_rgb = phase(gx_rgb, gy_rgb)
    @test (mag_rgb, gphase_rgb) == magnitude_phase(gx_rgb, gy_rgb)

    @test gx_rgb[1,SZ,1]   < 0.0   # white to black transition
    @test gx_rgb[1,2*SZ,1] > 0.0   # black to white transition
    @test gy_rgb[1,1,SZ]   < 0.0   # white to black transition
    @test gy_rgb[1,1,2*SZ] > 0.0   # black to white transition

    @test cos(gphase_rgb[1,SZ,1])   - (-1.0) < EPS   # increasing left  (=  pi   radians)
    @test cos(gphase_rgb[1,2*SZ,1]) -   1.0  < EPS   # increasing right (=   0   radians)
    @test sin(gphase_rgb[1,1,SZ])   -   1.0  < EPS   # increasing up    (=  pi/2 radians)
    @test sin(gphase_rgb[1,1,2*SZ]) - (-1.0) < EPS   # increasing down  (= -pi/2 radians)

    # Test that orientation is perpendicular to gradient
    orient_rgb = orientation(gx_rgb, gy_rgb)
    @test all((cos(gphase_rgb).*cos(orient_rgb) .+ sin(gphase_rgb).*sin(orient_rgb) .< EPS) |
                ((gphase_rgb .== 0.0) & (orient_rgb .== 0.0)))  # this part is where both are
                                                                # zero because there is no gradient
end

# Create an image with white along diagonals -2:2 and black elsewhere
m = zeros(Uint8, 20,20)
for i = -2:2; m[diagind(m,i)] = 0xff; end

m_xy = grayim(m')
m_yx = grayim(m)
m_yx["spatialorder"] = ["y","x"]

for method in ["sobel", "prewitt", "ando3", "ando4", "ando5", "ando4_sep", "ando5_sep"]
    ## Diagonal array

    (agx, agy) = imgradients(m, method)
    amag = magnitude(agx, agy)
    agphase = phase(agx, agy)
    @test (amag, agphase) == magnitude_phase(agx, agy)

    @test agx[7,9]  < 0.0   # white to black transition
    @test agx[10,8] > 0.0   # black to white transition
    @test agy[10,8] < 0.0   # white to black transition
    @test agy[7,9]  > 0.0   # black to white transition

    # Test direction of increasing gradient
    @test abs(agphase[10,8] -    pi/4 ) < EPS   # lower edge (increasing up-right  =   pi/4 radians)
    @test abs(agphase[7,9]  - (-3pi/4)) < EPS   # upper edge (increasing down-left = -3pi/4 radians)

    # Test that orientation is perpendicular to gradient
    aorient = orientation(agx, agy)
    @test all((cos(agphase).*cos(aorient) .+ sin(agphase).*sin(aorient) .< EPS) |
                ((agphase .== 0.0) & (aorient .== 0.0)))  # this part is where both are
                                                          # zero because there is no gradient

    ## Diagonal Image, row-major order

    (gx, gy) = imgradients(m_xy, method)
    mag = magnitude(gx, gy)
    gphase = phase(gx, gy)
    @test (mag, gphase) == magnitude_phase(gx, gy)

    @test gx[9,7]  < 0.0   # white to black transition
    @test gx[8,10] > 0.0   # black to white transition
    @test gy[8,10] < 0.0   # white to black transition
    @test gy[9,7]  > 0.0   # black to white transition

    # Test direction of increasing gradient
    @test abs(gphase[8,10] -    pi/4 ) < EPS   # lower edge (increasing up-right  =   pi/4 radians)
    @test abs(gphase[9,7]  - (-3pi/4)) < EPS   # upper edge (increasing down-left = -3pi/4 radians)

    # Test that orientation is perpendicular to gradient
    orient = orientation(gx, gy)
    @test all((cos(gphase).*cos(orient) .+ sin(gphase).*sin(orient) .< EPS) |
                ((gphase .== 0.0) & (orient .== 0.0)))  # this part is where both are
                                                        # zero because there is no gradient

    ## Diagonal Image, column-major order

    (gx, gy) = imgradients(m_yx, method)
    mag = magnitude(gx, gy)
    gphase = phase(gx, gy)
    @test (mag, gphase) == magnitude_phase(gx, gy)

    @test gx[7,9]  < 0.0   # white to black transition
    @test gx[10,8] > 0.0   # black to white transition
    @test gy[10,8] < 0.0   # white to black transition
    @test gy[7,9]  > 0.0   # black to white transition

    # Test direction of increasing gradient
    @test abs(gphase[10,8] -    pi/4 ) < EPS   # lower edge (increasing up-right  =   pi/4 radians)
    @test abs(gphase[7,9]  - (-3pi/4)) < EPS   # upper edge (increasing down-left = -3pi/4 radians)

    # Test that orientation is perpendicular to gradient
    orient = orientation(gx, gy)
    @test all((cos(gphase).*cos(orient) .+ sin(gphase).*sin(orient) .< EPS) |
                ((gphase .== 0.0) & (orient .== 0.0)))  # this part is where both are
                                                        # zero because there is no gradient
end

# Nonmaximal suppression

function thin_edges(img)
    # Get orientation
    gx,gy = imgradients(img)
    orient = phase(gx,gy)

    # Do NMS thinning
    thin_edges_nonmaxsup_subpix(img, orient, radius=1.35)
end

function nms_test_horiz_vert(img, which)
    ## which = :horizontal or :vertical

    # Do NMS thinning
    t,s = thin_edges(img)

    # Calc peak location by hand

    # Interpolate values 1.35 pixels left and right
    # Orientation is zero radians -> to the right
    v1 = 6 - 0.35   # slope on right is 1
    v2 = 5 - 0.35*2 # slope on left is 2
    c = 7.0         # peak value

    # solve v = a*r^2 + b*r + c
    a = (v1 + v2)/2 - c
    b = a + c - v2
    r = -b/2a

    @test abs(r - 1/6) < EPS

    # Location and value at peak
    peakloc = r*1.35 + 3
    peakval = a*r^2 + b*r + c

    transposed = spatialorder(img)[1] == "x"
    horizontal = which == :horizontal

    test_axis1 = transposed $ !horizontal

    @test test_axis1 ? all(t[:,[1,2,4,5]] .== 0) : all(t[[1,2,4,5],:] .== 0)
    @test test_axis1 ? all(t[:,3]   .== peakval) : all(t[3,:]   .== peakval)
    @test test_axis1 ? all(s[:,[1,2,4,5]] .== zero(Graphics.Point)) :
                         all(s[[1,2,4,5],:] .== zero(Graphics.Point))

    if transposed
        if which == :horizontal
            @test     [pt.x for pt in s[:,3]]  == [1:5;]
            @test all([pt.y for pt in s[:,3]] .== peakloc)
        else
            @test all([pt.x for pt in s[3,:]] .== peakloc)
            @test     [pt.y for pt in s[3,:]]  == [1:5;]
        end
    else
        if which == :horizontal
            @test     [pt.x for pt in s[3,:]]  == [1:5;]
            @test all([pt.y for pt in s[3,:]] .== peakloc)
        else
            @test all([pt.x for pt in s[:,3]] .== peakloc)
            @test     [pt.y for pt in s[:,3]]  == [1:5;]
        end
    end
end

# Test image: vertical edge
m = [3.0  5.0  7.0  6.0  5.0
     3.0  5.0  7.0  6.0  5.0
     3.0  5.0  7.0  6.0  5.0
     3.0  5.0  7.0  6.0  5.0
     3.0  5.0  7.0  6.0  5.0]

m_xy = grayim(m')
m_yx = grayim(m)
m_yx["spatialorder"] = ["y","x"]

nms_test_horiz_vert(m, :vertical)
nms_test_horiz_vert(m_xy, :vertical)
nms_test_horiz_vert(m_yx, :vertical)

# Test image: horizontal edge
m = m'
m_xy = grayim(m')
m_yx = grayim(m)
m_yx["spatialorder"] = ["y","x"]

nms_test_horiz_vert(m, :horizontal)
nms_test_horiz_vert(m_xy, :horizontal)
nms_test_horiz_vert(m_yx, :horizontal)

function nms_test_diagonal(img)
    # Do NMS thinning
    t,s = thin_edges(img)

    # Calc peak location by hand

    # Interpolate values 1.35 pixels up and left, down and right
    # using bilinear interpolation
    # Orientation is π/4 radians -> 45 degrees up
    fr = 1.35*cos(π/4)
    lower = (7 + fr*(6-7))
    upper = (6 + fr*(5-6))
    v1 = lower + fr*(upper-lower)

    lower = (7 + fr*(5-7))
    upper = (5 + fr*(3-5))
    v2 = lower + fr*(upper-lower)

    c = 7.0         # peak value

    # solve v = a*r^2 + b*r + c
    a = (v1 + v2)/2 - c
    b = a + c - v2
    r = -b/2a

    @test (r - 1/6) < EPS

    transposed = spatialorder(img)[1] == "x"

    # Location and value at peak

    x_peak_offset, y_peak_offset = r*fr, -r*fr
    peakval = a*r^2 + b*r + c

    @test all(diag(data(t))[2:4] .== peakval)  # Edge pixels aren't interpolated here
    @test all(t - diagm(diag(data(t))) .== 0)

    diag_s = copyproperties(s, diagm(diag(data(s))))
    @test s == diag_s

    @test all([pt.x for pt in diag(data(s))[2:4]] - ((2:4) + x_peak_offset) .< EPS)
    @test all([pt.y for pt in diag(data(s))[2:4]] - ((2:4) + y_peak_offset) .< EPS)

end

# Test image: diagonal edge
m = [7.0  6.0  5.0  0.0  0.0
     5.0  7.0  6.0  5.0  0.0
     3.0  5.0  7.0  6.0  5.0
     0.0  3.0  5.0  7.0  6.0
     0.0  0.0  3.0  5.0  7.0]

m_xy = grayim(m')
m_yx = grayim(m)
m_yx["spatialorder"] = ["y","x"]

nms_test_diagonal(m)
nms_test_diagonal(m_xy)
nms_test_diagonal(m_yx)


text:
#=
    Dummy controller
=#

mutable struct Dummy <: AbstractController
    decisions::NamedTuple
    Dummy() = new()
end

### Offline
function initialize_controller!(mg::Microgrid, controller::Dummy, ω::AbstractScenarios)
    # Preallocation
    preallocate!(mg, controller)
    return controller
end

### Online
function compute_operation_decisions!(h::Int64, y::Int64, s::Int64, mg::Microgrid, controller::Dummy)
    return controller
end

ext: jl


text:
using Optim

using SteinDistributions: SteinDistribution, gradlogdensity, logdensity, numdimensions

function getapproxmode(target::SteinDistribution;
                       x0::Array{Float64,1}=zeros(numdimensions(target)))
    function objective(x::Vector)
        -logdensity(target, x)
    end
    function gradient!(x::Vector, gradientres::Vector)
        gradient = gradlogdensity(target, x)
        for ii in 1:length(gradient)
            gradientres[ii] = -gradient[ii]
        end
    end
    res = optimize(objective,
                   gradient!,
                   x0,
                   method = LBFGS())
    res.minimizer
end


text:
module BuoyancyModels

export
    Buoyancy, BuoyancyTracer, SeawaterBuoyancy, buoyancy_perturbation,
    LinearEquationOfState, RoquetIdealizedNonlinearEquationOfState, TEOS10,
    ∂x_b, ∂y_b, ∂z_b, buoyancy_perturbation, x_dot_g_b, y_dot_g_b, z_dot_g_b,
    top_buoyancy_flux,
    buoyancy_frequency_squared,
    BuoyancyField

using Printf
using Oceananigans.Grids
using Oceananigans.Operators
using Oceananigans.BoundaryConditions: getbc

import SeawaterPolynomials: ρ′, thermal_expansion, haline_contraction

# Physical constants
# https://en.wikipedia.org/wiki/Gravitational_acceleration#Gravity_model_for_Earth (30 Oct 2019)
const g_Earth = 9.80665

"""
    AbstractBuoyancyModel{EOS}

Abstract supertype for buoyancy models.
"""
abstract type AbstractBuoyancyModel{EOS} end

"""
    AbstractEquationOfState

Abstract supertype for equations of state.
"""
abstract type AbstractEquationOfState end

function validate_buoyancy(buoyancy, tracers)
    req_tracers = required_tracers(buoyancy)

    all(tracer ∈ tracers for tracer in req_tracers) ||
        error("$(req_tracers) must be among the list of tracers to use $(typeof(buoyancy).name.wrapper)")

    return nothing
end

include("buoyancy.jl")
include("no_buoyancy.jl")
include("buoyancy_tracer.jl")
include("seawater_buoyancy.jl")
include("linear_equation_of_state.jl")
include("nonlinear_equation_of_state.jl")
include("g_dot_b.jl")
include("buoyancy_field.jl")

end # module


text:
include("$(pwd())/startup.jl")

#fp = "./ExampleFiles/STOFDATA/" # All files in this path will be processed

fp = "/media/wiebke/Extreme SSD/PSM_vs_PTR3/Data/apiTOFdata/CLOUD10/run1734_02/"

filefilterRegexp = r"\.h5$"
#rf = "./ExampleFiles/STOFDATA/2017-05-24_12h50m39_NH4.h5"  # The mass scale from this file defines the mass scale of all

rf = "$(fp)APi4_Data_2015.10.20-13h34m44s.h5"

#masslist = MasslistFunctions.loadMasslist("./ExampleFiles/MASSLISTS/exampleMasslistSTOF.csv")

masslist = MasslistFunctions.loadMasslist("$(fp)MassList_NO3-_AP_03-01-20.csv")

cr = [37 137]

# alternatively: use an auto generated masslist
# masslistMasses, masslistElements, masslistElementsMasses, masslistCompositions = createMassList(C=0:20, O=0:20, N=0:1, allowRadicals=false) #
(masslistMasses, masslistElements, masslistElementsMasses, masslistCompositions) = masslist
s = (masslistMasses.>0) .& ( masslistMasses.<600)
masslistMasses = masslistMasses[s]
masslistCompositions = masslistCompositions[s,:]

####################### END OF SETTINGS ###############################################################

####################### Processing sequence ###########################################################

correctMassScaleAndExtractSumSpec(
    fp,
    masslistMasses,
    masslistElements,
    masslistElementsMasses,
    masslistCompositions,
    rf,
    cr,
    filefilterRegexp=filefilterRegexp,
    onlyUseAverages = true,
    plotControlMass = true,
    recalibInterval = 300,
    resolution = 1500,
    firstNFiles=0,
    lastNFiles = 0
    )

baselineAndPeakshape(
    fp,
    peakshapeRegions=4,
    peakshapeRegionStretch=1,
    peakshapeQuantileValue = 0.2,
    peakfindingNoiseThresholdValue = 2,
    peakfindingSignalLimit=0.01
    )

mtrx = deconvolute(
    fp,
    calcTransposed = false
    )


text:
using PackageCompiler: bundle_artifacts, bundle_julia_libraries, create_pkg_context, create_sysimage
using Libdl

builddir = length(ARGS) > 0 ? ARGS[1] : "build"
project = joinpath(@__DIR__, "Clustering.jl")
sysimage_path = joinpath(builddir, "lib", "sysimage.$(Libdl.dlext)")

mkpath(builddir)

@info "Building Julia system image..." builddir=builddir project=project sysimage_path=sysimage_path

# From PackageCompiler.jl `create_app(...)`
ctx = create_pkg_context(project)
bundle_artifacts(ctx, builddir)
bundle_julia_libraries(builddir)

cp(joinpath(project, "Project.toml"), joinpath(builddir, "lib", "Project.toml"))

create_sysimage(
    :Clustering,
    sysimage_path = sysimage_path,
    project = project,
    # NOTE: Test only kmeans, others failing on Julia 1.5
    # TODO: Enable other tests
    precompile_execution_file = joinpath(project, "test", "kmeans.jl"),
    # precompile_execution_file = joinpath(project, "test", "runtests.jl"),
    # Optionally, to minimize the image size.
    # incremental = false,
    # filter_stdlibs = true,
    # Optionally, to target a specific CPU.
    # Use `julia -C help` to list CPU targets.
    # cpu_target = skylake
)

max_stars_repo_name: maxmouchet/julia-node-extension-demo

max_issues_repo_name: maxmouchet/julia-node-extension-demo

max_forks_repo_name: maxmouchet/julia-node-extension-demo


text:
# A type for keeping track of the current line number when printing Exprs
struct LineNumberIO <: IO
    io::IO
    linenos::Vector{Union{Missing,Int}} # source line number for each printed line of the Expr
    file::Symbol   # used to avoid confusion when we are in expanded macros
end

LineNumberIO(io::IO, line::Integer, file::Symbol) = LineNumberIO(io, Union{Missing,Int}[line], file)
LineNumberIO(io::IO, line::Integer, file::AbstractString) = LineNumberIO(io, line, Symbol(file))

const LNIO = Union{LineNumberIO, IOContext{LineNumberIO}}

# Instead of printing the source line number to `io.io`, associate it with the
# corresponding line of the printout
function Base.show_linenumber(io::LineNumberIO, line, file)
    if file == io.file
        # Count how many newlines we've encountered
        data = io.io.data
        nlines = count(isequal(UInt8('\n')), data) + 1 # TODO: O(N^2), optimize? See below
        # If there have been more printed lines than assigned line numbers, fill
        # with `missing`
        while nlines > length(io.linenos)
            push!(io.linenos, missing)
        end
        # Record this line number
        io.linenos[nlines] = line
    end
    return nothing
end
Base.show_linenumber(io::LNIO, line, file) = Base.show_linenumber(io.io, line, file)
Base.show_linenumber(io::LNIO, line, ::Nothing) = nothing
Base.show_linenumber(io::LNIO, line) = nothing

# TODO? intercept `\n` here and break the result up into lines at writing time?
Base.write(io::LineNumberIO, x::UInt8) = write(io.io, x)

# See docstring below
function expression_lines(method::Method)
    def = definition(method)
    if def === nothing
        # If the expression is not available, use the source text. This happens for methods in e.g., boot.jl.
        src, line1 = definition(String, method)
        methstrings = split(chomp(src), '\n')
        return Vector(range(Int(line1), length=length(methstrings))), line1, methstrings
    end
    def = unwrap(def)
    # We'll use the file in LineNumberNodes to make sure line numbers refer to the "outer"
    # method (and does not get confused by macros etc). Because of symlinks and non-normalized paths,
    # it's more reliable to grab the first LNN for the template filename than to use method.file.
    lnn = findline(def)
    mfile = lnn === nothing ? method.file : lnn.file
    buf = IOBuffer()
    io = LineNumberIO(buf, method.line, mfile) # deliberately using the in-method numbering
    print(io, def)
    seek(buf, 0)
    methstrings = readlines(buf)
    linenos = io.linenos
    while length(linenos) < length(methstrings)
        push!(linenos, missing)
    end
    if startswith(methstrings[1], ":(")
        # Chop off the Expr-quotes from the printing
        methstrings[1] = methstrings[1][3:end]
        methstrings[end] = methstrings[end][1:end-1]
    end
    # If it prints with `function`, adjust numbering for the signature line
    # Note that this works independently of whether it's written as an `=` method
    # in the source code (the contents of that line may not be the same but that's
    # true quite generally)
    startswith(methstrings[1], "function") && (linenos[1] -= 1)
    @assert issorted(skipmissing(linenos))
    # Strip out blank lines from the printed expression
    # These arise from the fact that we intercepted the printing of LineNumberNodes
    keeplinenos, keepstrings = Union{Missing,Int}[], String[]
    for (i, line) in enumerate(methstrings)
        if !all(isspace, line)
            push!(keepstrings, line)
            ln = linenos[i]
            if ismissing(ln) && i > 1
                # Line numbers get associated with the missing LineNumberNodes rather than
                # the succeeding expression line. Thus we look to the previous entry.
                # Deliberately go back only one entry.
                ln = linenos[i-1]
            end
            push!(keeplinenos, ln)
        end
    end
    linenos, methstrings = keeplinenos, keepstrings
    # Fill in missing lines where possible. There is no line info for things like
    # `end`, `catch` and similar; these are subsumed by the block structure.
    # So we try to assign line numbers to these missing elements.
    # However, source text like
    #   for i = 1:5 s += i end
    # and
    #   for i = 1:5 s += i
    #   end
    # and
    #   for i = 1:5
    #       s +=i
    #   end
    # all get printed in the latter format---so the matching will work only under
    # certain conditions.
    lastknown = 1
    for i = 2:length(linenos)-1
        if ismissing(linenos[i])
            if !ismissing(linenos[i+1])
                # Process the previous block of missing statements
                # If the printout has advanced by the same number of lines as the source,
                # we know what the answer must be.
                Δidx = i+1 - lastknown
                Δsrc = linenos[i+1] - linenos[lastknown]
                if Δsrc == Δidx
                    for j = lastknown+1:i
                        linenos[j] = linenos[j-1] + 1
                    end
                end
            end
        else
            lastknown = i
        end
    end
    _, line1 = whereis(method)
    return linenos, line1, keepstrings
end

"""
    linenos, line1, methstrings = expression_lines(frame)

Compute the source lines associated with each printed line of the expression
associated with the method executed in `frame`.
`methstrings` is a vector of strings, one per line of the expression.
`linenos` is a vector in 1-to-1 correspondence with `methstrings`,
 containing either the *compiled* line number or `missing`
if the line number is not available. `line1` contains the *actual* (current) line
number of the first line of the method body.
"""
function expression_lines(frame::Frame)
    m = scopeof(frame)
    mlinenos, line1, msrc = expression_lines(m)
    isdefined(m, :generator) || return mlinenos, line1, msrc
    # The rest of this is specific to generated functions.
    # Call the generator to get the expression. First we have to build up the arguments.
    g = m.generator
    gg = g.gen
    vars = JuliaInterpreter.locals(frame)
    ggargs = []
    # Static parameters come first
    for v in vars
        v.isparam || continue
        push!(ggargs, v.value)
    end
    # Slots are next. Naturally, the generator takes only their types, not their values.
    for v in vars
        v.isparam && continue
        push!(ggargs, JuliaInterpreter._Typeof(v.value))
    end
    # Call the generator
    def = gg(ggargs...)
    # `def` contains just the body. Wrap in a `function` to ensure proper indentation,
    # and then print it.
    def = Expr(:function, :(generatedtmp()), def)
    buf = IOBuffer()
    print(buf, def)  # there are no linenos, so no need for LineNumberIO
    seek(buf, 0)
    glines = readlines(buf)
    # Extract the signature line from `msrc`, and paste the body in
    gsrc = [msrc[1]; glines[2:end]]
    # Assign line numbers. Other than the first line, there *aren't* any, so use `missing`
    linenos = [mlinenos[1]; fill(missing, length(gsrc)-1)]
    return linenos, line1, gsrc
end

function show_code(term, frame, deflines, nlines)
    width = displaysize(term)[2]
    method = scopeof(frame)
    linenos, line1, showlines = deflines   # linenos is in "compiled" numbering, line1 in "current" numbering
    offset = line1 - method.line           # compiled + offset -> current
    known_linenos = skipmissing(linenos)
    nd = isempty(known_linenos) ? 0 : ndigits(offset + maximum(known_linenos))
    line = linenumber_unexpanded(frame)    # this is in "compiled" numbering
    # lineidx = searchsortedfirst(linenos, line)
    # Can't use searchsortedfirst with missing
    lineidx = 0
    for (i, l) in enumerate(linenos)
        if !ismissing(l) && l >= line
            lineidx = i
            break
        end
    end
    idxrange = max(1, lineidx-2):min(length(linenos), lineidx+2)
    iochar = IOBuffer()
    for idx in idxrange
        thisline, codestr = linenos[idx], showlines[idx]
        print(term, breakpoint_style(frame.framecode, thisline))
        if ismissing(thisline)
            print(term, " "^nd)
        else
            linestr = lpad(thisline + offset, nd)
            printstyled(term, linestr; color = thisline==line ? Base.warn_color() : :normal)
        end
        linestr = linetrunc(iochar, codestr, width-nd-3)
        print(term, "  ", linestr, '\n')
    end
    return length(idxrange)
end

# Limit output to a single line
function linetrunc(iochar::IO, linestr, width)
    nchars = 0
    for c in linestr
        if nchars == width-2
            print(iochar, '…')
            break
        else
            print(iochar, c)
        end
        nchars += 1
    end
    return String(take!(iochar))
end
linetrunc(linestr, width) = linetrunc(IOBuffer(), linestr, width)

function breakpoint_style(framecode, thisline)
    rng = coderange(framecode, thisline)
    style = ' '
    breakpoints = framecode.breakpoints
    for i in rng
        if isassigned(breakpoints, i)
            bps = breakpoints[i]
            if !bps.isactive
                if bps.condition === JuliaInterpreter.falsecondition  # removed
                else
                    style = style == ' ' ? 'd' : 'm'   # disabled
                end
            else
                if bps.condition === JuliaInterpreter.truecondition
                    style = style == ' ' ? 'b' : 'm'   # unconditional
                else
                    style = style == ' ' ? 'c' : 'm'   # conditional
                end
            end
        end
    end
    return style
end
breakpoint_style(framecode, ::Missing) = ' '

### Header display

function HeaderREPLs.print_header(io::IO, header::RebugHeader)
    if header.nlines != 0
        HeaderREPLs.clear_header_area(io, header)
    end
    iocount = IOBuffer()  # for counting lines
    for s in (io, iocount)
        if !isempty(header.warnmsg)
            printstyled(s, header.warnmsg, '\n'; color=Base.warn_color())
        end
        if !isempty(header.errmsg)
            printstyled(s, header.errmsg, '\n'; color=Base.error_color())
        end
        if header.current_method != dummymethod
            printstyled(s, header.current_method, '\n'; color=:light_magenta)
        end
        if header.uuid != dummyuuid
            data = stored[header.uuid]
            ds = displaysize(io)
            printer(args...) = printstyled(args..., '\n'; color=:light_blue)
            for (name, val) in zip(data.varnames, data.varvals)
                # Make sure each only spans one line
                if val === nothing
                    val = "nothing"
                end
                try
                    printf_maxsize(printer, s, "  ", name, " = ", val; maxlines=1, maxchars=ds[2]-1)
                catch # don't error just because a print method is borked
                    printstyled(s, "  ", name, " errors in its show method"; color=:red)
                end
            end
        end
    end
    header.nlines = count_display_lines(iocount, displaysize(io))
    header.warnmsg = ""
    header.errmsg = ""
    return nothing
end

function HeaderREPLs.print_header(io::IO, header::InterpretHeader)
    if header.nlines != 0
        HeaderREPLs.clear_header_area(io, header)
    end
    header.frame == nothing && return nothing
    frame, Δ = frameoffset(header.frame, header.leveloffset)
    header.leveloffset -= Δ
    frame === nothing && return nothing
    iocount = IOBuffer()  # for counting lines
    for s in (io, iocount)
        ds = displaysize(io)
        printer(args...) = printstyled(args..., '\n'; color=:light_blue)
        if !isempty(header.warnmsg)
            printstyled(s, header.warnmsg, '\n'; color=Base.warn_color())
        end
        if !isempty(header.errmsg)
            printstyled(s, header.errmsg, '\n'; color=Base.error_color())
        end
        indent = ""
        f = root(frame)
        while f !== nothing
            scope = scopeof(f)
            if f === frame
                printstyled(s, indent, scope, '\n'; color=:light_magenta, bold=true)
            else
                printstyled(s, indent, scope, '\n'; color=:light_magenta)
            end
            indent *= ' '
            f = f.callee
        end
        for (i, var) in enumerate(JuliaInterpreter.locals(frame))
            name, val = var.name, var.value
            name == Symbol("#self#") && (isa(val, Type) || sizeof(val) == 0) && continue
            name == Symbol("") && (name = "@_" * string(i))
            if val === nothing
                val = "nothing"
            end
            try
                printf_maxsize(printer, s, "  ", name, " = ", val; maxlines=1, maxchars=ds[2]-1)
            catch # don't error just because a print method is borked
                printstyled(s, "  ", name, " errors in its show method"; color=:red)
            end
        end
    end
    header.nlines = count_display_lines(iocount, displaysize(io))
    header.warnmsg = ""
    header.errmsg = ""
    return nothing
end

function frameoffset(frame, offset)
    while offset > 0
        cframe = frame.caller
        cframe === nothing && break
        frame = cframe
        offset -= 1
    end
    return frame, offset
end

function linenumber_unexpanded(frame)
    framecode, pc = frame.framecode, frame.pc
    scope = framecode.scope::Method
    codeloc = JuliaInterpreter.codelocation(framecode.src, pc)
    codeloc == 0 && return nothing
    lineinfo = framecode.src.linetable[codeloc]
    while lineinfo.file != scope.file && codeloc > 0
        codeloc -= 1
        lineinfo = framecode.src.linetable[codeloc]
    end
    return JuliaInterpreter.getline(lineinfo)
end


text:
using UnsteadyFlowSolvers

alphadef = ConstDef(4. *pi/180)

hdef = SinDef(0., 0.05, 3.93, 0.)

udef = ConstDef(1.)

full_kinem = KinemDef(alphadef, hdef, udef)

pvt = 0.25

geometry = "sd7003.dat"

surf = TwoDSurf(geometry, pvt, full_kinem)

curfield = TwoDFlowField()

dtstar = find_tstep(hdef)

t_tot = 5. *pi/hdef.k

nsteps =Int(round(t_tot/dtstar))+1

startflag = 0

writeflag = 1

writeInterval = t_tot/20.

#delvort = delSpalart(500, 12, 1e-5)
delvort = delNone()

mat, surf, curfield = lautat(surf, curfield, nsteps, dtstar,startflag, writeflag, writeInterval, delvort, wakerollup=1)

makeVortPlots2D()

makeForcePlots2D()

cleanWrite()

max_stars_repo_name: MaximeUofG/UNSflowTutorials

max_issues_repo_name: MaximeUofG/UNSflowTutorials

max_forks_repo_name: MaximeUofG/UNSflowTutorials


text:
datadir = joinpath(srcdir, "..", "data") # TODO: make OS agnostic

"""Load a well-known public regression dataset with nominal features."""
function load_boston()
    df = CSV.read(joinpath(datadir, "Boston.csv"), copycols=true,
                  categorical=true, allowmissing=:none)
    return SupervisedTask(verbosity=0, data=df,
                          target=:MedV,
                          ignore=[:Chas,],
                          is_probabilistic=false)
end

"""Load a reduced version of the well-known Ames Housing task,
having six numerical and six categorical features."""
function load_reduced_ames()
    df = CSV.read(joinpath(datadir, "reduced_ames.csv"), copycols=true,
                  categorical=true, allowmissing=:none)
    df[:target] = exp.(df[:target])
    # TODO: uncomment following after julia #29501 is resolved
#    df.OverallQual = categorical(df.OverallQual, ordered=true)
#    df[:GarageCars] = categorical(df[:GarageCars], ordered=true)
#    df[:YearBuilt] = categorical(df[:YearBuilt], ordered=true)
#    df[:YearRemodAdd] = categorical(df[:YearRemodAdd], ordered=true)
    return SupervisedTask(verbosity=0, data=df,
                          target=:target,
                          is_probabilistic=false)
end

"""Load the full version of the well-known Ames Housing task."""
function load_ames()
    df = CSV.read(joinpath(datadir, "ames.csv"), copycols=true,
                  categorical=true, allowmissing=:none)
    df[:target] = exp.(df[:target])
    return SupervisedTask(verbosity=0, data=df,
                          target=:target,
                          ignore=[:Id,],
                          is_probabilistic=false)
end

"""Load a well-known public classification task with nominal features."""
function load_iris()
    df = CSV.read(joinpath(datadir, "iris.csv"), pool=true, copycols=true,
                  categorical=true, allowmissing=:none)
    return SupervisedTask(verbosity=0, data=df,
                          target=:target,
                          is_probabilistic=false)
end

"""Load a well-known crab classification dataset with nominal features."""
function load_crabs()
    df = CSV.read(joinpath(datadir, "crabs.csv"), pool=true, copycols=true,
                  categorical=true, allowmissing=:none)
    return SupervisedTask(verbosity=0, data=df,
                          target=:sp,
                          ignore=[:sex, :index],
                          is_probabilistic=true)
end

"""Get some supervised data now!!"""
function datanow()
    X, y = X_and_y(load_boston())
    return (selectrows(X, 1:75), y[1:75])
end


text:
extrema(tg::TextGrid) = (minimum(x -> x.xmin, tg), maximum(x -> x.xmax, tg))

function reindex(interval::Interval, num::Int)
    interval_new = interval
    interval_new.index = num
    interval_new
end


text:
using StatisticalRethinking, CmdStan
#gr(size=(700,700));

ProjDir = rel_path("..", "scripts", "04")
cd(ProjDir)

howell1 = CSV.read(rel_path("..", "data", "Howell1.csv"), delim=';')
df = convert(DataFrame, howell1);

first(df, 5)

first(df, 5)

df2 = filter(row -> row[:age] >= 18, df);

m4_1 = "
  height ~ Normal(μ, σ) # likelihood
  μ ~ Normal(178,20) # prior
  σ ~ Uniform(0, 50) # prior
";

p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 3)
p[1] = density(df2[:height], xlim=(100,250), lab="All heights", xlab="height [cm]", ylab="density")

d1 = Normal(178, 20)
p[2] = plot(100:250, [pdf(d1, μ) for μ in 100:250], lab="Prior on mu")

d2 = Uniform(0, 50)
p[3] = plot(0:0.1:50, [pdf(d2, σ) for σ in 0:0.1:50], lab="Prior on sigma")

plot(p..., layout=(3,1))

sample_mu = rand(d1, 10000)
sample_sigma = rand(d2, 10000)
prior_height = [rand(Normal(sample_mu[i], sample_sigma[i]), 1)[1] for i in 1:10000]
df2 = DataFrame(mu = sample_mu, sigma=sample_sigma, prior_height=prior_height);
first(df2, 5)

density(prior_height, lab="prior_height")

!isfile(joinpath(ProjDir, "m4.1s.jls")) && include(joinpath(ProjDir, "m4.1s.jl"))

chn = deserialize(joinpath(ProjDir, "m4.1s.jls"))

MCMCChains.describe(chn)

density(chn)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
println("testing mesh")
# function test_mesh()
img = rand(100,100)
offset = [0,0]
dist = 20
mesh = create_mesh(img, offset, dist)
i1 = 100 % (dist*sin(pi/3)) / 2
@test mesh.src_nodes[1,:] ≈ [i1,0]
@test mesh.src_nodes[2,:] ≈ [i1,20]
@test mesh.src_nodes[3,:] ≈ [i1,40]
@test mesh.src_nodes[4,:] ≈ [i1,60]
@test mesh.src_nodes[5,:] ≈ [i1,80]
@test mesh.src_nodes[6,:] ≈ [i1,100]
@test mesh.src_nodes[7,:] ≈ [i1 + 20*(sin(pi/3)), 10]
@test maximum(mesh.src_nodes[:,1]) <= 100
@test maximum(mesh.src_nodes[:,2]) <= 100
@test size(mesh.src_nodes, 1) == 33
# test that src_nodes match dst_nodes, initially
function wrapper() # otherwise gives warning if a and b are already defined globally
	for (a, b) in zip(mesh.src_nodes, mesh.dst_nodes)
		@test a == b
	end
end
wrapper()
@test size(mesh.edges, 2) == 81


text:
"""
    _D_gain_method(parameters::Parameters1PL, latent::Latent1D)
"""
function _D_gain_method(parameters::Parameters1PL, latent::Latent1D)
    return _item_expected_information(parameters, latent)
end

"""
    _D_gain_method(parameters::Parameters2PL, latent::Latent1D)
"""
function _D_gain_method(parameters::Parameters2PL, latent::Latent1D)
    old_exp_info = copy(parameters.expected_information)
    return LinearAlgebra.det(old_exp_info + _item_expected_information(parameters, latent)) -
           LinearAlgebra.det(old_exp_info)
end

"""
    D_gain_method(item::AbstractItem, examinee::AbstractExaminee)

# Description
Computes the gain in the determinant of the expected information matrix for an item.

# Arguments
- **`item::AbstractItem`**: The item.
- **`examinee::AbstractExaminee`**: The examinee at which computing the information.

# Output
It returns a `Float64` scalar.

# References
__Ren H, van der Linden WJ, Diao Q. Continuous online item calibration: Parameter recovery and item calibration. Psychometrika. 2017;82:498–522. doi: 10.1007/s11336-017-9553-1.__
"""
function D_gain_method(item::AbstractItem, examinee::AbstractExaminee)
    return _D_gain_method(item.parameters, examinee.latent)
end


text:
## squared_exponential.jl : implementation of squared exponential covariance function

## SquaredExponential ##
struct SquaredExponential{T} <: IsotropicCovarianceStructure{T}
    λ::T
    σ::T
    p::T
end

"""
    SquaredExponential(λ, σ=1, p=2)

Create a squared exponential (Gaussian) covariance structure with correlation length `λ`, (optional) marginal standard deviation `σ` and (optional) `p`-norm.

# Examples
```jldoctest
julia> s1 = SquaredExponential(0.1)
Gaussian (λ=0.1, σ=1.0, p=2.0)

julia> s2 = SquaredExponential(0.1, σ=2.)
Gaussian (λ=0.1, σ=2.0, p=2.0)

```
"""
function SquaredExponential(λ::T where {T<:Real}; σ=1.0::T where {T<:Real}, p=2::T where {T<:Real})
    λ > 0 || throw(ArgumentError("correlation length λ of squared exponential covariance cannot be negative or zero"))
    σ > 0 || throw(ArgumentError("marginal standard deviation σ of squared exponential covariance cannot be negative or zero"))
    p >= 1 || throw(ArgumentError("in p-norm, p must be greater than or equal to 1"))
	isinf(p) && throw(ArgumentError("in p-norm, p cannot be infinity"))
    SquaredExponential{promote_type(typeof(λ),typeof(σ),typeof(p))}(promote(λ,σ,p)...)
end

"""
    Gaussian(λ, σ=1, p=2)

Create a Gaussian (squared exponential) covariance structure with correlation length `λ`, (optional) marginal standard deviation `σ` and (optional) `p`-norm.

# Examples
```jldoctest
julia> g1 = Gaussian(0.1)
Gaussian (λ=0.1, σ=1.0, p=2.0)

julia> g2 = Gaussian(0.1, σ=2.)
Gaussian (λ=0.1, σ=2.0, p=2.0)

```
"""
const Gaussian = SquaredExponential

# evaluate Gaussian covariance
function apply(s::SquaredExponential,x::T) where {T<:Real}
    exp(-(x/s.λ)^2)
end

show(io::IO, s::SquaredExponential) = print(io, "Gaussian (λ=$(s.λ), σ=$(s.σ), p=$(s.p))")


text:
using DelimitedFiles
using Dates
using PyPlot, Colors
using NetCDF
include("startup.jl")

fp = "/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/DatafromOtherGroups/DailyTrace/"
(dailyPTRTraces, dailyPTRHeader) = readdlm("$(fp)ptr3Dat_ppt.csv", ',', header = true)
(dailyNO3cimsTraces, dailyNO3cimsHeader) = readdlm("$(fp)C7H8Ox_NO3-CIMS_cm-3.csv", ',', header = true)
dailyNO3cimsTraces[:,2:end] = dailyNO3cimsTraces[:,2:end]
dailyPTRTraces[:,2:end] = dailyPTRTraces[:,2:end].*2.47e7
(FigaeroTracesPart, FigaeroHeaderPart) = readdlm("$(fp)FIGAERO_CIMS_C7_for_BAM.csv", ',', header = true)
(FigaeroTracesGas, FigaeroHeaderGas) = readdlm("$(fp)FIGAERO_CIMS_C7H8Ox_gasphase__for_BAM.csv", ',', header = true)
figaeroPart_time = Dates.DateTime.(FigaeroTracesPart[:,1], "dd-mm-yyyyTHH:MM:SS")
figaerogas_time = Dates.DateTime.(FigaeroTracesGas[:,1], "mm/dd/yyyy HH:MM:SS")

global Xpart = DatasetFunctions.nanmean(FigaeroTracesPart[( hour.(figaeroPart_time).== 0),2:end], 1)
for i in collect(1:23)
	x = DatasetFunctions.nanmean(FigaeroTracesPart[( hour.(figaeroPart_time).== i),2:end], 1)
	global Xpart = vcat(Xpart, x)
end

global Xgas = DatasetFunctions.nanmean(FigaeroTracesGas[( hour.(figaerogas_time).== 0),2:end], 1)
for i in collect(1:23)
	x = DatasetFunctions.nanmean(FigaeroTracesGas[( hour.(figaerogas_time).== i),2:end], 1)
	global Xgas = vcat(Xgas, x)
end

SRRfilename = "/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/DatafromOtherGroups/Flexpart/v03/cluster_series_v3.nc"
SRRutctimes = Dates.Hour.(ncread(SRRfilename, "releases")) .+ Dates.DateTime(2017,12,06)
SRR_DT = SRRutctimes .- Dates.Hour(4)
SRR_z0 = ncread(SRRfilename, "conc_all")[:,:,2][:,3]
mask = Dates.DateTime(2018,05,08,00) .< SRR_DT .< Dates.DateTime(2018,05,23,00)
SRR_z0 = SRR_z0[mask]
SRR_DT = SRR_DT[mask]

global z0 = DatasetFunctions.nanmean(SRR_z0[( hour.(SRR_DT).== 0)], 1)
for i in collect(1:23)
	z = DatasetFunctions.nanmean(SRR_z0[( hour.(SRR_DT).== i)], 1)
	global z0 = vcat(z0, z)
end

#
#open("/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/DatafromOtherGroups/DailyTrace/ptr3Dat.csv", "w") do io
#           writedlm(io, vcat(hcat(["hour"],permutedims(labels)), hcat(collect(0:23), X)))
#       end

figure()
subplot(311)
subplots_adjust(hspace=0.0)
PTRdats1 = dailyPTRTraces[:,2]./transpose(dailyPTRTraces[1,2])
semilogy(dailyPTRTraces[:,1], PTRdats1, "-")
PTR_Ioddats2 = DatasetFunctions.nanmean(hcat(dailyPTRTraces[:,3:5]./transpose(dailyPTRTraces[1,3:5]),Xgas[:,1]./Xgas[1,1]) ,2)
semilogy(dailyPTRTraces[:,1], PTR_Ioddats2, "-")
PTR_Iod_Nitrdats3 = DatasetFunctions.nanmean(hcat(dailyPTRTraces[:,6:end]./transpose(dailyPTRTraces[1,6:end]), Xgas[:,2]./Xgas[1,2], Xgas[:,4]./Xgas[1,4], dailyNO3cimsTraces[:,2:3]./transpose(dailyNO3cimsTraces[1,2:3])),2)
semilogy(dailyPTRTraces[:,1], PTR_Iod_Nitrdats3, "-")
#semilogy(dailyPTRTraces[:,1], Xgas[:,1]./Xgas[1,1], "--", color = "#9467bd")
#semilogy(dailyPTRTraces[:,1], Xgas[:,2]./Xgas[1,2], "--" , color = "#8c564b")
#semilogy(dailyPTRTraces[:,1], Xgas[:,4]./Xgas[1,4], "--", color = "#7f7f7f")
#semilogy(dailyNO3cimsTraces[:,1], dailyNO3cimsTraces[:,2]./dailyNO3cimsTraces[1,2], "-.", color = "#8c564b")
#semilogy(dailyNO3cimsTraces[:,1], dailyNO3cimsTraces[:,3]./dailyNO3cimsTraces[1,3], "-.", color = "#7f7f7f")
legend(vcat(["C7H8_ptr","mean(C7H8O1-O3_H+ and C7H8O4_I-)","mean(C7H8O4-C7H8O7_H+ and C7H8O5,7_I- and C7H8O5,7_NO3-)"], FigaeroHeaderGas[2:3], FigaeroHeaderGas[5], dailyNO3cimsHeader[2:end]))
ylabel("relative change compared to 00:00")
ax = gca()
setp(ax.get_xticklabels(),visible=false)

subplot(312, sharex = ax)
semilogy(InterpolationFunctions.averageSamples(dailyPTRTraces[:,1],2), InterpolationFunctions.averageSamples(Xpart,2; ignoreNaNs = true))
ylabel("particle mass [ug/m³]")
legend(FigaeroHeaderPart[2:end])
ax2 = gca()
setp(ax2.get_xticklabels(),visible=false)

subplot(313, sharex = ax)
plot(collect(0:23), z0, "-.")
xlabel("hour of the day")
ylabel("influence of surface layer [%]")
#ylim(0,1)

newFileDailyTraces = "/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/DatafromOtherGroups/DailyTrace/combined.csv"

# relative change gasphase:
header1 = ["hourOfDay" "C7H8_H+" "mean(C7H8O1-O3_H+ and C7H8O4_I-)" "mean(C7H8O4-C7H8O7_H+ and C7H8O5,7_I- and C7H8O5,7_NO3-)"]
data1 = hcat(dailyPTRTraces[:,1], PTRdats1, PTR_Ioddats2, PTR_Iod_Nitrdats3)
outfile1 = "/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/DatafromOtherGroups/DailyTrace/allGasPhase.csv"

# particle phase
header2 = permutedims(vcat(["hourOfDay"], FigaeroHeaderPart[2:end]))
data2 = hcat(InterpolationFunctions.averageSamples(dailyPTRTraces[:,1],2), InterpolationFunctions.averageSamples(Xpart,2; ignoreNaNs = true))
outfile2 = "/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/DatafromOtherGroups/DailyTrace/particlePhase.csv"

# z0 influence
header3 = ["hourOfDay" "influenceSurfaceLayer_[%]"]
data3 = hcat(collect(0:23), z0)
outfile3 = "/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/DatafromOtherGroups/DailyTrace/surfaceLayerInfluence.csv"

# writedlm :)
writedlm(outfile1, vcat(header1, data1))
writedlm(outfile2, vcat(header2, data2))
writedlm(outfile3, vcat(header3, data3))


text:
#
# This functions compute VdW interactions between two selections
#

#
# Computes the electrostatic interaction of a pair of atoms, given the distance
#

# Given the parameters for each atom

function vdwpair(d2 :: Float64, eps1 :: Float64, rmin1 :: Float64,
                                eps2 :: Float64, rmin2 :: Float64)
  epspair = sqrt( eps1*eps2 )
  rminpair6 = (rmin1 + rmin2)^6
  return vdwpair( d2, epspair, rminpair6 )
end

# If the pair parameters were computed already

function vdwpair(d2 :: Float64, epspair :: Float64, rminpair6 :: Float64 )
  d = sqrt(d2)
  p6 = rminpair6 / d2^3
  p12 = p6*p6
  vdwpair = epspair*( p12 - 2*p6 )
  return vdwpair
end

#
# For a trajectory: x, y, z coordinates read from a DCD are given
#

# First, the general function with all parameters

function vdw(atoms :: Vector{Atom}, sel1 :: Vector{Int}, sel2 :: Vector{Int},
             sides :: Vector{Float64},
             x :: Vector{Float32}, y :: Vector{Float32}, z :: Vector{Float32},
             usecutoff :: Bool, cutoff :: Float64, pbc :: Bool)
  n1 = length(sel1)
  n2 = length(sel2)
  cutoff2 = cutoff^2
  vdw = 0.
  for i in 1:n1
    center = [ x[sel1[i]], y[sel1[i]], z[sel1[i]] ]
    if pbc
      wrap!(sides,x,y,z,center=center,sel=sel2)
    end
    for j in 1:n2
      d2 = (x[sel1[i]] - x[sel2[j]])^2 + (y[sel1[i]] - y[sel2[j]])^2 + (z[sel1[i]] - z[sel2[j]])^2
      epspair = sqrt( atoms[i].eps*atoms[j].eps )
      rminpair6 = ( atoms[i].rmin + atoms[j].rmin )^6
      vdw = vdw + vdwpair(d2,epspair,rminpair6)
      if usecutoff
        if d2 > cutoff2
          vdw = vdw - vdwpair(cutoff2,epspair,rminpair6)
        end
      end
    end
  end
  return coulomb
end

# For a lazy user which does not want to setup the cutoff and the pbc, use defaults

vdw(atoms :: Vector{Atom}, sel1 :: Vector{Int}, sel2 :: Vector{Int},
    sides :: Vector{Float64},
    x :: Vector{Float32}, y :: Vector{Float32}, z :: Vector{Float32} ) =
  vdw(atoms, sel1, sel2, sides, x, y, z, true, 12., true)

# If the user does not provide the sides, than we will not use PBCs

vdw(atoms :: Vector{Atom}, sel1 :: Vector{Int}, sel2 :: Vector{Int},
        x :: Vector{Float32}, y :: Vector{Float32}, z :: Vector{Float32} ) =
  vdw(atoms, sel1, sel2, [0., 0., 0.], x, y, z, true, 12., false)

#
# For computing interaction energies from the coordinates available in the atoms array
#

function vdw2( atoms :: Vector{Atom}, sel1 :: Vector{Int}, sel2 :: Vector{Int},
               usecutoff :: Bool, cutoff :: Float64 )
  n1 = length(sel1)
  n2 = length(sel2)
  vdw = 0.
  cutoff2 = cutoff^2
  for isel1 in 1:n1
    i = sel1[isel1]
    for isel2 in 1:n2
      j = sel2[isel2]
      d2 = (atoms[i].coor[1] - atoms[j].coor[1])^2 +
           (atoms[i].coor[2] - atoms[j].coor[2])^2 +
           (atoms[i].coor[3] - atoms[j].coor[3])^2
      epspair = sqrt( atoms[i].eps*atoms[j].eps )
      rminpair6 = ( atoms[i].rmin + atoms[j].rmin )^6
      vdw = vdw + vdwpair(d2,epspair,rminpair6)
      if usecutoff
        if d2 > cutoff2
          vdw = vdw - vdwpair(cutoff2,epspair,rminpair6)
        end
      end
    end
  end
  return vdw
end

# Mapping to use a single function name

vdw(  atoms :: Vector{Atom}, sel1 :: Vector{Int}, sel2 :: Vector{Int}, usecutoff :: Bool, cutoff :: Float64 ) =
  vdw2( atoms, sel1, sel2, usecutoff, cutoff )

# If the cutoff is not set, use default

vdw( atoms :: Vector{Atom}, sel1 :: Vector{Int}, sel2 :: Vector{Int} ) =
  vdw2( atoms, sel1, sel2, true, 12. )


text:
"""
    _sort!(
        lines::Vector{String},
        lines_x_anchors::Vector{Float64},
        lines_yl_anchors::Vector{Float64},
        lines_yu_anchors::Vector{Float64},
    ) -> Vector{String}

Sort the lines vector using the arrays of anchors. Basic principle: if rectangles cross
by ordinate, sort them by abscissa.

# Arguments
- `lines::Vector{String}`: the lines vector of the highlight (text found in the rectangles
  of the highlight)
- `lines_x_anchors::Vector{Float64}`: the coordinate vector of the abscissa of the left
  side of the highlight rectangles
- `lines_yl_anchors::Vector{Float64}`: the coordinate vector of the ordinate of the
  lower-left corner of the highlight rectangles
- `lines_yu_anchors::Vector{Float64}`: the coordinate vector of the ordinate of the
  upper-left corner of the highlight rectangles

# Returns
- `Vector{String}`: the sorted lines vector of the highlight

# Example
```jldoctest; output = false
using PDFHighlights

lines = ["High", "high", "light"]
quad_x_anchors = [0.21, 0.15, 0.17]
quad_yl_anchors = [0.10, 0.10, 0.10]
quad_yu_anchors = [0.15, 0.12, 0.15]

PDFHighlights.Internal.PDF._sort!(
    lines,
    quad_x_anchors,
    quad_yl_anchors,
    quad_yu_anchors,
) == ["high", "light", "High"]

# output

true
```
"""
function _sort!(
    lines::Vector{String},
    lines_x_anchors::Vector{Float64},
    lines_yl_anchors::Vector{Float64},
    lines_yu_anchors::Vector{Float64},
)::Vector{String}

    start_index = 0
    previous_index = 1
    last_index = lastindex(lines)

    # Sort by `x`
    for index in 2:length(lines)

        previous_index = index - 1

        # Unpack the quad anchors
        yl1, yu1 = lines_yl_anchors[previous_index], lines_yu_anchors[previous_index]
        yl2, yu2 = lines_yl_anchors[index], lines_yu_anchors[index]

        # Check if the lines are crossing
        if (yl1 > yl2 && yl1 < yu2) ||
           (yu1 > yl2 && yl2 < yu2) ||
           (yl2 > yl1 && yl2 < yu1)

            start_index == 0 && (start_index = previous_index)

            # Sort the last chain of lines
            if index == last_index

                perm = sortperm(lines_x_anchors[start_index:index])
                permute!(@view(lines[start_index:index]), perm)
                permute!(@view(lines_yl_anchors[start_index:index]), perm)

                if start_index == 1
                    permute!(@view(lines_x_anchors[start_index:index]), perm)
                    permute!(@view(lines_yu_anchors[start_index:index]), perm)
                end

            end

        # Sort the last chain of lines
        elseif start_index != 0

            perm = sortperm(lines_x_anchors[start_index:previous_index])
            permute!(@view(lines[start_index:previous_index]), perm)

            if start_index == 1
                permute!(@view(lines_x_anchors[start_index:previous_index]), perm)
                permute!(@view(lines_yu_anchors[start_index:previous_index]), perm)
            end

            start_index = 0

        end

    end

    return lines

end

"""
    _sort!(
        highlights::Vector{String},
        comments::Vector{String},
        pages::Vector{Int32},
        highlights_x_anchors::Vector{Float64},
        highlights_yl_anchors::Vector{Float64},
        highlights_yu_anchors::Vector{Float64},
    ) -> Tuple{Vector{String}, Vector{String}}

Sort the highlights and comments vectors using the arrays of anchors and pages.
Basic principle: if highlights cross by ordinate, sort them by abscissa.

# Arguments
- `highlights::Vector{String}`: the highlights
- `comments::Vector{String}`: the comments
- `pages::Vector{Int32}`: the pages
- `highlights_x_anchors::Vector{Float64}`: the coordinate vector of the abscissa of the left
  side of the upper-left rectangle of the highlight
- `highlights_yl_anchors::Vector{Float64}`: the coordinate vector of the ordinate of the
  lower-left corner of the lower-left rectangle of the highlight
- `highlights_yu_anchors::Vector{Float64}`: the coordinate vector of the ordinate of the
  upper-left corner of the upper-left rectangle of the highlight

# Returns
- `Tuple{Vector{String}, Vector{String}}`: the sorted vectors of the highlights and comments

# Example
```jldoctest; output = false
using PDFHighlights

highlights = ["High", "high", "light"]
comments = ["Com", "com", "ment"]
pages = Int32[1, 1, 1]
highlights_x_anchors = [0.21, 0.15, 0.17]
highlights_yl_anchors = [0.10, 0.10, 0.10]
highlights_yu_anchors = [0.12, 0.15, 0.15]

PDFHighlights.Internal.PDF._sort!(
    highlights,
    comments,
    pages,
    highlights_x_anchors,
    highlights_yl_anchors,
    highlights_yu_anchors,
) ==
(
    ["high", "light", "High"],
    ["com", "ment", "Com"],
)

# output

true
```
"""
function _sort!(
    highlights::Vector{String},
    comments::Vector{String},
    pages::Vector{Int32},
    highlights_x_anchors::Vector{Float64},
    highlights_yl_anchors::Vector{Float64},
    highlights_yu_anchors::Vector{Float64},
)::Tuple{Vector{String}, Vector{String}}

    last_index = lastindex(highlights)

    start_index = 1
    previous_index = 0
    current_page = pages[1]

    # Sort by `y`
    for (index, page) in enumerate(pages)
        if page != current_page
            previous_index = index - 1
            if start_index != previous_index
                perm = sortperm(highlights_yu_anchors[start_index:previous_index])
                permute!(@view(highlights[start_index:previous_index]), perm)
                permute!(@view(comments[start_index:previous_index]), perm)
            end
            start_index = index
            current_page = page
        elseif index == last_index && pages[start_index] == pages[index]
            perm = sortperm(highlights_yu_anchors[start_index:index])
            permute!(@view(highlights[start_index:index]), perm)
            permute!(@view(comments[start_index:index]), perm)
        end
    end

    start_index = 0
    previous_index = 1

    # Sort by `x`
    for index in 2:length(highlights)
        previous_index = index - 1
        if pages[index] == pages[previous_index]

            # Unpack the quad anchors
            yl1 = highlights_yl_anchors[previous_index]
            yu1 = highlights_yu_anchors[previous_index]
            yl2 = highlights_yl_anchors[index]
            yu2 = highlights_yu_anchors[index]

            # Check if the highlights are crossing
            if (yl1 > yl2 && yl1 < yu2) ||
               (yu1 > yl2 && yl2 < yu2) ||
               (yl2 > yl1 && yl2 < yu1)

                start_index == 0 && (start_index = previous_index)

                # Sort the last chain of highlights
                if index == last_index
                    perm = sortperm(highlights_x_anchors[start_index:index])
                    permute!(@view(highlights[start_index:index]), perm)
                    permute!(@view(comments[start_index:index]), perm)
                end

            # Sort the current chain of highlights
            elseif start_index != 0

                perm = sortperm(highlights_x_anchors[start_index:previous_index])
                permute!(@view(highlights[start_index:previous_index]), perm)
                permute!(@view(comments[start_index:previous_index]), perm)

                start_index = 0

            end

        end
    end

    return highlights, comments

end