text:
################################################
#           jtwalk
#   Twalk in Julia.
#   Twalk by Andres Christen & Colin Fox
#   Julia implementation by Nicolas Kuschinski
#
################################################

#   This is able to run the twalk algorithm, but:
#     -does not estimate remaining time
#     -does not perform output analysis
#       - does not include a RWMH.
#   Run and onemove methods are exported.
#   And can be accessed independently.

# This is mostly just translating the python version of twalk,
# some optimization has been done for Julia:
# Type instability has been resolved everywhere except in Energy! SetupInitialValues! and onemove.
# SetupInitialValues! and Energy! do not matter because they are only used once.

# Type instability in onemove has been isolated to the invocation of the Supp and U functions.
# Profiling and analysis done by José Zubieta indicates that this type instability is not serious
# and that the performance hit is minimal

###############################################
# NOTE:
# This uses nothing outside of basic julia
# therefore no external libraries are required.
###############################################

module JTwalk
    export Twalk
    export jtwalk
    export Run!
    export onemove

#Auxilliary function:

    function SqrNorm(x::Array{Float64,1})
        sum(x .* x)
    end

    #Twalk structure
    mutable struct Twalk
        n::Int64
        U::Function
        Supp::Function
        t::Int64
        u::Float64
        up::Float64
        aw::Float64
        at::Float64
        pphi::Float64
        Output::Array{Float64,2}
        Output_u::Array{Float64,2}
        T::Int64
        Fw::Array{Float64,1}
        LikelihoodEnergy::Function
        PriorEnergy::Function
        ll_e::Float64
        prior_e::Float64
        nphi::Int64
        phi::Array{Bool,1}
        sigma::Float64
        x::Array{Float64,1}
        xp::Array{Float64,1}
    end

    function Energy!(T::Twalk, x::Array{Float64,1})
        T.ll_e=T.LikelihoodEnergy(x)
        T.prior_e=T.PriorEnergy(x)
        T.t*T.ll_e+T.prior_e
    end

    #Initializer All arguments are keywords
    function jtwalk(;n::Int64, U=x->x[1]^2/2, Supp=x->true, t=-1, u=0, up=0, w=x->0., ww=[0., 0.4918, 0.4918, 0.0082, 0.0082], aw=1.5, at=6., n1phi=4.)
        if t>0  #Penalized likelihood
            LikelihoodEnergy=u
            PriorEnergy=w
            selft=t
        else
            PriorEnergy=x->0
            LikelihoodEnergy=U
            selft=1
        end
        pphi=min(n,n1phi)/n ##Prob of choosing each parameter
        #
        #Create the object first & then modify it a little. Matches python version more closely.
        twalkinitial=Twalk(n, U, Supp, t, 0., up, aw, at, pphi, zeros(2,2), zeros(2,2), 1, cumsum(ww), LikelihoodEnergy, PriorEnergy, 0., 0., 0,[true],0.,[0],[0])
        #
        twalkinitial.U=x->Energy!(twalkinitial,x)
        #
        twalkinitial
    end

    function SetupInitialValues!(twalk::Twalk, x0::Array{Float64,1}, xp0::Array{Float64,1})
        #Check if x0 & xp0 are in the support
        (any(abs.(x0-xp0) .<=0)) && error("Not all entires of the initial values are different")
        (twalk.Supp(x0)) || error("x0 out of support")
        (twalk.Supp(xp0)) || error("xp0 out of support")
        u=twalk.U(x0)
        up=twalk.U(xp0)
        [true u up]
    end

    ###############################
    #Run the twalk

    function Run!(twalk::Twalk;T::Int64,x0::Array{Float64,1},xp0::Array{Float64,1},t=1)
        #Run the twalk. T=number of iterations, x0 & xp0 are start points.
        twalk.t=t
        #
        (rt, u, up)=SetupInitialValues!(twalk,x0,xp0)
        Bool(rt) || return 0
        x=x0 #last values used
        xp=xp0
        twalk.Output = zeros(T+1,twalk.n+1)
        twalk.Output_u=zeros(T+1,twalk.n+1)
        twalk.T = T+1
        kercall = zeros(6)
        #
        twalk.Output[1, 1:twalk.n] = copy(x)
        twalk.Output[1, twalk.n+1]=u
        twalk.Output_u[1]=twalk.ll_e
        #
        j1=1
        j=0
        #
        #Sampling
        for it in 1:T
            y, yp, ke, A, u_prop, up_prop = onemove(twalk, x, u, xp, up)
            kercall[ke]+=1
            kercall[6]+=1
            if rand()<A
                x=copy(y)
                u=u_prop
                xp=copy(yp)
                up=up_prop
            end
            #
            #To retrieve current values
            twalk.x=x
            twalk.xp=xp
            twalk.u=u
            twalk.up=up
            #
            twalk.Output[it+1,1:twalk.n]=copy(x)
            twalk.Output[it+1,twalk.n+1]=u
            twalk.Output_u[it+1]=twalk.ll_e
        end
        #
        return 1
    end

    ###############################
    #One move of twalk

    function onemove(twalk::Twalk, x::Array{Float64,1}, u::Float64, xp::Array{Float64,1}, up::Float64)
        #One move of the twalk
        n=twalk.n
        U=twalk.U
        Supp=twalk.Supp
        Fw=twalk.Fw
        ker=rand()
        ke=1
        u_prop::Float64=0.0
        up_prop::Float64=0.0
        A::Float64=0.0
        #
        #kernel nothing exchange x with xp, ~ used
        if ((0<=ker) & (ker<Fw[1]))
            ke=0
            y=copy(xp)
            up_prop=u
            yp=copy(x)
            u_prop=up
            A=1.0
        end
        #
        #The walk move
        if (Fw[1] <= ker < Fw[2])
            ke=1
            dir=rand()
            if(dir<0.5) #x as pivot
                yp=SimWalk(twalk,xp,x)
                y=copy(x)
                u_prop=u
                if (Supp(yp)& all(yp.!=y))
                    up_prop=U(yp)
                    A=exp(up-up_prop)
                else
                    up_prop=0.0
                    A=0.0
                end
            else #xp as pivot
                y=SimWalk(twalk,x,xp)
                yp=copy(xp)
                up_prop=up
                if(Supp(y) & (all(yp.!=y)))
                    u_prop=U(y)
                    A=exp(u-u_prop)
                else
                    u_prop=0.0
                    A=0.0
                end
            end
        end
        #
        #The traverse move
        if (Fw[2] <= ker < Fw[3])
            ke=2
            dir=rand()
            if(dir<0.5) #x as pivot
                beta=Simbeta(twalk)
                yp=SimTraverse(twalk,xp,x,beta)
                y=copy(x)
                u_prop=u
                if Supp(yp)
                    up_prop=U(yp)
                    (twalk.nphi==0) ? A=1.0 : A=(exp((up-up_prop) + (twalk.nphi-2)*log(beta)))
                else
                    up_prop=0.0
                    A=0.0
                end
            else #xp as pivot
                beta=Simbeta(twalk)
                y=SimTraverse(twalk,x,xp,beta)
                yp=copy(xp)
                up_prop=up
                if Supp(y)
                    u_prop=U(y)
                    (twalk.nphi==0) ? A=1.0 : A=(exp((u-u_prop) + (twalk.nphi-2)*log(beta)))
                else
                    u_prop=0.0
                    A=0.0
                end
            end
        end
        #
        #The blow move
        if (Fw[3] <= ker < Fw[4])
            ke=3
            dir=rand()
            if(dir<0.5) #x as pivot
                yp=SimBlow(twalk,xp,x)
                y=copy(x)
                u_prop=u
                if (Supp(yp)& all(yp.!=x))
                    up_prop=U(yp)
                    W1=GBlowU(twalk,yp,xp,x)
                    W2=GBlowU(twalk,xp,yp,x)
                    A=exp((up-up_prop)+(W1-W2))
                else
                    up_prop=0.0
                    A=0.0
                end
            else #xp as pivot
                y=SimBlow(twalk,x,xp)
                yp=copy(xp)
                up_prop=up
                if(Supp(y) & (all(y.!=xp)))
                    u_prop=U(y)
                    W1=GBlowU(twalk,y,x,xp)
                    W2=GBlowU(twalk,x,y,xp)
                    A=exp((u-u_prop)+(W1-W2))
                else
                    u_prop=0.0
                    A=0.0
                end
            end
        end
        #
        #The Hop move
        if (Fw[4] <= ker < Fw[5])
            ke=4
            dir=rand()
            if(dir<0.5) #x as pivot
                yp=SimHop(twalk,xp,x)
                y=copy(x)
                u_prop=u
                if (Supp(yp)& all(yp.!=x))
                    up_prop=U(yp)
                    W1=GHopU(twalk,yp,xp,x)
                    W2=GHopU(twalk,xp,yp,x)
                    A=exp((up-up_prop)+(W1-W2))
                else
                    up_prop=0.0
                    A=0.0
                end
            else #xp as pivot
                y=SimHop(twalk,x,xp)
                yp=copy(xp)
                up_prop=up
                if(Supp(y) & (all(y.!=xp)))
                    u_prop=U(y)
                    W1=GHopU(twalk,y,x,xp)
                    W2=GHopU(twalk,x,y,xp)
                    A=exp((u-u_prop)+(W1-W2))
                else
                    u_prop=0.0
                    A=0.0
                end
            end
        end
        return [y,yp,ke,A,u_prop,up_prop]
    end

    ###############################
    #Auxiliary functions for the kernels

    #used by the walk kernel
    function SimWalk(twalk::Twalk,x::Array{Float64,1},xp::Array{Float64,1})
        aw=twalk.aw
        n=twalk.n
        phi=(rand(n).<twalk.pphi)
        twalk.nphi=sum(phi)
        z=zeros(n)
        for i in 1:n
            if phi[i]
                u=rand()
                z[i]=(aw/(1+aw))*(aw*u^2+2*u-1)
            end
        end
        return x+((x-xp) .* z)
    end

    #used by the traverse kernel
    function Simbeta(twalk::Twalk)
        at=twalk.at
        if rand()<((at+1)/(2*at))
            return exp(1/(at+1)*log(rand()))
        else
            return exp(1/(1-at)*log(rand()))
        end
    end

    function SimTraverse(twalk, x::Array{Float64,1}, xp::Array{Float64,1}, beta::Float64)
        n=twalk.n
        phi=(rand(n).<twalk.pphi)
        twalk.nphi=sum(phi)
        rt=copy(x)
        for i in 1:n
            if phi[i]
                rt[i]=xp[i]+beta*(xp[i]-x[i])
            end
        end
        return rt
    end

    #used by the blow kernel
    function SimBlow(twalk::Twalk,x::Array{Float64,1},xp::Array{Float64,1})
        n=twalk.n
        twalk.phi=rand(n).<twalk.pphi
        twalk.nphi=sum(twalk.phi)
        twalk.sigma=maximum(twalk.phi .* abs.(xp-x))
        rt=copy(x)
        for i in 1:n
            if twalk.phi[i]
                rt[i]=xp[i]+twalk.sigma*randn()
            end
        end
        return rt
    end

    function GBlowU(twalk::Twalk,h,x::Array{Float64,1},xp::Array{Float64,1})
        nphi=twalk.nphi
        twalk.sigma=maximum(twalk.phi .* abs.(xp-x))
        if (nphi>0)
            return(nphi/2.0)*log(2*pi) + nphi*log(twalk.sigma)+0.5*SqrNorm(h-xp)/(twalk.sigma^2)
        else
            return 0.0
        end
    end

    #used by the hop kernel
    function SimHop(twalk::Twalk,x::Array{Float64,1},xp::Array{Float64,1})
        n=twalk.n
        twalk.phi=rand(n).<twalk.pphi
        twalk.nphi=sum(twalk.phi)
        twalk.sigma=maximum(twalk.phi .* abs.(xp-x))/3
        rt=copy(x)
        for i in 1:n
            if twalk.phi[i]
                rt[i]=x[i]+twalk.sigma*randn()
            end
        end
        return rt
    end

    function GHopU(twalk::Twalk,h,x::Array{Float64,1},xp::Array{Float64,1})
        nphi=twalk.nphi
        twalk.sigma=maximum(twalk.phi .* abs.(xp-x))/3
        if (nphi>0)
            return(nphi/2.0)*log(2*pi) + nphi*log(twalk.sigma)+0.5*SqrNorm(h-x)/(twalk.sigma^2)
        else
            return 0.0
        end
    end

end #module

max_stars_repo_name: tuerda/JTwalk

max_issues_repo_name: tuerda/JTwalk

max_forks_repo_name: tuerda/JTwalk


text:
module libclang
	using CBinding
	using Clang_jll

	c`$([
		Main.SYSROOT...,

		"-fparse-all-comments", "-fno-blocks",
		"-I$(Clang_jll.artifact_dir)/include",
		"-L$(dirname(Clang_jll.libclang_path))", "-lclang",
	])`

	const c"time_t" = Cvoid  # what is time_t?
	const c"size_t" = Csize_t

	c"""
		#include <clang-c/Documentation.h>
		#include <clang-c/Index.h>
	"""ji
end

ext: jl


text:
module lab5

    include("kepler_eqn.jl")
    include("calc_rv.jl")
    include("periodogram.jl")
    #include("periodogram_orig.jl")
end


text:
using CairoMakie, Random, Distributions
Random.seed!(1234)
b = Binomial(10, 0.85)
n = 500
function someSegments(; n = 50)
    Point2f0.(vec([[x, rand(b)] for i in 1:2, x = rand(n)]))
end
with_theme(theme_dark()) do
    fig, = linesegments(someSegments(; n = n); color = rand(n),
        colormap = :Spectral_11, linewidth = abs.(randn(n)),
        axis = (; xlabel = "x", ylabel = "y"),
        figure = (; resolution = (600,400)))
    save("linesegments.png", fig, px_per_unit = 2)
    display(fig)
end


text:
module TestMarkdown

using Test
using Markdown
using Comonicon.AST
using Comonicon
using Comonicon:
    JLArgument,
    JLOption,
    JLFlag,
    JLMD,
    JLMDFlag,
    JLMDOption,
    cast,
    cast_args,
    cast_flags,
    cast_options,
    default_name,
    get_version,
    split_leaf_command,
    split_docstring,
    read_arguments,
    read_intro,
    read_description,
    read_options,
    read_flags,
    split_hint,
    split_option

@testset "read_arguments" begin
    doc = Markdown.parse("""
        command <args>

    description of the command.

    # Intro

    asdmwds dasklsam xasdklqm dasdm, qwdjiolkasjdsa
    dasklmdas weqwlkjmdas kljnsadlksad qwlkdnasd
    dasklmdlqwoi, dasdasklmd qw,asd. dasdjklnmldqw.

    # Args

    - `arg1`: argument 1.
    - `arg2`: argument 2.
    - `arg3`: argument 3.
    - `arg4`: argument 4.
    """)

    @test read_description(doc) == "description of the command."
    @test read_intro(doc) ==
          "asdmwds dasklsam xasdklqm dasdm, qwdjiolkasjdsa " *
          "dasklmdas weqwlkjmdas kljnsadlksad qwlkdnasd " *
          "dasklmdlqwoi, dasdasklmd qw,asd. dasdjklnmldqw."

    args = read_arguments(doc)
    @test args["arg1"] == "argument 1."
    @test args["arg2"] == "argument 2."
    @test args["arg3"] == "argument 3."
    @test args["arg4"] == "argument 4."

    doc = Markdown.parse("""
    description of the command.

    # Arguments

    - `arg1`: argument 1.
    - `arg2`: argument 2.
    - `arg3`: argument 3.
    - `arg4`: argument 4.
    """)

    @test read_description(doc) == "description of the command."

    args = read_arguments(doc)
    @test args["arg1"] == "argument 1."
    @test args["arg2"] == "argument 2."
    @test args["arg3"] == "argument 3."
    @test args["arg4"] == "argument 4."
end

@testset "split_hint" begin
    @test split_hint("long") == ("long", nothing)
    @test split_hint("long <value>") == ("long", "value")
    @test split_hint("long=<value>") == ("long", "value")
    @test split_hint("s") == ("s", nothing)
    @test split_hint("s=<value>") == ("s", "value")
    @test split_hint("s <value>") == ("s", "value")
    @test_throws Meta.ParseError split_hint("s=s")
    @test_throws Meta.ParseError split_hint("s=s=s")
end

@testset "split_option" begin
    @test split_option("--long") == ("long", nothing, nothing)
    @test split_option("--long <value>") == ("long", nothing, "value")
    @test split_option("--long=<value>") == ("long", nothing, "value")
    @test split_option("--short, -s") == ("short", "s", nothing)
    @test split_option("--short, -s <value>") == ("short", "s", "value")
    @test split_option("--short, -s=<value>") == ("short", "s", "value")
    @test_throws Meta.ParseError split_option("--short, -ss, -st")
    @test_throws Meta.ParseError split_option("--short, -ss")
    @test_throws Meta.ParseError split_option("--short, -ss=<value>")
    @test_throws Meta.ParseError split_option("--short, -t=<value>")
end

@testset "read_options" begin
    doc = Markdown.parse("""
    description of the command.

    # Options

    - `--short, -s`: short option using default hint.
    - `--short-space, -s <value>`: short option using given hint.
    - `--short-assign, -s=<value>`: short option using given hint.
    - `--long`: long option using default hint.
    - `--long-space <value>`: long option using given hint.
    - `--long-assign=<value>`: long option using given hint.
    - `--short_underscore, -s <value>`: short option with underscore.
    """)

    options = read_options(doc)
    @test options["short"] == JLMDOption(nothing, "short option using default hint.", true)
    @test options["short-space"] == JLMDOption("value", "short option using given hint.", true)
    @test options["short-assign"] == JLMDOption("value", "short option using given hint.", true)
    @test options["long"] == JLMDOption(nothing, "long option using default hint.", false)
    @test options["long-space"] == JLMDOption("value", "long option using given hint.", false)
    @test options["long-assign"] == JLMDOption("value", "long option using given hint.", false)
    @test options["short-underscore"] == JLMDOption("value", "short option with underscore.", true)
end

@testset "read_flags" begin
    doc = Markdown.parse("""
    description of the command.

    # Flags

    - `--short, -s`: short flag.
    - `--short-space, -s`: short flag with dash.
    - `--long`: long flag.
    - `--long-space`: long flag with dash.
    """)

    flags = read_flags(doc)
    @test flags["short"] == JLMDFlag("short flag.", true)
    @test flags["short-space"] == JLMDFlag("short flag with dash.", true)
    @test flags["long"] == JLMDFlag("long flag.", false)
    @test flags["long-space"] == JLMDFlag("long flag with dash.", false)
end

@testset "split_docstring" begin
    content = Markdown.parse("""
    description of the command.

    # Args

    - `arg1`: argument 1.
    - `arg2`: argument 2.
    - `arg3`: argument 3.
    - `arg4`: argument 4.

    # Options

    - `--short, -s`: short option using default hint.
    - `--short-space, -s <value>`: short option using given hint.
    - `--short-assign, -s=<value>`: short option using given hint.
    - `--long`: long option using default hint.
    - `--long-space <value>`: long option using given hint.
    - `--long-assign=<value>`: long option using given hint.
    - `--short_underscore, -s <value>`: short option with underscore.

    # Flags

    - `--short, -s`: short flag.
    - `--short-space, -s`: short flag with dash.
    - `--long`: long flag.
    - `--long-space`: long flag with dash.
    """)

    doc = split_docstring(content)
    @test doc.desc == "description of the command."
    @test doc.arguments["arg1"] == "argument 1."
    @test doc.flags["long"] == JLMDFlag("long flag.", false)
    @test doc.options["long"] == JLMDOption(nothing, "long option using default hint.", false)
end

@testset "reverse order" begin
    content = Markdown.parse("""
    description of the command.

    # Options

    - `-o, --option=<value>`: some random option.
    - `-o,--option_space=<value>`: some random option.
    """)

    doc = split_docstring(content)

    @test doc.options["option"] == JLMDOption("value", "some random option.", true)
    @test doc.options["option-space"] == JLMDOption("value", "some random option.", true)

    content = Markdown.parse("""
    description of the command.

    # Options

    - `-o, option=<value>`: some random option.
    """)

    @test_throws ErrorException split_docstring(content)
end

@testset "hint with space" begin
    doc = Markdown.parse("""
    release a package.

    # Arguments

    - `version_spec`: version number you want to release. Can be a specific version, "current"
        or either of
    - `path`: path to the project you want to release.

    # Options

    - `-r,--registry <registry name>`: registry you want to register the package.
        If the package has not been registered, ion will try to register
        the package in the General registry. Or the user needs to specify
        the registry to register using this option.
    - `-b, --branch <branch name>`: branch you want to register.
    - `--note <release note>`: optional, release note you would like to specify.
    """)

    options = read_options(doc)
    @test options["note"].hint == "release note"
    @test options["registry"].hint == "registry name"
end

end # TestMarkdown


text:
using Plots

cList = [complex(-0.4,-0.6), complex(0,-1),
            complex(-0.12,-0.75), complex(-0.6,0),
                complex(-0.8,0.16), complex(-0.4,0.6)]

function JuliaSet(c, opnum)
    abst = zeros((601,601))
    k = (1 + sqrt(1 + (4 * abs(c)))) / 2
    inum = 0
    for i in -1.5:0.005:1.5
        inum += 1
        jnum = 0
        for j in -1.5:0.005:1.5
            jnum += 1
            z = complex(i,j)
            for n in 1:opnum
                z = (z^2) + c
                if abs(z) > k
                    abst[inum,jnum] = abs(z)
                    break
                # else
                #     abst[inum,jnum] = k
                end
            end
        end
    end
    return abst
end

# theme(:dark)

for C in 1:length(cList)
    data = JuliaSet(cList[C], 100)
    heatmap(data, c = cgrad(:solar, rev = false), legend = false, border=:none, title = "C = $(cList[C])")
    savefig("C:\\Users\\Yaghoub\\Documents\\GitHub\\ComputationalPhysics-Fall2021\\ProblemSet2\\Figs\\Q1\\JuliaSet$C.png")
end


text:
using InteractiveUtils
using LLVM.Interop
using Test

@testset "atomics" begin

coreptr(r::Base.RefValue) = reinterpret(Ptr{eltype(r)}, pointer_from_objref(r))
llvmptr(r::Base.RefValue) = reinterpret(Core.LLVMPtr{eltype(r),0}, pointer_from_objref(r))

# Testing that ordering can be const-prop'ed
atomic_pointerref_monotonic(ptr) = Interop.atomic_pointerref(ptr, :monotonic)
atomic_pointerset_monotonic(ptr, x) = Interop.atomic_pointerset(ptr, x, :monotonic)
atomic_pointermodify_monotonic(ptr, op::OP, x) where {OP} =
    Interop.atomic_pointermodify(ptr, op, x, :monotonic)
atomic_pointerswap_monotonic(ptr, x) = Interop.atomic_pointerswap(ptr, x, :monotonic)
atomic_pointerreplace_monotonic(ptr, expected, desired) =
    Interop.atomic_pointerreplace(ptr, expected, desired, :monotonic, :monotonic)

muladd1(x, y) = muladd(x, y, one(x))

MONOTONIC = :monotonic

@testset for T in
             [Int8, Int16, Int32, Int64, UInt8, UInt16, UInt32, UInt64, Float32, Float64]
    r1 = Ref(one(T))
    r2 = Ref(one(T))
    GC.@preserve r1 r2 begin
        p1 = llvmptr(r1)
        sp1 = coreptr(r1)
        sp2 = coreptr(r2)

        @test atomic_pointerref_monotonic(p1) === r1[]
        @test Interop.atomic_pointerref(p1, MONOTONIC) === r1[]
        @test Interop.atomic_pointerref(p1, Val(:monotonic)) === r1[]
        if VERSION ≥ v"1.7"
            @test Interop.atomic_pointerref(p1, Val(:monotonic)) ===
                  Core.Intrinsics.atomic_pointerref(sp1, :monotonic)
        end

        types = typeof((p1, Val(:sequentially_consistent)))
        ir = sprint(io -> code_llvm(io, Interop.atomic_pointerref, types))
        @test occursin(r"load atomic .* seq_cst"m, ir)

        ir = sprint(io -> code_llvm(io, atomic_pointerref_monotonic, Tuple{typeof(p1)}))
        @test occursin(r"load atomic .* monotonic"m, ir)
        @test !occursin(r"load atomic .* seq_cst"m, ir)

        val = r1[] + one(r1[])
        @test begin
            Interop.atomic_pointerset(p1, val, MONOTONIC)
            r1[]
        end === val

        val = r1[] + one(r1[])
        @test begin
            Interop.atomic_pointerset(p1, val, Val(:monotonic))
            r1[]
        end === val

        if VERSION ≥ v"1.7"
            val = r1[] + one(r1[])
            r2[] = r1[]
            @test begin
                Interop.atomic_pointerset(p1, val, Val(:monotonic))
                r1[]
            end === begin
                Core.Intrinsics.atomic_pointerset(sp2, val, :monotonic)
                r2[]
            end
        end

        types = typeof((p1, val, Val(:sequentially_consistent)))
        ir = sprint(io -> code_llvm(io, Interop.atomic_pointerset, types))
        @test occursin(r"store atomic .* seq_cst"m, ir)

        ir = sprint(io -> code_llvm(io, atomic_pointerset_monotonic, typeof((p1, val))))
        @test occursin(r"store atomic .* monotonic"m, ir)
        @test !occursin(r"store atomic .* seq_cst"m, ir)

        ops = if T <: AbstractFloat
            Any[Interop.right, +, -]
        else
            Any[op for (_, op, _) in Interop.binoptable]
        end
        push!(ops, muladd1)

        old = r1[]
        @testset for op in ops
            r1[] = old
            val = one(old)
            @test begin
                Interop.atomic_pointermodify(p1, op, val, MONOTONIC)
                r1[]
            end === op(old, val)

            r1[] = old
            val = one(old)
            @test begin
                Interop.atomic_pointermodify(p1, op, val, Val(:monotonic))
                r1[]
            end === op(old, val)

            if VERSION ≥ v"1.7"
                r1[] = r2[] = old
                val = one(old)
                @test begin
                    Interop.atomic_pointermodify(p1, op, val, Val(:monotonic))
                    r1[]
                end === begin
                    Core.Intrinsics.atomic_pointermodify(sp2, op, val, :monotonic)
                    r2[]
                end
            end

            types = typeof((p1, op, val, Val(:sequentially_consistent)))
            ir = sprint(io -> code_llvm(io, Interop.atomic_pointermodify, types))
            if op === muladd1
                @test occursin(r"cmpxchg .* seq_cst"m, ir)
            else
                @test occursin(r"atomicrmw .* seq_cst"m, ir)
            end

            types = typeof((p1, op, val))
            ir = sprint(io -> code_llvm(io, atomic_pointermodify_monotonic, types))
            if op === muladd1
                @test occursin(r"cmpxchg .* monotonic"m, ir)
                @test !occursin(r"cmpxchg .* seq_cst"m, ir)
            else
                @test occursin(r"atomicrmw .* monotonic"m, ir)
                @test !occursin(r"atomicrmw .* seq_cst"m, ir)
            end
        end

        r1[] = old
        val = one(old)
        @test (atomic_pointerswap_monotonic(p1, val), r1[]) === (old, val)

        r1[] = old
        val = one(old)
        @test (Interop.atomic_pointerswap(p1, val, MONOTONIC), r1[]) === (old, val)

        r1[] = old
        val = one(old)
        @test (Interop.atomic_pointerswap(p1, val, Val(:monotonic)), r1[]) === (old, val)

        if VERSION ≥ v"1.7"
            r1[] = r2[] = old
            val = one(old)
            @test (Interop.atomic_pointerswap(p1, val, Val(:monotonic)), r1[]) ===
                  (Core.Intrinsics.atomic_pointerswap(sp2, val, :monotonic), r2[])
        end

        types = typeof((p1, val, Val(:sequentially_consistent)))
        ir = sprint(io -> code_llvm(io, Interop.atomic_pointerswap, types))
        @test occursin(r"atomicrmw xchg .* seq_cst"m, ir)

        ir = sprint(io -> code_llvm(io, atomic_pointerswap_monotonic, typeof((p1, val))))
        @test occursin(r"atomicrmw xchg .* monotonic"m, ir)
        @test !occursin(r"atomicrmw xchg .* seq_cst"m, ir)

        r1[] = old
        val = old + one(old)
        @test (atomic_pointerreplace_monotonic(p1, old, val), r1[]) ===
              ((; old, success = true), val)
        @test (atomic_pointerreplace_monotonic(p1, old, val), r1[]) ===
              ((; old = val, success = false), val)

        r1[] = old
        val = old + one(old)
        @test (Interop.atomic_pointerreplace(p1, old, val, MONOTONIC, MONOTONIC), r1[]) ===
              ((; old, success = true), val)
        @test (Interop.atomic_pointerreplace(p1, old, val, MONOTONIC, MONOTONIC), r1[]) ===
              ((; old = val, success = false), val)

        r1[] = old
        val = old + one(old)
        @test (
            Interop.atomic_pointerreplace(p1, old, val, Val(:monotonic), Val(:monotonic)),
            r1[],
        ) === ((; old, success = true), val)
        @test (
            Interop.atomic_pointerreplace(p1, old, val, Val(:monotonic), Val(:monotonic)),
            r1[],
        ) === ((; old = val, success = false), val)

        types = typeof((p1, old, val, Val(:acquire_release), Val(:acquire)))
        ir = sprint(io -> code_llvm(io, Interop.atomic_pointerreplace, types))
        @test occursin(r"cmpxchg .* acq_rel acquire"m, ir)

        types = typeof((p1, old, val))
        ir = sprint(io -> code_llvm(io, atomic_pointerreplace_monotonic, types))
        @test occursin(r"cmpxchg .* monotonic monotonic"m, ir)
        @test !occursin(r"cmpxchg .* seq_cst seq_cst"m, ir)
    end
end

primitive type B64 64 end
B64(x) = reinterpret(B64, convert(Int64, x))

@testset "primitive type" begin
    ref = Ref(B64(123))
    ord = Val(:monotonic)
    GC.@preserve ref begin
        ptr = llvmptr(ref)
        @test Interop.atomic_pointerref(ptr, ord) === ref[]
        @test begin
            Interop.atomic_pointerset(ptr, B64(456), ord)
            ref[]
        end === B64(456)
        @test Interop.atomic_pointerswap(ptr, B64(789), ord) === B64(456)
        @test ref[] === B64(789)
        @test Interop.atomic_pointerreplace(ptr, B64(789), B64(123), ord, ord) ===
              (old = B64(789), success = true)
        @test ref[] === B64(123)
        @test Interop.atomic_pointerreplace(ptr, B64(789), B64(123), ord, ord) ===
              (old = B64(123), success = false)
    end
end

end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/AsterReader.jl/blob/master/LICENSE

using Documenter, AsterReader

makedocs(modules=[AsterReader],
         format = :html,
         sitename = "AsterReader.jl",
         pages = [
                  "Introduction" => "index.md",
                  "API" => "api.md"
                 ])


text:
# exnode.jl - ExNode, the building block of ExGraph
#
# There are several categories of ExNodes:
#
#  * :call     - single function call, e.g. ExNode{:call}(z = x + y)
#  * :bcast    - broadcasting, e.g. ExNode{:bcast}(y = exp.(x))
#  * :(=)      - assignment, e.g. ExNode{:(=)}(y = x)
#  * :input    - input variable
#  * :constant - constant, e.g. ExNode{:constant}(x = 42)
#  * :opaque   - unparsed expression that can contain any subexpression;
#                a graph with opaque nodes isn't valid for most tasks,
#                but it may be converted to normalized graph using `reparse(g)`

# exnode

mutable struct ExNode{C}  # C - category of node, e.g. :call, :=, etc.
    var::Union{Symbol, Expr}       # variable name, possibly with indices
    ex::Any                        # primitive expression that produces the var
    guards::Vector{Expr}           # guards, turning ex to 0 when false
    val::Any                       # example value
    meta::Dict                     # node metadata & optional parameters
end

function ExNode{C}(var::Union{Symbol,Expr}, ex::Any;
                   guards=[], val=nothing, meta=Dict()) where C
    return ExNode{C}(var, ex, guards, val, meta)
end

function ExNode{C}(full_ex::Expr; val=nothing) where C
    @assert full_ex.head == :(=)
    var = full_ex.args[1]
    ex = without_guards(full_ex.args[2])
    guards = find_guards(full_ex.args[2])
    return ExNode{C}(var, ex; guards=guards, val=val)
end

## accessors

getcategory(nd::ExNode{C}) where {C} = C

getvar(nd::ExNode) = nd.var
setvar!(nd::ExNode, var::Union{Symbol,Expr}) = (nd.var = var)

varname(nd::ExNode)::Symbol = isa(nd.var, Symbol) ? nd.var : nd.var.args[1]
varidxs(nd::ExNode)::Vector = isa(nd.var, Symbol) ? [] : nd.var.args[2:end]
varidxs(nd::ExNode{:input})::Vector = IDX_NAMES[1:ndims(getvalue(nd))]

getexpr(nd::ExNode) = nd.ex
setexpr!(nd::ExNode, ex::Any) = (nd.ex = ex)

getguards(nd::ExNode) = nd.guards
setguards!(nd::ExNode, guards::Vector{Expr}) = (nd.guards = guards)

getvalue(nd::ExNode) = nd.val
setvalue!(nd::ExNode, val) = (nd.val = val)

Base.copy(nd::ExNode{C}; category=C, var=nd.var, ex=nd.ex,
          guards=nd.guards, val=nd.val, meta=nd.meta) where {C} =
    ExNode{category}(var, ex, guards, val, meta)

## pseudo-accessors

function getexpr_kw(nd::ExNode)
    if haskey(nd.meta, :kw) && !isempty(nd.meta[:kw])
        ex = copy(getexpr(nd))
        insert!(ex.args, 2, make_kw_params(nd.meta[:kw]))
        return ex
    else
        return getexpr(nd)
    end
end

function setexpr_kw!(nd::Union{ExNode{:call}, ExNode{:ctor}}, ex)
    f = ex.args[1]
    args, kw = parse_call_args(ex)
    nd.ex = :($f($(args...)))
    nd.meta[:kw] = kw
end

setexpr_kw!(nd::ExNode, ex) = setexpr!(nd, ex)

## to_expr and friends

"""
Convert ExNode to a full expression, e.g. for vectorized notation:

    z = x + y

or for indexed notation:

    z[i] = x[i] + y[i]
"""
function to_expr(nd::ExNode)
    var = getvar(nd)
    ex = with_guards(getexpr(nd), getguards(nd))
    return :($var = $ex)
end

"""
Same as to_expr(ExNode), but includes keyword arguments if any
"""
function to_expr_kw(nd::ExNode)
    var = getvar(nd)
    ex = with_guards(getexpr_kw(nd), getguards(nd))
    return :($var = $ex)
end

"""
Convert ExNode to a fortmat compatible with Einsum.jl
"""
function to_einsum_expr(nd::ExNode)
    depwarn_eingraph(:to_einsum_expr)
    ex = getexpr(nd)
    v = varname(nd)
    assign_ex = Expr(:(:=), getvar(nd), ex)
    macrocall_ex = Expr(:macrocall, Symbol("@einsum"), assign_ex)
    return Expr(:block, macrocall_ex)
end

function to_einsum_expr(nd::ExNode, var_size)
    depwarn_eingraph(:to_einsum_expr)
    ex = getexpr(nd)
    v = varname(nd)
    init_ex = :($v = zeros($(var_size)))  # TODO: handle data types other than Float64
    assign_ex = Expr(:(=), getvar(nd), ex)
    macrocall_ex = Expr(:macrocall, Symbol("@einsum"), assign_ex)
    return Expr(:block, init_ex, macrocall_ex)
end

## dependencies

"""Get names of dependenices of this node"""
dependencies(nd::ExNode{:input}) = Symbol[]
dependencies(nd::ExNode{:constant}) = get_var_names(getexpr(nd))
dependencies(nd::ExNode{:(=)}) = get_var_names(getexpr(nd))
dependencies(nd::ExNode{:call}) = get_var_names(getexpr(nd))
dependencies(nd::ExNode{:bcast}) = get_var_names(getexpr(nd))
dependencies(nd::ExNode{:tuple}) = [split_indexed(dep)[1] for dep in getexpr(nd).args]
dependencies(nd::ExNode{:ref}) = getexpr(nd).args
dependencies(nd::ExNode{:opaque}) = get_var_names(getexpr(nd); rec=true)
dependencies(nd::ExNode{:field}) = [getexpr(nd).args[1]]
dependencies(nd::ExNode{:ctor}) =
    [getexpr(nd).args[2], values(get(nd.meta, :kw, Dict()))...]

## node utils

# varsize(nd::ExNode{:tuple}) = map(size, getvalue(nd))
# varsize(nd::ExNode) = size(getvalue(nd))

# buffer_expr(nd::ExNode{:tuple})

function Base.show(io::IO, nd::ExNode{C}) where C
    val = getvalue(nd)
    if val == nothing
        val = "nothing"
    elseif isa(val, AbstractArray)
        val = "<$(typeof(val))>"
    elseif isa(val, Tuple)
        val = ([isa(v, AbstractArray) ?  "<$(typeof(v))>" : v for v in val]...,)
    elseif isstruct(val)
        val = "$(typeof(getvalue(nd)))"
    end
    ex_str = "ExNode{$C}($(to_expr_kw(nd)) | $val)"
    print(io, ex_str)
end

isindexed(nd::ExNode) = any(isref, get_vars(to_expr(nd)))

function rewrite(nd::ExNode{C}, pat, rpat; rw_opts...) where C
    full_ex = to_expr(nd)
    new_full_ex = rewrite(full_ex, pat, rpat; rw_opts...)
    @assert new_full_ex.head == :(=)
    var, ex = new_full_ex.args[1:2]
    return copy(nd; category=:opaque, var=var, ex=ex, val=nothing)
end


text:
function update_global_vars!(inputs::SimulationInputs, x::AbstractArray)
    index = get_global_vars(inputs)[:ω_sys_index]
    index == 0 && return
    #TO DO: Make it general for cases when ω is not a state (droop)!
    get_global_vars(inputs)[:ω_sys] = x[index]
    return
end

function system!(out::Vector{<:Real}, dx, x, inputs::SimulationInputs, t::Float64)
    sys = get_system(inputs)
    I_injections_r = get_aux_arrays(inputs)[1]
    I_injections_i = get_aux_arrays(inputs)[2]
    injection_ode = get_aux_arrays(inputs)[3]
    branches_ode = get_aux_arrays(inputs)[4]

    #Index Setup
    bus_size = get_bus_count(inputs)
    bus_vars_count = 2 * bus_size
    bus_range = 1:bus_vars_count
    injection_start = get_injection_pointer(inputs)
    injection_count = 1
    branches_start = get_branches_pointer(inputs)
    branches_count = 1
    update_global_vars!(inputs, x)

    #Network quantities
    V_r = @view x[1:bus_size]
    V_i = @view x[(bus_size + 1):bus_vars_count]
    Sbase = PSY.get_base_power(sys)
    fill!(I_injections_r, 0.0)
    fill!(I_injections_i, 0.0)

    for d in get_injectors_data(inputs)
        dynamic_device = PSY.get_dynamic_injector(d)
        bus_n = PSY.get_number(PSY.get_bus(d))
        bus_ix = get_lookup(inputs)[bus_n]
        n_states = PSY.get_n_states(dynamic_device)
        ix_range = range(injection_start, length = n_states)
        ode_range = range(injection_count, length = n_states)
        injection_count = injection_count + n_states
        injection_start = injection_start + n_states
        device!(
            x,
            injection_ode,
            view(V_r, bus_ix),
            view(V_i, bus_ix),
            view(I_injections_r, bus_ix),
            view(I_injections_i, bus_ix),
            ix_range,
            ode_range,
            dynamic_device,
            inputs,
        )
        out[ix_range] = injection_ode[ode_range] - dx[ix_range]
    end

    for d in PSY.get_components(PSY.ElectricLoad, sys)
        bus_n = PSY.get_number(PSY.get_bus(d))
        bus_ix = get_lookup(inputs)[bus_n]
        device!(
            view(V_r, bus_ix),
            view(V_i, bus_ix),
            view(I_injections_r, bus_ix),
            view(I_injections_i, bus_ix),
            d,
            inputs,
        )
    end

    for d in PSY.get_components(PSY.Source, sys)
        bus_n = PSY.get_number(PSY.get_bus(d))
        bus_ix = get_lookup(inputs)[bus_n]
        device!(
            view(V_r, bus_ix),
            view(V_i, bus_ix),
            view(I_injections_r, bus_ix),
            view(I_injections_i, bus_ix),
            d,
            inputs,
        )
    end

    if get_dyn_lines(inputs)
        dyn_branches = PSY.get_components(PSY.DynamicBranch, sys)
        for br in dyn_branches
            arc = PSY.get_arc(br)
            n_states = PSY.get_n_states(br)
            from_bus_number = PSY.get_number(arc.from)
            to_bus_number = PSY.get_number(arc.to)
            bus_ix_from = get_lookup(inputs)[from_bus_number]
            bus_ix_to = get_lookup(inputs)[to_bus_number]
            ix_range = range(branches_start, length = n_states)
            ode_range = range(branches_count, length = n_states)
            branches_count = branches_count + n_states
            branch!(
                x,
                dx,
                branches_ode,
                #Get Voltage data
                view(V_r, bus_ix_from),
                view(V_i, bus_ix_from),
                view(V_r, bus_ix_to),
                view(V_i, bus_ix_to),
                #Get Current data
                view(I_injections_r, bus_ix_from),
                view(I_injections_i, bus_ix_from),
                view(I_injections_r, bus_ix_to),
                view(I_injections_i, bus_ix_to),
                ix_range,
                ode_range,
                br,
                inputs,
            )
            out[ix_range] = branches_ode[ode_range] - dx[ix_range]
        end
    end

    kirchoff_laws!(inputs, V_r, V_i, I_injections_r, I_injections_i, dx)
    out[bus_range] = get_aux_arrays(inputs)[6]
end


text:
# ------------------------------------------------------------------------------

PairIntTuples = Tuple{Tuple{Vararg{Int}},Tuple{Vararg{Int}}}

function det_minor_expansion_inner(
        m::MatElem{<: T},
        discarded::PairIntTuples,
        cache::Dict{PairIntTuples,T}
    ) where T <: RingElem
    n = size(m, 1);
    if length(discarded[1]) == n
        return 1;
    end
    if discarded in keys(cache)
        return cache[discarded];
    end
    result = 0;
    row = minimum(setdiff(Set(1:n), Set(discarded[1])));
    dis_rows = Tuple(sort([[i for i in discarded[1]]; row]));
    sign = 1;
    for col = 1:n
        if !(col in discarded[2])
            dis_cols = Tuple(sort([[i for i in discarded[2]]; col]));
            result += sign * m[row, col] * det_minor_expansion_inner(m, (dis_rows, dis_cols), cache);
            sign = -1 * sign;
        end
    end
    if length(discarded[1]) > 1
        cache[discarded] = result;
    end
    if length(discarded[1]) < 3
        @debug "Discarded: $discarded; $(Dates.now())"
        flush(stdout)
    end
    return result;
end

function det_minor_expansion(m::MatElem{T}) where T <: RingElem
    cache = Dict{PairIntTuples,T}();
    return det_minor_expansion_inner(m, (Tuple{}(), Tuple{}()), cache);
end

# ------------------------------------------------------------------------------

"""
    Bezout_matrix(f, g, var_elim)

Compute the Bezout matrix of two polynomials f, g with respect to var_elim

Inputs:
- `f` - first polynomial
- `g` - second polynomial
- `var_elim` - variable, of which f and g are considered as polynomials

Output:
- `M::MatrixElem` - The Bezout matrix
"""
function Bezout_matrix(f::P, g::P, var_elim::P) where P <: MPolyElem
    parent_ring = parent(f)
    deg_f = Nemo.degree(f, var_elim)
    deg_g = Nemo.degree(g, var_elim)
    n = max(deg_f, deg_g)
    coeffs_f = [coeff(f, [var_elim], [i]) for i in 0:n]
    coeffs_g = [coeff(g, [var_elim], [i]) for i in 0:n]
    GL = AbstractAlgebra.MatrixSpace(parent_ring, n, n)
    M = zero(GL)
    for i in 1:n
        for j in 1:n
            M[i, j] = sum(coeffs_f[j + k + 1] * coeffs_g[i - k] - coeffs_g[j + k + 1] * coeffs_f[i - k] for k in 0:min(i - 1, n - j))
        end
    end
    return M
end

# ------------------------------------------------------------------------------

"""
    Sylvester_matrix(f, g, var_elim)

Compute the Sylvester matrix of two polynomials f, g with respect to var_elim
Inputs:
- `f` - first polynomial
- `g` - second polynomial
- `var_elim` - variable, of which f and g are considered as polynomials
Output:
- `M::MatrixElem` - The Sylvester matrix
"""
function Sylvester_matrix(f::P, g::P, var_elim::P) where P <: MPolyElem
    parent_ring = parent(f)
    deg_f = Nemo.degree(f, var_elim)
    deg_g = Nemo.degree(g, var_elim)
    n = deg_f + deg_g
    GL = AbstractAlgebra.MatrixSpace(parent_ring, n, n)
    M = zero(GL)
    for i in 1:deg_f
        for j in 0:deg_g
            M[i, j + i] = coeff(g, [var_elim], [j])
        end
    end
    for i in 1:deg_g
        for j in 0:deg_f
            M[i + deg_f, j + i] = coeff(f, [var_elim], [j])
        end
    end

    return M
end

# ------------------------------------------------------------------------------

"""
    simplify_matrix(M)

Eliminate GCD of entries of every row and column

Input:
- `M::MatrixElem` - matrix to be simplified

Output:
- `M::MatrixElem` - Simplified matrix
- `extra_factors::Vector{AbstractAlgebra.MPolyElem}` - array of GCDs eliminated from `M`.
"""
function simplify_matrix(M::MatElem{P}) where P <: MPolyElem
    """
    An auxiliary function taking a list of coordinates of cells
    and dividing them by their gcd.
    Returns the gcd
    """
    function _simplify_range(coords::Array{Tuple{Int,Int},1})
        gcd_temp = M[coords[1]...]
        for c in coords[2:end]
            gcd_temp = gcd(gcd_temp, M[c...])
        end
        if gcd_temp != 1
            for c in coords
                M[c...] = divexact(M[c...], gcd_temp)
            end
        end
        return gcd_temp
    end

    extra_factors = Array{P,1}()
    rows_cols = Array{Array{Tuple{Int,Int},1},1}()
    # adding all rows
    for i in 1:nrows(M)
        push!(rows_cols, [(i, j) for j in 1:ncols(M)])
    end
    # adding all columns
    for j in 1:ncols(M)
        push!(rows_cols, [(i, j) for i in 1:nrows(M)])
    end
    # eliminating factors
    for r in rows_cols
        gcd_temp = _simplify_range(r)
        if gcd_temp != 1
            push!(extra_factors, gcd_temp)
        end
    end

    return M, extra_factors
end

# ------------------------------------------------------------------------------

# Definition: for an irreducible polynomial P in n variables, we will call
# an iterator Julia object *generic point generator* if it generates
# an infinite random sequence of points in C^n such that
#    1. P vanishes at each of the points
#    2. for every other irreducible polynomial, the probability of
#    vanishing at each of these points is zero
#
# Each returned point is dictionary from the variables to the values

# ------------------------------------------------------------------------------

abstract type PointGenerator{P} end

mutable struct ODEPointGenerator{P} <: PointGenerator{P}
    ode::ODE{P}
    big_ring::MPolyRing
    precision::Int
    cached_points::Array{Dict{P,<: FieldElem},1}
    number_type::Type

    function ODEPointGenerator{P}(ode::ODE{P}, big_ring::MPolyRing) where P <: MPolyElem
        prec = length(ode.x_vars) + 1
        number_type = typeof(one(base_ring(big_ring)))
        return new(ode, big_ring, prec, Array{Dict{P,number_type}}[], number_type)
    end
end

# ------------------------------------------------------------------------------

function Base.iterate(gpg::ODEPointGenerator{P}, i::Int=1) where P <: MPolyElem{<: FieldElem}
    if i > length(gpg.cached_points)
        @debug "Generating new point on the variety"
        sample_max = i * 50
        result = undef
        while true
            @debug "Preparing initial condition"
            flush(stdout)
            base_field = base_ring(gpg.big_ring)
            param_values = Dict{P,Int}(p => rand(1:sample_max) for p in gpg.ode.parameters)
            initial_conditions = Dict{P,Int}(x => rand(1:sample_max) for x in gpg.ode.x_vars)
            input_values = Dict{P,Array{Int,1}}(u => [rand(1:sample_max) for _ in 1:gpg.precision] for u in gpg.ode.u_vars)
            @debug "Computing a power series solution"
            flush(stdout)
            ps_solution = undef
            try
                ps_solution = power_series_solution(gpg.ode, param_values, initial_conditions, input_values, gpg.precision)
            catch e
                @debug "$e"
                flush(stdout)
                continue
            end

            @debug "Constructing the point"
            flush(stdout)
            result = Dict{P,gpg.number_type}(switch_ring(p, gpg.big_ring) => base_field(c) for (p, c) in param_values)
            for u in gpg.ode.u_vars
                for i in 0:(gpg.precision - 1)
                    result[str_to_var(var_to_str(u) * "_$i", gpg.big_ring)] = coeff(ps_solution[u], i) * factorial(i)
                end
            end
            for y in gpg.ode.y_vars
                for j in 0:(gpg.precision - 1)
                    result[str_to_var(var_to_str(y) * "_$j", gpg.big_ring)] = coeff(ps_solution[y], j) * factorial(j)
                end
            end
            for x in gpg.ode.x_vars
                result[switch_ring(x, gpg.big_ring)] = coeff(ps_solution[x], 0)
                result[str_to_var(var_to_str(x) * "_dot", gpg.big_ring)] = coeff(ps_solution[x], 1)
            end
            break
        end
        push!(gpg.cached_points, result)
    end
    return (gpg.cached_points[i], i + 1)
end

# ------------------------------------------------------------------------------

"""
    choose(polys, generic_point_generator)

Input:
- `polys` - an array of distinct irreducible polynomials in the same ring
- `generic_point_generator` - a generic point generator as described above for one of polys

Output:
- the polynomial that vanishes at the `generic_point_generator`
"""
function choose(polys::Array{P,1}, generic_point_generator) where P <: MPolyElem{<: FieldElem}
    vars = gens(parent(polys[1]))
    for p in generic_point_generator
        if length(polys) <= 1
            break
        end
        # get accounts for the fact that the big ring may contain some auxiliary variables, e.g. rand_proj_var
        point = [get(p, v, zero(base_ring(parent(polys[1])))) for v in vars]
        polys = filter(e -> (evaluate(e, point) == 0), polys)
        flush(stdout)
    end
    return polys[1]
end

# ------------------------------------------------------------------------------

"""
    eliminate_var(f, g, var_elim, generic_point_generator)

Eliminate variable from a pair of polynomials

Input:
- `f` and `g` - polynomials
- `var_elim` - variable to be eliminated
- `generic_point_generator` - a generic point generator object for the factor
      of the resultant of `f` and `g` of interest

Output:
- `polynomial` - the desired factor of the resultant of `f` and `g`
"""
function eliminate_var(f::P, g::P, var_elim::P, generic_point_generator) where P <: MPolyElem{<: FieldElem}
    # Linear comb
    while f != 0 && g != 0
        if Nemo.degree(f, var_elim) > Nemo.degree(g, var_elim)
            f, g = g, f
        end
        lf = coeff(f, [var_elim], [Nemo.degree(f, var_elim)])
        lg = coeff(g, [var_elim], [Nemo.degree(g, var_elim)])
        (flag, q) = divides(lg, lf)
        if flag
            @debug "\t Decreasing degree with linear combination $(Dates.now())"
            flush(stdout)
            g = g - q * f * var_elim^(Nemo.degree(g, var_elim) - Nemo.degree(f, var_elim))
        elseif (Nemo.degree(g, var_elim) == Nemo.degree(f, var_elim))
            (flag, q) = divides(lf, lg)
            if flag
                @debug "\t Decreasing degree with linear combination $(Dates.now())"
                flush(stdout)
                f = f - q * g
            else
                break
            end
        else
            break
        end
    end

    if f == 0 || g == 0
        return f + g
    end

    # Case (f(v^d), g(v^d)):
    list_deg = []
    for d in 0:Nemo.degree(f, var_elim)
        if coeff(f, [var_elim], [d]) != 0
            push!(list_deg, d)
        end
    end
    for d in 0:Nemo.degree(g, var_elim)
        if coeff(g, [var_elim], [d]) != 0
            push!(list_deg, d)
        end
    end
    gcd_deg = list_deg[1]
    for ele in list_deg[2:end]
        gcd_deg = gcd(gcd_deg, ele)
    end
    if gcd_deg > 1
        f = sum(coeff(f, [var_elim], [gcd_deg * i]) * (var_elim^i) for i in 0:(Nemo.degree(f, var_elim) ÷ gcd_deg))
        g = sum(coeff(g, [var_elim], [gcd_deg * i]) * (var_elim^i) for i in 0:(Nemo.degree(g, var_elim) ÷ gcd_deg))
    end

    resultant = undef
    matrix_factors = []
    if Nemo.degree(f, var_elim) == 0
        resultant = f
    else
        if Nemo.degree(f, var_elim) > 1
            @debug "Calculating Bezout Matrix $(Dates.now())"
            flush(stdout)
            M = Bezout_matrix(f, g, var_elim)
        else
            @debug "Calculating Sylvester matrix"
            flush(stdout)
            M = Sylvester_matrix(f, g, var_elim)
        end
        @debug "Simplifying the matrix $(Dates.now())"
        flush(stdout)
        M_simp, matrix_factors = simplify_matrix(M)
        @debug "Removed factors $(map(length, matrix_factors))"
        M_size = zero(Nemo.MatrixSpace(Nemo.ZZ, ncols(M_simp), ncols(M_simp)))
        for i in 1:ncols(M_simp)
            for j in 1:ncols(M_simp)
                M_size[i,j] = length(M_simp[i,j])
            end
        end
        @debug "\t Matrix size: \n $M_size"
        @debug "\t Computing determinant $(Dates.now())"
        flush(stdout)
        resultant = det_minor_expansion(M_simp)
    end
    # Step 4: Eliminate extra factors
    factors = fast_factor(resultant)
    for mfac in matrix_factors
        for fac in fast_factor(mfac)
            if !(fac in factors)
                push!(factors, fac)
            end
        end
    end
    res = choose(factors, generic_point_generator)
    for f in factors
        if f != res
            @debug "\t \t Size of extra factor: $(length(f)); $(Dates.now())"
            @debug "\t \t It is $f"
            flush(stdout)
        end
    end
    return res
end


text:
const SECRET_TOKEN = "c7c1784996e9592516a880b9cee2b635aff92d1dcd694c5862b2dc4fb3217fd4"

# Follow the instructions at https://github.com/randyzwitch/Twitter.jl to setup
# access for the Twitter library
const TWITTER_TOKENS = ["... your consumer key ...",
                        "... your consumer token ...",
                        "... your oauth token ...",
                        "... your oauth secret..."]

# In order to use the Aylien Sentiment Analysis API you need to setup a free account with them: http://docs.aylien.com/docs/api-key
const AYLIEN_TOKENS = Dict{Symbol,String}(:app_id => "... your app id ...", :api_key => "... your api key ...")

max_stars_repo_name: essenciary/genie-demo-tweet-stats

max_issues_repo_name: essenciary/genie-demo-tweet-stats

max_forks_repo_name: essenciary/genie-demo-tweet-stats


text:
using StaticArrays, ForwardDiff, SciMLBase
import SciMLBase: init, step!, isinplace, reinit!, u_modified!
import Base: show

export reinit!

#####################################################################################
#                            DiscreteIntegrator                                    #
#####################################################################################
mutable struct MinimalDiscreteIntegrator{IIP, S, D, F, P}
    f::F      # integrator f
    u::S      # integrator state
    t::Int    # integrator "time" (counter)
    dummy::S  # dummy, used only in the IIP version
    p::P      # parameter container, I don't know why
    t0::Int   # initial time (only for reinit!)
end
const MDI = MinimalDiscreteIntegrator
isinplace(::MDI{IIP}) where {IIP} = IIP
stateeltype(::MDI{IIP, S}) where {IIP, S} = eltype(S)
stateeltype(::MDI{IIP, S}) where {IIP, S<:Vector{<:AbstractArray{T}}} where {T} = T

function reinit!(integ::MDI, u = integ.u, Q0 = nothing; t0 = integ.t0)
    integ.u = u
    integ.dummy = copy(u)
    integ.t = t0
    if Q0 !== nothing
        set_deviations!(integ, Q0)
    end
    return
end

function (integ::MinimalDiscreteIntegrator)(t::Real)
    if t == integ.t
        return integ.u
    else
        error("Can't extrapolate discrete systems!")
    end
end

# Get state for parallel:
get_state(a::MDI{IIP, S}) where {IIP, S<:Vector{<:AbstractVector}} = a.u[1]
get_state(a::MDI{IIP, S}, k) where {IIP, S<:Vector{<:AbstractVector}} = a.u[k]
function set_state!(
    integ::MDI{Alg, S}, u::AbstractVector, k::Int = 1
    ) where {Alg, S<:Vector{<:AbstractVector}}
    integ.u[k] = u
end

# for autodiffed in-place version (where state is Matrix)
get_state(integ::MDI{Alg, S}) where {Alg, S<:AbstractMatrix} = integ.u[:, 1]
set_state!(integ::MDI{Alg, S}, u) where {Alg, S<:AbstractMatrix} = (integ.u[:, 1] .= u)
get_deviations(integ::MDI{Alg, S}) where {Alg, S<:AbstractMatrix} =
    @view integ.u[:, 2:end]
set_deviations!(integ::MDI{Alg, S}, Q) where {Alg, S<:AbstractMatrix} =
    (integ.u[:, 2:end] = Q)

#####################################################################################
#                                   Stepping                                        #
#####################################################################################
# IIP version
function step!(integ::MDI{true})
    # vector swap
    integ.dummy, integ.u = integ.u, integ.dummy
    integ.f(integ.u, integ.dummy, integ.p, integ.t)
    integ.t += 1
    return
end
function step!(integ::MDI{true}, N)
    for i in 1:N
        integ.dummy, integ.u = integ.u, integ.dummy
        integ.f(integ.u, integ.dummy, integ.p, integ.t)
        integ.t += 1
    end
    return
end

# OOP version
step!(integ::MDI{false}) =
(integ.u = integ.f(integ.u, integ.p, integ.t); integ.t +=1; nothing)
function step!(integ::MDI{false}, N)
    for i in 1:N
        integ.u = integ.f(integ.u, integ.p, integ.t)
        integ.t += 1
    end
    return
end

u_modified!(integ::MDI, ::Bool) = nothing

step!(integ::MDI, N, stop_at_tdt) = step!(integ, N)

#####################################################################################
#                           TangentDiscreteIntegrator                               #
#####################################################################################
# For discrete systems a special, super-performant version of tangent integrator
# can exist. We make this here.
mutable struct TangentDiscreteIntegrator{IIP, S, D, F, P, JAC, JM, WM}
    f::F            # integrator f
    u::S            # integrator state
    t::Int          # integrator "time" (counter)
    dummy::S        # dummy, used only in the IIP version
    p::P            # parameter container, I don't know why
    t0::Int         # initial time (only for reinit!)
    jacobian::JAC   # jacobian function
    J::JM           # jacobian matrix
    W::WM           # tangent vectors (in form of matrix)
    dummyW::WM      # dummy, only used in IIP version
end

const TDI = TangentDiscreteIntegrator
stateeltype(::TDI{IIP, S}) where {IIP, S} = eltype(S)

u_modified!(t::TDI, a) = nothing

# set_state is same as in standard, see dynamicalsystem.jl

get_deviations(t::TDI) = t.W
set_deviations!(t::TDI, Q) = (t.W = Q) # works for both Matrix/SMatrix Q
function set_deviations!(t::TDI, Q::LinearAlgebra.AbstractQ)
    # This method exists because there is a more efficient way to copy columns of Q
    # without transforming `Q` into a matrix
    # integ.u[:, 2:end] .= Matrix(Q) <- this is bad, it allocates!
    copyto!(t.W, I)
    lmul!(Q, t.W)
end

function reinit!(integ::TDI, u = integ.u, Q0 = nothing; t0 = integ.t0)
    set_state!(integ, u)
    if Q0 !== nothing
        set_deviations!(integ, Q0)
    end
    integ.t = t0
    return
end

#####################################################################################
#                                 Tangent Stepping                                  #
#####################################################################################
function step!(integ::TDI{true})
    integ.dummy, integ.u = integ.u, integ.dummy
    integ.dummyW, integ.W = integ.W, integ.dummyW

    integ.f(integ.u, integ.dummy, integ.p, integ.t)
    integ.jacobian(integ.J, integ.u, integ.p, integ.t)

    mul!(integ.W, integ.J, integ.dummyW)
    integ.t += 1
    return
end
function step!(integ::TDI{true}, N::Real)
    for i in 1:N
        integ.dummy, integ.u = integ.u, integ.dummy
        integ.dummyW, integ.W = integ.W, integ.dummyW

        integ.f(integ.u, integ.dummy, integ.p, integ.t)
        integ.jacobian(integ.J, integ.u, integ.p, integ.t)

        mul!(integ.W, integ.J, integ.dummyW)
        integ.t += 1
    end
    return
end

function step!(integ::TDI{false})
    integ.u = integ.f(integ.u, integ.p, integ.t)
    J = integ.jacobian(integ.u, integ.p, integ.t)
    integ.W = J*integ.W
    integ.t += 1
    return
end
function step!(integ::TDI{false}, N::Real)
    for i in 1:N
        integ.u = integ.f(integ.u, integ.p, integ.t)
        J = integ.jacobian(integ.u, integ.p, integ.t)
        integ.W = J*integ.W
        integ.t += 1
    end
    return
end

step!(integ::TDI, N, stop_at_tdt) = step!(integ, N)

#####################################################################################
#                            Creating Integrators                                   #
#####################################################################################
function integrator(ds::DDS{IIP, S, D, F, P}, u0 = ds.u0; kwargs...) where {IIP, S, D, F, P}
    return MinimalDiscreteIntegrator{IIP, S, D, F, P}(
    ds.f, S(u0), ds.t0, S(copy(u0)), ds.p, ds.t0)
end

function tangent_integrator(ds::DDS, k::Int = dimension(ds); kwargs...)
    return tangent_integrator(
    ds, orthonormal(dimension(ds), k); kwargs...)
end

function tangent_integrator(ds::DDS{IIP, S, D, F, P, JAC, JM}, Q0::AbstractMatrix;
    u0 = ds.u0, kwargs...) where {IIP, S, D, F, P, JAC, JM}

    Q = safe_matrix_type(Val{IIP}(), Q0)
    s = safe_state_type(Val{IIP}(), u0)
    size(Q)[2] > dimension(ds) && throw(ArgumentError(
    "It is not possible to evolve more tangent vectors than the system's dimension!"
    ))

    WM = typeof(Q)

    return TDI{IIP, S, D, F, P, JAC, JM, WM}(ds.f, s, ds.t0, deepcopy(s),
    ds.p, ds.t0, ds.jacobian::JAC, deepcopy(ds.J), Q, deepcopy(Q))
end

# Auto-diffed in-place version
function tangent_integrator(ds::DDS{true, S, D, F, P, JAC, JM, true},
    Q0::AbstractMatrix;
    u0 = ds.u0, kwargs...) where {S, D, F, P, JAC, JM}

    t0 = ds.t0
    R = D + length(Q0)
    k = size(Q0)[2]
    Q = safe_matrix_type(Val{true}(), Q0)
    u = safe_state_type(Val{true}(), u0)
    size(Q)[2] > dimension(ds) && throw(ArgumentError(
    "It is not possible to evolve more tangent vectors than the system's dimension!"
    ))

    tangentf = create_tangent_iad(ds.f, ds.J, u, ds.p, t0, Val{k}())

    TF = typeof(tangentf)

    s = hcat(u, Q)
    SS = typeof(s)

    return MDI{true, SS, R, TF, P}(tangentf, s, t0, deepcopy(s), ds.p, t0)
end

function parallel_integrator(ds::DDS{IIP, S, D}, states; kwargs...) where {IIP, S, D}
    peom, st = create_parallel(ds, states)
    F = typeof(peom); X = typeof(st); P = typeof(ds.p)
    return MDI{true, X, D, F, P}(peom, st, ds.t0, deepcopy(st), ds.p, ds.t0)
end

#####################################################################################
#                                 Trajectory                                        #
#####################################################################################
function trajectory(
        ds::DDS{IIP, S, D}, t, u = ds.u0;
        Δt = 1, Ttr = 0, save_idxs = nothing, kwargs...
    ) where {IIP, S, D}
    a = svector_access(save_idxs)
    trajectory(ds, t, u, Δt, Ttr, a)
end

function trajectory(ds::DDS{IIP, S, D}, t, u, Δt, Ttr, a)  where {IIP, S, D}
    Δt = round(Int, Δt)
    T = eltype(S)
    integ = integrator(ds, u)
    ti = ds.t0
    tvec = (ti+Ttr):Δt:(t+ti+Ttr)
    L = length(tvec)
    T = eltype(get_state(ds))
    X = isnothing(a) ? D : length(a)
    data = Vector{SVector{X, T}}(undef, L)
    Ttr != 0 && step!(integ, Ttr)
    data[1] = obtain_access(integ.u, a)
    for i in 2:L
        step!(integ, Δt)
        data[i] = SVector{X, T}(obtain_access(integ.u, a))
    end
    return Dataset(data)
end

# One-dimensional version → Vector
function trajectory(
        ds::DDS{false, S, 1}, t, u = ds.u0;
    	Δt = 1, Ttr = 0, kwargs...
    ) where {S}
    Δt = round(Int, Δt)
    ti = ds.t0
    tvec = (ti+Ttr):Δt:(t+ti+Ttr)
    L = length(tvec)
    integ = integrator(ds, u)
    data = Vector{S}(undef, L)
    Ttr != 0 && step!(integ, Ttr)
    data[1] = integ.u
    for i in 2:L
        step!(integ, Δt)
        data[i] = integ.u
    end
    return data
end

#####################################################################################
#                                Pretty-Printing                                    #
#####################################################################################
Base.summary(ds::MDI) =
"Discrete integrator (e.o.m.: $(eomstring(ds.f)))"
Base.summary(ds::MDI{true, S}) where {S<:Vector{<:AbstractArray}} =
"Discrete parallel integrator with $(length(ds.u)) states "*
"(e.o.m.: $(eomstring(ds.f)))"

function Base.show(io::IO, ds::MDI)
    ps = 3
    text = summary(ds)
    println(io, text)
    println(io, rpad(" t: ", ps)*"$(ds.t)")
    prefix = rpad(" u: ", ps)
    print(io, prefix); printlimited(io, get_state(ds)', Δx = length(prefix)); print(io, "\n")
end
function Base.show(io::IO, ds::MDI{true, S}) where {S<:Vector{<:AbstractArray}}
    ps = 3
    text = summary(ds)
    print(io, text*"\n",
    rpad(" t: ", ps)*"$(ds.t)\n",
    rpad(" states: ", ps)*"\n"
    )
    s = sprint(io -> show(IOContext(io, :limit=>true),
    MIME"text/plain"(), ds.u))
    s = join(split(s, '\n')[2:end], '\n')
    print(io, s)
end

Base.summary(ds::TDI) =
"Discrete tangent-space integrator "*
"(e.o.m.: $(eomstring(ds.f)), jacobian: $(eomstring(ds.jacobian)))"

function Base.show(io::IO, ds::TDI)
    ps = 3
    text = summary(ds)
    println(io, text)
    println(io, rpad(" t: ", ps)*"$(ds.t)")
    prefix = rpad(" u: ", ps)
    print(io, prefix); printlimited(io, get_state(ds)', Δx = length(prefix)); print(io, "\n")
    println(io, " deviation vectors:")
    s = sprint(io -> show(IOContext(io, :limit=>true),
    MIME"text/plain"(), get_deviations(ds)))
    s = join(split(s, '\n')[2:end], '\n')
    print(io, s)
end


text:
ChainRulesCore.@scalar_rule(xlogx(x::Real), (1 + log(x),))
ChainRulesCore.@scalar_rule(xlogy(x::Real, y::Real), (log(y), x / y,))
ChainRulesCore.@scalar_rule(xlog1py(x::Real, y::Real), (log1p(y), x / (1 + y),))

ChainRulesCore.@scalar_rule(logistic(x::Real), (Ω * (1 - Ω),))
ChainRulesCore.@scalar_rule(logit(x::Real), (inv(x * (1 - x)),))
ChainRulesCore.@scalar_rule(logcosh(x::Real), tanh(x))
ChainRulesCore.@scalar_rule(log1psq(x::Real), (2 * x / (1 + x^2),))
ChainRulesCore.@scalar_rule(log1pexp(x::Real), (logistic(x),))
ChainRulesCore.@scalar_rule(log1mexp(x::Real), (-exp(x - Ω),))
ChainRulesCore.@scalar_rule(log2mexp(x::Real), (-exp(x - Ω),))
ChainRulesCore.@scalar_rule(logexpm1(x::Real), (exp(x - Ω),))

ChainRulesCore.@scalar_rule(logaddexp(x::Real, y::Real), (exp(x - Ω), exp(y - Ω)))
ChainRulesCore.@scalar_rule(
    logsubexp(x::Real, y::Real),
    (x > y ? exp(x - Ω) : -exp(x - Ω), x > y ? -exp(y - Ω) : exp(y - Ω)),
)

function ChainRulesCore.frule((_, Δx), ::typeof(logsumexp), x::AbstractArray{<:Real}; dims=:)
    Ω = logsumexp(x; dims=dims)
    ΔΩ = sum(exp.(x .- Ω) .* Δx; dims=dims)
    return Ω, ΔΩ
end
function ChainRulesCore.rrule(::typeof(logsumexp), x::AbstractArray{<:Real}; dims=:)
    Ω = logsumexp(x; dims=dims)
    project_x = ChainRulesCore.ProjectTo(x)
    function logsumexp_pullback(Ω̄)
        x̄ = ChainRulesCore.InplaceableThunk(
            Δ -> Δ .+= Ω̄ .* exp.(x .- Ω),
            ChainRulesCore.@thunk(project_x(Ω̄ .* exp.(x .- Ω))),
        )
        return ChainRulesCore.NoTangent(), x̄
    end
    return Ω, logsumexp_pullback
end

# no rules for mutating functions currently:
# https://juliadiff.org/ChainRulesCore.jl/stable/writing_good_rules.html#Which-functions-need-rules?
function ChainRulesCore.frule((_, Δx), ::typeof(softmax), x::AbstractArray{<:Real}; dims=:)
    Ω = softmax(x; dims=dims)
    ΔΩ = if dims === (:)
        Ω .* (Δx .- LinearAlgebra.dot(Ω, Δx))
    else
        ΩΔx = Ω .* Δx
        ΩΔx .- Ω .* sum(ΩΔx; dims=dims)
    end
    return Ω, ΔΩ
end
function ChainRulesCore.rrule(::typeof(softmax), x::AbstractArray{<:Real}; dims=:)
    Ω = softmax(x; dims=dims)
    Ωcopy = copy(Ω)
    project_x = ChainRulesCore.ProjectTo(x)
    function softmax_pullback(Ω̄)
        x̄ = if dims === (:)
            ChainRulesCore.InplaceableThunk(
                Δ -> Δ .+= Ωcopy .* (Ω̄ .- LinearAlgebra.dot(Ωcopy, Ω̄)),
                ChainRulesCore.@thunk(project_x(Ωcopy .* (Ω̄ .- LinearAlgebra.dot(Ωcopy, Ω̄)))),
            )
        else
            ΩcopyΩ̄  = Ωcopy .* Ω̄
            ChainRulesCore.InplaceableThunk(
                Δ -> Δ .+= ΩcopyΩ̄  .- Ωcopy .* sum(ΩcopyΩ̄; dims=dims),
                ChainRulesCore.@thunk(project_x(ΩcopyΩ̄  .- Ωcopy .* sum(ΩcopyΩ̄; dims=dims))),
            )
        end
        return ChainRulesCore.NoTangent(), x̄
    end
    return Ω, softmax_pullback
end