text:
#Demo 16: Clock tests
#-------------------------------------------------------------------------------

using CMDimCircuits
CMDimCircuits.@using_CData()

#Get a demo display:
include(CMDimCircuits.demoplotcfgscript)
#Normally use something like:
#CMDimData.@includepkg EasyPlotInspect

#==Constants
===============================================================================#
vvst = cons(:a, labels = set(yaxis="Amplitude (V)", xaxis="Time (s)"))
Tnsvst = cons(:a, labels = set(yaxis="Period (ns)", xaxis="Time (s)"))
dutyvst = cons(:a, labels = set(yaxis="Duty Cycle (%)", xaxis="Time (s)"))
dutyoutvsmaxphi = cons(:a, labels = set(yaxis="Duty Cycle (%)", xaxis="Maximum Phase Shift"))
delayattr = cons(:a,
	line = set(style=:solid, width=3),
	glyph = set(shape=:x, size=2),
)

#==Input data
===============================================================================#
tbit = 1e-9 #Bit period
osr = 20 #samples per bit
ncycles = 100
noiseamp = 100e-3
ϕmax=2pi*.1

#==Computations
===============================================================================#
#Warning: Base functions not yet robust in the presence of noise.
tck = 2*tbit
t = DataF1(0:(tbit/osr):((ncycles)*tck))

#Generate parameter sweeps:
sweeplist = PSweep[
	PSweep("phimax", 2pi*[-.05, 0, .1])
]

#Generate data:
ϕ = fill(DataHR, sweeplist) do ϕmax
	return ϕmax*cos(2pi*t/tck)
end
#ϕ = ϕmax*cos(2pi*t/tck)

noise = (rand(length(t)) .- 0.5)*(2*noiseamp)
noise = DataF1(t.x, noise)
ck = sin(2pi*t/tck+ϕ)+noise

#cycle-to-cycle period:
Tc2c = delta(xcross(ck, allow=CrossType(:rise)))

#Duty:
duty = measduty(ck)

stats = measckstats(ck, tck=tck)
kv = sort(collect(keys(stats)))
for i in 1:length(stats)
	k, v = (kv[i], stats[kv[i]])
	if maximum(v) < 10e-9
		println("\n", k, "/1p: ", v/1e-12)
	else
		println("\n", k, ": ", v)
	end
end

ϕmax = parameter(ϕ, "phimax")
println("\n", "maximum duty-cycle phase shift: ", ϕmax)

#==Helper functions
===============================================================================#
#TODO: supply eyeparam() as a FoldedAxis or Plot constructor???
eyeparam(tbit; teye=1.5*tbit, tstart=0) = cons(:a,
	xfolded = set(tbit, xstart=tstart, xmax=teye),
	xaxis = set(min=0, max=teye), #Force limits on exact data range.
)

#==Generate plot
===============================================================================#
p1 = push!(cons(:plot, vvst, title="Clock"),
	cons(:wfrm, ck),
)
p2 = push!(cons(:plot, Tnsvst, title="Period (cycle-to-cycle)"),
	cons(:wfrm, Tc2c/1e-9),
)
p3 = push!(cons(:plot, eyeparam(tck, teye=1*tck), vvst, title="Clock Eye"),
	cons(:wfrm, ck),
)
p4 = push!(cons(:plot, dutyvst, title="Duty"),
	cons(:wfrm, duty*100),
)
p5 = push!(cons(:plot, dutyoutvsmaxphi, title="Duty"),
	cons(:wfrm, stats[:mean_duty]*100, delayattr, label="mean"),
	cons(:wfrm, stats[:min_duty]*100, delayattr, label="min"),
	cons(:wfrm, stats[:max_duty]*100, delayattr, label="max"),
)

pcoll = push!(cons(:plotcoll, title="Clock Stats Tests"), p1, p2, p3, p4, p5)
	pcoll.ncolumns = 2

#==Display results in pcoll
===============================================================================#
EasyPlot.displaygui(pcoll)

#==Return pcoll to user (call evalfile(...))
===============================================================================#
pcoll


text:
using BinaryBuilder

name = "IntelOpenMP"
version = v"2018.0.3"

sources_win32 = [
    "https://anaconda.org/intel/openmp/2018.0.3/download/win-32/openmp-2018.0.3-intel_0.tar.bz2" =>
    "86ed603332ed7b4004e8a474943468589b222ef16d0d9aaf3ebb4ceaf743a39d",
]

sources_win64 = [
    "https://anaconda.org/intel/openmp/2018.0.3/download/win-64/openmp-2018.0.3-intel_0.tar.bz2" =>
    "0aee3d9debb8b1c2bb9a202b780c2b2d2179e4cee9158f7d0ad46125cf6f3fa2",
]

sources_macos = [
    "https://anaconda.org/intel/openmp/2018.0.3/download/osx-64/openmp-2018.0.3-intel_0.tar.bz2" =>
    "110b94d5ff3c4df66fc89030c30ad42378da02817b3962f14cb5c268f9d94dae"
]

sources_linux32 = [
    "https://anaconda.org/intel/openmp/2018.0.3/download/linux-32/openmp-2018.0.3-intel_0.tar.bz2" =>
    "f06edc0c52337658fd4b780d0b5c704b0ffb1c156dced7f5038c1ebbda3d891b",
]

sources_linux64 = [
    "https://anaconda.org/intel/openmp/2018.0.3/download/linux-64/openmp-2018.0.3-intel_0.tar.bz2" =>
    "cae3ef59d900f12c723a3467e7122b559f0388c08c40c332da832131c024409b",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
if [[ ${target} == *mingw* ]]; then
    mv Library/bin ${libdir}
else
    mv lib ${libdir}
fi
install_license info/*.txt
"""

# The products that we will ensure are always built
products = [
    LibraryProduct(["libiomp5", "libiomp5md"], :libiomp),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Install first for win32, then win64.  This will accumulate files into `products` and also wrappers into the JLL package.
non_reg_ARGS = filter(arg -> arg != "--register", ARGS)

include("../../fancy_toys.jl")

if should_build_platform("i686-w64-mingw32")
    build_tarballs(non_reg_ARGS, name, version, sources_win32, script, [Platform("i686", "windows")], products, [])
end
if should_build_platform("x86_64-w64-mingw32")
    build_tarballs(non_reg_ARGS, name, version, sources_win64, script, [Platform("x86_64", "windows")], products, [])
end
if should_build_platform("x86_64-apple-darwin14")
    build_tarballs(non_reg_ARGS, name, version, sources_macos, script, [Platform("x86_64", "macos")], products, [])
end
if should_build_platform("i686-linux-gnu")
    build_tarballs(non_reg_ARGS, name, version, sources_linux32, script, [Platform("i686", "linux")], products, [])
end
if should_build_platform("x86_64-linux-gnu")
    build_tarballs(ARGS, name, version, sources_linux64, script, [Platform("x86_64", "linux")], products, [])
end

ext: jl

max_stars_repo_name: mortenpi/Yggdrasil

max_issues_repo_name: mortenpi/Yggdrasil

max_forks_repo_name: mortenpi/Yggdrasil


text:
# Test naive implementation of entropy calculations towards efficient calculation of entropy on a manifold

using DocStringExtensions
using KernelDensityEstimate

include(joinpath(dirname(@__FILE__), "circularEntropyUtils.jl"))

using Gadfly, Colors
using Distributions
using Random

using TransformUtils
using Optim

const TU = TransformUtils

using ApproxManifoldProducts

# some test points to work with
pts = 0.3*randn(30)
pts = [pts; TU.wrapRad.(0.6*randn(70) .- pi)]
shuffle!(pts)

## Construct circular KDE

pc = kde!_CircularNaiveCV(pts)
getBW(pc)[:,1]

plotKDECircular(pc)

## Create two more densities

pts2 = TU.wrapRad.(0.1*randn(100).-pi.+0.5)
pts3 = TU.wrapRad.(0.1*randn(100).+pi.-0.5)

pc2 = kde!_CircularNaiveCV(pts2)
pc3 = kde!_CircularNaiveCV(pts3)

plotKDECircular([pc2;pc3])

#


text:
"""
    CircularLayout

A circular layout has a single field, `radius`.
"""
struct CircularLayout
    radius::Float64
end

"""
    CircularLayout()

Create a circular layout with a radius of 1.
"""
CircularLayout() = CircularLayout(1.0)

function nodeangle(x, y)
    hyp = sqrt(x * x + y * y)
    θ = x < 0.0 ? π - asin(y / hyp) : asin(y / hyp) + 2π
    return θ
end

nodeangle(n::NodePosition) = nodeangle(n.x, n.y)

"""
    position!(LA::CircularLayout, L::Dict{K,NodePosition}, N::T) where {T <: AbstractEcologicalNetwork} where {K}

Nodes will be positioned at equal distances along a circle, and nodes that
are densely connected will be closer to one another. This is an efficient
way to represent modular networks.

#### References

McGuffin, M.J., 2012. Simple algorithms for network visualization: A tutorial.
Tsinghua Science and Technology 17, 383–398.
https://doi.org/10.1109/TST.2012.6297585
"""
function position!(
    LA::CircularLayout, L::Dict{K,NodePosition}, N::T
) where {T<:AbstractEcologicalNetwork} where {K}
    S = richness(N)
    Θ = Dict([s => nodeangle(L[s]) for s in species(N)])
    for (i, n1) in enumerate(species(N))
        θ = L[n1].r * 2π / S
        sx, sy = cos(θ), sin(θ)
        nei = Set{EcologicalNetworks._species_type(N)}[]
        if n1 ∈ species(N; dims=2)
            nei = union(nei, N[:, n1])
        end
        if n1 ∈ species(N; dims=1)
            nei = union(nei, N[n1, :])
        end
        for (j, n2) in enumerate(nei)
            θ2 = L[n2].r * 2π / S
            sx = sx + cos(θ2)
            sy = sy + sin(θ2)
        end
        Θ[n1] = mean(nodeangle(sx, sy))
    end
    m = minimum(collect(values(Θ)))
    M = maximum(collect(values(Θ)))
    for n in species(N)
        Θ[n] = (Θ[n] - m) / (M - m) * 2π + rand() * 0.001
    end
    v = collect(values(Θ))
    r = ordinalrank(v)
    R = Dict(zip(v, r))
    for n in species(N)
        i = 2 * R[Θ[n]] * π / S
        x, y = LA.radius * cos(i), LA.radius * sin(i)
        L[n] = NodePosition(x, y)
    end
end


text:
function exponent(x::Float64)
    xu = reinterpret(Uint64,x)
    int(xu >> 52) & 0x07ff - 1023
end

function exponent(x::Float32)
    xu = reinterpret(Uint32,x)
    int(xu >> 23) & 0x00ff -126
end

function setexponent(x::Float64,k::Int)
    xu = reinterpret(Uint64,x)
    ku = uint64((k + 1023) & 0x07ff) << 52
    xu = (xu & 0x800f_ffff_ffff_ffff) | ku
    reinterpret(Float64,xu)
end


text:
"""
Re-parameterization routine
that modifies the fields
`χvp` and `χρ` of an input seismic model
using two input vectors.

# Arguments

* `mod::Seismic` : to be updated
* `x1::Array{Float64,2}` : contrast of inverse bulk modulus
* `x2::Array{Float64,2}` : contrast of inverse density
* `attribvec:::Vector{Symbol}` : [:χKI, :χρI]
"""
function Seismic_reparameterize!(
	              mod::Seismic,
		      x,
		      attribvec::Vector{Symbol}=[:χKI, :χρI, :null]
		      )
	iszero(mod) ? error("mod cannot be zero") : nothing
	nznx=prod(length.(mod.mgrid))
	(length(x)≠(count(attribvec.≠ :null)*nznx)) &&  error("size x")
	if(attribvec == [:χKI, :χρI, :null])
		@inbounds for i in 1:nznx
			K=inv(χ(x[i],mod.ref.KI,-1))
			ρ=inv(χ(x[nznx+i],mod.ref.ρI,-1))
			mod.χvp[i]=χ(sqrt(K*inv(ρ)),mod.ref.vp,1)
			mod.χρ[i]=χ(ρ,mod.ref.ρ,1)
		end
	elseif(attribvec == [:χKI, :null, :null])
		@inbounds for i in 1:nznx
			K=inv(χ(x[i],mod.ref.KI,-1))
			ρ=χ(mod.χρ[i],mod.ref.ρ,-1)
			mod.χvp[i]=χ(sqrt(K*inv(ρ)),mod.ref.vp,1)
			mod.χρ[i]=χ(ρ,mod.ref.ρ,1)
		end
	elseif(attribvec == [:χvp, :χρ, :null])
		@inbounds for i in 1:nznx; mod.χvp[i]=x[i]; end
		@inbounds for i in 1:nznx; mod.χρ[i]=x[nznx+i]; end
	elseif(attribvec == [:null, :χρ, :null])
		@inbounds for i in 1:nznx; mod.χρ[i]=x[i]; end
	elseif(attribvec == [:χvp, :null, :null])
		@inbounds for i in 1:nznx; mod.χvp[i]=x[i]; end
	else
		error("invalid attribvec")
	end
	return mod
end

"""
Input perturbations corresponding to `parameterization`.
Output perturbations of KI and rhoI, i.e., without the contrast.
`mod` is just used for reference values.
"""
function pert_reparameterize!(δxout, δx, mod, parameterization)
	nznx=prod(length.(mod.mgrid))
	fill!(δxout, 0.0)
	if(parameterization == [:χKI, :χρI, :null])
		@inbounds for i in 1:nznx
			δxout[i]=δx[i]*mod.ref.KI
			δxout[nznx+i]=δx[nznx+i]*mod.ref.ρI
		end
	elseif(parameterization == [:χKI, :null, :null])
		@inbounds for i in 1:nznx
			δxout[i]=δx[i]*mod.ref.KI
		end
	elseif(parameterization == [:χvp, :χρ, :null])
		@inbounds for i in 1:nznx
			δxout[i]=(-1.0*δx[nznx+i]-2.0*δx[i])*mod.ref.KI
			δxout[nznx+i]=(-1.0*δx[nznx+i])*mod.ref.ρI
		end
	elseif(parameterization == [:χvp, :null, :null])
		@inbounds for i in 1:nznx
			δxout[i]=-2.0*δx[i]*mod.ref.KI
		end
	elseif(parameterization == [:null, :χρ, :null])
		@inbounds for i in 1:nznx
			δxout[i]=-1.0*δx[i]*mod.ref.ρI
		end
	end
end

#pert_of_KI(dvp, dρ, vp0, ρ0) = - 2. / (ρ0 * vp0^3) * dvp - 1 / (ρ0^2 * vp0^2) * dρ
#pert_of_ρI(dρ, vp0, ρ0) = -1. / (ρ0^2) * dρ


text:
############# Définir les différentes mesures applicables ######################

module mesures
    # Pas d'intervention
    const INTERVENTION = true ;
    # Tracabilité
    const TRACAGE = true ;
    # Distanciation sociale
    const DISTANCE = false ;
    # Port du masque
    const MASQUE = false ;
    # Confinement
    const CONFINEMENT = false ;
    # Vaccination
    const VACCIN = false ;
    # Tele-travail
    const TELETRAVAIL = false ;
    # Fermeture
    const FERMETURE = false ;
    # Cas contact
    const CAS_CONTACT = false ;
end

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

############################## Tracabilite #####################################

# Module contenant les paramètres liés à la tracabilité

module Trace
    # Seuil d'alerte (nombre de nouvelles contamination)
    const SEUIL_INCIDENCE = 90 ;
    # Délai, en jours, pour isoler l'individu contaminé
    const DELAI_TRACE = 2 ;
    # % de personnes testés
    const TAUX_TEST = 0.5 ;
    # Nombre de jours confinés
    const JOURS_CONF = 14 ;
end

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

############################ Distanciation #####################################

# Module contenant les paramètres liés à la distanciation sociale

module Distanciation
    # Augmentation de la distance moyenne entre deux individus, exprimée en m
    const DISTANCE_SUP = 2 ;

end

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

############################ Port du masque ####################################

# Module contenant les paramètres liés au port du masque

module PortMasque
    # Augmentation de la distance moyenne entre deux individus, exprimée en m
    const DISTANCE_SUP = 2 ;

end

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

############################## Confinement #####################################

# Module contenant les paramètres liés au confinement

module Confinement
    # Seuil d'alerte
    const SEUIL_INCIDENCE = 300 ;
    # Durée du confinement, exprimée en jours
    const DUREE = 40 ;
    # Teletravail
    const TELETRAVAIL = true ;
    # Taux de télétravail (actif)
    const TAUX_TELE = 0.2 ;
    # fermeture commerce non-essentiel
    const POI_NE = true ;
    # fermeture ecole
    const FERMETURE_ECOLE = true ;
end

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

####################### Fermeture des écoles ###################################

# Module contenant les paramètres liés à la fermeture des écoles

module Fermeture
    # Seuil d'alerte
    const SEUIL_ALERTE = 100 ;
    # Ecoles
    const ECOLES = true ;
    # POI non-essentiel
    const POI_NE = true ;
    # POI essentiel
    const POI_E = false ;
end

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

####################### Fermeture des écoles ###################################

# Module contenant les paramètres liés au télé-travail

module Teletravail
    # Seuil d'alerte
    const SEUIL_ALERTE = 100 ;
    # taux d'actif en télé-travail (par entreprise ?)
    const TAUX_TELE = 0.9 ;
end

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

########################## Calcul des indicateurs ##############################

function Calcul_R0(nb_transmi,nb_contact_cont)

    #fin = length(nb_transmi) ;
    d = Dict([(i,count(x->x.==i,nb_contact_cont)) for i in unique(nb_contact_cont)]) ;
    dv = collect(keys(d)) ;

    R0_local = zeros(Int32,length(dv)) ;

    for i=1:length(dv)

        Ind = findall(nb_contact_cont .== dv[i]) ;
        R0_local[i] = sum(nb_transmi[Ind]) ;

    end

    return(R0_local) ;

end

#function Calcul_Reff()
#end

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

######################## Prise en compte des mesures ###########################

function MiseEnPlaceMesure(Ncontam)

    if(scenario == 0) # Si pas de mesures
        Indice_mesure = 0 ;
    else # Si on a mis en place une mesure au debut de la simulation
        if(Ncontam > 1000) # si on a plus de 1000 personnes contaminées, alors
            Indice_mesure = 1 ;
        end

    end

    return(Indice_mesure)

end

function RetourNormal(cpt)

    if(scenario == 4)
        cpt = cpt - 1 ;
    end

end

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


text:
using Test

# Build the example docs
include("examples/make.jl")

# Test missing docs
include("missingdocs/make.jl")

# Primary @testset

# Error reporting.
println("="^50)
@info("The following errors are expected output.")
include(joinpath("errors", "make.jl"))
@info("END of expected error output.")
println("="^50)

@testset "Documenter" begin
    # Unit tests for module internals.
    include("utilities.jl")

    # DocChecks tests
    include("docchecks.jl")

    # NavNode tests.
    include("navnode.jl")

    # DocSystem unit tests.
    include("docsystem.jl")

    # DocTest unit tests.
    include("doctests/doctests.jl")

    # DOM Tests.
    include("dom.jl")

    # MDFlatten tests.
    include("mdflatten.jl")

    # HTMLWriter
    include("htmlwriter.jl")

    # Mock package docs.
    include("examples/tests.jl")

    # Documenter package docs with other formats.
    include("formats/markdown.jl")

    # A simple build outside of a Git repository
    include("nongit/tests.jl")
end

# Additional tests

## `Markdown.MD` to `DOM.Node` conversion tests.
module MarkdownToNode
    import Documenter.DocSystem
    import Documenter.Writers.HTMLWriter: mdconvert

    # Exhaustive Conversion from Markdown to Nodes.
    for mod in Base.Docs.modules
        for (binding, multidoc) in DocSystem.getmeta(mod)
            for (typesig, docstr) in multidoc.docs
                md = DocSystem.parsedoc(docstr)
                string(mdconvert(md))
            end
        end
    end
end

# Docstring signature syntax highlighting tests.
module HighlightSig
    using Test
    import Markdown
    import Documenter.Expanders: highlightsig!

    @testset "highlightsig!" begin
        s = """
                foo(bar::Baz)
            ---
                foo(bar::Baz)
            """
        original = Markdown.parse(s)
        md = Markdown.parse(s)
        highlightsig!(md)
        @test isempty(original.content[1].language)
        @test md.content[1].language == "julia"
        @test original.content[end].language == md.content[end].language

        s = """
            ```lang
             foo(bar::Baz)
            ```
            """
        original = Markdown.parse(s)
        md = Markdown.parse(s)
        highlightsig!(md)
        @test original == md
    end
end


text:
STUDENTNAME = "JuanPabloOssaZapata"

using Logging
using XLSX
using Profile

include("psmsdst.jl")

@info "Compiling solvers..."

solvers = Dict(
    # Does not use any parameters
    include("solvers/random.jl")(attempts=10000),

    # Does not use any parameters
    include("solvers/greedy.jl")(),

    # -- PARAMETERS --
    # attempts: The number of solution attempts with no improvement until the
    # method stops. In other words, if the method is unable to get a better
    # solution in this many attempts, it will stop.
    # hardness: Determines how biased the method will be towards greedyness.
    # Mathematically, the hardness parameter controls the β parameter of a
    # softmax function.
    include("solvers/softgreedy.jl")(attempts=40,hardness=5.5)
)

function main()

    allproblems = loadproblems()

    for (solvername,solver) in solvers

        shortname = split(solvername,r"\s")[1]

        filename = "PSMSDST_$(STUDENTNAME)_$(shortname).xlsx"

        XLSX.openxlsx(filename,mode="w") do xf

            newfile = true

            for (problemname, problem) in allproblems

                sheetname = "$problemname"

                if newfile
                    newfile = false
                    sheet = xf[1]
                    XLSX.rename!(sheet,sheetname)
                else
                    sheet = XLSX.addsheet!(xf,sheetname)
                end

                @info "Solving problem" solvername problemname

                tstart = time()

                solution = solver(problem)

                tend = time()

                elapsed = tend - tstart

                feasible, reason = isfeasible(problem,solution)

                if !feasible
                    @error "SOLUTION PROVIDED BY SOLVER IS NOT FEASIBLE" #=
                    =# solvername problemname reason
                end

                starttimes, finishtimes, cost, relative =
                    evaluatesolution(problem,solution)

                optimal = problem.lowerbound == cost

                @info "Problem solved" solvername problemname elapsed #=
                =# feasible cost problem.lowerbound relative optimal

                if optimal
                    @warn "OPTIMAL SOLUTION FOUND!"
                end

                for (machineid,jobs) in enumerate(solution)
                    # First cell is the number of jobs
                    sheet[machineid,1] = length(jobs)
                    # Next cells are the job ids
                    for (index,job) in enumerate(jobs)
                        sheet[machineid,index+1] = job
                    end
                    # Next cells are the job start times
                    for (index,job) in enumerate(jobs)
                        sheet[machineid,index+1+length(jobs)] = job
                    end
                    # Final cell is the finish time
                    sheet[machineid,2+2*length(jobs)] = finishtimes[machineid]
                end
                # First cell of last row is the total cost
                sheet[length(solution)+1,1] = cost
            end
        end
    end
end

function loadproblems()::Vector{Tuple{String,PSMSDSTProblem}}
    result = Vector{Tuple{String,PSMSDSTProblem}}()
    @info "Loading problems..."
    for (root,_,files) in walkdir("PSMSDSTdata/")
        for file in files
            push!(result,(file,parsepsmsdst(joinpath(root,file))))
        end
    end

    reg = r"PSMSDST(\d+).txt"

    sort!(result,by=x->parse(Int64,match(reg,x[1])[1]))

    @info "Problems loaded"
    return result
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end

max_stars_repo_name: mkualquiera/huepl

max_issues_repo_name: mkualquiera/huepl

max_forks_repo_name: mkualquiera/huepl


text:
import Stimulus
import Spiketrains
import JLD
function processTemporalEntropy(sptrains::Dict, bins::Array{Float64,1}, trials::Union(Array{Stimulus.Trial,1},Array{Stimulus.SaccadeTrial,1}),ncombos::Integer, word_size::Integer;skipMissing::Bool=false,alignment::Symbol=:target)
	cells = Spiketrains.sortCells(sptrains)
	tmin = minimum(bins)
	tmax = maximum(bins)
	trial_labels = Stimulus.getTrialLocationLabel(trials,:target)
    D = pmap( _cell -> begin
            fname = string(_cell,"temporalEntropy.jd")
            if isfile(fname)
                println("File already exists. Skipping...")
			else
                if skipMissing
                    println("Skipping missing data..")
                else
                    isfile(fname) && println("No data found. Computing...")
                    S = sptrains[_cell]
                    println("Processing cell $_cell")
					spikes = Spiketrains.getTrialRaster(S,trials, alignment,tmin,tmax)
					psth = Spiketrains.getTrialSpikeCount(spikes, bins)
					GTE = Entropy.GroupedTemporalEntropy(int(psth.counts),bins,trial_labels[1:psth.ntrials],word_size)
					JLD.save(fname, {"GTE"=>GTE})
				end
			end
			fname
		end,cells)
	return D
end

@doc """
	Process the spike trains specified in the dictionary ``sptrains`` by first binning them using ``bins`` and aligning them according to ``trials``. The optional keyword ``alignment`` specifies which trial event to align the spike trains to (defaults to :target). ``word_size`` specified the number of bins to conatenate before hashign the responses. The return value is an array of filenames, each containing the estimated entropy for each cell.
""" [ :returns => Array{String,1} ] ->
function process{T<:GroupedShannonEntropy}(::Type{T},sptrains::Dict, bins::Array{Float64,1}, trials::Union(Array{Stimulus.Trial,1},Array{Stimulus.SaccadeTrial,1}),ncombos::Integer, word_size::Integer;skipMissing::Bool=false,alignment::Symbol=:target)
	cells = Spiketrains.sortCells(sptrains)
	tmin = minimum(bins)
	tmax = maximum(bins)
	trial_labels = Stimulus.getTrialLocationLabel(trials,:target)
	fnames = Array(ASCIIString,length(cells))
    D = pmap( _cell -> begin
            fname = string(_cell,filename(T))
            if isfile(fname)
                println("File already exists. Skipping...")
			else
                if skipMissing
                    println("Skipping missing data..")
                else
                    isfile(fname) && println("No data found. Computing...")
                    S = sptrains[_cell]
                    println("Processing cell $_cell")
					spikes = Spiketrains.getTrialRaster(S,trials, alignment,tmin,tmax)
					psth = Spiketrains.getTrialSpikeCount(spikes, bins)
					GTE = T(int(psth.counts),bins,trial_labels[1:psth.ntrials],word_size)
					JLD.save(fname, {varname(T)=>GTE})
				end
			end
			fname
		end,cells)
		for (i,d) in enumerate(D)
			fnames[i] = d
		end
	return fnames
end


text:
module JSONRPC

import JSON, UUIDs

include("packagedef.jl")

end


text:
#=
P35 (**) Determine the prime factors of a given positive integer.
Construct a flat list containing the prime factors in ascending order.
Example:
* (prime-factors 315)
(3 3 5 7)
=#

SMALLEST_PRIME = 2

function next_prime(n)
    # Brute-force incremental check
    isprime(n) ? n : next_prime(n + 1)
end

function _prime_factors(n, x, pf)
    if isprime(n)
        vcat(pf, n)
    else
        if mod(n, x) == 0
            _prime_factors(div(n, x), SMALLEST_PRIME, vcat(pf, x))
        else
            _prime_factors(n, next_prime(x + 1), pf)
        end
    end
end

function prime_factors(n)
    # Starting with the smallest prime, 2
    isprime(n) ? n : _prime_factors(n, SMALLEST_PRIME, [])
end

@assert prime_factors(315) == [3, 3, 5, 7]
@assert prime_factors(360) == [2, 2, 2, 3, 3, 5]
@assert prime_factors(144) == [2, 2, 2, 2, 3, 3]

println("Tests passed: JL-35.jl")


text:
@testset "Storage data misspecification" begin
    # See https://discourse.julialang.org/t/how-to-use-test-warn/15557/5 about testing for warning throwing
    info_message = "The data doesn't include devices of type GenericBattery, consider changing the device models"
    device_model = DeviceModel(GenericBattery, BookKeeping)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5)
    @test_logs (:info, info_message) match_mode = :any mock_construct_device!(
        model,
        device_model,
    )
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys14)
    @test_logs (:info, info_message) match_mode = :any mock_construct_device!(
        model,
        device_model,
    )
end

@testset "Storage Basic Storage With DC - PF" begin
    device_model = DeviceModel(
        GenericBattery,
        BookKeeping;
        attributes=Dict{String, Any}("reservation" => false),
    )
    c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_bat)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 72, 0, 72, 72, 24, false)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Storage Basic Storage With AC - PF" begin
    device_model = DeviceModel(
        GenericBattery,
        BookKeeping;
        attributes=Dict{String, Any}("reservation" => false),
    )
    c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat")
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5_bat)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 96, 0, 96, 96, 24, false)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Storage with Reservation  & DC - PF" begin
    device_model = DeviceModel(GenericBattery, BookKeeping)
    c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_bat)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 96, 0, 72, 72, 24, true)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Storage with Reservation  & AC - PF" begin
    device_model = DeviceModel(GenericBattery, BookKeeping)
    c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat")
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5_bat)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 120, 0, 96, 96, 24, true)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Storage with BatteryAncillaryServices and Reservation DC - PF" begin
    device_model = DeviceModel(GenericBattery, BatteryAncillaryServices)
    c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_bat)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 96, 0, 72, 72, 24, true)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Storage with BatteryAncillaryServices and Reservation With AC - PF" begin
    device_model = DeviceModel(GenericBattery, BatteryAncillaryServices)
    c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat")
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5_bat)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 120, 0, 96, 96, 24, true)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "BatteryEMS with EnergyTarget with DC - PF" begin
    device_model = DeviceModel(BatteryEMS, EnergyTarget)
    c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat_ems")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_bat)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 144, 0, 72, 72, 48, true)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "BatteryEMS with EnergyTarget With AC - PF" begin
    device_model = DeviceModel(BatteryEMS, EnergyTarget)
    c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat_ems")
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5_bat)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 168, 0, 96, 96, 48, true)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "BatteryEMS with EnergyTarget Formulations (energy target - cases 1b-2b" begin
    template = ProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, BatteryEMS, EnergyTarget)
    set_device_model!(template, RenewableDispatch, RenewableFullDispatch)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    c_sys5 = PSB.build_system(PSITestSystems, "batt_test_case_b_sys")
    model =
        DecisionModel(EconomicDispatchProblem, template, c_sys5; optimizer=HiGHS_optimizer)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    moi_tests(model, false, 21, 0, 12, 9, 9, true)
    psi_checksolve_test(model, [MOI.OPTIMAL], 5811.0, 10.0)
end

@testset "BatteryEMS with EnergyTarget Formulations (energy target - cases 1c-2c" begin
    template = ProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, BatteryEMS, EnergyTarget)
    set_device_model!(template, RenewableDispatch, RenewableFullDispatch)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    c_sys5 = PSB.build_system(PSITestSystems, "batt_test_case_c_sys")

    model =
        DecisionModel(EconomicDispatchProblem, template, c_sys5; optimizer=HiGHS_optimizer)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    moi_tests(model, false, 21, 0, 12, 9, 9, true)
    psi_checksolve_test(model, [MOI.OPTIMAL], -63.0, 10.0)
end

@testset "BatteryEMS with EnergyTarget Formulations (energy target - cases 1d-2d" begin
    template = ProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, BatteryEMS, EnergyTarget)
    set_device_model!(template, RenewableDispatch, RenewableFullDispatch)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    c_sys5 = PSB.build_system(PSITestSystems, "batt_test_case_d_sys")

    model =
        DecisionModel(EconomicDispatchProblem, template, c_sys5; optimizer=HiGHS_optimizer)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    moi_tests(model, false, 28, 0, 16, 12, 12, true)
    psi_checksolve_test(model, [MOI.OPTIMAL], -11118.0, 10.0)
end

@testset "BatteryEMS with EnergyTarget Formulations (energy target - cases 1e-2e" begin
    template = ProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, BatteryEMS, EnergyTarget)
    set_device_model!(template, RenewableDispatch, RenewableFullDispatch)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    c_sys5 = PSB.build_system(PSITestSystems, "batt_test_case_e_sys")
    model =
        DecisionModel(EconomicDispatchProblem, template, c_sys5; optimizer=HiGHS_optimizer)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    moi_tests(model, false, 21, 0, 12, 9, 9, true)
    psi_checksolve_test(model, [MOI.OPTIMAL], 5547.0, 10.0)
end

@testset "BatteryEMS with EnergyTarget Formulations (energy target - cases 1f-2f" begin
    template = ProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, BatteryEMS, EnergyTarget)
    set_device_model!(template, RenewableDispatch, RenewableFullDispatch)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    c_sys5 = PSB.build_system(PSITestSystems, "batt_test_case_f_sys")

    model =
        DecisionModel(EconomicDispatchProblem, template, c_sys5; optimizer=HiGHS_optimizer)
    @test build!(model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    moi_tests(model, false, 21, 0, 12, 9, 9, true)
    psi_checksolve_test(model, [MOI.OPTIMAL], -1825.0, 10.0)
end

### Feedforward Test ###

@testset "Test EnergyTargetFeedforward to GenericBattery with BookKeeping model" begin
    device_model = DeviceModel(GenericBattery, BookKeeping)

    ff_et = EnergyTargetFeedforward(
        component_type=GenericBattery,
        source=EnergyVariable,
        affected_values=[EnergyVariable],
        target_period=12,
        penalty_cost=1e5,
    )

    PSI.attach_feedforward!(device_model, ff_et)
    sys = PSB.build_system(PSITestSystems, "c_sys5_bat")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, sys)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 120, 0, 72, 73, 24, true)
end

@testset "Test EnergyLimitFeedforward to GenericBattery with BookKeeping model" begin
    device_model = DeviceModel(GenericBattery, BookKeeping)

    ff_il = EnergyLimitFeedforward(
        component_type=GenericBattery,
        source=ActivePowerOutVariable,
        affected_values=[ActivePowerOutVariable],
        number_of_periods=12,
    )

    PSI.attach_feedforward!(device_model, ff_il)
    sys = PSB.build_system(PSITestSystems, "c_sys5_bat")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, sys)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 96, 0, 74, 72, 24, true)
end

@testset "Test EnergyTargetFeedforward to GenericBattery with BookKeeping model" begin
    device_model = DeviceModel(GenericBattery, BatteryAncillaryServices)

    ff_et = EnergyTargetFeedforward(
        component_type=GenericBattery,
        source=EnergyVariable,
        affected_values=[EnergyVariable],
        target_period=12,
        penalty_cost=1e5,
    )

    PSI.attach_feedforward!(device_model, ff_et)
    sys = PSB.build_system(PSITestSystems, "c_sys5_bat")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, sys)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 120, 0, 72, 73, 24, true)
end

@testset "Test EnergyLimitFeedforward to GenericBattery with BatteryAncillaryServices model" begin
    device_model = DeviceModel(GenericBattery, BatteryAncillaryServices)

    ff_il = EnergyLimitFeedforward(
        component_type=GenericBattery,
        source=ActivePowerOutVariable,
        affected_values=[ActivePowerOutVariable],
        number_of_periods=12,
    )

    PSI.attach_feedforward!(device_model, ff_il)
    sys = PSB.build_system(PSITestSystems, "c_sys5_bat")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, sys)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 96, 0, 74, 72, 24, true)
end

@testset "Test EnergyTargetFeedforward to GenericBattery with BookKeeping model" begin
    device_model = DeviceModel(BatteryEMS, BookKeeping)

    ff_et = EnergyTargetFeedforward(
        component_type=BatteryEMS,
        source=EnergyVariable,
        affected_values=[EnergyVariable],
        target_period=12,
        penalty_cost=1e5,
    )

    PSI.attach_feedforward!(device_model, ff_et)
    sys = PSB.build_system(PSITestSystems, "c_sys5_bat_ems")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, sys)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 120, 0, 72, 73, 24, true)
end

@testset "Test EnergyLimitFeedforward to BatteryEMS with BookKeeping model" begin
    device_model = DeviceModel(BatteryEMS, BookKeeping)

    ff_il = EnergyLimitFeedforward(
        component_type=BatteryEMS,
        source=ActivePowerOutVariable,
        affected_values=[ActivePowerOutVariable],
        number_of_periods=12,
    )

    PSI.attach_feedforward!(device_model, ff_il)
    sys = PSB.build_system(PSITestSystems, "c_sys5_bat_ems")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, sys)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 96, 0, 74, 72, 24, true)
end

@testset "Test EnergyTargetFeedforward to GenericBattery with BatteryAncillaryServices model" begin
    device_model = DeviceModel(BatteryEMS, BatteryAncillaryServices)

    ff_et = EnergyTargetFeedforward(
        component_type=BatteryEMS,
        source=EnergyVariable,
        affected_values=[EnergyVariable],
        target_period=12,
        penalty_cost=1e5,
    )

    PSI.attach_feedforward!(device_model, ff_et)
    sys = PSB.build_system(PSITestSystems, "c_sys5_bat_ems")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, sys)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 120, 0, 72, 73, 24, true)
end

@testset "Test EnergyLimitFeedforward to BatteryEMS with BatteryAncillaryServices model" begin
    device_model = DeviceModel(BatteryEMS, BatteryAncillaryServices)

    ff_il = EnergyLimitFeedforward(
        component_type=BatteryEMS,
        source=ActivePowerOutVariable,
        affected_values=[ActivePowerOutVariable],
        number_of_periods=12,
    )

    PSI.attach_feedforward!(device_model, ff_il)
    sys = PSB.build_system(PSITestSystems, "c_sys5_bat_ems")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, sys)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 96, 0, 74, 72, 24, true)
end


text:
module SNNPlot
using Statistics, Plots, GraphPlot

export rasterplot, voltageplot, smooth, rateplot, rArBplot, rArBgif

"""
    rasterplot(S::AbstractMatrix{<:Real}; timestep=1, kwargs...)

Generate a raster plot of neuronal spiking given the matrix `S[neuron, time]`
of neuronal spike trains, where each time-step in `C` corresponds
to `timestep=1` seconds.
"""
function rasterplot(S::AbstractMatrix{Bool}; timestep=1, kwargs...)
    spikes = findall(S)
    N = getindex.(spikes, 1)
    T = getindex.(spikes, 2) * timestep
    return scatter(T, N; leg=false, marker=:vline, ylab="Neuron", kwargs...)
end

"""
    voltageplot(V::AbstractArray{<:Real}; timestep=1, kwargs...)

Generate a plot of neuronal voltage over time from the matrix V[neuron, time],
where each column of V is separated by `timestep=1` seconds.
"""
function voltageplot(V::AbstractArray{<:Real}; timestep=1, kwargs...)
    return plot(range(0, length=size(V, 2), step=timestep), V'; ylab="Voltage (V)", leg=false, kwargs...)
end

"""
    smooth(S::AbstractArray{<:Real}; window=0, dt=1, timestep=1, average=true)

    Calculate the average of S for all rows (if `average==true`) or each
individual row (if `average==false`) over a window of size `window`
which slides in strides of `dt` columns, where each column corresponds
to `timestep` seconds. Multiplies the values by `tohz`, which by default converts
boolean spikes into firing rates.  Return the tuple of `smoothed[index, neuron]`
and `time[index]` where `index` is the number of strides, `smoothed`
contains the smoothed value at stride `index`, and `time` contains the
continuous timepoint of that stride in seconds.
"""
function smooth(S::AbstractArray{<:Real}; window=0, dt=1, timestep=1, average=true, tohz=1/(timestep*dt))
    if ndims(S) == 1
        S = reshape(S, 1, length(S))  # Reshape to 1 Neuron x length(S) timesteps
        average = false
    end

    N, T = size(S)
    bins = floor(Int, (T - window) / dt)
    if bins <= 0
        return [], []
    end

    R = Matrix{Float64}(undef, bins, average ? 1 : N)
    t = range(window*timestep/2, length=bins, step=timestep*dt)

    @inbounds for t in 0:bins-1
        if average
            R[t+1, 1] = mean(S[:,(dt*t + 1):(dt*t + window + 1)]) * tohz
        else
            for n in 1:N
                R[t+1, n] = mean(S[n, (dt*t + 1):(dt*t + window + 1)]) * tohz
            end
        end
    end

    return R, t
end

"""
    rateplot(S::AbstractArray{<:Real}; window=0, dt=1, timestep=1, average=true, tohz=1/(timestep*dt), kwargs...)

Generate a population rate plot of neuronal spiking given the matrix
`S[neuron, time]` of neuronal spike trains. Rates are averaged
over a window of `window=1` time-steps which slides in strides of
 `dt=1` time-steps, where a single time-step corresponds to `timestep=1` seconds.
 If `average` is `true`, then average over all neurons in each window. Otherwise,
 plot a separate line for each neuron.
"""
function rateplot(S::AbstractArray{<:Real}; window=0, dt=1, timestep=1, average=true, tohz=1/(timestep*dt), kwargs...)
    R, t = smooth(S; window=window, dt=dt, timestep=timestep, average=average, tohz=tohz)
    return plot(t, R; ylab="Firing Rate (Hz)", leg=false, kwargs...)
end

"""
    rArBplot(SA::AbstractArray{<:Real}, SB::AbstractArray{<:Real};
             window=0, dt=1, timestep=1, average=true, tohz=1/(timestep*dt), kwargs...)

Generate a plot of the rate of two neural populations over time, where
each axis corresponds to the rate of a population, and the line is a traversal
through this rate space over time.
"""
function rArBplot(SA::AbstractArray{<:Real}, SB::AbstractArray{<:Real}; window=0, dt=1, timestep=1, average=true, tohz=1/(timestep*dt), kwargs...)
    rA, _ = smooth(SA; window=window, dt=dt, timestep=timestep, average=average, tohz=tohz)
    rB, _ = smooth(SB; window=window, dt=dt, timestep=timestep, average=average, tohz=tohz)
    return plot(rA, rB; leg=false, kwargs...)
end

"""
    rArBgif(SA::AbstractArray{<:Real}, SB::AbstractArray{<:Real};
            window=0, dt=1, timestep=1, average=true, tohz=1/(timestep*dt), fps=1, tail=0, skip=1, kwargs...)

Generate a gif of the rate of two neural populations over time, where
each axis corresponds to the rate of a population, and the line is a traversal
through this rate space over time.
"""
function rArBgif(SA::AbstractArray{<:Real}, SB::AbstractArray{<:Real};
                 window=0, dt=1, timestep=1, average=true, tohz=1/(timestep*dt), fps=1, tail=0, skip=1, filename=nothing, kwargs...)
    rA, _ = smooth(SA; window=window, dt=dt, timestep=timestep, average=average, tohz=tohz)
    rB, _ = smooth(SB; window=window, dt=dt, timestep=timestep, average=average, tohz=tohz)

    anim = @animate for i in 1:size(rA, 1)
        scatter(rA[i,:], rB[i,:]; kwargs...)
        if !iszero(tail)
            range = max(1, i-tail):i
            plot!(rA[range,:], rB[range,:]; leg=false, kwargs...)
        end
    end every skip
    if isnothing(filename)
        return gif(anim; fps=fps)
    end
    return gif(anim, filename; fps=fps)
end

function snnplot(snn; kwargs...)
    gplot(snn.graph; kwargs...)
end

end


text:
using Distributions, ForwardDiff
using Random: seed!

seed!(1)
n = 1000
d = 10
X = rand(Normal(), n, d)
b = rand(Normal(), d)
ϵ = rand(Normal(0.0, 0.1), n)
Y = X * b + ϵ
obj(Y, X, b) = mean((Y - X * b) .^ 2)

epochs = 100

# momentum
θ_momentum = rand(Normal(), d)
opt = Momentum(η=1.0)

for i in 1:epochs
    g = ForwardDiff.gradient(θ -> obj(Y, X, θ), θ_momentum)

    δ = update(opt, g)
    global θ_momentum -= δ
end

# adagrad
θ_adagrad = rand(Normal(), d)
opt = Adagrad(η=1.0)

for i in 1:epochs
    g = ForwardDiff.gradient(θ -> obj(Y, X, θ), θ_adagrad)

    δ = update(opt, g)
    global θ_adagrad -= δ
end

# rmsprop
θ_rmsprop = rand(Normal(), d)
opt = RMSprop(η=0.1)

for i in 1:epochs
    g = ForwardDiff.gradient(θ -> obj(Y, X, θ), θ_rmsprop)

    δ = update(opt, g)
    global θ_rmsprop -= δ
end

# adam
θ_adam = rand(Normal(), d)
opt = Adam(α=1.0)

for i in 1:epochs
    g = ForwardDiff.gradient(θ -> obj(Y, X, θ), θ_adam)

    δ = update(opt, g)
    global θ_adam -= δ
end

# adamax
θ_adamax = rand(Normal(), d)
opt = Adamax(α=1.0)

for i in 1:epochs
    g = ForwardDiff.gradient(θ -> obj(Y, X, θ), θ_adamax)

    δ = update(opt, g)
    global θ_adamax -= δ
end

# nadam
θ_nadam = rand(Normal(), d)
opt = Nadam(η=1.0)

for i in 1:epochs
    g = ForwardDiff.gradient(θ -> obj(Y, X, θ), θ_nadam)

    δ = update(opt, g)
    global θ_nadam -= δ
end

## vanilla
θ_vanilla = rand(Normal(), d)
opt = VanillaGradDescent(η=0.001)

for i in 1:20*epochs
    g = ForwardDiff.gradient(θ -> obj(Y, X, θ), θ_vanilla)

    δ = update(opt, g)
    global θ_vanilla -= δ
end

## inversedecay
θ_invdec = rand(Normal(), d)
opt = Inversedecay(κ=0.9)

for i in 1:epochs
    g = ForwardDiff.gradient(θ -> obj(Y, X, θ), θ_invdec)

    δ = update(opt, g)
    global θ_invdec -= δ
end

##

@testset "Linear Regression" begin
    @test b ≈ θ_momentum atol=1e-1
    @test b ≈ θ_adagrad atol=1e-1
    @test b ≈ θ_adam atol=1e-1
    @test b ≈ θ_adamax atol=1e-1
    @test b ≈ θ_nadam atol=1e-1
    @test b ≈ θ_vanilla atol=1e-1
    @test b ≈ θ_invdec atol=1e-1
end


text:
using Test
using RedisGraph: Graph, Node, Edge, Path, addnode!, addedge!, commit, delete, query, getdatabase

function creategraph()
    db_conn = getdatabase()
    g = Graph("TestGraph", db_conn)
    return g
end

function simplerelation!(g::Graph)
    node1 = Node("FirstSimpleNode", Dict("IntProp" => 1, "StringProp" => "node prop", "BoolProp" => true))
    node2 = Node("SecondSimpleNode")
    edge = Edge("SimpleEdge", node1, node2, Dict("IntProp" => 1, "StringProp" => "node prop", "BoolProp" => false))

    addnode!(g, node1)
    addnode!(g, node2)
    addedge!(g, edge)
    res = commit(g)
end

function deletegraph!(g::Graph)
    delete(g)
end

@testset "RedisGraph.jl tests" begin
    @testset "check simple types" begin
        g = creategraph()
        try
            @test query(g, "RETURN null").results[1] === nothing
            @test query(g, "RETURN 2").results[1] == 2
            @test query(g, "RETURN 2.0").results[1] == 2.0
            @test query(g, "RETURN true").results[1] == true
            @test query(g, "RETURN [1, null, 'test', 3.0, false]").results[1] == [1, nothing, "test", 3.0, false]
            @test query(g, "RETURN {a: 1, b: null, c: 'test', d: 3.0, e: false}").results[1] == Dict("a" => 1, "b" => nothing, "c" => "test", "d" => 3.0, "e" => false)
        finally
            deletegraph!(g)
        end
    end
    @testset "check point data type" begin
        g = creategraph()
        try
            # this set is from https://github.com/RedisGraph/redisgraph-py
            expected_lat = 32.070794860
            expected_lon = 34.820751118
            actual = query(g, "RETURN point({latitude: 32.070794860, longitude: 34.820751118})").results[1]
            @test abs(actual["latitude"] - expected_lat) < 0.001
            @test abs(actual["longitude"] - expected_lon) < 0.001

            expected_lat = 32
            expected_lon = 34
            actual = query(g, "RETURN point({latitude: 32, longitude: 34.0})").results[1]
            @test abs(actual["latitude"] - expected_lat) < 0.001
            @test abs(actual["longitude"] - expected_lon) < 0.001
        finally
            deletegraph!(g)
        end
    end
    @testset "check simple relation" begin
        g = creategraph()
        try
            simplerelation!(g)

            q = query(g, "MATCH (n1)-[e]->(n2) RETURN n1, e, n2").results
            node1, edge, node2 = q[1:3]

            @test typeof(node1) == Node
            @test node1.label == "FirstSimpleNode"
            @test node1.properties == Dict("IntProp" => 1, "StringProp" => "node prop", "BoolProp" => true)
            @test typeof(edge) == Edge
            @test edge.relation == "SimpleEdge"
            @test edge.properties == Dict("IntProp" => 1, "StringProp" => "node prop", "BoolProp" => false)
            @test edge.src_node == node1
            @test edge.dest_node == node2
            @test typeof(node2) == Node
            @test node2.label == "SecondSimpleNode"
        finally
            deletegraph!(g)
        end
    end
    @testset "check path type" begin
        g = creategraph()
        try
            node1 = Node("FirstSimpleNode", Dict("IntProp" => 1, "StringProp" => "node prop", "BoolProp" => true))
            node2 = Node("SecondSimpleNode")
            edge = Edge("SimpleEdge", node1, node2, Dict("IntProp" => 1, "StringProp" => "node prop", "BoolProp" => false))

            addnode!(g, node1)
            addnode!(g, node2)
            addedge!(g, edge)
            res = commit(g)

            result_path = query(g, "MATCH p=(n1)-[e]->(n2) RETURN p").results[1]

            @test typeof(result_path) == Path
            @test result_path.nodes == Path([node1, node2], [edge]).nodes
            # it's not clear how to check edges in this case
            # @test result_path.edges == Path([node1, node2], [edge]).edges
        finally
            deletegraph!(g)
        end
    end
end


text: println("Hello, Julia!")

max_stars_repo_name: JuliaTagBot/JuliaBook-Samples

max_issues_repo_name: JuliaTagBot/JuliaBook-Samples

max_forks_repo_name: JuliaTagBot/JuliaBook-Samples


text:
using BinaryProvider

dependencies = [
    ("gmp", "6.1.2", "https://github.com/JuliaMath/GMPBuilder/releases/download/v6.1.2-2/build_GMP.v6.1.2.jl",
        "91b2197c4d2209985b2c156e04f3fc11e4beff633cd86bb0053bfefe71cdfba6"),

    ("mpfr", "4.0.1", "https://github.com/JuliaMath/MPFRBuilder/releases/download/v4.0.1-3/build_MPFR.v4.0.1.jl",
        "40f3d06abf45b571b23135a079d66827282e587be15d540874bb45dac163e2c7"),

    ("mpc", "1.1.0", "https://github.com/isuruf/MPCBuilder/releases/download/v1.1.0-2/build_MPC.v1.1.0.jl",
        "85cac0057832da9c9d965531e9a1bada7150032aea4dbead59ff76e95bbdc47f"),

    ("symengine", "0.3.0", "https://github.com/symengine/SymEngineBuilder/releases/download/v0.3.0-3/build_SymEngine.v0.3.0.jl",
        "c6122fbb9ef8198f413c645e66a251593d50a19414528d21f14d4a53ca5e299d"),
]

libdir = "lib"

if Sys.iswindows()
    libdir = "bin"
    dependencies[4]= ("symengine", "0.3.0", "https://github.com/symengine/SymEngineBuilder/releases/download/v0.3.0-2/build_SymEngine.v0.3.0.jl",
        "664b7df2b2e173625fa5742aa194e63392692489b73e6ba4005dcbf661093c9d")
end

prefix = joinpath(@__DIR__, "libsymengine-0.3")
downloads_dir = joinpath(@__DIR__, "downloads")

all_products = LibraryProduct[]

if !isdir(joinpath(prefix, libdir))
    mkpath(joinpath(prefix, libdir))
end

for (name, version, url, hash) in dependencies
    product = LibraryProduct(joinpath(prefix, libdir), String[string("lib", name)], Symbol(string("lib", name)))
    push!(all_products, product)

    if !satisfied(product)
        tmp_file = joinpath(downloads_dir, "build_$hash.jl")
        download_verify(url, hash, tmp_file)
        contents = read(tmp_file, String)
        tmp_prefix = joinpath(downloads_dir, string(name, "-", version))
        m = Module(:__anon__)
        Core.eval(m, quote
            using BinaryProvider
            function write_deps_file(path, products; verbose=true) end
            ARGS = [$tmp_prefix]
        end)
        Base.include_string(m, contents)
        for file in readdir(joinpath(tmp_prefix, libdir))
            cp(joinpath(tmp_prefix, libdir, file), joinpath(prefix, libdir, file), force=true)
        end
    end
    @info(locate(product))
end

for product in all_products
    locate(product, verbose=true)
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), all_products)


text:
module Logic

VAR_COUNTER = 0

struct Var{T}
    name::String
end

function Var{T}() where T
    global VAR_COUNTER

    VAR_COUNTER += 1
    Var{T}("_VAR$VAR_COUNTER")
end

Var(args...) = Var{Any}(args...)

Variable{T} = Union{T, Var{T}} where T

macro var(decl)
    name, type = decl isa Symbol ? (decl, Any) : decl.args
    quote
        $(esc(name)) = Var{$type}($(string(name)))
    end
end

abstract type Predicate end

macro predicate(decl)
    name, args = decl.args[1], decl.args[2:end]
    arg_names = [field isa Symbol ? field : field.args[1] for field in args]
    quote
        struct $(esc(name)) <: Predicate
            $(map(esc, args)...)
        end

        # params(p::$(esc(name))) = tuple($([:(p.$arg_name) for arg_name in arg_names]...)) #tuple((p.$(esc(arg_name)) for arg_name in arg_names)...)
    end
end

@generated function params(p::P) where P <: Predicate
    fields = fieldnames(P)
    :(
        tuple($([:(p.$name) for name in fields]...))
    )
end

nameof(::T) where T<:Predicate  = Symbol(T)

struct Conj{T<:NTuple{N, Predicate} where N}
    predicates::T
end
Conj(args...) = Conj(args)

struct Disj{T<:NTuple{N, Predicate} where N}
    predicates::T
end
Disj(args...) = Disj(args)

struct Clause
    head::Predicate
    body::Union{Predicate, Conj, Disj, Nothing}
end

Clause(head::Predicate) = Clause(head, nothing)

abstract type Result end

struct Success<:Result
    substitutions::Dict{Var, Any}
end

Success() = Success(Dict{Var, Any}())

Base.:(==)(s1::Success, s2::Success) = s1.substitutions == s2.substitutions

struct Failure<:Result
    val1
    val2
end

unify(v::Var{T}, x::Variable{T}) where T = Success(Dict(v => x))
unify(x::T, v::Var{T}) where T = unify(v, x)

unify(p1::P, p2::P) where P<:Predicate = intersection(unify(p1, p2) for (p1, p2) in zip(params(p1), params(p2)))

unify(p::Predicate, cl::Clause) = unify(p, cl.head)

# fallback
unify(x, y) = x == y ? Success() : Failure(x, y)

intersection(unifications) = reduce(intersection, unifications)
function intersection(unification1::Success, unification2::Success)
    u1, u2 = unification1.substitutions, unification2.substitutions
    merge_dict = Dict()
    common_vars = intersect(keys(u1), keys(u2))
    for var in common_vars
        v1, v2 = u1[var], u2[var]
        v, possible = unionof(v1, v2)
        !possible && return Failure(v1, v2)
        merge_dict[var] = v
    end
    Success(merge(u1, u2, merge_dict))
end
intersection(u1::Failure, u2) = u1
intersection(u1, u2::Failure) = u2

unionof(v1::Var{T}, v2::Variable{T}) where T = v2, true
unionof(v1::T, v2::Var{T}) where T = v1, true
unionof(v1, v2) = v1 == v2 ? (v1, true) : (nothing, false)

function query(p::Predicate, db::AbstractArray{Clause})
    matches = Iterators.filter(result -> results isa Success, ((clause, unify(p, clause)) for clause in db))
    isempty(matches) && return false
    for (clause, match) in matches
        result = reduce(intersection, (query(clause.body, db)))
    end
end

end


text:
using GLVisualize, GeometryTypes, Reactive, GLAbstraction

if !isdefined(:runtests)
    window = glscreen()
    timesignal = bounce(linspace(0,1,360))
end

function spiral(i, start_radius, offset)
	Point2f0(sin(i), cos(i)) * (start_radius + ((i/2pi)*offset))
end
# 2D particles
particle_data2D(i, N) = Point2f0[spiral(i+x, 3, 10) for x=1:N]
# stretch time a bit:
t = const_lift(*, timesignal, 30f0)

# the simplest of all, plain 2D particles.
# to make it a little more interesting, we animate the particles a bit!
particles = const_lift(particle_data2D, t, 256)

# create a visualisation with each particle being 15px wide
# if you omit the primitive, it defaults to a SimpleRectangle
vis = visualize(particles, scale=Vec2f0(15))
_view(vis, window, camera=:orthographic_pixel)

if !isdefined(:runtests)
    renderloop(window)
end


text:
import GeoStatsBase: fit, predict, status

# ----------------
# DUMMY ESTIMATOR
# ----------------

struct DummyEstimator end
struct FittedDummyEstimator end

fit(::DummyEstimator, data) = FittedDummyEstimator()
predict(::FittedDummyEstimator, var, pₒ) = (0, 1)
status(::FittedDummyEstimator) = true

# ------------------------
# DUMMY SIMULATION SOLVER
# ------------------------

import GeoStatsBase: solvesingle

@simsolver DummySimSolver begin end
function solvesingle(problem::SimulationProblem, covars::NamedTuple,
                     ::DummySimSolver, preproc)
  npts = nelements(domain(problem))
  mactypeof = Dict(name(v) => mactype(v) for v in variables(problem))
  reals = map(covars.names) do var
    V    = mactypeof[var]
    real = vcat(fill(zero(V), npts÷2), fill(one(V), npts÷2))
    var => real
  end
  Dict(reals)
end


text:
module MESolve

using DifferentialEquations, LinearAlgebra

# export dRho, dRho_vec, dCV, d_av
export dRho!, dCV!, d_av!
export me_solve_time_independent, me_solve_H_time_dependent,
								  me_solve_L_time_dependent,
								  me_solve_full_time_dependent
								  # me_solve_time_independent_vec

export CV_solve_time_independent, CV_solve_time_dependent
export create_HO, create_HO_network, create_HO_network!,
									 create_KPO_network,
									 create_KPO_network!,
									 create_Qubit,
									 create_Qubit_network,
									 create_Qubit_network!,
									 husimi_Q,
									 husimi_Q_pw

include("DerivativeFunctions.jl")
include("SolverFunctions.jl")
include("StateConstructorFunctions.jl")
include("CV_QI_Functions.jl")
include("CV_SolverFunctions.jl")

end


text:
using IterTools
using LinearAlgebra

function newton_update_theta(eta, eta_beta, theta)
    tensor_shape = size(eta)
    d = ndims(eta)
    idxs = []

    vecs_theta = []
    vecs_eta = []
    vecs_eta_beta = []
    for k = 1:d
        for i = 2:tensor_shape[k]
            idx = ntuple( j->j == k ? i : 1, d)
            push!(idxs, idx)
            push!(vecs_theta, theta[idx...])
            push!(vecs_eta, eta[idx...])
            push!(vecs_eta_beta, eta_beta[idx...])
        end
    end

    n_idxs = length(idxs)
    G = zeros(n_idxs,n_idxs)
    for gi = 1:n_idxs
        for gj = 1:n_idxs
            ix = ntuple( j->max(idxs[gi][j], idxs[gj][j]), d)
            G[gi,gj] += eta[ix...]
            G[gi,gj] -= eta[ idxs[gi]... ] * eta[ idxs[gj]...]
        end
    end

    # newton_update
    inv_G = inv(G)
    vecs_theta = vecs_theta - inv_G * (vecs_eta - vecs_eta_beta)

    for gi = 1:n_idxs
        theta[ idxs[gi]... ] = vecs_theta[gi]
    end

    return theta
end

max_stars_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_issues_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_forks_repo_name: gkazunii/Legendre-tucker-rank-reduction


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

"""
    Beam2D

Euler-Bernoulli beam for 2d problems.
"""
struct Beam2D <: FieldProblem end

FEMBase.get_unknown_field_name(::Problem{Beam2D}) = "displacement"

function FEMBase.assemble_elements!(::Problem{Beam2D}, ::Assembly,
                               elements::Vector{Element{B}}, ::Float64) where B

    for element in elements
        @info("Not doing anything useful right now (someone should implement).")
    end

    return nothing
end


text:
using HDF5, JLD2

function jld2h5(JLDfileName::String)
    checkExtensionIsJLD(JLDfileName)

    dict = load(JLDfileName)
    h5IO = createIOStreamForH5(JLDfileName)
    writeJLDDataToH5(h5IO,dict)
    close(h5IO)

    nothing
end

function checkExtensionIsJLD(fileName::String)
    nCharsInFile = length(fileName)
    fileExtension = fileName[nCharsInFile-4 : nCharsInFile]
    if fileExtension != ".jld2"
        error("input extension was '$(fileExtension)', not '.jld2'")
    end
end

function getFileBase(fileName::String)
    nCharsInFile = length(fileName)
    fileBase = fileName[1 : nCharsInFile-5]
    return fileBase
end

function createIOStreamForH5(fileName::String)
    fileBase = getFileBase(fileName)
    h5FileName = string(fileBase,".h5")
    h5IO = h5open(h5FileName, "w")
    return h5IO
end

function writeJLDDataToH5(h5IO,dict)
    for key in keys(dict)
        checkAndModifyDataIfSpecialCase(dict,key)
        write(h5IO, string(key), dict[key])
    end
end

function valueIsStepRangeLen(value)
    if typeof(value) == StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
        return true
    else
        return false
    end
end

function modifyForStepRangeLen!(dict,key)
    println("StepRangeLen converted to Vector")
    dict[key] = collect(dict[key])
    nothing
end

function valueIsArrayOfArrays(value)
    if typeof(value) == Array{Array{Float64,1},1}
        return true
    else
        return false
    end
end

function modifyForArrayOfArrays!(dict,key)
    n1 = length(dict[key])
    n2 = length(dict[key][1])
    newArr = zeros(Float64,n1,n2)
    for i=1:n1, j=1:n2
        newArr[i,j] = dict[key][i][j]
    end
    println("Array of Array converted to matrix")
    dict[key] = newArr
    nothing
end

function checkAndModifyDataIfSpecialCase(dict,key)
    if valueIsStepRangeLen(dict[key])
        modifyForStepRangeLen!(dict,key)
    elseif valueIsArrayOfArrays(dict[key])
        modifyForArrayOfArrays!(dict,key)
    end
    nothing
end


text:
@testset "HW2" begin
    using DMUStudent.HW2
    using POMDPs
    using DiscreteValueIteration
    using POMDPTesting
    using POMDPModelTools: render

    for n in 1
        m = UnresponsiveACASMDP(n)

        @test all(convert_s(ACASState, convert_s(Int, s, m), m) == s for s in states(m))
        @test has_consistent_transition_distributions(m)

        solver = ValueIterationSolver()
        policy = solve(solver, m)
        v = [value(policy, s) for s in states(m)]

        @test HW2.evaluate(v).score == n
        @test HW2.evaluate(v, "test1234@colorado.edu").score == n

        v[3] += 1.0
        @test HW2.evaluate(v).score == 0

        short = v[1:end-1]
        @test_throws ErrorException HW2.evaluate(short)

        s = rand(initialstate(m))
        sprint(show, render(m, s))
        io = IOBuffer()
        show(io, MIME("text/html"), render(m,s))

        v = convert(Vector{Int}, s)
        a = 0
        @test reward(m, v, a, v) == reward(m, s, a, s)
        @test support(transition(m, v, a)) == support(transition(m, s, a))
        @test isterminal(m, v) == isterminal(m, s)
        @test convert_s(Int, v, m) == convert_s(Int, s, m)
        @test stateindex(m, v) == stateindex(m, s)
    end

    medium = ones(1250*19^3)
    @test HW2.evaluate(medium).score == 0

    large = ones(1250*21^3)
    @test HW2.evaluate(large).score == 0
end


text:
include("Shape.jl")

a = Array(Float32, 1, 2)

function matches(s :: Shape, x)
    return s == shape(x)
end

@assert(matches(shape(a), Array(Float32, 1, 2)))
@assert(!matches(shape(a), Array(Int32, 1, 2)))
@assert(!matches(shape(a), Array(Float32, 1, 3)))
@assert(!matches(shape(a), Array(Float32, 1)))
@assert(!matches(shape(a), Array(Float32, 1, 2, 3)))

s = (Array(Float32, 1, 2), Array(Int32, 1, 3))

@assert(!matches(shape(a), s))
@assert(!matches(shape(s), a))
@assert(matches(shape(s), (Array(Float32, 1, 2), Array(Int32, 1, 3))))
@assert(!matches(shape(s), (Array(Float32, 1, 2),)))
@assert(!matches(shape(s), (Array(Float32, 1, 2), Array(Int32, 1, 3),
                            Array(Int32, 1, 4))))
@assert(!matches(shape(s), (Array(Float32, 1, 3), Array(Int32, 1, 3))))
@assert(!matches(shape(s), (Array(Float32, 1, 2), Array(Int32, 1, 2))))

x = 1.0

@assert(matches(shape(x), 2.0))
@assert(!matches(shape(x), a))
@assert(!matches(shape(x), s))
@assert(!matches(shape(x), 7))


text:
# obs_weight is a shortcut function for getting the relative likelihood of an observation without having to construct the observation distribution. Useful for particle filtering
# maintained by @zsunberg

"""
    obs_weight(pomdp, s, a, sp, o)

Return a weight proportional to the likelihood of receiving observation o from state sp (and a and s if they are present).

This is a useful shortcut for particle filtering so that the observation distribution does not have to be represented.
"""
obs_weight(p, s, a, sp, o) = pdf(observation(p, s, a, sp), o)


text:
using MetabolicFluxAnalysis

s1 = MetabolicFluxAnalysis.Species("s1",4)
s2 = MetabolicFluxAnalysis.Species("s2",0)
s3 = MetabolicFluxAnalysis.Species("s3",3)

# test parseReactand:
# "sA()" --> (1.0,"sA","")
# " 4.0 sB " --> (4.0,"sB","")
# " sC(a-b-X) " --> (1.0,"sC","a-b-X")
 (stoich1,sid1,atm1) = MetabolicFluxAnalysis.parseReactand(" sA() ")
 (stoich2,sid2,atm2) = MetabolicFluxAnalysis.parseReactand(" 4.0 sB ")
 (stoich3,sid3,atm3) = MetabolicFluxAnalysis.parseReactand("  sC(a-b-X) ")

@test stoich1==1.0
@test sid1=="sA"
@test atm1==""
@test stoich2==4.0
@test sid2=="sB"
@test atm2==""
@test stoich3==1.0
@test sid3=="sC"
@test atm3=="a-b-X"

f1 = MetabolicFluxAnalysis.parseFlux( "TestFlux", " A(1-2-3) + B(A-b-x) + B(f-g-j) + X(u) --> C(1-2-3-x-b) + D(A) + D(f) + X(g-j-u) " )
#MetabolicFluxAnalysis.sas(f1)

## benchmark basic functions
#function create_rand_network(n_fluxes::Int64, n_species::Int64 , avg_spec_per_flux:: )


text:
export LLVMType, issized, context, show

# subtypes are expected to have a 'ref::API.LLVMTypeRef' field
abstract type LLVMType end

Base.eltype(typ::LLVMType) = Any

Base.unsafe_convert(::Type{API.LLVMTypeRef}, typ::LLVMType) = typ.ref

const type_kinds = Vector{Type}(fill(Nothing, typemax(API.LLVMTypeKind)+1))
function identify(::Type{LLVMType}, ref::API.LLVMTypeRef)
    kind = API.LLVMGetTypeKind(ref)
    typ = @inbounds type_kinds[kind+1]
    typ === Nothing && error("Unknown type kind $kind")
    return typ
end
function register(T::Type{<:LLVMType}, kind::API.LLVMTypeKind)
    type_kinds[kind+1] = T
end

function refcheck(::Type{T}, ref::API.LLVMTypeRef) where T<:LLVMType
    ref==C_NULL && throw(UndefRefError())
    if Base.JLOptions().debug_level >= 2
        T′ = identify(LLVMType, ref)
        if T != T′
            error("invalid conversion of $T′ type reference to $T")
        end
    end
end

# Construct a concretely typed type object from an abstract type ref
function LLVMType(ref::API.LLVMTypeRef)
    ref == C_NULL && throw(UndefRefError())
    T = identify(LLVMType, ref)
    return T(ref)::LLVMType
end

issized(typ::LLVMType) =
    convert(Core.Bool, API.LLVMTypeIsSized(typ))
context(typ::LLVMType) = Context(API.LLVMGetTypeContext(typ))

function Base.show(io::IO, typ::LLVMType)
    output = unsafe_message(API.LLVMPrintTypeToString(typ))
    print(io, output)
end

Base.isempty(@nospecialize(T::LLVMType)) = false

## integer

export width

@checked struct IntegerType <: LLVMType
    ref::API.LLVMTypeRef
end
register(IntegerType, API.LLVMIntegerTypeKind)

for T in [:Int1, :Int8, :Int16, :Int32, :Int64, :Int128]
    jl_fname = Symbol(T, :Type)
    api_fname = Symbol(:LLVM, jl_fname)
    @eval begin
        $jl_fname(ctx::Context) =
            IntegerType(API.$(Symbol(api_fname, :InContext))(ctx))
    end
end

IntType(bits::Integer; ctx::Context) =
    IntegerType(API.LLVMIntTypeInContext(ctx, bits))

width(inttyp::IntegerType) = API.LLVMGetIntTypeWidth(inttyp)

## floating-point

# NOTE: this type doesn't exist in the LLVM API,
#       we add it for convenience of typechecking generic values (see execution.jl)
 abstract type FloatingPointType <: LLVMType end

for T in [:Half, :Float, :Double, :FP128, :X86_FP80, :PPC_FP128]
    CleanT = Symbol(replace(String(T), "_"=>""))    # only the type kind retains the underscore
    jl_fname = Symbol(CleanT, :Type)
    api_typename = Symbol(:LLVM, CleanT)
    api_fname = Symbol(:LLVM, jl_fname)
    enumkind = Symbol(:LLVM, T, :TypeKind)
    @eval begin
        @checked struct $api_typename <: FloatingPointType
            ref::API.LLVMTypeRef
        end
        register($api_typename, API.$enumkind)

        $jl_fname(ctx::Context) =
            $api_typename(API.$(Symbol(api_fname, :InContext))(ctx))
    end
end

## function types

export isvararg, return_type, parameters

@checked struct FunctionType <: LLVMType
    ref::API.LLVMTypeRef
end
register(FunctionType, API.LLVMFunctionTypeKind)

FunctionType(rettyp::LLVMType, params::Vector{<:LLVMType}=LLVMType[];
             vararg::Core.Bool=false) =
    FunctionType(API.LLVMFunctionType(rettyp, params,
                                      length(params), convert(Bool, vararg)))

isvararg(ft::FunctionType) =
    convert(Core.Bool, API.LLVMIsFunctionVarArg(ft))

return_type(ft::FunctionType) =
    LLVMType(API.LLVMGetReturnType(ft))

function parameters(ft::FunctionType)
    nparams = API.LLVMCountParamTypes(ft)
    params = Vector{API.LLVMTypeRef}(undef, nparams)
    API.LLVMGetParamTypes(ft, params)
    return LLVMType[LLVMType(param) for param in params]
end

## composite types

abstract type CompositeType <: LLVMType end

## sequential types

export addrspace

abstract type SequentialType <: CompositeType end

Base.eltype(typ::SequentialType) = LLVMType(API.LLVMGetElementType(typ))

@checked struct PointerType <: SequentialType
    ref::API.LLVMTypeRef
end
register(PointerType, API.LLVMPointerTypeKind)

function PointerType(eltyp::LLVMType, addrspace=0)
    return PointerType(API.LLVMPointerType(eltyp, addrspace))
end

addrspace(ptrtyp::PointerType) =
    API.LLVMGetPointerAddressSpace(ptrtyp)

@checked struct ArrayType <: SequentialType
    ref::API.LLVMTypeRef
end
register(ArrayType, API.LLVMArrayTypeKind)

function ArrayType(eltyp::LLVMType, count)
    return ArrayType(API.LLVMArrayType(eltyp, count))
end

Base.length(arrtyp::ArrayType) = API.LLVMGetArrayLength(arrtyp)

Base.isempty(@nospecialize(T::ArrayType)) = length(T) == 0 || isempty(eltype(T))

@checked struct VectorType <: SequentialType
    ref::API.LLVMTypeRef
end
register(VectorType, API.LLVMVectorTypeKind)

function VectorType(eltyp::LLVMType, count)
    return VectorType(API.LLVMVectorType(eltyp, count))
end

Base.size(vectyp::VectorType) = API.LLVMGetVectorSize(vectyp)

## structure types

export name, ispacked, isopaque, elements!

@checked struct StructType <: SequentialType
    ref::API.LLVMTypeRef
end
register(StructType, API.LLVMStructTypeKind)

function StructType(name::String; ctx::Context)
    return StructType(API.LLVMStructCreateNamed(ctx, name))
end

StructType(elems::Vector{<:LLVMType}; packed::Core.Bool=false, ctx::Context) =
    StructType(API.LLVMStructTypeInContext(ctx, elems, length(elems), convert(Bool, packed)))

name(structtyp::StructType) =
    unsafe_string(API.LLVMGetStructName(structtyp))
ispacked(structtyp::StructType) =
    convert(Core.Bool, API.LLVMIsPackedStruct(structtyp))
isopaque(structtyp::StructType) =
    convert(Core.Bool, API.LLVMIsOpaqueStruct(structtyp))

elements!(structtyp::StructType, elems::Vector{<:LLVMType}, packed::Core.Bool=false) =
    API.LLVMStructSetBody(structtyp, elems,
                          length(elems), convert(Bool, packed))

Base.isempty(@nospecialize(T::StructType)) =
    isempty(elements(T)) || all(isempty, elements(T))

# element iteration

export elements

struct StructTypeElementSet
    typ::StructType
end

elements(typ::StructType) = StructTypeElementSet(typ)

Base.eltype(::StructTypeElementSet) = LLVMType

function Base.getindex(iter::StructTypeElementSet, i)
    @boundscheck 1 <= i <= length(iter) || throw(BoundsError(iter, i))
    return LLVMType(API.LLVMStructGetTypeAtIndex(iter.typ, i-1))
end

function Base.iterate(iter::StructTypeElementSet, i=1)
    i >= length(iter) + 1 ? nothing : (iter[i], i+1)
end

Base.length(iter::StructTypeElementSet) = API.LLVMCountStructElementTypes(iter.typ)

Base.lastindex(iter::StructTypeElementSet) = length(iter)

# NOTE: optimized `collect`
function Base.collect(iter::StructTypeElementSet)
    elems = Vector{API.LLVMTypeRef}(undef, length(iter))
    API.LLVMGetStructElementTypes(iter.typ, elems)
    return LLVMType[LLVMType(elem) for elem in elems]
end

## other

@checked struct VoidType <: LLVMType
    ref::API.LLVMTypeRef
end
register(VoidType, API.LLVMVoidTypeKind)

VoidType(ctx::Context) = VoidType(API.LLVMVoidTypeInContext(ctx))

@checked struct LabelType <: LLVMType
    ref::API.LLVMTypeRef
end
register(LabelType, API.LLVMLabelTypeKind)

LabelType(ctx::Context) = LabelType(API.LLVMLabelTypeInContext(ctx))

@checked struct MetadataType <: LLVMType
    ref::API.LLVMTypeRef
end
register(MetadataType, API.LLVMMetadataTypeKind)

MetadataType(ctx::Context) = MetadataType(API.LLVMMetadataTypeInContext(ctx))

@checked struct TokenType <: LLVMType
    ref::API.LLVMTypeRef
end
register(TokenType, API.LLVMTokenTypeKind)

TokenType(ctx::Context) = TokenType(API.LLVMTokenTypeInContext(ctx))

## type iteration

export types

struct ContextTypeDict <: AbstractDict{String,LLVMType}
    ctx::Context
end

# FIXME: remove on LLVM 12
function LLVMGetTypeByName2(ctx::Context, name)
    Module("dummy"; ctx) do mod
        API.LLVMGetTypeByName(mod, name)
    end
end

types(ctx::Context) = ContextTypeDict(ctx)

function Base.haskey(iter::ContextTypeDict, name::String)
    return LLVMGetTypeByName2(iter.ctx, name) != C_NULL
end

function Base.getindex(iter::ContextTypeDict, name::String)
    objref = LLVMGetTypeByName2(iter.ctx, name)
    objref == C_NULL && throw(KeyError(name))
    return LLVMType(objref)
end


text:
using SISLES
using Util
using HDF5, JLD
using PyPlot

const PARAMFILE = "../Encounter/CorrAEMImpl/params/cor.txt"
const INITFILE = "rr_initial.txt"
const TRANFILE = "rr_transition.txt"
const NSAMPLES = 10000
const NTRANS = 50

aem = CorrAEM(PARAMFILE, INITFILE, NSAMPLES,
        TRANFILE, NTRANS)

function validate_initial(aem)

    params = aem.parameters

    marginal_dist = Array{Any}(params.n_initial)
    marginal_prob = Array{Any}(params.n_initial)

    for i = 1:params.n_initial
        marginal_dist[i] = vec(sum(params.N_initial[i], 2))
        marginal_prob[i] = marginal_dist[i] / sum(marginal_dist[i])
    end

    f = open(aem.initial_sample_filename, "r")

    readline(f)

    n_nmac = 0
    n_lines = 0

    sample_dist = Array{Any}(params.n_initial)
    sample_prob = Array{Any}(params.n_initial)

    for i = 1:params.n_initial
        sample_dist[i] = zeros(Int, params.r_initial[i])
    end

    for line = eachline(f) #(ix, line) in enumerate(eachline(f))
        values = float(split(chomp(line)))

        # 1       id
        # 2       A = [1:4]
        # 3       L = [1:5]
        # 4       chi = [90, 270]   : deg
        # 5       beta = [0:30:360] : deg
        # 6       C1 = [1, 2]
        # 7       C2 = [1, 2]
        # 8,9     v1, v2 = [50, 100, 200, 300, 400, 500, 600]   : kt
        # 10, 11  v1d, v2d = [-5, -2, -0.25, 0.25, 2, 5]        : kt/s
        # 12, 13  h1d, h2d = [-5000, -3000, -2000, -1000, -400, 400, 1000, 2000, 3000, 5000]    : ft/min
        # 14, 15  psi1d, psi2d = [-8, -6, -3.5, -1, -0.25, 0.25, 1, 3.5, 6, 8]  : deg/s
        # 16      hmd = [0, 0.0822896, 0.5, 1, 3]   : nm
        # 17      vmd = [0, 100, 200, 300, 400, 500, 600, 700, 800, 900, 6000]  : ft

        for i = [2, 3, 4, 6, 7]
            if !(Int(values[i]) in 1:params.r_initial[i-1])
                println("sample: ", chomp(line))
                println("i: ", i, ", value: ", Int(values[i]), ", boundary: ", params.r_initial[i-1])
                error("The value violates the boundary.")
            end
            #println("Updating sample_dist[", i-1, "][", Int(values[i]), "] += 1")
            sample_dist[i - 1][Int(values[i])] += 1
        end

        for i = [5, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]
            if !isempty(params.boundaries[i-1]) && (values[i] < params.boundaries[i-1][1] || values[i] > params.boundaries[i-1][end])
                println("sample: ", chomp(line))
                println("i: ", i, ", value: ", values[i], ", boundary: ", params.boundaries[i-1][end])
                error("The value violates the boundary.")
            end

            j = findfirst(x -> (x > values[i]), params.boundaries[i - 1])

            if j == 0 && values[i] == params.boundaries[i - 1][end]
                sample_dist[i - 1][params.r_initial[i - 1]] += 1
            else
                sample_dist[i - 1][j - 1] += 1
            end
        end

        hmd_ft = 6076.12 * values[16]
        vmd = values[17]

        if hmd_ft <= 500 && vmd <= 200
            n_nmac += 1
        end

        n_lines += 1
    end

    close(f)

    SimulationResult = load("result.jld", "data")
    labels, initial, AC1_trajectory, AC2_trajectory = SimulationResult

    println("n_nmac: ", n_nmac, " n_lines: ", n_lines)
    println()
    println("P(NMAC | enc, no TCAS) = ", n_nmac / n_lines, " ($n_nmac / $n_lines)")

    println()
end


text:
import Base.==

"""
Define the meaning of "==" operator
"""
function =={T<:AbstractMolecularContainer}(atom1::T, atom2::T)
    tol = 1e-5 # tolerance for floating point number comparison

    for field in fieldnames(atom1)
        value1 = getfield(atom1, field)
        value2 = getfield(atom2, field)

        if issubtype(typeof(value1), AbstractFloat)
            if abs(value1 - value2) > tol
                return false
            end

        elseif issubtype(typeof(value1), AbstractArray)
            if length(value1) != length(value2)
                return false
            end

            for i = 1:length(value1)
                if issubtype(typeof(value1[i]), AbstractFloat)
                    if abs(value1[i] - value2[i]) > tol
                        return false
                    else
                        continue
                    end
                else
                    if value1[i] != value2[i]
                        return false
                    else
                        continue
                    end
                end
            end

            return true

        else
            if value1 != value2
                return false
            end

        end
    end

    return true
end


text:
module LimitOfDetectionScript

using SynapseClient
using DISSEQT
using DISSEQT.SynapseTools
using DISSEQT.AnnotatedArrays
using Statistics
using DataFrames
using CSV
using JLD

function cachemetadata(syn,alignmentFolder,metadataID)
    metadataCacheFilename = "metadatacache.csv"
    metadata = DataFrame()
    if isfile(metadataCacheFilename)
        metadata = CSV.read(metadataCacheFilename; missingstrings=["","NA"])
    else
        # Get metadata (metadataID can be the synapse ID of a .csv file or of a folder containing .csv files).
        metadata = getsamplemetadata(syn, metadataID, :Rejected=>"No") # Always remove rejected samples

        # extend metadata
        appendswarmids!(syn,metadata,alignmentFolder,[:both,:forward,:reverse])
        appendalignedids!(syn,metadata,alignmentFolder)

        # download swarms and consensus sequences
        downloadswarms!(syn,metadata,cacheDir="MutantSwarms",[:both,:forward,:reverse])
        downloadconsensuses!(syn,metadata,cacheDir="MutantSwarms")

        CSV.write(metadataCacheFilename, metadata; missingstring="NA")
    end
    metadata
end

function main()
    syn = SynapseClient.login()

    # prohectFolder should point to MyProject folder (containing subfolders Metadata and Analysis)
    projectFolder   = "syn11639899" # VignuzziLabPublic/Projects/FitnessLandscapes
    analysisFolder  = getchildbyname(syn, projectFolder, "Analysis")
    alignmentFolder = getchildbyname(syn, analysisFolder, "Alignment")
    referenceFolder = getchildbyname(syn, alignmentFolder, "ReferenceGenomes")

    isdir("MutantSwarms") || mkdir("MutantSwarms")
    isdir("MutantSwarms/forward") || mkdir("MutantSwarms/forward")
    isdir("MutantSwarms/reverse") || mkdir("MutantSwarms/reverse")

    # set to true to upload files
    doUpload = true
    uploadName = "fitness landscapes invitro"

    metadataID = getchildbyname(syn, projectFolder, "Metadata", "fitness landscapes invitro.csv")
    metadata = cachemetadata(syn, alignmentFolder, metadataID)

    # get CVB3 feature list
    featureFileID = getchildbyname(syn,referenceFolder,"CVB3_features.csv")
    features = CSV.read(localpath(syn, featureFileID); missingstrings=["","NA"])

    positions = featurepositions(features, "ORF") # coding positions

    # load swarms
    swarms  = loadswarm(metadata[:SampleID],metadata[:SwarmPath],metadata[:ConsensusPath])
    swarmsF = loadswarm(metadata[:SampleID],metadata[:SwarmForwardPath],metadata[:ConsensusPath])
    swarmsR = loadswarm(metadata[:SampleID],metadata[:SwarmReversePath],metadata[:ConsensusPath])

    coverage = dropdims(swarms[:coverage]; dims=(1,2,3))
    meanCoverage = mean(coverage[positions,:]; dims=1)[:]

    # filter by mean read coverage at coding sites
    sampleMask = meanCoverage .>= 1000
    println("Removed ", count(.~sampleMask), " samples due to low coverage.")
    swarms  = swarms[:,:,:,:,sampleMask]
    swarmsF = swarmsF[:,:,:,:,sampleMask]
    swarmsR = swarmsR[:,:,:,:,sampleMask]
    metadata = metadata[sampleMask,:]

    println("Total number of samples used: ", count(sampleMask))

    uploads = UploadList()
    dependencies = unique(vcat(featureFileID, metadata[:SwarmID], metadata[:SwarmForwardID], metadata[:SwarmReverseID], metadata[:MetadataID]))

    lod = limitofdetection(swarms, swarmsF, swarmsR, groupID=metadata[:Run])
    lodDict = Dict( string(p[1])=>p[2] for p in pairs(convert(Dict,lod)) )
    lodfn = "limitofdetection.jld"
    save(lodfn, lodDict)
    markforupload!(uploads, "", lodfn, dependencies)

    doUpload && uploadfiles(syn, createchildfolder(syn, analysisFolder, "LimitOfDetection", uploadName), uploads, @__FILE__, "Limit of Detection")

    #swarms,swarmsF,swarmsR,lod,metadata
    nothing
end

end
LimitOfDetectionScript.main();


text:
# This script process the simulation data using non-linear least squares fitting
# for both the original Reference Region Model (NRRM) and the Constrained approach (CNRRM)

# Run this script by launching julia then entering: include("path/to/e01_2a_simProcessNLSQ.jl")

# Estimated runtimes:
# ~2.6 hours for NRRM
# ~1.2 hours for CNRRM
# So total run time is almost 4 hours
# Note: These runtimes were obtained using 4 threads.
#       If only a single core is used, then the runtime could be 4 times longer

# Pick the choice of parameters for the reference tissue parameters
# (This should match the choice from step e01_1)
refName = "refY"

# The temporal resolution (in seconds) to process over
listTRes = [1,5,10,15,30,60]
# The choice of Contrast-Noise Ratios
# (This does not change noise, rather it is only here to remind the code what the CNRs are. Should match values from simMaker.)
listCNR = collect(5:5:50)

# First, figure out the current file path
curFilePath = Base.source_path() # Path to current file
workDir = dirname(curFilePath) # Our working directory (directory of current file)
## If the above doesn't work, then user has to manually enter the location of the CRRM directory
## i.e. wordDir = "/path/to/CLRRM/Directory"

import Hwloc
topology = Hwloc.topology_load()
counts = Hwloc.histmap(topology)
nCores = counts[:Core]
if nprocs()<nCores
  addprocs(nCores-nprocs()+1)
end

## Load packages
using DCEMRI  # Needed for the levenberg-marquardt fitting
using MAT      # Needed for loading/saving .mat files

# Some required directories/files
matDir = workDir * "/data/simData/$refName/rawData"
warmupFile = workDir * "/data/QIBA-ToftsV6/QIBAv6-Mini4Jl-Warmup.mat"
auxCode = workDir * "/jlfiles/auxCode.jl"
refRegCode = workDir * "/jlfiles/refRegionFunction.jl"

outDir = workDir * "/data/simData/$refName/NRRM"

# Make the output directory if it doesn't already exist
if !(isdir(outDir))
  mkdir(outDir)
end

# Load pre-requisite code
include(auxCode)

## Warmup run to make sure all is well
println("==== Warmup Run ====")
results = fitdata(datafile=warmupFile, models=[2])
rm("output.mat") # Delete output since it isn't needed

# If warmup succeeds, then proceed by loading Reference Region Model
include(refRegCode)

# Define the reference tissue parameters
# Default is 'refY'
refKTrans = 0.1;
refVe = 0.1;
if (refName == "refWS")
  refKTrans = 0.07;
  refVe = 0.14;
elseif (refName == "refP")
  refKTrans = 0.137;
  refVe = 0.115;
end
refKep = refKTrans/refVe

# Get the names of the mat files containing simulated data
matFiles = readdir(matDir)
numFiles = length(matFiles)

println("")
println("===============================")
println("==== Start processing data ====")
# Loop through the simulated datasets
for q=1:numFiles
  # Load the simulated data for each CNR
  curCNR = listCNR[q]
  println("")
  println("------------")
  println("Processing CNR = $curCNR")
  curFile = matDir * "/Sim-CNR-$curCNR.mat"
  matData = matread(curFile)
  cnrInd = listCNR[q]
  for i=1:length(listTRes)
    # Appropriately downsample simulated data to obtain desired temporal resolutions
    curTRes = listTRes[i]
    println("")
    println("---")
    println("Temporal Resolution = $curTRes s")
    # Downsample the data
    Ct = downsample(matData["CtNoisy"],curTRes)
    Crr = downsample(vec(matData["CrrClean"]),curTRes)
    t = downsample(vec(matData["T"]),curTRes)
    # Build a dummy mask (ones everywhere)
    (nT, nV) = size(Ct)
    # Non-Linear Reference Region Model
    # (Unconstrained and Constrained are both run if estKepRR is set to zero)
    (pkParamsCN, residCN, medianKepRR, pkParamsN, residN, runtimeN, runtimeC) = fitCNRRM(Ct, Crr, t, estKepRR=0.0, doTime=true)
    # Output results
    results = Dict()
    results["pkParamsN"] = pkParamsN
    results["residN"] = residN
    results["pkParamsCN"] = pkParamsCN
    results["residCN"] = residCN
    results["runtimeN"] = runtimeN
    results["runtimeC"] = runtimeC
    results["medianKepRR"] = medianKepRR
    matwrite(outDir * "/Sim-CNR-$curCNR-TRes-$curTRes.mat", results)
  end
end
println("==== Finished ====")
println("==================")

max_stars_repo_name: MPUmri/CRRM

max_issues_repo_name: MPUmri/CRRM

max_forks_repo_name: MPUmri/CRRM


text:
function POMDPs.observation(pomdp::DoubleCovidPOMDP, s::DoubleCovidState, a::CovidAction, sp::DoubleCovidState)
    iszero(a.testing_prop) && return CovidPOMDPs.UninformedDist() # not getting any info; any transition equally likely
    tot_mean = 0.0
    tot_variance = 0.0
    p1 = pomdp.params
    p2 = s.params
    for (i,inf) in enumerate(s.I1)
        num_already_tested = sum(@view s.T1[:,i])
        num_tested = floor(Int,(inf-num_already_tested)*a.testing_prop)
        dist = Binomial(num_tested,p1.pos_test_probs[i])
        tot_mean += Statistics.mean(dist)
        tot_variance += Statistics.var(dist)
    end

    for (i,inf) in enumerate(s.I2)
        num_already_tested = sum(@view s.T2[:,i])
        num_tested = floor(Int,(inf-num_already_tested)*a.testing_prop)
        dist = Binomial(num_tested,p2.pos_test_probs[i])
        tot_mean += Statistics.mean(dist)
        tot_variance += Statistics.var(dist)
    end

    return Normal(tot_mean, sqrt(tot_variance))
end