text:
"""
    add_day(start, duration)

Calculate dates by adding `duration` (number of days) to the `start` date.  `duration` could be
decimal.

# Examples
```julia-repl
julia> start = Dates.DateTime(Dates.now())
2021-01-01T15:41:31.91

julia> later = add_day(start, 1.50)
2021-01-03T03:41:31.91
```
"""
function add_day(start::Dates.DateTime, duration)
    #d = start + Dates.Millisecond(Int(duration*86400*1000)) # convert days to ms

    d = start + Dates.Millisecond( # convert days to ms
            # Not sure which one is better
            #convert(Int64, round(duration*86400*1000, digits=4)) # 4 is from trial&error, no idea why
            Int64(trunc(duration*86400*1000))
        )
    return d
end

# PRT file has the following lines
#
#SECTION  Starting the simulation on 01-Feb-2019.
#SECTION  The simulation has reached 27-Feb-2024 1852 d.
#SECTION  The simulation has reached 28-Feb-2024 1853 d.
#SECTION  Simulation complete.

"""
    find_prt_start_date(io; to_datetime = true)

Find simulation start date from a PRT file.  Return `DateTime` if `to_datetime`
is true otherwise a string in the "d-u-Y" format such as "01-Feb-2019".

# Examples
```julia-repl
julia> start_datetime = find_prt_start_date(prt_file)
2019-02-01T00:00:00
```
"""
function find_prt_start_date(io::IOStream; to_datetime = true)
    for line in eachline(io)
        if !startswith(line, "SECTION  Starting the simulation on")
            continue
        end
        line = strip(line)
        #tokens = split(line, x -> x in [' ', '.'], keepempty=false)
        tokens = split(line, [' ', '.'], keepempty=false)

        startdate_token = tokens[6] # "01-Feb-2019"
        if to_datetime
            return Dates.DateTime(startdate_token, Dates.dateformat"d-u-Y")
        else
            return startdate_token
        end

    end

    nothing
end
function find_prt_start_date(path::String; to_datetime = true)
    open(path, "r") do io
        return find_prt_start_date(io, to_datetime=to_datetime)
    end
end

"""
    find_prt_current_time(io::IOStream, pos::Int64)
    find_prt_current_time(path::String, pos::Int64)

Find the current time in days in floating point number from the PRT file.
The function looks for a snippet like this.
```
LOG                    TIME  TSTEP       GOR      WCT      OPR  ....
                          d      d  MSCF/STB  STB/STB    STB/d  ....
          ----------------------------------------------------- ....
          SCT    ;   16.450  0.450      1.37    0.950   17.116  ....
```
"""
function find_prt_current_time(io::IOStream, pos::Int64)
    seek(io, pos) # re-position, used skip(io, pos) before

    # Find the following line to get TIME and TSTEP:
    #
    #LOG                    TIME  TSTEP       GOR      WCT      OPR      WPR      GPR      FPR      WIR      GIR  ITER  IMPL
    #                          d      d  MSCF/STB  STB/STB    STB/d    STB/d   MSCF/d      psi    STB/d   MSCF/d           %
    #          --------------------------------------------------------------------------------------------------------------
    #          SCT    ;   16.450  0.450      1.37    0.950   17.116  323.707   23.449  6102.16    0.000    0.000     6   100

    time_token  = nothing
    tstep_token = nothing
    code_token  = nothing
    gor_token   = nothing
    opr_token   = nothing
    wpr_token   = nothing
    gpr_token   = nothing
    for line in eachline(io)
        # Check if we are at the right place
        if !startswith(line, "LOG") continue end
        if length(line) < 34        continue end
        tokens = split(line, [' '], keepempty=false)
        if !((tokens[1] == "LOG") && (tokens[2] == "TIME") && (tokens[3] == "TSTEP"))
            continue
        end

        # skip 2 lines
        # "" == readline(io) when reach EOF
        if "" == readline(io) continue end
        if "" == readline(io) continue end

        #LOG                    TIME  TSTEP       GOR      WCT      OPR      WPR      GPR      FPR      WIR      GIR  ITER  IMPL
        #                  d      d  MSCF/STB  STB/STB    STB/d    STB/d   MSCF/d      psi    STB/d   MSCF/d           %
        # --------------------------------------------------------------------------------------------------------------
        # Read this line
        #          SCT    ;   16.450  0.450      1.37    0.950   17.116  323.707   23.449  6102.16    0.000    0.000     6   100
        line = readline(io)
        if "" == line continue end
        tokens = split(line, [' '], keepempty=false)
        code_token  = tokens[1] # Time step code, see IX pp 1315, Reporting node reference | Reporting property identifiers | Simulation engine properties | Simulation engine properties requiring special description | TIMESTEP_CONTROL_MODE
        time_token  = tokens[3] # 16.450 Days
        tstep_token = tokens[4] # 0.450 Days
        gor_token   = tokens[5]
        opr_token   = tokens[7]
        wpr_token   = tokens[8]
        gpr_token   = tokens[9]
    end

    if isnothing(time_token)
        return nothing
    end

    time = parse(Float64, time_token)
    tstep = parse(Float64, tstep_token)
    gor   = parse(Float64, gor_token)
    opr   = parse(Float64, opr_token)
    wpr   = parse(Float64, wpr_token)
    gpr   = parse(Float64, gpr_token)

    return (time=time, tstep=tstep, code=code_token, gor=gor, opr=opr, gpr=gpr, wpr=wpr)
end
function find_prt_current_time(path::String, pos::Int64)
    open(path, "r") do io
        return find_prt_current_time(io, pos)
    end
end

function find_prt_current_date(io::IOStream, pos::Int64; to_datetime = true)
    seek(io, pos) # re-position, used skip(io, pos) before

    # Possible SECTION lines:
    #
    #SECTION  Starting the simulation on 18-Jan-2020
    #SECTION  The simulation has reached 26-Jan-2021 374 d. ...
    #SECTION  Simulation complete.

    date_token = nothing
    for line in eachline(io)
        if !startswith(line, "SECTION")
            continue
        end
        if length(line) < 48
            continue
        end
        if startswith(line, "SECTION  The simulation has reached")
            tokens = split(line, [' '], keepempty = false)
            date_token = tokens[6]

            #days = convert(Day, datetime - start_datetime)
        end
    end

    if isnothing(date_token)
        #throw(ErrorException("No DATE found"))
        return nothing
    end

    if to_datetime
        return Dates.DateTime(date_token, Dates.dateformat"d-u-Y")
    else
        return date_token
    end
end
function find_prt_current_date(path::String, pos::Int64, to_datetime = true)
    open(path, "r") do io
        return find_prt_current_date(io, pos, to_datetime=to_datetime)
    end
end

"""
    find_schedule_end_date(io::IOStream; to_datetime=true)
    find_schedule_end_date(path::String; to_datetime = true)

Find the last date in the schedule.  Returns value in `DateTime` or `String`.
"""
function find_schedule_end_date(io::IOStream; to_datetime=true)
    date_token = nothing
    for line in eachline(io)
        # DATE "04-Feb-2021"
        if !startswith(line, "DATE")
            continue
        end

        line = strip(line)
        tokens = split(line, [' ', '"'], keepempty=false)
        date_token = tokens[2]
    end

    if isnothing(date_token)
        #throw(ErrorException("No DATE found"))
        return nothing
    end

    if to_datetime
        return Dates.DateTime(date_token, Dates.dateformat"d-u-Y")
    else
        return date_token
    end

end
function find_schedule_end_date(path::String; to_datetime = true)
    open(path, "r") do io
        return find_schedule_end_date(io, to_datetime=to_datetime)
    end
end

function read_hdf5_times(filename)
    h5open(filename, "r") do h5
        times = h5["general"]["time"]
        return times[:]
    end
end

# RSM file
# "1" marks the beginning of a section of data and starts at column 1.
# Data starts at column 2
#
# Use index to access headers.  There are 10 or less columns.  Only the first batch has YEARS column.
# Column indices are [2, 17, 32, 47, 62, 77, 92, 107, 122, 137]
# if the first blank is removed the indices are [1, 16, 31, 46, 61, 76, 91, 106, 121, 136]
# add 14 to reach the end of column so grep [2, 2+14].
# 1st row of header is unique and can be used as column names in DataFrame.
# 2nd row is the unit.
# 3rd and 4th row just call them 3rd and 4th rows.
#
# Use split to read body.
#
#
#1
# ------------------------------------------------------------------------------------------------------------------------------------------------------
# SUMMARY OF RUN Pad-X_51_19 :: Run Date and Time : XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
# ------------------------------------------------------------------------------------------------------------------------------------------------------
# TIME           YEARS          FPPW           FPPO           FPPG           FNQT           FNQR           FEIP           FWIT           FWIR
# DAYS           YEARS          PSIA           PSIA           PSIA           STB            STB/DAY        BTU            STB            STB/DAY
#                               FIELD          FIELD          FIELD          FIELD          FIELD          FIELD          FIELD          FIELD
#                               1              1              1                                            1
# ------------------------------------------------------------------------------------------------------------------------------------------------------
# 0              0              5316.56        6415.72        7839.82        0              0              0              0              0
# 1              0.00273785     5316.56        6415.62        7839.75        0              0              0              0              0
# 4              0.0109514      5316.68        6420.9         7843.84        0              0              0              1.20911        0.403038
#1
# ------------------------------------------------------------------------------------------------------------------------------------------------------
# SUMMARY OF RUN Pad-X_51_19 :: Run Date and Time : XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
# ------------------------------------------------------------------------------------------------------------------------------------------------------
# TIME           WOPR           WOPR           WOPT           WOPT           WOPT           WOPT
# DAYS           STB/DAY        STB/DAY        STB            STB            STB            STB
#                ADA-1763       ADA-1764       ADA-1761       ADA-1762       ADA-1763       ADA-1764
#
# ---------------------------------------------------------------------------------------------------------
# 0              0              0              0              0              0              0
# 1              0              0              0              0              0              0
# 4              0              0              0              0              0              0
"""
    read_rsm(io::IOStream)
    read_rsm(filename::String)

Parse RSM file and return (body=DataFrame, meta=DataFrame).

Parse
```
1
 ------------------------------------------------------------------------------------------------------------------------------------------------------
 SUMMARY OF RUN Pad-X_51_19 :: Run Date and Time : XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
 ------------------------------------------------------------------------------------------------------------------------------------------------------
 TIME           YEARS          FPPW           FPPO           FPPG           FNQT           FNQR           FEIP           FWIT           FWIR
 DAYS           YEARS          PSIA           PSIA           PSIA           STB            STB/DAY        BTU            STB            STB/DAY
                               FIELD          FIELD          FIELD          FIELD          FIELD          FIELD          FIELD          FIELD
                               1              1              1                                            1
 ------------------------------------------------------------------------------------------------------------------------------------------------------
 0              0              5316.56        6415.72        7839.82        0              0              0              0              0
 1              0.00273785     5316.56        6415.62        7839.75        0              0              0              0              0
 4              0.0109514      5316.68        6420.9         7843.84        0              0              0              1.20911        0.403038
1
 ------------------------------------------------------------------------------------------------------------------------------------------------------
 SUMMARY OF RUN Pad-X_51_19 :: Run Date and Time : XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
 ------------------------------------------------------------------------------------------------------------------------------------------------------
 TIME           WOPR           WOPR           WOPT           WOPT           WOPT           WOPT
 DAYS           STB/DAY        STB/DAY        STB            STB            STB            STB
                ADA-1763       ADA-1764       ADA-1761       ADA-1762       ADA-1763       ADA-1764

 ---------------------------------------------------------------------------------------------------------
 0              0              0              0              0              0              0
 1              0              0              0              0              0              0
 4              0              0              0              0              0              0
```

and return a named tuple (body=DataFrame, meta=DataFrame)

`body`
```
  Row │ TIME       YEARS       FPPW     FPPO     FPPG
      │ Float64    Float64     Float64  Float64  Float64
──────┼───────────────────────────────────────────────────
    1 │   0.0      0.0         5316.25  6415.77  7839.76
    2 │   1.0      0.00273785  5316.22  6415.66  7839.66
    3 │   1.35154  0.00370032  5322.81  6434.31  7855.26
    4 │   1.66463  0.0045575   5326.12  6446.29  7865.1
    5 │   2.30742  0.00631736  5331.27  6462.0   7878.14
```

`meta`
```
 Row │ column       1         2         3         4
     │ String       String    String    String    String
─────┼───────────────────────────────────────────────────────
   1 │ TIME         TIME      DAYS
   2 │ YEARS        YEARS     YEARS
   3 │ FPPW         FPPW      PSIA      FIELD     1
   4 │ FPPO         FPPO      PSIA      FIELD     1
 ...
 199 │ WOPR         WOPR      STB/DAY   ADA-1761
 200 │ WOPR_3       WOPR      STB/DAY   ADA-1762
 201 │ WOPR_1       WOPR      STB/DAY   ADA-1763
 202 │ WOPR_2       WOPR      STB/DAY   ADA-1764
```

 where
 * `column` is the unique column name in `body`
 * "1" is the first row in RSM header
 * "2" is the second row in RSM header
 * "3" is the third row in RSM header
 * "4" is the fourth row in RSM header
"""
function read_rsm(io::IOStream)
    line = readline(io) # 1

    # Read first SUMMARY
    nt = read_rsm_section(io) # (body = df_body, meta = df_meta)
    df_body = nt.body
    df_meta = nt.meta

    # Append the rest of SUMMARY
    while !eof(io)
        nt = read_rsm_section(io)

        # Must have the same number or rows
        @assert size(df_body)[1] == size(nt.body)[1] "SUMMARY sections has different number of rows"

        # skip the first column which is DAYS again
        df_body = hcat(df_body, nt.body[!,2:end], makeunique=true)
        df_meta = hcat(df_meta, nt.meta[!,2:end], makeunique=true)
    end

    # assert
    @assert names(nt.meta) == names(nt.body) "DataFrame generates different column names?"

    # Transpose meta table
    # https://stackoverflow.com/questions/37668312/transpose-of-julia-dataframe
    df = df_meta
    df_meta = DataFrame([[names(df)]; collect.(eachrow(df))], [:column; Symbol.(axes(df, 1))])

    return (body = df_body, meta = df_meta)
end
function read_rsm(filename)
    nt = nothing
    open(filename, "r") do io
        nt = read_rsm(io)
    end

    return nt
end

"""
    read_rsm_section(io::IOStream)

Reads a section of RSM file.

A section could end with "1"
```
 ------------------------------------------------------------------------------------------------------------------------------------------------------
 SUMMARY OF RUN Pad-X_51_19 :: Run Date and Time : XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
 ------------------------------------------------------------------------------------------------------------------------------------------------------
 TIME           YEARS          FPPW           FPPO           FPPG           FNQT           FNQR           FEIP           FWIT           FWIR
 DAYS           YEARS          PSIA           PSIA           PSIA           STB            STB/DAY        BTU            STB            STB/DAY
                               FIELD          FIELD          FIELD          FIELD          FIELD          FIELD          FIELD          FIELD
                               1              1              1                                            1
 ------------------------------------------------------------------------------------------------------------------------------------------------------
 0              0              5316.56        6415.72        7839.82        0              0              0              0              0
 1              0.00273785     5316.56        6415.62        7839.75        0              0              0              0              0
 4              0.0109514      5316.68        6420.9         7843.84        0              0              0              1.20911        0.403038
1
```

or `eof`

```
 ------------------------------------------------------------------------------------------------------------------------------------------------------
 SUMMARY OF RUN Pad-X_51_19 :: Run Date and Time : XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
 ------------------------------------------------------------------------------------------------------------------------------------------------------
 TIME           YEARS          FPPW           FPPO           FPPG           FNQT           FNQR           FEIP           FWIT           FWIR
 DAYS           YEARS          PSIA           PSIA           PSIA           STB            STB/DAY        BTU            STB            STB/DAY
                               FIELD          FIELD          FIELD          FIELD          FIELD          FIELD          FIELD          FIELD
                               1              1              1                                            1
 ------------------------------------------------------------------------------------------------------------------------------------------------------
 0              0              5316.56        6415.72        7839.82        0              0              0              0              0
 1              0.00273785     5316.56        6415.62        7839.75        0              0              0              0              0
 4              0.0109514      5316.68        6420.9         7843.84        0              0              0              1.20911        0.403038
```
"""
function read_rsm_section(io::IOStream)
    # starts right after "1"
    # ends with "1" or eof

    colinds = [2, 17, 32, 47, 62, 77, 92, 107, 122, 137]
    delta = 14
    delimiters = [' ']
    #delimiters = [' ', ',', '\t']

    readline(io) # ---------------...
    readline(io) # SUMMARY OF RUN ...
    readline(io) # ---------------...

    # Column name:  TIME           YEARS          FPPW           FPPO           FPPG           FNQT           FNQR           FEIP           FWIT           FWIR
    line = readline(io)
    tokens = split(line , delimiters, keepempty=false)
    col_count = length(tokens)
    row_1 = tokens

    # Create DataFrame
    df_body = DataFrame(fill(Float64[], col_count), row_1, makeunique=true)
    df_meta = DataFrame(fill(String[],  col_count), row_1, makeunique=true)

    push!(df_meta, row_1)

    # units:  DAYS           YEARS          PSIA           PSIA           PSIA           STB            STB/DAY        BTU            STB            STB/DAY
    line = readline(io)
    tokens = String[]
    for colind in colinds
        if colind > length(line)
            break
        end
        token = strip(line[colind:colind+delta])
        push!(tokens, token)
    end
    row_2 = tokens
    push!(df_meta, row_2)

    # 3rd row:
    line = readline(io)
    tokens = String[]
    for colind in colinds
        if colind > length(line)
            break
        end
        token = strip(line[colind:colind+delta])
        push!(tokens, token)
    end
    row_3 = tokens
    push!(df_meta, row_3)

    # 4th row:
    line = readline(io)
    tokens = String[]
    for colind in colinds
        if colind > length(line)
            break
        end
        token = strip(line[colind:colind+delta])
        push!(tokens, token)
    end
    row_4 = tokens
    push!(df_meta, row_4)

    readline(io) # ---------------...

    # Body
    while !eof(io)
        line = readline(io)

        if "1" == line
            break
        end
        if isempty(line)
            break
        end
        tokens = split(line, delimiters, keepempty=false)
        tokens = parse.(Float64, tokens)
        push!(df_body, tokens)
    end

    return (body = df_body, meta = df_meta)
end

"""

    find_column_name(df_meta, c1, c3, c4)

Filter rows of the meta table based on column 1, 3, 4.  Check if there is
1 and only 1 row, then return [1,1] that is "WOPR_3" in the example.

```
 Row │ column  1       2        3         4
     │ String  String  String   String    String
─────┼───────────────────────────────────────────
   1 │ WOPR_3  WOPR    STB/DAY  ADA-1762
```
"""
function find_column_name(df_meta, c1, c3, c4)
    df = filter(
        ["1", "3", "4"] => (x1,x3,x4) -> c1 == x1 && c3 == x3 && c4 == x4,
        df_meta
    )

    @assert size(df)[1] == 1 "None or multiple column names found"

    return df[1,1]
end

"""
    skip_grdecl_keyword_data!(io)

Skip the content of a keyword after the keyword has been read and determined
the content should be skipped.
"""
function skip_grdecl_keyword_data!(io)
    time_to_break = false

    for line in eachline(io)
        if time_to_break
            break
        end
        if endswith(line, '/')
            time_to_break = true
        end
    end
end

"""
    read_grdecl_keyword_data!(io::IO, values::Array{T}) where {T<:String}

Read keyword data as `String`.
"""
function read_grdecl_keyword_data!(io::IO, values::Array{T}) where {T<:String}
    # Read tokens
    tokens = String[]
    for line in eachline(io)
        line = strip(line)
        if isempty(line)           continue  end
        if startswith(line, "--")  continue  end

        append!(tokens, split(line))

        if endswith(line, '/')  break  end
    end

    values .= tokens

    return nothing
end

"""
    read_grdecl_keyword_data!(io::IO, values::Array{T}) where {T<:Real}

Read keyword data as subtypes of `Real`.
Read values of a keyword.  Notice the keyword should has been read before
calling this function.
"""
function read_grdecl_keyword_data!(io::IO, values::Array{T}) where {T<:Real}
    # Read tokens as strings
    tokens = Vector{String}()
    for line in eachline(io)

        line = strip(line)
        if isempty(line)           continue  end
        if startswith(line, "--")  continue  end

        append!(tokens, split(line))

        if endswith(line, '/')  break  end
    end

    # remove "/"
    tokens = tokens[1:end-1]

    # Parse tokens into Float64
    ind = 1
    for token in tokens
        count_value = split(token,"*")
        if length(count_value) == 2
            count = parse(Int64,    count_value[1])
            value = parse(T, count_value[2]) # eltype(T)?
            for x = 1:count
                values[ind] = value
                ind += 1
            end
        else
            values[ind] = parse(T, count_value[1]) # eltype(T)?
            ind += 1
        end
    end

    return nothing
end

"""
    read_grdecl_keyword_data!(io::IO, values::Array{T}, prop::String) where {T<:Real}

Read static data.
Read data of a given keyword as a subtype of `Real`.
"""
function read_grdecl_keyword_data!(io::IO, values::Array{T}, prop::String) where {T<:Real}
    for line in eachline(io)
        line = strip(line)

        if isempty(line)          continue end
        if startswith(line, "--") continue end
        tokens = split(line)
        if isempty(tokens)        continue end # will this every happen?

        keyword = tokens[1]
        if keyword != prop
            skip_grdecl_keyword_data!(io)
        else
            read_grdecl_keyword_data!(io, values)

            return true
        end
    end

    return false
end

"""

Read dynamic data (with timestamp)

# Arguments
- `values`: Pre-allocated array to save values to
- `prop`: Property name
- `date`: Date
"""
function read_grdecl_keyword_data!(io::IO, values::Array{T}, prop::String, date) where {T<:Real}
    for line in eachline(io)
        line = strip(line)
        if startswith(line, "--")
            continue
        end
        if isempty(line)
            continue
        end

        tokens = split(line)
        if isempty(tokens)
            continue
        end

        keyword = tokens[1]
        if keyword != prop
            skip_grdecl_keyword_data!(io)
        else
            # Read timestep
            line = readline(io)
            tokens = split(line, [' ', ','], keepempty=false)
            # last 4 tokens are time: Jan 07,2022 00:00:00
            timee = tokens[end]
            yearr = tokens[end-1]
            datee = tokens[end-2]
            monthh = tokens[end-3]

            #current_date = Date("$yearr-$monthh-$datee", dateformat"y-u-d") # 2022-Jan-07
            current_date = Date("$yearr-$monthh-$(datee)T$(timee)", dateformat"y-u-dTHH:MM:SS") # 2022-Jan-07
            #DateTime("2021-Jun-18T15:38:17", DateFormat("y-u-dTHH:MM:SS"))

            if current_date != date
                skip_grdecl_keyword_data!(io)
            else
                read_grdecl_keyword_data!(io, values)

                return true
            end
        end
    end

    return false
end


text:
#=
Copyright 2017 INSIGNEO Institute for in silico Medicine

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

# All the pre-processing and initialisation functions are herein contained.

# *function* __`projectPreamble`__
#
# ----------------------------------------------------------------------------
# Parameters:
# ---------------- -----------------------------------------------------------
# `project_name`   `::String` project name.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning:
# ----------------------------------------------------------------------------
# This function checks whether all the project files are present in the
# working directory. Then, the directory structure to save the output files
# is created. Eventually, `openBF` logo and version information are printed to
# screen.
# ----------------------------------------------------------------------------
#
function projectPreamble(project_name :: String, no_out :: Bool, no_inputs :: Bool)

  # `project_constants.jl`, `project_inlet.dat`, and `project.csv` file names
  # strings are created with `project_name` variable.
  p_const = join([project_name, "_constants.jl"])
  p_inlet = join([project_name, "_inlet.dat"])
  p_model = join([project_name, ".csv"])
  v_model = join([project_name, "_veins.csv"])
  # The existence of these three files is checked in the working directory.
  # If a file is missing an error is raised.
  if (isfile(p_const)) == false
    error("$p_const file missing")
  end

  if (isfile(p_inlet)) == false
    println("$p_inlet file missing")
  end

  if (isfile(p_model)) == false
    error("$p_model file missing")
  end
  # `r_folder` is the string used to create the directory containing the
  # results.
  r_folder = join([project_name, "_results"])
  # If the results directory does not exist, it is created.
  if (isdir(r_folder)) == false
    mkdir(r_folder)
  end
  # Project files are moved to the results directory.
  cp(p_const, join([r_folder, "/", p_const]), remove_destination=true)
  cp(p_inlet, join([r_folder, "/", p_inlet]), remove_destination=true)
  cp(p_model, join([r_folder, "/", p_model]), remove_destination=true)
  if (isfile(v_model)) == true
    cp(v_model, join([r_folder, "/", v_model]), remove_destination=true)
  end
  # All files in the working directory are deleted except those containing
  # the `project_name` string by the
  # [`cleanLibrary`](initialise.html#cleanLibrary) function.
  # cleanLibrary(project_name)
  # run(`rm clean_lib.sh`)
  # The working directory is changed to the results directory.
  cd(r_folder)
  # Bash scripts to handle I/O are written and saved in the working directory
  # by [`writeScripts`](initialise.html#writeScripts) function.
  # writeScripts(no_out, no_inputs)

end

# Two bash scripts are written in two `.sh` files.
#
# *   `appender.sh` is used to concatenate two files which names should be
#     given as script arguments. It can be used as
#
#         sh appender.sh file1 file2
#
# *   `cleaner.sh` is used at the end of the simulation to remove all the
#     temporary files from the working directory.
#
# These `.sh` are saved in the working directory and used through the `julia`
# command `run`.

# *function* __`writeScripts`__
#
# ----------------------------------------------------------------------------
# Functioning
# ----------------------------------------------------------------------------
# Two bash scripts are written in two `.sh` files.
# ----------------------------------------------------------------------------
#
# *   `appender.sh` is used to concatenate two files which names should be
#     given as script arguments. It can be used as
#
#         sh appender.sh file1 file2
#
# *   `cleaner.sh` is used at the end of the simulation to remove all the
#     temporary files from the working directory.
#
# These `.sh` are saved in the working directory and used through the `julia`
# command `run`.
#
function writeScripts(no_out :: Bool, no_inputs :: Bool)

  appsh = open("appender.sh", "w")
  write(appsh, "#!/bin/bash", "\n")
  write(appsh, "cat \$1 >> \$2")
  close(appsh)

  # clesh = open("cleaner.sh", "w")
  # write(clesh, "#!/bin/bash", "\n")
  # write(clesh, "rm *.temp\n")
  # write(clesh, "rm appender.sh")
  # if no_out == true
  #     write(clesh, "\nrm *.out")
  # end
  # if no_inputs == true
  #     write(clesh, "\nrm *.csv")
  #     write(clesh, "\nrm *.dat")
  #     write(clesh, "\nrm *.jl")
  # end
  # close(clesh)

end

# *function* __`cleanLibrary`__
#
# ----------------------------------------------------------------------------
# Parameters:
# ---------------- -----------------------------------------------------------
# `project_name`   `::String` project name.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning:
# ----------------------------------------------------------------------------
# This function writes and executes a bash scripts which delete all the
# files with the name not containing the string `project_name`.
# ----------------------------------------------------------------------------
#
function cleanLibrary(project_name)

  clibsh = open("clean_lib.sh", "w")
  write(clibsh, "#!/bin/bash", "\n")
  write(clibsh, "find . -type f -not -name '*$project_name*' | xargs rm -rf", "\n")
  write(clibsh, "rm -rf doc/")
  close(clibsh)

  #run(`sh clean_lib.sh`)

end

# *function* __`readModelData`__ $\rightarrow$ `model_matrix::Array{Any, 2}`
#
# ----------------------------------------------------------------------------
# Parameters:
# ---------------- -----------------------------------------------------------
# `model_data`     `.csv` file containing model data (see
#                  [tutorial.jl](../index.html#project)).
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning:
# ----------------------------------------------------------------------------
# `model_data.csv` is parsed by `readdlm` by using `,` as a separator. The
# first header row is discarded (`[2:end]`) as well as the last column
# (`[1:end-1]`) because it contains
# only blank spaces. The remaining is stored in matrix `m` and returned.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Returns:
# ---------------- -----------------------------------------------------------
# `m`              `::Array{Any, 2}` model matrix.
# ----------------------------------------------------------------------------
#
function readModelData(model_data)

  m = readdlm(model_data, ',')
  m = m[2:end, 1:end-1]

  return m
end

# *function* __`initialiseVessel`__ $\rightarrow$
# `vessel::`[`Vessel`](BTypes.html#Vessel)
#
# ----------------------------------------------------------------------------
# Parameters:
# ------------------- --------------------------------------------------------
# `m`                 `::Array{Any, 2}` model matrix row.
#
# `ID`                `::Int` topological ID of the vessel in the arterial
#                     tree.
#
# `h`                 [`::Heart`](BTypes.html#Heart) data structure
#                     containing inlet boundary conditions.
#
# `b`                 [`::Blood`](BTypes.html#Blood) data structure
#                     containing blood mechanical properties.
#
# `Pext`              `::Float` external pressure.
#
# `initial_pressure`  `::Float` initial pressure to be set along the entire
#                     vessel.
#
# `Ccfl`              `::Float` Courant-Friedrichs-Lewy number to [compute
#                     local $\Delta t$](godunov.html#calculateDeltaT).
# ----------------------------------------------------------------------------
#
function initialiseVessel(m :: Array{Any, 1}, ID :: Int64, h :: Heart,
                              b :: Blood,
                              initial_pressure :: Float64,
                              Ccfl :: Float64)

  # --------------------------------------------------------------------------
  # Functioning:
  # --------------------------------------------------------------------------
  # Model matrix row `m` is parsed  to retrieve vessel geometrical and
  # mechanical properties; conversion to `::Int` type is operated where
  # needed.
  # --------------------------------------------------------------------------
  vessel_name = m[1]
  sn = convert(Int, m[2])
  tn = convert(Int, m[3])
  rn = convert(Int, m[4])
  L = m[5]

  #M = convert(Int, m[6])
  M = maximum([5, convert(Int, m[6]), convert(Int, ceil(L*1e3))])

  Rp = m[7]
  Rd = m[8]

#   h0          =     m[7]

  E = m[9]
  Pext = convert(Float64, m[10])
  # Poisson's ratio `sigma` is set by default to 0.5 because vessel wall is
  # assumed to be incompressible.
  sigma = 0.5
  # $\Delta x$ for local numerical discretisation: lenght/number of nodes.
  #
  dx = L/M
  invDx = M/L
  halfDx = 0.5*dx

  # Unstressed cross-sectional area `A0` is computed from `R0`.
#   A0 = pi*R0*R0

#   Rp = 4e-3
#   Rd = 2e-3

#   Rp = r0
#   Rd = r0

  A0 = zeros(Float64, M)
  R0 = zeros(Float64, M)
  h0 = zeros(Float64, M)

  dA0dx = zeros(Float64, M)
  dTaudx = zeros(Float64, M)
  radius_slope = (Rd-Rp)/(M-1)
  ah = 0.2802
  bh = -5.053e2
  ch = 0.1324
  dh = -0.1114e2
  for i = 1:M
    R0[i] = radius_slope*(i-1)*dx + Rp
    h0[i] = R0[i] * (ah * exp(bh*R0[i]) + ch * exp(dh*R0[i])) #1e-3
    A0[i] = pi*R0[i]*R0[i]
    dA0dx[i] = 2*pi*R0[i]*radius_slope
    dTaudx[i] = sqrt(pi)*E*radius_slope*1.3*(h0[i]/R0[i]+R0[i]*(ah*bh*exp(bh*R0[i]) + ch*dh*exp(dh*R0[i])))
  end

#   println(dA0dx[end])

  # Elastic constants for trans-mural [pressure](converter.html#pressure)
  # (`beta`) and [wave speed](converter.html#waveSpeed) (`gamma`)
  # calculation:
  # $$
  #   \beta = \sqrt{\frac{\pi}{A_0}} \frac{h_0 E}{1 - \sigma^2}, \quad
  #   \gamma = \frac{\beta}{3 \rho R_0 \sqrt{\pi}} .
  # $$
  beta  = sqrt.(pi./A0) .* h0*E / (1 - sigma*sigma)
  gamma = beta ./ (3*b.rho*R0*sqrt(pi))

  gamma_ghost = zeros(Float64, M+2)
  gamma_ghost[2:M+1] = gamma
  gamma_ghost[1] = gamma[1]
  gamma_ghost[end] = gamma[end]

  # Initialise conservative (`A` and `Q`) and primitive (`u`, `c`, and `P`)
  # variables arrays with initial values. `A` is taken as the unstressed
  # cross-sectional area `A0` and `Q` is everywhere the `initial_flow`
  # specified in [`project_constants.jl`](../index.html#project_constants).
  # The longitudinal velocity `u` comes from the definition of `Q`
  # $$
  #     Q = u A .
  # $$
  # `c` and `P` are computed with `beta` and `gamma` by
  # [`pressure`](converter.html#pressure) and
  # [`waveSpeed`](converter.html#waveSpeed) functions.
  A = zeros(Float64, M)  + A0
  Q = zeros(Float64, M)  + h.initial_flow
  u = zeros(Float64, M)  + Q./A
  c = zeros(Float64, M)
  c = waveSpeed(A, gamma, c)
  P = zeros(Float64, M)
  P = pressure( A, A0, beta, Pext, P)
  # [Ghost cells](godunov.html) are initialised
  # as the internal nodes.
  U00A = A0[1]
  U01A = A0[2]
  UM1A = A0[M]
  UM2A = A0[M-1]

  U00Q = h.initial_flow
  U01Q = h.initial_flow
  UM1Q = h.initial_flow
  UM2Q = h.initial_flow
  # Forward and backward
  # [Riemann invariants](converter.html#riemannInvariants)
  # at the vessel outlet node.
  W1M0 = u[end] - 4*c[end]
  W2M0 = u[end] + 4*c[end]
  # Temporary files names are built by `join`ing the vessel `label` and a
  # string referring to the quantity stored, and a `.temp` extension. Same
  # procedure is followed for results files with `.out` extension. Data input
  # and output are handled by [`IOutils.jl`](IOutils.html) functions.
  temp_A_name = join((vessel_name,"_A.temp"))
  temp_Q_name = join((vessel_name,"_Q.temp"))
  temp_u_name = join((vessel_name,"_u.temp"))
  temp_c_name = join((vessel_name,"_c.temp"))
  temp_P_name = join((vessel_name,"_P.temp"))

  last_A_name = join((vessel_name,"_A.last"))
  last_Q_name = join((vessel_name,"_Q.last"))
  last_u_name = join((vessel_name,"_u.last"))
  last_c_name = join((vessel_name,"_c.last"))
  last_P_name = join((vessel_name,"_P.last"))

  out_A_name = join((vessel_name,"_A.out"))
  out_Q_name = join((vessel_name,"_Q.out"))
  out_u_name = join((vessel_name,"_u.out"))
  out_c_name = join((vessel_name,"_c.out"))
  out_P_name = join((vessel_name,"_P.out"))
  # Temporary data and output files are created with `open` command. Data from
  # temporary files to output files are
  # [transferred](IOutils.jl#transferTempToOut) by means of a bash script.
  # Thus, the `IOStream` to `.out` files must be `close`d. Writing in `.temp`
  # is handled within julia, hence `.temp` are left open, ready to be
  # written.
  temp_A = open(temp_A_name, "w")
  temp_Q = open(temp_Q_name, "w")
  temp_u = open(temp_u_name, "w")
  temp_c = open(temp_c_name, "w")
  temp_P = open(temp_P_name, "w")

  last_A = open(last_A_name, "w")
  last_Q = open(last_Q_name, "w")
  last_u = open(last_u_name, "w")
  last_c = open(last_c_name, "w")
  last_P = open(last_P_name, "w")

  out_A = open(out_A_name, "w")
  out_Q = open(out_Q_name, "w")
  out_u = open(out_u_name, "w")
  out_c = open(out_c_name, "w")
  out_P = open(out_P_name, "w")

  node2 = convert(Int, floor(M*0.25))
  node3 = convert(Int, floor(M*0.5))
  node4 = convert(Int, floor(M*0.75))

  close(last_A)
  close(last_Q)
  close(last_u)
  close(last_c)
  close(last_P)

  close(out_A)
  close(out_Q)
  close(out_u)
  close(out_c)
  close(out_P)
  # The outlet boundary condition is specified by the last columns in the
  # `project.csv` file.
  #
  # 9 columns means that only the reflection coefficient `Rt` has been
  # specified. Since the three elements windkessel parameters must be
  # specified anyway, they are set to zero.
  if length(m) == 11

    Rt = m[11]

    R1 = 0.
    R2 = 0.
    Cc = 0.
  # 10 columns means that the three element windkessel is specified with
  # Reymond (2009) notation. Peripheral resistance (`R2`) is supplied along
  # with the outlet impedance (`R1`). To decouple these two parameters, `R1`
  # is computed as the vessel characteristic impedance
  # $$
  #   R_1 = Z_c = \rho \frac{c}{A_0}.
  # $$
  elseif length(m) == 12

    R1 = b.rho*waveSpeed(A0[end], gamma[end])/A0[end]
    #R1 = 1e5 #b.rho * 13.3/(A0*(2*R0*1e3)^0.3)
    R2 = m[11] - R1
    Cc = m[12]

    Rt = 0.
  # 11 columns means that the three parameters are specified.
  elseif length(m) >= 13

    R1 = m[11]
    R2 = m[12]
    Cc = m[13]

    Rt = 0.

  end

#     elseif length(m) == 10

#     R1 = b.rho*waveSpeed(A0[end], gamma[end])/A0[end]
#     #R1 = 1e5 #b.rho * 13.3/(A0*(2*R0*1e3)^0.3)
#     R2 = m[ 9] - R1
#     Cc = m[10]

#     Rt = 0.
#   # 11 columns means that the three parameters are specified.
#   elseif length(m) >= 11

#     R1 = b.rho*waveSpeed(A0[end], gamma[end])/A0[end] #m[ 9]
#     R2 = m[ 9] +m[10]#- R1 #m[10]
#     Cc = m[11]

#     Rt = 0.

#   end

  # `Pcn` is the pressure through the peripheral compliance of the three
  # elements windkessel. It is set to zero to simulate the pressure at the
  # artery-vein interface.
  Pcn = 0.

  #Slope
  slope = zeros(Float64, M)

  # MUSCL arrays
  flux  = zeros(Float64, 2, M+2)
  uStar = zeros(Float64, 2, M+2)

  vA = zeros(Float64, M+2)
  vQ = zeros(Float64, M+2)

  dU = zeros(Float64, 2, M+2)

  slopesA = zeros(Float64, M+2)
  slopesQ = zeros(Float64, M+2)

  Al = zeros(Float64, M+2)
  Ar = zeros(Float64, M+2)

  Ql = zeros(Float64, M+2)
  Qr = zeros(Float64, M+2)

  Fl = zeros(Float64, 2, M+2)
  Fr = zeros(Float64, 2, M+2)
#
# ----------------------------------------------------------------------------
# Returns:
# -------------- -------------------------------------------------------------
# `vessel_data`  [`::Vessel`](BTypes.html#Vessel) data structure
#                containing vessel mechanical and geometrical properties.
# ----------------------------------------------------------------------------
  vessel_data = BTypes.Vessel(vessel_name,
                    ID, sn, tn, rn,
                    M,
                    dx, invDx, halfDx,
                    Ccfl,
                    beta, gamma, gamma_ghost,
#                     R0,
                    A0, dA0dx, dTaudx, Pext,
                    A, Q,
                    u, c, P,
                    W1M0, W2M0,
                    U00A, U00Q, U01A, U01Q,
                    UM1A, UM1Q, UM2A, UM2Q,
                    temp_P_name, temp_Q_name, temp_A_name,
                    temp_c_name, temp_u_name,
                    last_P_name, last_Q_name, last_A_name,
                    last_c_name, last_u_name,
                    out_P_name, out_Q_name, out_A_name,
                    out_c_name, out_u_name,
                    temp_P, temp_Q, temp_A, temp_c, temp_u,
                    last_P, last_Q, last_A, last_c, last_u,
                    node2, node3, node4,
                    Rt,
                    R1, R2, Cc,
                    Pcn,
                    slope,
                    flux, uStar, vA, vQ,
                    dU, slopesA, slopesQ,
                    Al, Ar, Ql, Qr, Fl, Fr)

  return vessel_data
end

# *function* __`loadInputData`__ $\rightarrow$ `::Array{Float, 2}`
#
# ----------------------------------------------------------------------------
# Parameters:
# ------------------- --------------------------------------------------------
# `project_name`      `project` name string.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning
# ----------------------------------------------------------------------------
# `project_inlet.dat` is parsed by `readdlm` function and its content is
# saved in a `in_data` matrix. Its first column contains time and the second
# column contains flow waveform data.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Returns:
# --------- ------------------------------------------------------------------
# `in_data` `::Array{Float, 2}` `[time, flow]` matrix for inlet boundary
# condition.
# ----------------------------------------------------------------------------
#
function loadInputData(project_name)
  in_data = readdlm(join([project_name, "_inlet.dat"]))
  return in_data
end

# *function* __`loadGlobalConstants`__ $\rightarrow$ `::Heart`,
# `::Blood`, `::Float`
#
# ----------------------------------------------------------------------------
# Parameters:
# --------------------- ------------------------------------------------------
# `project_name`        `project` name string.
#
# `inlet_BC_switch`     `::Int` inlet boundary condtion selector:
#
#                        1 $\rightarrow$ heaviside sine function
#
#                        2 $\rightarrow$ gaussian
#
#                        3 $\rightarrow$ from `project_inlet.dat`
#
# `cycles`              `::Int` number of cardiac cycles to be simulated.
#
# `rho`                 `::Float` blood density
#
# `mu`                  `::Float` blood dynamic viscosity
# ----------------------------------------------------------------------------
#
# --------------------------------------------------------------------------
# Functioning
# --------------------------------------------------------------------------
# This function load into [heart](BTypes.html#Heart) and
# [blood](BTypes.html#Blood) datastructures the parameters constant
# within all vessels.
# --------------------------------------------------------------------------
#
function loadGlobalConstants(project_name,
                             inlet_BC_switch :: Int64,
                             inlet_type :: String,
                             cycles :: Int64,
                             rho :: Float64, mu:: Float64, gamma_profile :: Int64)
  # `heart` data structure is filled depending on the inlet boundary condition
  # chosen by the user. When the inlet flow function is given
  # (`inlet_BC_swithc`=3), the flow waveform is imported by
  # [`loadInputData`](initialise.html#loadInputData) function. The
  # `cardiac_period` is taken as the last row in the time column of the input
  # file. `sys_T`, `initial_flow`, and `flow_amplitude` quantities will not
  # be used and they are set to zero.
  if inlet_BC_switch == 3
    input_data = loadInputData(project_name)

#     if input_data[1,2] < 1e-6
#       error("Inlet flow initial value < 1e-6")
#     end

    cardiac_period = input_data[end, 1]
    sys_T = 0.
    initial_flow = 0.
    flow_amplitude = 0.
  # In any other case, all the parameters are take from `project_constants.jl`
  # file. `input_data` array is initialised to zero and it will not be used.
  else
    include(join([project_name, "_constants.jl"]))
    input_data = zeros(Float64, 1,1)
  end
  # The total simulation time `total_time` is calculated as number of cycles
  # times the cardiac period in seconds.
  const total_time = cycles*cardiac_period

  #Initialise heart data structure
  heart_data = BTypes.Heart(inlet_BC_switch,
                            inlet_type,
                            cardiac_period, sys_T,
                            initial_flow, flow_amplitude,
                            input_data)
  # Kinematic viscosity `nu` and viscous resitance term `Cf` are calculated as
  # $$
  #   \nu = \frac{\mu}{\rho}, \quad C_f = 8\pi\nu.
  # $$
  # where `Cf` is computed by assuming a parabolic velocity profile
  # along the radial direction.
  const nu = mu/rho
  const Cf = 8*pi*nu

  #Initialise blood data structure
  blood_data = BTypes.Blood(mu, rho, Cf, gamma_profile)
  # --------------------------------------------------------------------------
  # Returns:
  # -------------- -----------------------------------------------------------
  # `heart_data`   `::Heart`
  #
  # `blood_data`   `::Blood`
  #
  # `total_time`   `::Float`
  # --------------------------------------------------------------------------
  return heart_data, blood_data, total_time
end

ext: jl


text:
import DataFrames: DataFrame, DataFrameRow, nrow, ncol, eltypes, mapcols
import Random: shuffle
import Statistics: mean, std
import CUDA: CuVector, CuMatrix

export num_examples, num_features,
    example, feature_values,
    iscomplete, isweighted, isfpdata, isbinarydata,
    num_chunks, chunks, eltype,
    add_sample_weights, split_sample_weights, get_weights,
    shuffle_examples, batch, batch_size, isbatched,
    threshold, soften, marginal_prob, bagging_dataset,
    to_gpu, to_cpu, isgpu, same_device,
    ll_per_example, bits_per_pixel

# Basic data-related utilities.
# These methods assume the following possible representations.
# 1/ Complete data
#  - `Bool` data:
#      * `DataFrame` of `(Cu)BitVector`
#  - `Float` data:
#      * `DataFrame` of `(Cu)Vector{Float}`
# 2/ Incomplete data:
#  - `Bool` data:
#      * `DataFrame` of `Vector{Union{Bool,Missing}}`
#      * `DataFrame` of `CuVector{Int8}` where typemax represents missing
#  - `Float` data:
#      * `DataFrame` of `Vector{Union{Float,Missing}}`
#      * `DataFrame` of `CuVector{Float}` where typmax represents missing

"""
    num_examples(df::DataFrame)

Number of examples in data
"""
num_examples(df::DataFrame) = nrow(df)
num_examples(df::Array{DataFrame}) = mapreduce(d -> num_examples(d), +, df)

"""
    num_features(df::DataFrame)

Number of features in the data
"""
num_features(df::DataFrame) = isweighted(df) ? ncol(df) - 1 : ncol(df)
num_features(df::Array{DataFrame}) = num_features(df[1])

"Get the ith example"
example(d,i) = vec(convert(Array, d[i,:]))

"Get the ith feature values"
feature_values(df::DataFrame, i) = df[!,i]

import Base: eltype # extend

"Find a type that can capture all feature values"
eltype(df::DataFrame) = reduce(typejoin, eltype.(eachcol(df)))

"Is the data complete (no missing values)?"
iscomplete(data::DataFrame) = begin
    if isweighted(data)
        all(iscomplete, eachcol(data)[1:end-1])
    else
        all(iscomplete, eachcol(data))
    end
end
iscomplete(data::DataFrameRow) = begin
    if isweighted(data)
        all(v -> (v isa Bool), data[1:end-1])
    else
        all(v -> (v isa Bool), data)
    end
end
iscomplete(::AbstractArray{Bool}) = true
iscomplete(::Array{<:AbstractFloat}) = true
iscomplete(x::Array{Union{Bool,Missing}}) = !any(ismissing, x)
iscomplete(x::Array{Union{<:AbstractFloat,Missing}}) = !any(ismissing, x)
iscomplete(x::Union{CuArray{<:Int8},CuArray{<:AbstractFloat}}) =
    !any(v -> v == typemax(v), x)
iscomplete(data::Array{DataFrame}) = all(d -> iscomplete(d), data)

"Is the dataset binary?"
isbinarydata(df::DataFrame) = begin
    if isweighted(df)
        all(t -> nonmissingtype(t) <: Union{Bool,UInt8}, eltype.(eachcol(df))[1:end-1])
    else
        all(t -> nonmissingtype(t) <: Union{Bool,UInt8}, eltype.(eachcol(df)))
    end
end
isbinarydata(df::Array{DataFrame}) =
    all(d -> isbinarydata(d), df)

"Is the dataset consisting of floating point data?"
isfpdata(df::DataFrame) =
    all(t -> nonmissingtype(t) <: AbstractFloat, eltype.(eachcol(df)))
isfpdata(df::Array{DataFrame}) =
    all(d -> isfpdata(d), df)

"For binary complete data, how many `UInt64` bit strings are needed to store one feature?"
num_chunks(d::DataFrame) = num_chunks(feature_values(d,1))
num_chunks(v::AbstractBitVector) = length(v.chunks)

"Retrieve chunks of bit vector"
chunks(v::AbstractBitVector) = v.chunks
chunks(df::DataFrame, i) = chunks(feature_values(df, i))

# WEIGHTED DATASETS

"Add sample weights by adding a named column `weight`
 at the end of the DataFrame."
add_sample_weights(df::DataFrame, weights::DataFrame) = begin
    if isweighted(df)
        df = split_sample_weights(df)[1]
    end
    df.weight = weights[!,1]

    df
end
add_sample_weights(df::DataFrame, weights::AbstractArray) = begin
    if isweighted(df)
        df = split_sample_weights(df)[1]
    end
    df.weight = collect(Iterators.flatten(weights))

    df
end
add_sample_weights(df::Array{DataFrame}, weights) = begin
    if (weights isa DataFrame) || (weights isa AbstractArray{F} where F <: AbstractFloat)
        weights = batch(weights, batch_size(df))
    end

    map(zip(df, weights)) do (d, dw)
        add_sample_weights(d, dw)
    end
end

"Split a weighted dataset into unweighted dataset
 and its corresponding weights."
split_sample_weights(df::DataFrame) = begin
    @assert isweighted(df) "`df` is not weighted."
    weights = df[!, end]
    df = df[!, 1:end-1]

    df, weights
end
split_sample_weights(df::Array{DataFrame}) = begin
    @assert isweighted(df) "`df` is not weighted."
    weights = map(df) do d
        d[!, end]
    end

    df = map(df) do d
        d[!, 1:end-1]
    end

    df, weights
end

"Get the weights from a weighted dataset."
get_weights(df::DataFrame) = begin
    @assert isweighted(df) "`df` is not weighted."
    df[!, end]
end
get_weights(df::Array{DataFrame}) = begin
    @assert isweighted(df) "`df` is not weighted."
    map(df) do d
        d[!, end]
    end
end

"Is the dataset weighted?"
isweighted(df::Union{DataFrame, DataFrameRow}) =
    cmp(names(df)[end], "weight") === 0
isweighted(df::Array{DataFrame}) =
    all(d -> isweighted(d), df)

# DATA TRANSFORMATIONS

"""
    shuffle_examples(df::DataFrame)

Shuffle the examples in the data
"""
shuffle_examples(data) = data[shuffle(axes(data, 1)), :]

"Create mini-batches"
function batch(data, batchsize=1024; shuffle::Bool = true)
    to_gpu_flag = false
    if isgpu(data)
        to_gpu_flag = true
        data = to_cpu(data)
    end

    if data isa Array{T} where T <: AbstractFloat
        n_examples = length(data)
    else
        n_examples = num_examples(data)
    end

    if shuffle
        data = shuffle_examples(data)
    end
    data = map(1:batchsize:n_examples) do start_index
        stop_index = min(start_index + batchsize - 1, n_examples)
        data[start_index:stop_index, :]
    end

    if to_gpu_flag
        data = to_cpu(data)
    end

    data
end

"Is the dataset batched?"
isbatched(data::DataFrame) = false
isbatched(data::Array{DataFrame}) = true

"Batch size of the dataset"
batch_size(data::Array{DataFrame}) = num_examples(data[1])

"Dataset bagging"
function bagging_dataset(data::DataFrame; num_bags::Integer = 1, frac_examples::AbstractFloat = 1.0,
                         batch_size::Integer = 0)
    if isweighted(data)
        data, weights = split_sample_weights(data)
    else
        weights = num_examples(data)
    end

    # Randomly draw samples with replacement
    function random_sample(weights::AbstractArray, num_examples::Integer)
        weights = cumsum(weights)
        cum_weights = weights[end]
        len = length(weights)
        map(1:num_examples) do idx
            randval = rand() * cum_weights
            i, j, k = 1, len, 0
            while i < j - 1
                k = (i + j) ÷ 2
                if randval < weights[k]
                    j = k
                else
                    i = k
                end
            end
            i == j ? i : (randval >= weights[i] ? j : i)
        end
    end
    function random_sample(weights::Integer, num_examples::Integer)
        map(1:num_examples) do idx
            rand(1:weights)
        end
    end

    # Returns an array of DataFrames where each DataFrame is randomly sampled from the
    # original dataset `data`.
    map(1:num_bags) do dataset_idx
        example_idxs = random_sample(weights, convert(UInt32, floor(num_examples(data) * frac_examples)))
        batch_size == 0 ? data[example_idxs, :] : batch(data[example_idxs, :], batch_size)
    end
end

"Threshold a numeric dataset making it binary"
threshold(train, valid, test) = threshold(train, valid, test, 0.05) # default threshold offset (used for MNIST)

function threshold(train::DataFrame, valid, test, offset)
    @assert isfpdata(train) "DataFrame to be thresholded contains non-numeric columns: $(eltypes(x))"
    train = convert(Matrix, train)
    valid = issomething(valid) ? convert(Matrix, valid) : nothing
    test = issomething(test) ? convert(Matrix, test) : nothing
    means = mean(train, dims=1)
    stds = std(train, dims=1)
    threshold_value = means .+ offset .* stds
    train = DataFrame(BitArray(train .> threshold_value))
    valid = issomething(valid) ? DataFrame(BitArray(valid .> threshold_value)) : nothing
    test = issomething(test) ? DataFrame(BitArray(test .> threshold_value)) : nothing
    return train, valid, test
end

"Turn binary data into floating point data close to 0 and 1."
soften(data::DataFrame, softness=0.05; scale_by_marginal=true, precision=Float32) = begin
    n_col = ncol(data)
    data_weighted = isweighted(data)

    marginals = marginal_prob(data; precision = precision)
    vs = AbstractVector[]
    col_idx = 1
    for v in eachcol(data)
        if data_weighted && col_idx == n_col
            fv = v
        elseif scale_by_marginal
            fv = v .* precision(1.0 - softness) .+ (precision(softness) * marginals[col_idx])
        else
            fv = v .* precision(1.0 - 2.0 * softness) .+ precision(softness)
        end
        push!(vs, fv === v ? copy(fv) : fv)

        col_idx += 1
    end
    DataFrame(vs, names(data), copycols = false)
end

"Compute the marginal prob of each feature in a binary dataset."
marginal_prob(data; precision=Float32) = begin
    @assert isbinarydata(data) "marginal_prob only support binary data."
    n_examples = num_examples(data)
    map(1:num_features(data)) do idx
        precision(sum(feature_values(data, idx))) / n_examples
    end
end

@inline _msk_end(df::DataFrame) = _msk_end(num_examples(df))

"Move data to the GPU"
to_gpu(d::Union{CuBitVector,CuArray}) = d
to_gpu(m::Array) = CuArray(m)
to_gpu(v::BitVector) =
    AbstractBitVector(to_gpu(chunks(v)), length(v))
to_gpu(v::Vector{Union{F,Missing}}) where F<:AbstractFloat =
    CuArray(T.(coalesce(v,typemax(T))))
to_gpu(v::Vector{Union{Bool,Missing}}) =
    CuArray(UInt8.(coalesce.(v,typemax(UInt8))))
to_gpu(df::DataFrame) = isgpu(df) ? df : mapcols(to_gpu, df)
to_gpu(df::Array{DataFrame}) = map(df) do d
    to_gpu(d)
end

"Move data to the CPU"
to_cpu(m::Union{BitVector,Array}) = m
to_cpu(m::CuArray) = Array(m)
to_cpu(cv::CuBitVector) =
    AbstractBitVector(to_cpu(cv.chunks), length(cv))
to_cpu(v::CuVector{T}) where T<:AbstractFloat =
    replace(Array(v), typemax(T) => missing)
to_cpu(v::CuVector{UInt8}) =
    convert(Vector{Union{Bool,Missing}},
        replace(Array(v), typemax(UInt8) => missing))
to_cpu(df::DataFrame) = mapcols(to_cpu, df)
to_cpu(df::Array{DataFrame}) = map(df) do d
    to_gpu(d)
end

"Check whether data resides on the GPU"
isgpu(::Union{Array, BitArray}) = false
isgpu(::Union{CuArray, CuBitVector}) = true
isgpu(df::DataFrame) = all(isgpu, eachcol(df))
isgpu(df::Array{DataFrame}) = all(isgpu, df)

"Ensure that `x` resides on the same device as `data`"
same_device(x, data) = isgpu(data) ? to_gpu(x) : to_cpu(x)

# LIKELIHOOD HELPERS

"Normalize the given log-likelihood by the number of examples in `data`"
ll_per_example(ll, data) = ll / num_examples(data)

"Normalize the given log-likelihood as bits per pixel in `data`"
bits_per_pixel(ll, data) = -(ll_per_example(ll, data)  / num_features(data)) / log(2)

"Computer the per-example log-likelihood of a fully factorized ML model on Bool data"
function fully_factorized_log_likelihood(data; pseudocount=0)
    @assert isbinarydata(data) && iscomplete(data) "This method requires complete binary data"
    m = convert(Matrix,data)
    counts = sum(m, dims=1)
    smoothed_counts = counts .+ pseudocount/2.0
    log_estimates = log.(smoothed_counts ./ (num_examples(data) + pseudocount))
    ll = sum(log_estimates .* counts)
    ll += sum(log.(1 .- exp.(log_estimates)) .* (num_examples(data) .- counts))
    ll_per_example(ll, data)
end


text:
const fdc1_2 = [0, 1/2]
const fdc1_4 = [0, 2/3, -1/12]
const fdc1_6 = [0, 3/4, -3/20,1/60 ]
const fdc1_8 = [0, 4/5, -1/5 ,4/105, -1/280]
const fdc2_2 = [-2, 1]
const fdc2_4 = [-5/2,    4/3,-1/12]
const fdc2_6 = [-49/18 , 3/2,-3/20, 1/90]
const fdc2_8 = [-205/72, 8/5,-1/5,  8/315, -1/560]
const fdc1 = [fdc1_2, fdc1_4, fdc1_6, fdc1_8]
const fdc2 = [fdc2_2, fdc2_4, fdc2_6, fdc2_8]
fdcoeff1 = Vector{OffsetVector{Float64, Vector{Float64}}}(undef, 4)
fdcoeff2 = Vector{OffsetVector{Float64, Vector{Float64}}}(undef, 4)
function init!(fdcoeff1, fdcoeff2)
	 for i in 1:4
		tmp1 = zeros(2*i + 1)
		tmp1[i + 1] = fdc1[i][1]
		tmp2 = zeros(2*i + 1)
		tmp2[i + 1] = fdc2[i][1]
		for k in 1:i
			tmp1[i + 1 + k] =  fdc1[i][k+1]
			tmp1[i + 1 - k] = -fdc1[i][k+1]
			tmp2[i + 1 + k] =  fdc2[i][k+1]
			tmp2[i + 1 - k] =  fdc2[i][k+1]
		end
		fdcoeff1[i] = OffsetVector(tmp1, -i:i)
		fdcoeff2[i] = OffsetVector(tmp2, -i:i)
	end
end
init!(fdcoeff1, fdcoeff2)
cycindex(i,n) = (i + n - 1)%n + 1
function dfdt(f::TVF, dh::T; order = 2, bc = :Zero)where {T, TVF<:AbstractVector{T}}
	@assert order in Set([2,4,6,8])
	@assert bc in Set([:Zero, :Periodic])
	n = length(f)
	m = div(order, 2)
	g = similar(f)
	if bc == :Zero
		for i in 1:n
			for k in 1:m
				g[i] = fdcoeff1[m][0]*f[i]
				if i + k ≤ n
					g[i] += fdcoeff1[m][ k]*f[i+k]
				end
				if 1 ≤ i - k
					g[i] += fdcoeff1[m][-k]*f[i-k]
				end
			end
		end
	else
		for i in 1:n
			g[i] = fdcoeff1[m][0]*f[i]
			for k in 1:m
				g[i] += fdcoeff1[m][ k]*f[cycindex(i + k, n)] + fdcoeff1[m][-k]*f[cycindex(i - k, n)]
			end
		end
	end
	g ./= dh
	return g
end
function d2fdt(f::TVF, dh::T; order = 2, bc = :Zero)where {T, TVF<:AbstractVector{T}}
	@assert order in Set([2,4,6,8])
	@assert bc in Set([:Zero, :Periodic])
	n = length(f)
	m = div(order, 2)
	g = similar(f)
	if bc == :Zero
		for i in 1:n
			for k in 1:m
				g[i] = fdcoeff2[m][0]*f[i]
				if i + k ≤ n
					g[i] += fdcoeff2[m][ k]*f[i+k]
				end
				if 1 ≤ i - k
					g[i] += fdcoeff2[m][-k]*f[i-k]
				end
			end
		end
	else
		for i in 1:n
			g[i] = fdcoeff2[m][0]*f[i]
			for k in 1:m
				g[i] += fdcoeff2[m][ k]*f[cycindex(i + k, n)] + fdcoeff2[m][-k]*f[cycindex(i - k, n)]
			end
		end
	end
	g ./= dh^2
	return g
end
function makefd1mat(n::Integer, dh, order = 2; bc = :Zero)
	@assert order in Set([2,4,6,8])
	@assert bc in Set([:Zero, :Periodic])
	mat = zeros(typeof(float(dh)), (n, n))
	m = div(order, 2)
	if bc == :Zero
		for i in 1:n
			mat[i,i] = fdcoeff1[m][0]
			for k in 1:m
				if i + k ≤ n
					mat[i,i + k] = fdcoeff1[m][ k]
				end
				if 1 ≤ i - k
					mat[i,i - k] = fdcoeff1[m][-k]
				end
			end
		end
	else
		for i in 1:n
			mat[i,i] = fdcoeff1[m][0]
			for k in 1:m
				mat[i,cycindex(i + k, n)] = fdcoeff1[m][ k]
				mat[i,cycindex(i - k, n)] = fdcoeff1[m][-k]
			end
		end
	end
	mat ./= dh
	return mat
end
function makefd2mat(n::Integer, dh, order = 2; bc = :Zero)
	@assert order in Set([2,4,6,8])
	@assert bc in Set([:Zero, :Periodic])
	mat = zeros(typeof(float(dh)), (n, n))
	m = div(order, 2)
	if bc == :Zero
		for i in 1:n
			mat[i,i] = fdcoeff2[m][0]
			for k in 1:m
				if i + k ≤ n
					mat[i,i + k] = fdcoeff2[m][ k]
				end
				if 1 ≤ i - k
					mat[i,i - k] = fdcoeff2[m][-k]
				end
			end
		end
	else
		for i in 1:n
			mat[i,i] = fdcoeff2[m][0]
			for k in 1:m
				mat[i,cycindex(i + k, n)] = fdcoeff2[m][ k]
				mat[i,cycindex(i - k, n)] = fdcoeff2[m][-k]
			end
		end
	end
	mat ./= dh^2
	return mat
end


text:
module Thrift

using Distributed

using Sockets
import Sockets: TCPServer, listen, accept

import Base: open, close, isopen, read, read!, write, flush, skip, show, copy!

export open, close, isopen, read, read!, write, flush, skip, listen, accept, show, copy!

# from base.jl
export TSTOP, TVOID, TBOOL, TBYTE, TI08, TDOUBLE, TI16, TI32, TI64, TSTRING, TUTF7, TSTRUCT, TMAP, TSET, TLIST, TUTF8, TUTF16
export TType, TProcessor, TTransport, TServerTransport, TServer, TProtocol
export writeMessageBegin, writeMessageEnd, writeStructBegin, writeStructEnd, writeFieldBegin, writeFieldEnd, writeFieldStop, writeMapBegin, writeMapEnd, writeListBegin, writeListEnd, writeSetBegin, writeSetEnd, writeBool, writeByte, writeI16, writeI32, writeI64, writeDouble, writeString
export readMessageBegin, readMessageEnd, readStructBegin, readStructEnd, readFieldBegin, readFieldEnd, readMapBegin, readMapEnd, readListBegin, readListEnd, readSetBegin, readSetEnd, readBool, readByte, readI16, readI32, readI64, readDouble, readString
export ApplicationExceptionType, MessageType, TException, TApplicationException
export ThriftMetaAttribs, ThriftMeta, meta
export isinitialized, set_field, set_field!, get_field, clear, has_field, fillunset, fillset, filled, isfilled, thriftbuild, enumstr

# from transports.jl
export TFramedTransport, TSASLClientTransport, TSocket, TServerSocket, TSocketBase, TMemoryTransport, TFileTransport
export TransportExceptionTypes, TTransportException

# from sasl.jl
export SASL_MECH_PLAIN, SASL_MECH_KERB, SASL_MECH_LDAP, SASLException

# from protocols.jl
export TBinaryProtocol, TCompactProtocol

# from processor.jl
export ThriftProcessor, ThriftHandler, process, handle, extend, distribute

# from server.jl
export TSimpleServer, TTaskServer, TProcessPoolServer, serve

# enable logging only during debugging
#using Logging
#const logger = Logging.configure(level=DEBUG)
##const logger = Logging.configure(filename="/tmp/thrift$(getpid()).log", level=DEBUG)
#macro logmsg(s)
#    quote
#        debug($(esc(s)))
#    end
#end
macro logmsg(s)
end

include("base.jl")
include("codec.jl")
include("sasl.jl")
include("transports.jl")
include("protocols.jl")
include("processor.jl")
include("server.jl")

end # module