text:
export BrickletLCD16x2Config
struct BrickletLCD16x2Config
    cursor::Bool
    blinking::Bool
end

export BrickletLCD16x2Identity
struct BrickletLCD16x2Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletLCD16x2
"""
16x2 character alphanumeric display with blue backlight
"""
mutable struct BrickletLCD16x2 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletLCD16x2(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_lcd_16x2")
        deviceInternal = package.BrickletLCD16x2(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export write_line
"""
    $(SIGNATURES)

Writes text to a specific line with a specific position.
The text can have a maximum of 16 characters.

For example: (0, 5, "Hello") will write *Hello* in the middle of the
first line of the display.

The display uses a special charset that includes all ASCII characters except
backslash and tilde. The LCD charset also includes several other non-ASCII characters, see
the `charset specification <https://github.com/Tinkerforge/lcd-16x2-bricklet/raw/master/datasheets/standard_charset.pdf>`__
for details. The Unicode example above shows how to specify non-ASCII characters
and how to translate from Unicode to the LCD charset.
"""
function write_line(device::BrickletLCD16x2, line, position, text)
    device.deviceInternal.write_line(line, position, text)
end

export clear_display
"""
    $(SIGNATURES)

Deletes all characters from the display.
"""
function clear_display(device::BrickletLCD16x2)
    device.deviceInternal.clear_display()
end

export backlight_on
"""
    $(SIGNATURES)

Turns the backlight on.
"""
function backlight_on(device::BrickletLCD16x2)
    device.deviceInternal.backlight_on()
end

export backlight_off
"""
    $(SIGNATURES)

Turns the backlight off.
"""
function backlight_off(device::BrickletLCD16x2)
    device.deviceInternal.backlight_off()
end

export is_backlight_on
"""
    $(SIGNATURES)

Returns *true* if the backlight is on and *false* otherwise.
"""
function is_backlight_on(device::BrickletLCD16x2)
    return device.deviceInternal.is_backlight_on()
end

export set_config
"""
    $(SIGNATURES)

Configures if the cursor (shown as "_") should be visible and if it
should be blinking (shown as a blinking block). The cursor position
is one character behind the the last text written with
:func:`Write Line`.
"""
function set_config(device::BrickletLCD16x2, cursor, blinking)
    device.deviceInternal.set_config(cursor, blinking)
end

export get_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Config`.
"""
function get_config(device::BrickletLCD16x2)
    return device.deviceInternal.get_config()
end

export is_button_pressed
"""
    $(SIGNATURES)

Returns *true* if the button is pressed.

If you want to react on button presses and releases it is recommended to use the
:cb:`Button Pressed` and :cb:`Button Released` callbacks.
"""
function is_button_pressed(device::BrickletLCD16x2, button)
    return device.deviceInternal.is_button_pressed(button)
end

export set_custom_character
"""
    $(SIGNATURES)

The LCD 16x2 Bricklet can store up to 8 custom characters. The characters
consist of 5x8 pixels and can be addressed with the index 0-7. To describe
the pixels, the first 5 bits of 8 bytes are used. For example, to make
a custom character "H", you should transfer the following:

* ``character[0] = 0b00010001`` (decimal value 17)
* ``character[1] = 0b00010001`` (decimal value 17)
* ``character[2] = 0b00010001`` (decimal value 17)
* ``character[3] = 0b00011111`` (decimal value 31)
* ``character[4] = 0b00010001`` (decimal value 17)
* ``character[5] = 0b00010001`` (decimal value 17)
* ``character[6] = 0b00010001`` (decimal value 17)
* ``character[7] = 0b00000000`` (decimal value 0)

The characters can later be written with :func:`Write Line` by using the
characters with the byte representation 8 ("\\\\x08" or "\\\\u0008") to 15
("\\\\x0F" or "\\\\u000F").

You can play around with the custom characters in Brick Viewer since
version 2.0.1.

Custom characters are stored by the LCD in RAM, so they have to be set
after each startup.

.. versionadded:: 2.0.1\$nbsp;(Plugin)
"""
function set_custom_character(device::BrickletLCD16x2, index, character)
    device.deviceInternal.set_custom_character(index, character)
end

export get_custom_character
"""
    $(SIGNATURES)

Returns the custom character for a given index, as set with
:func:`Set Custom Character`.

.. versionadded:: 2.0.1\$nbsp;(Plugin)
"""
function get_custom_character(device::BrickletLCD16x2, index)
    return device.deviceInternal.get_custom_character(index)
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletLCD16x2)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletLCD16x2, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
using Dash
using DashHtmlComponents

app = dash()

app.layout = html_div() do
    html_div("Hello Dash.jl testing", id="container")
end

run_server(app)

ext: jl


text:
struct RGB
    r::UInt8
    g::UInt8
    b::UInt8

    function RGB(r::Int64,g::Int64,b::Int64)
        @assert 0 <= r <= 255 || 0 <= g <= 255 || 0 <= b <= 255 "wrong input for a color. Float input should be between [0,1], Integer input should be between [0,255]"
        new(r,g,b)
    end
end

struct Color
    rgb::RGB
    opacity::Real
    function Color(rgb::RGB, opacity::Real=1)
        @assert 0 <= opacity <= 1 "wrong choice of opacity"
        new(rgb, opacity)
    end
end

const default_color = (white   = Color(RGB(255, 255, 255)),
                       black   = Color(RGB(0  , 0  , 0  )),
                       red     = Color(RGB(255 ,0   ,0  )),
                       lime    = Color(RGB(0  , 255 ,0  )),
                       blue    = Color(RGB(0  , 0   ,255)),
                       yellow  = Color(RGB(255, 255, 0  )),
                       magenta = Color(RGB(255 ,0   ,255)),
                       cyan    = Color(RGB(0   ,255 ,255)),
                       orange  = Color(RGB(255 ,128, 0  )),
                       gray    = Color(RGB(128, 128, 128)),
                       silver  = Color(RGB(191, 191, 191)),
                       pink    = Color(RGB(255, 128, 128)),
                       maroon  = Color(RGB(144, 53,  59 )),
                       olive   = Color(RGB(128, 128, 0  )),
                       green   = Color(RGB(0  , 128, 0  )),
                       purple  = Color(RGB(128, 0  , 128)),
                       teal    = Color(RGB(110, 142, 132)),
                       navy    = Color(RGB(26 , 71 , 111)),
                       forest  = Color(RGB(85 , 117, 47 )),
                       dkorg   = Color(RGB(227, 126, 0  )),
                       cranbrry= Color(RGB(193, 5  , 52 )),
                       lavender= Color(RGB(147, 141, 210)),
                       khaki   = Color(RGB(202, 194, 126)),
                       sienna  = Color(RGB(160, 82 , 45 )),
                       emidblue= Color(RGB(123, 146, 168)),
                       emerald = Color(RGB(45 , 109, 102)),
                       brown   = Color(RGB(156, 136, 71 )),
                       erose   = Color(RGB(191, 161, 156)),
                       gold    = Color(RGB(255, 210,   0)),
                       bluegray= Color(RGB(217, 230, 235)),
                       jl_blue = Color(RGB(64 , 99 , 216)),
                       jl_purp = Color(RGB(149, 88 , 178)),
                       jl_green= Color(RGB(56 , 152, 38 )),
                       jl_red  = Color(RGB(203, 60 , 51 )))

const stata_color = [:navy, :maroon, :forest, :dkorg, :teal, :cranbrry, :lavender, :khaki,
                     :sienna, :emidblue, :emerald, :brown, :erose, :gold, :bluegray]

const julia_color = [:jl_blue, :jl_green, :jl_red, :jl_purp]

COLOR_PALATTE = Iterators.cycle(broadcast(x->get(default_color,x,:black),julia_color))

changeopacity(color::Color,o::Real) = Color(color.rgb, o)

rbg_as_xcolor(rgb::RGB) = "rgb,255: red, $(rgb.r); green, $(rgb.g); blue, $(rgb.b)"
color_as_xcolor(color::Color) = rbg_as_xcolor(color.rgb)
get_default_color(color_name::Symbol) = color_name in keys(default_color) ? default_color[color_name] : throw(ArgumentError("$(color_name) is not defined."))
colorname_as_xcolor(color_name::Symbol) = color_as_xcolor(get_default_color(color_name))

color_as_option(color::Color) = PGFPlotsX.Options(:color => color_as_xcolor(color),
                                                  :opacity => color.opacity)
color_as_fill_option(color::Color) = PGFPlotsX.Options(:fill => color_as_xcolor(color),
                                                       Symbol("fill opacity") => color.opacity)
color_as_draw_option(color::Color) = PGFPlotsX.Options(:draw => color_as_xcolor(color),
                                                       Symbol("draw opacity") => color.opacity)

color_as_option(color_name::Symbol) = PGFPlotsX.Options(:color => colorname_as_xcolor(color_name),
                                                        :opacity => get_default_color(color_name).opacity)
color_as_fill_option(color_name::Symbol) = PGFPlotsX.Options(:fill => colorname_as_xcolor(color_name),
                                                             Symbol("fill opacity") => get_default_color(color_name).opacity)
color_as_draw_option(color_name::Symbol) = PGFPlotsX.Options(:draw => colorname_as_xcolor(color_name),
                                                             Symbol("draw opacity") => get_default_color(color_name).opacity)


text:
# TODO: Implement using adjacency lists.
# TODO: Implement using watched literals.
function unit_propagation!(fla::Formula, pa::PartialAssignment)::Bool
  while true
    l = find_unit_literal(fla, pa)
    if 0 == l
      break
    end

    pa[abs(l)] = sign(l)
    if any((cl) -> is_clause_falsified(cl, pa), fla)
      return false
    end
  end

  return !any((cl) -> is_clause_falsified(cl, pa), fla)
end


text:
using PyCall
enmap, curvedsky = pyimport("pixell.enmap"), pyimport("pixell.curvedsky")

##
using PowerSpectra
using Pixell
shape, wcs = fullsky_geometry(20 * Pixell.arcminute)

Dl = PowerSpectra.planck_theory_Dl()
channelindex(X) = findfirst(first(X), "TEB")

ells = collect(eachindex(Dl[:TT]))
Dlfac = ells .* (ells .+1 ) ./ (2π)

𝐂 = zeros(3, 3, length(ells))
for X in ("T", "E", "B"), Y in ("T", "E", "B")
    c₁, c₂ = channelindex(X), channelindex(Y)
    XY = Symbol(X * Y)
    if XY in keys(Dl)
        𝐂[c₁, c₂ ,:] .= parent(Dl[XY]) ./ Dlfac
        𝐂[ c₂, c₁ ,:] .= 𝐂[c₁, c₂ ,:]
    end
end

lmax = Pixell.getlmax(wcs)
alms = (Alm(lmax, lmax), Alm(lmax, lmax), Alm(lmax, lmax))
synalm!(𝐂, alms; lmin=2, lmax=lmax)

##
m = alm2map(alms[1], shape, wcs)
plot(m)

##
I, Q, U = alm2map(alms, shape, wcs)
aT, aE, aB = map2alm((I, Q, U); lmax=lmax)
ClTT = alm2cl(aT)
ClEE = alm2cl(aE)
ClTE = alm2cl(aT, aE)
ClBB = alm2cl(aB)

IQU = cat(I, Q, U, dims=3)
write_map("data/IQU.fits", IQU)

##
aT.lmax
##
using DelimitedFiles

open("data/TEB_alms_real.dat", "w") do io
    writedlm(io, [real.(a.alm) for a in alms])
end;
open("data/TEB_alms_imag.dat", "w") do io
    writedlm(io, [imag.(a.alm) for a in alms])
end;

##
pyI = pycall(enmap.read_map, PyObject, "data/IQU.fits", sel=PyObject(0:0))  # get WCS
newpymap = pycall(enmap.enmap, PyObject, pyI, pyI.wcs)
pyI = pycall(curvedsky.alm2map, PyObject, PyObject(alms[1].alm), newpymap, spin=0)

pyQU = pycall(enmap.read_map, PyObject, "data/IQU.fits", sel=PyObject(1:2))
newpymap = pycall(enmap.enmap, PyObject, pyQU, pyQU.wcs)
pyQUalm = [alms[2].alm  alms[3].alm]
pyQU = pycall(curvedsky.alm2map, PyObject, PyReverseDims(pyQUalm), newpymap, spin=2)

enmap.write_map("data/pyI.fits", pyI)
enmap.write_map("data/pyQU.fits", pyQU)

pyI = pycall(enmap.read_map, PyObject, "data/IQU.fits", sel=PyObject(0:0))
pyQU = pycall(enmap.read_map, PyObject, "data/IQU.fits", sel=PyObject(1:2))

# islice = py"lambda x : x[3,:,:]"
# quslice = py"lambda x : x[3,:,:]"

# pyI = pycall(islice, PyObject, pyIQU)

##
pyT = curvedsky.map2alm(pyI, lmax=lmax, spin=0)
pyEB = curvedsky.map2alm(pyQU, lmax=lmax, spin=2)
pyE = pyEB[1,:]
pyB = pyEB[2,:]

aT, aE, aB = map2alm((I, Q,U); lmax=lmax)
hp = pyimport("healpy")

pyClTT = (hp.alm2cl(pyT)')[:,1]
pyClTE = (hp.alm2cl(pyT, pyE)')[:,1]
pyClEE = hp.alm2cl(pyE)
pyClBB = hp.alm2cl(pyB)

##
using DelimitedFiles

data = [pyClTT pyClTE pyClEE pyClBB]
open("data/test_cls_IQU.txt", "w") do io
    writedlm(io, data)
end;


text:
using Random: MersenneTwister, randperm

# helper function that permutates the vertices of a graph
function shuffle_vertices(g::AbstractGraph, σ)
    G = typeof(g)
    gperm = G(nv(g))
    n = nv(g)
    for e in edges(g)
        add_edge!(gperm, σ[src(e)], σ[dst(e)])
    end
    return gperm
end

# all (up to isomorphism) 3-regular, connected graphs on 8 vertices
# list taken from https://en.wikipedia.org/wiki/Table_of_simple_cubic_graphs#4_vertices
cubic_graphs = [SimpleGraph(Edge.([(1,2), (1,3), (1,8), (2,3), (2,4), (3,4),
                                   (4,5), (5,6), (5,7), (6,7), (6,8), (7,8)])),
                SimpleGraph(Edge.([(1,2), (1,3), (1,8), (2,3), (2,5), (3,4),
                                   (4,5), (4,7), (5,6), (6,7), (6,8), (7,8)])),
                SimpleGraph(Edge.([(1,2), (1,3), (1,8), (2,3), (2,6), (3,4),
                                   (4,5), (4,7), (5,6), (5,8), (6,7), (7,8)])),
                grid([2,2,2]),
                SimpleGraph(Edge.([(1,2), (1,5), (1,8), (2,3), (2,6), (3,4),
                                   (3,7), (4,5), (4,8), (5,6), (6,7), (7,8)]))
               ]
# vertex permutated versions of cubic_graphs
cubic_graphs_perm = []
_rng = MersenneTwister(1)
for i = 1:length(cubic_graphs)
    push!(cubic_graphs_perm, shuffle_vertices(cubic_graphs[i], randperm(_rng, 8)))
end

@testset "Isomorphism" begin
    @test has_isomorph(grid([2,3]), grid([3,2]))

    # the cubic graphs should only be isomorph to themself
    # the same holds for subgraph isomorphism and induced subgraph isomorphism
    for i in 1:length(cubic_graphs), j in 1:length(cubic_graphs)
        @test (i == j) == has_isomorph(cubic_graphs[i], cubic_graphs[j])
        @test (i == j) == has_subgraphisomorph(cubic_graphs[i], cubic_graphs[j])
        @test (i == j) == has_induced_subgraphisomorph(cubic_graphs[i], cubic_graphs[j])
    end

    # the cubic graphs should only be isomorph a permutation of themself
    # the same holds for subgraph isomorphism and induced subgraph isomorphism
    for i in 1:length(cubic_graphs), j in 1:length(cubic_graphs_perm)
        @test (i == j) == has_isomorph(cubic_graphs[i], cubic_graphs_perm[j])
        @test (i == j) == has_subgraphisomorph(cubic_graphs[i], cubic_graphs_perm[j])
        @test (i == j) == has_induced_subgraphisomorph(cubic_graphs[i], cubic_graphs_perm[j])
    end

    # count_isomorph, count_subgraphisomorph and count_induced_subgraphisomorph are commutative
    for i in 1:length(cubic_graphs)
        g1 = cubic_graphs[i]
        g2 = cubic_graphs_perm[i]
        @test count_isomorph(g1, g1) == count_isomorph(g1, g2) == count_isomorph(g2, g1) == count_isomorph(g2,g2)
        @test count_subgraphisomorph(g1, g1) == count_subgraphisomorph(g1, g2) ==
                   count_subgraphisomorph(g2, g1) == count_subgraphisomorph(g2,g2)
        @test count_induced_subgraphisomorph(g1, g1) == count_subgraphisomorph(g1, g2) ==
                   count_subgraphisomorph(g2, g1) == count_subgraphisomorph(g2,g2)
    end

    for i in 1:length(cubic_graphs)
        g1 = cubic_graphs[i]
        g2 = cubic_graphs_perm[i]
        length(collect(all_isomorph(g1, g1))) == count_isomorph(g1, g2)
        length(collect(all_subgraphisomorph(g1, g1))) == count_isomorph(g1, g2)
        length(collect(all_induced_subgraphisomorph(g1, g1))) == count_isomorph(g1, g2)
    end

    # The frucht-graph is a 3-regular graph that has only one automorphism
    g = smallgraph(:frucht)
    perm = randperm(_rng, 12)
    g2 = shuffle_vertices(g, perm)
    perm_reconstructed = map(uv -> first(uv), first(all_isomorph(g2, g)))
    @test perm == perm_reconstructed

    # TODO better tests for vertex_relation and edge_relation
    vrel(u, v) = false
    erel(e1, e2) = false
    @test has_isomorph(complete_graph(4), complete_graph(4), vertex_relation=vrel) == false
    @test has_isomorph(complete_graph(4), complete_graph(4), edge_relation=erel) == false
    @test has_subgraphisomorph(complete_graph(4), complete_graph(3), vertex_relation=vrel) == false
    @test has_subgraphisomorph(complete_graph(4), complete_graph(3), edge_relation=erel) == false
    @test has_induced_subgraphisomorph(complete_graph(4), complete_graph(3), vertex_relation=vrel) == false
    @test has_induced_subgraphisomorph(complete_graph(4), complete_graph(3), edge_relation=erel) == false

    @test count_isomorph(complete_graph(4), complete_graph(4), vertex_relation=vrel) == 0
    @test count_isomorph(complete_graph(4), complete_graph(4), edge_relation=erel) == 0
    @test count_subgraphisomorph(complete_graph(4), complete_graph(3), vertex_relation=vrel) == 0
    @test count_subgraphisomorph(complete_graph(4), complete_graph(3), edge_relation=erel) == 0
    @test count_induced_subgraphisomorph(complete_graph(4), complete_graph(3), vertex_relation=vrel) == 0
    @test count_induced_subgraphisomorph(complete_graph(4), complete_graph(3), edge_relation=erel) == 0

    @test isempty(all_isomorph(complete_graph(4), complete_graph(4), vertex_relation=vrel))
    @test isempty(all_isomorph(complete_graph(4), complete_graph(4), edge_relation=erel))
    @test isempty(all_subgraphisomorph(complete_graph(4), complete_graph(3), vertex_relation=vrel))
    @test isempty(all_subgraphisomorph(complete_graph(4), complete_graph(3), edge_relation=erel))
    @test isempty(all_induced_subgraphisomorph(complete_graph(4), complete_graph(3), vertex_relation=vrel))
    @test isempty(all_induced_subgraphisomorph(complete_graph(4), complete_graph(3), edge_relation=erel))

    # some test for early returns if there is no isomorphism
    @test count_isomorph(complete_graph(4), cycle_graph(4)) == 0
    @test isempty(all_isomorph(complete_graph(4), cycle_graph(4)))
    @test count_subgraphisomorph(complete_graph(3), complete_graph(4)) == 0

    # this tests triggers the shortcut in the vf2 algorithm if the first graph is smaller than the second one
    @test has_isomorph(complete_graph(3), complete_graph(4)) == false

    # this test is mainly to cover a certain branch of vf2
    # TODO cover directed graphs better
    gd = SimpleDiGraph(Edge.([(2,1)]))
    @test has_isomorph(gd, gd, VF2()) == true

    # Test for correct handling of self-loops
    g1 = SimpleGraph(1)
    add_edge!(g1, 1, 1)
    g2 = SimpleGraph(1)
    @test has_isomorph(g1, g1) == true
    @test has_isomorph(g1, g2) == false
    @test has_isomorph(g2, g1) == false
    @test has_isomorph(g1, g1) == true

    @test has_induced_subgraphisomorph(g1, g1) == true
    @test has_induced_subgraphisomorph(g1, g2) == false
    @test has_induced_subgraphisomorph(g2, g1) == false
    @test has_induced_subgraphisomorph(g2, g2) == true

    @test has_subgraphisomorph(g1, g1) == true
    @test has_subgraphisomorph(g1, g2) == true
    @test has_subgraphisomorph(g2, g1) == false
    @test has_subgraphisomorph(g2, g2) == true
end


text:
function readvtu(filename)
    points = []
    tetrahedra = Array{Int64,1}[]

    open(filename) do fid
        while !eof(fid)
            line = readline(fid)
            #println(idx,"::",line)
            line=split(line)
            if line[1]=="<Points>"
                line=readline(fid)
                line=split(readline(fid))

                while line[1][1]!='<'
                    numbers=parse.(Float64,line)
                    append!(points,[numbers[1:3]])
                    if length(numbers)==6
                        append!(points,[numbers[4:6]])
                    end
                    line=split(readline(fid))
                end
            end

            if line[1]=="<Cells>"
                line=readline(fid)
                line=split(readline(fid))
                tet=[]
                while line[1][1]!='<'
                    numbers=parse.(Int,line)
                    while length(numbers)!=0 && length(tet)!=4
                        append!(tet,popfirst!(numbers))
                    end
                    append!(tetrahedra,[tet.+1])
                    tet=numbers[1:end]
                    if length(tet)==4
                        append!(tetrahedra,[tet.+1])
                        tet=[]
                    end
                    line=split(readline(fid))
                end
            end

        end
    end
    return points, tetrahedra
end

points, tetrahedra = readvtu("mesh.vtu")

import ProgressMeter

p = ProgressMeter.Progress(length(tetrahedra),desc="Assemble all triangles ", dt=1,
     barglyphs=ProgressMeter.BarGlyphs('|','█', ['▁' ,'▂' ,'▃' ,'▄' ,'▅' ,'▆', '▇'],' ','|',),
     barlen=20)
triangles=[]
simplices=[]
for (idx,tet) in enumerate(tetrahedra)
    append!(simplices,idx)

    for tri in (tet[[1,2,3]],tet[[1,2,4]],tet[[1,3,4]],tet[[2,3,4]])
        idx,ins=sort_smplx(triangles,tri)
        if ins
            insert!(triangles,idx,tri)
        else
            deleteat!(triangles,idx)
        end
    end
    ProgressMeter.next!(p)
end
# convert points

Points=zeros(Float64,3,length(points))
for (idx,pnt) in enumerate(points)
    Points[:,idx]=pnt
end
#
x_min=minimum(Points[1,:])
x_max=maximum(Points[1,:])
#
domains=Dict()
domain=Dict()
domain["simplices"]=simplices
domain["dimension"]=0x00000003
domains["Interior"]=domain
#
x_min=minimum(Points[1,:])
x_max=maximum(Points[1,:])
simplices=[]
for (tri_idx,tri) in enumerate(triangles)
    if all(abs.(Points[1,tri].-x_max).<=0.0001)
        insert_smplx!(simplices,tri_idx)
    end
end
domain=Dict()
domain["simplices"]=simplices
domain["dimension"]=0x00000002
domains["Outlet"]=domain
#
tri2tet=zeros(UInt32,length(triangles))
tri2tet[1]=0xffffffff

mesh=Mesh("mesh.vtu",Points,[],triangles,tetrahedra,domains,"mesh.vtu",tri2tet,1)
vtk_write("eth_diff", mesh, data_diff)


text:
using Plots
R = 0.0821
T = 298
isoterm(x) = R * T / x
v = 5 : -0.1 : 2.5
p = isoterm.(v)
plot(xlabel = "V[L]", ylabel = "P[atm]", xlims = xlims = (0, 7),
		ylims = ylims = (0, 14),
		xticks = 0 : 0.5 : xlims[2], yticks = 0 : 1 : ylims[2])
plot!(v, p;
		color = :red, label = "A: Trinn 1", arrow = :arrow, linewidth = 2)
plot!([v[end], v[end]], [p[end], p[1]];
		color = :blue, label = "A: Trinn 2", arrow = :arrow, linewidth = 2)
plot!([v[1], v[end]], [p[1], p[1]];
		color = :green, label = "B: Trinn 1", arrow = :arrow, linewidth = 2)
#


text:
using ContextualBandits
using Documenter

DocMeta.setdocmeta!(ContextualBandits, :DocTestSetup, :(using ContextualBandits); recursive=true)

makedocs(;
    modules=[ContextualBandits],
    authors="Douglas Corbin <dfcorbin98@gmail.com>",
    repo="https://github.com/dfcorbin/ContextualBandits.jl/blob/{commit}{path}#{line}",
    sitename="ContextualBandits.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://dfcorbin.github.io/ContextualBandits.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/dfcorbin/ContextualBandits.jl",
    devbranch="main",
)


text:
"""
Created on Thu Apr 10 2020

@author: Yoann Pradat

Functions for initializing W and H matrices in the NMF V = WH.
"""

"""
    initialize_nmf(V[, rank, init, dist, ϵ, seed])

Compute initial guess for the non-negative rank K approximation of V: V = WH.
# Arguments
- `V::Matrix{<:Real}`   : matrix, shape (n_features, n_samples)
- `rank::Integer`       : integer, number of components desired
- `init::Symbol`        : symbol,  Method used to initialize the procedure. Valid options are:
    + :default: :nndsvd if rank <= min(n_samples, n_features), otherwise :random.
    + :random: Nonnegative random matrices, scaled with: sqrt(mean(V) / rank)
    + :nndsvd: Nonnegative Double Singular Value Decomposition (NNDSVD) initialization
    (better for sparseness)
    + :nndsvda: NNDSVD with zeros filled with the average of V (better when sparsity is
    not desired)
    + :nndsvdar: NNDSVD with zeros filled with small random values (generally faster, less
    accurate alternative to NNDSVDa for when sparsity is not desired)
- `dist::Distribution` : the distribution from which random entries are drawn. Be aware
    that entries will be multiplied by sqrt(mean(V) / rank). Used for intializations with randomness
- `seed:Integer`       : fix the seed of RNG
- `ϵ::<:Real`          : threshold below which values are assigned to 0

# References
- C. Boutsidis, E. Gallopoulos: SVD based initialization: A head start for nonnegative
matrix factorization - Pattern Recognition, 2007
"""
function initialize_nmf(V::Matrix{T}; rank::Integer=2, init::Symbol=:default, dist::Distribution=truncated(Normal(0,1), 0, Inf), ϵ::Real=1e-6, seed::Integer=0) where T <: Real
    rng = Random.seed!(seed)

    if init == :default
        n_features, n_samples = size(V)
        if rank <= min(n_features, n_samples)
            init = :nndsvd
        else
            init = :random
        end
    end

    if init == :random
        W, H = _initialize_random(V, rank=rank, dist=dist, rng=rng)
    elseif init == :nndsvd
        W, H = _initialize_nndsvd(V, rank=rank)
        W[W .< ϵ] .= 0
        H[H .< ϵ] .= 0
    elseif init == :nndsvda
        W, H = _initialize_nndsvd(V, rank=rank)
        W[W .< ϵ] .= 0
        H[H .< ϵ] .= 0

        avg_V = mean(V)
        W[W .== 0] .= avg_V
        H[H .== 0] .= avg_V
    elseif init == :nndsvdar
        W, H = _initialize_nndsvd(V, rank=rank)
        W[W .< ϵ] .= 0
        H[H .< ϵ] .= 0

        avg_V = mean(V)
        W[W .== 0] .= avg_V * rand(rng, dist, size(W[W .== 0], 1)) ./ 100
        H[H .== 0] .= avg_V * rand(rng, dist, size(H[H .== 0], 1)) ./ 100

        if minimum(W) < 0 || minimum(H) < 0
            throw(ArgumentError("Please specify a distribution that does not produce negative values. Use `truncated()` if necessary"))
        end
    else
        throw(ArgumentError("""Invalid init parameter: choose one of :default, :random, :nndsvd, :nndsvda, :nndsvdar """))
    end
    W, H
end

"""
    _initialize_random(V[, rank, dist, rng])

Initialize W and H by random matrices of size (n_features, rank) and (rank, n_samples).
Entries are drawn from half-normal distribution i.e the absolute value of a normally distributed random
variable. The variance  is chosen to be the mean(V) / rank.
"""
function _initialize_random(V::Matrix{T}; rank, dist::Distribution, rng::AbstractRNG) where T <: Real
    n_features, n_samples = size(V)
    avg = sqrt(mean(V) / rank)

    W = avg * rand(rng, dist, n_features, rank)
    H = avg * rand(rng, dist, rank, n_samples)

    if minimum(W) < 0 || minimum(H) < 0
        throw(ArgumentError("Please specify a distribution that does not produce negative values. Use `truncated()` if necessary"))
    end

    W,H
end

"""
    _initialize_nndsvd(X; rank)

Initialize matrices of the decomposition X = WH using SVD.

# References
- C. Boutsidis, E. Gallopoulos: SVD based initialization: A head start for
nonnegative matrix factorization - Pattern Recognition, 2007
"""
function _initialize_nndsvd(X::Matrix{T}; rank::Integer=2) where T <: Real
    n_features, n_samples = size(X)
    W = zeros(T, n_features, rank)
    H = zeros(T, rank, n_samples)

    #### Using Perron-Frobenius theorem, one can prove that the leading singular
    #### vectors in SVD can be chosen to have non negative coordinates.
    #### These can therefore be used for initialization.
    U, S, Vt = svd(X, full=false)
    V = transpose(Vt)

    W[:, 1] = sqrt(S[1]) * abs.(U[:, 1])
    H[1, :] = sqrt(S[1]) * abs.(V[1, :])

    for j = 2:rank
        x, y = U[:, j], V[j, :]

        # extract positive and negative parts of column vectors
        x_p, y_p = max.(x, 0), max.(y, 0)
        x_n, y_n = abs.(min.(x, 0)), abs.(min.(y, 0))

        # and their norms
        x_p_nrm, y_p_nrm = norm(x_p, 2), norm(y_p, 2)
        x_n_nrm, y_n_nrm = norm(x_n, 2), norm(y_n, 2)

        # multiply positive and negative norms
        m_p, m_n = x_p_nrm * y_p_nrm, x_n_nrm * y_n_nrm

        # choose update
        if m_p > m_n
            u = x_p / x_p_nrm
            v = y_p / y_p_nrm
            σ = m_p
        else
            u = x_n / x_n_nrm
            v = y_n / y_n_nrm
            σ = m_n
        end

        λ = sqrt(S[j] * σ)
        W[:, j] = λ * u
        H[j, :] = λ * v
    end

    W, H
end


text:
using Documenter
using IntervalArithmetic, IntervalRootFinding

makedocs(
    modules = [IntervalRootFinding],
    doctest = true,
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "IntervalRootFinding.jl",
    authors = "David P. Sanders and Luis Benet",
    pages = Any[
        "Home" => "index.md",
        "`roots` interface" => "roots.md",
        "Internals" => "internals.md",
        "Bibliography" => "biblio.md",
        "API" => "api.md"
        ]
    )

deploydocs(
    repo = "github.com/JuliaIntervals/IntervalRootFinding.jl.git",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
module Geometry

export Point, Points, Size, BoundingBox, xy, yx

struct Point{T}
    x::T
    y::T
end

xy(point::Point) = (point.x, point.y)

yx(point::Point) = (point.y, point.x)

Base.ceil(point::Point) = Point(ceil(point.x),
                                ceil(point.y))

Base.round(point::Point) = Point(round(point.x),
                                 round(point.y))

Base.:+(a::Point, b::Point) = Point(a.x + b.x,
                                    a.y + b.y)

Base.:-(a::Point, b::Point) = Point(a.x - b.x,
                                    a.y - b.y)

Base.:*(a::Point, b::Point) = Point(a.x * b.x,
                                    a.y * b.y)

Base.:/(a::Point, b::Int64) = Point(a.x / b,
                                    a.y / b)

function Base.ceil(::Type{T}, point::Point{B}) where {T, B}
    return Point(ceil(T, point.x),
                 ceil(T, point.y))
end

function Base.round(::Type{T}, point::Point{B}) where {T, B}
    return Point(round(T, point.x),
                 round(T, point.y))
end

const Points{T} = Vector{Point{T}}

struct Size{T}
    x::T
    y::T
end

Size(point::Point{T}) where {T} = Size(point.x, point.y)

Point(size::Size{T}) where {T} = Point(size.x, size.y)

xy(size::Size) = (size.x, size.y)

yx(size::Size) = (size.y, size.x)

Base.ceil(size::Size) = Size(ceil(size.x),
                             ceil(size.y))

Base.:+(a::Size, b::Size) = Size(a.x + b.x,
                                 a.y + b.y)

Base.:-(a::Size, b::Size) = Size(a.x - b.x,
                                 a.y - b.y)

Base.:*(a::Size, b::Size) = Size(a.x * b.x,
                                 a.y * b.y)

Base.:*(a::Size, b::Number) = Size(a.x * b,
                                   a.y * b)

Base.:/(a::Size, b::Int64) = Size(a.x / b,
                                  a.y / b)

# TODO: bad decision?
Base.:*(a::Size, b::Point) = Point(a.x * b.x,
                                   a.y * b.y)
Base.:*(a::Point, b::Size) = Point(a.x * b.x,
                                   a.y * b.y)

struct BoundingBox
    x_min::Float64
    x_max::Float64
    y_min::Float64
    y_max::Float64
end

Base.:+(a::BoundingBox, b::BoundingBox) = BoundingBox(min(a.x_min, b.x_min),
                                                      max(a.x_max, b.x_max),
                                                      min(a.y_min, b.y_min),
                                                      max(a.y_max, b.y_max))

BoundingBox() = BoundingBox(0.0, 0.0, 0.0, 0.0)

BoundingBox(point::Point) = BoundingBox(point.x, point.x + 1,
                                        point.y, point.y + 1)

BoundingBox(points::Points) = reduce(+, BoundingBox.(points), init=BoundingBox(points[1]))

Size(box::BoundingBox) = Size(box.x_max - box.x_min,
                              box.y_max - box.y_min)

function Base.ceil(::Type{T}, size::Size{B}) where {T, B}
    return Size(ceil(T, size.x),
                ceil(T, size.y))
end

end

max_stars_repo_name: Tiendil/pcg-julia

max_issues_repo_name: Tiendil/pcg-julia

max_forks_repo_name: Tiendil/pcg-julia


text:
# support for nanoseconds in dates
using Dates

struct DateTimeN
    d::Date
    t::Time
end

str = "2019-10-23T12:01:15.123456789"

parseDateTimeN(str)
parseDateTimeN( "2019-10-23T12:01:15.230")

parseDateTimeN(str)

max_stars_repo_name: xiaodaigh/DataConvenience

max_issues_repo_name: xiaodaigh/DataConvenience

max_forks_repo_name: xiaodaigh/DataConvenience


text:
using Yao, YaoBase, Zygote, LuxurySparse, YaoFlux
import Zygote: Context

t = ArrayReg(bit"0") + ArrayReg(bit"1")
normalize!(t)

# function rotmat(x)
#     I = IMatrix{2, ComplexF64}()
#     return I * cos(x/2) - im * sin(x / 2) * X
# end

# # circuit(x, y, z) = chain(Rx(x), Ry(y), Rz(z))

# # function loss(x, y, z)
# #     r = zero_state(1)
# #     # U = mat(chain(Rx(ps[1]), Rz(ps[2]), Rx(ps[3])))
# #     U = mat(chain(Rx(x), Ry(y), Rz(z)))
# #     return abs(statevec(t)' * U * statevec(r))
# # end

# # Zygote.refresh()
# # gradient(loss, 0.1, 0.1, 0.1)

# function loss(x, y, z)
#     r = zero_state(1)
#     apply!(r, Rx(0.1))
#     return abs(statevec(t)' * statevec(r))
# end

# gradient(loss, 0.1)

# (loss(0.1 + 1e-8) - loss(0.1)) / 1e-8

# Zygote.refresh()

# gradient(loss, 0.1)

function mapply!(r::ArrayReg, rb::RotationGate)
    # v0 = copy(r.state)
    # apply!(r, rb.block)
    # NOTE: we should not change register's memory address,
    # or batch operations may fail
    # r.state .= -im*sin(rb.theta/2)*r.state + cos(rb.theta/2)*r.state
    t = cos(rb.theta/2)*r.state
    copyto!(r.state, t)
    return r
end

Zygote.refresh()
_, back = Zygote._forward(x->mapply!(zero_state(1), Rx(x)), 0.1)

_, back = Zygote._forward(0.1) do x
    r = mapply!(zero_state(1), Rx(x))
    fidelity(t, r)[]
end

# _, back = Zygote._forward(r->abs(statevec(t)' * statevec(r)), zero_state(1))
# back(1)

# _, back = Zygote._forward(0.1) do x
#     r = zeros(2)
#     t = cos(x/2) * r
#     copyto!(r, t)
#     return r
# end

# st = rand(2)
# _, back = Zygote._forward(0.1) do x
#     abs(sum(x * st))
# end

# back(1)

# Zygote.refresh()
st = rand(2)
f(x) = abs(sum(x * st))

(f(0.1im+1e-8im) - f(0.1im)) / 1e-8im
f'(0.1im)


text:
include("testsystem1.jl")
srand(11)
Y, X = sample(20, 100, M)
kalmanfilter(Y, M)

include("testsystem.jl")
srand(11)
Y, X = sample(20, 100, M)
kalmanfilter(Y, M)


text:
"""
# Quickstart for PopGen
Documentation: https://pdimens.github.io/PopGen.jl/

Motivational(?) quote: $(motivational_quote())

\nA few things things you can do to get started:

## Load in data
- `read_from(filename; kwargs...)`
- `genepop(infile; kwargs...)`  or similar file-specific importer
- use available `gulfsharks()` or `nancycats()` datasets

## Explore PopData
- `populations(PopData)` to view population information
- `loci(PopData)` to view locus names
- `samples(PopData)` to view sample names
- `missing(PopData, by = ...)` to view missing information

## Manipulate PopData
- `populations!(PopData, ...)` to rename populations
- `locations!(PopData, ...)` to add geographical coordinates
- `exclude!(PopData; kwargs...)` to selectively remove data

## Analyses
- `richness(PopData)` to calculate allelic richness
- `allele_avg(PopData)` to calculate average # of alleles
- `summary(PopData)` to calculate F-statistics, heterozygosity, etc.
- `hwe_test(PopData)` to test for Hardy-Weinberg Equilibrium
"""
module PopGen

#   o O       o O       o O       o O       o O
# o | | O   o | | O   o | | O   o | | O   o | | O
# | | | | O | | | | Dependencies| | | | O | | | | O
# O | | o   O | | o   O | | o   O | | o   O | | o
#   O o       O o       O o       O o       O o

using CSV, Distributions, DataFrames
using FileIO, JLD2
using MultipleTesting, Random, StatsBase
using GeneticVariation: VCF, BCF, header

#=
using Convex,
      ECOS,
      LinearAlgebra,
      ProgressMeter,
=#

#   o O       o O       o O       o O       o O
# o | | O   o | | O   o | | O   o | | O   o | | O
# | | | | O | | | |Include Files| | | | O | | | | O
# O | | o   O | | o   O | | o   O | | o   O | | o
#   O o       O o       O o       O o       O o

# the types
include("Types.jl")
# file io
include("io/ioUtils.jl")
include("io/Delimited.jl")
include("io/Genepop.jl")
include("io/Read.jl")
include("io/VariantCall.jl")
# example data
include("io/Datasets.jl")
# utility functions
include("Utils.jl")
include("Permutations.jl")
# manipulation and exploration
include("DataExploration.jl")
include("Manipulate.jl")
# allele frequency and heterozygosity functions
include("AlleleFreq.jl")
include("Heterozygosity.jl")
# summary information
include("SummaryInfo.jl")
#Analyses
include("HardyWeinberg.jl")
#include(PairwiseRelatedness.jl)  # not yet ready
#include(PlotRecipes.jl)  # not yet ready

end # module PopGen


text:
@testset "  PD                                                       " begin
    @test_nowarn MetidaNCA.pdimport(pddata, :time, :obs, :subj; bl = 0, th = NaN)
end


text:
using Phonetics
using Test
using DelimitedFiles
using StatsBase
using DataFrames
using Statistics
using Random
using DynamicAxisWarping
using Distances

@testset "Vowel space and density" begin

  f = Formants("sample_formants.txt")
  v = VowelSpace(f)
  @test area(v, unstandardize=true) ≈ 1920.3803415872135
  @test vdi(v, unstandardize=true) ≈ 766.8716362331788
end

@testset "Random formant generation" begin

  dat = generateFormants(10^7, cats=["uw"], gender=["w"])
  datμ = mean(Array(dat[:, 1:2]), dims=1)
  @test isapprox(datμ, [459.67 1105.52], atol=1)

  dat = generateFormants(10^7, cats=["aa"], gender=["w"])
  datμ = mean(Array(dat[:, 1:2]), dims=1)
  @test isapprox(datμ, [916.36 1525.83], atol=1)

  dat = generateFormants(30, cats=["iy", "aa"], gender=["w", "m"])
  @test dat.vowel == repeat(["iy", "aa"], inner=60)
  @test dat.gender == repeat(["w", "m"], inner=30, outer=2)
end

@testset "Normalization" begin

  dat = generateFormants(50)
  # Use gender as speaker for tests to give two separate speakers

  lobCustom = copy(dat)
  rename!(lobCustom, :gender => :speaker)

  lobCustom[lobCustom.speaker .== "m", :f1] = zscore(lobCustom[lobCustom.speaker .== "m", :f1])
  lobCustom[lobCustom.speaker .== "w", :f1] = zscore(lobCustom[lobCustom.speaker .== "w", :f1])

  lobCustom[lobCustom.speaker .== "m", :f2] = zscore(lobCustom[lobCustom.speaker .== "m", :f2])
  lobCustom[lobCustom.speaker .== "w", :f2] = zscore(lobCustom[lobCustom.speaker .== "w", :f2])

  lobFunc = lobanov(dat.f1, dat.f2, dat.vowel, dat.gender)
  # check for approximate equality because sum and related functions
  # called on the view genereated by `groupby` produces a different value
  # than sum called on the original array (due to how the view is generated
  #  and how order matters when performing sequential floating point operations).
  @test lobCustom.f1 ≈ lobFunc.f1
  @test lobCustom.f2 ≈ lobFunc.f2
  @test lobCustom.vowel == lobFunc.vowel
  @test lobCustom.speaker == lobFunc.speaker

  nICustom = copy(dat)
  rename!(nICustom, :gender => :speaker)

  mf1 = nICustom[nICustom.speaker .== "m", :f1]
  nICustom[nICustom.speaker .== "m", :f1] = log.(mf1) .- mean(log.(mf1))
  wf1 = nICustom[nICustom.speaker .== "w", :f1]
  nICustom[nICustom.speaker .== "w", :f1] = log.(wf1) .- mean(log.(wf1))

  mf2 = nICustom[nICustom.speaker .== "m", :f2]
  nICustom[nICustom.speaker .== "m", :f2] = log.(mf2) .- mean(log.(mf2))
  wf2 = nICustom[nICustom.speaker .== "w", :f2]
  nICustom[nICustom.speaker .== "w", :f2] = log.(wf2) .- mean(log.(wf2))

  nIFunc = neareyI(dat.f1, dat.f2, dat.vowel, dat.gender)
  @test nICustom.f1 ≈ nIFunc.f1
  @test nICustom.f2 ≈ nIFunc.f2
  @test nICustom.vowel == nIFunc.vowel
  @test nICustom.speaker == nIFunc.speaker

  nECustom = copy(dat)

  rename!(nECustom, :gender => :speaker)

  mf1 = nECustom[nECustom.speaker .== "m", :f1]
  mfall = [mf1; nECustom[nECustom.speaker .== "m", :f2]]
  nECustom[nECustom.speaker .== "m", :f1] = log.(mf1) .- mean(log.(mfall))

  wf1 = nECustom[nECustom.speaker .== "w", :f1]
  wfall = [wf1; nECustom[nECustom.speaker .== "w", :f2]]
  nECustom[nECustom.speaker .== "w", :f1] = log.(wf1) .- mean(log.(wfall))

  mf2 = nECustom[nECustom.speaker .== "m", :f2]
  nECustom[nECustom.speaker .== "m", :f2] = log.(mf2) .- mean(log.(mfall))

  wf2 = nECustom[nECustom.speaker .== "w", :f2]
  nECustom[nECustom.speaker .== "w", :f2] = log.(wf2) .- mean(log.(wfall))

  nEFunc = neareyE(dat.f1, dat.f2, dat.vowel, dat.gender)
  @test nECustom.f1 ≈ nEFunc.f1
  @test nECustom.f2 ≈ nEFunc.f2
  @test nECustom.vowel == nEFunc.vowel
  @test nECustom.speaker == nEFunc.speaker
end

@testset "Acoustic distance" begin

  rng = MersenneTwister(9)
  x = rand(rng, 1000)
  y = rand(rng, 3000)
  z = rand(rng, 10000)
  a = [Sound(x, 8000), Sound(y, 8000), Sound(z, 8000)]
  amfcc = Phonetics.sound2mfcc.(a)
  @test acdist(a[1], a[2]) ≈ dtw(amfcc[1], amfcc[2])[1]
  @test distinctiveness(a[1], a[2:3]) ≈ mean([acdist(a[1], a[2]), acdist(a[1], a[3])])
  dtwr = maximum(size(s, 2) for s in amfcc)
  @test avgseq(a) == dba(amfcc, DTW(dtwr, SqEuclidean()), init_center=amfcc[2], show_progress=false)[1]
end

const SAMPLE_CORPUS = [
["K", "AE1", "T"], # cat
["K", "AA1", "B"], # cob
["B", "AE1", "T"], # bat
["T", "AE1", "T", "S"], # tats
["M", "AA1", "R", "K"], # mark
["K", "AE1", "B"], # cab
]

@testset "Uniqueness point" begin

  ## Unique at second character
  desired = DataFrame(Query=[["K", "AA1", "B"]], UPT=[2])
  @test upt(SAMPLE_CORPUS, [["K", "AA1", "B"]]; inCorpus=true) == desired

  ## Unique at first character
  desired = DataFrame(Query=[["HH", "AE1", "T"]], UPT=[1])
  @test upt(SAMPLE_CORPUS, [["HH", "AE1", "T"]]; inCorpus=false) == desired

  ## Unique at the last character
  desired = DataFrame(Query=[["K", "AE1", "D"]], UPT=[3])
  @test upt(SAMPLE_CORPUS, [["K", "AE1", "D"]]; inCorpus=false) == desired

  ## Not unique
  desired = DataFrame(Query=[["T", "AE1", "T"]], UPT=[4])
  @test upt(SAMPLE_CORPUS, [["T", "AE1", "T"]]; inCorpus=false) == desired
end

@testset "Phonological neighborhood density" begin

  # Showing progress

  ## Neighbors with bat and cab
  desired = DataFrame(Query=[["K", "AE1", "T"]], PND=[2])
  @test pnd(SAMPLE_CORPUS, [["K", "AE1", "T"]]; progress=true) == desired

  ### Neighbors with cat and bat
  desired = DataFrame(Query=[["HH", "AE1", "T"]], PND=[2])
  @test pnd(SAMPLE_CORPUS, [["HH", "AE1", "T"]]; progress=true) == desired

  ### No neighbors
  desired = DataFrame(Query=[["L", "IH1", "V"]], PND=[0])
  @test pnd(SAMPLE_CORPUS, [["L", "IH1", "V"]]; progress=true) == desired

  ## Not showing progress

  ### Neighbors with bat and cab
  desired = DataFrame(Query=[["K", "AE1", "T"]], PND=[2])
  @test pnd(SAMPLE_CORPUS, [["K", "AE1", "T"]]; progress=false) == desired

  ### Neighbors with cat and bat
  desired = DataFrame(Query=[["HH", "AE1", "T"]], PND=[2])
  @test pnd(SAMPLE_CORPUS, [["HH", "AE1", "T"]]; progress=false) == desired

  ### No neighbors
  desired = DataFrame(Query=[["L", "IH1", "V"]], PND=[0])
  @test pnd(SAMPLE_CORPUS, [["L", "IH1", "V"]]; progress=false) == desired
end

@testset "Phonotactic probability" begin

  #=
  **Co-occurrence**

  20 total monophones
  K:      4
  AE1:    4
  T:      4
  S:      1
  M:      1
  AA1:    2
  B:      3
  R:      1

  ["K", "AE1", "T"], get counts
  [4, 4, 4], divide by total monophones
  [0.2, 0.2, 0.2], take product
  0.2^3 = 0.008
  =#
  freq = [1,1,1,1,1,1]
  p = prod([4,4,4] / 20)
  desired = DataFrame(Query=[["K", "AE1", "T"]], Probability=[p])
  prb = phnprb(SAMPLE_CORPUS, freq, [["K", "AE1", "T"]])
  @test prb == desired

  #=

  **Co-occurrence diphones**

  26 total diphones
  . K:    3
  K  AE1: 2
  AE1 T:  3
  T .:    2
  K AA1:  1
  AA1 B:  1
  B .:    2
  . B:    1
  B AE1:  1
  . T:    1
  T AE1:  1
  T S:    1
  S .:    1
  . M:    1
  M AA1:  1
  AA1 R:  1
  R K:    1
  K .:    1
  AE1 B:  1

  [(., K), (K, AE1), (AE1, T), (T .)]
  [3, 2, 3, 2]
  [3/26, 2/26, 3/26, 2/26]
  36/(26^4) = 7.9e-5
  =#
  p = prod([3,2,3,2]/26)
  desired = DataFrame(Query=[["K", "AE1", "T"]], Probability=[p])
  prb = phnprb(SAMPLE_CORPUS, freq, [["K", "AE1", "T"]]; nchar=2)
  @test prb == desired

  #=
  **Co-occurrence diphones, no padding**

  14 total diphones
  K  AE1: 2
  AE1 T:  3
  K AA1:  1
  AA1 B:  1
  B AE1:  1
  T AE1:  1
  T S:    1
  M AA1:  1
  AA1 R:  1
  R K:    1
  AE1 B:  1

  [(K, AE1), (AE1, T)]
  [2, 3]
  [3/14, 2/14] # take product
  6/(14^2) = 7.9e-5
  =#
  p = prod([3,2]/14)
  desired = DataFrame(Query=[["K", "AE1", "T"]], Probability=[p])
  prb = phnprb(SAMPLE_CORPUS, freq, [["K", "AE1", "T"]]; nchar=2, pad=false)
  @test prb == desired

  #=
  **Positional*

  20 total monophones
  K 1:    3
  AE1 2:  4
  AA1 2:  2
  B 3:    2
  B 1:    1
  T 3:    3
  T 1:    1
  S 4:    1
  M 1:    1
  R 3:    1
  K 4:    1

  ["K", "AE1", "T"], get counts
  [3, 4, 3], divide by total monophones
  [3/20, 4/20, 3/20], take product
  36/8000 = 0.0045
  =#
  p = prod([3,4,3]/20)
  desired = DataFrame(Query=[["K", "AE1", "T"]], Probability=[p])
  prb = phnprb(SAMPLE_CORPUS, freq, [["K", "AE1", "T"]]; positional=true)
  @test prb == desired

  # Phone not in corpus
  p = 0
  desired = DataFrame(Query=[["K", "AE1", "F"]], Probability=[p])
  @test phnprb(SAMPLE_CORPUS, freq, [["K", "AE1", "F"]]) == desired

  # With non-1 frequencies; doubling each frequency should produce same output
  freq = [2,2,2,2,2,2]
  p = prod([8,8,8] / 40)
  desired = DataFrame(Query=[["K", "AE1", "T"]], Probability=[p])
  prb = phnprb(SAMPLE_CORPUS, freq, [["K", "AE1", "T"]])
  @test prb == desired

end

@testset "Vantage point tree" begin

  tree = TextVPTree(SAMPLE_CORPUS, lev)

  # Test nneighbors
  res = collect(keys(nneighbors(tree, ["K", "AE1", "T"], 2)))
  res = sort(join.(res, " "))
  @test res == ["B AE1 T", "K AE1 B"]

  ## Test on point not inside of the tree
  res = collect(keys(nneighbors(tree, ["HH", "AE1", "T"], 2)))
  res = sort(join.(res, " "))
  @test res == ["B AE1 T", "K AE1 T"]

  # Test radius search
  res = radiusSearch(tree, ["K", "AE1", "T"], 1)
  res = sort(join.(res, " "))
  @test res == ["B AE1 T", "K AE1 B"]

  ## Test on point not inside of the tree
  res = radiusSearch(tree, ["HH", "AE1", "T"], 1)
  res = sort(join.(res, " "))
  @test res == ["B AE1 T", "K AE1 T"]
end


text:
# !! auto-generated code, do not edit !!
module trajectory_msgs
using RobotOSData.Messages
using RobotOSData.StdMsgs
using RobotOSData.CommonMsgs
struct MultiDOFJointTrajectoryPoint <: Readable
    transforms::Vector{geometry_msgs.Transform}
    velocities::Vector{geometry_msgs.Twist}
    accelerations::Vector{geometry_msgs.Twist}
    time_from_start::ROSDuration
end
struct MultiDOFJointTrajectory <: Readable
    header::Header
    joint_names::Vector{String}
    points::Vector{MultiDOFJointTrajectoryPoint}
end
struct JointTrajectoryPoint <: Readable
    positions::Vector{Float64}
    velocities::Vector{Float64}
    accelerations::Vector{Float64}
    effort::Vector{Float64}
    time_from_start::ROSDuration
end
struct JointTrajectory <: Readable
    header::Header
    joint_names::Vector{String}
    points::Vector{JointTrajectoryPoint}
end
end


text:
using Mimi

include("../../src/compute_scc.jl")

function setorup_param!(m::Model, param::Symbol, value)
    try
        set_param!(m, param, value)
    catch e
        update_param!(m, param, value)
    end
end

"""
Applies undiscounting factor to get the SCC, discounted to the emissions year instead of the base year.
"""
function undiscount_scc(m::Model, year::Int)
    df = m[:EquityWeighting, :df_utilitydiscountfactor_ann]
    consfocus0 = m[:GDP, :cons_percap_consumption_0][1]
    consfocus = m[:GDP, :cons_percap_consumption_ann][:, 1]
    emuc = m[:EquityWeighting, :emuc_utilityconvexity]
    sccii = getpageindexfromyear(year)

    return df[sccii] * ((consfocus[sccii] / consfocus0)^-emuc)
end
"""
compute_scc(m::Model = get_model(); year::Union{Int, Nothing} = nothing, eta::Union{Float64, Nothing} = nothing, prtp::Union{Float64, Nothing} = nothing)

Computes the social cost of CO2 for an emissions pulse in `year` for the provided Mimi-PAGE model.
If no model is provided, the default model from main_model.get_model() is used.
Discounting scheme can be specified by the `eta` and `prtp` parameters, which will update the values of emuc_utilitiyconvexity and ptp_timepreference in the model.
If no values are provided, the discount factors will be computed using the default PAGE values of emuc_utilitiyconvexity=1.1666666667 and ptp_timepreference=1.0333333333.
"""
function compute_scc(m::Model=get_model(); year::Union{Int,Nothing}=nothing, eta::Union{Float64,Nothing}=nothing, prtp::Union{Float64,Nothing}=nothing, pulse_size=75000.)
    year === nothing ? error("Must specify an emission year. Try `compute_scc(m, year=2020)`.") : nothing
    !(year in page_years) ? error("Cannot compute the scc for year $year, year must be within the model's time index $page_years.") : nothing

    eta == nothing ? nothing : setorup_param!(m, :emuc_utilityconvexity, eta)
    prtp == nothing ? nothing : setorup_param!(m, :ptp_timepreference, prtp * 100.)

    mm = get_marginal_model(m, year=year, pulse_size=pulse_size)   # Returns a marginal model that has already been run
    scc = mm[:EquityWeighting, :td_totaldiscountedimpacts_ann] / undiscount_scc(mm.base, year)

    return scc
end

"""
compute_scc_mm(m::Model = get_model(); year::Union{Int, Nothing} = nothing, eta::Union{Float64, Nothing} = nothing, prtp::Union{Float64, Nothing} = nothing)

Returns a NamedTuple (scc=scc, mm=mm) of the social cost of carbon and the MarginalModel used to compute it.
Computes the social cost of CO2 for an emissions pulse in `year` for the provided Mimi-PAGE model.
If no model is provided, the default model from main_model.get_model() is used.
Discounting scheme can be specified by the `eta` and `prtp` parameters, which will update the values of emuc_utilitiyconvexity and ptp_timepreference in the model.
If no values are provided, the discount factors will be computed using the default PAGE values of emuc_utilitiyconvexity=1.1666666667 and ptp_timepreference=1.0333333333.
"""
function compute_scc_mm(m::Model=get_model(); year::Union{Int,Nothing}=nothing, eta::Union{Float64,Nothing}=nothing, prtp::Union{Float64,Nothing}=nothing, pulse_size=75000.)
    year === nothing ? error("Must specify an emission year. Try `compute_scc(m, year=2020)`.") : nothing
    !(year in page_years) ? error("Cannot compute the scc for year $year, year must be within the model's time index $page_years.") : nothing

    eta == nothing ? nothing : setorup_param!(m, :emuc_utilityconvexity, eta)
    prtp == nothing ? nothing : setorup_param!(m, :ptp_timepreference, prtp * 100.)

    mm = get_marginal_model(m, year=year, pulse_size=pulse_size)   # Returns a marginal model that has already been run
    scc = mm[:EquityWeighting, :td_totaldiscountedimpacts_ann] / undiscount_scc(mm.base, year)
    scc_disaggregated = mm[:EquityWeighting, :addt_equityweightedimpact_discountedaggregated_ann] / undiscount_scc(mm.base, year)

    return (scc = scc, scc_disaggregated = scc_disaggregated, mm = mm)
end

"""
get_marginal_model(m::Model = get_model(); year::Union{Int, Nothing} = nothing)
Returns a Mimi MarginalModel where the provided m is the base model, and the marginal model has additional emissions of CO2 in year `year`.
If no Model m is provided, the default model from main_model.get_model() is used as the base model.
Note that the returned MarginalModel has already been run.
"""
function get_marginal_model(m::Model=get_model(); year::Union{Int,Nothing}=nothing, pulse_size=75000.)
    year === nothing ? error("Must specify an emission year. Try `get_marginal_model(m, year=2020)`.") : nothing
    !(year in page_years) ? error("Cannot add marginal emissions in $year, year must be within the model's time index $page_years.") : nothing

    mm = create_marginal_model(m, pulse_size)

    add_comp!(mm.modified, ExtraEmissions, :extra_emissions; after=:co2emissions)
    connect_param!(mm.modified, :extra_emissions => :e_globalCO2emissions, :co2emissions => :e_globalCO2emissions)
    set_param!(mm.modified, :extra_emissions, :pulse_size, pulse_size)
    set_param!(mm.modified, :extra_emissions, :pulse_year, year)

    connect_param!(mm.modified, :CO2Cycle => :e_globalCO2emissions, :extra_emissions => :e_globalCO2emissions_adjusted)

    run(mm)
    return mm
end

function compute_scc_mcs(m::Model, samplesize::Int; year::Union{Int,Nothing}=nothing, eta::Union{Float64,Nothing}=nothing, prtp::Union{Float64,Nothing}=nothing, pulse_size=75000.)
    # Setup of location of final results
    scc_results = zeros(samplesize)

    function mc_scc_calculation(sim_inst::SimulationInstance, trialnum::Int, ntimesteps::Int, ignore::Nothing)
        marginal = sim_inst.models[1]
        marg_damages = marginal[:EquityWeighting, :td_totaldiscountedimpacts_ann] / undiscount_scc(mm.base, year)
        scc_results[trialnum] = marg_damages
    end

    # get simulation
    mcs = getsim()

    # Setup models
    eta == nothing ? nothing : setorup_param!(m, :emuc_utilityconvexity, eta)
    prtp == nothing ? nothing : setorup_param!(m, :ptp_timepreference, prtp * 100.)

    mm = get_marginal_model(m, year=year, pulse_size=pulse_size)   # Returns a marginal model that has already been run

    # Run
    res = run(mcs, mm, samplesize; post_trial_func=mc_scc_calculation)

    scc_results
end


text:
@testset "Sentiment" begin
m = SentimentAnalyzer()

d = StringDocument("a very nice thing that everyone likes")

m(d) > 0.5

d = StringDocument("a horrible thing that everyone hates")

m(d) < 0.5
end


text:
using PDFIO, DataFrames, Dates, CSV

# A routine to extract data from the IOM Government pdf snapshots

#folder - folder where the snapshots are stored
folder = ""

#outFold - folder where the data output csvs should be written
outFold = ""

#read in the files in the folder
lst = readdir(folder)

#select only valid files for processing
prolst = lst[lst .>= "2021-07-29.pdf"]

#initialise variables
io = IOBuffer()
testingOutput = Vector[]
frontOutput = Vector[]

#loop through each file for processing
for pdname in prolst

    #open pdf file for reading
    p = pdDocOpen(folder*pdname)

        #routine to extract relevant data from page 8
        pdPageExtractText(io,pdDocGetPage(p, 8))
        data = String(take!(io))

        date =
            try
                Date(match(r"Default Date\s+(\d\d/\d\d/\d\d\d\d)",data).captures[1],"dd/mm/yyyy")
            catch e
                0
            end
        tests =
            try
                parse(Int64,replace(match(r"((?:\d|,)+)\s+TESTS",data).captures[1],","=>""))
            catch e
                -1
            end
        concluded =
            try
                parse(Int64,replace(match(r"TESTS(?:.|\R)+?((?:\d|,)+)\s+TESTS",data).captures[1],","=>""))
            catch e
                -1
            end
        rate =
            try
                parse(Float64,replace(match(r"((?:\d|\.)+)%\s+RATE",data).captures[1],"%"=>""))
            catch e
                -1
            end
        awaitresult =
            try
                parse(Int64,replace(match(r"((?:\d|,)+)\s+AWAITING RESULT",data).captures[1],","=>""))
            catch e
                -1
            end
        bookedtest =
            try
                parse(Int64,replace(match(r"((?:\d|,)+)\s+BOOKED TESTS",data).captures[1],","=>""))
            catch e
                -1
            end
        push!(testingOutput, [date, tests, concluded, rate, awaitresult, bookedtest])

        #routine to extract relevant data from page 1
        pdPageExtractText(io,pdDocGetPage(p, 1))
        data = String(take!(io))

        date =
        try
            Date(match(r"Data Last Refreshed\s+(\d\d/\d\d/\d\d\d\d)",data).captures[1],"dd/mm/yyyy")
        catch e
            0
        end

        topush = String[]
        for m in collect(eachmatch(r"(?:\w|\s)((?:\d|,)+)(?:\w|\s)",data))
            push!(topush, m.captures[1])
        end
        for n in length(topush):20
            push!(topush,"-1")
        end
        topush

        push!(frontOutput, [date, topush...])

    pdDocClose(p)

end

CSV.write(outFold*"front.csv",DataFrame(frontOutput))
CSV.write(outFold*"testing.csv",DataFrame(testingOutput))

max_stars_repo_name: dannywinrow/iom-govt-data

max_issues_repo_name: dannywinrow/iom-govt-data

max_forks_repo_name: dannywinrow/iom-govt-data


text:
module DataFrames

using Statistics, Printf, REPL
using Reexport, SortingAlgorithms, Compat, Unicode, PooledArrays
@reexport using Missings, InvertedIndices
using Base.Sort, Base.Order, Base.Iterators, Base.Threads
using TableTraits, IteratorInterfaceExtensions
import LinearAlgebra: norm
using Markdown
using PrettyTables

import DataAPI,
       DataAPI.All,
       DataAPI.Between,
       DataAPI.Cols,
       DataAPI.describe,
       Tables,
       Tables.columnindex,
       Future.copy!

export AbstractDataFrame,
       All,
       AsTable,
       Between,
       ByRow,
       Cols,
       DataFrame,
       DataFrameRow,
       GroupedDataFrame,
       SubDataFrame,
       Tables,
       allowmissing!,
       antijoin,
       by,
       columnindex,
       combine,
       completecases,
       crossjoin,
       describe,
       disallowmissing!,
       dropmissing!,
       dropmissing,
       flatten,
       groupby,
       groupindices,
       groupcols,
       innerjoin,
       insertcols!,
       leftjoin,
       mapcols,
       mapcols!,
       ncol,
       nonunique,
       nrow,
       order,
       outerjoin,
       PrettyTables,
       rename!,
       rename,
       repeat!,
       rightjoin,
       rownumber,
       select!,
       select,
       semijoin,
       stack,
       subset,
       subset!,
       transform,
       transform!,
       unique!,
       unstack,
       valuecols

# TODO: remove these exports in year 2021
export by, aggregate

if VERSION >= v"1.1.0-DEV.792"
    import Base.eachcol, Base.eachrow
else
    import Compat.eachcol, Compat.eachrow
    export eachcol, eachrow
end

if VERSION < v"1.2"
    export hasproperty
end

if isdefined(Base, :only)  # Introduced in 1.4.0
    import Base.only
else
    import Compat.only
    export only
end

if VERSION >= v"1.3"
    using Base.Threads: @spawn
else
    # This is the definition of @async in Base
    macro spawn(expr)
        thunk = esc(:(()->($expr)))
        var = esc(Base.sync_varname)
        quote
            local task = Task($thunk)
            if $(Expr(:isdefined, var))
                push!($var, task)
            end
            schedule(task)
        end
    end
end

include("other/utils.jl")
include("other/index.jl")

include("abstractdataframe/abstractdataframe.jl")
include("dataframe/dataframe.jl")
include("subdataframe/subdataframe.jl")
include("dataframerow/dataframerow.jl")
include("groupeddataframe/groupeddataframe.jl")
include("dataframerow/utils.jl")

include("other/broadcasting.jl")

include("abstractdataframe/selection.jl")
include("abstractdataframe/subset.jl")
include("abstractdataframe/iteration.jl")
include("abstractdataframe/reshape.jl")

include("join/composer.jl")
include("join/core.jl")

include("groupeddataframe/splitapplycombine.jl")
include("groupeddataframe/callprocessing.jl")
include("groupeddataframe/fastaggregates.jl")
include("groupeddataframe/complextransforms.jl")

include("abstractdataframe/prettytables.jl")
include("abstractdataframe/show.jl")
include("groupeddataframe/show.jl")
include("dataframerow/show.jl")
include("abstractdataframe/io.jl")

include("abstractdataframe/sort.jl")
include("dataframe/sort.jl")

include("deprecated.jl")

include("other/tables.jl")

include("other/precompile.jl")
precompile()

end # module DataFrames


text:
"""
    generate!(
        registry_packages::Vector{String},
        registry_name::String,
        registry_repo::String,
        [registry_description::String = "", writedir::String = pwd()];
        force::bool = true,
        METADATA_path::String = Pkg.Pkg2.dir("METADATA"),
        registry_dirs::Vector{String} = [joinpath(homedir(), ".julia", "registries", "General")],
        stdlib_dir::String = ""
     )

Generate a registry named `registry_name` with repository `registry_repo`
that contains the packages listed in `registry_packages`. These packages must be
valid METADATA packages in the METADATA directory at `METADATA_path`, but do not
need to be in the official Julia METADATA repository.

The packages in the new registry can have dependencies on packages in other
registries. By default the "General" registry is available to resolve dependencies,
any set of registries can be specified by setting `registry_dirs` to contain the
path for each registry.

Dependencies on the standard library can either be resolved by setting `stdlib_dir`
to be the path of the standard libraries in a git repository (i.e. in the julia
repository), or if `stdlib_dir` is empty then a cached set of UUIDs will be used.
"""
function generate! end
function generate!(
    pkgs::AbstractDict{<:AbstractString, Package},
    other_registries::AbstractDict{<:AbstractString, Package},
    args...; kwargs...
)
    check_packages(pkgs, other_registries)
    _generate!(pkgs, other_registries, args...; kwargs...)
end
function generate!(
    pkgs::AbstractDict{<:AbstractString, Package},
    other_registries::AbstractVector{<:AbstractDict{<:AbstractString, Package}},
    args...; kwargs...
)
    generate!(pkgs, merge(other_registries...), args...; kwargs...)
end
function generate!(
    pkgs::AbstractDict{<:AbstractString, Package},
    registry_name::AbstractString,
    args...;
    registry_dirs::AbstractVector{<:AbstractString} = [
        joinpath(homedir(), ".julia", "registries", "General")
    ],
    kwargs...
)
    generate!(
        pkgs, parse_registry.(registry_dirs), registry_name, args...; kwargs...
    )
end
function generate!(
    package_strings::AbstractVector{<:AbstractString},
    args...;
    METADATA_path::AbstractString = Pkg.Pkg2.dir("METADATA"),
    kwargs...
)
    isempty(package_strings) && throw(ArgumentError("Must specify packages for registry"))
    pkgs = load_packages_METADATA(METADATA_path, package_strings)
    generate!(pkgs, args...; kwargs...)
end

function _generate!(
    pkgs::AbstractDict{<:AbstractString, Package},
    other_registries::AbstractDict{<:AbstractString, Package},
    registry_name::AbstractString,
    registry_repo::AbstractString,
    registry_description::AbstractString = "",
    writedir::AbstractString = pwd();
    stdlib_dir::AbstractString = "",
    force::Bool = false,
)
    isdir(writedir) || throw(ArgumentError("$writedir does not exist"))
    prefix = joinpath(writedir, registry_name)
    if isdir(prefix)
        if ! isempty(readdir(prefix)) && ! force
            error("$prefix exists and is not empty, use 'force' to overwrite")
        end
    elseif ispath(prefix)
        if force
            mv(prefix, "$prefix.bak")
            mkdir(prefix)
        else
            error("$prefix is a file, use 'force' to overwrite")
        end
    else
        mkdir(prefix)
    end

    write_toml(prefix, "Registry") do io
        uuid = string(uuid5(uuid_registry, registry_repo))
        println(io, "name = ", repr(registry_name))
        println(io, "uuid = ", repr(uuid))
        println(io, "repo = ", repr(registry_repo))
        println(io, "\ndescription = \"\"\"")
        print(io, registry_description)
        println(io, "\"\"\"")
        println(io, "\n[packages]")
        for (pkg, p) in sort!(collect(pkgs), by=(p->p.uuid.value)∘last)
            bucket = string(uppercase(first(pkg)))
            path = joinpath(bucket, pkg)
            println(
                io,
                p.uuid, " = { name = ", repr(pkg), ", path = ", repr(path), " }"
            )
        end
    end

    buckets = Dict()
    for (pkg, p) in pkgs
        bucket = string(uppercase(first(pkg)))
        push!(get!(buckets, bucket, []), (pkg, p))
    end

    if isempty(stdlib_dir)
        for pkg in STDLIBS
            other_registries[pkg] = Package(
                UUID(STDLIB_UUID_STRS[pkg]),
                "https://github.com/JuliaLang/julia.git",
                Dict(
                    VersionNumber(0,7,0) => Version(
                        nothing, Dict{String,Require}()
                    )
                )
            )
        end
    else
        stdlib_uuids, stdlib_trees, stdlib_deps = stdlib_regdata(stdlib_dir)

        for pkg in STDLIBS
            tree = stdlib_trees[pkg]
            deps = Dict(
                dep => Require(VersionInterval()) for dep in stdlib_deps[pkg]
            )
            other_registries[pkg] = Package(
                UUID(stdlib_uuids[pkg]),
                "https://github.com/JuliaLang/julia.git",
                Dict(
                    VersionNumber(0,7,0,("DEV",),("r"*tree[1:8],)) =>
                    Version(tree, deps)
                ),
            )
        end
    end

    trees, stdlibs = gitmeta(pkgs)

    for (pkg, p) in pkgs
        uuid = string(p.uuid)
        haskey(stdlibs, uuid) || continue
        for (ver, v) in p.versions
            n = get(stdlibs[uuid], v.sha1, 0)
            n == 0 && continue
            for lib in STDLIBS
                if n & 1 != 0
                    v.requires[lib] = Require(VersionInterval())
                end
                n >>>= 1
            end
        end
    end

    for (bucket, b_pkgs) in buckets, (pkg, p) in b_pkgs
        haskey(stdlibs, pkg) && continue
        url = p.url
        uuid = string(p.uuid)
        startswith(url, "git://github.com") && (url = "https"*url[4:end])

        # Package.toml
        write_toml(prefix, bucket, pkg, "Package") do io
            println(io, "name = ", repr(pkg))
            println(io, "uuid = ", repr(uuid))
            println(io, "repo = ", repr(url))
        end

        # Versions.toml
        write_toml(prefix, bucket, pkg, "Versions") do io
            for (i, (ver, v)) in enumerate(sort!(collect(p.versions), by=first))
                i > 1 && println(io)
                println(io, "[", toml_key(string(ver)), "]")
                println(io, "git-tree-sha1 = ", repr(trees[uuid][v.sha1]))
            end
        end
        versions = sort!(collect(keys(p.versions)))

        function write_versions_data(f::Function, name::String; lt::Function=isless)
            data = Dict{VersionNumber,Dict{String,String}}()
            for (ver, v) in p.versions, (dep, d) in v.requires
                val = f(dep, d)
                val == nothing && continue
                haskey(data, ver) || (data[ver] = Dict{String,String}())
                # BinDeps injects a dependency on Libdl
                if name == "Deps" && dep == "BinDeps"
                    data[ver]["Libdl"] = "\"8f399da3-3557-5675-b5ff-fb832c97cbdb\""
                end
                data[ver][dep] = val
            end
            compressed = compress_versions_data(data, versions)
            !isempty(compressed) && write_toml(prefix, bucket, pkg, name) do io
                vers = unique(getindex.(compressed, 1))
                keys = sort!(unique(getindex.(compressed, 2)), lt=lt)
                what = (vers, keys)
                ord = (1, 2)
                for (i, x) in enumerate(what[ord[1]])
                    i > 1 && println(io)
                    println(io, "[", toml_key(x), "]")
                    for y in what[ord[2]]
                        for t in compressed
                            t[ord[1]] == x && t[ord[2]] == y || continue
                            println(io, toml_key(y), " = ", t[3])
                        end
                    end
                end
            end
        end

        # Deps.toml
        write_versions_data("Deps") do dep, d
            if dep == "julia"
                res = nothing
            else
                registry = haskey(pkgs, dep) ? pkgs : other_registries
                res = repr(string(registry[dep].uuid))
            end
            res
        end

        # Compat.toml
        write_versions_data("Compat", lt=packagelt) do dep, d
            if dep in STDLIBS
                res = nothing
            else
                registry = haskey(pkgs, dep) ? pkgs : other_registries
                res = versions_repr(compress_versions(
                    d.versions, collect(keys(registry[dep].versions))
                ))
            end
            res
        end
    end
end


text:
# AdaGrad + Analytic Gradient

invlogit(z::Real) = 1 / (1 + exp(-z))

type ApproxL2Logit <: StreamStat
    λ::Float64
    α::Float64
    β₀::Float64
    β::Vector{Float64}
    sum_gr₀_sq::Float64
    sum_gr_sq::Vector{Float64}
    n::Int
end

function ApproxL2Logit(p::Integer, λ::Real, α::Real = 0.1)
    return ApproxL2Logit(
        float64(λ),
        float64(α),
        0.0,
        fill(0.0, p),
        0.0,
        fill(0.0, p),
        0,
    )
end

function update!(stat::ApproxL2Logit, x::Vector{Float64}, y::Real)
    stat.n += 1

    y_pred = invlogit(stat.β₀ + dot(stat.β, x))

    ε = y - y_pred

    # Intercept
    gr₀ = ε - stat.λ * stat.β₀
    stat.sum_gr₀_sq += gr₀^2
    α₀ = stat.α / sqrt(stat.sum_gr₀_sq)
    stat.β₀ += α₀ * gr₀

    # Non-constant predictors
    for i in 1:length(x)
        grᵢ = ε * x[i] - stat.λ * stat.β[i]
        stat.sum_gr_sq[i] += grᵢ^2
        αᵢ = stat.α / sqrt(stat.sum_gr_sq[i])
        stat.β[i] += αᵢ * grᵢ
    end

    return
end

state(stat::ApproxL2Logit) = vcat(stat.β₀, stat.β)

nobs(stat::ApproxL2Logit) = stat.n

function Base.show(io::IO, stat::ApproxL2Logit)
    p = length(state(stat))
    n = nobs(stat)
    @printf(io, "Online L2-Regularized Logistic Regression\n")
    @printf(io, " * λ: %f\n", stat.λ)
    @printf(io, " * P: %d\n", p)
    @printf(io, " * N: %d\n", n)
    return
end

function Base.empty!(stat::ApproxL2Logit)
    stat.β₀ = 0
    fill!(stat.β, 0)
    stat.sum_gr₀_sq = 0
    fill!(stat.sum_gr_sq, 0)
    stat.n = 0
    return
end

# Base.merge


text:
using PyCall
using PrettyTables

#pydir = joinpath(dirname(pathof(FermiCG)), "python")
#pushfirst!(PyVector(pyimport("sys")."path"), pydir)
#ENV["PYTHON"] = Sys.which("python")
#print(ENV)

"""
	pyscf_do_scf(molecule::Molecule, conv_tol=1e-10)

Use PySCF to compute Hartree-Fock for a given molecule and basis set
and return a PYSCF mean field object
"""
function pyscf_do_scf(molecule::Molecule; conv_tol=1e-10, verbose=0)
    pyscf = pyimport("pyscf")
    pyscf.lib.num_threads(1)
    pymol = make_pyscf_mole(molecule)

    println(pymol.basis)
    #pymol.max_memory = 1000 # MB
    #pymol.symmetry = true
    mf = pyscf.scf.RHF(pymol).run(conv_tol=conv_tol, verbose=verbose)
    enu = mf.energy_nuc()
    println("MO Energies")
    display(mf.mo_energy)
    # print(np.linalg.eig(mf.get_fock())[0])

    # if pymol.symmetry == True:
    # 	from pyscf import symm
    # 	mo = symm.symmetrize_orb(mol, mf.mo_coeff)
    # 	osym = symm.label_orb_symm(mol, mol.irrep_name, mol.symm_orb, mo)
    # 	#symm.addons.symmetrize_space(mol, mo, s=None, check=True, tol=1e-07)
    # 	for i in range(len(osym)):
    # 		print("%4d %8s %16.8f"%(i+1,osym[i],mf.mo_energy[i]))

    return mf
end

"""
	make_pyscf_mole(molecule::Molecule)

Create a `pyscf.gto.Mole()` object
"""
function make_pyscf_mole(molecule::Molecule)
    pyscf = pyimport("pyscf")
    pyscf.lib.num_threads(1)
    pymol = pyscf.gto.Mole()
    pymol.basis = molecule.basis
    geomstr = ""
    for i in molecule.atoms
        geomstr = geomstr * string(i.symbol,", ", join(map(string, i.xyz), ", "),"\n")
    end
    pymol.atom = geomstr
    pymol.charge = molecule.charge
    pymol.spin = molecule.multiplicity-1
    pymol.build()
    return pymol
end

"""
	pyscf_write_molden(molecule::Molecule, C; filename="orbitals.molden")

# Arguments
- `molecule::Molecule`: Molecule object
- `C`: MO Coefficients
- `filename`: Filename to write to

Write MO coeffs `C` to a molden file for visualizing
"""
function pyscf_write_molden(molecule::Molecule, C; filename="orbitals.molden")
    pyscf = pyimport("pyscf")
    pyscf.lib.num_threads(1)
    molden = pyimport("pyscf.tools.molden")
    pymol = make_pyscf_mole(molecule)
    molden.from_mo(pymol, filename, C)
    return 1
end

"""
	pyscf_write_molden(mf; filename="orbitals.molden")

# Arguments
- `mf`: PySCF mean field object
- `filename`: Filename to write to

Write MO coeffs `C` to a molden file for visualizing
"""
function pyscf_write_molden(mf; filename="orbitals.molden")
    pyscf = pyimport("pyscf")
    tools = pyimport("pyscf.tools")
    molden = pyimport("tools.molden")
    molden.from_mo(mf.mol, filename, mf.mo_coeff)
    return 1
end

"""
	pyscf_build_1e(mol::Molecule)

build 1 electron integrals in AO basis
# Arguments
- `mol::Molecule`

returns a 2D matrix
"""
function pyscf_build_1e(mol::Molecule)

    pyscf = pyimport("pyscf")
    pyscf.lib.num_threads(1)
    #
    # get pyscf molecule type
    pymol = FermiCG.make_pyscf_mole(mol)
    return pyscf.scf.hf.get_hcore(pymol)
end

"""
	pyscf_build_eri(mol::Molecule, c1::Matrix, c2::Matrix, c3::Matrix, c4::Matrix)

build 2 electron integrals between different orbital spaces, (c1c2|c3c4)
# Arguments
- `mol::Molecule`
- `c1`: active space orbital MO coeffs for index 1
- `c2`: active space orbital MO coeffs for index 2
- `c3`: active space orbital MO coeffs for index 3
- `c4`: active space orbital MO coeffs for index 4

returns a 4D tensor
"""
function pyscf_build_eri(mol::Molecule, c1::Matrix, c2::Matrix, c3::Matrix, c4::Matrix)

    pyscf = pyimport("pyscf")
    pyscf.lib.num_threads(1)
    #
    # get pyscf molecule type
    pymol = FermiCG.make_pyscf_mole(mol)
    eri = pymol.intor("int2e")
    h2 = pyscf.ao2mo.incore.general(eri, (c1,c2,c3,c4))

    nm1 = size(c1)[2]
    nm2 = size(c2)[2]
    nm3 = size(c3)[2]
    nm4 = size(c4)[2]
    h2 = reshape(h2, (nm1, nm2, nm3, nm4))

    return h2
end

"""
	pyscf_get_jk(mol, density)

Build exchange matrix in AO basis

# Arguments
- `mol`: PySCF Molecule object
- `density`: 1rdm density matrix for

"""
function pyscf_get_jk(mol::Molecule, density)
	pyscf = pyimport("pyscf")
    pyscf.lib.num_threads(1)

    #
    # get pyscf molecule type
    pymol = FermiCG.make_pyscf_mole(mol)

	h0 = pyscf.gto.mole.energy_nuc(pymol)
	h  = pyscf.scf.hf.get_hcore(pymol)
    j, k = pyscf.scf.hf.get_jk(pymol, density, hermi=1)
    return h, j, k
end

"""
	pyscf_build_ints(mol, c_act, d1_embed)

build 1 and 2 electron integrals using a pyscf SCF object
# Arguments
- `mol`: PySCF Molecule object
- `c_act`: active space orbital MO coeffs
- `d1_embed`: 1rdm density matrix for the frozen part in the AO basis (e.g, doccs or frozen clusters)

returns an `InCoreInts` type
"""
function pyscf_build_ints(mol::Molecule, c_act, d1_embed)

	pyscf = pyimport("pyscf")
    pyscf.lib.num_threads(1)

	nact = size(c_act)[2]
	#mycas = pyscf.mcscf.CASSCF(mf, length(active), 0)
    #
    # get pyscf molecule type
    pymol = FermiCG.make_pyscf_mole(mol)

	h0 = pyscf.gto.mole.energy_nuc(pymol)
	h  = pyscf.scf.hf.get_hcore(pymol)
    j, k = pyscf.scf.hf.get_jk(pymol, d1_embed, hermi=1)

    # get core energy
    #h0 = tr(d1_embed * ( h + .5*j - .5*k))
    #mf = pyscf.scf.RHF(pymol)
    #h0 = pyscf.scf.hf.energy_elec(mf, d1_embed)[1]
    h0 += tr(d1_embed * ( h + .5*j - .25*k))

    # now rotate to MO basis
    h = c_act' * h * c_act
	j = c_act' * j * c_act;
	k = c_act' * k * c_act;
	h2 = pyscf.ao2mo.kernel(pymol, c_act, aosym="s4",compact=false)
	h2 = reshape(h2, (nact, nact, nact, nact))

	# The use of d1_embed only really makes sense if it has zero electrons in the
	# active space. Let's warn the user if that's not true
	S = pymol.intor("int1e_ovlp_sph")
	n_act = tr(S * d1_embed * S * c_act * c_act')
	if isapprox(abs(n_act),0,atol=1e-8) == false
		println(n_act)
		display(d1_embed)
		error(" I found embedded electrons in the active space?!")
	end

	#println(size(e2))
	#println(h0)

	h1 = h + j - .5*k;
	#display(h + j - .5*k)

	h = InCoreInts(h0, h1, h2);
	return h
end

#"""
#	pyscf_build_ints(mol, c_act)
#
#build 1 and 2 electron integrals using a pyscf SCF object
#active is list of orbital indices which are active
#
#returns an `InCoreInts` type
#"""
#function pyscf_build_ints(mol, c_act)
#
#	pyscf = pyimport("pyscf")
#
#	nact = size(c_act)[2]
#	#mycas = pyscf.mcscf.CASSCF(mf, length(active), 0)
#
#	h0 = pyscf.gto.mole.energy_nuc(mol)
#	h1 = c_act' * pyscf.scf.hf.get_hcore(mol) * c_act
#	h2 = pyscf.ao2mo.kernel(mol, c_act, aosym="s4",compact=false)
#	h2 = reshape(h2, (nact, nact, nact, nact))
#
#	# println(size(c_act))
#	h = InCoreInts(h0, h1, h2);
#	return h
#end

"""
	pyscf_fci(ham, na, nb; max_cycle=20, conv_tol=1e-8, nroots=1, verbose=1)

Use PySCF to compute Full CI
"""
function pyscf_fci(ham, na, nb; max_cycle=20, conv_tol=1e-8, nroots=1, verbose=1)
	# println(" Use PYSCF to compute FCI")
	pyscf = pyimport("pyscf")
    pyscf.lib.num_threads(1)
	fci = pyimport("pyscf.fci")
	cisolver = pyscf.fci.direct_spin1.FCI()
	cisolver.max_cycle = max_cycle
	cisolver.conv_tol = conv_tol
	nelec = na + nb
	norb = size(ham.h1)[1]
	efci, ci = cisolver.kernel(ham.h1, ham.h2, norb , (na,nb), ecore=0, nroots =nroots, verbose=100)
    #@printf(" Length of CI Vector: %i\n", length(ci[1]))
    #println(size(ci[1]))
	fci_dim = size(ci,1)*size(ci,2)
	# d1 = cisolver.make_rdm1(ci, norb, nelec)

	d1a,d1b = cisolver.make_rdm1s(ci, norb, (na,nb))

	d1,d2 = cisolver.make_rdm12(ci, norb, (na,nb))
	#@printf(" Energy2: %12.8f\n", FermiCG.compute_energy(ham.h0, ham.h1, ham.h2, d1a+d1b, d2))
	# print(" PYSCF 1RDM: ")
	F = eigen(d1)
	occs = F.values
	sum_n = sum(occs)
	# @printf(" Sum of diagonals = %12.8f\n", sum_n)
    if verbose > 1
        @printf(" Natural Orbital Occupations:\n")
        [@printf(" %4i %12.8f\n",i,occs[i]) for i in 1:size(occs)[1] ]
        @printf(" -----------------\n")
        @printf(" %4s %12.8f\n\n","sum",sum_n)
    end
    if verbose>1
        pretty_table(d1; formatters = ft_printf("%5.3f"), noheader=true)
    end
    if verbose>0
        @printf(" FCI:        %12.8f %12.8f \n", efci+ham.h0, efci)
    end

    return efci, d1a,d1b, d2, ci
end

function get_nuclear_rep(mol::Molecule)
    pyscf = pyimport("pyscf")
    pymol = FermiCG.make_pyscf_mole(mol)
    return pyscf.gto.mole.energy_nuc(pymol)
end

"""
	localize(C::Array{Float64,2},method::String, mf)

Localize the orbitals using method = `method`
"""
function localize(C::Array{Float64,2},method::String, mf)
    """
    mf is a pyscf scf object
    """
    pyscf = pyimport("pyscf")
    pyscf.lib.num_threads(1)
    pyscflo = pyimport("pyscf.lo")
    if lowercase(method) == "lowdin"
        Cl = mf.mol.intor("int1e_ovlp_sph")
        F = svd(Cl)
        # display(Cl - F.U * Diagonal(F.S) * F.Vt)
        # display(Cl - F.vectors * Diagonal(F.values) * F.vectors')
        return F.U * Diagonal(F.S.^(-.5)) * F.Vt
    elseif lowercase(method) == "boys"
        Cl = pyscflo.Boys(mf.mol, C).kernel(verbose=4)
        return Cl
    end
end

"""
	get_ovlp(mf)

Get overlap matrix from pyscf using mean-field object
"""
function get_ovlp(mf)
		return mf.mol.intor("int1e_ovlp_sph")
end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
using JuMP
using Ipopt
#using Hiop
using Printf
using DelimitedFiles
using SparseArrays
using acopf

function solve(opfmodel)
  optimize!(opfmodel)
  status = termination_status(opfmodel)
  if status != MOI.LOCALLY_SOLVED
      println("Could not solve the model to optimality.")
  end

  return opfmodel, status
end

function auglag_model(opfdata; max_iter=100, rho=10)
  Pg0,Qg0,Vm0,Va0 = acopf.initialPt_IPOPT(opfdata)
  lines = opfdata.lines; buses = opfdata.buses; generators = opfdata.generators; baseMVA = opfdata.baseMVA
  busIdx = opfdata.BusIdx; FromLines = opfdata.FromLines; ToLines = opfdata.ToLines; BusGeners = opfdata.BusGenerators;
  nbus=length(opfdata.buses); ngen=length(opfdata.generators)

  generators=opfdata.generators; buses=opfdata.buses; baseMVA = opfdata.baseMVA

  opfmodel = Model(optimizer_with_attributes(Ipopt.Optimizer, "max_iter" => max_iter))

  ncont = 1

  @variable(opfmodel, generators[i].Pmin <= Pg[i=1:ngen] <= generators[i].Pmax, start = Pg0[i])
  @variable(opfmodel, 0<=extra[i=1:ngen,0:ncont]<=0.05*generators[i].Pmax, start = 0.025*generators[i].Pmax)
  @variable(opfmodel, generators[i].Qmin <= Qg[i=1:ngen] <= generators[i].Qmax, start = Qg0[i])

  @variable(opfmodel, buses[i].Vmin <= Vm[i=1:nbus, c=0:ncont] <= buses[i].Vmax, start = Vm0[i])
  @variable(opfmodel, Va[i=1:nbus, c=0:ncont], start = Va0[i])

  #fix the voltage angle at the reference bus
  for c in 0:ncont
    set_lower_bound(Va[opfdata.bus_ref,c], buses[opfdata.bus_ref].Va)
    set_upper_bound(Va[opfdata.bus_ref,c], buses[opfdata.bus_ref].Va)
  end

  @variable(opfmodel, generators[g].Pmin <= Pgc[g=1:ngen,co=0:ncont] <= generators[g].Pmax)
  @NLexpression(opfmodel, post_gen[g=1:ngen,co=0:ncont], Pgc[g,co] - (Pg[g] + extra[g,co]))

  #@constraint(opfmodel, ex[i=1:ngen,co=0:ncont], generators[i].Pmin <= Pg[i] + extra[i,co] <= generators[i].Pmax)
  #@NLexpression(opfmodel, Pgc[g=1:ngen, co=0:ncont], Pg[g] + extra[g,co])
  @NLexpression(opfmodel, Qgc[g=1:ngen, co=0:ncont], Qg[g])
  zeroexpr = @NLexpression(opfmodel, 0)
  set_lower_bound(Pgc[1,1], 0)
  set_upper_bound(Pgc[1,1], 0)
  # Pgc[1,1] = zeroexpr
  # Pgc[2,2] = zeroexpr
  # Pgc[3,3] = zeroexpr
  Qgc[1,1] = zeroexpr
  # Qgc[2,2] = zeroexpr
  # Qgc[3,3] = zeroexpr
  # @show Pgc

  @NLexpression(opfmodel, objective, sum( generators[i].coeff[generators[i].n-2]*(baseMVA*(Pgc[i,c]))^2
			             +generators[i].coeff[generators[i].n-1]*(baseMVA*(Pgc[i,c]))
				     +generators[i].coeff[generators[i].n  ] for i=1:ngen, c=0:ncont) / (1+ncont))

  YffR,YffI,YttR,YttI,YftR,YftI,YtfR,YtfI,YshR,YshI = computeAdmitances(opfdata.lines, opfdata.buses, opfdata.baseMVA)
  lines=opfdata.lines
  busIdx=opfdata.BusIdx;FromLines=opfdata.FromLines; ToLines=opfdata.ToLines; BusGeners=opfdata.BusGenerators
  nline=length(lines)

  # power flow balance
  @NLexpression( opfmodel, real_power[b=1:nbus,co=0:ncont],
        ( sum( YffR[l] for l in FromLines[b]) + sum( YttR[l] for l in ToLines[b]) + YshR[b] ) * Vm[b,co]^2
        + sum(  Vm[b,co]*Vm[busIdx[lines[l].to],  co]*( YftR[l]*cos(Va[b,co]-Va[busIdx[lines[l].to],  co])
              + YftI[l]*sin(Va[b,co]-Va[busIdx[lines[l].to],co]  )) for l in FromLines[b] )
        + sum(  Vm[b,co]*Vm[busIdx[lines[l].from],co]*( YtfR[l]*cos(Va[b,co]-Va[busIdx[lines[l].from],co])
              + YtfI[l]*sin(Va[b,co]-Va[busIdx[lines[l].from],co])) for l in ToLines[b]   )
        -(sum(baseMVA*Pgc[g,co] for g in BusGeners[b]) - buses[b].Pd) / baseMVA      # Sbus part
        )

  @NLexpression( opfmodel, reactive_power[b=1:nbus,co=0:ncont],
        ( sum(-YffI[l] for l in FromLines[b]) + sum(-YttI[l] for l in ToLines[b]) - YshI[b] ) * Vm[b,co]^2
        + sum(  Vm[b,co]*Vm[busIdx[lines[l].to],co]  *(-YftI[l]*cos(Va[b,co]-Va[busIdx[lines[l].to],  co])
              + YftR[l]*sin(Va[b,co]-Va[busIdx[lines[l].to],co]  )) for l in FromLines[b] )
        + sum( Vm[b,co]*Vm[busIdx[lines[l].from],co] *(-YtfI[l]*cos(Va[b,co]-Va[busIdx[lines[l].from],co])
              + YtfR[l]*sin(Va[b,co]-Va[busIdx[lines[l].from],co])) for l in ToLines[b]   )
        -(sum(baseMVA*Qgc[g] for g in BusGeners[b]) - buses[b].Qd) / baseMVA      #Sbus part
       )

  # branch/lines flow limits
  limidx = findall(x -> x.rateA != 0 && x.rateA < 1.0e10, lines)
  nlinelim = length(limidx)
  Yff_abs2 = YffR.^2 .+ YffI.^2
  Yft_abs2 = YftR.^2 .+ YftI.^2
  Yre_fr   =  YffR.*YftR .+ YffI.*YftI
  Yim_fr   = -YffR.*YftI .+ YffI.*YftR
  flowmax = [(lines[l].rateA / baseMVA)^2 for l in limidx]

  @variable(opfmodel, Sfrom[l=1:nlinelim] >= 0, start=1e4)
  @variable(opfmodel, Sto[l=1:nlinelim] >= 0, start=1e4)

  @NLexpression(opfmodel, powerlimit_from[l=1:nlinelim,co=0:ncont],
    Vm[busIdx[lines[limidx[l]].from],co]^2 *
    ( Yff_abs2[limidx[l]]*Vm[busIdx[lines[limidx[l]].from],co]^2 + Yft_abs2[limidx[l]]*Vm[busIdx[lines[limidx[l]].to],co]^2
    + 2*Vm[busIdx[lines[limidx[l]].from],co]*Vm[busIdx[lines[limidx[l]].to],co]
              *(Yre_fr[limidx[l]]*cos(Va[busIdx[lines[limidx[l]].from],co]-Va[busIdx[lines[limidx[l]].to],co])-
                Yim_fr[limidx[l]]*sin(Va[busIdx[lines[limidx[l]].from],co]-Va[busIdx[lines[limidx[l]].to],co]))
    ) - flowmax[l] + Sfrom[l]
  )

  Ytf_abs2 = YtfR.^2 .+ YtfI.^2
  Ytt_abs2 = YttR.^2 .+ YttI.^2
  Yre_to   =  YtfR.*YttR .+ YtfI.*YttI
  Yim_to   = -YtfR.*YttI .+ YtfI.*YttR

  @NLexpression(opfmodel, powerlimit_to[l=1:nlinelim,co=0:ncont],
    Vm[busIdx[lines[limidx[l]].to],co]^2 *
    ( Ytf_abs2[limidx[l]]*Vm[busIdx[lines[limidx[l]].from],co]^2 + Ytt_abs2[limidx[l]]*Vm[busIdx[lines[limidx[l]].to],co]^2
    + 2*Vm[busIdx[lines[limidx[l]].from],co]*Vm[busIdx[lines[limidx[l]].to],co]
        *(Yre_to[limidx[l]]*cos(Va[busIdx[lines[limidx[l]].from],co]-Va[busIdx[lines[limidx[l]].to],co])-
          Yim_to[limidx[l]]*sin(Va[busIdx[lines[limidx[l]].from],co]-Va[busIdx[lines[limidx[l]].to],co]))
    ) - flowmax[l] + Sto[l]
  )

  # minimize active power
#  @NLobjective(opfmodel,
#		  Min,
#		  sum( generators[i].coeff[generators[i].n] +
#		       sum(generators[i].coeff[generators[i].n-k]*(baseMVA*Pg[i])^k for k=1:generators[i].n-1)
#		       for i=1:ngen)
#		 )

  lambda_post_pg = ones((ngen,ncont+1))
  lambda_pg = ones((nbus,ncont+1))
  lambda_qg = ones((nbus,ncont+1))
  lambda_sfrom = ones((nlinelim,ncont+1))
  lambda_sto = ones((nlinelim,ncont+1))

  @NLobjective(opfmodel, Min,
               objective +
               sum(lambda_post_pg[g,co+1]*post_gen[g,co] for g=1:ngen,co=0:ncont) +
               sum(lambda_pg[b,co+1]*real_power[b,co] for b=1:nbus,co=0:ncont) +
               sum(lambda_qg[b,co+1]*reactive_power[b,co] for b=1:nbus,co=0:ncont) +
               sum(lambda_sfrom[l,co+1]*powerlimit_from[l,co] for l=1:nlinelim,co=0:ncont) +
               sum(lambda_sto[l,co+1]*powerlimit_to[l,co] for l=1:nlinelim,co=0:ncont) +
               (rho/2)*sum((post_gen[g,co])^2 for g=1:ngen,co=0:ncont) +
               (rho/2)*sum((real_power[b,co])^2 for b=1:nbus,co=0:ncont) +
               (rho/2)*sum((reactive_power[b,co])^2 for b=1:nbus,co=0:ncont) +
               (rho/2)*sum((powerlimit_from[l,co])^2 for l=1:nlinelim,co=0:ncont) +
               (rho/2)*sum((powerlimit_to[l,co])^2 for l=1:nlinelim,co=0:ncont)
               )

  return opfmodel
end

function model(opfdata; max_iter=100)
  Pg0,Qg0,Vm0,Va0 = acopf.initialPt_IPOPT(opfdata)
  lines = opfdata.lines; buses = opfdata.buses; generators = opfdata.generators; baseMVA = opfdata.baseMVA
  busIdx = opfdata.BusIdx; FromLines = opfdata.FromLines; ToLines = opfdata.ToLines; BusGeners = opfdata.BusGenerators;
  nbus=length(opfdata.buses); ngen=length(opfdata.generators)

  generators=opfdata.generators; buses=opfdata.buses; baseMVA = opfdata.baseMVA

  opfmodel = Model(optimizer_with_attributes(Ipopt.Optimizer, "max_iter" => max_iter))

  ncont = 1

  @variable(opfmodel, generators[i].Pmin <= Pg[i=1:ngen] <= generators[i].Pmax, start = Pg0[i])
  @variable(opfmodel, 0<=extra[i=1:ngen,0:ncont]<=0.05*generators[i].Pmax, start = 0.025*generators[i].Pmax)
  @variable(opfmodel, generators[i].Qmin <= Qg[i=1:ngen] <= generators[i].Qmax, start = Qg0[i])

  @variable(opfmodel, buses[i].Vmin <= Vm[i=1:nbus, c=0:ncont] <= buses[i].Vmax, start = Vm0[i])
  @variable(opfmodel, Va[i=1:nbus, c=0:ncont], start = Va0[i])

  #fix the voltage angle at the reference bus
  for c in 0:ncont
    set_lower_bound(Va[opfdata.bus_ref,c], buses[opfdata.bus_ref].Va)
    set_upper_bound(Va[opfdata.bus_ref,c], buses[opfdata.bus_ref].Va)
  end

  @constraint(opfmodel, ex[i=1:ngen,co=0:ncont], generators[i].Pmin <= Pg[i] + extra[i,co] <= generators[i].Pmax)
  @NLexpression(opfmodel, Pgc[g=1:ngen, co=0:ncont], Pg[g] + extra[g,co])
  @NLexpression(opfmodel, Qgc[g=1:ngen, co=0:ncont], Qg[g])
  zeroexpr = @NLexpression(opfmodel, 0)
  Pgc[1,1] = zeroexpr
  # Pgc[2,2] = zeroexpr
  # Pgc[3,3] = zeroexpr
  Qgc[1,1] = zeroexpr
  # Qgc[2,2] = zeroexpr
  # Qgc[3,3] = zeroexpr
  # @show Pgc

  @NLobjective(opfmodel, Min, sum( generators[i].coeff[generators[i].n-2]*(baseMVA*(Pgc[i,c]))^2
			             +generators[i].coeff[generators[i].n-1]*(baseMVA*(Pgc[i,c]))
				     +generators[i].coeff[generators[i].n  ] for i=1:ngen, c=0:ncont) / (1+ncont))

  for co=0:ncont
    #branch admitances
    # if co==0
      YffR,YffI,YttR,YttI,YftR,YftI,YtfR,YtfI,YshR,YshI = computeAdmitances(opfdata.lines, opfdata.buses, opfdata.baseMVA)
      lines=opfdata.lines
      busIdx=opfdata.BusIdx;FromLines=opfdata.FromLines; ToLines=opfdata.ToLines; BusGeners=opfdata.BusGenerators
    # else
    #   opfm1=opf_loaddata(ARGS[1], sd.lines_off[co])
    #   YffR,YffI,YttR,YttI,YftR,YftI,YtfR,YtfI,YshR,YshI = computeAdmitances(opfm1.lines, opfm1.buses, opfm1.baseMVA)
    #   lines=opfm1.lines
    #   busIdx = opfm1.BusIdx; FromLines = opfm1.FromLines; ToLines = opfm1.ToLines; BusGeners = opfm1.BusGenerators
    # end
    nline=length(lines)

    # power flow balance
    for b in 1:nbus
      #real part
      @NLconstraint( opfmodel,
        ( sum( YffR[l] for l in FromLines[b]) + sum( YttR[l] for l in ToLines[b]) + YshR[b] ) * Vm[b,co]^2
        + sum(  Vm[b,co]*Vm[busIdx[lines[l].to],  co]*( YftR[l]*cos(Va[b,co]-Va[busIdx[lines[l].to],  co])
              + YftI[l]*sin(Va[b,co]-Va[busIdx[lines[l].to],co]  )) for l in FromLines[b] )
        + sum(  Vm[b,co]*Vm[busIdx[lines[l].from],co]*( YtfR[l]*cos(Va[b,co]-Va[busIdx[lines[l].from],co])
              + YtfI[l]*sin(Va[b,co]-Va[busIdx[lines[l].from],co])) for l in ToLines[b]   )
        -(sum(baseMVA*Pgc[g,co] for g in BusGeners[b]) - buses[b].Pd) / baseMVA      # Sbus part
        ==0)

      #imaginary part
      @NLconstraint( opfmodel,
        ( sum(-YffI[l] for l in FromLines[b]) + sum(-YttI[l] for l in ToLines[b]) - YshI[b] ) * Vm[b,co]^2
        + sum(  Vm[b,co]*Vm[busIdx[lines[l].to],co]  *(-YftI[l]*cos(Va[b,co]-Va[busIdx[lines[l].to],  co])
              + YftR[l]*sin(Va[b,co]-Va[busIdx[lines[l].to],co]  )) for l in FromLines[b] )
        + sum( Vm[b,co]*Vm[busIdx[lines[l].from],co] *(-YtfI[l]*cos(Va[b,co]-Va[busIdx[lines[l].from],co])
              + YtfR[l]*sin(Va[b,co]-Va[busIdx[lines[l].from],co])) for l in ToLines[b]   )
        -(sum(baseMVA*Qgc[g] for g in BusGeners[b]) - buses[b].Qd) / baseMVA      #Sbus part
        ==0)
    end # of for: power flow balance loop

    # branch/lines flow limits

    nlinelim=0
    for l in 1:nline
      if lines[l].rateA!=0 && lines[l].rateA<1.0e10
        nlinelim += 1
        flowmax=(lines[l].rateA/baseMVA)^2

        #branch apparent power limits (from bus)
        Yff_abs2=YffR[l]^2+YffI[l]^2; Yft_abs2=YftR[l]^2+YftI[l]^2
        Yre=YffR[l]*YftR[l]+YffI[l]*YftI[l]; Yim=-YffR[l]*YftI[l]+YffI[l]*YftR[l]
        @NLconstraint( opfmodel,
	  Vm[busIdx[lines[l].from],co]^2 *
	  ( Yff_abs2*Vm[busIdx[lines[l].from],co]^2 + Yft_abs2*Vm[busIdx[lines[l].to],co]^2
	  + 2*Vm[busIdx[lines[l].from],co]*Vm[busIdx[lines[l].to],co]
               *(Yre*cos(Va[busIdx[lines[l].from],co]-Va[busIdx[lines[l].to],co])-Yim*sin(Va[busIdx[lines[l].from],co]-Va[busIdx[lines[l].to],co]))
	  )
          - flowmax <=0)

        #branch apparent power limits (to bus)
        Ytf_abs2=YtfR[l]^2+YtfI[l]^2; Ytt_abs2=YttR[l]^2+YttI[l]^2
        Yre=YtfR[l]*YttR[l]+YtfI[l]*YttI[l]; Yim=-YtfR[l]*YttI[l]+YtfI[l]*YttR[l]
        @NLconstraint(
          opfmodel,
    	  Vm[busIdx[lines[l].to],co]^2 *
          ( Ytf_abs2*Vm[busIdx[lines[l].from],co]^2 + Ytt_abs2*Vm[busIdx[lines[l].to],co]^2
          + 2*Vm[busIdx[lines[l].from],co]*Vm[busIdx[lines[l].to],co]
               *(Yre*cos(Va[busIdx[lines[l].from],co]-Va[busIdx[lines[l].to],co])-Yim*sin(Va[busIdx[lines[l].from],co]-Va[busIdx[lines[l].to],co]))
          )
          - flowmax <=0)
      end
    end # of for: branch power limits

    @printf("Contingency %d -> Buses: %d  Lines: %d  Generators: %d\n", co, nbus, nline, ngen)
    println("     lines with limits:  ", nlinelim)

  end # of for: contingencies

  # minimize active power
#  @NLobjective(opfmodel,
#		  Min,
#		  sum( generators[i].coeff[generators[i].n] +
#		       sum(generators[i].coeff[generators[i].n-k]*(baseMVA*Pg[i])^k for k=1:generators[i].n-1)
#		       for i=1:ngen)
#		 )

  return opfmodel
end

  #######################################################

  #values = zeros(2*nbus+2*ngen)
  ## values[1:2*nbus+2*ngen] = readdlm("/sandbox/petra/work/installs/matpower5.1/vars2.txt")
  #values[1:2*nbus+2*ngen] = readdlm("/sandbox/petra/work/installs/matpower5.1/vars3_9241.txt")
  #d = JuMP.NLPEvaluator(opfmodel)
  #MathProgBase.initialize(d, [:Jac])

  #g = zeros(2*nbus+2*nlinelim)
  #MathProgBase.eval_g(d, g, values)
  #println("f=", MathProgBase.eval_f(d,values))

  #gmat=zeros(2*nbus+2*nlinelim)
  #gmat[1:end] = readdlm("/sandbox/petra/work/installs/matpower5.1/cons3_9241.txt")
  #println("diff: ", norm(gmat-g))

  #println(opfmodel)

  #############################################################

function scopf_outputAll(opfmodel, scopf_data)
  #shortcuts for compactness
  sd=scopf_data; opf_data=sd.opfdata
  lines = opf_data.lines; buses = opf_data.buses; generators = opf_data.generators; baseMVA = opf_data.baseMVA
  busIdx = opf_data.BusIdx; FromLines = opf_data.FromLines; ToLines = opf_data.ToLines; BusGeners = opf_data.BusGenerators;

  nbus  = length(buses); nline = length(lines); ngen  = length(generators)

  # OUTPUTING
  println("Objective value: ", objective_value(opfmodel), "USD/hr")
  VM=value.(getindex(opfmodel,:Vm)); VA=value.(getindex(opfmodel,:Va));
  PG=value.(getindex(opfmodel,:Pg)); QG=value.(getindex(opfmodel,:Qg));

  VM=VM[:,0]; VA=VA[:,0]; #base case

  EX=value.(getindex(opfmodel,:extra));
  EX=EX[:,0];

  # printing the first stage variables
  println("============================= BUSES ==================================")
  println("  Generator  |  extra ")   # |    P (MW)     Q (MVAr)")  #|         (load)   ")
  println("----------------------------------------------------------------------")
  for i in 1:ngen
      @printf("  %10d | %6.2f \n",generators[i].bus, EX[i])
  end
  println("\n")

  println("============================= BUSES ==================================")
  println("  BUS    Vm     Va   |   Pg (MW)    Qg(MVAr) ")   # |    P (MW)     Q (MVAr)")  #|         (load)   ")

  println("                     |     (generation)      ")
  println("----------------------------------------------------------------------")
  for i in 1:nbus
    @printf("%4d | %6.2f  %6.2f | %s  | \n",
	    buses[i].bus_i, VM[i], VA[i]*180/pi,
	    length(BusGeners[i])==0 ? "   --          --  " : @sprintf("%7.2f     %7.2f", baseMVA*PG[BusGeners[i][1]], baseMVA*QG[BusGeners[i][1]]))
  end
  println("\n")

  within=20 # percentage close to the limits

  nflowlim=0
  for l in 1:nline
    if lines[l].rateA!=0 && lines[l].rateA<1.0e10
      nflowlim += 1
    end
  end
  return

  if nflowlim>0
    println("Number of lines with flow limits: ", nflowlim)

    optvec=zeros(2*nbus+2*ngen)
    optvec[1:ngen]=PG
    optvec[ngen+1:2*ngen]=QG
    optvec[2*ngen+1:2*ngen+nbus]=VM
    optvec[2*ngen+nbus+1:2*ngen+2*nbus]=VA

    d = JuMP.NLPEvaluator(opfmodel)
    MOI.initialize(d, [:Jac])

    consRhs = zeros(2*nbus+2*nflowlim)
    MOI.eval_constraint(d, consRhs, optvec)

    #println(consRhs)

    @printf("================ Lines within %d %s of flow capacity ===================\n", within, "%")
    println("Line   From Bus    To Bus    At capacity")

    nlim=1
    for l in 1:nline
      if lines[l].rateA!=0 && lines[l].rateA<1.0e10
        flowmax=(lines[l].rateA/baseMVA)^2
        idx = 2*nbus+nlim

        if( (consRhs[idx]+flowmax)  >= (1-within/100)^2*flowmax )
          @printf("%3d      %3d      %3d        %5.3f%s\n", l, lines[l].from, lines[l].to, 100*sqrt((consRhs[idx]+flowmax)/flowmax), "%" )
          #@printf("%7.4f   %7.4f    %7.4f \n", consRhs[idx], consRhs[idx]+flowmax,  flowmax)
        end
        nlim += 1
      end
    end
  end

  return
end

function main(args)
  if length(args) < 1
    println("Usage: julia scopf.jl case_name")
    println("Cases are in 'data' directory: case9 case30 case57 case118 case300 case1354pegase case2383wp case2736sp case2737sop case2746wop case2869pegase case3012wp  case3120sp case3375wp case9241pegase")
    return
  end

  max_iter = 100
  opfdata = acopf.opf_loaddata(args)

  scopfmodel = auglag_model(opfdata; rho=10)
  opfmodel, status = solve(scopfmodel)

  d = NLPEvaluator(scopfmodel)
  MOI.initialize(d, [:Grad, :Jac, :Hess])
  op_model = backend(scopfmodel).optimizer.model
  inner = backend(scopfmodel).optimizer.model.inner
  nvar = length(inner.x)

  IJ = MOI.hessian_lagrangian_structure(d)
  Hv_tmp = zeros(length(IJ))
  MOI.eval_hessian_lagrangian(d, Hv_tmp, inner.x, 1.0, inner.mult_g)

  # Make it upper-triangular.
  I = zeros(length(IJ))
  J = zeros(length(IJ))
  for i=1:length(IJ)
    I[i] = IJ[i][1]
    J[i] = IJ[i][2]
    if I[i] > J[i]
      I[i], J[i] = J[i], I[i]
    end
  end

  # Merge duplicates.
  Ih, Jh, Vh = findnz(sparse(I, J, [Int[i] for i=1:length(I)], nvar, nvar, vcat))
  Hv = [sum(Hv_tmp[Vh[i]]) for i=1:length(Ih)]

  return scopfmodel, Hv

  # if status == MOI.LOCALLY_SOLVED
  #   scopf_outputAll(opfmodel,scopfdata)
  # end
end