text:
import Base: getindex, setindex!

abstract type  Absolute end
abstract type  Brake end
abstract type  ControlPosition end
abstract type  Info end
abstract type  Edges end
abstract type  Engine end
abstract type  Feedback end
abstract type  Hard end
abstract type  Home end
abstract type  Kind end
abstract type  Loft end
abstract type  Move end
abstract type  NameFRAM end
abstract type  NameStage end
abstract type  Position end
abstract type  Protection end
abstract type  Power end
abstract type  PID end
abstract type  Relative end
abstract type  Settings end
abstract type  SetZero end
abstract type  Soft end
abstract type  State end
abstract type  Status end
abstract type  Stop end

export Device

"""
    Device
Abstract super-class for devices.
"""
abstract type  Device end

const PrimType = Union{Bool, Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Int128, UInt128, Float16, Float32, Float64}

function _value_to_tuple(x::T) where T
    types = Tuple(T.types)
    Tuple([getfield(x, i) for i in eachindex(types)])
end

_value_to_tuple(x::PrimType) = (x,)

function ximc_crc16(data::AbstractVector{UInt8})
    crc = UInt16(0xffff)
    for byte in data
        crc = xor(crc,byte);
        for bitno in 0:7
            a = crc;
            carry_flag = a & UInt16(0x0001);
            crc = crc >>> 1;
            if (carry_flag == 1)
                crc = UInt16(xor(crc, 0xa001))
            end
        end
    end
    crc::UInt16
end

@assert ximc_crc16(UInt8.([75, 29, 243])) == 49992

ximcio_to_jl(x) = ltoh(x)
jl_to_ximcio(x) = htol(x)

mutable struct XIMCString{N}
    bytes::NTuple{N, UInt8}
end

Base.read(io::IO, ::Type{XIMCString{N}}) where N = XIMCString((read(io, N)...))
Base.write(io::IO, ::Type{XIMCString{N}}) where N = write(io, x.bytes...)

ximcio_to_jl(x::XIMCString) = x
jl_to_ximcio(x::XIMCString) = x

Base.convert(::Type{String}, x::XIMCString) = String([x.bytes...])
Base.convert(::Type{XIMCString}, x::String) = XIMCString((String(x).data...))

XIMCString(x::String) = convert(XIMCString, x)
String(x::XIMCString) = convert(String, x)

Base.show(io::IO, x::XIMCString) = print(io,"XIMCString(\"$(String(x))\")")

mutable struct XIMCReserved{N}
    bytes::NTuple{N, UInt8}
end

Base.read(io::IO, ::Type{XIMCReserved{N}}) where N = XIMCReserved( Tuple(read(io, N)) )
Base.write(io::IO, x::XIMCReserved) = write(io, x.bytes...)

Base.show(io::IO, x::XIMCReserved) = print(io,"XIMCReserved(...)")

function _ximc_send_cmd(io::IO, cmd_code::String, data::Vector{UInt8} = UInt8[])
    @assert(length(cmd_code) == 4)
    if length(data) > 0
        crc = ximc_crc16(data)
        write(io, cmd_code, data, jl_to_ximcio(crc))
        # info("checksum = $crc")
    else
        write(io, cmd_code)
    end
end

function _ximc_send_cmd(io::IO, cmd_code::String, x)
    data_buf = IOBuffer()
    println(map(jl_to_ximcio, _value_to_tuple(x))...)
    write(data_buf, map(jl_to_ximcio, _value_to_tuple(x))...)
    data = take!(data_buf)
    _ximc_send_cmd(io, cmd_code, data)
end

function _ximc_read_resp(io::IO, cmd_code::String, n_data_bytes::Integer)
    @assert(length(cmd_code) == 4)
    @assert(n_data_bytes >= 0)

    with_data = n_data_bytes > 0

    nbytes_total = n_data_bytes + 4 + (with_data ? 2 : 0)
    # info("Reading $nbytes_total bytes")
    resp = read(io, nbytes_total)
    @assert(length(resp) == nbytes_total)
    buf = IOBuffer(resp)
    # info("_ximc_read_resp: resp = $resp")
    recv_cmd_bytes = read(buf, 4)
    recv_cmd = String(recv_cmd_bytes)
    # info("_ximc_read_resp: recv_cmd = $recv_cmd")
    (recv_cmd != cmd_code) && error("Received $recv_cmd, expected $cmd_code")
    if with_data
        data = read(buf, n_data_bytes)
        crc = ximcio_to_jl(read(buf, UInt16))
        # info("_ximc_read_resp: data = $data")
        crc_expected = ximc_crc16(data)
        (crc != crc_expected) && error("Checksum error, got 0x$(hex(crc)), expected 0x$(hex(crc_expected))")
        data::Vector{UInt8}
    else
        UInt8[]
    end
end

function _ximc_read_resp(io::IO, cmd_code::String, ::Type{T}) ::T where T
    types = T.types
    # println(types)
    n_data_bytes = sum(map(sizeof, types))
    data = _ximc_read_resp(io::IO, cmd_code::String, n_data_bytes)
    data_buf = IOBuffer(data)
    # println(map(U -> ximcio_to_jl(read(data_buf, U)), types)...)
    # println("rrg")
    T(map(U -> ximcio_to_jl(read(data_buf, U)), types)...)
end

export XIMCMotor

mutable struct XIMCMotor <: Device
    io
    # io::Lockable{IO}
end

XIMCMotor(hostname::AbstractString, port::Integer) =
    XIMCMotor(connect(hostname, port))

# XIMCMotor(hostname::AbstractString, port::Integer) =
#     XIMCMotor(Lockable{IO}(connect(hostname, port)))
#

function ximc_cmd(motor::XIMCMotor, cmd_code::String)
    _ximc_send_cmd(motor.io, cmd_code)
    _ximc_read_resp(motor.io, cmd_code, 0)
    nothing
end
# ximc_cmd(motor::XIMCMotor, cmd_code::String) = map(motor.io) do io
#     _ximc_send_cmd(io, cmd_code)
#     _ximc_read_resp(io, cmd_code, 0)
#     nothing
# end

function ximc_cmd(motor::XIMCMotor, cmd_code::String, x)
    _ximc_send_cmd(motor.io, cmd_code, x)
    _ximc_read_resp(motor.io, cmd_code, 0)
nothing
end

# ximc_cmd(motor::XIMCMotor, cmd_code::String, x) = map(motor.io) do io
#     _ximc_send_cmd(io, cmd_code, x)
#     _ximc_read_resp(io, cmd_code, 0)
#     nothing
# end

# ximc_qry(motor::XIMCMotor, cmd_code::String, ::Type{T}) where T = map(motor.io) do io
#     _ximc_send_cmd(io, cmd_code)
#     _ximc_read_resp(io, cmd_code, T)
# end

function ximc_qry(motor::XIMCMotor, cmd_code::String, ::Type{T}) where T
    _ximc_send_cmd(motor.io, cmd_code)
    _ximc_read_resp(motor.io, cmd_code, T)
end

export XIMCDeviceInformation

mutable struct XIMCDeviceInformation
    manufacturer::XIMCString{4}
    manufacturer_id::XIMCString{2}
    product_description::XIMCString{8}
    hw_version_major::UInt8
    hw_version_minor::UInt8
    hw_release::UInt16
    reserved::XIMCReserved{12}
end

getindex(motor::XIMCMotor, ::Info) = ximc_qry(motor, "geti", XIMCDeviceInformation)

export XIMCDeviceStatus

mutable struct XIMCDeviceStatus
    move_state::UInt8
    move_command_state::UInt8
    power_state::UInt8
    encoder_state::UInt8
    winding_state::UInt8
    current_position::Int32
    u_current_position::Int16
    encoder_position::Int64
    current_speed::Int32
    u_current_speed::Int16
    engine_current::Int16
    supply_voltage::Int16
    usb_current::Int16
    usb_voltage::Int16
    current_temperature::Int16
    flags::UInt32
    gpio_flags::UInt32
    cmd_buf_free_space::UInt8
    reserved::XIMCReserved{4}
end

getindex(motor::XIMCMotor, ::Type{Status}) = ximc_qry(motor, "gets", XIMCDeviceStatus)

export XIMCPosition

mutable struct XIMCPosition
    pos::Int32
    u_pos::Int16
    enc_pos::Int64
    pos_flags::Int8
    reserved::XIMCReserved{5}
end

getindex(motor::XIMCMotor, ::Type{Position}) = ximc_qry(motor, "gpos", XIMCPosition)
setindex!(motor::XIMCMotor, x::XIMCPosition, ::Type{Position}) = ximc_cmd(motor, "spos", x)

const BORDER_IS_ENCODER = 0x1
const BORDER_STOP_LEFT = 0x2
const BORDER_STOP_RIGHT = 0x4
const BORDERS_SWAP_MISSET_DETECTION = 0x8

const ENDER_SWAP = 0x1
const ENDER_SW1_ACTIVE_LOW = 0x2
const ENDER_SW2_ACTIVE_LOW = 0x4

export XIMCEdgesSettings

mutable struct XIMCEdgesSettings
    border_flags::UInt8
    ender_flags::UInt8
    left_border::Int32
    uleft_border::Int16
    right_border::Int32
    uright_border::Int16
    reserved::XIMCReserved{6}
end

getindex(motor::XIMCMotor, ::Type{Edges}, ::Type{Settings}) = ximc_qry(motor, "geds", XIMCEdgesSettings)
setindex!(motor::XIMCMotor, x::XIMCEdgesSettings, ::Type{Edges}, ::Type{Settings}) = ximc_cmd(motor, "seds", x)

const HOME_DIR_FIRST = 0x01
const HOME_DIR_SECOND = 0x02
const HOME_MV_SEC_EN = 0x04
const HOME_HALF_MV = 0x08
const HOME_STOP_FIRST_BITS = 0x30
const HOME_STOP_FIRST_REV = 0x10
const HOME_STOP_FIRST_SYN = 0x20
const HOME_STOP_FIRST_LIM = 0x30
const HOME_STOP_SECOND_BITS = 0xC0
const HOME_STOP_SECOND_REV = 0x40
const HOME_STOP_SECOND_SYN = 0x80
const HOME_STOP_SECOND_LIM = 0xC0

export XIMCHomeSettings

mutable struct XIMCHomeSettings
    fast_home_spd::UInt32
    ufast_home_spd::UInt8
    slow_home_spd::UInt32
    uslow_home_spd::UInt8
    home_delta::Int32
    uhome_delta::Int16
    home_flags::UInt16
    reserved::XIMCReserved{9}
end

getindex(motor::XIMCMotor, ::Type{Home}, ::Type{Settings}) = ximc_qry(motor, "ghom", XIMCHomeSettings)
setindex!(motor::XIMCMotor, x::XIMCHomeSettings, ::Type{Home}, ::Type{Settings}) = ximc_cmd(motor, "shom", x)

const FEEDBACK_ENCODER = 0x1
const FEEDBACK_ENCODERHALL = 0x3
const FEEDBACK_EMF = 0x4
const FEEDBACK_NONE = 0x5

const FEEDBACK_ENC_REVERSE = 0x1
const FEEDBACK_HALL_REVERSE = 0x2

export XIMCFeedbackSettings
mutable struct XIMCFeedbackSettings
    ips::UInt16
    feedback_type::UInt8
    feedback_flags::UInt8
    hall_spr::UInt16
    hall_shift::Int8
    reserved::XIMCReserved{5}
end

getindex(motor::XIMCMotor, ::Type{Feedback}, ::Type{Settings}) = ximc_qry(motor, "gfbs", XIMCFeedbackSettings)
setindex!(motor::XIMCMotor, x::XIMCFeedbackSettings, ::Type{Feedback}, ::Type{Settings}) = ximc_cmd(motor, "sfbs", x)

export XIMCMoveSettings

mutable struct XIMCMoveSettings
    speed::UInt32
    u_speed::UInt8
    accel::UInt16
    decel::UInt16
    antiplay_speed::UInt32
    u_aptiplay_speed::UInt8
    reserved::XIMCReserved{10}
end

#=  03.04. modified integer sizes and types
type XIMCMoveSettings
    speed::UInt8                   +6
    u_speed::UInt8
    accel::Int32                   -4
    decel::Int16
    antiplay_speed::Int32
    u_aptiplay_speed::Int16        -2
    reserved::XIMCReserved{10}
end
=#

getindex(motor::XIMCMotor, ::Type{Move}, ::Type{Settings}) = ximc_qry(motor, "gmov", XIMCMoveSettings)
setindex!(motor::XIMCMotor, x::XIMCMoveSettings, ::Type{Move}, ::Type{Settings}) = ximc_cmd(motor, "smov", x)

const ENGINE_TYPE_NONE      = 0x00 # A value that shouldn't be used.
const ENGINE_TYPE_DC        = 0x01 # DC motor.
const ENGINE_TYPE_2DC       = 0x02 # 2 DC motors.
const ENGINE_TYPE_STEP      = 0x03 # Step motor.
const ENGINE_TYPE_TEST      = 0x04 # Duty cycle are fixed. Used only manufacturer.
const ENGINE_TYPE_BRUSHLESS = 0x05 # Brushless motor.

export XIMCEngineKind

mutable struct XIMCEngineKind
    engine_type::UInt8
    driver_type::UInt8
    reserved::XIMCReserved{6}
end

getindex(motor::XIMCMotor, ::Type{Engine}, ::Type{Kind}) = ximc_qry(motor, "gent", XIMCEngineKind)
#Remove comment to implement "sent"# setindex!(motor::XIMCMotor, x::XIMCEngineKind, ::Type{Engine}, ::Type{Kind}) = ximc_cmd(motor, "sent", x)

# Enginge Flags
const ENGINE_REVERSE        = 0x01 # Reverse positive direction of engine shaft rotation
const ENGINE_MAX_SPEED      = 0x04 # If set, maximum achievable speed with current settings is used as nominal speed
const ENGINE_ANTIPLAY       = 0x08 # Activate play compensation procedure
const ENGINE_ACCEL_ON       = 0x10 # Activate smooth acceleration-deceleration for movements
const ENGINE_LIMIT_VOLT     = 0x20 # Only for DC motor
const ENGINE_LIMIT_CURR     = 0x40 # Only for DC motor
const ENGINE_LIMIT_RPM      = 0x80 # Enable motor speed limit

#Microstep Modes - step size is 1/"FRAC_X" - ranges from full step mode to 1/256 step mode
const MICROSTEP_MODE_FULL       = 0x01
const MICROSTEP_MODE_FRAC_2     = 0x02
const MICROSTEP_MODE_FRAC_4     = 0x03
const MICROSTEP_MODE_FRAC_8     = 0x04
const MICROSTEP_MODE_FRAC_16    = 0x05
const MICROSTEP_MODE_FRAC_32    = 0x06
const MICROSTEP_MODE_FRAC_64    = 0x07
const MICROSTEP_MODE_FRAC_128   = 0x08
const MICROSTEP_MODE_FRAC_256   = 0x09

export XIMCEngineSettings

mutable struct XIMCEngineSettings
    nom_voltage::UInt16
    nom_current::UInt16
    nom_speed::UInt32
    u_nom_speed::UInt8
    engine_flags::UInt16
    antiplay::Int16
    microstep_mode::UInt8
    steps_per_rev::UInt16
    reserved::XIMCReserved{12}
end

getindex(motor::XIMCMotor, ::Type{Engine}, ::Type{Settings}) = ximc_qry(motor, "geng", XIMCEngineSettings)
setindex!(motor::XIMCMotor, x::XIMCEngineSettings, ::Type{Engine}, ::Type{Settings}) = ximc_cmd(motor, "seng", x)

(motor::XIMCMotor)(::Type{Stop}) = motor(Stop, Soft)
(motor::XIMCMotor)(::Type{Stop}, ::Type{Hard}) = ximc_cmd(motor, "stop")
(motor::XIMCMotor)(::Type{Stop}, ::Type{Soft}) = ximc_cmd(motor, "sstp")

(motor::XIMCMotor)(::Type{Loft}) = ximc_cmd(motor, "loft")  #Execute play compensation movement

#PowerFlags
const POWER_REDUCT_ENABLED  = 0x01
const POWER_OFF_ENABLED     = 0x02 # Power off after PowerOffDelay
const POWER_SMOOTH_CURRENT  = 0x04 # Smooth current ramp up/down during current_set_time

export XIMCPowerSettings

mutable struct XIMCPowerSettings
    hold_current::UInt8
    curr_reduct_delay::UInt16
    power_off_delay::UInt16
    current_set_time::UInt16
    power_flags::UInt8
    reserved::XIMCReserved{6}
end

getindex(motor::XIMCMotor, ::Type{Power}, ::Type{Settings}) = ximc_qry(motor, "gpwr", XIMCPowerSettings)
setindex!(motor::XIMCMotor, x::XIMCPowerSettings, ::Type{Power}, ::Type{Settings}) = ximc_cmd(motor, "spwr", x)

setindex!(motor::XIMCMotor, x::Bool, ::Type{Power}, ::Type{State}) =
    !x && ximc_cmd(motor, "pwof")

#Protection - Critical Parameter flags - enable conditions to enter ALARM state
const ALARM_ON_DRIVER_OVERHEATING   = 0x01
const LOW_UPWR_PROTECTION           = 0x02
const H_BRIDGE_ALERT                = 0x04
const ALARM_ON_BORDERS_SWAP_MISSET  = 0x08
const ALARM_FLAGS_STICKING          = 0x10 # If set, only STOP command can turn alarms to 0
const USB_BREAK_RECONNECT           = 0x20

export XIMCProtectionSettings

mutable struct XIMCProtectionSettings
    low_volt_pwr_off::UInt16
    critical_curr_pwr::UInt16
    critical_volt_pwr::UInt16
    critical_temp::UInt16
    critical_curr_usb::UInt16
    critical_volt_usb::UInt16
    minimum_volt_usb::UInt16
    flags::UInt8
    reserved::XIMCReserved{7}
end

getindex(motor::XIMCMotor, ::Type{Protection}, ::Type{Settings}) = ximc_qry(motor, "gsec", XIMCProtectionSettings)
setindex!(motor::XIMCMotor, x::XIMCProtectionSettings, ::Type{Protection}, ::Type{Settings}) = ximc_cmd(motor, "ssec", x)

export XIMCPIDSettings

mutable struct XIMCPIDSettings
    kp_u::UInt16
    ki_u::UInt16
    kd_u::UInt16
    reserved::XIMCReserved{36}
end

getindex(motor::XIMCMotor, ::Type{PID}, ::Type{Settings}) = ximc_qry(motor, "gpid", XIMCPIDSettings)
#Remove comment to implement "spid"# setindex!(motor::XIMCMotor, x::XIMCPIDSettings, ::Type{PID}, ::Type{Settings}) = ximc_cmd(motor, "spid", x)

#Brake flags
const BRAKE_ENABLED     = 0x01
const BRAKE_ENG_PWROFF  = 0x02   # Brake turns off power on step motor

export XIMCBrakeSettings

mutable struct XIMCBrakeSettings
    t1::UInt16              #time[ms] between motor pwr on and brake off
    t2::UInt16              #time[ms] between brake off and move ready
    t3::UInt16              #time[ms] between motor stop and brake on
    t4::UInt16              #time[ms] between brake on and motor pwr off
    brake_flags::UInt8
    reserved::XIMCReserved{10}
end

getindex(motor::XIMCMotor, ::Type{Brake}, ::Type{Settings}) = ximc_qry(motor, "gbrk", XIMCBrakeSettings)
setindex!(motor::XIMCMotor, x::XIMCBrakeSettings, ::Type{Brake}, ::Type{Settings}) = ximc_cmd(motor, "sbrk", x)

#Control Position Flags
const CTP_ENABLED           = 0x01
const CTP_BASE              = 0x02
const CTP_ALARM_ON_ERROR    = 0x04
const REV_SENS_INV          = 0x08
const CTP_ERROR_CORRECTION  = 0x10

export XIMCControlPositionSettings

struct XIMCControlPositionSettings
    ctp_min_error::UInt8
    ctp_flags::UInt8
    reserved::XIMCReserved{10}
end

getindex(motor::XIMCMotor, ::Type{ControlPosition}, ::Type{Settings}) = ximc_qry(motor, "gctp", XIMCControlPositionSettings)
setindex!(motor::XIMCMotor, x::XIMCControlPositionSettings, ::Type{ControlPosition}, ::Type{Settings}) = ximc_cmd(motor, "sctp", x)

#CtrlFlags
const EEPROM_PRECEDENCE = 0x01

export XIMCNameFRAMSettings

struct XIMCNameFRAMSettings
    controller_name::XIMCString{16}
    ctrl_flags::UInt8
    reserved::XIMCReserved{7}
end

getindex(motor::XIMCMotor, ::Type{NameFRAM}, ::Type{Settings}) = ximc_qry(motor, "gnmf", XIMCNameFRAMSettings)
setindex!(motor::XIMCMotor, x::XIMCNameFRAMSettings, ::Type{NameFRAM}, ::Type{Settings}) = ximc_cmd(motor, "snmf", x)

export XIMCNameStageSettings

struct XIMCNameStageSettings
    positioner_name::XIMCString{16}
    reservred::XIMCReserved{8}
end

getindex(motor::XIMCMotor, ::Type{NameStage}, ::Type{Settings}) = ximc_qry(motor, "gnme", XIMCNameStageSettings)
setindex!(motor::XIMCMotor, x::XIMCNameStageSettings, ::Type{NameStage}, ::Type{Settings}) = ximc_cmd(motor, "snme", x)

getindex(motor::XIMCMotor, ::Type{Position}, ::Type{SetZero}) =
    ximc_cmd(motor, "zero")

#(motor::XIMCMotor)(::Type{Position}, ::Type{SetZero}) =
#    ximc_cmd(motor, "zero")

struct XIMCCommandMove
    position::Int32
    u_position::Int16
    reserved::XIMCReserved{6}
end

XIMCCommandMove(p, u_p = 0) = XIMCCommandMove(p, u_p, XIMCReserved(Tuple(zeros(UInt8, 6))))

setindex!(motor::XIMCMotor, x::Real, ::Type{Move}, ::Type{Absolute}) =
    ximc_cmd(motor, "move", XIMCCommandMove(x))

setindex!(motor::XIMCMotor, x::Symbol, ::Type{Move}, ::Type{Absolute}) = begin
    if x == :home
        ximc_cmd(motor, "home")
    else
        error("Invalid target position $x")
    end
end

struct XIMCCommandMovr
    delta_position::Int32
    u_delta_position::Int16
    reserved::XIMCReserved{6}
end

XIMCCommandMovr(p, u_p = 0) = XIMCCommandMovr(p, u_p, XIMCReserved(Tuple(zeros(UInt8, 6))))

setindex!(motor::XIMCMotor, x::Real, ::Type{Move}, ::Type{Relative}) =
    ximc_cmd(motor, "movr", XIMCCommandMovr(x))


text:
using Serialization

const ATTACHMENT_BUFFER = Ref{Vector{Expr}}(Expr[])
const LAST_ATTACHMENT_IDX = Ref{Int}(0)

save_attachment(ex::Expr) = push!(ATTACHMENT_BUFFER[], ex)
get_attachments() = ATTACHMENT_BUFFER[]
clear_attachments!() = ATTACHMENT_BUFFER[] = Expr[]

function generate_temp_filename(prefix::AbstractString, suffix::AbstractString)
    LAST_ATTACHMENT_IDX[] = LAST_ATTACHMENT_IDX[] + 1
    attachment_idx = lpad(LAST_ATTACHMENT_IDX[], 8, '0')
    prefix * attachment_idx * suffix
end

function convert_to_uri(path::AbstractString)
    # add_file wants URI format in order to work in local mode
    # and it does seem to work unless we have three slashes at the statr
    if Sys.iswindows()
        "file:///" * replace(path, '\\' => '/')
    elseif startswith(path, "/")
        # absolute path
        "file://" * path
    else
        # relative path
        # "file:./" * path
        "file://" * abspath(path)
    end
end

function process_attachments(sc::SparkContext)
    for ex in get_attachments()
        temp_filename = generate_temp_filename("attached_", ".jl")
        path = joinpath(get_temp_dir(sc), temp_filename)
        open(path, "w") do io
            write(io, string(ex))
        end

        add_file(sc, convert_to_uri(path))

        # if it's `include` expression, also attach included file
        if isa(ex, Expr) && ex.head == :call && ex.args[1] == :include
            add_file(sc, convert_to_uri(ex.args[2]))
        end
    end
    clear_attachments!()
end

macro attach(ex)
    save_attachment(ex)
    esc(ex)
end

"Makes the value of data available on workers as symbol name"
function share_variable(sc::SparkContext, name::Symbol, data::Any)
    temp_filename = generate_temp_filename("sharevar_", ".jbin")
    path = joinpath(get_temp_dir(sc), temp_filename)
    open(path, "w") do io
        serialize(io, data)
    end
    add_file(sc, convert_to_uri(path))
    ex = quote
            using Serialization
            open($temp_filename, "r") do io
                global $name = deserialize(io)
            end
         end
    save_attachment(ex)
    data
end

"Makes the variable available on workers"
macro share(sc, var)
    q = Expr(:quote, var)
    :(share_variable($(esc(sc)), $q, $(esc(var))))
end

ext: jl


text:
using Utilities,StaticNets,JLD,MLBase
using PyCall; pygui(:qt); using PyPlot

## Load data
#whole period data file
loadFileName = "/home/Domenico/Dropbox/Dynamic_Networks/data/emid_data/juliaFiles/Weekly_eMid_Data_from_2009_06_22_to_2015_02_27.jld"
@load(loadFileName, AeMidWeekly_Active_T,banksIDs,inactiveBanks, YeMidWeekly_Active_T,weekInd,datesONeMid ,degsIO_T,strIO_T)

## Methods definitions

struct SnapSeqNetDirBin1
    obsT:: Array{<:Real,2}  # In degrees and out Degrees for each t . The last
                            # dimension is always time. 2NxT
    Par::Array{<:Real,2} # In and Out Parameters [InW_T;outW_T] for each t
end
SnapSeqNetDirBin1(obs::Array{<:Real,2}) = SnapSeqNetDirBin1(obs,zeros(obs))
function estimate(Model::SnapSeqNetDirBin1; degsIO_T::Array{<:Real,2}=Model.obsT,
                        targetErr::Real=1e-12,identPost=false,identIter=false,zeroDegFit="small")
   N2,T = size(degsIO_T);N = round(Int,N2/2)
   UnParT = zeros(2N,T)
   prog = 0
   for t = 1:T
       degs_t = degsIO_T[:,t]
       UnPar_t , Niter  = estimate(fooErgmDirBin1; degIO = degs_t , targetErr =  targetErr,identIter=identIter)
       identPost?UnParT[:,t] = identify!(fooErgmDirBin1,UnPar_t): UnParT[:,t] = UnPar_t
       round(t/T,2)>prog ? (prog=round(t/T,2);println((prog,Niter)) ):()
   end
   if uppercase(zeroDegFit) == "INF"
       # zero degrees have -Inf fitnesses
   elseif uppercase(zeroDegFit) == "10AS21"
       # such that the difference  fit(2) - fit(1) == fit(1)-fit(0) , where
       # fit(n) means the fitness associated with a node of deg n.
       # when more than one choice is possible (the fitness depends on both degrees)
       # than choose fit(2) and fit(1) that maximize the difference

   elseif uppercase(zeroDegFit) == "SMALL"
       # such that the probability of observing a link is very small
       UnParT[degsIO_T.==0] = -5
   end
 return UnParT
end
function bigNegParFun(Model::ErgmDirBin1,N::Int; smallProb = 1e-1)
    #define a number big enough to play the role of Inf for
    # purposes of sampling N(N-1) bernoulli rvs with prob 1/(1+exp(  bigNumb))
    # prob of sampling degree > 0 (N-1)/(1+exp(bigNumb)) < 1e-6
    bigNumb = - log((N-1)/smallProb - 1)
    return bigNumb
 end
## Application to Data
snapModData = SnapSeqNetDirBin1(degsIO_T)
singSnapEstNoIdent_T = estimate(snapModData)
singSnapEstIdent_T = estimate(snapModData; identPost = true,identIter=false)
singSnapEstIdentIter_T = estimate(snapModData; identPost = false,identIter=true)
N2,T = size(degsIO_T);N = round(Int,N2/2)

# plot(singSnapEstNoIdent_T[1:30,:]')
# plot(singSnapEstIdent_T[1:3,:]')
# plot(singSnapEstIdentIter_T[1:3,:]')
## check that equal degrees are associated to equal fitnesses.
using StatsBase
singSnapEst_T = singSnapEstIdent_T
fit,degs = [singSnapEst_T[1:N,1] singSnapEst_T[1+N:2N,1]],[degsIO_T[1:N,1] degsIO_T[1+N:2N,1]]
[unique(degs,1) unique(fit,1)]
# actually is a pair (InDeg,OutDeg) that is associated with a pair (InFit,OutFit). That is because
# there cannot be single loops. Hence the value of a fitness has to take into account
# that links cannot come from the node itself. So all things equal a node with
# (1,1) should have an higher  InFit than a (1,0) because it needs to attract more
# from the other ones
##
D = maximum(degsIO_T)
storeFitDistrPerDeg =  Matrix{Vector{Float64}}(D,T)
meansFit = Matrix{Float64}(D,T)
stdFit = Matrix{Float64}(D,T)
testSingSnapEst_T = singSnapEst_T
t=1;d=1
for t=1:T
    degsIOt = degsIO_T[:,t]
    degsUniIOt = sort(unique(degsIOt))
    degsUniIOt = degsUniIOt[2:end]
    println(length(degsUniIOt))
    for d=1:length(degsUniIOt)
        deg = degsUniIOt[d]
        indsDeg = degsIOt[1:N].== deg
        storeFitDistrPerDeg[d,t] = testSingSnapEst_T[1:N,t][indsDeg]
        meansFit[d,t] = mean(testSingSnapEst_T[1:N,t][indsDeg])
        stdFit[d,t] = std(testSingSnapEst_T[1:N,t][indsDeg])
    end

end

deg = 1:4
plot(meansFit[deg,:]',"-");plot(meansFit[deg,:]' + stdFit[deg,:]',"--");plot(meansFit[deg,:]' - stdFit[deg,:]',"--")
tmp = 1./sumSq(degsIO_T,1)
cor(meansFit[2,:],tmp)
mean(meansFit[4,:])
#plot(tmp)
# the mean value of the fitnesses associated to a given degree is highly correlated
# with the inverse of total number of links

## lets estimate some AR(1)

using AReg

function estManyAR(obs::Array{<:Real,2};p::Int=1)
    #obs is a matrix with a time series for each row (time is last dimension)
    D,T = size(obs)
    isConst = ([length(unique(obs[i,:])) for i=1:D ] .==1)
    parEst = zeros(D,p+1)
    arp = AReg.ARp(ones(p+1),.1,[])
    for d=1:D
        if ~isConst[d]
            parEst[d,:] = AReg.fitARp(obs[d,:],arp).estimates
        end
        #println(d)
    end
    return parEst,isConst
end
function oneStepForAR1(obs_t::Float64,AR1Par::Vector{<:Real})
    #given w and β AR1 parameters, compute the one step ahead forecast from obs_t
    obstp1 = AR1Par[1] + AR1Par[2] * obs_t
    return obstp1
end
#Estimate AR1 on the training sample
function forecastEvalAR1(allFit::Array{Float64,2})
    N2,T = size(allFit)
    Ttrain = round(Int,T/2)
    Ttest = T-Ttrain
    trainFit = allFit[:,1:Ttrain]
    testFit = allFit[:,1:Ttest]
    estAR1,isCon = estManyAR(trainFit)
    ## One step ahead forecasts for each fitness for each time in Test sample
    foreFit = zeros(testFit[:,1:end])
    for n=1:2N,t=1:Ttest-1
        foreFit[n,t+1] = oneStepForAR1(testFit[n,t],estAR1[n,:])
    end

    foreVals = 100 * ones(N*(N-1)*(Ttest-1))
    realVals = trues(N*(N-1)*(Ttest-1))
    noDiagInd = .!Bool.(eye(N))
    lastInd=1
    NmaxLinks = N*(N-1)
    Nobs = NmaxLinks*T
    for t=2:Ttest
        adjMat =  AeMidWeekly_Active_T[:,:,Ttrain+1:end][:,:,t]
        expMat = expMatrix(fooErgmDirBin1,foreFit[:,t])
        foreVals[lastInd:lastInd+NmaxLinks-1] = expMat[noDiagInd]
        realVals[lastInd:lastInd+NmaxLinks-1] = adjMat[noDiagInd]

        lastInd += NmaxLinks
    end
    return realVals,foreVals
end
##

forecastEvalAR1(allFit)
function rocCurve(realVals::BitArray{1},foreVals::Vector{Float64}; Th::Vector{Float64}=Vector(0:0.00001:1))
    rocOut = roc(realVals,foreVals,Th)
    Nth = length(Th)
    tpr = [rocOut[i].tp/rocOut[i].p for i=1:Nth]
    fpr = [rocOut[i].fp/rocOut[i].n for i=1:Nth]
    auc = sum(-diff(fpr).*tpr[1:end-1])
    plot(fpr,tpr)
    xlim([0,1])
    ylim([0,1])
    title("AUC =  $(round(auc,3))")
    grid()
    return tpr, fpr,auc
end
rocCurve(allFit::Array{Float64,2}) = (tmp = forecastEvalAR1(allFit); rocCurve(tmp[1],tmp[2]))
rocCurve(realVals,foreVals)

rocCurve(singSnapEstIdent_T)
rocCurve(singSnapEstNoIdent_T)
rocCurve(singSnapEstIdentIter_T)

#

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
using Test, Arrow, Tables, Dates

struct MapTable
    x
end
Tables.columnnames(x::MapTable) = propertynames(x.x)
Tables.getcolumn(x::MapTable, i::Int) = getfield(x.x, i)
Tables.getcolumn(x::MapTable, nm::Symbol) = getproperty(x.x, nm)
Tables.schema(x::MapTable) = Tables.Schema(propertynames(x.x), eltype.(getproperty(x.x, nm) for nm in propertynames(x.x)))
Tables.columns(x::MapTable) = x

@testset "Arrow" begin

# basic
t = (col1=Int64[1,2,3,4,5,6,7,8,9,10],)
io = IOBuffer()
Arrow.write(io, t)
seekstart(io)
tt = Arrow.Table(io)
@test length(tt) == length(t)
@test tt.col1 == t.col1
@test eltype(tt.col1) === Int64

# missing values
t = (col1=Union{Int64, Missing}[1,2,3,4,5,6,7,8,9,missing],)
io = IOBuffer()
Arrow.write(io, t)
seekstart(io)
tt = Arrow.Table(io)
@test length(tt) == length(t)
@test isequal(tt.col1, t.col1)
@test eltype(tt.col1) === Union{Int64, Missing}

# primitive types
t = (
    col1=[missing, missing, missing, missing],
    col2=Union{UInt8, Missing}[0, 1, 2, missing],
    col3=Union{UInt16, Missing}[0, 1, 2, missing],
    col4=Union{UInt32, Missing}[0, 1, 2, missing],
    col5=Union{UInt64, Missing}[0, 1, 2, missing],
    col6=Union{Int8, Missing}[0, 1, 2, missing],
    col7=Union{Int16, Missing}[0, 1, 2, missing],
    col8=Union{Int32, Missing}[0, 1, 2, missing],
    col9=Union{Int64, Missing}[0, 1, 2, missing],
    col10=Union{Float16, Missing}[0, 1, 2, missing],
    col11=Union{Float32, Missing}[0, 1, 2, missing],
    col12=Union{Float64, Missing}[0, 1, 2, missing],
    col13=[true, false, true, missing],
)
io = IOBuffer()
Arrow.write(io, t)
seekstart(io)
tt = Arrow.Table(io)
@test length(tt) == length(t)
@test all(isequal.(values(t), values(tt)))

t = (
    col14=[zero(Arrow.Decimal{Int32(2), Int32(2)}), zero(Arrow.Decimal{Int32(2), Int32(2)}), zero(Arrow.Decimal{Int32(2), Int32(2)}), missing],
    col15=[zero(Arrow.Date{Arrow.Meta.DateUnit.DAY, Int32}), zero(Arrow.Date{Arrow.Meta.DateUnit.DAY, Int32}), zero(Arrow.Date{Arrow.Meta.DateUnit.DAY, Int32}), missing],
    col16=[zero(Arrow.Time{Arrow.Meta.TimeUnit.SECOND, Int32}), zero(Arrow.Time{Arrow.Meta.TimeUnit.SECOND, Int32}), zero(Arrow.Time{Arrow.Meta.TimeUnit.SECOND, Int32}), missing],
    col17=[zero(Arrow.Timestamp{Arrow.Meta.TimeUnit.SECOND, nothing}), zero(Arrow.Timestamp{Arrow.Meta.TimeUnit.SECOND, nothing}), zero(Arrow.Timestamp{Arrow.Meta.TimeUnit.SECOND, nothing}), missing],
    col18=[zero(Arrow.Interval{Arrow.Meta.IntervalUnit.YEAR_MONTH, Int32}), zero(Arrow.Interval{Arrow.Meta.IntervalUnit.YEAR_MONTH, Int32}), zero(Arrow.Interval{Arrow.Meta.IntervalUnit.YEAR_MONTH, Int32}), missing],
    col19=[zero(Arrow.Duration{Arrow.Meta.TimeUnit.SECOND}), zero(Arrow.Duration{Arrow.Meta.TimeUnit.SECOND}), zero(Arrow.Duration{Arrow.Meta.TimeUnit.SECOND}), missing],
)
io = IOBuffer()
Arrow.write(io, t)
seekstart(io)
tt = Arrow.Table(io; convert=false)
@test length(tt) == length(t)
@test all(isequal.(values(t), values(tt)))

# list types
t = (
    col1=Union{String, Missing}["hey", "there", "sailor", missing],
    col2=Union{Vector{UInt8}, Missing}[b"hey", b"there", b"sailor", missing],
    col3=Union{Vector{Int64}, Missing}[Int64[1], Int64[2], Int64[3], missing],
    col4=Union{NTuple{2, Vector{Int64}},Missing}[(Int64[1], Int64[2]), missing, missing, (Int64[3], Int64[4])],
    col5=Union{NTuple{2, UInt8}, Missing}[(0x01, 0x02), (0x03, 0x04), missing, (0x05, 0x06)],
    col6=NamedTuple{(:a, :b), Tuple{Int64, String}}[(a=Int64(1), b="hey"), (a=Int64(2), b="there"), (a=Int64(3), b="sailor"), (a=Int64(4), b="jo-bob")],
)
io = IOBuffer()
Arrow.write(io, t)
seekstart(io)
tt = Arrow.Table(io)
@test length(tt) == length(t)
@test all(isequal.(values(t), values(tt)))

# unions
t = (
    col1=Arrow.DenseUnionVector([1, 2.0, 3, 4.0, missing]),
    col2=Arrow.SparseUnionVector([1, 2.0, 3, 4.0, missing]),
)
io = IOBuffer()
Arrow.write(io, t)
seekstart(io)
tt = Arrow.Table(io)
@test length(tt) == length(t)
@test all(isequal.(values(t), values(tt)))

# dict encodings
t = (
    col1=Arrow.DictEncode(NamedTuple{(:a, :b), Tuple{Int64, Union{String, Missing}}}[(a=Int64(1), b=missing), (a=Int64(1), b=missing), (a=Int64(3), b="sailor"), (a=Int64(4), b="jo-bob")]),
)
io = IOBuffer()
Arrow.write(io, t; debug=true)
seekstart(io)
tt = Arrow.Table(io; debug=true)
@test length(tt) == length(t)
@test all(isequal.(values(t), values(tt)))

# multiple record batches
if isdefined(Tables, :partitioner)
    xx = Tables.partitioner
else
    xx = Tuple
end
t = xx(((col1=Union{Int64, Missing}[1,2,3,4,5,6,7,8,9,missing],), (col1=Union{Int64, Missing}[1,2,3,4,5,6,7,8,9,missing],)))
io = IOBuffer()
Arrow.write(io, t)
seekstart(io)
tt = Arrow.Table(io)
@test length(tt) == 1
@test isequal(tt.col1, vcat([1,2,3,4,5,6,7,8,9,missing], [1,2,3,4,5,6,7,8,9,missing]))
@test eltype(tt.col1) === Union{Int64, Missing}

# auto-converting types
t = (
    col1=[Date(2001, 1, 2), Date(2010, 10, 10), Date(2020, 12, 1)],
    col2=[Time(1, 1, 2), Time(13, 10, 10), Time(22, 12, 1)],
    col3=[DateTime(2001, 1, 2), DateTime(2010, 10, 10), DateTime(2020, 12, 1)]
)
io = IOBuffer()
Arrow.write(io, t; debug=true)
seekstart(io)
tt = Arrow.Table(io; debug=true)
@test length(tt) == length(t)
@test all(isequal.(values(t), values(tt)))

# Map
t = MapTable((
    col1=Dict(Int32(1) => Float32(3.14)),
))
io = IOBuffer()
Arrow.write(io, t; debug=true)
seekstart(io)
tt = Arrow.Table(io; debug=true)
for (k, v) in tt.col1
    @test isequal(t.x.col1[k], v)
end

# file format
t = (
    col1=[missing, missing, missing, missing],
    col2=Union{UInt8, Missing}[0, 1, 2, missing],
    col3=Union{UInt16, Missing}[0, 1, 2, missing],
    col4=Union{UInt32, Missing}[0, 1, 2, missing],
    col5=Union{UInt64, Missing}[0, 1, 2, missing],
    col6=Union{Int8, Missing}[0, 1, 2, missing],
    col7=Union{Int16, Missing}[0, 1, 2, missing],
    col8=Union{Int32, Missing}[0, 1, 2, missing],
    col9=Union{Int64, Missing}[0, 1, 2, missing],
    col10=Union{Float16, Missing}[0, 1, 2, missing],
    col11=Union{Float32, Missing}[0, 1, 2, missing],
    col12=Union{Float64, Missing}[0, 1, 2, missing],
    col13=[true, false, true, missing],
)
io = IOBuffer()
Arrow.write(io, t; file=true)
seekstart(io)
tt = Arrow.Table(io)
@test length(tt) == length(t)
@test all(isequal.(values(t), values(tt)))

end


text:
"""
    About GLOQ.jl
"""
function about()
    println("    ...`.`.`.`````.`````...`.`..``...`.`........`....................................................---
    `.````````.```.`````````.`..``.```.......`................`.......................................--
    .`````.`.``..```````````````..`............``............`.......................................---
    `````````````````````````.`.........----..........................................................--
    ````````````````````````...-/++oo+oyhhhhyyyssooosososssoso++osss+:................................--
    ``````````````````````...-/shddddddddddddddddddmmdmdddddddddddddddy+........`....``...............--
    ``````````````````````.:syhddmmdmmmmmmmmmmmmmmmmmmmmmmmmdmmmddmddmdds......`.```.``...............--
    ````````````````````.-ohhhddddmmdmmmmmmmmmmmmmmmmmmmmmmdddmddddddmmddo-.....``.`.```...............-
    ```````````````````./sshdddmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmdho/:-....`````````...........-
    ``````````````````.+yyhddmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmdys+:..``````.````..........
    `````````````````./hhhhddmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmdddho/-.``````.```.`.......
    `````````````````:hhhhhddmdmmdddddmmmmmddddddmmmmmmmmmmmmmdmdmdddddddmmmmddmmdds:..``````````.......
    `````````````````ohhhhhddddddhhhhhdhhhhhhhhddddddddddddddhhhhhhhhhhhhhddmmddmmmmdo-.``````````......
    ````````````````-yhhhhhhddddysssssyssysssyyyyyyyyyyyyyyyyyyyyyyyyyyyyyhddmmdmmmmmmdy-``````````.....
    ``````````````./yhhhhhhhddhysoo++ooo+++o++oooosoooooooossssssssssssyyyyhdmmmdmmmmmmmy.``````````.`..
    `````````````-sdddhhhhdddhyo+++++++++++/+++++++++++++++o+ooooooossssyyyyhdmmmmmmmmmmd/```````````...
    ````````````-yddddhhhddddyo+++++++++//////////////+++/+++++++ooooossssyyyhdmmdmmmmmmho````````````..
    ````````````ydddddhhhdddhs+++++++++////////////////++/+++++++++oooossssyyhhddddmmmmmdo`````````````.
    ```````````.hdddddhhhhhhyo+++++//+//////////////////////++++++++ooosssssyyhhddhdmmmmd/````````````..
    ```````````-ddddddhhyyyys++++++////////////////////+////++++++++ooosssssyyhhhdhdmmmmd/`````````````.
    ` `````````.yddddddhhyyso++++++/////////////////////////+//+++++ooosssssyyhhhhhhdmmmdo`````````````.
    ````````````-sddddddhyyso++++++///////////////////////////+/++++oooosssyyyyyyhdhdmmmd/``````````````
    ````` ```````/hhdhhhyysoo++++++////////////////////+//////+++++oooooossyyyyyyhhhddmdd:``````````````
    ```` ````````/hhhhhyysso+////+/+///////////////////++////++++ooooooosssyyyyyyyhhddddy.``````````````
      ` ```` ````-hhyhysso++/++//+++/+//////////////////++//+++++oosssssyhhhhhyyyyyyhhhh+.``````````````
        ` ````` `.shhyssoo+++/++++++ooooooo++o+++++++++++o+++ooossyyyhhhdddddhhyyyyyyhhhy/``````````````
      ` ` ` `` `./+syyosso+++//++++ossssyyyyyyyyysssysooosssyyyyyyyhhhhhhyhhhhhyhyyyyhyhy.``````````````
         `  ` ```/o++osssso++//++oooooossyyyyyyyyyyyyysoosyhhhhhhhhyyyyyyyyyhhhyyyhysyyh/```````````````
          ` `  ` `+o+ooooooo+++sooosssssssyyhhhhyyyyoysooosyhhhhhhhhdmmdhddhhhhyyyhyyyyy`````` ` ```````
               ` `:+++o+oo+++ssyosssyyhyyhddhyyyysossyo//++shhyyyyssyyyyyyyhhyyssyyyyys:`` `````` ``````
              `  ``/++++o++///+s+oossoooooooossssooshy+++++oydhyyysssosssyyyyssssshyyyo`` ` ``  ``` ````
                   :+ooo++++///s//++++++++oooosso++yys++/++oshhyysssssssssssosssssyyys/`   ``    `` ````
                  `-+oo++++++//+o///++++++++++++++oso+++/++osyhhysoooooooooossssssysss-`  `        `````
                   `/++++o+++///+o+++///////////+os+++++++++osyhhhsoooooooosssysssyys+`           ` ` ``
                   `.///oo+++++///ooo+///:////++sso+++++++++ooyhhyssssssssyysssssyyss+`      ``      ```
                     ://+oo++++//////+oooooo+++ooo++++++///+ooyyyyysooooooooossssyyss+        `     ` ``
                     .//+o++++++/////////////+ooo++++oo++++oosyhyyyysoooooooosssyyyss+`             ` ``
                     `/+oo+++++++//////////++o+///+ossssoossyhhdhhyyyssooooosssyyyyss/`              ` `
                      ./oo+++++++/+/////+++o++////++oooossyyhhhhhyyyssssoooosssyyyyss-`              ` `
                      `-+oo+++++++++///+++++//+++++++++ooosssssyyyyyyssssssssssyyyyso`               ```
                       ``/oo+o+++++++++++++++++++++++++ooooooossssyyyyyysooosssyyyysy:.`              ``
                    ``.-+sooooo++++++++///+++ooo+++++++++++ooosssyyhhyssoooossyyyyyhmmhs:.`          ` `
              ```.-/oyhdddhsooo+++++++///+++oooosssssssssssyyyyyhhhyysoo++oosyyyhhhmmmmmmhs:.``       `
       `  ``.:+syhdddddmmmmhsooo+++++++//++////+////++++++oooooossssoooo+oosyyyhhhhmNNNmmmmmds+-````  ``
    ```.:+syhddddddddmmmmmmmhssooo++++++/+///////++++ooooooosssssssoossooosyyyhhhhhhNNNNNNNmmmmdyo:.````
    +syddddddddddmmmmmmmmmmhysssoo+++++++++///////+++ooooosooooooooosssssssyhhhhhhhhdNNNNNNNNmmmmmmhs+-`
    dddddddddddmdmmmmmmmmmhssoosssoooo++++++////////+++++oooo+oooooossssyyyyhhdhhhhhhmNNNNNNNNNmmmmmdddy
    dddddmmmmmmmmmmmmmmmmyssooosssoooooo++++++///////+++++++ooooooosssyyyhhhhdddhhhddmNNNNNNNNNNNNmmmmmm
    mdmmmmmmmmmmmmmmmmmdsssoooosysssooooooo+++++/+////++++++ooosoosyyyyhhhhddhdhhhddddNNNNNNNNNNNNNNNNNN
    mdmmmmmmmmmmmmmmmmdsssooooosyyssssssssoooo++++++++oooooooossssyyhhhdhdddhddhhdddddNNNNNNNNNNNNNNNNNN
    mmmmmmmmmmmmmmmmmhsoooooooosyyosyyysssssssoooooooosssosossyyyyhhhdddddddhdhhddddddNNNNNNNNNNNNNNNNNN
    mmmmmmmmmmmmmmmdyooooooo++osyysosyyyyyyyyyyyyyssyyysssyyyhhhhhdddddddddhhhhhddddddmNNNMNMNNNMNNNNMMM
    mmmmmmmmmmmmmmdsoooooooo+ooyyyyoossyyyyyyyyyyyhyhhhhhhhhhhddhdhdhhhddddhhhhhdddddddNNMNNMMNNMNNMNMMM
    mmmmmmmmmmmmdhsooo+oooo+++osyyyoooosyyyyyyyyyyyyyhhhhhhhhhhhhdhhhhdddddhhhhdddddddhNNMMNMNNNMNNMMMMM
    mmmmmmmmmmmhsooooooooo++++osssss++oossyyyyyhyhhhhhhhhhhhhhhhhhhhhdddddhhhhhhddddhhhmNNMNMMMNMMNMMMMM
    mmmmmmmmmdyoooooooo+++++++ossysyo++ooosssyyyyhhhhhhhhhhhhhhhhhhhdddddhhhhhhhhhhhhyyhNNNMMMMMMMMNMMMM
    mmmmmmmdhsooooooo+oo++++++ossysss++++ooosssyyyyyyyyhhhhhhhhyyhhddddhhhhhhhhhhhyyyyyymNMNMNMNMNNMMMMM
    mmmmmdyssoooooo+++++++++++osssssyso++++ooooossssssyyyyyyyyyyyhdhhhhyyyhhyyyyyyyyyyyyhNMNMNMNMMNNMMMM
    mmmdysoooooooooooooo+++++oosssssyhso+++oooooossssssssssssyyyyhysyyysyyhyyyyyyyyyyssssmNNNMNNNMNMNMMM
    mmmdhhhyyssssoooooooooooosysossssyhso++oooooooooooooooosssyhyssssysssyyssssssysssssssyNNMMMNNNMNMMNM
    mmmmmmmmmmmmdddhhhhhhhddmmmhssssssyys+++++ooooooooo+ooossyhysooossooyyysssssssssssoooodNNNNNMNNMNNMM
    mmmmmmmmmmmmmmmmmmmmmmmmmmddossosssyyo++++ooo++++++++oosyysooooosooosyssosoooooooooooosNNNNNNNNNNMMM
    mmmmmmmmmmmmmmmmmmmmmmmmmmhyoosssssysso+++++++++////++oossooooosdhysssssooooooooooooooosNNMMNNNNNNMM
    mmmmmmmmmmmmmmmmmmmmmmmmmmhysssosssyssso++///////////+oossysyhdmmmNNmddhhyyyyhhhhhhhyyssdMNMNNNMNMMM")
    return nothing
end


text:
@testset "reflectometryNN.jl" begin

    # Write your tests here.
end


text:
using Radixsort
using Test

@testset "Radixsort.jl" begin
    # Write your tests here.
end


text:
using Documenter, InvariantSets

makedocs(;
    modules=[InvariantSets],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/ueliwechsler/InvariantSets.jl/blob/{commit}{path}#L{line}",
    sitename="InvariantSets.jl",
    authors="Ueli Wechsler",
    assets=String[],
)

deploydocs(;
    repo="github.com/ueliwechsler/InvariantSets.jl",
)

max_stars_repo_name: mforets/InvariantSets

max_issues_repo_name: mforets/InvariantSets


text:
using BenchmarkTools, JetPackDSP, Jets, LinearAlgebra

const SUITE = BenchmarkGroup()

s = JetSpace(Float64, 50, 51, 52)
n = 5
A = JopConvolve(s, smoother=:rect, n=(n,n,n))
m = rand(domain(A))
d = rand(range(A))

SUITE["JopConvolve"] = BenchmarkGroup()
SUITE["JopConvolve"]["construct smoother"] = @benchmarkable JopConvolve($s; smoother=:rect, n=($n,$n,$n))
SUITE["JopConvolve"]["forward"] = @benchmarkable mul!($d, $A, $m)
SUITE["JopConvolve"]["adjoint"] = @benchmarkable mul!($m, ($A)', $d)

SUITE


text:
function mst_prim_test()
    A = load_matrix_network("airports")
    A = -A
    T = mst_prim_matrix(A)

    A = load_matrix_network("clr-24-1")
    A[2,3] = 9
    A[3,2] = 9
    T = mst_prim_matrix(A)
    Ttrueijv = [
     2     8     1     4     6     9     3     5     4     3     7     6     8     1     7     3
     1     1     2     3     3     3     4     4     5     6     6     7     7     8     8     9
     4     8     4     7     4     2     7     9     9     4     2     2     1     8     1     2 ]
    Ttrue = sparse(vec(Ttrueijv[1,:]),  vec(Ttrueijv[2,:]), vec(Ttrueijv[3,:]), 9,9)

    if nnz(T - Ttrue) != 0
        error("mst_prim failed")
    end

    A = load_matrix_network("clr-24-1")
    T1 = mst_prim_matrix(A)
    T2 = mst_prim_matrix(A,false,5)
    T1T2diff = sparse([2,1],[3,8],[8,-8],9,9)

    if nnz(triu(T1-T2) - T1T2diff) != 0
        error("mst_prim failed")
    end

    return true
end


text:
using Base.Test
using QuantumOptics.sparsematrix

# SparseMatrix = quantumoptics.sparsematrix.SparseMatrix
typealias SparseMatrix SparseMatrixCSC{Complex128}

# Set up test matrices
A = rand(Complex128, 5, 5)
A_sp = sparse(A)

B = eye(Complex128, 5)
B_sp = speye(Complex128, 5)

C = rand(Complex128, 3, 3)
C[2,:] = 0
C_sp = sparse(C)

R_sp = A_sp + B_sp
R = A + B

# Test arithmetic
@test_approx_eq 0. norm(full(R_sp) - R)
@test_approx_eq 0. norm(full(Complex128(0.5,0)*A_sp) - 0.5*A)
@test_approx_eq 0. norm(full(A_sp/2) - A/2)
@test_approx_eq 0. norm(full(A_sp*B_sp) - A*B)

# Test kronecker product
@test_approx_eq 0. norm(full(kron(A_sp, C_sp)) - kron(A, C))
@test_approx_eq 0. norm(full(kron(A_sp, B_sp)) - kron(A, B))


text:
"""
### print a Mera timestamp on the screen if the global variable: verbose_mode == true
```julia
function printtime(text::String="", verbose::Bool=verbose_mode)
```
"""
function printtime(text::String="", verbose::Bool=verbose_mode)
    if verbose
        printstyled( "[Mera]: $text",now(), "\n", bold=true, color=:normal)
        println()
    end
end

function printtablememory(data, verbose::Bool)
    if verbose
        arg_value, arg_unit = usedmemory(data, false)
        println("Memory used for data table :", arg_value, " ", arg_unit)
        println("-------------------------------------------------------")
        println()
    end
end

"""
### Get the memory that is used for an object in human-readable units
```julia
function usedmemory(object, verbose::Bool=true)
return value, unit
```
"""
function usedmemory(object, verbose::Bool=true)
    obj_value = Base.summarysize(object)

    return usedmemory(obj_value, verbose)
end

function usedmemory(obj_value::Real, verbose::Bool=true)

    value_buffer = obj_value
    value_unit = "Bytes"
    if obj_value > 1000.
        value_buffer = obj_value / 1024.
        value_unit = "KB"
        if value_buffer > 1000.
            value_buffer = value_buffer / 1024.
            value_unit = "MB"
            if value_buffer > 1000.
                value_buffer = value_buffer / 1024.
                value_unit = "GB"
                if value_buffer > 1000.
                    value_buffer = value_buffer / 1024.
                    value_unit = "TB"
                end
            end
        end
    end

    if verbose == true
        println("Memory used: ", round(value_buffer, digits=3), " ", value_unit)
    end

    return value_buffer, value_unit
end

"""
#### Print overview of the used storage per file type for a given timestep

```julia
function storageoverview(dataobject::InfoType; verbose::Bool=verbose_mode)
return dictionary in bytes
```
"""
function storageoverview(dataobject::InfoType; verbose::Bool=verbose_mode)
    # todo simplyfy to single function calls

    dictoutput = Dict()

    output = dataobject.output
    printstyled("Overview of the used disc space for output: [$output]\n", bold=true, color=:normal)
    printstyled("------------------------------------------------------\n", bold=true, color=:normal)
    #path = dataobject.path
    #fnames = createpath(output, path)
    #println(fnames)
    fnames = dataobject.fnames

    all_files = readdir(fnames.output)
    folder = filesize.(  fnames.output .* "/" .* all_files)
    folder_size = sum( folder )
    folder_mean = mean( folder )
    folder_value, folder_unit =  usedmemory(folder_size, false)
    folder_meanvalue, folder_meanunit =  usedmemory(folder_mean, false)
    println( "Folder:         ", round(folder_value,digits=2),   " ", folder_unit,  " \t<", round(folder_meanvalue, digits=2),  " ", folder_meanunit,">/file" )

    amr_files = all_files[ occursin.( "amr", all_files) ]
    amr = filesize.( fnames.output .* "/" .* amr_files )
    amr_size = sum( amr )
    amr_mean = mean( amr )
    amr_value, amr_unit =  usedmemory(amr_size, false)
    amr_meanvalue, amr_meanunit =  usedmemory(amr_mean, false)
    println( "AMR-Files:      ", round(amr_value, digits=2),     " ", amr_unit,     " \t<", round(amr_meanvalue, digits=2),     " ", amr_meanunit,">/file" )

    if dataobject.hydro
        hydro_files = all_files[ occursin.( "hydro", all_files) ]
        hydro = filesize.( fnames.output .* "/" .* hydro_files )
        hydro_size = sum( hydro )
        hydro_mean = mean( hydro )
        hydro_value, hydro_unit =  usedmemory(hydro_size, false)
        hydro_meanvalue, hydro_meanunit =  usedmemory(hydro_mean, false)
        println( "Hydro-Files:    ", round(hydro_value, digits=2),   " ", hydro_unit,   " \t<", round(hydro_meanvalue, digits=2),   " ", hydro_meanunit,">/file" )
    else
        hydro_size = 0.
    end

    if dataobject.gravity
        gravity_files = all_files[ occursin.( "grav", all_files) ]
        gravity = filesize.(fnames.output .* "/" .* gravity_files )
        gravity_size = sum( gravity )
        gravity_mean = mean( gravity )
        gravity_value, gravity_unit =  usedmemory(gravity_size, false)
        gravity_meanvalue, gravity_meanunit =  usedmemory(gravity_mean, false)
        println( "Gravity-Files:  ", round(gravity_value, digits=2), " ", gravity_unit, " \t<", round(gravity_meanvalue, digits=2), " ", gravity_meanunit,">/file" )
    else
        gravity_size = 0.
    end

    if dataobject.particles
        particle_files = all_files[ occursin.( "part", all_files) ]
        particle = filesize.( fnames.output .* "/" .* particle_files )
        particle_size = sum( particle )
        particle_mean = mean( particle )
        particle_value, particle_unit =  usedmemory(particle_size, false)
        particle_meanvalue, particle_meanunit =  usedmemory(particle_mean, false)
        println( "Particle-Files: ", round(particle_value, digits=2)," ", particle_unit," \t<", round(particle_meanvalue, digits=2)," ", particle_meanunit,">/file" )
    else
        particle_size = 0.
    end

    if dataobject.clumps
        clump_files = all_files[ occursin.( "clump", all_files) ]
        clump = filesize.( fnames.output .* "/" .* clump_files )
        clump_size = sum( clump )
        clump_mean = mean( clump )
        clump_value, clump_unit =  usedmemory(clump_size, false)
        clump_meanvalue, clump_meanunit =  usedmemory(clump_mean, false)
        println( "Clump-Files:    ", round(clump_value, digits=2),   " ", clump_unit,   " \t<", round(clump_meanvalue, digits=2),   " ", clump_meanunit,">/file" )
    else
        clump_size = 0.
    end

    if dataobject.rt
        rt_files = all_files[ occursin.( "rt", all_files) ]
        rt = filesize.( fnames.output .* "/" .* rt_files )
        rt_size = sum( rt )
        rt_mean = mean( rt )
        rt_value, rt_unit =  usedmemory(rt_size, false)
        rt_meanvalue, rt_meanunit =  usedmemory(rt_mean, false)
        println( "RT-Files:       ", round(rt_value, digits=2),   " ", rt_unit,   " \t<", round(rt_meanvalue, digits=2),   " ", rt_meanunit,">/file" )
    else
        rt_size = 0.
    end

    # todo: check for sink files
    if dataobject.sinks
        sink_files = all_files[ occursin.( "sink", all_files) ]
        sink = filesize.( fnames.output .* "/" .* sink_files )
        sink_size = sum( sink )
        sink_mean = mean( sink )
        sink_value, sink_unit =  usedmemory(sink_size, false)
        sink_meanvalue, sink_meanunit =  usedmemory(sink_mean, false)
        println( "Sink-Files:    ", round(sink_value, digits=2),   " ", sink_unit,   " \t<", round(sink_meanvalue, digits=2),   " ", sink_meanunit,">/file" )

    else
        sink_size = 0.
    end

    println()
    println()
    println("mtime: ", dataobject.mtime)
    println("ctime: ", dataobject.ctime)

    # prepare output
    dictoutput[:folder] = folder_size
    dictoutput[:amr] = amr_size
    dictoutput[:hydro] = hydro_size
    dictoutput[:gravity] = gravity_size
    dictoutput[:particle] = particle_size
    dictoutput[:clump] = clump_size
    dictoutput[:rt] = rt_size
    dictoutput[:sink] = sink_size

    return dictoutput

end

"""
### Get the number of cells and/or the CPUs per level
```julia
function overview_amr(dataobject::HydroDataType)
return a JuliaDB table
```
"""
function amroverview(dataobject::HydroDataType)

    checkforAMR(dataobject)

    # check if cpu column exists
    fn = propertynames(dataobject.data.columns)
    cpu_col = false
    Ncols = 2
    if  in(Symbol("cpu"), fn)
        cpu_col = true
        Ncols = 3
    end
    cells = zeros(Int, dataobject.lmax - dataobject.lmin + 1, Ncols)
    cellsize = zeros(Float64, dataobject.lmax - dataobject.lmin + 1,1)

    println("Counting...")
    @showprogress 1 "" for ilevel=dataobject.lmin:dataobject.lmax
        if cpu_col
         cpus_ilevel = length( unique( select( filter(p->p.level==ilevel, select(dataobject.data, (:level, :cpu) ) ), :cpu) ) )
         cells[Int(ilevel-dataobject.lmin+1),3] = cpus_ilevel
        end

        cells[Int(ilevel-dataobject.lmin+1),1] = ilevel
        cellsize[Int(ilevel-dataobject.lmin+1)] = dataobject.boxlen / 2^ilevel
    end

    cells_per_level = fit!(CountMap(Int), select(dataobject.data, (:level)) )
    #Nlevels = length(cells_per_level.value.keys)
    #for ilevel=1:(dataobject.lmax-dataobject.lmin)

        #if ilevel <= Nlevels
        for (ilevel,j) in enumerate(cells_per_level.value.keys)

                cells[j-dataobject.lmin+1,2] = cells_per_level.value.vals[ilevel]
            #else
            #    cells[ilevel,2] = 0.

        end
    #end

    if cpu_col
        amr_hydro_table = table(cells[:,1], cells[:,2], cellsize[:], cells[:,3], names=[:level, :cells, :cellsize, :cpus])
    else
        amr_hydro_table = table(cells[:,1], cells[:,2], cellsize[:], names=[:level, :cells, :cellsize])
    end

    return amr_hydro_table

end

"""
### Get the number of particles and/or the CPUs per level
```julia
function overview_amr(dataobject::PartDataType)
return a JuliaDB table
```
"""
function amroverview(dataobject::PartDataType)

    checkforAMR(dataobject)

    # check if cpu column exists
    fn = propertynames(dataobject.data.columns)
    cpu_col = false
    Ncols = 2
    if  in(Symbol("cpu"), fn)
        cpu_col = true
        Ncols = 3
    end

    parts = zeros(Int, dataobject.lmax - dataobject.lmin + 1, Ncols)

    part_tot = 0
    part_masstot = 0
    println("Counting...")
    @showprogress 1 "" for ilevel=dataobject.lmin:dataobject.lmax
        if cpu_col
            cpus_ilevel = length( unique( select( filter(p->p.level==ilevel, select(dataobject.data, (:level, :cpu) ) ), :cpu) ) )
            parts[Int(ilevel-dataobject.lmin+1),3] = cpus_ilevel
        end

        parts[Int(ilevel-dataobject.lmin+1),1] = ilevel

    end

    part_per_level = fit!(CountMap(Int32), select(dataobject.data, (:level)) )
    Nlevels = length(part_per_level.value.keys)
    for ilevel=1:Nlevels
        if ilevel <= Nlevels
            parts[ilevel,2] = part_per_level.value.vals[ilevel]
        else
            parts[ilevel,2] = 0.
        end
    end

    if cpu_col
        amr_part_table = table(parts[:,1], parts[:,2], parts[:,3], names=[:level, :particles, :cpus])
    else
        amr_part_table = table(parts[:,1], parts[:,2], names=[:level, :particles])
    end

    return amr_part_table

end

function checkforAMR(dataobject::DataSetType)
    if dataobject.lmax == dataobject.lmin
        error("[Mera]: Works only with AMR data!")
    end
end

"""
### Get the mass and min/max value of each variable in the database per level

```julia
function overview_data(dataobject::HydroDataType)
return a JuliaDB table
```
"""
function dataoverview(dataobject::HydroDataType)

    nvarh = dataobject.info.nvarh
    lmin = dataobject.lmin
    lmax = dataobject.lmax
    isamr = checkuniformgrid(dataobject, lmax)

    cells_tot = 0
    cells_masstot = 0
    density_var = :rho
    skip_vars = [:cpu, :level, :cx, :cy, :cz]

    if dataobject.info.descriptor.usehydro == true
        if haskey(dataobject.used_descriptors, 1)
            density_var = dataobject.used_descriptors[1]
        end
    end
    names_constr = [Symbol("level")]
    fn = propertynames(dataobject.data.columns)
    for i in fn
        if !in(i, skip_vars)
            if i == density_var
                append!(names_constr, [Symbol("mass")] )
            end
            append!(names_constr, [Symbol("$(i)_min")] )
            append!(names_constr, [Symbol("$(i)_max")] )

        end
    end

    cells = Array{Any,2}(undef, (dataobject.lmax - dataobject.lmin + 1,length(names_constr) ) )
    println("Calculating...")
    @showprogress 1 "" for ilevel=lmin:lmax
        cell_iterator = 1

        if isamr
            filtered_level = filter(p->p.level==ilevel, dataobject.data )
        else # if uniform grid
            filtered_level = dataobject.data
        end

        cells[Int(ilevel-lmin+1),cell_iterator] = ilevel
        cell_iterator= cell_iterator + 1

        for ifn in fn
            if !in(ifn, skip_vars)
                if ifn == density_var
                    cells_msum = sum(select(filtered_level , density_var)) * (dataobject.boxlen / 2^ilevel)^3
                    #todo: introduce humanize for mass
                    #cells_masstot = cells_masstot + cells_msum
                    cells[Int(ilevel-lmin+1),cell_iterator] = cells_msum
                    cell_iterator= cell_iterator + 1
                    if length(select(filtered_level, density_var)) != 0
                        rho_minmax = reduce((min, max), filtered_level, select=density_var)
                        rhomin= rho_minmax.min
                        rhomax= rho_minmax.max
                    else
                        rhomin= 0.
                        rhomax= 0.
                    end
                    cells[Int(ilevel-lmin+1),cell_iterator] = rhomin
                    cell_iterator= cell_iterator + 1
                    cells[Int(ilevel-lmin+1),cell_iterator] = rhomax
                    cell_iterator= cell_iterator + 1

                else
                    if length(select(filtered_level, ifn)) != 0
                        value_minmax = reduce((min, max), filtered_level, select=ifn)
                        valuemin = value_minmax.min
                        valuemax = value_minmax.max
                    else
                        valuemin = 0.
                        valuemax = 0.
                    end
                    cells[Int(ilevel-lmin+1),cell_iterator] = valuemin
                    cell_iterator= cell_iterator + 1
                    cells[Int(ilevel-lmin+1),cell_iterator] = valuemax
                    cell_iterator= cell_iterator + 1
                end
            end
        end

    end

    hydro_overview_table = table( [cells[:, i ] for i = 1:length(names_constr)]..., names=[names_constr...] )
    return hydro_overview_table
end

"""
### Get the extrema of each variable in the database

```julia
function overview_data(dataobject::ClumpDataType)
return a JuliaDB table
```
"""
function dataoverview(dataobject::ClumpDataType)
    fn = propertynames(dataobject.data.columns)
    s = Series(Extrema())
    Ncolumns = length(fn)
    values = Array{Any}(undef, Ncolumns+1, 2)
    values[1,1] = "min"
    values[1,2] = "max"
    for i = 1:Ncolumns
        a = reduce(s, dataobject.data; select = fn[i])
        values[i+1, 1] = a.stats[1].min
        values[i+1, 2] = a.stats[1].max
    end
    column_names = [Symbol("extrema")]
    append!(column_names, fn )

    clump_overview_table = table( [values[k, : ] for k = 1:(Ncolumns+1) ]...,names=collect(column_names) )
    return clump_overview_table
end

"""
### Get the min/max value of each variable in the database per level

```julia
function overviewdata(dataobject::PartDataType)
return a JuliaDB table
```
"""
function dataoverview(dataobject::PartDataType)

    lmin = dataobject.lmin
    lmax = dataobject.lmax
    isamr = checkuniformgrid(dataobject, lmax)

    parts_tot = 0
    parts_masstot = 0
    skip_vars = [:cpu, :level]
    fn = propertynames(dataobject.data.columns)
    #fn = convert(Array{String,1}, fn) #todo delte
    names_constr = [Symbol("level")]
    for i in fn
        if !in(i, skip_vars)
            append!(names_constr, [Symbol("$(i)_min")] )
            append!(names_constr, [Symbol("$(i)_max")] )

        end
    end

    parts = Array{Any,2}(undef, (dataobject.lmax - dataobject.lmin + 1,length(names_constr) ) )
    #@showprogress 1 "Searching..."
    for ilevel=lmin:lmax
        part_iterator = 1

        if isamr
            filtered_level = filter(p->p.level==ilevel, dataobject.data )
        else # if uniform grid
            filtered_level = dataobject.data
        end

        parts[Int(ilevel-lmin+1),part_iterator] = ilevel
        part_iterator= part_iterator + 1

        for ifn in fn
            if !in(ifn, skip_vars)
                if length(select(filtered_level, ifn)) != 0
                    value_minmax = reduce((min, max), filtered_level, select=ifn)
                    valuemin = value_minmax.min
                    valuemax = value_minmax.max
                else
                    valuemin = 0.
                    valuemax = 0.
                end
            parts[Int(ilevel-lmin+1),part_iterator] = valuemin
            part_iterator= part_iterator + 1
            parts[Int(ilevel-lmin+1),part_iterator] = valuemax
            part_iterator= part_iterator + 1
            end
        end

    end

    particle_overview_table = table( [parts[:, i ] for i = 1:length(names_constr)]..., names=[names_constr...] )
    return particle_overview_table
end

"""
#### Get the existing simulation snapshots in a given folder
- returns field `outputs` with Array{Int,1} containing the output-numbers of the existing simulations
- returns field `missing` with Array{Int,1} containing the output-numbers of empty simulation folders

```julia
checkoutputs(path::String="./")
return CheckOutputNumberType
```

#### Examples
```julia
# Example 1:
# look in current folder
julia> N = checkoutputs();
julia> N.outputs
julia> N.missing

# Example 2:
# look in given path
# without any keyword
julia>N = checkoutputs("simulation001");
```

"""
function checkoutputs(path::String="./")

    if path == "" || path == " " path="./" end
    folder = readdir(path)

    # filter "output_" - names
    ftrue = occursin.("output_", folder)
    folder = folder[ftrue]

    # create full path to supposed folders
    output_path = joinpath.(path,folder)

    # filter real folders
    folder = folder[isdir.(output_path)]

    # get existing output numbers
    missing_outputs = Int[]
    existing_outputs = Int[]
    Nfolders = length(folder)
    if Nfolders > 0
        Noutputs_string = [folder[x][8:end] for x in 1:Nfolders]
        Noutputs = parse.(Int, Noutputs_string)

        # find missing output numbers
        maxoutput = maximum(Noutputs)
        #expected_outputs = 1:maxoutput
        for Nout in Noutputs
            fnames = createpath(Nout, path)
            isinfofile = isfile(fnames.info)
            if isinfofile
                append!(existing_outputs, Nout)
            else
                append!(missing_outputs, Nout)
            end
        end
    end

    return CheckOutputNumberType(existing_outputs, missing_outputs)
end

"""
#### Get physical time in selected units
returns Float

```julia
gettime(output::Real; path::String="./", unit::Symbol=:standard)
gettime(dataobject::DataSetType; unit::Symbol=:standard)

return time
```

#### Arguments Function 1
##### Required:
- **`output`:** give the output-number of the simulation

##### Predefined/Optional Keywords:
- **`path`:** the path to the output folder relative to the current folder or absolute path
- **`unit`:** return the variable in given unit

#### Arguments Function 2
##### Required:
- **`dataobject`:** needs to be of type: "DataSetType"

##### Predefined/Optional Keywords:
- **`unit`:** return the variable in given unit

"""
function gettime(output::Real, path::String, unit::Symbol;)
    return gettime(output, path=path, unit=unit)
end

function gettime(output::Real, path::String; unit::Symbol=:standard)
    return gettime(output, path=path, unit=unit)
end

function gettime(output::Real; path::String="./", unit::Symbol=:standard)
    info = getinfo(output, path, verbose=false)
    return info.time * getunit(info, unit)
end

function gettime(dataobject::DataSetType, unit::Symbol;)
    return gettime(dataobject, unit=unit)
end

function gettime(dataobject::DataSetType; unit::Symbol=:standard)
    return dataobject.info.time * getunit(dataobject.info, unit)
end


text:
module TBIO

    using Formatting
    export readTB!, writeTB

    function readTB!(
        filename   :: AbstractString,
        txt_nchars :: Integer,
        arrs       :: AbstractArray;
        endianess  :: Symbol = :LITTLE,   # Small: 0x04030201,  Big: 0x01020304
        trim       :: Bool   = true,
    )

        if endianess != :LITTLE && endianess != :BIG
            throw(ErrorException("Unknown symbol: " * string(endianess)))
        end

        if txt_nchars < 0
            throw(ErrorException("txt_nchars cannot be negative."))
        end

        msg = nothing
        if isfile(filename)

            get_filesize = filesize(filename)
            expect_filesize = sum(length.(arrs)) * 8 + txt_nchars + 8
            if get_filesize == expect_filesize

                open(filename, "r") do io
                    msg = String(read(io, txt_nchars))

                    trim && (msg = strip(msg))

                    for i = 1:length(arrs)
                        read!(io, arrs[i])
                    end

                    if     endianess == :LITTLE && Base.ENDIAN_BOM == 0x01020304
                        for i = 1:length(arrs)
                            arrs[i][:] = ltoh.(arrs[i])
                        end
                    elseif endianess == :BIG && Base.ENDIAN_BOM == 0x04030201
                        for i = 1:length(arrs)
                            arrs[i][:] = ntoh.(arrs[i])
                        end
                    end

                    received_cs   = reinterpret(UInt64, read(io, 8))[1]
                    calculated_cs = calChecksum(arrs)

                    if received_cs != calculated_cs
                        msg = nothing
                        println(format("[readTB!] Checksum does not match. Received: {:X}, calculated: {:X}", received_cs, calculated_cs))
                    end

                end
            else
                println(format("[readTB!] Expecting filesize: {} bytes, but got {} bytes", expect_filesize, get_filesize))
            end
        else
            println(format("[readTB!] File {} does not exist.", filename))
        end

        return msg
    end

    function writeTB(
        filename   :: AbstractString,
        msg        :: AbstractString,
        txt_nchars :: Integer,
        arrs       :: AbstractArray{T};
        endianess  :: Symbol = :LITTLE,
    ) where T <: AbstractArray{Float64}

        if endianess != :LITTLE && endianess != :BIG
            throw(ErrorException("Unknown symbol: " * string(endianess)))
        end

        if length(msg) > txt_nchars
            throw(ErrorException("Message length exceeds txt_nchars."))
        end

        if txt_nchars < 0
            throw(ErrorException("txt_nchars cannot be negative."))
        end

        open(filename, "w") do io
            write(io, msg)

            for i = 1:(txt_nchars - length(msg))
                write(io, " ")
            end

            if     endianess == :LITTLE && Base.ENDIAN_BOM == 0x01020304
                for i = 1:length(arrs)
                    write(io, htol.(arrs[i]))
                end
                write(io, htol(calChecksum(arrs)))
            elseif endianess == :BIG && Base.ENDIAN_BOM == 0x04030201
                for i = 1:length(arrs)
                    write(io, hton.(arrs[i]))
                end
                write(io, hton(calChecksum(arrs)))
            else
                for i = 1:length(arrs)
                    write(io, arrs[i])
                end
                write(io, calChecksum(arrs))
            end

        end
    end

    # ror, rol :: bitwise circular shift to right / left
    # Source:
    # https://discourse.julialang.org/t/efficient-bit-rotate-functions-ror-rol-what-is-the-official-way-for-julia-v1/19062
    ror(x::UInt64, k::Int) = (x >>> (0x3f & k)) | (x << (0x3f & -k))
    rol(x::UInt64, k::Int) = (x >>> (0x3f & -k)) | (x << (0x3f & k))

    function calChecksum(
        arrs :: AbstractArray,
    )

        k = 0
        s = UInt64(0)
        for i = 1:length(arrs)
            arr = reinterpret(UInt64, arrs[i])

            for j = 1:length(arr)

                s = xor(s, rol(arr[j], k))
                #s = xor(s, arr[j])
                k = mod(k+1, 64)
            end

        end

#        println(format("{:X}", s))
        return s
    end
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
function test_indexing(AT)
    # TODO: more fine-grained allowscalar within test_indexing
    GPUArrays.@allowscalar @testset "indexing" begin
        for T in (Float32, Int32#=, SVector{3, Float32}=#)
            @testset "Indexing with $T" begin
                x = rand(T, 32)
                src = AT(x)
                for (i, xi) in enumerate(x)
                    @test src[i] == xi
                end
                @test Array(src[1:3]) == x[1:3]
                @test Array(src[3:end]) == x[3:end]
            end
            @testset "multi dim, sliced setindex" begin
                x = fill(AT{T}, T(0), (10, 10, 10, 10))
                y = rand(AT{T}, 5, 5, 10, 10)
                x[2:6, 2:6, :, :] = y
                x[2:6, 2:6, :, :] == y
           end

        end

        for T in (Float32, Int32)
            @testset "Indexing with $T" begin
                x = fill(zero(T), 7)
                src = AT(x)
                for i = 1:7
                    src[i] = i
                end
                @test Array(src) == T[1:7;]
                src[1:3] = T[77, 22, 11]
                @test Array(src[1:3]) == T[77, 22, 11]
                src[1] = T(0)
                src[2:end] = T(77)
                @test Array(src) == T[0, 77, 77, 77, 77, 77, 77]
            end
        end

        for T in (Float32, Int32)
            @testset "issue #42 with $T" begin
                Ac = rand(Float32, 2, 2)
                A = AT(Ac)
                @test A[1] == Ac[1]
                @test A[end] == Ac[end]
                @test A[1, 1] == Ac[1, 1]
            end
        end
        for T in (Float32, Int32)
            @testset "Colon() $T" begin
                Ac = rand(T, 10)
                A = AT(Ac)
                A[:] = T(1)
                @test all(x-> x == 1, A)
                A[:] = AT(Ac)
                @test Array(A) == Ac
            end
        end
    end
end


text:
export plot!, plot, scatter!, scatter,
       fill_between!,
       hist!, hist, bar!, bar,
       boxplot, heatmap, heatmap_ticks

"""
    plot_data(x)
    plot_data(x, ys...)

I/ Preprocess data for a 2D plot.
"""
function plot_data(
            x::AVM{<:CanMiss{<:Real}}
        )::NamedTuple

    isempty(x) && throw(
        ArgumentError("Cannot display empty vectors.")
    )
    hasmissing = Missing <: eltype(x) ||
                   any(isinf, x)      ||
                   any(isnan, x)
    return (
        data       = zip(1:size(x, 1), eachcol(x)...),
        hasmissing = hasmissing,
        nobj       = size(x, 2)
    )
end
# NOTE: these typechecks within the function body is to avoid clashes/ambiguity with plot_data(x)
function plot_data(
            x,
            ys...
        )::NamedTuple

    isempty(x) && throw(
        ArgumentError("Cannot display empty vectors.")
    )
    x isa AV{<:CanMiss{<:Real}} || throw(
        ArgumentError("x has un-handled type $(typeof(x))")
    )

    nobj       = 0
    hasmissing = Missing <: eltype(x) ||
                   any(isinf, x)      ||
                   any(isnan, x)
    for y ∈ ys
        y isa AVM{<:CanMiss{<:Real}} || throw(
            ArgumentError("y has un-handled type $(typeof(y))")
        )
        size(y, 1) == length(x) || throw(
            DimensionMismatch("y data must match x")
        )

        nobj      += size(y, 2)
        hasmissing |= Missing <: eltype(y) ||
                        any(isinf, y)      ||
                        any(isnan, y)
    end
    return (
        data       = zip(x, (view(y, :, j) for y ∈ ys for j ∈ axes(y, 2))...),
        hasmissing = hasmissing,
        nobj       = nobj
    )
end

"""
    fill_data(x, y1, y2)

I/ Preprocess data for a Fill2D.
"""
function fill_data(
            x::AVR,
            y1::Union{Real, AVR},
            y2::Union{Real,AVR}
        )::NamedTuple

    isempty(x) && throw(
        ArgumentError("Cannot display empty vectors.")
    )

    y1 isa AV || (y1 = fill(y1, length(x)))
    y2 isa AV || (y2 = fill(y2, length(x)))
    length(x) == length(y1) == length(y2) ||throw(
        DimensionMismatch("vectors x,y1,y2 must have matching lengths"))
    return (
        data = zip(x, y1, y2),
    )
end

"""
    hist_data(x)

I/ Preprocess data for a Hist2D.
Missing values are allowed but are not counted as observations.
"""
function hist_data(
            x::AV{<:CanMiss{<:Real}}
        )::NamedTuple

    isempty(x) && throw(
        ArgumentError("Cannot have a histogram with an empty vector.")
    )

    sx = skipmissing(x)
    isempty(sx) && throw(
        ArgumentError("Cannot have a histogram with only missing values")
    )
    return (
        data       = zip(x),
        hasmissing = (Missing <: eltype(x)),
        nobs       = sum(e -> 1, sx),
        range      = Float64.(extrema(sx))
    )
end

#################################################################################
#################################################################################

"""
    plot!(...)

E/ Add a plot. Keyword arguments can be passed to specify the linestyle(s),
label(s) and markerstyle(s).

## Example

    x = range(-2, 2, length=100)
    y = @. exp(-abs(x)+sin(x))
    plot(x, y,
         color  = "blue",
         lstyle = "--",
         marker = "o",
         lwidth = 0.05,
         label  = "First plot"
    )

"""
function plot!(
            x, ys...;
            #
            axes::Option{Axes2D} = nothing,
            overwrite::Bool      = false,
            o...
        )::DrawingHandle{<:Scatter2D}

    # are current axes empty?
    # -> if so don't do anything as the user may have pre-specified things
    #    like xlim etc and want the plot to appear with those
    # -> if it's not empty, reset the axes (will destroy xlim settings etc
    #    as it should)
    axes = check_axes(axes)
    if overwrite && !all(isempty, (axes.drawings, axes.objects))
        reset!(axes; parent=axes.parent)
    end

    # Form the scatter object
    pd      = plot_data(x, ys...)
    scatter = Scatter2D(pd.data, pd.hasmissing, pd.nobj)

    # Add properties & add scatter to axes
    set_properties!(scatter; o...)
    push!(axes.drawings, scatter)

    return DrawingHandle(scatter)
end

function plot!(
            f::Function,
            from::Real,
            to::Real;
            #
            length::Int = 100,
            o...
        )

    x = range(from, stop=to, length=length)
    plot!(x, f.(x); o...)
end

plot(a...; o...)     = plot!(a...;    overwrite=true, o...)
scatter!(a...; o...) = plot!(a...;    ls="none", marker=".", o...)
scatter(a...; o...)  = scatter!(a...; overwrite=true, o...)

"""
    fill_between!(...)

Add a fill plot between two lines. The arguments must not have missings but
`y1` and/or `y2` can be specified as single numbers (= horizontal line).
"""
function fill_between!(
            x, y1, y2;
            #
            axes::Option{Axes2D} = nothing,
            overwrite::Bool      = false,
            o...
        )::DrawingHandle{<:Fill2D}

    axes = check_axes(axes)
    if overwrite && !all(isempty, (axes.drawings, axes.objects))
        reset!(axes; parent=axes.parent)
    end

    # form the fill2d object
    fd   = fill_data(x, y1, y2)
    fill = Fill2D(fd)
    set_properties!(fill; o...)

    push!(axes.drawings, fill)
    return DrawingHandle(fill)
end

fill_between(a...; o...) = fill_between!(a...; overwrite=true, o...)

"""
    hist!(...)

Add a histogram of `x` on the current axes.
"""
function hist!(
            x;
            #
            axes::Option{Axes2D} = nothing,
            overwrite::Bool      = false,
            o...
        )::DrawingHandle{<:Hist2D}

    axes = check_axes(axes)
    if overwrite
        all(isempty, (axes.drawings, axes.objects)) || reset!(axes; parent=axes.parent)
    end

    hd   = hist_data(x)
    hist = Hist2D(hd.data, hd.hasmissing, hd.nobs, hd.range)
    set_properties!(hist; o...)

    push!(axes.drawings, hist)
    return DrawingHandle(hist)
end

hist(a...; o...)  = hist!(a...; overwrite=true, o...)

####
#### bar!, bar
####

"""
    bar!(...)

Add a bar plot.
"""
function bar!(
            x, ys...;
            #
            axes::Option{Axes2D} = nothing,
            overwrite::Bool      = false,
            o...
        )::DrawingHandle{<:Bar2D}

    axes = check_axes(axes)
    if overwrite
        all(isempty, (axes.drawings, axes.objects)) || reset!(axes; parent=axes.parent)
    end

    bd  = plot_data(x, ys...)
    bar = Bar2D(bd.data, bd.hasmissing, bd.nobj)
    set_properties!(bar; o...)

    push!(axes.drawings, bar)
    return DrawingHandle(bar)
end

bar(a...; o...) =  bar!(a...; overwrite=true, o...)

#################################################################################
#################################################################################
#
# Extra drawings that are not meant to be overlaid with anything
# else such as boxplot, heatmap, polarplot, pieplot
#

"""
    boxplot(...)

Erase previous drawings and add a boxplot. Missing values are allowed but not
Infinities or Nans.
"""
function boxplot(
            ys...;
            #
            axes::Option{Axes2D} = nothing,
            o...
        )::DrawingHandle{Boxplot}

    isempty(first(ys)) && throw(
        ArgumentError("Cannot display empty vectors.")
    )

    # always on fresh axes
    axes = check_axes(axes)
    reset!(axes; parent=axes.parent)

    # setting an empty struct first so that we can exploit the options
    # the actual data will be provided after analysis
    nobj = sum(size(y, 2) for y ∈ ys)
    bp = Boxplot(Matrix{Float64}(undef,0,0), nobj)
    set_properties!(bp; o...)

    # analyse data
    # 6 -> 1:wlow, 2:q25, 3:q50, 4:q75, 5:whigh, 6:mean
    stats    = Matrix{Float64}(undef, nobj, 6)
    outliers = Vector{Vector{Float64}}(undef, nobj)

    boxcounter = 1
    axmin      = Inf
    axmax      = -Inf

    for y ∈ ys
        for k ∈ Base.axes(y, 2)
            yk = collect(skipmissing(view(y, :, k)))
            if any(isnan, yk) || any(isinf, yk)
                throw(
                    ArgumentError("Inf or NaN values not allowed in boxplot.")
                )
            end

            q00, q25, q50, q75, q100 = quantile(yk, [.0, .25, .5, .75, 1.0])
            iqr = q75 - q25
            mu  = mean(yk)

            wrlength = bp.boxstyles[k].wrlength
            wlow  = q25 - wrlength * iqr
            whigh = q75 + wrlength * iqr
            if isinf(wrlength)
                wlow, whigh = q00, q100 # min/max values
            end

            stats[k, :] = [wlow, q25, q50, q75, whigh, mu]

            # outliers
            outliers[k] = filter(e->(e<wlow || whigh<e), yk)

            # Keep track of extremes to set axis limits afterwards
            axmin_ = min(q00, wlow)
            axmin_ -= 0.5abs(axmin_)
            axmin_ < axmin && (axmin = axmin_)
            axmax_ = max(q100, whigh)
            axmax_ += 0.5abs(axmax_)
            axmax_ > axmax && (axmax = axmax_)

            boxcounter += 1
        end
    end
    #
    bp.stats = stats
    push!(axes.drawings, bp)

    if bp.horiz # horizontal boxplot
        ylim(0, nobj+1)
        yticks(1:nobj)
        xlim(axmin, axmax)
        for k ∈ 1:nobj
            bp.boxstyles[k].oshow || continue
            nok = length(outliers[k])
            s = bp.boxstyles[k].omstyle # style of outliers
            if nok > 0
                scatter!(outliers[k], fill(k, nok); marker=s.marker, msize=s.msize, mcol=s.color)
            end
        end
    else # vertical boxplot
        xlim(0, nobj+1)
        xticks(1:nobj)
        ylim(axmin, axmax)
        for k ∈ 1:nobj
            bp.boxstyles[k].oshow || continue
            nok = length(outliers[k])
            s = bp.boxstyles[k].omstyle # style of outliers
            if nok > 0
                scatter!(fill(k, nok), outliers[k]; marker=s.marker, msize=s.msize, mcol=s.color)
            end
        end
    end
    return DrawingHandle(bp)
end

"""
    heatmap_ticks(ax)

Returns ticks position centered appropriately for a heatmap.

### Example

    heatmap(randn(2, 2))
    xticks(heatmap_ticks("x"), ["var1", "var2"]; angle=45)

"""
function heatmap_ticks(
            ax::String,
            axes::Axes2D
        )::Vector{F64}

    # it must necessarily be the first object since heatmap resets the axes
    (isempty(axes.drawings) || !(axes.drawings[1] isa Heatmap)) && throw(
        ArgumentError("No heatmap found.")
    )

    nrows, ncols = size(axes.drawings[1].data)
    ax_lc        = lowercase(ax)

    ax_lc ∈ ["x", "x2"] && return (collect(0:(ncols-1))    .+ 0.5) ./ ncols
    ax_lc ∈ ["y", "y2"] && return (collect((nrows-1):-1:0) .+ 0.5) ./ nrows

    throw(
        ArgumentError("""
            Unrecognised ax descriptor expected one of [x, x2, y, y2].
            """
        )
    )
end
heatmap_ticks(ax::String) = heatmap_ticks(ax, gca())

"""
    heatmap(X; opts...)

Creates a heatmap corresponding to matrix `X`. Doing this with matrices larger
than 100x100 with the GLE backend can be slow. As an indication, 500x500 takes
about 15-20 seconds on a standard laptop. Time scales with the number of elements
with, as a rough order of magnitude, around 10k elements per second.

You can provide ticks using `xticks(heatmap_ticks("x"), ["tick1", ...]; o...)`.
See also `heatmap_ticks`.

Note: handles missing values but not NaN or Inf.
"""
function heatmap(
            data::Matrix{<:CanMiss{Real}};
            #
            axes::Option{Axes2D} = nothing,
            o...
        )::DrawingHandle{Heatmap}

    nrows, ncols = size(data)
    min(nrows, ncols) <= 1000 || throw(
        ArgumentError(
            """
            The matrix is too large to be displayed as a heatmap.
            """
        )
    )

    # always on fresh axes
    axes = check_axes(axes)
    reset!(axes; parent=axes.parent)

    zmin, zmax = extrema(skipmissing(data))

    h = Heatmap(data=Matrix{Int}(undef, 0, 0), zmin=zmin, zmax=zmax)
    h.transpose = (nrows < ncols)
    set_properties!(h; o...)
    push!(axes.drawings, h)

    ncolors = length(h.cmap)

    minv, maxv = extrema(skipmissing(data))
    incr = (maxv - minv) / ncolors

    data2 = Matrix{Int}(undef, nrows, ncols)
    if iszero(incr)
        data2 .= ceil(ncols/2)
    else
        @inbounds for i=1:nrows, j=1:ncols
            dij = data[i,j]
            if ismissing(dij)
                data2[i,j] = ncolors + 1 # offset doesn't matter, anything ≥ 1 is fine
            else
                data2[i,j] = ceil(Integer, (dij - minv) / incr)
            end
        end
    end
    # NOTE: there may be a few 0 (min values), will be treated as 1 see add_sub_heatmap!

    # store the color assignments
    h.data = data2

    # set axis
    xlim(0,1); ylim(0,1)

    # add labels if provided to the heatmap object
    xplaces = heatmap_ticks("x")
    yplaces = heatmap_ticks("y")

    # used for sparse labelling in case ncols > 20
    stepx    = (ncols<=50)*10 + (50<ncols<=200)*50 + (200<ncols)*100
    lastx    = ncols - mod(ncols, stepx) + 1
    maskx    = collect(1:stepx:lastx) .- 1
    maskx[1] = 1
    stepy    = (nrows<=50)*10 + (50<nrows<=200)*50 + (200<nrows)*100
    lasty    = nrows - mod(nrows, stepy) + 1
    masky    = collect(1:stepy:lasty) .- 1
    masky[1] = 1

    if ncols <= 20
        xticks(xplaces, ["$i" for i ∈ 1:ncols]; axes=axes, fontsize=12)
    else
        xticks(xplaces[maskx], ["$i" for i ∈ maskx]; axes=axes, fontsize=12)
    end
    if nrows <= 20
        yticks(yplaces, ["$i" for i ∈ 1:nrows]; axes=axes, fontsize=12)
    else
        yticks(yplaces[masky], ["$i" for i ∈ masky]; axes=axes, fontsize=12)
    end
    x2ticks("off")
    y2ticks("off")
    return DrawingHandle(h)
end


text:
"""
        color_graph(g::LightGraphs.AbstractGraphs, ::AcyclicColoring)

Returns a coloring vector following the acyclic coloring rules (1) the coloring
corresponds to a distance-1 coloring, and (2) vertices in every cycle of the
graph are assigned at least three distinct colors. This variant of coloring is
called acyclic since every subgraph induced by vertices assigned any two colors
is a collection of trees—and hence is acyclic.

Reference: Gebremedhin AH, Manne F, Pothen A. **New Acyclic and Star Coloring Algorithms with Application to Computing Hessians**
"""
function color_graph(g::LightGraphs.AbstractGraph, ::AcyclicColoring)
    color = zeros(Int, nv(g))
    two_colored_forest = DisjointSets{Int}(())

    first_visit_to_tree = fill((0,0), ne(g))
    first_neighbor = fill((0,0), ne(g))

    forbidden_colors = zeros(Int, nv(g))

    for v in vertices(g)
        for w in outneighbors(g, v)
            if color[w] != 0
                forbidden_colors[color[w]] = v
            end
        end

        for w in outneighbors(g, v)
            if color[w] != 0
                for x in outneighbors(g, w)
                    if color[x] != 0
                        if forbidden_colors[color[x]] != v
                            prevent_cycle!(first_visit_to_tree,forbidden_colors,v, w, x, g, two_colored_forest,color)
                        end
                    end
                end
            end
        end

        color[v] = min_index(forbidden_colors, v)

        for w in outneighbors(g, v)
            if color[w] != 0
                grow_star!(two_colored_forest,first_neighbor,v, w, g, color)
            end
        end

        for w in outneighbors(g, v)
            if color[w] != 0
                for x in outneighbors(g, w)
                    if color[x] != 0 && x != v
                        if color[x] == color[v]
                            merge_trees!(two_colored_forest,v,w,x,g)
                        end
                    end
                end
            end
        end
    end
    return color
end

"""
        prevent_cycle!(first_visit_to_tree::AbstractVector{<:Tuple{Integer,Integer}},
                        forbidden_colors::AbstractVector{<:Integer},
                        v::Integer,
                        w::Integer,
                        x::Integer,
                        g::LightGraphs.AbstractGraph,
                        two_colored_forest::DisjointSets{<:Integer},
                        color::AbstractVector{<:Integer})

Subroutine to avoid generation of 2-colored cycle due to coloring of vertex v,
which is adjacent to vertices w and x in graph g. Disjoint set is used to store
the induced 2-colored subgraphs/trees where the id of set is an integer
representing an edge of graph 'g'
"""
function prevent_cycle!(first_visit_to_tree::AbstractVector{<:Tuple{Integer,Integer}},
                        forbidden_colors::AbstractVector{<:Integer},
                        v::Integer,
                        w::Integer,
                        x::Integer,
                        g::LightGraphs.AbstractGraph,
                        two_colored_forest::DisjointSets{<:Integer},
                        color::AbstractVector{<:Integer})
    e = find(w, x, g, two_colored_forest)
    p, q = first_visit_to_tree[e]

    if p != v
        first_visit_to_tree[e] = (v,w)
    elseif q != w
        forbidden_colors[color[x]] = v
    end
end

"""
        grow_star!(two_colored_forest::DisjointSets{<:Integer},
                    first_neighbor::AbstractVector{<: Tuple{Integer,Integer}},
                    v::Integer,
                    w::Integer,
                    g::LightGraphs.AbstractGraph,
                    color::AbstractVector{<:Integer})

Grow a 2-colored star after assigning a new color to the
previously uncolored vertex v, by comparing it with the adjacent vertex w.
Disjoint set is used to store stars in sets, which are identified through key
edges present in g.
"""
function grow_star!(two_colored_forest::DisjointSets{<:Integer},
                    first_neighbor::AbstractVector{<: Tuple{Integer,Integer}},
                    v::Integer,
                    w::Integer,
                    g::LightGraphs.AbstractGraph,
                    color::AbstractVector{<:Integer})
    insert_new_tree!(two_colored_forest,v,w,g)
    p, q = first_neighbor[color[w]]

    if p != v
        first_neighbor[color[w]] = (v,w)
    else
        e1 = find(v,w,g,two_colored_forest)
        e2 = find(p,q,g,two_colored_forest)
        union!(two_colored_forest, e1, e2)
    end
end

"""
        merge_trees!(two_colored_forest::DisjointSets{<:Integer},
                      v::Integer,
                      w::Integer,
                      x::Integer,
                      g::LightGraphs.AbstractGraph)

Subroutine to merge trees present in the disjoint set which have a
common edge.
"""
function merge_trees!(two_colored_forest::DisjointSets{<:Integer},
                      v::Integer,
                      w::Integer,
                      x::Integer,
                      g::LightGraphs.AbstractGraph)
    e1 = find(v,w,g,two_colored_forest)
    e2 = find(w,x,g,two_colored_forest)
    if e1 != e2
        union!(two_colored_forest, e1, e2)
    end
end

"""
        insert_new_tree!(two_colored_forest::DisjointSets{<:Integer},
                          v::Integer,
                          w::Integer,
                          g::LightGraphs.AbstractGraph)

creates a new singleton set in the disjoint set 'two_colored_forest' consisting
of the edge connecting v and w in the graph g
"""
function insert_new_tree!(two_colored_forest::DisjointSets{<:Integer},
                          v::Integer,
                          w::Integer,
                          g::LightGraphs.AbstractGraph)
    edge_index = find_edge_index(v,w,g)
    push!(two_colored_forest,edge_index)
end

"""
        min_index(forbidden_colors::AbstractVector{<:Integer}, v::Integer)

Returns min{i > 0 such that forbidden_colors[i] != v}
"""
function min_index(forbidden_colors::AbstractVector{<:Integer}, v::Integer)
    return findfirst(!isequal(v), forbidden_colors)
end

"""
        find(w::Integer,
             x::Integer,
             g::LightGraphs.AbstractGraph,
             two_colored_forest::DisjointSets{<:Integer})

Returns the root of the disjoint set to which the edge connecting vertices w and x
in the graph g belongs to
"""
function find(w::Integer,
              x::Integer,
              g::LightGraphs.AbstractGraph,
              two_colored_forest::DisjointSets{<:Integer})
    edge_index = find_edge_index(w, x, g)
    return find_root(two_colored_forest, edge_index)
end

"""
        find_edge(g::LightGraphs.AbstractGraph, v::Integer, w::Integer)

Returns an integer equivalent to the index of the edge connecting the vertices
v and w in the graph g
"""
function find_edge_index(v::Integer, w::Integer, g::LightGraphs.AbstractGraph)
    pos = 1
    for i in edges(g)

        if (src(i) == v && dst(i) == w) || (src(i) == w && dst(i) == v)
            return pos
        end
        pos = pos + 1
    end
    throw(ArgumentError("$v and $w are not connected in the graph"))
end


text:
using LaTeXUnicode
using Test

@testset "LaTeXUnicode.jl" begin
    # Write your tests here.
end


text:
stress_scale = 1.0e5
strain_scale = 1

include("nnutil.jl")

# H0 = constant(H1/stress_scale)
testtype = "NeuralNetwork2D"
force_scales = [50.0]
nntype = "linear"

# ! define H0
# Trained with nx, ny = 10, 5

# reset_default_graph()
H0 = Variable(zeros(3,3))
H0 = H0 + H0'
M = [1.0 1.0 0.0; 1.0 1.0 0.0; 0.0 0.0 1.0]
H0 = M .* H0
n_data = [100, 101, 102, 103, 104, 200, 201, 202, 203, 204]
porder = 2
# density 4.5*(1 - 0.25) + 3.2*0.25
fiber_fraction = 0.25
#todo
#fiber_fraction = 1.0
prop = Dict("name"=> testtype, "rho"=> 4.5*(1 - fiber_fraction) + 3.2*fiber_fraction, "nn"=>nn)

T = 200.0
NT = 200

# DNS computaional domain
fiber_size = 5
# nx_f, ny_f = 40*fiber_size, 20*fiber_size
nx_f, ny_f = 80*fiber_size, 40*fiber_size

# nx_f, ny_f = 12, 4

# homogenized computaional domain
# number of elements in each directions
nx, ny = 10, 5
# number of subelements in one element in each directions
sx_f, sy_f = div(nx_f,nx), div(ny_f,ny)

ndofs = 2
fine_to_coarse = zeros(Int64, ndofs*(nx*porder+1)*(ny*porder+1))
for idof = 1:ndofs
    for iy = 1:ny*porder+1
        for ix = 1:nx*porder+1
            fine_to_coarse[ix + (iy - 1)*(nx*porder+1) + (idof-1)*(nx*porder+1)*(ny*porder+1)] =
            1 + (ix - 1) * sx_f + (iy - 1) * (nx_f*porder + 1) * sy_f + (nx_f*porder + 1)*(ny_f*porder + 1)*(idof - 1)
        end
    end
end

# #todo only for first order
# function compute_fine_to_coarse_fext(tid)
#     @assert(porder == 1)
#     # Attention fix left
#     if div(tid,100)==1 # fix bottom
#         fine_to_coarse_fext = zeros(Int64, ndofs*(nx + 1)* ny)
#         for idof = 1:ndofs
#             for iy = 1:ny
#                 for ix = 1:nx + 1
#                     fine_to_coarse_fext[ix + (iy - 1)*(nx + 1) + (idof-1)*(nx + 1)*ny] =
#                      1 + (nx_f + 1) * (sy_f - 1)  + (iy - 1) * (nx_f + 1) * sy_f + (nx_f + 1)*ny_f*(idof - 1)
#                 end
#             end
#         end

#     elseif div(tid,100)==2 # fix left
#         fine_to_coarse_fext = zeros(Int64, ndofs*nx*(ny+1))
#         for idof = 1:ndofs
#             for iy = 1:ny+1
#                 for ix = 1:nx
#                     fine_to_coarse_fext[ix + (iy - 1)*(nx) + (idof-1)*(nx)*(ny+1)] = sx_f + (ix - 1) * sx_f + (iy - 1) * (nx_f) * sy_f + (nx_f)*(ny_f + 1)*(idof - 1)
#                 end
#             end
#         end
#     end
#     fine_to_coarse_fext
# end

function compute_loss(tid, force_scale)
    nodes, EBC, g, gt, FBC, fext, ft = BoundaryCondition(tid, nx, ny, porder; force_scale=force_scale)
    domain = Domain(nodes, elements, ndofs, EBC, g, FBC, fext)
    state = zeros(domain.neqs)
    ∂u = zeros(domain.neqs)
    globdat = GlobalData(state,zeros(domain.neqs), zeros(domain.neqs),∂u, domain.neqs, gt, ft)
    assembleMassMatrix!(globdat, domain)
    # full_state_history, full_fext_history = read_data("$(@__DIR__)/Data/order$porder/$(tid)_$(force_scale)_$(fiber_size).dat")
    full_state_history, full_fext_history = read_data("$(@__DIR__)/Data/order$porder/$(tid)_$(force_scale)_$(fiber_size).dat")

    #update state history and fext_history on the homogenized domain
    state_history = [x[fine_to_coarse] for x in full_state_history]

    fext_history = []
    setNeumannBoundary!(domain, FBC, fext)
    for i = 1:NT
        globdat.time = Δt*i
        updateDomainStateBoundary!(domain, globdat)
        push!(fext_history, domain.fext[:])
    end

    # domain.state = state_history[end]
    # visσ(domain)
    # error()
    sum_loss = DynamicMatLawLoss(domain, globdat, state_history, fext_history, nn,Δt)
end

Δt = T/NT

nodes, _, _, _, _, _, _ = BoundaryCondition(n_data[1], nx, ny, porder)
elements = []
for j = 1:ny
    for i = 1:nx
        n = (nx*porder+1)*(j-1)*porder + (i-1)porder+1
        #element (i,j)
        if porder == 1
            #   4 ---- 3
            #
            #   1 ---- 2

            elnodes = [n, n + 1, n + 1 + (nx + 1), n + (nx + 1)]
        elseif porder == 2
            #   4 --7-- 3
            #   8   9   6
            #   1 --5-- 2
            elnodes = [n, n + 2, n + 2 + 2*(2*nx+1),  n + 2*(2*nx+1), n+1, n + 2 + (2*nx+1), n + 1 + 2*(2*nx+1), n + (2*nx+1), n+1+(2*nx+1)]
        else
            error("polynomial order error, porder= ", porder)
        end
        coords = nodes[elnodes,:]
        push!(elements,SmallStrainContinuum(coords,elnodes, prop, 3))
    end
end

losses = Array{PyObject}(undef, length(n_data)*length(force_scales))
k = 1
for i in n_data
    global k
    for force_scale in force_scales
        losses[k] = compute_loss(i, force_scale)
        k += 1
    end
end

@show stress_scale^2
loss = sum(losses)

sess = tf.Session(); init(sess)
@info run(sess, loss)
# error()
for i = 1:2
    println("************************** Outer Iteration = $i ************************** ")
    BFGS!(sess, loss, 2000)
    ADCME.save(sess, "$(@__DIR__)/Data/linear_train.mat")
end

@info run(sess, H0)

@show [1335174.0968380707 326448.3267263398 0.0;
        326448.3267263398 1326879.2022994285 0.0;
        0.0 0.0 526955.763626241]


text:
# generate sample timeseries data
# a port of https://github.com/LeeDoYup/RobustSTL/blob/master/sample_generator.py
using StatsBase, Distributions

function getRandomChoice(total_length::Int, num_choice::Int;
                         replace = false)
    @assert !(replace == false && num_choice > total_length) "Draws without replacement must not exceed population size"
    sample(1:total_length, num_choice; replace = replace)
end

function getChange(change_type::String, change_level::Number)
    @assert change_type in ["random", "spike", "increase", "dip", "decrease"] "Not a valid change type"

    if change_type in ["random"]
        rand([-1, 1]) * change_level
    elseif change_type in ["dip", "decrease"]
        -1.0 * change_level
    else
        change_level
    end
end

function getSeason(season_n::Int, season_type::String, season_level::Number)
    @assert season_type in ["random", "stair"] "Not a valid season type"

    season = zeros(season_n)
    if season_type == "random"
        season = 2.0 * (rand(season_n) .- 0.5) * season_level
    else
        centre = season_n ÷ 2
        season[1:centre] .+= season_level
        season[centre + 1: end] .-= season_level
    end
    season
end

function genSeasons(total_n, season_n, season_num, season_type, season_level)
    season = getSeason(season_n, season_type, season_level)
    repeat(season, season_num)[1:total_n]
end

function genRemainder(total_n, noise_mean, noise_std)
    rand(Normal(noise_mean, noise_std), total_n)
end

function genAnomalies(total_n, anomaly_num, anomaly_type, anomaly_level)
    anomalyLocs = getRandomChoice(total_n, anomaly_num)
    anomalies = zeros(total_n)
    map(x -> anomalies[x] += getChange(anomaly_type, anomaly_level), anomalyLocs)
    anomalies
end

function genTrends(total_n, trend_change_num, trend_type, trend_level)
    trends = zeros(total_n)
    change_points = getRandomChoice(total_n, trend_change_num)
    map(x -> trends[x:end] .+= getChange(trend_type, trend_level), change_points)
    trends
end

"""
    gen_ts(total_n, season_n)

Generates a basic timeseries with trend, seasonality, noise and some anomalies.

**Parameters**:

- _total_n_: size of required dataset (e.g. 365, 52, 12, etc.).
- _season_n: number of seasonal points (e.g. 12 in case of a simulated monthly seasonality).

The following are optional parameters that can be passed to customise the individual timeseries components:

- _season_type_: set trajectory of seasonal component. `random` or `stair`. Default is `stair`.
- _season_level_: set amplitude of the seasonal effect (essentially acts as a static multiplier). Default is `1`.

- _trend_type_: set trajectory of trend component. `random` (default), `increase`, or `decrease`.
- _trend_level_: set amplitude of the trend component. Default is `3`
- _trend_change_num_: set how many instances of changes in the trend to be present. These are randomly distributed across the dataset. Default is `10`

- _anomaly_num_: set how many anomalies are to be randomly distributed across the dataset. Default is `6`.
- _anomaly_type_: set trajectory of the anomalies. `random` (default), `increase`, or `decrease`.
- _anomaly_level_: set amplitude of the anomalies. Default is `4`.

- _noise_mean_: mean of the normaly distributed noise component. Default is `0`.
- _noise_std_: standard deviation of the noise component. Default is `0.316`.

Returns a dataframe with individual components as well as the composite `sample`.
"""
function gen_ts(total_n::Number, season_n::Number;
    season_type = "stair",
    season_level = 1,
    trend_type = "random",
    trend_level = 3,
    trend_change_num = 10,
    anomaly_num = 6,
    anomaly_type = "random",
    anomaly_level = 4,
    noise_mean = 0,
    noise_std = 0.316)

    @assert total_n >= season_n "Dataset length must be equal or greater to seasonality"

    season_num = (total_n ÷ season_n) + 1

    trends = genTrends(total_n, trend_change_num, trend_type, trend_level)
    seasons = genSeasons(total_n, season_n, season_num, season_type, season_level)
    remainders = genRemainder(total_n, noise_mean, noise_std)
    anomalies = genAnomalies(total_n, anomaly_num, anomaly_type, anomaly_level)
    sample = trends + seasons + remainders + anomalies

    DataFrame(sample = sample,
            trends = trends,
            seasonality = seasons,
            remainder = remainders + anomalies)
end

"""
    gen_ts()

Generates a basic timeseries with trend, seasonality, noise and some anomalies.

Returns a dataframe with individual components as well as the composite `sample`.
"""
function gen_ts()
    gen_ts(365, 30)
end