text:
const XDMF_TRIANGLE = 4
const XDMF_QUAD = 5
const XDMF_QUADRATIC_TRIANGLE = 36
const XDMF_QUADRATIC_QUAD = 37
const XDMF_TETRA = 6
const XDMF_HEXAHEDRON = 9
const XDMF_QUADRATIC_TETRA = 38
const XDMF_QUADRATIC_HEXAHEDRON = 48

function vtk2xdmf(vtk_type)
    if vtk_type == VTK_TRIANGLE
        return XDMF_TRIANGLE
    elseif vtk_type == VTK_QUAD
        return XDMF_QUAD
    elseif vtk_type == VTK_QUADRATIC_TRIANGLE
        return XDMF_QUADRATIC_TRIANGLE
    elseif vtk_type == VTK_QUADRATIC_QUAD
        return XDMF_QUADRATIC_QUAD
    elseif vtk_type == VTK_TETRA
        return XDMF_TETRA
    elseif vtk_type == VTK_HEXAHEDRON
        return XDMF_HEXAHEDRON
    elseif vtk_type == VTK_QUADRATIC_TETRA
        return XDMF_QUADRATIC_TETRA
    elseif vtk_type == VTK_QUADRATIC_HEXAHEDRON
        return XDMF_QUADRATIC_HEXAHEDRON
    else
        error("Invalid VTK type.")
        return nothing
    end
end
#################################################################################
#                                    WRITE
#################################################################################
function write_xdmf(filename::String, mpt::MeshPartitionTree)
    # Check valid filename
    if !endswith(filename, ".xdmf")
        error("Invalid filename.")
    end

    material_groups = String[]
    for mesh in leaf_meshes(mpt)
        for group_name in keys(mesh.groups)
            if startswith(group_name, "Material") && group_name ∉ material_groups
                push!(material_groups, group_name)
            end
        end
    end
    sort!(material_groups)
    # h5 filename
    h5_filename = filename[1:end-4]*"h5"
    h5_file = h5open(h5_filename, "w")
    # XML
    xdoc = XMLDocument()
    try
        # Xdmf
        xroot = ElementNode("Xdmf")
        setroot!(xdoc, xroot)
        link!(xroot, AttributeNode("Version", "3.0"))
        # Domain
        xdomain = ElementNode("Domain")
        link!(xroot, xdomain)
        # Material names
        materials = [ mat[11:end] for mat in material_groups ]
        if materials != String[]
            xmaterials = ElementNode("Information")
            link!(xdomain, xmaterials)
            link!(xmaterials, AttributeNode("Name", "Materials"))
            link!(xmaterials, TextNode(join(materials, " ")))
        end
        _add_mesh_partition_xdmf!(xdomain, h5_filename, h5_file,
                                  mpt.partition_tree, mpt.leaf_meshes, material_groups)
    finally
        write(filename, xdoc)
        close(h5_file)
    end
    return nothing
end

function write_xdmf(filename::String, mesh::AbstractMesh)
    # Check valid filename
    if !endswith(filename, ".xdmf")
        error("Invalid filename.")
    end

    # h5
    h5_filename = filename[1:end-4]*"h5"
    h5_file = h5open(h5_filename, "w")

    # XML
    xdoc = XMLDocument()
    try
        # Xdmf
        xroot = ElementNode("Xdmf")
        setroot!(xdoc, xroot)
        link!(xroot, AttributeNode("Version", "3.0"))
        # Domain
        xdomain = ElementNode("Domain")
        link!(xroot, xdomain)

        material_groups = sort!(filter!(x->startswith(x, "Material"),
                                        collect(keys(mesh.groups))
                                       )
                               )
        materials = [ mat[11:end] for mat in material_groups ]

        # Material names
        if materials != String[]
            xmaterials = ElementNode("Information")
            link!(xdomain, xmaterials)
            link!(xmaterials, AttributeNode("Name", "Materials"))
            link!(xmaterials, TextNode(join(materials, " ")))
        end

    # Add uniform grid
    _add_uniform_grid_xdmf!(xdomain, h5_filename, h5_file, mesh, material_groups)

    finally
        write(filename, xdoc)
        close(h5_file)
    end
    return nothing
end

# Helper functions for write_xdmf
# -------------------------------------------------------------------------------------------------
function _add_uniform_grid_xdmf!(xml::EzXML.Node,
                                h5_filename::String,
                                h5_mesh::Union{HDF5.Group, HDF5.File},
                                mesh::AbstractMesh,
                                materials::Vector{String})
    # Grid
    xgrid = ElementNode("Grid")
    link!(xml, xgrid)
    link!(xgrid, AttributeNode("Name", mesh.name))
    link!(xgrid, AttributeNode("GridType", "Uniform"))

    # h5
    h5_group = create_group(h5_mesh, mesh.name)

    # Geometry
    _write_xdmf_geometry!(xgrid, h5_filename, h5_group, mesh)

    # Topology
    _write_xdmf_topology!(xgrid, h5_filename, h5_group, mesh)

    # Non-material groups
    if mesh.groups != Dict{String,BitSet}()
        _write_xdmf_groups!(xgrid, h5_filename, h5_group, mesh)
    end

    # Materials
    if materials != String[]
        _write_xdmf_materials!(xgrid, h5_filename, h5_group, mesh, materials)
    end
    return nothing
end

function _write_xdmf_geometry!(xml::EzXML.Node,
                              h5_filename::String,
                              h5_mesh::HDF5.Group,
                              mesh::AbstractMesh)
    verts = points(mesh)
    float_type = typeof(verts[1][1])
    float_precision = string(sizeof(float_type))
    nverts = length(verts)
    nverts_str = string(nverts)
    point_dim = length(verts[1])
    if point_dim == 2
        XYZ = "XY"
        dim = " 2"
    else
        XYZ= "XYZ"
        dim = " 3"
    end
    # Convert the points into an array
    vert_array = zeros(float_type, point_dim, nverts)
    for i in eachindex(verts)
        vert_array[:, i] = coordinates(verts[i])
    end
    # Geometry
    xgeom = ElementNode("Geometry")
    link!(xml, xgeom)
    link!(xgeom, AttributeNode("GeometryType", XYZ))
    # DataItem
    xdataitem = ElementNode("DataItem")
    link!(xgeom, xdataitem)
    link!(xdataitem, AttributeNode("DataType", "Float"))
    link!(xdataitem, AttributeNode("Dimensions", nverts_str*dim))
    link!(xdataitem, AttributeNode("Format", "HDF"))
    link!(xdataitem, AttributeNode("Precision", float_precision))
    h5_text_item = split(string(h5_mesh))[2]
    link!(xdataitem, TextNode(string(h5_filename, ":", h5_text_item, "/points")))
    # Write the h5
    h5_mesh["points", compress = 3] = vert_array
    return nothing
end

_nelements_xdmf(mesh::VolumeMesh) = length(mesh.types)
_nelements_xdmf(mesh::PolytopeVertexMesh) = length(mesh.polytopes)

_topo_length_xdmf(mesh::VolumeMesh) = length(mesh.types) + length(mesh.connectivity)
function _topo_length_xdmf(mesh::PolytopeVertexMesh)
    return mapreduce(x->length(vertices(x)), +, polytopes(mesh)) + length(mesh.polytopes)
end

_uint_type_xdmf(mesh::VolumeMesh{Dim,T,U}) where {Dim,T,U} = U
_uint_type_xdmf(mesh::PolytopeVertexMesh) = typeof(mesh.polytopes[1][1])

function _populate_topo_array_xdmf!(topo_array, mesh::VolumeMesh)
    topo_ctr = 1
    for i in eachindex(mesh.types)
        topo_array[topo_ctr] = vtk2xdmf(mesh.types[i])
        topo_ctr += 1
        len_element = points_in_vtk_type(mesh.types[i])
        offset = mesh.offsets[i]
        # adjust 1-based to 0-based indexing
        @. topo_array[topo_ctr:topo_ctr + len_element - 1] =
            mesh.connectivity[offset:offset + len_element - 1] - 1
        topo_ctr += len_element
    end
    return nothing
end

function _populate_topo_array_xdmf!(topo_array, mesh::PolytopeVertexMesh)
    ptopes = mesh.polytopes
    topo_ctr = 1
    for i in eachindex(ptopes)
        p = ptopes[i]
        topo_array[topo_ctr] = vtk2xdmf(vtk_type(typeof(p)))
        topo_ctr += 1
        len_element = length(vertices(p))
        # adjust 1-based to 0-based indexing
        topo_array[topo_ctr:topo_ctr + len_element - 1] .= vertices(p) .- 1
        topo_ctr += len_element
    end
    return nothing
end

function _write_xdmf_topology!(xml::EzXML.Node,
                              h5_filename::String,
                              h5_mesh::HDF5.Group,
                              mesh::AbstractMesh)
    nelements = _nelements_xdmf(mesh)
    nelements_str = string(nelements)
    topo_length = _topo_length_xdmf(mesh)
    uint_type = _uint_type_xdmf(mesh)
    uint_precision = string(sizeof(uint_type))
    topo_array = Vector{uint_type}(undef, topo_length)
    _populate_topo_array_xdmf!(topo_array, mesh)
    ndims = string(length(topo_array))
    # Topology
    xtopo = ElementNode("Topology")
    link!(xml, xtopo)
    link!(xtopo, AttributeNode("TopologyType", "Mixed"))
    link!(xtopo, AttributeNode("NumberOfElements", nelements_str))
    # DataItem
    xdataitem = ElementNode("DataItem")
    link!(xtopo, xdataitem)
    link!(xdataitem, AttributeNode("DataType", "UInt"))
    link!(xdataitem, AttributeNode("Dimensions", ndims))
    link!(xdataitem, AttributeNode("Format", "HDF"))
    link!(xdataitem, AttributeNode("Precision", uint_precision))
    h5_text_item = split(string(h5_mesh))[2]
    link!(xdataitem, TextNode(string(h5_filename, ":", h5_text_item, "/connectivity")))
    # Write the h5
    h5_mesh["connectivity", compress = 3] = topo_array
    return nothing
end

function _write_xdmf_materials!(xml::EzXML.Node,
                               h5_filename::String,
                               h5_mesh::HDF5.Group,
                               mesh::AbstractMesh,
                               materials::Vector{String})
    nelements = _nelements_xdmf(mesh)
    mat_id_array = fill(-1, nelements)
    mesh_group_names = keys(mesh.groups)
    for (mat_id, material_name) in enumerate(materials)
        if material_name in mesh_group_names
            for cell in mesh.groups[material_name]
                if mat_id_array[cell] === -1
                    mat_id_array[cell] = mat_id - 1 # 0 based indexing
                else
                    error("Mesh cell "*string(cell)*" has multiple materials assigned to it.")
                end
            end
        end
    end
    if any(x->x === -1, mat_id_array)
        error("Some mesh cells do not have a material.")
    end
    N = length(materials)
    uint_type = _select_uint_type(N)
    uint_precision= string(sizeof(uint_type))

    # MaterialID
    xmaterial = ElementNode("Attribute")
    link!(xml, xmaterial)
    link!(xmaterial, AttributeNode("Center", "Cell"))
    link!(xmaterial, AttributeNode("Name", "Material"))
    # DataItem
    xdataitem = ElementNode("DataItem")
    link!(xmaterial, xdataitem)
    link!(xdataitem, AttributeNode("DataType", "UInt"))
    link!(xdataitem, AttributeNode("Dimensions", string(nelements)))
    link!(xdataitem, AttributeNode("Format", "HDF"))
    link!(xdataitem, AttributeNode("Precision", uint_precision))
    h5_text_item = split(string(h5_mesh))[2]
    link!(xdataitem, TextNode(string(h5_filename, ":", h5_text_item, "/material")))
    # Write the h5
    h5_mesh["material", compress = 3] = mat_id_array
    return nothing
end

function _write_xdmf_groups!(xml::EzXML.Node,
                             h5_filename::String,
                             h5_mesh::HDF5.Group,
                             mesh::AbstractMesh)
    for set_name in keys(mesh.groups)
        if startswith(set_name, "Material")
            continue
        end
        grp = mesh.groups[set_name]
        N = maximum(grp)
        uint_type = _select_uint_type(N)
        uint_precision= string(sizeof(uint_type))

        id_array = collect(grp) .- 1
        # Set
        xset = ElementNode("Set")
        link!(xml, xset)
        link!(xset, AttributeNode("Name", set_name))
        link!(xset, AttributeNode("SetType", "Cell"))
        # DataItem
        xdataitem = ElementNode("DataItem")
        link!(xset, xdataitem)
        link!(xdataitem, AttributeNode("DataType", "UInt"))
        nelements = string(length(grp))
        link!(xdataitem, AttributeNode("Dimensions", nelements))
        link!(xdataitem, AttributeNode("Format", "HDF"))
        link!(xdataitem, AttributeNode("Precision", uint_precision))
        h5_text_item = split(string(h5_mesh))[2]
        link!(xdataitem, TextNode(string(h5_filename, ":", h5_text_item, "/", set_name)))
        # Write the h5
        h5_mesh[set_name, compress = 3] = id_array
    end
    return nothing
end

function _add_mesh_partition_xdmf!(xml::EzXML.Node,
                                   h5_filename::String,
                                   h5_mesh::Union{HDF5.Group, HDF5.File},
                                   node::Tree{Tuple{Int64, String}},
                                   leaf_meshes::Vector{<:AbstractMesh},
                                   materials::Vector{String})
    id = node.data[1]
    if id === 0 # Internal node
        name = node.data[2]
        # Grid
        xgrid = ElementNode("Grid")
        link!(xml, xgrid)
        link!(xgrid, AttributeNode("Name", name))
        link!(xgrid, AttributeNode("GridType", "Tree"))
        # h5_group
        h5_group = create_group(h5_mesh, name)
        for child in node.children
            _add_mesh_partition_xdmf!(xgrid, h5_filename, h5_group, child,
                                      leaf_meshes, materials)
        end
    else # Leaf node
        _add_uniform_grid_xdmf!(xml, h5_filename, h5_mesh,
                                leaf_meshes[id], materials)
    end
    return nothing
end

#################################################################################
#                                    READ
#################################################################################
#xdmf_read_error(x::String) = error("Error reading XDMF file.")
#function read_xdmf(path::String, ::Type{T}) where {T<:AbstractFloat}
#    xdoc = readxml(path)
#    xroot = root(xdoc)
#    nodename(xroot) != "Xdmf" && xdmf_read_error()
#    try
#        version = xroot["Version"]
#        version != "3.0" && xdmf_read_error()
#        xdomain = firstnode(xroot)
#        nodename(xdomain) != "Domain" && xdmf_read_error()
#        material_names = String[]
#        if 1 < countnodes(xdomain) && nodename(firstnode(xdomain)) == "Information"
#            append!(material_names,
#        end
#
#    finally
#
#    end
#end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
module ProxBTF

export pbtfProblem,pbtf,pbsrtfProblem,pbsrtf,visualize

using SparseArrays: SparseMatrixCSC,spdiagm
using Roots: find_zero, Bisection
using Lasso: fit, FusedLasso
using Gurobi: Env, GRBsetdblparam, Optimizer
using Convex: Variable, minimize, sumsquares, norm, solve!, evaluate
import LogDensityProblems.dimension,LogDensityProblems.capabilities
import LogDensityProblems.logdensity_and_gradient,LogDensityProblems.LogDensityOrder
using Parameters: @unpack
using Random: GLOBAL_RNG
using DynamicHMC: mcmc_with_warmup, default_warmup_stages, ProgressMeterReport, NoProgressReport
using DynamicHMC.Diagnostics: summarize_tree_statistics
using DataFrames: DataFrame
using MCMCChains: Chains,summarize,quantile
using Statistics: mean,norm,var
using StatsPlots: plot, scatter!, plot!

include("utils.jl")
include("projectors.jl")
include("model_pbtf.jl")
include("model_pbsrtf.jl")
include("mcmc.jl")
include("visualization.jl")
end # module

ext: jl


text:
# transform (x = :x, t = 2 * :x - :y)
#
# gormap(r -> (x = r.x, t = 2 * r.x + r.y))
#
# @transform rows body
# gormap(rows, row -> replace_vars(body))

# mutate begin ...  (x = :x, t = t) end
#
# gormap(r -> merge(r, (row -> ...)(r)))
#

# allow sequential assignment
#
# i.e. we might have
# begin
#  x = :y + 2
#  x = :x
#  t = 2x + :y
#  (t,)
# end
#

#
# 1. traverse body and replace all quoted symbols with gensyms
#    keeping track of replacements
#

"look up or define variable ´var´ in environment ´vars´"
function lookup(vars, var)
    if !haskey(vars, var)
	vars[var] = gensym()
    end
    vars[var]
end

function replace_vars(ex, vars) end

# default
replace_vars(ex, vars) = ex

# treat QuoteNode as regular quote
replace_vars(ex::QuoteNode, vars) = replace_vars(Meta.quot(ex.value), vars)

function replace_vars(ex::Expr, vars)
    if ex.head == :quote
	lookup(vars, Meta.quot(ex.args[1]))
    else
	Expr(ex.head, map(arg -> replace_vars(arg, vars), ex.args)...)
    end
end

#
# 2. create declaration block that assigns all vars from
#

# then we would first unpack row (but macro does not know which columns exist
#
# _x = row[:x]
# _y = row[:y]
#
# then execute body with transformed vars
# _x = _y + 2
# x = _x
# _t = 2x + _y
#
#
# then return merge(row, (:x = _x, :y = _y, :t = _t))

# default
replace_vars(ex) = ex

# treat QuoteNode as regular quote
replace_vars(ex::QuoteNode) = replace_vars(Meta.quot(ex.value))

function replace_vars(ex::Expr)
    if ex.head == :quote
        :(row[$(Meta.quot(ex.args[1]))])
    else
	Expr(ex.head, map(arg -> replace_vars(arg), ex.args)...)
    end
end


text:
f = open("jawiki-country.txt", "r")

for line in readlines(f)
  category_reg = r"\[\[Category:"
  if ismatch(category_reg, line)
    print(replace(line, r"\[\[Category:|(\|\*)*\]\].*", ""))
  end
end


text:
module Kinematics

using LinearAlgebra
using StaticArrays, ComponentArrays
using UnPack

using Flight.Plotting
using Flight.Attitude
using Flight.Geodesy

import Flight.Modeling: init_x0
import Flight.Plotting: plots

export AbstractKinematics, KinLTF, KinECEF
export VelX, PosData, VelData, KinData, KinInit
export init!, f_kin!, renormalize!

abstract type AbstractKinematics end

Base.@kwdef struct KinInit
    ω_lb_b::SVector{3, Float64} = zeros(SVector{3})
    v_eOb_b::SVector{3, Float64} = zeros(SVector{3})
    q_nb::RQuat = RQuat()
    Ob::Geographic{NVector,Ellipsoidal} = Geographic()
    Δx::Float64 = 0.0
    Δy::Float64 = 0.0
end

struct PosData
    q_nb::RQuat
    q_eb::RQuat
    e_nb::REuler
    q_en::RQuat
    n_e::NVector
    ϕ_λ::LatLon
    h_e::Altitude{Ellipsoidal}
    h_o::Altitude{Orthometric}
    Δxy::SVector{2,Float64}
end

struct VelData
    ω_eb_b::SVector{3,Float64}
    ω_el_n::SVector{3,Float64}
    ω_lb_b::SVector{3,Float64}
    ω_ie_b::SVector{3,Float64}
    ω_ib_b::SVector{3,Float64}
    v_eOb_b::SVector{3,Float64}
    v_eOb_n::SVector{3,Float64}
end

struct KinData
    pos::PosData
    vel::VelData
end

function KinData(init::KinInit)

    @unpack q_nb, Ob, ω_lb_b, v_eOb_b, Δx, Δy = init

    h_e = Ob.alt
    q_el = ltf(Ob)

    #arbitrarily initialize ψ_nl to 0
    q_nl = RQuat()
    q_lb = q_nb
    q_en = q_el
    q_eb = q_el ∘ q_lb

    n_e = NVector(q_el)

    v_eOb_n = q_nb(v_eOb_b)

    (R_N, R_E) = radii(Ob)
    ω_el_n = SVector{3}(
        v_eOb_n[2] / (R_E + Float64(h_e)),
        -v_eOb_n[1] / (R_N + Float64(h_e)),
        0.0)
    ω_el_b = q_nb'(ω_el_n)
    ω_eb_b = ω_el_b + ω_lb_b

    ω_el_l = q_nl'(ω_el_n)
    ω_el_b = q_lb'(ω_el_l)
    ω_lb_b = ω_eb_b - ω_el_b

    ω_ie_e = SVector{3}(0, 0, ω_ie)
    ω_ie_b = q_eb'(ω_ie_e)
    ω_ib_b = ω_ie_b + ω_eb_b

    pos = PosData(q_nb, q_eb, REuler(q_nb), q_en, n_e, LatLon(n_e), h_e,
        Altitude{Orthometric}(h_e, n_e), SVector{2}(Δx, Δy))

    vel = VelData(ω_eb_b, ω_el_n, ω_lb_b, ω_ie_b, ω_ib_b, v_eOb_b, v_eOb_n)

    return KinData(pos, vel)

end

KinData() = KinData(KinInit())

const VelXTemplate = ComponentVector(ω_eb_b = zeros(3), v_eOb_b = zeros(3))
const VelX{T, D} = ComponentVector{T, D, typeof(getaxes(VelXTemplate))} where {T, D}

function renormalize_block!(x, ε) #returns true if norm restored, false otherwise
    norm_x = norm(x)
    abs(norm_x - 1.0) > ε ? (x ./= norm_x; return true) : return false
end

######################### LTF Kinematics #########################

struct KinLTF <: AbstractKinematics end

const PosLTFXTemplate = ComponentVector(q_lb = zeros(4), q_el = zeros(4), Δx = 0.0, Δy = 0.0, h_e = 0.0)
const KinLTFXTemplate = ComponentVector(pos = PosLTFXTemplate, vel = VelXTemplate)
const PosLTFX{T, D} = ComponentVector{T, D, typeof(getaxes(PosLTFXTemplate))} where {T, D}
const KinLTFX{T, D} = ComponentVector{T, D, typeof(getaxes(KinLTFXTemplate))} where {T, D}

init_x0(::KinLTF, init::KinInit = KinInit()) = (x=similar(KinLTFXTemplate); init!(x, init); return x)

function init!(x::KinLTFX, init::KinInit)

    @unpack q_nb, Ob, ω_lb_b, v_eOb_b, Δx, Δy = init

    h_e = Ob.alt
    (R_N, R_E) = radii(Ob)
    v_eOb_n = q_nb * v_eOb_b
    ω_el_n = SVector{3}(
        v_eOb_n[2] / (R_E + Float64(h_e)),
        -v_eOb_n[1] / (R_N + Float64(h_e)),
        0.0)

    ω_el_b = q_nb'(ω_el_n)
    ω_eb_b = ω_el_b + ω_lb_b

    q_lb = q_nb #arbitrarily initialize ψ_nl to 1

    x.pos.q_lb .= q_lb[:]
    x.pos.q_el .= ltf(Ob)[:]
    x.pos.Δx = Δx
    x.pos.Δy = Δy
    x.pos.h_e = h_e
    x.vel.ω_eb_b .= ω_eb_b
    x.vel.v_eOb_b .= v_eOb_b

    return nothing

end

function f_kin!(ẋ_pos::PosLTFX, x::KinLTFX)

    q_lb = RQuat(x.pos.q_lb, normalization = false)
    q_el = RQuat(x.pos.q_el, normalization = false)
    ω_eb_b = SVector{3}(x.vel.ω_eb_b)
    v_eOb_b = SVector{3}(x.vel.v_eOb_b)
    h_e = Altitude{Ellipsoidal}(x.pos.h_e[1])

    n_e = NVector(q_el)
    ψ_nl = get_ψ_nl(q_el)
    q_nl = Rz(ψ_nl)
    q_nb = q_nl ∘ q_lb
    q_eb = q_el ∘ q_lb
    q_en = q_el ∘ q_nl'

    (R_N, R_E) = radii(n_e)
    v_eOb_n = q_nb(v_eOb_b)
    ω_el_n = SVector{3}(
        v_eOb_n[2] / (R_E + Float64(h_e)),
        -v_eOb_n[1] / (R_N + Float64(h_e)),
        0.0)

    ω_el_l = q_nl'(ω_el_n)
    ω_el_b = q_lb'(ω_el_l)
    ω_lb_b = ω_eb_b - ω_el_b

    ω_ie_e = SVector{3}(0, 0, ω_ie)
    ω_ie_b = q_eb'(ω_ie_e)
    ω_ib_b = ω_ie_b + ω_eb_b

    #update ẋ_pos
    ẋ_pos.q_lb .= dt(q_lb, ω_lb_b)
    ẋ_pos.q_el .= dt(q_el, ω_el_l)
    ẋ_pos.Δx = v_eOb_n[1]
    ẋ_pos.Δy = v_eOb_n[2]
    ẋ_pos.h_e = -v_eOb_n[3]

    #build output
    pos = PosData(q_nb, q_eb, REuler(q_nb), q_en, n_e, LatLon(n_e), h_e,
        Altitude{Orthometric}(h_e, n_e), SVector{2}(x.pos.Δx, x.pos.Δy))

    vel = VelData(ω_eb_b, ω_el_n, ω_lb_b, ω_ie_b, ω_ib_b, v_eOb_b, v_eOb_n)

    return KinData(pos, vel)

end

function renormalize!(x_kin::KinLTFX, ε = 1e-10)
    #we need both calls executed, so | must be used here instead of ||
    renormalize_block!(x_kin.pos.q_lb, ε) | renormalize_block!(x_kin.pos.q_el, ε)
end

######################### ECEF Kinematics #########################

struct KinECEF <: AbstractKinematics end

const PosECEFXTemplate = ComponentVector(q_eb = zeros(4), n_e = zeros(3), Δx = 0.0, Δy = 0.0, h_e = 0.0)
const KinECEFXTemplate = ComponentVector(pos = PosECEFXTemplate, vel = VelXTemplate)
const PosECEFX{T, D} = ComponentVector{T, D, typeof(getaxes(PosECEFXTemplate))} where {T, D}
const KinECEFX{T, D} = ComponentVector{T, D, typeof(getaxes(KinECEFXTemplate))} where {T, D}

init_x0(::KinECEF, init::KinInit = KinInit()) = (x=similar(KinECEFXTemplate); init!(x, init); return x)

function init!(x::KinECEFX, init::KinInit)

    @unpack q_nb, Ob, ω_lb_b, v_eOb_b, Δx, Δy = init

    n_e = Ob.l2d
    h_e = Ob.alt
    (R_N, R_E) = radii(Ob)
    v_eOb_n = q_nb * v_eOb_b
    ω_el_n = SVector{3}(
        v_eOb_n[2] / (R_E + Float64(h_e)),
        -v_eOb_n[1] / (R_N + Float64(h_e)),
        0.0)

    ω_el_b = q_nb'(ω_el_n)
    ω_eb_b = ω_el_b + ω_lb_b

    q_en = ltf(n_e)
    q_eb = q_en ∘ q_nb

    x.pos.q_eb .= q_eb[:]
    x.pos.n_e .= n_e[:]
    x.pos.Δx = Δx
    x.pos.Δy = Δy
    x.pos.h_e = h_e
    x.vel.ω_eb_b .= ω_eb_b
    x.vel.v_eOb_b .= v_eOb_b

end

function f_kin!(ẋ_pos::PosECEFX, x::KinECEFX)

    q_eb = RQuat(x.pos.q_eb, normalization = false)
    n_e = NVector(x.pos.n_e, normalization = false)
    ω_eb_b = SVector{3}(x.vel.ω_eb_b)
    v_eOb_b = SVector{3}(x.vel.v_eOb_b)
    h_e = Altitude{Ellipsoidal}(x.pos.h_e[1])

    q_en = ltf(n_e)
    q_nb = q_en' ∘ q_eb

    (R_N, R_E) = radii(n_e)
    v_eOb_n = q_nb(v_eOb_b)
    ω_el_n = SVector{3,Float64}(
        v_eOb_n[2] / (R_E + Float64(h_e)),
        -v_eOb_n[1] / (R_N + Float64(h_e)),
        0.0)

    ω_el_b = q_nb'(ω_el_n)
    ω_lb_b = ω_eb_b - ω_el_b

    ω_ie_e = SVector{3,Float64}(0, 0, ω_ie)
    ω_ie_b = q_eb'(ω_ie_e)
    ω_ib_b = ω_ie_b + ω_eb_b

    #update ẋ_pos
    ẋ_pos.q_eb .= dt(q_eb, ω_eb_b)
    ẋ_pos.n_e .= q_en(ω_el_n × SVector{3,Float64}(0,0,-1))
    ẋ_pos.Δx = v_eOb_n[1]
    ẋ_pos.Δy = v_eOb_n[2]
    ẋ_pos.h_e = -v_eOb_n[3]

    #build output
    pos = PosData(q_nb, q_eb, REuler(q_nb), q_en, n_e, LatLon(n_e), h_e,
        Altitude{Orthometric}(h_e, n_e), SVector{2}(x.pos.Δx, x.pos.Δy))

    vel = VelData(ω_eb_b, ω_el_n, ω_lb_b, ω_ie_b, ω_ib_b, v_eOb_b, v_eOb_n)

    return KinData(pos, vel)

end

function renormalize!(x_kin::KinECEFX, ε = 1e-10)
    #we need both calls executed, so | must be used here instead of ||
    renormalize_block!(x_kin.pos.q_eb, ε) | renormalize_block!(x_kin.pos.n_e, ε)
end

############################ Plotting ################################

function plots(t, data::AbstractVector{<:KinData}; mode, save_path, kwargs...)

    sa = StructArray(data)

    plots(t, sa.pos; mode, save_path, kwargs...)
    plots(t, sa.vel; mode, save_path, kwargs...)
end

function plots(t, data::AbstractVector{<:PosData}; mode, save_path, kwargs...)

    @unpack e_nb, ϕ_λ, h_e, h_o, Δxy = StructArray(data)

    pd = Dict{String, Plots.Plot}()

    splt_h_e = thplot(t, h_e; title = "", kwargs...)

    #remove the title added by the Altitude TH recipe
    splt_h = thplot(t, h_e; title = "", kwargs...)
             thplot!(t, h_o; title = "", kwargs...)

    #remove the title added by the LatLon TH recipe
    splt_latlon = thplot(t, ϕ_λ; title = "", th_split = :v, kwargs...)

    splt_xy = thplot(t, Δxy;
        label = [L"$\int v_{eb}^{x_n} dt$" L"$\int v_{eb}^{y_n} dt$"],
        ylabel = [L"$\Delta x\ (m)$" L"$\Delta y \ (m)$"],
        th_split = :h, link = :none, kwargs...)

    pd["01_e_nb"] = thplot(t, e_nb;
        plot_title = "Attitude (Airframe/NED)",
        rot_ref = "n", rot_target = "b",
        kwargs...)

    pd["02_Ob_geo"] = plot(splt_latlon, splt_h;
        layout = grid(1, 2, widths = [0.67, 0.33]),
        plot_title = "Position (WGS84)",
        kwargs..., plot_titlefontsize = 20) #override titlefontsize after kwargs

    pd["03_Ob_xyh"] = plot(splt_xy, splt_h_e;
        layout = grid(1, 2, widths = [0.67, 0.33]),
        plot_title = "Position (Local Cartesian)",
        kwargs..., plot_titlefontsize = 20) #override titlefontsize after kwargs

    save_plots(pd; save_path)

    #when we assemble a plot from multiple subplots, the plot_titlefontsize
    #attribute no longer works, and it is titlefontisze what determines the font
    #size of the overall figure title (which normally is used for subplots).
    #however, we can still override it specifically for this plot

    #debug mode plots:
    # wander angle

    #maybe add a Trajectory user recipe for Vectors of 3DVector so that i can
    #pass it a Vector series directly.
    #also trplot
    # Ob_xyh_voa = VectorOfArray(Ob_xyh)
    # plt_Ob_xyh_3D = plot(collect(view(Ob_xyh_voa,i,:) for i ∈ 1:3)...;
    #     camera = (30, 45))
    #aspect_ratio attribute does not work for 3d figures
    # savefig(plt_Ob_xyh_3D, joinpath(save_path, "Ob_xyh_3D.png"))
end

function plots(t, data::AbstractVector{<:VelData}; mode, save_path, kwargs...)

    @unpack v_eOb_b, v_eOb_n, ω_lb_b, ω_el_n = StructArray(data)

    pd = Dict{String, Plots.Plot}()

    pd["04_ω_lb_b"] = thplot(t, ω_lb_b;
        plot_title = "Angular Velocity (Airframe/LTF) [Airframe]",
        label = ["Roll Rate" "Pitch Rate" "Yaw Rate"],
        ylabel = [L"$p \ (rad/s)$" L"$q \ (rad/s)$" L"$r \ (rad/s)$"],
        th_split = :h,
        kwargs...)

    pd["05_ω_el_n"] = thplot(t, ω_el_n;
        plot_title = "Local Tangent Frame Transport Rate (LTF/ECEF) [NED]",
        ylabel = L"$\omega_{el}^{l} \ (rad/s)$",
        th_split = :h,
        kwargs...)

    pd["06_v_eOb_n"] = thplot(t, v_eOb_n;
        plot_title = "Velocity (Airframe/ECEF) [NED]",
        label = ["North" "East" "Down"],
        ylabel = [L"$v_{eb}^{N} \ (m/s)$" L"$v_{eb}^{E} \ (m/s)$" L"$v_{eb}^{D} \ (m/s)$"],
        th_split = :h,
        kwargs...)

    pd["07_v_eOb_b"] = thplot(t, v_eOb_b;
        plot_title = "Velocity (Airframe/ECEF) [Airframe]",
        ylabel = [L"$v_{eb}^{x_b} \ (m/s)$" L"$v_{eb}^{y_b} \ (m/s)$" L"$v_{eb}^{z_b} \ (m/s)$"],
        th_split = :h,
        kwargs...)

    save_plots(pd; save_path)

end

end #module


text:
# only use a function and no module to keep it simple
function scrape_heise(maxpages, ignore_year)
    host = "https://www.heise.de"
    hostname = split(host, ".")[2]
    println("host: $(host)")

    # collect
    articles = []
    i = 1
    while i <= maxpages
        println("calling site $(i)")
        req = get("$host/thema/K%C3%BCnstliche-Intelligenz/seite-$i")
        html = parsehtml(convert(String, req.data))
        relevant = Gumbo.HTMLNode
        if i == 1   # skip first list of articles ("Top Artikel")
            relevant = matchall(Selector(".themenseite-thema__article-list"),html.root)[2]
        else
            relevant = matchall(Selector(".themenseite-thema__article-list"),html.root)[1]
        end
        nodes = matchall(Selector(".a-article-teaser__link"),relevant)
        if (length(nodes) > 0)
            append!(articles, nodes)
        else
            break
        end
        i += 1
    end
    println("found $(min(i,maxpages)) sites with $(length(articles)) articles")

    # parse
    println("extracting data from html")
    out = []
    dates = OrderedDict{Int, Int}()
    for article in articles
        date = matchall(Selector("time"), article)[1].attributes["datetime"]
        date = split(date, "T")[1]  # yyyy-dd-mm
        year = parse(Int, split(date, "-")[1])  # yyyy
        if year == ignore_year
            continue
        end

        if haskey(dates, year)
            dates[year] += 1
        else
            dates[year] = 1
        end

        title = nodeText(matchall(Selector("span.a-article-teaser__title-text"), article)[1])
        title = strip(title)    # remove whitespaces
        link = host * matchall(Selector("a"), article)[1].attributes["href"]
        push!(out, [date, title, link])
    end

    # export
    println("writing data to csv")
    writecsv("output/$hostname.csv", out)

    # plot
    return dates
end


text:
#!/usr/bin/env julia
Pkg.test("GLMNet")   # force it to compile the dep
using InvariantCausal
using Base.Test
println("Starting tests")
@time begin
include("test_regression.jl")
include("test_search.jl")
end


text:
prefix(::Type{Float32}) = Symbol()
prefix(::Type{Float64}) = :d
prefix(::Type{Complex{Float32}}) = :c
prefix(::Type{Complex{Float64}}) = :z
prefix(::Any) = Symbol()

const spell_out_greek = Dict{String,String}(
    "α" => "alpha",
    "Α" => "capitol_alpha",
    "β" => "beta",
    "Β" => "capitol_beta",
    "π" => "pi",
    "Π" => "capitol_pi",
    "σ" => "sigma",
    "Σ" => "capitol_sigma",
    "λ" => "lambda",
    "Λ" => "capitol_lambda",
    "γ" => "gamma",
    "Γ" => "capitol_gamma",
    "δ" => "delta",
    "Δ" => "capitol_delta",
    "ϵ" => "epsilon_err",
    "Ε" => "capitol_epsilon_err",
    "ζ" => "zeta",
    "Ζ" => "capitol_zeta",
    "η" => "eta",
    "Η" => "capitol_eta",
    "θ" => "theta",
    "Θ" => "capitol_theta",
    "ι" => "iota",
    "Ι" => "capitol_iota",
    "κ" => "kappa",
    "Κ" => "capitol_kappa",
    "ν" => "nu",
    "ξ" => "xi",
    "Ξ" => "capitol_xi",
    "ρ" => "rho",
    "Ρ" => "capitol_rho",
    "τ" => "tau",
    "Τ" => "capitol_tau",
    "υ" => "upsilon",
    "Υ" => "capitol_upsilon",
    "ϕ" => "phi",
    "Φ" => "capitol_phi",
    "χ" => "chi",
    "Χ" => "capitol_chi",
    "ψ" => "psi",
    "Ψ" => "capitol_psi",
    "ω" => "omega",
    "Ω" => "capitol_omega",
    "ℯ" => "euler",
    "₀" => "0",
    "₁" => "1",
    "₂" => "2",
    "₃" => "3",
    "₄" => "4",
    "₅" => "5",
    "₆" => "6",
    "₇" => "7",
    "₈" => "8",
    "₉" => "9",
    "⁰" => "_zeroth",
    "¹" => "_first",
    "²" => "_second",
    "³" => "_third",
    "⁴" => "_fourth",
    "⁵" => "_fifth",
    "⁶" => "_sixth",
    "⁷" => "_seventh",
    "⁸" => "_eight",
    "⁹" => "_ninth",
    "ₐ" => "a",
    "ₑ" => "e",
    "ᵢ" => "i",
    "ₖ" => "k",
    "ⱼ" => "j",
    "ₗ" => "l",
    "ₘ" => "m",
    "ₙ" => "n",
    "ₒ" => "o",
    "ₚ" => "p",
    "ᵣ" => "r",
    "ₕ" => "h",
    "ₛ" => "s",
    "ₜ" => "t",
    "ᵤ" => "u",
    "ᵥ" => "v",
    "ₓ" => "x"
)
const greek_regex_match = r"α|Α|β|Β|π|Π|σ|Σ|λ|Λ|γ|Γ|δ|Δ|ϵ|Ε|ζ|Ζ|η|Η|θ|Θ|ι|Ι|κ|Κ|ν|ξ|Ξ|ρ|Ρ|τ|Τ|υ|Υ|ϕ|Φ|χ|Χ|ψ|Ψ|ω|Ω|ℯ|₀|₁|₂|₃|₄|₅|₆|₇|₈|₉|⁰|¹|²|³|⁴|⁵|⁶|⁷|⁸|⁹|ₐ|ₑ|ᵢ|ₖ|ⱼ|ₗ|ₘ|ₙ|ₒ|ₚ|ᵣ|ₕ|ₛ|ₜ|ᵤ|ᵥ|ₓ"

function remove_unicode(str)
    replace(str, greek_regex_match => m -> spell_out_greek[m])
end

# const doloop_regex = r" enddo\n\s+end"
for_to_do_loops(str) = replace(replace(str, "for" => "do"), r"    enddo\n\s+end" => "end do")
strip_doublebreaks(str) = replace(str, r"\n\s+\n" => "\n")
clip_ends(str) = str[8:end-4]
clean_fortran_string(str) = str |> strip_doublebreaks |> clip_ends |> for_to_do_loops |> remove_unicode

# function translate(expr::Expr, ::Type{StaticArrays.SArray{Tuple{N},Float64,1,N}})  where N
#     @show expr
#     Expr(:vect, expr.args[2].args...), 0
# end
function translate(func::Symbol, return_type)::Tuple{Symbol,Int}
    get(FortranIntrinsics, (func, return_type)) do
        if func == :+ || func == :- || func == :*
            return func, 0
        elseif func == :/ #What if inputs are integers?
            return func, 0
        elseif func == :^
            return :power, 0 #define a macro
        end
        Symbol(prefix(return_type), func), 0
    end
end

const FortranIntrinsics = Dict{Tuple{Symbol,DataType},Tuple{Symbol,Int}}(
    (:round, Float32) => (:aint, 4),
    (:round, Float64) => (:aint, 8),
    (:imag, Float32) => (:aimag, 0),
    (:besselj0, Float64) => (:bessel_j0, 0),
    (:besselj1, Float64) => (:bessel_j1, 0),
    (:besselj, Float64) => (:bessel_jn, 0),
    (:bessely0, Float64) => (:bessel_y0, 0),
    (:bessely1, Float64) => (:bessel_y1, 0),
    (:bessely, Float64) => (:bessel_yn, 0),
    (:complex, Float32) => (:complex, 0),
    (:complex, Float64) => (:complex, 0),
    (:conj, Float32) => (:conjg, 0),
    (:conj, Float64) => (:dconjg, 0),
    (:cot, Float32) => (:cotan, 0),
    (:cot, Float64) => (:dcotan, 0),
    (:cotd, Float32) => (:cotand, 0),
    (:cotd, Float64) => (:dcotand, 0),
    (:dot, Float32) => (:dot_product, 0),
    (:dot, Float64) => (:dot_product, 0),
    (:dot, Complex{Float32}) => (:dot_product, 0),
    (:dot, Complex{Float64}) => (:dot_product, 0),
    (:dot, Complex{Float64}) => (:dot_product, 0),
    (:eps, Float32) => (:epsilon, 0),
    (:eps, Float64) => (:epsilon, 0),
    (:Int, Int32) => (:int, 4),
    (:Int, Int64) => (:int, 8),
    (:length, Int32) => (:size, 4),
    (:length, Int64) => (:size, 8),
    (:log, Float32) => (:alog, 0),
    (:log10, Float32) => (:alog10, 0),
    (:lgamma, Float32) => (:algama, 0),
    (:lgamma, Float64) => (:dlgama, 0),
    (:max, Float64) => (:max, 0),
    (:maximum, Float64) => (:maxval, 0),
    (:indmax, Float64) => (:maxloc, 0),
    (:min, Float64) => (:min, 0),
    (:minimum, Float64) => (:minval, 0),
    (:indmin, Float64) => (:minloc, 0),
    (:mod, Float32) => (:modulus, 0),
    (:mod, Float64) => (:modulus, 0),
    (:precision, Int32) => (:precision, 0),
    (:sum, Int32) => (:sum, 0),
    (:sum, Int64) => (:sum, 0),
    (:sum, Float64) => (:sum, 0),
    (:sum, Complex{Float32}) => (:sum, 0),
    (:sum, Complex{Float64}) => (:sum, 0),
    (:size, Int32) => (:ubound, 3),#Not even possible?
    (:size, Int64) => (:ubound, 8)
)

const OrderSwap = Dict{Tuple{Symbol,DataType},Tuple{String,String}}(
    (:copysign, Int32) => ("sign(", ")"),
    (:copysign, Int64) => ("sign(", ")"),
    (:copysign, Float32) => ("sign(", ")"),
    (:copysign, Float64) => ("sign(", ")"),
)

const supported_types = Tuple{DataType,String}[(Float32,"real(4)"), (Float64,"real(8)"), (Int32,"integer(4)"), (Int64,"integer(8)"), (Complex{Float32},"complex(4)"), (Complex{Float64},"complex(8)")]

for (t, s) in supported_types
    @eval sized_matrix_type(::Type{<:SVector{N,$t}}) where N = $s * ", dimension($N)"
    @eval sized_matrix_type(::Type{<:SMatrix{M,N,$t}}) where {M,N} = $s * ", dimension($M,$N)"
    @eval sized_matrix_type(::Type{<:MMatrix{M,N,$t}}) where {M,N} = $s * ", dimension($M,$N)"
    @eval sized_matrix_type(::Type{<:StaticArrays.SArray{Tuple{M,N,O},$t}}) where {M,N,O} = $s * ", dimension($M,$N,$O)"
    @eval sized_matrix_type(::Type{<:StaticArrays.SArray{Tuple{M,N,O,P},$t}}) where {M,N,O,P} = $s * ", dimension($M,$N,$O,$P)"
    @eval sized_matrix_type(::Type{<:StaticArrays.SArray{Tuple{M,N,O,P,Q},$t}}) where {M,N,O,P,Q} = $s * ", dimension($M,$N,$O,$P,$Q)"
    @eval sized_matrix_type(::Type{<:StaticArrays.SArray{Tuple{M,N,O,P,Q,R},$t}}) where {M,N,O,P,Q,R} = $s * ", dimension($M,$N,$O,$P,$Q,$R)"
    @eval sized_matrix_type(::Type{<:StaticArrays.SArray{Tuple{M,N,O,P,Q,R,S},$t}}) where {M,N,O,P,Q,R,S} = $s * ", dimension($M,$N,$O,$P,$Q,$R,$S)"

    @eval sized_matrix_type(::Type{<:NTuple{N,$t}}) where N = $s * ", dimension($N)"
    # @eval sized_matrix_type(::Type{<:NTuple{N,$t}}) where N = $s * ", dimension($N)"

    # @eval sized_matrix_type(::Type{<:TriangularMatrices.AbstractTriangularMatrix{$t,N,L}}) where {N,L} = $s * ", dimension($L)"
    # @eval sized_matrix_type(::Type{<:TriangularMatrices.AbstractSymmetricMatrix{$t,N,L}}) where {N,L} = $s * ", dimension($L)"
    # @eval sized_matrix_type(::Type{<:TriangularMatrices.MutableRecursiveMatrix{$t,M,N,L}}) where {M,N,L} = $s * ", dimension($L)"
end

#     @eval sized_matrix_type(::Type{<:MutableRecursiveMatrix{$t,M,N,L}}) where {M,N,L} = $s * ", dimension($L)"
# end

function translate_type(t::DataType)
    local out::String
    if t == Float32
        out = "real(4)"
    elseif t == Float64
        out = "real(8)"
    elseif t == Int32
        out = "integer(4)"
    elseif t == Int64
        out = "integer(8)"
    elseif t <: NTuple
        out = sized_matrix_type(t)
    elseif t <: StaticArrays.SArray
        out = sized_matrix_type(t)
    elseif t <: StaticArrays.MArray
        out = sized_matrix_type(t)
    # elseif t <: TriangularMatrices.TriangularMatrix
    #     out = sized_matrix_type(t)
    # elseif t <: TriangularMatrices.SymmetricMMatrix
    #     out = sized_matrix_type(t)
    # elseif t <: TriangularMatrices.MutableRecursiveMatrix
    #     out = sized_matrix_type(t)
    else
        throw("Could not match type $t.")
    end
    out
end

function symbol_vector_to_string(v::Vector{Symbol})
    if length(v) == 1
        return string(v[1])
    elseif length(v) < 10
        return reduce((x,y) -> string(x)*", "*string(y), v)
    end
    out = string(v[1]) * ", "
    for i ∈ 2:length(v)-1
        out *= string(v[i]) * ", "
        i % 7 == 0 && (out *= "&\n            ")
    end
    out *= string(v[end])
    out
end

function add_fortran_types!(types::Dict{Union{Symbol,Expr},DataType}, arguments, d::Dict{DataType,Vector{Symbol}}, intent = "")
    out = ""
    for k ∈ keys(d)
        syms = d[k]
        append!(arguments, syms)
        out *= "    " * translate_type(k) * intent * " :: " * symbol_vector_to_string(syms) * "\n"
        for sym ∈ syms
            types[sym] = k
        end
    end
    out
end
function add_fortran_types(d::Dict{DataType,Vector{Symbol}}, intent = "")
    out = ""
    for k ∈ keys(d)
        syms = d[k]
        out *= "    " * translate_type(k) * intent * " :: " * symbol_vector_to_string(syms) * "\n"
    end
    out
end

function update_work_vars!(work_vars, types, x)
    if !haskey(work_vars, types[x.args[2]])
        types[x.args[1]] = types[x.args[2]]
        work_vars[types[x.args[2]]] = Symbol[x.args[1]]
    else
        type_vec = work_vars[types[x.args[2]]]
        x.args[1] ∉ type_vec && push!(type_vec, x.args[1])
    end
    nothing
end

function update_work_vars!(work_vars, types, assignee_sym, T_assigned, x)
    if !haskey(work_vars, T_assigned)
        types[assignee_sym] = T_assigned
        work_vars[T_assigned] = Symbol[ assignee_sym ]
    else
        type_vec = work_vars[T_assigned]
        assignee_sym ∉ type_vec && push!(type_vec, assignee_sym)
    end
    nothing
end

const broadcast_ops = (:.+, :.-, :.*, :./)
const translate_broadcast = Dict{Symbol,Symbol}(
    :.+ => :+,
    :.- => :-,
    :.* => :*,
    :./ => :/
)
function remove_broadcast(x)
    if isexpr(x, :call) && x.args[1] ∈ broadcast_ops
        x.args[1] = translate_broadcast[x.args[1]]
    elseif isexpr(x, :.=)
        x.head = :(=)
    end
    x
end

function check_stuff(x, types, work_vars, num_blocks)
    if isa(x, Expr)
        x.head == :line && return Symbol() # :! We want to skip line annotations
        if x.head == :call # if the expression is a call, translate it.
            # @show x
            if x.args[1] == :*

                # @show x

                if types[x.args[2]] <: AbstractMatrix && types[x.args[3]] <: AbstractMatrix
                    x.args[1] = :matmul
                    println("Shouldn't see this atm.")
                    types[x] = Test.Base.return_types(*, ntuple(i -> types[x.args[i+1]], length(x.args)-1))[1]
                elseif types[x.args[2]] <: AbstractMatrix && types[x.args[3]] <: AbstractVector
                    x.args[1] = :matmul
                elseif types[x.args[3]] <: AbstractMatrix && types[x.args[2]] <: AbstractVector#then second arg isn't
                    x.args[1] = :matmul
                elseif types[x.args[3]] <: AbstractVector && types[x.args[2]] <: AbstractVector#then second arg isn't
                    x.args[1] = :dot_product
                    types[x] = promote_type(eltype(types[x.args[2]]),eltype(types[x.args[3]]))
                #Hope this works. If not, gemv, double do loops, or a do loop of dot products?
                else
                    # @show x
                    types[x] = Test.Base.return_types(*, ntuple(i -> types[x.args[i+1]], length(x.args)-1))[1]
                end

                # @show x

                # x.args[1] == :dot_product && (x.args[2] = deepcopy(x.args[2].args[2]) )

                # @show x

                # @show x

            elseif isa(x.args[1], Expr) && x.args[1].args[1] == :SVector
                popfirst!(x.args)
                x.head = :vect
                x.args = deepcopy(x.args[1].args)
                types[x] = SVector{length(x.args),types[x.args[1]]}
            elseif isa(x.args[1], Expr) && x.args[1].args[1] == :NTuple
                popfirst!(x.args)
                x.head = :vect
                x.args = deepcopy(x.args[1].args)
                types[x] = NTuple{length(x.args),types[x.args[1]]}
            else
                local return_type::DataType
                try
                    return_type = get!(() -> Test.Base.return_types(eval(x.args[1]), ntuple(i -> begin
                        xarg = x.args[i+1]
                        (isa(xarg, Symbol) || isa(xarg, Expr)) ? types[xarg] : typeof(xarg)
                    end, length(x.args)-1))[1], types, x)
                catch err
                    @show x
                    @show types
                    @show ntuple(i -> x.args[i+1], length(x.args)-1)
                    rethrow(err)
                end
                # @show x, return_type
                x.args[1], kind = translate(x.args[1], return_type)
                kind != 0 && push!(x.args, kind)
                types[x] = return_type
            end

        elseif x.head == :tuple
            types[x] = NTuple{length(x.args), types[x.args[1]]}
        elseif x.head == :for
            push!(x.args[2].args, :enddo)
        elseif x.head == :ref
            types[x] = eltype(types[x.args[1]])
        elseif x.head == :(:)
            x.head = :call
            pushfirst!(x.args, :iterrange)
            types[x] = Int
        elseif x.head == Symbol("'") ##We're transposing. Lets not allocate a temporary...
            # Currently not supported.
            # x.head = :call
            # pushfirst!(x.args, :transpose)
            # types[x] = types[x.args[2]]
            x.head = deepcopy(x.args[1].head)
            x.args = deepcopy(x.args[1].args)
        elseif x.head == Symbol("=") && isa(x.args[1], Symbol)
            # @show x
            if x.args[2] isa Number
                types[x.args[1]] = typeof(x.args[2])
                update_work_vars!(work_vars, types, x.args[1], typeof(x.args[2]), x)
            else
                types[x.args[1]] = types[x.args[2]]
                update_work_vars!(work_vars, types, x)
            end
        elseif x.head == :+=
            x.head = Symbol("=")
            x.args[2] = :( $(x.args[1]) + ( $(x.args[2]) ) )
            types[x.args[2]] = types[x.args[1]]
            isa(x.args[1], Symbol) && update_work_vars!(work_vars, types, x)
        elseif x.head == :-=
            x.head = Symbol("=")
            x.args[2] = :( $(x.args[1]) - ( $(x.args[2]) ) )
            types[x.args[2]] = types[x.args[1]]
            isa(x.args[1], Symbol) && update_work_vars!(work_vars, types, x)
        elseif x.head == :*=
            x.head = Symbol("=")
            x.args[2] = :( $(x.args[1]) * ( $(x.args[2]) ) )
            types[x.args[2]] = types[x.args[1]]
            isa(x.args[1], Symbol) && update_work_vars!(work_vars, types, x)
        elseif x.head == :/=
            x.head = Symbol("=")
            x.args[2] = :( $(x.args[1]) / ( $(x.args[2]) ) )
            types[x.args[2]] = types[x.args[1]]
            isa(x.args[1], Symbol) && update_work_vars!(work_vars, types, x)
        elseif x == :π
            x = Symbol(Float64(π))
        elseif x == :ℯ
            x = Symbol(Float64(ℯ))
        elseif x.head == :macrocall
            x = x.args[2]
        elseif x.head == :block
            num_blocks[] += 1
        end
    end
    x
end

#define power(X, Y) ((X) ** (Y))

# Test.Base.return_type
function translate_expr(expr::Expr, outputs::Dict{DataType,Vector{Symbol}} = Dict{DataType,Vector{Symbol}}(), inputs::Dict{DataType,Vector{Symbol}} = Dict{DataType,Vector{Symbol}}(), size_inputs = Dict{DataType,Vector{Symbol}}(), inout = Dict{DataType,Vector{Symbol}}())

    arguments = Symbol[]
    types = Dict{Union{Symbol,Expr},DataType}()
    FortranCode = add_fortran_types!(types, arguments, size_inputs, ", intent(in)")
    FortranCode *= add_fortran_types!(types, arguments, inputs, ", intent(in)")
    FortranCode *= add_fortran_types!(types, arguments, inout, ", intent(inout)")
    FortranCode *= add_fortran_types!(types, arguments, outputs, ", intent(out)")

    work_vars = Dict{DataType,Vector{Symbol}}()

    num_blocks = Ref(0)

    q = MacroTools.postwalk(x -> check_stuff(x, types, work_vars, num_blocks), expr)
    q = MacroTools.postwalk(x -> isexpr(x, :ref) ? (x.head = :call; x) : remove_broadcast(x), q)

    if num_blocks[] > 1
        # squashed = MacroTools.postwalk(x -> isexpr(x,:block) && length(x.args) > 1 ? :($[x.args...]) : x, q)
        # squashed = MacroTools.striplines(MacroTools.flatten(q))
        # display(squashed)
        # # @show typeof(squashed)
        # finalq = MacroTools.@q begin end
        # for squashed_block ∈ squashed
        #     isa(squashed_block, AbstractArray) ? append!(finalq.args, squashed_block) : push!(finalq.args, squashed_block)
        # end
        finalq = MacroTools.striplines(MacroTools.flatten(q))
    else
        finalq = MacroTools.striplines(q)
    end
    FortranCode *= add_fortran_types(work_vars)

    FortranCode *= clean_fortran_string(string(finalq))
    FortranCode, arguments
end
function create_subroutine(name::String, expr::Expr; outputs::Dict{DataType,Vector{Symbol}} = Dict{DataType,Vector{Symbol}}(), inputs::Dict{DataType,Vector{Symbol}} = Dict{DataType,Vector{Symbol}}(), size_inputs = Dict{DataType,Vector{Symbol}}(), inout = Dict{DataType,Vector{Symbol}}())
    FortranCode, arguments = translate_expr(expr, outputs, inputs, size_inputs, inout)
    FortranCode = "subroutine $name($(symbol_vector_to_string(arguments)))\n" * FortranCode * "\nend subroutine $name"
    FortranCode
end

function create_module(module_name, body)
    """
    module $module_name

    #define iterrange(X,Y) (X),(Y)
    #define power(X,Y) (X) ** (Y)

    implicit none

    contains

    """ * body * """

    end module $module_name
    """
end


text:
function f()
         for file in readdir()
                         file2 = file * "2"
                         open(file2, "w") do io
                                 for line in eachline(file)
                                         println(io, split(line)[end])
                                 end
                         end
                         cp(file2, file; force=true)
                         rm(file2)
         end
end

cd(f, joinpath(pwd(), "usrdict", "sogou"))


text:
######################
# Instruction on functions for formulas:
# The formula should accept named arguments weight = kg, height = cm, and age = years
# The function should return total daily (24 hour) dosage in milligrams or else as must be specified in the rlabel argument to the app.
# In general, return values in floating point form.
######################

"""
    function amoxicillin_po(; age, weight, height)
Source: https://www.mayoclinic.org/drugs-supplements/penicillin-oral-route-injection-route-intravenous-route-intramuscular-route/proper-use/drg-20062334
For amoxicillin:
For bacterial infections:
For oral dosage forms (capsules, oral suspension, tablets, and chewable tablets):
Adults, teenagers, and children weighing more than 40 kilograms (kg) (88 pounds)
  250 to 500 milligrams (mg) every eight hours or 500 to 875 mg every twelve hours, depending on the type and severity of the infection.
Neonates and infants up to 3 months of age
  The usual dose is 15 mg per kg (6.8 mg per pound) of body weight or less every twelve hours.
Infants 3 months of age and older and children weighing up to 40 kg (88 lbs.)
  The usual dose is 6.7 to 13.3 mg per kg (3 to 6 mg per pound) of body weight every eight hours or
  12.5 to 22.5 mg per kg (5.7 to 10.2 mg per pound) of body weight every twelve hours.
"""
function amoxicillin_po(; age, weight, height)
    if age <= 3/12
        return weight * 30.0
    elseif age <= 16 && weight <= 40
        return weight * 37.0
    else
        return 1500.0
  end

"""
    function metformin(; age, weight, height)
Metformin in age < 7 years or nonobese is not generally used. Otherwise monotherapy dose starts at 1000 mg.
"""
function metformin_mg(; age, weight, height)
    if age < 7 || weight < 20
        return 0.0
    else
        return 1000.0
    end
end

"""
    function metoprolol_extended(; age, weight, height)

"""
function metoprolol_extended(; age, weight, height)
    if age < 18
        return min(weight, 12.5)
    else
        return 25.0
    end
end

"""
    levetiracetam_starting_po(; age, weight, height)

Source:    https://www.rxlist.com/keppra-drug.htm#indications
Initial treatment
Adults 16 Years Of Age And Older
Initiate treatment with a daily dose of 1000 mg/day
Pediatric Patients
1 Month To < 6 Months
Initiate treatment with a daily dose of 14 mg/kg in 2 divided doses
6 Months To < 16 Years
Initiate treatment with a daily dose of 20 mg/kg in 2 divided doses
"""
function levetiracetam_starting_po(; age, weight, height)
    if age <= 0.5
        return weight * 14.0
    elseif age <= 16
        return min(1000.0, weight * 20.0)
    else
        return 1000.0
    end
end

"""
    lamotrigine_starting_po(; age, weight, height)

Refer to tables at https://www.rxlist.com/lamictal-drug.htm#dosage
Note dose may change based on other medications ESPECIALLY VALPROATE
"""
function lamotrigine_starting_po(; age, weight, height)
    if age < 2
        return 0.0
    elseif age <= 12
        weight < 14 && return 1
        weight < 27 && return 2
        weight < 34 && return 4
        weight < 40 && return 5
        return min(25, weight * 0.3)
    else
        return 25
    end
end


text:
function BarycentricCoordinates(simplex1::AbstractArray{Float64, 2},
								simplex2::AbstractArray{Float64, 2},
                                    orientation1::Float64, orientation2::Float64,
                                    tolerance::Float64)
    n = size(simplex1, 1)

    numof1in2 = 0
    numof2in1 = 0

    convexexp1in2 = zeros(Float64, n+1, n+1)
    convexexp2in1 = zeros(Float64, n+1, n+1)

    ordered_vertices1 = collect(1:n+1)
    ordered_vertices2 = collect(1:n+1)

    ########################################
    # Simplex 1 with respect to simplex 2
    ########################################

    # Will contain the indices of the vertices of simplex1 that are contained in simplex2
    contained_vertices = zeros(Int, n + 1)

    for column1 = 1:n+1
        # Initialize convex expansion coefficient
        beta_min = 0.0

        for row2 = 1:n+1
            #tmp = copy(simplex2)
			tmp = similar(simplex2)
			tmp[:, :] = simplex2[:, :]
            tmp[:, row2] = simplex1[:, column1]

            beta = det([ones(1, n + 1); tmp]) / orientation2

            # If beta is close enough to zero, then the coefficient is set to zero
            if abs(beta) <= tolerance
                beta = 0.0
            end

            # If beta is close enough to one, then the coefficient is set to 1
            if abs(beta - 1) <= tolerance
                beta = 1.0
            end

            convexexp1in2[row2, column1] = beta

            if beta < beta_min
                beta_min = beta
            end
        end

        if beta_min >= 0
            numof1in2 += 1
            contained_vertices[numof1in2] = column1
        end
    end

    if numof1in2 > 0
        ordered_vertices1[1:numof1in2] = contained_vertices[1:numof1in2]

        if numof1in2 < n + 1
            ordered_vertices1[numof1in2+1:end] = complementary(ordered_vertices1[1:numof1in2], n + 1)
        end
    end

    ########################################
    # Simplex 2 with respect to simplex 1
    ########################################
    # Will contain the indices of the vertices of simplex1 that are contained in simplex2
    contained_vertices = zeros(Int, 1, n+1)

    for column2 = 1:n+1
        # Initialize convex expansion coefficient
        beta_min = 0

        for row1 = 1:n+1
            #tmp = copy(simplex1)
            tmp = similar(simplex1)
            tmp[:, :] = simplex1[:, :]
            tmp[:, row1] = simplex2[:, column2]

            beta = det([ones(1, n + 1); tmp]) / orientation1

            # If beta is close enough to zero, then the coefficient is set to zero
            if abs(beta) <= tolerance
                beta = 0
            end

            # If beta is close enough to one, then the coefficient is set to 1
            if abs(beta - 1) <= tolerance
                beta = 1
            end

            convexexp2in1[row1, column2] = beta

            if beta < beta_min
                beta_min = beta
            end
        end

        if beta_min >= 0
            numof2in1 = numof2in1 + 1
            contained_vertices[numof2in1] = column2
        end
    end

    if numof2in1 > 0
        ordered_vertices2[1:numof2in1] = contained_vertices[1:numof2in1]

        if numof2in1 < n + 1
            ordered_vertices2[numof2in1+1:end] = complementary(ordered_vertices2[1:numof2in1], n + 1)
        end
    end

    return convexexp1in2, convexexp2in1, ordered_vertices1, ordered_vertices2, numof1in2, numof2in1

end


text:
using JLD2, Base.Test

fn = joinpath(tempdir(),"test.jld")
f = jldopen(fn, "w")
write(f, "x1", 1.0+2.0im)
close(f)

f = jldopen(fn, "a")
write(f, "x2", 3.0+4.0im)
close(f)

f = jldopen(fn, "r")
@test read(f, "x1") === 1.0+2.0im
@test read(f, "x2") === 3.0+4.0im
close(f)
f = open(fn)
data = read(f, UInt8, filesize(fn))
close(f)

f = jldopen(fn, "w")
write(f, "x1", 1.0+2.0im)
write(f, "x2", 3.0+4.0im)
close(f)

f = open(fn)
@test read(f, UInt8, filesize(fn)) == data
close(f)


text:
using Test

@test 1+1 == 2


text:
using Snopt
using DelimitedFiles
using PyPlot
import ForwardDiff

using CSV
using DataFrames
# using YAML

# set up boundary constraint wrapper function
function boundary_wrapper(x, params)
    # include relevant globals
    params.boundary_vertices_a
    params.boundary_normals_a
    params.boundary_vertices_b
    params.boundary_normals_b
    params.boundary_vertices_c
    params.boundary_normals_c
    params.boundary_vertices_d
    params.boundary_normals_d
    params.boundary_vertices_e
    params.boundary_normals_e

    # get number of turbines
    nturbines = Int(length(x)/2)

    # extract x and y locations of turbines from design variables vector
    turbine_x = x[1:nturbines]
    turbine_y = x[nturbines+1:end]

    boundcon_a = ff.ray_trace_boundary(boundary_vertices_a, boundary_normals_a, turbine_x[1:31], turbine_y[1:31])
    boundcon_b = ff.ray_trace_boundary(boundary_vertices_b, boundary_normals_b, turbine_x[32:42], turbine_y[32:42])
    boundcon_c = ff.ray_trace_boundary(boundary_vertices_c, boundary_normals_c, turbine_x[43:58], turbine_y[43:58])
    boundcon_d = ff.ray_trace_boundary(boundary_vertices_d, boundary_normals_d, turbine_x[59:72], turbine_y[59:72])
    boundcon_e = ff.ray_trace_boundary(boundary_vertices_e, boundary_normals_e, turbine_x[73:81], turbine_y[73:81])

    # get and return boundary distances
    return [boundcon_a; boundcon_b; boundcon_c; boundcon_d; boundcon_e]
end

# set up spacing constraint wrapper function
function spacing_wrapper(x, params)
    # include relevant globals
    params.rotor_diameter

    # get number of turbines
    nturbines = Int(length(x)/2)

    # extract x and y locations of turbines from design variables vector
    turbine_x = x[1:nturbines]
    turbine_y = x[nturbines+1:end]

    spacecon_a = 2.0*rotor_diameter[1] .- ff.turbine_spacing(turbine_x[1:31], turbine_y[1:31])
    spacecon_b = 2.0*rotor_diameter[1] .- ff.turbine_spacing(turbine_x[32:42], turbine_y[32:42])
    spacecon_c = 2.0*rotor_diameter[1] .- ff.turbine_spacing(turbine_x[43:58], turbine_y[43:58])
    spacecon_d = 2.0*rotor_diameter[1] .- ff.turbine_spacing(turbine_x[59:72], turbine_y[59:72])
    spacecon_e = 2.0*rotor_diameter[1] .- ff.turbine_spacing(turbine_x[73:81], turbine_y[73:81])

    # get and return spacing distances
    return [spacecon_a; spacecon_b; spacecon_c; spacecon_d; spacecon_e]
end

# set up objective wrapper function
@everywhere function aep_wrapper(x, params)
    # include relevant globals
    turbine_z = params.turbine_z
    rotor_diameter = params.rotor_diameter
    hub_height = params.hub_height
    turbine_yaw = params.turbine_yaw
    ct_models = params.ct_models
    generator_efficiency = params.generator_efficiency
    cut_in_speed = params.cut_in_speed
    cut_out_speed = params.cut_out_speed
    rated_speed = params.rated_speed
    rated_power = params.rated_power
    windresource = params.windresource
    power_models = params.power_models
    model_set = params.model_set
    rotor_points_y = params.rotor_points_y
    rotor_points_z = params.rotor_points_z
    obj_scale = params.obj_scale

    # get number of turbines
    nturbines = Int(length(x)/2)

    # extract x and y locations of turbines from design variables vector
    turbine_x = x[1:nturbines]
    turbine_y = x[nturbines+1:end]

    # calculate AEP
    AEP = obj_scale*ff.calculate_aep(turbine_x, turbine_y, turbine_z, rotor_diameter,
                hub_height, turbine_yaw, ct_models, generator_efficiency, cut_in_speed,
                cut_out_speed, rated_speed, rated_power, windresource, power_models, model_set,
                rotor_sample_points_y=rotor_points_y,rotor_sample_points_z=rotor_points_z, hours_per_year=365.0*24.0)

    # return the objective as an array
    return [AEP]
end

# set up optimization problem wrapper function
function wind_farm_opt(x, params)
    it = params.it

    # calculate spacing constraint value and jacobian
    spacing_con = spacing_wrapper(x)
    ds_dx = ForwardDiff.jacobian(spacing_wrapper, x)

    # calculate boundary constraint and jacobian
    boundary_con = boundary_wrapper(x)
    db_dx = ForwardDiff.jacobian(boundary_wrapper, x)

    # # combine constraint jacobians into an overall array of constraint derivatives
    # for i = 1:prob.m
    #     for j = 1:prob.n
    #         values[(i-1)*prob.n+j] = [ds_dx; db_dx][i, j]
    #     end
    # end

    # combine constaint values and jacobians into overall constaint value and jacobian arrays
    c = [spacing_con; boundary_con]
    dcdx = [ds_dx; db_dx]

    # calculate the objective function and jacobian (negative sign in order to maximize AEP)
    AEP = -aep_wrapper(x)[1]
    dAEP_dx = -ForwardDiff.jacobian(aep_wrapper,x)

    it[1] += 1
    params.funcalls_AEP[it[1]] = it[1]
    params.iter_AEP[it[1]] = -AEP

    # set fail flag to false
    fail = false

    # return objective, constraint, and jacobian values
    return AEP, c, dAEP_dx, dcdx, fail
end

# set globals for iteration history
iter_AEP = zeros(Float64, 10000)
funcalls_AEP = zeros(Float64, 10000)

# import model set with wind farm and related details
include("./model_sets/model_set_7_ieacs4_reduced_wind_rose.jl")

# scale objective to be between 0 and 1
obj_scale = 1E-12

# set wind farm boundary parameters
boundary_vertices_a = [10363.8 6490.3; 9449.7 1602.2; 9387.0 1056.6; 9365.1 625.5; 9360.8 360.2; 9361.5 126.9; 9361.3 137.1; 7997.6 1457.9; 6098.3 3297.5;
    8450.3 6455.3; 8505.4 6422.3; 9133.0 6127.4; 9332.8 6072.6; 9544.2 6087.1; 9739.0 6171.2; 9894.9 6316.9; 10071.8 6552.5; 10106.9 6611.1]
boundary_normals_a = [0.9829601758936983 -0.1838186405319916; 0.9934614633172962 -0.11416795042154541; 0.9987121579438882 -0.050734855622757584;
    0.9998686751666075 -0.01620593781838486; 0.9999954987444023 0.0030004151269687495; -0.9998078216567232 -0.019604074934516894; -0.6957179389375846 -0.718315076718037;
    -0.6957275377423737 -0.7183057797532565; -0.8019887481131871 0.5973391397688945; 0.5138086803485797 0.8579047965820281; 0.4252760929807897 0.905063668886888;
    0.2645057513093967 0.9643841078762402; -0.0684295708121141 0.9976559496331737; -0.39636379138742883 0.9180935381958544; -0.6828023205475376 0.7306031693435896;
    -0.7996740386176392 0.6004343694034798; -0.8578802011411015 0.5138497450520954; 0.42552559023380465 0.9049463918134445]
boundary_vertices_b = [5588.4 3791.3; 4670.7 4680.2; 7274.9 7940.8; 7369.9 7896.2; 7455.1 7784.3; 7606.5 7713.0; 7638.9 7708.4; 8297.1 7398.9]
boundary_normals_b = [-0.6957460043611584 -0.7182878931288504; -0.7813688797257963 0.6240694462926818; 0.4249708760634733 0.9052070230051488; 0.7956275395848184 0.6057861159305391;
    0.4260560153872896 0.9046967844268629; 0.14056568619461773 0.9900713549359138; 0.4255255464063141 0.9049464124220882; 0.7996806883794807 -0.6004255129763556]
boundary_vertices_c = [3267.1 10100.6; 4164.1 9586.6; 5749.8 9068.6; 6054.7 8925.3; 1468.5 7781.7; 107.4 9100.0]
boundary_normals_c = [0.49718026396417986 0.8676472699919642; 0.31052117525343714 0.9505664625470563; 0.42535384615162936 0.9050271297392228; 0.24194817066179167 -0.9702891747893577;
    -0.6957228969594285 -0.7183102746351193; -0.30189947425802094 0.9533397649540959]
boundary_vertices_d = [6764.9 8399.7; 4176.8 5158.6; 2047.8 7220.7]
boundary_normals_d = [0.7814306689309158 -0.6239920749930895; -0.6957310325444781 -0.7183023947855072; -0.24248239299288069 0.9701558066045093]
boundary_vertices_e = [8953.7 11901.5; 7048.3 9531.5; 6127.7 9962.7; 4578.1 10464.9; 4524.1 10498.7]
boundary_normals_e = [0.7793586677376737 -0.6265780613955122; -0.4241667101838764 -0.9055841219742026; -0.30829751674447764 -0.9512899879475178; -0.5305632140423848 -0.847645371546978; -0.3019099610801309 0.9533364439695956]

# set globals for use in wrapper functions
struct params_struct{MS, AF, F, ACTM, WR, APM, AF2, AI}
    model_set::MS
    rotor_points_y::AF
    rotor_points_z::AF
    turbine_z::AF
    rotor_diameter::AF
    obj_scale::F
    hub_height::AF
    turbine_yaw::AF
    ct_models::ACTM
    generator_efficiency::AF
    cut_in_speed::AF
    cut_out_speed::AF
    rated_speed::AF
    rated_power::AF
    windresource::WR
    power_models::APM

    boundary_vertices_a::AF2
    boundary_normals_a::AF2
    boundary_vertices_b::AF2
    boundary_normals_b::AF2
    boundary_vertices_c::AF2
    boundary_normals_c::AF2
    boundary_vertices_d::AF2
    boundary_normals_d::AF2
    boundary_vertices_e::AF2
    boundary_normals_e::AF2
    iter_AEP::AF
    funcalls_AEP::AF
    it::AI
end

params = params_struct(model_set, rotor_points_y, rotor_points_z, turbine_z,
    rotor_diameter, obj_scale, hub_height, turbine_yaw, ct_models, generator_efficiency, cut_in_speed,
    cut_out_speed, rated_speed, rated_power, windresource, power_models, boundary_vertices_a, boundary_normals_a, boundary_vertices_b, boundary_normals_b,
    boundary_vertices_c, boundary_normals_c, boundary_vertices_d, boundary_normals_d, boundary_vertices_e,
    boundary_normals_e, iter_AEP, funcalls_AEP, [0])

# initialize design variable array
x_initial = [copy(turbine_x);copy(turbine_y)]
turbine_x_init = deepcopy(turbine_x)
turbine_y_init = deepcopy(turbine_y)
x = [copy(turbine_x);copy(turbine_y)]

# add boundary to plot
plt.clf()
plot(boundary_vertices_a[:,1],boundary_vertices_a[:,2], "b")
plot([boundary_vertices_a[1,1],boundary_vertices_a[end,1]],[boundary_vertices_a[1,2],boundary_vertices_a[end,2]], "b")
plot(boundary_vertices_b[:,1],boundary_vertices_b[:,2], "b")
plot([boundary_vertices_b[1,1],boundary_vertices_b[end,1]],[boundary_vertices_b[1,2],boundary_vertices_b[end,2]], "b")
plot(boundary_vertices_c[:,1],boundary_vertices_c[:,2], "b")
plot([boundary_vertices_c[1,1],boundary_vertices_c[end,1]],[boundary_vertices_c[1,2],boundary_vertices_c[end,2]], "b")
plot(boundary_vertices_d[:,1],boundary_vertices_d[:,2], "b")
plot([boundary_vertices_d[1,1],boundary_vertices_d[end,1]],[boundary_vertices_d[1,2],boundary_vertices_d[end,2]], "b")
plot(boundary_vertices_e[:,1],boundary_vertices_e[:,2], "b")
plot([boundary_vertices_e[1,1],boundary_vertices_e[end,1]],[boundary_vertices_e[1,2],boundary_vertices_e[end,2]], "b")

# for i = 1:length(boundary_vertices_a)
#     vx1 = boundary_vertices_a[i][1]
#     vy1 = boundary_vertices_a[i][2]
#     vx2 = boundary_vertices_a[i+1][1]
#     vy2 = boundary_vertices_a[i+1][2]
#     plot([vx1, ])
# end

# add initial turbine location to plot
for i = 1:length(turbine_x)
    plt.gcf().gca().add_artist(plt.Circle((turbine_x[i],turbine_y[i]), rotor_diameter[1]/2.0, fill=false,color="C0"))
end
plt.show()
# get number of design variables
n_designvariables = length(x_initial)

# get number of constraints
n_spacingconstraints = length(spacing_wrapper(x_initial, params))
n_boundaryconstraints = length(boundary_wrapper(x_initial, params))
n_constraints = n_spacingconstraints + n_boundaryconstraints

# set general lower and upper bounds for design variables
lb = ones(n_designvariables) * -Inf
ub = ones(n_designvariables) * Inf

# set lower and upper bounds for constraints
lb_g = ones(n_constraints) * -Inf
ub_g = zeros(n_constraints)

# generate wrapper function surrogates
spacing_wrapper(x) = spacing_wrapper(x, params)
aep_wrapper(x) = aep_wrapper(x, params)
boundary_wrapper(x) = boundary_wrapper(x, params)
wind_farm_opt(x) = wind_farm_opt(x, params)

# set up options for SNOPT
options = Dict{String, Any}()
options["Derivative option"] = 1
options["Verify level"] = 0
options["Major optimality tolerance"] = 1e-6
options["Major iteration limit"] = 1e6
options["Summary file"] = "./snopt_summarycs4.out"
options["Print file"] = "./snopt_printcs4.out"

# print initial objective value
println("nturbines: ", nturbines)
println("rotor diameter: ", rotor_diameter[1])
println("starting AEP value (MWh): ", aep_wrapper(x_initial)[1]*1E5)

# run and time optimization
t1 = time()
xopt, fopt, info = snopt(wind_farm_opt, x, lb, ub, options)
t2 = time()
clkt = t2-t2

# print optimization results
println("Finished in : ", clkt, " (s)")
println("info: ", info)
println("end AEP value (MWh): ", aep_wrapper(xopt)[1]*1E5)

# extract final turbine locations
turbine_x = copy(xopt[1:nturbines])
turbine_y = copy(xopt[nturbines+1:end])

# add final turbine locations to plot
for i = 1:length(turbine_x)
    plt.gcf().gca().add_artist(plt.Circle((turbine_x_init[i],turbine_y_init[i]), rotor_diameter[1]/2.0, fill=false,color="C2", linestyle="-"))
    plt.gcf().gca().add_artist(plt.Circle((turbine_x[i],turbine_y[i]), rotor_diameter[1]/2.0, fill=false,color="C1", linestyle="-"))
end

# add wind farm boundary to plot
plt.gcf().gca().plot([boundary_vertices_a[:,1];boundary_vertices_a[1,1]],[boundary_vertices_a[:,2];boundary_vertices_a[1,2]], color="C2")
plt.gcf().gca().plot([boundary_vertices_b[:,1];boundary_vertices_b[1,1]],[boundary_vertices_b[:,2];boundary_vertices_b[1,2]], color="C2")
plt.gcf().gca().plot([boundary_vertices_c[:,1];boundary_vertices_c[1,1]],[boundary_vertices_c[:,2];boundary_vertices_c[1,2]], color="C2")
plt.gcf().gca().plot([boundary_vertices_d[:,1];boundary_vertices_d[1,1]],[boundary_vertices_d[:,2];boundary_vertices_d[1,2]], color="C2")
plt.gcf().gca().plot([boundary_vertices_e[:,1];boundary_vertices_e[1,1]],[boundary_vertices_e[:,2];boundary_vertices_e[1,2]], color="C2")

# set up and show turbine layout plot
axis("square")
xlim(0, 11000)
ylim(-500, 13000)
plt.show()
savefig("snoptresults_turbinelayout.png")

# create obj value history plots
plt.clf()
plt.gcf().gca().plot(1:length(funcalls_AEP), funcalls_AEP)
title("Objective Value History")
xlabel("Function Call")
ylabel("AEP (TWh)")
savefig("snoptresults_plot_AEP_history_by_funcall.png")

plt.clf()
plt.gcf().gca().plot(1:length(iter_AEP), iter_AEP)
title("Objective Value at Each Iteration")
xlabel("Iteration")
ylabel("AEP (TWh)")
savefig("snoptresults_plot_AEP_history_by_iter.png")

# write text file with optimization summary info
io = open("snoptresults_summary.txt", "a")
write(io, "number of turbines: ")
writedlm(io, nturbines)
write(io, "rotor diameter: ")
writedlm(io, rotor_diameter[1])
write(io, "starting AEP value (MWh): ")
writedlm(io, aep_wrapper(x_initial)[1]*1E5)
write(io, "Finished in : ")
writedlm(io, clkt)
write(io, "end objective value: ")
writedlm(io, aep_wrapper(xopt)[1]*1E5)
write(io, "\noptimal turbine coordinates:\n")
writedlm(io, [turbine_x turbine_y])
write(io, "\neach function call value\n")
writedlm(io, funcalls_AEP)
close(io)

# write results to csv files
dataforcsv_xopt = DataFrame(xopt = xopt)
dataforcsv_funceval = DataFrame(function_value = funcalls_AEP)
dataforcsv_funceval_at_iteration = DataFrame(function_value = iter_AEP)
CSV.write("cs4_xopt.csv", dataforcsv_xopt)
CSV.write("cs4_functionvalue_log.csv", dataforcsv_funceval)
CSV.write("cs4_functionvalue_at_iteration_log.csv", dataforcsv_funceval_at_iteration)

# # write results to a YAML file (not complete)
# dataforyaml = [
#     "optimal_turbine_coordinates" => [turbine_x turbine_y]
#     "function_evaluation_log" => funcalls_AEP
#     "function_evaluation_at_iterations_log" => iter_AEP
#     ]
# YAML.write_file("cs4_results.yml", dataforyaml)


text: abstract type ProjectiveEntity end


text:
import Blosc

abstract type Compressor end
getCompressor(compdict::Dict) = getCompressor(compressortypes[compdict["id"]],compdict)
getCompressor(::Nothing) = NoCompressor()

struct BloscCompressor <: Compressor
    blocksize::Int
    clevel::Int
    cname::String
    shuffle::Bool
end

"""
    BloscCompressor(;blocksize=0, clevel=5, cname="lz4", shuffle=true)

Returns a `BloscCompressor` struct that can serve as a Zarr array compressor. Keyword arguments are:

* `clevel=5` the compression level, number between 0 (no compression) and 9 (max compression)
* `cname="lz4"` compressor name, can be one of `"blosclz"`, `"lz4"`, and `"lz4hc"`
* `shuffle=true` enables/disables bit-shuffling
"""
BloscCompressor(;blocksize=0, clevel=5, cname="lz4", shuffle=true) =
    BloscCompressor(blocksize, clevel, cname, shuffle)

function getCompressor(::Type{BloscCompressor}, d::Dict)
    BloscCompressor(d["blocksize"], d["clevel"], d["cname"], d["shuffle"] > 0)
end

zuncompress(a, r::AbstractArray, ::BloscCompressor) = copyto!(a, Blosc.decompress(Base.nonmissingtype(eltype(a)), r))

function zcompress(a, f::AbstractArray, c::BloscCompressor)
    Blosc.set_compressor(c.cname)
    r = Blosc.compress(a, level=c.clevel, shuffle=c.shuffle)
    empty!(f)
    append!(f, r)
end

JSON.lower(c::BloscCompressor) = Dict("id"=>"blosc", "cname"=>c.cname,
    "clevel"=>c.clevel, "shuffle"=>c.shuffle ? 1 : 0, "blocksize"=>c.blocksize)

"""
    NoCompressor()

Creates an object that can be passed to ZArray constructors without compression.
"""
struct NoCompressor <: Compressor end

function zuncompress(a, r::AbstractArray, ::NoCompressor)
  copyto!(a, reinterpret(eltype(a),r))
end

function zcompress(a, f::AbstractArray, ::NoCompressor)
  a2 = reinterpret(UInt8,a)
  empty!(f)
  append!(f, a2)
end

JSON.lower(::NoCompressor) = nothing

compressortypes = Dict("blosc"=>BloscCompressor, nothing=>NoCompressor)


text:
using NonUniformRandomVariateGeneration
using Distributions
using BenchmarkTools

function benchNative(α::Float64, β::Float64, N::Int64)
  m::Float64 = 0.0
  for i in 1:N
    m += sampleBeta(α, β)
  end
  return m / N
end

function benchDistributions(α::Float64, β::Float64, N::Int64)
  beta = Distributions.Beta(α, β)
  m::Float64 = 0.0
  for i in 1:N
    m += rand(beta)
  end
  return m / N
end

@btime benchNative(1.5, 2.5, 1024*1024*8)
@btime benchDistributions(1.5, 2.5, 1024*1024*8)


text:
function default_methods!(mt)
  mt["builtins.set"] = Set
  mt["builtins.bytes"] = () -> b""
  mt["builtins.list"] = Base.vect
  mt["collections.defaultdict"] = DataStructures.DefaultDict
  mt["codecs.encode"] = (s, c) -> codeunits(s)
  _setentry!(mt.head, mt["builtins"], "__builtin__")
  _setentry!(mt.head, mt["codecs"], "_codecs")
  mt["__julia__.Set"] = "builtins.set"
  mt["__julia__.Base.CodeUnits"] = "codecs.encode"
  mt["__julia__.__py__.bytes"] = "builtins.bytes"
  mt
end


text:
function part1(input)
    return sum(eval_line1.(readlines(input)))
end

function part2(input)
    return sum(eval_line2.(readlines(input)))
end

function eval_line1(line)
    expr = Meta.parse(replace(replace(line, "+" => "⊙"), "*" => "⊗"))
    return eval(expr)
end

⊙(x, y) = x + y
⊗(x, y) = x * y

function eval_line2(line)
    expr = Meta.parse(replace(replace(line, "+" => "⊠"), "*" => "⊞"))
    return eval(expr)
end

⊠(x, y) = x + y
⊞(x, y) = x * y


text:
function abc()
 println("in abc, evalon=$evalon")
end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
"""
Shapey value using Harsanyi dividends
"""
function fₛⁱ(G::Game, i::Int)
    ∑(Δₕ(G, T) / length(T) for T in subsets(G.N) if i ∈ T)
end

"""
Shapey value using permutations
"""
function fₛ(G::SimpleGame)
    N = length(G.N)
    Nσ = factorial(N)

    shapley = zeros(N)

    for n in 1:Nσ
        σ = nthperm(G.N, n)

        for (i, s) in enumerate(σ)
            shapley[s] += G.v(σ[1:i]) - G.v(σ[1:i - 1])
        end

    end

    return shapley ./ Nσ
end
fₛ(G::GraphGame) = fₛ(graphtoMyerson(G))

"""
Myerson solution, it is equivalent to the Shapley value of the Myerson game.
"""
μ(G::SimpleGame, L::Array{NTuple{2,Int}}) = fₛ(GraphGame(G.N, G.v, L))


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2020
#
# This is a Julia package for discrete event simulation
#

const _handle_exceptions = [true]

# ---------------------------------------------------------
# methods for active clocks
# ---------------------------------------------------------
tau(ac::ActiveClock) = tau(ac.clock)
_busy(ac::ActiveClock) = ac.clock.state == Busy()
sync!(ac::ActiveClock, clk::Clock) = sync!(ac.clock, clk)

delay!(ac::ActiveClock, args...) = delay!(ac.clock, args...)
wait!(ac::ActiveClock, args...; kwargs...) = wait!(ac.clock, args...; kwargs...)
now!(ac::ActiveClock, ex::A) where {A<:Action} = now!(ac.clock, ex)

# ---------------------------------------------------------
# state machine operations of active clocks
# ---------------------------------------------------------
step!(A::ActiveClock, ::Idle, ::Start) = ( A.clock.evcount = 0; A.clock.scount = 0)

function step!(A::ActiveClock, ::Idle, σ::Run)
    t0 = time_ns()
    _cycle!(A.clock, σ.duration, σ.sync)
    put!(A.back, Done(time_ns()-t0))
end

function step!(A::ActiveClock, ::Idle, σ::Finish)
    _finish!(A.clock, σ.tend)
    put!(A.back, Response((A.clock.evcount, A.clock.scount)))
end

step!(A::ActiveClock, ::Union{Idle, Busy}, ::Sync) = nothing

step!(A::ActiveClock, ::ClockState, ::Query) = put!(A.back, Response(A))

step!(A::ActiveClock, ::Union{Idle, Busy}, σ::Register) = _assign(A, σ.x, A.id)

function step!(A::ActiveClock, ::Union{Idle, Busy}, σ::Reset)
    m = A.master[]
    resetClock!(A.clock, m.Δt, t0=m.time, hard=σ.type, unit=m.unit)
    put!(A.back, Response(1))
end

step!(A::ActiveClock, q::ClockState, σ::ClockEvent) = error("transition q=$q, σ=$σ not implemented")

# -------------------------------------
# event loop for an active clock.
# -------------------------------------
function _activeClock(cmd::Channel, resp::Channel)
    info = take!(cmd) # get a pointer to the master clock and id
    ac = ActiveClock(Clock(), info.m, cmd, resp, info.id, threadid())
    ac.clock.id = info.id
    sf = Array{Base.StackTraces.StackFrame,1}[]
    exc = nothing
    ac.clock.state = Idle()
    sync!(ac.clock, ac.master[])

    while true
        σ = take!(cmd)
        if σ isa Stop
            break
        elseif σ isa Diag
            put!(resp, Response((exc, sf)))
        elseif _handle_exceptions[end]
            try
                step!(ac, ac.clock.state, σ)
            catch exc
                sf = stacktrace(catch_backtrace())
                @warn "clock $(ac.id), thread $(ac.thread) exception: $exc"
                put!(resp, Error(exc))  # send error to avoid master hangs
            end
        else
            step!(ac, ac.clock.state, σ)
        end
    end
    # stop task
    # close channel
end

# ---------------------------------------------------------
# starting and destroying active clocks
# ---------------------------------------------------------

# Start a task on each available thread (other than 1).
# - `f::Function`: function to start, has to take two channels as arguments,
# - `ch_size=256`: channel capacity for event transfer between clocks during
#                  each time step.
function _start_threads(f::F, ch_size=256)::Vector{ClockChannel} where {F<:Function}
    ac = ClockChannel[]
    for i in 1:nthreads()
        if i > 1
            ai = ClockChannel(Ref{Task}(),
                    Channel{ClockEvent}(ch_size),
                    Channel{ClockEvent}(ch_size),
                    i, false, 0)
            push!(ac, ai)
        end
    end
    @threads for i in 1:nthreads()
        i > 1 && (ai = ac[i-1]; ai.ref = Ref(@async f(ai.forth, ai.back)))
    end
    return ac
end

"""
    fork!(master::Clock)

Establish copies of a master clock (thread 1) on all parallel threads.
"""
function fork!(master::Clock)
    if (master.id == 0) && (threadid() == 1)
        if VERSION >= v"1.3"
            if nthreads() > 1
                if isempty(master.ac)
                    master.ac = _start_threads(_activeClock)           # startup steps
                    for i in eachindex(master.ac)                    # to all active clocks …
                        put!(master.ac[i].forth, Startup(Ref(master), i)) # send pointer and id
                    end
                else
                    println(stderr, "clock already has $(length(clk.ac)) active clocks!")
                end
            else
                println(stderr, "no parallel threads available!")
            end
        else
            println(stderr, "you have Julia $VERSION, threading is available for ≥ 1.3!")
        end
    else
        println(stderr, "only the master clock on thread 1 can be multiplied!")
    end
end

"""
    collapse!(master::Clock)

Transfer the schedules of the parallel clocks to master and them stop them.

!!! note
    If there are processes on other threads registered to parallel clocks,
    make sure that they aren't needed anymore before calling `collapse`. They
    are not transferred to and cannot be controlled by master.
"""
function collapse!(master::Clock)
    if (master.id == 0) && (threadid() == 1)
        for ac in master.ac
            put!(ac.forth, Query())
            c = take!(ac.back).x.clock
            for (ev, t) in pairs(c.sc.events)    # transfer events to master
                while any(i->i==t, values(master.sc.events)) # in case an event at that time exists
                    t = nextfloat(float(t))                  # increment scheduled time
                end
                master.sc.events[ev] = t
            end
            append!(master.sc.cevents, c.sc.cevents)
            append!(master.sc.samples, c.sc.samples)
            put!(ac.forth, Stop())
        end
        empty!(master.ac)
    else
        println(stderr, "only the master clock on thread 1 can be collapsed!")
    end
end

"""
    PClock(Δt::T=0.01; t0::U=0.0, unit::FreeUnits=NoUnits) where {T<:Number,U<:Number}

Setup a clock with parallel clocks on all available threads.

# Arguments

- `Δt::T=0.01`: time increment > 0. If given Δt ≤ 0, it gets set to 0.01.
- `t0::U=0.0`: start time for simulation,
- `unit::FreeUnits=NoUnits`: clock time unit for explicitl unit setting.

!!! note
    Processes on multiple threads are possible in Julia ≥ 1.3 and with
    [`JULIA_NUM_THREADS > 1`](https://docs.julialang.org/en/v1/manual/environment-variables/#JULIA_NUM_THREADS-1).
"""
function PClock(Δt::T=0.01; t0::U=0.0, unit::FreeUnits=NoUnits) where {T<:Number,U<:Number}
    Δt = Δt > 0 ? Δt : 0.01
    clk = Clock(Δt, t0=t0, unit=unit)
    fork!(clk)
    return clk
end

# ---------------------------------------------------------
# control of active clocks
# ---------------------------------------------------------
"""
```
pclock(clk::Clock, id::Int ) :: C where {C<:AbstractClock}
pclock(ac::ActiveClock, id::Int ) :: C where {C<:AbstractClock}
```
Get a parallel clock to a given clock.

# Arguments
- `master::Clock`: a master clock or
- `ac::ActiveClock`: an active clock,
- `id::Int=threadid()`: thread id, defaults to the caller's current thread.

# Returns
- the master `Clock` if id==0,
- a parallel `ActiveClock` else
"""
function pclock(clk::Clock, id::Int) ::AbstractClock
    if id == 0
        @assert clk.id == 0 "you cannot get master from a local clock!"
        return clk
    elseif id in eachindex(clk.ac)
        put!(clk.ac[id].forth, Query())
        return take!(clk.ac[id].back).x
    else
        println(stderr, "parallel clock $id not available!")
    end
end
function pclock(ac::ActiveClock, id::Int) ::AbstractClock
    if id == ac.clock.id
        return ac
    else
        return pclock(ac.master[], id)
    end
end

"""
    diagnose(master::Clock, id::Int)

Return the stacktrace from parallel clock id.
"""
function diagnose(master::Clock, id::Int)
    if id in eachindex(master.ac)
        if istaskfailed(master.ac[id].ref[])
            return master.ac[id].ref[]
        else
            while isready(master.ac[id].back)
                msg = take!(master.ac[id].back)
                println(msg)
            end
            put!(master.ac[id].forth, Diag())
            return take!(master.ac[id].back).x
        end
    else
        println(stderr, "parallel clock $id not available!")
    end
end

"""
    onthread(f::F, tid::Int) where {F<:Function}

Execute a function f on thread tid (to speed it up).

# Example
```jldoctest
julia> using DiscreteEvents, .Threads

julia> onthread(2) do; threadid(); end
2
```
"""
function onthread(f::F, tid::Int) where {F<:Function}
    t = Task(nothing)
    @assert tid in 1:nthreads() "thread $tid not available!"
    @threads for i in 1:nthreads()
        if i == tid
            t = @async f()
        end
    end
    fetch(t)
end


text:
####################################################################################################
# imgcodecs. Image file reading and writing
####################################################################################################

# 1. Reading and Writing Images

# imdecode: Reads an image from a buffer in memory
imdecode(buf, flags) = @cxx cv::imdecode(buf, flags)  # returs Mat
# buf – Input array or vector of bytes
# flags – The same flags as in imread()

# imencode: Encodes an image into a memory buffer
imencode(ext, img, buf) = @cxx cv::imencode(ext, img, buf)
# ext    – File extension that defines the output format
# img    – Image to be written
# buf    – Output buffer resized to fit the compressed image  => vector<uchar>&
# params – Format-specific parameters. See imwrite()   => const vector<int>&, default = vector<int>()

# imread: Read images from file
imread(filename::String, flags=IMREAD_UNCHANGED) = @cxx cv::imread(pointer(filename), flags)
# IMREAD_UNCHANGED  # 8bit, color or not
# IMREAD_GRAYSCALE  # 8bit, gray
# IMREAD_COLOR      # ?, color
# IMREAD_ANYDEPTH   # any depth, ?
# IMREAD_ANYCOLOR   # ?, any color
# IMREAD_LOAD_GDAL

# customized imread with QFileDialog::getOpenFileName
imread(filename = getOpenFileName(), flags=IMREAD_UNCHANGED) = @cxx cv::imread(pointer(filename), flags)
# invoke with imread()

# imwrite: Saves an image to a specified file
imwrite(filename::String, img) = @cxx cv::imwrite(pointer(filename), img)
imwrite(filename::String, img, params) = @cxx cv::imwrite(pointer(filename), img, params)
# optional: const vector<int>& params=vector<int>()
# filename – Name of the file
# image    – Image to be saved
# params   –
# IMWRITE_JPEG_QUALITY      Default value is 95 {0,100}
# IMWRITE_WEBP_QUALITY      Default value is 100 {1,100}
# IMWRITE_PNG_COMPRESSION   Default value is 3 {0,9}
# IMWRITE_PXM_BINARY        Default value is 1 {0,1}

# customized imread with QFileDialog::getOpenFileName
imwrite(img, filename = getSaveFileName()) = @cxx cv::imwrite(pointer(filename), img)


text:
# License for this file: MIT (expat)
# Copyright 2019, DLR Institute of System Dynamics and Control

"""
    module Simulate_PendulumDAE

Demonstrate how a pendulum described as index 3 DAE can be defined, simulated
and plotted with ModiaMath.

The model is a mass point attached via a rod to a revolute joint where "x" and "y"
coordinates are used as describing variables. The index 3 model is transformed to
an index 1 model using the Gear-Gupta-Leimkuhler transformation:

Starting equations (der(v) = dv/dt):

     der(x) = vx
     der(y) = vy
  m*der(vx) = -x*lambda
  m*der(vy) = -y*lambda-m*g
  x*x + y*y = L*L

Gear/Gupta/Leimkuhler transformation:
   # lambda = der(lambda_int)
   # mue    = der(mue_int)
   # g      = x*x + y*y - L*L
   # G      = dg/d[x;y] = [2*x 2*y]

   der([x,y]) = [vx,vy] - G'*mue
   0 = m*der(vx) + lambda*x
   0 = m*der(vy) + lambda*y + m*g
   0 = x*x + y*y - L*L
   0 = 2*x*vx + 2*y*vy

Modelica model: ModelicaReferenceModels.ODAEs.Pendulum
"""
module Simulate_PendulumDAE

using ModiaMath

@component PendulumDAE(;L=1.0, m=1.0, g=9.81, x_start=L/2.0, y_start=0.5, x_fixed=false) begin    # y_start = -sqrt(L*L - x_start*x_start)
    @assert(L > 0.0)
    @assert(m > 0.0)

    x      = RealScalar(           numericType=ModiaMath.XD_EXP,                  start=x_start, unit="m"    , fixed=x_fixed, nominal=1.0, info="Absolute x-position of mass")
    y      = RealScalar(           numericType=ModiaMath.XD_EXP,                  start=y_start, unit="m"    , fixed=false  , nominal=1.0, info="Absolute y-position of mass")
    der_x  = RealScalar("der(x)",  numericType=ModiaMath.DER_XD_EXP, integral=x,                 unit="m/s"  ,                             info="= der(x)")
    der_y  = RealScalar("der(y)",  numericType=ModiaMath.DER_XD_EXP, integral=y,                 unit="m/s"  ,                             info="= der(y)")

    vx     = RealScalar(           numericType=ModiaMath.XD_EXP,                  start=0.0,     unit="m/s"  , fixed=x_fixed, nominal=1.0, info="Absolute x-velocity of mass")
    vy     = RealScalar(           numericType=ModiaMath.XD_EXP,                  start=0.0,     unit="m/s"  , fixed=false  , nominal=1.0, info="Absolute y-velocity of mass")
    ax     = RealScalar("der(vx)", numericType=ModiaMath.DER_XD_EXP, integral=vx,                unit="m/s^2",                             info="= der(vx)")
    ay     = RealScalar("der(vy)", numericType=ModiaMath.DER_XD_EXP, integral=vy,                unit="m/s^2",                             info="= der(vy)")

    lambda        = RealScalar(numericType=ModiaMath.LAMBDA, unit="N"    , info="Cut-force in the rod")
    mue           = RealScalar(numericType=ModiaMath.MUE,    unit="1/s"  , info="Mue stabilizer (should be zero)")
    residue_xy    = RealScalar(numericType=ModiaMath.FC,     unit="m^2"  , info="Residue of position constraint")
    residue_vx_vy = RealScalar(numericType=ModiaMath.FC,     unit="m^2/s", info="Residue of velocity constraint")
end

function ModiaMath.computeVariables!(p::PendulumDAE, sim::ModiaMath.SimulationState)
    x      = p.x.value
    y      = p.y.value
    vx     = p.vx.value
    vy     = p.vy.value
    lambda = p.lambda.value
    mue    = p.mue.value

    p.der_x.value = vx + x * mue
	p.der_y.value = vy + y * mue
    p.ax.value    = lambda * x/p.m
    p.ay.value    = lambda * y/p.m + p.g

#	if ModiaMath.compute_der_fc(m)
#		_r[5] = x * _derx[1] + y * _derx[2]
#		_r[6] = x * dervx + y * dervy + vx*vx + vy*vy
#	else
	p.residue_xy.value    = (x * x + y * y - p.L * p.L) / 2.0
	p.residue_vx_vy.value = x * vx + y * vy
#    end
    return nothing
end

simulationModel = ModiaMath.SimulationModel(PendulumDAE(L=0.8, m=0.5); stopTime=5.0, tolerance=1e-6)
ModiaMath.print_ModelVariables(simulationModel)
result          = ModiaMath.simulate!(simulationModel, log=true)
ModiaMath.plot(result, [("x", "y"), ("vx", "vy"), "lambda", "mue"])

end


text:
using PATH
using SparseArrays
using Test

@testset "CheckLicense" begin
    @test PATH.c_api_License_SetString("bad_license") != 0
    @test PATH.c_api_Path_CheckLicense(1, 1) > 0
end

@testset "PathVersion" begin
    s = PATH.c_api_Path_Version()
    @test match(r"Path [0-9]\.[0-9]\.[0-9][0-9]", s) !== nothing
end

@testset "Example" begin
    M = convert(
        SparseArrays.SparseMatrixCSC{Cdouble, Cint},
        SparseArrays.sparse([
            0  0 -1 -1;
            0  0  1 -2;
            1 -1  2 -2;
            1  2 -2  4
        ])
    )
    status, z, info = PATH.solve_mcp(
        M,
        Float64[2, 2, -2, -6],
        fill(0.0, 4),
        fill(10.0, 4),
        [0.0, 0.0, 0.0, 0.0];
        output = "no",
    )
    @test status == PATH.MCP_Solved
    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])
end

@testset "Example LUSOL" begin
    M = convert(
        SparseArrays.SparseMatrixCSC{Cdouble, Cint},
        SparseArrays.sparse([
            0  0 -1 -1;
            0  0  1 -2;
            1 -1  2 -2;
            1  2 -2  4
        ])
    )
    status, z, info = PATH.solve_mcp(
        M,
        Float64[2, 2, -2, -6],
        fill(0.0, 4),
        fill(10.0, 4),
        [0.0, 0.0, 0.0, 0.0];
        output = "yes",
        # TODO(odow): when enabled, I get segfaults :(
        factorization_method = "blu_lusol",
        factorization_library_name = PATH.LUSOL_LIBRARY_PATH,
    )
    @test status == PATH.MCP_Solved
    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])
end

@testset "Example II" begin
    M = [
        0  0 -1 -1;
        0  0  1 -2;
        1 -1  2 -2;
        1  2 -2  4
    ]

    q = [2; 2; -2; -6]

    function F(n::Cint, x::Vector{Cdouble}, f::Vector{Cdouble})
        @assert n == length(x) == length(f) == 4
        f[1] = -x[3]^2 - x[4] + q[1]
        f[2] =  x[3]^3 - 2x[4]^2 + q[2]
        f[3] = x[1]^5 - x[2] + 2x[3] - 2x[4] + q[3]
        f[4] = x[1] + 2x[2]^3 - 2x[3] + 4x[4] + q[4]
        return Cint(0)
    end

    function J(
        n::Cint,
        nnz::Cint,
        x::Vector{Cdouble},
        col::Vector{Cint},
        len::Vector{Cint},
        row::Vector{Cint},
        data::Vector{Cdouble}
    )
        JAC = [
            0       0     -2x[3]      -1;
            0       0    3x[3]^2  -4x[4];
            5x[1]^4      -1          2      -2;
            1   6x[2]         -2       4
        ]
        @assert n == length(x) == length(col) == length(len) == 4
        @assert nnz == length(row) == length(data)
        i = 1
        for c in 1:n
            col[c] = i
            len[c] = 0
            for r in 1:n
                if !iszero(JAC[r, c])
                    data[i] = JAC[r, c]
                    row[i] = r
                    len[c] += 1
                    i += 1
                end
            end
        end
        return Cint(0)
    end

    status, z, info = PATH.solve_mcp(
        F,
        J,
        fill(0.0, 4),
        fill(10.0, 4),
        [1.0, 1.0, 1.0, 1.0];
        output = "yes"
    )
    @test status == PATH.MCP_Solved
    @test isapprox(z, [1.28475, 0.972916, 0.909376, 1.17304], atol=1e-4)
end


text:
# Implementation of the CouchDB CRUD API.
#
# This follows the example set by PouchDB of implementing all write operations via the
# CouchDB _bulk_docs end point, instead of the mess of a mixture of POST, PUT and DELETE.
# The advantage of this is much simpler code paths.
#
# There are some side effects of this, notably that some parameters that the other end points
# support isn't supported here.

import Base.isempty

"""
    immutable Database
      url
      name
      client

      Database(client::Client, name::AbstractString) =
        new(URI(client.url.scheme, client.url.host, client.url.port, "/\$name"), name, client)
    end

The Database immutable is a client connection tied to a specific remote DB. It is
normally not created directly, but via a call to `connectdb()`, or `createdb()`.

### Examples

    # Connect to existing DB. Does not verify it exists.
    db = connectdb(client; database="mydb")

    # Create a new db if it doesn't exist, otherwise connect
    db, created = createdb(client; database="mydb")
"""
immutable Database
  url
  name
  client

  Database(client::Client, name::AbstractString) =
    new(URI(client.url.scheme, client.url.host, client.url.port, "/$name"), name, client)
end

"""
    result = bulkdocs(db::Database; data=[], options=Dict())

Raw _bulk_docs.

This is a function primarily intended for internal use, but can
be used directly to create, update or delete documents in bulk,
so as to save on the HTTP overhead.

[API reference](http://docs.couchdb.org/en/1.6.1/api/database/bulk-api.html?#post--db-_bulk_docs)
"""
function bulkdocs(db::Database, data::AbstractArray; options=Dict())
  post_url = endpoint(db.url, "_bulk_docs")
  relax(post, post_url; json=Dict("docs" => data), cookies=db.client.cookies, query=options)
end

"""
    result = createdoc(db::Database, body::Dict)

Create one new document.

Note that this is implemented via the `_bulk_docs` endpoint, rather
than a `POST` to the `/{DB}`.

[API reference](http://docs.couchdb.org/en/1.6.1/api/database/bulk-api.html?#post--db-_bulk_docs)
"""
function createdoc(db::Database, body::Dict)
  if isempty(body)
    error("No body given")
  end
  response = bulkdocs(db, [body])
  response[1]
end

"""
    result = createdoc(db::Database, data::AbstractArray

Bulk create a set of new documents via the CouchDB `_bulk_docs` endpoint.

[API reference](http://docs.couchdb.org/en/1.6.1/api/database/bulk-api.html?#post--db-_bulk_docs)
"""
function createdoc(db::Database, data::AbstractArray)
  if isempty(data)
    error("No data given")
  end

  bulkdocs(db, data)
end

"""
    result = readdoc(db::Database, id::AbstractString;
      rev               = "",
      attachments       = false,
      att_encoding_info = false,
      atts_since        = [],
      open_revs         = [],
      conflicts         = false,
      deleted_conflicts = false,
      latest            = false,
      meta              = false,
      revs              = false,
      revs_info         = false)

Fetch a document by `id`.

For a description of the parameters, see reference below. To use the `open_revs` parameter as `all`, use

    result = readdoc(db, id; open_revs=["all"])

[API reference](http://docs.couchdb.org/en/1.6.1/api/document/common.html#get--db-docid)
"""
function readdoc(db::Database, id::AbstractString;
  rev               = "",
  attachments       = false,
  att_encoding_info = false,
  atts_since        = [],
  open_revs         = [],
  conflicts         = false,
  deleted_conflicts = false,
  latest            = false,
  meta              = false,
  revs              = false,
  revs_info         = false)

  query::Dict{AbstractString, Any} = Dict()
  headers::Dict{AbstractString, Any} = Dict("Accept" => "application/json")

  if attachments
    query["attachments"] = true
  end

  if att_encoding_info
    query["att_encoding_info"] = true
  end

  if conflicts
    query["conflicts"] = true
  end

  if deleted_conflicts
    query["deleted_conflicts"] = true
  end

  if latest
    query["latest"] = true
  end

  if meta
    query["meta"] = true
  end

  if revs
    query["revs"] = true
  end

  if revs_info
    query["revs_info"] = true
  end

  if rev != ""
    query["rev"] = rev
  end

  if length(atts_since) > 0
    query["atts_since"] = JSON.json(atts_since)
  end

  if length(open_revs) > 0
    if length(open_revs) == 1 && open_revs[1] == "all"
      query["open_revs"] = "all"
    else
      query["open_revs"] = JSON.json(open_revs)
    end
  end

  relax(get, endpoint(db.url, id); cookies=db.client.cookies, query=query, headers=headers)
end

"""
    result = updatedoc(db::Database; id::AbstractString=nothing, rev::AbstractString=nothing, body=Dict())

Update an existing document, creating a new revision.

Implemented via the _bulk_docs endpoint.

[API reference](http://docs.couchdb.org/en/1.6.1/api/database/bulk-api.html?#post--db-_bulk_docs)
"""
function updatedoc(db::Database; id::AbstractString=nothing, rev::AbstractString=nothing, body=Dict())
  response = bulkdocs(db, [merge(body, Dict("_id" => id, "_rev" => rev))])
  response[1]
end

"""
    result = deletedoc(db::Database; id::AbstractString=nothing, rev::AbstractString=nothing)

Delete a document revision. Implemented via the _bulk_docs endpoint:

[API reference](http://docs.couchdb.org/en/1.6.1/api/database/bulk-api.html?#post--db-_bulk_docs)
"""
function deletedoc(db::Database; id::AbstractString=nothing, rev::AbstractString=nothing)
  response = bulkdocs(db, [Dict("_id" => id, "_rev" => rev, "_deleted" => true)])
  response[1]
end


text:
using JuMP, Ipopt

function tcm_estimateMNL(N, assortments, transaction_counts)

	M = size(assortments, 1)

	total_trans_by_assortment = sum(transaction_counts, 2)

	m_MNL = Model(solver = IpoptSolver(print_level=1))

	@variable(m_MNL, -20 <= u[1:N-1] <= 20)

	@NLobjective(m_MNL, Max, sum(sum( transaction_counts[m,p] * u[p] for p in 1:N-1) for m in 1:M) - sum(total_trans_by_assortment[m]*log(1+sum(assortments[m,p]* exp(u[p]) for p in 1:N-1)) for m in 1:M) )

	start_time = time()
	solve(m_MNL)
	end_time = time()
	elapsed_time = end_time - start_time

	u_val = getvalue(u)

	obj = getobjectivevalue(m_MNL)

	return u_val, obj, elapsed_time
end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense