text:
# This file reimplements gmsh/tutorial/t2.geo in Julia. Comments focus on the new
# API functions used, compared to the ones introduced in t1.jl.

import gmsh

# nice shortcuts
model = gmsh.model
factory = model.geo

# If ARGS is passed, Gmsh will parse the commandline in the same way as the
# standalone Gmsh app.
gmsh.initialize(ARGS)

gmsh.option.setNumber("General.Terminal", 1)

model.add("t2")

# Copied from t1.jl...
lc = 1e-2
factory.addPoint(0, 0, 0, lc, 1)
factory.addPoint(.1, 0,  0, lc, 2)
factory.addPoint(.1, .3, 0, lc, 3)
factory.addPoint(0, .3, 0, lc, 4)
factory.addLine(1, 2, 1)
factory.addLine(3, 2, 2)
factory.addLine(3, 4, 3)
factory.addLine(4, 1, 4)
factory.addCurveLoop([4, 1, -2, 3], 1)
factory.addPlaneSurface([1], 1)
model.addPhysicalGroup(0, [1, 2], 1)
model.addPhysicalGroup(1, [1, 2], 2)
model.addPhysicalGroup(2, [1], 6)
model.setPhysicalName(2, 6, "My surface")
# ...end of copy

factory.addPoint(0, .4, 0, lc, 5)
factory.addLine(4, 5, 5)

# Geometrical transformations take a vector of pairs of integers as first
# argument, which contains the list of entities, represented by (dimension, tag)
# pairs. Here we thus translate point 3 (dimension=0, tag=3), by dx=-0.05, dy=0,
# dz=0.
factory.translate([(0, 3)], -0.05, 0, 0)

# The "Duplicata" functionality in .geo files is handled by
# factory.Copy(), which takes a vector of (dim, tag) pairs as input, and
# returns another vector of (dim, tag) pairs.

ov = factory.copy([(0, 3)])

factory.translate(ov, 0, 0.1, 0)

factory.addLine(3, ov[1][2], 7)
factory.addLine(ov[1][2], 5, 8)
factory.addCurveLoop([5,-8,-7,3], 10)
factory.addPlaneSurface([10], 11)

ov = factory.copy([(2, 1), (2, 11)])
factory.translate(ov, 0.12, 0, 0)

println("New surfaces ", ov[1][2], " and ", ov[2][2])

factory.addPoint(0., 0.3, 0.13, lc, 100)
factory.addPoint(0.08, 0.3, 0.1, lc, 101)
factory.addPoint(0.08, 0.4, 0.1, lc, 102)
factory.addPoint(0., 0.4, 0.13, lc, 103)

factory.addLine(4, 100, 110)
factory.addLine(3, 101, 111)
factory.addLine(6, 102, 112)
factory.addLine(5, 103, 113)
factory.addLine(103, 100, 114)
factory.addLine(100, 101, 115)
factory.addLine(101, 102, 116)
factory.addLine(102, 103, 117)

factory.addCurveLoop([115, -111, 3, 110], 118)
factory.addPlaneSurface([118], 119)
factory.addCurveLoop([111, 116, -112, -7], 120)
factory.addPlaneSurface([120], 121)
factory.addCurveLoop([112, 117, -113, -8], 122)
factory.addPlaneSurface([122], 123)
factory.addCurveLoop([114, -110, 5, 113], 124)
factory.addPlaneSurface([124], 125)
factory.addCurveLoop([115, 116, 117, 114], 126)
factory.addPlaneSurface([126], 127)

# The API to create surface loops ("shells") and volumes is similar to the
# one used to create curve loops and surfaces.
factory.addSurfaceLoop([127, 119, 121, 123, 125, 11], 128)
factory.addVolume([128], 129)

# Extrusion works as expected, by providing a vector of (dim, tag) pairs as
# input, the translation vector, and a vector of (dim, tag) pairs as output.
ov2 = factory.extrude([ov[1]], 0, 0, 0.12)

# Mesh sizes associated to geometrical points can be set by passing a vector of
# (dim, tag) pairs for the corresponding points.

factory.mesh.setSize([(0,103), (0,105), (0,109), (0,102), (0,28),
                      (0, 24), (0,6), (0,5)], lc * 3)

model.addPhysicalGroup(3, [129,130], 1)
model.setPhysicalName(3, 1, "The volume")

factory.synchronize()

model.mesh.generate(3)

gmsh.write("t2.msh")

gmsh.finalize()

max_stars_repo_name: Poofee/fastFEM

max_issues_repo_name: Christophe-Foyer/pyFEA

max_forks_repo_name: Christophe-Foyer/pyFEA


text:
module day23

using LinkedLists

wrap(i, len) = i == 0 ? len : i

function playcups(cups, iterations)
    len = length(cups)
    circle = LinkedList{Int}()
    append!(circle, cups)
    lookups = repeat([lastindex(circle)], len)
    for index in keys(circle)
        lookups[getindex(circle, index)] = index
    end

    for i in 1:iterations
        # Put current back at the end of the list
        currentlabel = popfirst!(circle)
        push!(circle, currentlabel)
        lookups[currentlabel] = lastindex(circle)
        # Remove next three
        removed = [popfirst!(circle) for i in 1:3]
        # Find destination
        destlabel = wrap(currentlabel - 1, len)
        while destlabel ∈ removed destlabel = wrap(destlabel - 1, len) end
        destindex = lookups[destlabel]
        # Replace three removed after dest
        for r in removed
            destindex = last(iterate(circle, destindex))
            destindex = insert!(circle, destindex, r)
            lookups[r] = destindex
        end
    end
    cups = [i for i in circle]
end

cups = parse.(Int, collect("247819356"))
cups = playcups(cups, 10)
println("Part 1: remaining = $(join(circshift(cups, -(only(indexin(1, cups)) - 1))[2:end]))")

# Part 2
bigcups = collect(1:1000000)
bigcups[1:9] = parse.(Int, collect("247819356"))
bigcups = playcups(bigcups, 10000000)
println("Part 2: score = $(prod(circshift(bigcups, -(only(indexin(1, bigcups)) - 1))[2:3]))")

end

ext: jl

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
################
# Comment Type #
################

type Comment <: GitHubType
    body::Nullable{GitHubString}
    path::Nullable{GitHubString}
    diff_hunk::Nullable{GitHubString}
    original_commit_id::Nullable{GitHubString}
    commit_id::Nullable{GitHubString}
    id::Nullable{Int}
    original_position::Nullable{Int}
    position::Nullable{Int}
    line::Nullable{Int}
    created_at::Nullable{Dates.DateTime}
    updated_at::Nullable{Dates.DateTime}
    url::Nullable{HttpCommon.URI}
    html_url::Nullable{HttpCommon.URI}
    issue_url::Nullable{HttpCommon.URI}
    pull_request_url::Nullable{HttpCommon.URI}
    user::Nullable{Owner}
end

Comment(data::Dict) = json2github(Comment, data)
Comment(id::Real) = Comment(Dict("id" => id))

namefield(comment::Comment) = comment.id

kind_err_str(kind) = ("Error building comment request: :$kind is not a valid kind of comment.\n"*
                      "The only valid comment kinds are: :issue, :review, :commit")
###############
# API Methods #
###############

function comment(repo, item, kind = :issue; options...)
    if (kind == :issue) || (kind == :pr)
        path = "/repos/$(name(repo))/issues/comments/$(name(item))"
    elseif kind == :review
        path = "/repos/$(name(repo))/pulls/comments/$(name(item))"
    elseif kind == :commit
        path = "/repos/$(name(repo))/comments/$(name(item))"
    else
        error(kind_err_str(kind))
    end
    return Comment(gh_get_json(path; options...))
end

function comments(repo, item, kind = :issue; options...)
    if (kind == :issue) || (kind == :pr)
        path = "/repos/$(name(repo))/issues/$(name(item))/comments"
    elseif kind == :review
        path = "/repos/$(name(repo))/pulls/$(name(item))/comments"
    elseif kind == :commit
        path = "/repos/$(name(repo))/commits/$(name(item))/comments"
    else
        error(kind_err_str(kind))
    end
    results, page_data = gh_get_paged_json(path; options...)
    return map(Comment, results), page_data
end

function create_comment(repo, item, kind = :issue; options...)
    if (kind == :issue) || (kind == :pr)
        path = "/repos/$(name(repo))/issues/$(name(item))/comments"
    elseif kind == :review
        path = "/repos/$(name(repo))/pulls/$(name(item))/comments"
    elseif kind == :commit
        path = "/repos/$(name(repo))/commits/$(name(item))/comments"
    else
        error(kind_err_str(kind))
    end
    return Comment(gh_post_json(path; options...))
end

function edit_comment(repo, item, kind = :issue; options...)
    if (kind == :issue) || (kind == :pr)
        path = "/repos/$(name(repo))/issues/comments/$(name(item))"
    elseif kind == :review
        path = "/repos/$(name(repo))/pulls/comments/$(name(item))"
    elseif kind == :commit
        path = "/repos/$(name(repo))/comments/$(name(item))"
    else
        error(kind_err_str(kind))
    end
    return Comment(gh_patch_json(path; options...))
end

function delete_comment(repo, item, kind = :issue; options...)
    if (kind == :issue) || (kind == :pr)
        path = "/repos/$(name(repo))/issues/comments/$(name(item))"
    elseif kind == :review
        path = "/repos/$(name(repo))/pulls/comments/$(name(item))"
    elseif kind == :commit
        path = "/repos/$(name(repo))/comments/$(name(item))"
    else
        error(kind_err_str(kind))
    end
    return gh_delete(path; options...)
end


text:
function wind_chill(T_a, v)
    T_wc = 0

    # Your code goes here!

    return T_wc
end

max_stars_repo_name: basimr/lovelace

max_issues_repo_name: basimr/project-lovelace

max_forks_repo_name: basimr/project-lovelace


text:
module TestSort
using DataArrays, Base.Test, Compat

dv1 = @data([9, 1, 8, NA, 3, 3, 7, NA])
dv2 = 1.0 * dv1
dv3 = DataArray(collect(1:8))
pdv1 = convert(PooledDataArray, dv1)

@test sortperm(dv1) == sortperm(dv2)
@test sortperm(dv1) == sortperm(pdv1)
@test isequal(sort(dv1), convert(DataArray, sort(dv1)))
@test isequal(sort(dv1), convert(DataArray, sort(pdv1)))

for T in (Float64, BigFloat)
    n = 1000
    na = bitrand(n)
    nna = sum(na)
    a = Array(T, n)
    ra = randn(n-nna)
    a[!na] = ra
    for da in (DataArray(a, na), PooledDataArray(a, na), (pda = PooledDataArray(a, na); setlevels!(pda, shuffle!(pda.pool))))
        @test isequal(sort(da), [DataArray(sort(dropna(da))); DataArray(T, nna)])
        @test isequal(sort(da; lt=(x,y)->isless(x,y)), [DataArray(sort(dropna(da))); DataArray(T, nna)])
        @test isequal(da[sortperm(da)], [DataArray(sort(dropna(da))); DataArray(T, nna)])
        @test isequal(sort(da, rev=true), [DataArray(T, nna); DataArray(sort(dropna(da), rev=true))])
        @test isequal(da[sortperm(da, rev=true)], [DataArray(T, nna); DataArray(sort(dropna(da), rev=true))])
    end
end
end


text:
# State transition function
function POMDPs.transition(mdp::VerticalCAS_MDP, s::stateType, ra::actType)
    h = s[1]; vown = s[2]; vint = s[3]; pra = s[4]; prev_resp = s[5]

    # Computation is faster when using vector of static size
    nextStates = MVector{18, stateType}(undef)
    nextProbs = @MVector(zeros(18))
    ind=1

    # Compute probability that new advisory will be followed
    prob_resp = 0.0
    if (ra==COC) .| (prev_resp .& (pra==ra))
        prob_resp=1.0
    elseif (pra==COC)
        prob_resp = 1.0/(1+5.)
    elseif sameSense(pra,ra)
        prob_resp = 1.0/(1+3.)
    else
        prob_resp = 1.0/(1+5.)
    end

    # Compute probabilities of next states using sigma point sampling
    next_resp = true
    next_pra  = ra
    ownProbs, ownAccels = mdp.accels[ra]
    intProbs, intAccels = mdp.accels[COC]
    for i = 1:3
        for j = 1:3
            next_h,next_vown,next_vint = dynamics(h,vown,vint,ownAccels[i],intAccels[j],ra,mdp)
            nextStates[ind] = (next_h,next_vown,next_vint,next_pra,next_resp)
            nextProbs[ind]  = prob_resp*ownProbs[i]*intProbs[j]
            ind+=1
        end
    end

    if prob_resp<1.0
        prob_notResp = 1.0-prob_resp
        next_resp=false
        ownProbs, ownAccels = mdp.accels[COC]
        for i = 1:3
            for j = 1:3
                next_h,next_vown,next_vint = dynamics(h,vown,vint,ownAccels[i],intAccels[j],COC,mdp)
                nextStates[ind] = (next_h,next_vown,next_vint,next_pra,next_resp)
                nextProbs[ind]  = prob_notResp*ownProbs[i]*intProbs[j]
                ind+=1
            end
        end
    end

    return SparseCat(nextStates,nextProbs)
end

# Dynamic equations
function dynamics(h::Float64,vown::Float64,vint::Float64,ownAccel::Float64, intAccel::Float64, ra::Int, mdp::VerticalCAS_MDP)
    vLow, vHigh = mdp.velRanges[ra]
    if (vLow >= vown) .| (vHigh <= vown)
        ownAccel = 0
    elseif vLow > vown + ownAccel
        ownAccel = vLow-vown
    elseif vHigh < vown + ownAccel
        ownAccel = vHigh-vown
    end

    vLow, vHigh = mdp.velRanges[COC]
    if (vLow >= vint) .| (vHigh <= vint)
        intAccel = 0
    elseif vLow > vint + intAccel
        intAccel = vLow-vint
    elseif vHigh < vint + intAccel
        intAccel = vHigh-vint
    end

    next_h = h-vown-0.5*ownAccel+vint+0.5*intAccel
    next_vown = vown+ownAccel
    next_vint = vint+intAccel
    return next_h,next_vown,next_vint
end

max_stars_repo_name: alisamiii/VerticalCAS

max_issues_repo_name: alisamiii/VerticalCAS

max_forks_repo_name: bethelhall/CoCoCAS


text:
δ=0.10; ρ=0.15; # Param
z= ρ + δ +.02     # Need: z > ρ + δ
I_SS = (z-ρ-δ)/(ρ+δ)
K_SS = I_SS/δ

if 1==1
    a_0(s;δ=δ) = δ*s
    min_s = eps(); max_s = 10.0*K_SS; n_s = 800;
    min_a = eps(); max_a = K_SS; n_a = 200;
    "states:"
    states = range(min_s, max_s, length=n_s)
    ceil_state(s) = states[searchsortedfirst(states, s)]
    "actions:"
    valid_actions()  = range(min_a, max_a, length=n_a)        # possible actions any state.
    valid_actions(s) = filter(>=(-(1-δ)*s), valid_actions())  # valid actions @ state=s
    "Transition:"
    μ(s,a;δ=δ)      = a +(1-δ)*s
    "reward:"
    r(s,a;z=z) = z*s -a -0.5*(a^2)
    "discount:"
    β = 1/(1+ρ);
    using QuickPOMDPs: QuickMDP           #QuickMDP()
    using POMDPModelTools: Deterministic
    m = QuickMDP(
        states     = states,
        actions    = valid_actions,
        transition = (s, a) -> Deterministic( ceil_state( μ(s,a) ) ),
        reward     = (s, a) -> r(s,a),
        discount   = β
    )
    # DiscreteValueIteration: Both work fast!
    using DiscreteValueIteration
    s1 = DiscreteValueIteration.SparseValueIterationSolver()
    s2 = DiscreteValueIteration.ValueIterationSolver()
    @time sol1 = solve(s1, m) #
    @time sol2 = solve(s2, m) #
    # value(sol1, states[2]), action(sol1, states[2])
    # value(sol2, states[2]), action(sol2, states[2])
    value(sol1, states[end])

    using Plots
    # Value
    plot(legend=:bottomright, title="Value Functions");
    #plot!(states[2:end], i->A0 + A1 * i, lab="closed form")
    plot!(states[2:end], i->value(sol1, i),   lab="sol1")
    plot!(states[2:end], i->value(sol2, i),   lab="sol2")

    # Simulation
    Tsim=150; s0=0.5*K_SS; sim1 = []; push!(sim1, s0);
    for tt in 1:Tsim
        s = sim1[tt]
        a = valid_actions()[sol1.policy[searchsortedfirst(states, s)]]
        sp = μ(s,a)
        sp = ceil_state(sp)
        push!(sim1, sp)
    end

    plot(legend=:bottomright, title="Simulation");
    #plot!(simcf,   lab="closed form")
    plot!(sim1,   lab="sol1")
end


text:
using HDF5

"""
`savedaqdata(h5, dev, X; kw...)`

Saves data read by a DAQ device in a HDF5 group

The group where the data will be stored is specified by the parameter `h5`.

The data is stored under the device name of the daq device (`daqdevname(dev)`).

Other parameters related to the data acquisition can be provided using the key word
arguments `kw`. They will be stored under the keyword name.

The user should ensure that the data being stored (including attributes) have types
compatible with the HDF5.jl package.

"""
function savedaqdata(h5, dev, X; kw...)

    devname = daqdevname(dev)
    h5[devname] = X
    d = h5[devname]

    for (k, v) in kw
        attributes(d)[string(k)] = v
    end
    return
end

"""
`savedaqconfig(h5, dev; kw...)`

Save device configuration. Stores information contained in the `conf::DAQConfig`
field of the daq device. Other information can be stored using the keyword parameters
where the keyword names are used as names in the HDF5 files.

The configuration is stored under the device name (`daqdevname(dev)`). User should make
sure that different devices do not have the same name.

The user should ensure that the data being stored ( have types
compatible with the HDF5.jl package.

"""
function savedaqconfig(h5, dev; kw...)

    devname = daqdevname(dev)
    g = create_group(h5, devname)
    g["devname"] = devname
    g["ip"] = daqdevip(dev)
    g["model"] = daqdevmodel(dev)
    g["sn"] = daqdevserialnum(dev)
    g["tag"] = daqdevtag(dev)

    g["iparameters_names"] = collect(keys(iparameters(dev)))
    g["iparameters"] = collect(values(iparameters(dev)))

    g["sparameters_names"] = collect(keys(sparameters(dev)))
    g["sparameters"] = collect(values(sparameters(dev)))

    g["fparameters_names"] = collect(keys(fparameters(dev)))
    g["fparameters"] = collect(values(fparameters(dev)))

    g["daqchannels"] = daqchannels(dev)

    for (k,v) in kw
        g[string(k)] = v
    end

    return
end


text:
function sqarea_triangle(a::T,b::T,c::T) where {T <: Real}
# How to compute (twice) area squared of triangle with
# high precision (Goldberg 1991):
a,b,c=reverse(sort([a,b,c]))
area = (a+(b+c))*(c-(a-b))*(c+(a-b))*(a+(b-c))
return area
end

# Shows the overlap area accuracy for the two formulae:

# 1). Mandel & Agol:

function circle_overlap(r,b)
kap1 = acos((1-r^2+b^2)/(2*b))
kap0 = acos((r^2+b^2-1)/(2*b*r))
akite = sqrt((4*b^2-(1+b^2-r^2)^2)/4)
area = r^2*kap0+kap1-akite

# 2). Luger & Agol:

kite_area2 = sqrt(sqarea_triangle(one(r),b,r))
# Angle of section for occultor:
kap0  = atan2(kite_area2,(r-1)*(r+1)+b^2)
# Angle of section for source:
kap1 = atan2(kite_area2,-((r-1)*(r+1)-b^2))
# Flux of visible uniform disk:
sn = kap1 + r^2*kap0 - 0.5*kite_area2

# Now at high precision:
r_big = big(r); b_big=big(b)
kite_area2_big = sqrt(sqarea_triangle(big(1.0),b_big,r_big))
# Angle of section for occultor:
kap0_big  = atan2(kite_area2,(r_big-1)*(r_big+1)+b_big^2)
# Angle of section for source:
kap1_big = atan2(kite_area2_big,-((r_big-1)*(r_big+1)-b_big^2))
# Flux of visible uniform disk:
sn_big = convert(Float64,kap1_big + r_big^2*kap0_big - 0.5*kite_area2_big)

return area,sn,sn_big
end

using PyPlot
fig,axes = subplots(2,2)

r = 1.0
nb = 1000
db = logspace(-15,-2,nb)
b = 1-r+db
area=zeros(nb); sn=zeros(nb); sn_big=zeros(nb)
for i=1:nb
  area[i],sn[i],sn_big[i] = circle_overlap(r,b[i])
end

ax = axes[1]
ax[:plot](log10.(db),log10.(abs.(pi*r^2-area)),label="r=1.0, MA(2002)",lw=3)
ax[:plot](log10.(db),log10.(abs.(pi*r^2-sn)),label="r=1.0, RA(2018)",".")
ax[:plot](log10.(db),log10.(abs.(pi*r^2-sn_big)),label="r=1.0, BigFloat",linestyle="--")
#ax[:set_xlabel](L"$\log_{10}(b-(1-r))$")
ax[:set_title ](L"$\log_{10}(\pi r^2$-Area of overlap)")
ax[:legend](loc="lower right",fontsize=8)
ax[:axis]([-15,-2,-22,-2])
t1 = linspace(0,2pi,100)
aspect = 1.4*24.0/17.0
x1 = -12.0+cos.(t1); y1=-5.0+aspect*sin.(t1)
x2 = -11.1+.1*cos.(t1); y2=-5.0+0.1*aspect*sin.(t1)
ax[:plot](x1,y1)
ax[:plot](x2,y2)
ax[:plot]([-16,-2],[1,1]*log10(pi*0.1^2/2^53),c="grey",linestyle="-.")

ax = axes[2]
ax[:plot](log10.(db),log10.(abs.(area-sn_big)),label="r=1.0, MA(2002)",lw=3)
ax[:plot](log10.(db),log10.(abs.(sn-sn_big)),label="r=1.0, RA(2018)",".")
ax[:set_xlabel](L"$\log_{10}(b-(1-r))$")
ax[:set_ylabel]("Log Abs(Error in area of overlap)")
ax[:legend](loc="center left",fontsize=8)
ax[:axis]([-15,-2,-25,-5])
ax[:plot]([-16,-2],[1,1]*log10(pi*0.1^2/2^53),c="grey",linestyle="-.")

b = 1+r-db
area=zeros(nb); sn=zeros(nb); sn_big=zeros(nb)
for i=1:nb
  area[i],sn[i],sn_big[i] = circle_overlap(r,b[i])
end

ax = axes[3]
mask = area .> 0.0
ax[:plot](log10.(db[mask]),log10.(area[mask]),label="r=1.0, MA(2002)",lw=3)
ax[:plot](log10.(db),log10.(sn),label="r=1.0, RA(2018)",".")
ax[:plot](log10.(db),log10.(sn_big),label="r=1.0, BigFloat",linestyle="--")
#ax[:set_xlabel](L"$\log_{10}(1+r-b)$")
ax[:set_title](L"$\log_{10}$ ( Area of overlap)")
ax[:legend](loc="lower right",fontsize=8)
ax[:axis]([-15,-2,-22,-2])
x1 = -12.0+cos.(t1); y1=-5.0+aspect*sin.(t1)
x2 = -10.9+0.1*cos.(t1); y2=-5.0+0.1*aspect*sin.(t1)
ax[:plot](x1,y1)
ax[:plot](x2,y2)

ax = axes[4]
ax[:plot](log10.(db),log10.(abs.(area-sn_big)),label="r=1.0, MA(2002)",lw=3)
ax[:plot](log10.(db),log10.(abs.(sn-sn_big)),label="r=1.0, RA(2018)",".")
ax[:set_xlabel](L"$\log_{10}(1+r-b)$")
#ax[:set_ylabel]("Log Abs(Error in area of overlap)")
ax[:legend](loc="center left",fontsize=8)
ax[:axis]([-15,-2,-25,-5])
savefig("area_of_overlap_r10.pdf", bbox_inches="tight")

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
function cart_to_cyl(x,y)
    return [sqrt(x*x + y*y), (atan(y, x) + 2π)%(2π)]
end
function cyl_to_cart(r,t)
    return [r*cos(t), r*sin(t)]
end
function dcart_to_cyl(x,y)
    r2 = x*x + y*y
    r = sqrt(r2)
    return [x/r y/r; -y/r2 x/r2]
end
function dcyl_to_cart(r,t)
    ct, st = cos(t), sin(t)
    return [ct -r*st;st r*ct]
end
function step(x, s, n)
    x_trj = zeros(3, n+1)
    x_trj[:,1] = x
    s0, s1, s2 = s[1], s[2], s[3]
    for i = 2:n+1
           x_trj[:,i] = next(x_trj[:,i-1],s)
    end
    return x_trj
end
function next(u, s)
    x, y, z = u
    s0, s1, s2 = s
    r, t = cart_to_cyl(x, y)
    r1 = s0 + (r - s0)/s1 + cos(t)/2
    t1 = (2*t + s2/4*sin(4*t)) % (2π)
    z1 = z/s1 + sin(t)/2
    x1, y1 = cyl_to_cart(r1,t1)
    return [x1, y1, z1]
end
function dstep(u::Array{Float64,1},s::Array{Float64,1})
    du = zeros(2,2)
    s0, s1, s2 = s
    x, y, z = u
    r, t = cart_to_cyl(x,y)
    r1 = s0 + (r - s0)/s1 + cos(t)/2
    t1 = (2*t + s2/4*sin(4*t)) % (2π)
    z1 = z/s1 + sin(t)/2

    drt1drt = [1/s1 -sin(t)/2;0 2 + s2*cos(4*t)]
    dz1dt = cos(t)/2
    dz1dz = 1/s1

    drtdxy = dcart_to_cyl(x,y)
    dtdx, dtdy = drtdxy[2,1], drtdxy[2,2]
    dxy1dxy = dcyl_to_cart(r1,t1)*drt1drt*drtdxy
    dz1 = [dz1dt*dtdx dz1dt*dtdy dz1dz]

    du = [dxy1dxy [0;0]]
    du = [du;dz1]
    return du
end
function pert(u::Array{Float64,1}, s::Array{Float64,1},
              p::Int64)
    x, y, z = u
    s0, s1, s2 = s
    r, t = cart_to_cyl(x,y)
    r1 = s0 + (r - s0)/s1 + cos(t)/2
    t1 = 2*t + s2/4*sin(4*t)
    z1 = z/s1 + sin(t)/2
    x1, y1 = cyl_to_cart(r1,t1)
    dxyz1drtz1 = [dcyl_to_cart(r1,t1) [0;0]]
    dxyz1drtz1 = [dxyz1drtz1; [0 0 1]]

    if p==1
        drtz1ds = [1-1/s1; 0; 0]
    elseif p==2
        drtz1ds = [-(r - s0)/s1^2; 0; -z/s1^2]
    elseif p==3
        drtz1ds = [0; (1/4)*sin(4*t); 0]
    else
        println("parameter perturbation is only defined
                for indices 1 through 3")
        return 0
    end
    return dxyz1drtz1*drtz1ds
end

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
export train_proj!
function train_proj!{T<:Number}(net::HopfieldNetwork, patterns::Array{T,2})
	nn,np = size(patterns)

	# calc synaptic weight (Projection learning rule)
	Qinv=pinv(patterns' * patterns/nn)
	for p in 1:nn
		for q in p+1:nn
			s = zero(T)
			for mu in 1:np
				for nu in 1:np
					s += patterns[p, mu] * Qinv[mu,nu] * conj(patterns[q, nu])
				end
			end
			s /= nn
			net.weight[p, q] = s
			net.weight[q, p] = conj(s)
		end
	end

	return
end


text:
# Run the block sparse benchmarks using the specified number of Julia threads and write the results
# By default, a verbose output is written to a log file
for n in 1:4
  run(`julia -t $n -e 'using ITensorBenchmarks; runbenchmarks(write_results = true, benchmarks = ["dmrg_1d_qns", "dmrg_2d_qns", "dmrg_2d_conserve_ky"], blocksparse_num_threads = Threads.nthreads())'`)
end


text:
# ------------------------------------
# calculate right hand side
# ------------------------------------
function setup_RHS(RHS, cellxmax, cellymax, cellzmax, E_adv_hat, F_adv_hat, G_adv_hat,
                    E_vis_hat, F_vis_hat, G_vis_hat, nval, volume, icell)
    Threads.@threads for l in 1:nval
        for k in 1+icell:cellzmax-icell
            for j in 1+icell:cellymax-icell
                for i in 1+icell:cellxmax-icell
                    RHS[i,j,k,l] = - (E_adv_hat[i+1,  j,  k, l] - E_vis_hat[i+1,  j,  k, l]) +
                                     (E_adv_hat[  i,  j,  k, l] - E_vis_hat[  i,  j,  k, l]) -
                                     (F_adv_hat[  i,j+1,  k, l] - F_vis_hat[  i,j+1,  k, l]) +
                                     (F_adv_hat[  i,  j,  k, l] - F_vis_hat[  i,  j,  k, l]) -
                                     (G_adv_hat[  i,  j,k+1, l] - G_vis_hat[  i,  j,k+1, l]) +
                                     (G_adv_hat[  i,  j,  k, l] - G_vis_hat[  i,  j,  k, l])
                end
            end
        end
    end

    #=
    println("RHS[20,50,:]")
    i = 50
    j = 50
    println(RHS[i,j,2])
    println((E_adv_hat[i+1,j,2] - E_adv_hat[i,j,2])/RHS[i,j,2])
    println((E_vis_hat[i+1,j,2] - E_vis_hat[i,j,2])/RHS[i,j,2])
    println((F_adv_hat[i,j+1,2] - F_adv_hat[i,j,2])/RHS[i,j,2])
    println((F_vis_hat[i,j+1,2] - F_vis_hat[i,j,2])/RHS[i,j,2])
    =#

    return RHS
end


text:
# Filter some noise, try to learn the coefficients
@testset "$method" for method in (LPCBurg(), LPCLevinson())
    coeffs = [1, .5, .3, .2]
    x = filt([1], coeffs, randn(50000))

    # Analyze the filtered noise, attempt to reconstruct the coefficients above
    ar, e = lpc(x[1000:end], length(coeffs)-1, method)

    @test all(abs.([1; ar] .- coeffs) .<= 1e-2)
end


text:
using Documenter, IndexedGraphs
makedocs(sitename="IndexedGraphs Documentation",
    # format = Documenter.HTML(prettyurls = false),
    pages = [
        "Home" => "index.md",
        "Graph types" => [
            "graph.md",
            "bidigraph.md",
            "digraph.md",
            "factor.md"
        ],
        "Reference" => "reference.md"
    ]
    )

deploydocs(
    repo = "github.com/stecrotti/IndexedGraphs.jl.git",
)


text:
function array_split(N, size, comm)

    counts = zeros(Int64, size)
    displs = zeros(Int64, size)

    if rank == 0
        counts_guess = Int64(div(N, size, RoundDown))
        Remainder    = Int64(N - counts_guess*size)
        counts[:]   .= counts_guess

        for i in 1:Remainder
            counts[i] += 1
        end

        displs[:] = cumsum(append!([0], counts))[1:size]
    end

    MPI.Bcast!(counts, 0, comm)
    MPI.Bcast!(displs, 0, comm)

    return counts, displs
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'dre'_quill

"""
    'dre'_quill(;kwargs...)

A Quill component.
ExampleComponent is an example component.
It takes a property, `label`, and
displays it.
It renders an input with the property `value`
which is editable by the user.
Keyword arguments:
- `id` (String; optional): The ID used to identify this component in Dash callbacks.
slateContent={SlateRTE.deserializeHTMLString(value)}
- `charCount` (Real; optional): The number of charaters in the editor (excl HTML)
- `maxLength` (Real; optional): The value displayed in the input.
- `modules` (Dict; optional): The toolbar options modules.
Should be {'toolbar':[list of options]}
- `value` (String; optional): The value displayed in the input.
"""
function 'dre'_quill(; kwargs...)
        available_props = Symbol[:id, :charCount, :maxLength, :modules, :value]
        wild_props = Symbol[]
        return Component("'dre'_quill", "Quill", "dash_quill", available_props, wild_props; kwargs...)
end

max_stars_repo_name: pandamodium/dash-quill

max_issues_repo_name: pandamodium/dash-quill

max_forks_repo_name: pandamodium/dash-quill


text:
using PerformanceTestTools
using PerformanceTestTools._Testing: XFailTestSet
using Test

@testset "expected pass" begin
    PerformanceTestTools.@include("__test_ir.jl")
    PerformanceTestTools.@include_foreach(
        "__test_a_eq_b.jl",
        [["A" => "1", "B" => "1"], ["A" => "2", "B" => "2", `--compile=min`]],
    )
    PerformanceTestTools.@include_foreach(
        "__test_a_eq_b.jl",
        [["A" => "1", "B" => "1"], ["A" => "2", "B" => "2", `--compile=min`]],
        parallel = true,
    )
end

@testset "expected failures" begin
    @info "=== Starting xfail tests ==="
    @testset for parallel in [false, true]
        @testset XFailTestSet "only envs" begin
            PerformanceTestTools.@include_foreach(
                "__test_a_eq_b.jl",
                [["A" => "1", "B" => "2"]],
                parallel = parallel,
            )
        end
    end
    @testset XFailTestSet "envs and opts" begin
        PerformanceTestTools.@include_foreach(
            "__test_a_eq_b.jl",
            [["A" => "1", "B" => "2", `--compile=min`]],
            parallel = true,
        )
    end
    @testset XFailTestSet "with in-process" begin
        withenv("A" => "1", "B" => "2") do
            PerformanceTestTools.@include_foreach(
                "__test_a_eq_b.jl",
                [nothing, [`--compile=min`]],
                parallel = true,
            )
        end
    end
    @info "=== Finished xfail tests ==="
end


text:
"""
    PyArray{R,N,T,M,L}(o) :: AbstractArray{T,N}

Interpret the Python array `o` as a Julia array.

The input may be anything supporting the buffer protocol or the numpy array interface.
This includes `bytes`, `bytearray`, `array.array`, `numpy.ndarray`, `pandas.Series`.

All type parameters are optional:
- `R` is the type of elements of the underlying buffer.
- `N` is the number of dimensions.
- `T` is the element type.
- `M` is true if the array is mutable.
- `L` is true if the array supports fast linear indexing.

There are alias types with names of the form `Py[Mutable/Immutable/][Linear/Cartesian/][Array/Vector/Matrix]`.
"""
mutable struct PyArray{R,N,T,M,L} <: AbstractArray{T,N}
    ref::PyRef
    ptr::Ptr{R}
    size::NTuple{N,Int}
    length::Int
    bytestrides::NTuple{N,Int}
    handle::Any
end
export PyArray

# aliases
for M in (true, false, missing)
    for L in (true, false, missing)
        for N in (1, 2, missing)
            name = Symbol("Py", M===missing ? "" : M ? "Mutable" : "Immutable", L===missing ? "" : L ? "Linear" : "Cartesian", N===missing ? "Array" : N==1 ? "Vector" : "Matrix")
            name == "PyArray" && continue
            lparams = filter(x->x!==nothing, (:R, N===missing ? :N : nothing, :T, M===missing ? :M : nothing, L===missing ? :L : nothing))
            rparams = (:R, N===missing ? :N : N, :T, M===missing ? :M : M, L===missing ? :L : L)
            @eval const $name{$(lparams...)} = PyArray{$(rparams...)}
            @eval export $name
        end
    end
end

ispyreftype(::Type{<:PyArray}) = true
pyptr(x::PyArray) = pyptr(x.ref)
Base.unsafe_convert(::Type{CPyPtr}, x::PyArray) = pyptr(x.ref)
C.PyObject_TryConvert__initial(o, ::Type{T}) where {T<:PyArray} =
    CTryConvertRule_trywrapref(o, T)

(::Type{A})(o; opts...) where {A<:PyArray} = begin
    ref = PyRef(o)
    info = pyarray_info(ref; opts...)
    R = pyarray_R(A, info)
    N = pyarray_N(A, info)
    T = pyarray_T(A, info, R)
    M = pyarray_M(A, info)
    L = pyarray_L(A, info, Val(N))
    size = pyarray_size(info, Val(N))
    PyArray{R,N,T,M,L}(
        ref,
        Ptr{R}(pyarray_ptr(info)),
        size,
        N == 0 ? 1 : prod(size),
        pyarray_bytestrides(info, Val(N)),
        pyarray_handle(info),
    )
end

pyarray_R(::Type{A}, info) where {R, A<:PyArray{R}} =
    !isa(R, Type) ? error("R must be a type, got R=$R") :
    !allocatedinline(R) ? error("source elements must be allocated inline, got R=$R") :
    sizeof(R) != pyarray_elsize(info) ? error("source elements must have size $(pyarray_elsize(info)), got R=$R") : R
pyarray_R(::Type{A}, info) where {A<:PyArray} = pyarray_R(info)::Type

pyarray_N(::Type{A}, info) where {N,A<:(PyArray{R,N} where R)} =
    !isa(N, Int) ? error("N must be an Int") :
    N != pyarray_N(info) ? error("source dimension is $(pyarray_N(info)), got N=$N") : N
pyarray_N(::Type{A}, info) where {A<:PyArray} = pyarray_N(info)::Int

pyarray_T(::Type{A}, info, ::Type{R}) where {T,A<:(PyArray{R,N,T} where {R,N}),R} =
    !isa(T, Type) ? error("T must be a type, got T=$T") : T
pyarray_T(::Type{A}, info, ::Type{R}) where {A<:PyArray,R} = pyarray_default_T(R)

pyarray_M(::Type{A}, info) where {M,A<:(PyArray{R,N,T,M} where {R,N,T})} =
    !isa(M, Bool) ? error("M must be a bool, got M=$M") :
    (M && !pyarray_M(info)) ? error("source is immutable, got M=$M") : M
pyarray_M(::Type{A}, info) where {A<:PyArray} = pyarray_M(info)

pyarray_L(::Type{A}, info, ::Val{N}) where {L,A<:(PyArray{R,N,T,M,L} where {R,N,T,M}),N} =
    !isa(L, Bool) ? error("L must be a bool, got L=$L") :
    (L && !pyarray_L(info, Val(N))) ? error("source is not linearly indexable, got L=$L") : L
pyarray_L(::Type{A}, info, ::Val{N}) where {A<:PyArray, N} = pyarray_L(info, Val(N))

abstract type PyArrayInfo end

pyarray_L(info::PyArrayInfo, ::Val{N}) where {N} = begin
    N ≤ 1 && return true
    strides = pyarray_bytestrides(info, Val(N))
    size = pyarray_size(info, Val(N))
    strides == size_to_fstrides(strides[1], size...)
end

struct PyArrayInfo_Buffer <: PyArrayInfo
    buf :: PyBuffer
end

pyarray_R(info::PyArrayInfo_Buffer) = info.buf.eltype
pyarray_N(info::PyArrayInfo_Buffer) = Int(info.buf.ndim)
pyarray_M(info::PyArrayInfo_Buffer) = !info.buf.readonly
pyarray_ptr(info::PyArrayInfo_Buffer) = info.buf.buf
pyarray_size(info::PyArrayInfo_Buffer, ::Val{N}) where {N} = begin
    p = info.buf.info[].shape
    ntuple(i->Int(unsafe_load(p,i)), Val(N))::NTuple{N,Int}
end
pyarray_bytestrides(info::PyArrayInfo_Buffer, ::Val{N}) where {N} = begin
    p = info.buf.info[].strides
    ntuple(i->Int(unsafe_load(p,i)), Val(N))::NTuple{N,Int}
end
pyarray_handle(info::PyArrayInfo_Buffer) = info.buf
pyarray_elsize(info::PyArrayInfo_Buffer) = info.buf.itemsize

struct PyArrayInfo_ArrayInterface <: PyArrayInfo
    dict :: PyObject
    ptr :: Ptr{Cvoid}
    mutable :: Bool
    handle :: Any
end

pyarray_R(info::PyArrayInfo_ArrayInterface) = begin
    # converts descr to something more Julian
    pyconvertdescr(x) = begin
        @py ```
        def convert(x):
            def fix(x):
                a = x[0]
                a = (a, a) if isinstance(a, str) else (a[0], a[1])
                b = x[1]
                c = x[2] if len(x)>2 else 1
                return (a, b, c)
            if x is None or isinstance(x, str):
                return x
            else:
                return [fix(y) for y in x]
        $(r::Union{Nothing,String,Vector{Tuple{Tuple{String,String}, PyObject, Int}}}) = convert($x)
        ```
        r isa Vector ? [(a, pyconvertdescr(b), c) for (a, b, c) in r] : r
    end
    typestr = @pyv `$(info.dict)["typestr"]`::String
    descr = pyconvertdescr(@pyv `$(info.dict).get("descr")`)
    pytypestrdescr_to_type(typestr, descr)
end
pyarray_N(info::PyArrayInfo_ArrayInterface) = @pyv `len($(info.dict)["shape"])`::Int
pyarray_M(info::PyArrayInfo_ArrayInterface) = info.mutable
pyarray_ptr(info::PyArrayInfo_ArrayInterface) = info.ptr
pyarray_size(info::PyArrayInfo_ArrayInterface, ::Val{N}) where {N} = @pyv `$(info.dict)["shape"]`::NTuple{N,Int}
pyarray_bytestrides(info::PyArrayInfo_ArrayInterface, ::Val{N}) where {N} = begin
    strides = @pyv `$(info.dict).get("strides")`::Union{Nothing,NTuple{N,Int}}
    if strides === nothing
        size = pyarray_size(info, Val(N))
        elsize = pyarray_elsize(info)
        size_to_cstrides(elsize, size...)
    else
        strides
    end
end
pyarray_handle(info::PyArrayInfo_ArrayInterface) = info.handle
pyarray_elsize(info::PyArrayInfo_ArrayInterface) = begin
    sz = @pya ```
    s = $(info.dict)["typestr"]
    c = s[1]
    if c == "O":
        ans = -1
    elif c == "m" or c == "M":
        ans = int(s[2:].split("[")[0])
    else:
        ans = int(s[2:])
    ```::Int
    sz < 0 ? sizeof(CPyPtr) : sz
end

pyarray_info(ref; buffer::Bool=true, array::Bool=true, copy::Bool=true) = begin
    # __array_interface__
    if array
        @py ```
        try:
            d = $ref.__array_interface__
            assert d["version"] == 3
            assert d.get("mask") is None
        except:
            d = None
        $(d::Union{Nothing,PyObject}) = d
        ```
        if d !== nothing
            @py ```
            data = $d.get("data", $ref)
            if isinstance(data, tuple):
                ptr = data[0]
                mut = not data[1]
                buf = None
            else:
                ptr = 0
                mut = False
                buf = data
            $(ptr::UInt) = ptr
            $(mut::Bool) = mut
            $(buf::Union{PyRef,Nothing}) = buf
            ```
            if buf === nothing
                return PyArrayInfo_ArrayInterface(d, Ptr{Cvoid}(ptr), mut, (ref, d))
            else
                b = PyBuffer(buf, C.PyBUF_RECORDS_RO)
                offset = @pyv `$d.get("offset", 0)`::Int
                return PyArrayInfo_ArrayInterface(d, b.buf+offset, !b.readonly, (ref, d, b))
            end
        end
    end
    # __array_struct__
    if array
        @py ```
        d = None
        try:
            d = $ref.__array_struct__
        except:
            pass
        $(d::Union{Nothing,PyObject}) = d
        ```
        d === nothing || error("not implemented: creating PyArray from __array_struct__")
    end
    # buffer protocol
    if buffer
        b = PyBuffer()
        err = C.PyObject_GetBuffer(ref, pointer(b.info), C.PyBUF_RECORDS_RO)
        if ism1(err)
            C.PyErr_Clear()
        else
            return PyArrayInfo_Buffer(b)
        end
    end
    # __array__
    if array && copy
        @py ```
        a = None
        try:
            a = $ref.__array__()
        except:
            pass
        $(a::Union{Nothing,PyRef}) = a
        ```
        a === nothing || return pyarray_info(a, buffer=buffer, array=array, copy=false)
    end
    error("this object cannot be interpreted as a strided array")
end

Base.size(x::PyArray) = x.size
Base.length(x::PyArray) = x.length
ismutablearray(::PyArray{R,N,T,M,L}) where {R,N,T,M,L} = M
Base.IndexStyle(::Type{PyArray{R,N,T,M,L}}) where {R,N,T,M,L} =
    L ? Base.IndexLinear() : Base.IndexCartesian()

Base.@propagate_inbounds Base.getindex(
    x::PyArray{R,N,T,M,L},
    i::Vararg{Int,N},
) where {R,N,T,M,L} = pyarray_getindex(x, i...)
Base.@propagate_inbounds Base.getindex(
    x::PyArray{R,N,T,M,true},
    i::Int
) where {R,N,T,M} = pyarray_getindex(x, i)
Base.@propagate_inbounds Base.getindex(
    x::PyArray{R,1,T,M,true},
    i::Int
) where {R,T,M} = pyarray_getindex(x, i)

Base.@propagate_inbounds pyarray_getindex(x::PyArray, i...) = begin
    @boundscheck checkbounds(x, i...)
    pyarray_load(eltype(x), x.ptr + pyarray_offset(x, i...))
end

Base.@propagate_inbounds Base.setindex!(
    x::PyArray{R,N,T,true,L},
    v,
    i::Vararg{Int,N},
) where {R,N,T,L} = pyarray_setindex!(x, v, i...)
Base.@propagate_inbounds Base.setindex!(
    x::PyArray{R,N,T,true,true},
    v,
    i::Int
) where {R,N,T} = pyarray_setindex!(x, v, i)
Base.@propagate_inbounds Base.setindex!(
    x::PyArray{R,1,T,true,true},
    v,
    i::Int
) where {R,T} = pyarray_setindex!(x, v, i)

Base.@propagate_inbounds pyarray_setindex!(x::PyArray, v, i...) = begin
    @boundscheck checkbounds(x, i...)
    pyarray_store!(x.ptr + pyarray_offset(x, i...), convert(eltype(x), v))
    x
end

pyarray_default_T(::Type{R}) where {R} = R
pyarray_default_T(::Type{C.PyObjectRef}) = PyObject

pyarray_load(::Type{T}, p::Ptr{T}) where {T} = unsafe_load(p)
pyarray_load(::Type{T}, p::Ptr{C.PyObjectRef}) where {T} = begin
    o = unsafe_load(p).ptr
    isnull(o) && throw(UndefRefError())
    ism1(C.PyObject_Convert(o, T)) && pythrow()
    takeresult(T)
end

pyarray_store!(p::Ptr{T}, v::T) where {T} = unsafe_store!(p, v)
pyarray_store!(p::Ptr{C.PyObjectRef}, v::T) where {T} = begin
    o = C.PyObject_From(v)
    isnull(o) && pythrow()
    C.Py_DecRef(unsafe_load(p).ptr)
    unsafe_store!(p, C.PyObjectRef(o))
end

pyarray_offset(x::PyArray{R,N,T,M,true}, i::Int) where {R,N,T,M} =
    N == 0 ? 0 : (i - 1) * x.bytestrides[1]

pyarray_offset(x::PyArray{R,1,T,M,true}, i::Int) where {R,T,M} = (i - 1) .* x.bytestrides[1]

pyarray_offset(x::PyArray{R,N}, i::Vararg{Int,N}) where {R,N} =
    sum((i .- 1) .* x.bytestrides)


text:
using FRGRealTime.MathFuns

function massflow1(k, m2, T, lambda, Nf=4)
    (1 / 6) * Nf * lambda * TF.Fb1(k, m2, T)
end

function massflow2(k, m2, T, lambda, Nf=4)
    (1 / 6) * (Nf + 2) * lambda * TF.Fb1(k, m2, T)
end

function lambdaflow1(k, m2, T, lambda)
    -(1 / 3) * lambda^2 * TF.Fb2(k, m2, T)
end

function lambdaflow2(k, m2, T, lambda, Nf=4)
    -(1 / 6) * lambda^2 * (Nf + 8) * TF.Fb2(k, m2, T)
end

function Imlambdaflow1(p0, k, m2, lambda, T, Imlambda, Relambda, ϵ=0.1)
    -(1 / 3) * (
        -Imlambda^2 * TF.ImFb2(p0, k, m2, T, ϵ) +
        2 * TF.ReFb2(p0, k, m2, T, ϵ) * Relambda * Imlambda +
        TF.ImFb2(p0, k, m2, T, ϵ) * Relambda^2
    )
end

function Relambdaflow1(p0, k, m2, lambda, T, Imlambda, Relambda, ϵ=0.1)
    -(1 / 3) * (
        -Imlambda^2 * TF.ReFb2(p0, k, m2, T, ϵ) -
        2 * TF.ImFb2(p0, k, m2, T, ϵ) * Relambda * Imlambda +
        TF.ReFb2(p0, k, m2, T, ϵ) * Relambda^2
    )
end

function Imlambdaflow2(p0, k, m2, lambda, T, Imlambda, Relambda, ϵ=0.1)
    -(2 / 3) * (
        -Imlambda^2 * TF.ImFb2(p0, k, m2, T, ϵ) +
        2 * TF.ReFb2(p0, k, m2, T, ϵ) * Relambda * Imlambda +
        TF.ImFb2(p0, k, m2, T, ϵ) * Relambda^2
    )
end

function Relambdaflow2(
    p0,
    k,
    m2,
    lambda,
    T,
    Imlambda,
    Relambda,
    ϵ=0.1,
    Nf=4,
)
    -(2 / 3) * (
        -Imlambda^2 * TF.ReFb2(p0, k, m2, T, ϵ) -
        2 * TF.ImFb2(p0, k, m2, T, ϵ) * Relambda * Imlambda +
        TF.ReFb2(p0, k, m2, T, ϵ) * Relambda^2
    ) - (1 / 6) * lambda^2 * (Nf + 4) * TF.Fb2(k, m2, T)
end

function Imlambdaflow3(p0, q0, k, m2, lambda, T, Imlambda, Relambda, ϵ=0.1)
    -(1 / 3) * (
        -Imlambda^2 * TF.ImFb2(p0 + q0, k, m2, T, ϵ) +
        2 *
        TF.ReFb2(p0 + q0, k, m2, T, ϵ) *
        Relambda *
        Imlambda +
        TF.ImFb2(p0 + q0, k, m2, T, ϵ) * Relambda^2
    ) -
    (1 / 3) * (
        -Imlambda^2 * TF.ImFb2(p0 - q0, k, m2, T, ϵ) +
        2 *
        TF.ReFb2(p0 - q0, k, m2, T, ϵ) *
        Relambda *
        Imlambda +
        TF.ImFb2(p0 - q0, k, m2, T, ϵ) * Relambda^2
    )
end

function Imlambdaflow3_ac(p0, q0, k, m2, lambda, T, Imlambda, Relambda, ϵ=0.1)
    -(1 / 3) * (
        -Imlambda^2 * TF.ImFb2(p0 + q0, k, m2, T, ϵ) +
        2 *
        TF.ReFb2_ac(p0 + q0, k, m2, T, ϵ) *
        Relambda *
        Imlambda +
        TF.ImFb2(p0 + q0, k, m2, T, ϵ) * Relambda^2
    ) -
    (1 / 3) * (
        -Imlambda^2 * TF.ImFb2(p0 - q0, k, m2, T, ϵ) +
        2 *
        TF.ReFb2_ac(p0 - q0, k, m2, T, ϵ) *
        Relambda *
        Imlambda +
        TF.ImFb2(p0 - q0, k, m2, T, ϵ) * Relambda^2
    )
end
function Relambdaflow3_ac(
    p0,
    q0,
    k,
    m2,
    lambda,
    T,
    Imlambda,
    Relambda,
    ϵ=0.1,
    Nf=4,
)
    -(1 / 3) * (
        -Imlambda^2 * TF.ReFb2_ac(p0 + q0, k, m2, T, ϵ) -
        2 *
        TF.ImFb2(p0 + q0, k, m2, T, ϵ) *
        Relambda *
        Imlambda + TF.ReFb2_ac(p0 + q0, k, m2, T, ϵ) * Relambda^2
    ) -
    (1 / 3) * (
        -Imlambda^2 * TF.ReFb2_ac(p0 - q0, k, m2, T, ϵ) -
        2 *
        TF.ImFb2(p0 - q0, k, m2, T, ϵ) *
        Relambda *
        Imlambda + TF.ReFb2_ac(p0 - q0, k, m2, T, ϵ) * Relambda^2
    ) - (1 / 6) * lambda^2 * (Nf + 4) * TF.Fb2(k, m2, T)
end

function Imlambdaflow4(p0, q0, k, m2, lambda, T, Imlambda, Relambda, ϵ = 0.1)
    -(1 / 3) * (TF.ImFb2(p0 + q0, k, m2, T, ϵ) * Relambda^2) -
    (1 / 3) * (TF.ImFb2(p0 - q0, k, m2, T, ϵ) * Relambda^2)
end

function Imlambdaflow4_ac(p0, q0, k, m2, lambda, T, Imlambda, Relambda, ϵ = 0.1)
    -(1 / 3) * (TF.ImFb2(p0 + q0, k, m2, T, ϵ) * Relambda^2) -
    (1 / 3) * (TF.ImFb2(p0 - q0, k, m2, T, ϵ) * Relambda^2)
end

function Relambdaflow4(
    p0,
    q0,
    k,
    m2,
    lambda,
    T,
    Imlambda,
    Relambda,
    ϵ = 0.1,
    Nf = 4,
)
    -(1 / 3) * (TF.ReFb2(p0 + q0, k, m2, T, ϵ) * Relambda^2) -
    (1 / 3) * (TF.ReFb2(p0 - q0, k, m2, T, ϵ) * Relambda^2) -
    (1 / 6) * lambda^2 * (Nf + 4) * TF.Fb2(k, m2, T)
end

function Relambdaflow4_ac(
    p0,
    q0,
    k,
    m2,
    lambda,
    T,
    Imlambda,
    Relambda,
    ϵ = 0.1,
    Nf = 4,
)
    -(1 / 3) * (TF.ReFb2_ac(p0 + q0, k, m2, T, ϵ) * Relambda^2) -
    (1 / 3) * (TF.ReFb2_ac(p0 - q0, k, m2, T, ϵ) * Relambda^2) -
    (1 / 6) * lambda^2 * (Nf + 4) * TF.Fb2(k, m2, T)
end

function twopointflow(p0, k, m2, T, lambda)
    -1 / 48 * (
        k^4 *
        csch(sqrt(k^2 + m2) / (2 * T))^2 *
        (
            lambda(k, abs(sqrt(k^2 + m2) - p0)) +
            lambda(k, abs(sqrt(k^2 + m2) + p0))
        ) *
        (sqrt(k^2 + m2) + T * sinh(sqrt(k^2 + m2) / T))
    ) / ((k^2 + m2)^(3 / 2) * pi^2 * T)
end

function twopointflow2(p0, k, m2, T, lambda)
    -1 / 24 * (
        k^4 *
        csch(sqrt(k^2 + m2) / (2 * T))^2 *
        lambda(k) *
        (sqrt(k^2 + m2) + T * sinh(sqrt(k^2 + m2) / T))
    ) / ((k^2 + m2)^(3 / 2) * pi^2 * T)
end

function twopointflow3(p0, k, m2, T, lambda)
    -1 / 24 * (
        k^4 *
        csch(sqrt(k^2 + m2) / (2 * T))^2 *
        lambda(k) *
        (sqrt(k^2 + m2) + T * sinh(sqrt(k^2 + m2) / T))
    ) / ((k^2 + m2)^(3 / 2) * pi^2 * T)
end

function RelambdaflowMktest(
    p0,
    q0,
    k,
    m2,
    lambda,
    T,
    Imlambda,
    RelambdaMk,
    ϵ = 0.1,
    Nf = 4,
)
    -(1 / 3) * (ReFb2Mkpositivetest(p0, k, m2, T, ϵ) * RelambdaMk^2)
end

function ReFb2Mkpositivetest(p0, k, m2, T, ϵ = 0.1)
    -1 / 48 * (k^4 * (ϵ - 4 * p0 * Eb(k, m2) + 8 * Eb(k, m2)^2)) / (
        pi^2 *
        Eb(k, m2)^3 *
        (p0 + 2 * Eb(k, m2)) *
        (p0^2 + ϵ + 4 * Eb(k, m2) * (-p0 + Eb(k, m2)))
    )
end


text:
# Functions for time-domain modeling and inversion using devito
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January, 2017
# Updated, December 2020, Mathias Louboutin, mlouboutin3@gatech.edu

#############################################################################
# Containers
include("Types/ModelStructure.jl")    # model container
include("Types/GeometryStructure.jl") # source or receiver setup, recording time and sampling
include("Types/OptionsStructure.jl")

#############################################################################
# Abstract vectors
include("Types/abstract.jl")
include("Types/broadcasting.jl")
include("Types/judiWavefield.jl") # dense RHS (wavefield)
include("Types/judiWeights.jl")    # Extended source weight vector
include("Types/judiVector.jl")    # Julia data container
include("Types/judiComposites.jl")    # A composite type to work with hcat/vcat of judi types

# Utility types
const SourceType{T} = Union{Vector{T}, judiMultiSourceVector{T}, PhysicalParameter{T}}
const dmType{T} = Union{Vector{T}, PhysicalParameter{T}}

#############################################################################
# Utils
include("Utils/auxiliaryFunctions.jl")

#############################################################################
# Linear operators
include("LinearOperators/basics.jl")
include("LinearOperators/lazy.jl")
include("LinearOperators/operators.jl")

#############################################################################
# PDE solvers
include("Modeling/distributed.jl") # Modeling functions utilities
include("Modeling/python_interface.jl")  # forward/adjoint linear/nonlinear modeling
include("Modeling/time_modeling_serial.jl")  # forward/adjoint linear/nonlinear modeling
include("Modeling/misfit_fg.jl")  # FWI/LSRTM objective function value and gradient
include("Modeling/twri_objective.jl")  # TWRI objective function value and gradient
include("Modeling/propagation.jl")

#############################################################################
# Preconditioners and optimization
include("Utils/seismic_preconditioners.jl")

#############################################################################
if VERSION>v"1.2"
  function (F::judiPropagator)(m::Model)
    Fl = deepcopy(F)
    Fl.model.n = m.n
    Fl.model.d = m.d
    Fl.model.o = m.o
    for (k, v) in m.params
      Fl.model.params[k] = v
    end
    Fl
  end

  function (F::judiPropagator)(;kwargs...)
    Fl = deepcopy(F)
    for (k, v) in kwargs
      k in keys(Fl.model.params) && Fl.model.params[k] .= v
    end
    Fl
  end
end


text:
"""
    `get_rew_rate(rewards, alpha, initialvalue)`

Use accumulate to update reward rate according to
a starting point 'initialvalue'
a sequence of outcomes 'rewards'
and a learning rate 'alpha'
"""

function get_rew_rate(rewards, alpha, initialvalue)
   return accumulate(rewards; init=initialvalue) do value, reward
       #=(1 - alpha) * value + alpha*reward=#
       value + alpha*(reward - value)
   end
end

"""
    `get_consecutive_omissions_rate(n, alpha, initialvalue)`

Calculate rew_rate update for n consecutive no rewards
"""
function get_consecutive_omissions_rate(n, alpha, initialvalue)
    initialvalue*(1-alpha)^n
end

"""
    `get_consecutive_rewards_rate(n, alpha, initialvalue)`

Calculate rew_rate update for n consecutive no rewards
"""

function get_consecutive_rewards_rate(n, alpha, initialvalue)
    initialvalue*((1-alpha)^n) + alpha*(ploynomial_alpha(1-alpha,n))
end

"""
    `ploynomial_alpha(n, alpha, initialvalue)`

    Helper function to calculate the value of the polynomial necessary to calculate consecutive rewards rate
"""
function ploynomial_alpha(oneminusalpha,n)
    n<1 && error("minimum one outcome to consider")
    sum(oneminusalpha.^(0:n-1))
end

function get_rew_rate(rewards,n, alpha, initialvalue)
    x = collect(zip(rewards,n))
    return accumulate(x; init=initialvalue) do value, x
        Bool(x[1]) ? get_consecutive_rewards_rate(x[2], alpha, value) : get_consecutive_omissions_rate(x[2], alpha, value)
    end
end

# rew = [0,0,1,1,0,0,0,1,1,1]
# alpha_val = 0.1
# initialvalue = 0.5
# get_rew_rate(rew, alpha_val, initialvalue)
# get_rew_rate([0,1,0,1],[2,2,3,3],alpha_val,initialvalue)


text:
using Sundials

##
## Example problem for IDA: 2D heat equation, serial, banded.
##
## This example solves a discretized 2D heat equation problem.  This
## version uses the band solver IDABand, and IDACalcIC.
##
## The DAE system solved is a spatial discretization of the PDE
##
##          du/dt = d^2u/dx^2 + d^2u/dy^2
##
## on the unit square. The boundary condition is u = 0 on all edges.
## Initial conditions are given by u = 48 x (1 - x) y (1 - y).  The
## PDE is treated with central differences on a uniform M x M
## grid. The values of u at the interior points satisfy ODEs, and
## equations u = 0 at the boundaries are appended, to form a DAE
## system of size N = M^2. Here M = 10.
##
## The system is solved with IDA using the banded linear system
## solver, half-bandwidths equal to M, and default difference-quotient
## Jacobian. For purposes of illustration, IDACalcIC is called to
## compute correct values at the boundary, given incorrect values as
## input initial guesses. The constraints u >= 0 are posed for all
## components.
##

MGRID = 10
NEQ = MGRID*MGRID

dx = 1.0/(MGRID - 1.0)
coeff = 1.0 / (dx * dx)
bval = 0.1
##
## heatres: heat equation system residual function This uses 5-point
## central differencing on the interior points, and includes algebraic
## equations for the boundary values.
##
## So for each interior point, the residual component has the form
##    r_i = u'_i - (central difference)_i
## while for each boundary point, it is r_i = u_i.
##

function heatres(t, u, up, r)
    r[:] = u ## Initialize r to u, to take care of boundary equations.

    ## Loop over interior points; set res = up - (central difference).
    for j = 2:(MGRID-2)
        offset = MGRID * j
        for i = 2:(MGRID-2)
            loc = offset + i
            r[loc] = up[loc] - coeff * (u[loc-1] + u[loc+1] +
                                        u[loc-MGRID] + u[loc+MGRID] -
                                        4.0 * u[loc])
        end
    end

    return Sundials.CV_SUCCESS
end

function initial()
    mm = MGRID
    mm1 = mm - 1

    u  = zeros(NEQ)
    id = ones(NEQ)

    ## initialize u on all grid points
    for j = 1:mm-1
        yfact = dx * j
        offset = mm*j
        for i = 1:mm-1
            xfact = dx * i
            loc = offset + i
            u[loc] = 48.0 * xfact * (1.0 - xfact) * yfact * (1.0 - yfact)
        end
    end

    up = zeros(NEQ)
    r  = zeros(NEQ)

    heatres(0.0, u, up, r)

    ## Copy -res into up to get correct interior initial up values.
    up[:] = -1.0 * r

    ## Finally, set values of u, up, and id at boundary points.
    for j = 1:mm-1
        offset = mm * j
        for i = 1:mm-1
            loc = offset + i;
            if j == 1 || j == mm1 || i == 1 || i == mm1
                u[loc] = bval
                up[loc] = 0
                id[loc] = 0
            end
        end
    end

    constraints = ones(NEQ)
    return (u,up,id,constraints)
end

function idabandsol(f::Function, y0::Vector{Float64}, yp0::Vector{Float64},
                    id::Vector{Float64}, constraints::Vector{Float64},
                    t::Vector{Float64};
                    reltol::Float64=1e-4, abstol::Float64=1e-6)
    neq = length(y0)
    mem = Sundials.IDACreate()
    Sundials.@checkflag Sundials.IDAInit(mem, Sundials.old_cfunction(
                                         Sundials.idasolfun, Cint,
                                         Tuple{Sundials.realtype, Sundials.N_Vector, Sundials.N_Vector, Sundials.N_Vector, Ref{Function}}),
                                         t[1], y0, yp0)
    Sundials.@checkflag Sundials.IDASetId(mem, id)
    Sundials.@checkflag Sundials.IDASetConstraints(mem, constraints)
    Sundials.@checkflag Sundials.IDASetUserData(mem, f)
    Sundials.@checkflag Sundials.IDASStolerances(mem, reltol, abstol)

    A = Sundials.SUNBandMatrix(neq, MGRID, MGRID,2MGRID)
    LS = Sundials.SUNBandLinearSolver(y0,A)
    Sundials.@checkflag Sundials.IDADlsSetLinearSolver(mem, LS, A)

    rtest = zeros(neq)
    Sundials.@checkflag Sundials.IDACalcIC(mem, Sundials.IDA_YA_YDP_INIT, t[2])
    yres = zeros(Float64, length(y0), length(t))
    ypres = zeros(Float64, length(y0), length(t))
    yres[:, 1] = y0
    ypres[:, 1] = yp0
    y = copy(y0)
    yp = copy(yp0)
    tout = [0.0]
    for k in 2:length(t)
        Sundials.@checkflag Sundials.IDASolve(mem, t[k], tout, y, yp, Sundials.IDA_NORMAL)
        yres[:, k] = y
        ypres[:, k] = yp
    end

    Sundials.SUNLinSolFree_Band(LS)
    Sundials.SUNMatDestroy_Band(A)

    return yres, ypres
end

nsteps = 10
tstep = 0.005
t = collect(0.0:tstep:(tstep*nsteps))
u0, up0, id, constraints = initial()

idabandsol(heatres, u0, up0, id, constraints, map(x -> x, t),
                         reltol = 0.0, abstol = 1e-3)


text:
module QPSParser

import Base: show

using SparseArrays

export  canonical,
        objective,
        inequalities,
        equalities,
        bounds,
        variables,
        name,
        parseqps,
        parse_sparseqps

include("mpsdesc.jl")
include("canonicaldesc.jl")
include("parseqps.jl")

end # module


text:
# code reflection entry-points

using InteractiveUtils

#
# code_* replacements
#

# NOTE: these functions replicate parts of the main compiler driver in order to generate
#       more compact code (i.e. without the run-time library) and/or to support generating
#       otherwise invalid code (e.g. with missing symbols).

"""
    code_llvm([io], f, types; optimize=true, cap::VersionNumber, kernel=false,
              optimize=true, raw=false, dump_module=false, strict=false)

Prints the device LLVM IR generated for the method matching the given generic function and
type signature to `io` which defaults to `stdout`.

The following keyword arguments are supported:

- `cap` which device to generate code for
- `kernel`: treat the function as an entry-point kernel
- `optimize`: determines if the code is optimized, which includes kernel-specific
  optimizations if `kernel` is true
- `raw`: return the raw IR including all metadata
- `dump_module`: display the entire module instead of just the function
- `strict`: verify generate code as early as possible

See also: [`@device_code_llvm`](@ref), [`InteractiveUtils.code_llvm`](@ref)
"""
function code_llvm(io::IO, @nospecialize(func), @nospecialize(types);
                   cap::VersionNumber=current_capability(), kernel::Bool=false,
                   optimize::Bool=true, raw::Bool=false,
                   dump_module::Bool=false, strict::Bool=false, kwargs...)
    tt = Base.to_tuple_type(types)
    job = CompilerJob(func, tt, cap, kernel; kwargs...)
    code_llvm(io, job; optimize=optimize,
              raw=raw, dump_module=dump_module, strict=strict)
end
function code_llvm(io::IO, job::CompilerJob; optimize::Bool=true, raw::Bool=false,
                   dump_module::Bool=false, strict::Bool=false)
    ir, entry = codegen(:llvm, job; optimize=optimize, strip=!raw, strict=strict)
    if dump_module
        show(io, ir)
    else
        show(io, entry)
    end
end
code_llvm(@nospecialize(func), @nospecialize(types); kwargs...) =
    code_llvm(stdout, func, types; kwargs...)

"""
    code_ptx([io], f, types; cap::VersionNumber, kernel=false, raw=false, strict=false)

Prints the PTX assembly generated for the method matching the given generic function and
type signature to `io` which defaults to `stdout`.

The following keyword arguments are supported:

- `cap` which device to generate code for
- `kernel`: treat the function as an entry-point kernel
- `raw`: return the raw code including all metadata
- `strict`: verify generate code as early as possible

See also: [`@device_code_ptx`](@ref)
"""
function code_ptx(io::IO, @nospecialize(func), @nospecialize(types);
                  cap::VersionNumber=current_capability(), kernel::Bool=false,
                  raw::Bool=false, strict::Bool=false, kwargs...)
    tt = Base.to_tuple_type(types)
    job = CompilerJob(func, tt, cap, kernel; kwargs...)
    code_ptx(io, job; raw=raw, strict=strict)
end
function code_ptx(io::IO, job::CompilerJob; raw::Bool=false, strict::Bool=false)
    asm, _ = codegen(:ptx, job; strip=!raw, strict=strict)
    print(io, asm)
end
code_ptx(@nospecialize(func), @nospecialize(types); kwargs...) =
    code_ptx(stdout, func, types; kwargs...)

"""
    code_sass([io], f, types, cap::VersionNumber)

Prints the SASS code generated for the method matching the given generic function and type
signature to `io` which defaults to `stdout`.

The following keyword arguments are supported:

- `cap` which device to generate code for
- `kernel`: treat the function as an entry-point kernel

See also: [`@device_code_sass`](@ref)
"""
function code_sass(io::IO, @nospecialize(func), @nospecialize(types);
                   cap::VersionNumber=current_capability(), kernel::Bool=true, kwargs...)
    tt = Base.to_tuple_type(types)
    job = CompilerJob(func, tt, cap, kernel; kwargs...)
    code_sass(io, job)
end
function code_sass(io::IO, job::CompilerJob)
    if !job.kernel
        error("Can only generate SASS code for kernel functions")
    end

    ptx, _ = codegen(:ptx, job)

    fn = tempname()
    gpu = "sm_$(job.cap.major)$(job.cap.minor)"
    # NOTE: this might not match what is being executed, due to the PTX->SASS conversion
    #       by the driver possibly not matching what `ptxas` (part of the toolkit) does.
    # TODO: see how `nvvp` extracts SASS code when doing PC sampling, and copy that.
    Base.run(`$(ptxas[]) --gpu-name $gpu --output-file $fn --input-as-string $ptx`)
    try
        cmd = `$(nvdisasm[]) --print-code --print-line-info $fn`
        for line in readlines(cmd)
            # nvdisasm output is pretty verbose;
            # perform some clean-up and make it look like @code_native
            line = replace(line, r"/\*[0-9a-f]{4}\*/" => "        ") # strip inst addr
            line = replace(line, r"^[ ]{30}" => "   ")               # reduce leading spaces
            line = replace(line, r"[\s+]//##" => ";")                # change line info tag
            line = replace(line, r"^\." => "\n.")                    # break before new BBs
            line = replace(line, r"; File \"(.+?)\", line (\d+)" => s"; Location \1:\2") # rename line info
            println(io, line)
        end
    finally
        rm(fn)
    end
end
code_sass(@nospecialize(func), @nospecialize(types); kwargs...) =
    code_sass(stdout, func, types; kwargs...)

#
# @device_code_* functions
#

export @device_code_lowered, @device_code_typed, @device_code_warntype,
       @device_code_llvm, @device_code_ptx, @device_code_sass,
       @device_code

function emit_hooked_compilation(inner_hook, ex...)
    user_code = ex[end]
    user_kwargs = ex[1:end-1]
    quote
        # wipe the compile cache to force recompilation
        empty!(CUDAnative.compilecache)

        local kernels = 0
        function outer_hook(job)
            kernels += 1
            $inner_hook(job; $(map(esc, user_kwargs)...))
        end

        if CUDAnative.compile_hook[] != nothing
            error("Chaining multiple @device_code calls is unsupported")
        end
        try
            CUDAnative.compile_hook[] = outer_hook
            $(esc(user_code))
        finally
            CUDAnative.compile_hook[] = nothing
        end

        if kernels == 0
            error("no kernels executed while evaluating the given expression")
        end

        nothing
    end
end

# NOTE: these hooks take both a `f` and an inner `f`, because of how `@cuda`/`_cuda` work:
#       kernels are automatically wrapper in a function returning nothing, for usability.
#
#       Julia-level reflection (lowered/typed/warntype) skips these wrapper, because we
#       can't do call-site inlining and the kernel wrapper would hide any meaningful code.
#
#       at the LLVM level, we inline everything so there's no need to hide the wrapper.

"""
    @device_code_lowered ex

Evaluates the expression `ex` and returns the result of
[`InteractiveUtils.code_lowered`](@ref) for every compiled CUDA kernel.

See also: [`InteractiveUtils.@code_lowered`](@ref)
"""
macro device_code_lowered(ex...)
    quote
        buf = Any[]
        function hook(job::CompilerJob)
            append!(buf, code_lowered(job.f, job.tt))
        end
        $(emit_hooked_compilation(:hook, ex...))
        buf
    end
end

"""
    @device_code_typed ex

Evaluates the expression `ex` and returns the result of
[`InteractiveUtils.code_typed`](@ref) for every compiled CUDA kernel.

See also: [`InteractiveUtils.@code_typed`](@ref)
"""
macro device_code_typed(ex...)
    quote
        buf = Any[]
        function hook(job::CompilerJob)
            if VERSION >= v"1.1.0"
                append!(buf, code_typed(job.f, job.tt, debuginfo=:source))
            else
                append!(buf, code_typed(job.f, job.tt))
            end
        end
        $(emit_hooked_compilation(:hook, ex...))
        buf
    end
end

"""
    @device_code_warntype [io::IO=stdout] ex

Evaluates the expression `ex` and prints the result of
[`InteractiveUtils.code_warntype`](@ref) to `io` for every compiled CUDA kernel.

See also: [`InteractiveUtils.@code_warntype`](@ref)
"""
macro device_code_warntype(ex...)
    function hook(job::CompilerJob; io::IO=stdout, kwargs...)
        code_warntype(io, job.f, job.tt; kwargs...)
    end
    emit_hooked_compilation(hook, ex...)
end

"""
    @device_code_llvm [io::IO=stdout, ...] ex

Evaluates the expression `ex` and prints the result of [`InteractiveUtils.code_llvm`](@ref)
to `io` for every compiled CUDA kernel. For other supported keywords, see
[`CUDAnative.code_llvm`](@ref).

See also: [`InteractiveUtils.@code_llvm`](@ref)
"""
macro device_code_llvm(ex...)
    hook(job::CompilerJob; io::IO=stdout, kwargs...) = code_llvm(io, job; kwargs...)
    emit_hooked_compilation(hook, ex...)
end

"""
    @device_code_ptx [io::IO=stdout, ...] ex

Evaluates the expression `ex` and prints the result of [`CUDAnative.code_ptx`](@ref) to `io`
for every compiled CUDA kernel. For other supported keywords, see
[`CUDAnative.code_ptx`](@ref).
"""
macro device_code_ptx(ex...)
    hook(job::CompilerJob; io::IO=stdout, kwargs...) = code_ptx(io, job; kwargs...)
    emit_hooked_compilation(hook, ex...)
end

"""
    @device_code_sass [io::IO=stdout, ...] ex

Evaluates the expression `ex` and prints the result of [`CUDAnative.code_sass`](@ref) to
`io` for every compiled CUDA kernel. For other supported keywords, see
[`CUDAnative.code_sass`](@ref).
"""
macro device_code_sass(ex...)
    hook(job::CompilerJob; io::IO=stdout, kwargs...) = code_sass(io, job; kwargs...)
    emit_hooked_compilation(hook, ex...)
end

"""
    @device_code dir::AbstractString=... [...] ex

Evaluates the expression `ex` and dumps all intermediate forms of code to the directory
`dir`.
"""
macro device_code(ex...)
    only(xs) = (@assert length(xs) == 1; first(xs))
    localUnique = 1
    function hook(job::CompilerJob; dir::AbstractString)
        name = something(job.name, nameof(job.f))
        fn = "$(name)_$(localUnique)"
        mkpath(dir)

        open(joinpath(dir, "$fn.lowered.jl"), "w") do io
            code = only(code_lowered(job.f, job.tt))
            println(io, code)
        end

        open(joinpath(dir, "$fn.typed.jl"), "w") do io
            if VERSION >= v"1.1.0"
                code = only(code_typed(job.f, job.tt, debuginfo=:source))
            else
                code = only(code_typed(job.f, job.tt))
            end
            println(io, code)
        end

        open(joinpath(dir, "$fn.unopt.ll"), "w") do io
            code_llvm(io, job; dump_module=true, raw=true, optimize=false)
        end

        open(joinpath(dir, "$fn.opt.ll"), "w") do io
            code_llvm(io, job; dump_module=true, raw=true)
        end

        open(joinpath(dir, "$fn.ptx"), "w") do io
            code_ptx(io, job)
        end

        open(joinpath(dir, "$fn.sass"), "w") do io
            code_sass(io, job)
        end

        localUnique += 1
    end
    emit_hooked_compilation(hook, ex...)
end


text:
#=
  From
  "Probabilistic logic programming and its applications"
  Luc De Raedt, Leuven
  https://www.youtube.com/watch?v=3lnVBqxjC88
  @ 3:44

  """
  Mike has a bag of marbles with 4 white, 8 blue, and
  6 red marbles. He pulls out one marble from the bag
  and it is red. What is the probability that the
  second marble he pulls out of the bag is white?

  The answer is 0.234931.
  """

  Cf ~/cplint/bag_of_marbles.pl
     ~/blog/bag_of_marbles.blog
     ~/psi/bag_of_marbles.psi
     ~/webppl/bag_of_marbles.wppl

  Also see bag_of_marbles2.jl for an alternative approach (using recursion).

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

#=
Distributions of variable draw1 (num:0)
2.00000 =>   46948  (0.469480)
3.00000 =>   29289  (0.292890)
1.00000 =>   23763  (0.237630)

White is 1 so the probability is 0.237630

=#

@model function bag_of_marbles()

    white = 1
    blue  = 2
    red   = 3
    colors = [white,blue,red]

    start = [4,8,6];

    draw0 ~ Categorical(simplex([start[white],start[blue],start[red]]))

    if draw0 == white
        draw1 ~ Categorical(simplex([start[white]-1,start[blue],start[red]]))
    elseif draw0==blue
        draw1 ~ Categorical(simplex([start[white],start[blue]-1,start[red]]))
    else
        # red
        draw1 ~ Categorical([start[white],start[blue],start[red]-1]|>simplex);
    end

    true ~ Dirac(draw0 == red)

    return draw1==white

end

model = bag_of_marbles()
num_chains = 4

# HH has problem with this!
# chains = sample(model, MH(), MCMCThreads(), 100_000, num_chains)
# chains = sample(model, MH(), 100_000)

# chains = sample(model, PG(15), MCMCThreads(), 10_000, num_chains)
chains = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chains)

# Note: IS don't generate chains the same way as MH, PG, and SMC!
# chains = sample(model, IS(), MCMCThreads(), 1000, num_chains)

display(chains)

show_var_dist_pct(chains,:draw0)
println("\ndraw1: The probability of drawing white is the probability for 1.0")
show_var_dist_pct(chains,:draw1)

# println("\ndraw1==white:")
# genq = generated_quantities(model,chains)
# show_var_dist_pct(genq,20)
# println("quantile:")
# println(quantile(vec(genq), [0.0, 0.025, 0.25, 0.5, 0.75, 0.975, 1.0]))

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
include("multigrid/utils.jl")
include("multigrid/poisson.jl")
include("multigrid/transfer.jl")
include("multigrid/smoothers.jl")
include("multigrid/multigrid.jl")
include("mgpcg.jl")

function solve_poisson_mgpcg(
    d2u::AbstractArray{<:AbstractFloat, N},
    mask::AbstractArray{Bool, 3},
    dx::NTuple{3, Real};
    kwargs...
) where {N}
    N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d"))
    return solve_poisson_mgpcg!(tzero(d2u), d2u, mask, dx; kwargs...)
end

function solve_poisson_mgpcg!(
    u::AbstractArray{<:AbstractFloat, N},
    d2u::AbstractArray{T, N},
    mask::AbstractArray{Bool, 3},
    dx::NTuple{3, Real};
    kwargs...
) where {T, N}
    N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d"))
    size(u) == size(d2u) || throw(DimensionMismatch())
    size(mask) == size(u)[1:3] || throw(DimensionMismatch())

    # crop to avoid unnecessary work. u[!mask] = 0
    sz = size(mask)
    Rc = crop_indices(mask)

    # pad for finite diffy stencil
    I1, I2 = first(Rc), last(Rc)
    Rc = CartesianIndices(ntuple(ndims(Rc)) do n
        max(1, I1[n]-1) : min(sz[n], I2[n]+1)
    end)

    # TODO: performance test proper cutoff
    szc = size(Rc)
    _crop = maximum(sz .- szc) > 10

    if !_crop
        mc = mask
        d2uc = d2u
        uc = u
    else
        mc = similar(mask, szc)
        mc = _tcopyto!(mc, @view(mask[Rc]))

        if N == 3
            d2uc = @view(d2u[Rc])
            uc = similar(u, szc)
            uc = _tcopyto!(uc, @view(u[Rc]))
        else
            d2uc = @view(d2u[Rc,:])
            uc = similar(u, (szc..., size(u, 4)))
            uc = _tcopyto!(uc, @view(u[Rc,:]))
        end
    end

    # solve
    A = Poisson2(T, mc, dx)
    uc = mgpcg!(uc, A, d2uc; kwargs...)

    if _crop
        if N == 3
            _tcopyto!(@view(u[Rc]), uc)
        else
            _tcopyto!(@view(u[Rc,:]), uc)
        end
    end

    return u
end

function solve_poisson_dct(
    d2u::AbstractArray{<:AbstractFloat, N},
    dx::NTuple{3, Real}
) where {N}
    N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d"))
    return solve_poisson_dct!(tcopy(d2u), dx)
end

function solve_poisson_dct!(
    u::AbstractArray{<:AbstractFloat, N},
    d2u::AbstractArray{<:AbstractFloat, N},
    dx::NTuple{3, Real}
) where {N}
    N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d"))
    size(u) == size(d2u) || throw(DimensionMismatch())
    return solve_poisson_dct!(_tcopyto!(u, d2u), dx)
end

function solve_poisson_dct!(
    d2u::AbstractArray{<:AbstractFloat, N},
    dx::NTuple{3, Real}
) where {N}
    N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d"))

    nx, ny, nz = size(d2u)[1:3]
    idx2 = inv(dx[1].*dx[1])
    idy2 = inv(dx[2].*dx[2])
    idz2 = inv(dx[3].*dx[3])

    u = d2u

    FFTW.set_num_threads(num_cores())
    P = plan_dct!(u, 1:3)
    iP = inv(P)

    u = P*u

    X = [2*(cospi(i)-1)*idx2 for i in range(0, step=1/nx, length=nx)]
    Y = [2*(cospi(j)-1)*idy2 for j in range(0, step=1/ny, length=ny)]
    Z = [2*(cospi(k)-1)*idz2 for k in range(0, step=1/nz, length=nz)]

    @inbounds for t in axes(u, 4)
        d2ût = @view(u[:,:,:,t])
        @batch for k in 1:nz
            for j in 1:ny
                for i in 1:nx
                    d2ût[i,j,k] *= inv(X[i] + Y[j] + Z[k])
                end
            end
        end
        # X = Y = Z = 0. Solution is not unique. Set mean to 0
        d2ût[1,1,1] = 0
    end

    # inverse dct
    return iP*u
end

function solve_poisson_fft(
    d2u::AbstractArray{<:AbstractFloat, N},
    dx::NTuple{3, Real}
) where {N}
    N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d"))
    return solve_poisson_fft!(similar(d2u), d2u, dx)
end

function solve_poisson_fft!(
    u::AbstractArray{<:AbstractFloat, N},
    d2u::AbstractArray{<:AbstractFloat, N},
    dx::NTuple{3, Real}
) where {N}
    N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d"))
    size(u) == size(d2u) || throw(DimensionMismatch())

    nx, ny, nz = size(d2u)[1:3]
    idx2 = inv(dx[1].*dx[1])
    idy2 = inv(dx[2].*dx[2])
    idz2 = inv(dx[3].*dx[3])

    _rfft = iseven(nx)
    FFTW.set_num_threads(num_cores())

    # FFTW's rfft is extremely slow with some odd lengths in the first dim
    if _rfft
        nx = nx>>1 + 1
        P = plan_rfft(d2u, 1:3)
        iP = inv(P)
        d2û = P*d2u
    else
        d2û = _tcopyto!(similar(d2u, complex(eltype(d2u))), d2u)
        P = plan_fft!(d2û, 1:3)
        iP = inv(P)
        d2û = P*d2û
    end

    X = [2*(cospi(i)-1)*idx2 for i in fftfreq(size(u, 1), 2)]
    Y = [2*(cospi(j)-1)*idy2 for j in fftfreq(size(u, 2), 2)]
    Z = [2*(cospi(k)-1)*idz2 for k in fftfreq(size(u, 3), 2)]

    @inbounds for t in axes(d2û, 4)
        d2ût = @view(d2û[:,:,:,t])
        @batch for k in 1:nz
            for j in 1:ny
                for i in 1:nx
                    d2ût[i,j,k] *= inv(X[i] + Y[j] + Z[k])
                end
            end
        end
        # X = Y = Z = 0. Solution is not unique. Set mean to 0
        d2ût[1,1,1] = 0
    end

    # inverse fft
    if _rfft
        u = mul!(u, iP, d2û)
    else
        d2û = iP*d2û
        u = _tcopyto!(real, u, d2û)
    end

    return u
end


text:
using CryptoTools
using Test
using Random

@testset "ChaCha20" begin
    plaintext = hex2bytes("00112233445566778899aabbccddeeff")
    key = hex2bytes("00"^32)
    nonce = hex2bytes("00"^12)
    cipher = ChaCha20(key, nonce)
    # @test encrypt(cipher, plaintext) == hex2bytes("04d5394f430044ccd4435310f690935a")
    @test decrypt(cipher, encrypt(cipher, plaintext)) == plaintext

    # https://mkyong.com/java/java-11-chacha20-stream-cipher-examples/
    plaintext = hex2bytes("4a617661202620436861436861323020656e6372797074696f6e206578616d706c652e")
    # key = hex2bytes("ee416df8b5154a4ac48f3930fcfa53ef7f677c8fd7cd093f1328eedfd831db1a")
    # nonce = hex2bytes("9806308f4d1732d2d39beaba")
    # cipher = ChaCha20(key, nonce)
    # ciphertext = hex2bytes("2149db2c32bf82f9e8dc0a709d8c15d5a22eb79d5f692e04f070d46cc7e264631f85e0")
    # @test encrypt(cipher, plaintext) == ciphertext

    # plaintext = hex2bytes("4a617661202620436861436861323020656e6372797074696f6e206578616d706c652e")
    # key = hex2bytes("f95fd5b41783595e41f4cbcd8dc26a782599184e97ccd768ac531aae729781d3")
    # nonce = hex2bytes("84133f2261ef44796e3669dc")
    # cipher = ChaCha20(key, nonce)
    # ciphertext = hex2bytes("7738807b409f3349dbefbeae988482e0e5959c35ee8f8ee8987357db459e10d7fb8c7e")
    # @test encrypt(cipher, plaintext) == ciphertext

    # plaintext = hex2bytes("000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000")
    # key = hex2bytes("0000000000000000000000000000000000000000000000000000000000000000")
    # nonce = hex2bytes("000000000000000000000000")
    # cipher = ChaCha20(key, nonce)
    # ciphertext = hex2bytes("76b8e0ada0f13d90405d6ae55386bd28bdd219b8a08ded1aa836efcc8b770dc7da41597c5157488d7724e03fb8d84a376a43b8f41518a11cc387b669")
    # @test encrypt(cipher, plaintext) == ciphertext

    # https://asecuritysite.com/encryption/salsa20
    # plaintext = base64decode("XKxPmA/tw9Px+ZtL40csmzDVZSPmMtFRI37JMJBIvak=")
    # key = base64decode("WTeEmfKonMWESy6xHvt7Eqi0WJQkUmlbNC19EXI8M3k=")
    # nonce = base64decode("4mgPYADJIPw=")
    # cipher = ChaCha20(key, nonce)
    # ciphertext = base64decode("kER+ZiYPMeLWxUihfPRS1VR+uA==")
    # @test encrypt(cipher, plaintext) == ciphertext
end


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

 function uniaxial_increment!(material, dstrain11,
        dt; dstrain=[dstrain11, -0.3*dstrain11, -0.3*dstrain11, 0.0, 0.0, 0.0],
        max_iter=50, norm_acc=1e-9)
    converged = false
    stress0 = tovoigt(material.variables.stress)
    for i=1:max_iter
        material.ddrivers.time = dt
        material.ddrivers.strain = fromvoigt(SymmetricTensor{2,3,Float64}, dstrain; offdiagscale=2.0)
        integrate_material!(material)
        stress = tovoigt(material.variables_new.stress)
        dstress = stress - stress0
        D = tovoigt(material.variables_new.jacobian)
        dstr = -D[2:end,2:end]\dstress[2:end]
        dstrain[2:end] .+= dstr
        norm(dstr) < norm_acc && (converged = true; break)
    end
    converged || error("No convergence in strain increment")
    return nothing
end


text:
using PaddedMatrices, VectorizationBase

function kcc(K, kc)
    num_k_iter = cld(K, kc)
    kreps_per_iter = ((div(K, num_k_iter)) + 3) & -4
    Krem = K - (num_k_iter - 1) * kreps_per_iter
    kreps_per_iter, Krem
end

using PaddedMatrices, VectorizationBase, ProgressMeter
function pick_K(M_c, r = 1.6, Koff=0)
    L2 = VectorizationBase.CACHE_SIZE[2] ÷ sizeof(Float64)

    Kbase = round(Int, (L2 * r) / M_c)
    Kbase + Koff
end
function pick_N(K_c, r)
    L3 = VectorizationBase.CACHE_SIZE[3] ÷ sizeof(Float64)
    round(Int, (L3 * r) / (K_c * PaddedMatrices.nᵣ)) * PaddedMatrices.nᵣ
end
function jmultpackab!(C, A, B, ::Val{M_c}, ::Val{K_c}, ::Val{N_c}) where {M_c, K_c, N_c}
    M, N = size(C); K = size(B,1)
    GC.@preserve C A B begin
        zc, za, zb = PaddedMatrices.zstridedpointer.((C,A,B))
        @elapsed(PaddedMatrices.jmultpackAB!(
            zc, za, zb, StaticInt{1}(), StaticInt{0}(), Val{M_c}(), Val{K_c}(), Val{N_c}(), M, K, N,
            PaddedMatrices.CloseOpen(0, VectorizationBase.NUM_CORES), Val{VectorizationBase.CACHE_COUNT[3]}()
        ))
    end
end

function bench_size(S, Cs, As, Bs, ::Val{M_c}, ::Val{K_c}, ::Val{N_c}) where {M_c, K_c, N_c}
    jmultpackab!(first(Cs), first(As), first(Bs), Val{M_c}(), Val{K_c}(), Val{N_c}()) # compile
    gflop = 0.0
    for sCAB ∈ zip(S,Cs,As,Bs)
        (s,C,A,B) = sCAB::Tuple{Int,Matrix{Float64},Matrix{Float64},Matrix{Float64}}
        # sleep(0.5)
        t = jmultpackab!(C, A, B, Val{M_c}(), Val{K_c}(), Val{N_c}())
        gf = 2e-9*s*s*s / t
        gflop += gf
    end
    gflop / length(S)
end
function search(Mc_range = 72:24:120, Kc_range = 1.5:0.25:1.7, Nc_range = 1.0:0.1:1.4)
    best = Ref((0,0,0,(0.0,0.0),-Inf))
    gflop_array = let S = round.(Int, exp.(range(log(1500), stop = log(10000), length = 100))),
        As = map(s -> rand(s,s), S), Bs = map(s -> rand(s,s), S), Cs = map(similar, As),
        iter_prod = Iterators.product(Mc_range, Kc_range, Nc_range), p = Progress(length(iter_prod)),
        best = best

        map(iter_prod) do P
            M_c, rk, rn = P
            gflops = bench_size(S, Cs, As, Bs, Val(M_c), Val(rk), Val(rn))
            b = best[]
            K_c = pick_K(M_c, rk)
            N_c = pick_N(K_c, rn)
            recent = (M_c, K_c, N_c, (rk, rn), gflops)
            bb = if last(b) > gflops
                b
            else
                best[] = recent
            end
            ProgressMeter.next!(p, showvalues = [(:Last, recent), (:Best, bb)])
            gflops
        end
    end
    gflop_array, best
end

search_range = (72:24:120, 0.52:0.02:0.56, 0.9:0.01:0.92)
gflop_array, best = search(search_range...); getindex.(search_range, Tuple(argmax(gflop_array)))

# julia> gflop_array, best = search(search_range...); getindex.(search_range, Tuple(argmax(gflop_array)))
# Progress: 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| Time: 0:23:06
#   Last:  (120, 557, 5301, (0.51, 0.91), 1468.106008702092)
#   Best:  (72, 819, 3600, (0.45, 0.91), 1563.9777269382319)
# (72, 0.45, 0.91)

# 72 with (0.4-0.44) x 0.9-9.92 looks good
# julia> permutedims(gflop_array, (2,3,1))
# 7×5×2 Array{Float64, 3}:
# [:, :, 1] =   # (96 x 0.40:0.01:0.46 x 0.9:0.005:0.92)
#  1566.74  1512.84  1521.73  1558.3   1559.93
#  1569.36  1563.85  1563.67  1559.56  1563.38
#  1566.34  1551.83  1562.32  1558.95  1561.91
#  1567.51  1561.19  1560.29  1557.78  1564.33
#  1567.02  1560.08  1547.32  1559.42  1565.47
#  1501.42  1550.91  1484.06  1554.05  1553.47
#  1503.18  1545.86  1549.18  1550.41  1548.85

# [:, :, 2] =   # (72 x 0.40:0.01:0.46 x 0.9:0.005:0.92)
#  1531.42  1521.71  1520.14  1520.35  1520.72
#  1539.86  1536.97  1532.17  1532.33  1536.38
#  1543.98  1541.24  1535.58  1540.67  1536.43
#  1550.4   1543.33  1542.23  1540.41  1548.76
#  1552.08  1544.08  1548.16  1547.86  1546.74
#  1548.65  1487.99  1544.13  1542.94  1541.12
#  1476.7   1512.51  1542.46  1537.65  1535.09

using JuMP
model = Model()
@variable(model, M_cf, Int)
@variable(model, K_c, Int)
@variable(model, N_cf, Int)
@constraints(
    model, begin
    1 <= M_cf <= 20 # M_c = M_cf * 24
    1 <= K_c <= 4000
    1 <= N_cf <= 2000 # N_c = N_cf * 9
    end
)
@objective(model, Max, gflops)

using MaBLAS, BenchmarkTools

function matmul_params(::Type{T}, mᵣ, nᵣ, mcratio, kcratio, ncratio) where {T}
    L₁, L₂, L₃ = MaBLAS.VectorizationBase.CACHE_SIZE .÷ sizeof(T)
    kc = round(Int, (kcratio * L₁ / nᵣ))
    mc = round(Int, (mcratio * L₂ / (kc*mᵣ))) * mᵣ
    nc = round(Int, (ncratio * L₃ / (kc*nᵣ))) * nᵣ
    mc, kc, nc
end

function search()
    kernels = [
        (16,12),
        (16,14),
        (24,9),
        (32,6),
        (40,5)
    ]
    Kcratio = 0.4:0.125:0.9
    Mcratio = 0.4:0.125:0.9
    Ncratio = 0.4:0.125:0.9
    matmul_sizes = (6:20) .^ 3
    GFLOPS = fill(-Inf, length(Mcratio), length(Kcratio), length(Ncratio), length(kernels), length(matmul_sizes))
    params = Array{Tuple{Tuple{Int,Int},NTuple{3,Int}}}(undef, length(Mcratio), length(Kcratio), length(Ncratio), length(kernels))
    BenchmarkTools.DEFAULT_PARAMETERS.seconds = timelim = 1e-2

    for (szi,sz) ∈ enumerate(matmul_sizes)
        C = Matrix{Float64}(undef, sz, sz);
        A = rand(sz, sz); B = rand(sz, sz);
        for (rri,(mr,nr)) ∈ enumerate(kernels), (nci,ncf) ∈ enumerate(Ncratio), (kci,kcf) ∈ enumerate(Kcratio), (mci,mcf) ∈ enumerate(Mcratio)
            mc, kc, nc = matmul_params(Float64, mr, nr, mcf, kcf, ncf)

            t1 = @elapsed MaBLAS.mul!(C, A, B; packing=true, cache_params=(cache_m=mc, cache_k=kc, cache_n=nc), kernel_params=(Val(mr), Val(nr)))
            t0 = if isone(szi) || t1 < 1e-3
                t2 = @belapsed MaBLAS.mul!($C, $A, $B; packing=true, cache_params=(cache_m=$mc, cache_k=$kc, cache_n=$nc), kernel_params=(Val($mr), Val($nr)))
                min(t1,t2)
            elseif 5t1 < timelim
                t2 = @elapsed MaBLAS.mul!(C, A, B; packing=true, cache_params=(cache_m=mc, cache_k=kc, cache_n=nc), kernel_params=(Val(mr), Val(nr)))
                min(t1,t2)
            else
                t1
            end
            GFLOPS[mci,kci,nci,rri,szi] = gflops = 2e-9*sz^3/t0
            isone(szi) && (params[mci,kci,nci,rri] = ((mr,nr), (mc,kc,nc)))
            @show sz, (mr,nr), (mc,kc,nc), gflops
        end
    end
    GFLOPS, params
end

#gflops, params = search()


text:
function respond_keyboard_event(E::Ensamble, M::SelectionMode, event)
    if event.key == Keyboard.i
        println("Insert Mode")
        E.mode = InsertMode()
    elseif event.key in (Keyboard.delete, Keyboard.d, Keyboard.x)
        to_be_deleted = Any[]
        for i in E.selected
            if i > 0
                push!(to_be_deleted, E.atoms[i])
            else
                push!(to_be_deleted, E.bonds[-i])
            end
        end

        for s in to_be_deleted
            delete!(E, s)
        end

        E.selected = Int[]
    end
    false
end

function respond_keyboard_event(E::Ensamble, M::InsertMode, event)
    if event.key == Keyboard.escape
        println("Selection Mode")
        E.mode = SelectionMode()
    end
    false
end


text:
@with_kw mutable struct LagrangeConstrainedPOMDP{S, A, O} <: POMDP{S,A,O}
    pomdp::POMDP{S,A,O}
    λ::Float32 = 0f0
    terminate_on_violation::Bool = false
    terminate_next::Bool = false
end

POMDPs.initialstate(mdp::LagrangeConstrainedPOMDP) = initialstate(mdp.pomdp)
POMDPs.initialobs(mdp::LagrangeConstrainedPOMDP, s) = initialobs(mdp.pomdp, s)

POMDPs.actions(mdp::LagrangeConstrainedPOMDP) = actions(mdp.pomdp)
POMDPs.actionindex(mdp::LagrangeConstrainedPOMDP, a) = actionindex(mdp.pomdp)

function POMDPs.isterminal(mdp::LagrangeConstrainedPOMDP, s)
    if mdp.terminate_next
        mdp.terminate_next = false
        return true
    else
        isterminal(mdp.pomdp, s)
    end
end
POMDPs.discount(mdp::LagrangeConstrainedPOMDP) = discount(mdp.pomdp)

POMDPs.observation(mdp::LagrangeConstrainedPOMDP, s) = observation(mdp.pomdp, s)

render(mdp::LagrangeConstrainedPOMDP, s, a = nothing; kwargs...) = render(mdp.pomdp, s, a; kwargs...)
render(mdp::LagrangeConstrainedPOMDP; kwargs...) = render(mdp.pomdp; kwargs...)

function POMDPs.gen(mdp::LagrangeConstrainedPOMDP, s, a, rng::AbstractRNG = Random.GLOBAL_RNG; info=Dict())
    sp, o, r = gen(mdp.pomdp, s, a, rng, info=info)
    if haskey(info, "cost") && info["cost"] >0
        r -= mdp.λ*info["cost"]
        if mdp.terminate_on_violation
            mdp.terminate_next = true
        end
    end
    return (sp=sp, o=o, r=r)
end

max_forks_repo_name: ancorso/POMDPGym


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct LCLFilter <: Filter
        lf::Float64
        rf::Float64
        cf::Float64
        lg::Float64
        rg::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
    end

Parameters of a LCL filter outside the converter, the states are in the grid's reference frame

# Arguments
- `lf::Float64`: Series inductance in p.u. of converter filter, validation range: `(0, nothing)`
- `rf::Float64`: Series resistance in p.u. of converter filter, validation range: `(0, nothing)`
- `cf::Float64`: Shunt capacitance in p.u. of converter filter, validation range: `(0, nothing)`
- `lg::Float64`: Series inductance in p.u. of converter filter to the grid, validation range: `(0, nothing)`
- `rg::Float64`: Series resistance in p.u. of converter filter to the grid, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states of the LCLFilter model are:
	ir_cnv: Real current out of the converter,
	ii_cnv: Imaginary current out of the converter,
	vr_filter: Real voltage at the filter's capacitor,
	vi_filter: Imaginary voltage at the filter's capacitor,
	ir_filter: Real current out of the filter,
	ii_filter: Imaginary current out of the filter
- `n_states::Int`: LCLFilter has 6 states
"""
mutable struct LCLFilter <: Filter
    "Series inductance in p.u. of converter filter"
    lf::Float64
    "Series resistance in p.u. of converter filter"
    rf::Float64
    "Shunt capacitance in p.u. of converter filter"
    cf::Float64
    "Series inductance in p.u. of converter filter to the grid"
    lg::Float64
    "Series resistance in p.u. of converter filter to the grid"
    rg::Float64
    ext::Dict{String, Any}
    "The states of the LCLFilter model are:
	ir_cnv: Real current out of the converter,
	ii_cnv: Imaginary current out of the converter,
	vr_filter: Real voltage at the filter's capacitor,
	vi_filter: Imaginary voltage at the filter's capacitor,
	ir_filter: Real current out of the filter,
	ii_filter: Imaginary current out of the filter"
    states::Vector{Symbol}
    "LCLFilter has 6 states"
    n_states::Int
end

function LCLFilter(lf, rf, cf, lg, rg, ext=Dict{String, Any}(), )
    LCLFilter(lf, rf, cf, lg, rg, ext, [:ir_cnv, :ii_cnv, :vr_filter, :vi_filter, :ir_filter, :ii_filter], 6, )
end

function LCLFilter(; lf, rf, cf, lg, rg, ext=Dict{String, Any}(), states=[:ir_cnv, :ii_cnv, :vr_filter, :vi_filter, :ir_filter, :ii_filter], n_states=6, )
    LCLFilter(lf, rf, cf, lg, rg, ext, states, n_states, )
end

# Constructor for demo purposes; non-functional.
function LCLFilter(::Nothing)
    LCLFilter(;
        lf=0,
        rf=0,
        cf=0,
        lg=0,
        rg=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`LCLFilter`](@ref) `lf`."""
get_lf(value::LCLFilter) = value.lf
"""Get [`LCLFilter`](@ref) `rf`."""
get_rf(value::LCLFilter) = value.rf
"""Get [`LCLFilter`](@ref) `cf`."""
get_cf(value::LCLFilter) = value.cf
"""Get [`LCLFilter`](@ref) `lg`."""
get_lg(value::LCLFilter) = value.lg
"""Get [`LCLFilter`](@ref) `rg`."""
get_rg(value::LCLFilter) = value.rg
"""Get [`LCLFilter`](@ref) `ext`."""
get_ext(value::LCLFilter) = value.ext
"""Get [`LCLFilter`](@ref) `states`."""
get_states(value::LCLFilter) = value.states
"""Get [`LCLFilter`](@ref) `n_states`."""
get_n_states(value::LCLFilter) = value.n_states

"""Set [`LCLFilter`](@ref) `lf`."""
set_lf!(value::LCLFilter, val) = value.lf = val
"""Set [`LCLFilter`](@ref) `rf`."""
set_rf!(value::LCLFilter, val) = value.rf = val
"""Set [`LCLFilter`](@ref) `cf`."""
set_cf!(value::LCLFilter, val) = value.cf = val
"""Set [`LCLFilter`](@ref) `lg`."""
set_lg!(value::LCLFilter, val) = value.lg = val
"""Set [`LCLFilter`](@ref) `rg`."""
set_rg!(value::LCLFilter, val) = value.rg = val
"""Set [`LCLFilter`](@ref) `ext`."""
set_ext!(value::LCLFilter, val) = value.ext = val


text:
function se2 end
function so2 end
function rot2 end
function trot2 end
function transl2 end

function se3 end
function so3 end
function rotx end
function roty end
function rotz end
function trotx end
function troty end
function trotz end
function transl3 end

function transl end


text:
using GeneralizedDynamicsFromData
using ModelingToolkit
using DataDrivenDiffEq
using LinearAlgebra
using Plots

#### Setup parameters of the ODE, UDE system and training
parameters_steadystate = Float32[0.1, 0.15]
parameters_oscillation = Float32[0.1, 0.6]

parameters = parameters_one_stable

u₀ = Float64[1.0, 1.0]
tspan =(0.0f0, 15.0f0)
s = 0.1

loss_weights = [0.2, 0.8]
max_iters = 3000

#### Setup functions for neural network and differential equations
UA = regression_model(
    2, # input_d
    1, # output_d
    16, # hidden_d
    1, # hidden_layers
    tanh, # nonlinearity
    Flux.glorot_normal # initialization
)

real_de!, univ_de = selkov(parameters, UA)

#### solve the real ODE and create a target time series
real_prob = ODEProblem(real_de!, u₀, tspan)
real_solution = solve(real_prob, Tsit5(), saveat = s)
y = Array(real_solution)

#### set up the UDE problem and training of the neural network
θ₀ = initial_params(UA)
prob_nn = ODEProblem(univ_de, u₀, tspan, θ₀)

first_guess = concrete_solve(prob_nn, Tsit5(), u₀, θ₀, saveat=real_solution.t)
predict(θ) = Array(concrete_solve(prob_nn, Tsit5(), u₀, θ, saveat=real_solution.t))
loss(θ) = polar_loss(θ, y, loss_weights, predict)

#### train the UDE with ADAM
callback = CallbackLog(T=Float64)
opt = Flux.Optimiser(WeightDecay(1e-4), ExpDecay(1e-1,0.5,100,1e-4), ADAM())
res = DiffEqFlux.sciml_train(loss, θ₀, opt, cb=callback, maxiters=max_iters)

trained_parameters = res.minimizer
ŷ = predict(trained_parameters)
ŷ_derivs = UA(ŷ, trained_parameters)

@variables u[1:2]
b = polynomial_basis(u,5)
basis = Basis(b, u)

opt = SR3(Float64(1e-2), Float64(1e-1))
λ = exp10.(-7:0.1:5)
g(x) = x[1] > 1 ? Inf : norm(x, 2)

Ψ = SINDy(ŷ,
          ŷ_derivs,
          basis,
          λ,
          opt,
          g = g,
          maxiter = 50000,
          normalize = true,
          denoise = true,
          convergence_error = Float64(1e-10))
println(Ψ)
print_equations(Ψ)

max_stars_repo_name: pnieters/GeneralizedDynamicsFromData

max_issues_repo_name: pnieters/GeneralizedDynamicsFromData

max_forks_repo_name: pnieters/GeneralizedDynamicsFromData


text:
# naiveunify.jl

## original algorithm

"""
find disagreement pair
"""
function disagreement(t1::Any, t2::Any)
 if t1==t2; return () end
 if !isa(t1,Expr) || !isa(t2,Expr)
   return (t1,t2)
 end
 if t1.args[1] != t2.args[1]; return (t1,t2) end

 for ix = 2:length(t1.args)
  a1 = t1.args[ix]
  a2 = t2.args[ix]
  Δ=disagreement(a1, a2)
  if Δ != (); return Δ end
 end
 return ()
end

### unify0(vars, t1, t2,
function unify0(vars, t1, t2, sigma)

end

### t1==t2 is resolved in unify
## t1!=t2 is assumed

function direction(vars::Vlist, t1::Any, t2::Any)
  if isvar(t1, vars)
    isa(t2, Expr) && loopcheck(t1, t2)
    return (t1, t2)
  end
  if isvar(t2, vars);
    isa(t1, Expr) && loopcheck(t2, t1)
    return (t2, t1)
  end
#@show typeof(t1), typeof(t2)
  throw(ICMP(t1,t2,:directionaa))
end

function unify1(vars::Vlist, t1::Expr, t2::Expr, σ::Tlist)::Tlist
  if t1.args[1] != t2.args[1]; throw(ICMP(t1,t2,:unify1ee)) end
  while true
    Δ=disagreement(t1, t2)
# one of Δ should be a Symbol(var)
    if Δ==(); break end
    Δa=direction(vars, Δ[1],Δ[2])
    σa = putasubst(vars,Δa[1],Δa[2],σ)
    t1 = apply(vars, t1, σa)
    t2 = apply(vars, t2, σa)
    σ = apply2(vars, σ, σa)
  end
  return σ
end

"""
unify 2 terms
"""
function unify(vars::Vlist, t1::Any, t2::Any)::Tlist
  σ = deepcopy(vars)
  if t1==t2; return σ end
  unify1(vars, t1, t2, σ)
end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
function saliency_map(network, input; coeffs = [-0.74, -0.44])
    """ Computes the gradient of coeffs' * y with respect to the input x
        Args:
            - network: flux network
            - input: input to get saliency map for
            ------------------
            - coeffs: coefficients to get output (linear combination of network outputs)
    """
    grad_function(x) = coeffs' * network(x)
    _, b = Flux.pullback(() -> grad_function(input), Flux.params(input))
    return b(1)[input]
end

function integrated_gradients(network, input, baseline; coeffs = [-0.74, -0.44], res = 20)
    """
    Args:
        - network: flux network
        - input: input to get integrated gradients for
        - baseline: IG baseline (should be same size as input)
        ----------------
        - coeffs: coefficients to get output (linear combination of network outputs)
        - number of samples used to approximate integral
    """
    grad_function(x) = coeffs' * network(x)
    αs = collect(range(0.0, stop = 1.0, length = res))
    attrs = zeros(length(input))
    for i = 1:res
        interp_inp = αs[i] .* input .+ (1 .- αs[i]) .* baseline
        attrs += saliency_map(network, interp_inp)
    end
    return attrs ./ res
end

function get_mask(attrs, x; threshold = 0.9, masked_opacity = 0.2)
    """ Returns a masked image that highlights the pixels with the highest
        absolute attributions
        Args:
            - attrs: vector of attributions (same size as x)
            - x: input associated with the attrs (vector)
            -------------
            - threshold: highlight pixels above this percentile
            - masked_opacity: opacity for non-highlighted pixels
    """
    perc = quantile(attrs, threshold)
    mask = ones(length(attrs))
    mask[attrs .< perc] .= masked_opacity
    masked_input = mask .* x
    return masked_input
end


text:
import Base.convert
using Dates

convert(::Type{Int}, v::SubString{String}) = parse(Int, v)
convert(::Type{Float64}, v::SubString{String}) = parse(Float64, v)
convert(::Type{Dates.Date}, v::String) = Date(v)

max_stars_repo_name: ky-s/DiaryApp

max_issues_repo_name: ky-s/DiaryApp

max_forks_repo_name: ky-s/DiaryApp


text:
@testset "Remove Vulkan prefixes" begin
    @test remove_prefix(s2) == s2_novk
    @test remove_prefix(c1) == c1_novk
    @test remove_prefix(c2) == c2_novk
end


text:
# Michelle Medina
# June 18, 2020
# Code that does linear interpolation of a function

export interpolate

function interpolate(p1::NTuple{3,T}, p2::NTuple{3,T}, p3::NTuple{3,T}, p4::NTuple{3,T}, p5::NTuple{2,T})::T where {T}
	x1, y1, F11 = p1
	x1, y2, F12 = p2
	x2, y1, F21 = p3
	x2, y2, F22 = p4
	xeval, yeval = p5
	dx = x2 - x1
	dy = y2 - y1
	Fxy = ((y2 - yeval)/dy)*(((x2 - xeval)/dx)*F11 + ((xeval - x1)/dx)*F21) + ((yeval - y1)/dy)*(((x2 - xeval)/dx)*F12 + ((xeval - x1)/dx)*F22)
	return Fxy
end

function interpolate(p1::NTuple{2,T}, p2::NTuple{2,T}, xeval::T)::T where {T}
	x1, y1 = p1
	x2, y2 = p2
	y = y1 + (xeval - x1)*((y2 - y1)/(x2 - x1))
    return y
end


text:
@defcomp TotalAdaptationCosts begin
    region = Index()
    # Total Adaptation Costs
    pop_population = Parameter(index=[time, region], unit="million person")
    ac_adaptationcosts_economic = Parameter(index=[time, region], unit="\$million")
    ac_adaptationcosts_noneconomic = Parameter(index=[time, region], unit="\$million")
    ac_adaptationcosts_sealevelrise = Parameter(index=[time, region], unit="\$million")

    act_adaptationcosts_total = Variable(index=[time, region], unit="\$million")
    act_percap_adaptationcosts = Variable(index=[time, region], unit="\$/person")

    function run_timestep(p, v, d, t)

        for r in d.region
            v.act_adaptationcosts_total[t,r] = p.ac_adaptationcosts_economic[t,r] + p.ac_adaptationcosts_sealevelrise[t,r] + p.ac_adaptationcosts_noneconomic[t,r]
            v.act_percap_adaptationcosts[t,r] = v.act_adaptationcosts_total[t,r] / p.pop_population[t,r]
        end
    end
end


text:
# # SVM with L^1 regularization
## Generate data for SVM classifier with L1 regularization.
using Random
Random.seed!(3);
n = 20;
m = 1000;
TEST = m;
DENSITY = 0.2;
beta_true = randn(n, 1);
idxs = randperm(n)[1:round(Int, (1 - DENSITY) * n)];
beta_true[idxs] .= 0
offset = 0;
sigma = 45;
X = 5 * randn(m, n);
Y = sign.(X * beta_true .+ offset .+ sigma * randn(m, 1));
X_test = 5 * randn(TEST, n);

#-

## Form SVM with L1 regularization problem.
using Convex, SCS, ECOS

beta = Variable(n);
v = Variable();
loss = sum(pos(1 - Y .* (X * beta - v)));
reg = norm(beta, 1);

## Compute a trade-off curve and record train and test error.
TRIALS = 100
train_error = zeros(TRIALS);
test_error = zeros(TRIALS);
lambda_vals = exp10.(range(-2, stop = 0, length = TRIALS);)
beta_vals = zeros(length(beta), TRIALS);
for i in 1:TRIALS
    lambda = lambda_vals[i]
    problem = minimize(loss / m + lambda * reg)
    solve!(problem, SCS.Optimizer; silent_solver = true)
    train_error[i] =
        sum(
            float(
                sign.(X * beta_true .+ offset) .!=
                sign.(evaluate(X * beta - v)),
            ),
        ) / m
    test_error[i] =
        sum(
            float(
                sign.(X_test * beta_true .+ offset) .!=
                sign.(evaluate(X_test * beta - v)),
            ),
        ) / TEST
    beta_vals[:, i] = evaluate(beta)
end

#-

# Plot the train and test error over the trade-off curve.
using Plots
plot(lambda_vals, train_error, label = "Train error");
plot!(lambda_vals, test_error, label = "Test error");
plot!(xscale = :log, yscale = :log, ylabel = "errors", xlabel = "lambda")

#-

# Plot the regularization path for beta.

plot()
for i in 1:n
    plot!(lambda_vals, vec(beta_vals[i, :]), label = "beta$i")
end
plot!(xscale = :log, ylabel = "betas", xlabel = "lambda")