text:
## globals

function init_console!(v, b, main_window)

    set_gtk_property!(b, :style_scheme, main_window.style_and_language_manager.main_style)

    highlight_matching_brackets(b, true)

    show_line_numbers!(v, false)
    auto_indent!(v, true)
    highlight_current_line!(v, true)
    set_gtk_property!(v, :wrap_mode, 1)
    #set_gtk_property!(v, :expand, true)

    set_gtk_property!(v, :tab_width, 4)
    set_gtk_property!(v, :insert_spaces_instead_of_tabs, true)

    set_gtk_property!(v, :margin_bottom, 10)

    style_css(v, style_provider(main_window))

end

function run()#__init__()

    global is_running = true #should probably use g_main_loop_is_running or something of the sort
    global default_settings = init_opt()
    global main_window = MainWindow()

    ## Console

    console_mng = ConsoleManager(main_window, GtkIDE)
    GtkREPL.init!(console_mng)
    init_console_commands()#insert console commands into GtkREPL

    ## Editor

    editor = Editor(main_window)

    search_window = SearchWindow(editor)
    init!(search_window)
    visible(search_window, false)

    init!(editor, search_window)

    upgrade_project()

    global project = Project(main_window, "default")

    pathCBox = PathComboBox(main_window)
    statusBar = GtkStatusbar()

    menubar = MainMenu(main_window)

    global sidepanel_ntbook = SidePanelManager(main_window)

    init!(main_window, editor, console_mng, pathCBox, statusBar, project, menubar, sidepanel_ntbook)
    GtkREPL.set_main_window(main_window)

    load(project)
    cd(project.path)
    load_tabs(editor, project)

    ## Ploting window

    global fig_ntbook = GtkNotebook()
    global _display = Immerse._display

    #FIXME need init!
    signal_connect(fig_ntbook_key_press_cb, fig_ntbook, "key-press-event", Cint, (Ptr{Gtk.GdkEvent}, ), false)
    signal_connect(fig_ntbook_switch_page_cb, fig_ntbook, "switch-page", Nothing, (Ptr{Gtk.GtkWidget}, Int32), false)

    ## completion window

    global completion_window = CompletionWindow(main_window)
    visible(completion_window, false)

    ## Main layout
    global mainPan = GtkPaned(:h)
    rightPan = GtkPaned(:v)

    main_window |>
        ((mainVbox = GtkBox(:v)) |>
            menubar |>
            (topBarBox = GtkBox(:h) |>
                (sidePanelButton = GtkButton("F1")) |>
                 pathCBox   |>
                (editorButton = GtkButton("F2"))
            ) |>
            (global sidePan = GtkPaned(:h)) |>
            statusBar
        )

    mainPan |>
        (rightPan |>
            #(canvas = GtkCanvas())  |>
            (fig_ntbook)  |>
            console_mng
        ) |>
        ((editorVBox = GtkBox(:v)) |>
            ((editorBox = GtkBox(:h)) |>
                editor |>
                editor.sourcemap
            ) |>
            search_window
        )

    @info 1
    sidePan |>
        sidepanel_ntbook |>
        mainPan

    # Console
    lang = main_window.style_and_language_manager.languageDefinitions[".jl"]
    RemoteGtkREPL.estalbish_connection(console_mng.port, 1, "GtkIDE")

    #@assert length(console_mng) == 1

    set_gtk_property!(statusBar, :margin, 2)
    push!(statusBar, "main", "Julia $VERSION")
    Gtk.G_.position(sidePan, 160)

    set_gtk_property!(editor, :vexpand, true)
    set_gtk_property!(editorBox, :expand, editor, true)
    set_gtk_property!(mainPan, :margin, 0)
    Gtk.G_.position(mainPan, 600)
    Gtk.G_.position(rightPan, 450)
    #-

    ## set some style

    nbtbookcss =
    "* tab {
        padding:0px;
        padding-left:6px;
        padding-right:1px;
        margin:0px;
        font-size:$(main_window.style_and_language_manager.fontsize-1)px;
    }"
    style_css(main_window.editor, nbtbookcss)
    style_css(main_window.console_manager, nbtbookcss)
    style_css(fig_ntbook, nbtbookcss)
    style_css(sidepanel_ntbook, nbtbookcss)

    set_gtk_property!(topBarBox, :hexpand, true)

    ################
    # Side Panels

    #FIXME globals
    global filespanel = FilesPanel(main_window)
    global workspacepanel = WorkspacePanel(main_window)
    global projectspanel = ProjectsPanel(main_window)

    init!(sidepanel_ntbook, [filespanel, workspacepanel, projectspanel], ["F", "W", "P"])

    ################
    ## Plots
    #GtkREPL.gadfly()

    sleep(0.01)
    figure()
    drawnow() = sleep(0.001)

    init!(pathCBox)#need on_path_change to be defined

    signal_connect(sidePanelButton_clicked_cb, sidePanelButton, "clicked", Nothing, (), false)
    signal_connect(editorButtonclicked_cb, editorButton, "clicked", Nothing, (), false)

    showall(main_window)
    visible(search_window, false)
    visible(sidepanel_ntbook, false)
    visible(editor.sourcemap, opt("Editor", "show_source_map"))

    ## starting task and such

    sleep(0.01)

    if false

        global stdout_ = stdout
        global stderr_ = stderr

        read_stdout, wr = redirect_stdout()
        #read_stderr, wre = redirect_stderr()

        function watch_stdout()
            @async GtkREPL.watch_stream(read_stdout, console)
        end
        function watch_stderr()
            @async GtkREPL.watch_stream(read_stderr, console)
        end

        global watch_stdout_task = watch_stdout()
        #global watch_stderr_task = watch_stderr()

       GtkREPL.init_stdout!(main_window.console_manager, watch_stdout_task, stdout_, stderr_)

        g_timeout_add(()->print_to_console(console), 10)
    end

    println("Warming up, hold on...")
    sleep(0.05)
    #@async logo()
    #new_prompt(console)

end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule pigpiod_if2_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("pigpiod_if2")
JLLWrappers.@generate_main_file("pigpiod_if2", UUID("50515f5d-e8ff-5b26-ab2b-0473a017b18c"))
end  # module pigpiod_if2_jll

ext: jl


text:
struct Rosenbrock23Tableau{T}
  c₃₂::T
  d::T
end

function Rosenbrock23Tableau(T)
  c₃₂ = convert(T,6 + sqrt(2))
  d = convert(T,1/(2+sqrt(2)))
  Rosenbrock23Tableau(c₃₂,d)
end

struct Rosenbrock32Tableau{T}
  c₃₂::T
  d::T
end

function Rosenbrock32Tableau(T)
  c₃₂ = convert(T,6 + sqrt(2))
  d = convert(T,1/(2+sqrt(2)))
  Rosenbrock32Tableau(c₃₂,d)
end

struct ROS3PTableau{T,T2}
  a21::T
  a31::T
  a32::T
  C21::T
  C31::T
  C32::T
  b1::T
  b2::T
  b3::T
  btilde1::T
  btilde2::T
  btilde3::T
  gamma::T2
  c2::T2
  c3::T2
  d1::T
  d2::T
  d3::T
end

function ROS3PTableau(T, T2)
  gamma = convert(T,1/2 + sqrt(3)/6)
  igamma = inv(gamma)
  a21 = convert(T,igamma)
  a31 = convert(T,igamma)
  a32 = convert(T,0)
  C21 = convert(T,-igamma^2)
  tmp = -igamma*(2 - (1/2)*igamma)
  C31 = -igamma*(1-tmp)
  C32 = tmp
  tmp = igamma*(2/3 - (1/6)*igamma)
  b1 = igamma*(1+tmp)
  b2 = tmp
  b3 = (1/3)*igamma
  # btilde1 = convert(T,2.113248654051871)
  # btilde2 = convert(T,1.000000000000000)
  # btilde3 = convert(T,0.4226497308103742)
  btilde1 = b1 - convert(T,2.113248654051871)
  btilde2 = b2 - convert(T,1.000000000000000)
  btilde3 = b3 - convert(T,0.4226497308103742)
  c2 = convert(T,1)
  c3 = convert(T,1)
  d1 = convert(T,0.7886751345948129)
  d2 = convert(T,-0.2113248654051871)
  d3 = convert(T,-1.077350269189626)
  ROS3PTableau(a21,a31,a32,C21,C31,C32,b1,b2,b3,btilde1,btilde2,btilde3,gamma,c2,c3,d1,d2,d3)
end

struct Rodas3Tableau{T,T2}
  a21::T
  a31::T
  a32::T
  a41::T
  a42::T
  a43::T
  C21::T
  C31::T
  C32::T
  C41::T
  C42::T
  C43::T
  b1::T
  b2::T
  b3::T
  b4::T
  btilde1::T
  btilde2::T
  btilde3::T
  btilde4::T
  gamma::T2
  c2::T2
  c3::T2
  d1::T
  d2::T
  d3::T
  d4::T
end

function Rodas3Tableau(T, T2)
  gamma = convert(T,1//2)
  a21 = convert(T,0)
  a31 = convert(T,2)
  a32 = convert(T,0)
  a41 = convert(T,2)
  a42 = convert(T,0)
  a43 = convert(T,1)
  C21 = convert(T,4)
  C31 = convert(T,1)
  C32 = convert(T,-1)
  C41 = convert(T,1)
  C42 = convert(T,-1)
  C43 = convert(T,-8//3)
  b1 = convert(T,2)
  b2 = convert(T,0)
  b3 = convert(T,1)
  b4 = convert(T,1)
  btilde1 = convert(T,0.0)
  btilde2 = convert(T,0.0)
  btilde3 = convert(T,0.0)
  btilde4 = convert(T,1.0)
  c2 = convert(T,0.0)
  c3 = convert(T,1.0)
  c4 = convert(T,1.0)
  d1 = convert(T,1//2)
  d2 = convert(T,3//2)
  d3 = convert(T,0)
  d4 = convert(T,0)
  Rodas3Tableau(a21,a31,a32,a41,a42,a43,C21,C31,C32,C41,C42,C43,b1,b2,b3,b4,btilde1,btilde2,btilde3,btilde4,gamma,c2,c3,d1,d2,d3,d4)
end

@ROS34PW(:tableau)

@Rosenbrock4(:tableau)

struct RodasTableau{T,T2}
  a21::T
  a31::T
  a32::T
  a41::T
  a42::T
  a43::T
  a51::T
  a52::T
  a53::T
  a54::T
  C21::T
  C31::T
  C32::T
  C41::T
  C42::T
  C43::T
  C51::T
  C52::T
  C53::T
  C54::T
  C61::T
  C62::T
  C63::T
  C64::T
  C65::T
  gamma::T
  c2::T2
  c3::T2
  c4::T2
  d1::T
  d2::T
  d3::T
  d4::T
  h21::T
  h22::T
  h23::T
  h24::T
  h25::T
  h31::T
  h32::T
  h33::T
  h34::T
  h35::T
end

function Rodas4Tableau(T, T2)
  gamma=convert(T,1//4)
  #BET2P=0.0317D0
  #BET3P=0.0635D0
  #BET4P=0.3438D0
  a21= convert(T,1.544000000000000)
  a31= convert(T,0.9466785280815826)
  a32= convert(T,0.2557011698983284)
  a41= convert(T,3.314825187068521)
  a42= convert(T,2.896124015972201)
  a43= convert(T,0.9986419139977817)
  a51= convert(T,1.221224509226641)
  a52= convert(T,6.019134481288629)
  a53= convert(T,12.53708332932087)
  a54=-convert(T,0.6878860361058950)
  C21=-convert(T,5.668800000000000)
  C31=-convert(T,2.430093356833875)
  C32=-convert(T,0.2063599157091915)
  C41=-convert(T,0.1073529058151375)
  C42=-convert(T,9.594562251023355)
  C43=-convert(T,20.47028614809616)
  C51= convert(T,7.496443313967647)
  C52=-convert(T,10.24680431464352)
  C53=-convert(T,33.99990352819905)
  C54= convert(T,11.70890893206160)
  C61= convert(T,8.083246795921522)
  C62=-convert(T,7.981132988064893)
  C63=-convert(T,31.52159432874371)
  C64= convert(T,16.31930543123136)
  C65=-convert(T,6.058818238834054)

  c2=convert(T2,0.386)
  c3=convert(T2,0.21)
  c4=convert(T2,0.63)

  d1= convert(T,0.2500000000000000)
  d2=-convert(T,0.1043000000000000)
  d3= convert(T,0.1035000000000000)
  d4=-convert(T,0.03620000000000023)

  h21= convert(T,10.12623508344586)
  h22=-convert(T,7.487995877610167)
  h23=-convert(T,34.80091861555747)
  h24=-convert(T,7.992771707568823)
  h25= convert(T,1.025137723295662)
  h31=-convert(T,0.6762803392801253)
  h32= convert(T,6.087714651680015)
  h33= convert(T,16.43084320892478)
  h34= convert(T,24.76722511418386)
  h35=-convert(T,6.594389125716872)

  RodasTableau(a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,
                    C21,C31,C32,C41,C42,C43,C51,C52,C53,C54,C61,C62,C63,C64,C65,
                    gamma,c2,c3,c4,d1,d2,d3,d4,
                    h21,h22,h23,h24,h25,h31,h32,h33,h34,h35)
end

function Rodas42Tableau(T, T2)
  gamma= convert(T,1//4)
  #BET2P=0.0317D0
  #BET3P=0.0047369D0
  #BET4P=0.3438D0
  a21= convert(T,1.402888400000000)
  a31= convert(T,0.6581212688557198)
  a32=-convert(T,1.320936088384301)
  a41= convert(T,7.131197445744498)
  a42= convert(T,16.02964143958207)
  a43=-convert(T,5.561572550509766)
  a51= convert(T,22.73885722420363)
  a52= convert(T,67.38147284535289)
  a53=-convert(T,31.21877493038560)
  a54= convert(T,0.7285641833203814)
  C21=-convert(T,5.104353600000000)
  C31=-convert(T,2.899967805418783)
  C32= convert(T,4.040399359702244)
  C41=-convert(T,32.64449927841361)
  C42=-convert(T,99.35311008728094)
  C43= convert(T,49.99119122405989)
  C51=-convert(T,76.46023087151691)
  C52=-convert(T,278.5942120829058)
  C53= convert(T,153.9294840910643)
  C54= convert(T,10.97101866258358)
  C61=-convert(T,76.29701586804983)
  C62=-convert(T,294.2795630511232)
  C63= convert(T,162.0029695867566)
  C64= convert(T,23.65166903095270)
  C65=-convert(T,7.652977706771382)
  c2=convert(T2,0.3507221)
  c3=convert(T2,0.2557041)
  c4=convert(T2,0.6817790)
  d1= convert(T,0.2500000000000000)
  d2=-convert(T,0.06902209999999998)
  d3=-convert(T,0.0009671999999999459)
  d4=-convert(T,0.08797900000000025)

  h21=-convert(T,38.71940424117216)
  h22=-convert(T,135.8025833007622)
  h23= convert(T,64.51068857505875)
  h24=-convert(T,4.192663174613162)
  h25=-convert(T,2.531932050335060)
  h31=-convert(T,14.99268484949843)
  h32=-convert(T,76.30242396627033)
  h33= convert(T,58.65928432851416)
  h34= convert(T,16.61359034616402)
  h35=-convert(T,0.6758691794084156)

  RodasTableau(a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,
                    C21,C31,C32,C41,C42,C43,C51,C52,C53,C54,C61,C62,C63,C64,C65,
                    gamma,c2,c3,c4,d1,d2,d3,d4,
                    h21,h22,h23,h24,h25,h31,h32,h33,h34,h35)
end

function Rodas4PTableau(T, T2)
  gamma = convert(T,1//4)
  #BET2P=0.D0
  #BET3P=c3*c3*(c3/6.d0-GAMMA/2.d0)/(GAMMA*GAMMA)
  #BET4P=0.3438D0
  a21= convert(T,3)
  a31= convert(T,1.831036793486759)
  a32= convert(T,0.4955183967433795)
  a41= convert(T,2.304376582692669)
  a42=-convert(T,0.05249275245743001)
  a43=-convert(T,1.176798761832782)
  a51=-convert(T,7.170454962423024)
  a52=-convert(T,4.741636671481785)
  a53=-convert(T,16.31002631330971)
  a54=-convert(T,1.062004044111401)
  C21=-convert(T,12)
  C31=-convert(T,8.791795173947035)
  C32=-convert(T,2.207865586973518)
  C41= convert(T,10.81793056857153)
  C42= convert(T,6.780270611428266)
  C43= convert(T,19.53485944642410)
  C51= convert(T,34.19095006749676)
  C52= convert(T,15.49671153725963)
  C53= convert(T,54.74760875964130)
  C54= convert(T,14.16005392148534)
  C61= convert(T,34.62605830930532)
  C62= convert(T,15.30084976114473)
  C63= convert(T,56.99955578662667)
  C64= convert(T,18.40807009793095)
  C65=-convert(T,5.714285714285717)
  c2=convert(T2,3*gamma)
  c3=convert(T2,0.21)
  c4=convert(T2,0.63)
  d1=convert(T, 0.2500000000000000)
  d2=convert(T,-0.5000000000000000)
  d3=convert(T,-0.0235040000000000)
  d4=convert(T,-0.0362000000000000)

  h21= convert(T,25.09876703708589)
  h22= convert(T,11.62013104361867)
  h23= convert(T,28.49148307714626)
  h24=-convert(T,5.664021568594133)
  h25= convert(T,0)
  h31= convert(T,1.638054557396973)
  h32=-convert(T,0.7373619806678748)
  h33= convert(T,8.477918219238990)
  h34= convert(T,15.99253148779520)
  h35=-convert(T,1.882352941176471)

  RodasTableau(a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,
                    C21,C31,C32,C41,C42,C43,C51,C52,C53,C54,C61,C62,C63,C64,C65,
                    gamma,c2,c3,c4,d1,d2,d3,d4,
                    h21,h22,h23,h24,h25,h31,h32,h33,h34,h35)
end

function Rodas4P2Tableau(T, T2)
  gamma = convert(T,1//4)
 a21= convert(T,3.000000000000000)
 a31= convert(T,0.906377755268814)
 a32=-convert(T,0.189707390391685)
 a41= convert(T,3.758617027739064)
 a42= convert(T,1.161741776019525)
 a43=-convert(T,0.849258085312803)
 a51= convert(T,7.089566927282776)
 a52= convert(T,4.573591406461604)
 a53=-convert(T,8.423496976860259)
 a54=-convert(T,0.959280113459775)

 C21= convert(T,-12.00000000000000)
 C31= convert(T,-6.354581592719008)
 C32= convert(T, 0.338972550544623)
 C41= convert(T,-8.575016317114033)
 C42= convert(T,-7.606483992117508)
 C43= convert(T,12.224997650124820)
 C51= convert(T,-5.888975457523102)
 C52= convert(T,-8.157396617841821)
 C53= convert(T,24.805546872612922)
 C54= convert(T,12.790401512796979)
 C61= convert(T,-4.408651676063871)
 C62= convert(T,-6.692003137674639)
 C63= convert(T,24.625568527593117)
 C64= convert(T,16.627521966636085)
 C65= convert(T,-5.714285714285718)

 c2=convert(T2,0.750000000000000)
 c3=convert(T2,0.321448134013046)
 c4=convert(T2,0.519745732277726)
 d1=convert(T, 0.250000000000000)
 d2=convert(T,-0.500000000000000)
 d3=convert(T,-0.189532918363016)
 d4=convert(T, 0.085612108792769)

 h21= convert(T, -5.323528268423303)
 h22= convert(T,-10.042123754867493)
 h23= convert(T, 17.175254928256965)
 h24= convert(T, -5.079931171878093)
 h25= convert(T, -0.016185991706112)
 h31= convert(T,  6.984505741529879)
 h32= convert(T,  6.914061169603662)
 h33= convert(T, -0.849178943070653)
 h34= convert(T, 18.104410789349338)
 h35= convert(T, -3.516963011559032)

  RodasTableau(a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,
                    C21,C31,C32,C41,C42,C43,C51,C52,C53,C54,C61,C62,C63,C64,C65,
                    gamma,c2,c3,c4,d1,d2,d3,d4,
                    h21,h22,h23,h24,h25,h31,h32,h33,h34,h35)
end

struct Rodas5Tableau{T,T2}
  a21::T
  a31::T
  a32::T
  a41::T
  a42::T
  a43::T
  a51::T
  a52::T
  a53::T
  a54::T
  a61::T
  a62::T
  a63::T
  a64::T
  a65::T
  C21::T
  C31::T
  C32::T
  C41::T
  C42::T
  C43::T
  C51::T
  C52::T
  C53::T
  C54::T
  C61::T
  C62::T
  C63::T
  C64::T
  C65::T
  C71::T
  C72::T
  C73::T
  C74::T
  C75::T
  C76::T
  C81::T
  C82::T
  C83::T
  C84::T
  C85::T
  C86::T
  C87::T
  gamma::T2
  d1::T
  d2::T
  d3::T
  d4::T
  d5::T
  c2::T2
  c3::T2
  c4::T2
  c5::T2
  h21::T
  h22::T
  h23::T
  h24::T
  h25::T
  h26::T
  h27::T
  h28::T
  h31::T
  h32::T
  h33::T
  h34::T
  h35::T
  h36::T
  h37::T
  h38::T
  h41::T
  h42::T
  h43::T
  h44::T
  h45::T
  h46::T
  h47::T
  h48::T
end

function Rodas5Tableau(T, T2)
  gamma = convert(T2,.19)
  a21 = convert(T,2.0)
  a31 = convert(T,3.040894194418781 )
  a32 = convert(T,1.041747909077569 )
  a41 = convert(T,2.576417536461461 )
  a42 = convert(T,1.622083060776640 )
  a43 = convert(T,-.9089668560264532)
  a51 = convert(T,2.760842080225597 )
  a52 = convert(T,1.446624659844071 )
  a53 = convert(T,-.3036980084553738)
  a54 = convert(T,.2877498600325443 )
  a61 = convert(T,-14.09640773051259)
  a62 = convert(T,6.925207756232704 )
  a63 = convert(T,-41.47510893210728)
  a64 = convert(T,2.343771018586405 )
  a65 = convert(T,24.13215229196062 )
  C21 = convert(T,-10.31323885133993)
  C31 = convert(T,-21.04823117650003)
  C32 = convert(T,-7.234992135176716)
  C41 = convert(T,32.22751541853323 )
  C42 = convert(T,-4.943732386540191)
  C43 = convert(T,19.44922031041879 )
  C51 = convert(T,-20.69865579590063)
  C52 = convert(T,-8.816374604402768)
  C53 = convert(T,1.260436877740897 )
  C54 = convert(T,-.7495647613787146)
  C61 = convert(T,-46.22004352711257)
  C62 = convert(T,-17.49534862857472)
  C63 = convert(T,-289.6389582892057)
  C64 = convert(T,93.60855400400906 )
  C65 = convert(T,318.3822534212147 )
  C71 = convert(T,34.20013733472935 )
  C72 = convert(T,-14.15535402717690)
  C73 = convert(T,57.82335640988400 )
  C74 = convert(T,25.83362985412365 )
  C75 = convert(T,1.408950972071624 )
  C76 = convert(T,-6.551835421242162)
  C81 = convert(T,42.57076742291101 )
  C82 = convert(T,-13.80770672017997)
  C83 = convert(T,93.98938432427124 )
  C84 = convert(T,18.77919633714503 )
  C85 = convert(T,-31.58359187223370)
  C86 = convert(T,-6.685968952921985)
  C87 = convert(T,-5.810979938412932)
  c2 = convert(T2,0.38              )
  c3 = convert(T2,0.3878509998321533)
  c4 = convert(T2,0.4839718937873840)
  c5 = convert(T2,0.4570477008819580)
  d1 = convert(T,gamma                 )
  d2 = convert(T,-0.1823079225333714636)
  d3 = convert(T,-0.319231832186874912 )
  d4 = convert(T, 0.3449828624725343   )
  d5 = convert(T,-0.377417564392089818 )

  h21= convert(T,27.354592673333357 )
  h22= convert(T,-6.925207756232857 )
  h23= convert(T,26.40037733258859 )
  h24= convert(T,0.5635230501052979 )
  h25= convert(T,-4.699151156849391 )
  h26= convert(T,-1.6008677469422725)
  h27= convert(T,-1.5306074446748028)
  h28= convert(T,-1.3929872940716344)

  h31= convert(T,44.19024239501722 )
  h32= convert(T,1.3677947663381929e-13)
  h33= convert(T,202.93261852171622 )
  h34= convert(T,-35.5669339789154 )
  h35= convert(T,-181.91095152160645 )
  h36= convert(T,3.4116351403665033 )
  h37= convert(T,2.5793540257308067 )
  h38= convert(T,2.2435122582734066)

  h41= convert(T,-44.0988150021747 )
  h42= convert(T,-5.755396159656812e-13)
  h43= convert(T,-181.26175034586677 )
  h44= convert(T,56.99302194811676 )
  h45= convert(T,183.21182741427398 )
  h46= convert(T,-7.480257918273637 )
  h47= convert(T,-5.792426076169686 )
  h48= convert(T,-5.32503859794143)

 # println("---Rodas5---")

  #=
  a71 = -14.09640773051259
  a72 = 6.925207756232704
  a73 = -41.47510893210728
  a74 = 2.343771018586405
  a75 = 24.13215229196062
  a76 = convert(T,1)
  a81 = -14.09640773051259
  a82 = 6.925207756232704
  a83 = -41.47510893210728
  a84 = 2.343771018586405
  a85 = 24.13215229196062
  a86 = convert(T,1)
  a87 = convert(T,1)
  b1 = -14.09640773051259
  b2 = 6.925207756232704
  b3 = -41.47510893210728
  b4 = 2.343771018586405
  b5 = 24.13215229196062
  b6 = convert(T,1)
  b7 = convert(T,1)
  b8 = convert(T,1)
  =#

  Rodas5Tableau(a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,
                      a61,a62,a63,a64,a65,
                      C21,C31,C32,C41,C42,C43,C51,C52,C53,C54,
                      C61,C62,C63,C64,C65,C71,C72,C73,C74,C75,C76,
                      C81,C82,C83,C84,C85,C86,C87,
                      gamma,d1,d2,d3,d4,d5,c2,c3,c4,c5,
                      h21,h22,h23,h24,h25,h26,h27,h28,h31,h32,h33,h34,h35,h36,h37,h38,h41,h42,h43,h44,h45,h46,h47,h48)
end

function Rodas5PTableau(T, T2)

gamma = convert(T2,0.21193756319429014)

a21 = convert(T,3.0 )
a31 = convert(T,2.849394379747939 )
a32 = convert(T,0.45842242204463923 )
a41 = convert(T,-6.954028509809101 )
a42 = convert(T,2.489845061869568 )
a43 = convert(T,-10.358996098473584 )
a51 = convert(T,2.8029986275628964 )
a52 = convert(T,0.5072464736228206 )
a53 = convert(T,-0.3988312541770524 )
a54 = convert(T,-0.04721187230404641 )
a61 = convert(T,-7.502846399306121 )
a62 = convert(T,2.561846144803919 )
a63 = convert(T,-11.627539656261098 )
a64 = convert(T,-0.18268767659942256 )
a65 = convert(T,0.030198172008377946 )

C21 = convert(T,-14.155112264123755 )
C31 = convert(T,-17.97296035885952 )
C32 = convert(T,-2.859693295451294 )
C41 = convert(T,147.12150275711716 )
C42 = convert(T,-1.41221402718213 )
C43 = convert(T,71.68940251302358 )
C51 = convert(T,165.43517024871676 )
C52 = convert(T,-0.4592823456491126 )
C53 = convert(T,42.90938336958603 )
C54 = convert(T,-5.961986721573306 )
C61 = convert(T,24.854864614690072 )
C62 = convert(T,-3.0009227002832186 )
C63 = convert(T,47.4931110020768 )
C64 = convert(T,5.5814197821558125 )
C65 = convert(T,-0.6610691825249471 )
C71 = convert(T,30.91273214028599 )
C72 = convert(T,-3.1208243349937974 )
C73 = convert(T,77.79954646070892 )
C74 = convert(T,34.28646028294783 )
C75 = convert(T,-19.097331116725623 )
C76 = convert(T,-28.087943162872662 )
C81 = convert(T,37.80277123390563 )
C82 = convert(T,-3.2571969029072276 )
C83 = convert(T,112.26918849496327 )
C84 = convert(T,66.9347231244047 )
C85 = convert(T,-40.06618937091002 )
C86 = convert(T,-54.66780262877968 )
C87 = convert(T,-9.48861652309627 )

c2 = convert(T2,0.6358126895828704)
c3 = convert(T2,0.4095798393397535)
c4 = convert(T2,0.9769306725060716)
c5 = convert(T2,0.4288403609558664)

d1 = convert(T,0.21193756319429014)
d2 = convert(T,-0.42387512638858027)
d3 = convert(T,-0.3384627126235924)
d4 = convert(T,1.8046452872882734)
d5 = convert(T,2.325825639765069)

h21= convert(T,25.948786856663858 )
h22= convert(T,-2.5579724845846235 )
h23= convert(T,10.433815404888879 )
h24= convert(T,-2.3679251022685204 )
h25= convert(T,0.524948541321073 )
h26= convert(T,1.1241088310450404 )
h27= convert(T,0.4272876194431874)
h28= convert(T,-0.17202221070155493)

h31= convert(T,-9.91568850695171 )
h32= convert(T,-0.9689944594115154 )
h33= convert(T,3.0438037242978453 )
h34= convert(T,-24.495224566215796 )
h35= convert(T,20.176138334709044 )
h36= convert(T,15.98066361424651 )
h37= convert(T,-6.789040303419874 )
h38= convert(T,-6.710236069923372)

h41= convert(T,11.419903575922262 )
h42= convert(T,2.8879645146136994 )
h43= convert(T,72.92137995996029 )
h44= convert(T,80.12511834622643 )
h45= convert(T,-52.072871366152654 )
h46= convert(T,-59.78993625266729 )
h47= convert(T,-0.15582684282751913)
h48= convert(T,4.883087185713722)

  Rodas5Tableau(a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,
                      a61,a62,a63,a64,a65,
                      C21,C31,C32,C41,C42,C43,C51,C52,C53,C54,
                      C61,C62,C63,C64,C65,C71,C72,C73,C74,C75,C76,
                      C81,C82,C83,C84,C85,C86,C87,
                      gamma,d1,d2,d3,d4,d5,c2,c3,c4,c5,
                      h21,h22,h23,h24,h25,h26,h27,h28,h31,h32,h33,h34,h35,h36,h37,h38,h41,h42,h43,h44,h45,h46,h47,h48)
end

@RosenbrockW6S4OS(:tableau)

#=
# alpha_ij
A = [0 0 0 0 0 0 0 0
     big"0.38" 0 0 0 0 0 0 0
     big"0.1899188971074152"    big"0.1979321027247381"  0 0 0 0 0 0
     big"0.1110729281178426"    big"0.5456026683145674"  big"-0.1727037026450261" 0 0 0 0 0
     big"0.2329444418850307"    big"0.025099380960713898" big"0.1443314046300300"  big"0.054672473406183418" 0 0 0 0
     big"-0.036201017843430883" big"4.208448872731939"   big"-7.549674427720996"  big"-0.2076823626400282" big"4.585108935472517" 0 0 0
     big"7.585261698003052"     big"-15.57426208319938"  big"-8.814406895608121"  big"1.534698996826085"   big"16.07870828397837" big"0.19" 0 0
     big"0.4646018839086969"    big"0"                   big"-1.720907508837576"  big"0.2910480220957973"  big"1.821778861539924" big"-0.046521258706842056" big"0.19" 0]
# bi
B = [big"0.464601884",0,big"-1.72090751",big"0.29104802",big"1.82177886",big"-0.02674488",big"-0.01977638",big"0.19"]

# Beta_ij

Beta = [big"0.19" 0 0 0 0 0 0 0
        big"0.0076920774666285364"  big"0.19" 0 0 0 0 0 0
        big"-0.058129718999580252"  big"-0.063251113355141360"  big"0.19" 0 0 0 0 0
        big"0.7075715596134048"     big"-0.5980299539145789"   big"0.5294131505610923"    big"0.19" 0 0 0 0
        big"-0.034975026573934865"  big"-0.1928476085817357"    big"0.089839586125126941"  big"0.027613185520411822" big"0.19" 0 0 0
        big"7.585261698003052"      big"-15.57426208319938"     big"-8.814406895608121"    big"1.534698996826085"    big"16.07870828397837"  big"0.19" 0 0
        big"0.4646018839086969"     0                           big"-1.720907508837576"    big"0.2910480220957973"   big"1.821778861539924"  big"-0.046521258706842056" big"0.19" 0
        big"0.4646018839086969"     0                           big"-1.720907508837576"    big"0.2910480220957973"   big"1.821778861539924"  big"-0.026744882930135193" big"-0.019776375776706864" big"0.19"]

Gamma = Beta - A
a = A*Gamma
m = B'*inv(Gamma) # = b_i
C = inv(Diagonal(diag(Gamma))) - inv(Gamma)
c = sum(A,2)
D = Beta*inv(Gamma)
d = sum(Gamma,2)

# Dense output
D_i = tanspose(D)*k_i
y1(Θ) = y₀*(1-Θ) + Θ*(y₁ + (Θ-1)*(D₂ + D₄ + (Θ + 1)*(D₃ + ΘD₄)))

# Determining coefficients
gamma = 0.19
c3 = 0.3878509998321533 == alpha3
c4 = 0.4839718937873840 == alpha4
c5 = 0.4570477008819580 == alpha5
beta3 = 6.8619167645278386e-2
beta4 = 0.8289547562599182
beta5 = 7.9630136489868164e-2
alpha64 = -0.2076823627400282
=#


text:
"""
Caching functionality for Genie.
"""
module Cache

import SHA, Logging
import Genie

export cachekey, withcache, @cachekey
export purge, purgeall

function init() :: Nothing
  @eval Genie.config.cache_storage === nothing && (Genie.config.cache_storage = :File)
  @eval Genie.config.cache_storage == :File && include(joinpath(@__DIR__, "cache_adapters", "FileCache.jl"))

  nothing
end

"""
    withcache(f::Function, key::Union{String,Symbol}, expiration::Int = Genie.config.cache_duration; dir = "", condition::Bool = true)

Executes the function `f` and stores the result into the cache for the duration (in seconds) of `expiration`. Next time the function is invoked,
if the cache has not expired, the cached result is returned skipping the function execution.
The optional `dir` param is used to designate the folder where the cache will be stored (within the configured cache folder).
If `condition` is `false` caching will be skipped.
"""
function withcache end

"""
    purge()

Removes the cache data stored under the `key` key.
"""
function purge end

"""
    purgeall()

Removes all cached data.
"""
function purgeall end

### PRIVATE ###

"""
    cachekey(args...) :: String

Computes a unique cache key based on `args`. Used to generate unique `key`s for storing data in cache.
"""
function cachekey(args...) :: String
  key = IOBuffer()
  for a in args
    print(key, string(a))
  end

  bytes2hex(SHA.sha1(String(take!(key))))
end

end


text:
push!(LOAD_PATH,"../src/")

using Documenter, EconTools

makedocs(modules = [EconTools])

deploydocs(
    deps   = Deps.pip("mkdocs", "pygments","python-markdown-math"),
    repo = "github.com/pmarg/EconTools.jl.git",
    julia  = "0.6"
)


text:
function ricker(f₀::Float64, Δt::Float64)
    @assert f₀ > 0 "Peak frequency f₀ must be positive!"
    @assert Δt > 0 "Time interval Δt must be positive!"

    # signal length
    length = 2.0 / f₀ / Δt

    # number of points in each direction
    n = ceil(Int, length / 2)

    # time axis
    t = Δt * (-n:n)

    # amplitude
    x = (π * f₀ * t) .^ 2
    A = (1 .- 2x) .* exp.(-x)
end

function d2t(lyr::Dict{Symbol,Vector{Float64}})
    vₚ = lyr[:vp]
    d = lyr[:dh]
    t = cumsum(2d ./ vₚ)
    pushfirst!(t,0)
end

function d2t(lyr::Vector{Dict{Symbol,Vector{Float64}}})
    stacked = stack(lyr)
    d2t(stacked)
end

function seqcounts(labels::Vector)
    # find sequence boundaries
    indices = findall(diff(labels) .!= 0)

    # get sequence of unique labels
    sequence = labels[indices]

    # get label counts
    push!(sequence, labels[end])
    pushfirst!(indices, 0)
    push!(indices, length(labels))
    counts = diff(indices)

    sequence, Int64.(counts)
end


text:
# pseudo-definitions to show how everything behaves
#
# throw(label, val) = # throw a value to a dynamically enclosing block
#
# function rethrow(val)
#     global current_exception = val
#     throw(current_handler(), current_exception)
# end
#
# rethrow() = rethrow(current_exception)
#
# function throw(val)
#     global catch_backtrace = backtrace()
#     rethrow(val)
# end

## native julia error handling ##

error(e::Exception) = throw(e)
error{E<:Exception}(::Type{E}) = throw(E())
error(s::String) = throw(ErrorException(s))
error(s...)      = throw(ErrorException(string(s...)))

macro unexpected()
    :(error("unexpected branch reached"))
end

rethrow() = ccall(:jl_rethrow, Void, ())::None
rethrow(e) = ccall(:jl_rethrow_other, Void, (Any,), e)::None
backtrace() = ccall(:jl_backtrace_from_here, Array{Ptr{Void},1}, ())
catch_backtrace() = ccall(:jl_get_backtrace, Array{Ptr{Void},1}, ())

## system error handling ##

errno() = ccall(:jl_errno, Cint, ())
errno(e::Integer) = ccall(:jl_set_errno, Void, (Cint,), e)
strerror(e::Integer) = bytestring(ccall(:strerror, Ptr{Uint8}, (Int32,), e))
strerror() = strerror(errno())
systemerror(p, b::Bool) = b ? throw(SystemError(string(p))) : nothing

## assertion functions and macros ##

assert(x) = x ? nothing : error("assertion failed")
macro assert(ex,msgs...)
    msg = isempty(msgs) ? ex : msgs[1]
    if !isempty(msgs) && isa(msg, Expr)
        # message is an expression needing evaluating
        msg = :(string("assertion failed: ", $(esc(msg))))
    elseif isdefined(Base,:string)
        msg = string("assertion failed: ", msg)
    else
        # string() might not be defined during bootstrap
        msg = :(string("assertion failed: ", $(Expr(:quote,msg))))
    end
    :($(esc(ex)) ? $(nothing) : error($msg))
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
if VERSION >= v"0.6.0-"
    using SpecialFunctions: besseli
end

# This file contains different window functions.
# The function getWindow returns a pair of window functions based on a string

function getWindow(window::Symbol)
    if window == :gauss
        return window_gauss, window_gauss_hat
    elseif window == :spline
        return window_spline, window_spline_hat
    elseif window == :kaiser_bessel_rev
        return window_kaiser_bessel_rev, window_kaiser_bessel_rev_hat
    else # default to kaiser_bessel
        return window_kaiser_bessel, window_kaiser_bessel_hat
    end
end

function window_kaiser_bessel(x,n,m,sigma)
    b = pi*(2-1/sigma)
    arg = m^2-n^2*x^2
    if abs(x) < m/n
        y = sinh(b*sqrt(arg))/sqrt(arg)/pi
    elseif abs(x) > m/n
        y = zero(x)
    else
        y = b/pi
    end
    return y
end

function window_kaiser_bessel_hat(k,n,m,sigma)
    b = pi*(2-1/sigma)
    return besseli(0,m*sqrt(b^2-(2*pi*k/n)^2))
end

function window_kaiser_bessel_rev(x,n,m,sigma)
    b = pi*(2-1/sigma)
    if abs(x) < m/n
        arg = m*b*sqrt(1-(n*x/m)^2)
        y = 0.5/m*besseli(0,arg)
    else
        y = zero(x)
    end
    return y
end

function window_kaiser_bessel_rev_hat(k,n,m,sigma)
    b = pi*(2-1/sigma)

    arg = sqrt(complex((2*pi*m*k/n)^2-(m*b)^2))
    return sinc(arg/pi)
end

function window_gauss(x,n,m,sigma)
    b = m / pi
    if abs(x) < m/n
        y = 1 / sqrt(pi*b) * exp(-(n*x)^2 / b)
    else
        y =  zero(x)
    end
    return y
end

function window_gauss_hat(k,n,m,sigma)
    b = m / pi
    return exp(-(pi*k/(n))^2 * b)
end

function cbspline(m,x)
    if m == 1
        if x>=0 && x<1
            y = one(x)
        else
            y = zero(x)
        end
    else
        y = x/(m-1)*cbspline(m-1,x) + (m-x)/(m-1)*cbspline(m-1,x-1)
    end
    return y
end

function window_spline(x,n,m,sigma)
    if abs(x) < m/n
        y = cbspline(2*m, n*x+m)
    else
        y = zero(x)
    end
    return y
end

function window_spline_hat(k,n,m,sigma)
    return (sinc(k/n))^(2*m)
end


text:
using Test, BellScenario

@testset "./src/bell_inequalities.jl" begin

using SignalingDimension

@testset "k_guessing_game()" begin
    @testset "edge cases small Y" begin
        G_K = k_guessing_game(1,1,1)

        @test G_K == ones(Int64,1,1)
        @test G_K.β == 1

        G_K = k_guessing_game(1,1,0)

        @test G_K == zeros(Int64,1,1)
        @test G_K.β == 0
    end

    @testset "edge cases for k" begin
        G_K = k_guessing_game(5,2,5)

        @test G_K == ones(Int64, 5, 1)
        @test G_K.β == 1

        G_K = k_guessing_game(5,2,0)

        @test G_K == zeros(Int64, 5, 1)
        @test G_K.β == 0
    end

    @testset "trivial cases" begin
        G_K = k_guessing_game(4,2,2)

        @test G_K == [1 1 1 0 0 0;1 0 0 1 1 0;0 1 0 1 0 1;0 0 1 0 1 1;]
        @test G_K.β == 5
    end

    @testset "LocalSignaling" begin
        G_K = k_guessing_game(LocalSignaling(6,4,2),2)

        @test G_K == [1 1 1 0 0 0;1 0 0 1 1 0;0 1 0 1 0 1;0 0 1 0 1 1;]
        @test G_K.β == 5
    end

    @testset "errors" begin
        @test_throws DomainError k_guessing_game(4,2,-1)
        @test_throws DomainError k_guessing_game(4,2,5)
        @test_throws DomainError k_guessing_game(LocalSignaling(5,4,2),2)
    end
end

@testset "ambiguous_guessing_game()" begin
    @testset "X, Y, d parameters" begin
        @test ambiguous_guessing_game(3,4,2,2) == [2 0 0;0 2 0;1 1 1;1 1 1]
    end

    @testset "varying k" begin
        scenario = LocalSignaling(3,4,2)

        @test ambiguous_guessing_game(scenario, 0) == [1 1 1;1 1 1;1 1 1;1 1 1]
        @test ambiguous_guessing_game(scenario, 1) == [2 0 0;1 1 1;1 1 1;1 1 1]
        @test ambiguous_guessing_game(scenario, 2) == [2 0 0;0 2 0;1 1 1;1 1 1]
        @test ambiguous_guessing_game(scenario, 3) == [2 0 0;0 2 0;0 0 2;1 1 1]
        @test ambiguous_guessing_game(scenario, 4) == [2 0 0;0 2 0;0 0 2;0 0 2]
    end

    @testset "d edge cases" begin
        BG = ambiguous_guessing_game(LocalSignaling(4,6,4), 4)
        @test BG == [
            1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1;1 1 1 1;1 1 1 1
        ]
        @test BG.β == 4

        BG = ambiguous_guessing_game(LocalSignaling(6,4,4), 3)
        @test BG == [
            3 0 0 0 0 0;0 3 0 0 0 0;0 0 3 0 0 0;1 1 1 1 1 1
        ]
        @test BG.β == 12
    end

    @testset "errors" begin
        @test_throws DomainError ambiguous_guessing_game(LocalSignaling(3,3,4),2)
        @test_throws DomainError ambiguous_guessing_game(LocalSignaling(4,5,2),7)
    end
end

@testset "maximum_likelihood_facet()" begin
    @testset "trivial case" begin
        game = maximum_likelihood_facet(3,2)

        @test game isa BellGame
        @test game == [1 0 0;0 1 0;0 0 1]
        @test game.β == 2
    end

    @testset "example case" begin
        game = maximum_likelihood_facet(6,4)

        @test game isa BellGame
        @test game == [
            1 0 0 0 0 0;
            0 1 0 0 0 0;
            0 0 1 0 0 0;
            0 0 0 1 0 0;
            0 0 0 0 1 0;
            0 0 0 0 0 1
        ]
        @test game.β == 4
    end

    @testset "errors" begin
        @test_throws DomainError maximum_likelihood_facet(2,2)
        @test_throws DomainError maximum_likelihood_facet(4,5)
        @test_throws DomainError maximum_likelihood_facet(4,1)
    end
end

@testset "anti_guessing_facet()" begin
    @testset "trivial case" begin
        game = anti_guessing_facet(4,2,3)

        @test game isa BellGame
        @test game == [0 1 1 0;1 0 1 0;1 1 0 0;0 0 0 1]
        @test game.β == 3
    end

    @testset "example case" begin
        game = anti_guessing_facet(6,3,4)

        @test game isa BellGame
        @test game == [
            0 1 1 1 0 0;
            1 0 1 1 0 0;
            1 1 0 1 0 0;
            1 1 1 0 0 0;
            0 0 0 0 1 0;
            0 0 0 0 0 1
        ]
        @test game.β == 5
    end

    @testset "errors" begin
        @test_throws DomainError anti_guessing_facet(3,2,3)
        @test_throws DomainError anti_guessing_facet(4,3,3)
        @test_throws DomainError anti_guessing_facet(4,1,3)
        @test_throws DomainError anti_guessing_facet(7,4,2)
        @test_throws DomainError anti_guessing_facet(7,4,2)
        @test_throws DomainError anti_guessing_facet(7,4,5)
    end
end

@testset "ambiguous_guessing_facet()" begin
    @testset "trivial case" begin
        game = ambiguous_guessing_facet(4,2)

        @test game isa BellGame
        @test game == [2 0 0;0 2 0;0 0 2;1 1 1]
        @test game.β == 4
    end

    @testset "example cases" begin
        game1 = ambiguous_guessing_facet(6,2)

        @test game1 isa BellGame
        @test game1 == [
            4 0 0 0 0;
            0 4 0 0 0;
            0 0 4 0 0;
            0 0 0 4 0;
            0 0 0 0 4;
            1 1 1 1 1
        ]
        @test game1.β == 8

        game2 = ambiguous_guessing_facet(6,4)
        @test game2 isa BellGame
        @test game2 == [
            2 0 0 0 0;
            0 2 0 0 0;
            0 0 2 0 0;
            0 0 0 2 0;
            0 0 0 0 2;
            1 1 1 1 1
        ]
        @test game2.β == 8
    end

    @testset "errors" begin
        @test_throws DomainError ambiguous_guessing_facet(3,2)
        @test_throws DomainError ambiguous_guessing_facet(5,1)
        @test_throws DomainError ambiguous_guessing_facet(5,4)
    end
end

@testset "k_guessing_facet()" begin
    @testset "maximum_likelihood_facet generalization" begin
        k = 1

        for N in 3:6
            d = N - k
            gen_err_game = k_guessing_facet(N, d, k)

            succ_game = maximum_likelihood_facet(N, d)

            @test gen_err_game isa BellGame
            @test gen_err_game == succ_game
            @test gen_err_game.β == succ_game.β
        end
    end

    @testset "errors" begin
        @test_throws DomainError k_guessing_facet(2,1,1)
        @test_throws DomainError k_guessing_facet(5,2,5)
        @test_throws DomainError k_guessing_facet(5,2,0)
        @test_throws DomainError k_guessing_facet(5,4,2)
        @test_throws DomainError k_guessing_facet(5,1,2)
    end

    @testset "simple non-trivial cases" begin
        @testset "6-2-4 polytope" begin
            gen_err_game = k_guessing_facet(4, 2, 2)

            @test gen_err_game isa BellGame
            @test gen_err_game == [
                1 1 1 0 0 0;
                1 0 0 1 1 0;
                0 1 0 1 0 1;
                0 0 1 0 1 1
            ]
            @test gen_err_game.β == 5
        end

        @testset "10-d-5 polytope" begin
            gen_err_game1 = k_guessing_facet(5, 3, 2)

            @test gen_err_game1 isa BellGame
            @test gen_err_game1 == [
                1 1 1 1 0 0 0 0 0 0;
                1 0 0 0 1 1 1 0 0 0;
                0 1 0 0 1 0 0 1 1 0;
                0 0 1 0 0 1 0 1 0 1;
                0 0 0 1 0 0 1 0 1 1;
            ]
            @test gen_err_game1.β == 9

            gen_err_game2 = k_guessing_facet(5,2,3)

            @test gen_err_game2 isa BellGame
            @test gen_err_game2 == [
                1 1 1 1 1 1 0 0 0 0;
                1 1 1 0 0 0 1 1 1 0;
                1 0 0 1 1 0 1 1 0 1;
                0 1 0 1 0 1 1 0 1 1;
                0 0 1 0 1 1 0 1 1 1;
            ]
            @test gen_err_game2.β == 9
        end
    end
end

@testset "non_negativity_facet()" begin
    @testset "simple 3x3 case" begin
        game = non_negativity_facet(3,3)

        @test game isa BellGame
        @test game == [1 0 0;1 0 0;0 0 0]
        @test game.β == 1
    end

    @testset "simple 5x5 case" begin
        game = non_negativity_facet(5,5)

        @test game isa BellGame
        @test game == [1 0 0 0 0;1 0 0 0 0;1 0 0 0 0;1 0 0 0 0;0 0 0 0 0]
        @test game.β == 1
    end

    @testset "DomainErrors" begin
        @test_throws DomainError non_negativity_facet(1, 5)
        @test_throws DomainError non_negativity_facet(5, 1)
    end
end

@testset "coarse_grained_input_ambiguous_guessing_facet()" begin
    @testset "simple 4x4 case" begin
        game = coarse_grained_input_ambiguous_guessing_facet(4,2)

        @test game isa BellGame
        @test game == [2 0 0 0;0 2 0 0;0 0 1 1;1 1 1 0]
        @test game.β == 4
    end

    @testset "simple 6x6 case" begin
        game = coarse_grained_input_ambiguous_guessing_facet(6,2)

        @test game isa BellGame
        @test game == [4 0 0 0 0 0;0 4 0 0 0 0;0 0 4 0 0 0;0 0 0 4 0 0;0 0 0 0 1 3;1 1 1 1 1 0]
        @test game.β == 8
    end

    @testset "errors" begin
        @test_throws DomainError coarse_grained_input_ambiguous_guessing_facet(6,5)
        @test_throws DomainError coarse_grained_input_ambiguous_guessing_facet(3,1)
    end
end

end


text:
# Julia compiler integration

## cache

using Core.Compiler: CodeInstance, MethodInstance

struct CodeCache
    dict::Dict{MethodInstance,Vector{CodeInstance}}

    CodeCache() = new(Dict{MethodInstance,Vector{CodeInstance}}())
end

function Base.show(io::IO, ::MIME"text/plain", cc::CodeCache)
    print(io, "CodeCache with $(mapreduce(length, +, values(cc.dict); init=0)) entries")
    if !isempty(cc.dict)
        print(io, ": ")
        for (mi, cis) in cc.dict
            println(io)
            print(io, "  ")
            show(io, mi)

            function worldstr(min_world, max_world)
                if min_world == typemax(UInt)
                    "empty world range"
                elseif max_world == typemax(UInt)
                    "worlds $(Int(min_world))+"
                else
                    "worlds $(Int(min_world)) to $(Int(max_world))"
                end
            end

            for (i,ci) in enumerate(cis)
                println(io)
                print(io, "    CodeInstance for ", worldstr(ci.min_world, ci.max_world))
            end
        end
    end
end

Base.empty!(cc::CodeCache) = empty!(cc.dict)

const GLOBAL_CI_CACHE = CodeCache()

## method invalidations

function Core.Compiler.setindex!(cache::CodeCache, ci::CodeInstance, mi::MethodInstance)
    # make sure the invalidation callback is attached to the method instance
    callback(mi, max_world) = invalidate(cache, mi, max_world)
    if !isdefined(mi, :callbacks)
        mi.callbacks = Any[callback]
    else
        if all(cb -> cb !== callback, mi.callbacks)
            push!(mi.callbacks, callback)
        end
    end

    cis = get!(cache.dict, mi, CodeInstance[])
    push!(cis, ci)
end

# invalidation (like invalidate_method_instance, but for our cache)
function invalidate(cache::CodeCache, replaced::MethodInstance, max_world, depth=0)
    cis = get(cache.dict, replaced, nothing)
    if cis === nothing
        return
    end
    for ci in cis
        if ci.max_world == ~0 % Csize_t
            @assert ci.min_world - 1 <= max_world "attempting to set illogical constraints"
            ci.max_world = max_world
        end
        @assert ci.max_world <= max_world
    end

    # recurse to all backedges to update their valid range also
    if isdefined(replaced, :backedges)
        backedges = replaced.backedges
        # Don't touch/empty backedges `invalidate_method_instance` in C will do that later
        # replaced.backedges = Any[]

        for mi in backedges
            invalidate(cache, mi, max_world, depth + 1)
        end
    end
end

## method overrides

@static if isdefined(Base.Experimental, Symbol("@overlay"))

# use an overlay method table

Base.Experimental.@MethodTable(GLOBAL_METHOD_TABLE)

else

# use an overlay world -- a special world that contains all method overrides

const GLOBAL_METHOD_TABLE = nothing

const override_world = typemax(Csize_t) - 1

struct WorldOverlayMethodTable <: Core.Compiler.MethodTableView
    world::UInt
end

function Core.Compiler.findall(@nospecialize(sig::Type{<:Tuple}), table::WorldOverlayMethodTable; limit::Int=typemax(Int))
    _min_val = Ref{UInt}(typemin(UInt))
    _max_val = Ref{UInt}(typemax(UInt))
    _ambig = Ref{Int32}(0)
    ms = Base._methods_by_ftype(sig, limit, override_world, false, _min_val, _max_val, _ambig)
    if ms === false
        return Core.Compiler.missing
    elseif isempty(ms)
        # no override, so look in the regular world
        _min_val[] = typemin(UInt)
        _max_val[] = typemax(UInt)
        ms = Base._methods_by_ftype(sig, limit, table.world, false, _min_val, _max_val, _ambig)
    else
        # HACK: inference doesn't like our override world
        _min_val[] = table.world
    end
    if ms === false
        return Core.Compiler.missing
    end
    return Core.Compiler.MethodLookupResult(ms::Vector{Any}, Core.Compiler.WorldRange(_min_val[], _max_val[]), _ambig[] != 0)
end

end

"""
    @override mt def

!!! warning

    On Julia 1.6, evaluation of the expression returned by this macro should be postponed
    until run time (i.e. don't just call this macro or return its returned value, but
    save it in a global expression and `eval` it during `__init__`, additionally guarded
    by a check to `ccall(:jl_generating_output, Cint, ()) != 0`).
"""
macro override(mt, ex)
    if isdefined(Base.Experimental, Symbol("@overlay"))
        esc(quote
            Base.Experimental.@overlay $mt $ex
        end)
    else
        quote
            world_counter = cglobal(:jl_world_counter, Csize_t)
            regular_world = unsafe_load(world_counter)

            $(Expr(:tryfinally, # don't introduce scope
                quote
                    unsafe_store!(world_counter, $(override_world-1))
                    $(esc(ex))
                end,
                quote
                    unsafe_store!(world_counter, regular_world)
                end
            ))
        end
    end
end

## interpreter

using Core.Compiler: AbstractInterpreter, InferenceResult, InferenceParams, InferenceState, OptimizationParams

struct GPUInterpreter <: AbstractInterpreter
    global_cache::CodeCache
    method_table::Union{Nothing,Core.MethodTable}

    # Cache of inference results for this particular interpreter
    local_cache::Vector{InferenceResult}
    # The world age we're working inside of
    world::UInt

    # Parameters for inference and optimization
    inf_params::InferenceParams
    opt_params::OptimizationParams

    function GPUInterpreter(cache::CodeCache, mt::Union{Nothing,Core.MethodTable}, world::UInt)
        @assert world <= Base.get_world_counter()

        return new(
            cache,
            mt,

            # Initially empty cache
            Vector{InferenceResult}(),

            # world age counter
            world,

            # parameters for inference and optimization
            InferenceParams(unoptimize_throw_blocks=false),
            OptimizationParams(unoptimize_throw_blocks=false),
        )
    end
end

Core.Compiler.InferenceParams(interp::GPUInterpreter) = interp.inf_params
Core.Compiler.OptimizationParams(interp::GPUInterpreter) = interp.opt_params
Core.Compiler.get_world_counter(interp::GPUInterpreter) = interp.world
Core.Compiler.get_inference_cache(interp::GPUInterpreter) = interp.local_cache
Core.Compiler.code_cache(interp::GPUInterpreter) = WorldView(interp.global_cache, interp.world)

# No need to do any locking since we're not putting our results into the runtime cache
Core.Compiler.lock_mi_inference(interp::GPUInterpreter, mi::MethodInstance) = nothing
Core.Compiler.unlock_mi_inference(interp::GPUInterpreter, mi::MethodInstance) = nothing

function Core.Compiler.add_remark!(interp::GPUInterpreter, sv::InferenceState, msg)
    @safe_debug "Inference remark during GPU compilation of $(sv.linfo): $msg"
end

Core.Compiler.may_optimize(interp::GPUInterpreter) = true
Core.Compiler.may_compress(interp::GPUInterpreter) = true
Core.Compiler.may_discard_trees(interp::GPUInterpreter) = true
if VERSION >= v"1.7.0-DEV.577"
Core.Compiler.verbose_stmt_info(interp::GPUInterpreter) = false
end

if isdefined(Base.Experimental, Symbol("@overlay"))
Core.Compiler.method_table(interp::GPUInterpreter, sv::InferenceState) =
    Core.Compiler.OverlayMethodTable(interp.world, interp.method_table)
else
Core.Compiler.method_table(interp::GPUInterpreter, sv::InferenceState) =
    WorldOverlayMethodTable(interp.world)
end

## world view of the cache

using Core.Compiler: WorldView

function Core.Compiler.haskey(wvc::WorldView{CodeCache}, mi::MethodInstance)
    Core.Compiler.get(wvc, mi, nothing) !== nothing
end

function Core.Compiler.get(wvc::WorldView{CodeCache}, mi::MethodInstance, default)
    # check the cache
    for ci in get!(wvc.cache.dict, mi, CodeInstance[])
        if ci.min_world <= wvc.worlds.min_world && wvc.worlds.max_world <= ci.max_world
            # TODO: if (code && (code == jl_nothing || jl_ir_flag_inferred((jl_array_t*)code)))
            src = if ci.inferred isa Vector{UInt8}
                ccall(:jl_uncompress_ir, Any, (Any, Ptr{Cvoid}, Any),
                       mi.def, C_NULL, ci.inferred)
            else
                ci.inferred
            end
            return ci
        end
    end

    return default
end

function Core.Compiler.getindex(wvc::WorldView{CodeCache}, mi::MethodInstance)
    r = Core.Compiler.get(wvc, mi, nothing)
    r === nothing && throw(KeyError(mi))
    return r::CodeInstance
end

function Core.Compiler.setindex!(wvc::WorldView{CodeCache}, ci::CodeInstance, mi::MethodInstance)
    src = if ci.inferred isa Vector{UInt8}
        ccall(:jl_uncompress_ir, Any, (Any, Ptr{Cvoid}, Any),
                mi.def, C_NULL, ci.inferred)
    else
        ci.inferred
    end
    Core.Compiler.setindex!(wvc.cache, ci, mi)
end

## codegen/inference integration

function ci_cache_populate(cache, mt, mi, min_world, max_world)
    interp = GPUInterpreter(cache, mt, min_world)
    src = Core.Compiler.typeinf_ext_toplevel(interp, mi)

    # inference populates the cache, so we don't need to jl_get_method_inferred
    wvc = WorldView(cache, min_world, max_world)
    @assert Core.Compiler.haskey(wvc, mi)

    # if src is rettyp_const, the codeinfo won't cache ci.inferred
    # (because it is normally not supposed to be used ever again).
    # to avoid the need to re-infer, set that field here.
    ci = Core.Compiler.getindex(wvc, mi)
    if ci !== nothing && ci.inferred === nothing
        ci.inferred = src
    end

    return ci
end

function ci_cache_lookup(cache, mi, min_world, max_world)
    wvc = WorldView(cache, min_world, max_world)
    return Core.Compiler.get(wvc, mi, nothing)
end

## interface

function compile_method_instance(@nospecialize(job::CompilerJob),
                                 method_instance::MethodInstance, world)
    # populate the cache
    cache = ci_cache(job)
    mt = method_table(job)
    if ci_cache_lookup(cache, method_instance, world, typemax(Cint)) === nothing
        ci_cache_populate(cache, mt, method_instance, world, typemax(Cint))
    end

    # set-up the compiler interface
    debug_info_kind = if Base.JLOptions().debug_level == 0
        LLVM.API.LLVMDebugEmissionKindNoDebug
    elseif Base.JLOptions().debug_level == 1
        LLVM.API.LLVMDebugEmissionKindLineTablesOnly
    elseif Base.JLOptions().debug_level >= 2
        LLVM.API.LLVMDebugEmissionKindFullDebug
    end
    if job.target isa PTXCompilerTarget && !job.target.debuginfo
        debug_info_kind = LLVM.API.LLVMDebugEmissionKindNoDebug
    end
    lookup_fun = (mi, min_world, max_world) -> ci_cache_lookup(cache, mi, min_world, max_world)
    lookup_cb = @cfunction($lookup_fun, Any, (Any, UInt, UInt))
    params = Base.CodegenParams(;
        track_allocations  = false,
        code_coverage      = false,
        prefer_specsig     = true,
        gnu_pubnames       = false,
        debug_info_kind    = Cint(debug_info_kind),
        lookup             = Base.unsafe_convert(Ptr{Nothing}, lookup_cb))

    # generate IR
    GC.@preserve lookup_cb begin
        native_code = ccall(:jl_create_native, Ptr{Cvoid},
                            (Vector{MethodInstance}, Base.CodegenParams, Cint),
                            [method_instance], params, #=extern policy=# 1)
        @assert native_code != C_NULL
        llvm_mod_ref = ccall(:jl_get_llvm_module, LLVM.API.LLVMModuleRef,
                            (Ptr{Cvoid},), native_code)
        @assert llvm_mod_ref != C_NULL
        llvm_mod = LLVM.Module(llvm_mod_ref)
    end

    # get the top-level code
    code = ci_cache_lookup(cache, method_instance, world, typemax(Cint))

    # get the top-level function index
    llvm_func_idx = Ref{Int32}(-1)
    llvm_specfunc_idx = Ref{Int32}(-1)
    ccall(:jl_get_function_id, Nothing,
          (Ptr{Cvoid}, Any, Ptr{Int32}, Ptr{Int32}),
          native_code, code, llvm_func_idx, llvm_specfunc_idx)
    @assert llvm_func_idx[] != -1
    @assert llvm_specfunc_idx[] != -1

    # get the top-level function)
    llvm_func_ref = ccall(:jl_get_llvm_function, LLVM.API.LLVMValueRef,
                          (Ptr{Cvoid}, UInt32), native_code, llvm_func_idx[]-1)
    @assert llvm_func_ref != C_NULL
    llvm_func = LLVM.Function(llvm_func_ref)
    llvm_specfunc_ref = ccall(:jl_get_llvm_function, LLVM.API.LLVMValueRef,
                              (Ptr{Cvoid}, UInt32), native_code, llvm_specfunc_idx[]-1)
    @assert llvm_specfunc_ref != C_NULL
    llvm_specfunc = LLVM.Function(llvm_specfunc_ref)

    # configure the module
    triple!(llvm_mod, llvm_triple(job.target))
    if julia_datalayout(job.target) !== nothing
        datalayout!(llvm_mod, julia_datalayout(job.target))
    end

    return llvm_specfunc, llvm_mod
end


text:
export get_removed_quantifier_form

"""
    get_removed_quantifier_form(formula::Union{Formula, Term})

Return `formula` without any quantifiers.
"""
function get_removed_quantifier_form(formula::Union{Formula, Term})
    @match formula begin
        f::AFormula || f::EFormula => get_removed_quantifier_form(f.formula)
        f::Conjunction || f::Disjunction => typeof(f)(
            get_removed_quantifier_form(f.formula1),
            get_removed_quantifier_form(f.formula2)
        )
        f::Negation => Negation(
            get_removed_quantifier_form(f.formula)
        )
        f::AtomicFormula || f::Function => typeof(f)(
            f.name,
            get_removed_quantifier_form.(f.arguments)
        )
        f::Variable => f
    end
end


text:
# exports
export create, Cuboid, show

mutable struct Cuboid <: Shape
    bounds::Array{Float64, 2}
end

function Base.show(io::IO, x::Cuboid)
    println(io, "Cuboid")
    for i in 1:size(x.bounds)[1]
        println(io, "Dimension $i: Low $(x.bounds[i, 1]) High $(x.bounds[i, 2])")
    end
end

function create(c::Cuboid; resolution=nothing, rand_=0.0, type::Int64=1)
    if isa(resolution, Nothing)
        N = [10 for i in 1:size(c.bounds)[1]]
    else
        N = Int64.(round.((c.bounds[:, 2] - c.bounds[:, 1])/resolution))
    end
    lattice = (c.bounds[:, 2] - c.bounds[:, 1]) ./ N
    mesh = zeros(3, prod(N))
    a = 1
    R = CartesianIndices(Tuple(N))
    for I in R
        I_ = Tuple(I)
        for i in eachindex(I_)
            mesh[i, a] = c.bounds[i, 1] - lattice[i]/2 + (I_[i] + rand_*randn())*lattice[i]
        end
        a += 1
    end
    return Dict(
        :x => mesh,
        :v => zeros(size(mesh)),
        :y => copy(mesh),
        :volume => ones(prod(N))*prod(lattice),
        :type => type*ones(Int64, prod(N))
        )
end

max_stars_repo_name: sajidmannan/PDMesh

max_issues_repo_name: sajidmannan/PDMesh

max_forks_repo_name: sajidmannan/PDMesh


text:
using POMDPGallery
using Test
# using Plots

# try
#     pyplot()
# catch ex
#     warn("PyPlot is not working. Attempting to install it.")
#     ENV["PYTHON"]=""
#     Pkg.build("PyCall")
#     Pkg.build("PyPlot")
# end

# allow failure for pyplot problems :(
@test run_scripts(allow_failure=["ContinuumWorld", "LaserTag"])
@test gen_readme("/tmp/test_README.md")

@test read("/tmp/test_README.md", String) == read(joinpath(dirname(pathof(POMDPGallery)), "..", "README.md"), String)


text:
__precompile__()

module Codecs

export encode, decode, Base64, Zlib, BCD

abstract type Codec end

function encode(codec::Type{T}, s::AbstractString) where {T <: Codec}
    encode(codec, convert(Vector{UInt8}, codeunits(String(s))))
end

function decode(codec::Type{T}, s::AbstractString) where {T <: Codec}
    decode(codec, convert(Vector{UInt8}, codeunits(String(s))))
end

@static if isdefined(Base, Symbol("@gc_preserve"))
    import Base: @gc_preserve
else
    # This is not the correct definition of `@gc_preserve` on 0.5 and 0.6.
    # It should be no worse than the old code in this package though.
    macro gc_preserve(args...)
        esc(args[end])
    end
end

# RFC3548/RFC4648 base64 codec

abstract type Base64 <: Codec end

const b64enc_tbl = UInt8[
    'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H',
    'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P',
    'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X',
    'Y', 'Z', 'a', 'b', 'c', 'd', 'e', 'f',
    'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n',
    'o', 'p', 'q', 'r', 's', 't', 'u', 'v',
    'w', 'x', 'y', 'z', '0', '1', '2', '3',
    '4', '5', '6', '7', '8', '9', '+', '/']

const base64_pad = UInt8('=')
const sentinel = typemax(UInt8)

const b64dec_tbl = fill(sentinel, 256)
for (val, ch) in enumerate(b64enc_tbl)
    b64dec_tbl[ch] = val - 1
end

function base64dec(c::UInt8)
    v = b64dec_tbl[c]
    if v == sentinel
        error("Invalid base64 symbol: $(char(c))")
    end
    v
end

function encode(::Type{Base64}, input::Vector{UInt8})
    n = length(input)
    if n == 0
        return Vector{UInt8}(undef, 0)
    end

    m = Int(4 * ceil(n / 3))
    output = Vector{UInt8}(undef, m)

    k = 1
    for ii = 1:3:length(input)-2
        u, v, w = input[ii], input[ii+1], input[ii+2]
        output[k]     = b64enc_tbl[  (u >> 2)                     + 1]
        output[k + 1] = b64enc_tbl[(((u << 4) | (v >> 4)) & 0x3f) + 1]
        output[k + 2] = b64enc_tbl[(((v << 2) | (w >> 6)) & 0x3f) + 1]
        output[k + 3] = b64enc_tbl[  (w                   & 0x3f) + 1]
        k += 4
    end

    if n % 3 == 1
        output[end - 3] = b64enc_tbl[ (input[end] >> 2)         + 1]
        output[end - 2] = b64enc_tbl[((input[end] << 4) & 0x3f) + 1]
        output[end - 1] = base64_pad
        output[end]     = base64_pad
    elseif n % 3 == 2
        output[end - 3] = b64enc_tbl[  (input[end - 1] >> 2)      + 1]
        output[end - 2] = b64enc_tbl[(((input[end - 1] << 4) |
                                       (input[end] >> 4)) & 0x3f) + 1]
        output[end - 1] = b64enc_tbl[ ((input[end] << 2) & 0x3f)  + 1]
        output[end] = base64_pad
    end

    output
end

function decode(::Type{Base64}, input::Vector{UInt8})
    n = length(input)
    if n == 0
        return Vector{UInt8}(undef, 0)
    end

    if n % 4 != 0
        warn("Length of Base64 input is not a multiple of four.")
    end

    m = 3 * div(n,  4)
    if input[end] == base64_pad
        m -= 1
    end
    if input[end - 1] == base64_pad
        m -= 1
    end
    output = Vector{UInt8}(undef, m)

    k = 1
    @inbounds for ii = 1:4:length(input)-4
        # This loop is performance-critical, so inline base64dec and consolidate into a single branch point for error-checking
        u = b64dec_tbl[input[ii]]
        v = b64dec_tbl[input[ii + 1]]
        w = b64dec_tbl[input[ii + 2]]
        z = b64dec_tbl[input[ii + 3]]
        if u | v | w | z == sentinel
            if u == sentinel
                error("Invalid base64 symbol: $(char(input[ii]))")
            elseif v == sentinel
                error("Invalid base64 symbol: $(char(input[ii + 1]))")
            elseif w == sentinel
                error("Invalid base64 symbol: $(char(input[ii + 2]))")
            else
                error("Invalid base64 symbol: $(char(input[ii + 3]))")
            end
        end
        output[k]     = (u << 2) | (v >> 4)
        output[k + 1] = (v << 4) | (w >> 2)
        output[k + 2] = (w << 6) | z
        k += 3
    end

    u = input[end-3]
    v = input[end-2]
    w = input[end-1]
    z = input[end]
    u = base64dec(u)
    v = base64dec(v)

    if w == base64_pad && z == base64_pad
        output[end] = (u << 2) | (v >> 4)
    elseif z == base64_pad
        w = base64dec(w)
        output[end - 1] = (u << 2) | (v >> 4)
        output[end]     = (v << 4) | (w >> 2)
    else
        w = base64dec(w)
        z = base64dec(z)
        output[end - 2] = (u << 2) | (v >> 4)
        output[end - 1] = (v << 4) | (w >> 2)
        output[end]     = (w << 6) | z
    end

    output
end

# Zlib/Gzip

abstract type Zlib <: Codec end

const Z_NO_FLUSH      = 0
const Z_PARTIAL_FLUSH = 1
const Z_SYNC_FLUSH    = 2
const Z_FULL_FLUSH    = 3
const Z_FINISH        = 4
const Z_BLOCK         = 5
const Z_TREES         = 6

const Z_OK            = 0
const Z_STREAM_END    = 1
const Z_NEED_DICT     = 2
const ZERRNO          = -1
const Z_STREAM_ERROR  = -2
const Z_DATA_ERROR    = -3
const Z_MEM_ERROR     = -4
const Z_BUF_ERROR     = -5
const Z_VERSION_ERROR = -6

if Sys.isunix()
    const libz = "libz"
elseif Sys.iswindows()
    const libz = "zlib1"
end

# The zlib z_stream structure.
mutable struct z_stream
    next_in::Ptr{UInt8}
    avail_in::Cuint
    total_in::Culong

    next_out::Ptr{UInt8}
    avail_out::Cuint
    total_out::Culong

    msg::Ptr{UInt8}
    state::Ptr{Nothing}

    zalloc::Ptr{Nothing}
    zfree::Ptr{Nothing}
    opaque::Ptr{Nothing}

    data_type::Cint
    adler::Culong
    reserved::Culong

    function z_stream()
        strm = new()
        strm.next_in   = C_NULL
        strm.avail_in  = 0
        strm.total_in  = 0
        strm.next_out  = C_NULL
        strm.avail_out = 0
        strm.total_out = 0
        strm.msg       = C_NULL
        strm.state     = C_NULL
        strm.zalloc    = C_NULL
        strm.zfree     = C_NULL
        strm.opaque    = C_NULL
        strm.data_type = 0
        strm.adler     = 0
        strm.reserved  = 0
        strm
    end
end

function zlib_version()
    unsafe_string(ccall((:zlibVersion, libz), Ptr{UInt8}, ()))
end

function encode(::Type{Zlib}, input::Vector{UInt8}, level::Integer)
    if !(1 <= level <= 9)
        error("Invalid zlib compression level.")
    end

    strm = z_stream()
    ret = ccall((:deflateInit_, libz),
                Int32, (Ref{z_stream}, Int32, Ptr{UInt8}, Int32),
                strm, level, zlib_version(), sizeof(z_stream))

    if ret != Z_OK
        error("Error initializing zlib deflate stream.")
    end

    strm.next_in = pointer(input)
    strm.avail_in = length(input)
    strm.total_in = length(input)
    output = Vector{UInt8}(undef, 0)
    outbuf = Vector{UInt8}(undef, 1024)

    @gc_preserve input outbuf while ret != Z_STREAM_END
        strm.avail_out = length(outbuf)
        strm.next_out = pointer(outbuf)
        flush = strm.avail_in == 0 ? Z_FINISH : Z_NO_FLUSH
        ret = ccall((:deflate, libz),
                    Int32, (Ref{z_stream}, Int32),
                    strm, flush)
        if ret != Z_OK && ret != Z_STREAM_END
            error("Error in zlib deflate stream ($(ret)).")
        end

        if length(outbuf) - strm.avail_out > 0
            append!(output, outbuf[1:(length(outbuf) - strm.avail_out)])
        end
    end

    ret = ccall((:deflateEnd, libz), Int32, (Ref{z_stream},), strm)
    if ret == Z_STREAM_ERROR
        error("Error: zlib deflate stream was prematurely freed.")
    end

    output
end

function encode(::Type{Zlib}, input::AbstractString, level::Integer)
    encode(Zlib, convert(Vector{UInt8}, codeunits(String(input))), level)
end

encode(::Type{Zlib}, input::Vector{UInt8}) = encode(Zlib, input, 9)

function decode(::Type{Zlib}, input::Vector{UInt8})
    strm = z_stream()
    ret = ccall((:inflateInit_, libz),
                Int32, (Ref{z_stream}, Ptr{UInt8}, Int32),
                strm, zlib_version(), sizeof(z_stream))

    if ret != Z_OK
        error("Error initializing zlib inflate stream.")
    end

    strm.next_in = pointer(input)
    strm.avail_in = length(input)
    strm.total_in = length(input)
    output = Vector{UInt8}(undef, 0)
    outbuf = Vector{UInt8}(undef, 1024)

    @gc_preserve input outbuf while ret != Z_STREAM_END
        strm.next_out = pointer(outbuf)
        strm.avail_out = length(outbuf)
        ret = ccall((:inflate, libz),
                    Int32, (Ref{z_stream}, Int32),
                    strm, Z_NO_FLUSH)
        if ret == Z_DATA_ERROR
            error("Error: input is not zlib compressed data.")
        elseif ret != Z_OK && ret != Z_STREAM_END
            error("Error in zlib inflate stream ($(ret)).")
        end

        if length(outbuf) - strm.avail_out > 0
            append!(output, outbuf[1:(length(outbuf) - strm.avail_out)])
        end
    end

    ret = ccall((:inflateEnd, libz), Int32, (Ref{z_stream},), strm)
    if ret == Z_STREAM_ERROR
        error("Error: zlib inflate stream was prematurely freed.")
    end

    output
end

# Packed binary-coded decimal (BCD) integers
# Two decimal digits per byte, each represented with 4 bits

abstract type BCD <: Codec end

function encode(::Type{BCD}, i::Integer, bigendian::Bool = false)
    if i < 0
        error("Integer must not be negative")
    end
    ndig = ndigits(i)
    ndig += isodd(ndig)
    v = digits(i; base=10, pad=ndig)
    nbytes = trunc(Integer, ndig/2)
    out = Vector{UInt8}(undef, nbytes)
    dj = 1-2*bigendian
    for i = 1:nbytes
        j = bigendian*length(v) + dj*2*(i-1) + 1-bigendian
        out[i] = (v[j]<<4) | v[j+dj]
    end
    return out
end

function decode(::Type{BCD}, v::Vector{UInt8}, bigendian::Bool = false)
    ret = 0
    if bigendian
        for i = 1:length(v)
            ret = 100*ret + 10 * Int(v[i]>>>4) + Int(v[i]&0x0f)
        end
    else
        for i = length(v):-1:1
            ret = 100*ret + Int(v[i]>>>4) + 10 * Int(v[i]&0x0f)
        end
    end
    return ret
end

end # module Codecs


text:
using MiniFB

export Window
export draw!, fillcolor!

import Base: show, close

"""
    Window
Floating window which we can render pixels into

```julia
w = Window("Hello World", 400, 500)

fillcolor!(w, 255, 0, 0)
draw!(w, Rect(10, 10, 200, 120))

renderloop(w)
```
"""
mutable struct Window
    width::Int
    height::Int
    buffer::Vector{UInt32}
    wndptr::Ptr{Cvoid}
    fillcolor::UInt32
end

"""
    Window(title, width, height)
Create a window on the screen of given proportions which we can draw in.
"""
function Window(title::AbstractString, width::Integer, height::Integer)
    wndptr = mfb_open_ex(title, width, height, MiniFB.WF_RESIZABLE)
    Window(width, height,
        zeros(UInt32, width*height),
        wndptr,
        mfb_rgb(255, 255, 255))
end

function Window(width, height)
    Window("Untitled", width, height)
end

function close(w::Window)
    mfb_close(w.wndptr)
end

function show(io::IO, wnd::Window)
    print(io, "Window($(wnd.width), $(wnd.height))")
end

function draw!(wnd::Window, r::Rect)
    x, y = r.x, r.y
    w, h = r.width, r.height

    for j in y:(y + h)
        x0 = x + j*wnd.width
        for i in x0:(x0 + w)
            wnd.buffer[i] = wnd.fillcolor
        end
    end
end

"""
    fillcolor!(w::Window, red, green, blue)
Specify color to be used by drawing commands
"""
function fillcolor!(w::Window, red, green, blue)
    w.fillcolor = mfb_rgb(red, green, blue)
end

quit = false


text:
using Test
using LogicCircuits
using ProbabilisticCircuits

# This tests are supposed to test queries on the circuits
@testset "Probability of Full Evidence" begin
    # Uses a PSDD with 4 variables, and tests 3 of the configurations to
    # match with python. Also tests all probabilities sum up to 1.

    EPS = 1e-7;
    prob_circuit = zoo_psdd("little_4var.psdd");
    @test prob_circuit isa Vector{<:ProbΔNode};

    flow_circuit = FlowΔ(prob_circuit, 16, Bool)
    @test flow_circuit isa Vector{<:FlowΔNode};

    # Step 1. Check Probabilities for 3 samples
    data = XData(Bool.([0 0 0 0; 0 1 1 0; 0 0 1 1]));
    true_prob = [0.07; 0.03; 0.13999999999999999]

    calc_prob = log_likelihood_per_instance(flow_circuit, data)
    calc_prob = exp.(calc_prob)

    for i = 1:3
        @test true_prob[i] ≈ calc_prob[i] atol= EPS;
    end

    # Step 2. Add up all probabilities and see if they add up to one
    N = 4;
    data_all = XData(generate_data_all(N))

    calc_prob_all = log_likelihood_per_instance(flow_circuit, data_all)
    calc_prob_all = exp.(calc_prob_all)
    sum_prob_all = sum(calc_prob_all)

    @test 1 ≈ sum_prob_all atol = EPS;
end

@testset "Probability of partial Evidence (marginals)" begin
    EPS = 1e-7;
    prob_circuit = zoo_psdd("little_4var.psdd");

    data = XData(
        Int8.([0 0 0 0; 0 1 1 0; 0 0 1 1;
                0 0 0 -1; -1 1 0 -1; -1 -1 -1 -1; 0 -1 -1 -1])
    );
    true_prob = [0.07; 0.03; 0.13999999999999999;
                    0.3499999999999; 0.1; 1.0; 0.8]

    opts = (compact⋀=false, compact⋁=false)
    flow_circuit = UpFlowΔ(prob_circuit, 16, Float64, opts)
    calc_prob = marginal_log_likelihood_per_instance(flow_circuit, data)
    calc_prob = exp.(calc_prob)

    for i = 1:length(true_prob)
        @test true_prob[i] ≈ calc_prob[i] atol= EPS;
    end

    # Now trying the other api without instantiating a flow circuit
    fc2, calc_prob2 = marginal_log_likelihood_per_instance(prob_circuit, data)
    calc_prob2 = exp.(calc_prob2)
    for i = 1:length(true_prob)
        @test true_prob[i] ≈ calc_prob2[i] atol= EPS;
    end

end

@testset "Marginal Pass Down" begin
    EPS = 1e-7;
    prob_circuit = zoo_psdd("little_4var.psdd");

    N = 4
    data_full = XData(Int8.(generate_data_all(N)))
    opts= (compact⋀=false, compact⋁=false)

    flow_circuit   = FlowΔ(prob_circuit, 16, Float64, opts)
    flow_circuit_marg = FlowΔ(prob_circuit, 16, Float64, opts)

    # Comparing with down pass with fully obeserved data
    pass_up_down(flow_circuit, data_full)
    marginal_pass_up_down(flow_circuit_marg, data_full)

    for (ind, node) in enumerate(flow_circuit)
        if node isa HasDownFlow
            @test all(  isapprox.(downflow(flow_circuit[ind]), downflow(flow_circuit_marg[ind]), atol = EPS) )
        end
    end

    # Validating one example with missing features done by hand
    data_partial = XData(Int8.([-1 1 -1 1]))
    flow_circuit_part  = FlowΔ(prob_circuit, 16, Float64, opts)
    ProbabilisticCircuits.marginal_pass_up_down(flow_circuit_part, data_partial)

    # (node index, correct down_flow_value)
    true_vals = [(1, 0.5),
                (2, 1.0),
                (3, 1/3),
                (4, 1.0),
                (5, 0.5),
                (6, 0.0),
                (7, 2/3),
                (8, 0.0),
                (9, 0.3333333333333),
                (10, 0.0),
                (11, 0.6666666666666),
                (12, 0.0),
                (13, 1.0),
                (14, 0.5),
                (15, 0.0),
                (16, 0.5),
                (17, 0.0),
                (18, 1.0),
                (19, 1.0),
                (20, 1.0)]

    for ind_val in true_vals
        @test downflow(flow_circuit_part[ind_val[1]])[1] ≈ ind_val[2] atol= EPS
    end

end

function test_mpe_brute_force(prob_circuit, evidence)
    EPS = 1e-9;
    result = MPE(prob_circuit, evidence);
    for idx = 1 : num_examples(evidence)
        marg = XData(generate_all(evidence.x[idx,:]));
        fc, lls = log_likelihood_per_instance(prob_circuit, marg);
        brute_mpe = marg.x[argmax(lls), :]

        # Compare and validate p(result[idx]) == p(brute_mpe)
        comp_data = XData(vcat(result[idx,:]',  brute_mpe'))
        fc2, lls2 = log_likelihood_per_instance(prob_circuit, comp_data);

        @test lls2[1] ≈ lls2[2] atol= EPS
    end
end

@testset "MPE Brute Force Test Small (4 var)" begin
    prob_circuit = zoo_psdd("little_4var.psdd");
    evidence = XData( Int8.( [-1 0 0 0;
                                0 -1 -1 0;
                                1 1 1 -1;
                                1 0 1 0;
                                -1 -1 -1 1;
                                -1 -1 -1 -1] ))

    test_mpe_brute_force(prob_circuit, evidence)

end

@testset "MPE Brute Force Test Big (15 var)" begin
    N = 15
    COUNT = 10

    prob_circuit = zoo_psdd("exp-D15-N1000-C4.psdd");
    evidence = XData(Int8.(rand( (-1,0,1), (COUNT, N) )))

    test_mpe_brute_force(prob_circuit, evidence)
end


text:
## Connection and Linearization

## A connection between polynomial systems `P` and `Q` is a function `C(m,n)` satisfying
## p_n = ∑C(m,n) q_m or  pᵢ = ∑ Cʲᵢqⱼ (if unicode is helpful)
##
## From this a polynomial in P = ∑ aᵢ⋅pᵢ can be represented in a Q through
##
##  ∑ᵢ aᵢ⋅pᵢ = ∑ᵢ aᵢ  ∑ⱼ Cʲᵢqⱼ  = ∑ᵢ (∑ⱼ (aᵢ⋅Cʲᵢ)) qⱼ or ∑ⱼ (∑ᵢ aᵢCʲᵢ) qⱼ
##
## The choice of which depends on how the connection coefficients are presented.
##
## A connection between P and Q defines a `convert` method `Base.convert(Q,  p::P)`
## Conversion can also be achieved through polynomial evaluation. For example
## if `x=variable(Q)` then `Base.convert(Q,p) = p(x)` will also convert `p` into
## the basis Q. This is helpful when Q is a standard basis polynomial.
##
## However, if Q is not a standard basis polynomial, it isn't helpful unless
## one can multiply monomomials within Q. For standard basis polynomials this
## is known, as  xⁱ⋅xʲ = xⁱ⁺ʲ.
##
## For general Q though, one needs a linearization: a function α with
##
##  qᵢqⱼ = ∑_k α(k,i,j) q_k
##
## with a given α, multiplication can be defined. The linearization can be defined
## through composition: If P is the standard basis, connect Q to P, multiply ini P, then connect
## to Q.
##
## In short, if there is a linearization then multiplication can be defined and conversion defined through
## polynomial evaluation
## If there is a connection defined then conversion is defined and multiplication can be defined through composition
##

## For `AbstractCCOP` polynomials, Koepf and Schmersau present recursive formulas for the connection coefficients
## between Q and the standard basis, the standard and Q, and between P and Q when σ and σ̂ are the same. (This is the case when a polynomial system is parameterized, like Laguerre, Gegenbauer, Bessel, and  Jacobi.)

ConvertibleTypes = Union{AbstractCOP, Polynomials.StandardBasisPolynomial,  FallingFactorial}
kn(::Type{P},  n) where{P <: Union{Polynomials.StandardBasisPolynomial,  FallingFactorial}} = one(eltype(one(P)))
k1k0(::Type{P},  n) where{P <: Union{Polynomials.StandardBasisPolynomial,  FallingFactorial}} = one(eltype(one(P)))
knk_1(::Type{P},  n) where{P <: Union{Polynomials.StandardBasisPolynomial,  FallingFactorial}} = one(eltype(one(P)))

## If  the connection  coefficients for  (P,Q) satisfy c₀⋅ C̃ⁱⱼ + c₁⋅C̃ⁱ⁺¹ⱼ  + c₂⋅*C̃ⁱ⁺²ⱼ = 0 for some computable
##  values c₀,c₁,c₂ then this  will implement  `convert(Q,p::P)`
function _convert_cop(::Type{Q}, p::P) where {P <: ConvertibleTypes,
                                               Q <: ConvertibleTypes}

    d = degree(p)
    T,S = eltype(Q), eltype(p)  #
    R = typeof(one(T) * p[0]/1)
    as = zeros(R, 1+d)

    λⱼⱼ = one(R)  * k0(P) / k0(Q)

    for j in  0:d

        λ = λⱼⱼ

        λⱼⱼ *= k1k0(P,j) / k1k0(Q,j)

        pⱼ = p[j]
        iszero(pⱼ) && continue

        C̃ʲⱼ = one(R)

        as[1+j] += pⱼ * (λ * C̃ʲⱼ)

        C̃ⁱ⁺²ⱼ, C̃ⁱ⁺¹ⱼ = zero(R), C̃ʲⱼ  # for recursion starting  with C̃ʲ⁺¹ⱼ, C̃ʲⱼ
        for i in j-1:-1:0

            λ *= k1k0(Q,i)

            # c₀⋅ C̃ⁱⱼ + c₁⋅C̃ⁱ⁺¹ⱼ  + c₂⋅*C̃ⁱ⁺²ⱼ) = 0
            c₀,c₁,c₂ = connection_m(P, Q, i, j)
            C̃ⁱⱼ = iszero(c₀) ? zero(R) : -one(R)*(c₁*C̃ⁱ⁺¹ⱼ + c₂*C̃ⁱ⁺²ⱼ) / c₀
            as[1+i] += pⱼ * (λ * C̃ⁱⱼ)

            C̃ⁱ⁺²ⱼ, C̃ⁱ⁺¹ⱼ = C̃ⁱ⁺¹ⱼ, C̃ⁱⱼ
        end
    end
    X = Polynomials.indeterminate(Q,p)
    ⟒(Q){R,X}(as)
end

# ## Use FastTransform for conversion, as  possible
# ## FastTransforms.kind2string.(0:15)
# # Legendre<--Chebyshev
# function _convert(::Type{Q}, p::P) where {
#     T <: AbstractFloat,
#     Q <: Union{Legendre, OrthonormalLegendre},
#     P <: Union{Chebyshev{T}, OrthonormalChebyshev{T}}}
#     ps = coeffs(p)
#     Q(cheb2leg(ps, normcheb=isorthonormal(P), normleg=isorthonormal(Q)))
# end

# # Chebyshev<--Legendre
# function _convert(::Type{Q}, p::P) where {
#     T <: AbstractFloat,
#     Q <: Union{Chebyshev, OrthonormalChebyshev},
#     P <: Union{Legendre{T}, OrthonormalLegendre{T}}
# }
#     ps =  coeffs(p)
#     Q(leg2cheb(ps, normleg=isorthonormal(P), normcheb=isorthonormal(Q)) )
# end

# # ultraspherical<--ultraspherical
# function _convert(::Type{Q}, p::P) where {
#     β, α,  T <: AbstractFloat,
#     Q <: Union{Gegenbauer{β}, OrthonormalGegenbauer{β}},
#     P <: Union{Gegenbauer{α, T}, OrthonormalGegenbauer{α, T}}
# }

#     ps =  coeffs(p)
#     Q( ultra2ultra(ps, α, β, norm1=isorthonormal(P), norm2=isorthonormal(Q)) )

# end

# # Jacobi<--Jacobi
# function _convert(::Type{Q}, p::P) where {
#     γ, δ, α, β,  T <: AbstractFloat,
#     Q <: Union{Jacobi{γ, δ}, OrthonormalJacobi{γ, δ}},
#     P <: Union{Jacobi{α, β, T}, OrthonormalJacobi{α, β, T}}
# }

#     ps =  coeffs(p)
#     Q( jac2jac(ps, α, β, γ, δ,  norm1=isorthonormal(P), norm2=isorthonormal(Q)) )

# end

# # Laguerre<--Laguerre
# function _convert(::Type{Q}, p::P) where {
#     α, β,  T <: AbstractFloat,
#     Q <: Union{Laguerre{β}, OrthonormalLaguerre{β}},
#     P <: Union{Laguerre{α, T}, OrthonormalLaguerre{α, T}}
# }

#     ps =  coeffs(p)
#     Q( lag2lag(ps, α, β,  norm1=isorthonormal(P), norm2=isorthonormal(Q)) )

# end

# # Jacobi<--ultraspherical
# function _convert(::Type{Q}, p::P) where {
#     γ, δ, α,  T <: AbstractFloat,
#     Q <: Union{Jacobi{γ, δ}, OrthonormalJacobi{γ, δ}},
#     P <: Union{Gegenbauer{α, T}, OrthonormalGegenbauer{α, T}}
# }

#     ps =  coeffs(p)
#     Q( ultra2jac(ps, α, γ, δ,  normultra=isorthonormal(P), normjac=isorthonormal(Q)) )

# end

# # ultraspherical<--Jacobi
# function _convert(::Type{Q}, p::P) where {
#     α, γ, δ,  T <: AbstractFloat,
#     Q <: Union{Gegenbauer{α, T}, OrthonormalGegenbauer{α, T}},
#     P <: Union{Jacobi{γ, δ}, OrthonormalJacobi{γ, δ}}
# }

#     ps =  coeffs(p)
#     Q( jac2ultra(ps,γ, δ, α, normjac=isorthonormal(P), normultra=isorthonormal(Q)) )

# end

# # Jacobi<--Chebyshev
# function _convert(::Type{Q}, p::P) where {
#     γ, δ,  T <: AbstractFloat,
#     Q <: Union{Jacobi{γ, δ}, OrthonormalJacobi{γ, δ}},
#     P <: Union{Chebyshev{T}, OrthonormalChebyshev{T}}
# }

#     ps =  coeffs(p)
#     Q( cheb2jac(ps,γ, δ, normcheb=isorthonormal(P), normjac=isorthonormal(Q)) )

# end

# # Chebyshev<--Jacobi
# function _convert(::Type{Q}, p::P) where {
#     γ, δ,  T <: AbstractFloat,
#     Q <: Union{Chebyshev, OrthonormalChebyshev},
#     P <: Union{Jacobi{γ, δ,T}, OrthonormalJacobi{γ, δ,T}}
# }

#     ps =  coeffs(p)
#     Q( jac2cheb(ps,γ, δ, normjac = isorthonormal(P), normcheb=isorthonormal(Q)) )

# end

# # ultraspherical<--Chebyshev
# function _convert(::Type{Q}, p::P) where {
#     α,  T <: AbstractFloat,
#     Q <: Union{Gegenbauer{α}, OrthonormalGegenbauer{α}},
#     P <: Union{Chebyshev{T}, OrthonormalChebyshev{T}}
# }

#     ps =  coeffs(p)
#     Q( cheb2ultra(ps,α, normcheb=isorthonormal(P), normultra=isorthonormal(Q)) )

# end

# # Chebyshev<--ultraspherical
# function _convert(::Type{Q}, p::P) where {
#     α,  T <: AbstractFloat,
#     Q <: Union{Chebyshev, OrthonormalChebyshev},
#     P <: Union{Gegenbauer{α,T}, OrthonormalGegenbauer{α,T}}
# }
#     ps =  coeffs(p)
#     Q( ultra2cheb(ps,α, normultra=isorthonormal(P), normcheb=isorthonormal(Q)) )

# end

##  Koepf and Schmersa Thm 2, p11 connection for P to Q when σ = σ̂
function connection_m(::Type{P},::Type{Q},m,n) where {
    P <: AbstractCCOP,
    Q <: AbstractCCOP}

    a,b,c,d,e = abcde(P)
    ā,b̄,c̄,d̄,ē = abcde(Q)

    (a == ā && b == b̄ && c == c̄) || throw(ArgumentError("This connection assume σ = σ̄"))
    a²,b²,c²,d²,d̄²,e²,ē²,m²,n² = a^2, b^2, c^2, d^2, d̄^2, e^2, ē^2, m^2, n^2

    c0 = -(m-n) ⋅ (a⋅m + d - a + a⋅n) ⋅ (d̄ + 2a⋅m) ⋅ (d̄ + a + 2a⋅m) ⋅ (d̄ + 3a + 2a⋅m)
    c0 *= (d̄ + 2a⋅m  + 2a)^2

    c1 = -d⋅b⋅n⋅d̄ + 2d⋅a⋅m²⋅b + d⋅b⋅d̄ + 2d⋅a⋅m⋅b + 2d⋅ē⋅n⋅a
    c1 += d⋅d̄⋅ē + 2d⋅d̄⋅b⋅m - m⋅b⋅d̄² - e⋅d̄² - 4a²⋅m²⋅e - m²⋅a⋅b⋅d̄ + b⋅n⋅d̄⋅a - 2e⋅d̄⋅a
    c1 += -4a²⋅m⋅e - 4e⋅d̄⋅a⋅m + 2m²⋅a²⋅ē + 2ē⋅a²⋅n² - 2ē⋅a²⋅n - m⋅a⋅b⋅d̄ + 2m⋅d̄⋅ē⋅a
    c1 += 2m⋅ē⋅a² - b⋅n²⋅d̄⋅a
    c1 *= (d̄ + 2a⋅m + 2a) ⋅ (m + 1) ⋅  (d̄ + a + 2a⋅m) ⋅ (d̄ + 3a + 2a⋅m)

    c2 =  -(d̄ + 2a⋅m) ⋅ (m + 1) ⋅ (-a⋅m - 2a + a⋅n - d̄ + d) ⋅ (a⋅m + a⋅n + a + d̄)
    c2a = a⋅b²⋅m² - 4a²⋅m²⋅c - 8a²⋅m⋅c + 2a⋅m⋅b² - 4a⋅d̄⋅m⋅c + m⋅b²⋅d̄ - 4a⋅d̄⋅c - a⋅ē² + a⋅b² - c⋅d̄²
    c2a += b⋅ē⋅d̄ - 4⋅a²⋅c + b²⋅d̄
    c2 *= c2a ⋅ (m + 2)

    return (c0, c1, c2)

end

##  Koepf and Schmersa Thm 4, p14 connection for P, Q=x^n
function connection_m(::Type{P},::Type{Q},m,n) where {
    P <: AbstractCCOP,
    Q <: Polynomials.StandardBasisPolynomial}

    a,b,c,d,e = abcde(P)

    c0 = (m-n) ⋅ (a⋅n + d - a + a⋅m)

    c1 = (m + 1) ⋅ (b⋅m + e)

    c2 = c ⋅ (m + 1) ⋅ (m + 2)

    return (c0, c1, c2)

end

##  Koepf and Schmersa Thm 5, p16 connection for P=x^n, Q
function connection_m(::Type{P},::Type{Q},m,n) where {
    P <: Polynomials.StandardBasisPolynomial,
    Q <: AbstractCCOP}

    ā,b̄,c̄,d̄,ē = abcde(Q)
    ā²,b̄²,c̄²,d̄²,ē²,m²,n² = ā^2, b̄^2, c̄^2, d̄^2, ē^2, m^2, n^2

    c0 =  (n - m) ⋅ (d̄ + 2ā⋅m) ⋅ (d̄ + 3ā + 2ā⋅m) ⋅ (d̄ + ā + 2ā⋅m) ⋅ (d̄ + 2ā⋅m + 2ā)^2

    c1 =  (d̄⋅ē + b̄⋅d̄ +2d̄⋅b̄⋅m + 2ā⋅m²⋅b̄ + 2ā⋅m⋅b̄ +2ē⋅ā⋅n - d̄⋅b̄⋅n) ⋅ (d̄ + 2ā⋅m + 2ā)
    c1 *= (m+1) ⋅ (d̄ + 3ā + 2ā⋅m) ⋅  (d̄ + ā + 2ā⋅m)

    c2 = -(m+2)
    c2a = -4ā²⋅c̄⋅m²
    c2a +=  ā⋅b̄²⋅m² + 2ā⋅b̄²⋅m - 4ā⋅c̄⋅m⋅d̄ - 8ā²⋅c̄⋅m + m⋅b̄²⋅d̄ - ā⋅ē^2 - d̄²⋅c̄ + b̄⋅ē⋅d̄ -4ā²⋅c̄
    c2a += -4ā⋅c̄⋅d̄ + ā⋅b̄² +b̄²⋅d̄
    c2 *= c2a ⋅ (ā⋅m + ā⋅n + ā + d̄) ⋅ (m + 1) ⋅ (d̄ + 2ā⋅m)

    return (c0, c1, c2)

end

## With an explicit connectioni
function connection(::Type{P}, q::Q) where
    {P <: Polynomials.AbstractPolynomial,
     Q<:Polynomials.AbstractPolynomial}

    n = degree(q)
    T = eltype(q)
    T = eltype(one(eltype(one(P)))*one(q)/1)
    cs = zeros(T, n+1)

    for k in 0:n
        for (i,val) in Connection{P,Q}(n, k)
            cs[1+k] = muladd(q[i], val, cs[1+k])
        end
    end
    X = Polynomials.indeterminate(P,q)
    ⟒(P){T,X}(cs)

end


text:
export husimi

using LinearAlgebra # for ⋅ (cdot)

struct Zero end
Base.:+(x, ::Zero) = x
Base.convert(::Type{Float64}, ::Zero) = 0.0
Base.convert(::Type{Any}, ::HusimiFunctions.Zero) = 0.0

const VSV{D} = Vector{SVector{D, Float64}}

"""
    husimi(ψ, r::Vector{<:SVector}, r₀::SVector, k₀::SVector, σ[, B]; thread)
Calculate the Husimi function of `ψ` at `(r₀, k₀)` with precision `σ`.
`ψ` is a vector whose values `ψ[i]` is the wavefunction value
at position `r[i]` (`r` is a vector of `SVector`s).
The keyword `thread` means to parallelize the operation over processor threads
(`true` by default).

Optionally calculate in a non-zero magnetic field `B`.
The units of `r, σ` are expected in nm, `k` in 1/nm and B in Tesla.

Notice that the magnetic transformation assumes Landau gauge in x orientation,
see https://arxiv.org/abs/1912.04622 for defining other gauges (by changing
the function `magneticphase`).
"""
function husimi(ψ, r, r₀::SVector, k₀::SVector, σ, B = Zero();
                thread = true)
    if thread
        abs2(coherent_projection_threaded(ψ, r, r₀, k₀, σ, B))
    else
        abs2(coherent_projection(ψ, r, r₀, k₀, σ, B))
    end
end

##############################################################
# Projection
##############################################################
# TODO: Check if dot is as optimized as  r[1]*r[1] + r[2]*r[2]
function coherent_projection_threaded(ψ, r, r₀, k₀, σ, B = Zero())
    D = length(r₀)
    sums = zeros(ComplexF64, Threads.nthreads())
    skip = (3*σ)^2
    denominator = (4*σ^2)

    @inbounds @fastmath Threads.@threads for i in 1:length(r)
        δr = r[i] - r₀
        rsq = δr⋅δr
        rsq > skip && continue
        space = exp(-rsq/denominator)
        phase = k₀⋅δr
        mphase = magneticphase(r[i], r₀, B)
        ss, cs = sincos(phase + mphase)
        sums[Threads.threadid()] += conj(ψ[i]) * space * complex(cs, ss)
    end
    return sum(sums)/(σ*sqrt(2π))^(D/2)
end

function coherent_projection(ψ, r, r₀, k₀, σ, B = Zero())
    D = length(r₀)
    sums = zero(ComplexF64)
    skip = (3*σ)^2
    denominator = (4*σ^2)

    @inbounds @fastmath for i in 1:length(r)
        δr = r[i] - r₀
        rsq = δr⋅δr
        rsq > skip && continue
        space = exp(-rsq/denominator)
        phase = k₀⋅δr
        mphase = magneticphase(r[i], r₀, B)
        ss, cs = sincos(phase + mphase)
        sums += conj(ψ[i]) * space * complex(cs, ss)
    end
    return sums/(σ*sqrt(2π))^(D/2)
end

const prefac = (10.0^(-2.0))/6.582 # e/ħ prefactor, when multiplying Tesla gives 1/nm^2

magneticphase(ri, r₀, B::Real) =
@inbounds prefac * B/2 * (ri[1]*ri[2] - r₀[1]*ri[2] + ri[1]*r₀[2])
magneticphase(ri, r₀, B::Zero) = Zero()


text:
using FormFactor
using Test

@testset "FormFactor.jl" begin
    @test  isapprox(mform(1,0,0,1,0,1.0,0.66),0.813233 ;atol=.000001)
    @test  isapprox(mform(2,0,0,2,0,2.0,0.66),-0.0168244 ;atol=.000001)
    @test  isapprox(mform(3,1,0,2,0,1.0,0.66),-0.136364im ;atol=.000001)
    @test  isapprox(mform(3,0,0,3,0,1.0,0.66),0.0433238 ;atol=.000001)
    @test  isapprox(mform(3,0,0,3,0,2.0,0.66),-0.00422855 ;atol=.000001)
    @test  isapprox(mform(4,3,2,4,3,2.0,0.66),-0.00334226 ;atol=.000001)
    @test  isapprox(mform(16,0,0,16,0,2.0,0.66),0.0000861344 ;atol=.000001)
    @test  isapprox(mform(16,0,0,16,0,1.0,0.66),0.0000898579 ;atol=.000001)
    @test  isapprox(mform(19,18,18,19,18,1.0,0.66),7.7825e-33 ;atol=.000001)
end