text:
# ===========================================================================
# VOLUME CONTROL BLOCKETTES
function blk_010!(io::IO, nb::Int64, v::Integer)
  sio = blk_string_read(io, nb, v)
  v > 1 && println("")
  close(sio)
  return nothing
end

function blk_011!(io::IO, nb::Int64, v::Integer)
  sio = blk_string_read(io, nb, v)
  if v > 1
    println("")
    if v > 2
      nsta = stream_int(sio, 3)
      println(" "^16, lpad("STA", 5), " | SEQ")
      println(" "^16, "------|-------")
      @inbounds for i = 1:nsta
        sta = String(fastread(sio, 5))
        seq = stream_int(sio, 6)
        println(" "^16, sta, " | ", seq)
      end
    end
  end
  close(sio)
  return nothing
end

function blk_012!(io::IO, nb::Int64, v::Integer)
  sio = blk_string_read(io, nb, v)
  if v > 1
    N = stream_int(sio, 4)
    println(", N_SPANS = ", N)
  end
  close(sio)
  return nothing
end

# ===========================================================================
# ABBREVIATION CONTROL HEADERS

# [30] is out of scope for SeisIO
function blk_030!(io::IO, nb::Int64, v::Integer)
  sio = blk_string_read(io, nb, v)
  if v > 1
    desc = string_field(sio)
    code = stream_int(sio, 4)
    fam = stream_int(sio, 3)
    println(", CODE ", code, ", DESC ", desc, ", FAM ", fam)
    if v > 2
      ND = stream_int(sio, 2)
      @inbounds for i = 1:ND
        desc = string_field(sio)
        println(" "^18, desc)
      end
    end
  end
  close(sio)
  return nothing
end

function blk_031!(io::IO, nb::Int64, v::Integer)
  sio = blk_string_read(io, nb, v)
  k = stream_int(sio, 4)
  class = fastread(sio)
  desc = string_field(sio)
  if v > 1
    units = stream_int(sio, 3)
    println(", CODE ", k, ", CLASS ", Char(class), ", DESC ", desc, ", UNITS ", units)
  end
  close(sio)
  comments[k] = replace(desc, "," => "-")

  return nothing
end

function blk_032!(io::IO, nb::Int64, v::Integer)
  sio = blk_string_read(io, nb, v)
  if v > 1
    code = stream_int(sio, 2)
    pub = string_field(sio)
    date = string_field(sio)
    name = string_field(sio)
    println(", CODE ", code, ", REFERENCE ", pub, ", DATE ", date, ", PUBLISHER ", name)
  end
  close(sio)

  return nothing
end

function blk_033!(io::IO, nb::Int64, v::Integer)
  sio = blk_string_read(io, nb, v)

  code = stream_int(sio, 3)
  desc = string_field(sio)
  close(sio)

  v > 1 && println(", CODE ", code, ", DESC ", desc)

  abbrev[code] = desc

  return nothing
end

function blk_034!(io::IO, nb::Int64, v::Integer)
  sio = blk_string_read(io, nb, v)

  code = stream_int(sio, 3)
  name = string_field(sio)
  if v > 1
    desc = string_field(sio)
    println(", CODE ", code, ", NAME ", name, ", DESC ", desc)
  end
  close(sio)

  units_lookup[code] = fix_units(name)

  return nothing
end

function blk_041!(io::IO, nb::Int64, v::Integer, units::Bool)
  sio = blk_string_read(io, nb, v)

  resp_lookup_key = stream_int(sio, 4)
  skip_string!(sio)
  symm_code = fastread(sio)
  uic = stream_int(sio, 3)
  uoc = stream_int(sio, 3)

  #= NF can be completely wrong here. Two issues:
  (1) The SEED manual quietly suggests writing only part of "symmetric" FIR
  filters to file; in which case NF is not the length of the vector we need. =#
  NF = stream_int(sio, 4)
  if symm_code == 0x42
    NF = div(NF, 2)+1
  elseif symm_code == 0x43
    NF = div(NF, 2)
  end

  #= (2) buried in the SEED manual, it's noted that the maximum size of
  a non-data blockette length is 9999 chars, and type [41] can exceed it.
  In this case, the sample files tell me that NF is dead wrong; workaround is
  to read as many FIR values as possible and fix with append! of next packet.

  We've already used 7 chars on blockette type & size; rest are in sio.
  =#
  p = fastpos(sio)
  NF_true = min(NF, div(9992-p, 14))

  if v > 1
    println(", KEY ", resp_lookup_key, ", SYMM ", Char(symm_code))
    if v > 2
      println(" "^16, "units in = ", units_lookup[uic])
      println(" "^16, "units out = ", units_lookup[uoc])
    end
    println(" "^16, "p = ", p, ", NF = ", NF, ", NF_true = ", NF_true)
  end
  #= This appears to fix it but results in some very odd behavior. =#

  F = Array{Float64, 1}(undef, NF_true)
  @inbounds for i = 1:NF_true
    fast_readbytes!(sio, BUF.hdr_old, 14)
    setindex!(F, buf_to_double(BUF.hdr_old, 14), i)
  end
  close(sio)

  if v > 2
    [println("F[", i, "] = ", F[i]) for i = 1:NF_true]
  end

  # Process to resps

  if haskey(responses, resp_lookup_key)
    R = responses[resp_lookup_key][1]
    append!(R.b, F)
  else
    resp = CoeffResp(F, Float64[])
    ui = units ? fix_units(units_lookup[uic]) : ""
    uo = units ? fix_units(units_lookup[uoc]) : ""
    responses[resp_lookup_key] = (resp, ui, uo)
  end
  return nothing
end

function blk_043!(io::IO, nb::Int64, v::Integer, units::Bool)
  sio = blk_string_read(io, nb, v)
  resp_lookup_key = stream_int(sio, 4)
  skip_string!(sio)
  fastskip(sio, 1)
  uic = stream_int(sio, 3)
  uoc = stream_int(sio, 3)
  fast_readbytes!(sio, BUF.hdr_old, 12)
  A0 = buf_to_double(BUF.hdr_old, 12)
  fast_readbytes!(sio, BUF.hdr_old, 12)
  F0 = buf_to_double(BUF.hdr_old, 12)
  NZ = stream_int(sio, 3)
  Z = Array{ComplexF64, 1}(undef, NZ)
  @inbounds for i = 1:NZ
    fast_readbytes!(sio, BUF.hdr_old, 12)
    rr = buf_to_double(BUF.hdr_old, 12)
    fast_readbytes!(sio, BUF.hdr_old, 12)
    ii = buf_to_double(BUF.hdr_old, 12)
    Z[i] = complex(rr, ii)
    fastskip(sio, 24)
  end
  NP = stream_int(sio, 3)
  P = Array{ComplexF64,1}(undef, NP)
  @inbounds for i = 1:NP
    fast_readbytes!(sio, BUF.hdr_old, 12)
    rr = buf_to_double(BUF.hdr_old, 12)
    fast_readbytes!(sio, BUF.hdr_old, 12)
    ii = buf_to_double(BUF.hdr_old, 12)
    P[i] = complex(rr, ii)
    fastskip(sio, 24)
  end
  close(sio)

  if v > 1
    println("A0 = ", string(A0), "F0 = ", string(F0))
    if v > 2
      println(" "^16, "NZ = ", NZ, ":")
      @inbounds for i = 1:NZ
        println(" "^16, Z[i])
      end
      println(" "^16, "NP = ", NP, ":")
      @inbounds for i = 1:NP
        println(" "^16, P[i])
      end
    end
  end

  # Process to resps
  R = PZResp64(a0 = A0, f0 = F0, z = Z, p = P)
  ui = units ? fix_units(units_lookup[uic]) : ""
  uo = units ? fix_units(units_lookup[uoc]) : ""
  responses[resp_lookup_key] = (R, ui, uo)

  return nothing
end

function blk_044!(io::IO, nb::Int64, v::Integer, units::Bool)
  sio = blk_string_read(io, nb, v)

  resp_lookup_key = stream_int(sio, 4)
  skip_string!(sio)
  fastskip(sio, 1)
  uic = stream_int(sio, 3)
  uoc = stream_int(sio, 3)
  NN = stream_int(sio, 4)
  N = Array{Float64,1}(undef, NN)
  @inbounds for i = 1:NN
    fast_readbytes!(sio, BUF.hdr_old, 12)
    N[i] = buf_to_double(BUF.hdr_old, 12)
    fastskip(sio, 12)
  end
  ND = stream_int(sio, 4)
  D = Array{Float64,1}(undef, ND)
  @inbounds for i = 1:ND
    fast_readbytes!(sio, BUF.hdr_old, 12)
    D[i] = buf_to_double(BUF.hdr_old, 12)
    fastskip(sio, 12)
  end
  close(sio)

  if v > 1
    println(", KEY ", resp_lookup_key)
    if v > 2
      println(" "^16, "units in code = ", uic)
      println(" "^16, "units out code = ", uoc)
      println(" "^16, "NN = ", NN, ":")
      @inbounds for i = 1:NN
        println(" "^16, N[i])
      end
      println(" "^16, "ND = ", ND, ":")
      @inbounds for i = 1:min(ND)
        println(" "^16, D[i])
      end
    end
  end

  if units
    responses[resp_lookup_key] = (CoeffResp(N, D),
                                  fix_units(units_lookup[uic]),
                                  fix_units(units_lookup[uoc])
                                  )
  else
    responses[resp_lookup_key] = (CoeffResp(N, D), "", "")
  end
  return nothing
end

function blk_047!(io::IO, nb::Int64, v::Integer)
  sio = blk_string_read(io, nb, v)
  resp_lookup_key = stream_int(sio, 4)
  skip_string!(sio)
  fast_readbytes!(sio, BUF.hdr_old, 10)
  fs = buf_to_double(BUF.hdr_old, 10)
  fac = stream_int(sio, 5)
  os = stream_int(sio, 5)
  fast_readbytes!(sio, BUF.hdr_old, 11)
  delay = buf_to_double(BUF.hdr_old, 11)
  fast_readbytes!(sio, BUF.hdr_old, 11)
  corr = buf_to_double(BUF.hdr_old, 11)
  close(sio)

  if v > 1
    println(", KEY ", resp_lookup_key)
    if v > 2
      println(" "^16, "fs = ", fs)
      println(" "^16, "decimation factor = ", fac)
      println(" "^16, "decimation offset = ", os)
      println(" "^16, "delay = ", delay)
      println(" "^16, "delay correction applied = ", corr)
    end
  end

  responses[resp_lookup_key] = Blk47(fs, delay, corr, fac, os)

  return nothing
end

function blk_048!(io::IO, nb::Int64, v::Integer)
  sio = blk_string_read(io, nb, v)

  resp_lookup_key = stream_int(sio, 4)
  skip_string!(sio)
  fast_readbytes!(sio, BUF.hdr_old, 12)
  gain = buf_to_double(BUF.hdr_old, 12)
  fast_readbytes!(sio, BUF.hdr_old, 12)
  fg = buf_to_double(BUF.hdr_old, 12)
  nv = stream_int(sio, 2)
  # channel histories are not in the scope of SeisIO
  @inbounds for i = 1:nv
    # fastskip(sio, 24)
    skip_string!(sio) # should cover the 24-Byte channel history
  end
  close(sio)

  if v > 1
    println(", KEY ", resp_lookup_key)
    if v > 2
      println(" "^16, "gain = ", gain, " (f = ", fg, " Hz)")
    end
  end

  responses[resp_lookup_key] = Blk48(gain, fg)

  return nothing
end

function blk_050(io::IO, nb::Int64, v::Integer)
  sio = blk_string_read(io, nb, v)

  # Station
  fill!(BUF.hdr, 0x00)
  p = pointer(BUF.hdr)
  fast_unsafe_read(sio, p, 5)

  #= I make the assumption here that channel coordinates are
  correctly set in blockette 52; if this is not true, then
  :loc will not be set.

  This shortcut makes support for multiplexing impossible, though
  blockette-50 multiplexing has never been encountered in my
  test files (or in any files of the ObsPy test suite.)

  Multiplexing would be annoying as one would need to define
  new SeisData objects S_subnet of length = n_subchans and use
  append!(S, S_subnet), rather than push!(S, C); indexing would
  be a mess. But it's doable in theory.
  =#
  fastskip(sio, 33)
  site_name = strip(string_field(sio))
  fastskip(sio, 9)
  ts = string_field(sio)
  te = string_field(sio)
  uc = fastread(sio)

  # Network
  p = pointer(BUF.hdr, 11)
  fast_unsafe_read(sio, p, 2)

  close(sio)

  if v  > 1
    println(", ID = ", String(copy(BUF.id)))
    if v > 2
      println(" "^16, "site name = ", site_name)
      println(" "^16, "start date = ", ts)
      println(" "^16, "end date = ", te)
      println(" "^16, "update code = ", Char(uc))
    end
  end

  return site_name
end

# ===========================================================================
# STATION CONTROL BLOCKETTES

# not necessary
function blk_051!(io::IO, nb::Int64, v::Integer)
  sio = blk_string_read(io, nb, v)
  if v > 1
     println("")
    if v > 2
      ts = string_field(sio)
      te = string_field(sio)
      k = stream_int(sio, 4)
      comment_level = stream_int(sio, 6)
      println(" "^16, "ts = ", ts)
      println(" "^16, "te = ", te)
      println(" "^16, "comment code key #", k)
      println(" "^16, "comment level = ", comment_level)
    end
  end
  close(sio)
  return nothing
end

function blk_052!(io::IO, nb::Int64, C::SeisChannel, ts_req::Int64, te_req::Int64, v::Integer)
  sio = blk_string_read(io, nb, v)

  # loc
  p = pointer(BUF.hdr, 6)
  fast_unsafe_read(sio, p, 2)

  # cha
  p = pointer(BUF.hdr, 8)
  fast_unsafe_read(sio, p, 3)
  fastskip(sio, 4)
  inst = stream_int(sio, 3)
  skip_string!(sio)
  units_code = stream_int(sio, 3)
  fastskip(sio, 3)

  # lat, lon, el, dep, az, inc
  fast_readbytes!(sio, BUF.hdr_old, 10)
  lat = buf_to_double(BUF.hdr_old, 10)
  fast_readbytes!(sio, BUF.hdr_old, 11)
  lon = buf_to_double(BUF.hdr_old, 11)
  fast_readbytes!(sio, BUF.hdr_old, 7)
  el = buf_to_double(BUF.hdr_old, 7)
  fast_readbytes!(sio, BUF.hdr_old, 5)
  dep = buf_to_double(BUF.hdr_old, 5)
  fast_readbytes!(sio, BUF.hdr_old, 5)
  az = buf_to_double(BUF.hdr_old, 5)
  fast_readbytes!(sio, BUF.hdr_old, 5)
  inc = 90.0 - buf_to_double(BUF.hdr_old, 5)
  fastskip(sio, 6)

  # fs
  fast_readbytes!(sio, BUF.hdr_old, 10)
  fs = buf_to_double(BUF.hdr_old, 10)

  # don't really need max. drift; skip_string passes it over
  skip_string!(sio)

  # ts, te
  ts = parse_resp_date(sio, BUF.uint16_buf)
  te = parse_resp_date(sio, BUF.uint16_buf)
  if te == -56504908800000000
    te = 19880899199000000
  end
  close(sio)

  if v > 1
    println(", ID = ", String(copy(BUF.hdr)), ", INST = ", inst)
    if v > 2
      println(" "^16, "lat = ", lat, ", lon = ", lon, ", z = ", el, ", dep = ", dep, ", θ = ", az, ", ϕ = ", inc)
      println(" "^16, "fs = ", fs)
      println(" "^16, "ts = ", u2d(div(ts, sμ)))
      println(" "^16, "te = ", u2d(div(te, sμ)))
    end
  end

  if ts ≤ te_req && te ≥ ts_req
    update_hdr!(BUF)
    C.id = getfield(BUF, :id_str)
    C.fs = fs
    C.loc = GeoLoc(lat = lat, lon = lon, el = el, dep = dep, az = az, inc = inc)
    C.resp = MultiStageResp(12)
    C.misc["ts"] = ts
    C.misc["te"] = te
    C.misc["timespan"] = string(u2d(div(ts, 1000000))) * " : " *  string(u2d(div(te, 1000000)))
    C.misc["inst"] = get(abbrev, inst, "")
    units = get(units_lookup, units_code, "")
    C.units = fix_units(units)
    skipping = false
  else
    if v > 1
      println(" "^16, "Skipping ", String(copy(BUF.hdr)), " (not in requested time range)")
    end
    skipping = true
  end
  return skipping
end

function blk_053(io::IO, nb::Int64, v::Integer, R::MultiStageResp, units::Bool)
  sio = blk_string_read(io, nb, v)
  tft = Char(fastread(sio))
  stage = stream_int(sio, 2)
  uic = stream_int(sio, 3)
  uoc = stream_int(sio, 3)
  fast_readbytes!(sio, BUF.hdr_old, 12)
  A0 = buf_to_double(BUF.hdr_old, 12)
  fast_readbytes!(sio, BUF.hdr_old, 12)
  F0 = buf_to_double(BUF.hdr_old, 12)
  NZ = stream_int(sio, 3)
  Z = Array{ComplexF64,1}(undef, NZ)
  @inbounds for i = 1:NZ
    fast_readbytes!(sio, BUF.hdr_old, 12)
    rr = buf_to_double(BUF.hdr_old, 12)
    fast_readbytes!(sio, BUF.hdr_old, 12)
    ii = buf_to_double(BUF.hdr_old, 12)
    Z[i] = complex(rr, ii)
    fastskip(sio, 24)
  end
  NP = stream_int(sio, 3)
  P = Array{ComplexF64,1}(undef, NP)
  @inbounds for i = 1:NP
    fast_readbytes!(sio, BUF.hdr_old, 12)
    rr = buf_to_double(BUF.hdr_old, 12)
    fast_readbytes!(sio, BUF.hdr_old, 12)
    ii = buf_to_double(BUF.hdr_old, 12)
    P[i] = complex(rr, ii)
    fastskip(sio, 24)
  end
  close(sio)

  if v > 1
    println(", TFT ", tft, ", STAGE ", stage)
    if v > 2
      println(" "^16, "units in code #", uic)
      println(" "^16, "units out code #", uoc)
      println(" "^16, "A0 = ", A0)
      println(" "^16, "F0 = ", F0)
      println(" "^16, "NZ = ", NZ, ":")
      @inbounds for i = 1:NZ
        println(" "^16, Z[i])
      end
      println(" "^16, "NP = ", NP, ":")
      @inbounds for i = 1:NP
        println(" "^16, P[i])
      end
    end
  end

  if stage > length(R.fs)
    append!(R, MultiStageResp(6))
  end
  resp = PZResp64(a0 = A0, f0 = F0, z = Z, p = P)
  ui = units ? fix_units(units_lookup[uic]) : ""
  uo = units ? fix_units(units_lookup[uoc]) : ""
  R.stage[stage] = resp
  R.i[stage] = ui
  R.o[stage] = uo
  return stage
end

function blk_054(io::IO, nb::Int64, v::Integer, R::MultiStageResp, units::Bool)
  sio = blk_string_read(io, nb, v)

  fastskip(sio, 1)
  stage = stream_int(sio, 2)
  uic = stream_int(sio, 3)
  uoc = stream_int(sio, 3)
  NN = stream_int(sio, 4)
  N = Array{Float64,1}(undef, NN)
  @inbounds for i = 1:NN
    fast_readbytes!(sio, BUF.hdr_old, 12)
    N[i] = buf_to_double(BUF.hdr_old, 12)
    fastskip(sio, 12)
  end
  ND = stream_int(sio, 4)
  D = Array{Float64,1}(undef, ND)
  @inbounds for i = 1:ND
    fast_readbytes!(sio, BUF.hdr_old, 12)
    D[i] = buf_to_double(BUF.hdr_old, 12)
    fastskip(sio, 12)
  end
  close(sio)

  if v > 1
    println(", STAGE ", stage)
    if v > 2
      println(" "^16, "units in code = ", uic)
      println(" "^16, "units out code = ", uoc)
      println(" "^16, "NN = ", NN, ":")
      @inbounds for i = 1:NN
        println(" "^16, N[i])
      end
      println(" "^16, "ND = ", ND, ":")
      @inbounds for i = 1:ND
        println(" "^16, D[i])
      end
    end
  end

  if stage > length(R.fs)
    append!(R, MultiStageResp(6))
  end
  R.stage[stage] = CoeffResp(N, D)
  if units
    @inbounds R.i[stage] = fix_units(units_lookup[uic])
    @inbounds R.o[stage] = fix_units(units_lookup[uoc])
  end
  return stage
end

function blk_057(io::IO, nb::Int64, v::Integer, R::MultiStageResp)
  sio = blk_string_read(io, nb, v)
  stage = stream_int(sio, 2)
  fast_readbytes!(sio, BUF.hdr_old, 10)
  fs = buf_to_double(BUF.hdr_old, 10)
  fac = stream_int(sio, 5)
  os = stream_int(sio, 5)
  fast_readbytes!(sio, BUF.hdr_old, 11)
  delay = buf_to_double(BUF.hdr_old, 11)
  fast_readbytes!(sio, BUF.hdr_old, 11)
  corr = buf_to_double(BUF.hdr_old, 11)
  close(sio)

  if v > 1
    println(", STAGE #", stage)
    if v > 2
      println(" "^16, "fs = ", fs)
      println(" "^16, "decimation factor = ", fac)
      println(" "^16, "decimation offset = ", os)
      println(" "^16, "delay = ", delay)
      println(" "^16, "delay correction applied = ", corr)
    end
  end

  @inbounds R.fs[stage] = fs
  @inbounds R.delay[stage] = delay
  @inbounds R.corr[stage] = corr
  @inbounds R.fac[stage] = fac
  @inbounds R.os[stage] = os
  return stage
end

function blk_058(io::IO, nb::Int64, v::Integer, C::SeisChannel)
  sio = blk_string_read(io, nb, v)
  stage = stream_int(sio, 2)
  fast_readbytes!(sio, BUF.hdr_old, 12)
  if stage == 0
    C.gain = buf_to_double(BUF.hdr_old, 12)
    close(sio)
    if v > 1
      println(", STAGE #", stage)
    end
    return stage
  else
    C.resp.gain[stage] = buf_to_double(BUF.hdr_old, 12)
    fast_readbytes!(sio, BUF.hdr_old, 12)
    C.resp.fg[stage] = buf_to_double(BUF.hdr_old, 12)

    if v > 1
      println(", STAGE #", stage)
      if v > 2
        println(" "^16, "gain = ", C.resp.gain[stage], " (f = ", C.resp.fg[stage], " Hz)")
      end
    end

    # station history is not in the scope of SeisIO
    close(sio)
    return stage
  end
end

# Not in scope of SeisIO
function blk_059!(io::IO, nb::Int64, v::Integer, C::SeisChannel, units::Bool)
  sio = blk_string_read(io, nb, v)
  v > 1 && println("")
  if units
    ts = parse_resp_date(sio, BUF.uint16_buf)
    te = parse_resp_date(sio, BUF.uint16_buf)
    if te == -56504908800000000
      te = 19880899199000000
    end
    k = stream_int(sio, 4)
    t_str = string("comment,", ts, ",", te, ",", comments[k])
    if v > 2
      println(u2d(div(ts, 1000000)), "–", u2d(div(te, 1000000)), ": ", comments[k])
    end
    note!(C, t_str)
  end
  close(sio)
  return nothing
end

# Assign dictionary elements with response info in the 41-49 blockettes
function blk_060(io::IO, nb::Int64, v::Integer, R::MultiStageResp)
  sio = blk_string_read(io, nb, v)

  nstg = stream_int(sio, 2)
  if v > 1
    println(", # STAGES = ", nstg)
  end

  for i = 1:nstg
    seq = stream_int(sio, 2)
    nr = stream_int(sio, 2)
    for j = 1:nr
      k = stream_int(sio, 4)
      rr = get(responses, k, "")
      if v > 2
        printstyled(" "^16, " assigning response to stage #", seq, ":\n", color=:green, bold=true)
        println(rr)
      end

      if isa(rr, Blk48)
        @inbounds R.gain[seq] = rr.gain
        @inbounds R.fg[seq] = rr.fg
      elseif isa(rr, Blk47)
        @inbounds R.fs[seq] = rr.fs
        @inbounds R.delay[seq] = rr.delay
        @inbounds R.corr[seq] = rr.corr
        @inbounds R.fac[seq] = rr.fac
        @inbounds R.os[seq] = rr.os
      elseif typeof(rr[1]) <: InstrumentResponse
        (length(R.fs) < seq) && append!(R, MultiStageResp(6))
        @inbounds R.stage[seq] = rr[1]
        @inbounds R.i[seq] = rr[2]
        @inbounds R.o[seq] = rr[3]
      end
    end
  end
  close(sio)
  return nstg
end

function blk_061(io::IO, nb::Int64, v::Integer, R::MultiStageResp, units::Bool)
  sio = blk_string_read(io, nb, v)

  stage = stream_int(sio, 2)
  skip_string!(sio)
  symm_code = Char(fastread(sio))
  uic = stream_int(sio, 3)
  uoc = stream_int(sio, 3)
  NF = stream_int(sio, 4)
  F = Array{Float64,1}(undef, NF)
  @inbounds for i = 1:NF
    fast_readbytes!(sio, BUF.hdr_old, 14)
    F[i] = buf_to_double(BUF.hdr_old, 14)
  end
  close(sio)

  if v > 1
    println(", STAGE #", stage, " , SYMM = ", Char(symm_code))
    if v > 2
      println(" "^16, "units in code #", uic)
      println(" "^16, "units out code #", uoc)
      println(" "^16, "NF = ", NF)
      @inbounds for i = 1:NF
        println(" "^16, F[i])
      end
    end
  end

  (length(R.fs) < stage) && append!(R, MultiStageResp(6))
  @inbounds R.stage[stage] = CoeffResp(F, Float64[])
  if units
    @inbounds R.i[stage] = fix_units(units_lookup[uic])
    @inbounds R.o[stage] = fix_units(units_lookup[uoc])
  end
  return stage
end


text:
using DiffEqFlux
using Pkg
using SafeTestsets
using Test

const GROUP = get(ENV, "GROUP", "All")
const is_APPVEYOR = ( Sys.iswindows() && haskey(ENV,"APPVEYOR") )
const is_CI = haskey(ENV,"CI")

function activate_gpu_env()
    Pkg.activate("gpu")
    Pkg.develop(PackageSpec(path=dirname(@__DIR__)))
    Pkg.instantiate()
end

@time begin
if GROUP == "All" || GROUP == "DiffEqFlux" || GROUP == "Layers"
    @safetestset "Layers Tests" begin include("layers.jl") end
    @safetestset "Fast Layers" begin include("fast_layers.jl") end
    @safetestset "Layers SciML Tests" begin include("layers_sciml.jl") end
    @safetestset "Layers SDE" begin include("layers_sde.jl") end
    @safetestset "Layers DDE" begin include("layers_dde.jl") end
    @safetestset "Collocation Regression" begin include("collocation_regression.jl") end
end

if GROUP == "All" || GROUP == "DiffEqFlux" || GROUP == "BasicNeuralDE"
    @safetestset "Neural DE Tests" begin include("neural_de.jl") end
    @safetestset "Augmented Neural DE Tests" begin include("augmented_nde.jl") end
    @safetestset "Neural Graph DE" begin include("neural_gde.jl") end
    @safetestset "Neural ODE MM Tests" begin include("neural_ode_mm.jl") end
    @safetestset "Fast Neural ODE Tests" begin include("fast_neural_ode.jl") end
    @safetestset "Tensor Product Layer" begin include("tensor_product_test.jl") end
    @safetestset "Spline Layer" begin include("spline_layer_test.jl") end
end

if GROUP == "All" || GROUP == "AdvancedNeuralDE"
    @safetestset "CNF Layer Tests" begin include("cnf_test.jl") end
    @safetestset "Newton Neural ODE Tests" begin include("newton_neural_ode.jl") end
    @safetestset "Neural Second Order ODE Tests" begin include("second_order_ode.jl") end
end

if GROUP == "All" || GROUP == "Integration"
    @safetestset "Partial Neural Tests" begin include("partial_neural.jl") end
    @safetestset "Size Handling in Adjoint Tests" begin include("size_handling_adjoint.jl") end
    @safetestset "odenet" begin include("odenet.jl") end
    @safetestset "GDP Regression Tests" begin include("gdp_regression_test.jl") end
    @safetestset "Stiff Nested AD Tests" begin include("stiff_nested_ad.jl") end
end
end

if !is_APPVEYOR && GROUP == "GPU"
    if is_CI
        activate_gpu_env()
    end
    @safetestset "odenet GPU" begin include("odenet_gpu.jl") end
    @safetestset "Neural DE GPU Tests" begin include("neural_de_gpu.jl") end
end

ext: jl


text:
function compute_residual_AH!(res::Array, sigma::Array,
                              gauge::Gauge, cache::HorizonCache,
                              sys::System{Outer})
    _, Nx, Ny = size(sys)

    Dx  = sys.Dx
    Dxx = sys.Dxx
    Dy  = sys.Dy
    Dyy = sys.Dyy

    # we assume here that the bulk functions and their u-derivatives have
    # already been interpolated to the uAH surface (u = 1/sigma)

    B1_uAH      = cache.bulkhorizon.B1_uAH
    B2_uAH      = cache.bulkhorizon.B2_uAH
    G_uAH       = cache.bulkhorizon.G_uAH
    S_uAH       = cache.bulkhorizon.S_uAH
    Fx_uAH      = cache.bulkhorizon.Fx_uAH
    Fy_uAH      = cache.bulkhorizon.Fy_uAH
    Sd_uAH      = cache.bulkhorizon.Sd_uAH

    Du_B1_uAH   = cache.bulkhorizon.Du_B1_uAH
    Du_B2_uAH   = cache.bulkhorizon.Du_B2_uAH
    Du_G_uAH    = cache.bulkhorizon.Du_G_uAH
    Du_S_uAH    = cache.bulkhorizon.Du_S_uAH
    Du_Fx_uAH   = cache.bulkhorizon.Du_Fx_uAH
    Du_Fy_uAH   = cache.bulkhorizon.Du_Fy_uAH
    Du_Sd_uAH   = cache.bulkhorizon.Du_Sd_uAH

    Duu_B1_uAH  = cache.bulkhorizon.Duu_B1_uAH
    Duu_B2_uAH  = cache.bulkhorizon.Duu_B2_uAH
    Duu_G_uAH   = cache.bulkhorizon.Duu_G_uAH
    Duu_S_uAH   = cache.bulkhorizon.Duu_S_uAH
    Duu_Fx_uAH  = cache.bulkhorizon.Duu_Fx_uAH
    Duu_Fy_uAH  = cache.bulkhorizon.Duu_Fy_uAH

    # compute AH equation residual
    @fastmath @inbounds Threads.@threads for j in 1:Ny
        @inbounds for i in 1:Nx
            uAH = 1 / sigma[1,i,j]
            u2  = uAH * uAH
            u3  = uAH * uAH * uAH
            u4  = uAH * uAH * uAH * uAH

            xi    = gauge.xi[1,i,j]
            xi_x  = Dx(gauge.xi, 1,i,j)
            xi_y  = Dy(gauge.xi, 1,i,j)
            xi_xx = Dxx(gauge.xi, 1,i,j)
            xi_yy = Dyy(gauge.xi, 1,i,j)
            xi_xy = Dx(Dy, gauge.xi, 1,i,j)

            sigma0    = sigma[1,i,j]
            sigma0_x  = Dx(sigma, 1,i,j)
            sigma0_y  = Dy(sigma, 1,i,j)
            sigma0_xx = Dxx(sigma, 1,i,j)
            sigma0_yy = Dyy(sigma, 1,i,j)
            sigma0_xy = Dx(Dy, sigma, 1,i,j)

            B1    = B1_uAH[1,i,j]
            B2    = B2_uAH[1,i,j]
            G     = G_uAH[1,i,j]
            S     = S_uAH[1,i,j]
            Fx    = Fx_uAH[1,i,j]
            Fy    = Fy_uAH[1,i,j]
            Sd    = Sd_uAH[1,i,j]

            # r derivatives

            B1p        = -u2 * Du_B1_uAH[1,i,j]
            B2p        = -u2 * Du_B2_uAH[1,i,j]
            Gp         = -u2 * Du_G_uAH[1,i,j]
            Sp         = -u2 * Du_S_uAH[1,i,j]
            Fxp        = -u2 * Du_Fx_uAH[1,i,j]
            Fyp        = -u2 * Du_Fy_uAH[1,i,j]
            Sdp        = -u2 * Du_Sd_uAH[1,i,j]

            B1pp       = 2*u3 * Du_B1_uAH[1,i,j]  + u4 * Duu_B1_uAH[1,i,j]
            B2pp       = 2*u3 * Du_B2_uAH[1,i,j]  + u4 * Duu_B2_uAH[1,i,j]
            Gpp        = 2*u3 * Du_G_uAH[1,i,j]   + u4 * Duu_G_uAH[1,i,j]
            Spp        = 2*u3 * Du_S_uAH[1,i,j]   + u4 * Duu_S_uAH[1,i,j]
            Fxpp       = 2*u3 * Du_Fx_uAH[1,i,j]  + u4 * Duu_Fx_uAH[1,i,j]
            Fypp       = 2*u3 * Du_Fy_uAH[1,i,j]  + u4 * Duu_Fy_uAH[1,i,j]

            # x and y derivatives

            B1_x    = Dx(B1_uAH,   1,i,j) - B1p * sigma0_x
            B2_x    = Dx(B2_uAH,   1,i,j) - B2p * sigma0_x
            G_x     = Dx(G_uAH,    1,i,j) -  Gp * sigma0_x
            S_x     = Dx(S_uAH,    1,i,j) -  Sp * sigma0_x
            Fx_x    = Dx(Fx_uAH,   1,i,j) - Fxp * sigma0_x
            Fy_x    = Dx(Fy_uAH,   1,i,j) - Fyp * sigma0_x
            Sd_x    = Dx(Sd_uAH,   1,i,j) - Sdp * sigma0_x

            B1_y    = Dy(B1_uAH,   1,i,j) - B1p * sigma0_y
            B2_y    = Dy(B2_uAH,   1,i,j) - B2p * sigma0_y
            G_y     = Dy(G_uAH,    1,i,j) -  Gp * sigma0_y
            S_y     = Dy(S_uAH,    1,i,j) -  Sp * sigma0_y
            Fx_y    = Dy(Fx_uAH,   1,i,j) - Fxp * sigma0_y
            Fy_y    = Dy(Fy_uAH,   1,i,j) - Fyp * sigma0_y
            Sd_y    = Dy(Sd_uAH,   1,i,j) - Sdp * sigma0_y

            B1p_x   = -u2 * Dx(Du_B1_uAH, 1,i,j) - B1pp * sigma0_x + 2 * u3 * sigma0_x * Du_B1_uAH[1,i,j]
            B2p_x   = -u2 * Dx(Du_B2_uAH, 1,i,j) - B2pp * sigma0_x + 2 * u3 * sigma0_x * Du_B2_uAH[1,i,j]
            Gp_x    = -u2 * Dx(Du_G_uAH,  1,i,j) -  Gpp * sigma0_x + 2 * u3 * sigma0_x *  Du_G_uAH[1,i,j]
            Sp_x    = -u2 * Dx(Du_S_uAH,  1,i,j) -  Spp * sigma0_x + 2 * u3 * sigma0_x *  Du_S_uAH[1,i,j]
            Fxp_x   = -u2 * Dx(Du_Fx_uAH, 1,i,j) - Fxpp * sigma0_x + 2 * u3 * sigma0_x * Du_Fx_uAH[1,i,j]
            Fyp_x   = -u2 * Dx(Du_Fy_uAH, 1,i,j) - Fypp * sigma0_x + 2 * u3 * sigma0_x * Du_Fy_uAH[1,i,j]

            B1p_y   = -u2 * Dy(Du_B1_uAH, 1,i,j) - B1pp * sigma0_y + 2 * u3 * sigma0_y * Du_B1_uAH[1,i,j]
            B2p_y   = -u2 * Dy(Du_B2_uAH, 1,i,j) - B2pp * sigma0_y + 2 * u3 * sigma0_y * Du_B2_uAH[1,i,j]
            Gp_y    = -u2 * Dy(Du_G_uAH,  1,i,j) -  Gpp * sigma0_y + 2 * u3 * sigma0_y *  Du_G_uAH[1,i,j]
            Sp_y    = -u2 * Dy(Du_S_uAH,  1,i,j) -  Spp * sigma0_y + 2 * u3 * sigma0_y *  Du_S_uAH[1,i,j]
            Fxp_y   = -u2 * Dy(Du_Fx_uAH, 1,i,j) - Fxpp * sigma0_y + 2 * u3 * sigma0_y * Du_Fx_uAH[1,i,j]
            Fyp_y   = -u2 * Dy(Du_Fy_uAH, 1,i,j) - Fypp * sigma0_y + 2 * u3 * sigma0_y * Du_Fy_uAH[1,i,j]

            vars = (
                sigma0, sigma0_x, sigma0_y, sigma0_xx, sigma0_yy, sigma0_xy,
                xi    , xi_x    , xi_y    , xi_xx    , xi_yy    , xi_xy,
                B1   , B2   , G   ,  S    , Fx    , Fy    , Sd ,
                B1p  , B2p  , Gp  ,  Sp   , Fxp   , Fyp   , Sdp,
                B1pp , B2pp , Gpp ,  Spp  , Fxpp  , Fypp  ,
                B1_x , B2_x , G_x ,  S_x  , Fx_x  , Fy_x  , Sd_x,
                B1_y , B2_y , G_y ,  S_y  , Fx_y  , Fy_y  , Sd_y,
            )

            res[1,i,j] = AH_eq_res(vars, sys.gridtype)
        end
    end

    res
end

function compute_coeffs_AH!(sigma::Array, gauge::Gauge, cache::HorizonCache,
                            sys::System{Outer})
    _, Nx, Ny = size(sys)

    Dx  = sys.Dx
    Dxx = sys.Dxx
    Dy  = sys.Dy
    Dyy = sys.Dyy

    B1_uAH      = cache.bulkhorizon.B1_uAH
    B2_uAH      = cache.bulkhorizon.B2_uAH
    G_uAH       = cache.bulkhorizon.G_uAH
    S_uAH       = cache.bulkhorizon.S_uAH
    Fx_uAH      = cache.bulkhorizon.Fx_uAH
    Fy_uAH      = cache.bulkhorizon.Fy_uAH
    Sd_uAH      = cache.bulkhorizon.Sd_uAH

    Du_B1_uAH   = cache.bulkhorizon.Du_B1_uAH
    Du_B2_uAH   = cache.bulkhorizon.Du_B2_uAH
    Du_G_uAH    = cache.bulkhorizon.Du_G_uAH
    Du_S_uAH    = cache.bulkhorizon.Du_S_uAH
    Du_Fx_uAH   = cache.bulkhorizon.Du_Fx_uAH
    Du_Fy_uAH   = cache.bulkhorizon.Du_Fy_uAH
    Du_Sd_uAH   = cache.bulkhorizon.Du_Sd_uAH

    Duu_B1_uAH  = cache.bulkhorizon.Duu_B1_uAH
    Duu_B2_uAH  = cache.bulkhorizon.Duu_B2_uAH
    Duu_G_uAH   = cache.bulkhorizon.Duu_G_uAH
    Duu_S_uAH   = cache.bulkhorizon.Duu_S_uAH
    Duu_Fx_uAH  = cache.bulkhorizon.Duu_Fx_uAH
    Duu_Fy_uAH  = cache.bulkhorizon.Duu_Fy_uAH

    axx         = cache.axx
    ayy         = cache.ayy
    axy         = cache.axy
    bx          = cache.bx
    by          = cache.by
    cc          = cache.cc
    #b_vec       = cache.b_vec

    ind2D  = LinearIndices(B1_uAH[1,:,:])

    # coefficients of the derivative operators
    @fastmath @inbounds Threads.@threads for j in 1:Ny
        @inbounds for i in 1:Nx
            idx   = ind2D[i,j]

            uAH   = 1 / sigma[1,i,j]
            u2    = uAH * uAH
            u3    = uAH * uAH * uAH
            u4    = uAH * uAH * uAH * uAH

            xi    = gauge.xi[1,i,j]
            xi_x  = Dx(gauge.xi, 1,i,j)
            xi_y  = Dy(gauge.xi, 1,i,j)
            xi_xx = Dxx(gauge.xi, 1,i,j)
            xi_yy = Dyy(gauge.xi, 1,i,j)
            xi_xy = Dx(Dy, gauge.xi, 1,i,j)

            sigma0    = sigma[1,i,j]
            sigma0_x  = Dx(sigma, 1,i,j)
            sigma0_y  = Dy(sigma, 1,i,j)
            sigma0_xx = Dxx(sigma, 1,i,j)
            sigma0_yy = Dyy(sigma, 1,i,j)
            sigma0_xy = Dx(Dy, sigma, 1,i,j)

            B1    = B1_uAH[1,i,j]
            B2    = B2_uAH[1,i,j]
            G     = G_uAH[1,i,j]
            S     = S_uAH[1,i,j]
            Fx    = Fx_uAH[1,i,j]
            Fy    = Fy_uAH[1,i,j]
            Sd    = Sd_uAH[1,i,j]

            # r derivatives

            B1p        = -u2 * Du_B1_uAH[1,i,j]
            B2p        = -u2 * Du_B2_uAH[1,i,j]
            Gp         = -u2 * Du_G_uAH[1,i,j]
            Sp         = -u2 * Du_S_uAH[1,i,j]
            Fxp        = -u2 * Du_Fx_uAH[1,i,j]
            Fyp        = -u2 * Du_Fy_uAH[1,i,j]
            Sdp        = -u2 * Du_Sd_uAH[1,i,j]

            B1pp       = 2*u3 * Du_B1_uAH[1,i,j]  + u4 * Duu_B1_uAH[1,i,j]
            B2pp       = 2*u3 * Du_B2_uAH[1,i,j]  + u4 * Duu_B2_uAH[1,i,j]
            Gpp        = 2*u3 * Du_G_uAH[1,i,j]   + u4 * Duu_G_uAH[1,i,j]
            Spp        = 2*u3 * Du_S_uAH[1,i,j]   + u4 * Duu_S_uAH[1,i,j]
            Fxpp       = 2*u3 * Du_Fx_uAH[1,i,j]  + u4 * Duu_Fx_uAH[1,i,j]
            Fypp       = 2*u3 * Du_Fy_uAH[1,i,j]  + u4 * Duu_Fy_uAH[1,i,j]

            # x and y derivatives

            B1_x    = Dx(B1_uAH,   1,i,j) - B1p * sigma0_x
            B2_x    = Dx(B2_uAH,   1,i,j) - B2p * sigma0_x
            G_x     = Dx(G_uAH,    1,i,j) -  Gp * sigma0_x
            S_x     = Dx(S_uAH,    1,i,j) -  Sp * sigma0_x
            Fx_x    = Dx(Fx_uAH,   1,i,j) - Fxp * sigma0_x
            Fy_x    = Dx(Fy_uAH,   1,i,j) - Fyp * sigma0_x
            Sd_x    = Dx(Sd_uAH,   1,i,j) - Sdp * sigma0_x

            B1_y    = Dy(B1_uAH,   1,i,j) - B1p * sigma0_y
            B2_y    = Dy(B2_uAH,   1,i,j) - B2p * sigma0_y
            G_y     = Dy(G_uAH,    1,i,j) -  Gp * sigma0_y
            S_y     = Dy(S_uAH,    1,i,j) -  Sp * sigma0_y
            Fx_y    = Dy(Fx_uAH,   1,i,j) - Fxp * sigma0_y
            Fy_y    = Dy(Fy_uAH,   1,i,j) - Fyp * sigma0_y
            Sd_y    = Dy(Sd_uAH,   1,i,j) - Sdp * sigma0_y

            B1p_x   = -u2 * Dx(Du_B1_uAH, 1,i,j) - B1pp * sigma0_x + 2 * u3 * sigma0_x * Du_B1_uAH[1,i,j]
            B2p_x   = -u2 * Dx(Du_B2_uAH, 1,i,j) - B2pp * sigma0_x + 2 * u3 * sigma0_x * Du_B2_uAH[1,i,j]
            Gp_x    = -u2 * Dx(Du_G_uAH,  1,i,j) -  Gpp * sigma0_x + 2 * u3 * sigma0_x *  Du_G_uAH[1,i,j]
            Sp_x    = -u2 * Dx(Du_S_uAH,  1,i,j) -  Spp * sigma0_x + 2 * u3 * sigma0_x *  Du_S_uAH[1,i,j]
            Fxp_x   = -u2 * Dx(Du_Fx_uAH, 1,i,j) - Fxpp * sigma0_x + 2 * u3 * sigma0_x * Du_Fx_uAH[1,i,j]
            Fyp_x   = -u2 * Dx(Du_Fy_uAH, 1,i,j) - Fypp * sigma0_x + 2 * u3 * sigma0_x * Du_Fy_uAH[1,i,j]

            B1p_y   = -u2 * Dy(Du_B1_uAH, 1,i,j) - B1pp * sigma0_y + 2 * u3 * sigma0_y * Du_B1_uAH[1,i,j]
            B2p_y   = -u2 * Dy(Du_B2_uAH, 1,i,j) - B2pp * sigma0_y + 2 * u3 * sigma0_y * Du_B2_uAH[1,i,j]
            Gp_y    = -u2 * Dy(Du_G_uAH,  1,i,j) -  Gpp * sigma0_y + 2 * u3 * sigma0_y *  Du_G_uAH[1,i,j]
            Sp_y    = -u2 * Dy(Du_S_uAH,  1,i,j) -  Spp * sigma0_y + 2 * u3 * sigma0_y *  Du_S_uAH[1,i,j]
            Fxp_y   = -u2 * Dy(Du_Fx_uAH, 1,i,j) - Fxpp * sigma0_y + 2 * u3 * sigma0_y * Du_Fx_uAH[1,i,j]
            Fyp_y   = -u2 * Dy(Du_Fy_uAH, 1,i,j) - Fypp * sigma0_y + 2 * u3 * sigma0_y * Du_Fy_uAH[1,i,j]

            vars = (
                sigma0, sigma0_x, sigma0_y, sigma0_xx, sigma0_yy, sigma0_xy,
                xi    , xi_x    , xi_y    , xi_xx    , xi_yy    , xi_xy,
                B1   , B2   , G   ,  S    , Fx    , Fy    , Sd ,
                B1p  , B2p  , Gp  ,  Sp   , Fxp   , Fyp   , Sdp,
                B1pp , B2pp , Gpp ,  Spp  , Fxpp  , Fypp  ,
                B1_x , B2_x , G_x ,  S_x  , Fx_x  , Fy_x  , Sd_x,
	        B1_y , B2_y , G_y ,  S_y  , Fx_y  , Fy_y  , Sd_y,
                B1p_x, B2p_x, Gp_x,  Sp_x , Fxp_x , Fyp_x ,
                B1p_y, B2p_y, Gp_y,  Sp_y , Fxp_y , Fyp_y ,
            )

            a11, a22, a12, b1, b2, c = AH_eq_coeff(vars, sys.gridtype)

            axx[idx]   = a11
            ayy[idx]   = a22
            axy[idx]   = a12
            bx[idx]    = b1
            by[idx]    = b2
            cc[idx]    = c
        end
    end

    nothing
end

#= Finding the Apparent Horizon

this is a 2D (non-linear) PDE of the type

  αxx f_xx + αyy f_yy + αxy f_xy + βxx (f_x)^2 + βyy (f_y)^2 + βxy (f_x) (f_y)
   + γx f_x + γy f_y + S = 0

we solve this equation with a Newton-Kantorovich method where, starting with a
guess, we solve the associated linear problem (using the same strategy as in
compute_xi_t!) to improve our guess until a sufficiently precise solution is
reached.
=#
function find_AH!(sigma::Array, bulkconstrain::BulkConstrained,
                  bulkevol::BulkEvolved, deriv::BulkDeriv, gauge::Gauge,
                  cache::HorizonCache, sys::System{Outer}, ahf::AHF)
    _, Nx, Ny = size(sys)
    bulk = Bulk(bulkevol, bulkconstrain)

    Du  = sys.Du
    Duu = sys.Duu
    Dx  = sys.Dx
    Dxx = sys.Dxx
    Dy  = sys.Dy
    Dyy = sys.Dyy

    interp = sys.uinterp

    B1_uAH      = cache.bulkhorizon.B1_uAH
    B2_uAH      = cache.bulkhorizon.B2_uAH
    G_uAH       = cache.bulkhorizon.G_uAH
    S_uAH       = cache.bulkhorizon.S_uAH
    Fx_uAH      = cache.bulkhorizon.Fx_uAH
    Fy_uAH      = cache.bulkhorizon.Fy_uAH
    Sd_uAH      = cache.bulkhorizon.Sd_uAH

    Du_B1_uAH   = cache.bulkhorizon.Du_B1_uAH
    Du_B2_uAH   = cache.bulkhorizon.Du_B2_uAH
    Du_G_uAH    = cache.bulkhorizon.Du_G_uAH
    Du_S_uAH    = cache.bulkhorizon.Du_S_uAH
    Du_Fx_uAH   = cache.bulkhorizon.Du_Fx_uAH
    Du_Fy_uAH   = cache.bulkhorizon.Du_Fy_uAH
    Du_Sd_uAH   = cache.bulkhorizon.Du_Sd_uAH

    Duu_B1_uAH  = cache.bulkhorizon.Duu_B1_uAH
    Duu_B2_uAH  = cache.bulkhorizon.Duu_B2_uAH
    Duu_G_uAH   = cache.bulkhorizon.Duu_G_uAH
    Duu_S_uAH   = cache.bulkhorizon.Duu_S_uAH
    Duu_Fx_uAH  = cache.bulkhorizon.Duu_Fx_uAH
    Duu_Fy_uAH  = cache.bulkhorizon.Duu_Fy_uAH

    axx         = cache.axx
    ayy         = cache.ayy
    axy         = cache.axy
    bx          = cache.bx
    by          = cache.by
    cc          = cache.cc
    b_vec       = cache.b_vec

    Dx_2D       = cache.Dx_2D
    Dxx_2D      = cache.Dxx_2D
    Dy_2D       = cache.Dy_2D
    Dyy_2D      = cache.Dyy_2D
    Dxy_2D      = cache.Dxy_2D
    _Dx_2D      = cache._Dx_2D
    _Dxx_2D     = cache._Dxx_2D
    _Dy_2D      = cache._Dy_2D
    _Dyy_2D     = cache._Dyy_2D
    _Dxy_2D     = cache._Dxy_2D

    # take all u-derivatives. we assume that they've not been computed beforehand
    @sync begin
        @spawn mul!(deriv.Du_B1,  Du,  bulk.B1)
        @spawn mul!(deriv.Du_B2,  Du,  bulk.B2)
        @spawn mul!(deriv.Du_G,   Du,  bulk.G)
        @spawn mul!(deriv.Du_S,   Du,  bulk.S)
        @spawn mul!(deriv.Du_Fx,  Du,  bulk.Fx)
        @spawn mul!(deriv.Du_Fy,  Du,  bulk.Fy)
        @spawn mul!(deriv.Du_Sd,  Du,  bulk.Sd)

        @spawn mul!(deriv.Duu_B1, Duu, bulk.B1)
        @spawn mul!(deriv.Duu_B2, Duu, bulk.B2)
        @spawn mul!(deriv.Duu_G,  Duu, bulk.G)
        @spawn mul!(deriv.Duu_S,  Duu, bulk.S)
        @spawn mul!(deriv.Duu_Fx, Duu, bulk.Fx)
        @spawn mul!(deriv.Duu_Fy, Duu, bulk.Fy)
        # @spawn mul!(deriv.Duu_Sd, Duu, bulk.Sd)
    end

    res = similar(sigma)
    f0  = similar(b_vec)

    println("INFO (AH): Looking for the apparent horizon...")
    println("    it \t max_res")
    # start relaxation method
    it = 0
    while true

        # interpolate bulk functions (and u-derivatives) to the 1/u = r = sigma surface
        @inbounds Threads.@threads for j in 1:Ny
            @inbounds for i in 1:Nx
                uAH = 1 / sigma[1,i,j]
                u2  = uAH * uAH
                u3  = uAH * uAH * uAH
                u4  = uAH * uAH * uAH * uAH

                B1_uAH[1,i,j]       = interp(view(bulk.B1,  :,i,j))(uAH)
                B2_uAH[1,i,j]       = interp(view(bulk.B2,  :,i,j))(uAH)
                G_uAH[1,i,j]        = interp(view(bulk.G,   :,i,j))(uAH)
                S_uAH[1,i,j]        = interp(view(bulk.S,   :,i,j))(uAH)
                Fx_uAH[1,i,j]       = interp(view(bulk.Fx,  :,i,j))(uAH)
                Fy_uAH[1,i,j]       = interp(view(bulk.Fy,  :,i,j))(uAH)
                Sd_uAH[1,i,j]       = interp(view(bulk.Sd,  :,i,j))(uAH)

                Du_B1_uAH[1,i,j]    = interp(view(deriv.Du_B1,  :,i,j))(uAH)
                Du_B2_uAH[1,i,j]    = interp(view(deriv.Du_B2,  :,i,j))(uAH)
                Du_G_uAH[1,i,j]     = interp(view(deriv.Du_G,   :,i,j))(uAH)
                Du_S_uAH[1,i,j]     = interp(view(deriv.Du_S,   :,i,j))(uAH)
                Du_Fx_uAH[1,i,j]    = interp(view(deriv.Du_Fx,  :,i,j))(uAH)
                Du_Fy_uAH[1,i,j]    = interp(view(deriv.Du_Fy,  :,i,j))(uAH)
                Du_Sd_uAH[1,i,j]    = interp(view(deriv.Du_Sd,  :,i,j))(uAH)

                Duu_B1_uAH[1,i,j]   = interp(view(deriv.Duu_B1,  :,i,j))(uAH)
                Duu_B2_uAH[1,i,j]   = interp(view(deriv.Duu_B2,  :,i,j))(uAH)
                Duu_G_uAH[1,i,j]    = interp(view(deriv.Duu_G,   :,i,j))(uAH)
                Duu_S_uAH[1,i,j]    = interp(view(deriv.Duu_S,   :,i,j))(uAH)
                Duu_Fx_uAH[1,i,j]   = interp(view(deriv.Duu_Fx,  :,i,j))(uAH)
                Duu_Fy_uAH[1,i,j]   = interp(view(deriv.Duu_Fy,  :,i,j))(uAH)
            end
        end

        compute_residual_AH!(res, sigma, gauge, cache, sys)
        max_res = maximum(abs.(res))
        println("    $it \t $max_res")

        if max_res < ahf.epsilon
            break
        end

        # compute axx, ayy, axy, bx, by and cc coefficients of the linearized
        # equation (stored in the cache struct)
        compute_coeffs_AH!(sigma, gauge, cache, sys)
        copyto!(cache.b_vec, res)

        if it >= ahf.itmax
            println("INFO (AH): maximum iteration reached.")
            println("INFO (AH): giving up")
            break
        end

        axxId = Diagonal(axx)
        ayyId = Diagonal(ayy)
        axyId = Diagonal(axy)
        bxId  = Diagonal(bx)
        byId  = Diagonal(by)
        ccId  = Diagonal(cc)

        # build the differential operators by multiplying with the coefficients
        # computed in the loop above (note that Dxx_2D, Dyy_2D, etc, are never
        # overwritten)
        mul!(_Dxx_2D, axxId, Dxx_2D)
        mul!(_Dyy_2D, ayyId, Dyy_2D)
        mul!(_Dxy_2D, axyId, Dxy_2D)
        mul!(_Dx_2D,  bxId,  Dx_2D)
        mul!(_Dy_2D,  byId,  Dy_2D)

        # build actual operator to be inverted
        A_mat = _Dxx_2D + _Dyy_2D + _Dxy_2D + _Dx_2D + _Dy_2D + ccId

        # solve system
        A_fact = lu(A_mat)
        ldiv!(f0, A_fact, b_vec)

        # update solution
        @inbounds for idx in eachindex(sigma)
            sigma[idx] -= f0[idx]
        end
        it += 1

        # check if sigma inside domain
        min_uAH = 1/maximum(sigma)
        max_uAH = 1/minimum(sigma)
        if ( min_uAH < sys.ucoord[1] || max_uAH > sys.ucoord[end] )
            println("INFO (AH): guess outside domain, min_uAH = $min_uAH, max_uAH = $max_uAH")
            println("INFO (AH): giving up")
            break
        end

    end # end relaxation loop

    nothing
end


text:
ENV["PYTHON"]="/usr/bin/python"
if Pkg.installed("PyPlot")==nothing
  Pkg.clone("PyPlot")
end
if isdir(Pkg.dir("jInv"))==false
  Pkg.clone("https://github.com/JuliaInv/jInv.jl","jInv")
end
if isdir(Pkg.dir("FactoredEikonalFastMarching"))==false
  Pkg.clone("https://github.com/JuliaInv/FactoredEikonalFastMarching.jl","FactoredEikonalFastMarching")
end
#Pkg.test("FactoredEikonalFastMarching")


text:
# sh test_benchmark.sh
using PyPlot
using DelimitedFiles

r = readdlm("data/benchmark.txt")
idx = sortperm(r[:,1])
r = r[idx,:]
figure(figsize=(10,4))
subplot(121)
bar(1:size(r,1), r[:,2], tick_label = Int64.(r[:,1]))
xlabel("Number of Processors")
ylabel("Time (seconds)")

subplot(122)
plot(r[:,1], r[1,2]./r[:,2], label="Speedup")
legend()
xlabel("Number of Processors")

savefig("benchmark.png")


text:
using AcFDTD
using Base.Test
using Base.Profile
srand(123)

geo = CuboidRoom(10,15,16,[10.;20.;30.;40.;50.;60.], FDTDEnv(0.1,IWB()))
Nt = 200
s_t = [randn(Nt,1) randn(Nt,1)]
S = randn(size(geo)...,Nt)
xs = [(5, 5, 5),(5,2,4)]
xr = [(5, 8, 7),(2,2,2)]
f = FDTD(geo)

p0ic = randn(10,15,16)
p1ic = randn(10,15,16)

println("test calls: \n")
println("----------      no IC")
p    = fdtd(xr, xs, Nt,geo, s_t)
p    = fdtd(xr, xs, Nt,  f, s_t)
p2 = copy(p)
fdtd!(p2, xr, xs, Nt,  f, s_t)
@time p    = fdtd(xr, xs, Nt,  f, s_t)
println("                       in-place call")
@time        fdtd!(p2,xr, xs, Nt,  f, s_t)

@test norm(p-p2) < 1e-8

println("----------         IC")
pic  = fdtd(xr, xs, Nt,geo, s_t, p0ic, p1ic)
ppic = fdtd(xr, xs, Nt,  f, s_t, p0ic, p1ic)
pic2 = copy(pic)
fdtd!(pic2, xr, xs, Nt,  f, s_t, p0ic, p1ic)
@time ppic = fdtd(xr, xs, Nt,  f, s_t, p0ic, p1ic)
println("                       in-place call")
@time        fdtd!(pic2, xr, xs, Nt,  f, s_t, p0ic, p1ic)

@test norm(pic-pic2) < 1e-8
@test norm(pic-ppic) < 1e-8

println("---------- full no IC")
P    = fdtd(    xs, Nt,geo, s_t)
P    = fdtd(    xs, Nt,  f, s_t)
P2 = copy(P)
fdtd!(P2,    xs, Nt,  f, s_t)
@time P    = fdtd(    xs, Nt,  f, s_t)
println("                       in-place call")
@time        fdtd!(P2, xs, Nt,  f, s_t)

@test vecnorm(P-P2) < 1e-8

println("---------- full    IC")
PIC  = fdtd(    xs, Nt,geo, s_t, p0ic, p1ic)
PIC  = fdtd(    xs, Nt,  f, s_t, p0ic, p1ic)
PIC2 = copy(PIC)
fdtd!(PIC2,    xs, Nt,  f, s_t, p0ic, p1ic)
@time PIC  = fdtd(    xs, Nt,  f, s_t, p0ic, p1ic)
println("                       in-place call")
@time      fdtd!(PIC2,xs, Nt,  f, s_t, p0ic, p1ic)

@test vecnorm(PIC-PIC2) < 1e-8

println("---------- source everywhere, no IC")
ps    = fdtd(    xr, Nt,geo, S)
ps    = fdtd(    xr, Nt,  f, S)
ps2 = copy(ps)
fdtd!(ps2, xr, Nt,  f, S)
@time ps    = fdtd(      xr, Nt,  f, S)
println("                       in-place call")
@time        fdtd!(ps2,  xr, Nt,  f, S)

println("---------- source everywhere, IC")
psic    = fdtd(    xr, Nt,geo, S, p0ic, p1ic)
psic    = fdtd(    xr, Nt,  f, S, p0ic, p1ic)
psic2 = copy(ps)
fdtd!(psic2, xr, Nt,  f, S, p0ic, p1ic)
@time psic    = fdtd(      xr, Nt,  f, S, p0ic, p1ic)
println("                       in-place call")
@time        fdtd!(psic2,  xr, Nt,  f, S, p0ic, p1ic)

@test vecnorm(psic-psic2) < 1e-8

println("---------- full no IC, source everywhere")
Ps    = fdtd(        Nt,geo, S)
Ps    = fdtd(        Nt,  f, S)
Ps2 = copy(Ps)
fdtd!(Ps2,    Nt,  f, S)
@time Ps    = fdtd(        Nt,  f, S)
println("                       in-place call")
@time        fdtd!(Ps2,    Nt,  f, S)

@test vecnorm(Ps-Ps2) < 1e-8

println("---------- full IC, source everywhere")
Psic    = fdtd(        Nt,geo, S, p0ic, p1ic)
Psic    = fdtd(        Nt,  f, S, p0ic, p1ic)
Psic2 = copy(Ps)
fdtd!(Psic2,    Nt,  f, S)
@time Psic    = fdtd(        Nt,  f, S, p0ic, p1ic)
println("                       in-place call")
@time        fdtd!(Psic2,    Nt,  f, S, p0ic, p1ic)

@test vecnorm(Psic-Psic2) < 1e-8

println("\ntest equivalence with matrix formulation. \n")

B = get_B(geo,Nt)
E = get_Exp(geo,Nt,xs)
F = get_Sel(geo,Nt,xr)
s = E*s_t[:]
s =[p0ic[:];p1ic[:];s[2*geo.Nx*geo.Ny*geo.Nz+1:end]]
S =[p0ic[:];p1ic[:];S[:]]
pm = F*(B\s)
println("sparse matrix      call")
@time pm = F*(B\s)
println("sparse matrix full call")
Pm = B\s
@time Pm = B\s

Pms = B\S

@test norm(PIC[:]-Pm[2*geo.Nx*geo.Ny*geo.Nz+1:end]) <1e-8
@test norm(pic[:]-pm[:]) <1e-8

@test norm(Psic[:]-Pms[2*geo.Nx*geo.Ny*geo.Nz+1:end]) <1e-8

println("\ntest resonances on cuboid room are correct \n")
Nt = 3000
X = 0.4
geo = CuboidRoom(10, 11, 12, 1e8*ones(6), FDTDEnv(X,IISO()))

using DSP
s_t = zeros(Nt)   #source signal
s_t[3] = 1
f2 = geo.env.Fs/2*0.14#cut-off frequency of source
filt!(s_t,digitalfilter(Bandpass(10,f2;fs = geo.env.Fs),Butterworth(15)),s_t)

xr = (2, 2, 2)
xs = (geo.Nx-1, geo.Ny-1, geo.Nz-1)
p_t2 = fdtd(xr, xs, Nt, geo, s_t)

f = linspace(0,geo.env.Fs/2,div(Nt,2)+1)
fres  = zeros(5,5,5)
for l = 0:4,m=0:4,i =0:4

	fres[l+1,m+1,i+1] = 343/(2*π)*sqrt(  (l*π/(geo.Lx))^2+(m*π/(geo.Ly))^2+(i*π/(geo.Lz))^2 )
end
fres = sort(fres[:])[2:7]

P = 10*log10(abs(rfft(p_t2)))

##using PyPlot
##figure()
##plot(f,P)
##plot(fres[:],5.* ones(length(fres[:])), "r*" )
##xlim([0,400])

@test norm(f[P.>5.5][1:6]-fres)<0.3 #test on resonance frequencies

println("\ntest acoustic delay is correct \n")
X = 0.4
Nt = 500
geo = CuboidRoom(20, 21, 24, 4*ones(6), FDTDEnv(X,IWB()))
s_t = zeros(Nt)   #source signal
s_t[3] = 1
f2 = geo.env.Fs/2*0.1#cut-off frequency of source
filt!(s_t,digitalfilter(Bandpass(10,f2;fs = geo.env.Fs),Butterworth(2)),s_t)

xs = (div(geo.Nx,2), div(geo.Ny,2)  , div(geo.Nz,2))
xr = (div(geo.Nx,2), div(geo.Ny,2)+4, div(geo.Nz,2))
d  = norm(X*[ xs[1]-xr[1],xs[2]-xr[2],xs[3]-xr[3]])

p_t2 = fdtd(xr, xs, Nt, geo, s_t)
t = 1/geo.env.Fs:1/geo.env.Fs:Nt/geo.env.Fs

##using PyPlot
##plot(t,p_t2)
##plot(d/343.+2*X/343, 0., "r*")#+2 due to initial conditions

#test delay
@test norm((indmax(xcorr(p_t2[:],s_t)[end-length(p_t2):end])-1)/geo.env.Fs-d/343.)<1e-3

max_stars_repo_name: nantonel/JuFdtd

max_issues_repo_name: nantonel/JuFdtd

max_forks_repo_name: nantonel/JuFdtd


text:
using Makie, FileIO
logo = load(joinpath(@__DIR__, "..", "docs", "src", "assets", "logo.png"))
# For rectangles + rotation

scatter(
    Point3f0[(1,0,0), (0,1,0), (0,0,1)],
    marker = [logo, logo, logo], markersize = reverse(size(logo) ./ 1500),
    # rotation around axis... can you Vec4f0(...) for a quaternion
    rotations = [(Vec3f0(0, 1, 0), 0.5pi), (Vec3f0(1, 0, 0), -0.5pi), (Vec3f0(0, 0, 1), -1.2pi)]
)

x = [heatmap(rand(100, 100)) for i = 1:3, j = 1:3]
grid(x)

# for more control:
scene = Scene(resolution = (500, 500))
mesh = GLNormalUVMesh(GeometryTypes.SimpleRectangle(0, 0, 1, 1))
# note, you can change the mesh.vertices to arbitrary values
Makie.mesh(mesh, color = logo, shading = false)

s2 = scatter(LinRange(0, 10, 10), rand(10))
s2 = scatter(LinRange(0, 1, 10), rand(10), rand(10))

scene = Scene()
x = range(0, stop = 6, length = 100)
s = surface!(scene, x, x, (x, y)-> sin(x) + cos(y), show_legend = true)
scene

scene = Scene()
r = range(-10, stop = 10, length = 512)
z = ((x, y)-> sin(x) + cos(y)).(r, r')
contour!(scene, r, r, z, levels = 5, color = :RdYlBu)
Makie.save(joinpath(homedir(), "Desktop", "contour.png"), scene)

scene = Scene()
x = range(0, stop = 6, length = 100)
s = heatmap!(scene, x, x, (x, y)-> sin(x) + cos(y), show_legend = true)
scene

using Makie
p = series(rand(10, 5));

# Makie.campixel!(scene)
# image!(scene, GLVisualize.loadasset("doge.png"), raw = true)

# Update limits when zooming with cam
# Update limits when subplot changes /new
# 3D  axis + limits
# remaining primitives ()

using GLVisualize, Makie
scene = Scene()
m = mesh!(scene, loadasset("cat.obj"))
scene

r = range(-10, stop = 10, length = 512)
z = ((x, y)-> sin(x) + cos(y)).(r, r')
scene = Scene()
# c = heatmap!(scene, r, r, z, levels = 5)
# s = scatter!(scene, rand(-10:10, 10), rand(-10:10, 10))
s2 = contour!(scene, r, r, z, levels = 5)
scene

scatter!(cscene, [1, 0, 1, 0], [1, 0, 0, 1],
    markersize = 0.02, color = :black
);

# a = text!(scene, "Hellooo", color = :black, textsize = 0.1, position = (0.5, 0.5))
# a = text!(scene, "Hellooo", color = :black, textsize = 0.1, position = ps)
tb = Makie.TextBuffer(scene)
start!(tb)
append!(tb, ["Hello there"], Point2f0[(0.5, 0.5)], textsize = 0.1, color = :black, rotation = 0.0, font = "default")
append!(tb, ["General Kenobi!"], Point2f0[(0.0, 0.0)])
finish!(tb)

scene = Scene()
lsb = Makie.LinesegmentBuffer(scene, Point2)
start!(lsb)
append!(lsb, Point2f0[(1, 1), (0, 0)], color = RGBAf0(0,0,0,1), linewidth = 1f0)
finish!(lsb)
scene

s = scatter!(scene, [0, 0, 1, 1], [0, 1.5, 0, 1.5])
s = lines!(scene, FRect(0, 0, 1, 1.5), color = :black, show_axis = true)
xy = range(0, stop = 2pi, length = 100)
f(x, y) = sin(x) + cos(y)
z = f.(xy, xy')
s1 = heatmap!(scene, 0.1 .. 0.9, 0.1 .. 0.44, z)
# TODO linewidth varies and comes out quite different compared to `lines`,
# because it is actually in defined in gradients of z (at least for opengl)
# Can user deal with this or does this need to be changed?
s2 = contour!(scene, 0.1 .. 0.9, 0.46 .. 0.9, z, linewidth = 0.1, fillrange = true)
s2 = contour!(scene, 0.1 .. 0.9, 0.91 .. 1.4, z, linewidth = 2)

# screen = Screen(scene)
scene = Scene()
a = text!(scene, "Hellooo", color = :black, textsize = 0.1, position = (0.5, 0.5))

b = scatter!(scene, rand(10), rand(10))
b = linesegments!(scene, rand(10), rand(10))
c = plot!(scene, rand(10), rand(10), color = :white)
d = meshscatter!(scene, rand(10), rand(10), rand(10), show_axis = true);
cam = cam2d!(scene)
cam.area[] = FRect(0, 0, normalize(widths(scene.px_area[])) * 3)
update_cam!(scene, cam)
scene

scene = Scene()
scene.px_area[] = IRect(0, 0, 1920, 1080)
cam = cam2d!(scene)
cam.area[] = FRect(0, 0, normalize(widths(scene.px_area[])) * 3)
update_cam!(scene, cam)
scatter!(scene, FRect(0, 0, 1, 1), scale_plot = false, linewidth = 5)
h = heatmap!(scene, range(0, stop = 1, length = 50), range(0, stop = 1, length = 50), rand(50, 50))
scene
# cam.rotationspeed[] = 0.1
# cam.pan_button[] = Mouse.right
# scene.events.window_dpi[]

# screenw = widths(scene.px_area[])
# camw = widths(scene.area[])
#
# screen_r = screenw ./ screenw[1]
# camw_r = camw ./ camw[1]
# r = (screen_r ./ camw_r)
# r = r ./ maximum(r)
#
# update_cam!(scene, FRect(minimum(scene.area[]), r .* camw))
using Makie

scene = Scene()
m = GLVisualize.loadasset("cat.obj")
mesh!(scene, m)
Makie.cam3d!(scene)
scene

using Makie, GeometryTypes
scene = Scene()
scene.px_area[] = IRect(0, 0, 1920, 1080)
cam = cam2d!(scene)
points = decompose(Point2f0, Circle(Point2f0(0), 500f0))
poly!(
    scene, points,
    color = :gray, linewidth = 10, linecolor = :black
)
scene

pol[:positions] = Circle(Point2f0(250), 500f0)
pol[:linewidth] = 2
# Optimized forms
y = poly([Circle(Point2f0(600+i, i), 50f0) for i = 1:150:800])
x = poly([Rectangle{Float32}(600+i, i, 100, 100) for i = 1:150:800], strokewidth = 10, strokecolor = :black)
x = linesegment([Point2f0(600+i, i) => Point2f0(i + 700, i + 100)

struct Test <: AbstractArray{Float32, 1} end
x = Test()
@which filter(x-> x > 0.5, x)
As[map(f, As)::AbstractArray{Bool}]

using Makie
scene = Scene()
t = text!(
    scene,
    ". This is an annotation!",
    position = (300, 200),
    align = (:center,  :center),
    textsize = 60,
    font = "Comic Sans MS",
    scale_plot = false
)
dlims = Makie.data_limits(t)[]
rect =  FRect(dlims[1], dlims[2] .- dlims[1])
lines!(scene, rect, scale_plot = false)

using Reactive

keys = (:position, :textsize, :font, :align, :rotation, :model)
x = t
args = (value(x.args[1]), value.(getindex.(x.attributes, keys))...)

positions, scale = Makie.layout_text(args...)
scatter!(scene, positions, markersize = 10, color = :red, scale_plot = false)
scatter!(scene, positions .+ scale, markersize = 10, color = :green, scale_plot = false)
rect = Makie.FRect2D(union(AABB(positions .+ scale), AABB(positions)))

ex = extrema(vcat(positions, positions .+ scale))
(last.(ex), first.(ex))

using Transpiler, GeometryTypes

import Transpiler.gli: texelFetch

function image_preprocess(img::AbstractMatrix{NTuple{4, T}}, uv, mini, maxi, kernel, brightness, contrast) where T
    @inbounds begin
        z = 1f0
        val = (z, z, z, 1f0)
        # filter
        for i in 1:3
            for j in 1:3
                tmp = img[Cint(i), Cint(j)] .* kernel[i, j]
                val = val .+ tmp
            end
        end
        # brightness
        val = val .+ brightness
        # contrast
        val = val .- 0.5f0
        val = (contrast .* val)
        val = val .+ 0.5f0
        val = val .- mini
        val = val ./ (maxi - mini)
        # min max
        return clamp.(val, 0f0, 1f0)
    end
end
decl = Transpiler.GLMethod((image_preprocess, (
    Matrix{NTuple{4, Float32}}, Tuple{Int, Int}, Float32, Float32, Mat3f0, Float32, Float32
)))

# image_preprocess(rand(50, 50), ())

source = Sugar.getsource!(decl);
println(source)

deps = reverse(collect(Sugar.dependencies!(decl, true)));
types = filter(Sugar.istype, deps)
Transpiler.is_fixedsize_array(Mat4f0)
funcs = filter(Sugar.isfunction, deps)
io = STDOUT
println(io, "// dependant type declarations")
for typ in types
    if !Sugar.isintrinsic(typ)
        println(io, Sugar.getsource!(typ))
    end
end
# defer printing for the function and types
io2 = GLIO(IOBuffer(), m)
for func in funcs
    if !Sugar.isintrinsic(func)
        println("//", func.signature)
        println(io, Sugar.getsource!(func))
    end
end


text:
# Resolve o Exemplo 6
# https://leonardosecchin.github.io/juliaopt_ex6/
# via WORHP

using JuMP, AmplNLWriter

P = Model(() -> AmplNLWriter.Optimizer("worhp_ampl"))

@variable(P, x[1:2])

@NLobjective(P, Min, (x[1] - 2.0)^2 + (x[2] - 1.0)^2)

@constraint(P, x[1] + x[2] - 2 <= 0)

@NLconstraint(P, x[1]^2 - x[2] <= 0)

optimize!(P)

max_stars_repo_name: leonardosecchin/temp

max_issues_repo_name: leonardosecchin/temp

max_forks_repo_name: leonardosecchin/temp


text:
using DefferedAcceptance
using Base.Test

const _deferred_acceptance = DefferedAcceptance.ikegami
const test_matrix = false

function mat2vecs{T<:Integer}(prefs::Matrix{T})
    return [prefs[1:findfirst(prefs[:, j], 0)-1, j] for j in 1:size(prefs, 2)]
end

function sort_matches!{T<:Integer}(matches::Vector{T}, indptr::Vector{T})
    for i in 1:length(indptr)-1
        sort!(view(matches, indptr[i]:indptr[i+1]-1))
    end
end

@testset "Testing deferred acceptance" begin
    matchings_one_to_one = Array{Dict, 1}()
    matchings_many_to_one = Array{Dict, 1}()

    # Test case: one-to-one
    m, n = 4, 3
    # Males' preference orders over females [1, 2, 3] and unmatched 0
    m_prefs = [1 3 2 3;
               2 1 3 1;
               3 2 1 2;
               0 0 0 0]
    # Females' preference orders over males [1, 2, 3, 4] and unmatched 0
    f_prefs = [3 1 3;
               1 2 0;
               2 3 2;
               4 4 1;
               0 0 4]
    # Unique stable matching
    m_matches_m_opt = m_matches_f_opt = [1, 2, 3, 0]
    f_matches_m_opt = f_matches_f_opt = [1, 2, 3]
    d = Dict(
        "m_prefs" => m_prefs,
        "f_prefs" => f_prefs,
        "m_matches_m_opt" => m_matches_m_opt,
        "f_matches_m_opt" => f_matches_m_opt,
        "m_matches_f_opt" => m_matches_f_opt,
        "f_matches_f_opt" => f_matches_f_opt,
    )
    push!(matchings_one_to_one, d)

    # Test case: one-to-one; from Roth and Sotomayor Example 2.9
    m, n = 5, 4
    m_prefs = [
        1, 2, 3, 4, 0,
        4, 2, 3, 1, 0,
        4, 3, 1, 2, 0,
        1, 4, 3, 2, 0,
        1, 2, 4, 0, 3,
    ]
    m_prefs = reshape(m_prefs, n+1, m)
    f_prefs = [
        2, 3, 1, 4, 5, 0,
        3, 1, 2, 4, 5, 0,
        5, 4, 1, 2, 3, 0,
        1, 4, 5, 2, 3, 0,
    ]
    f_prefs = reshape(f_prefs, m+1, n)
    m_matches_m_opt, f_matches_m_opt = [1, 2, 3, 4, 0], [1, 2, 3, 4]
    m_matches_f_opt, f_matches_f_opt = [4, 1, 2, 3, 0], [2, 3, 4, 1]
    d = Dict(
        "m_prefs" => m_prefs,
        "f_prefs" => f_prefs,
        "m_matches_m_opt" => m_matches_m_opt,
        "f_matches_m_opt" => f_matches_m_opt,
        "m_matches_f_opt" => m_matches_f_opt,
        "f_matches_f_opt" => f_matches_f_opt,
    )
    push!(matchings_one_to_one, d)

    # Test case: one-to-one; from Roth and Sotomayor Example 2.17
    m = n = 4
    m_prefs = [
        1, 2, 3, 4, 0,
        2, 1, 4, 3, 0,
        3, 4, 1, 2, 0,
        4, 3, 2, 1, 0,
    ]
    m_prefs = reshape(m_prefs, n+1, m)
    f_prefs = [
        4, 3, 2, 1, 0,
        3, 4, 1, 2, 0,
        2, 1, 4, 3, 0,
        1, 2, 3, 4, 0,
    ]
    f_prefs = reshape(f_prefs, m+1, n)
    m_matches_m_opt = f_matches_m_opt = [1, 2, 3, 4]
    m_matches_f_opt = f_matches_f_opt = [4, 3, 2, 1]
    d = Dict(
        "m_prefs" => m_prefs,
        "f_prefs" => f_prefs,
        "m_matches_m_opt" => m_matches_m_opt,
        "f_matches_m_opt" => f_matches_m_opt,
        "m_matches_f_opt" => m_matches_f_opt,
        "f_matches_f_opt" => f_matches_f_opt,
    )
    push!(matchings_one_to_one, d)

    # Test case: many-to-one
    # From Gusfield and Irving (1989, Section 1.6.5)

    m = 11  # Number of students
    n = 5   # Number of colleges

    # Students' preference orders over colleges 1, ..., 5 and unmatched 0
    s_prefs = [3, 1, 5, 4, 0, 2,
               1, 3, 4, 2, 5, 0,
               4, 5, 3, 1, 2, 0,
               3, 4, 1, 5, 0, 2,
               1, 4, 2, 0, 3, 5,
               4, 3, 2, 1, 5, 0,
               2, 5, 1, 3, 0, 4,
               1, 3, 2, 5, 4, 0,
               4, 1, 5, 0, 2, 3,
               3, 1, 5, 2, 4, 0,
               5, 4, 1, 3, 2, 0]
    s_prefs = reshape(s_prefs, n+1, m)

    # Colleges' preference orders over students 1, ..., 11 and unmatched 0
    c_prefs = [3, 7, 9, 11, 5, 4, 10, 8, 6, 1,
               2, 0, 5, 7, 10, 6, 8, 2, 3, 11,
               0, 1, 4, 9, 11, 6, 8, 3, 2, 4,
               7, 1, 10, 0, 5, 9, 10, 1, 2, 11,
               4, 9, 5, 3, 6, 8, 0, 7, 2, 4,
               10, 7, 6, 1, 8, 3, 11, 9, 0, 5]
    c_prefs = reshape(c_prefs, m+1, n)

    # Capacities for colleges
    caps = [4, 1, 3, 2, 1]
    indptr = [1, 5, 6, 9, 11, 12]

    # Optimal stable matchings

    # Student optimal
    s_matches_s_opt = [3, 1, 4, 3, 1, 3, 2, 1, 4, 1, 5]
    c_matches_s_opt = [2, 5, 8, 10, 7, 1, 4, 6, 3, 9, 11]

    # Matching the book claims to be hospital-(or college-)optimal
    # s_matches_c_opt = [4, 4, 3, 1, 1, 3, 2, 3, 1, 5, 1]
    # c_matches_c_opt = [4, 5, 9, 11, 7, 3, 6, 8, 1, 2, 10]

    # Matching actually hospital-(or college-)optimal
    s_matches_c_opt = [4, 5, 3, 1, 1, 3, 2, 3, 1, 4, 1]
    c_matches_c_opt = [4, 5, 9, 11, 7, 3, 6, 8, 1, 10, 2]

    d = Dict(
        "s_prefs" => s_prefs,
        "c_prefs" => c_prefs,
        "caps" => caps,
        "s_matches_s_opt" => s_matches_s_opt,
        "c_matches_s_opt" => c_matches_s_opt,
        "s_matches_c_opt" => s_matches_c_opt,
        "c_matches_c_opt" => c_matches_c_opt,
        "indptr" => indptr,
    )
    push!(matchings_many_to_one, d)

    # Test case: many-to-one; from Roth and Sotomayor Page 16
    m, n = 7, 5
    s_prefs = [
        5, 1, 0, 2, 3, 4,
        2, 5, 1, 0, 3, 4,
        3, 1, 0, 2, 4, 5,
        4, 1, 0, 2, 3, 5,
        1, 2, 0, 3, 4, 5,
        1, 3, 0, 2, 4, 5,
        1, 3, 4, 0, 2, 6,
    ]
    s_prefs = reshape(s_prefs, n+1, m)
    c_prefs = [
        1, 2, 3, 4, 5, 6, 7, 0,
        5, 2, 0, 1, 3, 4, 6, 7,
        6, 7, 3, 0, 1, 2, 4, 5,
        7, 4, 0, 1, 2, 3, 5, 6,
        2, 1, 0, 3, 4, 5, 6, 7,
    ]
    c_prefs = reshape(c_prefs, m+1, n)
    caps = [3, 1, 1, 1, 1]
    indptr = [1, 4, 5, 6, 7,  8]
    s_matches_s_opt = [5, 2, 3, 4, 1, 1, 1]
    c_matches_s_opt = [5, 6, 7, 2, 3, 4, 1]
    s_matches_c_opt = [1, 5, 1, 1, 2, 3, 4]
    c_matches_c_opt = [1, 3, 4, 5, 6, 7, 2]
    d = Dict(
        "s_prefs" => s_prefs,
        "c_prefs" => c_prefs,
        "caps" => caps,
        "s_matches_s_opt" => s_matches_s_opt,
        "c_matches_s_opt" => c_matches_s_opt,
        "s_matches_c_opt" => s_matches_c_opt,
        "c_matches_c_opt" => c_matches_c_opt,
        "indptr" => indptr,
    )
    push!(matchings_many_to_one, d)

    @testset "one-to-one: Vector of Vectors" begin
        for d in matchings_one_to_one
            # Convert Matrix to Vector{Vector}
            m_prefs, f_prefs = mat2vecs.([d["m_prefs"], d["f_prefs"]])

            # Male proposal
            m_matches, f_matches = _deferred_acceptance(m_prefs, f_prefs)
            @test m_matches == d["m_matches_m_opt"]
            @test f_matches == d["f_matches_m_opt"]

            # Female proposal
            f_matches, m_matches = _deferred_acceptance(f_prefs, m_prefs)
            @test m_matches == d["m_matches_f_opt"]
            @test f_matches == d["f_matches_f_opt"]
        end
    end

    @testset "many-to-one: Vector of Vectors" begin
        for d in matchings_many_to_one
            s_prefs, c_prefs = mat2vecs.(d["s_prefs"], d["c_prefs"])

            # Default (Student proposal)
            s_matches, c_matches, indptr =
                deferred_acceptance(s_prefs, c_prefs, d["caps"])
            sort_matches!(c_matches, indptr)
            @test s_matches == d["s_matches_s_opt"]
            @test c_matches == d["c_matches_s_opt"]
            @test indptr == d["indptr"]
        end
    end

    if test_matrix
        @testset "one-to-one: Matrix" begin
            for d in matchings_one_to_one
                m_prefs, f_prefs = d["m_prefs"], d["f_prefs"]

                # Male proposal
                m_matches, f_matches = _deferred_acceptance(m_prefs, f_prefs)
                @test m_matches == d["m_matches_m_opt"]
                @test f_matches == d["f_matches_m_opt"]

                # Female proposal
                f_matches, m_matches = _deferred_acceptance(f_prefs, m_prefs)
                @test m_matches == d["m_matches_f_opt"]
                @test f_matches == d["f_matches_f_opt"]
            end
        end

        @testset "many-to-one: Matrix" begin
            for d in matchings_many_to_one
                s_prefs, c_prefs = d["s_prefs"], d["c_prefs"]

                # Default (Student proposal)
                s_matches, c_matches, indptr =
                    deferred_acceptance(s_prefs, c_prefs, d["caps"])
                sort_matches!(c_matches, indptr)
                @test s_matches == d["s_matches_s_opt"]
                @test c_matches == d["c_matches_s_opt"]
                @test indptr == d["indptr"]
            end
        end
    end
end


text:
#!/usr/bin/env julia
# File              : Utils.jl
# Author            : Denis N Dumoulin <denis.dumoulin@uclouvain.be>
# Date              : 14.03.2019
# Last Modified Date: 08.12.2020

#Constant values for Gauss quadrature in triangles
const wgauss = [4.40470137440027E-04 2.73540341982354E-03 6.06288503309780E-03 7.32940618020361E-03 4.35981888065413E-03 2.54894897110758E-03 1.58294584351310E-02 3.50852038616632E-02 4.24144130415724E-02 2.52297599892837E-02 4.35986399676418E-03 2.70755855459964E-02 6.00116828031788E-02 7.25479695591883E-02 4.31543838150997E-02 2.80482254237839E-03 1.74184820315194E-02 3.86071952834637E-02 4.66721394461879E-02 2.77624230060559E-02 3.69367248237134E-04 2.29384093975323E-03 5.08418385425357E-03 6.14625683304730E-03 3.65603514489964E-03]

const xgauss = [5.65222820508010E-03 5.65222820508010E-03 5.65222820508010E-03 5.65222820508010E-03 5.65222820508010E-03 7.34303717426523E-02 7.34303717426523E-02 7.34303717426523E-02 7.34303717426523E-02 7.34303717426523E-02 2.84957404462558E-01 2.84957404462558E-01 2.84957404462558E-01 2.84957404462558E-01 2.84957404462558E-01 6.19482264084778E-01 6.19482264084778E-01 6.19482264084778E-01 6.19482264084778E-01 6.19482264084778E-01 9.15758083004698E-01 9.15758083004698E-01 9.15758083004698E-01 9.15758083004698E-01 9.15758083004698E-01]

const ygauss = [5.62028052139780E-03 7.30153265243790E-02 2.83346760183808E-01 6.15980808959171E-01 9.10582009338909E-01 5.23718298680677E-03 6.80383522483882E-02 2.64032876322051E-01 5.73993451145053E-01 8.48513626543325E-01 4.04158392633041E-03 5.25058436021453E-02 2.03756682104520E-01 4.42956206000591E-01 6.54806036556072E-01 2.15077307947324E-03 2.79415588029271E-02 1.08431346378371E-01 2.35723988569175E-01 3.48462192391012E-01 4.76154539290863E-04 6.18591528127868E-03 2.40053580139315E-02 5.21863734710916E-02 7.71452164162586E-02]

"""
    rotate(X, Y, θ)

Rotate frame by angle `θ` in the clockwise direction.
"""
function rotate(X, Y, θ)
	x = X*cos(θ) - Y*sin(θ)
	y = Y*cos(θ) + X*sin(θ)
	return [x, y]
end

"""
    norm2(x)

Squared norm of `x`.
"""
norm2(x...) = sum(x.^2)

"""
    simps(f, x)

Integrate function `f` using Simpson's rule.
"""
function simps(f, x)
	N = length(x) - 1
	N < 2 && return trapz(f, x)
    h = diff(x)

    result = 0.
    for i in 2:2:N
        hph = h[i] + h[i - 1]
        result += f[i]*(h[i]^3 + h[i-1]^3 + 3h[i]*h[i - 1]*hph)/(6h[i]*h[i-1])
        result += f[i-1]*(2h[i-1]^3 - h[i]^3 + 3h[i]*h[i-1]^2)/(6h[i-1]*hph)
        result += f[i+1]*(2h[i]^3 - h[i-1]^3 + 3h[i-1]*h[i]^2)/(6h[i]*hph)
	end

    if (N + 1) % 2 == 0
        result += f[N]*(2h[N-1]^2 + 3h[N-2]*h[N-1])/(6(h[N-2] + h[N-1]))
        result += f[N-1]*(h[N-1]^2 + 3h[N-1]*h[N-2])/(6h[N-2])
        result -= f[N-2]*h[N-1]^3/(6*h[N-2]*(h[N-2] + h[N-1]))
	end
	return result
end

"""
    trapz(f, x)

Integrate function `f` using trapeze rule.
"""
function trapz(y, x)
	n = length(y)
	length(y) == length(x) || error("Dimension mismatch between x and f(x).")

	r = zero(zero(eltype(x))*zero(eltype(y)))
	for i in 2:n
		r += (x[i] - x[i-1])*(y[i] + y[i-1])
	end
	.5r
end

"""
    backwarddifference(x, old_x, dt)

Compute the backward difference of nth order according to the length of `old_x`.
"""
function backwarddifference(x, old_x, dt)
	n = length(old_x)
	0 < n < 5 || error("5+th order Backward Difference undefined.")
	A = [-1 1 0 0 0;
		 -1.5 2 -.5 0 0;
		 -11/6 3 -1.5 1/3 0;
		 -25/12 4 -3 4/3 -.25]
	-sum(A[n, 1:n+1].*[x, old_x...])/dt
end

"""
    customfourthorder(x, old_x, dt)

Compute the 4th order time derivative at t using
	t-1/2
	t-3/2
	t-5/2
	t-7/2.
"""
function customfourthorder(x, old_x, dt)
	A = [.5 1.5 2.5 3.5;
		 1 9 25 49;
		 1 27 125 343;
		 1 81 625 2401]
	a = A\[1 0 0 0]'
	return sum([sum(a), -a...].*[x, old_x...])/dt
end

"""
    centereddifference(neighbours, dx)

Centered differences. Order depends on length of neighbours.
"""
function centereddifference(neighbours, dx)
	n = length(neighbours)
	if n == 2
		return [-.5 .5]*neighbours /dx
	elseif n == 4
		return [1/12 -2/3 2/3 -1/12]*neighbours /dx
	else
		error("6+th order Backward Difference undefined.")
	end
end

"""
    checkangle(θ)

Assert the continuity of panel angles by subtracting 2π if necessary.
"""
function checkangle(θ)
	θ > .5π && return θ-2π
	return θ
end

"""
    gausslegendretriangle(f, fx, fy)

Computes the integral of function `f(x,y)`
on the standard triangle abc:
a (0; 0)
b (1; 0)
c (0; 1)
using 25 points per dimension.
fx, fy are the mapping from the standard triangle to the real one.
Coefficients are extracted from Rathod et al (2004)
"Gauss Legendre quadrature over a triangle"
"""
function gausslegendretriangle(f, fx, fy)
	x = fx.(xgauss, ygauss)
	y = fy.(xgauss, ygauss)
	return (f.(x, y)*wgauss')[1]
end

"""
    get∫fdV(XY, f)

Computes the integral of `f` on the triangle created by the
set of points `XY`, `Y` using a Gauss-Legendre quadrature.
`f` can return a vector.
"""
function get∫fdV(XY1, XY2, XY3, f)
	J =	(XY3[1] - XY2[1])*(XY1[2] - XY2[2]) - (XY1[1] - XY2[1])*(XY3[2] - XY2[2])
	x(u, v)	= XY2[1] + u*(XY3[1] - XY2[1]) + v*(XY1[1] - XY2[1])
	y(u, v)	= XY2[2] + u*(XY3[2] - XY2[2]) + v*(XY1[2] - XY2[2])
	Is = gausslegendretriangle(f, x, y)
	return abs(J).*Is
end

"""
    cubicspline(z, X, U)

Computes the cubic spline at point `z ∈ [0, 1]`, build from
interpolation of points `(X, U)` as well as its derivative.
"""
function cubicspline(z, X, U)
	@assert length(X) == length(U) "X and U must have the same length!"
	@assert 0 ≤ z ≤ 1 "z must be on the interval [0, 1]"

	n = length(U)-1
	dv = 4ones(n+1)
	dv[1] = 2
	dv[end] = 2
	ev = ones(n)
	A = SymTridiagonal(dv, ev)
	x = map((a,b)->3(b-a), U, circshift(U, -2))
	x[1] = 3(U[2] - U[1])
	x[end] = 3(U[end] - U[end-1])
	b = (A\x)
	a = U[1:end-1]
	c = map((w, x, y, z)->3(x-w) - 2y - z, U, circshift(U, -1), b, circshift(b, -1))
	d = map((w, x, y, z)->2(w-x) + y + z, U, circshift(U, -1), b, circshift(b, -1))
	b = b[1:end-1]
	c = c[1:end-1]
	d = d[1:end-1]

	x = (z*ones(n+1) - X) ./ (circshift(X, -1) - X)
	i = findfirst(0 .≤ x[1:end] .≤ 1)
	return a[i] + b[i]*x[i] + c[i]*x[i]^2 + d[i]*x[i]^3

end

"""
    getbodymass(XYm, XYt, XYb)

Integrate the constant and uniform unit density on the body.
Integral computed with 25-point Gauss quadrature.
"""
function getbodymass(XYm, XYt, XYb)
	I12 = sum(map((XY1, XY2, XY3)->get∫fdV(XY1, XY2, XY3, (x, y)->1), XYm[1:end-1], XYm[2:end], XYt[2:end]))
	I22 = sum(map((XY1, XY2, XY3)->get∫fdV(XY1, XY2, XY3, (x, y)->1), XYm[1:end-1], XYm[2:end], XYb[2:end]))
	I32 = sum(map((XY1, XY2, XY3)->get∫fdV(XY1, XY2, XY3, (x, y)->1), XYm[1:end-1], XYt[1:end-1], XYt[2:end]))
	I42 = sum(map((XY1, XY2, XY3)->get∫fdV(XY1, XY2, XY3, (x, y)->1), XYm[1:end-1], XYb[1:end-1], XYb[2:end]))
	return (I12 + I22 + I32 + I42)
end

"""
    getinertia(XYm, XYt, XYb)

Return the moment of inertia of the body through 25-point Gauss quatradure on triangles.
"""
function getinertia(XYm, XYt, XYb)
	I1 = sum(map((XY1, XY2, XY3)->get∫fdV(XY1, XY2, XY3, (x, y)->norm2(x, y)), XYm[1:end-1], XYm[2:end], XYt[2:end]))
	I2 = sum(map((XY1, XY2, XY3)->get∫fdV(XY1, XY2, XY3, (x, y)->norm2(x, y)), XYm[1:end-1], XYm[2:end], XYb[2:end]))
	I3 = sum(map((XY1, XY2, XY3)->get∫fdV(XY1, XY2, XY3, (x, y)->norm2(x, y)), XYm[1:end-1], XYt[1:end-1], XYt[2:end]))
	I4 = sum(map((XY1, XY2, XY3)->get∫fdV(XY1, XY2, XY3, (x, y)->norm2(x, y)), XYm[1:end-1], XYb[1:end-1], XYb[2:end]))
	return I1 + I2 + I3 + I4
end

"""
    getcenterofmass(xm, ym, xb, yb, xt, yt)

Compute the center of mass of the body.
Integrals computed with 25-point Gauss quadrature.
"""
function getcenterofmass(xm, ym, xb, yb, xt, yt)
	XYm = map((x,y)->[x, y], xm, ym)
	XYt = map((x,y)->[x, y], xt, yt)
	XYb = map((x,y)->[x, y], xb, yb)
	I1 = sum(map((XY1, XY2, XY3)->get∫fdV(XY1, XY2, XY3, (x, y)->[x, y]), XYm[1:end-1], XYm[2:end], XYt[2:end]))
	I2 = sum(map((XY1, XY2, XY3)->get∫fdV(XY1, XY2, XY3, (x, y)->[x, y]), XYm[1:end-1], XYm[2:end], XYb[2:end]))
	I3 = sum(map((XY1, XY2, XY3)->get∫fdV(XY1, XY2, XY3, (x, y)->[x, y]), XYm[1:end-1], XYt[1:end-1], XYt[2:end]))
	I4 = sum(map((XY1, XY2, XY3)->get∫fdV(XY1, XY2, XY3, (x, y)->[x, y]), XYm[1:end-1], XYb[1:end-1], XYb[2:end]))
	mass = getbodymass(XYm, XYt, XYb)
	return  [((I1 .+ I2 .+ I3 .+ I4) ./ mass)..., mass]
end

"""
    getcenterofmassvelocity(xm, ym, xb, yb, xt, yt, ub, vb, ut, vt)

Compute the velocity of the center of mass of the body.
Integrals computed with bi-linear assumption on quadrilaterals.
"""
function getcenterofmassvelocity(xm, ym, xb, yb, xt, yt, ub, vb, ut, vt)
	XYm = map((x,y)->[x, y], xm, ym)
	XYt = map((x,y)->[x, y], xt, yt)
	XYb = map((x,y)->[x, y], xb, yb)
	diags1 = (XYb[2:end] - XYt[1:end-1])
	diags2 = (XYb[1:end-1] - XYt[2:end])
	Areas = .5abs.(map((xy1, xy2)->xy1[1]*xy2[2] - xy1[2]*xy2[1], diags1, diags2))
	I1 = sum(.25Areas.*(ut[1:end-1] .+ ut[2:end] .+ ub[1:end-1] .+ ub[2:end]))
	I2 = sum(.25Areas.*(vt[1:end-1] .+ vt[2:end] .+ vb[1:end-1] .+ vb[2:end]))
	return  [I1; I2] ./ getbodymass(XYm, XYt, XYb)
end

"""
    correctvelocity(xm, ym, xb, yb, xt, yt, ub, vb, ut, vt)

Compute the angular velocity (ratio between angular impulse and
scalar inertia moment).
Integrals computed with bi-linear assumption on quadrilaterals,
and 25-point Gauss quadrature.
"""
function correctvelocity(xm, ym, xb, yb, xt, yt, ub, vb, ut, vt)
	XYm = map((x,y)->[x, y], xm, ym)
	XYt = map((x,y)->[x, y], xt, yt)
	XYb = map((x,y)->[x, y], xb, yb)
	diags1 = (XYb[2:end] - XYt[1:end-1])
	diags2 = (XYb[1:end-1] - XYt[2:end])
	Areas = .5abs.(map((xy1, xy2)->xy1[1]*xy2[2] - xy1[2]*xy2[1], diags1, diags2))
	height1 = map((x, y, u, v)->x*v - u*y, xt, yt, ut, vt)
	height3 = map((x, y, u, v)->x*v - u*y, xb, yb, ub, vb)
	I = sum(.25Areas.*(height1[1:end-1] .+ height1[2:end] .+ height3[1:end-1] .+ height3[2:end]))
	J = getinertia(XYm, XYt, XYb)
	return I/J, J
end


text:
@define HeaderCoefficientsARK begin
    @HeaderCoefficientsRK
    r::Int
end

@define CoefficientsARK begin
    @CoefficientsRK
    α::Matrix{T}
    β::Vector{T}
end

@define CoefficientsPRK begin
    a::Matrix{T}
    c::Vector{T}
    α::Matrix{T}
end

"Holds the coefficients of an additive Runge-Kutta method."
struct CoefficientsARK{T} <: AbstractCoefficients{T}
    @HeaderCoefficientsARK
    @CoefficientsARK

    function CoefficientsARK{T}(name,o,s,r,a,b,c,α,β) where {T}
        @assert T <: Real
        @assert isa(name, Symbol)
        @assert isa(o, Integer)
        @assert isa(s, Integer)
        @assert isa(r, Integer)
        @assert s > 0 "Number of stages s must be > 0"
        @assert r > 0 "Number of stages r must be > 0"
        @assert s==size(a,1)==size(a,2)==size(α,1)==length(b)==length(c)
        @assert r==size(α,2)==length(β)
        new(name,o,s,r,a,b,c,α,β)
    end
end

function CoefficientsARK(name::Symbol, order::Int, a::Matrix{T}, b::Vector{T}, c::Vector{T}, α::Matrix{T}, β::Vector{T}) where {T}
    CoefficientsARK{T}(name, order, length(b), length(β), a, b, c, α, β)
end

Base.hash(tab::CoefficientsARK, h::UInt) = hash(tab.o, hash(tab.s, hash(tab.r, hash(tab.a, hash(tab.b, hash(tab.c, hash(tab.α, hash(tab.β, hash(:CoefficientsARK, h)))))))))

Base.:(==)(tab1::CoefficientsARK, tab2::CoefficientsARK) = (tab1.o == tab2.o
                                                         && tab1.s == tab2.s
                                                         && tab1.r == tab2.r
                                                         && tab1.a == tab2.a
                                                         && tab1.b == tab2.b
                                                         && tab1.c == tab2.c
                                                         && tab1.α == tab2.α
                                                         && tab1.β == tab2.β)

"Print additive Runge-Kutta coefficients."
function Base.show(io::IO, tab::CoefficientsARK)
    print(io, "Additive Runge-Kutta Coefficients ", tab.name, "with ", tab.s, " internal stages, ", tab.r, " projective stages and order ", tab.o)
    print(io, "  a = ", tab.a)
    print(io, "  b = ", tab.b)
    print(io, "  c = ", tab.c)
    print(io, "  α = ", tab.α)
    print(io, "  β = ", tab.β)
end

"Holds the coefficients of a projective Runge-Kutta method."
struct CoefficientsPRK{T} <: AbstractCoefficients{T}
    @HeaderCoefficientsARK
    @CoefficientsPRK

    function CoefficientsPRK{T}(name,o,s,r,a,c,α) where {T}
        @assert T <: Real
        @assert isa(name, Symbol)
        @assert isa(o, Integer)
        @assert isa(s, Integer)
        @assert isa(r, Integer)
        @assert s > 0 "Number of stages s must be > 0"
        @assert r > 0 "Number of stages r must be > 0"
        @assert r==size(a,1)==size(α,1)==size(α,2)==length(c)
        @assert s==size(a,2)
        new(name,o,s,r,a,c,α)
    end
end

function CoefficientsPRK(name::Symbol, order::Int, a::Matrix{T}, c::Vector{T}, α::Matrix{T}) where {T}
    CoefficientsPRK{T}(name, order, size(a,2), length(c), a, c, α)
end

Base.hash(tab::CoefficientsPRK, h::UInt) = hash(tab.o, hash(tab.s, hash(tab.r, hash(tab.a, hash(tab.c, hash(tab.α, hash(:CoefficientsPRK, h)))))))

Base.:(==)(tab1::CoefficientsPRK, tab2::CoefficientsPRK) = (tab1.o == tab2.o
                                                         && tab1.s == tab2.s
                                                         && tab1.r == tab2.r
                                                         && tab1.a == tab2.a
                                                         && tab1.c == tab2.c
                                                         && tab1.α == tab2.α)

"Print projective Runge-Kutta coefficients."
function Base.show(io::IO, tab::CoefficientsPRK)
    print(io, "Projective Runge-Kutta Coefficients ", tab.name, "with ", tab.s, " internal stages, ", tab.r, " projective stages and order ", tab.o)
    print(io, "  a = ", tab.a)
    print(io, "  c = ", tab.c)
    print(io, "  α = ", tab.α)
end

"Holds the multiplier Runge-Kutta coefficients."
struct CoefficientsMRK{T}
    name::Symbol
    r::Int
    b::Vector{T}
    c::Vector{T}

    function CoefficientsMRK{T}(name,r,b,c) where {T}
        @assert T <: Real
        @assert isa(name, Symbol)
        @assert isa(r, Integer)
        @assert r > 0 "Number of stages r must be > 0"
        @assert r==length(b)==length(c)
        new(name,r,b,c)
    end
end

function CoefficientsMRK(name::Symbol, b::Vector{T}, c::Vector{T}) where {T}
    CoefficientsMRK{T}(name, length(c), b, c)
end

Base.hash(tab::CoefficientsMRK, h::UInt) = hash(tab.r, hash(tab.b, hash(tab.c, hash(:CoefficientsMRK, h))))

Base.:(==)(tab1::CoefficientsMRK, tab2::CoefficientsMRK) = (tab1.r == tab2.r
                                                         && tab1.b == tab2.b
                                                         && tab1.c == tab2.c)

"Print multiplier Runge-Kutta coefficients."
function Base.show(io::IO, tab::CoefficientsMRK)
    print(io, "Multiplier Runge-Kutta coefficients ", tab.name, "with ", tab.r, " projective stages")
    print(io, "  b = ", tab.b)
    print(io, "  c = ", tab.c)
end


text:
# BFS simple algorithm using abhinavmehndiratta SuiteSparseGraphBLAS Julia wrapper
# (can be found at https://github.com/abhinavmehndiratta/SuiteSparseGraphBLAS.jl)

using SuiteSparseGraphBLAS, GraphBLASInterface

"""
    function GraphBLAS_bfs!(A, v, source)

Compute BFS on the graph represented by `A` starting from `source`
using GraphBLAS. The output is written in v.

# Arguments
- `A::GrB_Matrix{T}`: the matrix representing input graph
- `v::GrB_Vector{U}, U<:Integer`: the vector that will contain BFS output
- `source::ZeroBasedIndex`: the source from which BFS starts

# Examples
```julia-repl

julia> GrB_init(GrB_NONBLOCKING);

julia> dim = 100;

julia> A = GrB_Matrix(Int64, dim, dim);

julia> GraphBLAS_sparse_square_matrix!(A, dim);

julia> v = GrB_Vector(Int64, dim);

julia> GraphBLAS_bfs!(A, v, ZeroBasedIndex(0));

julia> dropzeros!(v);

julia> GrB_wait();
```
"""
function GraphBLAS_bfs!(
    A::GrB_Matrix{T},
    v::GrB_Vector{U},
    source::ZeroBasedIndex
    ) where {T, U<:Integer}

    n = GrB_Matrix_ncols(A)

    GrB_Matrix_nrows(A) == n || throw(DimensionMismatch("Input matrix must be square"))
    GrB_Vector_size(v) == n || throw(DimensionMismatch("Input matrix and vector dimension mismatch"))

    # Create and initialize GrB_Descriptor
    desc_rc = GrB_Descriptor(GrB_OUTP => GrB_REPLACE, GrB_MASK => GrB_SCMP)

    # Define the semiring that will be used
    semiring = GxB_LOR_LAND_BOOL

    # Make input vector dense
    GrB_Vector_assign(v, GrB_NULL, GrB_NULL, 0, GrB_ALL, n, GrB_NULL)

    # Create a boolean vector; set q(source) = true
    q = GrB_Vector{Bool}()
    GrB_Vector_new(q, GrB_BOOL, n)

    GrB_Vector_setElement(q, true, source)

    # BFS algorithm
    for  level = 1:n
        # v<q> = level
        GrB_Vector_assign(v, q, GrB_NULL, level, GrB_ALL, n, GrB_NULL)

        if (GrB_Vector_nvals(q) == 0)
            break
        end

        # q<!v> = q ||.&& A; updates the frontier
        GrB_vxm(q, v, GrB_NULL, GxB_LOR_LAND_BOOL, q, A, desc_rc)
    end

    # free workspace
    GrB_Semiring_free(semiring)
    GrB_Descriptor_free(desc_rc)
    GrB_Vector_free(q)

    # Force pending operations to complete and return
    GrB_wait()

    return (GrB_SUCCESS)
end


text:
@testset "Free Memory Blocks Tests" begin

    model = ALifeBenchmark.TierraModel(
        ALifeBenchmark.TierrianInstruction[
            ALifeBenchmark.NOP_0()
        ]
    )

    # test size of allocated memory

    memory_allocated::UInt64 = 1

    @test model.free_blocks == [
        ALifeBenchmark.FreeMemoryBlock(1, 2^16-1)
    ]
    @test model.used_memory == memory_allocated

    size::UInt64 = 0
    result::Int64 = 0

    for _ in 1:10_000
        # first: allocate some memory

        size = rand(1:10000)

        result = ALifeBenchmark.allocate_free_memory!(model, size)

        if result != -1
            memory_allocated += size
        end

        @test model.used_memory == memory_allocated

        # second: try to clear some memory

        size = rand(1:10000)

        memory_allocated -= ALifeBenchmark.clear_memory!(model, rand(0:(2^16-size-1)), size)
        @test model.used_memory == memory_allocated

        # test free blocks

        free_block_sum = 0
        for free_block in model.free_blocks
            free_block_sum += free_block.length
        end
        @test free_block_sum == ALifeBenchmark.get_free_memory_size(model)
    end

    # test if free blocks are exclusive and non-neighboring

    free = falses(2^16)

    non_exclusive = false
    no_neighbors = false

    for free_block in model.free_blocks
        for i in free_block.start_address:(free_block.start_address - 1 + free_block.length)
            if free[i]
                non_exclusive = true
            end
            free[i] = true
        end

        if 0 < free_block.start_address
            if free[free_block.start_address - 1]
                println(free_block)
                no_neighbors = true
            end
        end
        if free_block.start_address + free_block.length < 2^16
            if free[free_block.start_address + free_block.length]
                println(free_block)
                no_neighbors = true
            end
        end
    end

    @test !non_exclusive
    @test !no_neighbors
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
###############################################################################
#
# Dynamically change number of workers
#
###############################################################################
"""
    dynamic_workers!(nthread::Int)

Change the number of workers according to CPU_THREADS, given
- `nthread` Number of thread required
"""
function dynamic_workers!(nthread::Int)
    MaxThreads = Sys.CPU_THREADS;

    # if        no worker yet, and nthread <= MaxThreads, hire nthread
    # elseif    no worker yet, but nthread > MaxThreads, hire MaxThreads
    # elseif    some workers already, and nthread <= MaxThreads, hire more
    # elseif    some workers already, and nthread > MaxThreads, hire more
    # else      workers is more than expected, remove the extra
    if (length(workers())==1) && (workers()[1]==1) && (nthread<=MaxThreads)
        addprocs(nthread, exeflags="--project");
    elseif (length(workers())==1) && (workers()[1]==1) && (nthread>MaxThreads)
        addprocs(MaxThreads, exeflags="--project");
    elseif length(workers())<nthread && (nthread<=MaxThreads)
        addprocs(nthread-length(workers()), exeflags="--project");
    elseif length(workers())<MaxThreads && (nthread>MaxThreads)
        addprocs(MaxThreads-length(workers()), exeflags="--project");
    else
        _to_remove = workers()[(nthread+1):end];
        rmprocs(_to_remove...);
    end

    # load the module into each worker
    @everywhere Base.MainInclude.eval(using GriddingMachineOld);

    return nothing
end


text:
using Distributed: pmap, addprocs, nworkers, @everywhere
@everywhere using Infiltrator
const NUM_DTW = 6
const NWORKERS = 3

if nworkers() < NWORKERS
    new_workers = NWORKERS - nworkers() + 1
    println("adding $new_workers workers")
    # number of threads must be set vie env variable...
    ENV["JULIA_NUM_THREADS"] = 2 * NUM_DTW
    # note the `--project=.`
    addprocs(new_workers; exeflags="--project=.")
end

@everywhere begin
    using ProgressMeter: @showprogress
    using PyCall
    using Statistics: mean
    using Printf: @sprintf
    include("alignment_fede.jl")

    const DATASET_LEN = 0.05
    const MIN_NOTES = 0
    const MAX_NOTES = 2000
    const MAX_DURATION = 120
    const STEP_WEIGHTS = 13
    const NWEIGHT = $NUM_DTW * STEP_WEIGHTS
    const TOT_WEIGHTS = NWEIGHT + 2
end

@everywhere py"""
import sys
sys.path.insert(0, ".")
import skopt
import numpy as np
from alignment import fede_tuning
from alignment.evaluate_midi2midi import MissingExtraManager, Result
from alignment.asmd.asmd import asmd, dataset_utils
from skopt.space.space import Integer
"""

@everywhere const full_dataset = py"asmd.Dataset()"
@everywhere const dataset, _ = py"dataset_utils.choice"(full_dataset,
                                                        p=[DATASET_LEN, 1 - DATASET_LEN],
                                                        random_state=1992)

@everywhere function objective(hparams)
    α = convert(Int64, hparams[1])
    β = convert(Int64, hparams[2])
    # python was major mode, here each column is a set of weights, instead
    sp_weights = reshape(convert(Vector{Float64}, hparams[3:end]), STEP_WEIGHTS, :)
    # display("text/plain", sp_weights)

    L = convert(Int64, py"len"(dataset))
    println("Size of dataset: $L")
    res = Vector{Float64}(undef, 0)
    @showprogress 1 "Evaluating..." for i in 0:(L - 1)
        mem = py"MissingExtraManager"(dataset, i)
        if mem.num_notes > MAX_NOTES || mem.num_notes < MIN_NOTES
            continue
        else
            push!(res, eval_fede(mem, α, β, sp_weights))
        end
    end
    f1m = mean(res)
    println(@sprintf("Result: %.2e", f1m))
    println("N. of tests are $(length(res))")
    return 1 - f1m
end

@everywhere @views function eval_fede(mem, α, β, sp_weights)::Float64
    score = mem.get_score()
    perfm = mem.get_perfm()
    fede_match, _, _ = FedeAlignment.get_matching_notes(score, perfm, α, β, sp_weights, 0.0)
    if length(fede_match) > 0
        #transpose
        fede_match = permutedims(fede_match)
        mask = trues(size(fede_match)...)
        for i in 1:2
            # compute argmin and argmax
            m = argmin(fede_match[i, :])
            M = argmax(fede_match[i, :])
            # remove their column
            mask[i, m] = false
            mask[i, M] = false
        end
        fede_match = hcat(fede_match[1, mask[1, :]], fede_match[2, mask[2, :]])
    end
    fede_res = mem.evaluate(fede_match .- 2, "fede", 0.0)
    f1m = convert(Float64, fede_res.match_fmeasure)
    return f1m
end

struct EarlyStop
    range::Float64
    patience::Int64
end

function (self::EarlyStop)(res::Dict)
    m = minimum(res["func_vals"])
    vals = res["func_vals"] .- m
    if count(<(self.range), vals) > self.patience
        println("Early-stopping!")
        return true
    else
        return false
    end
end

function callbacks(res::Dict; early_stop::Union{Nothing,EarlyStop}=nothing)::Bool
    # early stopping
    if early_stop !== nothing
        return early_stop(res)
    end
    return false
end

function log(res, nworkers)
    println("----------------")
    println("""Total number of calls: $(length(res["func_vals"]))""")
    println("Values obtained: ")
    display("text/plain", res["func_vals"][(end - nworkers + 1):end])
    return println(@sprintf("Best val so far: %.2e", res["fun"]))
end

function tell!(optimizer, xx, yy, nworkers)
    # xx and yy are Julia vectors, that can be translated to numpy arrays or
    # lists; python lists are translated to Julia arrays. Moreover,
    # sometimes only the outer vector is translated into a list, while the
    # inners remain arrays, or the opposite.
    # so, let's do everything in python
    py"""
    xx = [x.tolist() if type(x) is np.ndarray else x for x in $xx]
    yy = [y.tolist() if type(y) is np.ndarray else y for y in $yy]
    res = $optimizer.tell(xx, yy)
    """
    res = py"res"
    log(res, nworkers)
    return res
end

"""
Perform `n_iter * n_jobs` calls and returns the sklearn ResultObject
"""
function optimize(optimizer::PyObject, fn::Function, n_iter; nworkers=NWORKERS,
                  early_stop::Union{Nothing,EarlyStop}=nothing,
                  x0::Vector{Vector{T}}=[]) where {T}

    for x in x0
        println("Testing provided initial point x=")
        println(x)
        y = fn(x)

        res = tell!(optimizer, [x], [y], 1)
    end
    for n in 1:n_iter
        println("----------------")
        println("Iteration number $n")

        if nworkers == 1
            xx = [optimizer.ask(n_points=nworkers)[1, :]]
        else
            xx = optimizer.ask(n_points=nworkers)
        end
        println("Testing x=")
        display("text/plain", xx)
        if nworkers == 1
            yy = [fn(xx[1])]
        else
            yy = pmap(fn, eachrow(xx))
        end

        res = tell!(optimizer, xx, yy, nworkers)
        if callbacks(res; early_stop=early_stop)
            break
        end
    end

    return optimizer.get_result()

end

function main()
#! format: off
    # x0 = [
    #     [
    #         20, 13,
    #         1, 3,
    #         1, 2,
    #         1, 1,
    #         1, 1, 3,
    #         1, 1, 2,
    #         2, 1,
    #         1, 3, 1,
    #         1, 2, 1,
    #         3, 1,
    #         1,

    #         1, 2,
    #         1, 3,
    #         1, 2,
    #         1, 2, 3,
    #         1, 2, 2,
    #         1, 1,
    #         1, 3, 2,
    #         1, 2, 1,
    #         2, 1,
    #         1,

    #         1, 2,
    #         1, 1,
    #         1, 3,
    #         1, 1, 3,
    #         1, 2, 2,
    #         1, 1,
    #         1, 3, 2,
    #         1, 2, 1,
    #         2, 1,
    #         1,

    #         2, 1,
    #         2, 1,
    #         2, 2,
    #         2, 1, 1,
    #         1, 2, 1,
    #         1, 2,
    #         1, 1, 2,
    #         1, 2, 1,
    #         1, 1,
    #         1,

    #         1, 3,
    #         1, 3,
    #         2, 1,
    #         2, 1, 3,
    #         1, 2, 1,
    #         1, 2,
    #         1, 1, 2,
    #         1, 1, 1,
    #         2, 1,
    #         1,

    #     ],
    # ]
    x0 = [
        [
            18, 9,
            1, 3, 1, 0,
            1, 2, 1, 0,
            1, 1, 1, 0,
            1,

            1, 1, 1, 0,
            1, 1, 1, 0,
            1, 1, 2, 0,
            1,

            1, 1, 1, 0,
            1, 2, 1, 0,
            1, 1, 1, 0,
            1,

            2, 1, 1, 0,
            1, 1, 1, 0,
            1, 1, 1, 0,
            1,

            1, 1, 1, 0,
            1, 1, 1, 0,
            1, 1, 2, 0,
            2,

            1, 1, 1, 0,
            1, 2, 1, 0,
            1, 1, 1, 0,
            2,
        ],
    ]
    optimizer = py"skopt.optimizer.Optimizer"(
        vcat([
            py"Integer"(10, 30, name="alpha", transform="normalize"),
            py"Integer"(5, 20, name="beta", transform="normalize")
        ], [
            py"Integer"(0, 7, name="step$i", transform="normalize") for i in 1:NWEIGHT
        ]);
        # base_estimator=SKDNGO(False),
        # base_estimator=py"skopt.learning.RandomForestRegressor"(n_estimators=10),
        # base_estimator="GBRT",
        base_estimator=py"skopt.learning.ExtraTreesRegressor"(
            n_estimators=1, min_samples_leaf=5, warm_start=true),
        # base_estimator="ET",
        # n_calls=4000,
        n_initial_points=100, random_state=1750,
        acq_func="EI",
        acq_func_kwargs=Dict("xi" => 0.1),
        # kappa=0.7,
        # kappa=BOFactor(update_kappa),
        acq_optimizer="sampling",
        # n_points=10^4,
        # verbose=true,
        # callback=[
        #     # py"skopt.callbacks.CheckpointSaver"("fede_model_dngo.pkl"),
        #     py"skopt.callbacks.DeltaYStopper"(0.01, 40),
        # ],
        # x0=x0,
        n_jobs=-1)
#! format: on
    res = optimize(optimizer, objective, 2000; early_stop=EarlyStop(0.01, 40), x0=x0)

    println("----------------")
    println(@sprintf("Best value: %.2e", res["fun"]))
    return println("""Point: $(res["x"])""")
end


text:
@inline function ExaTronKernels.daxpy(n::Int,da::Float64,
                       dx::ROCDeviceArray{Float64,1},incx::Int,
                       dy::ROCDeviceArray{Float64,1},incy::Int)
    tx = workitemIdx().x

    if tx <= n
        @inbounds dy[tx] = dy[tx] + da*dx[tx]
    end
    AMDGPU.sync_workgroup()

    return
end