text:
# This file is part of the ValidatedNumerics.jl package; MIT licensed

using ValidatedNumerics
using FactCheck

setprecision(Interval, Float64)

a = @interval(0.1, 1.1)
b = @interval(0.9, 2.0)
c = @interval(0.25, 4.0)

facts("Consistency tests") do

    @fact isa( @interval(1,2), Interval ) --> true
    @fact isa( @interval(0.1), Interval ) --> true
    @fact isa( zero(b), Interval ) --> true

    @fact zero(b) --> 0.0
    @fact zero(b) == zero(typeof(b)) --> true
    @fact one(a) --> 1.0
    @fact one(a) == one(typeof(a)) --> true
    @fact one(a) --> big(1.0)
    @fact a == b --> false
    @fact a != b --> true

    @fact a --> Interval(a.lo, a.hi)
    @fact @interval(1, Inf) --> Interval(1.0, Inf)
    @fact @interval(-Inf, 1) --> Interval(-Inf, 1.0)
    @fact @biginterval(1, Inf) --> Interval{BigFloat}(1.0, Inf)
    @fact @biginterval(-Inf, 1) --> Interval{BigFloat}(-Inf, 1.0)
    @fact @interval(-Inf, Inf) --> entireinterval(Float64)
    @fact emptyinterval(Rational{Int}) --> ∅

    @fact 1 == zero(a)+one(b) --> true
    @fact Interval(0,1) + emptyinterval(a) --> emptyinterval(a)
    @fact @interval(0.25) - one(c)/4 --> zero(c)
    @fact emptyinterval(a) - Interval(0,1) --> emptyinterval(a)
    @fact Interval(0,1) - emptyinterval(a) --> emptyinterval(a)
    @fact a*b --> Interval(a.lo*b.lo, a.hi*b.hi)
    @fact Interval(0,1) * emptyinterval(a) --> emptyinterval(a)
    @fact a * Interval(0) --> zero(a)
end

facts("inv") do

    @fact inv( zero(a) ) --> emptyinterval()
    @fact inv( @interval(0, 1) ) --> Interval(1, Inf)
    @fact inv( @interval(1, Inf) ) --> Interval(0, 1)
    @fact inv(c) --> c
    @fact one(b)/b --> inv(b)
    @fact a/emptyinterval(a) --> emptyinterval(a)
    @fact emptyinterval(a)/a --> emptyinterval(a)
    @fact inv(@interval(-4.0,0.0)) --> @interval(-Inf, -0.25)
    @fact inv(@interval(0.0,4.0)) --> @interval(0.25, Inf)
    @fact inv(@interval(-4.0,4.0)) --> entireinterval(Float64)
    @fact @interval(0)/@interval(0) --> emptyinterval()
    @fact typeof(emptyinterval()) --> Interval{Float64}
end

facts("fma consistency") do
    @fact fma(emptyinterval(), a, b) --> emptyinterval()
    @fact fma(entireinterval(), zero(a), b) --> b
    @fact fma(entireinterval(), one(a), b) --> entireinterval()
    @fact fma(zero(a), entireinterval(), b) --> b
    @fact fma(one(a), entireinterval(), b) --> entireinterval()
    @fact fma(a, zero(a), c) --> c
    @fact fma(Interval(1//2), Interval(1//3), Interval(1//12)) --> Interval(3//12)
end

facts("∈ tests") do
    @fact Inf ∈ entireinterval() --> false
    @fact 0.1 ∈ @interval(0.1) --> true
    @fact 0.1 in @interval(0.1) --> true
    @fact -Inf ∈ entireinterval() --> false
    @fact Inf ∈ entireinterval() --> false
end

facts("Inclusion tests") do
    @fact b ⊆ c --> true
    @fact emptyinterval(c) ⊆ c --> true
    @fact c ⊆ emptyinterval(c) --> false
    @fact interior(b,c) --> true
    @fact b ⪽ emptyinterval(b) --> false
    @fact emptyinterval(c) ⪽ c --> true
    @fact emptyinterval(c) ⪽ emptyinterval(c) --> true
    @fact isdisjoint(a, @interval(2.1)) --> true
    @fact isdisjoint(a, b) --> false
    @fact isdisjoint(emptyinterval(a), a) --> true
    @fact isdisjoint(emptyinterval(), emptyinterval()) --> true

end

facts("Comparison tests") do
    @fact ValidatedNumerics.islessprime(a.lo, b.lo) --> a.lo < b.lo
    @fact ValidatedNumerics.islessprime(Inf, Inf) --> true
    @fact ∅ <= ∅ --> true
    @fact Interval(1.0,2.0) <= ∅ --> false
    @fact Interval(-Inf,Inf) <= Interval(-Inf,Inf) --> true
    @fact Interval(-0.0,2.0) ≤ Interval(-Inf,Inf) --> false
    @fact precedes(∅,∅) --> true
    @fact precedes(Interval(3.0,4.0),∅) --> true
    @fact precedes(Interval(0.0,2.0),Interval(-Inf,Inf)) --> false
    @fact precedes(Interval(1.0,3.0),Interval(3.0,4.0)) --> true
    @fact strictprecedes(Interval(3.0,4.0),∅) --> true
    @fact strictprecedes(Interval(-3.0,-1.0),Interval(-1.0,0.0)) --> false
    @fact iscommon(emptyinterval()) --> false
    @fact iscommon(entireinterval()) --> false
    @fact iscommon(a) --> true
    @fact isunbounded(emptyinterval()) --> false
    @fact isunbounded(entireinterval()) --> true
    @fact isunbounded(Interval(-Inf, 0.0)) --> true
    @fact isunbounded(Interval(0.0, Inf)) --> true
    @fact isunbounded(a) --> false
end

facts("Intersection tests") do
    @fact emptyinterval() --> Interval(Inf, -Inf)
    @fact a ∩ @interval(-1) --> emptyinterval(a)
    @fact isempty(a ∩ @interval(-1) ) --> true
    @fact isempty(a) --> false
    @fact emptyinterval(a) == a --> false
    @fact emptyinterval() == emptyinterval() --> true

    @fact intersect(a, hull(a,b)) --> a
    @fact union(a,b) --> Interval(a.lo, b.hi)
end

facts("Special interval tests") do

    @fact entireinterval(Float64) --> Interval(-Inf, Inf)
    @fact isentire(entireinterval(a)) --> true
    @fact isentire(Interval(-Inf, Inf)) --> true
    @fact isentire(a) --> false
    @fact Interval(-Inf, Inf) ⪽ Interval(-Inf, Inf) --> true

    @fact nai(a) == nai(a) --> false
    @fact nai(a) === nai(a) --> true
    @fact nai(Float64) === Interval(NaN) --> true
    @fact isnan(nai(BigFloat).lo) --> true
    @fact isnai(nai()) --> true
    @fact isnai(a) --> false

    @fact infimum(a) == a.lo --> true
    @fact supremum(a) == a.hi --> true
    @fact infimum(emptyinterval(a)) --> Inf
    @fact supremum(emptyinterval(a)) --> -Inf
    @fact infimum(entireinterval(a)) --> -Inf
    @fact supremum(entireinterval(a)) --> Inf
    @fact isnan(supremum(nai(BigFloat))) --> true
end

facts("mid etc.") do

    @fact mid( Interval(1//2) ) --> 1//2
    @fact diam( Interval(1//2) ) --> 0//1
    @fact diam( @interval(1//10) ) --> eps(0.1)
    @fact diam( @interval(0.1) ) --> eps(0.1)
    @fact isnan(diam(emptyinterval())) --> true
    @fact mig(@interval(-2,2)) --> BigFloat(0.0)
    @fact mig( Interval(1//2) ) --> 1//2
    @fact isnan(mig(emptyinterval())) --> true
    @fact mag(-b) --> b.hi
    @fact mag( Interval(1//2) ) --> 1//2
    @fact isnan(mag(emptyinterval())) --> true
    @fact diam(a) --> 1.0000000000000002
end

facts("cancelplus tests") do

    x = Interval(-2.0, 4.440892098500622e-16)
    y = Interval(-4.440892098500624e-16, 2.0)
    @fact cancelminus(x, y) --> entireinterval(Float64)
    @fact cancelplus(x, y) --> entireinterval(Float64)
    x = Interval(-big(1.0), eps(big(1.0))/4)
    y = Interval(-eps(big(1.0))/2, big(1.0))
    @fact cancelminus(x, y) --> entireinterval(BigFloat)
    @fact cancelplus(x, y) --> entireinterval(BigFloat)
    x = Interval(-big(1.0), eps(big(1.0))/2)
    y = Interval(-eps(big(1.0))/2, big(1.0))
    @fact cancelminus(x, y) ⊆ Interval(-one(BigFloat), one(BigFloat)) --> true
    @fact cancelplus(x, y) --> Interval(zero(BigFloat), zero(BigFloat))
    @fact cancelminus(emptyinterval(), emptyinterval()) --> emptyinterval()
    @fact cancelplus(emptyinterval(), emptyinterval()) --> emptyinterval()
    @fact cancelminus(emptyinterval(), Interval(0.0, 5.0)) --> emptyinterval()
    @fact cancelplus(emptyinterval(), Interval(0.0, 5.0)) --> emptyinterval()
    @fact cancelminus(entireinterval(), Interval(0.0, 5.0)) --> entireinterval()
    @fact cancelplus(entireinterval(), Interval(0.0, 5.0)) --> entireinterval()
    @fact cancelminus(Interval(5.0), Interval(-Inf, 0.0)) --> entireinterval()
    @fact cancelplus(Interval(5.0), Interval(-Inf, 0.0)) --> entireinterval()
    @fact cancelminus(Interval(0.0, 5.0), emptyinterval()) --> entireinterval()
    @fact cancelplus(Interval(0.0, 5.0), emptyinterval()) --> entireinterval()
    @fact cancelminus(Interval(0.0), Interval(0.0, 1.0)) --> entireinterval()
    @fact cancelplus(Interval(0.0), Interval(0.0, 1.0)) --> entireinterval()
    @fact cancelminus(Interval(0.0), Interval(1.0)) --> Interval(-1.0)
    @fact cancelplus(Interval(0.0), Interval(1.0)) --> Interval(1.0)
    @fact cancelminus(Interval(-5.0, 0.0), Interval(0.0, 5.0)) --> Interval(-5.0)
    @fact cancelplus(Interval(-5.0, 0.0), Interval(0.0, 5.0)) --> Interval(0.0)
end

facts("mid and radius") do

    # NOTE: By some strange reason radius is not recognized here
    @fact ValidatedNumerics.radius(Interval(-1//10,1//10)) -->
        diam(Interval(-1//10,1//10))/2
    @fact isnan(ValidatedNumerics.radius(emptyinterval())) --> true
    @fact mid(c) == 2.125 --> true
    @fact isnan(mid(emptyinterval())) --> true
    @fact mid(entireinterval()) == 0.0 --> true
    @fact isnan(mid(nai())) --> true
    # In v0.3 it corresponds to AssertionError
    @fact_throws ArgumentError nai(Interval(1//2))
end

facts("abs, min, max, sign") do

    @fact abs(entireinterval()) --> Interval(0.0, Inf)
    @fact abs(emptyinterval()) --> emptyinterval()
    @fact abs(Interval(-3.0,1.0)) --> Interval(0.0, 3.0)
    @fact abs(Interval(-3.0,-1.0)) --> Interval(1.0, 3.0)
    @fact min(entireinterval(), Interval(3.0,4.0)) --> Interval(-Inf, 4.0)
    @fact min(emptyinterval(), Interval(3.0,4.0)) --> emptyinterval()
    @fact min(Interval(-3.0,1.0), Interval(3.0,4.0)) --> Interval(-3.0, 1.0)
    @fact min(Interval(-3.0,-1.0), Interval(3.0,4.0)) --> Interval(-3.0, -1.0)
    @fact max(entireinterval(), Interval(3.0,4.0)) --> Interval(3.0, Inf)
    @fact max(emptyinterval(), Interval(3.0,4.0)) --> emptyinterval()
    @fact max(Interval(-3.0,1.0), Interval(3.0,4.0)) --> Interval(3.0, 4.0)
    @fact max(Interval(-3.0,-1.0), Interval(3.0,4.0)) --> Interval(3.0, 4.0)
    @fact sign(entireinterval()) --> Interval(-1.0, 1.0)
    @fact sign(emptyinterval()) --> emptyinterval()
    @fact sign(Interval(-3.0,1.0)) --> Interval(-1.0, 1.0)
    @fact sign(Interval(-3.0,-1.0)) --> Interval(-1.0, -1.0)

    @fact log(@interval(-2,5)) --> @interval(-Inf,log(5.0))

    # Test putting functions in @interval:
    @fact @interval(sin(0.1) + cos(0.2)) --> sin(@interval(0.1)) + cos(@interval(0.2))

    f(x) = 2x
    @fact @interval(f(0.1)) --> f(@interval(0.1))

    # midpoint-radius representation
    a = @interval(0.1)
    midpoint, radius = midpoint_radius(a)

    @fact interval_from_midpoint_radius(midpoint, radius) -->
        Interval(0.09999999999999999, 0.10000000000000002)

end

facts("Precision tests") do
    setprecision(Interval, 100)
    @fact precision(Interval) == (BigFloat, 100) --> true

    setprecision(Interval, Float64)
    @fact precision(Interval) == (Float64, 100) --> true

    a = @interval(0.1, 0.3)

    b = setprecision(Interval, 64) do
        @interval(0.1, 0.3)
    end

    @fact b ⊆ a --> true

    @fact precision(Interval) == (Float64, 100) --> true

end

facts("Interval rounding tests") do
    setrounding(Interval, :wide)
    @fact rounding(Interval) == :wide --> true

    @fact_throws ArgumentError setrounding(Interval, :hello)

    setrounding(Interval, :narrow)
    @fact rounding(Interval) == :narrow --> true

end

facts("Interval power of an interval") do
    a = @interval(1, 2)
    b = @interval(3, 4)

    @fact a^b --> @interval(1, 16)
    @fact a^@interval(0.5, 1) --> a
    @fact a^@interval(0.3, 0.5) --> @interval(1, sqrt(2))

    @fact b^@interval(0.3) == Interval(1.3903891703159093, 1.5157165665103982) --> true

end


text:
using nevanlinna_ac
using Test

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

ext: jl


text:
using RecurrenceAnalysis, DelimitedFiles, Statistics

# Measure the times (in ms) of evaluating an expression n times
macro measuretime(ex, n)
    quote
        # Train the expression and get the result
        result = $(esc(ex))
        t = zeros($n)
        for i in 1:$n
            t[i] = 1000*(@elapsed $(esc(ex)))
        end
        t, result
    end
end

# Function that will be measured
function fun_rqa(x,metric)
    xe = embed(x,3,6)
    rmat = RecurrenceMatrix(xe,1.2,metric=metric)
    rqa(rmat,theiler=1)
end

# Analyse 12 series from 250 to 3000 points
# (With variable metric)
function benchmark(metric)
    m = readdlm("rossler.txt")
    for r=1:12
        x=m[1:250r,2r-1]
        tt, res = @measuretime fun_rqa(x,metric) 5
        t = median(tt)
        # Write table of results
        open("benchmark_rqa_julia_$metric.txt","a") do f
            write(f,"$r\t$t\t")
            for k in ["RR","DET","L","Lmax","ENTR","LAM","TT"]
                write(f, "$(res[k])\t")
            end
            write(f,"\n")
        end
    end
end
# Do it with max and euclidean norms
benchmark("euclidean")
benchmark("max")


text:
module PortableGameNotation

import Base.show, Base.length, Base.println

using Dates
using Printf

export readpgn, writepgn, Game, event, site, date, round, white, black, result,
  whiteelo, blackelo, eventdate, eco, movetext, plycount, length, movestring,
  headerstring, repr, intresult, whiteev, blackev, whitescore, blackscore,
  whiteperfelo, blackperfelo, isdecisive

struct Move
  san::String
  nag::UInt8
  comment::String
end
Move(san::String, comment::String) = Move(san, 0, comment)
struct Game
  header::Dict{String, String}
  moves::Array{Move,1}
end

RESULT_HASH = Dict{String,Int}("1-0" => 1, "1/2-1/2" => 0, "0-1" => -1, "*" => 0)
REQUIRED_TAGS = ["Event", "Site", "Date", "Round", "White", "Black", "Result"]
DEFAULT_HASH = Dict("Event"=>"","Site"=>"","Date"=>"","Round"=>"","White"=>"",
  "Black"=>"","Result"=>"")

const SIDE_WHITE = 0
const SIDE_BLACK = 1

const STATE_UNKNOWN = 0
const STATE_MOVE_NUMBER = 1  # currently reading a move number
const STATE_MOVE = 2         # currently reading a move
const STATE_COMMENT = 3      # currently reading a comment
const STATE_NAG = 4          # numeric annotation glyph
const STATE_RESULT = 5

move_scratch = Array{Move}(undef, 600)

"""
readpgn(filename; [skim=false, verbose=false])

Read games from PGN file `filename` and returns an array of `Game`
objects.

If `skim` is set to true, movetexts will be ignored, making the read
process faser and less memory intensive for cases when only header
information is required.
"""
function readpgn(pgnfilename::String; verbose=true, skim=false, prealloc=true)
  # This is an alternative constructor that performs a state-based
  # single scan of the move text to parse into Move objects, hopefully
  # as quickly and efficiently as possible
  # skim=true will skip move parsing and just return a Game object
  # with headers as quickly as possible for cases in which you don't
  # care about the moves
  # TODO:
  #   - deal with precomments
  #   - deal with NAGS
  #   - deal with ; comments
  #   - implement skim
  # COMMENTS
  # Comments are inserted by either a ; (a comment that continues
  # to the end of the line) or a { (which continues until a matching }).
  # Comments do not nest.
  f = open(pgnfilename,"r")
  lines = readlines(f)
  close(f)
  if prealloc
    games = Vector{Game}(undef, div(filesize(pgnfilename), 500)) # guess
  else
    games = Vector{Game}()
  end
  ngamesread = 0
  newgame = false
  current_header = DEFAULT_HASH
  current_movetext = ""
  for (j, l) in enumerate(lines)
    if length(l) == 0
      if length(current_movetext) > 0 && length(current_header) > 0 # finished a game
        newgame = true
      end
    elseif l[1] == '['
      k, v = parse_header_line(l)
      current_header[k] = v
    else
      current_movetext *= l
    end
    if newgame || j == length(lines)
      moves = parse_movetext(current_movetext)
      g = Game(current_header, moves)
      ngamesread += 1
      print("\r", ngamesread)
      if ngamesread <= length(games)
        games[ngamesread] = g
      else
        println("PUSHING")
        push!(games, g)
      end
      current_movetext = ""
      current_header = DEFAULT_HASH
      newgame = false
    end
  end
  return games[1:ngamesread]
end

function parse_header_line(line::String)::Tuple{String,String}
  # parse a PGN header line and return a (key,val) tuple
  k1 = findfirst("[", line)[1] + 1
  k2 = findnext(" ", line, k1)[1] - 1
  v1 = findnext("\"", line, k2)[1] + 1
  v2 = findnext("\"", line, v1)[1] - 1
  return (line[k1:k2], line[v1:v2])
end

function parse_movetext(txt::String)::Array{Move,1}
  # parse a line of movetext
  # assumes <= 100 moves per line
  nmoves = 0
  state = STATE_UNKNOWN
  current_move = ""
  current_nag = 0
  current_comment = ""
  k = 1
  #try
  while k < length(txt)
    if state == STATE_UNKNOWN
      if isdigit(txt[k])
        state = STATE_MOVE_NUMBER
      elseif txt[k] == '{'
        state = STATE_COMMENT
      elseif txt[k] == '$'
        state = STATE_NAG
      elseif txt[k] != ' '
        state = STATE_MOVE
      else  # skip everything else
        k += 1
      end
    elseif state == STATE_MOVE_NUMBER
      k += 1
      while true
        if isdigit(txt[k]) || txt[k] == '.'
          k += 1
        elseif txt[k] == '-' || txt[k] == '/'
          state = STATE_RESULT
          break
        else
          state = STATE_MOVE  # replace with ischar?
          break
        end
        if k > length(txt)
          state = STATE_RESULT
          break
        end
      end
    elseif state == STATE_COMMENT
      k += 1
      comment_depth = 1
      startidx = k
      while comment_depth > 0
        if txt[k] == '{'
          comment_depth += 1
        elseif txt[k] == '}'
          comment_depth -= 1
        end
        k += 1
      end
      current_comment = txt[startidx:k-2]
      state = STATE_UNKNOWN
    elseif state == STATE_MOVE
      if current_move != ""
        nmoves += 1
        move_scratch[nmoves] = Move(current_move, current_nag, current_comment)
        current_comment = ""
      end
      startidx = k
      while true
        if txt[k] == ' '
          current_move = txt[startidx:k-1]
          break
        elseif k >= length(txt)
          current_move = txt[startidx:k]
          break
        else
          k += 1
        end
      end
      state = STATE_UNKNOWN
    elseif state == STATE_NAG
      k += 1
      startidx = k
      while txt[k] != ' '
        k += 1
        current_nag = parse(Int, txt[startidx:k])
      end
    elseif state == STATE_RESULT
      break
    else
      @info "UNHANDLED STATE" state txt k moves
    end
  end
  #catch BoundsError
    #@info "BoundsError" length(txt) k txt
    #exit(1)
  #end
  if current_move != ""
    nmoves += 1
    move_scratch[nmoves] = Move(current_move, current_nag, current_comment)
  end
  return move_scratch[1:nmoves]
end

"""
headerstring(g)

Format the headers of game `g` according to the PGN specification.
"""
function headerstring(g::Game)
  s = String[]
  # tags required by standard must be printed first and in order
  for t in REQUIRED_TAGS
    push!(s, "[$t \"$(g.header[t])\"]\n")
  end
  for k in keys(g.header)
    if !(k in REQUIRED_TAGS)
      push!(s, "[$k \"$(g.header[k])\"]\n")
    end
  end
  join(s, "")
end

"""
movestring(g)

Format the move text of game `g` with wrapping and linebreaks.
"""
function movestring(g::Game; line=80, comments=true)
  output = String[]
  linelen = 0
  ply = 0
  for m in g.moves
    if ply % 2 == 0
      str ="$(div(ply,2)+1)."
      linelen += length(str)
      if linelen >= line
        push!(output, "\n")
        linelen = length(str)
      end
      push!(output, str)
    end
    linelen += length(m.san) + 1
    if linelen >= line
        push!(output, "\n")
        linelen = length(m.san) + 1
    end
    push!(output, "$(m.san) ") # TODO: avoid string interp?
    if comments && m.comment != ""
      linelen += length(m.comment) + 1
      if linelen >= line
        push!(output, "\n")
        linelen = length(m.comment) + 1
      else
        push!(output, " ")
      # TODO: wrap comments
      end
      push!(output, "{$(m.comment)} ")
    end
    ply += 1
  end
  join(output, "")
end

resultstring(g::Game) = g.header["Result"]

function Base.repr(mime, m::Move)
  if m.side == SIDE_WHITE && m.comment == ""
    return "$(m.number)."*m.san
  elseif m.side == SIDE_WHITE && m.comment != ""
    return "$(M.number)."*m.san*" "*m.comment
  elseif m.comment == ""
    return m.san
  else
    return m.san*" "*m.comment
  end
end
function Base.show(io::IO, g::Game)
  println(io, white(g), " - ", black(g), ", ", site(g), " ", Dates.year(date(g)))
end
function Base.println(g::Game)
    Base.println(headerstring(g))
    Base.print(movestring(g), " ")
    Base.println(resultstring(g), "\n")
end

plycount(g::Game) = length(g.moves)
"""
length(g)

Number of moves in the game `g`.
"""
length(g::Game) = div(plycount(g), 2) + 1

function validate(g::Game)
  for t in REQUIRED_TAGS
    if !(t in keys(g.header))
      return false
    end
  end
  return true
end

function query(g::Game, key::String, default="?")
  try
    return g.header[key]
  catch KeyError
    return default
  end
end

function intquery(g::Game, key::String, default=0)
  s = query(g, key)
  try
    t = parse(Int, s)
  catch ArgumentError
    return default
  end
end

function datequery(g::Game, key::String)
  q = split(query(g, key),'.')
  if length(q) < 3
    return Dates.Date(1970,1,1)
  else
    y, m, d = q[1], q[2], q[3]
  end
  if occursin("?", y)
    return Dates.Date(1970,1,1)
  elseif occursin("?", m)
    return Dates.Date(parse(Int,y))
  elseif occursin("?", d)
    return Dates.Date(parse(Int,y), parse(Int,m))
  else
    return Dates.Date(parse(Int,y), parse(Int,m), parse(Int,d))
  end
end

"""
white(g)

Name of white player in game `g`
"""
white(g::Game) = query(g, "White")
"""
black(g)

Name of black player in game `g`
"""
black(g::Game) = query(g, "Black")
"""
date(g)

Date game `g` was played.
"""
date(g::Game) = datequery(g, "Date")
"""
site(g)

Location where game `g` was played.
"""
site(g::Game) = query(g, "Site")
"""
event(g)

Name of event where game `g` was played.
"""
event(g::Game) = query(g, "Event")
"""
result(g)

Result of game `g`. Default is unknown, "*".
"""
result(g::Game) = query(g, "Result", "*")
"""
whiteelo(g)

Elo rating of white player in game `g`. Default is 0.
"""
whiteelo(g::Game) = intquery(g, "WhiteElo")
"""
blackelo(g)

Elo rating of black player in game `g`. Default is 0.
"""
blackelo(g::Game) = intquery(g, "BlackElo")
"""
eco(g)

ECO code for opening in the game `g`.
"""
eco(g::Game) = query(g, "ECO")
"""
eventdate(g)

Date event started in which game `g` was played.
"""
eventdate(g::Game) = datequery(g, "EventDate")
"""
plycount(g)

Ply count of game `g`.
"""
plycount(g::Game) = intquery(g, "PlyCount")

movetext(g::Game) = g.movetext

intresult(g::Game) = RESULT_HASH[query(g, "Result", "1/2-1/2")]

"""
whiteev(g)

Expected score of the white player based on Elo rating in game `g`.

"""
whiteev(g::Game,evslope=400,ev0=0) = ev0 + 1. / (1. + 10^((blackelo(g)-whiteelo(g)) / evslope))
"""
blackev(g)

Expected score of the black player based on Elo rating in game `g`.
"""
blackev(g::Game,evslope=400,ev0=0) = ev0 + 1. / (1. + 10^((whiteelo(g)-blackelo(g)) / evslope))
"""
whitescore(g)

Score for white in game `g`, based on the result.
"""
whitescore(g::Game) = 0.5*(intresult(g) + 1)
"""
blackscore(g)

Score for black in game `g`, based on the result.
"""
blackscore(g::Game) = 0.5*(1 - intresult(g))
"""
whiteperfelo(g)

Performance rating for white based on the result in game `g`.
"""
whiteperfelo(g::Game) =  intresult(g)*400 + blackelo(g)
"""
blackperfelo(g)

Performance rating for black based on the result in game `g`.
"""
blackperfelo(g::Game) = -intresult(g)*400 + whiteelo(g)

"""
isdecisive(g)

Boolean test of whether game `g` had a decisive result.
"""
isdecisive(g::Game) = intresult(g) != 0

const STATE_HEADER = 0
const STATE_MOVES = 1
const STATE_NEWGAME = 2

isblank(line) = all(isspace, line)

"""
sortpgnfile(filename)

Sort the games in a PGN file. (Implementation incomplete.)
"""
function sortpgnfile(pgnfilename, outfile)
  data = readpgn(pgnfilename)
  datasorted=sort(data, by=cpsort)
  f = open(outfile, "w")
  for d in datasorted
    write(f, repr(d))
  end
  close(f)
end

"""
browsepgn(filename)

Browse the games in a PGN file, one by one.
"""
function browsepgn(pgnfilename)
  data = readpgn(pgnfilename)
  datasorted = sort(data, by=cpsort)
  n = 1
  while true
    println("=========================== GAME $n ===========================")
    gamestats(datasorted[n])
    key = read(STDIN,Char)
    if key == 'q'
      break
    elseif key == 'j'
      n += 1
    elseif key == 'k'
      n -= 1
    elseif key == 'J'
      n += 10
    elseif key == 'K'
      n -= 10
    end
    n = n % length(datasorted)  # wrap around
  end
end

function main()
  games = readpgn(ARGS[1], verbose=true)
  println("Read $(length(games)) games.")
  #for g in games
  #  show(g)
  #end
end

PROGRAM_FILE == "PortableGameNotation.jl" && main()
end


text:
##############################################################################
##
## The lightest type that can (i) print table (ii) predict etc
##
##############################################################################

abstract AbstractRegressionResult <: RegressionModel

# fields
coef(x::AbstractRegressionResult) = x.coef
coefnames(x::AbstractRegressionResult) = x.coefnames
vcov(x::AbstractRegressionResult) = x.vcov
nobs(x::AbstractRegressionResult) = x.nobs
df_residual(x::AbstractRegressionResult) = x.df_residual
function confint(x::AbstractRegressionResult)
    scale = quantile(TDist(x.df_residual), 1 - (1-0.95)/2)
    se = stderr(x)
    hcat(x.coef -  scale * se, x.coef + scale * se)
end

# predict, residuals, modelresponse
function predict(x::AbstractRegressionResult, df::AbstractDataFrame)
    rf = deepcopy(x.formula)
    secondstage!(rf)

    newTerms = dropresponse(Terms(rf))
    mf = ModelFrame(newTerms, df)
    newX = ModelMatrix(mf).m

    out = DataArray(Float64, size(df, 1))
    out[mf.msng] = newX * x.coef
end

function residuals(x::AbstractRegressionResult, df::AbstractDataFrame)
    rf = deepcopy(x.formula)
    secondstage!(rf)

    mf = ModelFrame(Terms(rf), df)
    newX = ModelMatrix(mf).m
    out = DataArray(Float64, size(df, 1))
    out[mf.msng] = model_response(mf) -  newX * x.coef
end

function model_response(x::AbstractRegressionResult, df::AbstractDataFrame)
    rf = deepcopy(x.formula)
    secondstage!(rf)
    mf = ModelFrame(Terms(rf), df)
    model_response(mf)
end

# Display Results
function title(::AbstractRegressionResult)
    error("function title has no general method for AbstractRegressionResult")
end

function top(::AbstractRegressionResult)
    error("function top has no general method for AbstractRegressionResult")
end

function coeftable(x::AbstractRegressionResult)
    ctitle = title(x)
    ctop = top(x)
    cc = coef(x)
    se = stderr(x)
    coefnms = coefnames(x)
    conf_int = confint(x)
    # put (intercept) last
    if !isempty(coefnms) && ((coefnms[1] == Symbol("(Intercept)")) || (coefnms[1] == "(Intercept)"))
        newindex = vcat(2:length(cc), 1)
        cc = cc[newindex]
        se = se[newindex]
        conf_int = conf_int[newindex, :]
        coefnms = coefnms[newindex]
    end
    tt = cc ./ se
    CoefTable2(
        hcat(cc, se, tt, ccdf(FDist(1, df_residual(x)), abs2(tt)), conf_int[:, 1], conf_int[:, 2]),
        ["Estimate","Std.Error","t value", "Pr(>|t|)", "Lower 95%", "Upper 95%" ],
        ["$(coefnms[i])" for i = 1:length(cc)], 4, ctitle, ctop)
end

function Base.show(io::IO, x::AbstractRegressionResult)
    show(io, coeftable(x))
end

## Coeftalble2 is a modified Coeftable allowing for a top String matrix displayed before the coefficients.
## Pull request: https://github.com/JuliaStats/StatsBase.jl/pull/119

type CoefTable2
    mat::Matrix
    colnms::Vector
    rownms::Vector
    pvalcol::Integer
    title::AbstractString
    top::Matrix{AbstractString}
    function CoefTable2(mat::Matrix,colnms::Vector,rownms::Vector,pvalcol::Int=0,
                        title::AbstractString = "", top::Matrix = Any[])
        nr,nc = size(mat)
        0 <= pvalcol <= nc || error("pvalcol = $pvalcol should be in 0,...,$nc]")
        length(colnms) in [0,nc] || error("colnms should have length 0 or $nc")
        length(rownms) in [0,nr] || error("rownms should have length 0 or $nr")
        length(top) == 0 || (size(top, 2) == 2 || error("top should have 2 columns"))
        new(mat,colnms,rownms,pvalcol, title, top)
    end
end

## format numbers in the p-value column
function format_scientific(pv::Number)
    return @sprintf("%.3f", pv)
end

function show(io::IO, ct::CoefTable2)
    mat = ct.mat; nr,nc = size(mat); rownms = ct.rownms; colnms = ct.colnms;
    pvc = ct.pvalcol; title = ct.title;   top = ct.top
    if length(rownms) == 0
        rownms = AbstractString[lpad("[$i]",floor(Integer, log10(nr))+3) for i in 1:nr]
    end
    rnwidth = max(4,maximum([length(nm) for nm in rownms]) + 1)
    rownms = [rpad(nm,rnwidth) for nm in rownms]
    widths = [length(cn)::Int for cn in colnms]
    str = [sprint(showcompact,mat[i,j]) for i in 1:nr, j in 1:nc]
    if pvc != 0                         # format the p-values column
        for i in 1:nr
            str[i, pvc] = format_scientific(mat[i, pvc])
        end
    end
    for j in 1:nc
        for i in 1:nr
            lij = length(str[i, j])
            if lij > widths[j]
                widths[j] = lij
            end
        end
    end
    widths .+= 1
    totalwidth = sum(widths) + rnwidth
    if length(title) > 0
        halfwidth = div(totalwidth - length(title), 2)
        println(io, " " ^ halfwidth * string(title) * " " ^ halfwidth)
    end
    if length(top) > 0
        for i in 1:size(top, 1)
            top[i, 1] = top[i, 1] * ":"
        end
        println(io, "=" ^totalwidth)
        halfwidth = div(totalwidth, 2) - 1
        interwidth = 2 +  mod(totalwidth, 2)
        for i in 1:(div(size(top, 1) - 1, 2)+1)
            print(io, top[2*i-1, 1])
            print(io, lpad(top[2*i-1, 2], halfwidth - length(top[2*i-1, 1])))
            print(io, " " ^interwidth)
            if size(top, 1) >= 2*i
                print(io, top[2*i, 1])
                print(io, lpad(top[2*i, 2], halfwidth - length(top[2*i, 1])))
            end
            println(io)
        end
    end
    println("=" ^totalwidth)
    println(io," " ^ rnwidth *
            join([lpad(string(colnms[i]), widths[i]) for i = 1:nc], ""))
    println("-" ^totalwidth)
    for i in 1:nr
        print(io, rownms[i])
        for j in 1:nc
            print(io, lpad(str[i,j],widths[j]))
        end
        println(io)
    end
    println("=" ^totalwidth)
end

##############################################################################
##
## Subtypes of Regression Result
##
##############################################################################

type RegressionResult <: AbstractRegressionResult
    coef::Vector{Float64}   # Vector of coefficients
    vcov::Matrix{Float64}   # Covariance matrix

    esample::BitVector      # Is the row of the original dataframe part of the estimation sample?
    augmentdf::DataFrame

    coefnames::Vector       # Name of coefficients
    yname::Symbol           # Name of dependent variable
    formula::Formula        # Original formula

    nobs::Int64             # Number of observations
    df_residual::Int64      # degrees of freedoms

    r2::Float64             # R squared
    r2_a::Float64           # R squared adjusted
    F::Float64              # F statistics
    p::Float64              # p value for the F statistics
end
title(::RegressionResult) =  "Linear Model"
top(x::RegressionResult) = [
            "Number of obs" sprint(showcompact, nobs(x));
            "Degrees of freedom" sprint(showcompact, nobs(x) - df_residual(x));
            "R2" format_scientific(x.r2);
            "R2 Adjusted" format_scientific(x.r2_a);
            "F Statistic" sprint(showcompact, x.F);
            "p-value" format_scientific(x.p);
            ]

type RegressionResultIV <: AbstractRegressionResult
    coef::Vector{Float64}   # Vector of coefficients
    vcov::Matrix{Float64}   # Covariance matrix

    esample::BitVector      # Is the row of the original dataframe part of the estimation sample?
    augmentdf::DataFrame

    coefnames::Vector       # Name of coefficients
    yname::Symbol           # Name of dependent variable
    formula::Formula        # Original formula

    nobs::Int64             # Number of observations
    df_residual::Int64      # degrees of freedoms

    r2::Float64             # R squared
    r2_a::Float64           # R squared adjusted
    F::Float64              # F statistics
    p::Float64              # p value for the F statistics

    F_kp::Float64           # First Stage F statistics KP
    p_kp::Float64           # First Stage p value KP
end

title(::RegressionResultIV) = "IV Model"
top(x::RegressionResultIV) = [
            "Number of obs" sprint(showcompact, nobs(x));
            "Degrees of freedom" sprint(showcompact, nobs(x) - df_residual(x));
            "R2" format_scientific(x.r2);
            "R2 Adjusted" format_scientific(x.r2_a);
            "F-Statistic" sprint(showcompact, x.F);
            "p-value" format_scientific(x.p);
            "First Stage F-stat (KP)" sprint(showcompact, x.F_kp);
            "First State p-val (KP)" format_scientific(x.p_kp);
            ]

type RegressionResultFE <: AbstractRegressionResult
    coef::Vector{Float64}   # Vector of coefficients
    vcov::Matrix{Float64}   # Covariance matrix

    esample::BitVector      # Is the row of the original dataframe part of the estimation sample?
    augmentdf::DataFrame

    coefnames::Vector       # Name of coefficients
    yname::Symbol           # Name of dependent variable
    formula::Formula        # Original formula

    nobs::Int64             # Number of observations
    df_residual::Int64      # degrees of freedoms

    r2::Float64             # R squared
    r2_a::Float64           # R squared adjusted
    r2_within::Float64
    F::Float64              # F statistics
    p::Float64              # p value for the F statistics

    iterations::Int         # Number of iterations
    converged::Bool         # Has the demeaning algorithm converged?
end
function predict(::RegressionResultFE, ::AbstractDataFrame)
    error("predict is not defined by default for fixed effect models.  Run reg with the the option save = true")
end
function residuals(::RegressionResultFE, ::AbstractDataFrame)
    error("predict is not defined by default for fixed effect models.  Run reg with the the option save = true")
end
title(::RegressionResultFE) = "Fixed Effect Model"
top(x::RegressionResultFE) = [
            "Number of obs" sprint(showcompact, nobs(x));
            "Degrees of freedom" sprint(showcompact, nobs(x) - df_residual(x));
            "R2" format_scientific(x.r2);
            "R2 within" format_scientific(x.r2_within);
            "F-Statistic" sprint(showcompact, x.F);
            "p-value" format_scientific(x.p);
            "Iterations" sprint(showcompact, x.iterations);
            "Converged" sprint(showcompact, x.converged)
            ]

type RegressionResultFEIV <: AbstractRegressionResult
    coef::Vector{Float64}   # Vector of coefficients
    vcov::Matrix{Float64}   # Covariance matrix

    esample::BitVector      # Is the row of the original dataframe part of the estimation sample?
    augmentdf::DataFrame

    coefnames::Vector       # Name of coefficients
    yname::Symbol           # Name of dependent variable
    formula::Formula        # Original formula

    nobs::Int64             # Number of observations
    df_residual::Int64      # degrees of freedoms

    r2::Float64             # R squared
    r2_a::Float64           # R squared adjusted
    r2_within::Float64
    F::Float64              # F statistics
    p::Float64              # p value for the F statistics

    F_kp::Float64           # First Stage F statistics KP
    p_kp::Float64           # First Stage p value KP

    iterations::Int         # Number of iterations
    converged::Bool         # Has the demeaning algorithm converged?
end
function predict(::RegressionResultFEIV, ::AbstractDataFrame)
    error("predict is not defined for fixed effect models. Run reg with the the option savefe = true")
end
function residuals(::RegressionResultFEIV, ::AbstractDataFrame)
    error("residuals is not defined for fixed effect models. Use the function partial_out")
end
title(::RegressionResultFEIV) = "Fixed effect IV Model"
top(x::RegressionResultFEIV) = [
            "Number of obs" sprint(showcompact, nobs(x));
            "Degrees of freedom" sprint(showcompact, nobs(x) - df_residual(x));
            "R2" format_scientific(x.r2);
            "R2 within" format_scientific(x.r2_within);
            "F Statistic" sprint(showcompact, x.F);
            "p-value" format_scientific(x.p);
            "First Stage F-stat (KP)" sprint(showcompact, x.F_kp);
            "First State p-val (KP)" format_scientific(x.p_kp);
            "Iterations" sprint(showcompact, x.iterations);
            "Converged" sprint(showcompact, x.converged)
            ]


text:
# make your scripts automatically re-activate your project
cd(@__DIR__)
using Pkg; Pkg.activate("../..")

# load packages
using ml_for_pt
using DelimitedFiles
using LaTeXStrings
using Plots
using Random
ENV["GKSwstype"]="nul"

# set linear system size and corresponding data folder
L = 6
# L = 8
data_folder = "../../data/mbl_bose_hubbard/L="*string(L)*"/"

# set path to save folder
save_folder = "./results/L="*string(L)*"/"

# critical parameter values for many-body localization transition from literature
p_crit1 = 4.0
p_crit2 = 7.0

# define parameter ranges
p_min = 0.1
p_max = 20.0
n_points = 200
p_range = collect(LinRange(p_min, p_max, n_points))
dp = p_range[2]-p_range[1]
p_range_LBC = collect(p_min-dp/2:dp:p_max+dp/2)

# import data

# probabilities associated with each basis state (Fock basis) at p=p_1
prob = vcat(readdlm(data_folder*"probs_W="*string(1)*".txt")...)

# iterate over entire parameter range
probss = [prob]
for index in 2:201-1
  probb = vcat(readdlm(data_folder*"probs_W="*string(index)*".txt")...)
  push!(probss, probb)
end
samples = collect(1:length(probss[1]))
n_samples = length(samples)

# construct probability distribution from data
# distr(x,p) gives probability to draw the sample x at parameter value p
# samples correspond to enumeration of unique basis states with integers
function distribution(sample, p, p_range, probss)
  p_indx = Int(round((p-p_range[1])/(p_range[2]-p_range[1])))+1
  return probss[p_indx][Int(sample)]
end
distr = (x, p) -> distribution(x, p, p_range ,probss)

# supervised learning using analytical expression
# returns optimal predictions, indicator, and loss
pred_opt_SL, indicator_opt_SL, loss_opt_SL = MLP.get_indicators_SL_analytical(samples, distr, p_range, dp, p_min, p_max)

# prediction-based method using analytical expression
# returns optimal predictions, indicator, and loss
pred_opt_PBM, indicator_opt_PBM, loss_opt_PBM = MLP.get_indicators_PBM_analytical(samples, distr, p_range, dp)

# learning by confusion using analytical expression
# returns optimal indicator and loss
indicator_opt_LBC, loss_opt_LBC = MLP.get_indicators_LBC_analytical(samples, distr, p_range, p_range_LBC)

# save output in save_folder
open(save_folder*"p_range_SL.txt", "w") do io
  writedlm(io, p_range)
end

open(save_folder*"p_range_PBM.txt", "w") do io
  writedlm(io, p_range)
end

open(save_folder*"p_range_LBC.txt", "w") do io
  writedlm(io, p_range_LBC)
end

open(save_folder*"pred_opt_SL.txt", "w") do io
  writedlm(io, pred_opt_SL)
end

open(save_folder*"indicator_opt_SL.txt", "w") do io
  writedlm(io, indicator_opt_SL)
end

open(save_folder*"loss_opt_SL.txt", "w") do io
  writedlm(io, [loss_opt_SL])
end

open(save_folder*"pred_opt_PBM.txt", "w") do io
  writedlm(io, pred_opt_PBM)
end

open(save_folder*"indicator_opt_PBM.txt", "w") do io
  writedlm(io, indicator_opt_PBM)
end

open(save_folder*"loss_opt_PBM.txt", "w") do io
  writedlm(io, [loss_opt_PBM])
end

open(save_folder*"indicator_opt_LBC.txt", "w") do io
  writedlm(io, indicator_opt_LBC)
end

open(save_folder*"loss_opt_LBC.txt", "w") do io
  writedlm(io, loss_opt_LBC)
end

# plot and save results in save_folder
plt = plot(p_range, pred_opt_SL, dpi=300, c="black", label=false)
vline!([p_crit1, p_crit2], label=L"p_{c}", c="red")
xlabel!(L"$p$")
ylabel!(L"$\hat{y}^{\mathrm{opt}}_{\mathrm{SL}}$")
savefig(save_folder*"pred_opt_SL.png")

plt = plot(p_range[2:end-1], indicator_opt_SL, dpi=300, c="black", label=false)
vline!([p_crit1, p_crit2], label=L"p_{c}", c="red")
xlabel!(L"$p$")
ylabel!(L"$I^{\mathrm{opt}}_{\mathrm{SL}}$")
savefig(save_folder*"indicator_opt_SL.png")

plt = plot(p_range, pred_opt_PBM, dpi=300, c="black", label=false)
vline!([p_crit1, p_crit2], label=L"p_{c}", c="red")
xlabel!(L"$p$")
ylabel!(L"$\hat{y}^{\mathrm{opt}}_{\mathrm{PBM}}$")
savefig(save_folder*"pred_opt_PBM.png")

plt = plot(p_range[2:end-1], indicator_opt_PBM, dpi=300, c="black", label=false)
vline!([p_crit1, p_crit2], label=L"p_{c}", c="red")
xlabel!(L"$p$")
ylabel!(L"$I^{\mathrm{opt}}_{\mathrm{PBM}}$")
savefig(save_folder*"indicator_opt_PBM.png")

plt = plot(p_range_LBC, indicator_opt_LBC, dpi=300, c="black", label=false)
vline!([p_crit1, p_crit2], label=L"p_{c}", c="red")
xlabel!(L"$p$")
ylabel!(L"$I^{\mathrm{opt}}_{\mathrm{LBC}}$")
savefig(save_folder*"indicator_opt_LBC.png")

plt = plot(p_range_LBC, loss_opt_LBC, dpi=300, c="black", label=false)
vline!([p_crit1, p_crit2], label=L"p_{c}", c="red")
xlabel!(L"$p$")
ylabel!(L"$\mathcal{L}^{\mathrm{opt}}_{\mathrm{LBC}}$")
savefig(save_folder*"loss_opt_LBC.png")


text:
@testset "Pricing European Option" begin
    pricing_date = Date(2020, 5, 14)
    call = FinancialDSL.Core.european_call(:PETR4, 20.0BRL, Date(2020, 5, 19))
    put = FinancialDSL.Core.european_put(:PETR4, 20.0BRL, Date(2020, 5, 19))

    @testset "internal DailyDatesRange" begin
        dd = InterestRates.DailyDatesRange(pricing_date, call.maturity, InterestRates.BDays252(FinancialDSL.CAL_BRL))
        @test length(dd) == 4
        dd[1] == Date(2020, 5, 14)
        dd[2] == Date(2020, 5, 15)
        dd[3] == Date(2020, 5, 18)
        dd[4] == Date(2020, 5, 19)
        @test dd[end] == call.maturity
        @test InterestRates.yearfractionvalue(dd) ≈ 1/252
    end

    currency_to_curves_map = Dict( "onshore" => Dict( :BRL => :PRE, :USD => :cpUSD, :PETR4 => Symbol("PETR4 DIVIDEND YIELD") ))
    static_model = FinancialDSL.Core.StaticHedgingModel(BRL, currency_to_curves_map)
    binomial_daily_model = FinancialDSL.Core.BinomialModelDaily(static_model, FinancialDSL.Core.Stock(:PETR4), InterestRates.BDays252(BusinessDays.BRSettlement()))
    attr = FinancialDSL.Core.ContractAttributes("riskfree_curves" => "onshore", "carry_type" => "none")
    empty_provider = FinancialDSL.MarketData.EmptyMarketDataProvider()

#=
PETR4 (underlying price)    N/A 20
PRE 15/05/2020  0,9
PRE 18/05/2020  0,8
PRE 19/05/2020  0,7
underlying carry    15/05/2020  0,95
underlying carry    18/05/2020  0,93
underlying carry    19/05/2020  0,91
underlying volatility   N/A 200%
=#

    scenario_fixed = FinancialDSL.Core.FixedScenario()
    scenario_fixed[FinancialDSL.Core.Stock(:PETR4)] = 20.0BRL
    scenario_fixed[FinancialDSL.Core.DiscountFactor(:PRE, Date(2020, 5, 15))] = 0.9
    scenario_fixed[FinancialDSL.Core.DiscountFactor(:PRE, Date(2020, 5, 18))] = 0.8
    scenario_fixed[FinancialDSL.Core.DiscountFactor(:PRE, Date(2020, 5, 19))] = 0.7
    scenario_fixed[FinancialDSL.Core.DiscountFactor(Symbol("PETR4 DIVIDEND YIELD"), Date(2020, 5, 15))] = 0.95
    scenario_fixed[FinancialDSL.Core.DiscountFactor(Symbol("PETR4 DIVIDEND YIELD"), Date(2020, 5, 18))] = 0.93
    scenario_fixed[FinancialDSL.Core.DiscountFactor(Symbol("PETR4 DIVIDEND YIELD"), Date(2020, 5, 19))] = 0.91
    scenario_fixed[FinancialDSL.Core.Volatility(FinancialDSL.Core.Stock(:PETR4))] = 2.00 # 200%

    call_pricer = FinancialDSL.Core.compile_pricer(empty_provider, pricing_date, binomial_daily_model, call, attr)
    put_pricer = FinancialDSL.Core.compile_pricer(empty_provider, pricing_date, binomial_daily_model, put, attr)

    # pricing
    call_p = FinancialDSL.Core.price(call_pricer, scenario_fixed)
    @test call_p ≈ 4.312249551607606
    put_p = FinancialDSL.Core.price(put_pricer, scenario_fixed)
    @test put_p ≈ 0.11224955160760348

    native_call_pricer = FinancialDSL.Core.compile_pricer(empty_provider, pricing_date, binomial_daily_model, call, attr, compiler=:native)
    native_put_pricer = FinancialDSL.Core.compile_pricer(empty_provider, pricing_date, binomial_daily_model, put, attr, compiler=:native)
    @test call_p ≈ FinancialDSL.Core.price(native_call_pricer, scenario_fixed)
    @test put_p ≈ FinancialDSL.Core.price(native_put_pricer, scenario_fixed)

    # exposures
    call_exposures = FinancialDSL.Core.exposures(FinancialDSL.Core.DeltaGammaApproxExposuresMethod(), call_pricer, scenario_fixed)
    put_exposures = FinancialDSL.Core.exposures(FinancialDSL.Core.DeltaGammaApproxExposuresMethod(), put_pricer, scenario_fixed)

    #println("European Call Exposures")
    #for (k, v) in call_exposures
    #    println()
    #    println("$k => $v")
    #    println()
    #end

    native_call_exposures = FinancialDSL.Core.exposures(FinancialDSL.Core.DeltaGammaApproxExposuresMethod(), native_call_pricer, scenario_fixed)
    native_put_exposures = FinancialDSL.Core.exposures(FinancialDSL.Core.DeltaGammaApproxExposuresMethod(), native_put_pricer, scenario_fixed)

    #println("european native_call_exposures = $native_call_exposures")
end

@testset "Pricing American Option" begin
    pricing_date = Date(2020, 5, 14)
    call = FinancialDSL.Core.american_call(:PETR4, 20.0BRL, Date(2020, 5, 19))
    put = FinancialDSL.Core.american_put(:PETR4, 20.0BRL, Date(2020, 5, 19))

    currency_to_curves_map = Dict( "onshore" => Dict( :BRL => :PRE, :USD => :cpUSD, :PETR4 => Symbol("PETR4 DIVIDEND YIELD") ))
    static_model = FinancialDSL.Core.StaticHedgingModel(BRL, currency_to_curves_map)
    binomial_daily_model = FinancialDSL.Core.BinomialModelDaily(static_model, FinancialDSL.Core.Stock(:PETR4), InterestRates.BDays252(BusinessDays.BRSettlement()))
    attr = FinancialDSL.Core.ContractAttributes("riskfree_curves" => "onshore", "carry_type" => "none")
    empty_provider = FinancialDSL.MarketData.EmptyMarketDataProvider()

#=
PETR4 (underlying price)    N/A 20
PRE 15/05/2020  0,9
PRE 18/05/2020  0,8
PRE 19/05/2020  0,7
underlying carry    15/05/2020  0,95
underlying carry    18/05/2020  0,93
underlying carry    19/05/2020  0,91
underlying volatility   N/A 200%
=#

    scenario_fixed = FinancialDSL.Core.FixedScenario()
    scenario_fixed[FinancialDSL.Core.Stock(:PETR4)] = 20.0BRL
    scenario_fixed[FinancialDSL.Core.DiscountFactor(:PRE, Date(2020, 5, 15))] = 0.9
    scenario_fixed[FinancialDSL.Core.DiscountFactor(:PRE, Date(2020, 5, 18))] = 0.8
    scenario_fixed[FinancialDSL.Core.DiscountFactor(:PRE, Date(2020, 5, 19))] = 0.7
    scenario_fixed[FinancialDSL.Core.DiscountFactor(Symbol("PETR4 DIVIDEND YIELD"), Date(2020, 5, 15))] = 0.95
    scenario_fixed[FinancialDSL.Core.DiscountFactor(Symbol("PETR4 DIVIDEND YIELD"), Date(2020, 5, 18))] = 0.93
    scenario_fixed[FinancialDSL.Core.DiscountFactor(Symbol("PETR4 DIVIDEND YIELD"), Date(2020, 5, 19))] = 0.91
    scenario_fixed[FinancialDSL.Core.Volatility(FinancialDSL.Core.Stock(:PETR4))] = 2.00 # 200%

    call_pricer = FinancialDSL.Core.compile_pricer(empty_provider, pricing_date, binomial_daily_model, call, attr)
    put_pricer = FinancialDSL.Core.compile_pricer(empty_provider, pricing_date, binomial_daily_model, put, attr)

    # pricing
    call_p = FinancialDSL.Core.price(call_pricer, scenario_fixed)
    @test call_p ≈ 4.312249551607606
    put_p = FinancialDSL.Core.price(put_pricer, scenario_fixed)
    @test put_p ≈ 0.6732693664447074

    native_call_pricer = FinancialDSL.Core.compile_pricer(empty_provider, pricing_date, binomial_daily_model, call, attr, compiler=:native)
    native_put_pricer = FinancialDSL.Core.compile_pricer(empty_provider, pricing_date, binomial_daily_model, put, attr, compiler=:native)
    @test call_p ≈ FinancialDSL.Core.price(native_call_pricer, scenario_fixed)
    @test put_p ≈ FinancialDSL.Core.price(native_put_pricer, scenario_fixed)

    # exposures
    call_exposures = FinancialDSL.Core.exposures(FinancialDSL.Core.DeltaGammaApproxExposuresMethod(), call_pricer, scenario_fixed)
    put_exposures = FinancialDSL.Core.exposures(FinancialDSL.Core.DeltaGammaApproxExposuresMethod(), put_pricer, scenario_fixed)

    #println("american call_exposures = $call_exposures")

    native_call_exposures = FinancialDSL.Core.exposures(FinancialDSL.Core.DeltaGammaApproxExposuresMethod(), native_call_pricer, scenario_fixed)
    native_put_exposures = FinancialDSL.Core.exposures(FinancialDSL.Core.DeltaGammaApproxExposuresMethod(), native_put_pricer, scenario_fixed)

    #println("american native_call_exposures = $native_call_exposures")
end

@testset "Pricing ZCB" begin
    pricing_date = Date(2020, 5, 14)
    zcb = FinancialDSL.Core.WhenAt(Date(2020, 5, 19), FinancialDSL.Core.Amount(1000.0BRL))

    currency_to_curves_map = Dict( "onshore" => Dict( :BRL => :PRE, :USD => :cpUSD, :PETR4 => Symbol("PETR4 DIVIDEND YIELD") ))
    static_model = FinancialDSL.Core.StaticHedgingModel(BRL, currency_to_curves_map)
    binomial_daily_model = FinancialDSL.Core.BinomialModelDaily(static_model, FinancialDSL.Core.Stock(:PETR4), InterestRates.BDays252(BusinessDays.BRSettlement()))
    attr = FinancialDSL.Core.ContractAttributes("riskfree_curves" => "onshore", "carry_type" => "none")
    empty_provider = FinancialDSL.MarketData.EmptyMarketDataProvider()

    scenario_fixed = FinancialDSL.Core.FixedScenario()
    scenario_fixed[FinancialDSL.Core.Stock(:PETR4)] = 20.0BRL
    scenario_fixed[FinancialDSL.Core.DiscountFactor(:PRE, Date(2020, 5, 15))] = 0.9
    scenario_fixed[FinancialDSL.Core.DiscountFactor(:PRE, Date(2020, 5, 18))] = 0.8
    scenario_fixed[FinancialDSL.Core.DiscountFactor(:PRE, Date(2020, 5, 19))] = 0.7
    scenario_fixed[FinancialDSL.Core.DiscountFactor(Symbol("PETR4 DIVIDEND YIELD"), Date(2020, 5, 15))] = 0.95
    scenario_fixed[FinancialDSL.Core.DiscountFactor(Symbol("PETR4 DIVIDEND YIELD"), Date(2020, 5, 18))] = 0.93
    scenario_fixed[FinancialDSL.Core.DiscountFactor(Symbol("PETR4 DIVIDEND YIELD"), Date(2020, 5, 19))] = 0.91
    scenario_fixed[FinancialDSL.Core.Volatility(FinancialDSL.Core.Stock(:PETR4))] = 2.00 # 200%

    pricer = FinancialDSL.Core.compile_pricer(empty_provider, pricing_date, binomial_daily_model, zcb, attr)
    p = FinancialDSL.Core.price(pricer, scenario_fixed)
    @test p ≈ 1000.0 * 0.7

    zcb_exposures = FinancialDSL.Core.exposures(FinancialDSL.Core.DeltaGammaApproxExposuresMethod(), pricer, scenario_fixed)
    @test length(zcb_exposures) == 1
    @test zcb_exposures[FinancialDSL.Core.DiscountFactor(:PRE, Date(2020, 5, 19))] ≈ 1000.0 * 0.7
end


text:
using DataGraphs
using Dictionaries
using Graphs

vdata = dictionary([1 => "V1", 2 => "V2", 3 => "V3", 4 => "V4"])
edata = dictionary([Edge(1, 2) => "E12", Edge(2, 3) => "E23", Edge(3, 4) => "E34"])
dg = DataGraph(g, vdata, edata)

@show dg[1] == "V1"
@show dg[2] == "V2"
@show dg[3] == "V3"
@show dg[4] == "V4"

@show dg[1 => 2] == "E12"
@show dg[2 => 3] == "E23"
@show dg[3 => 4] == "E34"

@show DataGraph(g) isa DataGraph{Any,Any}
@show DataGraph{String}(g) isa DataGraph{String,Any}
@show DataGraph{<:Any,String}(g) isa DataGraph{Any,String}

vdata = [1 => "V1", 2 => "V2", 3 => "V3", 4 => "V4"]
edata = [(1, 2) => "E12", (2 => 3) => "E23", Edge(3, 4) => "E34"]
dg = DataGraph(g, vdata, edata)

@show dg[1] == "V1"
@show dg[2] == "V2"
@show dg[3] == "V3"
@show dg[4] == "V4"

@show dg[1 => 2] == "E12"
@show dg[2 => 3] == "E23"
@show dg[3 => 4] == "E34"

vdata = ["V1", "V2", "V3", "V4"]
edata = [(1, 2) => "E12", (2 => 3) => "E23", Edge(3, 4) => "E34"]
dg = DataGraph(g, vdata, edata)

@show dg[1] == "V1"
@show dg[2] == "V2"
@show dg[3] == "V3"
@show dg[4] == "V4"

@show dg[1 => 2] == "E12"
@show dg[2 => 3] == "E23"
@show dg[3 => 4] == "E34"

vdata = ["V1", "V2", "V3", "V4"]
edata = ["E12", "E23", "E34"]
dg = DataGraph(g, vdata, edata)

@show dg[1] == "V1"
@show dg[2] == "V2"
@show dg[3] == "V3"
@show dg[4] == "V4"

@show dg[1 => 2] == "E12"
@show dg[2 => 3] == "E23"
@show dg[3 => 4] == "E34"

vdata = ["V1", "V2", "V3", "V4"]
edata = [Edge(1, 2) => "E12", Edge(2, 3) => "E23"]
dg = DataGraph(g; vertex_data=vdata, edge_data=edata)

@show dg[1] == "V1"
@show dg[2] == "V2"
@show dg[3] == "V3"
@show dg[4] == "V4"

@show dg[1 => 2] == "E12"
@show dg[2 => 3] == "E23"
@show !isassigned(dg, 3 => 4)

vdata = ["V1", "V2", "V3", "V4"]
dg = DataGraph(g; vertex_data=vdata)

@show dg[1] == "V1"
@show dg[2] == "V2"
@show dg[3] == "V3"
@show dg[4] == "V4"

@show !isassigned(dg, 1 => 2)
@show !isassigned(dg, 2 => 3)
@show !isassigned(dg, 3 => 4)


text:
include("geneticAlgorithm.jl")

genesLength = 6
populationSize = 10

repeatSize = 3
selectSize = 3
crossoverPoint = 4
mutationPercentage = 0.2

population = generatePopulation(populationSize,genesLength)
calculatePopulationFitness!(population)
printPopulation(population)

popGen, repeatCounter, population = geneticAlgorithm(population,
                                                repeatSize,selectSize,
                                                crossoverPoint,mutationPercentage)

printPopulation(population)
print("Ukupan broj generacija: $(popGen)\n")
print("Ukupan broj ponavljanja: $(repeatCounter)\n")

max_stars_repo_name: pikacxe/Primenjeni-Algoritmi

max_issues_repo_name: pikacxe/Primenjeni-Algoritmi

max_forks_repo_name: pikacxe/Primenjeni-Algoritmi


text:
import NonlinearEquations
using OrdinaryDiffEq
import PyPlot
import SparseArrays

#solve u_t + c * u_x = 0 on a periodic 1d domain
@NonlinearEquations.equations exclude=(c,) function t1d(u, c)
	dx = 1 / length(u)
	NonlinearEquations.setnumequations(length(u))
	if c > 0#use an upwind discretization of u_x
		NonlinearEquations.addterm(1, -c * (u[1] - u[end]) / dx)
		for i = 2:length(u)
			NonlinearEquations.addterm(i, -c * (u[i] - u[i - 1]) / dx)
		end
	else
		for i = 1:length(u) - 1
			NonlinearEquations.addterm(i, -c * (u[i + 1] - u[i]) / dx)
		end
		NonlinearEquations.addterm(length(u), -c * (u[1] - u[end]) / dx)
	end
end

for c in [exp(1), -pi]
	n = 1000
	p = nothing
	tspan = [0.0, 1 / abs(c)]
	u0 = map(i->i < n / 8 || i > 7 * n / 8 ? 1.0 : 0.0, 1:n) + map(x->exp(-(x - 0.5) ^ 2 * 64), range(0, 1; length=n))
	f(u, p, t) = t1d_residuals(u, c)
	f_u(u, p, t) = t1d_u(u, c)
	odef = ODEFunction(f; jac=f_u, jac_prototype=f_u(u0, p, 0.0))
	prob = ODEProblem(odef, u0, tspan, p)
	@time u_explicit = solve(prob, Euler(); dt=0.1 / (abs(c) * n), u0=prob.u0, p=prob.p)
	@time u_implicit = solve(prob, ImplicitEuler(); u0=prob.u0, p=prob.p, abstol=1e-4, reltol=1e-4)

	fig, ax = PyPlot.subplots()
	ax.plot(range(0, 1; length=length(u0)), u0, label="init")
	ax.plot(range(0, 1; length=length(u0)), u_explicit(tspan[end]), label="explicit")
	ax.plot(range(0, 1; length=length(u0)), u_implicit(tspan[end]), label="implicit")
	ax.legend()
	display(fig)
	PyPlot.println()
	PyPlot.close(fig)
end


text:
module WhereTheWaterFlows

using StaticArrays, Requires

export waterflows, fill_dem, catchments,
    plotarea

const I11 = CartesianIndex(1,1)
const I22 = CartesianIndex(2,2)
#const NaNflow = -99
const NOFLOW = 5 # direction number indicating no flow.  Use a constant to better keep track.

"""
Direction numbers.  E.g. dirnums[1,1] will return the number
corresponding to the direction top-left.

Note, I use the conversion that the x-axis corresponds to the row of
the matrix and the y-axis the columns.  To print them in this "normal"
coordinate system use `showme`
"""
const dirnums = SMatrix{3,3}(reverse([ 7      8 9
                                       4 NOFLOW 6
                                       1      2 3]',
                                     dims=2))
# "As dirnums but with no-flow encoded as NaNflow."
# const dirnums_nan = SMatrix{3,3}(reverse([ 7       8 9
#                                            4 NaNflow 6
#                                            1       2 3]',
#                                          dims=2))

"Translation from dirnums to CartesianIndex"
const cartesian = SMatrix{3,3}(reverse(permutedims([CartesianIndex(-1,1)  CartesianIndex(0,1)  CartesianIndex(1,1)
                                                    CartesianIndex(-1,0)  CartesianIndex(0,0)  CartesianIndex(1,0)
                                                    CartesianIndex(-1,-1) CartesianIndex(0,-1) CartesianIndex(1,-1)]),
                                       dims=2))
"""
Show an array with its axes oriented such that they correspond
to x and y direction.
"""
showme(ar) = (display(reverse(ar',dims=1)); println(" "))

"Translate a CartesianIndex, giving the offset, into a direction number"
ind2dir(ind::CartesianIndex) = dirnums[ind + I22]
# ind2dir_nan(ind::CartesianIndex) = dirnums_nan[ind + I22]

"Translate a D8 direction number into a CartesianIndex"
dir2ind(dir) = cartesian[dir]
# "Translate a D8 direction number (with NaNflow at center) into a CartesianIndex"
# dir2ind_nan(dir) = dir==0 ? CartesianIndex(0,0) : dir2ind(dir)

"Translate a D8 direction number into a 2D vector."
dir2vec(dir) = [dir2ind(dir).I...]

"""
Tests whether a cell `J` with flowdir `dirJ` flows into cell `I`.
"""
flowsinto(J::CartesianIndex, dirJ, I::CartesianIndex) = ind2dir(I-J) == dirJ

"Return CartesianIndices corresponding to the 8 neighbors and the point itself"
iterate_D9(I, Iend, I1=I11) = max(I1, I-I1):min(Iend, I+I1)
function iterate_D9(I, ar::AbstractMatrix)
    R = CartesianIndices(size(ar))
    I1, Iend = first(R), last(R)
    return max(I1, I-I1):min(Iend, I+I1)
end

"""
    on_outer_boundary(ar, I::CartesianIndex)

Check whether on outer boundary of ar
"""
function on_outer_boundary(ar, I::CartesianIndex)
    i,j = I.I
    iend, jend = size(ar)
    if i==1 || i==iend || j==1 || j==jend
        return true
    else
        return false
    end
end

"""
    d8dir_feature(dem, bnd_as_pits)

D8 directions of a DEM and drainage features.

Elevations with NaN map to dir==NOFLOW, i.e. just like pits.
However, they are not treated in the pits-filling function `drainpits`.

The argument `bnd_as_pits` determines whether neighboring cells
flow out of the domain or into NaN-cells (`bnd_as_pits==true`)
or not.

Return
- dir  - direction, encoded as `dirnums`
- nout - number of outflow cells of a cell (0 or 1)
- nin  - number of inflow cells of a cell (0-8)
- pits - location of pits as a `Vector{CartesianIndex{2}}` (sorted)
"""
function d8dir_feature(dem, bnd_as_pits)
    # outputs
    diro = zeros(Int8, size(dem))
    nout = falses(size(dem))
    nin = zeros(Int8, size(dem))
    pits = CartesianIndex{2}[]

    R = CartesianIndices(size(dem))
    Iend = last(R)

    # get dir for all points
    for I in R
        # make pits on boundary if bnd_as_pits is set
        if bnd_as_pits && on_outer_boundary(dem,I)
            # make it a pit
            diro[I] = NOFLOW
            continue
        end

        ele = dem[I] # keeps track of lowest elevation of all 9 cells
        dir = NOFLOW
        if isnan(ele)
            # just mark as NOFLOW
        else
            for J in iterate_D9(I, Iend)
                I==J && continue
                ele2 = dem[J]
                if isnan(ele2)
                    if bnd_as_pits
                        # flow into first found NaN-cell
                        dir = ind2dir(J-I)
                        break
                    else
                        # ignore NaN-Cell
                        continue
                    end
                elseif ele > ele2
                    # lower elevation found, adjust dir
                    ele = ele2
                    dir = ind2dir(J-I)
                end
            end
        end
        diro[I] = dir
    end
    # flow features
    for I in R
        for J in iterate_D9(I, Iend)
            J==I && continue
            nin[I] += flowsinto(J, diro[J], I)
        end
        if diro[I]==NOFLOW
            nout[I] = false
            if !isnan(dem[I])
                push!(pits, I)
            elseif nin[I]>0
                # mark NaN points only as pits if something is flowing into them
                push!(pits, I)
            end
        else
            nout[I] = true
        end
    end

    return diro, nout, nin, pits
end

"""
    waterflows(dem, cellarea=ones(size(dem));
               maxitr=1000, calc_streamlength=true, drain_pits=true, bnd_as_pits=false)

Does the water flow routing according the D8 algorithm.

kwargs:
- drain_pits -- whether to route through pits
- maxiter -- maximum iterations of the algorithm
- calc_streamlength -- whether to calculate stream length
- bnd_as_pits -- whether the domain boundary an NaNs should be pits,
                 i.e. adjacent cells can drain into them,
                 or whether to ignore them

TODO: if bnd_as_pits rout water along boundary edges first. This would probably
substantially reduce the number of catchments, as currently every boundary point
is a pit and thus a catchment (if bnd_as_pits==true).

Returns
- area -- upslope area
- stream-length -- length of stream to the farthest source
- dir -- flow direction at each location
- nout -- whether the point has outlflow.  I.e. nout[I]==0 --> I is a pit
- nin -- number of inflow cells
- pits -- location of pits as Vector{CartesianIndex{2}}
- c -- catchment map
- bnds -- boundaries between catchments.  The boundary to the exterior/NaNs is not in here.
"""
function waterflows(dem, cellarea=ones(size(dem));
                     maxitr=1000, calc_streamlength=true, drain_pits=true, bnd_as_pits=false)
    area, slen, dir, nout, nin, pits = _waterflows(d8dir_feature(dem, bnd_as_pits)..., cellarea; maxitr=maxitr, calc_streamlength=calc_streamlength)
    c, bnds = catchments(dir, pits, bnd_as_pits)
    if drain_pits
        dir, nin, nout, pits, c, bnds = drainpits(dem, dir, nin, nout, pits, (c,bnds))
        area, slen = _waterflows(dir, nout, nin, pits)
    end
    #area[isnan.(dem)] .= NaN
    return area, slen, dir, nout, nin, pits, c, bnds
end
# this function does the actual routing
function _waterflows(dir, nout, nin, pits, cellarea=ones(size(dir)); maxitr=1000, calc_streamlength=true)
    area = copy(cellarea)
    tmp = convert(Matrix{Int}, nin)
    tmp2 = calc_streamlength ? copy(tmp) : tmp

    for counter = 1:maxitr
        n = 0
        for R in CartesianIndices(size(dir))
            # Consider only points with less than `counter` upstream points
            if tmp[R]==0
                tmp[R] = -counter # done with it
                if calc_streamlength
                    tmp2[R] = -counter # done with it
                end
                d = dir[R]
                if d!=NOFLOW
                    receiver = R + dir2ind(d)
                    area[receiver] += area[R]
                    tmp2[receiver] -= 1
                    tmp2[receiver] < 0 && error("This should not happen!")
                end
                n +=1
            end
        end
        if n==0
            #@show counter
            break
        end
        if counter==maxitr
            error("Maximum number of iterations reached in `_waterflows`: $counter")
        end
        copyto!(tmp, tmp2)
    end
    return area, -tmp, dir, nout, nin, pits
end

"""
    catchments(dir, pits, bnd_as_pits)

Calculate catchments from
- dir
- pits

Return: catchments Matrix{Int}.  Value==0 corresponds to NaNs in the DEM
which are not pits (i.e. where no water flows into).
"""
function catchments(dir, pits, bnd_as_pits)
    c = zeros(Int, size(dir))
    np = length(pits)
    # recursively traverse the drainage tree in up-flow direction,
    # starting at all pits
    for (n, pit) in enumerate(pits)
        _catchments!(n, c, dir, pit)
    end

    return c, make_boundaries(c, 1:np)
end
function _catchments!(n, c, dir, ij)
    c[ij] = n
    # proc upstream points
    for IJ in iterate_D9(ij, c)
        ij==IJ && continue
        if flowsinto(IJ, dir[IJ], ij)
            _catchments!(n, c, dir, IJ)
        end
    end
end

"""
    make_boundaries(catchments, colors, bnd_as_pits)

Make vectors of boundary points.  Note that points along the
edge of the domain as well as points bordering NaN-cells with no
inflow do not count as boundaries.

TODO: this is brute force...
"""
function make_boundaries(catchments, colors)
    bnds = [CartesianIndex[] for c in colors]
    for R in CartesianIndices(size(catchments))
        c = catchments[R]
        c==0 && continue # don't find boundaries for c==0 (NaNs with no inflow)
        bnd = bnds[c]
        for I in iterate_D9(R, catchments)
            co = catchments[I]
            if co!=c && co!=0
                push!(bnd, R)
                break
            end
        end
    end#
    bnds
end

"""
Checks that all points in `bnds[color]` are indeed on the
boundary; removes them otherwise.
"""
function _prune_boundary!(bnds, catchments::Matrix, color)
    del = Int[]
    # loop over all boundary cells
    for (i,P) in enumerate(bnds[color])
        keep = false
        # check cells around it
        for PP in iterate_D9(P, catchments)
            # if one is of different color, keep it.
            co = catchments[PP]
            if co!=color && co!=0
                keep = true
                break
            end
        end
        if !keep
            push!(del, i)
        end
    end
    deleteat!(bnds[color], del)
    return nothing
end

"""
     drainpits(dem, dir, nin, nout, pits, (c, bnds)=catchments(dir, pits);
               maxitr=100)

Return an update direction field which drains (interior) pits.
This is done by reversing the flow connecting the lowest boundary point
to the pits.

There needs to be a decision on how to treat the outer boundary and boundaries
to NaN-cells.  Possibilities:
- do not treat boundary cells as pits but do treat pits at the boundary as terminal.
- treat boundary cells as pits, i.e. flow reaching such a cell will vanish.  Again
  such cells would be terminal.

What to do if there are no terminal pits in the DEM?  Fill to the uppermost pit?  Or take
the lowermost as terminal?

Returns new dir, nin, nout, pits (sorted), c.

TODO: this is the performance bottleneck.
"""
function drainpits(dem, dir, nin, nout, pits, (c, bnds)=catchments(dir, pits);
                   maxitr=100)
    dir_ = copy(dir)
    dir2 = copy(dir)
    nin_ = copy(nin)
    nout_ = copy(nout)
    pits_ = copy(pits)
    c_ = copy(c)
    bnds_ = deepcopy(bnds)

    Iend = CartesianIndex(size(dem))

    pits_to_keep = trues(length(pits_))

    no_drainage_across_boundary = false
    # iterate until all interior pits are removed
    for i=1:maxitr
        n_removed = 0
        for (color, P) in enumerate(pits_)
            copyto!(dir2, dir_) # TODO hack...

            # Already removed pit, skip
            P==CartesianIndex(-1,-1) && continue
            # Don't process pits on the DEM boundary, because there water disappears.
            (P.I[1]==1 || P.I[2]==1 || P.I[1]==size(dir,1) || P.I[2]==size(dir,2)) && continue
            # Don't process pits which are NaNs, because there water disappears.
            # See also bnd_as_pits.
            isnan(dem[P]) && continue
            # If there are no more boundaries left, stop.  This should only occur when
            # bnd_as_pits==false and when there are only interior pits.
            if all(isempty.(bnds_))
                no_drainage_across_boundary = true # set flag to warn later
                break
            end

            ## Debug plotting
            # cls = [c=>sum(c_.==c) for c=1:6]
            # @show color, size(bnds_[color])
            # @show cls
            # imshow(Array(c_'), origin="lower")
            # if color==1
            #     colorbar()
            # end
            # readline(stdin)
            # if color==6
            #     imshow(Array(c_'), origin="lower")
            # end

            # find point on catchment boundary with minimum elevation
            Imin = bnds_[color][findmin(getindex.(Ref(dem), bnds_[color]))[2]]
            @assert c_[Imin]==c_[P] # Something is amiss if the found minimum is the pit!
            # make the outflow and find the next catchment
            min_ = Inf
            target = Imin
            for J in iterate_D9(Imin, dem)
                Imin==J && continue # do not look at the point itself
                c_[J] == color && continue # J is in Imin's catchment
                if dem[J] < min_
                    min_ = dem[J]
                    target = J
                end
            end
            if target==Imin
                # this means that point is on a NaN boundary -> don't process
                continue
            end
            # do the flow across the boundary
            _flow_from_to!(Imin, target, dir_, nin_, nout_)

            # reverse directions on path going from Imin to P
            if Imin!=P
                P1 = Imin
                while true
                    # get downstream point
                    P2 = dir2ind(dir2[P1]) + P1 # note usage of `dir2`
                    # The need for dir2 could be removed by first
                    # traversing the flow path without modifying it.
                    # And then do the modification in a second traversal.
                    if P2==P1
                        # With dir2, this should be fixed now.
                        println("This should not happen!")
                        @show color, Imin, P, P1
                        break
                    end
                    # reverse
                    _flow_from_to!(P2, P1, dir_, nin_, nout_)

                    P2==P && break # reached the pit
                    P1 = P2
                end
            end

            # remove from list of pits
            pits_to_keep[color] = false
            pits_[color] = CartesianIndex(-1,-1)

            # update catchments
            othercolor = c_[target]
            c_[c_.==color] .= othercolor
            append!(bnds_[othercolor], bnds_[color])
            empty!(bnds_[color])
            _prune_boundary!(bnds_, c_, othercolor)
            n_removed +=1
        end
        n_removed==0 && break
        if i==maxitr
            error("Maximum number of iterations reached in `drainpits`: $i")
        end
    end

    if no_drainage_across_boundary
        @warn """Water cannot leave this DEM!  Instead it drains into a random interior pit.
                 Consider setting `bnd_as_pits=true`."""
    end

    # remove removed pits
    pits_ = pits_[pits_.!=Ref(CartesianIndex(-1, -1))]
    # redo colors
    cols = Dict([c=>i for (i,c) in enumerate(unique(c_))])
    for i in eachindex(c_)
        c_[i] = cols[c_[i]]
    end

    return dir_, nin_, nout_, sort(pits_), c_, bnds_
end

"""
Update dir, nin, and nout such that flow at P1 is now
from P1 to P2.

It can potentially modify `dir, nin, nout` at three locations
- P1: dir, nout, nin
- P2: nin
  - if allow_reversion==true, then P2's dir, nout can also be modified.
- P3 (previous receiver cell of P1): nin
"""
function _flow_from_to!(P1, P2, dir, nin, nout, allow_reversion=false)
    # already right
    ind2dir(P2-P1)==dir[P1] && return nothing

    # otherwise update
    P3 = P1 + dir2ind(dir[P1]) # previous receiver cell
    if P3!=P1
        nin[P3] -= 1
    end
    dir[P1] = ind2dir(P2-P1)
    nin[P2] += 1
    nout[P1] = true # definitely not a pit anymore

    # if flow from P2 was into P1, then make P2 a pit
    # (otherwise dir will be inconsistent)
    if (dir2ind(dir[P2]) + P2 == P1)
        if allow_reversion
            dir[P2] = ind2dir(CartesianIndex(0,0))
            nout[P2] = false
            nin[P1] -= 1
        else
            error("Flow direction in P2 would need to be reversed but not allowed (set option `allow_reversion`).")
        end
    end
    return nothing
end

"""
    fill_dem(dem, pits, dir)

Fill the pits (aka sinks) of a DEM (apply this after applying "drainpits",
which is the default). Returns the filled DEM.

Note, this is not needed as pre-processing step to use `upstream` area.

This uses a tree traversal to fill the DEM. It does it depth-first (as it
is easier) which may lead to a stack overflow on a large DEM.
"""
function fill_dem(dem, pits, dir; small=0.0)
    dem = copy(dem)
    for pit in pits
        npts = _fill_ij!(0.0, dem, pit, dir, small)
    end
    return dem
end

function _fill_ij!(ele, dem, ij, dir, small)
    if ele > dem[ij]
        dem[ij] = ele
        ele += small
    else
        ele = dem[ij]
    end
    # proc upstream points
    for IJ in iterate_D9(ij, dem)
        ij==IJ && continue
        if flowsinto(IJ, dir[IJ], ij)
            _fill_ij!(ele, dem, IJ, dir, small)
        end
    end
end

## Plotting
function __init__()
    @require PyPlot="d330b81b-6aea-500a-939a-2ce795aea3ee" include("plotting.jl")
end

end # module


text:
using ContinuousBenchmarks

s = Ref(0)

res = @benchmarkd("Example", for i in 1:100 s[] += i end)

@show s[]

LOG_DATA = Dict(
    "name" => res[1],
    "engine" => "None",
    "time" => res[3],
    "mem" => res[4],
)