text:
@block SimonDanisch [interaction, record_events] begin
    @cell "Interaction with Mouse" [interactive, camera, scatter, lines, marker, record] begin
        using LinearAlgebra
        scene = Scene(raw = true, camera = cam2d!, resolution = (500, 500))
        r = range(0, stop = 3, length = 4)
        the_time = Node(time())
        last_open = false
        @async while true
            global last_open
            the_time[] = time()
            # this is a bit awkward, since the isopen(scene) is false
            # as long as the scene isn't displayed
            last_open && !isopen(scene) && break
            last_open = isopen(scene)
            sleep(1/25)
        end
        pos = lift(scene.events.mouseposition, the_time) do mpos, t
            map(LinRange(0, 2pi, 60)) do i
                circle = Point2f0(sin(i), cos(i))
                mouse = to_world(scene, Point2f0(mpos))
                secondary = (sin((i * 10f0) + t) * 0.09) * normalize(circle)
                (secondary .+ circle) .+ mouse
            end
        end
        lines!(scene, pos)
        p1 = scene[end]
        p2 = scatter!(
            scene,
            pos, markersize = 0.1f0,
            marker = :star5,
            color = p1[:color],
        )[end]
        center!(scene)
        t = Theme(raw = true, camera = campixel!)
        b1 = button(t, "color")
        b2 = button(t, "marker")
        msize = slider(t, 0.1:0.01:0.5)
        on(b1[end][:clicks]) do c
            p1[:color] = rand(RGBAf0)
        end
        markers = ('π', '😹', '⚃', '◑', '▼')
        on(b2[end][:clicks]) do c
            p2[:marker] = markers[rand(1:5)]
        end
        on(msize[end][:value]) do val
            p2[:markersize] = val
        end
        RecordEvents(hbox(
            vbox(b1, b2, msize),
            scene
        ), @replace_with_a_path)
    end

    @cell "Mouse Picking" [scatter, heatmap, interactive] begin
        img = rand(100, 100)
        scene = Scene(resolution = (500, 500))
        heatmap!(scene, img, scale_plot = false)
        clicks = Node(Point2f0[(0,0)])
        on(scene.events.mousebuttons) do buttons
           if ispressed(scene, Mouse.left)
               pos = to_world(scene, Point2f0(scene.events.mouseposition[]))
               push!(clicks, push!(clicks[], pos))
           end
           return
        end
        scatter!(scene, clicks, color = :red, marker = '+', markersize = 10)
        RecordEvents(scene, @replace_with_a_path)
    end

    @cell "sliders" [scatter, slider, vbox] begin
        s1 = slider(LinRange(0.01, 1, 100), raw = true, camera = campixel!, start = 0.3)
        s2 = slider(LinRange(-2pi, 2pi, 100), raw = true, camera = campixel!)
        data = lift(s2[end][:value]) do v
            map(LinRange(0, 2pi, 100)) do x
                4f0 .* Point2f0(sin(x) + (sin(x * v) .* 0.1), cos(x) + (cos(x * v) .* 0.1))
            end
        end
        p = scatter(data, markersize = s1[end][:value])

        RecordEvents(
            hbox(p, vbox(s1, s2), parent = Scene(resolution = (500, 500))),
            @replace_with_a_path
        )
    end

    @cell "Mouse Hover" [lines, hover, lift, poly, translate, text, popup, on] begin
        using Colors, Observables
        r = range(0, stop=5pi, length=100)
        scene = Scene(resolution = (500, 500))
        lines!(scene, r, sin.(r), linewidth = 3)
        lineplot = scene[end]
        visible = node(:visible, false)
        poprect = lift(scene.events.mouseposition) do mp
            FRect((mp .+ 5), 250, 40)
        end
        textpos = lift(scene.events.mouseposition) do mp
            Vec3f0((mp .+ 5 .+ (250/2, 40 / 2))..., 120)
        end
        popup = poly!(campixel(scene), poprect, raw = true, color = :white, strokewidth = 2, strokecolor = :black, visible = visible)
        rect = popup[end]
        translate!(rect, Vec3f0(0, 0, 100))
        text!(popup, "( 0.000,  0.000)", textsize = 30, position = textpos, color = :darkred, align = (:center, :center), raw = true, visible = visible)
        text_field = popup[end]
        scene
        x = Node(false)
        on(scene.events.mouseposition) do event
            plot, idx = mouse_selection(scene)
            if plot == lineplot && idx > 0
                visible[] = true
                text_field[1] = sprint(io-> print(io, round.(Float64.(Tuple(lineplot[1][][idx])), digits = 3)))
            else
                visible[] = false
            end
            return
        end
        RecordEvents(scene, @replace_with_a_path)
    end

    @cell "GUI for exploring Lorenz equation" [vbox, hbox, meshscatter, slider, textslider, colorswatch] begin
        using Colors
        using AbstractPlotting: textslider, colorswatch
        s1, a = textslider(0f0:50f0, "a", start = 13)
        s2, b = textslider(-20f0:20f0, "b", start = 10)
        s3, c = textslider(0f0:20f0, "c", start = 2)
        s4, d = textslider(range(0.0, stop = 0.02, length = 100), "d", start = 0.01)
        s5, scales = textslider(range(0.01, stop = 0.5, length = 100), "scale", start = 0.1)
        s6, colorsw, pop = colorswatch()

        function lorenz(t0, a, b, c, h)
            Point3f0(
                t0[1] + h * a * (t0[2] - t0[1]),
                t0[2] + h * (t0[1] * (b - t0[3]) - t0[2]),
                t0[3] + h * (t0[1] * t0[2] - c * t0[3]),
            )
        end
        # step through the `time`
        function lorenz(array::Vector, a = 5.0 ,b = 2.0, c = 6.0, d = 0.01)
            t0 = Point3f0(0.1, 0, 0)
            for i = eachindex(array)
                t0 = lorenz(t0, a,b,c,d)
                array[i] = t0
            end
            array
        end
        n1, n2 = 18, 30
        N = n1*n2
        args_n = (a, b, c, d)
        v0 = lorenz(zeros(Point3f0, N), to_value.(args_n)...)
        positions = lift(lorenz, Node(v0), args_n...)
        rotations = lift(diff, positions)
        rotations = lift(x-> push!(x, x[end]), rotations)

        mesh_scene = meshscatter(
            positions,
            markersize = scales, rotation = rotations,
            intensity = collect(range(0f0, stop = 1f0, length = length(positions[]))),
            color = colorsw
        )
        parent = Scene(resolution = (1000, 800))
        vbox(
            hbox(s1, s2, s3, s4, s5, s6),
            mesh_scene, parent = parent
        )
        RecordEvents(parent, @replace_with_a_path)
    end

    @cell "Edit Polygon" [poly, node, on, events] begin
        points = node(:poly, Point2f0[(0, 0), (0.5, 0.5), (1.0, 0.0)])
        scene = Scene(resolution = (500, 500))
        poly!(scene, points, strokewidth = 2, strokecolor = :black, color = :skyblue2, show_axis = false, scale_plot = false)
        scatter!(points, color = :white, strokewidth = 10, markersize = 0.05, strokecolor = :black, raw = true)
        pplot = scene[end]
        push!(points[], Point2f0(0.6, -0.3))
        points[] = points[]
        function add_move!(scene, points, pplot)
            idx = Ref(0); dragstart = Ref(false); startpos = Base.RefValue(Point2f0(0))
            on(events(scene).mousedrag) do drag
                if ispressed(scene, Mouse.left)
                    if drag == Mouse.down
                        plot, _idx = mouse_selection(scene)
                        if plot == pplot
                            idx[] = _idx; dragstart[] = true
                            startpos[] = to_world(scene, Point2f0(scene.events.mouseposition[]))
                        end
                    elseif drag == Mouse.pressed && dragstart[] && checkbounds(Bool, points[], idx[])
                        pos = to_world(scene, Point2f0(scene.events.mouseposition[]))
                        points[][idx[]] = pos
                        points[] = points[]
                    end
                else
                    dragstart[] = false
                end
                return
            end
        end

        function add_remove_add!(scene, points, pplot)
            on(events(scene).mousebuttons) do but
                if ispressed(but, Mouse.left) && ispressed(scene, Keyboard.left_control)
                    pos = to_world(scene, Point2f0(events(scene).mouseposition[]))
                    push!(points[], pos)
                    points[] = points[]
                elseif ispressed(but, Mouse.right)
                    plot, idx = mouse_selection(scene)
                    if plot == pplot && checkbounds(Bool, points[], idx)
                        deleteat!(points[], idx)
                        points[] = points[]
                    end
                end
                return
            end
        end
        add_move!(scene, points, pplot)
        add_remove_add!(scene, points, pplot)
        center!(scene)
        RecordEvents(scene, @replace_with_a_path)
    end

    @cell "Add and change points" [heatmap, on, ispressed, to_world, scatter, center] begin
        using LinearAlgebra
        img = rand(100, 100)
        scene = Scene(scale_plot = false, resolution = (500, 500))
        heatmap!(scene, img)
        clicks = Node(Point2f0[(0, 0)])
        blues = Node(Point2f0[])
        on(scene.events.mousebuttons) do buttons
            if ispressed(scene, Mouse.left)
                pos = to_world(scene, Point2f0(scene.events.mouseposition[]))
                found = -1
                c = clicks[]
                for i in 1:length(c)
                   if norm(pos - c[i]) < 1
                       found = i
                   end
                end
                if found >= 1
                    blues[] = push!(blues[], pos)
                    deleteat!(clicks[], found)
                else
                    push!(clicks[], pos)
                end
                clicks[] = clicks[]
           end
           return
        end
        t = Theme(markersize = 10, raw = true)
        scatter!(scene, t, clicks, color = :red, marker = '+')
        red_clicks = scene[end]
        scatter!(scene, t, blues, color = :blue, marker = 'o')
        center!(scene)
        RecordEvents(scene, @replace_with_a_path)
    end
    @cell "Orbit Diagram" [scatter, slider] begin
        # example by @datseris
        using Observables
        growth(🐇, 🥕) = 🐇 * 🥕 * (1.0 - 🐇)
        function orbitdiagram(growth, r1, r2, n = 500, a = zeros(1000, n); T = 1000)
            rs = range(r1, stop = r2, length = 1000)
            for (j, r) in enumerate(rs)
                x = 0.5
                for _ in 1:T; x = growth(x, r); end
                for i in 1:n
                    x = growth(x, r)
                    @inbounds a[j, i] = x
                end
            end
            rs, a
        end
        r1 = slider(0:0.001:4, raw = true, camera=campixel!, start = 0.0)
        r2 = slider(0:0.001:4, raw = true, camera=campixel!, start = 4)
        n1 = 500; n2 = 1000; a = zeros(n2, n1)
        positions = Vector{Point2f0}(undef, n1 * n2)
        r1node, r2node = r1[end][:value], r2[end][:value]
        r1r2 = async_latest(lift(tuple, r1node, r2node))
        pos = lift(r1r2) do (r1, r2,)
            global a
            rs, a = orbitdiagram(growth, r1, r2, size(a, 2), a)
            dim = size(a, 2)
            for (i, r) in enumerate(rs)
                positions[((i-1)*dim) + 1 : (i*dim)] .= Point2f0.(r, view(a, i, :))
            end
            positions
        end
        p = scatter(
            pos, markersize = 0.006, color = (:black, 0.2),
            show_axis = false
        )
        onany(pos) do pos
            # faster to give the boundingbox ourselves, since otherwise we'll need to
            # loop over pos!
            #AbstractPlotting.update_limits!(p) # <- would calculate bb
            AbstractPlotting.update_limits!(p, FRect(r1node[], 0, r2node[] - r1node[], 1))
            AbstractPlotting.update!(p)
        end
        scene = hbox(
            p,
            vbox(r1, r2)
        )
        RecordEvents(scene, @replace_with_a_path)
    end

    @cell "Robot Arm" [slider, interactive, linesegments, vbox] begin

        using AbstractPlotting: Mesh, Scene, LineSegments, translate!, rotate!, vbox, hbox, qrotation, mesh!
        using GeometryTypes: HyperRectangle, Vec3f0, Point3f0, Sphere
        using StaticArrays: SVector
        using AbstractPlotting: textslider
        using Observables: on

        """
            example by @pbouffard from JuliaPlots/Makie.jl#307
            https://github.com/pbouffard/miniature-garbanzo/
        """

        function triad!(scene, len; translation = (0f0,0f0,0f0), show_axis = false)
            ret = linesegments!(
                scene, [
                    Point3f0(0,0,0) => Point3f0(len,0,0),
                    Point3f0(0,0,0) => Point3f0(0,len,0),
                    Point3f0(0,0,0) => Point3f0(0,0,len)
                ],
                color = [:red, :green, :blue],
                linewidth = 3, show_axis = show_axis
            )[end]
            translate!(ret, translation)
            return ret
        end
        # Joint vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
        mutable struct Joint
            scene::Scene
            triad::LineSegments
            # link::Mesh
            angle::Float32
            axis::Vec3f0
            offset::Vec3f0
        end

        s = Scene()

        function Joint(s::Scene)
            newscene = Scene(s)
            triad = triad!(newscene, 1)
            Joint(newscene, triad, 0f0, (0, 1, 0), (0, 0, 0))
        end

        function Joint(j::Joint; offset::Point3f0=(0,0,0), axis=(0, 1, 0), angle=0)
            jnew = Joint(j.scene)
            translate!(jnew.scene, j.offset)
            linesegments!(jnew.scene, [Point3f0(0) => offset], linewidth=4, color=:magenta, show_axis=false)
            jnew.axis = axis
            jnew.offset = offset
            setangle!(jnew, angle)
            return jnew
        end

        function setangle!(j::Joint, angle::Real)
            j.angle = angle
            rotate!(j.scene, qrotation(j.axis, angle))
        end

        # Joint ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

        joints = Vector{Joint}()
        links = Float32[5, 5]
        #triad!(s, 10; show_axis=true)

        push!(joints, Joint(s))
        joints[1].axis = (0,0,1) # first joint is yaw
        joints[1].offset = (0, 0, 1)
        push!(joints, Joint(joints[end]; offset=Point3f0(3,0,0), axis=(0,1,0), angle=-pi/4)) # Pitch
        push!(joints, Joint(joints[end]; offset=Point3f0(3,0,0), axis=(0,1,0), angle=pi/2)) # Pitch
        push!(joints, Joint(joints[end]; offset=Point3f0(1,0,0), axis=(0,1,0), angle=-pi/4)) # Pitch
        push!(joints, Joint(joints[end]; offset=Point3f0(1,0,0), axis=(0,0,1))) # Yaw
        push!(joints, Joint(joints[end]; offset=Point3f0(0,0,0), axis=(1,0,0))) # Roll

        sliders = []
        vals = []
        for i = 1:length(joints)
            slider, val = textslider(-180.0:1.0:180.0, "Joint $(i)", start=rad2deg(joints[i].angle))
            push!(sliders, slider)
            push!(vals, val)
            on(val) do x
                setangle!(joints[i], deg2rad(x))
            end
        end

        # Add sphere to end effector:
        mesh!(joints[end].scene, Sphere(Point3f0(0.5, 0, 0), 0.25f0), color=:cyan, show_axis=false)
        center!(s)
        RecordEvents(vbox(hbox(sliders...), s), @replace_with_a_path)
    end

    @cell "Earth & Ships" [slider, interactive, lines, mesh, vbox] begin

        using AbstractPlotting: textslider
        using GeometryTypes, FileIO
        using LinearAlgebra

        """
            example by @pbouffard from JuliaPlots/Makie.jl#307
            https://github.com/pbouffard/miniature-garbanzo/
        """

        const rearth_km = 6371.1f0
        const mearth_kg = 5.972e24
        const rmoon_km = 1738.1f0
        const mmoon_kg = 7.34e22
        const rship_km = 500f0
        const mship_kg = 1000.0
        const dbetween_km = 378_000.0f0
        const radius_mult = 1.0f0
        const timestep = 1.0f0

        mutable struct Ship
            mass_kg::Float32
            position_m::Vec3f0
            velocity_mps::Vec3f0
            color::Symbol
            mesh::Mesh
        end

        function moveto!(ship::Ship, (x, y, z))
            ship.position_m = Vec3f0(x, y, z)
            translate!(ship.mesh, x, y, z)
            ship.position_m, ship.velocity_mps
        end

        function orbit(r; steps=80)
            return [Point3f0(r*cos(x), r*sin(x), 0) for x in range(0, stop=2pi, length=steps)]
        end

        function loadordownload(localfname, url)
            return isfile(localfname) ? load(localfname) : load(download(url, localfname))
        end

        # http://corysimon.github.io/articles/uniformdistn-on-sphere/
        function makestars(n)
            return map(1:n) do i
                v = [0, 0, 0]  # initialize so we go into the while loop
                while norm(v) < .0001
                    v = randn(Point3f0)
                end
                v = v / norm(v)  # normalize to unit norm
                v
            end
        end

        function makecontrols()
            orbit_slider, orbit_speed = textslider(0f0:10f0, "Speed", start=1.0)
            scale_slider, scale = textslider(1f0:20f0, "Scale", start=10.0)
            return orbit_slider, orbit_speed, scale_slider, scale
        end

        function makeships(scene, N)
            return map(1:N) do i
                sm = mesh!(scene, Sphere(Point3f0(0), rship_km*radius_mult), color=:green, show_axis=false)[end]
                ship = Ship(mship_kg, [0, 0, 0], [0, 0, 0], :green, sm)
                moveto!(ship, (100000f0 * randn(Float32), dbetween_km/2 + randn(Float32) * 50000f0, 50000f0*randn(Float32)))
                ship.velocity_mps = [40000*randn(Float32), -1000*randn(Float32), 1000*randn(Float32)]
                ship
            end
        end

        function makeobjects(scene)
            earthfname = "bluemarble-2048.png"
            earthurl = "https://svs.gsfc.nasa.gov/vis/a000000/a002900/a002915/" * earthfname
            earthbitmap = loadordownload(earthfname, earthurl)
            earth = mesh!(scene, Sphere(Point3f0(0), rearth_km*radius_mult), color=earthbitmap)[end]

            moonfname = "phases.0001_print.jpg"
            moonurl = "https://svs.gsfc.nasa.gov/vis/a000000/a004600/a004675/" * moonfname
            moonbitmap = loadordownload(moonfname, moonurl)
            moon = mesh!(scene, Sphere(Point3f0(0), rmoon_km*radius_mult), color=moonbitmap)[end]
            translate!(moon, Vec3f0(0, dbetween_km, 0))
            orb = lines!(scene, orbit(dbetween_km), color=:gray)
            # stars = scatter!(1000000*makestars(1000), color=:white, markersize=2000)

            return earth, moon #, stars#, ships
        end

        function spin(scene, orbit_speed, scale, moon, earth, ships)
            center!(scene) # update far / near automatically.
            # could also be done manually:
            # cam.near[] = 100
            # cam.far[] = 1.5*dbetween_km
            update_cam!(scene, Vec3f0(1.3032e5, 7.12119e5, 3.60022e5), Vec3f0(0, 0, 0))
            θ = 0.0
            # wait untill window is open
            while !isopen(scene)
                sleep(0.1)
            end
            while isopen(scene)
                scale!(moon, (scale[], scale[], scale[]))
                scale!(earth, (scale[], scale[], scale[]))
                for ship in ships
                    scale!(ship.mesh, (scale[], scale[], scale[]))
                end
                translate!(moon, Vec3f0(-dbetween_km*sin(θ/28), dbetween_km*cos(θ/28), 0))
                rotate!(earth, Vec3f0(0, 0, 1), θ)
                rotate!(moon, Vec3f0(0, 0, 1), π/2 + θ/28)
                sleep(0.01)
                θ += 0.05 * orbit_speed[]
                timestep = 10 * orbit_speed[]
                gravity(timestep, ships, earth, moon)
            end
        end

        function gravity(timestep, ships, earth, moon)
            for ship in ships
                rship = ship.position_m #translation(ship)[]
                rearth = translation(earth)[]
                rmoon = translation(moon)[]
                rship_earth = (rearth - rship) * 1e3 # convert to m
                rship_moon = (rmoon - rship) * 1e3 # convert to m
                G = 6.67408e-11 # m^3⋅kg^-2⋅s^-1
                nrship_earth = max(0.1, norm(rship_earth))
                nrship_moon = max(0.1, norm(rship_moon))
                fearth = G * (mship_kg * mearth_kg)/(nrship_earth^2) * rship_earth/nrship_earth
                fmoon = G * (mship_kg * mmoon_kg)/(nrship_moon^2) * rship_earth/nrship_earth
                ftot = fearth + fmoon
                ship.velocity_mps += timestep * ftot
                drship = timestep * ship.velocity_mps / 1e3
                rship = rship + drship
                moveto!(ship, rship)
            end
        end
        universe = Scene(backgroundcolor=:black, show_axis=false, resolution=(2048,1024))
        ships = makeships(universe, 600)
        earth, moon = makeobjects(universe);
        orbit_slider, orbit_speed, scale_slider, scale = makecontrols();
        controls = (orbit_slider, scale_slider);
        scene = hbox(vbox(controls...), universe)
        universe.center = false
        task = @async spin(universe, orbit_speed, scale, moon, earth, ships)
        scene.center = false
        RecordEvents(scene, @replace_with_a_path)
    end

    @cell "Interactive Differential Equation" [lines, slider, animated, record, diffeq, interactive] begin

        using DifferentialEquations, ParameterizedFunctions
        import AbstractPlotting: textslider
        lorenz = @ode_def Lorenz begin           # define the system
            dx = σ * (y - x)
            dy = x * (ρ - z) - y
            dz = x * y - β*z
        end σ ρ β

        u0 = [1.0,0.0,0.0]                       # initial conditions
        tspan0 = (0.0,100.0)                      # initial timespan
        p0 = [10.0,28.0,8/3]                      # initial parameters
        prob = ODEProblem(lorenz, u0, tspan0, p0)  # define the problem

        ## setup sliders and plotting

        "The order of magnitude to range between."
        OME = 8

        sσ, oσ = textslider(exp10.(-OME:0.001:OME), "σ", start = p0[1]);

        sρ, oρ = textslider(exp10.(-OME:0.001:OME), "ρ", start = p0[2]);

        sβ, oβ = textslider(exp10.(-OME:0.001:OME), "β", start = p0[3]);

        st, ot = textslider(exp10.(-OME:0.001:OME), "tₘₐₓ", start = tspan0[end]);

        sr, or = textslider(100:10000, "resolution", start = 2000);

        trange = lift(ot, or) do tmax, resolution
            LinRange(0.0, tmax, resolution)
        end

        data = lift(oσ, oρ, oβ, trange) do σ, ρ, β, ts
            Point3f0.(
                solve(
                    remake(
                        prob;
                        p = [σ, ρ, β],
                        tspan = (ts[1], ts[end])
                    )
                )(ts).u
            )  # change to fit the dimensionality - maybe even return 2 arrays, or a set of `Point2`s...
        end
        parent = Scene(resolution = (1000, 500))
        three = lines(data, linewidth = 1, transparency = true, color = ("#fe4a49", 0.6))
        meshscatter!(three, data, markersize = 0.1, color = :gray)
        scene = vbox(hbox(sσ, sρ, sβ, st, sr), three, parent = parent)
        RecordEvents(scene, @replace_with_a_path)
    end
end

function record_example_events()
    eval_examples(:record_events, :scatter, :slider, :interactive, outputfile = event_path) do example, value
        record_events(value.scene, value.path) do
            wait(value.scene)
        end
    end
end
function record_example(title = "Orbit Diagram")
    set_theme!(resolution = (500, 500))
    idx = findfirst(x-> x.title == title, database)
    entry = database[idx]
    value = MakieGallery.eval_example(entry, outputfile = MakieGallery.event_path)
    last = AbstractPlotting.use_display[]
    AbstractPlotting.inline!(false)
    record_events(value.scene, value.path) do
        wait(value.scene)
    end
    AbstractPlotting.use_display[] = last
end
# record_example("Interactive Differential Equation")


text:
module YahooFinance

using Dates
using DataFrames
using TimeSeries
using HTTP
using JSON

export
    get_symbols

include("get_symbols.jl")

end

ext: jl


text:
function solve()
    x = parse(Int, readline())
    x ≥ 30 ? "Yes" : "No"
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
module MyModules
# CourErr.jl functions
export load_data, load_data_select, gM, mystars
# MarkovChain.jl functions
export Par, Param, warm_start, jump, gchain
# Parallel_mod.jl functions
export shared_fill_chain!, fill_chain_sem!, myrange_sem
export shared_hM_fill!, hM_chunk!, myrange, objfun_par, objfun_bbo_par
export myrange_y, fill_table!, shared_fill_table!
# Testing.jl functions
export objfun, objfun_bbo, my_dist, RP_Cournot_test, acc_rate
export RP_Cournot_test_par

# Testing alternative transition kernel and alternative acceptance rule
export objfun_par_alt, objfun_bbo_par_alt, hM_chunk_alt!, RP_Cournot_test_par_alt
export gchain_algo, jump_algo, fill_chain_sem_algo!, fill_chain_sem_alt!, gchain_alt

include("CourErr.jl")
include("MarkovChain.jl")
include("Parallel_mod.jl")
include("Testing.jl")

end  # module MyModules

max_stars_repo_name: hans-mtz/SCR

max_issues_repo_name: hans-mtz/SCR

max_forks_repo_name: hans-mtz/SCR


text:
using Turing, Random, Test
using Turing: @VarName
using Distributions: Normal
using StatsFuns

dir = splitdir(splitdir(pathof(Turing))[1])[1]
include(dir*"/test/test_utils/AllUtils.jl")

@turing_testset "util.jl" begin
    i = 1
    @test @VarName(s)[1:end-1] == (:s,())
    @test @VarName(x[1,2][1+5][45][3][i])[1:end-1] ==
        (:x,("[1, 2]","[6]","[45]","[3]","[1]"))
    @test StatsFuns.logistic(1.1) == 1.0 / (exp(-1.1) + 1.0)
    @test StatsFuns.logit(0.3) ≈ -0.8472978603872036 atol=1e-9
    @test !isnan(StatsFuns.logit(1.0))
    @test isinf(StatsFuns.logit(1.0))
    @test !isnan(StatsFuns.logit(0.0))
    @test isinf(StatsFuns.logit(0.0))

    @test Turing.Core.data(Dict(:y => 1, :x => 2, :z => 1), [:x, :y]) ==
        (x = 2, y = 1)
end


text:
using CompScienceMeshes, BEAST
Γ = readmesh(joinpath(@__DIR__,"sphere2.in"))
@show numcells(Γ)
X = raviartthomas(Γ)

Δt, Nt = 0.6, 200
T = timebasisshiftedlagrange(Δt, Nt, 3)
U = timebasisdelta(Δt, Nt)

V = X ⊗ T
W = X ⊗ U

duration = 20 * Δt
delay = 1.5 * duration
amplitude = 1.0
gaussian = creategaussian(duration, delay, amplitude)

direction, polarisation = ẑ, x̂
E = BEAST.planewave(polarisation, direction, derive(gaussian), 1.0)

SL = TDMaxwell3D.singlelayer(speedoflight=1.0, numdiffs=1)

@hilbertspace j
@hilbertspace j′
tdefie = @discretise SL[j′,j] == -1.0E[j′]   j∈V  j′∈W
xefie = solve(tdefie)

Xefie, Δω, ω0 = fouriertransform(xefie, Δt, 0.0, 2)
ω = collect(ω0 .+ (0:Nt-1)*Δω)
_, i1 = findmin(abs.(ω.-1.0))

ω1 = ω[i1]
ue = Xefie[:,i1] / fouriertransform(gaussian)(ω1)


text:
module SDLParameters

using Distributions
using ShaleDrillingLikelihood
using Dates
using SparseArrays
using LinearAlgebra
using Random

using StatsModels: @formula
using CategoricalArrays: CategoricalVector
using MixedModels

using ShaleDrillingLikelihood: simulate,
    cost, extend, revenue, scrap,
    lrmean, lrvar, lrstd,
    zero_out_small_probs,
    extra_parm,
    LeasedProblem,
    idx_produce_ψ,
    idx_drill_ψ,
    idx_produce_g,
    idx_drill_g,
    idx_produce_t,
    idx_drill_t,
    idx_drill_D,
    idx_produce_D,
    total_wells_drilled,
    DataDynamicDrill,
    num_choices,
    _num_x,
    DrillingRevenueUnconstrained

using Base: product
using LinearAlgebra: checksquare
using ArgParse
using Base.Meta: parse
import ArgParse: parse_item # SDL

import ShaleDrillingLikelihood: ichars, zchars,
    DataRoyalty, DataProduce, DataDrillPrimitive, DataDrill

include("thetas.jl")
include("model-components.jl")
include("make-models.jl")
include("grids.jl")
include("data-import.jl")
include("argparse.jl")
include("argparse-simulations.jl")

end


text:
# ---
# title: 1091. Shortest Path in Binary Matrix
# id: problem1091
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Breadth-first Search
# link: <https://leetcode.com/problems/shortest-path-in-binary-matrix/description/>
# hidden: true
# ---
#
# In an N by N square grid, each cell is either empty (0) or blocked (1).
#
# A  _clear  path from top-left to bottom-right_ has length `k` if and only if
# it is composed of cells `C_1, C_2, ..., C_k` such that:
#
#   * Adjacent cells `C_i` and `C_{i+1}` are connected 8-directionally (ie., they are different and share an edge or corner)
#   * `C_1` is at location `(0, 0)` (ie. has value `grid[0][0]`)
#   * `C_k` is at location `(N-1, N-1)` (ie. has value `grid[N-1][N-1]`)
#   * If `C_i` is located at `(r, c)`, then `grid[r][c]` is empty (ie. `grid[r][c] == 0`).
#
# Return the length of the shortest such clear path from top-left to bottom-
# right.  If such a path does not exist, return -1.
#
#
#
# **Example 1:**
#
#
#
#     Input: [[0,1],[1,0]]
#     ![](https://assets.leetcode.com/uploads/2019/08/04/example1_1.png)
#
#     Output: 2
#     ![](https://assets.leetcode.com/uploads/2019/08/04/example1_2.png)
#
#
# **Example 2:**
#
#
#
#     Input: [[0,0,0],[1,1,0],[1,1,0]]
#     ![](https://assets.leetcode.com/uploads/2019/08/04/example2_1.png)
#
#     Output: 4
#     ![](https://assets.leetcode.com/uploads/2019/08/04/example2_2.png)
#
#
#
#
# **Note:**
#
#   1. `1 <= grid.length == grid[0].length <= 100`
#   2. `grid[r][c]` is `0` or `1`
#
#
## @lc code=start
using LeetCode

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


text:
function frontendfilterdesign(filter_type_num, zfunc, L, fig_num,
                                savefig_flag, plot_output_flag::Bool,
                                max_iters::Int, verbose_flag::Bool,
                                debug_flag::Bool)

    config = getdefaultoptparameters(L)
    config.max_iters = max_iters
    config.verbose_flag = verbose_flag
    config.plot_flag = debug_flag

    return frontendfilterdesign(filter_type_num, zfunc, L, fig_num,
                                    savefig_flag, plot_output_flag,
                                    config)

end

function frontendfilterdesign(filter_type_num, zfunc, L, fig_num,
                                savefig_flag, plot_output_flag::Bool)

    config = getdefaultoptparameters(L)

    return frontendfilterdesign(filter_type_num, zfunc, L, fig_num,
                                    savefig_flag, plot_output_flag,
                                    config)

end

function frontendfilterdesign(filter_type_num, zfunc, L, fig_num,
                                savefig_flag, plot_output_flag::Bool,
                                opt_parameters_file_name::String)

    config = getoptparameters(L,opt_parameters_file_name)

    return frontendfilterdesign(filter_type_num, zfunc, L, fig_num,
                                    savefig_flag, plot_output_flag,
                                    config)
end

function frontendfilterdesign(filter_type_num, zfunc, L, fig_num,
                                savefig_flag, plot_output_flag,
                                config::MiniMaxConfigType)

    𝓧 = collect( Chebyshev2ndnode(k,L+1) for k = 0:L+1 )

    ##
    wfunc = getweightfunc(filter_type_num)
    dfunc = xx->zfunc(acos(xx))

    𝓧_selection_string = "Chebyshev"
    candidate_multiple = 1

    @time 𝓧_new, h_history,
            iters_ran, status_msg, fig_num = designfilter(  config,
                                                            L,
                                                            candidate_multiple,
                                                            wfunc,
                                                            dfunc,
                                                            𝓧_selection_string,
                                                            fig_num)

    #
    println("Completed mini-max algorithm. Status: ", status_msg)

    println("Start to package filter solution.")
    @time h, fig_num = packagefiltersolution(𝓧_new, h_history, wfunc, zfunc, dfunc,
                                        filter_type_num, fig_num, savefig_flag, plot_output_flag)
    println("Finished packaging filter solution.")

    return h, 𝓧_new, fig_num
end

function frontendfilterdesignfast(filter_type_num, zfunc, L, fig_num,
                                savefig_flag, plot_output_flag,
                                config::MiniMaxConfigType)

    𝓧_tmp = SharedArray{Float64}(L+2)
    @sync @distributed for k = 0:L+1
        𝓧_tmp[k+1] = convert(Float64, Chebyshev2ndnode(k,L+1))
    end
    𝓧 = convert(Vector{BigFloat}, 𝓧_tmp)

    ##
    wfunc = getweightfunc(filter_type_num)
    dfunc = xx->zfunc(acos(xx))

    𝓧_selection_string = "Chebyshev"
    candidate_multiple = 1

    println("Begin mini-max algorithm.")
    @time 𝓧_new, h_history,
            iters_ran, status_msg, fig_num = designfilter(  config,
                                                            L,
                                                            candidate_multiple,
                                                            wfunc,
                                                            dfunc,
                                                            𝓧_selection_string,
                                                            fig_num)

    #

    println("Completed mini-max algorithm. Status: ", status_msg)

    println("Start to package filter solution.")
    @time h = packagefiltersolutionfast(𝓧_new, h_history, wfunc, zfunc, dfunc,
                                        filter_type_num)
    println("Finished packaging filter solution.")

    return h, 𝓧_new, fig_num
end

function visualizeoptimzationsolution(xq, 𝓧, f, fig_num)

    𝑤_new = get𝑤scaled(𝓧)
    h_new = geth(𝑤_new, constantfunc, f, 𝓧)

    f_itp = xx->getp(xx, 𝑤_new, h_new, constantfunc, f, 𝓧)
    f_itp_xq = collect( f_itp(xq[i]) for i = 1:length(xq) )
    f_itp_𝓧 = collect( f_itp(𝓧[i]) for i = 1:length(𝓧) )

    𝑒 = xx->(f(xx)-f_itp(xx))
    𝑒_xq = collect( 𝑒(xq[i]) for i = 1:length(xq) )
    𝑒_𝓧 = collect( 𝑒(𝓧[i]) for i = 1:length(𝓧) )

    # PyPlot.figure(fig_num)
    # fig_num += 1
    # PyPlot.plot(xq, f.(xq), "--", label = "target")
    # PyPlot.plot(xq, f_itp_xq, label = "solution")
    # PyPlot.plot(𝓧, f_itp_𝓧, ".")
    # PyPlot.title("f vs. solution")
    # PyPlot.legend()
    #
    # PyPlot.figure(fig_num)
    # fig_num += 1
    # PyPlot.plot(xq, 𝑒_xq)
    # PyPlot.plot(𝓧, 𝑒_𝓧, ".")
    # PyPlot.title("error")
    # PyPlot.legend()

    return f_itp_xq, f_itp_𝓧, 𝑒_xq, 𝑒_𝓧, fig_num
end

function frontendoptimization(targetfunc, L,
            opt_parameters_file_name::String = "0")

    𝓧 = collect( Chebyshev2ndnode(k,L+1) for k = 0:L+1 )

    config = getdefaultoptparameters(L)
    if opt_parameters_file_name != "0"
        config = getoptparameters(L,opt_parameters_file_name)
    end

    𝓧_new, h_history, iters_ran, status_msg, fig_num_unused = fitminimax(  constantfunc,
                                                                targetfunc,
                                                                𝓧,
                                                                1,
                                                                config)
    return 𝓧_new
end

max_stars_repo_name: RoyCCWang/FIRRemez

max_issues_repo_name: RoyCCWang/FIRRemez

max_forks_repo_name: RoyCCWang/FIRRemez


text:
using Test
using Random, Printf
using GRUtils

Random.seed!(111)
examplesdir = joinpath(dirname(pathof(GRUtils)), "../examples/docstrings")

GRUtils.GR.inline("pdf")

functionlist = (
    ("Line plots", ("plot", "oplot", "stair", "plot3", "polar")),
    ("Scatter plots", ("scatter", "scatter3")),
    ("Stem plots", ("stem",)),
    ("Bar plots", ("barplot",)),
    ("Vector fields", ("quiver", "quiver3")),
    ("Histograms", ("histogram", "polarhistogram", "hexbin")),
    ("Contour plots", ("contour", "contourf", "tricont")),
    ("Surface plots", ("surface", "trisurf", "wireframe")),
    ("Volume rendering", ("volume",)),
    ("Heatmaps", ("heatmap", "polarheatmap", "shade")),
    ("Images", ("imshow",)),
    ("Isosurfaces", ("isosurface",)),
    ("Text", ("annotations",)),
    ("Attributes", ("aspectratio", "ticks")),
    ("Control", ("subplot",)),
    ("Colors", ("background", "colorscheme", "colormap"))
)

# GR.gr3 had an issue in some systems
# win = Sys.iswindows()
# in_gr3 = ("surface", "volume", "shade", "isosurface")

#Plotting functions
@testset "$(functionlist[g][1])" for g in 1:length(functionlist)
    funs = functionlist[g][2]
    for f in funs
        include(joinpath(examplesdir, "$f.jl"))
        @test true
    end
end

include("extratests.jl")

file_path = ENV["GKS_FILEPATH"]
@test isfile(file_path)
GRUtils.GR.reset()
rm(file_path)


text:
#computes ms that appears in the part of the loglikelihood that depends on model parameters
function ms(Ns::Vector{Float64}, λ::Vector{Float64}) #Ns are counts of variant s throut the various rounds
	m = 0.0
	@inbounds @simd for t in 1:length(Ns)-1
		m += Ns[t+1]/λ[t]
	end
	return m
end

#computes Ms that appears in the part of the loglikelihood that depends on model parameters
function Ms(Ns::Vector{Float64})
	M = 0.0
	@inbounds @simd for t in 1:length(Ns)-1
		M += Ns[t]
	end
	return M
end

#computes the minimum value that λ(t) can assume consistent with the data
function λmin(N::Array{Float64,3},t::Int)
	(S, F, T) = size(N)
	@assert t > 0 && t < T
	l = -Inf
	@inbounds for s in 1:S
		if N[s,1,t] == 0
			continue
		end
		l = (l > N[s,1,t+1]/N[s,1,t]) ? l : N[s,1,t+1]/N[s,1,t]
	end
	return l
end

#computes the value of the term of the loglikelihood that depends only on model parameters
function Likelihood_fields(model::Model, data)
	L = 0.0
	@inbounds @simd for s in 1:data.S
		Hs = Energy(model.epistasis, model.x, data.variants[s])
		m = ms(data.counts[s,1,:], model.λ)
		M = Ms(data.counts[s,1,:])
		L += -M*log_oneplusexp(Hs) + (M-m)*Hs
	end
	return L
end

#computes the gradient with respect to model parameters of the loglikelihood
function gradient_Likelihood_fields(model::Model, data)
	grad = zeros(length(model.x))
	@inbounds for s in 1:data.S
		p = pb(model.epistasis, model.x, data.variants[s])
		m = ms(data.counts[s,1,:], model.λ)
		M = Ms(data.counts[s,1,:])
		g = m - p*M

		@simd for i in 1:data.L
			index = field_index(model.epistasis, :h, data.variants[s].sequence[i], i)
			grad[index] += g
		end
		for i in 1:data.L-1
			@simd for j in i+1:data.L
				index = field_index(model.epistasis, :J, data.variants[s].sequence[i], data.variants[s].sequence[j], i, j)
				grad[index] += g
			end
		end
		grad[length(model.x)] += g
	end
	return grad
end

#computes the part of the loglikelihood that depends on λ's'
function Likelihood_lambda(model::Model, data)
	L = 0.0
	@inbounds H = [Energy(model.epistasis, model.x, data.variants[s]) for s in 1:data.S]
	@inbounds for t in 1:data.T-1
		@simd for s in 1:data.S
			L += log_binom(data.counts[s,1,t], data.counts[s,1,t+1]/model.λ[t]) - H[s]*data.counts[s,1,t+1]/model.λ[t]
		end
		L -= data.S*log(model.λ[t])
	end
	return L
end

#computes the gradient of the loglikelihood with respect to λ's'
function gradient_Likelihood_lambda(model::Model, data)
	grad = zeros(model.T-1)
	@inbounds H = [Energy(model.epistasis, model.x, data.variants[s]) for s in 1:data.S]
	@inbounds for t in 1:model.T-1
		@simd for s in 1:data.S
			grad[t] += ( xlogy(data.counts[s,1,t+1]/model.λ[t]^2, data.counts[s,1,t+1]/model.λ[t])
						- xlogy(data.counts[s,1,t+1]/model.λ[t]^2, data.counts[s,1,t]-data.counts[s,1,t+1]/model.λ[t])
						+ H[s]*data.counts[s,1,t+1]/model.λ[t]^2)
		end
		grad[t] -= data.S/model.λ[t]
	end
	return grad
end

#computes the value of the total loglikelihood
function Likelihood(model::Model, data)
	L = 0.0
	@inbounds H = [Energy(model.epistasis, model.x, data.variants[s]) for s in 1:data.S]
	@inbounds for t in 1:model.T-1
		@simd for s in 1:data.S
			L += ( log_binom(data.counts[s,1,t], data.counts[s,1,t+1]/model.λ[t])
					- H[s]*data.counts[s,1,t+1]/model.λ[t] - log_oneplusexp(-H[s])*data.counts[s,1,t] )
		end
		L -= data.S*log(model.λ[t])
	end
	return L
end

#computes the gradient of the loglikelihood with respect to all the parameters that must be learnt
#first the parameters of the model and then the λ's'
function gradient_Likelihood(model::Model, data)
	gf = gradient_Likelihood_fields(model, data)
	gl = gradient_Likelihood_lambda(model, data)
	return vcat(gf, gl)
end


text:
module RobustAdaptiveMetropolisSampler

using LinearAlgebra, Random, Distributions, PDMats, ProgressMeter

export RAM_sample

# The following methods cover different ways to pass in a co-variance matrix
function RAM_sample(logtarget, x0::AbstractVector{<:Number}, M0::AbstractMatrix{<:Real}, n::Int; kwargs...)
    return RAM_sample(logtarget, x0, PDMat(M0), n; kwargs...)
end

function RAM_sample(logtarget, x0::AbstractVector{<:Number}, M0::AbstractVector{<:Real}, n::Int; kwargs...)
    return RAM_sample(logtarget, x0, PDiagMat(abs2.(M0)), n; kwargs...)
end

function RAM_sample(logtarget, x0::AbstractVector{<:Number}, M0::Diagonal{<:Real}, n::Int; kwargs...)
    return RAM_sample(logtarget, x0, PDiagMat(diag(M0)), n; kwargs...)
end

function RAM_sample(logtarget, x0::AbstractVector{<:Number}, M0::Real, n::Int; kwargs...)
    return RAM_sample(logtarget, x0, ScalMat(length(x0), abs2(M0)), n; kwargs...)
end

# Actual sampling code

function RAM_sample(logtarget, x0::AbstractVector{<:Number}, M0::AbstractPDMat, n::Int; opt_α=0.234, γ=2/3, q=Normal(), show_progress::Bool=true)
    length(x0) == size(M0, 1) || error("Covariance matrix M0 must match size of x0.")
    n > 0 || error("n must be larger than 0.")
    0 < opt_α < 1 || error("opt_α must be between 0 and 1.")
    0.5 < γ <= 1 || error("γ must be between 0.5 and 1.")

    x = copy(x0)
    s = cholesky(M0)

    d = length(x0)

    u = zeros(d)
    y = zeros(d)

    log_probability_x = logtarget(x0)

    output_chain = Matrix{Float64}(undef, n, d)

    stats_accepted_values = 0

    progress_meter = Progress(n, show_progress ? 0.1 : Inf)

    # This is a pre-allocated vector used in the loop
    scaled_proposal_vector = Vector{Float64}(undef, d)

    for i in 1:n
        # Step R1
        rand!(q, u)

        y[:] .= x .+ mul!(scaled_proposal_vector, s.L, u)

        # Step R2
        log_probability_y = logtarget(y)

        α = min(1, exp(log_probability_y - log_probability_x))

        if α > rand()
            stats_accepted_values += 1
            x, y = y, x
            log_probability_x = log_probability_y
        end

        # Step R3

        # This is taken from the second paragraph of section 5
        η = min(1, d * i^-γ)

        # Compute the new covariance matrix
        M = s.L * (I + η * (α-opt_α) * (u * u') / norm(u)^2 ) * s.L'

        # The paper has a proof that M is symmetric, so we declare that fact
        # to work around numerical rounding errors
        s = cholesky(Symmetric(M))

        output_chain[i, :] .= x

        next!(progress_meter; showvalues = [(:acceptance_rate,stats_accepted_values/i)])
    end

    return (chain=output_chain, acceptance_rate=stats_accepted_values/n, M=s.L * s.L')
end

end


text:
@marginalrule NormalMeanVariance(:out_μ_v) (m_out::NormalDistributionsFamily, m_μ::PointMass, m_v::PointMass) = begin
    return (out = prod(ProdAnalytical(), NormalMeanVariance(mean(m_μ), mean(m_v)), m_out), μ = m_μ, v = m_v)
end

@marginalrule NormalMeanVariance(:out_μ_v) (m_out::PointMass, m_μ::NormalDistributionsFamily, m_v::PointMass) = begin
    return (out = m_out, μ = prod(ProdAnalytical(), m_μ, NormalMeanVariance(mean(m_out), mean(m_v))), v = m_v)
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

run a CBLIB instance
cblib_dir (defined in JuMP_test.jl) is the directory of CBLIB files
=#

struct CBLIBJuMP{T <: Real} <: ExampleInstanceJuMP{T}
    name::String # filename of CBLIB instance
end

function build(inst::CBLIBJuMP{T}) where {T <: Float64}
    # read in model from file
    model = JuMP.read_from_file(joinpath(cblib_dir, inst.name * ".cbf.gz"))

    # delete integer constraints
    JuMP.delete.(model, JuMP.all_constraints(model, JuMP.VariableRef, MOI.Integer))

    return model
end


text:
using Test, CUDA, Random
import NNlib
using CUDA.CUDNN:
    cudnnPoolingForward,
    cudnnPoolingForward!,
    cudnnPoolingBackward,
    cudnnGetPoolingNdForwardOutputDim,
    cudnnPoolingDescriptor,
        cudnnPoolingDescriptor_t,
        cudnnCreatePoolingDescriptor,
        cudnnSetPoolingNdDescriptor,
        cudnnDestroyPoolingDescriptor,
    cudnnPoolingMode_t,
        CUDNN_POOLING_MAX,                           # 0,
        CUDNN_POOLING_AVERAGE_COUNT_INCLUDE_PADDING, # 1, /* count for average includes padded values */
        CUDNN_POOLING_AVERAGE_COUNT_EXCLUDE_PADDING, # 2, /* count for average does not include padded values */
        CUDNN_POOLING_MAX_DETERMINISTIC,             # 3
    cudnnNanPropagation_t,
        CUDNN_NOT_PROPAGATE_NAN, # 0
        CUDNN_PROPAGATE_NAN,     # 1
    cudnnTensorFormat_t,
        CUDNN_TENSOR_NCHW,        # 0, /* row major (wStride = 1, hStride = w) */
        CUDNN_TENSOR_NHWC,        # 1, /* feature maps interleaved ( cStride = 1 )*/
        CUDNN_TENSOR_NCHW_VECT_C, # 2, /* each image point is vector of element of C, vector length in data type */
    pooldims

@testset "cudnn/pooling" begin
    function pooltest(;
                      mode = CUDNN_POOLING_MAX,
                      nanOpt = CUDNN_NOT_PROPAGATE_NAN,
                      window = 2,
                      padding = 0,
                      stride = window,
                      format = CUDNN_TENSOR_NCHW,
                      dataType = Float32,
                      alpha = 1,
                      beta = 0)
        ax = randn(dataType,12,6,4,2)
        N = ndims(ax)
        window = expand(Val(N-2), window)
        stride = expand(Val(N-2), stride)
        padding = expand(Val(N-2), padding)
        pdims = NNlib.PoolDims(ax, window; padding = padding, stride = stride)
        #=
        if mode == CUDNN_POOLING_AVERAGE_COUNT_EXCLUDE_PADDING
            @warn "Pool mode=$mode not yet implemented in NNlib, using INCLUDE instead. See https://github.com/FluxML/NNlib.jl/issues/218" maxlog=1
        end
        if mode == CUDNN_POOLING_MAX_DETERMINISTIC
            @warn "Pool mode=$mode not yet implemented in NNlib, using MAX instead." maxlog=1
        end
        if nanOpt == CUDNN_NOT_PROPAGATE_NAN
            @warn "Pool nanOpt=$nanOpt not yet implemented in NNlib, using PROPAGATE instead. See https://github.com/FluxML/NNlib.jl/issues/218" maxlog=1
        end
        =#
        ay1 = (mode == CUDNN_POOLING_MAX ? NNlib.maxpool(ax, pdims) :
               mode == CUDNN_POOLING_AVERAGE_COUNT_INCLUDE_PADDING ? NNlib.meanpool(ax, pdims) :
               mode == CUDNN_POOLING_AVERAGE_COUNT_EXCLUDE_PADDING ? NNlib.meanpool(ax, pdims) :
               mode == CUDNN_POOLING_MAX_DETERMINISTIC ? NNlib.maxpool(ax, pdims) :
               error("mode=$mode is not supported."))
        ay1 = alpha * ay1
        ay  = randn!(similar(ay1))
        ay2 = ay1 .+ beta * ay
        d = cudnnPoolingDescriptor(mode, nanOpt, Cint(max(2,ndims(ax)-2)), pooldims(window,size(ax)), pooldims(padding,size(ax)), pooldims(stride,size(ax)))
        nhwc(a) = permutedims(a,(3,1,2,4))
        if format === CUDNN_TENSOR_NCHW
            cx, cy = CuArray.((ax, ay))
        else
            cx, cy = CuArray.(nhwc.((ax,ay)))
            ay1, ay2 = nhwc.((ay1, ay2))
        end
        @test ay1 ≈ cudnnPoolingForward(cx; mode, nanOpt, window, padding, stride, format, alpha) |> Array
        @test ay1 ≈ cudnnPoolingForward(cx, d; format, alpha) |> Array
        @test ay2 ≈ cudnnPoolingForward!(copy(cy), cx; mode, nanOpt, window, padding, stride, format, alpha, beta) |> Array
        @test ay2 ≈ cudnnPoolingForward!(copy(cy), cx, d; format, alpha, beta) |> Array
    end

    expand(::Val{N}, i::NTuple{N}) where {N} = i
    expand(::Val{N}, i::Integer) where {N} = ntuple(_ -> i, N)

    pooltest()
    pooltest(mode = CUDNN_POOLING_AVERAGE_COUNT_INCLUDE_PADDING)
    pooltest(mode = CUDNN_POOLING_AVERAGE_COUNT_EXCLUDE_PADDING)
    pooltest(mode = CUDNN_POOLING_MAX_DETERMINISTIC)
    pooltest(nanOpt = CUDNN_PROPAGATE_NAN)
    pooltest(window = 3)
    pooltest(padding = 1)
    pooltest(stride = 1)
    pooltest(format = CUDNN_TENSOR_NHWC)
    pooltest(dataType = Float16)
    pooltest(alpha = 2)
    pooltest(beta = 2)
end


text:
module mfasttoptest1
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexBeams.CrossSectionModule: CrossSectionRectangle
using FinEtoolsFlexBeams.MeshFrameMemberModule: frame_member, merge_members
using FinEtoolsFlexBeams.FEMMCorotBeamModule
using FinEtoolsFlexBeams.FEMMCorotBeamModule: FEMMCorotBeam
stiffness = FEMMCorotBeamModule.stiffness
mass = FEMMCorotBeamModule.mass
distribloads_global = FEMMCorotBeamModule.distribloads_global
restoringforce = FEMMCorotBeamModule.restoringforce
gyroscopic = FEMMCorotBeamModule.gyroscopic
using FinEtoolsFlexBeams.RotUtilModule: initial_Rfield, linear_update_rotation_field!, update_rotation_field!
using LinearAlgebra: dot
using Arpack
using LinearAlgebra
using SparseArrays
using FinEtoolsFlexBeams.VisUtilModule: plot_points, plot_nodes, plot_midline, render, plot_space_box, plot_solid, space_aspectratio
using PlotlyJS
using JSON
using Test
function test()
    # Parameters:
    E = 71240.0;#MPa
    nu = 0.31;# Poisson ratio
    rho = 2.7e-9;
    Width = 60; Height = 60; Length = 4*Width;; # cross-sectional dimensions
    Mass = Width*Height*Length*rho;
    IL =  (1/12)*Mass*(Width^2+Height^2);
    Ib =  (1/12)*Mass*(Length^2+Height^2)+Mass*(Length/2)^2;
    Ih =  (1/12)*Mass*(Width^2+Length^2)+Mass*(Length/2)^2;
    Omega0 = 313*pi;
    R0 = rotmat3([0.05 0 0]);
    g = 9.81e3; # 9.81 meters per second squared
    q = [0,0,-g*Width*Height*rho];
    utol = 1e-3;
    maxit = 12
    dt = min(2*pi/norm(Omega0)/10, 0.005);
    tend = 0.8;
    ng = 1/2; nb = 1/4*(1/2+ng)^2;
    # Choose the mass formulation:
    mass_type=1;

    # Cross-sectional properties
    cs = CrossSectionRectangle(s -> Width, s -> Width, s -> [1.0, 0.0, 0.0])

    # Select the number of elements per leg.
    spin_vector = R0*[0, 0, Omega0];
    X=[0 0 0;    reshape(R0*[0,0,Length], 1, 3)];
    n=8;
    tolerance=Length/n/100;
    members = []
    push!(members, frame_member(X, n, cs))
    fens, fes = merge_members(members; tolerance = tolerance);

    # Material properties
    material = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    box = fill(0.0, 6)
    initbox!(box, X[1,:])
    supportn = selectnode(fens; box = box, inflate = tolerance)
    initbox!(box, X[2,:])
    tipn = selectnode(fens; box = box, inflate = tolerance)
    for i in [1, 2, 3]
        setebc!(dchi, supportn, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Additional fields
    stepdchi = deepcopy(dchi)
    u1 = deepcopy(u0)
    v0 = deepcopy(dchi) # need the numbering of the degrees of freedom
    for i in 1:size(v0.values, 1)
        v0.values[i, 4:6] = spin_vector
    end
    v1 = deepcopy(dchi)
    a1 = deepcopy(dchi) # zero out the Acceleration
    a0 = deepcopy(dchi) # zero out the Acceleration
    Rfield1 = deepcopy(Rfield0)
    v0v = gathersysvec(v0)
    a0v = gathersysvec(a0)
    vpv = gathersysvec(v0);
    dchipv = gathersysvec(dchi);
    stepdchiv = gathersysvec(dchi);
    rhs = gathersysvec(dchi);
    TMPv = deepcopy(rhs)
    utol = 1e-13*dchi.nfreedofs;

    # tbox = plot_space_box([[-1.1*Width -1.1*Width 0]; [1.1*Width 1.1*Width 1.1*Length]])
    # tshape0 = plot_solid(fens, fes; x = geom0.values, u = 0.0.*dchi.values[:, 1:3], R = Rfield0.values, facecolor = "rgb(125, 155, 125)", opacity = 0.3);
    # plots = cat(tbox,  tshape0; dims = 1)
    # pl = render(plots)
    # sleep(3.5)

    tipx = Float64[]
    tipy = Float64[]
    push!(tipx, X[2,1])
    push!(tipy, X[2,2])
    # tbox = scatter(;x=[0.0, 0.06], y=[-0.06, 0.02], mode="markers")
    # plots = cat(tbox, scatter(;x=tipx./Length, y=tipy./Length, mode="markers+lines"); dims = 1)
    # layout = Layout(width=500, height=500)
    # pl = plot(plots, layout)
    # display(pl)

    femm = FEMMCorotBeam(IntegDomain(fes, GaussRule(1, 2)), material)
    fi = ForceIntensity(q);

    t = 0.0; #
    step = 0;
    while (t <= tend)
        t = t + dt;
        # println("Time $(t)"); # pause
        # Initialization
        applyebc!(dchi) # Apply boundary conditions
        u1.values[:] = u0.values[:]; # guess
        Rfield1.values[:] = Rfield0.values[:]; # guess
        stepdchi.values[:] .= 0.0;# Total increment in current step
        a1.values[:] = -(1/nb/dt)*v0.values[:] -(1/2-nb)/nb*a0.values[:];
        v1.values[:] = v0.values[:] + dt*((1-ng)*a0.values[:] + ng*a1.values[:]);
        gathersysvec!(v0, v0v)
        gathersysvec!(a0, a0v)
        dchipv = dt*v0v + (dt^2/2*(1-2*nb))*a0v
        vpv = v0v +(dt*(1-ng))*a0v;

        iter = 1;
        while true
            F = distribloads_global(femm, geom0, u1, Rfield1, dchi, fi)
            Fr = restoringforce(femm, geom0, u1, Rfield1, dchi);       # Internal forces
            @. rhs = F + Fr;
            K = stiffness(femm, geom0, u1, Rfield1, dchi);
            M = mass(femm, geom0, u1, Rfield1, dchi);
            G = gyroscopic(femm, geom0, u1, Rfield1, v1, dchi);
            gathersysvec!(stepdchi, stepdchiv)
            @. TMPv = ((-1/(nb*dt^2))*stepdchiv+(1/(nb*dt^2))*dchipv)
            rhs .+= M*TMPv
            @. TMPv = ((-ng/nb/dt)*stepdchiv+(ng/nb/dt)*dchipv - vpv)
            rhs .+= G*TMPv;
            dchi = scattersysvec!(dchi, (K+(ng/nb/dt)*G+(1/(nb*dt^2))*M)\rhs); # Disp. incr
            u1.values[:] += (dchi.values[:,1:3])[:];   # increment displacement
            stepdchi.values[:] += dchi.values[:]
            v1.values[:] += (ng/nb/dt)*dchi.values[:];
            a1.values[:] += (1/nb/dt^2)*dchi.values[:];
            update_rotation_field!(Rfield1, dchi)
            if maximum(abs.(dchi.values[:])) < utol# convergence check
                break;
            end
            if (iter > maxit)# bailout for failed convergence
                error("Possible failed convergence");
            end
            iter += 1;
        end
        u0.values[:] = u1.values[:];       # update the displacement
        Rfield0.values[:] = Rfield1.values[:]; # update the rotations
        v0.values[:] = v1.values[:];       # update the velocities
        a0.values[:] = a1.values[:];       # update the accelerations

        # tshape1 = plot_solid(fens, fes; x = geom0.values, u = u1.values, R = Rfield1.values, facecolor = "rgb(125, 15, 15)");
        # plots = cat(tbox,  tshape0,  tshape1; dims = 1)
        # react!(pl, plots, pl.plot.layout)

        if (mod(step,50)==0)
            push!(tipx, X[2,1]+u1.values[tipn[1], 1])
            push!(tipy, X[2,2]+u1.values[tipn[1], 2])
            # plots = cat(tbox, scatter(;x=tipx./Length, y=tipy./Length, mode="markers+lines"); dims = 1)
            # react!(pl, plots, pl.plot.layout)
            # sleep(0.01)
        end

        step=step+1;
    end

    # # Visualize vibration modes
    # scattersysvec!(dchi, v[:, 1])
    # update_rotation_field!(Rfield0, dchi)
    # plots = cat(plot_nodes(fens),
    #     plot_solid(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values);
    #     dims = 1)
    # render(plots; aspectratio = space_aspectratio(fens.xyz))
    reftipx = [0.0, 7.436031085824281e-7, 0.11777019275816547, 0.8203709506659435, 2.241864039281919, 3.9679895248072063, 5.3454217754301885, 5.971372979336259, 5.994274006124595, 5.998838219173172, 6.562336269473291, 7.816659235667364, 9.346281101346545, 10.499101141896105, 10.875849911112233, 10.628929261212697, 10.345322488144156, 10.602025422687818, 11.529299245601559, 12.710964257026657, 13.496913529329909, 13.492271339298975, 12.854269845420365, 12.172880227196368, 12.025322964890227, 12.540504127222404, 13.301744975038586]
    reftipy = [-11.995000624962799, -11.99507329138783, -12.34434321937684, -13.092699071921164, -13.606709453527161, -13.384277801594491, -12.433899377302394, -11.269034682594993, -10.521805428897146, -10.467440913151172, -10.818091285108137, -10.941387824929802, -10.338997425035933, -9.024387391352876, -7.51482956095552, -6.443555337060203, -6.084606620045228, -6.148711597765479, -6.004002866512018, -5.155146068191021, -3.6197117072489498, -1.9179359415572481, -0.6834506841073384, -0.18854220120395482, -0.1419789099985156, 0.08825549466116911, 0.9953416189668296]
    @test norm(reftipx - tipx) / norm(reftipx) <= 1.0e-3
    @test norm(reftipy - tipy) / norm(reftipy) <= 1.0e-3
    return true
end # test
end
using .mfasttoptest1
mfasttoptest1.test()


text:
"""
    mutable struct TensorParser

Define the `TensorParser` object which include the fields:
- `preprocessors::Vector{Any}`: a list of all preprocessing steps
- `contractiontreebuilder::Any`: determine a contraction tree for a contraction involving
                                    multiple tensors
- `contractiontreesorter::Any`: transforms the contraction expression into an expression of
                                    nested binary contractions using the tree output from
                                    the contractiontreebuilder
- `postprocessors::Vector{Any}`: a list of all postprocessing steps
"""
mutable struct TensorParser
    preprocessors::Vector{Any}
    contractiontreebuilder::Any
    contractiontreesorter::Any
    postprocessors::Vector{Any}
    function TensorParser()
        preprocessors = [normalizeindices, expandconj, nconindexcompletion,
                            extracttensorobjects]
        contractiontreebuilder = defaulttreebuilder
        contractiontreesorter = defaulttreesorter
        postprocessors = [_flatten, removelinenumbernode, addtensoroperations]
        return new(preprocessors, contractiontreebuilder, contractiontreesorter,
                    postprocessors)
    end
end

"""
    (parser::TensorParser)(ex::Expr)

Do the tensor contractions with the steps defined in `TensorParser`.
"""
function (parser::TensorParser)(ex::Expr)
    if ex.head == :function
        return Expr(:function, ex.args[1], parser(ex.args[2]))
    end
    for p in parser.preprocessors
        ex = p(ex)::Expr
    end
    treebuilder = parser.contractiontreebuilder
    treesorter = parser.contractiontreesorter
    ex = processcontractions(ex, treebuilder, treesorter)::Expr
    ex = tensorify(ex)::Expr
    for p in parser.postprocessors
        ex = p(ex)::Expr
    end
    return ex
end

"""
    defaulttreebuilder(network)

The input `network` is a list of the untraced indices of each tensorexpression in a tensor
contraction expression, and the `length(network)` is the number of tensorexpressions to be
contracted. Return `tree` in the form e.g. `Any[Any[Any[Any[1, 2], 3], 4], 5]` if
`length(network) == 5`.
"""
function defaulttreebuilder(network)
    if isnconstyle(network)
        tree = ncontree(network)
    else
        tree = Any[1,2]
        for k = 3:length(network)
            tree = Any[tree, k]
        end
    end
    return tree
end

"""
    defaulttreesorter(args, tree)

Sort the sequencee of the contractions as given by the input `tree` (default in the form
`Any[Any[Any[Any[1, 2], 3], 4], 5]` if `length(network) == 5`) by changing the
positions of the tensor expressions to be contracted. The input `args` is a list of the
tensor expressions to be contracted. The returned expression has the form e.g. (default)
`(((A[a,b]*B[c,d])*C[a,c]))*D[b,f])`. This default tree contracts the tensors in the
expression two-by-two from left to right.
"""
function defaulttreesorter(args, tree)
    if isa(tree, Int)
        return args[tree]
    else
        return Expr(:call, :*,
                    defaulttreesorter(args, tree[1]), defaulttreesorter(args, tree[2]))
    end
end

"""
    processcontractions(ex::Expr, treebuilder, treesorter)

Sort the contractions based on the `treebuilder` and `treesorter` if the number of tensor
expressions to be contracted is larger than two.
"""
function processcontractions(ex::Expr, treebuilder, treesorter)
    if ex.head == :macrocall && ex.args[1] == Symbol("@notensor")
        return ex
    end
    ex = Expr(ex.head, map(e->processcontractions(e, treebuilder, treesorter), ex.args)...)
    if istensorcontraction(ex) && length(ex.args) > 3
        args = ex.args[2:end]
        network = map(getindices, args)
        for a in getallindices(ex)
            count(a in n for n in network) <= 2 ||
                throw(ArgumentError("invalid tensor contraction: $ex"))
        end
        tree = treebuilder(network)
        ex = treesorter(args, tree)
    end
    return ex
end
processcontractions(ex, treebuilder, treesorter) = ex # if `ex` is not an Expr do nothing

"""
    tensorify(ex::Expr)

Functions for parsing and processing tensor expressions. Change the tensor expressions to
the actual functions like `add!`, `trace!` and `contract!` they represent.
"""
function tensorify(ex::Expr)
    if ex.head == :macrocall && ex.args[1] == Symbol("@notensor")
        return ex.args[3]
        # > dump(:(@notensor A))
        # > Expr
        #     head: Symbol macrocall
        #     args: Array{Any}((3,))
        #       1: Symbol @notensor
        #       2: LineNumberNode
        #         line: Int64 1
        #         file: Symbol REPL[29]
        #       3: Symbol A
    end
    # assignment case
    if isassignment(ex) || isdefinition(ex)
        lhs, rhs = getlhs(ex), getrhs(ex)
        if isa(rhs, Expr) && rhs.head == :call && rhs.args[1] == :throw
            return rhs
        end
        if istensor(lhs) && istensorexpr(rhs)
            indices = getindices(rhs)
            if hastraceindices(lhs)
                err = "left hand side of an assignment should have unique indices: $lhs"
                return :(throw(IndexError($err)))
            end
            dst, leftind, rightind = decomposetensor(lhs)
            if Set(vcat(leftind,rightind)) != Set(indices)
                err = "non-matching indices between left and right hand side: $ex"
                return :(throw(IndexError($err)))
            end
            if isassignment(ex)
                if ex.head == :(=)
                    return instantiate(dst, false, rhs, true, leftind, rightind)
                elseif ex.head == :(+=)
                    return instantiate(dst, true, rhs, 1, leftind, rightind)
                else # if ex.head == :(-=)
                    return instantiate(dst, true, rhs, -1, leftind, rightind)
                end
            else # if isdefinition(ex)
                return Expr(:(=), dst, instantiate(nothing, false, rhs, true, leftind,
                                                    rightind, false))
            end
        elseif isassignment(ex) && isscalarexpr(lhs)
            if istensorexpr(rhs) && isempty(getindices(rhs))
                return Expr(ex.head, instantiate_scalar(lhs), Expr(:call, :scalar,
                                    instantiate(nothing, false, rhs, true, [], [], true)))
            elseif isscalarexpr(rhs)
                return Expr(ex.head, instantiate_scalar(lhs), instantiate_scalar(rhs))
            end
        else
            return ex # likely an error
        end
    end
    if ex.head == :block
        return Expr(ex.head, map(tensorify, ex.args)...)
    end
    if ex.head == :for
        return Expr(ex.head, ex.args[1], tensorify(ex.args[2]))
    end
    if ex.head == :function
        return Expr(ex.head, ex.args[1], tensorify(ex.args[2]))
    end
    # constructions of the form: a = @tensor ..., where a is a scalar
    if isscalarexpr(ex)
        return instantiate_scalar(ex)
    end
    if istensorexpr(ex)
        if !isempty(getindices(ex))
            err = "cannot evaluate $ex to a scalar: uncontracted indices"
            return :(throw(IndexError($err)))
        end
        return Expr(:call, :scalar, instantiate(nothing, false, ex, true, [], [], true))
    end
    error("invalid syntax in @tensor macro: $ex")
end
tensorify(ex) = ex


text:
#abstract type State end
#abstract type FiniteBasisState <: State end

add!(v::FiniteBasisState, i) = set!(v, toint(v)+i)
zero!(v::FiniteBasisState) = set!(v, 0)
@inline config(v) = v
rand!(v::State) = rand!(GLOBAL_RNG, v)
flipat!(v::State, i) = flipat!(GLOBAL_RNG, v, i)
flipat_fast!(v::State, i) = flipat_fast!(GLOBAL_RNG, v, i)

init_lut!(s::State, net::NeuralNetwork) = nothing
update_lookup!(s::State, net::NeuralNetwork) = nothing

export NAryState, DoubleState, BinaryState
export local_dimension, spacedimension
export nsites, toint, index, index_to_int, flipped, row, col, config
export add!, zero!
export setat!, set!, set_index!, rand!


text:
"""
    CoulombDerivative(D, Z, ℓ)

Helper operator wrapping the `Derivative` operator in the case of a
Coulomb potential of charge `Z` and centrifugal potential
``ℓ(ℓ+1)/2r²``. Its main use is to correctly apply boundary conditions
at ``r=0``, when materializing derivatives.
"""
struct CoulombDerivative{T,Der<:Derivative{T},U<:Number} <: LazyQuasiMatrix{T}
    D::Der
    Z::U
    ℓ::Int
end

axes(D::CoulombDerivative) = axes(D.D)
==(a::CoulombDerivative, b::CoulombDerivative) = a.D == b.D && a.Z == b.Z && a.ℓ == b.ℓ
copy(D::CoulombDerivative) = CoulombDerivative(copy(D.D), D.Z, D.ℓ)

Base.show(io::IO, D::CoulombDerivative) =
    write(io, "CoulombDerivative($(D.D), Z = $(D.Z), ℓ = $(D.ℓ))")

export CoulombDerivative


text:
#
# Copyright (c) 2015, Intel Corporation
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met:
#
# * Redistributions of source code must retain the above copyright
#      notice, this list of conditions and the following disclaimer.
# * Redistributions in binary form must reproduce the above
#      copyright notice, this list of conditions and the following
#      disclaimer in the documentation and/or other materials provided
#      with the distribution.
# * Neither the name of Intel Corporation nor the names of its
#      contributors may be used to endorse or promote products
#      derived from this software without specific prior written
#      permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
# COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.

#*******************************************************************
# NAME:    Pipeline
#
# PURPOSE: This program tests the efficiency with which point-to-point
#          synchronization can be carried out. It does so by executing
#          a pipelined algorithm on an m*n grid. The first array dimension
#          is distributed among the threads (stripwise decomposition).
#
# USAGE:   The program takes as input the
#          dimensions of the grid, and the number of iterations on the grid
#
#                <progname> <iterations> <m> <n>
#
#          The output consists of diagnostics to make sure the
#          algorithm worked, and of timing statistics.
#
# FUNCTIONS CALLED:
#
#          Other than standard C functions, the following
#          functions are used in this program:
#
# HISTORY: - Written by Rob Van der Wijngaart, February 2009.
#          - Converted to Python by Jeff Hammond, February 2016.
#          - Converted to Julia by Jeff Hammond, June 2016.
#
# *******************************************************************

function iterate_over_grid!(grid, m, n)
    for j = 2:n
        for i = 2:m
            @inbounds grid[i,j] = grid[i-1,j] + grid[i,j-1] - grid[i-1,j-1]
        end
    end
end

function main()
    # ********************************************************************
    # read and test input parameters
    # ********************************************************************

    println("Parallel Research Kernels version ") #, PRKVERSION)
    println("Julia pipeline execution on 2D grid")

    if length(ARGS) != 3
        println("argument count = ", length(ARGS))
        println("Usage: ./synch_p2p <# iterations> <first array dimension> <second array dimension>")
        exit(1)
    end

    argv = map(x->parse(Int64,x),ARGS)

    iterations = argv[1]
    if iterations < 1
        println("ERROR: iterations must be >= 1")
        exit(2)
    end

    m = argv[2]
    if m < 1
        println("ERROR: array dimension must be >= 1")
        exit(3)
    end

    n = argv[3]
    if n < 1
        println("ERROR: array dimension must be >= 1")
        exit(4)
    end

    println("Grid sizes               = ", m, ",", n)
    println("Number of iterations     = ", iterations)

    grid = zeros(Float64,m,n)
    grid[1,1:n] = collect(Float64,0:n-1)
    grid[1:m,1] = collect(Float64,0:m-1)

    # precompile hot functions to smooth performance measurement
    precompile(iterate_over_grid!, (Array{Float64, 2}, Int64, Int64))

    t0 = time_ns()

    for k in 0:iterations
        if k==0
            t0 = time_ns()
        end
        iterate_over_grid!(grid, m, n)

        # copy top right corner value to bottom left corner to create dependency
        grid[1,1] = -grid[m,n]
    end
    t1 = time_ns()

    # convert time from nanoseconds to seconds
    pipeline_time = (t1 - t0) * 1.e-9

    # ********************************************************************
    # ** Analyze and output results.
    # ********************************************************************

    epsilon=1.e-8

    # verify correctness, using top right value
    corner_val = Float64((iterations+1)*(n+m-2))
    if (abs(grid[m,n] - corner_val)/corner_val) < epsilon
        println("Solution validates")
        avgtime = pipeline_time/iterations
        println("Rate (MFlops/s): ", 1.e-6*2*(m)*(n)/avgtime, " Avg time (s): ",avgtime)
    else
        println("ERROR: checksum ", grid[m,n], " does not match verification value ", corner_val)
        exit(9)
    end
end

main()

max_stars_repo_name: hattom/Kernels

max_issues_repo_name: hattom/Kernels

max_forks_repo_name: hattom/Kernels


text:
module PL

using CEnum

include("libpvcam_common.jl")
include("libpvcam_api.jl")

end


text:
function localmax(s)
	Y = zeros(size(s))
	I = {}
	for i=2:length(s)-1
		if s[i]>s[i-1] && s[i]>=s[i+1]
			Y[i] = 1
			push!(I,i)
		end
	end
	return Y,I
end

max_stars_repo_name: barzan/dbseer

max_issues_repo_name: barzan/dbseer

max_forks_repo_name: barzan/dbseer


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
"""
    big!(expr)

replaces `Int` and `Float64` values in `expr` by big ones and returns the result.
"""
big!(expr) = expr
big!(expr::Int) = big(expr)
big!(expr::Float64) = big(expr)
function big!(expr::Expr)
    for i in eachindex(expr.args)
        expr.args[i] = big!(expr.args[i])
    end
    expr
end

"""
    @big(expr)
    @big expr

replaces `Int` and `Float64` values in `expr` by big ones and executes the result.
"""
macro big(expr) big!(expr) end

# %%
@doc big!

# %%
@doc @big

# %%
1.0π

# %%
@big 1.0π

# %%
for n in 20:24
    println("factorial(", n, ") = ", factorial(n))
end

# %%
@big for n in 20:24
    println("factorial(", n, ") = ", factorial(n))
end

# %%
@big @eval begin
    n = 57
    c = Tuple(1/factorial(k) for k in 0:n)
end
evalpoly(1, c)

# %%
@big exp(1)

# %%
typeof(c)

# %%

max_stars_repo_path: 0017/@big macro.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0017/@big macro.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0017/@big macro.jl

max_forks_repo_name: genkuroki/public


text:
##### backup.jl #####

# In this file you'll find backup and reverse backup strategies for the
# Genetic Algorithm.

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

export backup, reverse_backup

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

"""
    backup(f ::IOStream, gene ::IntegerGene)

Writes the current state of a `IntegerGene` structure to a buffer.
"""
function backup(f ::IOStream, gene ::IntegerGene)
    write(f, 'I')
    write(f, Int8(length(gene.value)))
    for i in gene.value
        write(f, i)
    end
    write(f, Float64(gene.lbound), Float64(gene.ubound))
    write(f, Int8(-length(gene.name)), gene.name)
    return nothing
end

"""
    backup(f ::IOStream, gene ::FloatGene)

Writes the current state of a `FloatGene` structure to a buffer.
"""
function backup(f ::IOStream, gene ::FloatGene)
    write(f, 'F')
    write(f, gene.m)
    write(f, Int8(length(gene.value)))
    for i in [:value, :range, :lbound, :ubound]
        for j in getfield(gene, i)
            write(f, j)
        end
    end
    for i in gene.name
        l = Int8(-length(i))
        write(f, l, i)
    end

    return nothing
end

"""
    backup(f ::IOStream, gene ::BinaryGene)

Writes the current state of a `BinaryGene` structure to a buffer.
"""
function backup(f ::IOStream, gene ::BinaryGene)
    write(f, 'B')
    write(f, gene.value)
    write(f, Int8(-length(gene.name)), gene.name)
    return nothing
end

"""
    backup( ngens ::Int64              ,
            tgens ::Int64              ,
            chrom ::Vector{Individual} ,
            file  ::AbstractString     )

Writes number of generations `ngens`, total number of generations `tgens` and the population `chrom` into file `file`. Always writes to folder `backup-files` that is created, if nonexistent, inside the `ga` function.
"""
function backup( ngens ::Int64              ,
                 tgens ::Int64              ,
                 pop   ::Vector{Individual} ,
                 file  ::AbstractString     )
    file = "backup-files/$file"
    psize = length(pop   )
    gsize = length(pop[1])
    open(file, "w") do f
        write(f, ngens, tgens, psize, gsize)
        for i in pop
            for j in i
                backup(f, j)
            end
        end
    end
    return nothing
end

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

"""
    reverse_backup(file ::AbstractString)

Reads backup file `file` and saves the number of generations run, the total number of generations supposed to run and the population into variables for later continuing the Genetic Algorithm.
"""
function reverse_backup(file ::AbstractString)
    f = open(file, "r")

    ngens    = read(f, Int64)
    tgens    = read(f, Int64)
    popsize  = read(f, Int64)
    genesize = read(f, Int64)

    population = Vector{Individual}(undef, popsize)
    for p in 1:popsize
        population[p] = Individual(undef, genesize)
        for g in 1:genesize
            id = read(f, Char)
            if id == 'I'
                nvals   = read(f, Int8)
                bit_vec = BitVector(undef, nvals)
                readbytes!(f, reinterpret(UInt8, bit_vec))
                lb      = read(f, Float64)
                ub      = read(f, Float64)
                strsize = -read(f, Int8)
                name    = Vector{UInt8}(undef, strsize)
                readbytes!(f, name)
                population[p][g] =
                    IntegerGene(bit_vec, lb, ub, String(name))
            elseif id == 'F'
                m      = read(f, Int64)
                nvals  = read(f, Int8 )
                names  = Vector{String }(undef, nvals)
                values = Vector{Float64}(undef, nvals)
                ranges = Vector{Float64}(undef, nvals)
                lb     = Vector{Float64}(undef, nvals)
                ub     = Vector{Float64}(undef, nvals)
                readbytes!(f, reinterpret(UInt8, values))
                readbytes!(f, reinterpret(UInt8, ranges))
                readbytes!(f, reinterpret(UInt8, lb    ))
                readbytes!(f, reinterpret(UInt8, ub    ))
                for i in 1:nvals
                    strsize = -read(f, Int8)
                    name = Vector{UInt8}(undef, strsize)
                    readbytes!(f, name)
                    names[i] = String(name)
                end
                population[p][g] =
                    FloatGene(values, ranges, names, m, lb, ub)
            elseif id == 'B'
                value = read(f, Bool)
                strsize = -read(f, Int8)
                name = Vector{UInt8}(undef, strsize)
                readbytes!(f, name)
                population[p][g] = BinaryGene(value, String(name))
            end
        end
    end

    close(f)
    return ngens, tgens, population
end

"""
    reverse_backup(files ::Vector{<:AbstractString})

This should be used only for backup files of a parallel run, since each worker writes its own backup file.

Reads backup files `files` and returns the number of generations run in the slowest worker, the total number of generations supposed to be run and the entire population.
"""
function reverse_backup(files ::Vector{<:AbstractString})
    tgens = Vector{Int64     }(undef, length(files))
    gens  = Vector{Int64     }(undef, length(files))
    pop   = Vector{Individual}(undef, 0            )

    for (i,f) in enumerate(files)
        ngen, tgen, popul = reverse_backup(f)
        gens[i] = ngen
        tgens[i] = tgen
        append!(pop, popul)
    end

    min_gens = findmin( gens)[1]
    tot_gens = findmax(tgens)[1]

    return min_gens, tot_gens, pop
end


text:
using TypeClasses

# Monoid instances
# ================

# this is standard Applicative combine implementation
# there is no other sensible definition for combine and hence if this does fail, it fails correctly
TypeClasses.combine(a::Callable, b::Callable) = mapn(combine, a, b)

# we do not implement the same for `orelse`, as `orelse` is usually defined on the container level, but there is no sensible default implementation.

# it is not possible to implement neutral, as we do not know the element type without executing the function

# Monad instances
# ===============

# there is no definition for Base.foreach, as a callable is not runnable without knowing the arguments
TypeClasses.map(f, g::Callable) = Callable((args...; kwargs...) -> f(g(args...; kwargs...)))

TypeClasses.pure(::Type{<:Callable}, a) = (args...; kwargs...) -> a
TypeClasses.ap(f::Callable, g::Callable) = Callable((args...; kwargs...) -> f(args...; kwargs...)(g(args...; kwargs...)))

# we cannot overload this generically, because `Base.map(f, ::Vector...)` would get overwritten as well (even without warning surprisingly)
TypeClasses.map(f, a::Callable, b::Callable, more::Callable...) = mapn(f, a, b, more...)

# we don't use convert, but directly use function call,
# which should give readable errors that a something is not callable, and is a bit more flexible
TypeClasses.flatmap(f, g::Callable) = Callable((args...; kwargs...) -> f(g(args...; kwargs...))(args...; kwargs...))

# FlipTypes instance
# ==================

# there cannot be any flip_types for functions


text:
function RLCore.Experiment(
    ::Val{:rlpyt},
    ::Val{:PPO},
    ::Val{:Atari},
    name::AbstractString;
    save_dir = nothing,
    seed = 123,
)
    @warn "Currently setting the `seed` will not guarantee the reproducibility. The instability seems to be caused by the `CrossCor` layer when calculating gradient."
    rng = StableRNG(seed)
    if isnothing(save_dir)
        t = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
        save_dir = joinpath(pwd(), "checkpoints", "rlpyt_PPO_Atari_$(name)_$(t)")
    end

    lg = TBLogger(joinpath(save_dir, "tb_log"), min_level = Logging.Info)

    N_ENV = 32
    UPDATE_FREQ = 64
    N_FRAMES = 4
    STATE_SIZE = (80, 104)
    env = MultiThreadEnv([
        atari_env_factory(
            name,
            STATE_SIZE,
            N_FRAMES;
            repeat_action_probability = 0,
            seed = seed + i,
        ) for i in 1:N_ENV
    ])
    N_ACTIONS = length(action_space(env[1]))
    INIT_CLIP_RANGE = 0.1f0
    INIT_LEARNING_RATE = 1e-3

    init = orthogonal(rng)

    # share model
    model = Chain(
        x -> x ./ 255,
        CrossCor((8, 8), N_FRAMES => 32, relu; stride = 4, pad = 0, init = init),
        CrossCor((4, 4), 32 => 64, relu; stride = 2, pad = 1, init = init),
        x -> reshape(x, :, size(x)[end]),
        Dense(6912, 512, relu; initW = init),
    )

    agent = Agent(
        policy = RandomStartPolicy(
            num_rand_start = 1000,
            random_policy = RandomPolicy(action_space(env); rng = rng),
            policy = PPOPolicy(
                approximator = ActorCritic(
                    actor = Chain(model, Dense(512, N_ACTIONS; initW = init)),
                    critic = Chain(model, Dense(512, 1; initW = init)),
                    optimizer = ADAM(INIT_LEARNING_RATE),  # decrease learning rate with a hook
                ) |> gpu,
                γ = 0.99f0,
                λ = 0.98f0,
                clip_range = INIT_CLIP_RANGE,  # decrease with a hook
                max_grad_norm = 1.0f0,
                n_microbatches = 4,
                n_epochs = 4,
                actor_loss_weight = 1.0f0,
                critic_loss_weight = 0.5f0,
                entropy_loss_weight = 0.01f0,
                rng = rng,
                update_freq = UPDATE_FREQ,
            ),
        ),
        trajectory = PPOTrajectory(;
            capacity = UPDATE_FREQ,
            state = Array{Float32,4} => (STATE_SIZE..., N_FRAMES, N_ENV),
            action = Vector{Int} => (N_ENV,),
            reward = Vector{Float32} => (N_ENV,),
            terminal = Vector{Bool} => (N_ENV,),
            action_log_prob = Vector{Float32} => (N_ENV,),
        ),
    )

    N_TRAINING_STEPS = 50_000_000 ÷ N_ENV
    EVALUATION_FREQ = N_TRAINING_STEPS ÷ 100
    MAX_EPISODE_STEPS_EVAL = 27_000
    N_CHECKPOINTS = 3
    stop_condition = StopAfterStep(N_TRAINING_STEPS)

    total_batch_reward_per_episode = TotalBatchRewardPerEpisode(N_ENV)
    batch_steps_per_episode = BatchStepsPerEpisode(N_ENV)
    evaluation_result = []

    hook = ComposedHook(
        total_batch_reward_per_episode,
        batch_steps_per_episode,
        DoEveryNStep(UPDATE_FREQ) do t, agent, env
            p = agent.policy.policy
            with_logger(lg) do
                @info "training" loss = mean(p.loss) actor_loss = mean(p.actor_loss) critic_loss =
                    mean(p.critic_loss) entropy_loss = mean(p.entropy_loss) norm =
                    mean(p.norm) log_step_increment = UPDATE_FREQ
            end
        end,
        DoEveryNStep(UPDATE_FREQ) do t, agent, env
            decay = (N_TRAINING_STEPS - t) / N_TRAINING_STEPS
            agent.policy.policy.approximator.optimizer.eta = INIT_LEARNING_RATE * decay
            agent.policy.policy.clip_range = INIT_CLIP_RANGE * Float32(decay)
        end,
        DoEveryNStep() do t, agent, env
            with_logger(lg) do
                rewards = [
                    total_batch_reward_per_episode.rewards[i][end]
                    for i in 1:length(env) if is_terminated(env[i])
                ]
                if length(rewards) > 0
                    @info "training" rewards = mean(rewards) log_step_increment = 0
                end
                steps = [
                    batch_steps_per_episode.steps[i][end]
                    for i in 1:length(env) if is_terminated(env[i])
                ]
                if length(steps) > 0
                    @info "training" steps = mean(steps) log_step_increment = 0
                end
            end
        end,
        DoEveryNStep(EVALUATION_FREQ) do t, agent, env
            @info "evaluating agent at $t step..."
            # switch to GreedyExplorer?
            h = TotalBatchRewardPerEpisode(N_ENV)
            s = @elapsed run(
                agent.policy,
                MultiThreadEnv([
                    atari_env_factory(
                        name,
                        STATE_SIZE,
                        N_FRAMES,
                        MAX_EPISODE_STEPS_EVAL;
                        repeat_action_probability = 0,
                        seed = seed + t + i,
                    ) for i in 1:4
                ]),
                StopAfterStep(27_000; is_show_progress = false),
                h,
            )
            res = (avg_score = mean(Iterators.flatten(h.rewards)),)
            push!(evaluation_result, res)

            @info "finished evaluating agent in $s seconds" avg_score = res.avg_score
            with_logger(lg) do
                @info "evaluating" avg_score = res.avg_score log_step_increment = 0
            end

            policy = cpu(agent.policy)
            mkdir(joinpath(save_dir, string(t)))
            BSON.@save joinpath(save_dir, string(t), "policy.bson") policy
            BSON.@save joinpath(save_dir, string(t), "stats.bson") total_batch_reward_per_episode evaluation_result

            # only keep recent 3 checkpoints
            old_checkpoint_folder =
                joinpath(save_dir, string(t - EVALUATION_FREQ * N_CHECKPOINTS))
            if isdir(old_checkpoint_folder)
                rm(old_checkpoint_folder; force = true, recursive = true)
            end
        end,
    )

    description = """
    # Play Atari($name) with PPO
    The configurations use here should be aligned with [atari_ff_ppo.py](https://github.com/astooke/rlpyt/blob/master/rlpyt/experiments/configs/atari/pg/atari_ff_ppo.py) in rlpyt.
    You can also view the tensorboard logs with `tensorboard --logdir $(joinpath(save_dir, "tb_log"))`
    """

    Experiment(agent, env, stop_condition, hook, description)
end


text:
function HyperRectangle(square::Square{T}) where {T}
    HyperRectangle{2,T}(square.lowercorner, square.dimensions)
end

function HyperRectangle(cube::Cuboid{T}) where {T}
    HyperRectangle{3,T}(cube.lowercorner, cube.dimensions)
end

function HyperRectangle(sphere::Sphere{T}) where {T}
    HyperRectangle{3,T}(fill(-sphere.radius,3), fill(sphere.radius*2,3))
end

function HyperRectangle(p::Cylinder{T}) where {T}
    HyperRectangle{3,T}(Vec(-p.radius,-p.radius,0), Vec(p.radius*2,p.radius*2,p.height))
end

function HyperRectangle(p::Circle{T}) where {T}
    HyperRectangle{2,T}(Vec(-p.radius,-p.radius), Vec(p.radius*2,p.radius*2))
end

function HyperRectangle(p::Piping{T}) where {T}
    maxx, maxy, maxz = typemin(Float64), typemin(Float64), typemin(Float64)
    minx, miny, minz = typemax(Float64), typemax(Float64), typemax(Float64)
    for pt in p.points
        newx, newy, newz = pt
        maxx = max(newx,maxx)
        maxy = max(newy,maxy)
        maxz = max(newz,maxz)
        minx = min(newx,minx)
        miny = min(newy,miny)
        minz = min(newz,minz)
    end
    minv = Vec(minx,miny,minz) .- p.radius
    maxv = Vec(maxx,maxy,maxz) .+ p.radius
    HyperRectangle(minv, maxv-minv)
end

function HyperRectangle(m::MapContainer)
    transform(m.map, HyperRectangle(m.primitive))
end

function HyperRectangle(csg::CSGUnion)
    union(HyperRectangle(csg.left), HyperRectangle(csg.right))
end

function HyperRectangle(csg::RadiusedCSGUnion)
    union(HyperRectangle(csg.left), HyperRectangle(csg.right))
end

function HyperRectangle(csg::CSGIntersect)
    intersect(HyperRectangle(csg.left), HyperRectangle(csg.right))
end

function HyperRectangle(csg::CSGDiff)
    diff(HyperRectangle(csg.left), HyperRectangle(csg.right))
end

function HyperRectangle(s::Shell)
    HyperRectangle(s.primitive)
end

function HyperRectangle(s::LinearExtrude)
    h = HyperRectangle(s.primitive)
    HyperRectangle(Vec(h.origin[1], h.origin[2], 0), Vec(h.widths[1], h.widths[2], s.distance))
end

"""
Perform a union between two HyperRectangles.
"""
function union(h1::HyperRectangle, h2::HyperRectangle)
    m = min.(minimum(h1), minimum(h2))
    mm = max.(maximum(h1), maximum(h2))
    HyperRectangle(m, mm - m)
end

"""
Perform a difference between two HyperRectangles.
"""
diff(h1::HyperRectangle, h2::HyperRectangle) = h1

"""
Perform a intersection between two HyperRectangles.
"""
function intersect(h1::HyperRectangle, h2::HyperRectangle)
    m = max.(minimum(h1), minimum(h2))
    mm = min.(maximum(h1), maximum(h2))
    HyperRectangle(m, mm - m)
end


text:
module ScatteredInterpolation

using Distances, NearestNeighbors, Combinatorics, LinearAlgebra

export interpolate,
    evaluate

abstract type ScatteredInterpolant end
abstract type InterpolationMethod end

include("./rbf.jl")
include("./idw.jl")
include("./nearestNeighbor.jl")

# Fallback method for the case of just one point
function evaluate(itp::ScatteredInterpolant, points::AbstractArray{<:Real, 1})

    # pairwise requires the points array to be 2-d.
    n = length(points)
    points = reshape(points, n, 1)

    evaluate(itp, points)
end

"""
    interpolate(method, points, samples; metric = Euclidean(), returnRBFmatrix = false; smooth = false)

Create an interpolation of the data in `samples` sampled at the locations defined in
`points` based on the interpolation method `method`. `metric` is any of the metrics defined
by the `Distances` package. The RBF matrix used for solving the weights can be returned with
the boolean `returnRBFmatrix`. Note that this option is only valid for RadialBasisFunction
interpolations.

`points` should be an ``n×k`` matrix, where ``n`` is dimension of the sampled space and
``k`` is the number of points. This means that each column in the matrix defines one point.

`samples` is an ``k×m`` array, where ``k`` is the number of sampled points (same as for
`points`) and ``m`` is the dimension of the sampled data.

The RadialBasisFunction interpolation supports the use of unique RBF functions and widths
for each sampled point by supplying `method` with a vector of interpolation methods
of length ``k``.

The RadialBasisFunction interpolation also supports smoothing of the data points using
ridge regression. All points can be smoothed equally supplying a scalar value,
alternatively each point can be smoothed independently by supplying a vector of smoothing
values. Note that it is no longer interpolating when using smoothing.

The returned `ScatteredInterpolant` object can be passed to `evaluate` to interpolate the
data to new points.
"""
function interpolate end

"""
    evaluate(itp, points)

Evaluate an interpolation object `itp` at the locations defined in `points`.

`points` should be an ``n×k`` matrix, where ``n`` is dimension of the sampled space and
``k`` is the number of points. This means that each column in the matrix defines one point.
"""
function evaluate end
end


text:
using GLWindow, GLAbstraction, ImmutableArrays, GLFW, React, ModernGL, Quaternions, GLText, GLPlot, Images

window 	= createwindow("Mesh Display", 1000, 1000, debugging = false) # debugging just works on linux and windows
cam 	= PerspectiveCamera(window.inputs, Vec3(1,1,1), Vec3(0))

initplotting()

function textwithoffset(start::Vec3, text, advance, lineHeight)
	resulttext 	= Vec1[]
	offset 		= Vec3[start]
	for elem in text
		offset1 = last(offset)
		if elem == '\t'
        	offset[end] = Vec3(offset1[1] + advance * 3, offset1[2], 0f0)
        elseif  elem == ' '
        	offset[end] = Vec3(offset1[1] + advance, offset1[2], 0f0)
        elseif elem == '\r' || elem == '\n'
        	offset[end] = Vec3(start[1], offset1[2] - lineHeight, 0f0)
        else
        	glchar = float32(elem)
        	if glchar > 256
        		glchar = float32(0)
        	end
			push!(offset, Vec3(offset1[1] + advance, offset1[2], 0f0))
			push!(resulttext, Vector1(glchar))
        end
	end
	offset, resulttext
end

text = "ööähh Hallo, was ist denn hier los?#üa+#üsadi0ß24 \nakjskljd ="
font = GLFont("VeraMono")
offset, ctext 	= textwithoffset(Vec3(0), text, font.props[1]...)
println(offset)

parameters 		= [(GL_TEXTURE_WRAP_S,  GL_CLAMP_TO_EDGE),(GL_TEXTURE_MIN_FILTER, GL_NEAREST)]
texttex 		= Texture(ctext, parameters=parameters)
offset 			= Texture(offset, parameters=parameters)

view = [
  "GLSL_EXTENSIONS"     	=> "#extension GL_ARB_draw_instanced : enable",
  "offset_calculation"     => "texelFetch(offset, index, 0).rgb;",
]

data = merge([
	:index_offset		=> convert(GLint, 0),
    :rotation 			=> Vec4(1,0,0,0),
    :text 				=> texttex,
    :offset 			=> offset,
    :scale 				=> Vec2(1/500),
    :color 				=> Vec4(0,0,0,1),
    :backgroundcolor 	=> Vec4(0),
    :projectionview 	=> cam.projectionview
], font.data)

program = TemplateProgram("textShader.vert", "textShader.frag", view, data)
obj = instancedobject(data, program, length(text))
#prerender!(obj, glEnable, GL_DEPTH_TEST)
obj2 = toopengl(Texture("VeraMono.bmp"))

glClearColor(1,1,1,0)
while !GLFW.WindowShouldClose(window.glfwWindow)
  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)

  render(GRID)
  render(obj)
  #render(obj2)

  yield() # this is needed for react to work
  GLFW.SwapBuffers(window.glfwWindow)
  GLFW.PollEvents()
end
GLFW.Terminate()