text:
using SeriesIterators
using Test

# Test: Ramanujan's famous series to calculate 1/π
# https://en.wikipedia.org/wiki/Srinivasa_Ramanujan#Mathematical_achievements
ramagen(k) = big(1103.0) * 2 * sqrt(big(2)) / 9801, big(1.0)
function ramagen(k, (p,p1))
    k = oftype(p1, k)
    p1 *= (1-3/4k)*(1-2/4k)*(1-1/4k)/99^4
    p = p1 * (1103 + 26390*k) * 2 * sqrt(oftype(p1, 2)) / 9801
    p, p1
end

@testset "SeriesIterators.jl" begin

@testset "term iterator" begin
    setprecision(BigFloat, 256)
    t = iterm(ramagen, 0:5)
    @test sum(abs, (cumsum(t) .- big(1)/pi ) .* exp10.(8:8:48)) < 10
    @test inv(cumsum(t)[end]) ≈ pi
end

@testset "bernoulli" begin
    @test bernoulli(0) == 1
    @test bernoulli(1) == -1//2
    @test bernoulli(2) == 1//6
    @test bernoulli(3) == 0
    @test bernoulli.(4:2:20) == [
      -1//30
       1//42
      -1//30
       5//66
    -691//2730
       7//6
   -3617//510
   43867//798
   -174611//330]
    @test bernoulli(200) == -498384049428333414764928632140399662108495887457206674968055822617263669621523687568865802302210999132601412697613279391058654527145340515840099290478026350382802884371712359337984274122861159800280019110197888555893671151//1366530

end

@testset "bernoulli(Int128)" begin
    bernoulli128 = SeriesIterators.gen_bernoulli(Int128)
    @test bernoulli128(10) == 5 // 66
    @test bernoulli128(58)  == 84483613348880041862046775994036021 // 354
    @test_throws OverflowError bernoulli128(60)
end

@testset "bernoulli(Int16)" begin
    bernoulli16 = SeriesIterators.gen_bernoulli(Int16)
    @test bernoulli16(10) == 5 // 66
    @test bernoulli16(14) == 7 // 6
    @test_throws OverflowError bernoulli16(16)
end

@testset "cumulative sum iterator" begin
    itr = -10:2:10
    @test collect(icumsum(itr)) == cumsum(itr)
end

end


text:
export generate_JLD2

function generate_JLD2(fname,props_dics,Ts,Ps,fluid;include_std=true)

    pair_PT_table=[(P=P,T=T) for P in Ps, T in Ts]
    pair_PT=[pair_PT_table...]

    len_PT=length(pair_PT)
    jump_T=length(Ts)
    size_PT=size(pair_PT_table)

    props_keys=keys(props_dics)

    props_lines=Dict{Union{Symbol,String},Array{Float64}}(
        [k => zeros(Float64,len_PT) for k in props_keys]
    )

    props_tabs=Dict{Union{Symbol,String},Array{Float64}}(
        [k => zeros(Float64,size_PT) for k in props_keys]
    )

    for k in props_keys
        for (i, P) in enumerate(Ps)
            for (j, T) in enumerate(Ts)
                prop_calc = props_dics[k](P, T, fluid)

                props_tabs[k][i, j] = prop_calc

                props_lines[k][(i-1)*jump_T+j] = prop_calc

            end
        end
    end

    jldopen(fname,"w") do file
        file["fluid"]=fluid

        file["inputs/Ts"]=Ts
        file["inputs/Ps"]=Ps
        file["inputs/pairs"]=pair_PT
        file["inputs/table"]=pair_PT_table

        file["inputs/len_data"]=len_PT
        file["props"]=[k for k in props_keys]

        file["data/linear"]=props_lines
        file["data/tabular"]=props_tabs

    end

    return true

end

ext: jl


text:
using JLD
using BenchmarkTools
using Dates
using Base.Sys
using CellListsBenchmarks
using Plots

gettime(x) = getfield(x, :time)
gettime_nogc(x) = getfield(x, :time) - getfield(x, :gctime)

filepath_serial = "output/cell_list_constructor_serial_2021-04-09T08:28:43+0300/n1000-d2-r0.1.jld"
filepath_parallel = "output/cell_list_constructor_parallel_2021-04-09T08:28:47+0300/n1000-d2-r0.1.jld"

ls = load(filepath_serial)["benchmarks"]
lp = load(filepath_parallel)["benchmarks"]

d = getfield(ls[1], :d)
r = getfield(ls[1], :r)
ns = getfield.(ls, :n)
ts = getfield.(ls, :trials)
tp = getfield.(lp, :trials)

ms = [median.(t) for t in ts]
mp = [median.(t) for t in tp]

means = [mean(gettime.(v)) for v in ms]
meanp = [mean(gettime.(v)) for v in mp]

p1 = plot(legend=:topleft)
plot!(p1, ns, means, markershape=:circle, markersize=2, linestyle=:dash, label="serial")
plot!(p1, ns, meanp, markershape=:circle, markersize=2, linestyle=:dash, label="parallel")

p2 = plot(legend=:bottomright)
plot!(p2, ns, means./meanp, markershape=:circle, markersize=2, linestyle=:dash, label="serial/parallel")

directory = joinpath("figures", "constructor")
if !isdir(directory)
    mkpath(directory)
end

savefig(p1, joinpath(directory, "d$d-r$r-n$(ns[end]).svg"))
savefig(p2, joinpath(directory, "d$d-r$r-n$(ns[end])-ratio.svg"))

max_stars_repo_name: jaantollander/cell-lists-benchmarks

max_issues_repo_name: jaantollander/cell-lists-benchmarks

max_forks_repo_name: jaantollander/cell-lists-benchmarks


text:
using Revise
using Dice
using Dice: num_flips, num_nodes, to_dice_ir

@macroexpand @dice

code = @dice begin
    function uniform(b::Int, w::Int) # b is the bits for uniform, w is the bitwidth
        x = Vector(undef, b)
        for i = b:-1:1
            x[i] = flip(0.5)
        end
        return add_bits(DistInt(x), w - b)
    end
    uniform(2, 2)
end

# BDD analysis
bdd = compile(code)
num_flips(bdd)
num_nodes(bdd)
infer(code, :bdd)

# IR analysis
# to_dice_ir(code)
# has_dice_binary() && rundice(code)
# has_dice_binary() && infer(code, :ocaml)


text:
include("./util.jl")
import LinearAlgebra
using Graphs: SimpleGraph, is_connected, add_edge!, a_star

struct Orientation
    rotation_matrix::Array{Int,2}
    translation::Array{Int,1}
end

struct BeaconPosition
    rel_position::Array{Int,1}
end

struct Scanner
    index::Int
    detected_beacons::Array{BeaconPosition}
    abs_pos::Array{Int,1}
end

@enum Dim X Y Z

function rot(num_turns::Int, dim::Dim)::Array{Int,2}
    num_turns = num_turns % 4
    if num_turns == 0
        return [
            1 0 0
            0 1 0
            0 0 1
        ]
    end
    if dim == X
        if num_turns == 1
            [
                1 0 0
                0 0 -1
                0 1 0
            ]
        elseif num_turns == 2
            [
                1 0 0
                0 -1 0
                0 0 -1
            ]
        elseif num_turns == 3
            [
                1 0 0
                0 0 1
                0 -1 0
            ]
        end
    elseif dim == Y
        if num_turns == 1
            [
                0 0 1
                0 1 0
                -1 0 0
            ]
        elseif num_turns == 2
            [
                -1 0 0
                0 1 0
                0 0 -1
            ]
        elseif num_turns == 3
            [
                0 0 -1
                0 1 0
                1 0 0
            ]
        end
    elseif dim == Z
        if num_turns == 1
            [
                0 -1 0
                1 0 0
                0 0 1
            ]
        elseif num_turns == 2
            [
                -1 0 0
                0 -1 0
                0 0 1
            ]
        elseif num_turns == 3
            [
                0 1 0
                -1 0 0
                0 0 1
            ]
        end
    end
end

function apply(orient::Orientation, scanner::Scanner)::Scanner
    Scanner(
        scanner.index,
        map(bp -> BeaconPosition((orient.rotation_matrix * bp.rel_position) + orient.translation), scanner.detected_beacons),
        (orient.rotation_matrix * scanner.abs_pos) + orient.translation
    )
end

function apply_inv(orient::Orientation, scanner::Scanner)::Scanner
    Scanner(
        scanner.index,
        map(bp -> BeaconPosition(inv(orient.rotation_matrix) * (bp.rel_position - orient.translation)), scanner.detected_beacons),
        inv(orient.rotation_matrix) * (scanner.abs_pos - orient.translation)
    )
end

function parse_input()::Array{Scanner}
    lines = util.input_lines(19)
    scanners = []
    new_scanner_r = r"--- scanner (\d+) ---"
    triple_r = r"(-?\d+),(-?\d+),(-?\d+)"
    curr_scanner = -1
    curr_scanner_triples = []
    for line in lines
        m = match(new_scanner_r, line)
        if !isnothing(m)
            if curr_scanner != -1
                push!(scanners, Scanner(curr_scanner, curr_scanner_triples, [0, 0, 0]))
            end
            curr_scanner = parse(Int, m.captures[1]) + 1
            curr_scanner_triples = []
        end
        m = match(triple_r, line)
        if !isnothing(m)
            push!(
                curr_scanner_triples,
                BeaconPosition([
                    parse(Int, m.captures[1]),
                    parse(Int, m.captures[2]), parse(Int, m.captures[3])
                ])
            )
        end
    end
    push!(scanners, Scanner(curr_scanner, curr_scanner_triples, [0, 0, 0]))
    scanners
end

function check_align(s0::Scanner, s1::Scanner, orient::Orientation)::Bool
    bs0 = map(it -> it.rel_position, s0.detected_beacons)
    bs1 = map(it -> it.rel_position, apply(orient, s1).detected_beacons)
    length(intersect(bs0, bs1)) >= 12
end

function try_align(s0::Scanner, s1::Scanner)::Union{Orientation,Nothing}
    for x_turns = 0:3
        for y_turns = 0:3
            for z_turns = 0:3
                rot_mat = rot(x_turns, X) * rot(y_turns, Y) * rot(z_turns, Z)
                for b0 in s0.detected_beacons
                    for b1 in s1.detected_beacons
                        offset = b0.rel_position - (rot_mat * b1.rel_position)
                        orient = Orientation(rot_mat, offset)
                        if check_align(s0, s1, orient)
                            return orient
                        end
                    end
                end
            end
        end
    end
    nothing
end

function align(scanners::Array{Scanner})::Dict{Tuple{Int,Int},Orientation}
    alignments = Dict{Tuple{Int,Int},Orientation}()
    for s0 in scanners
        for s1 in scanners
            if s1.index <= s0.index
                continue
            end
            # println("Aligning $(s0.index) to $(s1.index)")
            alignment = try_align(s0, s1)
            if !isnothing(alignment)
                # println(alignment)
                alignments[(s0.index, s1.index)] = alignment
            end
        end
    end
    alignments
end

function transform_to_1(scanners::Array{Scanner}, alignment::Dict{Tuple{Int,Int},Orientation})::Array{Scanner}
    g = SimpleGraph(length(scanners))
    for (a, b) in keys(alignment)
        add_edge!(g, a, b)
    end
    @assert is_connected(g)
    output = []
    for scanner in scanners
        curr_transformed = scanner
        path_to_1 = a_star(g, scanner.index, 1)
        for e in path_to_1
            if haskey(alignment, (e.src, e.dst))
                curr_transformed = apply_inv(alignment[(e.src, e.dst)], curr_transformed)
            else
                curr_transformed = apply(alignment[(e.dst, e.src)], curr_transformed)
            end
        end
        push!(output, curr_transformed)
    end
    output
end

function extract_beacons(scanners::Array{Scanner})::Set{Array{Int,1}}
    Set([pos.rel_position for scan in scanners for pos in scan.detected_beacons])
end

function max_manhattan_distance(scanners::Array{Scanner})::Int
    best = 0
    for s0 in scanners
        for s1 in scanners
            best = max(best, LinearAlgebra.norm(s0.abs_pos - s1.abs_pos, 1))
        end
    end
    best
end

scanners = parse_input()
alignment = align(scanners)
transformed = transform_to_1(scanners, alignment)
beacons = extract_beacons(transformed)
println("Part 1: $(length(beacons))")
dist = max_manhattan_distance(transformed)
println("Part 2: $dist")


text:
import Base: getindex, haskey, rand, iterate

# ------------ Core ------------ #

Base.iterate(s::Symbol) = s

# Special calls recognized by tracer.
rand(addr::A, d::Distribution{T}) where {A <: Address, T} = error("(rand) call with address $addr evaluated outside of the tracer.\nThis normally occurs because you're not matching dispatch correctly, or you've forgotten to tell the tracer to recurse into a call site (wrap it with rand).")
rand(addr::A, fn::Function, args...) where A <: Address = error("(rand) call with address $addr evaluated outside of the tracer.\nThis normally occurs because you're not matching dispatch correctly, or you've forgotten to tell the tracer to recurse into a call site (wrap it with rand).")
rand(addr::A, fn::Function, args::Tuple, ret_score::Function) where A <: Address = error("(rand) call with address $addr evaluated outside of the tracer.\nThis normally occurs because you're not matching dispatch correctly, or you've forgotten to tell the tracer to recurse into a call site (wrap it with rand).")
learnable(addr::A) where {A <: Address} = error("(learnable) call with address $addr evaluated outside of the tracer.\nThis normally occurs because you're not matching dispatch correctly.")
plate(addr::A, args...) where A <: Address = error("(plate) call with address $addr evaluated outside of the tracer.\nThis normally occurs because you're not matching the dispatch correctly.")
markov(addr::A, args...) where A <: Address = error("(markov) call with address $addr evaluated outside of the tracer.\nThis normally occurs because you're not matching the dispatch correctly.")
cond(addr::A, args...) where A <: Address = error("(cond) call with address $addr evaluated outside of the tracer.\nThis normally occurs because you're not matching the dispatch correctly.")
factor(args...) = args

# Generic abstract types..
abstract type RecordSite end
abstract type CallSite <: RecordSite end
abstract type LearnableSite <: RecordSite end

# ------------ Choice sites - rand calls with distributions. ------------ #

struct ChoiceSite{T} <: RecordSite
    score::Float64
    val::T
end
get_score(cs::ChoiceSite) = cs.score
getindex(cs::ChoiceSite, addr::Address) = nothing
haskey(cs::ChoiceSite, addr::Address) = false
get_ret(cs::ChoiceSite) = cs.val

# ------------ Site with a user-declared parameter ------------ #

struct ParameterSite{T} <: LearnableSite
    val::T
end

# ------------ includes - traces + call sites ------------ #

abstract type Trace end

include("traces/hierarchical.jl")
include("traces/vectorized.jl")
include("traces/conditional.jl")

# ------------ Pretty printing ------------ #

function Base.display(call::C;
                      fields::Array{Symbol, 1} = [:val],
                      show_full = false) where C <: CallSite
    println("  __________________________________\n")
    println("               Playback\n")
    println(" type : $C\n")
    map(fieldnames(C)) do f
        val = getfield(call, f)
        typeof(val) <: Dict{Address, ChoiceSite} && begin
            vals = collect(val)
            if length(vals) > 5 && !show_full
                map(vals[1:5]) do (k, v)
                    println(" $(k)")
                    map(fieldnames(ChoiceSite)) do nm
                        !(nm in fields) && return
                        println("          $(nm)  = $(getfield(v, nm))")
                    end
                    println("")
                end
                println("                  ...\n")
                println("  __________________________________\n")
                return
            else
                map(vals) do (k, v)
                    println(" $(k)")
                    map(fieldnames(ChoiceSite)) do nm
                        !(nm in fields) && return
                        println("          $(nm)  = $(getfield(v, nm))")
                    end
                    println("")
                end
                println("  __________________________________\n")
                return
            end
        end
        typeof(val) <: Real && begin
            println(" $(f) : $(val)\n")
            return
        end
        println(" $(f) : $(typeof(val))\n")
    end
    println("  __________________________________\n")
end

import Base.collect
function collect(tr::Trace)
    addrs = Any[]
    chd = Dict{Any, Any}()
    meta = Dict()
    collect!(addrs, chd, tr, meta)
    return addrs, chd, meta
end

function Base.display(tr::Trace;
                      show_values = true,
                      show_types = false)
    println("  __________________________________\n")
    println("               Addresses\n")
    addrs, chd, meta = collect(tr)
    if show_values
        for a in addrs
            if haskey(meta, a)
                println("$(meta[a]) $(a) = $(chd[a])")
            else
                println("$(a) = $(chd[a])")
            end
        end
    elseif show_types
        for a in addrs
            println(" $(a) = $(typeof(chd[a]))")
        end
    elseif show_types && show_values
        for a in addrs
            println(" $(a) = $(chd[a]) : $(typeof(chd[a]))")
        end
    else
        for a in addrs
            println(" $(a)")
        end
    end
    println("  __________________________________\n")
end

# ------------ Documentation ------------ #

@doc(
"""
```julia
struct ChoiceSite{T} <: RecordSite
    score::Float64
    val::T
end
```
A record of a random sample at an addressed `rand` call with a `Distribution` instance. Keeps the value of the random sample and the `logpdf` score.
""", ChoiceSite)

@doc(
"""
```julia
abstract type CallSite <: RecordSite end
```
Abstract base type of all call sites.
""", CallSite)

@doc(
"""
```julia
abstract type Trace end
```
Abstract base type of all traces.
""", Trace)

@doc(
"""
```julia
struct HierarchicalTrace <: Trace
    calls::Dict{Address, CallSite}
    choices::Dict{Address, ChoiceSite}
    params::Dict{Address, LearnableSite}
    end
end
```
Execution trace which allows the tracking of randomness metadata in a function call.
""", HierarchicalTrace)

@doc(
"""
```julia
struct VectorizedTrace{C <: RecordSite} <: Trace
    subrecords::Vector{C}
    params::Dict{Address, LearnableSite}
end

```
Structured execution trace for `markov` and `plate` calls. The dependency structure interpretation of the `subrecords` vector depends on the call. For `markov`, the structure is Markovian. For `plate`, each element is drawn IID from the program or distribution provided to `plate`.
""", VectorizedTrace)

@doc(
"""
```julia
struct HierarchicalCallSite{J, K} <: CallSite
    trace::HierarchicalTrace
    score::Float64
    fn::Function
    args::J
    ret::K
end
```
A record of a black-box call (e.g. no special tracer language features). Records the `fn` and `args` for the call, as well as the `ret` return value.
""", HierarchicalCallSite)

@doc(
"""
```julia
struct VectorizedCallSite{F <: Function, C <: RecordSite, J, K} <: CallSite
    trace::VectorizedTrace{C}
    score::Float64
    fn::Function
    args::J
    ret::Vector{K}
end
```
A record of a call site using the special `plate` and `markov` tracer language features. Informs the tracer that the call conforms to a special pattern of randomness dependency, which allows the storing of `Trace` instances sequentially in a vector.
""", VectorizedCallSite)


text:
# ---
# title: 1095. Find in Mountain Array
# id: problem1095
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Binary Search
# link: <https://leetcode.com/problems/find-in-mountain-array/description/>
# hidden: true
# ---
#
# _(This problem is an   **interactive problem**.)_
#
# You may recall that an array `A` is a _mountain array_ if and only if:
#
#   * `A.length >= 3`
#   * There exists some `i` with `0 < i < A.length - 1` such that:
#     * `A[0] < A[1] < ... A[i-1] < A[i]`
#     * `A[i] > A[i+1] > ... > A[A.length - 1]`
#
# Given a mountain array `mountainArr`, return the **minimum**  `index` such
# that `mountainArr.get(index) == target`.  If such an `index` doesn't exist,
# return `-1`.
#
# **You can 't access the mountain array directly.**  You may only access the
# array using a `MountainArray` interface:
#
#   * `MountainArray.get(k)` returns the element of the array at index `k` (0-indexed).
#   * `MountainArray.length()` returns the length of the array.
#
# Submissions making more than `100` calls to `MountainArray.get` will be judged
# _Wrong Answer_.  Also, any solutions that attempt to circumvent the judge will
# result in disqualification.
#
#
#
# **Example 1:**
#
#
#
#     Input: array = [1,2,3,4,5,3,1], target = 3
#     Output: 2
#     Explanation: 3 exists in the array, at index=2 and index=5. Return the minimum index, which is 2.
#
# **Example 2:**
#
#
#
#     Input: array = [0,1,2,4,2,1], target = 3
#     Output: -1
#     Explanation: 3 does not exist in the array, so we return -1.
#
#
#
#
# **Constraints:**
#
#   * `3 <= mountain_arr.length() <= 10000`
#   * `0 <= target <= 10^9`
#   * `0 <= mountain_arr.get(index) <= 10^9`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module MuxMiddleware

using Mux

export hostbranch, regexbranch

"""
    hostbranch(host::AbstractString, app...)

Branches to `app...` when the request's hostname matches `host` exactly.
"""
hostbranch(host::AbstractString, app...) = branch(app...) do req
    req_host_all = first(filter(hdr->hdr[1]=="Host", req[:headers]))[2]
    # FIXME: Strip protocol?
    req_host = first(split(req_host_all, ":"))
    return host == req_host
end

"""
    regexbranch(regex::Regex, app...)

Branches to `app...` when `regex` matches the request path. The path will
always start with a '/'. If the regex matches, the RegexMatch result will be
stored in the request's `:regex_match` field.
"""
regexbranch(regex::Regex, app...) = branch(app...) do req
    matches = match(regex, "/" * join(req[:path], "/"))
    matches === nothing && return false
    req[:regex_match] = matches
    return true
end

end # module


text:
# # Supervised 3D Mesh Reconstruction
#
# ## Problem Description:
# We are given an initial source shape (sphere in this case) and we want to deform
# this source shape to fit target shape (dolphin in this case). For this
# demonstration, we will be using Triangle Mesh for the representation of source
# and target shape.
#
# Triangle Mesh has two main components, vertices and faces. Deformation of source
# shape to fit target shape can be achieved by offsetting source's vertices to fit
# target surface. Also, the number of vertices and faces is not equal in source
# and target shape.
#
# !!!note
#     For visualization purpose we will require to install Makie and
#     compatible backend (GLMakie or WGLMakie). To install it simply run
#     `] add AbstractPlotting GLMakie` in the julia prompt.

using Flux3D, Zygote, Flux, FileIO, Statistics
using Makie, GLMakie, WGLMakie, GeometricalPredicates, GeometryBasics

GLMakie.activate!()
set_theme!(show_axis = false)

# ## Downloading obj file of sphere and dolphin

download(
    "https://github.com/nirmalsuthar/public_files/raw/master/dolphin.obj",
    "dolphin.obj",
)
download("https://github.com/nirmalsuthar/public_files/raw/master/sphere.obj", "sphere.obj")

# ## Loading Triangle Mesh
# Triangle Mesh is handled by TriMesh in Flux3D. TriMesh also supports batched
# format, namely padded, packed and list, which allow us to use fast batched
# operations. We can load TriMesh with load_trimesh function
# (supports `obj`, `stl`, `ply`, `off` and `2DM`)

dolphin = load_trimesh("dolphin.obj")
src = load_trimesh("sphere.obj")

# ## Preprocessing data
# Preprocessing tgt (dolphin), such that its mean is zero and also scale it
# according to the bounding box of src (sphere), So that src can converge at
# greater speed.

tgt = deepcopy(dolphin)
verts = get_verts_packed(tgt)
center = mean(verts, dims = 2)
verts = verts .- center
scale = maximum(abs.(verts))
verts = verts ./ scale
tgt._verts_packed = verts

# ## Visualizing TriMesh
# We will use `visualize` function for visualizing TriMesh. This function uses
# **Makie** for plotting. In fact, we can also visualize PointCloud using this
# function, which makes this function handy dealing with different 3D format.
vbox(visualize(src), visualize(tgt))

# ```@raw html
# <p align="center">
#     <img width=480 height=270 src="../../assets/fitmesh_initial.png">
# </p>
# ```

# ## Defining loss objective
# Starting from the src mesh, we will deform src mesh by offsetting its vertices
# (by offset array), such that new deformed mesh is close to target mesh.
# Therefore, our loss function will optimize the offset array. We will be using
# the following metrics to define loss objective:
#
# * `chamfer_distance` - the distance between the deformed mesh and target mesh, which is calculated by taking randomly 5000 points from the surface of each mesh and calculating chamfer_distance between these two pointcloud.
# * `laplacian_loss` - also known as Laplacian smoothing will act as a regularizer.
# * `edge_loss` - this will minimize edges length in deformed mesh, also act as a regularizer.

function loss_dolphin(x::AbstractArray, src::TriMesh, tgt::TriMesh)
    src = Flux3D.offset(src, x)
    loss1 = chamfer_distance(src, tgt, 5000)
    loss2 = laplacian_loss(src)
    loss3 = edge_loss(src)
    return loss1 + 0.1 * loss2 + loss3
end
# ## Defining learning rate and optimizer

lr = 1.0
opt = Flux.Optimise.Momentum(lr, 0.9)
_offset = zeros(Float32, size(get_verts_packed(src))...)

# ## Using GPU for fast training [**Optional**]
# We can convert the TriMesh structure to GPU or CPU using`gpu` and `cpu`
# function which is exactly the same syntax as Flux.

tgt = tgt |> gpu
src = src |> gpu
_offset = _offset |> gpu

# ## Optimizing the offset array
# We first initialize offset array as zeros, hence deformed mesh is equivalent to
# src mesh (sphere). Next, we calculate loss using this offset array and we
# compute derivatives wrt. offset array and finally optimize the array.

@info("Training...")
θ = Zygote.Params([_offset])
for itr = 1:2001
    gs = gradient(θ) do
        loss_dolphin(_offset, src, tgt)
    end
    Flux.update!(opt, _offset, gs[_offset])
    if (itr % 10 == 1)
        loss = loss_dolphin(_offset, src, tgt)
        @show itr, loss
        save("src_$(itr).png", visualize(Flux3D.offset(src, _offset)))
    end
end

anim = @animate for i ∈ 1:8
    Plots.plot(load("src_$(1+250*(i-1)).png"), showaxis = false)
end
gif(anim, "src_deform.gif", fps = 2)

# ```@raw html
# <p align="center">

#     <img width=256 height=256 src="../../assets/fitmesh_anim.gif">
# </p>
# ```
# ## Postprocessing the predicted mesh
# We create a new TriMesh by offsetting src by final offset array and scale up
# the final_mesh by the same scaling factor we scale down tgt, such that
# final_mesh has similar bounding box as dolphin mesh.

final_mesh = Flux3D.offset(src, _offset)
final_mesh = Flux3D.scale!(final_mesh, scale)

# ## Saving the final_mesh
# Flux3D provide IO function `save_trimesh`  to save TriMesh
# (supports `obj`, `stl`, `ply`, `off` and `2DM`)

save_trimesh("results/final_mesh.off", final_mesh)
save("results/final_mesh.png", visualize(final_mesh))
Flux3D.AbstractPlotting.vbox(visualize(final_mesh), visualize(dolphin))

# ```@raw html
# <p align="center">
#     <img width=480 height=270 src="../../assets/fitmesh_final.png">
# </p>
# ```

# ## Finally..
# * Look into the other examples in `examples/`
# * Read more the TriMesh in TriMesh section in documentation and
# Metrics/Transforms section for manipulating TriMesh and computing standard
# metrics.


text:
using ENN.Neurons, GLMakie

neuron1 = Neuron((((((),[:E1,:E2,:E3,:E4],:E),((),[:E1,:E2],:F),((),[:G1,:G2],:G),((),[:H1,:H2],:H)),[:E1,:E2],:B),((((),Symbol[:E1],:I),),Symbol[],:C),(((((((),[:K1,:K2],:K),((),[:L1,:L2],:L),((),[:A1,:A2],:M))),Symbol[],:J),),Symbol[:H1],:D)),[:A1,:A2],:A)
neuron2 = Neuron((((((),[:E1,:E2,:E3,:E4],:E),((),[:E1,:E2],:F),((),[:G1,:G2],:G),((),[:H1,:H2],:H)),[:E1,:E2],:B),((((),Symbol[:E1],:I),),Symbol[],:C),(((((((),[:K1,:K2],:K),((),[:L1,:L2],:L),((),[:A1,:A2],:M))),Symbol[],:J),),Symbol[:H1],:D)),[:A1,:A2],:A)

net = NeuralNetwork(
    Dict(:n1=>neuron1,:n2=>neuron2),
    Dict(:in1=>[(:n2, :E2, :excitatory)]),
    Dict(:n2=>:out1),
    Dict(:n1=>[(:n2, :E1, :excitatory)],:n2=>[(:n1, :E1, :excitatory),(:n1, :H1, :excitatory)])
)

f = Figure()
ax = Axis(f[1, 1], aspect=DataAspect())
res=netplot!(ax, net, linewidth=2)
display(f)


text:
# make detection more robust
# 2. sub-grid detection based on fitting, rather than accurate zero finding
# rationale: core location isn't a good proxy for vortex location.
# vortex is a bulk flow excitation, whereas any small noise will shift zero - without shifting bulk

#---
using Plots, Revise, VortexDistributions

# unitary?
v0 = [.2 .4 1]
w0 = PointVortex(v0)
@test rawData(w0) == v0

# @test foundNear(1)

# no vortices
Nx = 400; Ny = 400
Lx = 200;
Ly = 200;
x = LinRange(-Lx / 2, Ly / 2, Nx)
y = LinRange(-Ly / 2, Ly / 2, Ny)

psi0 = one.(x * y') |> complex;
psi = Torus(psi0, x, y);
vortfound = findvortices(psi)
@test vortfound == nothing

#--- debugging functions
import VortexDistributions: findvortices, findvortices_grid, findvortices_interp, findvortices_jumps

function findvortices(psi::Field,interp=true)
    if interp
        return findvortices_interp(psi)
    else
        return findvortices_grid(psi)
    end
end

function findvortices_grid(psi::Torus;shift=true)
    vort = findvortices_jumps(psi,shift=shift)
    rawvort = rawData(vort)
    vort = sortslices(rawvort,dims=1)
    return PointVortex(vort)
end

function findvortices_jumps(psi::Field;shift=true)
    @unpack x,y,ψ = psi
    phase = angle.(ψ)

    # Compute all nearest neighbour jumps
    diffx = phasejumps(phase,1); diffy = phasejumps(phase,2)

    # Compute all plaquette loop integrals with in-place memory recycling
    circshift!(phase,diffx,(0,1))
    diffx .-= phase; diffx .-= diffy
    circshift!(phase,diffy,(1,0))
    diffx .+= phase

    # Find all windings
    ixp,iyp,vp = findwhere(diffx .> 0.0)
    xp = x[ixp]; yp = y[iyp]; np = length(vp)
    ixn,iyn,vn = findwhere(diffx .< 0.0)
    xn = x[ixn]; yn = y[iyn]; nn = length(vn)

    if shift
        dx = x[2]-x[1]; dy = y[2] - y[1]
        xp .-= dx/2; yp .-= dy/2; xn .-= dx/2; yn .-= dy/2
    end

    vortices = [xn yn -vn; xp yp vp]

    return PointVortex(vortices)
end

function findvortices_interp(psi::Field)
    vort = findvortices_grid(psi,shift=true)

        return vort
    end
end

# function findvortices_interp(psi::Field)
#     vort = findvortices_grid(psi,shift=true)
#     vort = removeedgevortices(vort,psi)
#     #TODO: allow for interp with periodic data (here edges are stripped)
#
#     for (j,vortex) in enumerate(vort)
#         try
#         vortz,psiz = corezoom(vortex,psi)
#         vortz,psiz = corezoom(vortz,psiz)
#         vortz,psiz = corezoom(vortz,psiz)
#         vort[j] = vortz
#         catch nothing
#         end
#         return vort
#     end
# end

#--- one vortex
Nx = 400; Ny = 400
Lx = 200;
Ly = 200;
x = LinRange(-Lx / 2, Ly / 2, Nx)
y = LinRange(-Ly / 2, Ly / 2, Ny)

psi0 = one.(x * y') |> complex;
psi = Torus(psi0, x, y);
vort = randVortex(1, psi)
vortex!(psi, vort)

vort = PointVortex(vort)

vortfound = findvortices(psi)

vfdata = rawData(vortfound)

vdata = rawData(vort)

heatmap(x,y,abs2.(psi.ψ),aspect_ratio=1)
heatmap(x,y,angle.(psi.ψ),aspect_ratio=1)

# Let's try fitting a Gaussian to the density and compare with root finding
# Do this with and without local grid noise

#--- local field (within 2 healing lengths)
x1 = x[@. abs(x -vort[1].xv)<2]
y1 = y[@. abs(x -vort[1].yv)<2]

psi1 = psi.ψ[(@. abs(x -vort[1].xv)<2) , @. abs(x -vort[1].yv)<2]

heatmap(x1,y1,abs2.(psi1),aspect_ratio=1)
heatmap(x1,y1,angle.(psi1),aspect_ratio=1)

m = abs2.(psi1)
m = maximum(m) .- m
heatmap(x1,y1,m,aspect_ratio=1)

using DataFitting

gaussian(x, y, p1, p2, p3,p4) = @.  exp(-((x-p1)^2  +  (y-p2)^2)/p3^2)*p4

dom = CartesianDomain(x1,y1)

data = Measures(m) #Measures(d + randn(rng, size(d)), 1.)

model1 = Model(:comp1 => FuncWrap(gaussian, vort[1].xv, vort[1].yv, 1.,1.))
prepare!(model1, dom, :comp1)
@time result = fit!(model1, data)

#--- add local noise
phi = Torus(copy(psi.ψ),x,y)
phi.ψ .+= (randn()+im*randn())*0.2

vort = findvortices_interp(phi)
x1 = x[@. abs(x -vort[1].xv)<2]
y1 = y[@. abs(x -vort[1].yv)<2]

phi1 = phi.ψ[(@. abs(x -vort[1].xv)<2) , @. abs(x -vort[1].yv)<2]

heatmap(x1,y1,abs2.(phi1),aspect_ratio=1)
heatmap(x1,y1,angle.(phi1),aspect_ratio=1)
heatmap(x,y,angle.(phi.ψ),aspect_ratio=1)

m = abs2.(phi1)
m = maximum(m) .- m
heatmap(x1,y1,m,aspect_ratio=1)

using DataFitting

gaussian(x, y, p1, p2, p3,p4) = @.  exp(-((x-p1)^2  +  (y-p2)^2)/p3^2)*p4

dom = CartesianDomain(x1,y1)

data = Measures(m) #Measures(d + randn(rng, size(d)), 1.)

model1 = Model(:comp1 => FuncWrap(gaussian, vort[1].xv, vort[1].yv, 1.,1.))
prepare!(model1, dom, :comp1)
@time result = fit!(model1, data)

vfoundnoise = findvortices(phi) |> rawData


text:
export viz_policy

function viz_policy(alpha::Matrix{Float64}, neuralNetworkPath::AbstractString, batch_size::Int64=500)

    nnet = NNet(neuralNetworkPath);
    grid  = RectangleGrid(pas,taus,sis,sos,psis,thetas,ranges)
    grid2 = RectangleGrid(thetas,ranges)

    policy = read_policy(ACTIONS, alpha)
    c = RGB{U8}(1.,1.,1.) # white
    e = RGB{U8}(.0,.0,.5) # pink
    a = RGB{U8}(.0,.600,.0) # green
    d = RGB{U8}(.5,.5,.5) # grey
    b = RGB{U8}(.7,.9,.0) # neon green
    colors = [a; b; c; d; e]

    @manipulate for psi_int  = convert(Array{Int32,1},round(rad2deg(psis))),
        v_own = sos,
        v_int = sis,
        tau in taus,
        pa = pas,
        zoom = [4, 3, 2, 1],
        nbin = [100,150,200,250],
        Interp = [false,true],
        Belief = [false,true],
        beliefProb = [0.333,0.111,0.01],

        deltaR   = [40, 400, 4000, 0],
        deltaTh = [5.0, 30, 60, 0],
        deltaPsi = [5, 30, 60,90, 0],
        deltaV   = [10, 100, 200, 0],
        worst    = [false,true]

        deltaTheta = deltaTh*pi/180.0
        #Load table with the inputs needed to plot the heat map
        if Belief
            numBelief = 9
            inputsNet= zeros(nbin*nbin*numBelief,7)
            ind = 1
            for i=linspace(round(Int,-1*RANGEMAX/zoom),round(Int,RANGEMAX/zoom),nbin)
                for j=linspace(round(Int,-1*RANGEMAX/zoom),round(Int,RANGEMAX/zoom),nbin)
                    r = sqrt(i^2+j^2)
                    th = atan2(j,i)
                    bel = belief_states(r,th,psi_int,v_own,v_int,tau,pa,deltaR,deltaPsi,deltaV,deltaTheta,true)
                    inputsNet[ind:ind+8,:] = bel
                    ind = ind+numBelief
                end
            end
        else
            numBelief = 1
            inputsNet= zeros(nbin*nbin,7)
            ind = 1
            for i=linspace(round(Int,-1*RANGEMAX/zoom),round(Int,RANGEMAX/zoom),nbin)
                for j=linspace(round(Int,-1*RANGEMAX/zoom),round(Int,RANGEMAX/zoom),nbin)
                    r = sqrt(i^2+j^2)
                    th = atan2(j,i)
                    inputsNet[ind,:] = [r,th,deg2rad(psi_int),v_own,v_int,tau,pasTrue[pa]];
                    ind = ind+1
                end
            end
        end

        #Calculate all of the Q values from the input array
        q_nnet = zeros(nbin*nbin*numBelief,5);
        ind = 1

        while ind+batch_size<nbin*nbin*numBelief
            input = inputsNet[ind:(ind+batch_size-1),:]'
            output = evaluate_network_multiple(nnet,input)
            q_nnet = [q_nnet[1:(ind-1),:];output';q_nnet[ind+batch_size:end,:]]
            ind=ind+batch_size
        end
        input = inputsNet[ind:end,:]'
        output = evaluate_network_multiple(nnet,input)
        q_nnet = [q_nnet[1:(ind-1),:];output']

        ind = 1
        # Q Table Heat Map
        function get_heat1(x::Float64, y::Float64)
            r = sqrt(x^2+y^2)
            th = atan2(y,x)
            bel = belief_states(r,th,psi_int,v_own,v_int,tau,pa,deltaR,deltaPsi,deltaV,deltaTheta,false)
            qvals = evaluate(policy, get_belief(bel[1,end:-1:1][:],grid,Interp))
            if Belief
                if !worst
                    qvals*=beliefProb
                end
                for i=2:9
                    temp = evaluate(policy, get_belief(bel[i,end:-1:1][:],grid,Interp))
                    if worst
                        if minimum(temp)>minimum(qvals)
                            qvals = temp
                        end
                    else
                        qvals += temp*(1.0-beliefProb)/(numBelief-1.0)

                    end
                end
            end
            return rad2deg(ACTIONS[indmin(qvals)])
       end # function get_heat1

        #Neural Net Heat Map
       function get_heat2(x::Float64, y::Float64)
           r = sqrt(x^2+y^2)
           th = atan2(y,x)
           qvals = q_nnet[ind,:]
            if !worst
                qvals*=beliefProb
            end
           if Belief
               for i = 1:8
                    qvalTemp = q_nnet[ind+i,:]
                    if worst
                        if minimum(qvalTemp)>minimum(qvals)
                            qvals = qvalTemp
                        end
                    else
                        qvals+=qvalTemp*(1.0-beliefProb)/(numBelief-1.0)
                    end
               end
           end

           ind +=numBelief
           return rad2deg(ACTIONS[indmin(qvals)])
       end # function get_heat2

        if Belief
            g = GroupPlot(2, 2, groupStyle = "horizontal sep=3cm, vertical sep = 3cm")
            Belief = false
            push!(g, Axis([
                Plots.Image(get_heat1, (round(Int,-1*RANGEMAX/zoom), round(Int,RANGEMAX/zoom)),
                           (round(Int,-1*RANGEMAX/zoom), round(Int,RANGEMAX/zoom)),
                           zmin = -3, zmax = 3,
                           xbins = nbin, ybins = nbin,
                colormap = ColorMaps.RGBArray(colors), colorbar=false),
               Plots.Node(L">", 0, 0, style=string("font=\\Huge")),
               Plots.Node(L">", 53000/zoom, 53000/zoom, style=string("rotate=", psi_int, ",font=\\Huge"))
                    ], width="10cm", height="10cm", xlabel="x (ft)", ylabel="y (ft)", title="Nominal Q Table action"))

           push!(g, Axis([
               Plots.Image(get_heat2, (round(Int,-1*RANGEMAX/zoom), round(Int,RANGEMAX/zoom)),
                           (round(Int,-1*RANGEMAX/zoom), round(Int,RANGEMAX/zoom)),
                           zmin = -3, zmax = 3,
                           xbins = nbin, ybins = nbin,
                           colormap = ColorMaps.RGBArray(colors)),
               Plots.Node(L">", 0, 0, style=string("font=\\Huge")),
               Plots.Node(L">", 53000/zoom, 53000/zoom, style=string("rotate=", psi_int, ",font=\\Huge"))
                    ], width="10cm", height="10cm", xlabel="x (ft)", ylabel="y (ft)", title="Nominal Neural Net action"))

            Belief = true
            ind = 1
            push!(g, Axis([
                Plots.Image(get_heat1, (round(Int,-1*RANGEMAX/zoom), round(Int,RANGEMAX/zoom)),
                           (round(Int,-1*RANGEMAX/zoom), round(Int,RANGEMAX/zoom)),
                           zmin = -3, zmax = 3,
                           xbins = nbin, ybins = nbin,
                colormap = ColorMaps.RGBArray(colors), colorbar=false),
               Plots.Node(L">", 0, 0, style=string("font=\\Huge")),
               Plots.Node(L">", 53000/zoom, 53000/zoom, style=string("rotate=", psi_int, ",font=\\Huge"))
                    ], width="10cm", height="10cm", xlabel="x (ft)", ylabel="y (ft)", title="Belief Q Table action"))

            push!(g, Axis([
               Plots.Image(get_heat2, (round(Int,-1*RANGEMAX/zoom), round(Int,RANGEMAX/zoom)),
                           (round(Int,-1*RANGEMAX/zoom), round(Int,RANGEMAX/zoom)),
                           zmin = -3, zmax = 3,
                           xbins = nbin, ybins = nbin,
                           colormap = ColorMaps.RGBArray(colors)),
               Plots.Node(L">", 0, 0, style=string("font=\\Huge")),
               Plots.Node(L">", 53000/zoom, 53000/zoom, style=string("rotate=", psi_int, ",font=\\Huge"))
                    ], width="10cm", height="10cm", xlabel="x (ft)", ylabel="y (ft)", title="Belief Neural Net action"))

        else
            g = GroupPlot(2, 1, groupStyle = "horizontal sep=3cm")
            push!(g, Axis([
                Plots.Image(get_heat1, (round(Int,-1*RANGEMAX/zoom), round(Int,RANGEMAX/zoom)),
                           (round(Int,-1*RANGEMAX/zoom), round(Int,RANGEMAX/zoom)),
                           zmin = -3, zmax = 3,
                           xbins = nbin, ybins = nbin,
                colormap = ColorMaps.RGBArray(colors), colorbar=false),
               Plots.Node(L">", 0, 0, style=string("font=\\Huge")),
               Plots.Node(L">", 53000/zoom, 53000/zoom, style=string("rotate=", psi_int, ",font=\\Huge"))
                ], width="10cm", height="10cm", xlabel="x (ft)", ylabel="y (ft)", title="Q Table action"))

            push!(g, Axis([
               Plots.Image(get_heat2, (round(Int,-1*RANGEMAX/zoom), round(Int,RANGEMAX/zoom)),
                           (round(Int,-1*RANGEMAX/zoom), round(Int,RANGEMAX/zoom)),
                           zmin = -3, zmax = 3,
                           xbins = nbin, ybins = nbin,
                           colormap = ColorMaps.RGBArray(colors)),
               Plots.Node(L">", 0, 0, style=string("font=\\Huge")),
               Plots.Node(L">", 53000/zoom, 53000/zoom, style=string("rotate=", psi_int, ",font=\\Huge"))
                ], width="10cm", height="10cm", xlabel="x (ft)", ylabel="y (ft)", title="Neural Net action"))
        end
        g
    end # for p_int, v0, v1, pa, ta
end # function viz_pairwise_policy

function viz_policy(neuralNetworkPath::AbstractString, batch_size=500)

    nnet = NNet(neuralNetworkPath);

    c = RGB{U8}(1.,1.,1.) # white
    e = RGB{U8}(.0,.0,.5) # pink
    a = RGB{U8}(.0,.600,.0) # green
    d = RGB{U8}(.5,.5,.5) # grey
    b = RGB{U8}(.7,.9,.0) # neon green
    colors =[a; b; c; d; e]

    @manipulate for psi_int  = convert(Array{Int32,1},round(rad2deg(psis))),
        v_own = sos,
        v_int = sis,
        tau in taus,
        pa = pas,
        zoom = [4, 3, 2, 1],
        nbin = [100,150,200,250],
        Belief = [false,true],
        beliefProb = [0.333,0.111,0.01],

        deltaR   = [40, 400, 4000, 0],
        deltaTh  = [5.0, 30, 60, 0],
        deltaPsi = [5, 30, 60,90, 0],
        deltaV   = [10, 100, 200, 0],
        worst    = [false,true]

        deltaTheta = deltaTh*pi/180.0

        #Load table with the inputs needed to plot the heat map
        if Belief
            numBelief = 9
            inputsNet= zeros(nbin*nbin*numBelief,7)
            ind = 1
            for i=linspace(round(Int,-1*RANGEMAX/zoom),round(Int,RANGEMAX/zoom),nbin)
                for j=linspace(round(Int,-1*RANGEMAX/zoom),round(Int,RANGEMAX/zoom),nbin)
                    r = sqrt(i^2+j^2)
                    th = atan2(j,i)
                    bel = belief_states(r,th,psi_int,v_own,v_int,tau,pa,deltaR,deltaPsi,deltaV,deltaTheta,true)
                    inputsNet[ind:ind+8,:] = bel
                    ind = ind+numBelief
                end
            end
        else
            numBelief = 1
            inputsNet= zeros(nbin*nbin,7)
            ind = 1
            for i=linspace(round(Int,-1*RANGEMAX/zoom),round(Int,RANGEMAX/zoom),nbin)
                for j=linspace(round(Int,-1*RANGEMAX/zoom),round(Int,RANGEMAX/zoom),nbin)
                    r = sqrt(i^2+j^2)
                    th = atan2(j,i)
                    inputsNet[ind,:] = [r,th,deg2rad(psi_int),v_own,v_int,tau,pasTrue[pa]];
                    ind = ind+1
                end
            end
        end

        #Calculate all of the Q values from the input array
        q_nnet = zeros(nbin*nbin*numBelief,5);
        ind = 1

        while ind+batch_size<nbin*nbin*numBelief
            input = inputsNet[ind:(ind+batch_size-1),:]'
            output = evaluate_network_multiple(nnet,input)
            q_nnet = [q_nnet[1:(ind-1),:];output';q_nnet[ind+batch_size:end,:]]
            ind=ind+batch_size
        end
        input = inputsNet[ind:end,:]'
        output = evaluate_network_multiple(nnet,input)
        q_nnet = [q_nnet[1:(ind-1),:];output']

        ind = 1
        #Neural Net Heat Map
       function get_heat2(x::Float64, y::Float64)
           r = sqrt(x^2+y^2)
           th = atan2(y,x)
           qvals = q_nnet[ind,:]
            if !worst
                qvals*=beliefProb
            end
           if Belief
               for i = 1:8
                    qvalTemp = q_nnet[ind+i,:]
                    if worst
                        if minimum(qvalTemp)>minimum(qvals)
                            qvals = qvalTemp
                        end
                    else
                        qvals+=qvalTemp*(1.0-beliefProb)/(numBelief-1.0)
                    end
               end
           end

           ind +=numBelief
           return rad2deg(ACTIONS[indmin(qvals)])
       end # function get_heat2

        if Belief
            g = GroupPlot(2, 1, groupStyle = "horizontal sep=3cm, vertical sep = 3cm")
            Belief = false

           push!(g, Axis([
               Plots.Image(get_heat2, (round(Int,-1*RANGEMAX/zoom), round(Int,RANGEMAX/zoom)),
                           (round(Int,-1*RANGEMAX/zoom), round(Int,RANGEMAX/zoom)),
                           zmin = -3, zmax = 3,
                           xbins = nbin, ybins = nbin,
                           colormap = ColorMaps.RGBArray(colors)),
               Plots.Node(L">", 0, 0, style=string("font=\\Huge")),
               Plots.Node(L">", 53000/zoom, 53000/zoom, style=string("rotate=", psi_int, ",font=\\Huge"))
                    ], width="10cm", height="10cm", xlabel="x (ft)", ylabel="y (ft)", title="Nominal Neural Net action"))

            Belief = true
            ind = 1
            push!(g, Axis([
               Plots.Image(get_heat2, (round(Int,-1*RANGEMAX/zoom), round(Int,RANGEMAX/zoom)),
                           (round(Int,-1*RANGEMAX/zoom), round(Int,RANGEMAX/zoom)),
                           zmin = -3, zmax = 3,
                           xbins = nbin, ybins = nbin,
                           colormap = ColorMaps.RGBArray(colors)),
               Plots.Node(L">", 0, 0, style=string("font=\\Huge")),
               Plots.Node(L">", 53000/zoom, 53000/zoom, style=string("rotate=", psi_int, ",font=\\Huge"))
                    ], width="10cm", height="10cm", xlabel="x (ft)", ylabel="y (ft)", title="Belief Neural Net action"))

        else
            g = GroupPlot(1, 1, groupStyle = "horizontal sep=3cm")
            push!(g, Axis([
               Plots.Image(get_heat2, (round(Int,-1*RANGEMAX/zoom), round(Int,RANGEMAX/zoom)),
                           (round(Int,-1*RANGEMAX/zoom), round(Int,RANGEMAX/zoom)),
                           zmin = -3, zmax = 3,
                           xbins = nbin, ybins = nbin,
                           colormap = ColorMaps.RGBArray(colors)),
               Plots.Node(L">", 0, 0, style=string("font=\\Huge")),
               Plots.Node(L">", 53000/zoom, 53000/zoom, style=string("rotate=", psi_int, ",font=\\Huge"))
                ], width="10cm", height="10cm", xlabel="x (ft)", ylabel="y (ft)", title="Neural Net action"))
        end
        g
    end # for p_int, v0, v1, pa, ta, etc
end # function viz_policy


text:
__precompile__()
module ComplexSystems
    using DataFrames
    using Distributions
    using GLM
    using Statistics
    using LibGEOS

    export distance,freeze,look_around,random_walk,walker_set,dla
    export neighboors,initialize_population,friends_opinion,avg_friends_opinion
    export convince,convince_2_sites,convince_3_sites,sznajd_one_site,sznajd_two_site,sznajd_three_site
    export stupid_bots, smart_bots, cambridge_analytica, gen_box, make_grid
    export inside_box, euler, heun, lv, lot3, my_integrator
    export evolve, paste_string, lsystem,mult_evol_lsystem,turtle_step,update_angle

    # An AbstractType with an empty property
    mutable struct box_2d
        squares::Polygon
        empty::Bool
    end

    # walker type, with a position and if it is a wanderer
    mutable struct walker
        x::Float64
        y::Float64
        loner::Bool
    end

    include("diffusion_limited_aggregation.jl")
    include("dynamical_system.jl")
    include("l_system.jl")
    include("box_counting.jl")
    include("CA.jl")
end # module


text:
module PowerMonitor

import battery_cli_jll: battery

Base.@kwdef mutable struct PowerStatus
    pow_stat::Symbol = :none
    perc::Float64 = 0
end

function status(;show_raw = false)
    pow_stat = nothing
    perc = 0
    iob = IOBuffer()
    good = battery() do battery_cmd
        success(pipeline(`$(battery_cmd)`, iob))
    end
    !good && return PowerStatus(:nobatteries, 100)
    lines = split(String(take!(iob)), "\n")
    st = lines[1] # just use first battery for now
    show_raw && @info lines
    pow_stat = parse_charge_status(st)
    perc = parse_perc(st)
    return PowerStatus(pow_stat, perc)
end

function parse_charge_status(st::AbstractString)::Union{Symbol,Nothing}
    occursin("Charging", st) && return :charging
    occursin("Full", st) && return :full
    occursin("Discharging", st) && return :discharging
    occursin("Unknown", st) && return :unknown
    error("Cannot find charge st in: $st")
end
function parse_perc(st::AbstractString)::Float64
    parse(Float64, split(split(st, ",")[2], "%")[1])
end

### AUTOMATION
const timers = Timer[]

function automate(src_map::Dict{Symbol, Function}; interval = 1)
    Timer((timer)->src_map[status().pow_stat](), 0, interval = interval)
end

function autoprecomp_notbattery()
    status().pow_stat == :nobatteries && error("This system reports that no batteries are present")
    t = automate(
        Dict(   :nobatteries => ()->ENV["JULIA_PKG_PRECOMPILE_AUTO"]=1,
                :charging =>    ()->ENV["JULIA_PKG_PRECOMPILE_AUTO"]=1,
                :full =>        ()->ENV["JULIA_PKG_PRECOMPILE_AUTO"]=1,
                :discharging => ()->ENV["JULIA_PKG_PRECOMPILE_AUTO"]=0,
                :unknown =>     ()->ENV["JULIA_PKG_PRECOMPILE_AUTO"]=0,
            );
            interval = 1
        )
    push!(timers, t)
    nothing
end

function stop_automation()
    close.(timers)
    empty!(timers)
    nothing
end

end # module


text:
import Optim.optimize

export frank_wolfe, FrankWolfeParams

type FrankWolfeParams<:OptParams
	maxiters::Int
	abstol::AbstractFloat
	reltol::AbstractFloat
	stepsize::StepSizeRule
end
stop(p::FrankWolfeParams, UB, LB) = (UB - LB < p.abstol) || ((UB - LB)/max(abs(LB), abs(UB)) < p.reltol) ? true : false
FrankWolfeParams(;maxiters=50, abstol=1e-4, reltol=1e-2, stepsize=DecreasingStepSize(2,1)) =
	FrankWolfeParams(maxiters, abstol, reltol, stepsize)

function frank_wolfe(x::Array{Float64,2}, # starting point
	objective, grad_objective,
	constraint, # evaluates constraint
	delta::AbstractFloat, # bound on constraint function
	min_lin_st_constraint, # solves min_{c(x)<=delta} g \dot x (as a function of delta and g)
	params::FrankWolfeParams = FrankWolfeParams();
	verbose = false,
	B = -Inf) # stopping condition

	# initialize
	objval = objective(x)
	if verbose
		@printf("%10s%12s%12s\n", "iter", "obj", "constr")
		@printf("%10d%12.4e%12.4e\n", 0, objval, constraint(x))
	end

	for k=1:params.maxiters
		g = grad_objective(x)

		# step 1
		# tilde_x is solution to min_{c(theta)<=delta} grad(theta_old) \dot theta
		tilde_x = min_lin_st_constraint(g, delta)

		# step 2

		# x is solution to minimize_alpha o(a*x_old + (1-a)*tilde_x)
		# f(a) = objective((1-a)*x + a*tilde_x)
		# g!(a,g) = (g[1] = dot(grad_objective((1-a)*x + a*tilde_x), - x + tilde_x); g)
		# a = optimize(f, g!, [.5], method = :l_bfgs)x, x_old = , x

		# fixed stepsize
		a = step(params.stepsize)

		# bisection
		# f(a) = dot(grad_objective((1-a)*x + a*tilde_x), tilde_x - x)
		# a = find_zero(f, 0, 1/sqrt(k), tol=1e-2, maxiters=10)

		x_old = copy(x)
		x = (1-a)*x + a*tilde_x
		# the following line segfaults:
		# x, x_old = (1-a)*x + a*tilde_x, x
		objval, objval_old = objective(x), objval
		cx = constraint(x)
		if verbose && k%1==0
			@printf("%10d%12.4e%12.4e\n", k, objval, cx)
		end

		# stopping condition
		tilde_B = objval_old + dot(g, tilde_x - x_old)
		B = max(B, tilde_B)
		stop(params, objval, B) && break
	end
	return x
end

max_stars_repo_name: madeleineudell/FirstOrderOptimization