text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

"""
    SourceLocation(file_name::String, line_num::Int64, col_num::Int64)

A specific position in a source file.

## Arguments

- `file_name::String` : the name of the file, or the empty string if there is no
  file associated with this location (e.g., because the source code was provided as
  a memory stream, or through a network connection)

- `line_num::Int64` : number of the line (starting from 1)

- `col_num::Int64` : number of the column (starting from 1)
"""
mutable struct SourceLocation
    file_name::String
    line_num::Int64
    col_num::Int64
    SourceLocation(fn::String, ln::Int64 = 0, cn::Int64 = 0) = new(fn, ln, cn)
end

copy(location::SourceLocation) = SourceLocation(location.file_name, location.line_num, location.col_num)

"""
     @enum KeywordEnum

Enumeration for all the possible keywords recognized by the lexer:
```ditaa
|:-----------------:|:-----------------:|:----------------------:|
| NEW = 1           | PIGMENT = 20      | TRANSFORMATION = 40    |
| FLOAT = 2         | UNIFORM = 21      | IDENTITY = 41          |
| STRING = 3        | CHECKERED = 22    | TRANSLATION = 42       |
| VECTOR = 4        | IMAGE = 23        | ROTATION_X = 43        |
| COLOR = 5         |                   | ROTATION_Y = 44        |
| MATERIAL = 6      |                   | ROTATION_Z = 45        |
| POINTLIGHT = 7    |                   | SCALING = 46           |
|                   |                   |                        |
|:-----------------:|:-----------------:|:----------------------:|
| BRDFS = 10        | CAMERA = 30       | BOOL = 50              |
| DIFFUSE = 11      | ORTHOGONAL = 31   | TRUE = 51              |
| SPECULAR = 12     | PERSPECTIVE = 32  | FALSE = 52             |
|                   |                   |                        |
|                   |                   |                        |
|:-----------------:|:-----------------:|:----------------------:|
| PLANE = 61        | PRINT = 71        |                        |
| SPHERE = 62       | ASSERT = 72       |                        |
| CUBE = 63         |                   |                        |
| TRIANGLE = 64     |                   |                        |
|                   |                   |                        |
|:-----------------:|:-----------------:|:----------------------:|
```
"""
@enum KeywordEnum begin
    NEW = 1
    FLOAT = 2
    STRING = 3
    VECTOR = 4
    COLOR = 5
    MATERIAL = 6
    POINTLIGHT = 7

    BRDFS = 10
    DIFFUSE = 11
    SPECULAR = 12

    PIGMENT = 20
    UNIFORM = 21
    CHECKERED = 22
    IMAGE = 23

    CAMERA = 30
    ORTHOGONAL = 31
    PERSPECTIVE = 32

    TRANSFORMATION = 40
    IDENTITY = 41
    TRANSLATION = 42
    ROTATION_X = 43
    ROTATION_Y = 44
    ROTATION_Z = 45
    SCALING = 46

    BOOL = 50
    TRUE = 51
    FALSE = 52

    PLANE = 61
    SPHERE = 62
    CUBE = 63
    TRIANGLE = 64

    PRINT = 71
    ASSERT = 72
end

KEYWORDS = Dict{String, KeywordEnum}(
    "NEW" => NEW,
    "FLOAT" => FLOAT,
    "STRING" => STRING,
    "VECTOR" => VECTOR,
    "COLOR" => COLOR,
    "MATERIAL" => MATERIAL,
    "POINTLIGHT" => POINTLIGHT,

    "BRDF" => BRDFS,
    "DIFFUSE" => DIFFUSE,
    "SPECULAR" => SPECULAR,

    "PIGMENT" => PIGMENT,
    "UNIFORM" => UNIFORM,
    "CHECKERED" => CHECKERED,
    "IMAGE" => IMAGE,

    "CAMERA" => CAMERA,
    "ORTHOGONAL" => ORTHOGONAL,
    "PERSPECTIVE" => PERSPECTIVE,

    "TRANSFORMATION" => TRANSFORMATION,
    "IDENTITY" => IDENTITY,
    "TRANSLATION" => TRANSLATION,
    "ROTATION_X" => ROTATION_X,
    "ROTATION_Y" => ROTATION_Y,
    "ROTATION_Z" => ROTATION_Z,
    "SCALING" => SCALING,

    "BOOL" => BOOL,
    "TRUE" => TRUE,
    "FALSE" => FALSE,

    "PLANE" => PLANE,
    "SPHERE" => SPHERE,
    "CUBE" => CUBE,
    "TRIANGLE" => TRIANGLE,

    "PRINT" => PRINT,
    "ASSERT" => ASSERT,
)

"""
    KeywordToken(keyword::KeywordEnum)

A token containing a keyword of the Photorealistic Object Applications Language.

See also: [`KeywordEnum`](@ref)
"""
struct KeywordToken
    keyword::KeywordEnum
end

"""
    IdentifierToken(identifier::String)

A token containing an identifier, i.e. a name of a variable.
"""
struct IdentifierToken
    identifier::String
end

"""
     StringToken(string::String)

A token containing a literal string, i.e. a sentence placed inside
two  double quotes ("...") symbols.
"""
struct StringToken
    string::String
end

"""
    LiteralNumberToken(number::Float64)

A token containing a literal number.
"""
struct LiteralNumberToken
    number::Float64
end

"""
    SymbolToken(symbol::String)

A token containing a recognised symbol by the Photorealistic Object
Applications Language.
"""
struct SymbolToken
    symbol::String
end

"""
    StopToken()

A token signalling the end of a file.
"""
struct StopToken
end

"""
    Token(
          location::SourceLocation,
          value::Union{
               KeywordToken,
               IdentifierToken,
               StringToken,
               LiteralNumberToken,
               SymbolToken,
               StopToken}
          )

A lexical token, used when parsing a scene file.

## Arguments

- `location::SourceLocation`: location of the last char read

- `value` : one of the basic 6 token types:
  - [`KeywordToken`](@ref)
  - [`IdentifierToken`](@ref)
  - [`StringToken`](@ref)
  - [`LiteralNumberToken`](@ref)
  - [`SymbolToken`](@ref)
  - [`StopToken`](@ref)

See also: [`SourceLocation`](@ref)
"""
struct Token
     location::SourceLocation
     value::Union{KeywordToken, IdentifierToken, StringToken, LiteralNumberToken, SymbolToken,StopToken}
end

copy(token::KeywordToken) = KeywordToken(token.keyword)
copy(token::IdentifierToken) = IdentifierToken(token.identifier)
copy(token::StringToken) = StringToken(token.string)
copy(token::LiteralNumberToken) = LiteralNumberToken(token.number)
copy(token::SymbolToken) = SymbolToken(token.symbol)
copy(token::StopToken) = StopToken()
copy(token::Token) = Token(token.location, copy(token.value))

"""
     GrammarError <: Exception(
          location::SourceLocation
          message::String
     )

An error found by the lexer/parser while reading a scene file.

## Arguments

- `location::SourceLocation` : location of the last char read

- `message::String` : a user-frendly error message

See also: [`SourceLocation`](@ref)
"""
struct GrammarError <: Exception
    location::SourceLocation
    message::String
end


text:
using Gtk
import GeometryTypes, Cairo, Colors, Reactive, Chipmunk, DelimitedFiles
const gt = GeometryTypes
const ca = Cairo
const co = Colors
const re = Reactive
const df = DelimitedFiles

"""
    ClickObject

Struct for keeping the clicked object on screen and additional rendering
information.
The object is accessed by reference, so changes are visible the multiple owners.
"""
mutable struct ClickObject
    kind::Symbol
    object::Union{Object{<:Number}, Nothing}
end

mutable struct FunctionWrapper
    fun::Function
end

"""
get_world_to_screen_coeffs(width::Integer, height::Integer,
    field_size::AbstractVector{<:AbstractFloat},
    border_strip::AbstractVector{<:AbstractFloat})

Get the coefficints for world to screen and screen to world coordiate
transformation.

# Examples

```jldoctest
julia> get_world_to_screen_coeffs(900, 1500, [9.0, 6.0], [1.0, 1.0])
(81.81818181818181, [0.0, 422.727])
```
"""
function get_world_to_screen_coeffs(width::Integer, height::Integer,
    field_size::AbstractVector{<:AbstractFloat},
    border_strip::AbstractVector{<:AbstractFloat})

    @inbounds field_aspect_ratio =
        (field_size[2] + 2*border_strip[2]) /
        (field_size[1] + 2*border_strip[1])

    if width*field_aspect_ratio > height
        width_offset = (width*field_aspect_ratio - height)/2.0
        height_offset = 0.0
        @inbounds multiplier = height/(field_size[2] + 2*border_strip[2])#8.0
    else
        height_offset = (height - width*field_aspect_ratio)/2.0
        width_offset = 0.0
        @inbounds multiplier = width/(field_size[1] + 2*border_strip[1])#12.0
    end
    multiplier, [width_offset, height_offset]
end

"""
world_to_screen(x::AbstractVector{<:AbstractFloat},
    multiplier::AbstractFloat, offset::AbstractVector{<:AbstractFloat},
    xy_inversion_correction::AbstractVector{<:Real}=[0, 0],
    border_strip::AbstractVector{<:AbstractFloat})

Converts world to screen coordinates by respecting the field origin.
The origin of worldcoordinates is the lower left corner of the field.
Respects the inversion of coordinate axis, if the width or the hieght is set to
the `xy_inversion_correction` array ([width, height]). Axis, which should not be
inverted are set to zero. For example: [0, height] for an inversion of the y-Axis,
what is quite common for the most gui-frameworks.

# Examples

```jldoctest
julia> world_to_screen([1.0, 1.0], 75.0, [0.0, 450.0], [1.0, 1.0])
2-element Array{Float64,1}:
 150.0
 600.0
```
"""
@inline function world_to_screen(x::AbstractVector{<:AbstractFloat},
    center_offset::AbstractVector{<:AbstractFloat},
    border_strip::AbstractVector{<:AbstractFloat},
    multiplier::AbstractFloat,
    offset::AbstractVector{<:AbstractFloat},
    xy_inversion_correction::AbstractVector{<:Real}=[0, 0])

    abs.(xy_inversion_correction .- ((x .+ border_strip .+ center_offset).*multiplier .+ offset))
end

"""
screen_to_world(x::AbstractVector{<:AbstractFloat},
    multiplier::AbstractFloat, offset::AbstractVector{<:AbstractFloat},
    xy_inversion_correction::AbstractVector{<:Real}=[0, 0],
    border_strip::AbstractVector{<:AbstractFloat})

Converts screen to worl coordinates by respecting the field origin.
The origin of worldcoordinates is the lower left corner of the field.
Respects the inversion of coordinate axis, if the width or the hieght is set to
the `xy_inversion_correction` array ([width, height]). Axis, which should not be
inverted are set to zero. For example: [0, height] for an inversion of the y-Axis,
what is quite common for the most gui-frameworks.

# Examples

```jldoctest
julia> screen_to_world([112.5, 525.0], 75.0, [0.0, 450.0], [1.0, 1.0])
2-element Array{Float64,1}:
 0.5
 0.0
```
"""
@inline function screen_to_world(x::AbstractVector{<:AbstractFloat},
    center_offset::AbstractVector{<:AbstractFloat},
    border_strip::AbstractVector{<:AbstractFloat},
    multiplier::AbstractFloat,
    offset::AbstractVector{<:AbstractFloat},
    xy_inversion_correction::AbstractVector{<:Real}=[0, 0])

    ((abs.(xy_inversion_correction .- x) .- offset)./multiplier .- border_strip .- center_offset)
end

"""
    draw_player!(ctx, pos::Robot{<:Number}, radius::Real, col::co.RGBA{<:Number}, multiplier::Number, wts::Function)

Draws a Robot as a circle with radius `get_size(robot)` on its position on the screen. Adds a isosceles triangle for the body orientation with its head heaing in the direction of the angle.
"""
function draw_player!(ctx, robot::Robot{<:Number}, radius::Real, col::co.RGBA{<:Number}, multiplier::Number, wts::Function)
    arc_open = robot.hfov/2 # β/2
    add_r = radius/3
    x, y = wts(robot)
    angle = robot.body_orientation
    arrow_a = wts([cos(angle), sin(angle)] .* (radius + add_r) .+ robot)
    arrow_b = wts([cos(angle + arc_open), sin(angle + arc_open)] .* radius .+ robot)
    arrow_c = wts([cos(angle - arc_open), sin(angle - arc_open)] .* radius .+ robot)
    # ca.circle(ctx, x, y, radius*multiplier)
    # ca.set_source_rgba(ctx, col.r, col.g, col.b, col.alpha)

    sketch(ctx, GraphicLine([arrow_a, arrow_b, arrow_c, arrow_a], fill=true, fillcolor=[col.r, col.g, col.b, col.alpha]))
    draw_circle!(ctx, robot, radius, col, multiplier, wts)
end

"""
    draw_circle!(ctx, pos::AbstractVector{<:AbstractFloat}, radius::Real,
        col::co.RGBA{<:AbstractFloat}, multiplier::AbstractFloat, wts::Function)

# Attributes
- `ctx`: grphics context for drawing
- `pos::AbstractVector{<:Number}`: the postion for drawing the circle in world coordinates
- `radius:Real`: the radius of the circle in world
- `col::co.RGBA{<:Number}`: color in rgba (0..1) format
- `multiplier::Number, offset::AbstractVector{<:Number}`: render coefficients, returned by the function `get_world_to_screen_coeffs`
- `wts::Function`: function for converting world coordinates to screen coordinates (must have exact one argument of type AbstractVector{<:Real})
(for drawing on canvas, the y axis is inverted respectivly to coordinatesystems used in math)
"""
function draw_circle!(ctx, pos::AbstractVector{<:Number}, radius::Real, col::co.RGBA{<:Number}, multiplier::Number, wts::Function)
    x, y = wts(pos)
    ca.circle(ctx, x, y, radius*multiplier)
    ca.set_source_rgba(ctx, col.r, col.g, col.b, col.alpha)
    ca.fill(ctx)
end

"""
    draw_and_save_pathlines!(ctx, cur_position::AbstractVector{<:AbstractFloat},
        pathline::AbstractVector{<:AbstractVector{<:AbstractFloat}},
        col::co.RGBA{<:AbstractFloat}, multiplier::AbstractFloat,
        offset::AbstractVector{<:AbstractFloat}, inv_y::Function,
        atol::Real=0.01, size::Real=1.0)

# Attributes
- `ctx`: graphics context for drawing
- `cur_position::AbstractVector{<:Number}`: the current position of the tracked object
- `pathline::AbstractVector{<:AbstractVector{<:Number}}`: vector for storing ll path line points
- `col::co.RGBA{<:Number}`: color in rgba (0..1) format
- `multiplier::Number, offset::AbstractVector{<:Number}`: render coefficients, returned by the function `get_world_to_screen_coeffs`
- `wts::Function`: function for converting world coordinates to screen coordinates (must have exact one argument of type AbstractVector{<:Real})
(for drawing on canvas, the y axis is inverted respectively to coordinatesystems used in math)
- `atol::Real=0.01`: tolerance for storing points to the path line, in m - default: only if `cur_postition` is more than 1 cm away from the last added position
- `size::Real=0.01`: size of the drawn points (drawn as circle) in m
"""
function draw_and_save_pathlines!(ctx,
    cur_position::AbstractVector{<:AbstractFloat},
    pathline::AbstractVector{<:AbstractVector{<:AbstractFloat}},
    col::co.RGBA{<:AbstractFloat}, multiplier::AbstractFloat,
    wts::Function, atol::Real=0.012, size::Real=0.01)

    if isempty(pathline) || !isapprox(pathline[end], cur_position, atol=atol)
        push!(pathline, copy(cur_position))
        if length(pathline) == 10_000
            @warn "More than 10_000 points saved for path lines."
        end
    end
    for line in pathline
        draw_circle!(ctx, line, size, col, multiplier, wts)
    end
end

using Images

mutable struct App
    config::Config
    background_cache::Dict{String, ca.CairoSurface{co.ColorTypes.RGB24}}

    main_window::GtkWindow
    root::Gtk.GtkWidget
    canvas::GtkCanvas
    status_bar::GtkStatusbar
    progress_bar::GtkProgressBar
    ball_check_button::GtkCheckButton

    myPlayer_pathline_signal::AbstractVector{re.Signal{Bool}}
    opponent_pathline_signal::AbstractVector{re.Signal{Bool}}
    obstacle_pathline_signal::AbstractVector{re.Signal{Bool}}

    center_offset::AbstractVector{<:Number} #TODO: move to config
    border_strip::AbstractVector{<:Number} #TODO: move to config

    App() = new()
end

"""
    Create a canvas for a football robot simulator app.
    Here the canvas is used to display the simulation of a strategy.
"""
function create_canvas(app::App)
    ws = app.config.ws
    field_size = app.config.field_size

    canvas = GtkCanvas()
    set_gtk_property!(canvas, :expand, true)
    set_gtk_property!(canvas, :width_request, 390)

    # Initialize arrays for pathlines
    # TODO: Move to Linked List from DataStructurres.jl for performance
    myPlayer_pathlines = fill(Vector{Float64}[], length(ws.myPlayers))
    opponent_pathlines = fill(Vector{Float64}[], length(ws.opponents))
    obstacle_pathlines = fill(Vector{Float64}[], length(ws.obstacles))

    # Register functions for deleting pathlines
    del_off_pathlines = [map(!, sig) for sig in app.myPlayer_pathline_signal]
    foreach(enumerate(del_off_pathlines)) do element
        idx, signal = element
        re.foreach(signal) do sig
            empty!(myPlayer_pathlines[idx])
            draw(canvas)
        end
    end
    del_opp_pathlines = [map(!, sig) for sig in app.opponent_pathline_signal]
    foreach(enumerate(del_opp_pathlines)) do element
        idx, signal = element
        re.foreach(signal) do sig
            empty!(opponent_pathlines[idx])
            draw(canvas)
        end
    end
    del_obs_pathlines = [map(!, sig) for sig in app.obstacle_pathline_signal]
    foreach(enumerate(del_obs_pathlines)) do element
        idx, signal = element
        re.foreach(signal) do sig
            empty!(obstacle_pathlines[idx])
            draw(canvas)
        end
    end

    @guarded draw(canvas) do widget
        # Initialization
        ctx = Gtk.getgc(canvas)
        h = Gtk.height(canvas)
        w = Gtk.width(canvas)
        ca.set_source_rgb(ctx, 1.0, 1.0, 1.0)
        ca.rectangle(ctx, 0.0, 0.0, w, h)
        ca.fill(ctx)
        ca.paint(ctx)

        multiplier, offset = get_world_to_screen_coeffs(w, h, field_size, app.border_strip)
        wts(pos::AbstractVector{<:Real}) = world_to_screen(pos, app.center_offset, app.border_strip, multiplier, offset, [0, h])

        # Draw all config.drawables, which are images and shall be below the other stuff
        for (name, (pos, fun)) in app.config.drawables
            if startswith(name, "_")
                @inbounds screen_tl_pos = wts(pos[1])
                @inbounds screent_br_pos = wts(pos[2])
                # screen_size = wts(pos[2], multiplier, offset) .- screen_tl_pos
                screen_size = screent_br_pos .- screen_tl_pos
                preimage, redraw = fun()
                if !haskey(app.background_cache, name) || redraw
                    @inbounds image = @view preimage[:, end:-1:1]
                    im_enlarged = imresize(image, trunc.(Int, screen_size)...)
                    if any(isnan.(im_enlarged))
                        @warn "resized image: isnan"
                        im_enlarged[isnan.(im_enlarged)] = co.ColorTypes.RGB(0.0, 0.0, 0.0)
                    end
                    app.background_cache[name] = ca.CairoImageSurface(co.RGB24.(im_enlarged))
                end
                ca.save(ctx)
                ca.rectangle(ctx, screen_tl_pos..., screen_size...)
                ca.clip(ctx)
                ca.new_path(ctx)
                ca.set_source(ctx, app.background_cache[name], screen_tl_pos...)
                ca.paint(ctx)
                ca.restore(ctx)
            end
        end

        # Draw field border
        ca.set_line_width(ctx,2)
        @inbounds tl_rect = wts([-field_size[1]/2, field_size[2]/2])
        size_rect = world_to_screen(field_size./2, app.center_offset, app.border_strip, multiplier, offset) .- tl_rect
        @inbounds ca.rectangle(ctx, tl_rect[1], tl_rect[2], size_rect[1], size_rect[2])
        ca.set_source_rgb(ctx, 0.0, 0.0, 0.0)
        ca.stroke(ctx)

        # Draw pathlines
        for idx in 1:length(ws.myPlayers)
            if re.value(app.myPlayer_pathline_signal[idx])
                @inbounds draw_and_save_pathlines!(ctx, ws.myPlayers[idx],
                    myPlayer_pathlines[idx], co.RGBA(0.1, 0.4, 1.0, 1.0),
                    multiplier, wts)
            end
        end
        for idx in 1:length(ws.opponents)
            if re.value(app.opponent_pathline_signal[idx])
                @inbounds draw_and_save_pathlines!(ctx, ws.opponents[idx],
                    opponent_pathlines[idx], co.RGBA(0.1, 0.1, 0.1, 1.0),
                    multiplier, wts)
            end
        end
        for idx in 1:length(ws.obstacles)
            if re.value(app.obstacle_pathline_signal[idx])
                @inbounds draw_and_save_pathlines!(ctx, ws.obstacles[idx],
                    obstacle_pathlines[idx], co.RGBA(ws.obstacles[idx].color...), multiplier, wts)
            end
        end

        # Draw player and ball
        for (idx, myPlayer) in enumerate(ws.myPlayers)
            draw_player!(ctx, myPlayer, get_size(myPlayer),
                co.RGBA(0.1, 0.4, 1.0, 1.0), multiplier, wts)
            sketch(ctx, GraphicText(wts(myPlayer),
                string(idx), textsize=15.0))
        end
        for (idx, opponent) in enumerate(ws.opponents)
            draw_player!(ctx, opponent, get_size(opponent),
                co.RGBA(0.1, 0.1, 0.1, 1.0), multiplier, wts)
            sketch(ctx, GraphicText(wts(opponent),
                string(idx), textsize=15.0))
        end
        for (idx, obstacle) in enumerate(ws.obstacles)
            !(typeof(obstacle.shape) <: Chipmunk.CircleShape) && continue
            draw_circle!(ctx, obstacle, get_size(obstacle),
                co.RGBA(obstacle.color...), multiplier, wts)
            sketch(ctx, GraphicText(wts(obstacle),
                string(idx), textsize=get_size(obstacle)*180))
        end
        # ws.isballwhere && draw_circle!(ctx, ws.ball, ws.ball_size, co.RGBA(1.0, 0.1, 0.1, 1.0),
            # multiplier, wts)

        # Draw config.drawables from GameState
        for (name, (pos, fun)) in app.config.drawables
            screen_pos = wts.(pos)
            if startswith(name, "_")
                continue
            end

            # If the element is too small, there will be an error
            @inbounds if length(screen_pos) == 2 && screen_pos[1] ≈ screen_pos[2]
                continue
            end
            sketch(ctx, fun(screen_pos))
        end
    end

    clicked_object = ClickObject(:nothing, nothing)
    stwWrap = FunctionWrapper(() -> nothing)

    function click_player!(clicked_object::ClickObject, player::Object{<:Number}, event_pos::gt.Point2)
        if is_inside(event_pos, gt.Circle(gt.Point2(player), get_size(player)))
            clicked_object.kind = :position
            clicked_object.object = player
            return true
        end
        norm_x, norm_y = event_pos .- player
        if typeof(player) <: Robot && get_size(player) < distance(event_pos, player) ≤ get_size(player)*1.3 # HACK: Magic number
            angle_diff = abs(player.body_orientation - atan(norm_y, norm_x))
            if angle_diff < π/4 || angle_diff > 15π/8
                clicked_object.kind = :angle
                clicked_object.object = player
                return true
            end
        end
        clicked_object.kind = :nothing
        false
    end

    canvas.mouse.button1press = @guarded (widget, event) -> begin
        pos = gt.Point2(stwWrap.fun([event.x, event.y]))
        for player in vcat(ws.myPlayers, ws.opponents, ws.obstacles)
            click_player!(clicked_object, player, pos) && return
        end
        return
    end

    function set_pos!(clicked_object::ClickObject, widget, event, kind::Symbol)
        pos = stwWrap.fun([event.x, event.y])
        if clicked_object.kind == :angle
            x, y = pos .- clicked_object.object
            angle = atan(y, x)
            clicked_object.object.body_orientation = angle
        elseif clicked_object.kind == :position
            set_position!(clicked_object.object, physic(ws, clicked_object.object, pos))
        else
            return
        end
        get_strategy(app.config).run_strategies(app.config, 0, kind, clicked_object.object)
        draw(widget)
    end

    canvas.mouse.button1motion = @guarded (widget, event) -> begin
        clicked_object.kind == :nothing && return
        # empty!(app.background_cache) # HACK: only no new drawing of cached
        set_pos!(clicked_object, widget, event, :motion)
    end

    canvas.mouse.button1release = @guarded (widget, event) -> begin
        clicked_object.kind == :nothing && return
        empty!(app.background_cache)
        set_pos!(clicked_object, widget, event, :release)
        clicked_object.kind = :nothing
    end

    canvas.mouse.button3release = @guarded (widget, event) -> begin
        pos = gt.Point2(stwWrap.fun([event.x, event.y]))
        if length(ws.myPlayers) ≥ 1
            for (idx, off) in enumerate(ws.myPlayers)
                if click_player!(clicked_object, off, pos)
                    remove_object!(app.config, ws.myPlayers[idx])
                    deleteat!(ws.myPlayers, idx)
                    # ws.player = vcat(ws.myPlayers, ws.opponents)
                    # ws.all = vcat(ws.myPlayers, ws.opponents, ws.obstacles)
                    break
                end
            end
        end
        if clicked_object.kind == :nothing && length(ws.opponents) ≥ 1
            for (idx, opp) in enumerate(ws.opponents)
                if click_player!(clicked_object, opp, pos)
                    remove_object!(app.config, ws.opponents[idx])
                    deleteat!(ws.opponents, idx)
                    # ws.player = vcat(ws.myPlayers, ws.opponents)
                    # ws.all = vcat(ws.myPlayers, ws.opponents, ws.obstacles)
                    break
                end
            end
        end
        if clicked_object.kind == :nothing && length(ws.obstacles) ≥ 1
            for (idx, obs) in enumerate(ws.obstacles)
                if click_player!(clicked_object, obs, pos)
                    remove_object!(app.config, ws.obstacles[idx])
                    deleteat!(ws.obstacles, idx)
                    # ws.player = vcat(ws.myPlayers, ws.opponents)
                    # ws.all = vcat(ws.myPlayers, ws.opponents, ws.obstacles)
                    break
                end
            end
        end
        clicked_object.kind == :nothing && return
        empty!(app.background_cache)
        clicked_object.kind = :nothing
        recreate_root!(app::App)
    end

    canvas.mouse.motion = @guarded (widget, event) -> begin
        x, y = stwWrap.fun([event.x, event.y])
        push!(app.status_bar, 0,
            @sprintf("Position: %.2f, %.2f - Screen: %.2f, %.2f",
                x, y, event.x, event.y))
    end

    # signal_connect(canvas_frame, "resize") do widget, width, height
    canvas.resize = @guarded (widget) -> begin
        h = Gtk.height(widget)
        w = Gtk.width(widget)
        multiplier, offset = get_world_to_screen_coeffs(w, h, field_size, app.border_strip)
        stwWrap.fun = x -> screen_to_world(x, app.center_offset, app.border_strip, multiplier, offset, [0, h])
    end

    return canvas
end

"""
    Create a strategy menu for a football robot simulator app.
"""
function create_strategy_menu(app::App)
    print_grid = GtkGrid()

    for (idx, p) in enumerate(app.config.printables)
        name, sig = p
        l1 = GtkLabel(string(name,": "))
        l2 = GtkLabel(string(re.value(sig)))
        foreach(sig) do s
            set_gtk_property!(l2, :label, s)
        end
        print_grid[1, idx] = l1
        print_grid[2, idx] = l2
    end

    manipulator_menu = GtkBox(:v)
    scrolled_manip_menu = GtkScrolledWindow(manipulator_menu)
    set_gtk_property!(scrolled_manip_menu, :vexpand, true)

    # Register all manipulators from the strategy
    for (name, widget) in app.config.manipulators
        widget, _, fun = widget
        if typeof(widget) <: Gtk.GtkRadioButtonGroupLeaf && length(widget) < 2
            continue
        end
        signal_connect(widget, "clicked") do w
            fun(w)
            empty!(app.background_cache)
            @async draw(app.canvas)
            return
        end
        if typeof(widget) <: Gtk.GtkRadioButtonLeaf
            continue
        end
        push!(manipulator_menu, GtkLabel(name))
        push!(manipulator_menu, widget)
    end

    result = GtkBox(:v)
    set_gtk_property!(result, :margin, 5)
    set_gtk_property!(result, :width_request, 150)
    push!(result, print_grid)
    push!(result, scrolled_manip_menu)

    return result
end

"""
    Create a general menu for a football robot simulator app.
"""
function create_general_menu(app::App)
    config = app.config
    ws = config.ws

    # Create stratgies combobox
    strategies_combobox = GtkComboBoxText()
    strategies_list = collect(keys(config.strategies_table))
    for s in strategies_list
        push!(strategies_combobox, s)
    end
    set_gtk_property!(strategies_combobox, :active,
        argmax(strategies_list.==config.active_strategy) - 1)

    # # Save worldstate for resetting
    # initial_ws = deepcopy(ws)

    # Create control elements
    config_filechooser_button = GtkButton("Load config")
    play_button = GtkToggleButton("Play Strategy")
    fps_label = GtkLabel("Frames per s")
    speed_scalebutton = GtkScale(false, 1:60)
    # velocity_scale_label = GtkLabel("Velocity scale in m/s")
    speed_adjustment = GtkAdjustment(speed_scalebutton)
    set_gtk_property!(speed_adjustment, :value, 30)
    # velocity_scalebutton = GtkScale(false, 0.05:0.05:1.0)
    # velocity_adjustment = GtkAdjustment(velocity_scalebutton)
    # set_gtk_property!(velocity_adjustment, :value, 0.1)
    step_button = GtkButton("Step")
    # reset_button = GtkButton("Reset")
    addmyPlayer_button = GtkButton("Add myPlayer")
    addopponent_button = GtkButton("Add opponent")
    pathlines_label = GtkLabel("Pathlines")

    set_gtk_property!(fps_label, :halign, Gtk.GConstants.GtkAlign.GTK_ALIGN_START)
    # set_gtk_property!(velocity_scale_label,
    #     :halign, Gtk.GConstants.GtkAlign.GTK_ALIGN_START)
    set_gtk_property!(pathlines_label,
        :halign, Gtk.GConstants.GtkAlign.GTK_ALIGN_START)

    # Create menu_box
    menu_box = GtkBox(:v)
    set_gtk_property!(menu_box, :margin, 5)
    set_gtk_property!(menu_box, :spacing, 5)

    push!(menu_box, strategies_combobox)
    push!(menu_box, config_filechooser_button)
    push!(menu_box, play_button)
    push!(menu_box, fps_label)
    push!(menu_box, speed_scalebutton)
    # push!(menu_box, velocity_scale_label)
    # push!(menu_box, velocity_scalebutton)
    push!(menu_box, step_button)
    # push!(menu_box, reset_button)
    push!(menu_box, addmyPlayer_button)
    push!(menu_box, addopponent_button)
    push!(menu_box, pathlines_label)

    pathline_menu_box = GtkBox(:v)
    set_gtk_property!(pathline_menu_box, :vexpand, true)

    # Register checkboxes for pathlines and emit signals
    app.myPlayer_pathline_signal = fill(re.Signal(false), length(ws.myPlayers))
    myPlayer_check_buttons = [GtkCheckButton("myPlayer $idx")
    for idx in 1:length(ws.myPlayers)]

    [signal_connect(off_check_button, "clicked") do widget
        push!(app.myPlayer_pathline_signal[idx],
            get_gtk_property(widget, :active, Bool))
        end
        for (idx, off_check_button) in enumerate(myPlayer_check_buttons)]
    append!(pathline_menu_box, myPlayer_check_buttons)

    app.opponent_pathline_signal = fill(re.Signal(false), length(ws.opponents))
    opponent_check_buttons =
        [GtkCheckButton("Opponent $idx") for idx in 1:length(ws.opponents)]
    [signal_connect(opp_check_button, "clicked") do widget
        push!(app.opponent_pathline_signal[idx],
            get_gtk_property(widget, :active, Bool))
        end
        for (idx, opp_check_button) in enumerate(opponent_check_buttons)]
    append!(pathline_menu_box, opponent_check_buttons)

    app.obstacle_pathline_signal = fill(re.Signal(false), length(ws.obstacles))
    obstacle_check_buttons =
        [GtkCheckButton("Obstacle $idx") for idx in 1:length(ws.obstacles)]
    [signal_connect(obs_check_button, "clicked") do widget
        push!(app.obstacle_pathline_signal[idx],
            get_gtk_property(widget, :active, Bool))
        end
        for (idx, obs_check_button) in enumerate(obstacle_check_buttons)]
    append!(pathline_menu_box, obstacle_check_buttons)

    push!(menu_box, GtkScrolledWindow(pathline_menu_box))

    signal_connect(config_filechooser_button, "clicked") do widget
        config_file_path = open_dialog("Choose config file.", app.main_window)
        if typeof(config_file_path)<:AbstractArray || config_file_path == ""
            return
        end
        empty!(app.background_cache)
        app.config = get_config(String["--file", config_file_path])
        recreate_root!(app::App)
    end

    signal_connect(strategies_combobox, "changed") do widget, others...
        idx = get_gtk_property(strategies_combobox, "active", Int)
        config.active_strategy = Gtk.bytestring(
            GAccessor.active_text(strategies_combobox))
        empty!(app.background_cache)
        recreate_root!(app::App)
    end

    # Perform animation
    animate = re.Signal(false)
    signal_connect(play_button, "clicked") do widget
        active = get_gtk_property(widget, :active, Bool)
        push!(animate, active)
        set_gtk_property!(widget, :label, active ? "Stop Strategy" : "Play Strategy")
    end
    animation_speed = re.Signal(get_gtk_property(speed_adjustment, :value, Int))
    signal_connect(speed_adjustment, "value-changed") do widget
        push!(animation_speed, get_gtk_property(widget, :value, Int))
    end
    animation = re.fpswhen(animate, animation_speed)

    function one_step(dt::Number)
        (dt ≤ 0.0 || dt > 2/get_gtk_property(speed_adjustment, :value, Float64)) && return
        # dt ≤ 0.0 && return
        update_game!(ws, dt, app.config)
        # physic_step(config, dt)
        @async draw(app.canvas)
        return
    end
    signal_connect(step_button, "clicked") do widget
        one_step(1/get_gtk_property(speed_adjustment, :value, Float64))
    end
    re.foreach(one_step, animation, init=nothing)

    # # Reset_button world
    # signal_connect(reset_button, "clicked") do widget
    #     empty!(app.background_cache)
    #     app.config.ws = initial_ws
    #     recreate_root!(app)
    # end

    # addmyPlayer_button world
    signal_connect(addmyPlayer_button, "clicked") do widget
        config.maximum_player_per_side ≤ length(config.ws.myPlayers) && return
        empty!(app.background_cache)
        add_myPlayer!(config)
        recreate_root!(app)
        # draw(app.canvas)
    end

    # addopponent_button world
    signal_connect(addopponent_button, "clicked") do widget
        config.maximum_player_per_side ≤ length(config.ws.opponents) && return
        empty!(app.background_cache)
        add_opponent!(config)
        recreate_root!(app)
    end

    # Initialize everything
    get_strategy(config).init_strategies(config)
    # get_strategy(config).init_strategies(config, x=reshape(df.readdlm("run01.txt"), (12)))
    # get_strategy(config).run_strategies(config)

    return menu_box
end

"""
Create and fill the root container of the window.
"""
function create_root!(app::App)

    empty!(app.config.drawables)
    empty!(app.config.printables)
    empty!(app.config.manipulators)

    left_menu = create_general_menu(app)
    app.canvas = create_canvas(app)
    right_menu = create_strategy_menu(app)

    # Setup main view
    main_view = GtkBox(:h)
    push!(main_view, left_menu)
    push!(main_view, app.canvas)
    push!(main_view, right_menu)

    # Setup status bar
    app.status_bar = GtkStatusbar()

    # app.progress_bar = ProgressBar()

    # function set_progressbar(fraction::Real)
    #     @debug fraction
    #     set_gtk_property!(app.progress_bar, :fraction, fraction)
    # end

    # app.config.set_progress_bar = set_progressbar

    # Setup root
    root = GtkBox(:v)
    push!(root, main_view)
    push!(root, app.status_bar)
    # push!(root, app.progress_bar)

    root
end

"""
Recreates the window root after changes.
"""
function recreate_root!(app::App)
    destroy(app.root)
    app.root = create_root!(app)
    push!(app.main_window, app.root)
    showall(app.main_window)
    # visible(app.ball_check_button, app.config.ws.isballwhere)
end

"""
    Create and start a new football robot simulator app.
"""
function start_view(config::Config)
    app = App()
    app.config = config
    app.border_strip = [0.7, 0.7] # HACK: magic numer
    app.center_offset = config.field_size./2 # HACK: magic numer
    app.background_cache = Dict{String, Any}()

    # Setup main window
    app.main_window = GtkWindow("Football Robot Simulator", 1440, 900) # HACK: magic numer
    app.root = create_root!(app)

    push!(app.main_window, app.root)

    showall(app.main_window)

    on_destroy(windowLeaf::Gtk.GtkWindowLeaf) = Gtk.gtk_quit()
    signal_connect(on_destroy, app.main_window, :destroy)

    # visible(app.ball_check_button, config.ws.isballwhere)

    Gtk.gtk_main()
end

ext: jl


text:
include("log.jl")
mutable struct Route{T}
    path::String
    page::T
    function Route(path::String = "", page::Function = http -> "")
        new{Function}(path, page)
    end
    function Route(path::String, page::Page)
        new{Page}(path, page)
    end
    function Route(path::String, page::FormComponent)
        new{FormComponent}(path, page)
    end
    function Route(path::String, page::Component)
        new{typeof(page)}(path, page)
    end
end

function route_from_dir(dir::String)
    dirs = readdir(dir)
    routes::Vector{String} = []
    for directory in dirs
        if isfile("$dir/" * directory)
            push!(routes, "$dir/$directory")
        else
            if ~(directory in routes)
                newread = dir * "/$directory"
                newrs = route_from_dir(newread)
                [push!(routes, r) for r in newrs]
            end
        end
    end
    rts::Vector{Route} = []
    for directory in routes
        if isfile("$dir/" * directory)
            push!(rts, Route("/$directory", file("$dir/" * directory)))
        end
    end
    rts
end

mutable struct ServerTemplate
    ip::String
    port::Integer
    routes::Vector{Route}
    logger::Logger
    remove::Function
    add::Function
    start::Function
    public::String
    function ServerTemplate(ip::String = "127.0.0.1", port::Int64 = 8001,
        routes::Vector{Route} = Vector{Route}(); logger::Logger = Logger(),
        public::String = "public")
        add, remove, start = funcdefs(routes, ip, port, logger, public)
        new(ip, port, routes, logger, remove, add, start, public)
    end
end

function funcdefs(routes::AbstractVector, ip::String, port::Integer,
    logger::Logger, public::String)
    add(r::Route{Function}) = push!(routes, r)
    add(r::Route{Component}) = push!(routes, r)
    add(r::Route{Page}) = begin push!(routes, r)
        for comp in r.page.components
            if typeof(comp) != Function
                if typeof(comp) <: FormComponent
                    push!(routes, Route(comp.action, fn(comp.onAction)))
                end
            end
        end
    end
    add(r::Route{FormComponent}) = begin push!(routes, r);
        push!(routes, Route(r.page.action, fn(r.page.onAction)))
    end
    remove(i::Int64) = deleteat!(routes, i)
    start() = _start(routes, ip, port, logger, public)
    return(add, remove, start)
end

function _start(routes::AbstractVector, ip::String, port::Integer,
    logger::Logger, public::String)
    public_rs = route_from_dir(public)
    routes = vcat(routes, public_rs)
    server = Sockets.listen(Sockets.InetAddr(parse(IPAddr, ip), port))
    logger.log(1, "Toolips Server starting on port " * string(port))
    routefunc = generate_router(routes, server, logger)
    @async HTTP.listen(routefunc, ip, port; server = server)
    logger.log(2, "Successfully started server on port " * string(port))
    logger.log(1,
    "You may visit it now at http://" * string(ip) * ":" * string(port))
    return(server)
end
function generate_router(routes::AbstractVector, server, logger::Logger)
    route_paths = Dict([route.path => route.page for route in routes])
    # CORE routing server lies here.
    # - Router itself is merely a function that gets called with the http
    #  stream. This trickles down the line all the way to the interface methods.
    routeserver = function serve(http)
    HTTP.setheader(http, "Content-Type" => "text/html")
    fullpath = http.message.target
    # Checks for argument data, because this is not in the route.
    if contains(http.message.target, '?')
         fullpath = split(http.message.target, '?')[1]
    end

     if fullpath in keys(route_paths)
         if typeof(route_paths[fullpath]) == Function
             write(http, route_paths[fullpath](http))
         else
             write(http, route_paths[fullpath].f(http))
         end
     else
         if typeof(route_paths["404"]) != Page
             write(http, route_paths["404"](http))
         else
            write(http, route_paths["404"].f(http))
        end
     end

 end # serve()
    return(routeserver)
end
function stop!(x::Any)
    close(x)
end


text:
#!/usr/bin/env julia

# License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

__precompile__()

module FidelityUMI

export
    ReferenceSequences,
    Coverage,
    Mutation,
    filter_variants,
    call_variants,
    check_operations,
    check_quality,
    collect_variants,
    count_mut_coverages!,
    coverage_counter!,
    dedup_mutations,
    get_read_shift,
    populate_variants!,
    translate_mutations,
    compare_seq,
    mapp_mutations,
    parse_chr,
    stringify,
    conf_coverage!,
    transvers,
    MutationsCluster,
    mapp_single_aa_subs

import BioAlignments: BioAlignments, BAM, Operation
import BioAlignments: OP_MATCH, OP_INSERT, OP_DELETE, OP_SKIP
import BioAlignments: OP_SOFT_CLIP, OP_HARD_CLIP, OP_PAD
import BioAlignments: OP_SEQ_MATCH, OP_SEQ_MISMATCH, OP_BACK, OP_START
import BioSequences: DNASequence, FASTA, ReferenceSequence, BioSequences
import BioSequences: sequence, seqname, RNASequence, translate
import AutoHashEquals: @auto_hash_equals
import CSV: CSV, eachrow
import DataFrames: DataFrame

type  ReferenceSequences
    file::String
    chromosomes::Array{String,1}
    chrom_lengths::Dict{String,Int64}
    sequences::Dict{String,BioSequences.ReferenceSequence}

    function ReferenceSequences(file::String)

        sequences::Dict{String,BioSequences.ReferenceSequence} = Dict()
        chromosomes::Array{String,1} = []
        chrom_lengths::Dict{String,Int64} = Dict()
        ext = split(file, ".")[end]

        if ext == "fasta" || ext == "fa" || ext == "fna"
            reader = FASTA.Reader(open(file,"r"))

        elseif ext == "2bit"
            reader = TwoBit.Reader(open(file, "r"))

        else
            println(STDERR, "ERROR: Reference FILE extention is not recognised: $ext")
            ext(1)
        end

        for record in reader
            sequences[seqname(record)] = convert(DNASequence,replace(string(sequence(record)),r"[^ATGCN]","N"))
            push!(chromosomes,seqname(record))
            chrom_lengths[seqname(record)] = length(sequences[seqname(record)])
        end

        return new(file, chromosomes, chrom_lengths, sequences)::ReferenceSequences
    end
end

type  Coverage
    chromosomes::Array{String,1}
    chrom_lengths::Dict{String,Int64}
    coverage::Dict{String,Array{UInt32,1}}

    function Coverage(reference::ReferenceSequences)
        chromosomes::Array{String,1} = reference.chromosomes
        chrom_lengths::Dict{String,Int64} = reference.chrom_lengths
        coverage::Dict{String,Array{UInt32,1}} = Dict()

        for chr in chromosomes
            coverage[chr] = zeros(UInt32, chrom_lengths[chr])
        end

        return new(chromosomes, chrom_lengths, coverage)
    end
end

@auto_hash_equals immutable Mutation
    class::Char
    chrom::String
    mStart::Int64
    mEnd::Int64
    reference::String
    variant::String
end

@auto_hash_equals mutable struct MutationsCluster
    mutation_count::Dict{Mutation,Int64}
    cluster_size::Int64
    R1_length::Int64
    R2_length::Int64
    total_mutations::Int64
    R1_coverage::Array{UInt16,1}
    R2_coverage::Array{UInt16,1}
    c_start::Int64
    c_end::Int64
end

function coverage_counter!(Coverage::Coverage, bamF::String; mappingQuality::Int64=0)
    reader = open(BAM.Reader, bamF)
    record = BAM.Record()
    ct::Int64 = 0
    ct2::Int64 = 0
    while !eof(reader)
        ct += 1
        read!(reader, record)

        if BAM.ismapped(record)
            refName::String = BAM.refname(record)
            flag = BAM.flag(record)
            mpq = BAM.mappingquality(record)
            # these conditions should match collect_variants() conditions
            if mpq >= mappingQuality && (flag&4 == 0) && (flag&256 == 0) && (flag&2048 == 0) && (refName in Coverage.chromosomes)   # filter out only mapped and primary alignments
                Coverage.coverage[refName][BAM.position(record):BAM.rightposition(record)]+=UInt32(1)
                ct2 += 1
            end
        end
    end
    close(reader)
end

function call_variants(record::BioAlignments.BAM.Record; quality::Int64=0, verbose::Bool=false)

    mutations::Array{Mutation,1} = []
    leftposition = BAM.position(record)
    md_tag = uppercase(record["MD"]::String)
    record_cigar = BAM.cigar(record)
    rfname = BAM.refname(record)
    read_seq = BAM.sequence(record)
    a_operations, o_lengths = BAM.cigar_rle(record)
    variants, md_positions = split_md_tag(md_tag)
    qualities = BAM.quality(record)

    if !(check_operations(a_operations))
        println(STDOUT, "INFO: CALLING VARIANTS, IN RECORD:")
        println(STDOUT, record)
        println()
        return Mutation[]
    end

    if verbose
        println(STDOUT, "INFO: CALLING VARIANTS, MD_TAG: ", md_tag)
        println(STDOUT, "INFO: CALLING VARIANTS, PARSED VARIANTS AND POSITIONS: ", variants, ":", md_positions)
        println(STDOUT, "INFO: CALLING VARIANTS, CIGAR RECORD: ", record_cigar)
        println(STDOUT, "INFO: CALLING VARIANTS, ALIGNMENT OPERATIONS: ", a_operations)
        println(STDOUT, "INFO: CALLING VARIANTS, LENGTHS OF OPERATIONS: ", o_lengths)
    end

    # Find variants and deletions in the read.
    ct = 0
    for i in variants
        ct += 1
        ref_shift = get_ref_shift(variants, md_positions, ct)
        read_shift = get_read_shift(a_operations, o_lengths, variants, md_positions, ct)

        if verbose
            println(STDOUT, "INFO: CALL VARIANTS, REF_SHIFT: ", ref_shift)
            println(STDOUT, "INFO: CALL VARIANTS, READ_SHIFT: ", read_shift)
        end

        if i[1] == '^'
                class = 'D'
            ix = 2
            s2 = length(i[2:end]) - 1
        else
            class = 'M'
            ix = 1
            s2 = 0
        end

        var_start_pos = leftposition + ref_shift - s2 - 1
        var_end_pos = leftposition + ref_shift - 1

        if verbose
            println(STDOUT, "INFO: CALL VARIANTS, READ LEFT MAPPING POSITION: ", leftposition)
            println(STDOUT, "INFO: CALL VARIANTS, CALCULATED SHIFT in REFERENCE: ", ref_shift)
            println(STDOUT, "INFO: CALL VARIANTS, VARIANT START POSITION IN REFERENCE: ", var_start_pos)
            println(STDOUT, "INFO: CALL VARIANTS, VARIANT END POSITION IN REFERENCE: ", var_end_pos)
            println(STDOUT, "INFO: CALL VARIANTS, QUALITY MAPPING START POSITION: ", read_shift)
            println(STDOUT, "INFO: CALL VARIANTS, QUALITY MAPPING end POSITION: ", read_shift + s2)
        end

        if i[1] == '^'
            if qualities[read_shift] >= quality
                push!(mutations, Mutation(class, rfname, var_start_pos, var_end_pos, i[ix:end], "*"))
            else
                if verbose
                    println(STDOUT, "SKIPED BY QUALITY: ", mean(qualities[read_shift:read_shift + s2]))
                end
            end
        else
            if mean(qualities[read_shift:read_shift]) >= quality
                v = read_seq[read_shift:read_shift]
                push!(mutations, Mutation(class, rfname, var_start_pos, var_end_pos, i[ix:end], convert(String, v)))
            else
                if verbose
                    println(STDOUT, "INFO: SKIPED BY QUALITY: ", mean(qualities[read_shift:read_shift]), " ", qualities[read_shift:read_shift])
                end
            end
        end
    end

    # find insertions in the read
    ct = 0
    read_pos = 0
    ref_pos = 0

    for i in a_operations
        ct += 1

        if i == BioAlignments.Operation(OP_MATCH)
            ref_pos += o_lengths[ct]
            read_pos += o_lengths[ct]

        elseif i == BioAlignments.Operation(OP_DELETE)
            ref_pos += o_lengths[ct]

        elseif i == BioAlignments.Operation(OP_SOFT_CLIP)
            read_pos += o_lengths[ct]

        elseif i == BioAlignments.Operation(OP_INSERT)
            in_start = read_pos + 1
            read_pos += o_lengths[ct]
            in_end = read_pos
            v = read_seq[in_start:in_end]
            var_start_pos = ref_pos + leftposition - 1
            var_end_pos = var_start_pos + 1

            if mean(qualities[in_start:in_end]) >= quality
                push!(mutations, Mutation('I', rfname, var_start_pos, var_end_pos, "-", v))
            end
        end
    end

    if verbose
        println(STDOUT, "INFO: VARIANTS CALCULATED: ", mutations)
        println(STDOUT, "INFO: FROM THIS RECORD: ", record)
        println()
    end

    return mutations
end

function check_operations(op::Array{BioAlignments.Operation,1})

    supported_operations = [BioAlignments.Operation(OP_MATCH),
                            BioAlignments.Operation(OP_INSERT),
                            BioAlignments.Operation(OP_SOFT_CLIP),
                            BioAlignments.Operation(OP_DELETE),
                            BioAlignments.Operation(OP_HARD_CLIP)]

    for i in op
        if !(i in supported_operations)
            println(STDOUT, "WARNING: FOUND NOT SUPPORTED ALIGMENT OPERATION: $i")
            return false
        end
    end

    return true
end

function split_md_tag(md::String)
    md_split::Array{String, 1} = split(replace(md, r"[\\^]*[ACGT]+", x -> "_"*x*"_"), "_")
    variants = String[]
    positions = Int64[]

    for i in md_split
        # finding numbers
        if ismatch(r"[0-9]+", i)
             pos = parse(Int64, i)
             push!(positions, pos)
        else
            push!(variants, i)
        end
    end

    return variants, positions
end

function get_read_shift(cigar_o::Array{BioAlignments.Operation,1},
                        lengths::Array{Int64,1},
                        variants::Array{String,1},
                        md_positions::Array{Int64,1},
                        pos::Int64
                        )
    shift = 0
    op_shitfs = 0
    ct = 0

    for i in variants
        ct += 1
        if pos >= ct && i[1] != '^'
            shift += md_positions[ct] + 1
            op_shitfs += md_positions[ct] + 1
        elseif pos >= ct
            shift += md_positions[ct]
            op_shitfs += md_positions[ct]
        end
    end

    ct = 0
    op_matches = 0
    for i in cigar_o
        ct += 1
        if i == BioAlignments.Operation(OP_MATCH)
            op_matches += lengths[ct]

        elseif i == BioAlignments.Operation(OP_SOFT_CLIP)
            shift += lengths[ct]

        elseif i == BioAlignments.Operation(OP_INSERT)
            shift += lengths[ct]
        end

        if  op_shitfs < op_matches
            break
        end
    end

    return shift
end

function get_ref_shift(variants::Array{String,1}, md_positions::Array{Int64,1}, pos::Int64)

    ct = 0
    shift = 0
    for i in md_positions
        ct += 1
        if pos >= ct
            shift += i
            if variants[ct][1] == '^'
                shift += length(variants[ct][2:end])
            else
                shift += 1
            end
        end
    end

    return shift
end

function collect_variants(bamfile::String, Coverage::Coverage, mappingq::Int64, mut_n_read::Int64, q_limit::Int64,
                          umi_s::String, umi_p::Int64, use_grouped::Bool, group_tag::String;
                          verbose::Bool=false, algo::String="umi")

    if !(algo in ["umi", "position", "paired", "naive"])
        println(STDERR, "ERROR: Selected --algorithm does not exist: $algo")
        println(STDERR, "ERROR: Options are: umi position paired naive")
        exit(1)
    end

    # collect coverages for each mutation from all clusters
    # TODO: Change mutcoverage structure to remove chromosome as key since Mutation has chromosome and is unique.
    mut_coverages = Dict{String,Dict{Mutation,Int64}}()
    # key is unique ID based either on mapping position or UMI + mapping position
    vars = Dict{String,MutationsCluster}()

    error_distr = Int64[]
    read_distr = Int64[]
    reader = open(BAM.Reader, bamfile)
    record = BAM.Record()
    record2 = BAM.Record()
    total_bases = 0
    total_reads = 0
    filtered_reads = 0
    filtered_R1 = 0
    filtered_R2 = 0

    while !eof(reader)
        read!(reader, record)

        try
            read!(reader, record2)
        catch
            break
        end

        if BAM.ismapped(record) && BAM.ismapped(record2)
            passed1 = false
            passed2 = false
            flag1 = BAM.flag(record)
            flag2 = BAM.flag(record2)
            tempname::String = BAM.tempname(record)
            tempname2::String = BAM.tempname(record2)
            rfname::String = BAM.refname(record)
            rfname2::String = BAM.refname(record2)
            mpq = BAM.mappingquality(record)
            mpq2 = BAM.mappingquality(record2)
            rpos = BAM.rightposition(record)
            lpos = BAM.position(record)
            rpos2 = BAM.rightposition(record2)
            lpos2 = BAM.position(record2)
            seq_len = BAM.seqlength(record)
            seq_len2 = BAM.seqlength(record2)
            mutations::Array{Mutation,1} = call_variants(record; quality=q_limit, verbose=verbose)
            mutations2::Array{Mutation,1} = call_variants(record2; quality=q_limit, verbose=verbose)
            l1 = length(mutations) > 0
            l2 = length(mutations2) > 0

            # This section is for counting mutations at positions
            # Conditions must match coverage_counter!() conditions
            if  mpq >= mappingq && flag1&4 == 0 && flag1&256 == 0 && flag1&2048 == 0 && l1 <= mut_n_read
                total_reads += 1
                passed1 = true
                if l1
                    # count before deduplication of overlaped paired regions
                    count_mut_coverages!(mut_coverages, mutations, rfname; verbose=verbose)
                    error_dist_over_length!(error_distr, mutations, lpos, seq_len)
                    read_dist!(read_distr, seq_len)
                end
            end

            if mpq2 >= mappingq && flag2&4 == 0 && flag2&256 == 0 && flag2&2048 == 0 && l2 <= mut_n_read
                total_reads += 1
                passed2 = true
                if l2
                    # count before deduplication of overlaped paired regions
                    count_mut_coverages!(mut_coverages, mutations2, rfname2; verbose=verbose)
                    error_dist_over_length!(error_distr, mutations2, lpos2, seq_len2)
                    read_dist!(read_distr, seq_len2)
                end
            end

            # only if first is R1 and second is R2
            # only if both mapped and primary alignments
            if tempname == tempname2 && passed1 && passed2 && rfname == rfname2
                if use_grouped && algo == "umi"
                    try
                        if flag1&64 == 64
                            UMI = uppercase(record[group_tag]::String)
                        else
                            UMI = uppercase(record2[group_tag]::String)
                        end
                    catch e
                        println(STDERR, "WARNING: GROUP TAG WAS NOT FOUND: $group_tag")
                        println(STDERR, "TRYING UMI FROM READ NAME ....")
                        UMI = split(tempname, umi_s)[umi_p]
                    end
                elseif algo == "umi"
                    UMI = split(tempname, umi_s)[umi_p]
                end

                spos = minimum([lpos, lpos2, rpos, rpos2])
                epos = maximum([lpos, lpos2, rpos, rpos2])

                if algo == "umi"
                    key = rfname*":"*string(spos)*":"*string(epos)*":"*UMI
                elseif algo == "position"
                    key = rfname*":"*string(spos)*":"*string(epos)
                else
                    key = rfname*":"*string(spos)*":"*string(epos)*":"*tempname
                end

                if algo == "paired"
                    # duplicated mutations found in R1 and R2 overlaping reads
                    variants = intersect(mutations, mutations2)
                else
                    # unique mutations from R1 and R2
                    variants = union(mutations, mutations2)
                end

                # if array is not empty after base quality filtering.

                if flag1&16 == 0
                    seq_len = rpos - lpos + 1
                    seq_len2 = rpos2 - lpos2 + 1
                else
                    seq_len2 = rpos - lpos + 1
                    seq_len = rpos2 - lpos2 + 1
                end

                total_bases += seq_len
                total_bases += seq_len2

                if haskey(vars, key)
                    populate_variants!(vars, variants, key, seq_len,
                                       seq_len2, spos, epos;
                                       verbose=verbose)
                else
                    Coverage.coverage[rfname][lpos:rpos] += 1
                    Coverage.coverage[rfname][lpos2:rpos2] += 1
                    # length = epos - spos + 1
                    populate_variants!(vars, variants, key, seq_len,
                                       seq_len2, spos, epos;
                                       verbose=verbose)
                end
            else
                if tempname != tempname2
                    println(STDOUT, "WARNING: READ NAMES DOES NOT MATCH. Is bam sorted by name?")
                    filtered_reads += 2
                end

                if rfname2 != rfname
                    println(STDOUT, "WARNING: READS ARE ON DIFF CHROMOSOMES [$tempname]")
                    filtered_reads += 2
                end

                if !(passed1)
                    filtered_R1 += 1
                end

                if !(passed2)
                    filtered_R2 += 1
                end
            end
        end
    end

    if verbose
        println(STDOUT, "INFO: COVERAGES OF MUTATIONS: ")
        for i in mut_coverages
            println(STDOUT, i)
        end
        println(STDOUT, "INFO: END OF COVERAGES OF MUTATIONS: ")
    end

    error_distr = error_distr./read_distr

    return vars, mut_coverages, total_bases, total_reads, filtered_reads, filtered_R1, filtered_R2, error_distr
end

function count_mut_coverages!(mut_coverages::Dict{String,Dict{Mutation,Int64}},
                              mutations::Array{Mutation,1}, refname::String;
                              verbose::Bool=false
                              )

    if verbose
        println(STDOUT, "INFO: COUNT MUTATION COVERAGES, input: $mut_coverages")
    end

    for mutation in mutations
        if verbose
            println(STDOUT, "INFO: COUNT MUTATION COVERAGES, MUTATION: $mutation")
        end

        # collecting coverages of mutations per reference
        if haskey(mut_coverages, refname)
            if haskey(mut_coverages[refname], mutation)
                mut_coverages[refname][mutation] += 1
            else
                mut_coverages[refname][mutation] = 1
            end
        else
            mut_coverages[refname] = Dict(mutation => 1)
        end

        if verbose
            println(STDOUT, "INFO: COUNT MUTATION COVERAGES, MUTATION COVERAGE: ", mut_coverages[refname][mutation])
        end

    end

    return mut_coverages
end

function populate_variants!(vars::Dict{String,MutationsCluster},
                            mutations::Array{Mutation, 1},
                            id::String, seq_len::Int64,
                            seq_len2::Int64,
                            spos::Int64, epos::Int64;
                            verbose::Bool=false)

    if !(haskey(vars, id))
        if verbose
            println(STDOUT, "INFO: POPULATE VARIANTS, init new cluster w/o mutations: R1_len:$seq_len R2_len:$seq_len2 spos:$spos epos:$epos, id:$id")
        end
        vars[id] = MutationsCluster(
            Dict(),
            0,
            seq_len,
            seq_len2,
            0,
            zeros(UInt16, seq_len),
            zeros(UInt16, seq_len2),
            spos,
            epos
        )
    end

    for mutation in mutations
        if verbose
            println(STDOUT, "INFO: POPULATE VARIANTS, adding $mutation $seq_len $seq_len2 $spos $epos to cluster $id", vars[id] )
        end
        if haskey(vars[id].mutation_count, mutation)
            vars[id].mutation_count[mutation] += 1
        else
            vars[id].mutation_count[mutation] = 1
            vars[id].total_mutations += 1
        end
    end

    if vars[id].R1_length < seq_len
        if verbose
            println(STDOUT, "INFO: POPULATE VARIANTS, R1_len increased: ", vars[id].R1_length, " $seq_len")
            println(STDOUT, "INFO: POPULATE VARIANTS, extending R1_coverage: ", length(vars[id].R1_coverage))
        end
        append!(vars[id].R1_coverage, zeros(seq_len-vars[id].R1_length))
        if verbose
            println(STDOUT, "INFO: POPULATE VARIANTS, extended R1_coverage: ", length(vars[id].R1_coverage))
        end
        vars[id].R1_length = seq_len
    end

    if vars[id].R2_length < seq_len2
        if verbose
            println(STDOUT, "INFO: POPULATE VARIANTS, R2_len increased: ", vars[id].R1_length, " $seq_len")
            println(STDOUT, "INFO: POPULATE VARIANTS, extending R1_coverage: ", length(vars[id].R1_coverage))
        end
        append!(vars[id].R2_coverage, zeros(seq_len2-vars[id].R2_length))
        if verbose
            println(STDOUT, "INFO: POPULATE VARIANTS, extended R1_coverage: ", length(vars[id].R1_coverage))
        end
        vars[id].R2_length = seq_len2
    end

    if verbose
        println(STDOUT, "INFO: POPULATE VARIANTS, incrementing R1_coverage: ", vars[id].R1_coverage)
    end

    vars[id].R1_coverage[1:seq_len] += 1

    if verbose
        println(STDOUT, "INFO: POPULATE VARIANTS, incremented R1_coverage: ", vars[id].R1_coverage)
        println(STDOUT, "INFO: POPULATE VARIANTS, incrementing R2_coverage: ", vars[id].R2_coverage)
    end

    vars[id].R2_coverage[1:seq_len2] += 1

    if verbose
        println(STDOUT, "INFO: POPULATE VARIANTS, incremented R2_coverage: ", vars[id].R2_coverage)
        println(STDOUT, "INFO: POPULATE VARIANTS, incrementing cluster size: ", vars[id].cluster_size)
    end

    vars[id].cluster_size += 1

    if verbose
        println(STDOUT, "INFO: POPULATE VARIANTS, incremented cluster size: ", vars[id].cluster_size)
    end

    return vars
end

function filter_variants(vars::Dict{String,MutationsCluster},
                     mut_coverages::Dict{String,Dict{Mutation,Int64}},
                     cvCoverage::Coverage, referenceSruct::ReferenceSequences,
                     cluster_size::Int64, ref_error::Float64; algo::String="umi",
                     freq_in_clust::Float64=1.0, verbose::Bool=false)

    if !(algo in ["umi", "position", "paired", "naive"])
        println(STDERR, "ERROR: Selected --algorithm does not exist: $algo")
        println(STDERR, "ERROR: Options are: umi position paired naive")
        exit(1)
    end

    true_var_count = 0
    err_distr = Float32[]
    clust_size_distr = Int64[]
    clusts_with_vars = 0
    total_clusters = 0
    found_insertions = 0
    found_deletions = 0
    found_mutations = 0
    found_ns = 0
    passing_vars = Dict{String,MutationsCluster}()
    confident_bp = 0
    total_bp = 0
    conf_coverage = Coverage(referenceSruct)

    for id in keys(vars)
        # parse key and get position of mutation
        chromosome = parse_chr(id)
        read_position = parse(Int64, split(id, ":")[2])
        found_var = false
        total_clusters += 1

        if verbose
            println("INFO: CALCULATE VARIANTS, CLUSTER ID: ", id)
            println("INFO: CALCULATE VARIANTS, CLUSTER: ", vars[id])
        end

        push!(clust_size_distr, vars[id].cluster_size)
        if vars[id].cluster_size >= cluster_size

            # ad cluster to confident coverage
            try
                start = vars[id].c_start
                end_ = vars[id].c_start + vars[id].R1_length - 1
                conf_coverage.coverage[chromosome][start:end_] += 1
                start2 = vars[id].c_end - vars[id].R2_length + 1
                end_2 = vars[id].c_end

                #    -------
                #               --------
                if start2 > end_
                    conf_coverage.coverage[chromosome][start2:end_2] += 1
                #    ---------
                # ---------
                elseif end_ > end_2 && start2 < start
                    conf_coverage.coverage[chromosome][start2:start-1] += 1
                #    -------
                #        --------
                elseif end_ > start2
                    conf_coverage.coverage[chromosome][end_+1:end_2] += 1
                end
            catch error
                println(STDERR, "ERROR: $error")
                println(STDERR, "ERROR: ", vars[id])
            end

            for variant in keys(vars[id].mutation_count)
                var_count = vars[id].mutation_count[variant]
                spos = variant.mStart
                epos = variant.mEnd

                # Filtering non confident variants in the cluster
                if var_count/vars[id].cluster_size >= freq_in_clust
                    mut_cov = mut_coverages[chromosome][variant]
                    coverage = mean(cvCoverage.coverage[chromosome][spos:epos])

                    # checking if mutation coverage is not equal to coverage
                    if mut_cov/coverage < ref_error

                        if verbose
                            println("INFO: CALCULATE VARIANTS, CLUSTER ID: ", id)
                            println("INFO: CALCULATE VARIANTS, READ_POSITION: ", read_position)
                            println("INFO: CALCULATE VARIANTS, COVERAGE: ", coverage)
                            println("INFO: CALCULATE VARIANTS, MUT_COVERAGE: ", mut_cov)
                            println("INFO: CALCULATE VARIANTS, ERROR_COUNT: ", var_count)
                        end

                        found_var = true
                        true_var_count += 1

                        if variant.class == 'M'
                            found_mutations += 1
                        elseif variant.class == 'I'
                            found_insertions += 1
                        elseif variant.class == 'D'
                            found_deletions += 1
                        elseif variant.class == 'N'
                            found_ns += 1
                        end

                        if haskey(passing_vars, id)
                            passing_vars[id].mutation_count[variant] = 1
                            passing_vars[id].total_mutations += 1
                        else
                            passing_vars[id] = MutationsCluster(
                            Dict(variant => 1),
                            vars[id].cluster_size,
                            vars[id].R1_length,
                            vars[id].R2_length,
                            1,
                            vars[id].R1_coverage,
                            vars[id].R2_coverage,
                            vars[id].c_start,
                            vars[id].c_end
                            )
                        end

                        if verbose
                            println(STDOUT, "INFO: CALCULATE VARIANTS, PASSING MUTATION: $variant")
                        end
                    end

                else
                    if verbose
                        println("INFO: CALCULATE VARIANTS, CLUSTER ID: ", id)
                        println("INFO: CALCULATE VARIANTS, SKIPPING VAR: ", variant)
                    end

                end

                # collect error frequencies for distribution ploting.
                if var_count != 0
                    push!(err_distr, round(var_count/vars[id].cluster_size, 2))
                end
            end

            if algo == "umi" || algo == "position" || algo == "naive"
                # adding only confident bases, equal to cluster coverage.
                conf_bp = count(x -> x >= vars[id].cluster_size, vars[id].R1_coverage)
                t_bp = conf_bp * vars[id].cluster_size
            elseif algo == "paired"
                # adding only confident bases from overlaping pairs
                if vars[id].c_start + vars[id].R1_length > vars[id].c_end - vars[id].R2_length && vars[id].total_mutations > 0
                    conf_bp = vars[id].c_start + vars[id].R1_length - (vars[id].c_end - vars[id].R2_length) + 1
                    t_bp = vars[id].R1_length + vars[id].R2_length
                else
                    conf_bp = 0
                    t_bp = 0
                end
            end

            if verbose
                println(STDOUT, "INFO: CALCULATE VARIANTS, counted confident bp for R1 $id: $conf_bp")
                println(STDOUT, "INFO: CALCULATE VARIANTS, counted total cluster bp for R1 $id: $t_bp")
            end

            confident_bp += conf_bp
            total_bp += t_bp

            if algo == "umi" || algo == "position" || algo == "naive"
                conf_bp = count(x -> x >= vars[id].cluster_size, vars[id].R2_coverage)
                t_bp = conf_bp * vars[id].cluster_size

                if verbose
                    println(STDOUT, "INFO: CALCULATE VARIANTS, counted confident bp for R2 $id: $conf_bp")
                    println(STDOUT, "INFO: CALCULATE VARIANTS, counted total cluster bp for R2 $id: $t_bp")
                end

                confident_bp += conf_bp
                total_bp += t_bp
            end

        else
            if verbose
                println("INFO: CALCULATE VARIANTS, CLUSTER ID: ", id)
                println("INFO: CALCULATE VARIANTS, SKIPPING CLUSTER: ", vars[id])
            end
        end

        if found_var
            clusts_with_vars += 1
        end
    end

    return (
        true_var_count,
        err_distr,
        clust_size_distr,
        found_mutations,
        found_insertions,
        found_deletions,
        found_ns,
        total_clusters,
        confident_bp,
        clusts_with_vars,
        passing_vars,
        total_bp,
        conf_coverage
    )
end

function error_dist_over_length!(error_distr::Array{Int64, 1}, mutations::Array{Mutation,1}, lpos::Int64, read_length::Int64)

    if read_length > length(error_distr)
        for i in 1:(read_length-length(error_distr))
            push!(error_distr, 0)
        end
    end

    # find if any deletions

    l = 0
    # TODO: mutations should be sorted.
    for mutation in mutations
        if mutation.class == 'D'
            l += length(mutation.reference)
        end
        position_in_read = mutation.mStart - lpos - l + 1
        error_distr[position_in_read] += 1
    end

    return error_distr
end

function read_dist!(read_distr::Array{Int64, 1}, read_length::Int64)

    if read_length > length(read_distr)
        for i in 1:(read_length-length(read_distr))
            push!(read_distr, 0)
        end
    end

    read_distr .+= 1

    return read_distr
end

"""
    Translates collection of mutations to collection of aminoacid subtitutions.
    For genetic codes check BioSequences.ncbi_trans_table.
"""
function translate_mutations(referenceSruct::ReferenceSequences, genesCSV::String,
                             vars::Dict{String,MutationsCluster};
                             genetic_code::Int64=11, verbose::Bool=false)

    # get AA sequence from reference
    # loop mutations, get codon, translate.

    if verbose
        println(STDOUT, "INFO: TRANSLATE MUTATIONS, Loading BED file. $genesCSV")
    end
    df = CSV.read(genesCSV, delim=',')
    proteins = BioSequences.AminoAcidSequence[]
    genes_seqs = BioSequences.BioSequence{BioSequences.DNAAlphabet{4}}[]
    chromosomes = String[]
    gene_sidx = Int64[]
    gene_eidx = Int64[]
    strands = String[]
    names = String[]
    frame_shifts = 0
    # chr => variant_clusters[variant_cluster[variants]]
    aa_substitutions = Dict{String,Array{Array{String,1},1}}()

    if verbose
        println(STDOUT, "INFO: TRANSLATE MUTATIONS, Parsing records.")
    end

    for row in eachrow(df)
        if verbose
            println(STDOUT, "INFO: TRANSLATE MUTATIONS, Parsing recod: $row")
        end

        chr = row[1]
        strand = row[5]
        gene_start = row[2]
        gene_end = row[3]
        gene_name = row[4]

        if verbose
            println(STDOUT, "INFO: TRANSLATE MUTATIONS, strand: $strand")
            println(STDOUT, "INFO: TRANSLATE MUTATIONS, chromosome: $chr")
            println(STDOUT, "INFO: TRANSLATE MUTATIONS, gene start: $gene_start")
            println(STDOUT, "INFO: TRANSLATE MUTATIONS, gene end: $gene_end")
            println(STDOUT, "INFO: TRANSLATE MUTATIONS, gene name: $gene_name")
            println(STDOUT, "INFO: TRANSLATE MUTATIONS, loading gene sequence from reference.")
        end

        dna_sequence = DNASequence(referenceSruct.sequences[chr][gene_start:gene_end])
        if verbose
            println(STDOUT, dna_sequence)
            println(STDOUT, "INFO: TRANSLATE MUTATIONS, loading gene sequence from reference. [DONE]")
            println(STDOUT, "INFO: TRANSLATE MUTATIONS, convert DNA sequence to RNA.")
        end

        push!(genes_seqs, dna_sequence)

        if strand == "-"
            if verbose
                println(STDOUT, "INFO: TRANSLATE MUTATIONS, geting reverse complement.")
            end
            BioSequences.reverse_complement!(dna_sequence)
            if verbose
                println(STDOUT, "INFO: TRANSLATE MUTATIONS, geting reverse complement. [DONE]")
            end
        end

        rna_sequence = convert(RNASequence, dna_sequence)
        if verbose
            println(STDOUT, rna_sequence)
            println(STDOUT, "INFO: TRANSLATE MUTATIONS, convert DNA sequence to RNA. [DONE]")
        end

        push!(gene_sidx, gene_start)
        push!(gene_eidx, gene_end)
        push!(chromosomes, chr)
        push!(strands, strand)
        push!(names, gene_name)

        if verbose
            println(STDOUT, "INFO: TRANSLATE MUTATIONS, translating to protein sequence.")
        end
        protein = translate(rna_sequence, code=BioSequences.ncbi_trans_table[genetic_code], allow_ambiguous_codons=true)
        if verbose
            println(STDOUT, "INFO: TRANSLATE MUTATIONS, CHR: $chr; STRAND: $strand; START: $gene_start; END: $gene_end; PROTEIN: $protein")
        end

        push!(proteins, protein)
    end

    if verbose
        println(STDOUT, "INFO: TRANSLATE MUTATIONS, slelected genes chromosomes: $chromosomes")
        println(STDOUT, "INFO: TRANSLATE MUTATIONS, slelected genes starts: $gene_sidx")
        println(STDOUT, "INFO: TRANSLATE MUTATIONS, slelected genes ends: $gene_eidx")
        println(STDOUT, "INFO: TRANSLATE MUTATIONS, slelected genes strands: $strands")
    end

    # mutaion is Chromosome:ReferencePositionMutation
    for id in keys(vars)
        chr = parse_chr(id)

        if verbose
            println(STDOUT, "INFO: TRANSLATE MUTATIONS, cluster of mutations at $chr: ", vars[id])
        end
        perm_seq = DNASequence()
        gene_idx = 0
        found_first = false
        found_gene = false
        found_indel = false
        indel_ct = 0

        for var in keys(vars[id].mutation_count)
            mut_pos = var.mStart
            # search for gene matching cluster range
            if verbose
                println(STDOUT, "INFO: TRANSLATE MUTATIONS, Searching if var $var is matching any gene from BED.")
            end

            chr_idxs = find(x -> x == chr, chromosomes)
            gene_starts = find(x -> x <= mut_pos, gene_sidx)
            gene_ends = find(x -> x >= mut_pos, gene_eidx)
            common1 = intersect(chr_idxs, gene_starts)
            gene_idx2 = intersect(common1, gene_ends)

            if verbose
                println(STDOUT, "INFO: TRANSLATE MUTATIONS, Matching chromosomes indexes: $chr_idxs")
                println(STDOUT, "INFO: TRANSLATE MUTATIONS, Matching gene starts: $gene_starts")
                println(STDOUT, "INFO: TRANSLATE MUTATIONS, Matching gene ends indexes: $gene_ends")
                println(STDOUT, "INFO: TRANSLATE MUTATIONS, Matching both chromosomes and gene starts: $common1")
                println(STDOUT, "INFO: TRANSLATE MUTATIONS, Matching all: $gene_idx2")
            end

            if length(gene_idx2) < 1
                if verbose
                    println(STDOUT, "INFO: TRANSLATE MUTATIONS, $var is not in your gene list therefore will not be translated.")
                end
                continue
            elseif length(gene_idx2) > 1
                if verbose
                    println(STDOUT, "WARNING: TRANSLATE MUTATIONS, $var is matching more than one gene. Does ranges overlap? This will not be handled!")
                end
                continue
            else
                found_gene = true
                gene_idx = gene_idx2
                if verbose
                    println(STDOUT, "INFO: TRANSLATE MUTATIONS, $var found in ", names[gene_idx][1])
                end
            end

            # get gene sequence
            if !(found_first)
                perm_seq = copy(genes_seqs[gene_idx[1]])
                found_first = true
            end

            # do permutations
            if found_gene
                if var.class != 'D' && var.class != 'I'
                    substitution = DNASequence(DNASequence(string(var.variant)))
                    mut_gene_pos = mut_pos - gene_sidx[gene_idx] + 1
                    perm_seq[mut_gene_pos] = substitution
                else
                    if verbose
                        println(STDOUT, "INFO: TRANSLATE MUTATIONS, indel found: ", var)
                        println(STDOUT, "INFO: TRANSLATE MUTATIONS, skipping this cluster: $id")
                    end
                    frame_shifts += 1
                    found_indel = true
                    break
                end
            end
        end

        if found_gene && !(found_indel)
            if strands[gene_idx] == "-"
                BioSequences.reverse_complement!(perm_seq)
            end

            mut_rna = convert(RNASequence, perm_seq)

            mut_protein = translate(mut_rna, code=BioSequences.ncbi_trans_table[genetic_code], allow_ambiguous_codons=true)

            # Comparing WT vs mutated protein and extracting
            protein = proteins[gene_idx][1]
            if verbose
                println(STDOUT, "INFO: TRANSLATE MUTATIONS, mutated sequence: $perm_seq")
                println(STDOUT, "INFO: TRANSLATE MUTATIONS, wt sequence:      ", convert(String, genes_seqs[gene_idx[1]]))
                println(STDOUT, "INFO: TRANSLATE MUTATIONS, mutated protein sequence: $mut_protein")
                println(STDOUT, "INFO: TRANSLATE MUTATIONS, wt protein sequence:      $protein")
            end
            aa_subs = compare_seq(protein, mut_protein; verbose=verbose)

            if haskey(aa_substitutions, chr)
                push!(aa_substitutions[chr], aa_subs)
            else
                aa_substitutions[chr] = [aa_subs]
            end
        end
    end

    return aa_substitutions, frame_shifts
end

function compare_seq(seqa::Any, seqb::Any; verbose::Bool=false)

    muts = String[]
    ix = 0
    for i in seqa
        ix += 1
        if i == seqb[ix]
            continue
        else
            mut = string(i)*string(ix)*string(seqb[ix])
            push!(muts, mut)
        end
    end

    if verbose
        println(STDOUT, "INFO: COMPARE SEQ, found translated mutations: $muts")
    end

    return muts
end

function mapp_mutations(vars::Dict{String, MutationsCluster}; cluster::Bool=true)

    if cluster
        mapp_count = Dict{Array{Mutation,1},Int64}()
    else
        mapp_count = Dict{Mutation,Int64}()
    end

    for id in keys(vars)
        mutations = collect(keys(vars[id].mutation_count))
        if cluster
            if mutations in keys(mapp_count)
                mapp_count[mutations] += 1
            else
                mapp_count[mutations] = 1
            end
        else
            for mutation in mutations
                if mutation in keys(mapp_count)
                    mapp_count[mutation] += 1
                else
                    mapp_count[mutation] = 1
                end
            end
        end
    end

    return mapp_count
end

function parse_chr(id::String)

    return String(split(id,":")[1])
end

function stringify(mutation::Mutation)

    mutation_str = string(mutation.reference)
    mutation_str = mutation_str * string(mutation.mStart)
    mutation_str = mutation_str * string(mutation.variant)

    return mutation_str
end

function conf_coverage!(Coverage::Coverage, vars::Dict{String,MutationsCluster})

    for id in keys(vars)
        chr = parse_chr(id)
        try
            start = vars[id].c_start
            end_ = vars[id].c_start + vars[id].R1_length - 1
            Coverage.coverage[chr][start:end_] += 1
            start2 = vars[id].c_end - vars[id].R2_length + 1
            end_2 = vars[id].c_end

            #    -------
            #               --------
            if start2 > end_
                Coverage.coverage[chr][start2:end_2] += 1
            #    ---------
            # ---------
            elseif end_ > end_2 && start2 < start
                Coverage.coverage[chr][start2:start-1] += 1
            #    -------
            #        --------
            elseif end_ > start2
                Coverage.coverage[chr][end_+1:end_2] += 1
            end
        catch error
            println(STDERR, "ERROR: $error")
            println(STDERR, "ERROR: ", vars[id])
        end
    end

    return Coverage
end

function transvers(vars::Dict{String,MutationsCluster})

    types = Dict(
        "AG"=>0,
        "GA"=>0,
        "TC"=>0,
        "CT"=>0,
        "AC"=>0,
        "CA"=>0,
        "GT"=>0,
        "TG"=>0,
        "AT"=>0,
        "TA"=>0,
        "CG"=>0,
        "GC"=>0
        )
    trans = 0
    vers = 0
    for id in keys(vars)
        for var in keys(vars[id].mutation_count)
            if var.class == 'M'
                if var.reference == "A" && var.variant == "G"
                    trans += 1
                    types["AG"] += 1
                elseif var.reference == "G" && var.variant == "A"
                    trans += 1
                    types["GA"] += 1
                elseif var.reference == "T" && var.variant == "C"
                    trans += 1
                    types["TC"] += 1
                elseif var.reference == "C" && var.variant == "T"
                    trans += 1
                    types["CT"] += 1
                elseif var.reference == "A" && var.variant == "C"
                    vers += 1
                    types["AC"] += 1
                elseif var.reference == "C" && var.variant == "A"
                    vers += 1
                    types["CA"] += 1
                elseif var.reference == "G" && var.variant == "T"
                    vers += 1
                    types["GT"] += 1
                elseif var.reference == "T" && var.variant == "G"
                    vers += 1
                    types["TG"] += 1
                elseif var.reference == "A" && var.variant == "T"
                    vers += 1
                    types["AT"] += 1
                elseif var.reference == "T" && var.variant == "A"
                    vers += 1
                    types["TA"] += 1
                elseif var.reference == "C" && var.variant == "G"
                    vers += 1
                    types["CG"] += 1
                elseif var.reference == "G" && var.variant == "C"
                    vers += 1
                    types["GC"] += 1
                end
            end
        end
    end

    return trans, vers, types
end

function mapp_single_aa_subs(arrays::Array{Array{String,1}})

    aa_subs = Dict{String, UInt64}()
    for array in arrays
        for aa_sub in array
            if haskey(aa_subs, aa_sub)
                aa_subs[aa_sub] += 1
            else
                aa_subs[aa_sub] = 1
            end
        end
    end

    return aa_subs
end

end  # module

max_stars_repo_name: KarolisM/rnd-it-molbio-fidelity-umi

max_issues_repo_name: KarolisM/rnd-it-molbio-fidelity-umi

max_forks_repo_name: KarolisM/rnd-it-molbio-fidelity-umi


text:
using BenchmarkTools

function explode(str::String)
    newstr = str
    while true
        a = [str[i:i] for i in eachindex(str)]
        isdig(s) = isdigit(s[begin])
        function int(a, i)
            return (j = findfirst(!isdig, a[(i + 1):end]);
            parse(Int, join(a[i:(j == nothing ? length(a) : j)])))
        end
        nesting = accumulate(+, [
            if c == "["
                1
            elseif c == "]"
                -1
            else
                0
            end for c in a
        ])
        n = findfirst(>(4), nesting)
        n == nothing && return str
        isdigits = [isdig(s) for s in a]
        nums = [
            if !isdigits[i]
                -1
            elseif isdigits[i - 1]
                -2
            else
                int(join(a[i:end]), 1)
            end for i in eachindex(a)
        ]
        n = findfirst(==(first(findmax(nesting))), nesting)
        x, y = nums[n + 1], nums[n + findfirst(==(-1), nums[(n + 1):end]) + 1]
        if (prev = findlast(>(-1), nums[begin:(n - 1)])) != nothing
            nums[prev] += x
        end
        rbrack = findfirst(==("]"), a[(n + 1):end]) + n
        if (succ = findfirst(>(-1), nums[(rbrack + 1):end])) != nothing
            nums[rbrack + succ] += y
        end
        rbpos = findfirst(==(-1), nums[(n + 4):end]) + n + 3
        newstr =
            join([
                if nums[i] == -2
                    ""
                elseif nums[i] == -1
                    a[i]
                else
                    string(nums[i])
                end for i in eachindex(a[begin:(n - 1)])
            ]) *
            "0" *
            join([
                if nums[rbpos + i] == -2
                    ""
                elseif nums[rbpos + i] == -1
                    a[rbpos + i]
                else
                    string(nums[rbpos + i])
                end for i in eachindex(a[(rbpos + 1):end])
            ])
        newstr == str && return newstr
        str = newstr
    end
end

function snailsplit(s)
    newstr = s
    while true
        if occursin(r"\d\d+", s)
            newstr = replace(
                s,
                r"\d\d+" =>
                    (x) ->
                        "[" *
                        string(parse(Int, x) ÷ 2) *
                        ',' *
                        string((parse(Int, x) + 1) ÷ 2) *
                        "]";
                count=1,
            )
            newstr = explode(newstr)
        end
        newstr == s && return newstr
        s = newstr
    end
end

function snailadd(s1, s2)
    s3 = "[" * s1 * "," * s2 * "]"
    while true
        newstring = snailsplit(explode(s3))
        newstring == s3 && break
        s3 = newstring
    end
    return s3
end

listadd(a) = reduce(snailadd, a)

function magnitude(arr::Vector)
    fir, las = first(arr), last(arr)
    return 3 * (fir isa Integer ? fir : magnitude(fir)) +
           2 * (las isa Integer ? las : magnitude(las))
end

magnitude(s::String) = magnitude(eval(Meta.parse(s)))
snailcheck(list) = magnitude(listadd(list))

function day18()
    lines = strip.(readlines("AoCdata/AoC_2021_day18.txt"))
    part = [0, 0]
    part[1] = snailcheck(lines)
    checksums = Int[]
    for x in lines, y in lines
        x == y && continue
        push!(checksums, snailcheck([x, y]))
    end
    part[2] = maximum(checksums)
    return part
end

part = day18()
println("Part 1: ", part[1])
println("Part 2: ", part[2])

@btime day18()


text:
#= ForwardDiff is removed from dependencies in this version.
import MechanicalSketch
import MechanicalSketch: empty_figure, background, sethue, O, WI, HE, EM, FS, finish
import MechanicalSketch: PALETTE, settext, setline
import MechanicalSketch: dimension_aligned, move, do_action, noise, rk4_steps!
import MechanicalSketch: generate_complex_potential_source, generate_complex_potential_vortex
import MechanicalSketch: @import_expand, Quantity, @layer, x_y_iterators_at_pixels
import MechanicalSketch: place_image, draw_legend, set_scale_sketch
import MechanicalSketch: ∙, ∇_rectangle, SVector, normalize_datarange, draw_streamlines
import Interpolations: interpolate, Gridded, Linear, Flat, extrapolate
import MechanicalSketch: Greys_9, BinLegend, circle

let
empty_figure(filename = joinpath(@__DIR__, "test_27.png"));

if !@isdefined m²
    @import_expand(m, s)
end

include("test_functions_24.jl")
# Reuse the flow field from test_23.jl, a matrix of complex velocities with an element per pixel.
A = flowfield_23();

# Duration to trace the streamlines forward and back
global const DU_2 = 2.0s
include("test_functions_27.jl")
# Make a matrix containing one complex quantity per pixel.
# In test_23 we displayed this using colour for direction
# and luminosity for magnitude.
physwidth = 10.0m
height_relative_width = 0.4
physheight = physwidth * height_relative_width
screen_width_frac = 2 / 3
set_scale_sketch(physwidth, round(Int, screen_width_frac * WI))

str = "Noise has maximum wavelength $CUTOFF_23 · $DU_2 = $(CUTOFF_23 * DU_2). Experiments with noise spectrum."
settext(str, O + (-WI/2 + EM, -0.5HE + 2EM), markup = true)

# We are going to use the velocity vector field in a lot of calculations,
# and interpolate between the calculated pixel values
xs, ys = x_y_iterators_at_pixels(;physwidth, physheight)
fxy_inter = interpolate((xs, reverse(ys)), map( cmplx -> (real(cmplx), imag(cmplx)), transpose(A)[ : , end:-1:1]), Gridded(Linear()));
fxy = extrapolate(fxy_inter, Flat());

global const OU_27 = O + (0.0, -0.25HE + EM)
global const NO_3 = noisepic_3(A, physheight, physwidth)
legend = BinLegend(;maxlegend = 1.0, noofbins = 128,
                       colorscheme = reverse(Greys_9), name = Symbol("Value{Float64}"))
upleftpoint, lowrightpoint = place_image(OU_27, legend.(NO_3))
draw_legend(lowrightpoint + (EM, 0) + (0.0m, PHYSHEIGHT_23), legend)

# And make a function that interpolates between pixels:
nxy_inter = interpolate((xs, reverse(ys)), transpose(NO_3)[ : , end:-1:1], Gridded(Linear()));
nxy = extrapolate(nxy_inter, Flat());
nxy(0.1m, 0.02m)

# Centre of the bottom figure
OB = O + (0.0, + 0.25HE + 0.5EM )

# Number of sampling points along our streamline in each direction (effectively one less than this)
global const NS_3 = 20

# Step length (duration)
h = DU_2 / (NS_3 -1 )

@time M = convolute_image_3(xs, reverse(ys), fxy, nxy, h, CUTOFF_23) # 36.019 s (10 allocations: 4.85 MiB)
                                                  # 58.803552 seconds (62.34 M allocations: 2.247 GiB, 0.56% gc time)
                                                  # 57.021068 seconds (62.33 M allocations: 2.247 GiB, 0.59% gc time)
                                                  # 59.596110 seconds (59.04 M allocations: 2.002 GiB, 0.58% gc time)

botlegend = BinLegend(;maxlegend = maximum(M), noofbins = 128,
                       colorscheme = reverse(Greys_9), name = Symbol("Value{Float64}"))
upleftpoint, lowrightpoint = place_image(OB, botlegend.(M) )
draw_legend(lowrightpoint + (EM, 0) + (0.0m, PHYSHEIGHT_23), botlegend)

@layer begin
    sethue(PALETTE[1])
    setline(8)
    draw_streamlines(OB, fxy)
    setline(1)
end

dimension_aligned(OB + (-physwidth / 2, physheight / 2), OB + (physwidth / 2, physheight / 2))
dimension_aligned(OB + (-physwidth / 2, - physheight / 2 ),  OB +  (-physwidth / 2, physheight / 2 ))

str = " "
settext(str, O + (-WI/2 + EM, EM) , markup = true)
str = "First implementation Finite Impulse Response filter, Blackman window."
settext(str, O + (-WI/2 + EM, 2 * EM) , markup = true)

finish()
set_scale_sketch()
nothing
end
=#


text:
"""
    QuasiKron(A, B, C...)

is a lazy representation of a Kronecker product of the quasi-arrays A, B, C.... Unlike
standard kronecker product, we also kronecker the axes, using tuples.
"""

struct QuasiKron{T, N, Args<:Tuple} <: AbstractQuasiArray{T, N}
    args::Args
end

QuasiKron{T,N}(a...) where {T,N} = QuasiKron{T,N,typeof(a)}(a)
QuasiKron{T}(a::AbstractQuasiOrVector...) where T = QuasiKron{mapreduce(eltype,promote_type,a),1}(a...)
QuasiKron{T}(a::AbstractQuasiOrMatrix...) where T = QuasiKron{mapreduce(eltype,promote_type,a),2}(a...)
QuasiKron(a...) = QuasiKron{mapreduce(eltype,promote_type,a)}(a...)

quasikron(a::Inclusion...) = Inclusion(ProductDomain(map(domain,a)...))
quasikron(a...) = QuasiKron(a...)

_kronaxes(::Tuple{}...) = ()
_kronaxes(a::Tuple...) = (quasikron(map(first,a)...), _kronaxes(map(tail,a)...)...)
axes(K::QuasiKron) = _kronaxes(map(axes,K.args)...)

function getindex(K::QuasiKron{<:Any,1}, i::Tuple)
    @boundscheck checkbounds(K, i)
    prod(getindex.(K.args,i))
end


text:
using Dates

c=Cube()

@testset "Access single variable" begin
d = subsetcube(c,variable="air_temperature_2m",lon=(10,11),lat=(51,50),
                time=(Date("2002-01-01"),Date("2008-12-31")))

d.data.v.indices==(18:21, 17:20, 93:414)
@test d.axes[1].values==10.125:0.25:10.875
@test d.axes[2].values==50.875:-0.25:50.125
@test d.axes[1] == RangeAxis("lon", 10.125:0.25:10.875)
@test d.axes[2] == RangeAxis("lat", 50.875:-0.25:50.125)
end

@testset "Access multiple variables" begin
d2 = subsetcube(c,variable=["air_temperature_2m","gross_primary_productivity"],lon=(10,11),lat=(50,51),
                time=(Date("2002-01-01"),Date("2008-12-31")))

  @test d2.axes[4].values == ["air_temperature_2m", "gross_primary_productivity"]
  @test d2.data.arrays[1].v.indices ==(18:21, 17:20, 93:414)
  @test d2.axes[1].values==10.125:0.25:10.875
  @test d2.axes[2].values==50.875:-0.25:50.125
  @test first(d2.axes[3].values) == Dates.Date(2002,1,5)
  @test last(d2.axes[3].values) == Dates.Date(2008, 12, 30)
end

@testset "Test values in MemCube" begin
d = subsetcube(c,variable="air_temperature_2m",lon=(10,11),lat=(51,50),
                time=(Date("2002-01-01"),Date("2008-12-31")))
d2 = subsetcube(c,variable=["air_temperature_2m","gross_primary_productivity"],lon=(10,11),lat=(50,51),
                time=(Date("2002-01-01"),Date("2008-12-31")))
data1=readcubedata(d)
data2=readcubedata(d2)

@test size(data1.data)==(4,4,322)
@test size(data2.data)==(4,4,322,2)

@test isapprox(data1.data[1,1,1:10],Float32[267.9917, 269.9631, 276.71036, 280.88998,
  280.90665, 277.02243, 274.5466, 276.919, 279.96243, 279.42276])

@test isapprox(data2.data[1,1,1:10,1],Float32[267.9917, 269.9631, 276.71036, 280.88998,
  280.90665, 277.02243, 274.5466, 276.919, 279.96243, 279.42276])

@test caxes(data1)[1:2]==CubeAxis[RangeAxis("lon",10.125:0.25:10.875),RangeAxis("lat",50.875:-0.25:50.125)]

tax = caxes(data1)[3]
@test YAXArrays.Cubes.Axes.axsym(tax)==:time
@test tax.values[1] == Date(2002,1,5)
@test tax.values[end] == Date(2008,12,30)
@test length(tax.values) == 7*46

@test caxes(data2)[[1,2,4]]==CubeAxis[RangeAxis("lon",10.125:0.25:10.875),RangeAxis("lat",50.875:-0.25:50.125),CategoricalAxis("Variable",["air_temperature_2m","gross_primary_productivity"])]

tax = caxes(data2)[3]
@test YAXArrays.Cubes.Axes.axsym(tax)==:time
@test tax.values[1] == Date(2002,1,5)
@test tax.values[end] == Date(2008,12,30)
@test length(tax.values) == 7*46

end

@testset "Coordinate extraction" begin
d = subsetcube(c,variable="air_temperature_2m",lon=(10,11),lat=(51,50),
                time=(Date("2002-01-01"),Date("2008-12-31")))
data1=readcubedata(d)
# Test reading of coordinate list
ll=[10.1 50.2;10.5 51.1;10.8 51.1]
llcube = readcubedata(extractLonLats(data1,ll))
@test llcube.data[1,:]==data1.data[1,4,:]
@test llcube.data[2,:]==data1.data[3,1,:]
@test llcube.data[3,:]==data1.data[4,1,:]
end

@testset "Accessing regions" begin
#Test access datacube by region
d3 = subsetcube(c,variable="gross_primary_productivity",region="Austria",time=Date("2005-01-01"))
@test d3.axes==[RangeAxis("lon",9.625:0.25:14.875),RangeAxis("lat",48.875:-0.25:47.375)]
end

using DiskArrayTools: DiskArrayStack

@testset "Saving and loading cubes" begin
  d = subsetcube(c,variable="air_temperature_2m",lon=(10,31),lat=(51,50),
                  time=(Date("2002-01-01"),Date("2008-12-31")))
  data1=readcubedata(d)
  #Test saving cubes
  dire=tempname()
  YAXdir(dire)
  savecube(data1,"mySavedCube")

  data3=readcubedata(loadcube("mySavedCube"))
  @test data1.axes==data3.axes
  @test data1.data==data3.data

  # Test loadOrGenerate macro
  d=subsetcube(c,time=Date(2001)..Date(2005),lon=(10,11),lat=(50,51),variable=["gross_primary_productivity","net_ecosystem_exchange"])

  rmcube("Anomalies")
  @loadOrGenerate danom=>"Anomalies" begin
      danom = removeMSC(d)
  end

  @test danom.data isa Array

  @loadOrGenerate danom=>"Anomalies" begin
      error("This should never execute")
  end;
  @test danom.data isa DiskArrayStack

  zp = tempname()
  savecube(danom, zp, overwrite=true)

  @test danom.data isa DiskArrayStack

  danom=readcubedata(danom)
  danom2=readcubedata(loadcube(zp))

  @test danom.axes==danom2.axes
  @test all(map(isequal,danom.data,danom2.data))

  ncf = string(tempname(),".nc")
  savecube(danom,ncf,backend=:netcdf)

  using NetCDF, Dates
  #Test exportcube
  @test ncread(ncf,"lon") == 10.125:0.25:10.875
  @test ncread(ncf,"lat") == 50.875:-0.25:50.125
  @test ncgetatt(ncf,"time","units") == "days since 1980-01-01"
  @test getAxis("time",danom).values .- DateTime(1980) == Millisecond.(Day.(ncread(ncf,"time")))

  anc = replace(ncread(ncf,"gross_primary_productivity")[:,:,:],-9999.0=>missing)
  @test all(isequal.(anc, danom.data[:,:,:,1]))
  neear = replace(danom.data[:,:,:,2],missing=>-9999.0)
  @test all(isequal.(ncread(ncf,"net_ecosystem_exchange")[:,:,:],neear))
  cnc = Cube(ncf)
  @test all(cnc.data[:,:,:,:] .== danom.data[:,:,:,:])
end