text:
include("graphs.jl")
using DelimitedFiles

function incidencematrix(graph::Graph)
	n = numnodes(graph)
	m = numedges(graph)
	E = zeros(Float64, n, m)
	for i=1:n
		v = graph.nodelist[i]
		for j in v.outedges
			E[i,j] = 1.0
		end
		for j in v.inedges
			E[i,j] = -1.0
		end
	end
	return E
end

function incidencematrixtonetwork(E::Matrix{Float64})
	n = size(E,1)
	m = size(E,2)
	net = Graph(Node[], Edge[], 1, 1)
	for i=1:n
		addnode!(net)
	end
	for j=1:m
		origin = 0
		target = 0
		for i=1:n
			if E[i,j] < -0.1
				origin = i
			end
			if E[i,j] > 0.1
				target = i
			end
		end
		addedge!(net, origin, target)
	end
	return net
end

function getdegrees(E::Matrix{Float64})
	n = size(E,1)
	m = size(E,2)
	degrees = zeros(Float64, n)
	for i=1:n
		for j=1:m
			degrees[i] += abs(E[i,j])
		end
	end
	return degrees
end

function adjacencymatrix(E::Matrix{Float64})
	return Matrix(Diagonal(getdegrees(E))) - laplacianmatrix(E)
end

function laplacianmatrix(E::Matrix{Float64})
	return E*transpose(E)
end

function weightedlaplacianmatrix(E::Matrix{Float64}, w::Vector{Float64})
	return E*Matrix(Diagonal(w))*transpose(E)
end

function wattsstrogatz!(graph::Graph, E::Matrix{Float64}, q::Float64)
	n = size(E,1)
	m = size(E,2)
	A = adjacencymatrix(E)
	for j=1:m
		z = rand()
		if z < q
			e = graph.edgelist[j]
			s = e.origin
			t = e.target
			while true
				possible_t = rand(1:n)
				if (possible_t ≠ s) && (A[s, possible_t] < 0.1)
					A[s, possible_t] = 1.0
					A[possible_t, s] = 1.0
					A[s, t] = 0.0
					A[t, s] = 0.0
					E[t, j] = 0.0
					E[possible_t, j] = -1.0
					break
				end
			end
		end
	end
end

function scalefreenetwork(n::Int64, m::Int64, m0::Int64)
	net = Graph(Node[], Edge[], 1, 1)
	for i=1:4
		addnode!(net)
	end
	addedge!(net, 1, 2)
	addedge!(net, 1, 3)
	addedge!(net, 3, 4)
	addedge!(net, 1, 4)

	while numnodes(net) < n
		addnode!(net)
		v = net.nodelist[end]
		i = v.nodeindex
		while length(v.outedges) < m0
			j = rand(1:i-1)
			w = net.nodelist[j]
			k = length(w.outedges) + length(w.inedges)
			if rand() < Float64(k)/(2.0*numedges(net))
				addedge!(net, i, j)
			end
		end
	end

	E = incidencematrix(net)
	return net, E
end

function smallworldnetwork(n::Int64, k::Int64, q::Float64)
	net = Graph(Node[], Edge[], 1, 1)
	for i=1:n
		addnode!(net)
	end
	if k ≥ n-1
		println("complete")
		for i=1:(n-1)
			for j=(i+1):n
				addedge!(net, i, j)
			end
		end
	else
		for i=1:n
			for j=1:div(k,2)
				addedge!(net, i, mod1((i+j),n))
			end
		end
	end
	E = incidencematrix(net)
	wattsstrogatz!(net, E, q)
	net = incidencematrixtonetwork(E)
	return net, E
end

function austriannetwork()
	E = readdlm("austria_incidence.txt")
	net = incidencematrixtonetwork(E)
	return net, E
end

function testnetwork()
	net = Graph(Node[], Edge[], 1, 1)
	for i=1:6
		addnode!(net)
	end
	addedge!(net, 1, 2)
	addedge!(net, 1, 3)
	addedge!(net, 2, 3)
	addedge!(net, 3, 4)
	addedge!(net, 4, 5)
	addedge!(net, 4, 6)
	addedge!(net, 5, 6)
	E = incidencematrix(net)
	X = zeros(Float64, 7)
	X .+= 1.0
	P = zeros(Float64, 6)
	P[1] = 1.0
	P[2] = -0.5
	P[6] = -0.5
	return net, E, X, P
end

function sourcesinklocs(ns::Int64, nd::Int64, n::Int64)
	sources = Array{Int64, 1}([])
	sinks = Array{Int64, 1}([])
	for i=1:ns
		while true
			z = rand(1:n)
			if z ∉ sources
				push!(sources, z)
				break
			end
		end
	end
	for i=1:nd
		while true
			z = rand(1:n)
			if (z ∉ sinks) & (z ∉ sources)
				push!(sinks, z)
				break
			end
		end
	end
	return sources, sinks
end

function sourcesinklocsmod(net::Graph, ns::Int64, nd::Int64, n::Int64)
	sources = Array{Int64, 1}([])
	sinks = Array{Int64, 1}([])
	if ns == 1
		while true
			z = rand(1:n)
			if length(net.nodelist[z].outedges) + length(net.nodelist[z].inedges) > 3
				push!(sources, z)
				break
			end
		end
	end
	if nd == 1
		while true
			z = rand(1:n)
			if length(net.nodelist[z].outedges) + length(net.nodelist[z].inedges) > 3
				push!(sinks, z)
				break
			end
		end
	end
	if ns > 1
		for i=1:ns
			while true
				z = rand(1:n)
				if (z ∉ sinks) & (z ∉ sources)
					push!(sources, z)
					break
				end
			end
		end
	end
	if nd > 1
		for i=1:nd
			while true
				z = rand(1:n)
				if (z ∉ sinks) & (z ∉ sources)
					push!(sinks, z)
					break
				end
			end
		end
	end
	return sources, sinks
end

function sourcesinkwellmixed(ns::Int64, nd::Int64, n::Int64)
	sources = Array{Int64, 1}([])
	sinks = Array{Int64, 1}([])
	for i=1:2:n
		push!(sources, i)
	end
	for i=1:n
		if i ∉ sources
			push!(sinks, i)
		end
	end
	return sources, sinks
end

function sourcesinkhalfhalf(ns::Int64, nd::Int64, n::Int64)
	sources = Array{Int64, 1}([])
	sinks = Array{Int64, 1}([])
	for i=1:div(n,2)
		push!(sources, i)
	end
	for i=1:n
		if i ∉ sources
			push!(sinks, i)
		end
	end
	return sources, sinks
end

function sourcessinkclockface(n::Int64, pos::Int64)
	sources = Array{Int64, 1}([])
	sinks = Array{Int64, 1}([])
	push!(sources,1)
	push!(sources, pos)
	for i=2:n
		if i ≠ pos
			push!(sinks, i)
		end
	end
	return sources, sinks
end

function sourcesinkvector(sources::Vector{Int64}, sinks::Vector{Int64}, n::Int64, Ptot::Float64)
	P = zeros(Float64, n)
	sourcenum = length(sources)
	sinknum = length(sinks)
	for i ∈ sources
		P[i] = Ptot/Float64(sourcenum)
	end
	for i ∈ sinks
		P[i] = -Ptot/Float64(sinknum)
	end
	return P
end

function makenoise(sourcenum::Int64, sinknum::Int64)
	d = Normal(0.0, 0.4)
	sourcenoise = rand(d, sourcenum)
	sinknoise = rand(d, sinknum)
	return sourcenoise, sinknoise
end

function makegammasourcesinkvector(sources::Vector{Int64}, sinks::Vector{Int64}, n::Int64, Ptot::Float64)
	P = P = zeros(Float64, n)
	sourcenum = length(sources)
	sinknum = length(sinks)
	d = Gamma(2.0, 0.5)
	runningsum = 0.0
	for i ∈ sources
		P[i] = rand(d)
		runningsum += P[i]
	end
	alpha = Ptot/runningsum
	for i ∈ sources
		P[i] *= alpha
	end
	runningsum = 0.0
	for i ∈ sinks
		P[i] = -rand(d)
		runningsum += abs(P[i])
	end
	alpha = Ptot/runningsum
	for i ∈ sinks
		P[i] *= alpha
	end
	return P
end

function noisysourcesinkvector(sources::Vector{Int64}, sinks::Vector{Int64}, n::Int64, Ptot::Float64)
	P = zeros(Float64, n)
	sourcenum = length(sources)
	sinknum = length(sinks)
	sourcenoise, sinknoise = makenoise(sourcenum, sinknum)
	runningsum = 0.0
	for (i, val) ∈ enumerate(sources)
		P[val] = (Ptot/Float64(sourcenum))*(1.0 + sourcenoise[i])
		runningsum += P[val]
	end
	delta = Ptot - runningsum
	for val ∈ sources
		P[val] += delta/sourcenum
	end
	runningsum = 0.0
	for (i, val) ∈ enumerate(sinks)
		P[val] = -(Ptot/Float64(sinknum))*(1.0 + sinknoise[i])
		runningsum += P[val]
	end
	delta = -Ptot - runningsum
	for val ∈ sinks
		P[val] += delta/sinknum
	end
	return P
end

max_stars_repo_name: o-smith/GridResilience

max_issues_repo_name: o-smith/GridResilience

max_forks_repo_name: o-smith/GridResilience


text:
module DataGetters

    HEAD = 0
    BODY = 1
    DIR = 2
    FOOD = 3
    SCORE = 4
    HAS_FOOD = 5
    ALIVE = 6
    HUNGER = 7
    HIST = 8
    POSITIVE_RESPONSE = 100.0
    NEGATIVE_RESPONSE = -20.0
    NEUTRAL_RESPONSE = 1.0
    BORDER_THRESHOLD = 10.0
    HUNGER_FACTOR = 0.01
    TILE_SIZE = 20
    TERRAIN_DIMS = (30, 30)
    HEIGHT = TERRAIN_DIMS[1] * TILE_SIZE
    WIDTH = TERRAIN_DIMS[2] * TILE_SIZE
    function data_simple(state)
        xs = []
        fx, fy = state[FOOD]
        dir = state[DIR]
        hx, hy = state[HEAD]
        fx, fy = state[FOOD]
        xs = vcat(xs, [hx, hy])
        xs = vcat(xs, [fx, fy])
        goal = NEUTRAL_RESPONSE
        dist = abs(fx - hx) + abs(fy - hy)
        append!(xs, dist)
        goal -= HUNGER_FACTOR * state[HUNGER]
        if !state[HAS_FOOD]
            goal = POSITIVE_RESPONSE
        end

        dist_wall_x = min(abs(hx - WIDTH), hx)
        dist_wall_y = min(abs(hy - HEIGHT), hy)
        dist_wall_food_x = min(abs(fx - WIDTH), fx)
        dist_wall_food_y = min(abs(fy - HEIGHT), fy)
        xs = vcat(xs, [dist_wall_x, dist_wall_y])
        xs = vcat(xs, [dist_wall_food_x, dist_wall_food_y])

        append!(xs, state[DIR])
        append!(xs, size(state[BODY]))
        append!(xs, state[HUNGER])
        return (xs, [goal])
    end

    function data_sensory(state, border_look_ahead=5)
        dir = state[DIR]
        fx, fy = state[FOOD]
        hx, hy = state[HEAD]
        body = state[BODY]
        data = []
        for alpha in dir-90:90:dir+90
            sensed = [0.0, 0.0, 0.0]
            for j in 0:(max(WIDTH, HEIGHT)/TILE_SIZE)
                mul = 1.0
                x = hx + j * TILE_SIZE * cos(deg2rad(alpha))
                y = hy + j * TILE_SIZE * cos(deg2rad(alpha))
                if isapprox(x, fx) && isapprox(y, fy)
                    sensed[1] = 1 / (j + 1)
                end
                if ((0 <= x < WIDTH) && (0 <= y < HEIGHT))
                    sensed[2] = 1 / (j + 1)
                end
                if j > 0
                    for part in body
                        if isapprox(x, part[1]) && isapprox(y, part[2])
                            sensed[3] = 1 / (j + 1)
                        end
                    end
                end
            end
            append!(data, sensed)

        end
        #Computing how many neighbouring tiles from the head are occupied
        #by the snake's body.
        radius = 3
        for r in 1:radius
            closeness = 0
            for tx in (hx - r * TILE_SIZE):TILE_SIZE:(hx + r * TILE_SIZE)
                for ty in (hy - r * TILE_SIZE):TILE_SIZE:(hy + r * TILE_SIZE)
                    for (px, py) in body[2:end] #2 to not take the head into account
                        if isapprox(px, tx) && isapprox(py, ty)
                            closeness += 1
                        end
                    end
                end
            end
            append!(data, closeness)
        end
        xs, y = data_simple(state)
        data = vcat(data, xs)
        return (vcat(data, state[HIST]), y)

    end
    export data_sensory
    export data_simple
end

ext: jl

max_stars_repo_name: evrardco/neural-snake

max_issues_repo_name: evrardco/neural-snake

max_forks_repo_name: evrardco/neural-snake


text:
import Base: display, redisplay

immutable InlineDisplay <: Display end

# supported MIME types for inline display in IPython, in descending order
# of preference (descending "richness")
const ipy_mime = [ "application/vnd.vegalite.v2+json", "text/html", "text/latex", "image/svg+xml", "image/png", "image/jpeg", "text/plain", "text/markdown", "application/javascript" ]

# need special handling for showing a string as a textmime
# type, since in that case the string is assumed to be
# raw data unless it is text/plain
israwtext(::MIME, x::AbstractString) = true
israwtext(::MIME"text/plain", x::AbstractString) = false
israwtext(::MIME, x) = false

# convert x to a string of type mime, making sure to use an
# IOContext that tells the underlying show function to limit output
function limitstringmime(mime::MIME, x)
    buf = IOBuffer()
    if istextmime(mime)
        if israwtext(mime, x)
            return String(x)
        else
            show(IOContext(buf, :limit=>true, :color=>true), mime, x)
        end
    else
        b64 = Base64EncodePipe(buf)
        if isa(x, Vector{UInt8})
            write(b64, x) # x assumed to be raw binary data
        else
            show(IOContext(b64, :limit=>true, :color=>true), mime, x)
        end
        close(b64)
    end
    return String(take!(buf))
end

for mime in ipy_mime
    @eval begin
        function display(d::InlineDisplay, ::MIME{Symbol($mime)}, x)
            send_ipython(publish[],
                         msg_pub(execute_msg, "display_data",
                                 Dict(
                                  "metadata" => metadata(x), # optional
                                  "data" => Dict($mime => limitstringmime(MIME($mime), x)))))
        end
        displayable(d::InlineDisplay, ::MIME{Symbol($mime)}) = true
    end
end

# deal with annoying application/x-latex == text/latex synonyms
display(d::InlineDisplay, m::MIME"application/x-latex", x) = display(d, MIME("text/latex"), limitstringmime(m, x))

# deal with annoying text/javascript == application/javascript synonyms
display(d::InlineDisplay, m::MIME"text/javascript", x) = display(d, MIME("application/javascript"), limitstringmime(m, x))

# If the user explicitly calls display("foo/bar", x), we send
# the display message, also sending text/plain for text data.
displayable(d::InlineDisplay, M::MIME) = istextmime(M)
function display(d::InlineDisplay, M::MIME, x)
    sx = limitstringmime(M, x)
    d = Dict(string(M) => sx)
    if istextmime(M)
        d["text/plain"] = sx # directly show text data, e.g. text/csv
    end
    send_ipython(publish[],
                 msg_pub(execute_msg, "display_data",
                         Dict("metadata" => metadata(x), # optional
                              "data" => d)))
end

# override display to send IPython a dictionary of all supported
# output types, so that IPython can choose what to display.
function display(d::InlineDisplay, x)
    undisplay(x) # dequeue previous redisplay(x)
    send_ipython(publish[],
                 msg_pub(execute_msg, "display_data",
                         Dict("metadata" => metadata(x), # optional
                              "data" => display_dict(x))))
end

# we overload redisplay(d, x) to add x to a queue of objects to display,
# with the actual display occuring when display() is called or when
# an input cell has finished executing.

function redisplay(d::InlineDisplay, x)
    if !in(x,displayqueue)
        push!(displayqueue, x)
    end
end

function display()
    q = copy(displayqueue)
    empty!(displayqueue) # so that undisplay in display(x) is no-op
    for x in q
        display(x)
    end
end


text:
#   Copyright 2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
Unit tests to test methods in `jltest/utils.jl`
"""

# --- Imports

using Logging
using Test

# --- Tests

@testset "@warn message tests" begin
    @warn(
        "Single line @warn message test",
        _module = nothing,
        _file = nothing,
        _group = nothing
    )

    @warn """
          Multi-line @warn message test.
          Second line.
          Third line.
          """
end
@testset "@info message tests" begin
    @info(
        "Single line @info message test",
        _module = nothing,
        _file = nothing,
        _group = nothing
    )

    @info """
          Multi-line @info message test.
          Second line.
          Third line.
          """
end

@testset "@debug message tests" begin
    @debug(
        "Single line @debug message test",
        _module = nothing,
        _file = nothing,
        _group = nothing
    )

    @debug """
           Multi-line @debug message test.
           Second line.
           Third line.
           """
end


text:
using Random: Float64
using LinearAlgebra
using Statistics
using Random
using Printf

# lattice parameters
struct Lattice
    Ns::Int64
    nup::Int64
    ndn::Int64
    t::Float64
    U::Float64
    KMatrix::Array{Float64,2}
end

# simulation parameters
struct QMC
    block_size::Int64
    sample_size::Int64
    cluster_size::Int64
    Δτ::Float64
    L::Int64
    γ::Float64
    μ::Float64
    KMatrixExp::Array{Float64,2}
end

const FloatType = Union{Float64, ComplexF64}

include("./2D/MatrixGenerator.jl")
include("./2D/Recursion.jl")
include("./2D/Calculations.jl")
include("./2D/Test.jl")

function get_inputs()
    lines = readlines("./InputFermiHubbard2D")

    # basic parameters
    β = parse(Float64, split(lines[1])[2])
    num_sitesx = parse(Int64, split(lines[2])[2])
    num_sitesy = parse(Int64, split(lines[3])[2])
    num_spinup = parse(Int64, split(lines[4])[2])
    num_spindn = parse(Int64, split(lines[5])[2])
    num_sites = num_sitesx * num_sitesy

    Δτ = parse(Float64, split(lines[6])[2])
    t = parse(Float64, split(lines[7])[2])
    U = parse(Float64, split(lines[8])[2])
    num_blocks = convert(Int64, parse(Float64, split(lines[9])[2]))
    num_samples = convert(Int64, parse(Float64, split(lines[10])[2]))

    # optimization parameters
    stab_interval = parse(Int64, split(lines[13])[2])

    # constants for calculations
    L = convert(Int64, round(β / Δτ))
    γ = atanh(sqrt(tanh(Δτ * U / 4)))
    mu = convert(Float64, parse(Float64, split(lines[16])[2]))
    KMatrix = kinetic_matrix_generator(num_sitesx, num_sitesy, t)

    lattice = Lattice(num_sites, num_spinup, num_spindn, t, U, KMatrix)
    qmc = QMC(num_blocks, num_samples, stab_interval,
    Δτ, L, γ, mu, exp(-KMatrix * Δτ / 2))

    lattice, qmc
end

function main()
    lattice, qmc = get_inputs()

    BlockEk_Array = Vector{Float64}()
    BlockEp_Array = Vector{Float64}()
    BlockE_Array = Vector{Float64}()
    BlockEkVar_Array = Vector{Float64}()
    BlockEpVar_Array = Vector{Float64}()
    BlockEVar_Array = Vector{Float64}()
    BlockSgn_Array = Vector{Float64}()

    filename = string("output_Ns", lattice.Ns, "_N", lattice.nup + lattice.ndn,
                    "_U", lattice.U, "_beta", qmc.Δτ * qmc.L)
    outputfile = open(filename, "a+")

    @printf(outputfile, "       Canonical Ensemble AFQMC for 2D Hubbard Model     \n\n")
    @printf(outputfile, "Number of Electrons: %d\n", lattice.nup + lattice.ndn)
    @printf(outputfile, "Number of Blocks: %d\n", qmc.block_size)
    @printf(outputfile, "Length of Img Timestep: %.4f\n", qmc.Δτ)
    @printf(outputfile, "Stablization Interval: %d\n\n", qmc.cluster_size)

    for i = 1 : qmc.block_size
        BlockE, BlockEVar, BlockEk, BlockEkVar, BlockEp, BlockEpVar, avgsgn = mc_integration(lattice, qmc)
        BlockE = BlockE / avgsgn
        BlockEk = BlockEk / avgsgn
        BlockEp = BlockEp / avgsgn

        @printf(outputfile, "Block # %d out of %d\n", i, qmc.block_size)
        @printf(outputfile, "--------------------\n")
        @printf(outputfile, "Average Sign: %.3f\n", avgsgn)
        @printf(outputfile, "Ek            Ep            Etot\n")
        @printf(outputfile, "%.5f     %.5f       %.5f\n", BlockEk, BlockEp, BlockE)
        @printf(outputfile, "Error:\n")
        @printf(outputfile, "%.5f     %.5f       %.5f\n\n", sqrt(BlockEkVar), sqrt(BlockEpVar), sqrt(BlockEVar))

        push!(BlockEk_Array, BlockEk)
        push!(BlockEkVar_Array, BlockEkVar)
        push!(BlockEp_Array, BlockEp)
        push!(BlockEpVar_Array, BlockEpVar)
        push!(BlockE_Array, BlockE)
        push!(BlockEVar_Array, BlockEVar)
        push!(BlockSgn_Array, avgsgn)
    end

    EkMean = mean(BlockEk_Array)
    EkVar = sum(BlockEkVar_Array * (qmc.sample_size - 1) .+
    qmc.sample_size * (BlockEk_Array .- EkMean) .^ 2) / (qmc.block_size * qmc.sample_size - 1)
    EpMean = mean(BlockEp_Array)
    EpVar = sum(BlockEpVar_Array * (qmc.sample_size - 1) .+
    qmc.sample_size * (BlockEp_Array .- EpMean) .^ 2) / (qmc.block_size * qmc.sample_size - 1)
    EMean = mean(BlockE_Array)
    EVar = sum(BlockEVar_Array * (qmc.sample_size - 1) .+
    qmc.sample_size * (BlockE_Array .- EMean) .^ 2) / (qmc.block_size * qmc.sample_size - 1)

    @printf(outputfile, "--------------------\n")
    @printf(outputfile, "Results for this run:\n")
    @printf(outputfile, "Average Sign: %.3f\n", mean(BlockSgn_Array))
    @printf(outputfile, "Ek            Ep            Etot\n")
    @printf(outputfile, "%.5f     %.5f       %.5f\n", EkMean, EpMean, EMean)
    @printf(outputfile, "Error:\n")
    @printf(outputfile, "%.5f     %.5f       %.5f\n\n", sqrt(EkVar / qmc.block_size), sqrt(EpVar / qmc.block_size), sqrt(EVar / qmc.block_size))

    close(outputfile)

end

#Random.seed!(1234)

main()
#main_gc()

max_stars_repo_name: TongSericus/CanonicalAFQMC

max_issues_repo_name: TongSericus/CanonicalAFQMC

max_forks_repo_name: TongSericus/CanonicalAFQMC


text:
# Native execution support

export @cuda, nearest_warpsize

using Base.Iterators: filter

#
# Auxiliary
#

# Determine which type to pre-convert objects to for use on a CUDA device.
#
# The resulting object type will be used as a starting point to determine the final
# specialization and argument types (there might be other conversions, eg. factoring in the
# ABI). This is different from `cconvert` in that we don't know which type to convert to.
function convert_type(t)
    # NOTE: this conversion was originally intended to be a user-extensible interface,
    #       a la cconvert (look for cudaconvert in f1e592e61d6898869b918331e3e625292f4c8cab).
    #
    #       however, the generated function behind @cuda isn't allowed to call overloaded
    #       functions (only pure ones), and also won't be able to see functions defined
    #       after the generated function's body (see JuliaLang/julia#19942).

    # Pointer handling
    if t <: DevicePtr
        return Ptr{t.parameters...}
    elseif t <: Ptr
        throw(InexactError())
    end

    # Array types
    if t <: CuArray
        return CuDeviceArray{t.parameters...}
    end

    return t
end

"""
Convert the arguments to a kernel function to their CUDA representation, and figure out what
types to specialize the kernel function for and how to actually pass those objects.
"""
function convert_arguments(args, types)
    argtypes = DataType[types...]
    argexprs = Union{Expr,Symbol}[args...]

    # convert types to their CUDA representation
    for i in 1:length(argexprs)
        t = argtypes[i]
        ct = convert_type(t)
        if ct != t
            argtypes[i] = ct
            if ct <: Ptr
                argexprs[i] = :( Base.unsafe_convert($ct, $(argexprs[i])) )
            else
                argexprs[i] = :( convert($ct, $(argexprs[i])) )
            end
        end
    end

    # figure out how to codegen and pass these types
    codegen_types, call_types = Array{DataType}(length(argtypes)), Array{Type}(length(argtypes))
    for i in 1:length(argexprs)
        codegen_types[i], call_types[i] = actual_types(argtypes[i])
    end

    # NOTE: DevicePtr's should have disappeared after this point

    return argexprs, (codegen_types...), (call_types...)
end

# NOTE: keep this in sync with jl_is_bitstype in julia.h
isbitstype(dt::DataType) =
    !dt.mutable && dt.layout != C_NULL && nfields(dt) == 0 && sizeof(dt) > 0

"""
Determine the actual types of an object, that is, 1) the type that needs to be used to
specialize (compile) the kernel function, and 2) the type which an object needs to be
converted to before passing it to a kernel.

These two types can differ, eg. when passing a bitstype that doesn't fit in a register, in
which case we'll be specializing a function that directly uses that bitstype (and the
codegen ABI will figure out it needs to be passed by ref, ie. a pointer), but we do need to
allocate memory and pass an actual pointer since we perform the call ourselves.
"""
function actual_types(argtype::DataType)
    if argtype.layout != C_NULL && Base.datatype_pointerfree(argtype)
        # pointerfree objects with a layout can be used on the GPU
        cgtype = argtype
        # but the ABI might require them to be passed by pointer
        if isbitstype(argtype)
            calltype = argtype
        else
            calltype = Ptr{argtype}
        end
    else
        error("don't know how to handle argument of type $argtype")
    end

    # special-case args which don't appear in the generated code
    # (but we still need to specialize for)
    if !argtype.mutable && sizeof(argtype) == 0
        # ghost type, ignored by the compiler
        calltype = Base.Bottom
    end

    return cgtype::Type, calltype::Type
end

function emit_allocations(args, codegen_types, call_types)
    # if we're generating code for a given type, but passing a pointer to that type instead,
    # this is indicative of needing to upload the value to GPU memory
    kernel_allocations = Expr(:block)
    for i in 1:length(args)
        if call_types[i] == Ptr{codegen_types[i]}
            @gensym dev_arg
            alloc = quote
                $dev_arg = Mem.alloc($(codegen_types[i]))
                # TODO: we're never freeing this (use refcounted single-value array?)
                Mem.upload($dev_arg, $(args[i]))
            end
            append!(kernel_allocations.args, alloc.args)
            args[i] = dev_arg
        end
    end

    return kernel_allocations, args
end

function emit_cudacall(func, dims, shmem, stream, types, args)
    # TODO: can we handle non-isbits types?
    all(t -> isbits(t) && sizeof(t) > 0, types) ||
        error("can only pass bitstypes of size > 0 to CUDA kernels")
    any(t -> sizeof(t) > 8, types) &&
        error("cannot pass objects that don't fit in registers to CUDA functions")

    return quote
        Profile.@launch begin
            cudacall($func, $dims[1], $dims[2], $shmem, $stream, Tuple{$(types...)}, $(args...))
        end
    end
end

#
# @cuda macro
#

"""
    @cuda (gridDim::CuDim, blockDim::CuDim, [shmem::Int], [stream::CuStream]) func(args...)

High-level interface for calling functions on a GPU, queues a kernel launch on the current
context. The `gridDim` and `blockDim` arguments represent the launch configuration, the
optional `shmem` parameter specifies how much bytes of dynamic shared memory should be
allocated (defaulting to 0), while the optional `stream` parameter indicates on which stream
the launch should be scheduled.

The `func` argument should be a valid Julia function. It will be compiled to a CUDA function
upon first use, and to a certain extent arguments will be converted and managed
automatically. Finally, a call to `cudacall` is performed, scheduling the compiled function
for execution on the GPU.
"""
macro cuda(config::Expr, callexpr::Expr)
    # sanity checks
    if config.head != :tuple || !(2 <= length(config.args) <= 4)
        throw(ArgumentError("first argument to @cuda should be a tuple (gridDim, blockDim, [shmem], [stream])"))
    end
    if callexpr.head != :call
        throw(ArgumentError("second argument to @cuda should be a function call"))
    end

    # handle optional arguments and forward the call
    # NOTE: we duplicate the CUDAdrv's default values of these arguments,
    #       because the kwarg version of `cudacall` is too slow
    stream = length(config.args)==4 ? esc(pop!(config.args)) : :(CuDefaultStream())
    shmem  = length(config.args)==3 ? esc(pop!(config.args)) : :(0)
    dims = esc(config)
    return :(generated_cuda($dims, $shmem, $stream, $(map(esc, callexpr.args)...)))
end

# Compile and execute a CUDA kernel from a Julia function
const func_cache = Dict{UInt, CuFunction}()
@generated function generated_cuda{F<:Core.Function,N}(dims::Tuple{CuDim, CuDim}, shmem, stream,
                                                       func::F, args::Vararg{Any,N})
    arg_exprs = [:( args[$i] ) for i in 1:N]
    arg_exprs, codegen_types, call_types = convert_arguments(arg_exprs, args)

    kernel_allocations, arg_exprs = emit_allocations(arg_exprs, codegen_types, call_types)

    # compile the function, once
    @gensym cuda_fun
    precomp_key = hash(tuple(func, codegen_types...))  # precomputable part of the key
    kernel_compilation = quote
        ctx = CuCurrentContext()
        key = hash(($precomp_key, ctx))
        if (haskey(func_cache, key))
            $cuda_fun = func_cache[key]
        else
            $cuda_fun, _ = cufunction(device(ctx), func, $codegen_types)
            func_cache[key] = $cuda_fun
        end
    end

    # filter out non-concrete args
    concrete = map(t->t!=Base.Bottom, call_types)
    call_types = map(x->x[2], filter(x->x[1], zip(concrete, call_types)))
    arg_exprs  = map(x->x[2], filter(x->x[1], zip(concrete, arg_exprs)))

    kernel_call = emit_cudacall(cuda_fun, :(dims), :(shmem), :(stream),
                                call_types, arg_exprs)

    quote
        Base.@_inline_meta
        $kernel_allocations
        $kernel_compilation
        $kernel_call
    end
end

"""
Return the nearest number of threads that is a multiple of the warp size of a device:

    nearest_warpsize(dev::CuDevice, threads::Integer)

This is a common requirement, eg. when using shuffle intrinsics.
"""
function nearest_warpsize(dev::CuDevice, threads::Integer)
    ws = CUDAdrv.warpsize(dev)
    return threads + (ws - threads % ws) % ws
end


text:
#
#--------------------------------------#
export pcg
#--------------------------------------#
"""
 Preconditioned conjugate gradient

 We can't just use any packaged iterative solver for the PCG because
 we need to overwrite not just A*x, A'*x ops for Ax=b which is
 complicated enough since x,b  are multidimensional arrays (would have
 to allocate memory in converting and reconverting arrays to vecs), but
 also <u,v>, the dot product. because of the SEM construction, the dot
 product isnt just sum(u.*v) , but it's sum(u.*v.* SEM_mult)

"""
function pcg(b,opA
            ;opM=x->x           # preconditioner
            ,mult=ones(size(b)) # SEM mult
            ,ifv=false          # verbose flag
            ,tol=1e-8           # tolerance
            ,maxiter=length(b)) # maximum number of iterations

n = length(b);

x   = zero(b) # zero out initial condition
Ax  = zero(b)
ra  = b - Ax;
ha  = zero(b)
hp  = zero(b)
hpp = zero(b)
rp  = zero(b)
rpp = zero(b)
u   = zero(b)
k   = 0

while(norm(ra,Inf) > tol)
ha = opM * ra # preconditinoer
#println("PCG iter: ",k,", res: ",norm(ra,2));
if(k==maxiter) println("warning: res:",norm(ra,Inf)); return x; end;
k  += 1;
hpp = hp;
rpp = rp;
hp  = ha;
rp  = ra;
t   = sum(rp.*hp.*mult);
if(k==1)
    u = copy(hp);
else
    u = hp+(t/sum(rpp.*hpp.*mult))*u;
end
Au = opA * u # operator
a = t / sum(u.*Au.*mult);
x = x + a * u;
ra = rp - a * Au;
end

if(ifv) println("PCG iter: ",k,", res: ",norm(ra,2)); end

return x
end
#--------------------------------------#
export pcg!
#--------------------------------------#
function pcg!(x,b,opA
             ;opM=x->x           # preconditioner
             ,mult=ones(size(b)) # SEM mult
             ,ifv=false          # verbose flag
             ,tol=1e-8           # tolerance
             ,maxiter=length(b)) # maximum number of iterations

    x .= pcg(b,opA
            ;opM=opM
            ,mult=mult
            ,ifv=ifv
            ,tol=tol
            ,maxiter=maxiter)

return
end
#--------------------------------------#


text:
using LogProbs
using Test

p = LogProb(0.2)
q = LogProb(0.5)
r = LogProb(0.1)

@test p * q ≈ r
@test float(p+q) ≈ 0.7
@test log(p) == -1.6094379124341003
@test information(q) ≈ 1
@test r * LogProb(2) ≈ p
@test q - p ≈ LogProb(.3)
@test p / q ≈ LogProb(.4)
@test p + zero(LogProb) ≈ p
@test p *  one(LogProb) ≈ p
@test p < q
@test LogProb(0) < p

@test one(LogProb).log == 0
@test zero(LogProb).log == -Inf

@test string(p) == "LogProb(0.2)"

@test let d = Dict(LogProb(x) => x for x in [0.2, 0.5, 0.1])
    (d[p], d[q], d[r]) == (0.2, 0.5, 0.1)
end

@test let s = rand(LogProb)
    s.log == log(s)
end


text:
"""
Read a text file consistent entirely of floating point numbers separated
by single spaces or newlines.

I saw that this was about 1.5-2x faster than CSV.jl, which was the fastest
of the existing read methods.
"""

@inline mytryparse(::Type{Float64}, s::Vector{UInt8}, pos::Int64, len::Int64) =
            ccall(:jl_try_substrtod, Nullable{Float64}, (Ptr{UInt8},Csize_t,Csize_t), s, pos, len)

@inline function myparse(::Type{Float64}, s::Vector{UInt8}, pos::Int64, last::Int64)
    result = mytryparse(Float64, s, pos-1, last-pos+1)
    if isnull(result)
        throw(ArgumentError("cannot parse $(repr(s)) as $Float64"))
    end
    if VERSION < v"0.6"
      # when introduced, the function of unsafe_get is
      # unsafe_get(x::Nullable) = x.value
      # unsafe_get(x) = x
      return result.value
    else
      return unsafe_get(result)
    end
end

@inline function myparse(::Type{Int32}, s::Vector{UInt8}, pos::Int64, last::Int64)
    val = myparse(Float64, s, pos, last)
    return convert(Int32, val)
end

"""
Find the next space in the buffer a, starting the search at position = pos
and continuing until we find a space or hit len. If we can't find a space,
return -1, otherwise, return the index of the space.
"""
@inline function next_space(a::Vector{UInt8}, pos, len)
    @inbounds for i = pos:len
        if a[i] == UInt8(' ') || a[i] == UInt8('\n') || a[i] == UInt8('\t')
            return i
        end
    end
    return -1
end

@inline function next_non_space(a::Vector{UInt8}, pos, len)
    @inbounds for i = pos:len
        if a[i] != UInt8(' ') && a[i] != UInt8('\n') && a[i] != UInt8('\t')
            return i
        end
    end
    return -1
end

@inline function next_space_with_leading(a::Vector{UInt8}, pos, len)
    skip = next_non_space(a, pos, len)
    skip != -1 && return next_space(a, skip, len)
    return skip
end

@inline function all_spaces(a::Vector{UInt8}, pos, len)
    @inbounds for i = pos:len
        if a[i] != UInt8(' ') && a[i] != UInt8('\n') && a[i] != UInt8('\t')
            return false
        end
    end
    return true
end

mydebug = false

function tseq2(io::IO, a; maxbuf::Int=2^10, ntok::Int=-1)
    buf = Vector{UInt8}(maxbuf) # 64k bytes
    curtok = 0
    curbuf = 1

    tokinc = 1
    if ntok == -1
        tokinc = 0
        ntok = 1
    end

    while !eof(io) && curtok < ntok
        x = read(io, UInt8)
        if curbuf > 1
            # then we are processing a token
            if x == UInt8(' ') || x == UInt8('\n') || x == UInt8('\t')
                # there is a token from buf[1:curbuf]
                #@show "token", String(buf[1:curbuf-1])
                push!(a, myparse(Float64, buf, 1, curbuf-1))
                curtok += tokinc
                curbuf = 1
            elseif curbuf > maxbuf
                throw(ArgumentError("current token exceeded maxbuf=$(maxbuf)"))
            else
                buf[curbuf] = x
                curbuf += 1
            end
        else
            # otherwise we are handling spaces
            if x == UInt8(' ') || x == UInt8('\n') || x == UInt8('\t')
                continue
            else
                buf[curbuf] = x
                curbuf += 1 # we are in a token now!
            end
        end
    end
    if eof(io) && curtok < ntok && curbuf > 1
        # @show "token", String(buf[1:curbuf-1])
        push!(a, myparse(Float64, buf, 1, curbuf-1))
    end
    return a
end

function tseq(io::IO, a; maxbuf::Int=2^10)
    buf = Vector{UInt8}(maxbuf) # 64k bytes
    buf2 = Vector{UInt8}(maxbuf) # 64k bytes
    eof = false
    nb = readbytes!(io, buf)
    if nb == 0
        return a # just return
    end
    (nb < maxbuf) && (eof = true)

    curspace = next_non_space(buf, 1, nb)

    @inbounds while nb >= 0 && curspace >= 0
        nextspace = next_space_with_leading(buf, curspace+1, nb)
        if nextspace >= 0
            mydebug && @show "parsing", String(buf[curspace:nextspace])
            push!(a, myparse(Float64, buf, curspace, nextspace))
        elseif curspace == 1 && eof == true
            push!(a, myparse(Float64, buf, curspace, nb))
            # We are at eof and just processed the last of the buffer
        elseif curspace == 1
            @show curspace, nextspace
            throw(ArgumentError("current value exceeded maxbuf=$(maxbuf)"))
        else
            # we didn't see a space, that means we need to read more buffer
            # move things to beginning
            mydebug && print("Refilling: ", replace(String(buf),"\n","@"), " ", curspace, " ", String(buf[curspace:nb]), "\n")
            copy!(buf, 1, buf, curspace, nb - curspace + 1)
            #print("     Move: ", replace(String(buf),"\n","@"), "\n")
            bufstart = nb - curspace + 2
            #buffree = maxbuf - bufstart + 1
            curspace = 1

            # Try reading
            #nread = readbytes!(io, @view(buf[bufstart:end]))
            nread = readbytes!(io, buf2, maxbuf-bufstart+1)
            copy!(buf, bufstart, buf2, 1, maxbuf-bufstart+1)
            (nread < maxbuf-bufstart+1) && (eof = true)

            mydebug && print("   Reload: ", replace(String(buf),"\n","@"), "\n")

            mydebug && @show nread

            if nread == 0
                # We couldn't read any more, that means we are at the end!
                #print("    Final: ", replace(String(buf),"\n","@"), "\n")
                mydebug && @show "final", String(buf[curspace:bufstart-1])
                if bufstart > 2 && !all_spaces(buf, curspace, bufstart-1)
                    push!(a, myparse(Float64, buf, curspace, bufstart-1))
                    return a
                end
            else
                nb = bufstart + nread - 1
                mydebug && @show String(buf[curspace:nb])
                continue
            end

        end
        curspace = nextspace
    end
    return a
end

function read_file_to_float64_array(filename::AbstractString)
	a = zeros(0)
	open(filename, "r") do fh
		return tseq(fh, a)
	end
end

function readarrayt(io; kwargs...)
    line = readline(io) # read the first line
    a = zeros(0) # allocate the output
    tseq(IOBuffer(line), a; kwargs...) # parse the first line
    ncols = length(a) # this is the number of columns
    tseq(io, a; kwargs...) # read the rest of the array
    nrows, nrem = divrem(length(a),ncols)
    if nrem != 0
        throw(ArgumentError("first line of file had $(ncols) entries," *
            " but file has $(length(a)) for $(nrows) rows and $(nrem) left"))
    else
        return reshape(a, ncols, nrows)
    end
end
readarray(io; kwargs...) = readarrayt(io; kwargs...)'

#=
function read_number_stream(io, record::Tuple{Type}, nrec; all=true)

end
=#


text:
export @stateMonitor_str
macro stateMonitor_str(index)
  # initialize recording
  str=""
  dim=length(eval(parse(index)))
  for k=1:length(vars)
    str*="net.$(vars[k])Record=Array{Float64,2}($dim,0)\n"
  end
  print(str)
  map(x->eval(parse(x)),split(str,"\n"))
  # construct function
  str="function record!(net::Network$timestr)\n"
  for k=1:length(vars)
    str*="net.$(vars[k])Record=[net.$(vars[k])Record net.$(vars[k])[$index]]\n"
  end
  str*="end\n"
  print(str)
  eval(parse(str))
end


text:
## Local potentials. Can be provided from external potentials, or from `model.atoms`.

# a local potential term. Must have the field `potential_values`, storing the
# potential in real space on the grid. If the potential is different in the α and β
# components then it should be a 4d-array with the last axis running over the
# two spin components.
abstract type TermLocalPotential <: Term end

@timing "ene_ops: local" function ene_ops(term::TermLocalPotential,
                                          basis::PlaneWaveBasis{T}, ψ, occ; kwargs...) where {T}
    potview(data, spin) = ndims(data) == 4 ? (@view data[:, :, :, spin]) : data
    ops = [RealSpaceMultiplication(basis, kpt, potview(term.potential_values, kpt.spin))
           for kpt in basis.kpoints]
    if :ρ in keys(kwargs)
        E = sum(total_density(kwargs[:ρ]) .* term.potential_values) * basis.dvol
    else
        E = T(Inf)
    end

    (E=E, ops=ops)
end

## External potentials

struct TermExternal <: TermLocalPotential
    potential_values::AbstractArray
end

"""
External potential from an analytic function `V` (in cartesian coordinates).
No low-pass filtering is performed.
"""
struct ExternalFromReal
    potential::Function
end

function (external::ExternalFromReal)(basis::PlaneWaveBasis{T}) where {T}
    pot_real = external.potential.(r_vectors_cart(basis))
    TermExternal(convert_dual.(T, pot_real))
end

"""
External potential from the (unnormalized) Fourier coefficients `V(G)`
G is passed in cartesian coordinates
"""
struct ExternalFromFourier
    potential::Function
end
function (external::ExternalFromFourier)(basis::PlaneWaveBasis{T}) where {T}
    unit_cell_volume = basis.model.unit_cell_volume
    pot_fourier = map(G_vectors_cart(basis)) do G
        convert_dual(complex(T), external.potential(G) / sqrt(unit_cell_volume))
    end
    TermExternal(G_to_r(basis, pot_fourier))
end

## Atomic local potential

struct TermAtomicLocal{AT} <: TermLocalPotential
    potential_values::AT
end

"""
Atomic local potential defined by `model.atoms`.
"""
struct AtomicLocal end
function (::AtomicLocal)(basis::PlaneWaveBasis{T}) where {T}
    model = basis.model

    # pot_fourier is <e_G|V|e_G'> expanded in a basis of e_{G-G'}
    # Since V is a sum of radial functions located at atomic
    # positions, this involves a form factor (`local_potential_fourier`)
    # and a structure factor e^{-i Gr}

    pot_fourier = map(G_vectors(basis)) do G
        pot = sum(model.atom_groups) do group
            element = model.atoms[first(group)]
            form_factor::T = local_potential_fourier(element, norm(model.recip_lattice * G))
            form_factor * sum(cis2pi(-dot(G, r)) for r in @view model.positions[group])
        end
        pot / sqrt(model.unit_cell_volume)
    end

    pot_real = G_to_r(basis, pot_fourier)
    TermAtomicLocal(pot_real)
end

@timing "forces: local" function compute_forces(::TermAtomicLocal, basis::PlaneWaveBasis{TT},
                                                ψ, occupation; ρ, kwargs...) where TT
    T = promote_type(TT, real(eltype(ψ[1])))
    model = basis.model
    recip_lattice = model.recip_lattice
    ρ_fourier = r_to_G(basis, total_density(ρ))

    # energy = sum of form_factor(G) * struct_factor(G) * rho(G)
    # where struct_factor(G) = cis(-2π G⋅r)
    forces = [zero(Vec3{T}) for _ in 1:length(model.positions)]
    for group in model.atom_groups
        element = model.atoms[first(group)]
        form_factors = [Complex{T}(local_potential_fourier(element, norm(G)))
                        for G in G_vectors_cart(basis)]
        for idx in group
            r = model.positions[idx]
            forces[idx] = _force_local_internal(basis, ρ_fourier, form_factors, r)
        end
    end
    forces
end

# function barrier to work around various type instabilities
function _force_local_internal(basis, ρ_fourier, form_factors, r)
    T = real(eltype(ρ_fourier))
    f = zero(Vec3{T})
    for (iG, G) in enumerate(G_vectors(basis))
        f -= real(conj(ρ_fourier[iG])
                  .* form_factors[iG]
                  .* cis2pi(-dot(G, r))
                  .* (-2T(π)) .* G .* im
                  ./ sqrt(basis.model.unit_cell_volume))
    end
    f
end


text:
#set the path where the modules are located
push!(LOAD_PATH, "./")

#A examle function ready to be called from java using Ju4ja
function greetings(name::String)
    str="Hi, $name"
    println(str)
    return str
end

# g=greetings("test")
# println(g)

#using Ju4ja and the example solver for linear programming
using Ju4ja
using RereDmlLpSolver

#start the Ju4ja server as a new coroutine
@async begin
    startServer()
end


text:
module Poisson

using TimerOutputs
using IterativeSolvers: cg
using Preconditioners: AMGPreconditioner, SmoothedAggregation
using SparseArrays

using Gridap
import Gridap: ∇

export poisson

# Define manufactured functions
u(x) = x[1] + x[2]
∇u(x) = VectorValue(1.0,1.0,0.0)
∇(::typeof(u)) = ∇u
f(x) = 0.0

# Define forms of the problem
a(v,u) = ∇(v)*∇(u)
l(v) = v*f

# Define norms to measure the error
l2(u) = inner(u,u)
h1(u) = a(u,u) + l2(u)

"""
n: number of elements per direction
"""
function poisson(n::Integer)

  reset_timer!()

  domain = (0,1,0,1,0,1)
  partition = (n,n,n)
  @timeit "model" model = CartesianDiscreteModel(domain,partition)

  order = 1
  @timeit "V" V = TestFESpace(
    model=model,dirichlet_tags="boundary", conformity=:H1,
    reffe=:Lagrangian, order=order, valuetype=Float64)

  @timeit "U" U = TrialFESpace(V,u)

  @timeit "trian" trian = Triangulation(model)

  degree = 2
  @timeit "quad" quad = CellQuadrature(trian,degree)

  @timeit "t_Ω" t_Ω = AffineFETerm(a,l,trian,quad)

  @timeit "op" op = AffineFEOperator(SparseMatrixCSC{Float64,Int32},V,U,t_Ω)

  A = get_matrix(op)
  b = get_vector(op)

  @timeit "AMG setup" p = AMGPreconditioner{SmoothedAggregation}(A)

  @timeit "PCG"  x = cg(A,b,verbose=true,Pl=p)

  @timeit "uh" uh = FEFunction(U,x)

  @timeit "e" e = u - uh

  @timeit "el2" el2 = sqrt(sum( integrate(l2(e),trian,quad) ))
  @timeit "eh1" eh1 = sqrt(sum( integrate(h1(e),trian,quad) ))

  @show el2
  @show eh1

  print_timer()

end

end # module


text:
module DEW

using DifferentialEquations

struct Solution
    t
    x
    orig
end

function solveODE(helper_function, initial_conditions, interval, parameters=[])
    initial_conditions = convert.(Float64, initial_conditions)
    interval = convert.(Float64, interval)
    interval = tuple(interval...)

    prob = ODEProblem(helper_function, initial_conditions, interval, parameters)
    sol = solve(prob)

    solution = Solution(sol.t, sol[:,:]', sol)
    return solution
end

export solveODE

end # module

max_stars_repo_name: Balta-zar/DifferentialEquationsWrapper

max_issues_repo_name: Balta-zar/DifferentialEquationsWrapper

max_forks_repo_name: Balta-zar/DifferentialEquationsWrapper


text:
module My
using Parameters
using Random

@with_kw mutable struct Config
  char = (skip='?',less='>',more='<',num='$',klass='!')
  str  = (skip="?",)
  some = (max=512,step=.5, cohen=.3, trivial=1.05)
  seed = 1
end

THE = Config()
Random.seed!(THE.seed)

# -------------------------------------------
same(s) = s
int(x)  = floor(Int,x)
any(a)  = a[ int(size(a) * rand()) + 1 ]

function say(i)
  s,pre="$(typeof(i)){",""
  for f in sort!([x for x in fieldnames(typeof(i))
                     if !("$x"[1] == '_')])
    g = getfield(i,f)
    s = s * pre * "$f=$g"
    pre=", "
  end
  print(s * "}")
end

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

adds!(init=[],i=Some) = incs!(i(),init, 1)
subs!(init=[],i=Some) = incs!(i(),init,-1)
add!(i,x)             = inc!( i  ,   x, 1)
sub!(i,x)             = inc!( i  ,   x,-1)

# my cols can do:
#    incs!, inc!, statel, fresh, mid, var
# my cols know about:
#    w,pos,txt,w,key,n
incs!(i,init=[],w=1) = begin [inc!(i,x,w) for x in init]; i end

function inc!(i,x,w=1)
  y=i.key(x)
  if y != THE.str.skip
    stale(i)
    i.n += w
    inc1!(i, y,w) end
end

# -------------------------------------------
@with_kw mutable struct Num
  pos=0; txt=""; w=1; key=same; n=0;
  lo=10^32; hi=-1*10^32; mu=0; m2=0; sd=nothing end

mid(i::Num)   = i.mu
stale(i::Num) = i.sd = nothing

function var(i::Num)
  if i.sd == nothing
    i.sd = i.n < 2 ? 0 : (i.m2 / (i.n - 1 + 10^-32))^0.5 end
  i.sd
end

function inc1!(i::Num,x)
  i.lo  = min(i.lo, x)
  i.hi  = max(i.hi, x)
  d     = x - i.mu
  i.mu += d / i.n
  i.m2 += d * (x - i.mu)
end

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

@with_kw mutable struct Some
  pos=0; txt=""; w=1; key=same; n=0;
  all=[]; max=THE.some.max ;tidy=false end

p(i::Some,n)   = begin fresh(i); i.all[int(n*length(i.all))+1] end
stale(i::Some) = i.tidy=false

function fresh(i::Some)
  if !i.tidy
    sort!(i.all)
    i.tidy=true end  end

has(i::Some,n) = begin fresh(i); i.all[n] end
mid(i::Some, lo=1, hi=length(i.all)) = has(i,int(lo+(hi-lo)*.5))

function var(i::Some,lo=1,hi=length(i.all))
  fresh(i)
  n10 = int(lo+(hi-lo)*.1) + 1
  n90 = int(lo+(hi-lo)*.9) + 1
  (i.all[n90] - i.all[n10])/2.7
end

function inc1!(i::Some, x,w=1)
  m = length(i.all)
  if m < i.max
    push!(i.all,x)
  elseif rand() < m/i.n
    i.all[ int(m*rand()) + 1 ] = x end
end

"If i.all is broken at the points listed in `a`
between `lo` and `hi`, what is the expected value?"
function xpect(i::Some,a,lo=1,hi=length(i.all))
  e1(x,y) = (y-x+1)/(hi-lo+1)*var(i,x,y)
  e,m = 0,lo
  for n in a
    e += e1(m,n)
    m = n+1
  end
  e + e1(m,hi)
end

div(i::Some) = begin fresh(i); div(i.all,i.key) end

# -------------------------------------------
@with_kw mutable struct Range
  lo=0; hi=0; _all=[]; start=0; stop=0; w=0; _kids=[] end

Base.show(io::IO, i::Range) = say(i)

"assumes lst is sorted"
function div(lst::Array,key=same)
  the = THE.some
  x(z)           = key(lst[int(z)])
  val(y,z,p=0.5) = x(y+(z-y)*p)
  var(y,z)       = (val(y,z,0.9) - val(y,z,0.1))/2.7
  function xchop(lo,hi,out=nothing)
    best = var(lo,hi)
    for j = lo+step:hi-step
      now, after = x(j), x(j+1)
      if now != after
        if after - start > epsilon
	  if stop - now > epsilon
            if abs(val(lo,j) - val(j+1,hi)) > epsilon
	      n1,n2 = j-lo+1, hi-j
              here  = (var(lo,j)*n1 + var(j+1,hi)*n2)/(n1+n2)
              if here*the.trivial < best
                best,out = here,j end end end end end end
    return out
  end
  function xchops(lo,hi,ranges, cut = chop(lo,hi))
    if cut == nothing
      push!(ranges, Range(lo=x(lo), hi=x(hi),
                          _all=lst[lo:hi],start=lo,stop=hi))
    else
      xchops(lo,    cut, ranges)
      xchops(cut+1, hi,  ranges) end
  end
  #----------------------------------------------
  n                 = length(lst)
  epsilon           = var(1,n) * the.cohen
  step, start, stop = int(n^the.step)-1, x(1), x(n)
  xchops(1,n,[])
end

function chops(lo,hi,ranges,chop)
    cut = chop(lo,hi)
    if cut == nothing
      push!(ranges, Range(lo=x(lo), hi=x(hi),
                          _all=lst[lo:hi],start=lo,stop=hi))
    else
      ychops(lo,    cut, ranges)
      ychops(cut+1, hi,  ranges) end
  end
end
#
#function unite(rs, y=same,better= <, yis=Num)
#  the = THE.some
#  all(x=yis(key=y),a=[])= begin [incs!(x,r._all) for r in a]; x end
#  function ychop(lo,hi,best,rs,out=nothing)
#    left = yis(key=y)
#    for j in lo:hi-1
#      l= all(x=left,[rs[j]])
#      rall(a=rs[j+1:hi])
#      now = (var(l)*l.n + var(r)*r.n)/(l.n + r.n)
#      if better(now*the.trivial, best)
#        best,out = now,j end end
#    out
#  end
#  f = (start,stop) -> ychop(start,stopr,)
#  chop(1,length(rs),[], var(all(ranges)))
#end
#
# -------------------------------------------
@with_kw mutable struct Sym
  pos=0; txt=""; w=1; key=same; n=0;
  seen=Dict();  mode=nothing; ent=nothing;  end

mid(i::Sym) = begin i.fresh(); i.mode end
var(i::Sym) = begin i.fresh(); i.ent  end

stale(i::Sym) = i.mode,i.ent = nothing,nothing
function fresh(i::Sym)
  if i.mode == nothing
    i.ent, most = 0,0
    for (k,n) in i.seen
      p = n/i.n
      i.ent -= p*log(2,p)
      if n > most most,i.mode = n.k end end end
 end

function inc1!(i::Sym,x,w=1)
  new = w + (haskey(i.seen, x) ? i.seen[x] : 0)
  i.seen[x] = max(new,0)
end

# --------------------------------------------
norm(i::Sym, x) = x
norm(i::Some,x) = begin fresh(i); (x-i.all[1])/(i.all[end]-i.all[1]) end

difference(i::Sym, x,y) = x==THE.string.skip ? 1 : x == y
function difference(i::Some,x,y, no = THE.string.skip)
  d(a,b) = begin a= norm(i,a); b= a<0.5 ? 1 : 0; abs(a-b) end
  if     x==no && y==no 1
  elseif x==no          d(y,x)
  elseif y==no          d(x,y)
  else                  abs(norm(i,x) - norm(i,y)) end
end

# -------------------------------------------
@with_kw struct Lines file; src=open(file) end

"Define an iterator that returns a comma-seperated file, one
 record at a time without loading the whole file into memory."
function Base.iterate(it::Lines, (n,want)=(1,[]))
  "Split on comma, coerce strings to numbers or strings, as approriate."
  coerce(s)  = map(coerce1, split(s,","))
  coerce1(s) = ((x = tryparse(Float64,s))==nothing) ? s : x

  "Coerce strings. If first row, check what columns we should use.
   Only return those columns."
  function cols(a)
    if n == 1
      want = [i for (i,s) in enumerate(a) if !('?' in s)] end
    [a[i] for i in want]
  end

  "Delete comments and whitespace. Lines ending in
   ',' are joined to the next. Skip empty lines."
  function row(txt="")
    while true
      if eof(it.src) return txt end
      new = readline(it.src)
      new = replace(new, r"([ \t\n]|#.*)"=>"")
      if sizeof(new) != 0
        txt *= new
        if txt[end] != ','
	  return txt end end end end

  new = row()
  if sizeof(new) > 0
    (n, cols(coerce(new))) , (n+1,want) end
 end

#--------------------------------------------
id=0

@with_kw mutable struct Tbl
  rows=[]; cols=Cols() end

@with_kw mutable struct Row
  cells=[]; cooked=[]; id=global id+= 1
end

say(Row())
say(Row())
@with_kw mutable struct Cols
  x = (all=[], nums=[], syms=[])
  y = (all=[], nums=[], syms=[], goals=[])
  klass=""
  all  = []; nums = []; syms = []; end

function table(file::String)
  t=Tbl()
  for (n,a) in Lines(file=file)
    n==1 ? head!(t,a) : row!(t,a) end
  t
end

function row!(i::Tbl,a)
  [add!(c,a[c.pos]) for c in i.cols.all]
  push!(i.rows, Row(cells=a) )
end

head!(i::Tbl,a) = [head!(i.cols,n,x) for (n,x) in enumerate(a)]

function head!(i::Cols, n,txt)
  the = THE.char
  goalp()  = the.less in txt || the.more in txt
  nump()   = the.num  in txt || goalp()
  yp()     = klassp() || goalp()
  klassp() = the.klass in txt
  x = nump() ? Some : Sym
  y = x(pos=n, txt=txt)
  if klassp() i.klass = y end
  if goalp()  push!(i.y.goals, y) end
  if nump()
    push!(i.nums,y); push!(yp() ? i.y.nums : i.x.nums, y)
  else
    push!(i.syms,y); push!(yp() ? i.y.syms : i.x.syms, y)
  end
  push!(yp()  ? i.y.all : i.x.all, y)
  push!(i.all, y)
end

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

function tbl1(f="data/auto.csv")
  t = table(f)
  println("n ",length(t.rows))
  for col in t.cols.x.nums
   println(div(col)) #println(var(col)," ",col.all)
  end
end

function nums(f="data/auto.csv")
  t = table(f)
  #println(t.rows[end].cells)
  for num in t.cols.x.nums
    d=div(num)
    println(num.txt, " ",length(d))
    println(d)
  end
end

function sym1()
  s=Sym()
  [add!(s,x) for x in "aaaabbc"]
end

function Lines1(f="data/weather.csv")
  m=1
  print(m)
  for (n,tmp) in Lines(file= f)
     m += sizeof(tmp)  #println(n," ",tmp)
     if mod(n,1000) == 0 println(n,":",m) end
  end
  print(m)
end

function num1(x)
  if x<0.3 return 0.1 end
  if x<0.7 return 0.8 end
  return 0.9
 end

function numbers1(s=Some())
  [add!(s,num1(rand())) for i in 1:100]
  println([has(s,i) for i in div(s)])
end

function numbers2(n=2, s=Some())
  [add!(s,rand()^0.5) for i in 1:10^n]
  println([(i,has(s,i)) for i in div(s)])
end

#some1()
#sym1()
#@time tbl1("data/xomo10000.csv")
@time tbl1("data/weather.csv")
#@time nums("data/xomo10000.csv")
#numbers1()
#@time numbers2(3)
end

max_stars_repo_name: timm/misc

max_issues_repo_name: timm/misc

max_forks_repo_name: timm/misc


text:
using DashBootstrapComponents

badge = dbc_button(
    [
        "Notifications",
        dbc_badge("4", color = "light", text_color = "primary", className = "ms-1"),
    ],
    color = "primary",
);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
using Test
using ForwardDiff
using IndirectTrajOpt

# Initialize parameters
ps = initCR3BPIndirectWithSTMParams("Low Thrust 10 CR3BP")

# Initial condition vector
y0 = [-0.019488511458668, -0.016033479812051, 0.0,
       8.918881923678198, -4.081793688818725, 0.0, 1.0,
       15.616017, 32.875896, -0.094522,
       -0.101606, 0.044791, -0.000150, 0.133266]

# Initial condition vector with STM
z0 = [y0; zeros(196)]
for i in 1:14
    z0[14 + i + (i-1)*14] = 1.0
end

ps.ϵ = 1.0
for utype in [0,1,2]
    ps.utype = utype
    # Evaluate eoms with Jacobian
    dz0 = zeros(210)
    IndirectTrajOpt.cr3bpEomIndirectWithSTM!(dz0, z0, ps, 0.0)

    # Evaluate Jacobian of eoms with ForwardDiff
    dy0 = zeros(14)
    jac = ForwardDiff.jacobian(
        (y,x)->IndirectTrajOpt.cr3bpEomIndirect!(y,x,ps,0.0),
        dy0, y0)
    ps.utype
    # Evaluate jacobian difference
    jacDiff = abs.(reshape(dz0[15:end], (14,14)) .- jac)

    # Tests
    for i in 1:196
        @test jacDiff[i] < 1e-9
    end
end

ps.ϵ = 0.0
for utype in [0,2]
    ps.utype = utype
    # Evaluate eoms with Jacobian
    dz0 = zeros(210)
    IndirectTrajOpt.cr3bpEomIndirectWithSTM!(dz0, z0, ps, 0.0)

    # Evaluate Jacobian of eoms with ForwardDiff
    dy0 = zeros(14)
    jac = ForwardDiff.jacobian(
        (y,x)->IndirectTrajOpt.cr3bpEomIndirect!(y,x,ps,0.0),
        dy0, y0)
    ps.utype
    # Evaluate jacobian difference
    jacDiff = abs.(reshape(dz0[15:end], (14,14)) .- jac)

    # Tests
    for i in 1:196
        @test jacDiff[i] < 1e-9
    end
end

# Define function for test
function cr3bpAnalyticJacTestFunc(x, tspan, ps)
    out, extra = cr3bpOptIntegrate(x, tspan, ps)
    return out
end

# Integrate state and co-states with STM
for ϵ in [1.0, 0.5, 0.0]
    ps.ϵ = ϵ
    tspan = (0.0, 8.6404*24*3600/ps.crp.TU)
    zf = cr3bpOptWithSTMIntegrate(z0, tspan, ps)

    # Compute STM with ForwardDiff
    stm = ForwardDiff.jacobian(
        (x)->cr3bpAnalyticJacTestFunc(x, tspan, ps), y0)

    # Evaluate STM difference
    jacDiff = abs.(reshape(zf[15:end], (14,14)) .- stm)
    display(jacDiff)

    # Tests
    for i in 1:196
        @test jacDiff[i] < 1e-2
    end
end


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

using CompilerTools
using Base.Test

## Tests for CompilerTools.LivenessAnalysis

function test_liveness_3(x::Int64, y::Int64, z::Int64)
	 A = random(x, y)
	 B = random(y, z)
	 C = zeros(x, z)
	 for i = 1:x
	     for j = 1:y
	     	 A[i,j] += 1;
	     end
	     for j = 1:z
	     	 for k = 1:y
		     C[i,j] += A[i,k] * B[k,j] + 1
		 end
	     end
	  end
end

ast_lv_2 = code_lowered(test_liveness_3, (Int64,Int64,Int64))[1]
cfg_2 = CompilerTools.CFGs.from_lambda(ast_lv_2) :: CompilerTools.CFGs.CFG

#CompilerTools.Loops.set_debug_level(3)

all_loops = CompilerTools.Loops.compute_dom_loops(cfg_2)

#println((all_loops.dom_dict))
@test (length(all_loops.loops) == 4)
@test (CompilerTools.Loops.isInLoop(all_loops, -1) == false)
# This is hard-coding a basic block number that is open to change by Julia.
#@test (CompilerTools.Loops.isInLoop(all_loops, 6) == true)


text:
using Test, Suppressor, DataFrames
using HMCResearchRandomWalks.ViewOut
using HMCResearchRandomWalks.ViewOut:
    isdfnull, isseednull, ispathnull, isfilenull

@testset "initialized CURRENT" begin
    @test isdfnull()
    @test isseednull()
    @test ispathnull()
    @test isfilenull()
end

@testset "current_df" begin
    path_err = "No DataFrame loaded. Run the getdata_...(...) method first."
    test_err = @capture_err begin
        current_df()
    end
    @test occursin(path_err, test_err)

    test_df = DataFrame(test_title = [1,2,3])
    current_df(test_df)
    @test test_df === current_df()
end

@testset "current_seed" begin
    path_err = "No seed loaded. Run the getdata_...(...[,seed]) method first."
    test_err = @capture_err begin
        current_seed()
    end
    @test occursin(path_err, test_err)

    test_seed = 4321
    current_seed(test_seed)
    @test "4321" === current_seed()
end

@testset "current_path" begin
    path_err = "No path loaded. Run the getdata_...(...) method first"
    test_err = @capture_err begin
        current_path()
    end
    @test occursin(path_err, test_err)

    test_path = "out/noflowdata/"
    current_path(test_path)
    @test test_path === current_path()
end

@testset "current_file" begin
    file_err = "No file loaded. Run savedata() or savetofolder() first."
    test_err = @capture_err begin
        current_file()
    end
    @test occursin(file_err, test_err)

    test_file = "C:/Users/username/.julia/dev/HMCResearchRandomWalks/out/test.csv"
    current_file(test_file)
    @test test_file === current_file()
end

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
using JLD2, Random, LinearAlgebra, Statistics, CSV, DataFrames, FreqTables, Distributions, Kronecker, JLD, GLM, HTTP, CategoricalArrays, TexTables

# Worked With Ahmed

function q1()
# 1,A
Random.seed!(1234)
A = rand(Uniform(-5,10),10,7)
B = rand(Normal(-2,15),10,7)
C = [A[1:5,1:5] B[1:5,6:7]]
D = (A .<= 0) .*A
# 1,B
size(A)
# 1,C
unique(D)
# 1,D
E = reshape(B,70,1)
E2 = vec(B)
# 1,E
F = reshape([A B], 10, 7, 2)
# 1,F
F = permutedims(F, [3, 1, 2])
#  1,G
G = kronecker(B,C)
# 1,H
save("matrixpractice.jld", "A", A, "B", B, "C", C, "D",
 D, "E", E, "F", F, "G", G )
 # 1,I
save("firstmatrix.jld", "A", A, "B", B, "C", C, "D", D)
# 1,J
CSV.write("Cmatrix.csv", DataFrame(C,:auto))
# 1,K
CSV.write("Dmatrix.dat", DataFrame(D,:auto), delim = ' ')
# 1,L

return A, B, C, D

end

A, B, C, D = q1()

function q2(A, B, C)

#2,A
AB = [A[i,j]*B[i,j] for i in 1:10, j in 1:7]
#2,B
Cprime = [C[i,j] for i in 1:5, j in 1:7 if -5 <= C[i,j]<= 5]

C[-5 .<= C .<= 5]
#2,C
N, K, T = 15169, 6, 5

X = zeros(N,K,T);

x5 = rand(Binomial(20,0.6),N)
x6 = rand(Binomial(20,0.5),N)

    for k in 1:T
    X[:,1,k] = ones(N)

    if rand(1)[1] <= 75*(6-k)/5
            X[:,2,k] = ones(N)
    end

    X[:,3,k] = rand(Normal(15 + k - 1, 5*(k-1)),N)

    X[:,4,k] = rand(Normal(pi*(6-k)/3, 1/exp(1)),N)

    # Columns 1, 5, and 6 remain stationary over time
    X[:,5, k] = x5
    X[:,6, k] = x6
end

#2,D

β = zeros(K,T);
β[1,:] = [1:0.25:2;];
β[2,:] = [log(j) for j in 1:T];
β[3,:] = [-sqrt(j) for j in 1:T];
β[4,:] = [exp(j)-exp(j+1) for j in 1:T];
β[5,:] = [j for j in 1:T];
β[6,:] = [j/3 for j in 1:T];

β

#2,E

Y = [X[:,:,t] * β[:,t] + rand(Normal(0,0.36),N,1) for t in 1:T]

end
q2(A, B, C)

function q3()
        #3,A

       # cd("C:/Users/wat20/OneDrive/Documents/ProblemSets/PS1-julia-intro")
        df = CSV.read("nlsw88.csv", DataFrame)

        #is missing (df)
        #[print(string(":", x, ", ")) for x in names(df)]

        #3,B
        NotMarrPer = 1 - sum(df[:,"married"])/size(df)[1];
        GradPer =sum(df[:,"collgrad"])/size(df)[1];

        println("Percent not married = ", round(NotMarrPer, digits = 4) * 100,
        "\nCollege graduate percentage = ", round(GradPer, digits = 4)* 100)

        #3,C)
        #note: tabulate command is undefined
        # source of freq. tables:
        # https://stackoverflow.com/questions/34654489/rs-table-function-in-julia-for-dataframes

        prop(freqtable(df.race))

        # d)

        intrq = (describe(df,:q75)[:,2] - describe(df,:q25)[:,2])
        summarystats = [describe(df,:nmissing, :mean, :median, :std, :min, :max, :nunique) intrq]
        rename(summarystats, :x1=> :itrq)

        #3,E
        ## Cross tab between industry and occupation
        freqtable(df.industry, df.occupation)

        #3,F
        # source: https://dataframes.juliadata.org/stable/man/split_apply_combine/

        qq = df[:, [:race, :occupation, :wage]];
        gdf = groupby(qq, [:race, :occupation]);
        combine(gdf, :wage => mean)
end

q3()

#4
fmat = load("firstmatrix.jld")
#keys(fmat)
#fmat["B"]

function q4()

#4A-E
function matrixops(A, B)
    #element by element product of inputs, A and b
    #input matrices product
    #input matrices sum
            if(size(A) != size(B))
            println("inputs must have the same size")
            else
                return println(
                    "A .* B = ", A .* B,
                    "\n\nA' * B = ", A' * B,
                    "\n\nA + B = ", A + B)
            end
        end

matrixops(A, B)

#4,F
# we get the error message back matrixops(C, D)

# g)

    matrixops(
    convert(Array,df.ttl_exp ),
    convert(Array,df.wage)
    )
end

q4()


text:
#=

  From
  https://mathematica.stackexchange.com/questions/128945/problems-encoding-a-bayesian-network-with-just-five-nodes-using-probabilitydistr

  Example is from "page 53 in Probabilistic Graphical Models (2009), by Daphne Koller and Neir Friedman:"

  """
  The network has five nodes (random variables):

  Difficulty of a class taken by a student (0 = easy, 1 = hard)
  Intelligence of the student (0 = low, 1 = high)
  Grade achieved by the student (1 = A, 2 = B, 3 = C)
  SAT score of the student (0 = low, 1 = high)
  Letter of recommendation by the teacher (0 = False, 1 = True)
  We would like to use this network to do probabilistic inference (causal or evidential) like:
  "What is the probability of the student achieving an A, given that he is intelligent?"
  """

  For observation/evidence: intelligence = intelligence_high

  Summary Statistics
    parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
        Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

    difficulty    1.4083    0.4915     0.0049    0.0101   2853.3277    1.0002      526.7358
  intelligence    2.0000    0.0000     0.0000    0.0000         NaN       NaN           NaN
         grade    1.3519    0.6449     0.0064    0.0105   3016.3692    1.0005      556.8339
        letter    0.7597    0.4273     0.0043    0.0072   2982.1758    1.0003      550.5216
           sat    1.1925    0.3943     0.0039    0.0072   3001.3436    0.9999      554.0601

  Distributions of variable grade
  a          =>    7420  (0.742000)
  b          =>    1641  (0.164100)
  c          =>     939  (0.093900)

  Distributions of variable intelligence
  high       =>   10000  (1.000000)

  Distributions of variable difficulty
  easy       =>    5917  (0.591700)
  hard       =>    4083  (0.408300)

  Distributions of variable letter (num:0)
  1.00000 =>    7597  (0.759700)
  0.00000 =>    2403  (0.240300)

  Distributions of variable sat
  high       =>    8075  (0.807500)
  low        =>    1925  (0.192500)

  Cf ~/blog/intelligence_test.blog
     ~/webppl/intelligence_test.wppl

=#

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

@model function intelligence_test()
    difficulty_easy = 1
    difficulty_hard = 2

    intelligence_low = 1
    intelligence_high = 2
    difficulty ~ Categorical([0.6,0.4]) # [difficulty_easy,difficulty_hard]

    intelligence ~ Categorical([0.7,0.3]) # [intelligence_low, intelligence_high]

    grade_a = 1
    grade_b = 2
    grade_c = 3
    grades = [grade_a,grade_b,grade_c]
    grade ~
        (intelligence == intelligence_low && difficulty==difficulty_easy)  ? Categorical([0.3 , 0.4 , 0.3 ]) :
        (intelligence == intelligence_low && difficulty==difficulty_hard)  ? Categorical([0.05, 0.25, 0.7 ]) :
        (intelligence == intelligence_high && difficulty==difficulty_easy) ? Categorical([0.9 , 0.08, 0.02]) :
        (intelligence == intelligence_high && difficulty==difficulty_hard) ? Categorical([0.5 , 0.3 , 0.2 ]) : Dirac(grade_c);

    letter ~
        (grade == grade_a) ? flip(0.9) :
        (grade == grade_b) ? flip(0.6) :
        (grade == grade_c) ? flip(0.01) : Dirac(false)

    sat ~
        (intelligence == intelligence_high) ? Categorical([0.8,0.2]) : # [sat_high,sat_low]}
        Categorical([0.05,0.95])

    # What is the probability of the student achieving an A, given that he is intelligent?
    # true ~ Dirac(sat == sat_high);
    # true ~ Dirac(letter == false);
    true ~ Dirac(intelligence == intelligence_high);

end

model = intelligence_test()

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 100_000)
# chns = sample(model, PG(5), 10_000)
chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

show_var_dist_pct(chns, :grade,["a","b","c"])
show_var_dist_pct(chns, :intelligence,["low","high"])
show_var_dist_pct(chns, :difficulty,["easy","hard"])
show_var_dist_pct(chns, :letter)
show_var_dist_pct(chns, :sat,["high","low"])

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
module TestSugar

import AtomixCUDA

using Atomix: @atomic
using CUDA
using CUDA: @allowscalar
using Test

using ..Utils: cuda

# Not implemented:
#=
function test_get_set()
    A = CUDA.ones(Int, 3)
    cuda() do
        GC.@preserve A begin
            x = @atomic A[begin]
            @atomic A[end-2] = -x
        end
    end
    @test collect(A) == [-1, 1, 1]
end
=#

function test_inc()
    A = CUDA.ones(Int, 3)
    cuda() do
        GC.@preserve A begin
            @atomic A[begin] += 1
        end
    end
    @test collect(A) == [2, 1, 1]
end

end  # module