text:
using DataFrames
showln(x) = (show(x); println())
# TODO: needs more links to docs.

# A DataFrame is an in-memory database
df = DataFrame(A = [1, 2], B = [e, pi], C = ["xx", "xy"])
showln(df)
#> 2x3 DataFrame
#> |-------|---|---------|------|
#> | Row # | A | B       | C    |
#> | 1     | 1 | 2.71828 | "xx" |
#> | 2     | 2 | 3.14159 | "xy" |

# The columns of a DataFrame can be indexed using numbers or names
showln(df[1])
#> [1,2]
showln(df[:A])
#> [1,2]

showln(df[2])
#> [2.718281828459045,3.141592653589793]
showln(df[:B])
#> [2.718281828459045,3.141592653589793]

showln(df[3])
#> ASCIIString["xx","xy"]
showln(df[:C])
#> ASCIIString["xx","xy"]

# The rows of a DataFrame can be indexed only by using numbers
showln(df[1, :])
#> 1x3 DataFrame
#> |-------|---|---------|------|
#> | Row # | A | B       | C    |
#> | 1     | 1 | 2.71828 | "xx" |
showln(df[1:2, :])
#> 2x3 DataFrame
#> |-------|---|---------|------|
#> | Row # | A | B       | C    |
#> | 1     | 1 | 2.71828 | "xx" |
#> | 2     | 2 | 3.14159 | "xy" |

# importing data into DataFrames
# ------------------------------

# DataFrames can be loaded from CSV files using readtable()
iris = readtable("iris.csv")

# the iris dataset (and plenty of others) is also available from
using RDatasets
iris = dataset("datasets","iris")

# you can directly import your own R .rda dataframe with
# mydf = DataFrame(read_rda("path/to/your/df.rda")["name_of_df"]), e.g.
diamonds = DataFrame(read_rda(joinpath(Pkg.dir("RDatasets"),"data","ggplot2","diamonds.rda"))["diamonds"])

# showing DataFrames
# ------------------

# Check the names and element types of the columns of our new DataFrame
showln(names(iris))
#> [:SepalLength,:SepalWidth,:PetalLength,:PetalWidth,:Species]
showln(eltypes(iris))
#> Type[Float64,Float64,Float64,Float64,UTF8String]

# Subset the DataFrame to only include rows for one species
showln(iris[iris[:Species] .== "setosa", :])
#> 50x5 DataFrame
#> |-------|-------------|------------|-------------|------------|----------|
#> | Row # | SepalLength | SepalWidth | PetalLength | PetalWidth | Species  |
#> | 1     | 5.1         | 3.5        | 1.4         | 0.2        | "setosa" |
#> | 2     | 4.9         | 3.0        | 1.4         | 0.2        | "setosa" |
#> | 3     | 4.7         | 3.2        | 1.3         | 0.2        | "setosa" |
#> | 4     | 4.6         | 3.1        | 1.5         | 0.2        | "setosa" |
#> | 5     | 5.0         | 3.6        | 1.4         | 0.2        | "setosa" |
#> | 6     | 5.4         | 3.9        | 1.7         | 0.4        | "setosa" |
#> | 7     | 4.6         | 3.4        | 1.4         | 0.3        | "setosa" |
#> | 8     | 5.0         | 3.4        | 1.5         | 0.2        | "setosa" |
#> | 9     | 4.4         | 2.9        | 1.4         | 0.2        | "setosa" |
#> ⋮
#> | 41    | 5.0         | 3.5        | 1.3         | 0.3        | "setosa" |
#> | 42    | 4.5         | 2.3        | 1.3         | 0.3        | "setosa" |
#> | 43    | 4.4         | 3.2        | 1.3         | 0.2        | "setosa" |
#> | 44    | 5.0         | 3.5        | 1.6         | 0.6        | "setosa" |
#> | 45    | 5.1         | 3.8        | 1.9         | 0.4        | "setosa" |
#> | 46    | 4.8         | 3.0        | 1.4         | 0.3        | "setosa" |
#> | 47    | 5.1         | 3.8        | 1.6         | 0.2        | "setosa" |
#> | 48    | 4.6         | 3.2        | 1.4         | 0.2        | "setosa" |
#> | 49    | 5.3         | 3.7        | 1.5         | 0.2        | "setosa" |
#> | 50    | 5.0         | 3.3        | 1.4         | 0.2        | "setosa" |

# Count the number of rows for each species
showln(by(iris, :Species, df -> size(df, 1)))
#> 3x2 DataFrame
#> |-------|--------------|----|
#> | Row # | Species      | x1 |
#> | 1     | "setosa"     | 50 |
#> | 2     | "versicolor" | 50 |
#> | 3     | "virginica"  | 50 |

# Discretize entire columns at a time
iris[:SepalLength] = iround(iris[:SepalLength])
iris[:SepalWidth] = iround(iris[:SepalWidth])

# Tabulate data according to discretized columns to see "clusters"
tabulated = by(
    iris,
    [:Species, :SepalLength, :SepalWidth],
    df -> size(df, 1)
)
showln(tabulated)
#> 17x4 DataFrame
#> |-------|--------------|-------------|------------|----|
#> | Row # | Species      | SepalLength | SepalWidth | x1 |
#> | 1     | "setosa"     | 4           | 3          | 4  |
#> | 2     | "setosa"     | 5           | 2          | 1  |
#> | 3     | "setosa"     | 5           | 3          | 23 |
#> | 4     | "setosa"     | 5           | 4          | 17 |
#> | 5     | "setosa"     | 6           | 4          | 5  |
#> | 6     | "versicolor" | 5           | 2          | 3  |
#> | 7     | "versicolor" | 5           | 3          | 3  |
#> | 8     | "versicolor" | 6           | 2          | 6  |
#> | 9     | "versicolor" | 6           | 3          | 29 |
#> | 10    | "versicolor" | 7           | 3          | 9  |
#> | 11    | "virginica"  | 5           | 3          | 1  |
#> | 12    | "virginica"  | 6           | 2          | 1  |
#> | 13    | "virginica"  | 6           | 3          | 22 |
#> | 14    | "virginica"  | 7           | 3          | 19 |
#> | 15    | "virginica"  | 7           | 4          | 1  |
#> | 16    | "virginica"  | 8           | 3          | 4  |
#> | 17    | "virginica"  | 8           | 4          | 2  |

# you can setup a grouped dataframe like this
gdf = groupby(iris,[:Species, :SepalLength, :SepalWidth])

# and then iterate over it
for idf in gdf
	println(size(idf,1))
end

# Adding/Removing columns
# -----------------------

# insert!(df::DataFrame,index::Int64,item::AbstractArray{T,1},name::Symbol)
# insert random numbers at col 5:
insert!(iris, 5, rand(nrow(iris)), :randCol)

# remove it
delete!(iris, :randCol)

# Same for rows
# push!(df::DataFrame,iterable)
# push an iterable as last row onto a dataframe:
push!(iris, [4, 2.01, -1.5, 100.1, "alienSpecies"])

# show last few rows
tail(iris)

# delete row with "alienSpecies"
deleterows!(iris, findin(iris[:Species], ["alienSpecies"]))
tail(iris)

# Using DataFramesMeta.jl
# -----------------------

using DataFramesMeta

# subset with DataFramesMeta
# @ix gets rows/colums
@ix(iris, :SepalWidth .> 4.0, [:SepalWidth, :Species])

# @where gets rows
 @where(iris, (:PetalWidth .> 1) & (:SepalWidth .< 3.0))

# create a new column with SepalWidth in 5 categories
# very similar to R cut() and R transform()
iris = transform( iris, widthBin = cut(iris[:SepalWidth],5))

# DataFramesMeta also has a convenient macro syntax,
# no need to re-specify the dataframe in the epxression!
iris = @transform( iris, widthBin = cut(:SepalWidth,5))

ext: jl

max_stars_repo_name: floswald/JuliaByExample

max_issues_repo_name: floswald/JuliaByExample

max_forks_repo_name: floswald/JuliaByExample


text:
#!/usr/bin/env julia

# LD_LIBRARY_PATH=/usr/lib/llvm-3.7/lib

using Clang.wrap_c

ctx = wrap_c.init(
        headers=["/usr/include/postgresql/libpq-fe.h"],
        common_file="common.jl",
        output_file="output.jl",
        header_library="libpq",
        header_wrapped=(header, cursorname) -> startswith(cursorname, "/usr/include/postgresql/"))
wrap_c.run(ctx)

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
using LightGraphs, Compose, GraphPlot, Cairo, Fontconfig

include("./parsers.jl")
####################################################################################
####################  Draw graph with edges and allocations  #######################
####################################################################################

function drawGraph(a, L, Tf, edges, num_nodes)
    edgelabels = []
    graph = LightGraphs.DiGraph(num_nodes)
    for e in 1:size(edges,1)
        LightGraphs.add_edge!(graph, edges[e][1], edges[e][2])
        s = 0
        for f in 1:size(a,1)
            for t in 1:size(a,2)
                s += a[f,t] * L[Tf[f][t],e]
            end
        end
        push!(edgelabels, s)
    end
    nodelabel = collect(1:num_nodes)
    Compose.draw(PNG("./graph.png", 1000, 1000), gplot(graph, edgelabelc="white", EDGELABELSIZE = 15.0, NODELABELSIZE=20.0, nodelabel=nodelabel, edgelabeldisty=0.5, edgelabel=edgelabels, EDGELINEWIDTH=3.0, arrowlengthfrac=.04))
end

function drawGraph(topology; outdir="./")
    edges, capacity, probabilities, nodes = readTopology(topology)
    graph = LightGraphs.DiGraph(length(nodes))
    for e in 1:length(edges)
        LightGraphs.add_edge!(graph, edges[e][1], edges[e][2])
    end
    nodelabel = collect(1:length(nodes))
    Compose.draw(PNG("$(outdir)/graph.png", 1000, 1000), gplot(graph, edgelabelc="white", EDGELABELSIZE = 15.0, NODELABELSIZE=20.0, nodelabel=nodelabel, edgelabeldisty=0.5, EDGELINEWIDTH=3.0, arrowlengthfrac=.04))
end


text:
using CompScienceMeshes

h = 1/6
rect = meshrectangle(1.0, 1.0, h, 3)
hole = meshrectangle(1/3, 1/3, h, 3)
translate!(hole, point(1/3, 1/3, 0))

pred = overlap_gpredicate(hole)
all_faces = submesh(c->!pred(chart(rect,c)), rect)

all_edges = skeleton(all_faces, 1)
all_verts = skeleton(all_faces, 0)

bnd_edges = boundary(all_faces)
bnd_verts = skeleton(bnd_edges, 0)

not(f) = x -> !f(x...)

onbnd1 = overlap_tpredicate(bnd_edges)
onbnd0 = overlap_tpredicate(bnd_verts)

interior_edges = submesh(not(onbnd1), all_edges)
interior_verts = submesh(not(onbnd0), all_verts)

@show numcells(interior_edges)
@show numcells(interior_verts)

D0 = connectivity(interior_verts, interior_edges)
D1 = connectivity(interior_edges, all_faces)

nullity(A) = size(A,2) - rank(A')
genus = nullity(full(D1)) - rank(full(D0))


text:
function stop_order(order, id, config, closed_orders, bot_path, client)
    #
    ## Cancel order
    result = try client.cancel_order(
        symbol=order["symbol"],
        orderId=order["orderId"])
    catch
        nothing
    end

    #
    ## Cancel order did not work, try it next round again
    if result == nothing
        return nothing
    end

    #
    ## Order was partially executed
    if parse(Float64, result["executedQty"]) != 0.0
        id_str = split(id, "_")[end]
        if order["side"] == "SELL"
            closed_orders[id_str]["SELL-" * lpad(rand(1:1:100000), 6, "0")] = deepcopy(result);
        else
            closed_orders[id_str]["EXTRA-" * lpad(rand(100:1:100000), 6, "0")] = deepcopy(result);
        end
        open(joinpath(bot_path, "closed_orders.json"), "w") do f
            JSON.print(f, closed_orders, 4)
        end
    end

    #
    ## Create new temp order
    new = Dict()
    new["price"] = order["price"]
    new["origQty"] = string(parse(Float64, order["origQty"]) - parse(Float64, order["executedQty"]))
    new["executedQty"] = order["executedQty"]
    new["symbol"] = order["symbol"]
    new["status"] = "stopped"
    new["side"] = order["side"]
    new["type"] = order["type"]
    return new
end
export stop_order

function revive_order(order, id, config, closed_orders, bot_path, client)
    qty = floor(parse(Float64, order["origQty"]), digits=config["precision"])
    if order["side"] == "SELL"
        result = try client.order_limit_sell(
            symbol=order["symbol"],
            quantity=qty,
            price=order["price"]);
            catch
            nothing
        end
    else
        result = try client.order_limit_buy(
            symbol=order["symbol"],
            quantity=qty,
            price=order["price"]);
            catch
            nothing
        end
    end

    return result
end
export revive_order

function check_stopped_orders(k, v, config, open_orders, closed_orders, bot_path, client, usr)
    @info("Stopping order")
    order = Pinguu.stop_order(v, k, config, closed_orders, bot_path, client)
    open(joinpath(bot_path, "trigger.json"), "w") do f
        JSON.print(f, Dict(), 4)
    end
    if order !== nothing
        #
        ## Change v with temp order
        open_orders[order["side"]][k] = order
        open(joinpath(bot_path, "open_orders.json"), "w") do f
            JSON.print(f, open_orders, 4)
        end
        return
    end
    @info("Cancelling order failed. Will try again soon.")
end
export check_stopped_orders

function check_orders_to_revive(k, v, config, open_orders, closed_orders, bot_path, client, usr)
    order = Pinguu.revive_order(v, k, config, closed_orders, bot_path, client)
    open(joinpath(bot_path, "trigger.json"), "w") do f
        JSON.print(f, Dict(), 4)
    end
    if order !== nothing
        #
        ## Change v with temp order
        open_orders[order["side"]][k] = order
        open(joinpath(bot_path, "open_orders.json"), "w") do f
            JSON.print(f, open_orders, 4)
        end
        v = deepcopy(order)
        return
    end
    @info("Reviving order failed. Will try again soon.")
end
export check_orders_to_revive


text:
using MacroTools

# Forwarding iteration

immutable SubIter{I,S}
  iter::I
  state::S
end

macro iter(ex)
  @capture(ex, x_::T_ -> it_) || error("Use @iter x::T -> y ...")
  @capture(it, $x.f_) &&
    return :(@forward $(esc(T)).$f Base.start, Base.next, Base.done)
  quote
    @inline function Base.start($x::$T)
      it = $it
      SubIter(it, Base.start(it))
    end
    @inline function Base.next(::$T, sub::SubIter)
      next, state = Base.next(sub.iter, sub.state)
      next, SubIter(sub.iter, state)
    end
    @inline function Base.done(::$T, sub::SubIter)
      Base.done(sub.iter, sub.state)
    end
  end |> esc
end


text:
"""
DirectImages.jl

A package of useful functions for interacting with FITS files.

# Available Functions

- `DirectImage(array)`
- `headers(img)`
- `axes(img)`
- `origin(img)`
- `size(img)`
- `copy(img)`
- `arraydata(img)` Gives underlying array

(All other Julia array functions as well)

- `readfits`
- `writefits`
- `readheader`  Only read headers. `readfits` gives both data and headers.

- `contrast`

"""
module DirectImages

include("Type.jl")
include("Metrics.jl")
include("IO.jl")
include("Display.jl")
include("Utilities.jl")
include("Stacking.jl")
include("precompile.jl")

end

max_stars_repo_name: sefffal/SpaceImages


text:
module MyPackage

using ForwardDiff

greet() = print("Hello World!")

include("my_func.jl")

export my_func, derivative_of_my_func

end # module


text:
@doc raw"""
    ellipk(m)

Computes Complete Elliptic Integral of 1st kind ``K(m)`` for parameter ``m`` given by

```math
\operatorname{ellipk}(m)
= K(m)
= \int_0^{ \frac{\pi}{2} } \frac{1}{\sqrt{1 - m \sin^2 \theta}} \, \mathrm{d}\theta
\quad \text{for} \quad m \in \left( -\infty, 1 \right] \, .
```

External links: [DLMF](https://dlmf.nist.gov/19.2.8), [Wikipedia](https://en.wikipedia.org/wiki/Elliptic_integral#Notational_variants).

See also: [`ellipe(m)`](@ref SpecialFunctions.ellipe).

# Arguments
- `m`: parameter ``m``, restricted to the domain ``(-\infty,1]``, is related to the elliptic modulus ``k`` by ``k^2=m`` and to the modular
    angle ``\alpha`` by ``k=\sin \alpha``.

# Implementation
Using piecewise approximation polynomial as given in
> 'Fast Computation of Complete Elliptic Integrals and Jacobian Elliptic Functions',
> Fukushima, Toshio. (2014). F09-FastEI. Celest Mech Dyn Astr,
> DOI 10.1007/s10569-009-9228-z,
> <https://pdfs.semanticscholar.org/8112/c1f56e833476b61fc54d41e194c962fbe647.pdf>

For ``m<0``, followed by
> Fukushima, Toshio. (2014).
> 'Precise, compact, and fast computation of complete elliptic integrals by piecewise
> minimax rational function approximation'.
> Journal of Computational and Applied Mathematics. 282.
> DOI 10.13140/2.1.1946.6245.,
> <https://www.researchgate.net/publication/267330394>
As suggested in this paper, the domain is restricted to ``(-\infty,1]``.
"""
function ellipk(m::Float64)
    flag_is_m_neg = false
    if m < 0.0
        x               = m / (m-1)         #dealing with negative args
        flag_is_m_neg   = true
    elseif m >= 0.0
        x               = m
        flag_is_m_neg   = false
    end

    if x == 0.0
        return pi/2

    elseif x == 1.0
        return Inf

    elseif x > 1.0
        throw(DomainError(m, "`m` must lie between -Inf and 1 ---- Domain: (-Inf,1.0]"))

    elseif 0.0 <= x < 0.1   #Table 2 from paper
        t = x-0.05
        t = @horner(t,
            1.591003453790792180 , 0.416000743991786912 , 0.245791514264103415 ,
            0.179481482914906162 , 0.144556057087555150 , 0.123200993312427711 ,
            0.108938811574293531 , 0.098853409871592910 , 0.091439629201749751 ,
            0.085842591595413900 , 0.081541118718303215)

    elseif 0.1 <= x < 0.2  #Table 3
        t = x-0.15
        t = @horner(t ,
            1.635256732264579992 , 0.471190626148732291 , 0.309728410831499587 ,
            0.252208311773135699 , 0.226725623219684650 , 0.215774446729585976 ,
            0.213108771877348910 , 0.216029124605188282 , 0.223255831633057896 ,
            0.234180501294209925 , 0.248557682972264071 , 0.266363809892617521)

    elseif 0.2 <= x < 0.3 #Table 4
        t = x-0.25
        t = @horner(t ,
            1.685750354812596043 , 0.541731848613280329 , 0.401524438390690257 ,
            0.369642473420889090 , 0.376060715354583645 , 0.405235887085125919 ,
            0.453294381753999079 , 0.520518947651184205 , 0.609426039204995055 ,
            0.724263522282908870 , 0.871013847709812357 , 1.057652872753547036)

    elseif 0.3 <= x < 0.4 #Table 5
        t = x-0.35
        t = @horner(t ,
            1.744350597225613243 , 0.634864275371935304 , 0.539842564164445538 ,
            0.571892705193787391 , 0.670295136265406100 , 0.832586590010977199 ,
            1.073857448247933265 , 1.422091460675497751 , 1.920387183402304829 ,
            2.632552548331654201 , 3.652109747319039160 , 5.115867135558865806 ,
            7.224080007363877411)

    elseif 0.4 <= x < 0.5 #Table 6
        t = x-0.45
        t = @horner(t,
            1.813883936816982644 , 0.763163245700557246 , 0.761928605321595831 ,
            0.951074653668427927 , 1.315180671703161215 , 1.928560693477410941 ,
            2.937509342531378755 , 4.594894405442878062 , 7.330071221881720772 ,
            11.87151259742530180 , 19.45851374822937738 , 32.20638657246426863 ,
            53.73749198700554656 , 90.27388602940998849)

    elseif 0.5 <= x < 0.6 #Table 7
        t = x-0.55
        t = @horner(t ,
            1.898924910271553526 , 0.950521794618244435 , 1.151077589959015808 ,
            1.750239106986300540 , 2.952676812636875180 , 5.285800396121450889 ,
            9.832485716659979747 , 18.78714868327559562 , 36.61468615273698145 ,
            72.45292395127771801 , 145.1079577347069102 , 293.4786396308497026 ,
            598.3851815055010179 , 1228.420013075863451 , 2536.529755382764488)
    elseif 0.6 <= x < 0.7 #Table 8
        t = x-0.65
        t = @horner(t ,
            2.007598398424376302 , 1.248457231212347337 , 1.926234657076479729 ,
            3.751289640087587680 , 8.119944554932045802 , 18.66572130873555361 ,
            44.60392484291437063 , 109.5092054309498377 , 274.2779548232413480 ,
            697.5598008606326163 , 1795.716014500247129 , 4668.381716790389910 ,
            12235.76246813664335 , 32290.17809718320818 , 85713.07608195964685 ,
            228672.1890493117096 , 612757.2711915852774)

    elseif 0.7 <= x < 0.8 #Table 9
        t = x-0.75
        t = @horner(t,
            2.156515647499643235 , 1.791805641849463243 , 3.826751287465713147 ,
            10.38672468363797208 , 31.40331405468070290 , 100.9237039498695416 ,
            337.3268282632272897 , 1158.707930567827917 , 4060.990742193632092 ,
            14454.00184034344795 , 52076.66107599404803 , 189493.6591462156887 ,
            695184.5762413896145 , 2567994.048255284686 , 9541921.966748386322 ,
            35634927.44218076174 , 133669298.4612040871 , 503352186.6866284541 ,
            1901975729.538660119 , 7208915015.330103756)

    elseif 0.8 <= x < 0.85 #Table 10
        t = x-0.825
        t = @horner(t ,
            2.318122621712510589 , 2.616920150291232841 , 7.897935075731355823 ,
            30.50239715446672327 , 131.4869365523528456 , 602.9847637356491617 ,
            2877.024617809972641 , 14110.51991915180325 , 70621.44088156540229 ,
            358977.2665825309926 , 1847238.263723971684 , 9600515.416049214109 ,
            50307677.08502366879 , 265444188.6527127967 , 1408862325.028702687 ,
            7515687935.373774627)

    elseif 0.85 <= x < 0.9 #Table 11
        t = x-0.875
        t = @horner(t,
            2.473596173751343912 , 3.727624244118099310 , 15.60739303554930496 ,
            84.12850842805887747 , 506.9818197040613935 , 3252.277058145123644 ,
            21713.24241957434256 , 149037.0451890932766 , 1043999.331089990839 ,
            7427974.817042038995 , 53503839.67558661151 , 389249886.9948708474 ,
            2855288351.100810619 , 21090077038.76684053 , 156699833947.7902014 ,
            1170222242422.439893 , 8777948323668.937971 , 66101242752484.95041 ,
            499488053713388.7989 , 37859743397240299.20)

    elseif x >= 0.9
        td  = 1-x
        td1 = td-0.05
        qd  = @horner(td,
            0.0, (1.0/16.0), (1.0/32.0), (21.0/1024.0), (31.0/2048.0), (6257.0/524288.0),
            (10293.0/1048576.0), (279025.0/33554432.0), (483127.0/67108864.0),
            (435506703.0/68719476736.0), (776957575.0/137438953472.0) ,
            (22417045555.0/4398046511104.0) , (40784671953.0/8796093022208.0) ,
            (9569130097211.0/2251799813685248.0) , (17652604545791.0/4503599627370496.0))

        kdm = @horner(td1 ,
            1.591003453790792180 , 0.416000743991786912 , 0.245791514264103415 ,
            0.179481482914906162 , 0.144556057087555150 , 0.123200993312427711 ,
            0.108938811574293531 , 0.098853409871592910 , 0.091439629201749751 ,
            0.085842591595413900 , 0.081541118718303215)
        km  = -Base.log(qd) * (kdm/pi)
        t   = km
    end

    if flag_is_m_neg
        ans = t / sqrt(1.0-m)
        return ans
    else
        return t
    end
end

@doc raw"""
    ellipe(m)

Computes Complete Elliptic Integral of 2nd kind ``E(m)`` for parameter ``m`` given by

```math
\operatorname{ellipe}(m)
= E(m)
= \int_0^{ \frac{\pi}{2} } \sqrt{1 - m \sin^2 \theta} \, \mathrm{d}\theta
\quad \text{for} \quad m \in \left( -\infty, 1 \right] \, .
```

External links: [DLMF](https://dlmf.nist.gov/19.2.8), [Wikipedia](https://en.wikipedia.org/wiki/Elliptic_integral#Complete_elliptic_integral_of_the_second_kind).

See also: [`ellipk(m)`](@ref SpecialFunctions.ellipk).

# Arguments
- `m`: parameter ``m``, restricted to the domain ``(-\infty,1]``, is related to the elliptic modulus ``k`` by ``k^2=m`` and to the modular
    angle ``\alpha`` by ``k=\sin \alpha``.

# Implementation
Using piecewise approximation polynomial as given in
> 'Fast Computation of Complete Elliptic Integrals and Jacobian Elliptic Functions',
> Fukushima, Toshio. (2014). F09-FastEI. Celest Mech Dyn Astr,
> DOI 10.1007/s10569-009-9228-z,
> <https://pdfs.semanticscholar.org/8112/c1f56e833476b61fc54d41e194c962fbe647.pdf>

For ``m<0``, followed by
> Fukushima, Toshio. (2014).
> 'Precise, compact, and fast computation of complete elliptic integrals by piecewise
> minimax rational function approximation'.
> Journal of Computational and Applied Mathematics. 282.
> DOI 10.13140/2.1.1946.6245.,
> <https://www.researchgate.net/publication/267330394>
As suggested in this paper, the domain is restricted to ``(-\infty,1]``.
"""
function ellipe(m::Float64)
    flag_is_m_neg = false
    if m < 0.0
        x               = m / (m-1)         #dealing with negative args
        flag_is_m_neg   = true
    elseif m >= 0.0
        x               = m
        flag_is_m_neg   = false
    end

    if x == 0.0
        return pi/2
    elseif x == 1.0
        return 1.0

    elseif x > 1.0
        throw(DomainError(m,"`m` must lie between -Inf and 1 ---- Domain : (-inf,1.0]"))

    elseif 0.0 <= x < 0.1   #Table 2 from paper
        t = x-0.05
        t = @horner(t  ,
            +1.550973351780472328 , -0.400301020103198524 , -0.078498619442941939 ,
            -0.034318853117591992 , -0.019718043317365499 , -0.013059507731993309 ,
            -0.009442372874146547 , -0.007246728512402157 , -0.005807424012956090 ,
            -0.004809187786009338)

    elseif 0.1 <= x < 0.2  #Table 3
        t = x-0.15
        t = @horner(t ,
            +1.510121832092819728 , -0.417116333905867549 , -0.090123820404774569 ,
            -0.043729944019084312 , -0.027965493064761785 , -0.020644781177568105 ,
            -0.016650786739707238 , -0.014261960828842520 , -0.012759847429264803 ,
            -0.011799303775587354 , -0.011197445703074968)

    elseif 0.2 <= x < 0.3 #Table 4
        t = x-0.25
        t = @horner(t ,
            +1.467462209339427155 , -0.436576290946337775 , -0.105155557666942554 ,
            -0.057371843593241730 , -0.041391627727340220 , -0.034527728505280841 ,
            -0.031495443512532783 , -0.030527000890325277 , -0.030916984019238900 ,
            -0.032371395314758122 , -0.034789960386404158)

    elseif 0.3 <= x < 0.4 #Table 5
        t = x-0.35
        t = @horner(t,
            +1.422691133490879171 , -0.459513519621048674 , -0.125250539822061878,
            -0.078138545094409477 , -0.064714278472050002 , -0.062084339131730311,
            -0.065197032815572477 , -0.072793895362578779 , -0.084959075171781003,
            -0.102539850131045997 , -0.127053585157696036 , -0.160791120691274606)

    elseif 0.4 <= x < 0.5 #Table 6
        t = x-0.45
        t = @horner(t  ,
            +1.375401971871116291 , -0.487202183273184837 , -0.153311701348540228 ,
            -0.111849444917027833 , -0.108840952523135768 , -0.122954223120269076 ,
            -0.152217163962035047 , -0.200495323642697339 , -0.276174333067751758 ,
            -0.393513114304375851 , -0.575754406027879147 , -0.860523235727239756 ,
            -1.308833205758540162)

    elseif 0.5 <= x < 0.6 #Table 7
        t = x-0.55
        t = @horner(t  ,
            +1.325024497958230082 , -0.521727647557566767 , -0.194906430482126213 ,
            -0.171623726822011264 , -0.202754652926419141 , -0.278798953118534762 ,
            -0.420698457281005762 , -0.675948400853106021 , -1.136343121839229244 ,
            -1.976721143954398261 , -3.531696773095722506 , -6.446753640156048150 ,
            -11.97703130208884026)
    elseif 0.6 <= x < 0.7 #Table 8
        t = x-0.65
        t = @horner(t,
            +1.270707479650149744 , -0.566839168287866583 , -0.262160793432492598 ,
            -0.292244173533077419 , -0.440397840850423189 , -0.774947641381397458 ,
            -1.498870837987561088 , -3.089708310445186667 , -6.667595903381001064 ,
            -14.89436036517319078 , -34.18120574251449024 , -80.15895841905397306 ,
            -191.3489480762984920 , -463.5938853480342030 , -1137.380822169360061)

    elseif 0.7 <= x < 0.8 #Table 9
        t = x-0.75
        t = @horner(t,
            +1.211056027568459525 , -0.630306413287455807 , -0.387166409520669145 ,
            -0.592278235311934603 , -1.237555584513049844 , -3.032056661745247199 ,
            -8.181688221573590762 , -23.55507217389693250 , -71.04099935893064956 ,
            -221.8796853192349888 , -712.1364793277635425 , -2336.125331440396407 ,
            -7801.945954775964673 , -26448.19586059191933 , -90799.48341621365251 ,
            -315126.0406449163424 , -1104011.344311591159)

    elseif 0.8 <= x < 0.85 #Table 10
        t = x-0.825
        t = @horner(t,
            +1.161307152196282836 , -0.701100284555289548 , -0.580551474465437362 ,
            -1.243693061077786614 , -3.679383613496634879 , -12.81590924337895775 ,
            -49.25672530759985272 , -202.1818735434090269 , -869.8602699308701437 ,
            -3877.005847313289571 , -17761.70710170939814 , -83182.69029154232061 ,
            -396650.4505013548170 , -1920033.413682634405)

    elseif 0.85 <= x < 0.9 #Table 11
        t = x-0.875
        t =  @horner(t,
            +1.124617325119752213 , -0.770845056360909542 , -0.844794053644911362 ,
            -2.490097309450394453 , -10.23971741154384360 , -49.74900546551479866 ,
            -267.0986675195705196 , -1532.665883825229947 , -9222.313478526091951 ,
            -57502.51612140314030 , -368596.1167416106063 , -2415611.088701091428 ,
            -16120097.81581656797 , -109209938.5203089915 , -749380758.1942496220 ,
            -5198725846.725541393 , -36409256888.12139973)

    elseif x >= 0.9
        td  = 1-x
        td1 = td-0.05

        kdm = @horner(td1 ,
            1.591003453790792180 , 0.416000743991786912 , 0.245791514264103415 ,
            0.179481482914906162 , 0.144556057087555150 , 0.123200993312427711 ,
            0.108938811574293531 , 0.098853409871592910 , 0.091439629201749751 ,
            0.085842591595413900 , 0.081541118718303215)
        edm = @horner(td1  ,
            +1.550973351780472328 , -0.400301020103198524 , -0.078498619442941939 ,
            -0.034318853117591992 , -0.019718043317365499 , -0.013059507731993309 ,
            -0.009442372874146547 , -0.007246728512402157 , -0.005807424012956090 ,
            -0.004809187786009338)
        hdm = kdm - edm
        km  = ellipk(Float64(x))
        #em =  km + (pi/2 - km*edm)/kdm
        em  = (pi/2 + km*hdm) / kdm   #to avoid precision loss near 1
        t   = em
    end
    if flag_is_m_neg
        return t * sqrt(1.0-m)
    else
        return t
    end
end

for f in (:ellipk,:ellipe)
    @eval begin
        ($f)(x::Float16)        = Float16(($f)(Float64(x)))
	    ($f)(x::Float32)        = Float32(($f)(Float64(x)))
        ($f)(x::Real)           = ($f)(float(x))
	    ($f)(x::AbstractFloat)  = throw(MethodError($f, (x, "")))
    end
end


text:
module QuickHullTests

using QuickHull
using Base.Test
using FixedSizeArrays

# test data
getpoints(N) = map(Vec, rand(N), rand(N))
getpoints(xs,ys) = map(Vec, xs, ys)
getxys(N::Int) = rand(N), rand(N)
getxys{T<:Vec{2}}(points::Vector{T}) = map(p->p[1],points), map(p->p[2],points)

function test_property(property)
    points = getpoints(50)
    hull = qhull(points)

    property(points, hull)

    xys = getxys(50)
    xyhull = qhull(xys...)

    property(getpoints(xys...), getpoints(xyhull...))

    nothing
end

# property: all points in hull are unique
test_property() do points, hull
    @test length(hull) == length(unique(hull))
end

# property: no point is to the right of any line through adjacent hull members
#           this property assumes counter-clockwise orientation of the hull,
#           and tests also for the pair of last-first points
test_property() do points, hull
    for (a,b) in zip(hull[1:end], vcat(hull[2:end], hull[1]))
        for q in points
            # skip points in the hull
            (q == a || q == b) && continue
            # q is to the right if |(b-a)x(q-a)| < 0 with implied z-coordinate = 0
            @test (b[1]-a[1])*(q[2]-a[2])-(b[2]-a[2])*(q[1]-a[1]) > 0
        end
    end
end

end


text:
using eFEMpart
include("AdvDiff_example_param.jl")

function ad_example(p::AD_Para)
  square = [p.xmin,p.xmax,p.ymin,p.ymax]
  mesh = squareMesh(square,p.Nintervals,p.basisorder)

  varname      = ""
  OperatorType = :AdvDiff2D

  dNodes = Dirichlet(:left,:top)
  nNodes = Neumann(:right,:bottom)

  dBCf   = Dirichlet((x,y) -> (x==p.xmin ? p.dirVal : 0.0))
  ff     = Forcing((x,y) -> 0.0)

  # parameter
  κ = 1/p.Pe
  windX(x,y) =  -y;  windY(x,y) = x
  wx = [windX.(mesh.xy[i].x,mesh.xy[i].y) for i=1:length(mesh.xy)]
  wy = [windY.(mesh.xy[i].x,mesh.xy[i].y) for i=1:length(mesh.xy)]
  param  = AdvDiffParam(wx,wy,κ)

  Nodes = [dNodes]
  bcfun = [dBCf,ff]

  prob = Problem(mesh,Nodes,bcfun,varname,OperatorType)
  sol = solve(prob,mesh,param)

  # save solution
  vtkname = Path(p.solFolder,p.solFile)
  sd = ScalarData(sol.u)
  sn = ScalarNames(p.solName)
  vd = VectorData([wx,wy])
  vn = VectorNames("velocity")
  vtksave(mesh,sd,sn,vd,vn,vtkname)

  return nothing
end

function main()
  p = AD_Para()
  ad_example(p)
end

main()

max_stars_repo_name: pseastham/eFEMpart


text:
module Evolutionary
using Random
    export Strategy, strategy, inverse, mutationwrapper,
           # ES mutations
           isotropic, anisotropic, isotropicSigma, anisotropicSigma,
           # GA mutations
           flip, domainrange, inversion, insertion, swap2, scramble, shifting,
           # ES recombinations
           average, marriage, averageSigma1, averageSigmaN,
           # GA recombinations
           singlepoint, twopoint, uniform,
           discrete, waverage, intermediate, line,
           pmx, ox1, cx, ox2, pos,
           # GA selections
           ranklinear, uniformranking, roulette, sus, tournament, truncation,
           # Optimization methods
           es, cmaes, ga

    const Strategy = Dict{Symbol,Any}
    const Individual = Union{Vector, Matrix, Function, Nothing}

    # Wrapping function for strategy
    function strategy(; kwargs...)
        result = Dict{Symbol,Any}()
        for (k, v) in kwargs
            result[k] = v
        end
        return result
    end

    # Inverse function for reversing optimization direction
    function inverseFunc(f::Function)
        function fitnessFunc(x::T) where {T <: Vector}
            return 1.0/(f(x)+eps())
        end
        return fitnessFunc
    end

    # Obtain individual
    function getIndividual(init::Individual, N::Int)
        if isa(init, Vector)
            @assert length(init) == N "Dimensionality of initial population must be $(N)"
            individual = init
        elseif isa(init, Matrix)
            @assert size(init, 1) == N "Dimensionality of initial population must be $(N)"
            populationSize = size(init, 2)
            individual = init[:, 1]
        elseif isa(init, Function) # Creation function
            individual = init(N)
        else
            individual = ones(N)
        end
        return  individual
    end

    # Collecting interim values
    function keep(interim, v, vv, col)
        if interim
            if !haskey(col, v)
                col[v] = typeof(vv)[]
            end
            push!(col[v], vv)
        end
    end

    # ES & GA recombination functions
    include("recombinations.jl")

    # ES & GA mutation functions
    include("mutations.jl")

    # GA selection functions
    include("selections.jl")

    # Evolution Strategy
    include("es.jl")
    include("cmaes.jl")

    # Genetic Algorithms
    include("ga.jl")

end


text:
#=
# 3D Poisson, Dirichlet bc
# Uses Dendro imported as a custom gen target
=#
if !@isdefined(Femshop)
    include("../Femshop.jl");
    using .Femshop
end
init_femshop("dendrotest");

useLog("dendrotestlog", level=3)

# default values (max_depth=6, wavelet_tol = 0.1, partition_tol = 0.3, solve_tol = 1e-6, max_iters = 100)
generateFor("target_dendro_cg.jl", params=(8, 0.001, 0.3, 0.000001, 500))

domain(3)
functionSpace(order=2)

u = variable("u")
testSymbol("v")

boundary(u, 1, DIRICHLET, "0")

# Write the weak form
coefficient("k1", "1.1")
coefficient("k2", "0.1")
coefficient("k3", "0.1")
coefficient("k4", "0.1")
coefficient("k5", "0.1")
f = coefficient("f", "-1.1*14*pi*pi*sin(3*pi*x)*sin(2*pi*y)*sin(pi*z) + 0.1*3*pi*cos(3*pi*x)*sin(2*pi*y)*sin(pi*z) + 0.1*2*pi*sin(3*pi*x)*cos(2*pi*y)*sin(pi*z) + 0.1*pi*sin(3*pi*x)*sin(2*pi*y)*cos(pi*z) - 0.1*sin(3*pi*x)*sin(2*pi*y)*sin(pi*z)")
weakForm(u, "-k1*dot(grad(u),grad(v)) + k2*deriv(u,1)*v + k3*deriv(u,2)*v + k4*deriv(u,3)*v - k5*u*v - f*v")

build_octree_with(f);
solve(u);

finalize_femshop()

max_stars_repo_name: paralab/femshop

max_issues_repo_name: paralab/femshop

max_forks_repo_name: paralab/femshop


text:
"""
Creates lossy converter model between AC and DC grid, assuming U_i is approximatley 1 numerically

```
pconv_ac[i] + pconv_dc[i] == a + b*pconv_ac
```
"""
function constraint_converter_losses(pm::_PM.AbstractDCPModel, n::Int,  i::Int, a, b, c, plmax)
    pconv_ac = _PM.var(pm, n, :pconv_ac, i)
    pconv_dc = _PM.var(pm, n, :pconv_dc, i)
    v = 1 #pu, assumption to approximate current
    cm_conv_ac = pconv_ac/v # can actually be negative, not a very nice model...
    if pm.setting["conv_losses_mp"] == true
        # _PM.con(pm, n, :conv_loss)[i] = JuMP.@constraint(pm.model, pconv_ac + pconv_dc == a + b*cm_conv_ac )
        # JuMP.@constraint(pm.model, pconv_ac + pconv_dc == a + b*cm_conv_ac )
        JuMP.@constraint(pm.model, pconv_ac + pconv_dc == 0 )

    else
        # _PM.con(pm, n, :conv_loss)[i] = JuMP.@constraint(pm.model, pconv_ac + pconv_dc >= a + b*cm_conv_ac )
        # _PM.con(pm, n, :conv_loss_aux)[i] = JuMP.@constraint(pm.model, pconv_ac + pconv_dc >= a - b*cm_conv_ac )
        # _PM.con(pm, n, :conv_loss_plmax)[i] = JuMP.@constraint(pm.model, pconv_ac + pconv_dc <= plmax)
        JuMP.@constraint(pm.model, pconv_ac + pconv_dc >= - a + b*cm_conv_ac )
        JuMP.@constraint(pm.model, pconv_ac + pconv_dc >= a - b*cm_conv_ac )
        JuMP.@constraint(pm.model, pconv_ac + pconv_dc <= plmax)
    end
end
"""
Converter transformer constraints

```
p_tf_fr == -btf*(v^2)/tm*(va-vaf)
p_tf_to == -btf*(v^2)/tm*(vaf-va)
```
"""
function constraint_conv_transformer(pm::_PM.AbstractDCPModel, n::Int,  i::Int, rtf, xtf, acbus, tm, transformer)
    ptf_fr = _PM.var(pm, n, :pconv_tf_fr, i)
    ptf_to = _PM.var(pm, n, :pconv_tf_to, i)
    # display(ptf_fr)
    # display(ptf_to)
    vaf = _PM.var(pm, n, :vaf, i)
    va = _PM.var(pm, n, :va, acbus)

    if transformer
        btf = imag(1/(im*xtf)) # classic DC approach to obtain susceptance form
        v = 1 # pu, assumption DC approximation
        # _PM.con(pm, n, :conv_tf_p_fr)[i] = JuMP.@constraint(pm.model, ptf_fr == -btf*(v^2)/tm*(va-vaf))
        # _PM.con(pm, n, :conv_tf_p_to)[i] = JuMP.@constraint(pm.model, ptf_to == -btf*(v^2)/tm*(vaf-va))
        JuMP.@constraint(pm.model, ptf_fr == -btf*(v^2)/tm*(va-vaf))
        JuMP.@constraint(pm.model, ptf_to == -btf*(v^2)/tm*(vaf-va))
    else
        # _PM.con(pm, n, :conv_tf_p_fr)[i] = JuMP.@constraint(pm.model, va == vaf)
        # _PM.con(pm, n, :conv_tf_p_to)[i] = JuMP.@constraint(pm.model, ptf_fr + ptf_to  == 0)
        JuMP.@constraint(pm.model, va == vaf)
        JuMP.@constraint(pm.model, ptf_fr + ptf_to  == 0)
    end
end
"""
Converter reactor constraints

```
p_pr_fr == -bc*(v^2)*(vaf-vac)
pconv_ac == -bc*(v^2)*(vac-vaf)
```
"""
function constraint_conv_reactor(pm::_PM.AbstractDCPModel, n::Int,  i::Int, rc, xc, reactor)
    ppr_fr = _PM.var(pm, n, :pconv_pr_fr, i)
    pconv_ac = _PM.var(pm, n, :pconv_ac, i)
    vaf = _PM.var(pm, n, :vaf, i)
    vac = _PM.var(pm, n, :vac, i)

    if reactor
        bc = imag(1/(im*xc))
        v = 1 # pu, assumption DC approximation
        # _PM.con(pm, n, :conv_pr_p)[i] = JuMP.@constraint(pm.model, ppr_fr == -bc*(v^2)*(vaf-vac))
        # _PM.con(pm, n, :conv_pr_p_to)[i] = JuMP.@constraint(pm.model, pconv_ac == -bc*(v^2)*(vac-vaf))
        JuMP.@constraint(pm.model, ppr_fr == -bc*(v^2)*(vaf-vac))
        JuMP.@constraint(pm.model, pconv_ac == -bc*(v^2)*(vac-vaf))
    else
        # _PM.con(pm, n, :conv_pr_p)[i] =  JuMP.@constraint(pm.model, vac == vaf)
        # _PM.con(pm, n, :conv_pr_p_to)[i] = JuMP.@constraint(pm.model, ppr_fr + pconv_ac  == 0)
        JuMP.@constraint(pm.model, vac == vaf)
        JuMP.@constraint(pm.model, ppr_fr + pconv_ac  == 0)
    end
end
"""
Converter filter constraints (no active power losses)
```
p_pr_fr + p_tf_to == 0
```
"""
function constraint_conv_filter(pm::_PM.AbstractDCPModel, n::Int,  i::Int, bv, filter)
    ppr_fr = _PM.var(pm, n, :pconv_pr_fr, i)
    ptf_to = _PM.var(pm, n, :pconv_tf_to, i)

    # _PM.con(pm, n, :conv_kcl_p)[i] = JuMP.@constraint(pm.model,   ppr_fr + ptf_to == 0 )
    JuMP.@constraint(pm.model,   ppr_fr + ptf_to == 0 )
end
"""
Converter current constraint (not applicable)
```
```
"""
function constraint_converter_current(pm::_PM.AbstractDCPModel, n::Int,  i::Int, Umax, Imax)
    # not used
end
function variable_dc_converter(pm::_PM.AbstractDCPModel; kwargs...)
    variable_converter_active_power(pm; kwargs...)
    variable_dcside_power(pm; kwargs...)
    variable_converter_filter_voltage(pm; kwargs...)
    variable_converter_internal_voltage(pm; kwargs...)
    variable_converter_to_grid_active_power(pm; kwargs...)

    variable_conv_transformer_active_power_to(pm; kwargs...)
    variable_conv_reactor_active_power_from(pm; kwargs...)
end

function variable_converter_filter_voltage(pm::_PM.AbstractDCPModel; kwargs...)
    variable_converter_filter_voltage_angle(pm; kwargs...)
end

function variable_converter_internal_voltage(pm::_PM.AbstractDCPModel; kwargs...)
    variable_converter_internal_voltage_angle(pm; kwargs...)
end
"""
Converter reactive power setpoint constraint (PF only, not applicable)
```
```
"""
function constraint_reactive_conv_setpoint(pm::_PM.AbstractDCPModel, n::Int,  i, qconv)
end
"""
Converter firing angle constraint (not applicable)
```
```
"""
function constraint_conv_firing_angle(pm::_PM.AbstractDCPModel, n::Int,  i::Int, S, P1, Q1, P2, Q2)
end
"""
Converter droop constraint (not applicable)
```
```
"""
function constraint_dc_droop_control(pm::_PM.AbstractDCPModel, n::Int,  i::Int, busdc_i, vref_dc, pref_dc, k_droop)
end
######################## TNEP Constraints #################
"""
Creates lossy converter model between AC and DC grid, assuming U_i is approximatley 1 numerically

```
pconv_ac[i] + pconv_dc[i] == a + b*pconv_ac
```
"""
function constraint_converter_losses_ne(pm::_PM.AbstractDCPModel, n::Int, i::Int, a, b, c, plmax)
    pconv_ac = _PM.var(pm, n, :pconv_ac_ne)[i]
    pconv_dc = _PM.var(pm, n, :pconv_dc_ne)[i]
    z = _PM.var(pm, n, :conv_ne)[i]
    v = 1 #pu, assumption to approximate current
    cm_conv_ac = pconv_ac/v # can actually be negative, not a very nice model...
    # binary to omit the no load losses, power omitted via binary constraint in varible bounds
    if pm.setting["conv_losses_mp"] == true
        # JuMP.@constraint(pm.model, pconv_ac + pconv_dc == a * z + b*cm_conv_ac )
        JuMP.@constraint(pm.model, pconv_ac + pconv_dc == 0 )

    else
        JuMP.@constraint(pm.model, pconv_ac + pconv_dc >= a * z + b*cm_conv_ac )
        JuMP.@constraint(pm.model, pconv_ac + pconv_dc >= a * z - b*cm_conv_ac )
        JuMP.@constraint(pm.model,  pconv_ac + pconv_dc <= plmax)
    end
     #pm.con[:nw][n][:conv_loss_ne][i] = @constraint(pm.model, pconv_ac + pconv_dc == a * z )
end

function constraint_conv_transformer_ne(pm::_PM.AbstractDCPModel, n::Int, i::Int, rtf, xtf, acbus, tm, transformer)
    ptf_fr = _PM.var(pm, n, :pconv_tf_fr_ne, i)
    ptf_to = _PM.var(pm, n, :pconv_tf_to_ne, i)
    #filter voltage angle
    vaf = _PM.var(pm, n, :vaf_ne, i)
    #acbus voltage angle
    va = _PM.var(pm, n, :va, acbus)
    va_du = _PM.var(pm, n, :va_du, i)
    JuMP.set_upper_bound(va, pi/2)
    JuMP.set_lower_bound(va, -pi/2)
    z = _PM.var(pm, n, :conv_ne, i)

    if transformer
        btf = imag(1/(im*xtf)) # classic DC approach to obtain susceptance form
        v = 1 # pu, assumption DC approximation
        JuMP.@constraint(pm.model, ptf_fr == -btf*(v^2)/tm*(va-vaf))
        JuMP.@constraint(pm.model, ptf_to == -btf*(v^2)/tm*(vaf-va))
    else
        JuMP.@constraint(pm.model, va_du == vaf)
        JuMP.@constraint(pm.model, ptf_fr + ptf_to  == 0)
    end
    relaxation_variable_on_off(pm.model, va, va_du, z)
    _IM.relaxation_equality_on_off(pm.model, va, va_du, z)
end
#
function constraint_conv_reactor_ne(pm::_PM.AbstractDCPModel, n::Int, i::Int, rc, xc, reactor)
    ppr_fr = _PM.var(pm, n, :pconv_pr_fr_ne, i)
    pconv = _PM.var(pm, n, :pconv_ac_ne, i)
    #filter voltage angle
    vaf = _PM.var(pm, n, :vaf_ne, i)
    #converter voltage angle
    vac = _PM.var(pm, n, :vac_ne, i)
    z = _PM.var(pm, n, :conv_ne)[i]

    if reactor
        bc = imag(1/(im*xc))
        v = 1 # pu, assumption DC approximation
        JuMP.@constraint(pm.model, ppr_fr == -bc*(v^2)*(vaf-vac))
        JuMP.@constraint(pm.model, pconv == -bc*(v^2)*(vac-vaf))
    else
        JuMP.@constraint(pm.model, vac == vaf)
        JuMP.@constraint(pm.model, ppr_fr + pconv == 0)
    end
end
#
function constraint_conv_filter_ne(pm::_PM.AbstractDCPModel, n::Int, i::Int, bv, filter)
    ppr_fr = _PM.var(pm, n, :pconv_pr_fr_ne, i)
    ptf_to = _PM.var(pm, n, :pconv_tf_to_ne, i)
    JuMP.@constraint(pm.model,   ppr_fr + ptf_to == 0 )
end
#
#
function constraint_converter_current_ne(pm::_PM.AbstractDCPModel, n::Int, i::Int, Umax, Imax)
     # not used
 end
#
#
function variable_dc_converter_ne(pm::_PM.AbstractDCPModel; kwargs...)
    variable_converter_ne(pm; kwargs...)
    variable_converter_active_power_ne(pm; kwargs...)
    variable_dcside_power_ne(pm; kwargs...)
    variable_converter_filter_voltage_ne(pm; kwargs...)
    variable_converter_internal_voltage_ne(pm; kwargs...)
    variable_converter_to_grid_active_power_ne(pm; kwargs...)

    variable_conv_transformer_active_power_to_ne(pm; kwargs...)
    variable_conv_reactor_active_power_from_ne(pm; kwargs...)
end

function variable_converter_filter_voltage_ne(pm::_PM.AbstractDCPModel; kwargs...)
    variable_converter_filter_voltage_angle_ne(pm; kwargs...)
end
#
#
function variable_converter_internal_voltage_ne(pm::_PM.AbstractDCPModel; kwargs...)
    variable_converter_internal_voltage_angle_ne(pm; kwargs...)
end

function variable_voltage_slack(pm::_PM.AbstractDCPModel; nw::Int=pm.cnw, bounded::Bool = true, report::Bool=false)
    va_ne = _PM.var(pm, nw)[:va_du] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_va_du",
    lower_bound = -2*pi,
    upper_bound = 2*pi,
    start = 0,
    )
    report && _IM.sol_component_value(pm, nw, :convdc_ne, :va, _PM.ids(pm, nw, :convdc_ne), va_ne)
end


text:
module MULTEM

using Configurations: @option, from_dict, to_dict, to_toml

import Configurations: toml_convert
import FieldMetadata: @description, description
import Base: convert

export from_dict, to_dict, to_toml
export MULTEMInput
export System, Theory, Simulation, Specimen, Instrument
# export ElectronPhononInteraction, ElectronSpecimenInteraction, IncidentWave, Illumination, Aberrations, DefocusSpreadFunction, ZeroDefocusReference, ObjectiveLens, SourceSpreadFunction, CondenserLens, ScanningArea, Detector, Rotation, Thickness, Amorphous, Crystal, XYSampling, Atom, OutputRegion, HCI, EELS, PET, EFTEM
# export Precision, Device, ElectronSpecimenInteractionModel, PotentialType, OperationMode, PhononModel, ThicknessType, PotentialSlicing, SimulationType, IncidentWaveType, IlluminationModel, IlluminationIncoherenceMode, ScanningType, ZeroDefocusType, DetectorType

@enum Precision begin
    eP_Float = 1
    eP_Double = 2
end

@enum Device begin
    eD_CPU = 1
    eD_GPU = 2
end

@enum ElectronSpecimenInteractionModel begin
    eESIM_Multislice = 1
    eESIM_Phase_Object = 2
    eESIM_Weak_Phase_Object = 3
end

@enum PotentialType begin
    ePT_Doyle_0_4 = 1
    ePT_Peng_0_4 = 2
    ePT_Peng_0_12 = 3
    ePT_Kirkland_0_12 = 4
    ePT_Weickenmeier_0_12 = 5
    ePT_Lobato_0_12 = 6
end

@enum OperationMode begin
    eOM_Normal = 1
    eOM_Advanced = 2
end

@enum PhononModel begin
    ePM_Still_Atom = 1
    ePM_Absorptive = 2
    ePM_Frozen_Phonon = 3
end

@enum ThicknessType begin
    eTT_Whole_Spec = 1
    eTT_Through_Thick = 2
    eTT_Through_Slices = 3
end

@enum PotentialSlicing begin
    ePS_Planes = 1
    ePS_dz_Proj = 2
    ePS_dz_Sub = 3
    ePS_Auto = 4
end

@enum SimulationType begin
    eST_STEM = 11
    eST_ISTEM = 12
    eST_CBED = 21
    eST_CBEI = 22

    eST_ED = 31
    eST_HRTEM = 32
    eST_PED = 41
    eST_HCTEM = 42

    eST_EWFS = 51
    eST_EWRS = 52
    eST_EELS = 61
    eST_EFTEM = 62

    eST_ProbeFS = 71
    eST_ProbeRS = 72
    eST_PPFS = 81
    eST_PPRS = 82

    eST_TFFS = 91
    eST_TFRS = 92
end

@enum IncidentWaveType begin
    eIWT_Plane_Wave = 1
    eIWT_Convergent_Wave = 2
    eIWT_User_Define_Wave = 3
    eIWT_Auto_Wave = 4
end

@enum IlluminationModel begin
    eIM_Coherente_Mode = 1
    eIM_Partial_Coherente_Mode = 2
    eIM_Transmission_Cross_Coefficient = 3
    eIM_Numerical_Integration = 4
end

@enum IlluminationIncoherenceMode begin
    eIIM_Temporal_And_Spatial_Incoherence = 1
    eIIM_Temporal_Incoherence = 2
    eIIM_Spatial_Incoherence = 3
end

@enum ScanningType begin
    eST_Line = 1
    eST_Area = 2
end

@enum ZeroDefocusType begin
    eZDT_First = 1
    eZDT_Middle = 2
    eZDT_Last = 3
    eZDT_User_Define = 4
end

@enum DetectorType begin
    eDT_Circular = 1
    eDT_Radial = 2
    eDT_Matrix = 3
end

@enum CenterType begin
    eCT_Geometric = 1
    eCT_User = 2
end

MULTEMEnum = Union{Precision, Device, ElectronSpecimenInteractionModel, PotentialType, OperationMode, PhononModel, ThicknessType, PotentialSlicing, SimulationType, IncidentWaveType, IlluminationModel, IlluminationIncoherenceMode, ScanningType, ZeroDefocusType, DetectorType}
# toml_convert(::Type, x::T) where T<:MULTEMEnum = string(x)
toml_convert(::Type, x::T) where T<:MULTEMEnum = Int(x)
convert(::Type{T}, x::Int) where T<:MULTEMEnum = T(x)

@description @option struct ElectronPhononInteraction
    model::PhononModel = ePM_Still_Atom | ""
    coherence_contribution::Bool = false | ""
    single_phonon_configuration::Bool = false | ""
    number_of_configuration::Int =
        1 |
        "???true: specific phonon configuration, false: number of frozen phonon configurations"
    dimensions::Int = 110 | "phonon dimensions (xyz)"
    seed::Int = 300183 | "Random seed(frozen phonon)"
end

@option struct ElectronSpecimenInteraction
    reverse_multislice::Bool = false
    model::ElectronSpecimenInteractionModel = eESIM_Multislice
    memory_size::Int = 0
    operation_mode::OperationMode = eOM_Normal
    potential_type::PotentialType = ePT_Lobato_0_12
end

"EM Theory"
@option struct Theory
    electron_phonon_interaction::ElectronPhononInteraction = ElectronPhononInteraction()
    electron_specimen_interaction::ElectronSpecimenInteraction =
        ElectronSpecimenInteraction()
end

@description @option struct IncidentWave
    acceleration_voltage::Real = 300.0 | "Acceleration Voltage (keV)"
    theta::Real = 0.0 | "Polar angle (°)"
    phi::Real = 0.0 | "Azimuthal angle (°)"
    type::IncidentWaveType = eIWT_Auto_Wave | ""
    psi::Real = 0.0 | ""
    x::Real = 0.0 | ""
    y::Real = 0.0 | ""
end

@option struct Illumination
    model::IlluminationModel = eIM_Transmission_Cross_Coefficient
    incoherence_mode::IlluminationIncoherenceMode = eIIM_Temporal_And_Spatial_Incoherence
end

@option struct Aberrations
    "defocus"
    c_10::Real = 14.0312
    c_12::Real = 0.0
    phi_12::Real = 0.0

    c_21::Real = 0.0
    phi_21::Real = 0.0
    c_23::Real = 0.0
    phi_23::Real = 0.0

    c_30::Real = 0.001
    c_32::Real = 0.0
    phi_32::Real = 0.0
    c_34::Real = 0.0
    phi_34::Real = 0.0

    c_41::Real = 0.0
    phi_41::Real = 0.0
    c_43::Real = 0.0
    phi_43::Real = 0.0
    c_45::Real = 0.0
    phi_45::Real = 0.0

    c_50::Real = 0.0
    c_52::Real = 0.0
    phi_52::Real = 0.0
    c_54::Real = 0.0
    phi_54::Real = 0.0
    c_56::Real = 0.0
    phi_56::Real = 0.0

    inner_aperture_angle::Real = 0.0
    outer_aperture_angle::Real = 21.0

end

@option struct ObjectiveLensDefocusSpreadFunction
    sigma::Real = 32.0
    integration_npoints::Int = 10
end

@option struct CondenserLensDefocusSpreadFunction
    # `a` and `beta` only for condenser lens???
    a::Real = 1.0
    beta::Real = 1.0

    sigma::Real = 32.0
    integration_npoints::Int = 10
end

@option struct SourceSpreadFunction
    a::Real = 1.0
    sigma::Real = 0.072
    beta::Real = 1.0
    radial_integration_npoints::Int = 4
    azimuthal_integration_npoints::Int = 4
end

@option struct ZeroDefocusReference
    type::ZeroDefocusType = eZDT_First
    plane::Real = 0.0
end

@option struct ObjectiveLens
    vortex_momentum::Real = 0.0
    aberrations::Aberrations = Aberrations()
    defocus_spread_function::ObjectiveLensDefocusSpreadFunction = ObjectiveLensDefocusSpreadFunction()
    zero_defocus_reference::ZeroDefocusReference = ZeroDefocusReference()
end

@option struct CondenserLens
    vortex_momentum::Real = 0.0
    aberrations::Aberrations = Aberrations()
    defocus_spread_function::CondenserLensDefocusSpreadFunction = CondenserLensDefocusSpreadFunction()
    source_spread_function::SourceSpreadFunction = SourceSpreadFunction()
    zero_defocus_reference::ZeroDefocusReference = ZeroDefocusReference()
end

@option struct ScanningArea
    type::ScanningType = eST_Area
    periodic_boundary_conditions::Bool = false
    square_pixel::Bool = true
    sampling_npoints::Int = 10
    x0::Real = 0.0
    y0::Real = 0.0
    xe::Real = 4.078
    ye::Real = 4.078
end

@option struct Detector
    type::DetectorType = eDT_Circular
    inner_angle::Real = 60.0
    outer_angle::Real = 180.0
    radial_sensitivity::Real = 0.0
    matrix_sensitivity::Real = 0.0
end

"Instrument"
@option struct Instrument
    incident_wave::IncidentWave = IncidentWave()
    illumination::Illumination = Illumination()
    objective_lens::ObjectiveLens = ObjectiveLens()
    scanning_area::ScanningArea = ScanningArea()
    condenser_lens::CondenserLens = CondenserLens()
    detector::Detector = Detector()
end

@option struct Rotation
    theta::Real = 0.0
    center_type::Real = 1
    center_p::AbstractVector{Real} = [0.0, 0.0, 0.0]
    u0::AbstractVector{Real} = [0.0, 0.0, 1.0]
end

@option struct Thickness
    thick_array::AbstractVector{Real} = [0.0]
    type::ThicknessType = eTT_Whole_Spec
end

@option struct Amorphous
    z0::Real = 0.0
    ze::Real = 0.0
    dz::Real = 2.0
end

@option struct Crystal
    na::Real = 1.0
    nb::Real = 1.0
    nc::Real = 1.0
    a::Real = 0.0
    b::Real = 0.0
    c::Real = 0.0
    x0::Real = 0.0
    y0::Real = 0.0
end

@option struct XYSampling
    nx::Int = 256
    ny::Int = 256
    band_width_limit::Real = 0.0
end

@option struct Atom
    Z::Int = 0
    x::Real = 0.0
    y::Real = 0.0
    z::Real = 0.0
    rms::Real = 0.085
    occupancy::Real = 1.0
    area::Int = 0
end

"Specimen"
@option struct Specimen
    rotation::Rotation = Rotation()
    thickness::Thickness = Thickness()
    amorphous::Amorphous = Amorphous()
    lx::Real = 10.0
    ly::Real = 10.0
    lz::Real = 10.0
    dz::Real = 0.25
    crystal::Crystal = Crystal()
    xy_sampling::XYSampling = XYSampling()
    potential_slicing::PotentialSlicing = ePS_Planes
    atoms::AbstractVector{Atom} = Vector{Atom}[]
end

@option struct OutputRegion
    ix_0::Real = 0.0
    iy_0::Real = 0.0
    ix_e::Real = 0.0
    iy_e::Real = 0.0
end

@option struct HCI
    theta::Real = 3.0
    nrotations::Real = 360.0
end

@option struct EELS
    atomic_number::Real = 79
    collection_angle::Real = 100.0
    m_selection::Real = 3.0
    E_loss::Real = 80.0
    channelling_type::Real = 1
end

@option struct PET
    theta::Real = 3.0
    nrotations::Real = 360.0
end

@option struct EFTEM
    atomic_number::Real = 79
    collection_angle::Real = 100.0
    m_selection::Real = 3.0
    E_loss::Real = 80.0
    channelling_type::Real = 1
end

"Simulation"
@option struct Simulation
    output_region::OutputRegion = OutputRegion()
    hci::HCI = HCI()
    eels::EELS = EELS()
    pet::PET = PET()
    eftem::EFTEM = EFTEM()
    type::SimulationType = eST_EWRS
end

"System"
@option struct System
    cpu_nthreads::Int = 1
    cpu_ncores::Int = 1
    device::Device = eD_GPU
    gpu_device::Int = 0
    precision::Precision = eP_Float
end

"MULTEM Input"
@option struct MULTEMInput
    "system"
    system::System = System()
    theory::Theory = Theory()
    simulation::Simulation = Simulation()
    specimen::Specimen = Specimen()
    instrument::Instrument = Instrument()
end

end

# @enum CondenserLensZeroDefocusReferenceType begin
#     eCLZDT_First = 1
#     eCLZDT_User_Define = 4
# end

# @enum ObjectiveLensZeroDefocusReferenceType begin
#     eOLZDT_First = 1
#     eOLZDT_Middle = 2
#     eOLZDT_Last = 3
#     eOLZDT_User_Define = 4
# end


text:
using StatsBase: sample

function select_col_sampling_strategy(colsample)
     # a sample of the features
    if 0 < colsample < 1
        col_sampling_bytree_strategy =
            function (features, args...; kwargs...)
                sample(features, floor(Int, length(features)*colsample), replace=false)
            end
    elseif colsample == 1
        col_sampling_bytree_strategy = (features, args...; kwargs...) -> features
    else
        error("colsample must be within [0, 1)")
    end

    col_sampling_bytree_strategy
end


text:
@inline function Flux(Cell,s,z,Def,j_tf,D,res0)
"""
Flux Transfer Function

Flux(Cell,s,z,Def)

"""

if Def == "Pos"
   Electrode = Cell.Pos #Electrode Length
else
   Electrode = Cell.Neg #Electrode Length
end

κ_eff = Cell.Const.κ*Electrode.ϵ_e^Electrode.κ_brug #Effective Electrolyte Conductivity
σ_eff = Electrode.σ*Electrode.ϵ_s^Electrode.σ_brug #Effective Electrode Conductivity

#Defining SOC
θ = Cell.Const.SOC * (Electrode.θ_100-Electrode.θ_0) + Electrode.θ_0

#Prepare for j0
ce0 = Cell.Const.ce0
cs_max = Electrode.cs_max
cs0 = cs_max*θ
α = Electrode.α

#Current Flux Density
if Cell.Const.CellTyp == "Doyle_94"
    κ = Electrode.k_norm/Electrode.cs_max/ce0^(1-α)
    j0 = κ*(ce0*(cs_max-cs0))^(1-α)*cs0^α
else
   j0 = Electrode.k_norm*(Cell.Const.ce0*cs0*(Electrode.cs_max-cs0))^(1-Electrode.α)
end

#Resistances
Rtot = R*Cell.Const.T /(j0*F^2) + Electrode.RFilm
#Rtot = R*Cell.Const.T/(j0*Cell.Const.CC_A*F) + Electrode.RFilm

#∂Uocp_Def
∂Uocp_elc = Cell.Const.∂Uocp(Def,θ)/cs_max

#Condensing Variable
ν = @. Electrode.L*sqrt((Electrode.as/σ_eff+Electrode.as/κ_eff)/(Rtot+∂Uocp_elc*(Electrode.Rs/(F*Electrode.Ds))*(tanh(Electrode.β)/(tanh(Electrode.β)-Electrode.β))))
ν_∞ = @. Electrode.L*sqrt(Electrode.as*((1/κ_eff)+(1/σ_eff))/(Rtot))

#Transfer Function
j_tf .= @. ν*(σ_eff*cosh(ν*z)+κ_eff*cosh(ν*(z-1)))/(Electrode.as*F*Electrode.L*Cell.Const.CC_A*(κ_eff+σ_eff)*sinh(ν))
D .= @. ν_∞*(σ_eff*cosh(ν_∞*z)+κ_eff*cosh(ν_∞*(z-1)))/(Electrode.as*F*Electrode.L*Cell.Const.CC_A*(κ_eff+σ_eff)*sinh(ν_∞))
zero_tf =ones(size(z,1))*1/(Cell.Const.CC_A*Electrode.as*F*Electrode.L)
j_tf[:,findall(s.==0)] .= zero_tf[:,findall(s.==0)]
res0 .= zeros(length(z))

if Def == "Pos"
   j_tf .= -j_tf
   D .= -D
end

end

max_stars_repo_name: BradyPlanden/LiBRA

max_forks_repo_name: BradyPlanden/LiBRA


text:
import Base.Libc: strerror

mean(x::Number, y::Number) = (x + y) / 2

const ESC = '\x1B'

# Encoded terminal modes from SSH.jl
const NCCS = Sys.islinux() ? 32 : 20
const tcflag_t = Sys.islinux() ? Cuint : Culong
const speed_t = tcflag_t

mutable struct termios
    c_iflag::tcflag_t
    c_oflag::tcflag_t
    c_cflag::tcflag_t
    c_lflag::tcflag_t
    @static if Sys.islinux()
        c_line::UInt8
    end
    c_cc::NTuple{NCCS, UInt8}
    c_uispeed::speed_t
    c_ospeed::speed_t
end
TERM = Ref{termios}(
                    @static if Sys.islinux() termios(
                            0,
                            0,
                            0,
                            0,
                            0,
                            Tuple([0 for _ in 1:NCCS]),
                            0,
                            0
                           )
                    else
                                             termios(
                            0,
                            0,
                            0,
                            0,
                            Tuple([0 for _ in 1:NCCS]),
                            0,
                            0
                           )
                    end
                   );
RESTORE = deepcopy(TERM);

const ICANON = Sys.islinux() ? 0o0000002 : 0x00000100
const ECHO = Sys.islinux()   ? 0o0000010 : 0x00000008
const TCSANOW = 0
const O_RDWR = Base.Filesystem.JL_O_RDWR
const O_NOCTTY = Base.Filesystem.JL_O_NOCTTY
const OS_HANDLE = Base.OS_HANDLE

function getColor(t::Terminals.TTYTerminal)
    fd = ccall(:jl_uv_file_handle, Base.OS_HANDLE, (Ptr{Cvoid},), stdin.handle)
    systemerror("tcgetattr",
                ccall(:tcgetattr, Cint, (Cint, Ptr{Cvoid}), fd, TERM) == -1)
    systemerror("tcgetattr",
                ccall(:tcgetattr, Cint, (Cint, Ptr{Cvoid}), fd, RESTORE) == -1)
    TERM[].c_lflag &= ~(ICANON|ECHO)
    systemerror("tcsetattr",
                ccall(:tcsetattr, Cint, (Cint, Cint, Ptr{Cvoid}), fd, TCSANOW, TERM) == -1
               )
    write(stdout, "\e]11;?\a")
    buf = UInt8[]
    while true
        ch = read(stdin, 1)[1]
        if Char(ch) == ';'
            break
        end
    end
    @assert Char(read(stdin, 1)[1]) == 'r'
    @assert Char(read(stdin, 1)[1]) == 'g'
    @assert Char(read(stdin, 1)[1]) == 'b'
    @assert Char(read(stdin, 1)[1]) == ':'
    r = parse(UInt8, "0x$(String(read(stdin, 2)))")
    @assert Char(read(stdin, 1)[1]) == '/'
    g = parse(UInt8, "0x$(String(read(stdin, 2)))")
    @assert Char(read(stdin, 1)[1]) == '/'
    b = parse(UInt8, "0x$(String(read(stdin, 2)))")
    systemerror("tcsetattr",
                ccall(:tcsetattr, Cint, (Cint, Cint, Ptr{Cvoid}), fd, TCSANOW, RESTORE) == -1
               )

    return r, g, b
end

function getXY(t::Terminals.TTYTerminal)
    fd = ccall(:jl_uv_file_handle, Base.OS_HANDLE, (Ptr{Cvoid},), stdin.handle)
    systemerror("tcgetattr",
                ccall(:tcgetattr, Cint, (Cint, Ptr{Cvoid}), fd, TERM) == -1)
    systemerror("tcgetattr",
                ccall(:tcgetattr, Cint, (Cint, Ptr{Cvoid}), fd, RESTORE) == -1)
    TERM[].c_lflag &= ~(ICANON|ECHO)
    systemerror("tcsetattr",
                ccall(:tcsetattr, Cint, (Cint, Cint, Ptr{Cvoid}), fd, TCSANOW, TERM) == -1
               )
    write(stdout, "$(Terminals.CSI)6n")
    buf = UInt8[]
    while true
        ch = read(stdin, 1)[1]
        if Char(ch) == 'R'
            break
        end
        push!(buf, ch)
    end
    systemerror("tcsetattr",
                ccall(:tcsetattr, Cint, (Cint, Cint, Ptr{Cvoid}), fd, TCSANOW, RESTORE) == -1
               )

    r, c = split(String(buf[3:end]), ";")
    return parse(Int, c), parse(Int, r)
end

abstract type Terminal end
abstract type JPEG end

size(::Type{Terminal}) = Base.displaysize(stdout) |> reverse

function size(::Type{JPEG}, filename)
    fhandle = open(filename)
    seek(fhandle, 0) # Read 0xff next
    size = 2
    ftype = 0
    while ! ( 0xc0 <= ftype <= 0xcf )
        read(fhandle, size)
        byte = read(fhandle, 1)
        while byte[1] == 0xff
            byte = read(fhandle, 1)
        end
        ftype = byte[1]
        size = read(fhandle, 2)[2] - 2
    end
    read(fhandle, 1)
    h1, h2, w1, w2 = read(fhandle, 4)
    return Int(UInt16(w1)<<8 + w2), Int(UInt16(h1)<<8 + h2)
end

canvassize() = size(Terminal)
pos() = pos(TERMINAL)
getXY() = getXY(TERMINAL)
getColor() = getColor(TERMINAL)
getX() = getX(TERMINAL)
getY() = getY(TERMINAL)
getW() = canvassize()[1]
getH() = canvassize()[2]
cmove(x, y) = cmove(TERMINAL, x, y)
clear() = clear(TERMINAL)
cmove_bottom() = cmove(1, getH() - 2)

pos(t::Terminals.TTYTerminal) = (getX(t), getY(t))

function getX(t::Terminals.TTYTerminal)
    x, y = getXY(t)
    return x
end

function getY(t::Terminals.TTYTerminal)
    x, y = getXY(t)
    return y
end

@eval cmove(t::Terminals.TTYTerminal, x::Int, y::Int) = print("$(Terminals.CSI)$(y);$(x)H")
@eval function clear(t::Terminals.TTYTerminal)
    r = getH()
    print("$(Terminals.CSI)$(r);1H")
    print("$(Terminals.CSI)2J")
end

wrap(s) = wrap(s, round(Int, getW() * 3 / 4))

function wrap(s, w::Int)
    length(s) < w && return String[s]
    i = findprev(" ", s, w)
    i = (i == nothing) ? w : i[1]
    first, remaining = s[1:i], s[i+1:end]
    first = first
    remaining = remaining
    return vcat(String[first], wrap(remaining, w))
end


text:
using TrajectoryOptimization
using Distributions
using Test
using BenchmarkTools
using LinearAlgebra
using Random
using StaticArrays
using SparseArrays
using Ipopt
using Logging
using ForwardDiff
const TO = TrajectoryOptimization

@testset "Logging" begin
    include("logger_tests.jl")
end

@testset "Dynamics" begin
    include("models_test.jl")
    include("dynamics_constraints.jl")
end

@testset "Full Solves" begin
    include("car_tests.jl")
    TEST_TIME = false  # don't test timing results
    include("benchmark_solves.jl")
end

@testset "Solvers" begin
    include("solver_options.jl")
end

@testset "Constraints" begin
    include("constraint_tests.jl")
end

# @testset "Rotations" begin
#     include("rotations_tests.jl")
#     include("retraction_maps.jl")
# end

@testset "Costs" begin
    include("cost_tests.jl")
end

@testset "Controllers" begin
    include("controllers_test.jl")
end


text:
# The old-gdemo model.
@model function gdemo(x, y)
  s ~ InverseGamma(2, 3)
  m ~ Normal(0, sqrt(s))
  x ~ Normal(m, sqrt(s))
  y ~ Normal(m, sqrt(s))
  return s, m
end

@model function gdemo_d()
  s ~ InverseGamma(2, 3)
  m ~ Normal(0, sqrt(s))
  1.5 ~ Normal(m, sqrt(s))
  2.0 ~ Normal(m, sqrt(s))
  return s, m
end

gdemo_default = gdemo_d()

@model function MoGtest(D)
    mu1 ~ Normal(1, 1)
    mu2 ~ Normal(4, 1)
    z1 ~ Categorical(2)
    if z1 == 1
        D[1] ~ Normal(mu1, 1)
    else
        D[1] ~ Normal(mu2, 1)
    end
    z2 ~ Categorical(2)
    if z2 == 1
        D[2] ~ Normal(mu1, 1)
    else
        D[2] ~ Normal(mu2, 1)
    end
    z3 ~ Categorical(2)
    if z3 == 1
        D[3] ~ Normal(mu1, 1)
    else
        D[3] ~ Normal(mu2, 1)
    end
    z4 ~ Categorical(2)
    if z4 == 1
        D[4] ~ Normal(mu1, 1)
    else
        D[4] ~ Normal(mu2, 1)
    end
    z1, z2, z3, z4, mu1, mu2
end

MoGtest_default = MoGtest([1.0 1.0 4.0 4.0])

# Declare empty model to make the Sampler constructor work.
@model empty_model() = x = 1


text:
# This example shows the pedestrain and crosswalk scenario -- nominal and failure scenario
using AdversarialDriving
using POMDPs, POMDPPolicies, POMDPSimulators

sut_agent = BlinkerVehicleAgent(get_ped_vehicle(id=1, s=5., v=15.), TIDM(ped_TIDM_template, noisy_observations = true))
adv_ped = NoisyPedestrianAgent(get_pedestrian(id=2, s=7., v=2.0), AdversarialPedestrian())
mdp = AdversarialDrivingMDP(sut_agent, [adv_ped], ped_roadway, 0.1)

null_action = Disturbance[PedestrianControl()]
noisy_action = Disturbance[PedestrianControl(noise = Noise((-10.,0.), -2))]

# Nominal Behavior
hist = POMDPSimulators.simulate(HistoryRecorder(), mdp, FunctionPolicy((s) -> null_action))
scenes_to_gif(state_hist(hist), mdp.roadway, "ped_crosswalk_nominal.gif", others = [crosswalk])

# Behavior with noise
hist = POMDPSimulators.simulate(HistoryRecorder(), mdp, FunctionPolicy((s) -> noisy_action))
scenes_to_gif(state_hist(hist), mdp.roadway, "ped_crosswalk_failure.gif", others = [crosswalk])


text:
import SMC
import Distributions
using DataFrames
include("hmm_serialization.jl")
include("schema.jl")

@everywhere begin
    using SMC
    using Distributions
    include("smc_samplers.jl")
    include("../aide.jl")
end

function generate_aide_estimates(hmm::HiddenMarkovModel,
                                 observations::Vector{Int},
                                 num_particles_list::Vector{Int},
                                 num_metainference_list::Vector{Int},
                                 num_replicates::Int)
    data = DataFrame()
    num_particles_column = Int[]
    num_metainference_column = Int[]
    proposal_name_column = String[]
    aide_estimate_column = Float64[]
    aide_stderr_column = Float64[]
    gold_standard_name_column = String[]

    exact_sampler = HMMExactSampler(hmm, observations)
    gold_standard_sampler = make_optimal_proposal_smc_sampler(hmm, observations, 1000)

    for num_particles in num_particles_list
        for num_metainference in num_metainference_list
            println("generating data for num_partices=$num_particles, num_metainference=$num_metainference...")

            prior_smc_sampler = make_prior_proposal_smc_sampler(hmm, observations, num_particles)
            kls = aide(exact_sampler, prior_smc_sampler, 1, num_metainference, num_replicates, num_replicates)
            push!(num_particles_column, num_particles)
            push!(num_metainference_column, num_metainference)
            push!(proposal_name_column, PRIOR_PROPOSAL_NAME)
            push!(aide_estimate_column, mean(kls))
            push!(aide_stderr_column, std(kls)/sqrt(length(kls)))
            push!(gold_standard_name_column, EXACT_GOLD_STANDARD)

            prior_smc_sampler = make_prior_proposal_smc_sampler(hmm, observations, num_particles)
            kls = aide(gold_standard_sampler, prior_smc_sampler, 1, num_metainference, num_replicates, num_replicates)
            push!(num_particles_column, num_particles)
            push!(num_metainference_column, num_metainference)
            push!(proposal_name_column, PRIOR_PROPOSAL_NAME)
            push!(aide_estimate_column, mean(kls))
            push!(aide_stderr_column, std(kls)/sqrt(length(kls)))
            push!(gold_standard_name_column, APPROXIMATE_GOLD_STANDARD)

            optimal_smc_sampler = make_optimal_proposal_smc_sampler(hmm, observations, num_particles)
            kls = aide(exact_sampler, optimal_smc_sampler, 1, num_metainference, num_replicates, num_replicates)
            push!(num_particles_column, num_particles)
            push!(num_metainference_column, num_metainference)
            push!(proposal_name_column, OPTIMAL_PROPOSAL_NAME)
            push!(aide_estimate_column, mean(kls))
            push!(aide_stderr_column, std(kls)/sqrt(length(kls)))
            push!(gold_standard_name_column, EXACT_GOLD_STANDARD)

            optimal_smc_sampler = make_optimal_proposal_smc_sampler(hmm, observations, num_particles)
            kls = aide(gold_standard_sampler, optimal_smc_sampler, 1, num_metainference, num_replicates, num_replicates)
            push!(num_particles_column, num_particles)
            push!(num_metainference_column, num_metainference)
            push!(proposal_name_column, OPTIMAL_PROPOSAL_NAME)
            push!(aide_estimate_column, mean(kls))
            push!(aide_stderr_column, std(kls)/sqrt(length(kls)))
            push!(gold_standard_name_column, APPROXIMATE_GOLD_STANDARD)
        end
    end
    data[COL_NUM_PARTICLES] = num_particles_column
    data[COL_NUM_METAINFERENCE] = num_metainference_column
    data[COL_PROPOSAL_NAME] = proposal_name_column
    data[COL_AIDE_ESTIMATE] = aide_estimate_column
    data[COL_AIDE_STDERR] = aide_stderr_column
    data[COL_GOLD_STANDARD_NAME] = gold_standard_name_column
    return data
end

# do experiment
data_dir = "data"
plot_dir = "plots"

# load HMM and observations
hmm = load_hmm("$data_dir/hmm.json")
num_states = num_states(hmm)
num_obs = num_observations(hmm)
observations = load_observations("$data_dir/observations.json")
num_steps = length(observations)

# do AIDE experiment, save data to CSV file
num_particles_list = [1, 3, 10, 30, 100]
num_metainference_list = [1, 100]
num_replicates = 100
aide_estimates = generate_aide_estimates(hmm, observations, num_particles_list,
                                         num_metainference_list, num_replicates)
writetable("$data_dir/aide_estimates.csv", aide_estimates)
println("done!")


text:
"""
    MinTimeGlobal <: AbstractSequencingAlgorithm

Algorithm that decides the sequence based on minimizing the time from the current
stop to the next stop.
"""
struct MinTimeGlobal{F} <: AbstractSequencingAlgorithm
    minfunc::F
end

MinTimeGlobal() = MinTimeGlobal(_min_length_path_exhaustive)

function solve(alg::MinTimeGlobal, route::AbstractRoute, times::TravelTimes; initial_stop=station(route))
    T = times.travel_times
    return solve(alg, route, T, initial_stop=initial_stop)
end

function solve(alg::MinTimeGlobal, route::AbstractRoute, T::AbstractMatrix;
               initial_stop=station(route))

    S = stops(route)
    num_stops = length(S)

    # salida
    route_id = name(route)
    seq = Vector{BareStop}()
    out = Sequence(route_id, seq)

    # para cada parada, encontrar la siguiente con el menor tiempo global
    if length(S) == 1
        push!(out, BareStop(name(S[1])))
        return out
    end

    _, result = alg.minfunc(T, S)
    for r in result
        push!(out, BareStop(name(r)))
    end
    return out
end


text:
using UnPack, OrderedCollections, ForwardDiff

# TODO
# 1. Create model object
# 2. Run code in MATLAB and generate comparison test output for all functions, say at the deterministic steady state
# 3. Test ForwardDiff works properly for all cases, for both in place and non in place models
# 4. Label equations for μ, ξ so creating Γ₅ and Γ₆ makes more sense
#=function RiskAdjustedLinearization(μ::M, Λ::L, Σ::S, ξ::X, ccgf::CF,
                                   z::AbstractVector{T}, y::AbstractVector{T}, Ψ::AbstractMatrix{T},
                                   Nε::Int = -1) where {T <: Number, M <: Function, L <: Function,
                                                        S <: Function, X <: Function, CF <: Function}
end=#
mutable struct WachterDisasterRisk{T <: Real}
    μₐ::T
    σₐ::T
    ν::T
    δ::T
    ρₚ::T
    pp::T
    ϕₚ::T
    ρ::T
    γ::T
    β::T
end

function WachterDisasterRisk(; μₐ::T = .0252 / 4., σₐ::T = .02 / sqrt(4.), ν::T = .3, δ::T = 0., ρₚ::T = .08^(1. / 4.), pp::T = .0355 / 4.,
                             ϕₚ::T = .0114 / 4. / (.02 / sqrt(4.)) / sqrt(.0355 / 4.), ρ::T = 2.0, γ::T = 3.0,
                             β::T = exp(-.012 / 4.)) where {T <: Real}
    return WachterDisasterRisk{T}(μₐ, σₐ, ν, δ, ρₚ, pp, ϕₚ, ρ, γ, β)
end

function inplace_wachter_disaster_risk(m::WachterDisasterRisk{T}) where {T <: Real}
    @unpack μₐ, σₐ, ν, δ, ρₚ, pp, ϕₚ, ρ, γ, β = m

    @assert ρ != 1. # Forcing ρ to be non-unit for this example

    S  = OrderedDict{Symbol, Int}(:p => 1,  :εc => 2, :εξ => 3) # State variables
    J  = OrderedDict{Symbol, Int}(:vc => 1, :xc => 2, :rf => 3) # Jump variables
    SH = OrderedDict{Symbol, Int}(:εₚ => 1, :εc => 2, :εξ => 3) # Exogenous shocks
    Nz = length(S)
    Ny = length(J)
    Nε = length(SH)

    function μ(F, z, y) # the equation labeling could be wrong
        F_type    = eltype(F)
        F[S[:p]]  = (1 - ρₚ) * pp + ρₚ * z[S[:p]]
        F[S[:εc]] = zero(F_type)
        F[S[:εξ]] = zero(F_type)
    end

    function ξ(F, z, y) # the equation labeling could be wrong
        F[J[:vc]] = log(β) - γ * μₐ + γ * ν * z[S[:p]] - (ρ - γ) * y[J[:xc]] + y[J[:rf]]
        F[J[:xc]] = log(1. - β + β * exp((1. - ρ) * y[J[:xc]])) - (1. - ρ) * y[J[:vc]]
        F[J[:rf]] = (1. - γ) * (μₐ - ν * z[S[:p]] - y[J[:xc]])
    end

    Λ = zeros(T, Nz, Ny)

    function Σ(F, z)
        F_type = eltype(F)
        F[SH[:εₚ], SH[:εₚ]] = sqrt(z[S[:p]]) * ϕₚ * σₐ
        F[SH[:εc], SH[:εc]] = one(F_type)
        F[SH[:εξ], SH[:εξ]] = one(F_type)
    end

    function ccgf(F, α, z)
        F .= .5 .* α[:, 1].^2 + .5 * α[:, 2].^2 + (exp.(α[:, 3] + α[:, 3].^2 .* δ^2 ./ 2.) .- 1. - α[:, 3]) * z[S[:p]]
    end

    Γ₅ = zeros(T, Ny, Nz)
    Γ₅[J[:vc], S[:εc]] = (-γ * σₐ)
    Γ₅[J[:vc], S[:εξ]] = (γ * ν)
    Γ₅[J[:rf], S[:εc]] = (1. - γ) * σₐ
    Γ₅[J[:rf], S[:εξ]] = -(1. - γ) * ν

    Γ₆ = zeros(T, Ny, Ny)
    Γ₆[J[:vc], J[:vc]] = (ρ - γ)
    Γ₆[J[:rf], J[:vc]] = (1. - γ)

    z = [pp, 0., 0.]
    xc_sss = log((1. - β) / (exp((1. - ρ) * (ν * pp - μₐ)) - β)) / (1. - ρ)
    vc_sss = xc_sss + ν * pp - μₐ
    y = [vc_sss, xc_sss, -log(β) + γ * (μₐ - ν * pp) - (ρ - γ) * (vc_sss - xc_sss)]
    Ψ = zeros(T, Ny, Nz)
    return RiskAdjustedLinearization(μ, Λ, Σ, ξ, Γ₅, Γ₆, ccgf, z, y, Ψ, Nε)
end

function outofplace_wachter_disaster_risk(m::WachterDisasterRisk{T}) where {T}
    @unpack μₐ, σₐ, ν, δ, ρₚ, pp, ϕₚ, ρ, γ, β = m

    @assert ρ != 1. # Forcing ρ to be non-unit for this example

    S  = OrderedDict{Symbol, Int}(:p => 1,  :εc => 2, :εξ => 3) # State variables
    J  = OrderedDict{Symbol, Int}(:vc => 1, :xc => 2, :rf => 3) # Jump variables
    SH = OrderedDict{Symbol, Int}(:εₚ => 1, :εc => 2, :εξ => 3) # Exogenous shocks
    Nz = length(S)
    Ny = length(J)
    Nε = length(SH)

    function μ(z, y) # the equation labeling could be wrong
        return [(1 - ρₚ) * pp + ρₚ * z[S[:p]], 0., 0.]
    end

    function ξ(z, y) # the equation labeling could be wrong
        F = RiskAdjustedLinearizations.dualvector(y, z)

        F[J[:vc]] = log(β) - γ * μₐ + γ * ν * z[S[:p]] - (ρ - γ) * y[J[:xc]] + y[J[:rf]]
        F[J[:xc]] = log(1. - β + β * exp((1. - ρ) * y[J[:xc]])) - (1. - ρ) * y[J[:vc]]
        F[J[:rf]] = (1. - γ) * (μₐ - ν * z[S[:p]] - y[J[:xc]])
        return F
    end

    Λ = zeros(T, Nz, Nz)

    function Σ(z)
        F = zeros(eltype(z), Nz, Nz)
        F[SH[:εₚ], SH[:εₚ]] = sqrt(z[S[:p]]) * ϕₚ * σₐ
        F[SH[:εc], SH[:εc]] = 1.
        F[SH[:εξ], SH[:εξ]] = 1.
        return F
    end

    ccgf(α, z) = .5 * α[:, 1].^2 + .5 * α[:, 2].^2 + (exp.(α[:, 3] + α[:, 3].^2 * δ.^2 ./ 2.) .- 1. - α[:, 3]) * z[S[:p]]

    Γ₅ = zeros(T, Ny, Nz)
    Γ₅[J[:vc], S[:εc]] = (-γ * σₐ)
    Γ₅[J[:vc], S[:εξ]] = (γ * ν)
    Γ₅[J[:rf], S[:εc]] = (1. - γ) * σₐ
    Γ₅[J[:rf], S[:εξ]] = -(1. - γ) * ν

    Γ₆ = zeros(T, Ny, Ny)
    Γ₆[J[:vc], J[:vc]] = (ρ - γ)
    Γ₆[J[:rf], J[:vc]] = (1. - γ)

    z = [pp, 0., 0.]
    xc_sss = log((1. - β) / (exp((1. - ρ) * (ν * pp - μₐ)) - β)) / (1. - ρ)
    vc_sss = xc_sss + ν * pp - μₐ
    y = [vc_sss, xc_sss, -log(β) + γ * (μₐ - ν * pp) - (ρ - γ) * (vc_sss - xc_sss)]
    Ψ = zeros(T, Ny, Nz)
    return RiskAdjustedLinearization(μ, Λ, Σ, ξ, Γ₅, Γ₆, ccgf, z, y, Ψ, Nε)
end


text:
import MacroTools: splitdef, combinedef, isexpr, postwalk

# XXX: Proper errors
function __kernel(expr)
    def = splitdef(expr)
    name = def[:name]
    args = def[:args]

    constargs = Array{Bool}(undef, length(args))
    for (i, arg) in enumerate(args)
        if isexpr(arg, :macrocall)
            if arg.args[1] === Symbol("@Const")
                # arg.args[2] is a LineInfo node
                args[i] = arg.args[3] # strip @Const
                constargs[i] = true
                continue
            end
        end
        constargs[i] = false
    end

    # create two functions
    # 1. GPU function
    # 2. CPU function with work-group loops inserted
    #
    # Without the deepcopy we might accidentially modify expr shared between CPU and GPU
    def_cpu = deepcopy(def)
    def_gpu = deepcopy(def)

    def_cpu[:name] = cpu_name = Symbol(:cpu_, name)
    def_gpu[:name] = gpu_name = Symbol(:gpu_, name)

    transform_cpu!(def_cpu, constargs)
    transform_gpu!(def_gpu, constargs)

    cpu_function = combinedef(def_cpu)
    gpu_function = combinedef(def_gpu)

    # create constructor functions
    constructors = quote
        if !@isdefined($name)
            $name(dev::$Device) = $name(dev, $DynamicSize(), $DynamicSize())
            $name(dev::$Device, size) = $name(dev, $StaticSize(size), $DynamicSize())
            $name(dev::$Device, size, range) = $name(dev, $StaticSize(size), $StaticSize(range))
            function $name(::Device, ::S, ::NDRange) where {Device<:$CPU, S<:$_Size, NDRange<:$_Size}
                return $Kernel{Device, S, NDRange, typeof($cpu_name)}($cpu_name)
            end
            function $name(::Device, ::S, ::NDRange) where {Device<:$GPU, S<:$_Size, NDRange<:$_Size}
                return $Kernel{Device, S, NDRange, typeof($gpu_name)}($gpu_name)
            end
        end
    end

    return Expr(:block, esc(cpu_function), esc(gpu_function), esc(constructors))
end

# The easy case, transform the function for GPU execution
# - mark constant arguments by applying `constify`.
function transform_gpu!(def, constargs)
    let_constargs = Expr[]
    for (i, arg) in enumerate(def[:args])
        if constargs[i]
            push!(let_constargs, :($arg = $constify($arg)))
        end
    end

    body = quote
        if $__validindex()
            $(def[:body])
        end
        return nothing
    end
    def[:body] = Expr(:let,
        Expr(:block, let_constargs...),
        body,
    )
end

# The hard case, transform the function for CPU execution
# - mark constant arguments by applying `constify`.
# - insert aliasscope markers
# - insert implied loop bodys
#   - handle indicies
#   - hoist workgroup definitions
#   - hoist uniform variables
function transform_cpu!(def, constargs)
    let_constargs = Expr[]
    for (i, arg) in enumerate(def[:args])
        if constargs[i]
            push!(let_constargs, :($arg = $constify($arg)))
        end
    end
    new_stmts = Expr[]
    body = MacroTools.flatten(def[:body])
    push!(new_stmts, Expr(:aliasscope))
    append!(new_stmts, split(body.args))
    push!(new_stmts, Expr(:popaliasscope))
    push!(new_stmts, :(return nothing))
    def[:body] = Expr(:let,
        Expr(:block, let_constargs...),
        Expr(:block, new_stmts...)
    )
end

struct WorkgroupLoop
    indicies :: Vector{Any}
    stmts :: Vector{Any}
    allocations :: Vector{Any}
    private :: Vector{Any}
end

is_sync(expr) = @capture(expr, @synchronize() | @synchronize(a_))

function is_scope_construct(expr::Expr)
    expr.head === :block # ||
    # expr.head === :let
end

function find_sync(stmt)
    result = false
    postwalk(stmt) do expr
        result |= is_sync(expr)
        expr
    end
    result
end

# TODO proper handling of LineInfo
function split(stmts,
               indicies = Any[], private=Any[])
    # 1. Split the code into blocks separated by `@synchronize`
    # 2. Aggregate `@index` expressions
    # 3. Hoist allocations
    # 4. Hoist uniforms

    current     = Any[]
    allocations = Any[]
    new_stmts   = Any[]
    for stmt in stmts
        has_sync = find_sync(stmt)
        if has_sync
            loop = WorkgroupLoop(deepcopy(indicies), current, allocations, deepcopy(private))
            push!(new_stmts, emit(loop))
            allocations = Any[]
            current     = Any[]
            is_sync(stmt) && continue

            # Recurse into scope constructs
            # TODO: This currently implements hard scoping
            #       probably need to implemet soft scoping
            #       by not deepcopying the environment.
            recurse(x) = x
            function recurse(expr::Expr)
                expr = unblock(expr)
                if is_scope_construct(expr) && any(is_sync, expr.args)
                    new_args = unblock(split(expr.args, deepcopy(indicies), deepcopy(private)))
                    return Expr(expr.head, new_args...)
                else
                    return Expr(expr.head, map(recurse, expr.args)...)
                end
            end
            push!(new_stmts, recurse(stmt))
            continue
        end

        if @capture(stmt, @uniform x_)
            push!(allocations, stmt)
            continue
        elseif @capture(stmt, lhs_ = rhs_ | (vs__, lhs_ = rhs_))
            if @capture(rhs, @index(args__))
                push!(indicies, stmt)
                continue
            elseif @capture(rhs, @localmem(args__) | @uniform(args__) )
                push!(allocations, stmt)
                continue
            elseif @capture(rhs, @private(args__))
                push!(allocations, stmt)
                push!(private, lhs)
                continue
            end
        end

        push!(current, stmt)
    end

    # everything since the last `@synchronize`
    if !isempty(current)
        loop = WorkgroupLoop(deepcopy(indicies), current, allocations, deepcopy(private))
        push!(new_stmts, emit(loop))
    end
    return new_stmts
end

function emit(loop)
    idx = gensym(:I)
    for stmt in loop.indicies
        # splice index into the i = @index(Cartesian, $idx)
        @assert stmt.head === :(=)
        rhs = stmt.args[2]
        push!(rhs.args, idx)
    end
    stmts = Any[]
    append!(stmts, loop.allocations)
    # don't emit empty loops
    if !(isempty(loop.stmts) || all(s->s isa LineNumberNode, loop.stmts))
        body = Expr(:block, loop.stmts...)
        body = postwalk(body) do expr
            if @capture(expr, A_[i__])
                if A in loop.private
                    return :($A[$(i...), $(idx).I...])
                end
            end
            return expr
        end
        loopexpr = quote
            for $idx in $__workitems_iterspace()
                $__validindex($idx) || continue
                $(loop.indicies...)
                $(unblock(body))
            end
        end
        push!(stmts, loopexpr)
    end

    return unblock(Expr(:block, stmts...))
end


text:
module IO

using StaticArrays
using Printf
using RecipesBase
using OrderedCollections
using Requires # for conditional loading of vtkIO

using WavePropBase
using WavePropBase.Utils
using WavePropBase.Geometry
using WavePropBase.Trees
using WavePropBase.Integration
using WavePropBase.Interpolation
using WavePropBase.Mesh

WavePropBase.@import_interface

include("plotsIO.jl")

function __init__()
    # if WriteVTK is available, include vtkIO
    @require WriteVTK="64499a7a-5c06-52f2-abe2-ccb03c286192" begin
        @info "including vtkIO.jl from WavePropBase/IO"
        include("vtkIO.jl")
    end
end

end # module

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase