text:
#запрос данных из бокса
function getData(baseIP::IPv4,port::Union{String,Int64},param::Dict)
#baseURI = $localIP:$port
    strParams = []
    for key in keys(param)
        push!(strParams,String(key)*"="*param[key]*"&")
    end

    strParams = join(strParams)

    if !haskey(param,"to") && !haskey(param,"count") #если не было указано конца, то читаем всё
        strParams = strParams*"all&"
    end

    if !haskey(param,"from")
        strParams = strParams*"from=0&"
    end
    # println("http://$baseIP:$port/apibox/getData?$strParams")
    res = HTTP.request("GET", "http://$baseIP:$port/api/getData?$strParams")
    data = res.body |> base64decode
    index = param["index"]

    dataType = getDataType(baseIP,port,param["dataName"],param["index"]) #узнаем тип данных
    convertedData = unpack_vec(data, dataType)
    # @info convertedData
    return convertedData
end

#запрос данных из в интервале как бы
function getStructData(baseIP::IPv4,port::Union{String,Int64},dataName::String, fields::String,index="0",from="0",to="",count="")
    param = Dict()
    param["dataName"] = dataName
    param["index"] = index
    param["from"] = from
    param["to"] = to
    param["count"] = count
    param["fields"] = fields

    data = getStructData(baseIP, port, param)
    return data
end
function getStructData(baseIP::IPv4,port::Union{String,Int64},param::Dict)
#baseURI = $localIP:$port
    strParams = []
    for key in keys(param)
        push!(strParams,String(key)*"="*param[key]*"&")
    end

    strParams = join(strParams)

    if !haskey(param,"to") && !haskey(param,"count")
        strParams = strParams*"all&"
    end

    if !haskey(param,"from")
        strParams = strParams*"from=0&"
    end

    println("http://$baseIP:$port/api/getStructData?$strParams")
    res = HTTP.request("GET", "http://$baseIP:$port/api/getStructData?$strParams")
    data = res.body |> base64decode
    # @info data
    fields = param["fields"]
    allFields = (split(fields,",")...,)
    #alldataType = [] #собираем типы данных, чтобы работать с ними gvg: ЗАБУДЬ об этом синтаксисе: []
    alldataType = map(allFields) do f
        getDataType(baseIP,port,string(f))
    end
    # for f in allFields
    #     push!(alldataType, getDataType(baseIP,port,string(f))) #узнаем тип данных
    # end

    convertedData = unpack_vec(data, alldataType...)
    @info alldataType
    return convertedData
end

function getData(baseIP::IPv4,port::Union{String,Int64},dataName::String,index=0,from=0,to=[],count=[])

    addToRequest = ""
    if isempty(to) && isempty(count) #если не было указано конца, то читаем всё
        addToRequest = "&all"
    else
        if ~isempty(to)
            addToRequest = addToRequest*"&to=$to"
        end
        if ~isempty(count)
            addToRequest = addToRequest*"&count=$count"
        end
        if ~isempty(from)
            addToRequest = addToRequest*"&from=$from"
        end
    end
    res = HTTP.request("GET", "http://$baseIP:$port/api/getData?dataName=$dataName&index=$index$addToRequest")
    data = res.body |> base64decode

    dataType = getDataType(baseIP,port,dataName, index) #узнаем тип данных
    convertedData = unpack_vec(data, dataType)
    return convertedData
end

function getCountDataInInterval(baseIP::IPv4,port::Union{String,Int64}, param::Dict)
    if haskey(param,"index")
        ind = param["index"]
    else
        ind = 0
    end

    if haskey(param,"from")
        from = param["from"]
    else
        from = 0
    end

    if haskey(param,"to")
        to = param["to"]
    else
        to = []
    end
    if haskey(param,"count")
        count = param["count"]
    else
        count = []
    end

    countData = getCountDataInInterval(baseIP,port,param["dataName"],ind,from,to,count)
end
#определяем, сколько точек попало в заданный диапазон
function getCountDataInInterval(baseIP::IPv4,port::Union{String,Int64},dataName::String,index=0,from=0,to=[],count=[])
    addToRequest = ""
    if isempty(to) && isempty(count) #если не было указано конца, то читаем всё
        addToRequest = "&all"
    else
        if ~isempty(to)
            addToRequest = addToRequest*"&to=$to"
        end
        if ~isempty(count)
            addToRequest = addToRequest*"&count=$count"
        end
    end
    res = HTTP.request("GET", "http://$baseIP:$port/apibox/getCountDataInInterval?dataName=$dataName&index=$index$addToRequest")
    data = String(res.body)
    count = parse(Int32, data)
    if count<0 count=0 end #если данные не найдены, то возвращает 0

    # @info data
    # convertedData = reinterpret(Int32, base64decode(data)) |> collect
    return count

end
#запрос атрибутов из бокса
#преобразуются в словарь, так как записаны в формате JSON
function getAttr(baseIP::IPv4,port::Union{String,Int64},dataName::String)
#dataName = /Mark/Reo_0/breathsegpar
    res = HTTP.request("GET", "http://$baseIP:$port/api/getAttributes?dataName=$dataName")
    data = String(res.body)
    attrs = JSON.parse(data)
    return attrs["attrs"]
end

#получение имени данного
function getTag(baseIP::IPv4,port::Union{String,Int64},dataName::String,index=0)
    # println("http://$baseIP:$port/apibox/getData?getTag=$dataName&index=$index")
    res = HTTP.request("GET", "http://$baseIP:$port/apibox/getTag?dataName=$dataName&index=$index")
    tag = String(res.body)

    return  tag
end

#получение количества детей
function getChildsCount(baseIP::IPv4,port::Union{String,Int64},dataName::String,index::Union{String,Int64})
    # println("http://$baseIP:$port/apibox/getData?getTag=$dataName&index=$index")
    res = HTTP.request("GET", "http://$baseIP:$port/apibox/getChildsCount?getTag=$dataName&index=$index")
    count = String(res.body)

    return  count
end

#структура файла
function getDataTree(localIP::IPv4,port::Union{String,Int64})
    res = HTTP.request("GET", "http://$localIP:$port/api/getDataTree")
    # @info res
    out = String(res.body)
    parse(out)
    return tree
end

#определение типа данных на основании EntityInfo
function getDataType(baseIP::IPv4,port::Union{String,Int64},dataName::String,index=0)
    #проверяем, не составное ли имя
    dnparts = split(dataName,"_")
    if length(dnparts)>1
        index = dnparts[2]
        dataName =  dnparts[1]
    end
    res = HTTP.request("GET", "http://$baseIP:$port/apibox/getEntityInfo?dataName=$dataName&index=$index")
    info = res.body
    info = String(info)
    if occursin("type: float",info)
        dataType = Float32
    elseif occursin("type: long",info)
        dataType = Int32
    elseif occursin("type: BYTE",info)
        dataType = UInt8
    elseif occursin("type: WORD",info)
        dataType = UInt16
    elseif occursin("type: char",info)
        dataType = Int8
    elseif occursin("type: Time",info)
        dataType = DateTime
    elseif occursin("type: bin16",info)
        dataType = BitArray #но тут я не уверена gvg: Int16 ??
    elseif occursin("type: object",info)
        dataType = [] # лучше сразу выдать ошибку
    else
        println("Неизвестный тип "*info)
    end
    return dataType
end


text:
@cell "Figure and Subplots" begin
    fig, _ = scatter(randn(100, 2), color = :red)
    scatter(fig[1, 2], randn(100, 2), color = :blue)
    scatter(fig[2, 1:2], randn(100, 2), color = :green)
    scatter(fig[1:2, 3][1:2, 1], randn(100, 2), color = :black)
    scatter(fig[1:2, 3][3, 1], randn(100, 2), color = :gray)
    fig
end

@cell "Figure with Layoutables" begin
    fig = Figure(resolution = (900, 900))
    ax, sc = scatter(fig[1, 1][1, 1], randn(100, 2), axis = (;title = "Random Dots", xlabel = "Time"))
    sc2 = scatter!(ax, randn(100, 2) .+ 2, color = :red)
    ll = fig[1, 1][1, 2] = Legend(fig, [sc, sc2], ["Scatter", "Other"])
    lines(fig[2, 1:2][1, 3][1, 1], 0..3, sin ∘ exp, axis = (;title = "Exponential Sine"))
    heatmap(fig[2, 1:2][1, 1], randn(30, 30))
    heatmap(fig[2, 1:2][1, 2], randn(30, 30), colormap = :grays)
    lines!(fig[2, 1:2][1, 2], cumsum(rand(30)), color = :red, linewidth = 10)
    surface(fig[1, 2], collect(1.0:40), collect(1.0:40), (x, y) -> 10 * cos(x) * sin(y))
    fig[2, 1:2][2, :] = Colorbar(fig, vertical = false,
        height = 20, ticklabelalign = (:center, :top), flipaxis = false)
    fig[3, :] = Menu(fig, options = ["A", "B", "C"])
    lt = fig[0, :] = Label(fig, "Figure Demo")
    fig[5, :] = Textbox(fig)
    fig
end

ext: jl


text:
println("pid: $(getpid())"); flush(stdout)
include("sample/util.jl") # compress_simdat!

using BSON
using Cytof5

if length(ARGS) == 0
  SEED = 0
  RESULTS_DIR = "results/test/vb-sim-paper/$(SEED)/"
  SIMDAT_PATH = "../sim_study/simdata/kills-flowsom/N5000/K10/1/simdat.bson"
  K_VB = 30
  BATCHSIZE = 100
  NITERS = 5000
else
  SEED = parse(Int, ARGS[1])
  RESULTS_DIR = ARGS[2]
  SIMDAT_PATH = ARGS[3]
  K_VB = parse(Int, ARGS[4])
  BATCHSIZE = parse(Int, ARGS[5])
  NITERS = 20000
end
mkpath(RESULTS_DIR)

# For BSON
using Flux, Distributions

# Load sim data
# SIMDAT_PATH = "../sim_study/simdata/kills-flowsom/N500/K5/90/simdat.bson"
# SIMDAT_PATH = "../sim_study/simdata/kills-flowsom/N5000/K10/1/simdat.bson"
simdat = BSON.load(SIMDAT_PATH)[:simdat]
simdat[:y] = Vector{Matrix{Float64}}(simdat[:y])

# Generate model constnats
c = Cytof5.VB.Constants(y=simdat[:y], K=K_VB, L=Dict(false=>5, true=>5),
                        yQuantiles=[.0, .25, .5], pBounds=[.05, .8, .05],
                        use_stickbreak=false, tau=1e-5)
c.priors.eps = Beta(1, 99)
c.priors.sig2 = LogNormal(log(.3), .2)

println("seed: $SEED")
# Fit model
out = Cytof5.VB.fit(y=simdat[:y], niters=NITERS, batchsize=BATCHSIZE, c=c,
                    nsave=30, seed=SEED, flushOutput=true)

# Save results
out[:simdat] = compress_simdat!(simdat)
BSON.bson("$(RESULTS_DIR)/output.bson", out)

# Post process
include("post_process_defs_pyplot.jl")
out = nothing
post_process("$(RESULTS_DIR)/output.bson")


text:
# # Overlays

#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/overlays.ipynb)

# Overlays are useful to display additional information on top of a given driving scene.
# For example, they can be used to display the ID of the vehicle, their speed, or other useful information.

# More generally, renderable objects that are passed to the render function
# are rendered in the order they are passed in.
# For example, if we pass a collection consisting of a roadway and a scene,
# the roadway is rendered first, and the scene is rendered on top of the roadway.

# This behavior leaves the possibility to create overlays as well as "underlays"
# that may be drawn under the vehicles.

# In this tutorial, we will use the same stadium scenario as shown in the camera example.

using AutomotiveDrivingModels
using AutoViz
AutoViz.colortheme["background"] = colorant"white"; # hide
AutoViz.set_render_mode(:fancy)

nlanes = 2
nveh = 8
nticks = 100
timestep = 0.1

roadway = gen_stadium_roadway(nlanes)

scene = Scene([
    Entity(
        VehicleState(
            Frenet(roadway[LaneTag(1,rand(1:nlanes))], 10.0*i),  # position
            roadway, 4.0 + 2.0rand()  # speed
        ),
        VehicleDef(), i
    ) for i in 1:nveh
])

models = Dict((
    i => LatLonSeparableDriver(ProportionalLaneTracker(), IntelligentDriverModel())
    for i in 1:nveh
))
set_desired_speed!.(values(models), 8.0 .+ 4.0rand(nveh))

scenes = simulate(scene, roadway, models, nticks, timestep)
#md nothing # hide

# The following code will help us to animate the simulation using `Reel`

using Reel

function animate(roadway, scenes, overlays, camera=TargetFollowCamera(3, zoom=15.))
    animation = roll(fps=1.0/timestep, duration=nticks*timestep) do t, dt
        i = Int(floor(t/dt)) + 1
        update_camera!(camera, scenes[i])
        renderables = [
            roadway, scenes[i], overlays[i]...
        ]
        render(renderables, camera=camera)
    end
    return animation
end
nothing # hide

# The `animate` function defined above takes an `overlays` object which holds a list of
# overlays for every frame in `scenes`.

#md # !!! note
#     Note on compatibility:
#     The interface for rendering overlays has changed in version 0.8 of AutoViz.
#     In older versions of AutoViz, overlays were rendered by passing in
#     the current scene and roadway explicitly.
#     Since AutoViz 0.8, every renderable object must have the same interface.
#     As a result, an overlay that requires the `roadway` or `scene` for rendering,
#     must store those objects as internal variables.
#
#     (`AutoViz.RenderableOverlay` is an overlay that can be used as a wrapper around older overlays
#     which rely on the `scene` and `roadway` objects.
#     However, it is recommended that overlays implement the new interface for `add_renderable!`
#     since `AutoViz.RenderableOverlay` may be removed in future versions.)

# ## IDOverlay

# Show IDs for all vehicles in the scene
overlays = [[
    IDOverlay(scene=scene, x_off=-2, y_off=1)
] for scene in scenes]
animation = animate(roadway, scenes, overlays)
#md write("stadium_id_overlay.gif", animation) # hide
#md # ![stadium with ID overlay](stadium_id_overlay.gif)

# ## CarFollowingStatsOverlay

# Diplays info about the front neighbor of a car
# such as the difference in velocity and the relative distance.
# Show statistics for vehicle 3:

# ### TODO: overlay name is not very descriptive
overlays = [[CarFollowingStatsOverlay(
    scene=scene, roadway=roadway, target_id=3, font_size=20, color=colorant"black"
)] for scene in scenes]
animation = animate(roadway, scenes, overlays)
#md write("stadium_car_stats_overlay.gif", animation) # hide
#md # ![stadium with car stats overlay](stadium_car_stats_overlay.gif)

# ## LineToFrontOverlay

# Show the line to the front vehicle for vehicle 3
overlays = [[
    LineToFrontOverlay(scene=scene, roadway=roadway, target_id=3)
] for scene in scenes]
animation = animate(roadway, scenes, overlays)
#md write("stadium_line_front_overlay.gif", animation) # hide
#md # ![stadium with ID overlay](stadium_line_front_overlay.gif)

# ## TextOverlay

# The `TextOverlay` displays some text at the desired location.
# By default, the coordinates are in camera pixels.

overlays = [[TextOverlay(
    text=["Overlays are nice!", "frame $(i)/$(length(scenes))"],
    font_size=30, pos=VecE2(50.0, 100.0), color=colorant"black",
    coordinate_system=:camera_pixels
)] for i in 1:length(scenes)]
animation = animate(roadway, scenes, overlays)
#md write("stadium_text_overlay.gif", animation) # hide
#md # ![stadium with text overlay](stadium_text_overlay.gif)

# setting the `coordinate_system` keyword to `:scene` places the text
# in scene coordinates
overlays = [[TextOverlay(
    text=["Scene Coordinates!"],
    font_size=10+round(20cos(.01i)),
    pos=VecE2(2i, 40+round(20*sin(.002i))),
    color=colorant"black", coordinate_system=:scene
)] for i in 1:length(scenes)]
animation = animate(roadway, scenes, overlays)
#md write("stadium_moving_text_overlay.gif", animation) # hide
#md # ![stadium with moving text overlay](stadium_moving_text_overlay.gif)

# ## HistogramOverlay
#
# Display a bar at the specified position `pos`.
# The bar has dimensions `width` and `height`
# and is filled up to a given proportion of its height.
# The fill proportion is set using `val`, which should be a number between 0 and 1.
# If it is 0, the bar is not filled, if it is 1 it is filled to the top.
# The default units are in the camera frame.

max_speed = 14.0
overlays = [[HistogramOverlay(
    pos=VecE2(40.0, 100.0), val=get_by_id(scene, 1).state.v/max_speed, label="veh1 speed",
    fill_color=colorant"blue", line_color=colorant"black",
    font_size=24, width=10, height=50
)] for scene in scenes]
animation = animate(roadway, scenes, overlays)
#md write("stadium_histogram_overlay.gif", animation) # hide
#md # ![stadium with moving text overlay](stadium_histogram_overlay.gif)

# ## NeighborsOverlay

# Draws a line between a vehicle and its neighbors.

overlays = [[
    NeighborsOverlay(
        scene=scene, roadway=roadway, target_id=3,
        textparams=TextParams(color=colorant"black")
    )
] for (i, scene) in enumerate(scenes)]
animation = animate(roadway, scenes, overlays)
#md write("neighbors_overlay.gif", animation) # hide
#md # ![stadium with neighbors overlay](neighbors_overlay.gif)

# ## Custom Overlays

# ### Lane Highlighting

# Just like roadways and vehicles, overlays are renderable objects.
# In order to make a custom overlay type renderable, one needs to implement
# the `add_renderable!(::RenderModel, ::YourOverlay)` function.

# To define an overlay which highlights a lane, let's start by defining a
# custom type `LaneOverlay`

struct LaneOverlay
    lane::Lane
    roadway::Roadway
    color::Colorant
end

function AutoViz.add_renderable!(rendermodel::RenderModel, overlay::LaneOverlay)
    add_renderable!(rendermodel, overlay.lane, overlay.roadway, color_asphalt=overlay.color)
    return rendermodel
end

# We can now initiate a `LaneOverlay` object and pass it to the `render` function
lane_overlay = LaneOverlay(roadway[LaneTag(1,1)], roadway, RGBA(0.0,0.0,1.0,0.5))
camera = TargetFollowCamera(3, zoom=15.)
update_camera!(camera, scene)
snapshot = render([roadway, lane_overlay, scene], camera=camera)
#md write("stadium_lane_overlay.svg", snapshot)
#md # ![lane overlay](stadium_lane_overlay.svg)

# ### Concentric Rectangle Overlay

# In the following, we will show a more advanced example in which direct
# manipulation of the `CairoContext` is required.
# We will create an overlay which draws concentric rectangles around
# a vehicle with ID `target_id`

struct ConcentricRectOverlay
    target_id
    scene::Scene
    n::Int64
    color::Colorant
end

# We implement a new rendering function for concentric rectangles

using Cairo

function render_concentric_rects(
    ctx::CairoContext,
    n::Int,
    x::Real, y::Real,
    w::Real, l::Real,
    yaw::Real,
    color::Colorant
)
    save(ctx)
    translate(ctx, x, y)
    rotate(ctx, yaw)

    for i in 1:n
        AutoViz.render_round_rect(
            ctx,
            0, 0, l+2(i-1), w+2(i-1), 1., (i-.5),
            RGBA(color.r, color.g, color.b, 1-i/n),
            true, false,
        )
    end

    restore(ctx)
end

# The `add_renderable` function adds the function defined above to the
# instruction stack of the render model via `add_instruction!`.
# A series of common rendering instructions are already implemented as part
# of the `AutoViz.jl` framework.

function AutoViz.add_renderable!(rendermodel::RenderModel, overlay::ConcentricRectOverlay)
    veh = get_by_id(overlay.scene, overlay.target_id)
    x, y, yaw = posg(veh)
    w, l = AutomotiveDrivingModels.width(veh.def), length(veh.def)

    add_instruction!(
        rendermodel, render_concentric_rects,
        (overlay.n, x, y, w, l, yaw, overlay.color)
    )
    return rendermodel
end

concentric_rect_overlay = ConcentricRectOverlay(3, scene, 5, colorant"green")
snapshot = render([roadway, concentric_rect_overlay, scene], camera=camera)
#md write("stadium_concentric_rect_overlay.svg", snapshot)
#md # ![lane overlay](stadium_concentric_rect_overlay.svg)

nothing

max_stars_repo_name: tawheeler/AutoViz

max_issues_repo_name: tawheeler/AutoViz

max_forks_repo_name: tawheeler/AutoViz


text:
module VariableFrequencyOPF

#### Variable Frequency AC Optimal Power Flow ####

# Developed by David Sehloff (@dsehloff) and Line Roald (@lroald)

using PowerModels
using PowerModelsSecurityConstrained: parse_con_file
using Ipopt
using JuMP
using CSV
using DataFrames
using Statistics
using JSON
using Combinatorics
using CategoricalArrays
using Plots
using Plots.PlotMeasures
using StatsPlots
# using AmplNLWriter
import Unicode
# import GR
import MathOptInterface

include("plot/bar-plots.jl")
include("plot/line-plots.jl")
include("plot/multi_plot_new.jl")

include("core/multifrequency-opf.jl")
include("core/variables.jl")
include("core/constraints.jl")
include("core/data-import.jl")
include("core/start-vals.jl")

include("app/control-multi-folder.jl")
include("app/fixf-indirPQ.jl")
include("app/frequency-sweep.jl")
include("app/hvdc-multi-folder.jl")
include("app/modify-network.jl")
include("app/run-models.jl")
include("app/benchmarking.jl")

include("util/collect-results.jl")
include("util/dict-util.jl")
include("util/multi-folder.jl")
include("util/no-loss-override.jl")
include("util/param-sweep.jl")
include("util/time_series_opf.jl")
include("util/time_parse.jl")

end


text:
using MatrixPolynomials, LinearAlgebra
using Test

x = variable("x")

@testset "matrix operations" begin
    A = [x x-1;1 0]
    @test Matrix(I, size(A))*A == A
    @test A+A == 2*A
    @test A+A == [2x 2x-2;2 0]
    @test A*A == A^2
    @test A^2 == [x^2+x-1 x^2-x;x x-1]
    @test iszero(A-A)

    @test x*A == [x^2 x^2-x;x 0]
    @test A*x == x*A
    @test A+2 == [x+2 x+1;3 2]
    @test 2+A == A+2
    @test A-1 == [x-1 x-2;0 -1]
    @test 1-A == [1-x 2-x;0 1]

    B = [1 1;1 1]
    @test x+B == [x+1 x+1; x+1 x+1]
    @test B+x == x+B
    @test x^2*B == [x^2 x^2;x^2 x^2]
    @test B*x^2 == x^2*B

    @test deg(A) == 1
    @test deg(A^2) == 2
    @test lc(A) == [1 1;0 0]
    @test coeff(A, deg(A)) == lc(A)
    @test coeff(A, 0) == [0 -1;1 0]
    @test rev(A) == [1 -x+1;x 0]

    A1, B1 = toPencil(A)
    @test A1 == -coeff(A, 0)
    @test B1 == coeff(A, 1)

    A = [x^2 1 0;0 0 0;0 0 x]
    A1, B1 = toPencil(A)
    @test A1 == [0 0 0 0 -1 0; 0 0 0 0 0 0; 0 0 -1 0 0 0; 1 0 0 0 0 0; 0 1 0 0 0 0; 0 0 1 0 0 0]
    @test B1 == [1 0 0 0 0 0; 0 0 0 0 0 0; 0 0 0 0 0 0; 0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 1]
end


text:
# sample to show robot module's sensing
module RobotSensingSample
  # include external modules
  include(joinpath(split(@__FILE__, "src")[1], "src/model/world/puddle_world.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/model/agent/agent.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/model/robot/differential_wheeled_robot/differential_wheeled_robot.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/model/map/map.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/model/object/object.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/model/sensor/sensor.jl"))

  # main function
  # show world's graph
  function main(delta_time=0.1, end_time=30; is_test=false)
    # save path of output
    path = "src/model/sensor/robot_sensing_sample.gif"

    # world instance
    world = PuddleWorld(-5.0, 5.0, -5.0, 5.0,
                        delta_time=delta_time,
                        end_time=end_time,
                        is_test=is_test,
                        save_path=path)

    # map including landmarks
    map = Map()
    add_object(map, Object(2.0, -2.0, id=1))
    add_object(map, Object(-1.0, -3.0, id=2))
    add_object(map, Object(3.0, 3.0, id=3))
    append(world, map)

    # robot instance
    robot = DifferentialWheeledRobot([-2.0, -1.0, pi/5*6], 0.2, "black",
                                     Agent(0.2, 10.0/180*pi), delta_time,
                                     sensor=Sensor(map))
    append(world, robot)

    # draw graph
    draw(world)
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
using Cambrian
import YAML

cfg = get_config("test.yaml")

@testset "Random fitness" begin
    e = Evolution{FloatIndividual}(cfg)

    evaluate(e::Evolution) = random_evaluate(e)
    populate(e::Evolution) = oneplus_populate(e)

    @test length(e.population) == cfg.n_population
    for i in e.population
        @test all(i.fitness .== -Inf)
    end

    evaluate(e)
    fits = [i.fitness[1] for i in e.population]
    evaluate(e)
    # random fitness, all values should change
    for i in eachindex(e.population)
        @test fits[i] != e.population[i].fitness[1]
    end

    fits = copy([i.fitness[:] for i in e.population])
    step!(e)
    # random fitness, all values should change
    for i in eachindex(e.population)
        @test fits[i] != e.population[i].fitness[1]
    end
end


text:
using Logging
using TensorBoardLogger

# something like tf.keras.metrics.Mean
# a good reference https://github.com/apache/incubator-mxnet/blob/master/julia/src/metric.jl
mutable struct MeanMetric{T}
    # FIXME syntax for T
    sum::T
    n::Int
    # CAUTION for T, you must define Base.convert(::T, v::Int64)
    MeanMetric{T}() where {T} = new(0, 0)
end
function add!(m::MeanMetric{T}, v) where {T}
    m.sum += v
    m.n += 1
end
_get(m::MeanMetric) = if m.n != 0 m.sum / m.n else m.sum / 1 end
function get!(m::MeanMetric)
    res = _get(m)
    reset!(m)
    res
end
function reset!(m::MeanMetric{T}) where {T}
    m.sum = 0
    m.n = 0
end

struct GraphMetric
    # All float64, because when I average, the / operator, it will be float
    nnz::Float64
    nny::Float64
    tpr::Float64
    fpr::Float64
    fdr::Float64
    shd::Float64
    prec::Float64
    recall::Float64
    # GraphMetric() = new(0,0,0,0,0,0)
end

function show(io::IO, gm::GraphMetric)
    prec = round(gm.prec, digits=2)
    recall = round(gm.recall, digits=2)
    shd = round(gm.shd, digits=2)
    print(io, "GraphMetric: prec: $(prec), recall: $(recall), shd: $(shd)")
end

# FIXME this only custom @show and REPL, but does not affect logger like @info
# function Base.show(io::IO, ::MIME"text/plain", m::GraphMetric)
#     for fname in fieldnames(typeof(m))
#         println(io, "$fname = $(getfield(m, fname))")
#     end
# end

# from https://discourse.julialang.org/t/get-fieldnames-and-values-of-struct-as-namedtuple/8991
to_named_tuple(p) = (; (v=>getfield(p, v) for v in fieldnames(typeof(p)))...)

# I should have used "value type", but that seems to be less performant. I'm
# setting it to 0 for whatever v.
Base.convert(::Type{GraphMetric}, v::Int64) = GraphMetric(0,0,0,0,0,0,0,0)

import Base.+
function +(a::GraphMetric, b::GraphMetric)
    # FIXME can I automate this?
    GraphMetric(a.nnz + b.nnz,
                a.nny + b.nny,
                a.tpr + b.tpr,
                a.fpr + b.fpr,
                a.fdr + b.fdr,
                a.shd + b.shd,
                a.prec + b.prec,
                a.recall + b.recall)
end
import Base./
function /(a::GraphMetric, n::Int64)
    GraphMetric(a.nnz / n,
                a.nny / n,
                a.tpr / n,
                a.fpr / n,
                a.fdr / n,
                a.shd / n,
                a.prec / n,
                a.recall / n)
end

# FIXME overwrite TBLogger
import TensorBoardLogger.log_value
function log_value(logger::TBLogger, name::AbstractString, value::GraphMetric; step=nothing)
    # FIXME loop through all fields, and log value for each of them
    fn = propertynames(value)

    # for f=fn
    #     prop = getfield(value, f)
    #     log_value(logger, name*"/$f", prop, step=step)
    # end

    # FIXME it turns out the order is random, but at least the prefix works
    #
    # I want to control the order of the fields

    # log_value(logger, name*"/tpr", value.tpr, step=step)
    # log_value(logger, name*"/fpr", value.fpr, step=step)
    # log_value(logger, name*"/fdr", value.fdr, step=step)
    log_value(logger, name*"/prec", value.prec, step=step)
    log_value(logger, name*"/recall", value.recall, step=step)

    # tfboard has good support to log scale something near 0
    # log_value(logger, name*"/1-prec", 1-value.prec, step=step)
    # log_value(logger, name*"/1-recall", 1-value.recall, step=step)

    # use a different prefix
    log_value(logger, name*"/v/shd", value.shd, step=step)

    # FIXME I want to plot nnz and nny in the same plot
    # log_value(logger, name*"/v/nnz", value.nnz, step=step)
    # log_value(logger, name*"/v/nny", value.nny, step=step)
end

function test()
    m = MeanMetric()

    m = MeanMetric{Float64}()
    add!(m, 1.3)
    get!(m)
    reset!(m)

    m = MeanMetric{GraphMetric}()
    add!(m, GraphMetric(1,1,1,1,3,1,0,0))
    reset!(m)
end

max_stars_repo_name: lihebi/DAG-EQ

max_issues_repo_name: lihebi/DAG-EQ

max_forks_repo_name: lihebi/DAG-EQ


text:
using Test, BlockDiagonalFactors
using LinearAlgebra, SparseArrays, SuiteSparse

@testset "Testing BlockDiagonalFactors" begin
    @testset "full" begin
        A = rand(3, 3) + I
        B = rand(4, 4) + I
        C = rand(5, 5) + I
        Ms = [A, B, C]
        indices = [1, 1, 2, 1, 3, 1, 3, 2, 1]
        BDF = factorize(Ms, indices)
        @test BDF isa BlockFactors
        m = BDF.m
        n = BDF.n
        BDM = Matrix(blockdiag(map(sparse, Ms[indices])...))
        @test (m,n) == size(BDM)
        b = rand(m)
        x_BDM = BDM \ b
        @testset "backslash" begin
            x_BDF = BDF \ b
            @test x_BDF ≈ x_BDM
            x_BDF = BDF \ (im*b)
            @test x_BDF ≈ (im*x_BDM)
        end
        @testset "ldiv!(A, b)" begin
            x_BDF2 = copy(b)
            ldiv!(BDF, x_BDF2)
            @test x_BDF2 ≈ x_BDM
        end
        @testset "ldiv!(x, A, b)" begin
            x_BDF3 = similar(b)
            ldiv!(x_BDF3, BDF, b)
            @test x_BDF3 ≈ x_BDM
        end
        @testset "adjoint" begin
            @test (x_BDM' * (BDF \ b)) ≈ (b' * (BDF' \ x_BDM))'
        end
    end
    @testset "sparse" begin
        A = sprand(10, 10, 0.5) + I
        B = sprand(20, 20, 0.5) + I
        C = sprand(30, 30, 0.5) + I
        Ms = [A, B, C]
        indices = [1, 1, 2, 1, 3, 1, 3, 2, 1]
        BDF = factorize(Ms, indices)
        @test BDF isa SparseBlockFactors
        m = BDF.m
        n = BDF.n
        BDM = blockdiag(Ms[indices]...)
        @test (m,n) == size(BDM)
        b = rand(m)
        x_BDM = BDM \ b
        @testset "backslash" begin
            x_BDF = BDF \ b
            @test x_BDF ≈ x_BDM
            x_BDF = BDF \ (im*b)
            @test x_BDF ≈ (im*x_BDM)
        end
        @testset "ldiv!(A, b)" begin
            x_BDF2 = copy(b)
            ldiv!(BDF, x_BDF2)
            @test x_BDF2 ≈ x_BDM
        end
        @testset "ldiv!(x, A, b)" begin
            x_BDF3 = similar(b)
            ldiv!(x_BDF3, BDF, b)
            @test x_BDF3 ≈ x_BDM
        end
        @testset "adjoint" begin
            @test (x_BDM' * (BDF \ b)) ≈ (b' * (BDF' \ x_BDM))'
        end
    end
end


text:
function calc_forces(surf::TwoDSurf, vels::Vector{Float64})

    # First term in eqn (2.30) Ramesh et al. in coefficient form
    cnc = 2*pi*((surf.kinem.u + vels[1])*cos(surf.kinem.alpha)/surf.uref + (surf.kinem.hdot - vels[2])*sin(surf.kinem.alpha)/surf.uref)*(surf.a0[1] + surf.aterm[1]/2.)

    # Second term in eqn (2.30) Ramesh et al. in coefficient form
    cnnc = 2*pi*(3*surf.c*surf.a0dot[1]/(4*surf.uref) + surf.c*surf.adot[1]/(4*surf.uref) + surf.c*surf.adot[2]/(8*surf.uref))

    # Suction force given in eqn (2.31) Ramesh et al.
    cs = 2*pi*surf.a0[1]*surf.a0[1]

    #The components of normal force and moment from induced velocities are calulcated in dimensional units and nondimensionalized later
    nonl=0
    nonl_m=0
    for ib = 1:surf.ndiv-1
        nonl = nonl + (surf.uind[ib]*cos(surf.kinem.alpha) - surf.wind[ib]*sin(surf.kinem.alpha))*surf.bv[ib].s
        nonl_m = nonl_m + (surf.uind[ib]*cos(surf.kinem.alpha) - surf.wind[ib]*sin(surf.kinem.alpha))*surf.x[ib]*surf.bv[ib].s
    end
    nonl = nonl*2. /(surf.uref*surf.uref*surf.c)
    nonl_m = nonl_m*2. /(surf.uref*surf.uref*surf.c*surf.c)

    # Normal force coefficient
    cn = cnc + cnnc + nonl

    # Lift and drag coefficients
    cl = cn*cos(surf.kinem.alpha) + cs*sin(surf.kinem.alpha)
    cd = cn*sin(surf.kinem.alpha)-cs*cos(surf.kinem.alpha)

    #Pitching moment is clockwise or nose up positive
    cm = cn*surf.pvt - 2*pi*(((surf.kinem.u + vels[1])*cos(surf.kinem.alpha)/surf.uref + (surf.kinem.hdot - vels[2])*sin(surf.kinem.alpha)/surf.uref)*(surf.a0[1]/4. + surf.aterm[1]/4. - surf.aterm[2]/8.) + (surf.c/surf.uref)*(7. *surf.a0dot[1]/16. + 3. *surf.adot[1]/16. + surf.adot[2]/16. - surf.adot[3]/64.)) - nonl_m
    return cl, cd, cm
end

function writeStamp(dirname::String, t::Float64, surf::TwoDSurf, curfield::TwoDFlowField)
    dirvec = readdir()
    if dirname in dirvec
        rm(dirname, recursive=true)
    end
    mkdir(dirname)
    cd(dirname)

    f = open("timeKinem", "w")
    Serialization.serialize(f, ["#time \t", "alpha (deg) \t", "h/c \t", "u/uref \t"])
    DelimitedFiles.writedlm(f, [t, surf.kinem.alpha, surf.kinem.h, surf.kinem.u]')
    close(f)

    f = open("FourierCoeffs", "w")
    Serialization.serialize(f, ["#Fourier coeffs (0-n) \t", "d/dt (Fourier coeffs) \n"])
    matfour = zeros(surf.naterm+1, 2)
    matfour[:,1] = [surf.a0[1];surf.aterm[:]]
    matfour[:,2] = [surf.a0dot[1];surf.adot[:]]
    DelimitedFiles.writedlm(f, matfour)
    close(f)

    # f = open("forces", "w")
    # write(f, ["#cl \t", "cd \t", "cm \t", "Gamma \t", "cn \t", "cs \t", "cnc \t",
    #           "cnnc \t", "nonl \t", "cm_n \t", "cm_pvt \t", "nonl_m \n"])
    # cl, cd, cm, gamma, cn, cs, cnc, cncc, nonl, cm_n, cm_pvt, nonl_m = calc_forces_more(surf)
    # writedlm(f, [cl, cd, cm, gamma, cn, cs, cnc, cnnc, nonl, cm_n, cm_pvt, nonl_m])
    # close(f)

    f = open("tev", "w")
    Serialization.serialize(f, ["#strength \t", "x-position \t", "z-position \n"])
    tevmat = zeros(length(curfield.tev), 3)
    for i = 1:length(curfield.tev)
        tevmat[i,:] = [curfield.tev[i].s curfield.tev[i].x curfield.tev[i].z]
    end
    DelimitedFiles.writedlm(f, tevmat)
    close(f)

    f = open("lev", "w")
    Serialization.serialize(f, ["#strength \t", "x-position \t", "z-position \n"])
    levmat = zeros(length(curfield.lev), 3)
    for i = 1:length(curfield.lev)
        levmat[i,:] = [curfield.lev[i].s curfield.lev[i].x curfield.lev[i].z]
    end
    DelimitedFiles.writedlm(f, levmat)
    close(f)

    f = open("boundv", "w")
    Serialization.serialize(f, ["#strength \t", "x-position \t", "z-position \n"])
    bvmat = zeros(length(surf.bv), 3)
    for i = 1:length(surf.bv)
        bvmat[i,:] = [surf.bv[i].s surf.bv[i].x surf.bv[i].z]
    end
    DelimitedFiles.writedlm(f, bvmat)
    close(f)
    cd("..")
end

function writeStamp(dirname::String, t::Float64, surf::Vector{TwoDSurf}, curfield::TwoDFlowField)
    dirvec = readdir()
    if dirname in dirvec
        rm(dirname, recursive=true)
    end
    mkdir(dirname)
    cd(dirname)

    nsurf = length(surf)
    for i = 1:nsurf
        f = open("timeKinem-$i", "w")
        Serialization.Serialization.serialize(f, ["#time \t", "alpha (deg) \t", "h/c \t", "u/uref \t"])
        DelimitedFiles.writedlm(f, [t, surf[i].kinem.alpha, surf[i].kinem.h, surf[i].kinem.u]')
        close(f)

        f = open("FourierCoeffs-$i", "w")
        Serialization.serialize(f, ["#Fourier coeffs (0-n) \t", "d/dt (Fourier coeffs) \n"])
        matfour = zeros(surf[i].naterm+1, 2)
        matfour[:,1] = [surf[i].a0[1];surf[i].aterm[:]]
        matfour[:,2] = [surf[i].a0dot[1];surf[i].adot[:]]
        DelimitedFiles.writedlm(f, matfour)
        close(f)

        # f = open("forces", "w")
        # write(f, ["#cl \t", "cd \t", "cm \t", "Gamma \t", "cn \t", "cs \t", "cnc \t",
        #           "cnnc \t", "nonl \t", "cm_n \t", "cm_pvt \t", "nonl_m \n"])
        # cl, cd, cm, gamma, cn, cs, cnc, cncc, nonl, cm_n, cm_pvt, nonl_m = calc_forces_more(surf)
        # writedlm(f, [cl, cd, cm, gamma, cn, cs, cnc, cnnc, nonl, cm_n, cm_pvt, nonl_m])
        # close(f)
    end

    f = open("tev", "w")
    Serialization.serialize(f, ["#strength \t", "x-position \t", "z-position \n"])
    tevmat = zeros(length(curfield.tev), 3)
    for i = 1:length(curfield.tev)
        tevmat[i,:] = [curfield.tev[i].s curfield.tev[i].x curfield.tev[i].z]
    end
    DelimitedFiles.writedlm(f, tevmat)
    close(f)

    f = open("lev", "w")
    Serialization.serialize(f, ["#strength \t", "x-position \t", "z-position \n"])
    levmat = zeros(length(curfield.lev), 3)
    levcount = 0
    for i = 1:length(curfield.lev)
        if curfield.lev[i].vc != 0.
            levcount += 1
            levmat[levcount,:] = [curfield.lev[i].s curfield.lev[i].x curfield.lev[i].z]
        end
    end

    levmat = levmat[1:levcount,:]
    DelimitedFiles.writedlm(f, levmat)
    close(f)

    for is = 1:nsurf
        f = open("boundv-$is", "w")
        Serialization.serialize(f, ["#strength \t", "x-position \t", "z-position \n"])
        bvmat = zeros(length(surf[is].bv), 3)
        for i = 1:length(surf[is].bv)
            bvmat[i,:] = [surf[is].bv[i].s surf[is].bv[i].x surf[is].bv[i].z]
        end
        DelimitedFiles.writedlm(f, bvmat)
        close(f)
    end

    cd("..")
end

function calc_forces(surf::TwoDSurfThick, int_wax_prev :: Vector{Float64}, dt :: Float64)
    l_x = zeros(surf.ndiv)
    t_z = zeros(surf.ndiv)
    l_z = zeros(surf.ndiv)
    t_x = zeros(surf.ndiv)
    ws_x = zeros(surf.ndiv)
    ws_z = zeros(surf.ndiv)
    wa_x = zeros(surf.ndiv)
    wa_z = zeros(surf.ndiv)

    #All terms involving l_x, l_z (having the singularity) are integrated in theta to eliminate the singularity

    for i = 1:surf.ndiv
        l_x[i] = surf.a0[1]*(1+cos(surf.theta[i]))
        l_z[i] = -surf.a0[1]*sin(surf.theta[i])

        ws_x[i] = 0.5*(surf.uind_u[i] + surf.uind_l[i])
        wa_x[i] = 0.5*(surf.uind_u[i] - surf.uind_l[i])

        ws_z[i] = 0.5*(surf.wind_u[i] - surf.wind_l[i])
        wa_z[i] = 0.5*(surf.wind_u[i] + surf.wind_l[i])

        for ia = 1:surf.naterm
            l_x[i] += surf.aterm[ia]*sin(ia*surf.theta[i])*sin(surf.theta[i])
            t_x[i] -= surf.bterm[ia]*cos(ia*surf.theta[i])
            l_z[i] += surf.aterm[ia]*cos(ia*surf.theta[i])*sin(surf.theta[i])
            t_z[i] += surf.bterm[ia]*sin(ia*surf.theta[i])
        end
        l_x[i] = l_x[i]*surf.uref*surf.c/2.
        t_x[i] = t_x[i]*surf.uref
        l_z[i] = l_z[i]*surf.uref*surf.c/2.
        t_z[i] = t_z[i]*surf.uref
    end

    int_lxtx = simpleTrapz(l_x.*t_x, surf.theta)
    int_lxwsx = simpleTrapz(l_x.*ws_x, surf.theta)
    int_txwax = simpleTrapz(t_x.*wa_x, surf.x)
    int_wsxwax = simpleTrapz(ws_x.*wa_x, surf.x)

    int_lztz = simpleTrapz(l_z.*t_z, surf.theta)
    int_lzwsz = simpleTrapz(l_z.*ws_z, surf.theta)
    int_tzwaz = simpleTrapz(t_z.*wa_z, surf.x)
    int_wszwaz = simpleTrapz(ws_z.*wa_z, surf.x)

    cnc1 = (4. /(surf.uref*surf.uref*surf.c))*(int_lxtx + int_lxwsx + int_txwax +
    int_wsxwax + int_lztz + int_lzwsz + int_tzwaz + int_wszwaz)

    int_wax_kinem = simpleTrapz(wa_x.*(surf.kinem.u*cos(surf.kinem.alpha) +
    surf.kinem.hdot*sin(surf.kinem.alpha) .- surf.kinem.alphadot*surf.cam), surf.x)
    int_lx_etac = simpleTrapz(l_x*surf.kinem.alphadot.*surf.cam, surf.theta)

    cnc2 = 2*pi*(surf.kinem.u*cos(surf.kinem.alpha)/surf.uref +
    surf.kinem.hdot*sin(surf.kinem.alpha)/surf.uref)*(surf.a0[1] + surf.aterm[1]/2.) +
    (4. /(surf.uref*surf.uref*surf.c))*(int_wax_kinem - int_lx_etac)

    int_tx_etat = simpleTrapz((t_x + wa_x)*surf.kinem.alphadot.*surf.thick, surf.x)

    int_wsz_kinem = simpleTrapz(ws_z.*(surf.kinem.u*sin(surf.kinem.alpha)/surf.uref -
    surf.kinem.hdot*cos(surf.kinem.alpha) .+ surf.kinem.alphadot*(surf.x .- surf.pvt*surf.c)), surf.x)

    cnc3 = pi*surf.bterm[1]*(surf.kinem.u*sin(surf.kinem.alpha)/surf.uref -
    surf.kinem.hdot*cos(surf.kinem.alpha)/surf.uref -
    surf.kinem.alphadot*surf.c*(0.5 - surf.pvt)/surf.uref) -
    surf.kinem.alphadot*surf.c*pi*surf.bterm[2]/(4*surf.uref) +
    (4. /(surf.uref*surf.uref*surf.c))*(int_wsz_kinem - int_tx_etat)

    int_wax = zeros(surf.ndiv)
    d_int_wax = zeros(surf.ndiv)

    for i = 1:surf.ndiv
        int_wax[i] = simpleTrapz(wa_x[1:i], surf.x[1:i])

        d_int_wax[i] = (int_wax[i] - int_wax_prev[i])/dt
    end

    cnnc = (2. *pi*surf.c/surf.uref)*(3. *surf.a0dot[1]/4. + surf.adot[1]/4. + surf.adot[2]/8.) +
    (4. /(surf.uref*surf.uref*surf.c))*simpleTrapz(d_int_wax, surf.x)

    # Suction force given in eqn (2.31) Ramesh et al.
    cs = 2*pi*surf.a0[1]*surf.a0[1]

    # Normal force coefficient
    cn = cnc1 + cnc2 + cnc3 + cnnc

    # Lift and drag coefficients
    cl = cn*cos(surf.kinem.alpha) + cs*sin(surf.kinem.alpha)
    cd = cn*sin(surf.kinem.alpha)-cs*cos(surf.kinem.alpha)

    #Pitching moment is clockwise or nose up positive

    return cnc1, cnc2, cnc3, cnnc, cn, cs, cl, cd
end

function calc_q_cp(surf::TwoDSurf, vels::Vector{Float64})

    q_in_u = zeros(surf.ndiv)
    q_in_l = zeros(surf.ndiv)
    q_out_u = zeros(surf.ndiv)
    q_out_l = zeros(surf.ndiv)
    p_in = zeros(surf.ndiv)
    p_out = zeros(surf.ndiv)
    gam = zeros(surf.ndiv)
    gamint = zeros(surf.ndiv)
    p_com = zeros(surf.ndiv)
    gammod = zeros(surf.ndiv)
    q_com_u = zeros(surf.ndiv)
    q_com_l = zeros(surf.ndiv)

    rho = 0.016

    for i = 1:surf.ndiv
        udash = surf.kinem.u*cos(surf.kinem.alpha) + surf.kinem.hdot*sin(surf.kinem.alpha) + surf.uind[i]*cos(surf.kinem.alpha) - surf.wind[i]*sin(surf.kinem.alpha)

        q_in_u[i] = (sqrt(surf.c)*surf.uref*surf.a0[1] + sqrt(surf.x[i])*udash)/(sqrt(surf.x[i]+rho*surf.c/2))
        q_in_l[i] = (sqrt(surf.c)*surf.uref*surf.a0[1] - sqrt(surf.x[i])*udash)/(sqrt(surf.x[i]+rho*surf.c/2))
        p_in[i] = 8*surf.uref*sqrt(surf.c)*surf.a0[1]*sqrt(surf.x[i])*udash/(rho*surf.c + 2*surf.x[i]) + 4*surf.uref*sqrt(surf.c)*surf.a0dot[1]*sqrt(surf.x[i])
    end
    for i = 2:surf.ndiv
        udash = surf.kinem.u*cos(surf.kinem.alpha) + surf.kinem.hdot*sin(surf.kinem.alpha) + surf.uind[i]*cos(surf.kinem.alpha) - surf.wind[i]*sin(surf.kinem.alpha)
        gam[i] = surf.a0[1]*cot(surf.theta[i]/2)
        gamint[i] = surf.a0dot[1]*(1. + cos(surf.theta[i]))
        for n = 1:surf.naterm
            gam[i] += surf.aterm[n]*sin(n*surf.theta[i])
            gammod[i] += surf.aterm[n]*sin(n*surf.theta[i])
            if n == 1
                gamint[i] += surf.adot[1]*(surf.theta[i]/2 - sin(2*surf.theta[i])/4)
            else
                gamint[i] += surf.adot[n]/2*(sin((n-1)*surf.theta[i])/(n-1) - sin((n+1)*surf.theta[i])/(n+1))
            end
        end
        gam[i] = gam[i]*surf.uref
        gammod[i] = gammod[i]*surf.uref
        gamint[i] = gamint[i]*surf.uref*surf.c

        q_out_u[i] = udash + gam[i]
        q_out_l[i] = udash - gam[i]
        p_out[i] = 2*udash*gam[i] + gamint[i]
        p_com[i] = 2*udash*surf.uref*surf.a0[1]*(cos(surf.theta[i]/2)-1)/sin(surf.theta[i]/2) + 2*udash*gammod[i] + gamint[i] + 8*surf.uref*udash*surf.c*surf.a0[1]*sin(surf.theta[i]/2)/(rho*surf.c + 2*surf.c*(sin(surf.theta[i]/2))^2)
        q_com_u[i] = surf.uref*surf.a0[1]*(cos(surf.theta[i]/2) - 1)/sin(surf.theta[i]/2) + surf.uref*gammod[i] + (sqrt(surf.c)*surf.uref*surf.a0[1] + sqrt(surf.x[i])*udash)/sqrt(surf.x[i] + rho*surf.c/2)
        q_com_l[i] = -surf.uref*surf.a0[1]*(cos(surf.theta[i]/2) - 1)/sin(surf.theta[i]/2) - surf.uref*gammod[i] + (-sqrt(surf.c)*surf.uref*surf.a0[1] + sqrt(surf.x[i])*udash)/sqrt(surf.x[i] + rho*surf.c/2)
    end

    q_com_u[1] = surf.uref*surf.a0[1]/sqrt(0.5*rho)
    q_com_l[1] = surf.uref*surf.a0[1]/sqrt(0.5*rho)

    return q_in_l, q_in_u, q_out_l, q_out_u, p_in, p_out, p_com, q_com_u, q_com_l
end

function calc_edgeVel(surf::TwoDSurf, vels::Vector{Float64})

    gammod = zeros(surf.ndiv)
    q_com_u = zeros(surf.ndiv)
    q_com_l = zeros(surf.ndiv)

    q_com_u[1] = surf.uref*surf.a0[1]/sqrt(0.5*surf.rho)
    q_com_l[1] = surf.uref*surf.a0[1]/sqrt(0.5*surf.rho)

    for i = 2:surf.ndiv
        udash = (surf.kinem.u + vels[1])*cos(surf.kinem.alpha) + (surf.kinem.hdot - vels[2])*sin(surf.kinem.alpha) + surf.uind[i]*cos(surf.kinem.alpha) - surf.wind[i]*sin(surf.kinem.alpha)

        for n = 1:surf.naterm
            gammod[i] += surf.aterm[n]*sin(n*surf.theta[i])
        end
        gammod[i] = gammod[i]*surf.uref

        q_com_u[i] = surf.uref*surf.a0[1]*(cos(surf.theta[i]/2) - 1)/sin(surf.theta[i]/2) + surf.uref*gammod[i] + (sqrt(surf.c)*surf.uref*surf.a0[1] + sqrt(surf.x[i])*udash)/sqrt(surf.x[i] + surf.rho*surf.c/2)
        q_com_l[i] = -surf.uref*surf.a0[1]*(cos(surf.theta[i]/2) - 1)/sin(surf.theta[i]/2) - surf.uref*gammod[i] + (-sqrt(surf.c)*surf.uref*surf.a0[1] + sqrt(surf.x[i])*udash)/sqrt(surf.x[i] + surf.rho*surf.c/2)
    end

    return q_com_u, q_com_l
end

function calc_q_cp(surf::TwoDSurfThick, vels::Vector{Float64})
    uc = zeros(surf.ndiv)
    ut = zeros(surf.ndiv)
    wt = zeros(surf.ndiv)
    wc = zeros(surf.ndiv)
    for ib = 1:surf.ndiv
        uc[ib] = sqrt(2. /(1-cos(surf.theta[ib]) + surf.rho/surf.c))*(surf.a0[1]*cos(surf.theta[ib]/2))
        ut[ib] = 0.
        for ia = 1:surf.naterm
            uc[ib] += sqrt(2. /(1-cos(surf.theta[ib]) + surf.rho/surf.c))*surf.aterm[ia]*sin(ia*surf.theta[ib])*sin(surf.theta[ib]/2)
            ut[ib] -= sqrt(surf.x[ib]./(surf.x[ib] + 0.5*surf.rho))*surf.bterm[ia]*cos(ia*surf.theta[ib])
        end
        uc[ib] = uc[ib]*surf.uref
        ut[ib] = ut[ib]*surf.uref
    end

    # w1t[1] = 2*w1t[2] - w1t[3]
    # w1t[surf.ndiv] = 2*w1t[surf.ndiv-1] - w1t[surf.ndiv-2]
    # u1t[1] = 2*u1t[2] - u1t[3]
    # u1t[surf.ndiv] = 2*u1t[surf.ndiv-1] - u1t[surf.ndiv-2]
    #w1c[1] = 2*w1c[2] - w1c[3]
    #w1c[surf.ndiv] = 2*w1c[surf.ndiv-1] - w1c[surf.ndiv-2]
    #u1c[1] = 2*u1c[2] - u1c[3]
    #u1c[surf.ndiv] = 2*u1c[surf.ndiv-1] - u1c[surf.ndiv-2]

    #lamb = mean(surf.cam_slope[1:5])

    vu = zeros(surf.ndiv)
    vl = zeros(surf.ndiv)

    for i = 1:surf.ndiv

        wlz = 0.5*(surf.wind_u[i]*cos(surf.kinem.alpha) + surf.uind_u[i]*sin(surf.kinem.alpha) +
                   surf.wind_l[i]*cos(surf.kinem.alpha) + surf.uind_l[i]*sin(surf.kinem.alpha))

        wtz = 0.5*(surf.wind_u[i]*cos(surf.kinem.alpha) + surf.uind_u[i]*sin(surf.kinem.alpha) -
                   surf.wind_l[i]*cos(surf.kinem.alpha) - surf.uind_l[i]*sin(surf.kinem.alpha))

        wtx = 0.5*(surf.uind_u[i]*cos(surf.kinem.alpha) - surf.wind_u[i]*sin(surf.kinem.alpha) +
                   surf.uind_l[i]*cos(surf.kinem.alpha) - surf.wind_l[i]*sin(surf.kinem.alpha))

        wlx = 0.5*(surf.uind_u[i]*cos(surf.kinem.alpha) - surf.wind_u[i]*sin(surf.kinem.alpha) -
                   surf.uind_l[i]*cos(surf.kinem.alpha) + surf.wind_l[i]*sin(surf.kinem.alpha))

        u_u = uc[i] + ut[i] + sqrt(surf.x[i]./(surf.x[i] + 0.5*surf.rho))*(surf.uind_u[i] + (surf.kinem.u + vels[1])*cos(surf.kinem.alpha) + (surf.kinem.hdot - vels[2])*sin(surf.kinem.alpha) - surf.kinem.alphadot*(surf.cam[i] + surf.thick[i]))

        w_u = 0#uc[i]*(surf.cam_slope[i] + surf.thick_slope[i]) + ut[i]*(surf.thick_slope[i] + surf.cam_slope[i]) +
            sqrt(surf.x[i]./(surf.x[i] + 0.5*surf.rho))*((surf.thick_slope[i] + surf.cam_slope[i])*((surf.kinem.u + vels[1])*cos(surf.kinem.alpha) +
                                                       (surf.kinem.hdot - vels[2])*sin(surf.kinem.alpha) -
                                                       surf.kinem.alphadot*surf.cam[i] + wtx)
        + (surf.cam_slope[i] + surf.thick_slope[i])*(wlx - surf.kinem.alphadot*surf.thick[i]))

        u_l = -uc[i] + ut[i] + sqrt(surf.x[i]./(surf.x[i] + 0.5*surf.rho))*(surf.uind_l[i] + (surf.kinem.u + vels[1])*cos(surf.kinem.alpha) +
            (surf.kinem.hdot - vels[2])*sin(surf.kinem.alpha) - surf.kinem.alphadot*(surf.cam[i] - surf.thick[i]))

        w_l = 0#uc[i]*(-surf.cam_slope[i] + surf.thick_slope[i]) + ut[i]*(-surf.thick_slope[i] + surf.cam_slope[i]) +
            sqrt(surf.x[i]./(surf.x[i] + 0.5*surf.rho))*((-surf.thick_slope[i] + surf.cam_slope[i])*((surf.kinem.u + vels[1])*cos(surf.kinem.alpha) +
                                                       (surf.kinem.hdot - vels[2])*sin(surf.kinem.alpha) -
                                                       surf.kinem.alphadot*surf.cam[i] + wtx)
        + (-surf.cam_slope[i] + surf.thick_slope[i])*(wlx - surf.kinem.alphadot*surf.thick[i]))

        vu[i] = sqrt(u_u^2 + w_u^2)
        vl[i] = sqrt(u_l^2 + w_l^2)

    end
    return vu, vl
end

# function calc_q_cp(surf::TwoDSurfThick, vels::Vector{Float64})
#      u1c = zeros(surf.ndiv)
#      u1t = zeros(surf.ndiv)
#      w1t = zeros(surf.ndiv)
#      w1c = zeros(surf.ndiv)
#     for ib = 2:surf.ndiv-1
#         u1c[ib] = surf.a0[1]*(1 + cos(surf.theta[ib]))/sin(surf.theta[ib])
#         w1t[ib] = 0. #surf.b0[1]*(1 + cos(surf.theta[ib]))/sin(surf.theta[ib])
#         #w1t[ib] += surf.b0[2]*(1 - cos(surf.theta[ib]))/sin(surf.theta[ib])
#         w1c[ib] = -surf.a0[1]
#         u1t[ib] = 0.#surf.b0[1] - surf.b0[2]
#         for ia = 1:surf.naterm
#             u1c[ib] += surf.aterm[ia]*sin(ia*surf.theta[ib])
#             u1t[ib] -= surf.bterm[ia]*cos(ia*surf.theta[ib])
#             w1t[ib] += surf.bterm[ia]*sin(ia*surf.theta[ib])
#             w1c[ib] += surf.aterm[ia]*cos(ia*surf.theta[ib])
#         end
#         u1c[ib] = u1c[ib]*surf.uref
#         u1t[ib] = u1t[ib]*surf.uref
#         w1t[ib] = w1t[ib]*surf.uref
#         w1c[ib] = w1c[ib]*surf.uref
#     end

#     # w1t[1] = 2*w1t[2] - w1t[3]
#     # w1t[surf.ndiv] = 2*w1t[surf.ndiv-1] - w1t[surf.ndiv-2]
#     # u1t[1] = 2*u1t[2] - u1t[3]
#     # u1t[surf.ndiv] = 2*u1t[surf.ndiv-1] - u1t[surf.ndiv-2]
#     w1c[1] = 2*w1c[2] - w1c[3]
#     #w1c[surf.ndiv] = 2*w1c[surf.ndiv-1] - w1c[surf.ndiv-2]
#     u1c[1] = 2*u1c[2] - u1c[3]
#     #u1c[surf.ndiv] = 2*u1c[surf.ndiv-1] - u1c[surf.ndiv-2]

#     for ib in [1, surf.ndiv]
#         w1t[ib] = 0. #surf.b0[1]*(1 + cos(surf.theta[ib]))/sin(surf.theta[ib])
#         #w1t[ib] += surf.b0[2]*(1 - cos(surf.theta[ib]))/sin(surf.theta[ib])
#         w1c[ib] = -surf.a0[1]
#         u1t[ib] = 0.#surf.b0[1] - surf.b0[2]
#         for ia = 1:surf.naterm
#             u1c[ib] += surf.aterm[ia]*sin(ia*surf.theta[ib])
#             u1t[ib] -= surf.bterm[ia]*cos(ia*surf.theta[ib])
#             w1t[ib] += surf.bterm[ia]*sin(ia*surf.theta[ib])
#             w1c[ib] += surf.aterm[ia]*cos(ia*surf.theta[ib])
#         end
#         u1t[ib] = u1t[ib]*surf.uref
#         w1t[ib] = w1t[ib]*surf.uref
#         w1c[ib] = w1c[ib]*surf.uref
#     end

#     lamb = mean(surf.cam_slope[1:5])
#     rho = surf.rho
#     vu = zeros(surf.ndiv)
#     vl = zeros(surf.ndiv)
#     vu_van = zeros(surf.ndiv)
#     vu2_van = zeros(surf.ndiv)
#     vl_van = zeros(surf.ndiv)
#     vl2_van = zeros(surf.ndiv)

#     for i = 1:surf.ndiv
#         u_u = u1c[i] + u1t[i] + surf.uind_u[i] + (surf.kinem.u + vels[1])*cos(surf.kinem.alpha) +
#             (surf.kinem.hdot - vels[2])*sin(surf.kinem.alpha) - surf.kinem.alphadot*(surf.cam[i] + surf.thick[i])

#         w_u = w1c[i] + w1t[i] + surf.wind_u[i] + (surf.kinem.u + vels[1])*sin(surf.kinem.alpha) -
#             (surf.kinem.hdot - vels[2])*cos(surf.kinem.alpha) + surf.kinem.alphadot*(surf.x[i] - surf.pvt*surf.c)

#         vu[i] = sqrt(u_u^2 + w_u^2)

#         u_l = -u1c[i] + u1t[i] + surf.uind_l[i] + (surf.kinem.u + vels[1])*cos(surf.kinem.alpha) +
#         (surf.kinem.hdot - vels[2])*sin(surf.kinem.alpha) - surf.kinem.alphadot*(surf.cam[i] - surf.thick[i])

#         w_l = w1c[i] - w1t[i] + surf.wind_l[i] + (surf.kinem.u + vels[1])*sin(surf.kinem.alpha) -
#         (surf.kinem.hdot - vels[2])*cos(surf.kinem.alpha) + surf.kinem.alphadot*(surf.x[i] - surf.pvt*surf.c)

#         vl[i] = sqrt(u_l^2 + w_l^2)

#         vu_van[i] = sqrt((surf.x[i]+lamb*sqrt(2*rho*surf.x[i]))/(surf.x[i] + lamb*sqrt(2*rho*surf.x[i]) + rho/2.))*(vu[i]/surf.uref+rho/(4*surf.x[i]))
#         #vu2_van[i] = sqrt((surf.x[i]+lamb*sqrt(2*rho*surf.x[i]))/(surf.x[i] + lamb*sqrt(2*rho*surf.x[i]) + rho/2.))*(vu[i]/surf.uref)
#         vu2_van[i] = 1. /sqrt(1. + surf.thick_slope[i]^2)*(vu[i]/surf.uref)
#         vl_van[i] = sqrt((surf.x[i]-lamb*sqrt(2*rho*surf.x[i]))/(surf.x[i] - lamb*sqrt(2*rho*surf.x[i]) + rho/2.))*(vl[i]/surf.uref+rho/(4*surf.x[i]))
#         vl2_van[i] = 1. /sqrt(1. + surf.thick_slope[i]^2)*(vl[i]/surf.uref)
#     end
#     return vu, vl, vu_van, vu2_van, vl_van, vl2_van
# end

function writeStamp(dirname::String, t::Float64, surf::TwoDSurfThick, curfield::TwoDFlowField)
    dirvec = readdir()
    if dirname in dirvec
        rm(dirname, recursive=true)
    end
    mkdir(dirname)
    cd(dirname)

    f = open("timeKinem", "w")
    Serialization.serialize(f, ["#time \t", "alpha (deg) \t", "h/c \t", "u/uref \t"])
    DelimitedFiles.writedlm(f, [t, surf.kinem.alpha, surf.kinem.h, surf.kinem.u]')
    close(f)

    f = open("FourierCoeffsA", "w")
    Serialization.serialize(f, ["#Fourier coeffs (0-n) \t", "d/dt (Fourier coeffs) \n"])
    mat = zeros(surf.naterm+1, 2)
    mat[:,1] = [surf.a0[1];surf.aterm[:]]
    mat[1:4,2] = [surf.a0dot[1];surf.adot[:]]
    DelimitedFiles.writedlm(f, mat)
    close(f)

    f = open("FourierCoeffsB", "w")
    Serialization.serialize(f, ["#Fourier coeffs (0-n) \t", "d/dt (Fourier coeffs) \n"])
    mat = zeros(surf.naterm, 1)
    mat[:,1] = surf.bterm[:]
    DelimitedFiles.writedlm(f, mat)
    close(f)

    # f = open("forces", "w")
    # write(f, ["#cl \t", "cd \t", "cm \t", "Gamma \t", "cn \t", "cs \t", "cnc \t",
    #           "cnnc \t", "nonl \t", "cm_n \t", "cm_pvt \t", "nonl_m \n"])
    # cl, cd, cm, gamma, cn, cs, cnc, cncc, nonl, cm_n, cm_pvt, nonl_m = calc_forces_more(surf)
    # writedlm(f, [cl, cd, cm, gamma, cn, cs, cnc, cnnc, nonl, cm_n, cm_pvt, nonl_m])
    # close(f)

    f = open("tev", "w")
    Serialization.serialize(f, ["#strength \t", "x-position \t", "z-position \n"])
    tevmat = zeros(length(curfield.tev), 3)
    for i = 1:length(curfield.tev)
        tevmat[i,:] = [curfield.tev[i].s curfield.tev[i].x curfield.tev[i].z]
    end
    DelimitedFiles.writedlm(f, tevmat)
    close(f)

    f = open("lev", "w")
    Serialization.serialize(f, ["#strength \t", "x-position \t", "z-position \n"])
    levmat = zeros(length(curfield.lev), 3)
    for i = 1:length(curfield.lev)
        levmat[i,:] = [curfield.lev[i].s curfield.lev[i].x curfield.lev[i].z]
    end
    DelimitedFiles.writedlm(f, levmat)
    close(f)

    f = open("boundv", "w")
    Serialization.serialize(f, ["#strength \t", "x-position \t", "z-position \n"])
    bvmat = zeros(length(surf.bv), 3)
    for i = 1:length(surf.bv)
        bvmat[i,:] = [surf.bv[i].s surf.bv[i].x surf.bv[i].z]
    end
    DelimitedFiles.writedlm(f, bvmat)
    close(f)

    f = open("src", "w")
    Serialization.serialize(f, ["#strength \t", "x-position \t", "z-position \n"])
    bvmat = zeros(length(surf.src), 3)
    for i = 1:length(surf.src)
        bvmat[i,:] = [surf.src[i].s surf.src[i].x surf.src[i].z]
    end
    DelimitedFiles.writedlm(f, bvmat)
    close(f)

    cd("..")
end

function getEndCycle(mat, k)

    T = pi/k

    end_cycle = mat[end,1]

    #Find number of cycles
    ncyc = 0
    for i = 1:1000
        t_l = real(i)*T
        if t_l > end_cycle
            break
        end
        ncyc = ncyc + 1
    end

    start_t = real(ncyc-1)*T
    end_t = real(ncyc)*T
    start_ind = argmin(abs.(mat[:,1] .- start_t))
    end_ind = argmin(abs.(mat[:,1] .- end_t))

    nlast = end_ind - start_ind + 1

    newmat = zeros(nlast, 8)
    newmat[:,1] = (mat[start_ind:end_ind,1] .- start_t)/T
    for i = 2:7
        newmat[:,i] = mat[start_ind:end_ind,i]
    end

    f = open("resultsSummaryEndCycle", "w")
    Serialization.serialize(f, ["#t/T \t", "alpha (deg) \t", "h/c \t", "u/uref \t", "A0 \t", "Cl \t", "Cd \t", "Cm \n"])
    DelimitedFiles.writedlm(f, newmat)
    close(f)

    return newmat
end


text:
__precompile__()

module MeshCat

using GeometryTypes
using CoordinateTransformations
using Rotations: rotation_between
using ZMQ

using Colors: Color, Colorant, RGB, RGBA, alpha
using StaticArrays: StaticVector, SVector, SDiagonal
using GeometryTypes: raw
using Parameters: @with_kw
using Base.Random: UUID, uuid1
using URIParser: escape
using DocStringExtensions

import Base: delete!, length
import MsgPack: pack, Ext
import GeometryTypes: origin, radius

export Visualizer,
       ViewerWindow,
       IJuliaCell,
       HyperEllipsoid,
       HyperCylinder,
       PointCloud,
       PointsMaterial,
       MeshLambertMaterial,
       MeshBasicMaterial,
       MeshPhongMaterial,
       Texture,
       PngImage,
       Mesh,
       Points,
       setobject!,
       settransform!,
       delete!,
       url

include("servers/trees.jl")
include("servers/zmqserver.jl")
using .ZMQServer: ZMQWebSocketBridge, zmq_url, web_url, VIEWER_ROOT

include("geometry.jl")
include("objects.jl")
include("commands.jl")
include("lowering.jl")
include("msgpack.jl")
include("visualizer.jl")
include("ijulia.jl")

end


text:
function get_schemas end
function get_tables end
function get_table_oid end
function get_pks end
function get_fks end
function get_columns_types end
function analyse_db_schema end
function check_if_table_is_partition_of_another_table end
function check_if_table_or_partition_exists end
function get_custom_types end


text:
using Test
import Pluto
import Pluto: update_run!, WorkspaceManager, ClientSession, ServerSession, Notebook, Cell

@testset "Rich output" begin

    🍭 = ServerSession()
    🍭.options.evaluation.workspace_use_distributed = false
    fakeclient = ClientSession(:fake, nothing)
    🍭.connected_clients[fakeclient.id] = fakeclient

    @testset "Tree viewer" begin
        @testset "Basics" begin
            notebook = Notebook([
                    Cell("[1,1,[1]]"),
                    Cell("Dict(:a => [:b, :c])"),
                    Cell("[3, Dict()]"),
                    Cell("[4,[3, Dict()]]"),
                    Cell("[5, missing, 5]"),
                    Cell("[]"),
                    Cell("(7,7)"),
                    Cell("(a=8,b=[8])"),
                    Cell("Ref(9)"),
                    Cell("Vector(undef, 10)"),
                    Cell("struct Eleven x end"),
                    Cell("Eleven(12)"),

                    Cell("struct Amazing{T} x::T end"),
                    Cell("Amazing(14)"),

                    Cell("1:15"),

                    Cell("""begin
                        mutable struct A x; y end
                        a = A(16, 0)
                        a.y = [1, (2, A(3, (r=4, t=(5 => Dict(6 => Ref(a))))))]
                        a
                    end""")
                ])
            fakeclient.connected_notebook = notebook

            update_run!(🍭, notebook, notebook.cells)

            @test notebook.cells[1].repr_mime isa MIME"application/vnd.pluto.tree+object"
            @test notebook.cells[2].repr_mime isa MIME"application/vnd.pluto.tree+object"
            @test notebook.cells[3].repr_mime isa MIME"application/vnd.pluto.tree+object"
            @test notebook.cells[4].repr_mime isa MIME"application/vnd.pluto.tree+object"
            @test notebook.cells[5].repr_mime isa MIME"application/vnd.pluto.tree+object"
            @test notebook.cells[6].repr_mime isa MIME"application/vnd.pluto.tree+object"
            @test notebook.cells[7].repr_mime isa MIME"application/vnd.pluto.tree+object"
            @test notebook.cells[8].repr_mime isa MIME"application/vnd.pluto.tree+object"
            @test notebook.cells[9].repr_mime isa MIME"application/vnd.pluto.tree+object"
            @test notebook.cells[10].repr_mime isa MIME"application/vnd.pluto.tree+object"
            @test notebook.cells[1].output_repr isa Dict
            @test notebook.cells[2].output_repr isa Dict
            @test notebook.cells[3].output_repr isa Dict
            @test notebook.cells[4].output_repr isa Dict
            @test notebook.cells[5].output_repr isa Dict
            @test notebook.cells[6].output_repr isa Dict
            @test notebook.cells[7].output_repr isa Dict
            @test notebook.cells[8].output_repr isa Dict
            @test notebook.cells[9].output_repr isa Dict
            @test notebook.cells[10].output_repr isa Dict

            @test notebook.cells[12].repr_mime isa MIME"application/vnd.pluto.tree+object"
            @test notebook.cells[12].output_repr isa Dict
            @test notebook.cells[14].repr_mime isa MIME"application/vnd.pluto.tree+object"
            @test notebook.cells[14].output_repr isa Dict

            @test notebook.cells[15].repr_mime isa MIME"text/plain"

            @test notebook.cells[16].repr_mime isa MIME"application/vnd.pluto.tree+object"
            @test notebook.cells[16].output_repr isa Dict
            @test occursin("circular", notebook.cells[16].output_repr |> string)

            WorkspaceManager.unmake_workspace((🍭, notebook))
        end

        @testset "Overloaded Base.show" begin

            notebook = Notebook([
                Cell("""begin
                    struct A x end
                    Base.show(io::IO, ::MIME"image/png", x::Vector{A}) = print(io, "1")
                    [A(1), A(1)]
                end"""),

                Cell("""begin
                    struct B x end
                    Base.show(io::IO, ::MIME"text/html", x::B) = print(io, "2")
                    B(2)
                end"""),

                Cell("""begin
                    struct C x end
                    Base.show(io::IO, ::MIME"text/plain", x::C) = print(io, "3")
                    C(3)
                end"""),
            ])
            fakeclient.connected_notebook = notebook

            update_run!(🍭, notebook, notebook.cells)

            @test notebook.cells[1].repr_mime isa MIME"image/png"
            @test notebook.cells[1].output_repr == codeunits("1")
            @test notebook.cells[2].repr_mime isa MIME"text/html"
            @test notebook.cells[2].output_repr == "2"
            @test notebook.cells[3].repr_mime isa MIME"text/plain"
            @test notebook.cells[3].output_repr == "3"

            WorkspaceManager.unmake_workspace((🍭, notebook))
        end

        @testset "Special arrays" begin

            notebook = Notebook([
                Cell("using OffsetArrays"),
                Cell("OffsetArray(zeros(3), 20:22)"),
            ])
            fakeclient.connected_notebook = notebook

            update_run!(🍭, notebook, notebook.cells)

            @test notebook.cells[2].repr_mime isa MIME"application/vnd.pluto.tree+object"
            s = string(notebook.cells[2].output_repr)
            @test occursin("OffsetArray", s)
            @test occursin("21", s)
            if VERSION >= v"1.3"
                # once in the prefix, once as index
                @test count("22", s) >= 2
            end

            WorkspaceManager.unmake_workspace((🍭, notebook))
        end
    end

    @testset "Table viewer" begin
        notebook = Notebook([
                Cell("using DataFrames, Tables"),
                Cell("DataFrame()"),
                Cell("DataFrame(:a => [])"),
                Cell("DataFrame(:a => [1,2,3], :b => [999, 5, 6])"),
                Cell("DataFrame(rand(20,20))"),
                Cell("DataFrame(rand(2000,20))"),
                Cell("DataFrame(rand(20,2000))"),
                Cell("@view DataFrame(rand(100,3))[:, 2:2]"),
                Cell("@view DataFrame(rand(3,100))[2:2, :]"),
                Cell("DataFrame"),
                Cell("Tables.table(rand(11,11))"),
                Cell("Tables.table(rand(120,120))"),
            ])
        fakeclient.connected_notebook = notebook

        update_run!(🍭, notebook, notebook.cells)

        @test notebook.cells[2].repr_mime isa MIME"application/vnd.pluto.table+object"
        @test notebook.cells[3].repr_mime isa MIME"application/vnd.pluto.table+object"
        @test notebook.cells[4].repr_mime isa MIME"application/vnd.pluto.table+object"
        @test notebook.cells[5].repr_mime isa MIME"application/vnd.pluto.table+object"
        @test notebook.cells[6].repr_mime isa MIME"application/vnd.pluto.table+object"
        @test notebook.cells[7].repr_mime isa MIME"application/vnd.pluto.table+object"
        @test notebook.cells[8].repr_mime isa MIME"application/vnd.pluto.table+object"
        @test notebook.cells[9].repr_mime isa MIME"application/vnd.pluto.table+object"
        @test notebook.cells[11].repr_mime isa MIME"application/vnd.pluto.table+object"
        @test notebook.cells[12].repr_mime isa MIME"application/vnd.pluto.table+object"
        @test notebook.cells[2].output_repr isa Dict
        @test notebook.cells[3].output_repr isa Dict
        @test notebook.cells[4].output_repr isa Dict
        @test notebook.cells[5].output_repr isa Dict
        @test notebook.cells[6].output_repr isa Dict
        @test notebook.cells[7].output_repr isa Dict
        @test notebook.cells[8].output_repr isa Dict
        @test notebook.cells[9].output_repr isa Dict
        @test notebook.cells[11].output_repr isa Dict
        @test notebook.cells[12].output_repr isa Dict

        @test notebook.cells[10].repr_mime isa MIME"text/plain"
        @test notebook.cells[10].errored == false

        # to see if we truncated correctly, we convert the output to string and check how big it is
        # because we don't want to test too specifically
        roughsize(x) = length(string(x))

        smallsize = roughsize(notebook.cells[5])
        manyrowssize = roughsize(notebook.cells[6])
        manycolssize = roughsize(notebook.cells[7])
        @test manyrowssize < 50 * smallsize
        @test manycolssize < 50 * smallsize

        # TODO: test lazy loading more rows/cols

        WorkspaceManager.unmake_workspace((🍭, notebook))
    end

    begin
        escape_me = "16 \\ \" ' / \b \f \n \r \t 💩 \x10 \$"
        notebook = Notebook([
            Cell("a\\"),
            Cell("1 = 2"),

            Cell("b = 3.0\nb = 3"),
            Cell("\n# uhm\n\nc = 4\n\n# wowie \n\n"),
            Cell("d = 5;"),
            Cell("e = 6; f = 6"),
            Cell("g = 7; h = 7;"),
            Cell("\n\n0 + 8; 0 + 8;\n\n\n"),
            Cell("0 + 9; 9;\n\n\n"),
            Cell("0 + 10;\n10;"),
            Cell("0 + 11;\n11"),
        ])
        fakeclient.connected_notebook = notebook

        @testset "Strange code"  begin
            update_run!(🍭, notebook, notebook.cells[1])
            update_run!(🍭, notebook, notebook.cells[2])
            @test notebook.cells[1].errored == true
            @test notebook.cells[2].errored == true

        end
        @testset "Mutliple expressions & semicolon"  begin

            update_run!(🍭, notebook, notebook.cells[3:end])
            @test occursinerror("syntax: extra token after", notebook.cells[3])

            @test notebook.cells[4].errored == false
            @test notebook.cells[4].output_repr == "4"
            @test notebook.cells[4].rootassignee == :c

            @test notebook.cells[5].errored == false
            @test notebook.cells[5].output_repr == ""
            @test notebook.cells[5].rootassignee === nothing

            @test notebook.cells[6].errored == false
            @test notebook.cells[6].output_repr == "6"
            @test notebook.cells[6].rootassignee === nothing

            @test notebook.cells[7].errored == false
            @test notebook.cells[7].output_repr == ""
            @test notebook.cells[7].rootassignee === nothing

            @test notebook.cells[8].errored == false
            @test notebook.cells[8].output_repr == ""

            @test notebook.cells[9].errored == false
            @test notebook.cells[9].output_repr == ""

            @test occursinerror("syntax: extra token after", notebook.cells[10])

            @test occursinerror("syntax: extra token after", notebook.cells[11])
        end

        WorkspaceManager.unmake_workspace((🍭, notebook))
    end

    @testset "Stack traces" begin
        escape_me = "16 \\ \" ' / \b \f \n \r \t 💩 \x10 \$"

        codes = [
            "sqrt(-1)",
            "let\n\nsqrt(-2)\nend",
            "\"Something very exciting!\"\nfunction w(x)\n\tsqrt(x)\nend",
            "w(-4)",
            "error(" * sprint(Base.print_quoted, escape_me) * ")",
        ]

        notebook1 = Notebook([
            Cell(code)
            for (i, code) in enumerate(codes)
        ])

        # create struct to disable the function-generating optimization
        notebook2 = Notebook([
            Cell("struct S$(i) end; $code")
            for (i, code) in enumerate(codes)
        ])

        @testset "$(wrapped ? "With" : "Without") function wrapping" for wrapped in [false, true]
            notebook = wrapped ? notebook1 : notebook2

            fakeclient.connected_notebook = notebook

            @test_nowarn update_run!(🍭, notebook, notebook.cells[1:5])

            @test occursinerror("DomainError", notebook.cells[1])
            let
                st = notebook.cells[1].output_repr
                @test length(st[:stacktrace]) == 4 # check in REPL
                if Pluto.can_insert_filename
                    @test st[:stacktrace][4][:line] == 1
                    @test occursin(notebook.cells[1].cell_id |> string, st[:stacktrace][4][:file])
                    @test occursin(notebook.path |> basename, st[:stacktrace][4][:file])
                else
                    @test_broken false
                end
            end

            @test occursinerror("DomainError", notebook.cells[2])
            let
                st = notebook.cells[2].output_repr
                @test length(st[:stacktrace]) == 4
                if Pluto.can_insert_filename
                    @test st[:stacktrace][4][:line] == 3
                    @test occursin(notebook.cells[2].cell_id |> string, st[:stacktrace][4][:file])
                    @test occursin(notebook.path |> basename, st[:stacktrace][4][:file])
                else
                    @test_broken false
                end
            end

            @test occursinerror("DomainError", notebook.cells[4])
            let
                st = notebook.cells[4].output_repr
                @test length(st[:stacktrace]) == 5

                if Pluto.can_insert_filename
                    @test st[:stacktrace][4][:line] == 3
                    @test occursin(notebook.cells[3].cell_id |> string, st[:stacktrace][4][:file])
                    @test occursin(notebook.path |> basename, st[:stacktrace][4][:file])

                    @test st[:stacktrace][5][:line] == 1
                    @test occursin(notebook.cells[4].cell_id |> string, st[:stacktrace][5][:file])
                    @test occursin(notebook.path |> basename, st[:stacktrace][5][:file])
                else
                    @test_broken false
                end
            end

            let
                st = notebook.cells[5].output_repr
                @test occursin(escape_me, st[:msg])
            end

            WorkspaceManager.unmake_workspace((🍭, notebook))
        end

    end

end


text:
module Diracoperators
    import ..Gaugefields:GaugeFields,GaugeFields_1d
    import ..Fermionfields:FermionFields,WilsonFermion,Wx!,Wdagx!,WdagWx!,
                            StaggeredFermion,set_wing_fermi!,Dx!,add!,clear!,Dxplus!
    import ..Clover:Make_CloverFμν
    import ..Actions:FermiActionParam_Wilson,FermiActionParam_WilsonClover,
                FermiActionParam_Staggered
    using LinearAlgebra

    abstract type Dirac_operator
    end

    abstract type Adjoint_Dirac_operator
    end

    abstract type DdagD_operator
    end

    function Dirac_operator(U::Array{T,1},x,fparam) where  T <: GaugeFields
        if typeof(fparam) == FermiActionParam_Wilson
            W = Wilson_operator(U,x,fparam)
        elseif typeof(fparam) == FermiActionParam_WilsonClover
            W = WilsonClover_operator(U,x,fparam)
        elseif typeof(fparam) == FermiActionParam_Staggered
            W = Staggered_operator(U,x,fparam)
        else
            error("unknown fparam: fparam is ",typeof(fparam))
        end
    end

    struct Wilson_operator{T} <: Dirac_operator  where  T <: GaugeFields
        U::Array{T,1}
        _temporal_fermi::Array{WilsonFermion,1}

        function Wilson_operator(U::Array{T,1},x,fparam) where  T <: GaugeFields
            num = 5
            _temporal_fermi = Array{WilsonFermion,1}(undef,num)
            for i=1:num
                _temporal_fermi[i] = similar(x)
            end
            return new{eltype(U)}(U,_temporal_fermi)
        end
    end

    struct DdagD_Wilson_operator <: DdagD_operator
        dirac::Wilson_operator
        function DdagD_Wilson_operator(U::Array{T,1},x,fparam) where  T <: GaugeFields
            return new(Wilson_operator(U,x,fparam))
        end
    end

    struct WilsonClover_operator{T} <: Dirac_operator  where  T <: GaugeFields
        U::Array{T,1}
        _temporal_fermi::Array{WilsonFermion,1}
        CloverFμν::Array{ComplexF64,4}

        function WilsonClover_operator(U::Array{T,1},x,fparam) where  T <: GaugeFields
            num = 5
            _temporal_fermi = Array{WilsonFermion,1}(undef,num)
            for i=1:num
                _temporal_fermi[i] = similar(x)
            end

            NV = U[1].NV
            NX = U[1].NX
            NY = U[1].NY
            NZ = U[1].NZ
            NT = U[1].NT
            NC = U[1].NC

            temp1 = GaugeFields_1d(NC,NX,NY,NZ,NT)
            temps = Array{typeof(temp1),1}(undef,num)
            for i=1:num
                temps[i] = deepcopy(temp1)
            end

            CloverFμν = Make_CloverFμν(fparam,U,temps)

            return new{eltype(U)}(U,_temporal_fermi,CloverFμν)
        end
    end

    struct DdagD_WilsonClover_operator <: DdagD_operator
        dirac::WilsonClover_operator
        function DdagD_WilsonClover_operator(U::Array{T,1},x,fparam) where  T <: GaugeFields
            return new(WilsonClover_operator(U,x,fparam))
        end
    end

    struct Staggered_operator{T} <: Dirac_operator  where  T <: GaugeFields
        U::Array{T,1}
        _temporal_fermi::Array{StaggeredFermion,1}
        mass::Float64

        function Staggered_operator(U::Array{T,1},x,fparam) where  T <: GaugeFields
            num = 9
            _temporal_fermi = Array{StaggeredFermion,1}(undef,num)
            for i=1:num
                _temporal_fermi[i] = similar(x)
            end
            return new{eltype(U)}(U,_temporal_fermi,fparam.mass)
        end
    end

    struct DdagD_Staggered_operator{perturb} <: DdagD_operator
        dirac::Staggered_operator
        t::Float64
        function DdagD_Staggered_operator(U::Array{T,1},x,fparam;t=1) where  T <: GaugeFields
            if t == 1
                perturb = false
            else
                perturb = true
            end

            return new{perturb}(Staggered_operator(U,x,fparam),t)
        end

        function DdagD_Staggered_operator(A::DdagD_Staggered_operator,t)
            if t == 1
                perturb = false
            else
                perturb = true
            end
            return new{perturb}(A.dirac,t)
        end
    end

    struct DdagDND_Staggered_operator <: DdagD_operator #nondiagonal part of DdagD
        dirac::Staggered_operator

        function DdagDND_Staggered_operator(A::DdagD_Staggered_operator)
            return new(A.dirac)
        end
    end

    function Base.size(A::Staggered_operator)
        NX = A.U[1].NX
        NY = A.U[1].NY
        NZ = A.U[1].NZ
        NT = A.U[1].NT
        NC = A.U[1].NC
        return NX*NY*NZ*NT*NC,NX*NY*NZ*NT*NC
    end

    function Base.size(A::DdagD_operator)
        return size(A.dirac)
    end

    function DdagD_operator(U::Array{T,1},x,fparam) where  T <: GaugeFields
        if typeof(fparam) == FermiActionParam_Wilson
            W = DdagD_Wilson_operator(U,x,fparam)
        elseif typeof(fparam) == FermiActionParam_WilsonClover
            W = DdagD_WilsonClover_operator(U,x,fparam)
        elseif typeof(fparam) == FermiActionParam_Staggered
            W = DdagD_Staggered_operator(U,x,fparam)
        else
            error("unknown fparam: fparam is ",typeof(fparam))
        end
    end

    struct Adjoint_Wilson_operator <: Adjoint_Dirac_operator
        parent::Wilson_operator
    end

    struct Adjoint_WilsonClover_operator <: Adjoint_Dirac_operator
        parent::WilsonClover_operator
    end

    struct Adjoint_Staggered_operator <: Adjoint_Dirac_operator
        parent::Staggered_operator
    end

    function Base.adjoint(A::Wilson_operator)
        Adjoint_Wilson_operator(A)
    end

    function Base.adjoint(A::WilsonClover_operator)
        Adjoint_WilsonClover_operator(A)
    end

    function Base.adjoint(A::Staggered_operator)
        Adjoint_Staggered_operator(A)
    end

    Base.adjoint(A::Adjoint_Dirac_operator) = A.parent

    function LinearAlgebra.mul!(y::WilsonFermion,A::Wilson_operator,x::WilsonFermion) #y = A*x
        Wx!(y,A.U,x,A._temporal_fermi)
        return
    end

    function LinearAlgebra.mul!(y::WilsonFermion,A::WilsonClover_operator,x::WilsonFermion) #y = A*x
        Wx!(y,A.U,x,A._temporal_fermi,A.CloverFμν)
        return
    end

    function LinearAlgebra.mul!(y::FermionFields,A::T,x::FermionFields)  where T <: DdagD_operator #y = A*x
        temp = A.dirac._temporal_fermi[5]
        mul!(temp,A.dirac,x)
        mul!(y,A.dirac',temp)

        return
    end

    function LinearAlgebra.mul!(y::FermionFields,A::DdagD_Staggered_operator{perturb},x::FermionFields)  where perturb
        temps = A.dirac._temporal_fermi
        temp = temps[5]
        temp2 = temps[6]
        Dx!(temp,A.dirac.U,x,temps)
        Dx!(temp2,A.dirac.U,temp,temps)
        clear!(y)

        if perturb
            add!(y,A.dirac.mass^2,x,-A.t,temp2)

            for ν=1:4
                Dxplus!(temp,ν,A.dirac.U,x,temps)
                Dxplus!(temp2,ν,A.dirac.U,temp,temps)

                add!(y,1,x,-(1-A.t),temp2)
            end

            set_wing_fermi!(y)
        else

            add!(y,A.dirac.mass^2,x,-1,temp2)
            set_wing_fermi!(y)
        end

        return
    end

    function LinearAlgebra.mul!(y::FermionFields,A::DdagDND_Staggered_operator,x::FermionFields)  where perturb
        temps = A.dirac._temporal_fermi
        temp = temps[5]
        temp2 = temps[6]
        Dx!(temp,A.dirac.U,x,temps)
        Dx!(temp2,A.dirac.U,temp,temps)
        clear!(y)
        add!(y,1,x,-1,temp2)

        for ν=1:4
            Dxplus!(temp,ν,A.dirac.U,x,temps)
            Dxplus!(temp2,ν,A.dirac.U,temp,temps)

            add!(y,1,x,1,temp2)
        end
        set_wing_fermi!(y)

        return
    end

    function LinearAlgebra.mul!(y::FermionFields,A::DdagD_operator,x::FermionFields,indices) #y = A*x
        mul!(y,A,x)
        return
    end

    function LinearAlgebra.mul!(y::FermionFields,A::DdagD_Staggered_operator,x::FermionFields,indices) #y = A*x
        WdagWx!(y,A.dirac.U,x,A.dirac._temporal_fermi,A.dirac.mass,indices)
        return
    end

    function LinearAlgebra.mul!(y::StaggeredFermion,A::Staggered_operator,x::StaggeredFermion) #y = A*x
        temps = A._temporal_fermi
        temp = temps[4]
        Dx!(temp,A.U,x,temps)
        clear!(y)
        add!(y,A.mass,x,1,temp)
        set_wing_fermi!(y)
        return
    end

    function Base.:*(A::Dirac_operator,x::FermionFields)
        y = similar(x)
        mul!(y,A,x)
        return y
    end

    function Base.:*(A::DdagD_operator,x::FermionFields)
        y = similar(x)
        mul!(y,A,x)
        return y
    end

    function LinearAlgebra.mul!(y::WilsonFermion,A::Adjoint_Wilson_operator,x::WilsonFermion) #y = A*x
        Wdagx!(y,A.parent.U,x,A.parent._temporal_fermi)
    end

    function LinearAlgebra.mul!(y::WilsonFermion,A::Adjoint_WilsonClover_operator,x::WilsonFermion) #y = A*x
        Wdagx!(y,A.parent.U,x,A.parent._temporal_fermi,A.parent.CloverFμν)
    end

    function LinearAlgebra.mul!(y::StaggeredFermion,A::Adjoint_Staggered_operator,x::StaggeredFermion) #y = A*x
        temps = A.parent._temporal_fermi
        temp = temps[4]
        Dx!(temp,A.parent.U,x,temps)
        clear!(y)
        add!(y,A.parent.mass,x,-1,temp)
        set_wing_fermi!(y)
        return
    end

    function Base.:*(A::Adjoint_Dirac_operator,x::WilsonFermion)
        y = similar(x)
        mul!(y,A,x)
        return y
    end

end


text:
doc"""
    Arcsine(a,b)

The *Arcsine distribution* has probability density function

$f(x) = \frac{1}{\pi \sqrt{(x - a) (b - x)}}, \quad x \in [a, b]$

```julia
Arcsine()        # Arcsine distribution with support [0, 1]
Arcsine(b)       # Arcsine distribution with support [0, b]
Arcsine(a, b)    # Arcsine distribution with support [a, b]

params(d)        # Get the parameters, i.e. (a, b)
minimum(d)       # Get the lower bound, i.e. a
maximum(d)       # Get the upper bound, i.e. b
location(d)      # Get the left bound, i.e. a
scale(d)         # Get the span of the support, i.e. b - a
```

External links

* [Arcsine distribution on Wikipedia](http://en.wikipedia.org/wiki/Arcsine_distribution)

"""

immutable Arcsine{T<:Real} <: ContinuousUnivariateDistribution
    a::T
    b::T

    Arcsine(a::T, b::T) = (@check_args(Arcsine, a < b); new(a, b))
end

Arcsine{T<:Real}(a::T, b::T) = Arcsine{T}(a, b)
Arcsine(a::Real, b::Real) = Arcsine(promote(a, b)...)
Arcsine(a::Integer, b::Integer) = Arcsine(Float64(a), Float64(b))
Arcsine(b::Real) = Arcsine(0.0, b)
Arcsine() = Arcsine(0.0, 1.0)

@distr_support Arcsine d.a d.b

#### Conversions
function convert{T<:Real}(::Type{Arcsine{T}}, a::Real, b::Real)
    Arcsine(T(a), T(b))
end
function convert{T <: Real, S <: Real}(::Type{Arcsine{T}}, d::Arcsine{S})
    Arcsine(T(d.a), T(d.b))
end

### Parameters

location(d::Arcsine) = d.a
scale(d::Arcsine) = d.b - d.a
params(d::Arcsine) = (d.a, d.b)
@inline partype{T<:Real}(d::Arcsine{T}) = T

### Statistics

mean(d::Arcsine) = (d.a + d.b) / 2
median(d::Arcsine) = mean(d)
mode(d::Arcsine) = d.a
modes(d::Arcsine) = [d.a, d.b]

var(d::Arcsine) = abs2(d.b - d.a) / 8
skewness{T<:Real}(d::Arcsine{T}) = zero(T)
kurtosis{T<:Real}(d::Arcsine{T}) = -T(3/2)

entropy(d::Arcsine) = -0.24156447527049044469 + log(scale(d))

### Evaluation

function pdf(d::Arcsine, x::Real)
    insupport(d, x) ? one(d.a) / (π * sqrt((x - d.a) * (d.b - x))) : zero(d.a)
end

function logpdf{T<:Real}(d::Arcsine{T}, x::Real)
    insupport(d, x) ? -(logπ + log((x - d.a) * (d.b - x))/2) : -T(Inf)
end

cdf{T<:Real}(d::Arcsine{T}, x::Real) = x < d.a ? zero(T) :
                              x > d.b ? one(T) :
                              0.636619772367581343 * asin(sqrt((x - d.a) / (d.b - d.a)))

quantile(d::Arcsine, p::Real) = location(d) + abs2(sin(halfπ * p)) * scale(d)

### Sampling

rand(d::Arcsine) = quantile(d, rand())


text:
using Pitaya
using Test

const dir = joinpath(dirname(pathof(Pitaya)), "..", "test")

@testset "Pitaya module" begin

    for f in ["test_utils.jl",
              "test_network.jl",
              "test_api.jl"]
        file = joinpath(dir, f)
        println("Running $f tests...")
        if isfile(file)
            include(file)
        else
            @show readdir(dirname(file))
        end
    end
end

max_stars_repo_name: naustica/Pitaya

max_issues_repo_name: naustica/Pitaya

max_forks_repo_name: naustica/Pitaya