text:
"""
### Cutout sub-regions of the data base of DataSetType
- select shape of a shell-region
- select size of a region (w/o intersecting cells)
- give the spatial center (with units) of the data relative to the full box
- relate the coordinates to a direction (x,y,z)
- inverse the selected region

```julia
shellregion(dataobject::DataSetType, shape::Symbol=:cylinder;
            radius::Array{<:Real,1}=[0.,0.],  # cylinder, sphere;
            height::Real=0.,                  # cylinder
            direction::Symbol=:z,                # cylinder

            center::Array{<:Any,1}=[0., 0., 0.],   # all
            range_unit::Symbol=:standard,  # all
            cell::Bool=true,                        # hydro and gravity
            inverse::Bool=false,                    # all
            verbose::Bool=verbose_mode)             # all
```

#### Arguments
##### Required:
- **`dataobject`:** needs to be of type: "DataSetType"
- **`shape`:** select between region shapes: :cuboid, :cylinder/:disc, :sphere
##### Predefined/Optional Keywords:
**For cylindrical shell-region, related to a given center:**
- **`radius`:** the inner and outer radius of the shell in units given by argument `range_unit` and relative to the given `center`
- **`height`:** the hight above and below a plane [-height, height] in units given by argument `range_unit` and relative to the given `center`
- **`direction`:** todo

**For spherical shell-region, related to a given center:**
- **`radius`:** the inner and outer radius of the shell in units given by argument `range_unit` and relative to the given `center`

**Keywords related to all region shapes**
- **`range_unit`:** the units of the given ranges: :standard (code units), :Mpc, :kpc, :pc, :mpc, :ly, :au , :km, :cm (of typye Symbol) ..etc. ; see for defined length-scales viewfields(info.scale)
- **`center`:** in units given by argument `range_unit`; by default [0., 0., 0.]; the box-center can be selected by e.g. [:bc], [:boxcenter], [value, :bc, :bc], etc..
- **`inverse`:** inverse the region selection = get the data outside of the region
- **`cell`:** take intersecting cells of the region boarder into account (true) or only the cells-centers within the selected region (false)
- **`verbose`:** print timestamp, selected vars and ranges on screen; default: set by the variable `verbose_mode`
"""
function shellregion(dataobject::DataSetType, shape::Symbol=:cylinder;
            radius::Array{<:Real,1}=[0.,0.],  # cylinder, sphere;
            height::Real=0.,                  # cylinder
            direction::Symbol=:z,                # cylinder

            center::Array{<:Any,1}=[0., 0., 0.],   # all
            range_unit::Symbol=:standard,  # all
            cell::Bool=true,                        # hydro and gravity
            inverse::Bool=false,                    # all
            verbose::Bool=verbose_mode)             # all

    # subregion = wrapper over all subregion shell functions

    if shape == :cylinder || shape == :disc
        if typeof(dataobject) == HydroDataType
            return shellregioncylinder(dataobject,
                                        radius=radius,
                                        height=height,
                                        center=center,
                                        range_unit=range_unit,
                                        direction=direction,
                                        cell=cell,
                                        inverse=inverse,
                                        verbose=verbose)
        else
            return shellregioncylinder(dataobject,
                                        radius=radius,
                                        height=height,
                                        center=center,
                                        range_unit=range_unit,
                                        direction=direction,
                                        inverse=inverse,
                                        verbose=verbose)
        end

    elseif shape == :sphere
        if typeof(dataobject) == HydroDataType
            return shellregionsphere(  dataobject,
                                        radius=radius,
                                        center=center,
                                        range_unit=range_unit,
                                        cell=cell,
                                        inverse=inverse,
                                        verbose=verbose)
        else
            return shellregionsphere(  dataobject,
                                        radius=radius,
                                        center=center,
                                        range_unit=range_unit,
                                        inverse=inverse,
                                        verbose=verbose)

        end

    end

end


text:
module Hotwire

export Thermistor, Resistor, AbstractResistor, temperature, resistance
export CTASensor, Wire, resistor, optemperature, overtemp, overheat_ratio, gain
export AbstractThermalAnemometer
export CCASensor, current
export tempcorr
export calibr_curve, CalibrCurve, ExtrapolateFit
export HWKingFit, BlendFit

export AbstractProbe, AbstractProbe1d, AbstractProbe2d, AbstractProbe3d
export HWCable, HWSupport, HWBridge
export Probe1d, Probe2d, Probe3d

include("resistor.jl")
include("wire.jl")
include("tempcorr.jl")
include("probes.jl")
include("calibr.jl")
end # module

ext: jl


text:
###############################
### WUE metrics calculation ###
###############################

#' Water-Use Efficiency Metrics
#'
#' Calculation of various water use efficiency (WUE) metrics.
#'
#' - data      DataFrame or matrix containing all required variables
#' - GPP       Gross primary productivity (umol CO2 m-2 s-1)
#' - NEE       Net ecosystem exchange (umol CO2 m-2 s-1)
#' - LE        Latent heat flux (W m-2)
#' - VPD       Vapor pressure deficit (kPa)
#' - Tair      Air temperature (deg C)
#' - constants Cmol - molar mass of carbon (kg mol-1)
#'                  umol2mol - conversion micromole (umol) to mole (mol)
#'                  kg2g - conversion kilogram (kg) to gram (g)
#'
#' # Details
 the following metrics are calculated:
#'
#'          Water-use efficiency (WUE):
#'
#'            ``WUE = GPP / ET``
#'
#'          Water-use efficiency based on NEE (WUE_NEE):
#'
#'            ``WUE_NEE = NEE / ET``
#'
#'          Inherent water-use efficiency (IWUE; Beer et al. 2009):
#'
#'            ``IWUE = (GPP * VPD) / ET``
#'
#'          Underlying water-use efficiency (uWUE; Zhou et al. 2014):
#'
#'            ``uWUE= (GPP * sqrt(VPD)) / ET``
#'
#'          All metrics are calculated based on the median of all values. E_g.
#'          WUE = median(GPP/ET,na_rm=true)
#'
#' # Value
 a named vector with the following elements:
#'         - WUE: Water-use efficiency (gC (kg H20)-1)
#'         - WUE_NEE: Water-use efficiency based on NEE (gC (kg H20)-1)
#'         - IWUE: Inherent water-use efficiency (gC kPa (kg H20)-1)
#'         - uWUE: Underlying water-use efficiency (gC kPa^0.5 (kg H20)-1)
#'
#' # Note
#' Units for VPD can also be hPa. Units change accordingly.
#'       WUE_NEE is calculated based on the absolute value of NEE (the sign convention does not matter here).
#'
#' #References
#' Beer, C., et al., 2009: Temporal and among-site variability of inherent
#'             water use efficiency at the ecosystem level. Global Biogeochemical Cycles 23, GB2018.
#'
#'             Zhou, S., et al., 2014: The effect of vapor pressure deficit on water
#'             use efficiency at the sub-daily time scale. Geophysical Research Letters 41.
#'
#' #See also
#' [`stomatal_slope`](@ref) for a measure of intrinsic WUE
#'
#' ```@example; output = false
#' ```
#' ## filter data for dry periods and daytime at DE-Tha in June 2014
#' DE_Tha_Jun_2014_2 = filter_data(DE_Tha_Jun_2014,quality_control=false,
#'                                  vars_qc=c("Tair","precip","VPD","H","LE"),
#'                                  filter_growseas=false,filter_precip=true,
#'                                  filter_vars=c("Tair","PPFD","ustar"),
#'                                  filter_vals_min=c(5,200,0.2),
#'                                  filter_vals_max=c(missing,missing,missing),NA_as_invalid=true,
#'                                  quality_ext="_qc",good_quality=c(0,1),
#'                                  missing_qc_as_bad=true,GPP="GPP",doy="doy",
#'                                  year="year",tGPP=0.5,ws=15,min_int=5,precip="precip",
#'                                  tprecip=0.1,precip_hours=24,records_per_hour=2)
#'
#' ## calculate WUE metrics in the filtered periods
#' WUE_metrics(DE_Tha_Jun_2014_2)
#'
#' @importFrom stats median
"""
"""
function WUE_metrics(data,GPP="GPP",NEE="NEE",LE="LE",VPD="VPD",Tair="Tair",
                        constants=BigleafConstants())

  check_input(data,list(GPP,NEE,LE,VPD,Tair))

  ET  = LE_to_ET(LE,Tair)                 # kg H2O m-2 s-1
  GPP = (GPP * constants.umol2mol] * constants[:Cmol]) * constants[:kg2g  # gC m-2 s-1
  NEE = (NEE * constants.umol2mol] * constants[:Cmol]) * constants[:kg2g  # gC m-2 s-1

  WUE     = median(GPP/ET,na_rm=true)
  WUE_NEE = median(abs(NEE)/ET,na_rm=true)
  IWUE    = median((GPP*VPD)/ET,na_rm=true)
  uWUE    = median((GPP*sqrt(VPD))/ET,na_rm=true)

  return(c(WUE=WUE,WUE_NEE=WUE_NEE,IWUE=IWUE,uWUE=uWUE))
end


text:
# Test for parallel external sort
#
# Call this file using
#
#     $ julia -p 2 ./pes-test.jl

using Base.Test

include("pes.jl")
include("helper.jl")

@testset "Parallel External Sort Tests" begin

    ############################ File with 12 integers ############################

    create_file(12)

    # Basic test
    parallel_external_sort(2, 4)
    @test check_sort()

end

max_stars_repo_name: rgaiacs/external-sort

max_issues_repo_name: rgaiacs/external-sort

max_forks_repo_name: rgaiacs/external-sort


text:
@testset "eval" begin
    include("../scripts/eval.jl")
    args_common = ["num_instances=1", "time_limit=1", "root=./local/"]
    yaml_file = "../scripts/config/eval/point2d.yaml"
    main(yaml_file, args_common...)
end


text:
using BenchmarkTools

using QOI
using PNGFiles
using TestImages
using FixedPointNumbers
using ColorTypes

struct TrialData
    load_time::Float64
    save_time::Float64
    size::Float64
end

struct BenchmarkResult
    file::String
    dim::Tuple{Int, Int}
    size::Int # bytes
    png_data::TrialData
    qoi_data::TrialData
end

function save_png_qoi(img)
    t_png = tempname()
    PNGFiles.save(t_png, img)
    t_qoi = tempname()
    QOI.qoi_encode(t_qoi, img)
    return t_png, t_qoi
end

function run_benchmarks(png_images)
    results = BenchmarkResult[]

    minimum_trial(trial) = minimum(trial.times)

    for (i, png_name) in enumerate(png_images)
        @info "Benchmarking $png_name, $i / $(length(png_images))"
        png_image = testimage(png_name)

        # QOI only supports 24-bit colors so convert to that
        T = eltype(png_image)
        if T <: TransparentColor
            png_image = RGBA{N0f8}.(png_image)
        else
            png_image = RGB{N0f8}.(png_image)
        end

        # Load
        f_png, f_qoi = save_png_qoi(png_image)
        b_png_load = @benchmark PNGFiles.load($f_png)
        b_qoi_load = @benchmark QOI.qoi_decode($f_qoi)

        # Save
        t_png = tempname()
        b_png_save = @benchmark PNGFiles.save($t_png, $png_image)
        t_qoi = tempname()
        b_qoi_save = @benchmark QOI.qoi_encode($t_qoi, $png_image)

        # Sizes
        orig_size = sizeof(png_image)
        png_size = filesize(t_png)
        qoi_size = filesize(t_qoi)

        # Sanity check implementation
        @assert QOI.qoi_decode(f_qoi) == PNGFiles.load(f_png)
        @assert QOI.qoi_decode(t_qoi) == PNGFiles.load(t_png)

        png_data = TrialData(minimum_trial(b_png_load), minimum_trial(b_png_save), png_size)
        qoi_data = TrialData(minimum_trial(b_qoi_load), minimum_trial(b_qoi_save), qoi_size)

        result = BenchmarkResult(png_name, size(png_image), orig_size, png_data, qoi_data)
        push!(results, result)
    end

    return results
end

function render_table(io::IO, results::Vector{BenchmarkResult})
    for result in results
        println(io, "`", result.file, "` ", result.dim[1], "x", result.dim[2], " ", round(result.size / 1024; digits=2), "kB")
        println(io)

        mega_pixels = result.size / (1024^2)

        png_dec_ms = round(result.png_data.load_time / 1e6; digits=2)
        qoi_dec_ms = round(result.qoi_data.load_time / 1e6; digits=2)
        png_enc_ms = round(result.png_data.save_time / 1e6; digits=2)
        qoi_enc_ms = round(result.qoi_data.save_time / 1e6; digits=2)
        png_dec_mpps = round(mega_pixels / (result.png_data.load_time/1e9); digits=2)
        qoi_dec_mpps = round(mega_pixels / (result.qoi_data.load_time/1e9); digits=2)
        png_enc_mpps = round(mega_pixels / (result.png_data.save_time/1e9); digits=2)
        qoi_enc_mpps = round(mega_pixels / (result.qoi_data.save_time/1e9); digits=2)
        png_size_kb = round(result.png_data.size / 1024; digits=2)
        qoi_size_kb = round(result.qoi_data.size / 1024; digits=2)
        png_rate = round(result.png_data.size / result.size * 100; digits=1)
        qoi_rate = round(result.qoi_data.size / result.size * 100; digits=1)

        print(io, """
        |           | decode ms | encode ms | decode mpps | encode mpps | size kb | rate |
        |:-----------|----------:|----------:|------------:|------------:|--------:|-----:|
        |  PNGFiles  |$png_dec_ms|$png_enc_ms|$png_dec_mpps|$png_enc_mpps|$png_size_kb |$png_rate%|
        |  QOI       |$qoi_dec_ms|$qoi_enc_ms|$qoi_dec_mpps|$qoi_enc_mpps|$qoi_size_kb |$qoi_rate%|

        ----------------------

        """)
    end
end

png_images = filter(endswith(".png"), TestImages.remotefiles)

results = run_benchmarks(png_images)
render_table(stdout, results)


text:
# v0.6

using Dates

function easter(year::Int)::Date
    a = rem(year, 19)
    b, c = divrem(year, 100)
    d = rem(19a + b - div(b, 4) - div(b - div(b + 8, 25) + 1, 3) + 15, 30)
    e = rem(32 + 2 * rem(b, 4) + 2 * div(c, 4) - d - rem(c, 4), 7)
    f = d + e - 7 * div(a + 11d + 22e, 451) + 114
    month, day = divrem(f, 31)
    day += 1
    return Date(year, month, day)
end

function holiday_values(year::Int)::Dict{String,Date}
    offsets = Dict{String,Day}("Easter" => Day(0), "Ascension" => Day(39),
    "Pentecost" => Day(49), "Trinity" => Day(56), "Corpus" => Day(60))
    easterdate = easter(year)
    rst = Dict{String,Date}(holiday => easterdate + days for (holiday, days) in offsets)
    return rst
end

function holiday2str(year::Int)::String
    function holiday2str(holiday::String, date::Date)::String
        dayname, daynumb, month = dayabbr(date), day(date), monthabbr(date)
        return @sprintf "%s: %s %2i %s" holiday dayname daynumb month
    end
    cal = holiday_values(year)
    rst = join(collect(holiday2str(hday, cal[hday]) for hday in
    ("Easter", "Ascension", "Pentecost", "Trinity", "Corpus")), ", ")
    return "$year -> " * rst
end

println("\nChristian holidays, related to Easter, for each centennial from 400 to 2100 CE:")
for yr in 400:100:2200
    println(holiday2str(yr))
end

println("\nChristian holidays, related to Easter, for years from 2010 to 2020 CE:")
for yr in 2010:2020
    println(holiday2str(yr))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using ErlPort.Exceptions
using ErlPort.ErlTerms.Decode

include("Testing.jl")

# XXX Do we need to test REFERENCE_EXT, PORT_EXT, PID_EXT, NEW_REFERENCE_EXT, FUN_EXT, NEW_FUN_EXT, EXPORT_EXT

# data for use by tests
badsizecompdata = vcat(b"\x83P\0\0\0\x16\x78\xda\xcb\x66\x10\x49\xc1\2\0\x5d",
                       b"\x60\x08\x50")
compdata1 = b"\x83P\0\0\0\x17\x78\xda\xcb\x66\x10\x49\xc1\2\0\x5d\x60\x08\x50"
compdata2 = vcat(b"\x83P\0\0\0\x17\x78\xda\xcb\x66\x10\x49\xc1\2\0\x5d\x60",
                 b"\x08\x50tail")

tuplepart = b"d\0\1ad\0\1bd\0\1cd\0\1dd\0\1ed\0\1fd\0\1gd\0\1hd\0\1id\0\1jd\0\1kd\0\1ld\0\1md\0\1nd\0\1od\0\1p"
largetuple1 = vcat(b"i\0\0\1\0",
                   tuplepart, tuplepart, tuplepart, tuplepart, tuplepart, tuplepart,
                   tuplepart, tuplepart, tuplepart, tuplepart, tuplepart, tuplepart,
                   tuplepart, tuplepart, tuplepart, tuplepart)

smallbigintmax = vcat(b"n\xff\0", fill(uint8(255), 255))
largebigintmin = vcat(b"o\0\0\1\0\0", fill(uint8(0), 255), b"\1")

# tests for supporting functions
testcase() do
    @test size1unpack(b"\x0a")             == 0x0a
    @test size1unpack(b"\xaa")             == 0xaa
    @test size2unpack(b"\x0a\x0b")         == 0x0a0b
    @test size2unpack(b"\xaa\xbb")         == 0xaabb
    @test size4unpack(b"\x0a\x0b\x0c\x0d") == 0x0a0b0c0d
    @test size4unpack(b"\xaa\xbb\xcc\xdd") == 0xaabbccdd

    expected = [107,0,20,100,100,100,100,100,100,100,100,100,100,
                100,100,100,100,100,100,100,100,100,100]
    @test decompressterm(compdata1) == expected
end

# basic decode errors
testcase() do
    @test_throws UnsupportedType decode("")
    @test_throws IncompleteData decode(b"")
    @test_throws UnknownProtocolVersion decode(b"\0")
    @test_throws IncompleteData decode(b"\x83")
end

# decode empty list (that is in fact a nil) (NIL_EXT)
testcase() do
    @test_throws IncompleteData decodenil(b"")
    #@test_throws IncompleteData decodenil(b"\0") # XXX decode functions don't check term tag
    @test decode(b"\x83j") == ([], b"")
    @test decodenil(b"j") == ([], b"")
    @test decodenil(b"jtail") == ([], b"tail")
end

# decode atoms (ATOM_EXT)
testcase() do
    @test_throws IncompleteData decode(b"\x83d")
    @test_throws IncompleteData decode(b"\x83d\0")
    @test_throws IncompleteData decode(b"\x83d\0\1")
    @test decode(b"\x83d\0\0") == (symbol(""), [])
    @test decode(b"\x83d\0\0tail") == (symbol(""), b"tail")
    @test decode(b"\x83d\0\4test") == (:test, b"")
    @test decode(b"\x83d\0\4testtail") == (:test, b"tail")
    @test decode(b"\x83d\0\4true") == (true, b"")
    @test decode(b"\x83d\0\5false") == (false, b"")
    @test decode(b"\x83d\0\x09undefined") == (nothing, b"")
    @test decodeatom(b"d\0\0") == (symbol(""), [])
    @test decodeatom(b"d\0\0tail") == (symbol(""), b"tail")
    @test decodeatom(b"d\0\4test") == (:test, b"")
    @test decodeatom(b"d\0\4testtail") == (:test, b"tail")
    @test decodeatom(b"d\0\4true") == (true, b"")
    @test decodeatom(b"d\0\5false") == (false, b"")
    @test decodeatom(b"d\0\x09undefined") == (nothing, b"")
end

# decode predefined atoms (ATOM_CACHE_REF)
testcase() do
    #atomcachereftag = uint8(82)
end

# decode small atoms (SMALL_ATOM_EXT)
testcase() do
    #smallatomtag = uint8(115)
end

# decode UTF-8 atoms (ATOM_UTF8_EXT)
testcase() do
    #atomutf8tag = uint8(118)
end

# decode small, UTF-8 atoms (SMALL_ATOM_UTF8_EXT)
testcase() do
    #smallatomutf8tag = uint8(119)
end

# decode string list (STRING_EXT)
testcase() do
    @test_throws IncompleteData decodestring(b"")
    @test_throws IncompleteData decodestring(b"\0")
    @test_throws IncompleteData decodestring(b"\0\0")
    @test decode(b"\x83k\0\1\1") == (b"\1", b"")
    @test decode(b"\x83k\0\6string") == (b"string", b"")
    @test decodestring(b"k\0\0") == (b"", b"")
    @test decodestring(b"k\0\0tail") == (b"", b"tail")
    @test decodestring(b"k\0\6string") == (b"string", b"")
    @test decodestring(b"k\0\x1f\xc3\xa1rv\xc3\xadzt\xc5\xb1r\xc5\x91\x20t\xc3\xbck\xc3\xb6rf\xc3\xbar\xc3\xb3g\xc3\xa9p") == (b"árvíztűrő tükörfúrógép", b"")
end

# decode list (LIST_EXT)
testcase() do
    @test decode(b"\x83l\0\0\0\4a\1d\0\1aa\3d\0\x09undefinedj") == ([1,:a,3,nothing], b"")
    @test decodelist(b"l\0\0\0\4a\1d\0\1aa\3d\0\x09undefinedj") == ([1,:a,3,nothing], b"")
    @test decodelist(b"l\0\0\0\4a\1d\0\1aa\3d\0\x09undefinedjtail") == ([1,:a,3,nothing], b"tail")
end

# decode improper list (LIST_EXT continued)
# XXX how should we handle improper lists in Julia?

# decode small tuple (SMALL_TUPLE_EXT)
testcase() do
    @test_throws IncompleteData decodesmalltuple(b"")
    @test_throws IncompleteData decodesmalltuple(b"\0")
    @test_throws IncompleteData decodesmalltuple(b"h\1")
    @test decode(b"\x83h\2a\1a\2") == ((1,2), b"")
    @test decode(b"\x83h\3a\x0aa\x14a\x1e") == ((10,20,30), b"")
    @test decode(b"\x83h\4d\0\1ad\0\1bd\0\1cd\0\1d") == ((:a,:b,:c,:d), b"")
    @test decode(b"\x83h\2a\1a\2tail") == ((1,2), b"tail")
    @test decodesmalltuple(b"h\0") == ((), b"")
    @test decodesmalltuple(b"h\0tail") == ((), b"tail")
    # it's a real tuple that has only 1 element
    @test decodesmalltuple(b"h\1a\1") == ((1,), b"")
    @test decodesmalltuple(b"h\2a\1a\2tail") == ((1,2), b"tail")
end

# decode large tuple (LARGE_TUPLE_EXT)
testcase() do
    (lt1, tail1) = decode(vcat(b"\x83", largetuple1, b"tail"))
    @test tail1 == b"tail"

    (lt2, tail2) = decodelargetuple(largetuple1)
    @test (length(lt2), tail2) == (256, b"")
    @test lt2[1]   == :a
    @test lt2[256] == :p
end

# decode opaque object

# decode small integer (SMALL_INTEGER_EXT)
testcase() do
    @test_throws IncompleteData decode(b"\x83a")
    @test decode(b"\x83a\0") == (0, b"")
    @test decode(b"\x83a\0tail") == (0, b"tail")
    @test decode(b"\x83a\xff") == (255, b"")
    @test decode(b"\x83a\xfftail") == (255, b"tail")
    @test_throws IncompleteData decodesmallint(b"a")
    @test decodesmallint(b"a\0") == (0, b"")
    @test decodesmallint(b"a\0tail") == (0, b"tail")
    @test decodesmallint(b"a\xff") == (255, b"")
    @test decodesmallint(b"a\xfftail") == (255, b"tail")
end

# decode integer (INTEGER_EXT)
testcase() do
    @test_throws IncompleteData decode(b"\x83b")
    @test_throws IncompleteData decode(b"\x83b\0")
    @test_throws IncompleteData decode(b"\x83b\0\0")
    @test_throws IncompleteData decode(b"\x83b\0\0\0")
    @test decode(b"\x83b\0\0\0\0") == (0, [])
    @test decode(b"\x83b\0\0\0\0tail") == (0, b"tail")
    @test decode(b"\x83b\x7f\xff\xff\xff") == (2147483647, [])
    @test decode(b"\x83b\x7f\xff\xff\xfftail") == (2147483647, b"tail")
    @test decode(b"\x83b\xff\xff\xff\xff") == (-1, [])
    @test decode(b"\x83b\xff\xff\xff\xfftail") == (-1, b"tail")
    @test_throws IncompleteData decodeint(b"b")
    @test_throws IncompleteData decodeint(b"b\0")
    @test_throws IncompleteData decodeint(b"b\0\0")
    @test_throws IncompleteData decodeint(b"b\0\0\0")
    @test decodeint(b"b\0\0\0\0") == (0, [])
    @test decodeint(b"b\0\0\0\0tail") == (0, b"tail")
    @test decodeint(b"b\x7f\xff\xff\xff") == (2147483647, [])
    @test decodeint(b"b\x7f\xff\xff\xfftail") == (2147483647, b"tail")
    @test decodeint(b"b\xff\xff\xff\xff") == (-1, [])
    @test decodeint(b"b\xff\xff\xff\xfftail") == (-1, b"tail")
end

# decode binary (BINARY_EXT)
testcase() do
    @test_throws IncompleteData decode(b"\x83m")
    @test_throws IncompleteData decode(b"\x83m\0")
    @test_throws IncompleteData decode(b"\x83m\0\0")
    @test_throws IncompleteData decode(b"\x83m\0\0\0")
    @test_throws IncompleteData decode(b"\x83m\0\0\0\1")
    @test decode(b"\x83m\0\0\0\0") == (b"", b"")
    @test decode(b"\x83m\0\0\0\0tail") == (b"", b"tail")
    @test decode(b"\x83m\0\0\0\4data") == (b"data", b"")
    @test decode(b"\x83m\0\0\0\4datatail") == (b"data", b"tail")
    @test decodebin(b"m\0\0\0\0") == (b"", b"")
    @test decodebin(b"m\0\0\0\0tail") == (b"", b"tail")
    @test decodebin(b"m\0\0\0\4data") == (b"data", b"")
    @test decodebin(b"m\0\0\0\4datatail") == (b"data", b"tail")
    @test_throws IncompleteData decodebin(b"m")
    @test_throws IncompleteData decodebin(b"m\0")
    @test_throws IncompleteData decodebin(b"m\0\0")
    @test_throws IncompleteData decodebin(b"m\0\0\0")
    @test_throws IncompleteData decodebin(b"m\0\0\0\1")
    @test decodebin(b"m\0\0\0\x1f\xc3\xa1rv\xc3\xadzt\xc5\xb1r\xc5\x91\x20t\xc3\xbck\xc3\xb6rf\xc3\xbar\xc3\xb3g\xc3\xa9p") == (b"árvíztűrő tükörfúrógép", b"")
end

# decode bitstring (BIT_BINARY_EXT)
testcase() do
    # XXX follow-up
    #@test_throws IncompleteData decode(b"\x83M")
    #@test_throws IncompleteData decode(b"\x83M\0")
    #@test_throws IncompleteData decode(b"\x83M\0\0")
    #@test_throws IncompleteData decode(b"\x83M\0\0\0")
    #@test_throws IncompleteData decode(b"\x83M\0\0\0\0")
    #@test_throws IncompleteData decode(b"\x83M\0\0\0\1\0")
    #@test_throws IncompleteData decode(b"\x83M\0\0\0\1\0")
    #@test decode(b"\x83M\0\0\0\0\0") == (b"", b"")
    #@test decode(b"\x83M\0\0\0\0\3") == (b"", b"")
    #@test decode(b"\x83M\0\0\0\0\5tail") == (b"", b"tail")
    #@test decode(b"\x83M\0\0\0\1\1\x80") == (b"\1", b"")
    #@test decode(b"\x83M\0\0\0\1\1\x80tail") == (b"\1", b"tail")
end

# decode float (FLOAT_EXT)
testcase() do
    # XXX that is represented by its sprintf'ed string...
end

# decode new float (NEW_FLOAT_EXT)
testcase() do
    @test_throws IncompleteData decode(b"\x83F")
    @test_throws IncompleteData decode(b"\x83F\0")
    @test_throws IncompleteData decode(b"\x83F\0\0")
    @test_throws IncompleteData decode(b"\x83F\0\0\0")
    @test_throws IncompleteData decode(b"\x83F\0\0\0\0")
    @test_throws IncompleteData decode(b"\x83F\0\0\0\0\0")
    @test_throws IncompleteData decode(b"\x83F\0\0\0\0\0\0")
    @test_throws IncompleteData decode(b"\x83F\0\0\0\0\0\0\0")
    @test decode(b"\x83F\0\0\0\0\0\0\0\0") == (0.0, b"")
    @test decode(b"\x83F\0\0\0\0\0\0\0\0tail") == (0.0, b"tail")
    @test decode(b"\x83F?\xf8\0\0\0\0\0\0") == (1.5, b"")
    @test decode(b"\x83F?\xf8\0\0\0\0\0\0tail") == (1.5, b"tail")
    @test decodenewfloat(b"F\0\0\0\0\0\0\0\0") == (0.0, b"")
    @test decodenewfloat(b"F\0\0\0\0\0\0\0\0tail") == (0.0, b"tail")
    @test decodenewfloat(b"F?\xf8\0\0\0\0\0\0") == (1.5, b"")
    @test decodenewfloat(b"F?\xf8\0\0\0\0\0\0tail") == (1.5, b"tail")
    @test_throws IncompleteData decodenewfloat(b"F")
    @test_throws IncompleteData decodenewfloat(b"F\0")
    @test_throws IncompleteData decodenewfloat(b"F\0\0")
    @test_throws IncompleteData decodenewfloat(b"F\0\0\0")
    @test_throws IncompleteData decodenewfloat(b"F\0\0\0\0")
    @test_throws IncompleteData decodenewfloat(b"F\0\0\0\0\0")
    @test_throws IncompleteData decodenewfloat(b"F\0\0\0\0\0\0")
    @test_throws IncompleteData decodenewfloat(b"F\0\0\0\0\0\0\0")
end

# decode small big integer (SMALL_BIG_EXT)
testcase() do
    @test_throws IncompleteData decode(b"\x83n")
    @test_throws IncompleteData decode(b"\x83n\0")
    @test_throws IncompleteData decode(b"\x83n\1\0")
    @test_throws IncompleteData decode(b"\x83n\2\0\0")
    @test decode(b"\x83n\1\0\0") == (0, b"")
    @test decode(b"\x83n\1\0\0tail") == (0, b"tail")

    @test decodesmallbigint(b"n\1\0\1") == (1, b"")
    @test decodesmallbigint(b"n\1\1\1") == (-1, b"")
    @test decodesmallbigint(b"n\2\0\1\2") == (513, b"")
    @test decodesmallbigint(smallbigintmax) == (256^255-1, b"")
end

# decode large big integer (LARGE_BIG_EXT)
testcase() do
    @test_throws IncompleteData decode(b"\x83o")
    @test_throws IncompleteData decode(b"\x83o\0")
    @test_throws IncompleteData decode(b"\x83o\0\0")
    @test_throws IncompleteData decode(b"\x83o\0\0\0")
    @test_throws IncompleteData decode(b"\x83o\0\0\0\0")
    @test_throws IncompleteData decode(b"\x83o\0\0\0\1\0")
    @test_throws IncompleteData decode(b"\x83o\0\0\0\2\0\0")
    @test decode(b"\x83o\0\0\0\1\0\0") == (0, b"")
    @test decode(b"\x83o\0\0\0\1\0\0tail") == (0, b"tail")

    @test decodelargebigint(b"o\0\0\0\1\0\1") == (1, b"")
    @test decodelargebigint(b"o\0\0\0\1\1\1") == (-1, b"")
    @test decodelargebigint(b"o\0\0\0\2\1\1\2") == (-513, b"")
    @test decodelargebigint(largebigintmin) == (256^255, b"")
end

# decode map (MAP_EXT)

# decode compressed term
testcase() do
    @test_throws IncompleteData decode(b"\x83P")
    @test_throws IncompleteData decode(b"\x83P\0")
    @test_throws IncompleteData decode(b"\x83P\0\0")
    @test_throws IncompleteData decode(b"\x83P\0\0\0")
    @test_throws IncompleteData decode(b"\x83P\0\0\0\0")
    @test_throws InvalidCompressedTag decode(badsizecompdata)
    @test decode(compdata1) == ([100,100,100,100,100,100,100,100,100,100,
                                 100,100,100,100, 100,100,100,100,100,100],
                                Uint8[])
    # XXX the following test fails because the Zlib library for Julia doesn't
    # provide a flush-like mechanism that the Python zlib library does
    # @test decode(compdata2) == ([100,100,100,100,100,100,100,100,100,100,
    #                              100,100,100,100, 100,100,100,100,100,100],
    #                             b"tail")
end


text:
@testset "Test cov_to_cor" begin

    cov1 = [2.74 1.15 0.86 1.31;
            1.15 1.67 0.82 1.32;
            0.86 0.82 2.12 1.05;
            1.31 1.32 1.05 2.93]

    cor1, unc1 = EFTfitter.cov_to_cor(cov1)

    @test cor1 ≈[1.0 0.5376062677719541 0.35682497116489126 0.46234078289360164;
                0.5376062677719541 1.0 0.4358004620316949 0.5967355144471027;
                0.35682497116489126 0.4358004620316949 1.0 0.4212962160310593;
                0.46234078289360164 0.5967355144471027 0.4212962160310593 1.0]

    @test unc1 ≈ [1.6552945357246849, 1.2922847983320085, 1.4560219778561037, 1.711724276862369]

    # test posdef warning
    cov2 = [2.74 -1.15 0.86 1.31;
            -1.15 1.67 0.82 1.32;
            0.86 0.82 2.12 1.05;
            1.31 1.32 1.05 2.93]
    @test_logs (:warn, "The covariance matrix [2.74 -1.15 0.86 1.31; -1.15 1.67 0.82 1.32; 0.86 0.82 2.12 1.05; 1.31 1.32 1.05 2.93] is not positive definite!") cor2, unc2 = EFTfitter.cov_to_cor(cov2)

end

@testset "Test all_observables_equal" begin
        estimator(params) = params.τ
        estimator2(params) = 2*params.τ

        parameters = BAT.NamedTupleDist(τ = 8..14, )

        covariance = [2.74 1.15 0.86 1.31;
                      1.15 1.67 0.82 1.32;
                      0.86 0.82 2.12 1.05;
                      1.31 1.32 1.05 2.93]

        corr, unc = EFTfitter.cov_to_cor(covariance)

        # all observables the same
        measurements = (
            τ1 = Measurement(estimator,  9.5, uncertainties = (stat=unc[1],) ),
            τ2 = Measurement(estimator, 11.9, uncertainties = (stat=unc[2],) ),
            τ3 = Measurement(estimator, 11.1, uncertainties = (stat=unc[3],) ),
            τ4 = Measurement(estimator,  8.9, uncertainties = (stat=unc[4],) ),
        )

        # not all observables the same
        measurements2 = (
            τ1 = Measurement(estimator,  9.5, uncertainties = (stat=unc[1],) ),
            τ2 = Measurement(estimator, 11.9, uncertainties = (stat=unc[2],) ),
            τ3 = Measurement(estimator2, 11.1, uncertainties = (stat=unc[3],) ),
            τ4 = Measurement(estimator,  8.9, uncertainties = (stat=unc[4],) ),
        )

        # all observables the same
        measurements3 = (
            τ1 = Measurement(estimator,  9.5, uncertainties = (stat=unc[1],) ),
            τ2 = Measurement(estimator, 11.9, uncertainties = (stat=unc[2],) ),
            τ3 = Measurement(estimator2, 11.1, uncertainties = (stat=unc[3],), active=false),
            τ4 = Measurement(estimator,  8.9, uncertainties = (stat=unc[4],) ),
        )

        correlations = (stat = Correlation(corr), )
        correlations2 = (stat = NoCorrelation(), )
        correlations3 = (stat = NoCorrelation(), )

        m = EFTfitterModel(parameters, measurements, correlations)
        m2 = EFTfitterModel(parameters, measurements2, correlations2)
        m3 = EFTfitterModel(parameters, measurements3, correlations3)

        @test  EFTfitter.all_observables_equal(m) == true
        @test  EFTfitter.all_observables_equal(m2) == false
        @test  EFTfitter.all_observables_equal(m3) == true
end

@testset "Test BLUE" begin
    estimator(params) = params.τ

    parameters = BAT.NamedTupleDist(τ = 8..14, )

    covariance = [2.74 1.15 0.86 1.31;
                  1.15 1.67 0.82 1.32;
                  0.86 0.82 2.12 1.05;
                  1.31 1.32 1.05 2.93]

    corr, unc = EFTfitter.cov_to_cor(covariance)

    measurements = (
        τ1 = Measurement(estimator,  9.5, uncertainties = (stat=unc[1],) ),
        τ2 = Measurement(estimator, 11.9, uncertainties = (stat=unc[2],) ),
        τ3 = Measurement(estimator, 11.1, uncertainties = (stat=unc[3],) ),
        τ4 = Measurement(estimator,  8.9, uncertainties = (stat=unc[4],) ),
    )

    correlations = (stat = Correlation(corr), )

    m = EFTfitterModel(parameters, measurements, correlations)

    blue = EFTfitter.BLUE(m)
    @test blue.value ≈ 11.1598305 rtol=0.005
    @test blue.unc ≈ 1.134037 rtol=0.005
    @test blue.weights ≈ [0.1450748, 0.46957738, 0.34729705, 0.0380508] rtol=0.005
end


text:
module ASTModels
export
  replace_subexpr!, substitute!, FloatX, cast, traverse_noconflict_expr_functions, expr, @expr,
  Model, default_input_fields, default_output_fields, merge_models, TypeWithModel,
  astfunc, @astfunc, clone, TestTypeWithModel, traverse_noconflict_methods_sig, create_function_file,
  traverse_ASTleaves

using Continuables
using DataStructures
include("./util.jl")

## expr functionalities ---------------------------------------------------------------------------------------------------------------------------

# when bringing anything else than Expr into an Expr, we preserve the possibility to cast the respective object
cast(x) = x  # identity

expr() = Expr(:_)  # empty expression

function expr(e)
  if isa(e, Symbol)
    :($e::Any)
  elseif !isa(e, Expr)
    :(cast($e))
  else
    e
  end
end

# expr(es...) = expr.(es)  # single case is already taken over by expr(e)

"""
@expr a b c
will expand to
```
a = :(a::Any)
b = :(b::Any)
c = :(c::Any)
```
creating new basic, identifyable Expressions
"""
macro expr(ss...)
  es = map(ss) do s
    # expr should be this package's `expr`
    rhs = QuoteNode(expr(s)) # we want this to stay expression

    if isa(s, Symbol)
      name = s
    elseif s.head == :(::) && isa(s.args[1], Symbol)
      name = s.args[1]
    else
      error("Can only handle blank symbols or additional ::Type checks.")
    end

    :($name = $rhs)
  end
  esc(Expr(:block, es..., nothing))  # hygiene
end

"""
creates a function expr which mimics a given function however with new types
"""
function create_expr_function(s::Symbol, types, m::Module=Base)
  vars = [Symbol(t===Expr ? :e : :o, i) for (i, t) in enumerate(types)]

  if (string(s)[1] ∈ ['.', '@']  # either is dot syntax (deprecated I thought, but still occurs, method seems to not exist) or a macro
    || eval(m, :(method_exists($s, $types)))  # or method already exists in the given package (maybe better to just execute this within this scope? would it be the global scope?)
    || (string(s) != function_name(eval(m, s)))) # or maybe this is a reference to something else, i.e. the symbol is not the same as the string transformation
    return
  end

  vars_typed = [Expr(:(::), x, t) for (x, t) in zip(vars, types)]
  vars_casted = [t === Expr ? x : :(expr($x)) for (x, t) in zip(vars, types)]

  f_head = Expr(:call, :($s), vars_typed...)
  f_body = Expr(:call, :Expr, :(:call), QuoteNode(s), vars_casted...)
  Expr(:function, f_head, f_body)
end

"""
continuable over function expressions for Expr types
"""
traverse_noconflict_expr_functions(mod::Module=Base, exclude=[]) = cont -> begin
  exclude = Set(exclude)
  for s in names(mod)
    if s in exclude
      continue
    end
      f = eval(mod, s)
      if isa(f, Function)
        traverse_noconflict_methods_sig(f) do types
          expr = create_expr_function(s, types, mod)
          if expr !== nothing
            cont(expr)
          end
        end
      end
  end
end

traverse_noconflict_expr_functions(cont, mod::Module=Base, exclude=[]) = traverse_noconflict_expr_functions(mod, exclude)(cont)

function create_function_file(path::String, mod::Module=Base)
  open(path, "w") do f
    write(f, "importall " * string(mod) * "\n")
    traverse_noconflict_expr_functions(Base) do expr
      write(f, string(expr) * "\n")
    end
  end
end
# csuppress(eval, Union{StackOverflowError, MethodError}))
# [:show, :haskey, :push!, :display, :setindex!]

# run for base (now everything seems to work fine, no exceptions anylonger):
# importall Base
# traverse_noconflict_expr_functions(Base)(eval)

# build using create_function_file()
include("./base_expr.jl")  # much faster than computing the functions itself

## General helpers for altering AST ---------------------------------------------------------------------------------------------------------

"""
Replace within an Expression on the toppest level only
"""
function replace_exprargs!(expr::Expr, old, new)
  for (i, subexpr) in enumerate(expr.args)
    if subexpr == old
      expr.args[i] = new
    end
  end
end

"""
Replace recursively within an Expression
"""
function replace_subexpr!(expr::Expr, old, new)
  replace_exprargs!(expr, old, new)
  for subexpr in expr.args
    replace_subexpr!(subexpr, old, new)
  end
end

replace_subexpr!(s, old, new) = ()  # there is no way to replace within anything else than expr

"""
Changes the original Expression to the substitute Expression.
Caution! this literally makes the old Expression to the substituted one. All references to ``ori`` will point to the contents of ``sub`` afterwards.
"""
function substitute!(ori::Expr, sub::Expr)
  old = Expr(ori.head, ori.args...)
  ori.head = sub.head
  ori.args = sub.args
  # ori.typ = sub.typ  # typ is not supported and might be deleted in future julia versions
  old
end

## models ---------------------------------------------------------------------------------------------------------------------------

typealias Model VectorAssociative{String, Expr}
# typealias Model Dict{String, Vector{Expr}}

# input field convention:
# outputs and anything starting with out is an outputfield (this is due to there normally being a lot of input things and only a few output things)
function merge_models(m::Model, others::Model...)
  d = merge_concat(m, others...)

  inputlike, outputlike = Vector{Expr}(), Vector{Expr}()
  for k ∈ keys(d)
    if startswith(k, "out")
      append!(outputlike, d[k])
    else
      append!(inputlike, d[k])
    end
  end
  inputlike, outputlike = Set(inputlike), Set(outputlike)

  for k ∈ keys(d)
    if startswith(k, "out")
      d[k] = [x for x ∈ d[k] if x ∉ inputlike]
    else
      d[k] = [x for x ∈ d[k] if x ∉ outputlike]
    end
    isempty(d[k]) && delete!(d, k)
  end
  d
end

# bringing merge_models to Struct with internal dictionary

abstract TypeWithModel
type TestTypeWithModel <: TypeWithModel
  dict::Model
  #... add you favourite specific attributes and your favourite constructor
end

import Base.convert
convert(::Type{Model}, x::TypeWithModel) = x.dict
merge_models(models::TypeWithModel...) = merge_models([convert(Model, m) for m in models]...)  # we somehow have to manually call convert.

function clone(model::Model)
  model2 = copy(model)

  # copy all parameters
  for k in keys(model2)
    if startswith(k, "param")
      model2[k] = copy.(model2[k])
    end
  end

  # collect everything which is not output
  oldinputs = collect(Base.flatten([model[k] for k in keys(model) if !startswith(k, "out")]))
  newinputs = collect(Base.flatten([model2[k] for k in keys(model) if !startswith(k, "out")]))

  # new output variables
  for k in keys(model2)
    if startswith(k, "out")
      model2[k] = clone(oldinputs, model2[k], newinputs)
    end
  end
  model2
end

## building function out of astmodels ---------------------------------------------------------------------------

# functions unfortunately don't work, as ``eval`` is going to execute the code within ASTModels, however we want it to be executed in the calling evironment
# still we need can use them for clone, where we only need a anonymous function

# """ build a function out of AST references """
function astfunc{T<:Union{Tuple, Array, Expr}}(func_name::Symbol, args, body::T, mod::Module=Main, types=[extract_type_info(a) for a in args])  # named version
  # we need to have Expressions
  args = [args...]  # for .() broadcasting
  if T <: Tuple
    body = Expr(:tuple, body...)
  elseif T <: Array
    body = Expr(:vect, body...)
  end

  new_args = [Expr(:(::), Symbol("p", i), t) for (i, t) in enumerate(types)]
  old_symbols = substitute!.(args, new_args)  # Dict(zip(args, symbols))  ??
  f_def = Expr(:call, func_name, new_args...)
  eval(mod, Expr(:(=), f_def, body))  # evaluate function
  new_args = substitute!.(args, old_symbols)
  eval(mod, Symbol(func_name))
end

# astfunc{T<:Union{Tuple, Array, Expr}}(func_name::Symbol, args, body::T, mod::Module=Main, types=[extract_type_info(a) for a in args]) = astfunc(func_name, args, [extract_type_info(a) for a in args], body, mod)

function astfunc{T<:Union{Tuple, Array, Expr}}(args, body::T, mod::Module=Main, types=[extract_type_info(a) for a in args])  # default to interpret body as iterable
  # we need to have Expressions
  args = [args...]  # for .() broadcasting
  if T <: Tuple
    body = Expr(:tuple, body...)
  elseif T <: Array
    body = Expr(:vect, body...)
  end

  new_args = [Expr(:(::), Symbol("p", i), t) for (i, t) in enumerate(types)]
  old_symbols = substitute!.(args, new_args)
  f_def = Expr(:tuple, new_args...)
  f = eval(mod, Expr(:(->), f_def, body))  # evaluate function

  new_args = substitute!.(args, old_symbols)
  f
end

# astfunc{T<:Union{Tuple, Array, Expr}}(args, body::T, mod::Module=Main) = astfunc(args, , body, mod)

function _eval_args_body(args_plain, body, mod::Module=Main)
  args = [eval(mod, ASTModels.extract_symbol(a)) for a in args_plain]
  if body.head ∈ [:tuple, :vect]
    # body.args .= eval(mod, body.args)  # in Julia 0.6 I think this should work
    for (i, a) in enumerate(body.args)
      body.args[i] = eval(mod, a)  # eval all subexpressions
    end
  else
    body = eval(mod, body)
  end

  types = Vector(length(args))
  for (i, (a, b)) in enumerate(zip(args_plain, args))
    if isa(a, Expr) && a.head == :(::)  # prefer overwritten types
      types[i] = a.args[2]
      if ASTModels.extract_type_info(b) != :Any
        warn("Overwrote Type Information for " * string(b))
      end
    else
      types[i] = ASTModels.extract_type_info(b)
    end
  end

  args, body, types
end

""" macro to build a function out of AST references """
macro astfunc(e)  # CAUTION: the body always has to be a single expression
  if e.head == :(=)
    f = e.args[1]
    f_name = QuoteNode(f.args[1])  # we need f_name to be a symbol if inserted below, not a function
    args = f.args[2:end]
    body = QuoteNode(ASTModels.ast_minify(e.args[2]))

    expr = quote
      mod = current_module()
      args, body, types = ASTModels._eval_args_body($args, $body, mod)
      ASTModels.astfunc($f_name, args, body, mod, types)
    end

  elseif e.head == :(->)
    if e.args[1].head == :tuple
      args = e.args[1].args
    else
      args = (e.args[1],)  # single symbol -> tuple of single symbol
    end
    body = QuoteNode(ASTModels.ast_minify(e.args[2]))

    expr = quote
      mod = current_module()
      args, body, types = ASTModels._eval_args_body($args, $body, mod)
      ASTModels.astfunc(args, body, mod, types)
    end

  else
    error("Use this macro for simple function syntax only.")
  end

  expr  # we don't want to escape things as all variables are indeed for internal purposes
end

"""
clone a given graph with optional new inputs (otherwise old inputs are preserved)
"""
function clone(inputs, outputs; new_inputs=[])
  f = astfunc(inputs, outputs)
  if length(new_inputs) == 0
    new_inputs = inputs  # use old inputs
  end
  f(new_inputs...)
end

clone(outputs) = clone(collect(traverse_ASTleaves(outputs)), outputs)  # new_inputs make no sense if you don't know yet what the inputs are

traverse_ASTleaves(l::Union{Array,Tuple}) = cont -> begin
  for a in l
    if isa(a, Expr)
      traverse_ASTleaves(a)(cont)
    end
  end
end

traverse_ASTleaves(expr::Expr) = cont -> begin
  if expr.head==:(::) && isa(expr.args[1], Symbol)
    # for simplicity we regard every symbol complemented with type check as a leave
    cont(expr)

  elseif expr.head ∈ [:(=), :function]
    # skip over assignment part, as this naturally includes symbols
    traverse_ASTleaves(expr.args[2:end])(cont)

  else
    # travers everything
    traverse_ASTleaves(expr.args)(cont)
  end
end
traverse_ASTleaves(cont, expr::Expr) = traverse_ASTleaves(expr)(cont)

## working with ASTModels --------------------------------------------------------------------------------------

# here we see one drawback of ASTModels: the explicit dependency structure captured normally in order of execution easily is lost here
# we have to reconstruct the dependencies by traversing the AST
# as this is only needed at construction time and not runtime, it is still okay to do so

# this one seems very efficient and short written
function sort_by_dependency(list_of_expr::Vector{Expr})
  l = pointer_from_objref.(list_of_expr)
  deps_in!(list_of_expr, Set(l), Set(l))
end
sort_by_dependency(cont, list_of_expr::Vector{Expr}) = sort_by_dependency(list_of_expr)(cont)

function deps_in!{T<:Any}(cont, es::Vector{T}, pointers!, pointers)
  for e ∈ es
    !isa(e, Expr) && continue

    if pointer_from_objref(e) ∈ pointers!  # go into everything we haven't seen yet
      setdiff!(pointers!, [pointer_from_objref(e)])
      deps_in!(cont, e.args[2:end], pointers!, pointers)  # the first element of e.args is the name of the function or co, i.e. not needed
      cont(e)
    elseif pointer_from_objref(e) ∈ pointers  # don't traverse things again which we already traversed
      continue   # the first element of e.args is the name of the function or co, i.e. not needed
    else
      deps_in!(cont, e.args[2:end], pointers!, pointers)
    end
  end
end
deps_in!{T<:Any}(es::Vector{T}, pointers!, pointers) = cont -> deps_in!(cont, es, pointers!, pointers)

end # module


text:
function specificity(tar, nontar, p)
    sum(nontar .< p) / length(nontar)
end

function sensitivity(tar, nontar, p)
    sum(tar .>= p) / length(tar)
end


text:
using Tricks
using InteractiveUtils
using Test

struct Iterable end; struct NonIterable end;
iterableness_static(::Type{T}) where T = static_hasmethod(iterate, Tuple{T}) ? Iterable() : NonIterable()

struct Foo end

# if code has no calls then it must be fully static
has_no_calls(ir) = all(stmt->!Meta.isexpr(stmt, :call), ir)

VERSION >= v"1.3" && @testset "static_hasmethod" begin
    @testset "positive: $(typeof(data))" for data in (
        "abc", [1,2,3], (2,3), ones(4,10,2), 'a',  1:100
    )
        T = typeof(data)
        @test iterableness_static(T) === Iterable()
        code_typed = (@code_typed iterableness_static(T))

        @test code_typed[2] === Iterable  # return type
        @test has_no_calls(code_typed[1].code)
    end

    @testset "negative: $(typeof(data))" for data in (
        :a, rand, Int
    )
        T = typeof(data)
        @test iterableness_static(T) === NonIterable()
        code_typed = (@code_typed iterableness_static(T))

        @test code_typed[2] === NonIterable  # return type
        @test has_no_calls(code_typed[1].code)
    end

    @testset "add method" begin
        @test iterableness_static(Foo) === NonIterable()

        Base.iterate(::Foo) = ("Foo", nothing);
        Base.iterate(::Foo, ::Nothing) = nothing;
        Base.length(::Foo) = 1;
        @test collect(Foo()) == ["Foo"]

        @test iterableness_static(Foo) === Iterable()
    end

    @testset "delete method" begin
        @test iterableness_static(Foo) === Iterable()
        meth = first(methods(iterate, Tuple{Foo}))
        Base.delete_method(meth)
        @test_throws MethodError collect(Foo())

        @test iterableness_static(Foo) === NonIterable()
    end

    @testset "abstract type args" begin
        # https://github.com/oxinabox/Tricks.jl/issues/14

        goo(x::Integer) = 1
        @assert !hasmethod(goo, Tuple{Real})  # the behaviour we want to match
        @test !static_hasmethod(goo, Tuple{Real})

        goo(x::Number) = 2
        @assert hasmethod(goo, Tuple{Real})   # Now it _is_ covered.
        @test static_hasmethod(goo, Tuple{Real})   # Now it _is_ covered.
    end
end

@testset "compat_hasmethod" begin
    @static if VERSION < v"1.3"
        @test compat_hasmethod == hasmethod
    else
        @test compat_hasmethod == static_hasmethod
    end
end

module Bar
h(::Int) = 1

end
using .Bar

VERSION >= v"1.3" && @testset "static_methods" begin
    # behavour
    f(x) = x + 1
    @test (length ∘ collect ∘ static_methods)(f) == 1
    f(::Int) = 1
    @test (length ∘ collect ∘ static_methods)(f) == 2

    g(::Int) = 1
    @test (length ∘ collect ∘ static_methods)(g) == 1
    g(x) = x+1
    @test (length ∘ collect ∘ static_methods)(g) == 2

    @test (length ∘ collect ∘ static_methods)(Bar.h) == 1
    Bar.h(x) = x
    @test (length ∘ collect ∘ static_methods)(Bar.h) == 2

    # Code Generation
    code_typed = (@code_typed static_methods(f))
    @test code_typed[2] === Base.MethodList  # return type
    @test has_no_calls(code_typed[1].code)

    @testset "delete method" begin
        i(::Int) = 1
        @test (length ∘ collect ∘ static_methods)(i) == 1
        i(x) = x+1
        @test (length ∘ collect ∘ static_methods)(i) == 2

        Base.delete_method((first ∘ methods)(i))
        @test (length ∘ collect ∘ static_methods)(i) == 1

        Base.delete_method((first ∘ methods)(i))
        @test (length ∘ collect ∘ static_methods)(i) == 0

        i(x) = x+1
        @test (length ∘ collect ∘ static_methods)(i) == 1
    end
end
@testset "type signatures" begin
    j(x) = x
    j(x, y, z...) = x + y + length(y)
    @test (length ∘ collect ∘ static_methods)(j) == 2
    @test static_method_count(j) == 2

    @test (length ∘ collect ∘ static_methods)(j, Tuple{Int}) == 1
    @test static_method_count(j, Tuple{Int,Int}) == 1
    @test (length ∘ collect ∘ static_methods)(j, Tuple{Int,Int}) == 1
    @test static_method_count(j, Tuple{Int,Int}) == 1
    @test (length ∘ collect ∘ static_methods)(j, Tuple{Int,Int,Int}) == 1
    @test static_method_count(j, Tuple{Int,Int,Int}) == 1

    # Varargs
    @test (length ∘ collect ∘ static_methods)(j, Tuple{Int, Int, Vararg{Int}}) == 1
    @test static_method_count(j, Tuple{Int, Int, Vararg{Int}}) == 1

    @test (length ∘ collect ∘ static_methods)(j, Tuple{Vararg{Int}}) == 2
    @test static_method_count(j, Tuple{Vararg{Int}}) == 2

    @test (length ∘ collect ∘ static_methods)(j, Tuple{Any, Vararg{Int}}) == 2
    @test static_method_count(j, Tuple{Any, Vararg{Int}}) == 2

    j(x, y::Int, z...) = x + y + length(y)

    @test (length ∘ collect ∘ static_methods)(j, Tuple{Int, Int, Vararg{Int}}) == 1
    @test static_method_count(j, Tuple{Int, Int, Vararg{Int}}) == 1
    @test (length ∘ collect ∘ static_methods)(j, Tuple{Int, Any, Vararg{Int}}) == 2
    @test static_method_count(j, Tuple{Int, Any, Vararg{Int}}) == 2

    @test (length ∘ collect ∘ static_methods)(j, Tuple{Vararg{Any}}) == 3
    @test static_method_count(j, Tuple{Vararg{Any}}) == 3

    @test (length ∘ collect ∘ static_methods)(j, Tuple{Any, Vararg{Any}}) == 3
    @test static_method_count(j, Tuple{Any, Vararg{Any}}) == 3

    @testset "delete method" begin

        while length(collect(methods(j))) > 0
            Base.delete_method((first ∘ methods)(j))

            T = Tuple{Any, Vararg{Any}}
            @test (length ∘ collect ∘ static_methods)(j, T) == length(collect(methods(j)))
            @test static_method_count(j, T) == length(collect(methods(j)))
        end
    end
end

# Redefine module Bar
module Bar
h(::Int) = 1
end
VERSION >= v"1.3" && @testset "static_method_count" begin
    # behavour
    f(x) = x + 1
    @test static_method_count(f) == 1
    f(::Int) = 1
    @test static_method_count(f) == 2

    g(::Int) = 1
    @test static_method_count(g) == 1
    g(x) = x+1
    @test static_method_count(g) == 2

    @test static_method_count(Bar.h) == 1
    Bar.h(x) = x
    @test static_method_count(Bar.h) == 2

    # Code Generation
    code_typed = (@code_typed static_method_count(f))
    @test code_typed[2] === Int  # return type
    @test has_no_calls(code_typed[1].code)

    @testset "delete method" begin
        i(::Int) = 1
        @test static_method_count(i) == 1
        i(x) = x+1
        @test static_method_count(i) == 2

        Base.delete_method((first ∘ methods)(i))
        @test static_method_count(i) == 1

        Base.delete_method((first ∘ methods)(i))
        @test static_method_count(i) == 0

        i(x) = x+1
        @test static_method_count(i) == 1
    end
end

VERSION >= v"1.3" && @testset "closures" begin
    make_add_n(n) = x->x+n
    func = make_add_n(2)

    # Add a 0-arg method:
    (::typeof(func))() = func(0)

    @assert func(1) == 3
    @assert func() == 2

    @testset "static_hasmethod" begin
        @assert hasmethod(func, Tuple{Int}) == true
        @test static_hasmethod(func, Tuple{Int}) == true

        @assert hasmethod(func, Tuple{}) == true
        @test static_hasmethod(func, Tuple{}) == true
    end

    @testset "static_methods" begin
        @test collect(methods(func, Tuple{Int})) ==
                collect(static_methods(func, Tuple{Int}))
        @test length(static_methods(func, Tuple{Int})) == 1

        @test collect(methods(func, Tuple{})) ==
                collect(static_methods(func, Tuple{}))
        @test length(static_methods(func, Tuple{})) == 1

        @test collect(methods(func, Tuple{Int,Int,Int})) ==
                collect(static_methods(func, Tuple{Int,Int,Int}))
        @test length(static_methods(func, Tuple{Int,Int,Int})) == 0
    end

    @testset "static_method_count" begin
        @test length(collect(methods(func, Tuple{Int}))) ==
                static_method_count(func, Tuple{Int}) == 1

        @test length(collect(methods(func, Tuple{}))) ==
                static_method_count(func, Tuple{}) == 1

        @test length(collect(methods(func, Tuple{Int,Int,Int}))) ==
                static_method_count(func, Tuple{Int,Int,Int}) == 0
    end
end

@testset "static_field____" begin
    function foo(data)
        names = static_fieldnames(typeof(data))
        map(name -> getproperty(data, name), names)
    end
    @test (@inferred foo(:a => 1)) == (:a, 1)

    bar(::Type{T}) where {T} = Val{static_fieldcount(T)}()
    @test @inferred(bar(Complex{Int})) == Val(2)

    baz(::Type{T}) where {T} = static_fieldtypes(T)
    @test @inferred(baz(Complex{Int})) == (Int, Int)
end


text:
function ulmfit_datadep_register()
    register(DataDep(
        "WikiText-103",
        """
        WikiText Long Term Dependency Language Modeling Dataset
        This is a language modelling dataset under Creative Commons Attribution-ShareAlike License. I contains over 100
        million tokens.
        """,
        "https://s3.amazonaws.com/research.metamind.io/wikitext/wikitext-103-v1.zip",
        post_fetch_method = function (file)
            unpack(file)
            dir = "wikitext-103"
            files = readdir(dir)
            mv.(joinpath.(dir, files), files)
            rm(dir)
            # some preprocessing
            corpus = read(open("wiki.train.tokens", "r"), String)
            corpus = split(lowercase(corpus), '\n')
            deleteat!(corpus, findall(x -> isequal(x, "")||isequal(x, " ")||(isequal(x[1:2], " =")&&isequal(x[prevind(x, lastindex(x), 1):end], "= ")), corpus))
            corpus = strip.(corpus)
            corpus = join(corpus, ' ')
            open("wiki.train.tokens", "w") do ftrain
                write(ftrain, corpus)
            end #do
        end # post_fetch_method
    ))

    #Pretrained ULMFiT language model weights
    register(DataDep(
        "Pretrained ULMFiT Language Model",
        """
        The pretrained Language Model weights trained on WikiText-103 corpus will be downloaded.
        These weights can be used to fine-tuning steps in ULMFiT classifier.
        This lanaguage model is originally trained by authors of ULMFiT paper Jeremy Howard and Sebastian Ruder.
        http://files.fast.ai/models/wt103/
        This contains the weights links for all original models by authors.
        """,
        "https://github.com/JuliaText/TextAnalysis.jl/releases/download/v0.6.0/ulmfit_lm_en.bson.tar.gz",
        "c71f580868fab7cc279cbceb0977ca26d252439a5cea604b2eebf6bf687127fa",
        post_fetch_method = function(fn)
            unpack(fn)
        end
    ))

    #Weights for ULMFiT binary sentiment classifier
    register(DataDep(
        "ULMFiT Sentiment Classifier",
        """
        Weights for the binary sentiment classifier, trained on IMDB movie review dataset,
        will be downloaded.
        """,
        "https://github.com/JuliaText/TextAnalysis.jl/releases/download/v0.6.0/ulmfit_sentiment_en.bson.tar.gz",# link
        post_fetch_method = function(fn)
            unpack(fn)
        end
    ))
end


text:
using JSON

@testset "Ext - JSON" begin
    rec = DefaultRecord("Logger.example", "info", 20, "blah")

    fmt = DictFormatter(JSON.json)
    result = Memento.format(fmt, rec)
    @test isa(JSON.parse(result), Dict)

    for key in [:date, :name, :level, :lookup, :stacktrace, :msg]
        @test occursin(string(key), result)
    end

    @test occursin("blah", result)

    aliases = Dict(
        :logger => :name,
        :level => :level,
        :timestamp => :date,
        :location => :lookup,
        :message => :msg,
        :process_id => :pid,
    )

    fmt2 = DictFormatter(aliases, JSON.json)
    result = Memento.format(fmt2, rec)

    @test isa(JSON.parse(result), Dict)

    for key in [:location, :message, :timestamp, :process_id]
        @test occursin(string(key), result)
        @test !occursin(string(aliases[key]), result)
    end

    @test occursin("blah", result)
end


text:
#=== wavefield initialization===#
#=== elastic 2d ===#
function initelwf(nDZ::Int64, nHX::Int64, ext::Int64, iflag::Int64)

        if iflag == 1 # free surface
                return elwf2d(zeros(Float64,(nDZ+ext),(nHX-1+2*ext)),
                              zeros(Float64,(nDZ-1+ext),(nHX+2*ext)),
                              zeros(Float64,(nDZ+ext),(nHX+2*ext)),
                              zeros(Float64,(nDZ+ext),(nHX+2*ext)),
                              zeros(Float64,(nDZ-1+ext),(nHX-1+2*ext)))
        elseif iflag == 2 # unlimited medium
                return elwf2d(zeros(Float64,(nDZ+2*ext),(nHX-1+2*ext)),
                              zeros(Float64,(nDZ-1+2*ext),(nHX+2*ext)),
                              zeros(Float64,(nDZ+2*ext),(nHX+2*ext)),
                              zeros(Float64,(nDZ+2*ext),(nHX+2*ext)),
                              zeros(Float64,(nDZ-1+2*ext),(nHX-1+2*ext)))
        else error("Please choose from 1 or 2. 1: free surface; 2: unlimited medium")
        end
end

#=== acoustic 2d ===#
function initacwf(nDZ::Int64, nHX::Int64, ext::Int64, iflag::Int64)

    # This function initialize the wavefield for forward modeling
    # It works for 2D cases

        if iflag == 1 # free surface
                return acwf2d(
                        zeros(Float64,(nDZ+ext),(nHX-1+2*ext)),
                        zeros(Float64,(nDZ-1+ext),(nHX+2*ext)),
                        zeros(Float64,(nDZ+ext),(nHX+2*ext)))
        elseif iflag == 2 # unlimited medium
                return acwf2d(
                        zeros(Float64,(nDZ+2*ext),(nHX-1+2*ext)),
                        zeros(Float64,(nDZ-1+2*ext),(nHX+2*ext)),
                        zeros(Float64,(nDZ+2*ext),(nHX+2*ext)))
        else error("Please choose from 1 or 2. 1: free surface; 2: unlimited medium")
        end
end

#=== elastic 3d ===#
function initelwf(nDZ::Int64, nHX::Int64, nHY::Int64, ext::Int64, iflag::Int64)

        if iflag == 1 # free surface
                return elwf3d(zeros(Float64,(nDZ+ext),(nHX-1+2*ext),(nHY+2*ext)),
                                 zeros(Float64,(nDZ+ext),(nHX+2*ext),(nHY-1+2*ext)),
                                 zeros(Float64,(nDZ-1+ext),(nHX+2*ext),(nHY+2*ext)),
                                 zeros(Float64,(nDZ+ext),(nHX+2*ext),(nHY+2*ext)),
                                 zeros(Float64,(nDZ+ext),(nHX+2*ext),(nHY+2*ext)),
                                 zeros(Float64,(nDZ+ext),(nHX+2*ext),(nHY+2*ext)),
                                 zeros(Float64,(nDZ-1+ext),(nHX-1+2*ext),(nHY+2*ext)),
                                 zeros(Float64,(nDZ-1+ext),(nHX+2*ext),(nHY-1+2*ext)),
                                 zeros(Float64,(nDZ+ext),(nHX-1+2*ext),(nHY-1+2*ext)))
        elseif iflag == 2 # unlimited medium
                return elwf3d(zeros(Float64,(nDZ+2*ext),(nHX-1+2*ext),(nHY+2*ext)),
                                 zeros(Float64,(nDZ+2*ext),(nHX+2*ext),(nHY-1+2*ext)),
                                 zeros(Float64,(nDZ-1+2*ext),(nHX+2*ext),(nHY+2*ext)),
                                 zeros(Float64,(nDZ+2*ext),(nHX+2*ext),(nHY+2*ext)),
                                 zeros(Float64,(nDZ+2*ext),(nHX+2*ext),(nHY+2*ext)),
                                 zeros(Float64,(nDZ+2*ext),(nHX+2*ext),(nHY+2*ext)),
                                 zeros(Float64,(nDZ-1+2*ext),(nHX-1+2*ext),(nHY+2*ext)),
                                 zeros(Float64,(nDZ-1+2*ext),(nHX+2*ext),(nHY-1+2*ext)),
                                 zeros(Float64,(nDZ+2*ext),(nHX-1+2*ext),(nHY-1+2*ext)))
        else error("Please choose from 1 or 2. 1: free surface; 2: unlimited medium")
        end
end

#=== acoustic 3d ===#
function initacwf(nDZ::Int64, nHX::Int64, nHY::Int64, ext::Int64, iflag::Int64)

    # This function initialize the wavefield for forward modeling
    # It works for 2D cases

        if iflag == 1 # free surface
                return acwf3d(zeros(Float64,(nDZ+ext),(nHX-1+2*ext),(nHY+2*ext)),
                                 zeros(Float64,(nDZ+ext),(nHX+2*ext),(nHY-1+2*ext)),
                                 zeros(Float64,(nDZ-1+ext),(nHX+2*ext),(nHY+2*ext)),
                                 zeros(Float64,(nDZ+ext),(nHX+2*ext),(nHY+2*ext)))
        elseif iflag == 2 # unlimited medium

                return acwf3d(zeros(Float64,(nDZ+2*ext),(nHX-1+2*ext),(nHY+2*ext)),
                                 zeros(Float64,(nDZ+2*ext),(nHX+2*ext),(nHY-1+2*ext)),
                                 zeros(Float64,(nDZ-1+2*ext),(nHX+2*ext),(nHY+2*ext)),
                                 zeros(Float64,(nDZ+2*ext),(nHX+2*ext),(nHY+2*ext)))

        else error("Please choose from 1 or 2. 1: free surface; 2: unlimited medium")
        end
end

#=== size of wavefield initialization ===#
#=== elastic 2d ===#
function initnelwf(nDZ::Int64, nHX::Int64, ext::Int64, iflag::Int64)

  if iflag == 1 # free surface
    return nelwf2d([nDZ+ext,nHX+2*ext],
                          [nDZ-1+ext,nHX-1+2*ext],
                          [nDZ+ext,nHX-1+2*ext],
                          [nDZ-1+ext,nHX+2*ext],
                          [nDZ,nHX],
                          [nDZ-1,nHX-1],
                          [nDZ,nHX-1],
                          [nDZ-1,nHX])

  elseif iflag == 2 # unlimited surface
    return nelwf2d([nDZ+2*ext,nHX+2*ext],
                          [nDZ-1+2*ext,nHX-1+2*ext],
                          [nDZ+2*ext,nHX-1+2*ext],
                          [nDZ-1+2*ext,nHX+2*ext],
                          [nDZ,nHX],
                          [nDZ-1,nHX-1],
                          [nDZ,nHX-1],
                          [nDZ-1,nHX])
  else error("Please choose from 1 or 2. 1: free surface; 2: unlimited medium")
  end

end

#=== acoustic 2d ===#
function initnacwf(nDZ::Int64, nHX::Int64, ext::Int64, iflag::Int64)

  if iflag == 1 # free surface
    return nacwf2d([nDZ+ext,nHX+2*ext],
                          [nDZ+ext,nHX-1+2*ext],
                          [nDZ-1+ext,nHX+2*ext],
                          [nDZ,nHX],
                          [nDZ,nHX-1],
                          [nDZ-1,nHX])

  elseif iflag == 2 # unlimited surface
    return nacwf2d([nDZ+2*ext,nHX+2*ext],
                          [nDZ+2*ext,nHX-1+2*ext],
                          [nDZ-1+2*ext,nHX+2*ext],
                          [nDZ,nHX],
                          [nDZ,nHX-1],
                          [nDZ-1,nHX])
  else error("Please choose from 1 or 2. 1: free surface; 2: unlimited medium")
  end

end

#=== elastic 3d ===#
function initnelwf(nDZ::Int64, nHX::Int64, nHY::Int64, ext::Int64, iflag::Int64)

  if iflag == 1 # free surface
    return nelwf3d([nDZ-1+ext,nHX+2*ext,nHY+2*ext],
                   [nDZ+ext,nHX-1+2*ext,nHY+2*ext],
                   [nDZ+ext,nHX+2*ext,nHY-1+2*ext],
                   [nDZ+ext,nHX+2*ext,nHY+2*ext],
                   [nDZ-1+ext,nHX-1+2*ext,nHY+2*ext],
                   [nDZ+ext,nHX-1+2*ext,nHY-1+2*ext],
                   [nDZ-1+ext,nHX+2*ext,nHY-1+2*ext],
                   [nDZ-1,nHX,nHY],
                   [nDZ,nHX-1,nHY],
                   [nDZ,nHX,nHY-1],
                   [nDZ,nHX,nHY],
                   [nDZ-1,nHX-1,nHY],
                   [nDZ,nHX-1,nHY-1],
                   [nDZ-1,nHX,nHY-1])

  elseif iflag == 2 # unlimited surface
    return nelwf3d([nDZ-1+2*ext,nHX+2*ext,nHY+2*ext],
                   [nDZ+2*ext,nHX-1+2*ext,nHY+2*ext],
                   [nDZ+2*ext,nHX+2*ext,nHY-1+2*ext],
                   [nDZ+2*ext,nHX+2*ext,nHY+2*ext],
                   [nDZ-1+2*ext,nHX-1+2*ext,nHY+2*ext],
                   [nDZ+2*ext,nHX-1+2*ext,nHY-1+2*ext],
                   [nDZ-1+2*ext,nHX+2*ext,nHY-1+2*ext],
                   [nDZ-1,nHX,nHY],
                   [nDZ,nHX-1,nHY],
                   [nDZ,nHX,nHY-1],
                   [nDZ,nHX,nHY],
                   [nDZ-1,nHX-1,nHY],
                   [nDZ,nHX-1,nHY-1],
                   [nDZ-1,nHX,nHY-1])
  else error("Please choose from 1 or 2. 1: free surface; 2: unlimited medium")
  end

end

#=== acoustic 3d ===#
function initnacwf(nDZ::Int64, nHX::Int64, nHY::Int64, ext::Int64, iflag::Int64)

  if iflag == 1 # free surface
    return nacwf3d([nDZ-1+ext,nHX+2*ext,nHY+2*ext],
                   [nDZ+ext,nHX-1+2*ext,nHY+2*ext],
                   [nDZ+ext,nHX+2*ext,nHY-1+2*ext],
                   [nDZ+ext,nHX+2*ext,nHY+2*ext],
                   [nDZ-1,nHX,nHY],
                   [nDZ,nHX-1,nHY],
                   [nDZ,nHX,nHY-1],
                   [nDZ,nHX,nHY])

  elseif iflag == 2 # unlimited surface
    return nacwf3d([nDZ-1+2*ext,nHX+2*ext,nHY+2*ext],
                   [nDZ+2*ext,nHX-1+2*ext,nHY+2*ext],
                   [nDZ+2*ext,nHX+2*ext,nHY-1+2*ext],
                   [nDZ+2*ext,nHX+2*ext,nHY+2*ext],
                   [nDZ-1,nHX,nHY],
                   [nDZ,nHX-1,nHY],
                   [nDZ,nHX,nHY-1],
                   [nDZ,nHX,nHY])
  else error("Please choose from 1 or 2. 1: free surface; 2: unlimited medium")
  end

end