text:
module SynapseUploadUtils

export FolderInfo,
       listfiles,
       confirmupload,
       uploadfolder,
       copyfolder

using SynapseClient
import SynapseClient: AbstractEntity, Project, Folder, File, Activity

using ArgParse
using Printf

struct FolderInfo
	path
	name
	files::Array{AbstractString,1}
	folders::Array{FolderInfo,1}
end
FolderInfo(path,name) = FolderInfo(path,name,[],[])

function describefolder!(res::Dict{AbstractString,Int}, fi::FolderInfo)
	totalSize = 0
	nbrFolders = 1
	for file in fi.files
		k = splitext(file)[2]
		res[k] = get(res,k,0)+1
		totalSize += filesize(joinpath(fi.path,file))
	end
	for folder in fi.folders
		nf,ts = describefolder!(res,folder)
		nbrFolders += nf
		totalSize += ts
	end
	nbrFolders, totalSize
end

function describefolder(fi::FolderInfo)
	res = Dict{AbstractString,Int}()
	nbrFolders, totalSize = describefolder!(res,fi)
	nbrFolders, totalSize, res
end

function listfiles(path::AbstractString; excludeFiles=[], excludeFolders=[])
	splitPath = split(path, ['/','\\'])
	while isempty(splitPath[end])
		pop!(splitPath)
	end
	isempty(splitPath) && error("\"$path\" is not a valid path")
	name = splitPath[end]

	folderInfo = FolderInfo(path, name)

	for file in readdir(path)
		isempty(file) && continue
		file[1] == '.' && continue

		fullFilePath = joinpath(path, file)
		if isdir(fullFilePath)
			any(r->ismatch(r,file), excludeFolders) && continue
			# println("folder: $file")
			push!(folderInfo.folders, listfiles(fullFilePath, excludeFiles=excludeFiles, excludeFolders=excludeFolders))
		else
			any(r->ismatch(r,file), excludeFiles) && continue
			filesize(fullFilePath)==0 && continue # Synapse doesn't allow uploading of empty files
			# println("file: $file")
			push!(folderInfo.files, file)
		end
	end

	folderInfo
end

function getchildbyname(syn, parentID::AbstractString, child::AbstractString)::String
    # if 0 results, return "",
    # if 1 result, return it
    # if 2 results, error

    res = ""
    for c in getchildren(syn, parentID)
        if get(c,"name",nothing) == child
            @assert isempty(res) "Unexpected error, multiple children with the same name."
            res = c["id"]
        end
    end
    res
end
getchildbyname(syn, parent::AbstractEntity, child::AbstractString) = getchildbyname(syn, parent.id,child)
getchildbyname(syn, parent, args::AbstractString...) = getchildbyname(syn,getchildbyname(syn, parent, args[1]),args[2:end]...)

function retrystore(syn::Synapse, args...; kwargs...)
	delays = Int[5, 30, 60]
	i = 0 # how many times did we fail?
	while true
		try
			return store(syn, args...; kwargs...)
		catch ex
			i += 1
			if i>length(delays)
				println("Synapse store failed. Aborting.")
				rethrow(ex)
			end

			println("Synapse store failed. Waiting $(delays[i])s before trying again.")
			sleep(delays[i])
		end
	end
end

function askforconfirmation(str::AbstractString)
	while true
		println(str, " (y/n)")
		response = strip(readline())
		lowercase(response) in ["y","yes"] && return true
		lowercase(response) in ["n","no"] && return false
	end
end

function nbrbytes2string(x::Integer)
	x < 1024 && return string(x, " bytes")
	y = Float64(x)/1024
	y < 1024 && return @sprintf("%.1f KB", y)
	y /= 1024
	y < 1024 && return @sprintf("%.1f MB", y)
	y /= 1024
	y < 1024 && return @sprintf("%.1f GB", y)
	y /= 1024
	@sprintf("%.1f TB", y)
end

function fullsynapsepath(syn::Synapse, id::AbstractString)
	try
		entity = get(syn, id, downloadFile=false)
		name = entity.name

		typeof(entity) <: Folder || return name # i.e. go upwards until we find the parent project
		return string(fullsynapsepath(syn,entity.parentId), '/', name)
	catch
	end
	return "[UNKNOWN]"
end

function confirmupload(syn::Synapse, parentFolderID::AbstractString, fi::FolderInfo; askForConfirmation=true)
	synapsePath = fullsynapsepath(syn, parentFolderID)

	child = getchildbyname(syn, parentFolderID, fi.name)
	if askForConfirmation && !isempty(child)
		askforconfirmation("Folder \"$(fi.name)\" already exists in \"$synapsePath\", continue?") || return false
	end

	nbrFolders, totalSize, desc = describefolder(fi)

	println("--- Summary for $(fi.name) ---")
	println("Number of folders: $nbrFolders")
	println("Number of files:")
	for (k,v) in desc
		println("\t$k: $v")
	end
	println("Total size: ", nbrbytes2string(totalSize))

	if askForConfirmation
		askforconfirmation("Upload to \"$synapsePath\"?") || return false
	end
	true
end

function _uploadfolder(syn::Synapse, parentID::AbstractString, fi::FolderInfo, exec::AbstractString)
	folder = getchildbyname(syn, parentID, fi.name)
	if isempty(folder)
		# create folder
		folder = Folder(fi.name, parent=parentID)
		act = Activity(name="Uploaded folder")
		isempty(exec) || executed(act,exec)
		folder = retrystore(syn,folder,activity=act)
	end

	annot = getannotations(syn, folder)
	annot.uploading = "In progress"
	setannotations(syn, folder, annot)

	for filename in fi.files
		# TODO: set contentType?
		file = File(path=joinpath(fi.path,filename), name=filename, parent=folder)
		act = Activity(name="Uploaded file")
		isempty(exec) || executed(act,exec)
		file = retrystore(syn,file,activity=act)
	end

	for subfolder in fi.folders
		_uploadfolder(syn,folder,subfolder,exec)
	end

	annot = getannotations(syn, folder)
	annot.uploading = "Finished"
	setannotations(syn, folder, annot)
end

_uploadfolder(syn::Synapse, parent::Folder, fi::FolderInfo, exec::AbstractString) = _uploadfolder(syn, parent.id, fi, exec)

function uploadfolder(syn::Synapse, parentFolderID::AbstractString, fi::FolderInfo; executed="")
	_uploadfolder(syn,parentFolderID,fi,executed)
end

# Upload folder
# 	1. Identify all folders and files recursively.
# 		Ignore hidden files/folders.
# 		Ignore files/folders starting with ".".
# 	2. Ask for confirmation to continue if folder already exists in Synapse.
# 	3. Show summary of what will be uploaded:
# 		Name of folder to be uploaded.
# 		Nbr of folders.
# 		Nbr of files with different file endings.
# 		Total size.
# 	4. Ask for confirmation to continue.
# 	5. Create the folder specified in "Your Project/Your Folder".
# 	6. Add annotation that upload is in progress
# 	7. Upload files.
# 		Create subfolders as needed. (Thus, if the upload stops, we will see how far it got.)
# 	8. Add annotation that upload finished.

function parse_uploadfolder_commandline(ARGS)
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--yes", "-y"
            help = "Upload without asking for confirmation"
            action = :store_true
        "--parent", "-p"
        	help = "Parent folder (synapse ID)"
            arg_type = UTF8String
        	default = UTF8String("syn6177609")
        "--excludefile", "-e"
            help = "Exclude files matching pattern"
            arg_type = Regex
            action = :append_arg
        "--excludefolder", "-x"
            help = "Exclude folders matching pattern"
            arg_type = Regex
            action = :append_arg
        "folders"
            help = "Folders to upload"
            nargs = '*'
            arg_type = UTF8String
            #action = :append_arg
            required = true
    end

    return parse_args(ARGS, s)
end

function uploadfolder(ARGS)
	parsedArgs = parse_uploadfolder_commandline(ARGS)
	askForConfirmation = !parsedArgs["yes"]
	sources = parsedArgs["folders"]

	excludeFiles = parsedArgs["excludefile"]
	excludeFolders = parsedArgs["excludefolder"]

	parentFolderID = parsedArgs["parent"]

	map!(abspath,sources)

	syn = SynapseClient.login()

	# prepare
	folders = Array{FolderInfo,1}(length(sources))
	map!( s->listfiles(s,excludeFiles=excludeFiles,excludeFolders=excludeFolders), folders, sources)

	# check that it is ok to upload each folder
	for fi in folders
		confirmupload(syn, parentFolderID, fi, askForConfirmation=askForConfirmation) || exit(0) # error("User abort")
	end

	# upload each folder
	for fi in folders
		uploadfolder(syn, parentFolderID, fi, executed="https://github.com/rasmushenningsson/SynapseUpload.jl/blob/master/synapseupload.jl")
	end
end

struct FileCopyInfo
	sourceID::AbstractString
	name::AbstractString
	# some folder info is needed for flatten=false
end

function listsynapsefiles!(syn::Synapse, list::Vector{FileCopyInfo}, parentID::AbstractString; recursive=true)
    for c in getchildren(syn, parentID)
		id = c["id"]
		entity = get(syn, id, downloadFile=false)

		if typeof(entity) == File
			push!(list, FileCopyInfo(id, entity.name))
		elseif typeof(entity) == Folder && recursive
			listsynapsefiles!(syn, list, id, recursive=recursive)
		end
    end
    list

	# results = chunkedquery(syn, "select id from entity where entity.parentId=='$parentID'")
	# for r in results
	# 	id = r["entity.id"]
	# 	entity = get(syn, id, downloadFile=false)

	# 	if typeof(entity) == File
	# 		push!(list, FileCopyInfo(id, entity["name"]))
	# 	elseif typeof(entity) == Folder && recursive
	# 		listsynapsefiles!(syn, list, id, recursive=recursive)
	# 	end
	# end
	# list
end

function listsynapsefiles(syn::Synapse,parentID::AbstractString; recursive=true)
	list = Array{FileCopyInfo,1}()
	listsynapsefiles!(syn, list, parentID, recursive=recursive)
end

function matchany(names::Array{T}, patterns::Array{Regex,1}) where {T<:AbstractString}
	Bool[any(pattern->match(pattern,name)!=nothing, patterns) for name in names]
end

function performcopy(syn::Synapse, list::Array{FileCopyInfo,1}, destinationID::AbstractString, exec::AbstractString)
	profile = getuserprofile(syn)

	for f in list
		file = get(syn, f.sourceID, downloadFile=false)
		@assert typeof(file)==File

		act = Activity(name="Copied file", used=file)
		isempty(exec) || executed(act,exec)

		local newFile

		if profile.ownerId == file.createdBy # no need to copy if the file has the same creator
	        newFile = File(name=file.name, parentId=destinationID)
	        newFile.dataFileHandleId = file.dataFileHandleId
		else # fallback to copy
            file = get(syn,f.sourceID,downloadFile=true) # download file
            newFile = File(path=file.path, name=file.name, parentId=destinationID)
		end
		retrystore(syn,newFile,activity=act)
	end
end

function copyfolder(syn::Synapse, folderName::AbstractString, parentID::AbstractString,
	                destinationParentID::AbstractString, patterns::Array{Regex,1};
	                destinationName::AbstractString="", flatten=true, askForConfirmation=true,
	                localScript="", externalScript="")
	flatten == true || error("Only flatten=true supported at the moment.")
	!isempty(localScript) && !isempty(externalScript) && error("Only one of localScript and externalScript can be specified.")

	if askForConfirmation && isempty(localScript) && isempty(externalScript)
		askforconfirmation("Please consider specifying localScript or externalScript, continue anyway?") || return false
	end

	destinationName = isempty(destinationName) ? folderName : destinationName

	synapsePath = fullsynapsepath(syn, destinationParentID)
	localScriptName = "copy_$(destinationName).jl" # only used if localScript is set

	if askForConfirmation
		destinationID = getchildbyname(syn, destinationParentID, destinationName)
		if !isempty(destinationID)
			askforconfirmation("Folder \"$(destinationName)\" already exists in \"$synapsePath\", continue?") || return false
		end
	end

	sourceFolderID = getchildbyname(syn,parentID,folderName)
	sourceFolder = get(syn,sourceFolderID,downloadFile=false)
	@assert typeof(sourceFolder)==Folder "Source should be a folder."

	list = listsynapsefiles(syn, sourceFolderID, recursive=true)
	mask = matchany(map(x->x.name,list), patterns)
	skipped = list[.!mask]
	list = list[mask] # only keep files with names matching at least one of the patterns

	if askForConfirmation
		println("--- Summary for $destinationName ---")
		println("-Skipped-")
		for f in skipped
			println("\t$(f.name)")
		end

		println("-Uploading-")
		for f in list
			println("\t$(f.name)")
		end
		isempty(localScript) || (fn = splitdir(localScript)[2]; println("Local script \"$fn\" will be uploaded to \"$synapsePath/$localScriptName\"."))
		isempty(externalScript) || println("External script: \"$externalScript\" is used for provenance.")
		askforconfirmation("Copy files to \"$synapsePath\"?") || return false
	end

	if flatten
		names = map(x->x.name,list)
		length(names) != length(unique(names)) && error("Could not flatten file structure due to name collisions.")
	end

	# create destination folder
	destFolder = Folder(name=destinationName,parentId=destinationParentID)
	destFolder = retrystore(syn, destFolder)

	# setup executed
	executed = ""
	if !isempty(localScript)
		# upload script and set executed as synapse id

		scriptFile = File(path=localScript, name=localScriptName, parentId=destinationParentID)

		# Had to remove this line since "FileNameOverride" is deprecated. It doesn't break anything, but if the script is downloaded, it will use the old filename which is annoying.
		#scriptFile["fileNameOverride"] = localScriptName # standardize the script file name (i.e. use same name as for uploaded folder)

		scriptFile = retrystore(syn, scriptFile)
		executed = scriptFile.id
	elseif !isempty(externalScript) # just refer to external script
		executed = externalScript
	end

	performcopy(syn,list,destFolder.id, executed) # copy files
	true
end

copyfolder(syn::Synapse, folderName::AbstractString, parentID::AbstractString, destinationParentID::AbstractString, pattern::Regex; kwargs...) =
	copyfolder(syn, folderName, parentID, destinationParentID, [pattern]; kwargs...)
copyfolder(syn::Synapse, folderName::AbstractString, parentID::AbstractString, destinationParentID::AbstractString, patterns::Array{AbstractString,1}; kwargs...) =
	copyfolder(syn, folderName, parentID, destinationParentID, map(Regex,pattern); kwargs...)
copyfolder(syn::Synapse, folderName::AbstractString, parentID::AbstractString, destinationParentID::AbstractString, pattern::AbstractString; kwargs...) =
	copyfolder(syn, folderName, parentID, destinationParentID, Regex(pattern); kwargs...)

end


text:
# Everything dealing with TWO classifier objects

function mpm_classifier(data1, data2; burn=1000, thin=50, d1=100.0, d2=100.0, kappa=10.0, usepriors=true)
    assert(size(data1,2) == size(data1,2))
    D = size(data1, 2)
    cov = eye(D,D) .* 0.1
    mu = ones(D)
    prior = MPM.MPMPrior(D=D, kappa=kappa, S=eye(D)*0.05*10)

    # Class 1
    start = MPM.MPMParams(mu, cov,
        clamp(log(data1./d1)',-8.0,Inf)) #lam
    obj_a = MPM.MPMSampler(prior, data1, deepcopy(start), d1)
    obj_a.usepriors = usepriors
    #mymh_a = MPM.MHRecord(obj_a,burn=burn,thin=thin)
    mymh_a = MPM.AMWGRecord(obj_a, blocks, burn=burn,thin=thin)

    # Class 2
    start.lam = clamp(log(data2./d2)',-8.0,Inf) # lam
    obj_b = MPM.MPMSampler(prior, data2, deepcopy(start), d2)
    obj_b.usepriors = usepriors
    #mymh_b = MPM.MHRecord(obj_b,burn=burn,thin=thin)
    mymh_b = MPM.AMWGRecord(obj_b, blocks, burn=burn,thin=thin)

    OBC.BinaryClassifier(obj_a, obj_b, mymh_a, mymh_b)
end

function sample!(cls::OBC.BinaryClassifier, iters=10000; verbose=false)
    t0 = time()
    OBC.sample!(cls.mcmc1, iters, verbose=verbose)
    t1 = time()
    OBC.sample!(cls.mcmc2, iters, verbose=verbose)
    t2 = time()
    return t1-t0, t2-t1
end

import Base: vcat

function vcat(xs::OBC.BinaryClassifier...)
    x = xs[1]
    for y in xs[2:end]
        append!(x.mcmc1.db, y.mcmc1.db)
        append!(x.mcmc2.db, y.mcmc2.db)
    end
    return x
end

function gelman_rubin(samplers::Vector{OBC.BinaryClassifier})
    return [gelman_rubin([x.mcmc1 for x in samplers]), gelman_rubin([x.mcmc1 for x in samplers])]
end

function acceptance_rates(cls::OBC.BinaryClassifier)
    mc1,mc2 = cls.mcmc1, cls.mcmc2
    accs = Dict{String,Any}()
    accs["1"] = cls.mcmc1.block_accept
    accs["2"] = cls.mcmc2.block_accept
    return accs
end

############################################################
# BEE Moment calculations
############################################################

function bee_moments(cls::OBC.BinaryClassifier; dmean=10.0, numlam=20, numpts=20)
    bee_moments(cls, (dmean,dmean), numlam=numlam, numpts=numpts)
end

function bee_moments(cls::OBC.BinaryClassifier, dmeans::NTuple{2,Float64}; numlam=20, numpts=20)
    # FIXME This is only valid for c=0.5
    dmean1,dmean2 = dmeans

    db1 = cls.mcmc1.db
    db2 = cls.mcmc2.db
    assert(length(db1) == length(db2))
    dblen = length(db1)
    dim = length(db1[1].mu)
    points = Array(Int, dim, numpts*2)
    lamsx1 = Array(Float64, dim, numpts)
    lamsx2 = Array(Float64, dim, numpts)
    lams1 = Array(Float64, dim, numlam)
    lams2 = Array(Float64, dim, numlam)
    llams1 = Array(Float64, dim, numlam)
    llams2 = Array(Float64, dim, numlam)
    rlams1 = Array(Float64, dim, numlam)
    rlams2 = Array(Float64, dim, numlam)

    beem1 = beem2 = 0.0

    for i in 1:dblen
        dbpass = 0.0
        #dbpass = 0.0
        curr1 = db1[i]
        curr2 = db2[i]
        # Generate x values (from lams or from new lams?)
        rand!(MultivariateNormal(curr1.mu, curr1.sigma), lamsx1)
        rand!(MultivariateNormal(curr2.mu, curr2.sigma), lamsx2)

        if any(lamsx1 .> 35) || any(lamsx2 .> 35)
            #warn("Encountered lambda value greater than 35, skipping iteration")
            continue
        end

        for i=1:dim, j=1:numpts
            points[i,j] = rand(Poisson(dmean1*exp(lamsx1[i,j])))
            points[i,j+numpts] = rand(Poisson(dmean2*exp(lamsx2[i,j])))
        end

        g1 = calc_g(points, db1, numlam, dmean=dmean1)
        g2 = calc_g(points, db2, numlam, dmean=dmean2)
        g1 = clamp(g1, -10_000_000.0, Inf)
        g2 = clamp(g2, -10_000_000.0, Inf)

        # Now compute p(y|x,\theta)
        rand!(MultivariateNormal(curr1.mu, curr1.sigma), lams1)
        rand!(MultivariateNormal(curr2.mu, curr2.sigma), lams2)
        for k=1:numlam*dim
            rlams1[k] = dmean1 * exp(lams1[k])
            llams1[k] = log(dmean1) + lams1[k]
            rlams2[k] = dmean2 * exp(lams2[k])
            llams2[k] = log(dmean2) + lams2[k]
        end
        for j in 1:numpts*2 # each point
            accumlam1 = 0.0
            accumlam2 = 0.0
            #accumlam1 = -Inf
            #accumlam2 = -Inf
            for s in 1:numlam # each lambda
                accumD1 = 0.0
                accumD2 = 0.0
                for d in 1:dim # each dimension
                    accumD1 += points[d,j]*llams1[d,s] - rlams1[d,s] - lgamma(points[d,j]+1)
                    accumD2 += points[d,j]*llams2[d,s] - rlams2[d,s] - lgamma(points[d,j]+1)
                end
                accumlam1 += exp(accumD1)
                accumlam2 += exp(accumD2)
                #accumlam1 = accumlam1 == -Inf ? accumD1 : logsum(accumlam1, accumD1)
                #accumlam2 = accumlam2 == -Inf ? accumD2 : logsum(accumlam2, accumD2)
            end
            temp = g1[j] < g2[j] ? accumlam1 : accumlam2 #FIXME c!=0.5
            z = accumlam1 + accumlam2 #FIXME c!=0.5
            if z != 0.0
                dbpass += temp/z #FIXME c!=0.5
            end
            #z = logsum(accumlam1, accumlam2)
            #@show exp(temp-z)
            #dbpass = dbpass == -Inf ? temp-z : logsum(dbpass,temp-z) #FIXME c!=0.5
            #dbpass += exp(temp-z)
        end
        temp = dbpass/numpts/2 #FIXME c!=0.5 AND CHECK 4 constant
        beem1 += temp
        beem2 += temp^2
    end
    beem1 /= dblen
    beem2 /= dblen
    return beem1, beem2
end

function bee_e_mc(cls::OBC.BinaryClassifier; dmean=10.0, numpts=100)
    bee_e_mc(cls, (dmean,dmean), numpts=numpts)
end

function bee_e_mc(cls::OBC.BinaryClassifier, dmeans::NTuple{2,Float64}; numpts=100)
    #FIXME only c=0.5
    #Generate points from each class
    dmean1,dmean2 = dmeans
    pts1 = gen_posterior_points(numpts, dmean1, cls.mcmc1.db)
    pts2 = gen_posterior_points(numpts, dmean2, cls.mcmc2.db)
    acc_numpts = size(pts1,2) + size(pts2,2) # requested != generated
    points = hcat(pts1,pts2)
    g0 = calc_g(points, cls.mcmc1.db, 20, dmean=dmean1)
    g1 = calc_g(points, cls.mcmc2.db, 20, dmean=dmean2)
    g0 = clamp(g0, -10_000_000.0, Inf)
    g1 = clamp(g1, -10_000_000.0, Inf)
    z = mapslices(logsum, hcat(g0,g1), 2)
    res = exp(logsum(min(g0,g1) .- z .- log(acc_numpts)))
    return res
end

bee_moments_2sample(cls::OBC.BinaryClassifier; dmean=10.0, numpts=100) = bee_moments_2sample(cls, (dmean,dmean), numpts=numpts)

function bee_moments_2sample(cls::OBC.BinaryClassifier, dmeans::NTuple{2,Float64}; numpts=100)
    #FIXME only c=0.5
    #Generate points from each class
    dmean1,dmean2 = dmeans

    @assert length(cls.mcmc1.db) == length(cls.mcmc2.db)
    dbsize = length(cls.mcmc1.db)
    bee = 0.0
    bee2 = 0.0
    numpts = div(numpts,2)
    for i=1:dbsize

        pts1 = gen_points(numpts*2, dmean1, cls.mcmc1.db[i]) # x and z
        pts2 = gen_points(numpts*2, dmean2, cls.mcmc2.db[i]) # x and z

        allpts = hcat(pts1,pts2)

        g1 = clamp(calc_g(allpts, [cls.mcmc1.db[i]], 20, dmean=dmean1), -1e6, Inf)
        g2 = clamp(calc_g(allpts, [cls.mcmc2.db[i]], 20, dmean=dmean2), -1e6, Inf)

        g1eff = clamp(calc_g(allpts, cls.mcmc1.db, 20, dmean=dmean1), -1e6, Inf)
        g2eff = clamp(calc_g(allpts, cls.mcmc2.db, 20, dmean=dmean2), -1e6, Inf)

        # Could inline and optimize calc_g here

        g1x = g1[[1:numpts,numpts*2+1:numpts*3]]
        g1z = g1[[numpts+1:numpts*2,numpts*3+1:end]]
        g2x = g2[[1:numpts,numpts*2+1:numpts*3]]
        g2z = g2[[numpts+1:numpts*2,numpts*3+1:end]]

        g1xeff = g1eff[[1:numpts,numpts*2+1:numpts*3]]
        g1zeff = g1eff[[numpts+1:numpts*2,numpts*3+1:end]]
        g2xeff = g2eff[[1:numpts,numpts*2+1:numpts*3]]
        g2zeff = g2eff[[numpts+1:numpts*2,numpts*3+1:end]]

        #agree = (g1x .< g2x) .== (g1z .< g2z)
        labelx = (g1xeff .< g2xeff)
        labelz = (g1zeff .< g2zeff)

        zx = Float64[logsum(g1x[j],g2x[j]) for j=1:numpts*2]
        zz = Float64[logsum(g1z[j],g2z[j]) for j=1:numpts*2]

        lpyx = min(g1x,g2x) .- zx
        lpyz = min(g1z,g2z) .- zz

        both1 = g1x .+ g1z .- zx .- zz
        both2 = g2x .+ g2z .- zx .- zz

        cross = min(both1,both2)

        ## BEE calculation
        temp = exp(logsum(vcat(lpyx,lpyz)))
        #@show exp(minimum(lpyx))
        #@show exp(minimum(lpyz))
        @show hist(exp(lpyx))
        @show extrema(exp(lpyx))
        @show temp/numpts/2
        #@show sum(exp(lpyx))/numpts/2
        @show size(lpyx)
        @show numpts
        #@show exp(lpyx[1:3])
        #@show exp(lpyz[1:3])
        bee += temp

        for j=1:numpts
            if labelx[j]==labelz[j]
                bee2 += exp(cross[j])
            end
        end
    end
    bee /= (numpts*4 * dbsize)
    bee2 /= (numpts*2 * dbsize) # FIXME not sure about this 2 constant

    return bee, bee2, bee2 - bee^2
end

function predict(db0, db1, points; dmean=10.0)
    g0 = calc_g(points, db0, 20, dmean=dmean)
    g1 = calc_g(points, db1, 20, dmean=dmean)
    return vec((g0 .- g1) .< 0) * 1.0
end

function error_points(cls::OBC.BinaryClassifier, points, labels; dmean=10.0)
    return sum(abs(predict(cls.mcmc1.db, cls.mcmc2.db, points, dmean=dmean) - labels))/size(labels,1)
end

ext: jl


text:
using EGM96
using Base.Test

tests = ["end_to_end"]
for t in tests
  include("$(t).jl")
end


text:
"""
    Rule(idxs::Tuple{Vararg{Int}}, char_attr::Tuple{Vararg{Char}},
    is_pure::Bool, num_group::Int, num_non_group::Int, occurances::Int)

Struct to store relevant information about a CA.
"""
struct Rule
    idxs::Tuple{Vararg{Int}}
    char_attr::Tuple{Vararg{Char}}
    is_pure::Bool
    num_group::Int
    num_non_group::Int
    occurances::Int
end

"""
    Node(CAs::Array{Rule,1}, taxa_label::String="")

Struct to store a node (is recursive).
"""
struct Node
    CAs::Array{Rule,1}
    children::Array{Node,1}
    taxa_label::String

    function Node(CAs ; label="")
        return new(CAs, Array{Node,1}(), label)
    end
end

"""
    remove_from_tree!(tree_tokens::Vector{String}, taxa_to_remove::Union{Array{String,1},Bool})

Takes a tree in Newick format, removes a specific taxa from the tree.

# Arguments
- `tree_tokens::Vector{String}`: the tree in Newick format, tokenized.
- `taxa_to_remove::Union{Array{String,1},Bool}`: the taxa that will be removed.
"""
function remove_from_tree!(tree_tokens::Vector{String}, taxa_to_remove::Union{Array{String,1},Bool})

    tokens_removed = 0

    old_tree_tokens = deepcopy(tree_tokens)

    # Iterate through tokens in tree
    for (idx,char) in enumerate(old_tree_tokens)

        # If the current character is the taxa to remove
        if (char in taxa_to_remove) || (char == "$taxa_to_remove")

            # Remove preceeding comma and taxa number in tree
            if old_tree_tokens[idx-1] == ","
                deleteat!(tree_tokens, idx-1-tokens_removed)
                deleteat!(tree_tokens, idx-1-tokens_removed)
                tokens_removed += 2

            # Remove succeeding comma and taxa number in tree
            elseif old_tree_tokens[idx+1] == ","
                deleteat!(tree_tokens, idx-tokens_removed)
                deleteat!(tree_tokens, idx-tokens_removed)
                tokens_removed += 2

            # Remove parenthesis, comma, and taxa number in tree
            elseif (old_tree_tokens[idx-1] == "(") && (old_tree_tokens[idx+1] == ")")
                if old_tree_tokens[idx-2] == ","
                    deleteat!(tree_tokens, idx-2-tokens_removed)
                    deleteat!(tree_tokens, idx-2-tokens_removed)
                    deleteat!(tree_tokens, idx-2-tokens_removed)
                    deleteat!(tree_tokens, idx-2-tokens_removed)
                elseif old_tree_tokens[idx+2] == ","
                    deleteat!(tree_tokens, idx-1-tokens_removed)
                    deleteat!(tree_tokens, idx-1-tokens_removed)
                    deleteat!(tree_tokens, idx-1-tokens_removed)
                    deleteat!(tree_tokens, idx-1-tokens_removed)
                end
                tokens_removed += 4
            else
                print("$char \n")
            end
        end
    end

    # Return new tree tokens with taxa removed
    return tree_tokens
end

"""
    get_nodes(tree::String ; taxa_to_remove::Union{Array{String,1},Bool}=false)

Takes a tree in Newick format, returns an internal representation of the tree.

# Arguments
- `tree::String`: the tree in Newick format.
- `taxa_to_remove::Union{Array{String,1},Bool}=false`: the taxa that will be removed (if applicable).
"""
function get_nodes(tree::String ; taxa_to_remove::Union{Array{String,1},Bool}=false)

    # Initialize variables
    num_nodes = 1
    curr_node = 1
    nodes = Array{Dict{String,Any}}(undef, 0)

    # Split tree string into tokens
    tree_tokens = [untokenize(token) for token in collect(tokenize(tree))]

    # Remove taxa from tree if applicable
    if !(taxa_to_remove == false)
        tree_tokens = remove_from_tree!(tree_tokens, taxa_to_remove)
    end

    # Iterate through each token
    for (idx, char) in enumerate(tree_tokens)

        # Start of a new node
        if char == "("
            push!(nodes, Dict{String,Any}("Groups" => 0, "Taxa" => Array{String}(undef, 0), "Descendents" => Array{Int}(undef, 0)))
            paren = 1

            # Iterate through the rest of the tree tokens
            for char in tree_tokens[idx+1:end]

                # Start of a new group
                if char == "("
                    paren += 1
                    curr_node += 1
                    if paren == 2
                        push!(nodes[num_nodes]["Descendents"], curr_node)
                    end

                # End of a group
                elseif char == ")"
                    paren -= 1
                    if paren == 0
                        nodes[num_nodes]["Groups"] += 1
                    end

                # End of a group
                elseif char == "," && paren == 1
                    nodes[num_nodes]["Groups"] += 1

                # Taxa to be added to node
                elseif char != "," && char != " "
                    push!(nodes[num_nodes]["Taxa"], char)
                end

                # End of node
                if paren == 0
                    break
                end
            end
            num_nodes += 1
            curr_node = num_nodes
        end
    end
    return nodes
end

"""
    parse_tree(file_path::String; taxa_to_remove::Union{Array{String,1},Bool}=false)

Takes a Nexus file for a tree, returns an internal representation of that tree (and other relevant information).

# Arguments
- `file_path::String`: file path to the Nexus file.
- `taxa_to_remove::Union{Array{String,1},Bool}=false`: the taxa that will be removed (if applicable).
"""
function parse_tree(file_path::String; taxa_to_remove::Union{Array{String,1},Bool}=false)

    # Initialize variables
    character_labels = Dict{String,String}()
    taxa_labels = Dict{String,String}()
    labels = false
    characters = false
    file = open(file_path)
    nodes = Array{Dict{String,Any}}(undef, 0)
    title_count = 0
    name = "NoName"

    # Iterate through each line of the file
    for line in eachline(file)

        line = String(strip(replace(line, r"\[(.*?)\]" => "")))
        if String(strip(line)) == ""
            continue
        end

        if occursin("TITLE ", line)
            title_count += 1
            if title_count == 2
                name = untokenize(collect(tokenize(line))[3])
            end
        end

        # Get nodes for the tree
        if occursin("TREE ", line) && occursin("=", line)
            clean_line = String(strip(replace(line, r":[-+]?[0-9]*\.?[0-9]*" => "")))

            nodes = get_nodes(clean_line,taxa_to_remove = taxa_to_remove)
            break
        end

        # Match the character sequences with their labels
        if characters
            if line == ";"
                characters = false
            else
                split_line = split(line)
                if !(split_line[1] in taxa_to_remove)
                    character_labels[split_line[1]] = split_line[2]
                end
            end
        end

        # Match the taxa with their labels
        if labels
            label_tokens = collect(tokenize(line))
            temp_labels = Array{String}(undef, 0)
            for token in label_tokens
                if token.kind != collect(tokenize(" "))[1].kind
                    push!(temp_labels, untokenize(token))
                    if length(temp_labels) == 2
                        break
                    end
                end
            end
            if !(temp_labels[2] in taxa_to_remove)
                taxa_labels[temp_labels[1]] = temp_labels[2]
            end
            if occursin(";", line)
                labels = false
            end
        end

        # Differentiate lines
        if occursin("TRANSLATE", line)
            labels = true
        end
        if String(strip(line)) == "MATRIX"
            characters = true
        end
    end
    close(file)
    return nodes, taxa_labels, character_labels, name
end

"""
    add_nodes!(tree::Node,sPu::Array{Dict{String,Any}},sPr::Array{Dict{String,Any}},
    cPu::Array{Dict{String,Any}},cPr::Array{Dict{String,Any}},taxa_labels::Dict{String,String},
    character_labels::Dict{String,String},nodes::Array{Dict{String,Any}},node_num::Int64;complex::Bool=true)

Takes a tree (Node), adds all the CA's from the entire tree into the internal representation.

# Arguments
- `tree::Node`: the tree represented as a Node.
- `sPu::Array{Dict{String,Any}}`: an array of simple pure rules.
- `sPr::Array{Dict{String,Any}}`: an array of simple private rules.
- `cPu::Array{Dict{String,Any}}`: an array of complex pure rules.
- `cPr::Array{Dict{String,Any}}`: an array of complex private rules.
- `taxa_labels::Dict{String,String}`: a mapping of the taxa labels to the character labels.
- `character_labels::Dict{String,String}`: a mapping of the character labels to the corresponding sequences.
- `nodes::Array{Dict{String,Any}}`: an array of nodes.
- `node_num::Int64`: the current node number.
- `complex::Bool=true`: indicates whether complex rules should be calculated
- `protein::Bool=false`: indicates whether dataset is a protein (or nucleotide)
"""
function add_nodes!(tree::Node,sPu::Array{Dict{String,Any}},sPr::Array{Dict{String,Any}},cPu::Array{Dict{String,Any}},cPr::Array{Dict{String,Any}},taxa_labels::Dict{String,String},character_labels::Dict{String,String},nodes::Array{Dict{String,Any}},node_num::Int64;complex::Bool=true,protein::Bool=false)

    # Iterate through each group at a node
    for (group_idx, taxa) in enumerate(get_group_taxa_at_node(nodes, node_num))

        # Initialize varibales
        num_group = length(taxa)
        num_non_group = length(nodes[node_num]["Taxa"]) - num_group
        num_nodes = length(nodes)

        # If terminal node, add taxa label
        if num_group == 1
            push!(tree.children, Node([],label=taxa_labels[sPr[group_idx]["Group_Taxa"][1]]))
        else
            push!(tree.children, Node([]))
        end

        # Add all the Simple Private CA's for the node
        for (idx, info) in sPr[group_idx]["sPr"]
            for CA in info
                push!(tree.children[group_idx].CAs, Rule(tuple(idx), tuple(CA[1]), false, num_group, num_non_group, CA[2]))
            end
        end

        # Add all the Simple Pure CA's for the node
        for (idx, info) in sPu[group_idx]["sPu"]
            push!(tree.children[group_idx].CAs, Rule(tuple(idx), tuple(info[1]), true, num_group, num_non_group, info[2]))
        end

        if complex
            # Add all the Complex Private CA's for the node
            for (idx, info) in cPr[group_idx]["cPr"]
                for CA in info
                    push!(tree.children[group_idx].CAs, Rule(idx, CA[1], false, num_group, num_non_group, CA[2]))
                end
            end

            # Add all the Complex Pure CA's for the node
            for (idx, info) in cPu[group_idx]["cPu"]
                push!(tree.children[group_idx].CAs, Rule(idx, info[1], true, num_group, num_non_group, info[2]))
            end
        end

        # Get descendent's CA's and recur on descendent to add next layer to tree
        if num_group > 1
            new_node_num = nodes[node_num]["Descendents"][group_idx]
            new_sPu, new_sPr = get_sPu_and_sPr(nodes, new_node_num, taxa_labels, character_labels, protein=protein)
            if complex
                new_cPu, new_cPr = get_cPu_and_cPr(nodes, new_node_num, taxa_labels, character_labels, new_sPu, new_sPr)
            else
                new_cPu, new_cPr = cPu, cPr
            end
            #print("$node_num \n")
            add_nodes!(tree.children[group_idx],new_sPu,new_sPr,new_cPu,new_cPr,taxa_labels,character_labels,nodes,new_node_num,complex=complex,protein=protein)
        end
    end
end


text:
__precompile__()

module FortranFiles

export FortranFile, rewind
export RECMRK4B, RECMRK8B
export FString, trimstring, trimlen, trim
export FortranFilesError
export @fread

include("types.jl")
include("exceptions.jl")
include("file.jl")
include("withoutsubrecords.jl")
include("withsubrecords.jl")
include("fixedlengthrecords.jl")
include("string.jl")
include("read.jl")
include("write.jl")
include("macro.jl")

end


text:
using OmegaCore
using Test

@testset "named tuple" begin
  x = (a = 3, b = 4, c = 12)
  x_ = OmegaCore.Util.rmkey(x, Val{:a})
  @test :b ∈ keys(x_)
  @test :c ∈ keys(x_)
  @test :a ∉ keys(x_)
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
module CutsceneNode

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Cutscene Node" => Ahorn.EntityPlacement(
        Maple.CutsceneNode
    )
)

function Ahorn.selection(entity::Maple.CutsceneNode)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle[Ahorn.Rectangle(x - 12, y - 12, 24, 24)]
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.CutsceneNode, room::Maple.Room) = Ahorn.drawImage(ctx, Ahorn.Assets.cutsceneNode, -12, -12)

end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text:
using SpecialFunctions, Nemo, StatsFuns, Memoization
import Nemo.binom, Nemo.gamma

"Session-wide precision"
const prec = 5000 ## Increase for better precision

"Converter from real to arbitrary precision real (Arb)"
const RR = RealField(prec)

"Converter from real to arbitrary precision complex (Acb)"
const CC = ComplexField(prec)

binom(n::Int64, k::Int64, r::Nemo.ArbField) = binom(convert(UInt64, n),convert(UInt64, k), r)

# Making sure the factorial function uses the Nemo version of the gamma and returns an Arb
gamma(x::Int64) = Nemo.gamma(RR(x))

# Useful typed constants
"1 with required arbitrary precision"
const arb_1 = RR(1)
"0 with required arbitrary precision"
const arb_0 = RR(0)

"""
  unsigned_Stirling1(n::Integer, k::Integer)

Computation of unsigned Stirling numbers of the first kind, using Memoization and a recursive formula.

# Examples
```julia-repl
julia> GibbsTypePriors.unsigned_Stirling1(10, 5)
269325
```
"""
@memoize function unsigned_Stirling1(n::Integer, k::Integer)
  # special cases
  if k<0 || n<0
    throw(DomainError())
  end
  if k>n
    return big(0)
  end
  if n==0  # and, by logic, k==0
    return big(1)
  end
  if k==0  # and, by logic, n>0
    return big(0)
  end
  # end of special cases, invoke recursion
  return unsigned_Stirling1(n-1,k-1) + (n-1)*unsigned_Stirling1(n-1,k)
end

@memoize function unsigned_Stirling2(n::Integer,k::Integer)
  # special cases
  if k<0 || n<0
    throw(DomainError())
  end
  if k>n
    return big(0)
  end
  if n==0  # and, by logic, k==0
    return big(1)
  end
  if k==0  # and, by logic, n>0
    return big(0)
  end
  # end of special cases, invoke recursion
  return  (k)*unsigned_Stirling2(n-1,k) + unsigned_Stirling2(n-1,k-1)
end

"""
  has_reasonable_precision(arb_num)

Check whether arb_num has at least roughly the double precision (maybe replace by 64 by 53).

# Examples
```julia-repl
julia> GibbsTypePriors.Cnk_rec(6, 5, 0.5)
0.234375
```
"""
function has_reasonable_precision(arb_num)
    return accuracy_bits(arb_num) ≥ 53
end

import Nemo.risingfac
"""
  risingfac(r, n)

Computes the rising factorial.

# Examples
```julia-repl
julia> GibbsTypePriors.Cnk_rec(6, 5, 0.5)
0.234375
```
"""
function risingfac(r, n)
  if r == 0
    return 0
  else
    return prod(r + i for i in 0:(n-1))
  end
end
function risingfac(r::arb, n)
  if r == 0
    return arb_0
  else
    return prod(r + i for i in 0:(n-1))
  end
end

max_stars_repo_name: konkam/GibbsTypePriors

max_issues_repo_name: konkam/GibbsTypePriors

max_forks_repo_name: konkam/GibbsTypePriors


text:
#####
##### `mean`
#####

_denom(x, dims::Colon) = length(x)
_denom(x, dims::Integer) = size(x, dims)
_denom(x, dims) = mapreduce(i->size(x, i), Base.mul_prod, unique(dims), init=1)

# TODO: We have `mean(f, x; dims)` as of 1.3.0-DEV.36
# https://github.com/JuliaDiff/ChainRules.jl/issues/85
function rrule(::typeof(mean), x::AbstractArray{<:Real}; dims=:)
    y_sum, sum_pullback = rrule(sum, x; dims=dims)
    n = _denom(x, dims)
    function mean_pullback(ȳ)
        ∂x = Thunk() do
            _, ∂sum_x = sum_pullback(ȳ)
            extern(∂sum_x) / n
        end
        return (NO_FIELDS, ∂x)
    end
    return y_sum / n, mean_pullback
end

function rrule(::typeof(mean), f, x::AbstractArray{<:Real})
    y_sum, sum_pullback = rrule(sum, f, x)
    n = _denom(x, :)
    function mean_pullback(ȳ)
        ∂x = Thunk() do
            _, _, ∂sum_x = sum_pullback(ȳ)
            extern(∂sum_x) / n
        end
        return (NO_FIELDS, DoesNotExist(), ∂x)
    end
    return y_sum / n, mean_pullback
end


text:
cd(@__DIR__); include("setups/grid7x3.jl")
pyplot(dpi = 200)

## (a) eigenvectors by nondecreasing eigenvalue ordering
plot(layout = Plots.grid(3, 7))
for i in 1:N
    heatmap!(reshape(𝚽[:, i], (Nx, Ny))', c = :viridis, cbar = false,
                clims = (-0.4,0.4), frame = :none, ratio = 1, ylim = [0, Ny + 1],
                title = latexstring("\\phi_{", i-1, "}"), titlefont = 12,
                subplot = i)
end
plt = current()
savefig(plt, "../figs/grid7x3_evsp_title.png")

## (b) eigenvectors by natural frequency ordering
plot(layout = Plots.grid(3, 7))
for i in 1:N
    k = grid2eig_ind[i]
    heatmap!(reshape(𝚽[:,k], (Nx, Ny))', c = :viridis, cbar = false,
                clims = (-0.4,0.4), frame = :none, ratio = 1, ylim = [0, Ny + 1],
                title = latexstring("\\varphi_{", string(eig2dct[k,1]),
                ",", string(eig2dct[k,2]), "}"), titlefont = 12, subplot = i)
end
plt = current()
savefig(plt, "../figs/grid7x3_dct_title2.png")


text:
#=

@authors: John Stachurski
Date: Thu Aug 21 11:09:30 EST 2014

Provides functions for working with and visualizing scalar ARMA processes.
Ported from Python module quantecon.arma, which was written by Doc-Jin Jang,
Jerry Choi, Thomas Sargent and John Stachurski

References
----------

https://lectures.quantecon.org/jl/arma.html

=#

doc"""
Represents a scalar ARMA(p, q) process

If ``\phi`` and ``\theta`` are scalars, then the model is
understood to be

```math
    X_t = \phi X_{t-1} + \epsilon_t + \theta \epsilon_{t-1}
```

where ``\epsilon_t`` is a white noise process with standard
deviation sigma.

If ``\phi`` and ``\theta`` are arrays or sequences,
then the interpretation is the ARMA(p, q) model

```math
    X_t = \phi_1 X_{t-1} + ... + \phi_p X_{t-p} +
    \epsilon_t + \theta_1 \epsilon_{t-1} + \ldots  +
    \theta_q \epsilon_{t-q}
```

where

* ``\phi = (\phi_1, \phi_2, \ldots , \phi_p)``
* ``\theta = (\theta_1, \theta_2, \ldots , \theta_q)``
* ``\sigma`` is a scalar, the standard deviation of the white noise

##### Fields

 - `phi::Vector` : AR parameters ``\phi_1, \ldots, \phi_p``
 - `theta::Vector` : MA parameters ``\theta_1, \ldots, \theta_q``
 - `p::Integer` : Number of AR coefficients
 - `q::Integer` : Number of MA coefficients
 - `sigma::Real` : Standard deviation of white noise
 - `ma_poly::Vector` : MA polynomial --- filtering representatoin
 - `ar_poly::Vector` : AR polynomial --- filtering representation

##### Examples

```julia
using QuantEcon
phi = 0.5
theta = [0.0, -0.8]
sigma = 1.0
lp = ARMA(phi, theta, sigma)
require(joinpath(dirname(@__FILE__),"..", "examples", "arma_plots.jl"))
quad_plot(lp)
```
"""
mutable struct ARMA
    phi::Vector      # AR parameters phi_1, ..., phi_p
    theta::Vector    # MA parameters theta_1, ..., theta_q
    p::Integer       # Number of AR coefficients
    q::Integer       # Number of MA coefficients
    sigma::Real      # Variance of white noise
    ma_poly::Vector  # MA polynomial --- filtering representatoin
    ar_poly::Vector  # AR polynomial --- filtering representation
end

# constructors to coerce phi/theta to vectors
ARMA(phi::Real, theta::Real, sigma::Real) = ARMA([phi;], [theta;], sigma)
ARMA(phi::Real, theta::Vector, sigma::Real) = ARMA([phi;], theta, sigma)
ARMA(phi::Vector, theta::Real, sigma::Real) = ARMA(phi, [theta;], sigma)

function ARMA(phi::AbstractVector, theta::AbstractVector=[0.0], sigma::Real=1.0)
    # == Record dimensions == #
    p = length(phi)
    q = length(theta)

    # == Build filtering representation of polynomials == #
    ma_poly = [1.0; theta]
    ar_poly = [1.0; -phi]
    return ARMA(phi, theta, p, q, sigma, ma_poly, ar_poly)
end

doc"""
Compute the spectral density function.

The spectral density is the discrete time Fourier transform of the
autocovariance function. In particular,

```math
    f(w) = \sum_k \gamma(k) \exp(-ikw)
```

where ``\gamma`` is the autocovariance function and the sum is over
the set of all integers.

##### Arguments

- `arma::ARMA`: Instance of `ARMA` type
- `;two_pi::Bool(true)`: Compute the spectral density function over ``[0, \pi]``
  if false and ``[0, 2 \pi]`` otherwise.
- `;res(1200)` : If `res` is a scalar then the spectral density is computed at
  `res` frequencies evenly spaced around the unit circle, but if `res` is an array
  then the function computes the response at the frequencies given by the array

##### Returns
- `w::Vector{Float64}`: The normalized frequencies at which h was computed, in
  radians/sample
- `spect::Vector{Float64}` : The frequency response
"""
function spectral_density(arma::ARMA; res=1200, two_pi::Bool=true)
    # Compute the spectral density associated with ARMA process arma
    wmax = two_pi ? 2pi : pi
    w = linspace(0, wmax, res)
    tf = TFFilter(reverse(arma.ma_poly), reverse(arma.ar_poly))
    h = freqz(tf, w)
    spect = arma.sigma.^2 .* abs.(h).^2
    return w, spect
end

"""
Compute the autocovariance function from the ARMA parameters
over the integers range(`num_autocov`) using the spectral density
and the inverse Fourier transform.

##### Arguments

- `arma::ARMA`: Instance of `ARMA` type
- `;num_autocov::Integer(16)` : The number of autocovariances to calculate

"""
function autocovariance(arma::ARMA; num_autocov::Integer=16)
    # Compute the autocovariance function associated with ARMA process arma
    # Computation is via the spectral density and inverse FFT
    (w, spect) = spectral_density(arma)
    acov = real(Base.ifft(spect))
    # num_autocov should be <= len(acov) / 2
    return acov[1:num_autocov]
end

doc"""
Get the impulse response corresponding to our model.

##### Arguments

- `arma::ARMA`: Instance of `ARMA` type
- `;impulse_length::Integer(30)`: Length of horizon for calcluating impulse reponse. Must be at least as long as the `p` fields of `arma`

##### Returns

- `psi::Vector{Float64}`: `psi[j]` is the response at lag j of the impulse
  response. We take `psi[1]` as unity.

"""
function impulse_response(arma::ARMA; impulse_length=30)
    # Compute the impulse response function associated with ARMA process arma
    err_msg = "Impulse length must be greater than number of AR coefficients"
    @assert impulse_length >= arma.p err_msg
    # == Pad theta with zeros at the end == #
    theta = [arma.theta; zeros(impulse_length - arma.q)]
    psi_zero = 1.0
    psi = Array{Float64}(impulse_length)
    for j = 1:impulse_length
        psi[j] = theta[j]
        for i = 1:min(j, arma.p)
            psi[j] += arma.phi[i] * (j-i > 0 ? psi[j-i] : psi_zero)
        end
    end
    return [psi_zero; psi[1:end-1]]
end

"""
Compute a simulated sample path assuming Gaussian shocks.

##### Arguments

- `arma::ARMA`: Instance of `ARMA` type
- `;ts_length::Integer(90)`: Length of simulation
- `;impulse_length::Integer(30)`: Horizon for calculating impulse response
  (see also docstring for `impulse_response`)

##### Returns

- `X::Vector{Float64}`: Simulation of the ARMA model `arma`

"""
function simulation(arma::ARMA; ts_length=90, impulse_length=30)
    # Simulate the ARMA process arma assuing Gaussian shocks
    J = impulse_length
    T = ts_length
    psi = impulse_response(arma, impulse_length=impulse_length)
    epsilon = arma.sigma * randn(T + J)
    X = Array{Float64}(T)
    for t=1:T
        X[t] = dot(epsilon[t:J+t-1], psi)
    end
    return X
end


text:
using AdalmPluto
using Documenter

if haskey(ENV, "DOCSARGS")
    for arg in split(ENV["DOCSARGS"])
        (arg in ARGS) || push!(ARGS, arg)
    end
end

makedocs(;
    modules=[AdalmPluto, libIIO_jl],
    authors="JuliaTelecom and contributors",
    repo="https://github.com/JuliaTelecom/AdalmPluto.jl/blob/{commit}{path}#L{line}",
    sitename="AdalmPluto.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "AdalmPluto" => "adalmpluto.md",
        "libIIO_jl" => Any[
            "libiio/scan.md",
            "libiio/toplevel.md",
            "libiio/context.md",
            "libiio/device.md",
            "libiio/channel.md",
            "libiio/buffer.md",
            "libiio/debug.md",
            "libiio/structures.md",
            "libiio/helpers.md"
        ],
    ],
)

deploydocs(
    repo = "github.com/JuliaTelecom/AdalmPluto.jl.git",
    target = "build",
    push_preview = true,
    devbranch = "main",
)


text:
import JSON
export BrainFlowModelParams

@enum BrainFlowMetrics begin

    MINDFULNESS = 0
    RESTFULNESS = 1
    USER_DEFINED = 2

end

MetricType = Union{BrainFlowMetrics, Integer}

@enum BrainFlowClassifiers begin

    DEFAULT_CLASSIFIER = 0
    DYN_LIB_CLASSIFIER = 1
    ONNX_CLASSIFIER = 2

end

ClassifierType = Union{BrainFlowClassifiers, Integer}

mutable struct BrainFlowModelParams
    metric::MetricType
    classifier::ClassifierType
    file::String
    other_info::String
    output_name::String
    max_array_size::Int32

    function BrainFlowModelParams(metric_::MetricType, classifier_::ClassifierType)
        new(metric_, classifier_, "", "", "", 8192)
    end

end

function JSON.json(params::BrainFlowModelParams)
    d = Dict(
        "metric" => Integer(params.metric),
        "classifier" => Integer(params.classifier),
        "file" => params.file,
        "other_info" => params.other_info,
        "output_name" => params.output_name,
        "max_array_size" => params.max_array_size,
        )
    return JSON.json(d)
end

@brainflow_rethrow function prepare(params::BrainFlowModelParams)
    input_json = JSON.json(params)
    ccall((:prepare, ML_MODULE_INTERFACE), Cint, (Ptr{UInt8},), input_json)
    return
end

@brainflow_rethrow function release(params::BrainFlowModelParams)
    input_json = JSON.json(params)
    ccall((:release, ML_MODULE_INTERFACE), Cint, (Ptr{UInt8},), input_json)
    return
end

@brainflow_rethrow function predict(data, params::BrainFlowModelParams)
    input_json = JSON.json(params)
    val = Vector{Float64}(undef, params.max_array_size)
    val_len = Vector{Float64}(undef, 1)
    ccall((:predict, ML_MODULE_INTERFACE), Cint, (Ptr{Float64}, Cint, Ptr{Float64}, Ptr{Cint}, Ptr{UInt8}),
        data, length(data), val, val_len, input_json)
    value = val[1:val_len[1]]
    return value
end

@brainflow_rethrow function release_all()
    ccall((:release_all, ML_MODULE_INTERFACE), Cint, ())
end


text:
# 2520 is the smallest number that can be divided by each of the numbers from 1
# to 10 without any remainder.
#
# What is the smallest positive number that is evenly divisible by all of the
# numbers from 1 to 20?

using ProjectEulerSolutions

# Find the number of prime factors for each number between 2 and n.  Store the
# count of each prime factor, and then take the max number of that prime factor
# for each number.  Multiply all the factors (exponentiated by their max number)
# to ensure the result is evenly divisible.
function p005solution(n::Integer=2)::Integer
    counts = zeros(Integer, n)
    for i in 2:n
        subcounts = Dict{Integer, Integer}()
        for j in primefactors(i)
            subcounts[j] = get(subcounts, j, 0) + 1
        end
        for k in subcounts
            @inbounds counts[k[1]] = max(counts[k[1]], k[2])
        end
    end

    result = 1
    for i in 2:n
        @inbounds if counts[i] > 0
            result = result * (i ^ counts[i])
        end
    end

    return result
end

p005 = Problems.Problem(p005solution)

Problems.benchmark(p005, 20)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
@testset "PolyharmonicSpline" begin
    x = RandomVariable.(Uniform(-π, π), [:x1, :x2, :x3])
    a = Parameter(7, :a)
    b = Parameter(0.1, :b)

    inputs = [x; a; b]

    ishigami = Model(
        df ->
            sin.(df.x1) .+ df.a .* sin.(df.x2) .^ 2 .+ df.b .* (df.x3 .^ 4) .* sin.(df.x1),
        :y,
    )

    Random.seed!(8128)

    data = sample(inputs, 100)
    evaluate!(ishigami, data)

    Random.seed!()
    spline = PolyharmonicSpline(data, 2, :y)

    splinedata = select(data, Not(:y))
    evaluate!(spline, splinedata)

    mse = mean((data.y .- splinedata.y) .^ 2)
    @test isapprox(mse, 0; atol=eps(Float64))
end

max_stars_repo_name: FriesischScott/julia-uq

max_issues_repo_name: FriesischScott/julia-uq

max_forks_repo_name: FriesischScott/julia-uq


text:
using JetPack, Jets, Test

@testset "restriction, real to real, 1D" begin
    A = JopRestriction(JetSpace(Float64,512), [1:4:512;])
    lhs, rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
    @test lhs ≈ rhs
    m = rand(domain(A))
    @test A*m ≈ m[1:4:512]
    d = m[1:4:512]
    m = zeros(domain(A))
    m[1:4:512] .= d
    @test A'*d ≈ m
end

@testset "restriction, real to real, 2D" begin
    indices = [ (i,j) for i=[1,4], j=[1,4] ]
    A = JopRestriction(JetSpace(Float64,4,4), vec(indices))
    m = [11.0 12.0 13.0 14.0 ; 21.0 22.0 23.0 24.0 ; 31.0 32.0 33.0 34.0 ; 41.0 42.0 43.0 44.0]
    @test A*m ≈ [11.0, 41.0, 14, 44.0]
    lhs, rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
    @test lhs ≈ rhs
end


text:
"""
    struct GOCNullImpedance_withtransformer <: AbstractLinkLabel

Structure descendant of AbstractLinkLabel

# Fields
- `link::Link`
- `id::String`
- `susceptance::Float64`: coefficient bc
- `transfo_ratio::Float64`: ratio τ
- `transfo_phase::Float64`: phase θ
- `power_magnitude_max::Float64`: Smax

"""
struct GOCNullImpedance_withtransformer <: AbstractLinkLabel
    orig_id::Int64
    dest_id::Int64
    id::Any
    susceptance::Float64 #bc
    transfo_ratio::Float64 #τ
    transfo_phase::Float64 #θ
    power_magnitude_max::Float64 #Smax
end

"""
    create_vars!(element::T, link::Link, elemid::String, elem_formulation::Symbol, link_vars::SortedDict{String, Variable}, scenario::String) where T <: GOCNullImpedance_withtransformer

Create voltage variables and power variables for origin and destination of `link` in `link_vars`.
Return nothing.
"""
function create_vars!(element::T, link::Link, elemid::String, elem_formulation::Symbol, link_vars::SortedDict{String, Variable}, scenario::String) where T <: GOCNullImpedance_withtransformer
    origid, destid = link.orig, link.dest
    link_vars["Volt_orig"] = Variable(variable_name("VOLT", origid, "", scenario), Complex)
    link_vars["Volt_dest"] = Variable(variable_name("VOLT", destid, "", scenario), Complex)

    link_vars[elemid*"_Sorig"] = Variable(variable_name("Sorig", origid, "", scenario), Complex)
    # link_vars[elemid*"_Sdest"] = Variable(variable_name("Sdest", destid, "", scenario), Complex)
    #elem_formulation modified ?
    return
end

"""
    Sorig(element::T, link::Link, elemid::String, elem_formulation::Symbol, link_vars::SortedDict{String, Variable}) where T<:GOCNullImpedance_withtransformer

Return power variable Sorig * baseMVA.
"""
function Sorig(element::T, link::Link, elemid::String, elem_formulation::Symbol, link_vars::SortedDict{String, Variable}) where T<:GOCNullImpedance_withtransformer
    return get_baseMVA(link.orig)*link_vars[elemid*"_Sorig"]
end

"""
    Sdest(element::T, link::Link, elemid::String, elem_formulation::Symbol, link_vars::SortedDict{String, Variable}) where T<:GOCNullImpedance_withtransformer

Return power variable Sdest * baseMVA.
"""
function Sdest(element::T, link::Link, elemid::String, elem_formulation::Symbol, link_vars::SortedDict{String, Variable}) where T<:GOCNullImpedance_withtransformer
   return get_baseMVA(link.dest)*(link_vars[elemid*"_Sorig"]- im * element.susceptance * abs2(link_vars["Volt_dest"]))
end

## 3. Constraints creation
function constraint(element::T, link::Link, elemid::String, elem_formulation::Symbol, link_vars::SortedDict{String, Variable}, scenario::String, OPFpbs::OPFProblems) where T <:GOCNullImpedance_withtransformer
    cstrs = SortedDict{String, Constraint}()

    #Smax constraints
    Sor = Sorig(element, link, elemid, elem_formulation, link_vars)
    Sde = Sdest(element, link, elemid, elem_formulation, link_vars)
    Smax = element.power_magnitude_max

    cstrs[get_Smax_orig_cstrname()] = (abs2(real(Sor)) + abs2(imag(Sor))) << Smax^2
    cstrs[get_Smax_dest_cstrname()] = (abs2(real(Sde)) + abs2(imag(Sde))) << Smax^2

    cstrs[get_Smax_orig_cstrname()].precond = :sqrt
    cstrs[get_Smax_dest_cstrname()].precond = :sqrt

    ## constraint linking Vorig and Vdest
    τ = element.transfo_ratio
    θ = element.transfo_phase
    cstrs[get_NullImpVolt_cstrname()] = Polynomial(link_vars["Volt_orig"] - link_vars["Volt_dest"] * 1/τ * exp(-im*θ))==0

    return cstrs
end


text:
module StatsAPI

#    pairwise(f, x[, y])
#
# Return a matrix holding the result of applying `f` to all possible pairs
# of entries in iterators `x` and `y`, and return it.
#
# This generic function is owned by StatsBase.jl, which is the sole provider
# of the default definition.
function pairwise end

#    pairwise!(f, dest::AbstractMatrix, x[, y])
#
# Store in matrix `dest` the result of applying `f` to all possible pairs
# of entries in iterators `x` and `y`, and return it.
#
# This generic function is owned by StatsBase.jl, which is the sole provider
# of the default definition.
function pairwise! end

end # module


text:
using FractionalDiffEq

K = 1
α = 0.5
dx = pi/2
dt = 0.5
n = 2
xStart = 0;
xEnd = pi
x = collect(0:dx:xEnd)
t = collect(0:dt:n)
u0t = 0
uendt = 0
u0(x) = sin(x)

U=solve(α, dx, dt, xStart, xEnd, n, K, u0t, uendt, u0, FiniteDiffIm())

using Plots

plotlyjs()
plot(x, t, U, st=:surface)


text:
        # Error maps (maps err to (NAME, Description))

const error_codes = @compat Dict(
 0   =>	("GNUTLS_E_SUCCESS","Success."),
-3	 => ("GNUTLS_E_UNKNOWN_COMPRESSION_ALGORITHM","Could not negotiate a supported compression method."),
-6	 => ("GNUTLS_E_UNKNOWN_CIPHER_TYPE","The cipher type is unsupported."),
-7	 => ("GNUTLS_E_LARGE_PACKET","The transmitted packet is too large (EMSGSIZE)."),
-8	 => ("GNUTLS_E_UNSUPPORTED_VERSION_PACKET","A record packet with illegal version was received."),
-9	 => ("GNUTLS_E_UNEXPECTED_PACKET_LENGTH","A TLS packet with unexpected length was received."),
-10	 => ("GNUTLS_E_INVALID_SESSION","The specified session has been invalidated for some reason."),
-12	 => ("GNUTLS_E_FATAL_ALERT_RECEIVED","A TLS fatal alert has been received."),
-15	 => ("GNUTLS_E_UNEXPECTED_PACKET","An unexpected TLS packet was received."),
-16	 => ("GNUTLS_E_WARNING_ALERT_RECEIVED","A TLS warning alert has been received."),
-18	 => ("GNUTLS_E_ERROR_IN_FINISHED_PACKET","An error was encountered at the TLS Finished packet calculation."),
-19	 => ("GNUTLS_E_UNEXPECTED_HANDSHAKE_PACKET","An unexpected TLS handshake packet was received."),
-21	 => ("GNUTLS_E_UNKNOWN_CIPHER_SUITE","Could not negotiate a supported cipher suite."),
-22	 => ("GNUTLS_E_UNWANTED_ALGORITHM","An algorithm that is not enabled was negotiated."),
-23	 => ("GNUTLS_E_MPI_SCAN_FAILED","The scanning of a large integer has failed."),
-24	 => ("GNUTLS_E_DECRYPTION_FAILED","Decryption has failed."),
-25	 => ("GNUTLS_E_MEMORY_ERROR","Internal error in memory allocation."),
-26	 => ("GNUTLS_E_DECOMPRESSION_FAILED","Decompression of the TLS record packet has failed."),
-27	 => ("GNUTLS_E_COMPRESSION_FAILED","Compression of the TLS record packet has failed."),
-28	 => ("GNUTLS_E_AGAIN","Resource temporarily unavailable, try again."),
-29	 => ("GNUTLS_E_EXPIRED","The requested session has expired."),
-30	 => ("GNUTLS_E_DB_ERROR","Error in Database backend."),
-31	 => ("GNUTLS_E_SRP_PWD_ERROR","Error in password file."),
-32	 => ("GNUTLS_E_INSUFFICIENT_CREDENTIALS","Insufficient credentials for that request."),
-33	 => ("GNUTLS_E_HASH_FAILED","Hashing has failed."),
-34	 => ("GNUTLS_E_BASE64_DECODING_ERROR","Base64 decoding error."),
-35	 => ("GNUTLS_E_MPI_PRINT_FAILED","Could not export a large integer."),
-37	 => ("GNUTLS_E_REHANDSHAKE","Rehandshake was requested by the peer."),
-38	 => ("GNUTLS_E_GOT_APPLICATION_DATA","TLS Application data were received, while expecting handshake data."),
-39	 => ("GNUTLS_E_RECORD_LIMIT_REACHED","The upper limit of record packet sequence numbers has been reached. Wow!"),
-40	 => ("GNUTLS_E_ENCRYPTION_FAILED","Encryption has failed."),
-43	 => ("GNUTLS_E_CERTIFICATE_ERROR","Error in the certificate."),
-44	 => ("GNUTLS_E_PK_ENCRYPTION_FAILED","Public key encryption has failed."),
-45	 => ("GNUTLS_E_PK_DECRYPTION_FAILED","Public key decryption has failed."),
-46	 => ("GNUTLS_E_PK_SIGN_FAILED","Public key signing has failed."),
-47	 => ("GNUTLS_E_X509_UNSUPPORTED_CRITICAL_EXTENSION","Unsupported critical extension in X.509 certificate."),
-48	 => ("GNUTLS_E_KEY_USAGE_VIOLATION","Key usage violation in certificate has been detected."),
-49	 => ("GNUTLS_E_NO_CERTIFICATE_FOUND","No certificate was found."),
-50	 => ("GNUTLS_E_INVALID_REQUEST","The request is invalid."),
-51	 => ("GNUTLS_E_SHORT_MEMORY_BUFFER","The given memory buffer is too short to hold parameters."),
-52	 => ("GNUTLS_E_INTERRUPTED","Function was interrupted."),
-53	 => ("GNUTLS_E_PUSH_ERROR","Error in the push function."),
-54	 => ("GNUTLS_E_PULL_ERROR","Error in the pull function."),
-55	 => ("GNUTLS_E_RECEIVED_ILLEGAL_PARAMETER","An illegal parameter has been received."),
-56	 => ("GNUTLS_E_REQUESTED_DATA_NOT_AVAILABLE","The requested data were not available."),
-57	 => ("GNUTLS_E_PKCS1_WRONG_PAD","Wrong padding in PKCS1 packet."),
-58	 => ("GNUTLS_E_RECEIVED_ILLEGAL_EXTENSION","An illegal TLS extension was received."),
-59	 => ("GNUTLS_E_INTERNAL_ERROR","GnuTLS internal error."),
-60	 => ("GNUTLS_E_CERTIFICATE_KEY_MISMATCH","The certificate and the given key do not match."),
-61	 => ("GNUTLS_E_UNSUPPORTED_CERTIFICATE_TYPE","The certificate type is not supported."),
-62	 => ("GNUTLS_E_X509_UNKNOWN_SAN","Unknown Subject Alternative name in X.509 certificate."),
-63	 => ("GNUTLS_E_DH_PRIME_UNACCEPTABLE","The Diffie-Hellman prime sent by the server is not acceptable (not long enough)."),
-64	 => ("GNUTLS_E_FILE_ERROR","Error while reading file."),
-67	 => ("GNUTLS_E_ASN1_ELEMENT_NOT_FOUND","ASN1 parser: Element was not found."),
-68	 => ("GNUTLS_E_ASN1_IDENTIFIER_NOT_FOUND","ASN1 parser: Identifier was not found"),
-69	 => ("GNUTLS_E_ASN1_DER_ERROR","ASN1 parser: Error in DER parsing."),
-70	 => ("GNUTLS_E_ASN1_VALUE_NOT_FOUND","ASN1 parser: Value was not found."),
-71	 => ("GNUTLS_E_ASN1_GENERIC_ERROR","ASN1 parser: Generic parsing error."),
-72	 => ("GNUTLS_E_ASN1_VALUE_NOT_VALID","ASN1 parser: Value is not valid."),
-73	 => ("GNUTLS_E_ASN1_TAG_ERROR","ASN1 parser: Error in TAG."),
-74	 => ("GNUTLS_E_ASN1_TAG_IMPLICIT","ASN1 parser: error in implicit tag"),
-75	 => ("GNUTLS_E_ASN1_TYPE_ANY_ERROR","ASN1 parser: Error in type ’ANY’."),
-76	 => ("GNUTLS_E_ASN1_SYNTAX_ERROR","ASN1 parser: Syntax error."),
-77	 => ("GNUTLS_E_ASN1_DER_OVERFLOW","ASN1 parser: Overflow in DER parsing."),
-78	 => ("GNUTLS_E_TOO_MANY_EMPTY_PACKETS","Too many empty record packets have been received."),
-79	 => ("GNUTLS_E_OPENPGP_UID_REVOKED","The OpenPGP User ID is revoked."),
-80	 => ("GNUTLS_E_UNKNOWN_PK_ALGORITHM","An unknown public key algorithm was encountered."),
-81	 => ("GNUTLS_E_TOO_MANY_HANDSHAKE_PACKETS","Too many handshake packets have been received."),
-84	 => ("GNUTLS_E_NO_TEMPORARY_RSA_PARAMS","No temporary RSA parameters were found."),
-86	 => ("GNUTLS_E_NO_COMPRESSION_ALGORITHMS","No supported compression algorithms have been found."),
-87	 => ("GNUTLS_E_NO_CIPHER_SUITES","No supported cipher suites have been found."),
-88	 => ("GNUTLS_E_OPENPGP_GETKEY_FAILED","Could not get OpenPGP key."),
-89	 => ("GNUTLS_E_PK_SIG_VERIFY_FAILED","Public key signature verification has failed."),
-90	 => ("GNUTLS_E_ILLEGAL_SRP_USERNAME","The SRP username supplied is illegal."),
-91	 => ("GNUTLS_E_SRP_PWD_PARSING_ERROR","Parsing error in password file."),
-93	 => ("GNUTLS_E_NO_TEMPORARY_DH_PARAMS","No temporary DH parameters were found."),
-94	 => ("GNUTLS_E_OPENPGP_FINGERPRINT_UNSUPPORTED","The OpenPGP fingerprint is not supported."),
-95	 => ("GNUTLS_E_X509_UNSUPPORTED_ATTRIBUTE","The certificate has unsupported attributes."),
-96	 => ("GNUTLS_E_UNKNOWN_HASH_ALGORITHM","The hash algorithm is unknown."),
-97	 => ("GNUTLS_E_UNKNOWN_PKCS_CONTENT_TYPE","The PKCS structure’s content type is unknown."),
-98	 => ("GNUTLS_E_UNKNOWN_PKCS_BAG_TYPE","The PKCS structure’s bag type is unknown."),
-99	 => ("GNUTLS_E_INVALID_PASSWORD","The given password contains invalid characters."),
-100	 => ("GNUTLS_E_MAC_VERIFY_FAILED","The Message Authentication Code verification failed."),
-101	 => ("GNUTLS_E_CONSTRAINT_ERROR","Some constraint limits were reached."),
-102	 => ("GNUTLS_E_WARNING_IA_IPHF_RECEIVED","Received a TLS/IA Intermediate Phase Finished message"),
-103	 => ("GNUTLS_E_WARNING_IA_FPHF_RECEIVED","Received a TLS/IA Final Phase Finished message"),
-104	 => ("GNUTLS_E_IA_VERIFY_FAILED","Verifying TLS/IA phase checksum failed"),
-105	 => ("GNUTLS_E_UNKNOWN_ALGORITHM","The specified algorithm or protocol is unknown."),
-106	 => ("GNUTLS_E_UNSUPPORTED_SIGNATURE_ALGORITHM","The signature algorithm is not supported."),
-107	 => ("GNUTLS_E_SAFE_RENEGOTIATION_FAILED","Safe renegotiation failed."),
-108	 => ("GNUTLS_E_UNSAFE_RENEGOTIATION_DENIED","Unsafe renegotiation denied."),
-109	 => ("GNUTLS_E_UNKNOWN_SRP_USERNAME","The SRP username supplied is unknown."),
-110	 => ("GNUTLS_E_PREMATURE_TERMINATION","The TLS connection was non-properly terminated."),
-201	 => ("GNUTLS_E_BASE64_ENCODING_ERROR","Base64 encoding error."),
-202	 => ("GNUTLS_E_INCOMPATIBLE_GCRYPT_LIBRARY","The crypto library version is too old."),
-203	 => ("GNUTLS_E_INCOMPATIBLE_LIBTASN1_LIBRARY","The tasn1 library version is too old."),
-204	 => ("GNUTLS_E_OPENPGP_KEYRING_ERROR","Error loading the keyring."),
-205	 => ("GNUTLS_E_X509_UNSUPPORTED_OID","The OID is not supported."),
-206	 => ("GNUTLS_E_RANDOM_FAILED","Failed to acquire random data."),
-207	 => ("GNUTLS_E_BASE64_UNEXPECTED_HEADER_ERROR","Base64 unexpected header error."),
-208	 => ("GNUTLS_E_OPENPGP_SUBKEY_ERROR","Could not find OpenPGP subkey."),
-209	 => ("GNUTLS_E_CRYPTO_ALREADY_REGISTERED","There is already a crypto algorithm with lower priority."),
-210	 => ("GNUTLS_E_HANDSHAKE_TOO_LARGE","The handshake data size is too large."),
-211	 => ("GNUTLS_E_CRYPTODEV_IOCTL_ERROR","Error interfacing with /dev/crypto"),
-212	 => ("GNUTLS_E_CRYPTODEV_DEVICE_ERROR","Error opening /dev/crypto"),
-213	 => ("GNUTLS_E_CHANNEL_BINDING_NOT_AVAILABLE","Channel binding data not available"),
-214	 => ("GNUTLS_E_BAD_COOKIE","The cookie was bad."),
-215	 => ("GNUTLS_E_OPENPGP_PREFERRED_KEY_ERROR","The OpenPGP key has not a preferred key set."),
-216	 => ("GNUTLS_E_INCOMPAT_DSA_KEY_WITH_TLS_PROTOCOL","The given DSA key is incompatible with the selected TLS protocol."),
-292	 => ("GNUTLS_E_HEARTBEAT_PONG_RECEIVED","A heartbeat pong message was received."),
-293	 => ("GNUTLS_E_HEARTBEAT_PING_RECEIVED","A heartbeat ping message was received."),
-300	 => ("GNUTLS_E_PKCS11_ERROR","PKCS #11 error."),
-301	 => ("GNUTLS_E_PKCS11_LOAD_ERROR","PKCS #11 initialization error."),
-302	 => ("GNUTLS_E_PARSING_ERROR","Error in parsing."),
-303	 => ("GNUTLS_E_PKCS11_PIN_ERROR","Error in provided PIN."),
-305	 => ("GNUTLS_E_PKCS11_SLOT_ERROR","PKCS #11 error in slot"),
-306	 => ("GNUTLS_E_LOCKING_ERROR","Thread locking error"),
-307	 => ("GNUTLS_E_PKCS11_ATTRIBUTE_ERROR","PKCS #11 error in attribute"),
-308	 => ("GNUTLS_E_PKCS11_DEVICE_ERROR","PKCS #11 error in device"),
-309	 => ("GNUTLS_E_PKCS11_DATA_ERROR","PKCS #11 error in data"),
-310	 => ("GNUTLS_E_PKCS11_UNSUPPORTED_FEATURE_ERROR","PKCS #11 unsupported feature"),
-311	 => ("GNUTLS_E_PKCS11_KEY_ERROR","PKCS #11 error in key"),
-312	 => ("GNUTLS_E_PKCS11_PIN_EXPIRED","PKCS #11 PIN expired"),
-313	 => ("GNUTLS_E_PKCS11_PIN_LOCKED","PKCS #11 PIN locked"),
-314	 => ("GNUTLS_E_PKCS11_SESSION_ERROR","PKCS #11 error in session"),
-315	 => ("GNUTLS_E_PKCS11_SIGNATURE_ERROR","PKCS #11 error in signature"),
-316	 => ("GNUTLS_E_PKCS11_TOKEN_ERROR","PKCS #11 error in token"),
-317	 => ("GNUTLS_E_PKCS11_USER_ERROR","PKCS #11 user error"),
-318	 => ("GNUTLS_E_CRYPTO_INIT_FAILED","The initialization of crypto backend has failed."),
-319	 => ("GNUTLS_E_TIMEDOUT","The operation timed out"),
-320	 => ("GNUTLS_E_USER_ERROR","The operation was cancelled due to user error"),
-321	 => ("GNUTLS_E_ECC_NO_SUPPORTED_CURVES","No supported ECC curves were found"),
-322	 => ("GNUTLS_E_ECC_UNSUPPORTED_CURVE","The curve is unsupported"),
-323	 => ("GNUTLS_E_PKCS11_REQUESTED_OBJECT_NOT_AVAILBLE","The requested PKCS #11 object is not available"),
-324	 => ("GNUTLS_E_CERTIFICATE_LIST_UNSORTED","The provided X.509 certificate list is not sorted (in subject to issuer order)"),
-325	 => ("GNUTLS_E_ILLEGAL_PARAMETER","An illegal parameter was found."),
-326	 => ("GNUTLS_E_NO_PRIORITIES_WERE_SET","No or insufficient priorities were set."),
-327	 => ("GNUTLS_E_X509_UNSUPPORTED_EXTENSION","Unsupported extension in X.509 certificate."),
-328	 => ("GNUTLS_E_SESSION_EOF","Peer has terminated the connection"),
-329	 => ("GNUTLS_E_TPM_ERROR","TPM error."),
-330	 => ("GNUTLS_E_TPM_KEY_PASSWORD_ERROR","Error in provided password for key to be loaded in TPM."),
-331	 => ("GNUTLS_E_TPM_SRK_PASSWORD_ERROR","Error in provided SRK password for TPM."),
-332	 => ("GNUTLS_E_TPM_SESSION_ERROR","Cannot initialize a session with the TPM."),
-333	 => ("GNUTLS_E_TPM_KEY_NOT_FOUND","TPM key was not found in persistent storage."),
-334	 => ("GNUTLS_E_TPM_UNINITIALIZED","TPM is not initialized."),
-340	 => ("GNUTLS_E_NO_CERTIFICATE_STATUS","There is no certificate status (OCSP)."),
-341	 => ("GNUTLS_E_OCSP_RESPONSE_ERROR","The OCSP response is invalid"),
-342	 => ("GNUTLS_E_RANDOM_DEVICE_ERROR","Error in the system’s randomness device."),
-343	 => ("GNUTLS_E_AUTH_ERROR","Could not authenticate peer."),
-344	 => ("GNUTLS_E_NO_APPLICATION_PROTOCOL","No common application protocol could be negotiated."),
-1250    => ("GNUTLS_E_UNIMPLEMENTED_FEATURE", "Unimplemented Feature")
)


text:
module LineCount
    export linecount

    function linecount(filename)
        count = 0
        for line in eachline(filename)
            count += 1
        end
        count
    end
end

max_stars_repo_name: planetnoob/LearnJulia

max_issues_repo_name: planetnoob/LearnJulia

max_forks_repo_name: planetnoob/LearnJulia


text:
using DiffEqDevTools, Test

@test stability_region(constructDormandPrince6(), initial_guess=-3.5) ≈ -3.95413 rtol=1e-3
@test stability_region(constructTsitourasPapakostas6(), initial_guess=-3.5) ≈ -3.95413 rtol=1e-3
@test stability_region(constructRadauIIA5(), initial_guess=12.) ≈ 11.84 rtol=1e-2


text:
using LiBRA, Plots

#---------- Cell Definition -----------------#
Cell = Construct("LG_M50") #Alternative "Doyle_94"
Init_SOC = 0.75

#---------- DRA Loop -----------------#
@inline function DRA_loop(Cell, Init_SOC)
    A = B = C = D = tuple()

    for Temp in 5.0:10.0:45.0
        #Arrhenius
        Cell.Const.T = 273.15+Temp
        Arr_Factor = (1/Cell.Const.T_ref-1/Cell.Const.T)/R

        #Set Cell Constants
        Cell.Const.SOC = Init_SOC
        Cell.Const.κ = Cell.Const.κf(Cell.Const.ce0)*exp(Cell.Const.Ea_κ*Arr_Factor)
        Cell.RA.Nfft = Cell.RA.Nfft!(Cell.RA.Fs, Cell.RA.Tlen)
        Cell.RA.f = Cell.RA.f!(Cell.RA.Nfft)
        Cell.RA.s = Cell.RA.s!(Cell.RA.Fs,Cell.RA.Nfft,Cell.RA.f)

        #DRA
        A_DRA, B_DRA, C_DRA, D_DRA = DRA(Cell,Cell.RA.s,Cell.RA.f)

        #Flatten output into tuples
        A = flatten_(A,A_DRA)
        B = flatten_(B,B_DRA)
        C = flatten_(C,C_DRA)
        D = flatten_(D,D_DRA)
    end
return A, B, C, D
end

#---------- Generate Model -----------------#
A, B, C, D = DRA_loop(Cell, Init_SOC)

#---------- Sim Loop -----------------#
function Sim_loop(Cell, Init_SOC)
    CellV = Ce = jNeg = jPos = RtotNeg = RtotPos = η0 = ηL = η_neg = η_pos = ϕ_ẽ1 = ϕ_ẽ2 = Uocp_Neg = Uocp_Pos = ϕ_e = Cse_Neg = Cse_Pos = tuple()

    #Set Experiment
    i = Int64(1/Cell.RA.SamplingT) #Sampling Frequency
    Iapp = [ones(1)*0.; ones(10*i)*4.8181; ones(40*i)*0.; ones(10*i)*-3.613; ones(40*i)*0.] #1C HPPC Experiment Current Profile
    Tk = ones(size(Iapp))*298.15 #Cell Temperature
    tDra = 0:(1.0/i):((length(Iapp)-1)/i)

    #Simulate Model
    CellV_, Ce_, jNeg, jPos, RtotNeg, RtotPos, η0, ηL, η_neg, η_pos, ϕ_ẽ1, ϕ_ẽ2, Uocp_Neg_, Uocp_Pos_, ϕ_e, Cse_Neg_, Cse_Pos_ = Sim_Model(Cell,Iapp,Tk,Init_SOC,A,B,C,D)

    #Flatten output into tuples
    CellV = flatten_(CellV,CellV_)
    Ce = flatten_(Ce,Ce_)
    Cse_Neg = flatten_(Cse_Neg,Cse_Neg_)
    Cse_Pos = flatten_(Cse_Pos,Cse_Pos_)
    Uocp_Neg = flatten_(Uocp_Neg,Uocp_Neg_)
    Uocp_Pos = flatten_(Uocp_Pos,Uocp_Pos_)

return CellV, Ce, jNeg, jPos, RtotNeg, RtotPos, η0, ηL, η_neg, η_pos, ϕ_ẽ1, ϕ_ẽ2, Uocp_Neg, Uocp_Pos, ϕ_e, Cse_Neg, Cse_Pos, tDra
end

#---------- Simulate Model -----------------#
CellV, Ce, jNeg, jPos, RtotNeg, RtotPos, η0, ηL, η_neg, η_pos, ϕ_ẽ1, ϕ_ẽ2, Uocp_Neg, Uocp_Pos, ϕ_e, Cse_Neg, Cse_Pos, tDra = Sim_loop(Cell, Init_SOC)

#----------- Plotting ---------------------------#
plotly()
display(plot(tDra[1:end-1],CellV[1][1:end-1], legend=:topright,color=:blue,bottom_margin=5Plots.mm, left_margin = 5Plots.mm, right_margin = 15Plots.mm, ylabel = "Terminal Voltage [V]", xlabel = "Time [s]"))
display(plot(tDra[1:end-1],Ce[1][1:end-1,[1,2,19,20]], legend=:topright,bottom_margin=5Plots.mm, left_margin = 5Plots.mm, right_margin = 15Plots.mm, ylabel = "Electrolyte Concen. [mol/m^3]", xlabel = "Time [s]"))
display(plot(tDra[1:end-1],Cse_Pos[1][1:end-1,:], legend=:topright,bottom_margin=5Plots.mm, left_margin = 5Plots.mm, right_margin = 15Plots.mm, ylabel = "Pos. Electrode Concen. [mol/m^3]", xlabel = "Time [s]"))
display(plot(tDra[1:end-1],Cse_Neg[1][1:end-1,:], legend=:topright,bottom_margin=5Plots.mm, left_margin = 5Plots.mm, right_margin = 15Plots.mm, ylabel = "Neg. Electrode Concen. [mol/m^3]", xlabel = "Time [s]"))

max_stars_repo_name: BradyPlanden/LIBRA

max_issues_repo_name: BradyPlanden/LIBRA

max_forks_repo_name: BradyPlanden/LIBRA


text:
#include("ER_graphs.jl")

using Memoize

@memoize function psi(PG::Pair_ER,k,d1,d2)
    exp(PG.λ*PG.s)*PG.s^k*(1-PG.s)^(d1+d2-2*k)/(PG.λ^k)
end

@memoize function psi_log(PG::Pair_ER,k,d1,d2)
    PG.λ*PG.s+log(PG.s^k*(1-PG.s)^(d1+d2-2*k)/(PG.λ^k))
end

@memoize function psi(λ,s,k,d1,d2)
    exp(λ*s)*s^k*(1-s)^(d1+d2-2*k)/(λ^k)
end

@memoize function test_psi(k::Int,d1::Int,d2::Int)
    deg_min = min(d1,d2)+1
    factorial(d1-k)*factorial(d2-k)*factorial(k)/(factorial(d1)*factorial(d2)*deg_min)
end

function overlap(PG,max_index1,max_index2)
    sum([PG.perm_G2G1[i[2]]==i[1] for i in max_index1]), sum([PG.perm_G1G2[i[1]]==i[2] for i in max_index2])
end

"""
    eval_M(PG,M)

Inputs: a pair of networks `PG` and a matrix `M` of scores where M_{iu} is
the score of pairing vertex i in G1 to vertex u in G2.
Creates the mapping from vertices in G1 to vertices in G2 by taking
the argmax on each row of `M`, then compute the overlap with the true permutation
(i.e. count the number of vertices in G1 matched correctly).
Returns also the sum of the max scores on each row.
Does a similar computation on the columns of `M`.
"""
function eval_M(PG,M)
    vmax2, max_index2 = findmax(M,dims=2)
    vmax1, max_index1 = findmax(M,dims=1)
    ov1, ov2 = overlap(PG,max_index1, max_index2)
    ov1, ov2, sum(vmax1), sum(vmax2)
end

"""
    match_edges(PG,map1,map2)

Inputs: a pair of networks `PG`, `map1` and `map2` such that map1[i] = u
for a vertex i in G1 and u in G2 and similarly map2[v] = j.
"""
function match_edges(PG,map1,map2)
    sum(has_edge(PG.N2.G,map1[src(e)],map1[dst(e)]) for e in edges(PG.N1.G)), sum(has_edge(PG.N1.G,map2[src(e)],map2[dst(e)]) for e in edges(PG.N2.G))
end

function eval_edges(PG,M)
    _, max_index2 = findmax(M,dims=2)
    _, max_index1 = findmax(M,dims=1)
    map1 = [c[2] for c in max_index2]
    map2 = [c[1] for c in max_index1]
    match_edges(PG,map1,map2)
end

function log_sum_exp(a::Float64, b::Float64)
    if a > b
        return a + log1p(exp(b-a))
    else
        return b + log1p(exp(a-b))
    end
end


text:
@testset "SHIP Basis" begin

##

@info("-------- TEST 🚢  BASIS ---------")
using SHIPs, JuLIP, BenchmarkTools, LinearAlgebra, Test, Random, StaticArrays
using SHIPs: eval_basis!, eval_basis, PolyCutoff1s, PolyCutoff2s
using JuLIP.MLIPs: IPSuperBasis
using JuLIP.Testing: print_tf
using Printf

function randR()
   R = rand(JVecF) .- 0.5
   return (0.9 + 2 * rand()) * R/norm(R)
end
randR(N) = [ randR() for n=1:N ], zeros(Int16, N)
randR(N, syms) = randR(N)[1], rand( Int16.(atomic_number.(syms)), N )
function randiso()
   K = @SMatrix rand(3,3)
   K = K - K'
   Q = rand([-1,1]) * exp(K)
end
function randiso(Rs)
   Q = randiso()
   return [ Q * R for R in shuffle(Rs) ]
end

##

trans = PolyTransform(2, 1.0)
cutf = PolyCutoff2s(2, 0.5, 3.0)

ship2 = SHIPBasis(SparseSHIP(2, 15; wL = 2.0), trans, cutf)
ship3 = SHIPBasis(SparseSHIP(3, 13; wL = 2.0), trans, cutf)
ship4 = SHIPBasis(SparseSHIP(4, 10; wL = 1.5), trans, cutf)
ship5 = SHIPBasis(SparseSHIP(5,  8; wL = 1.5), trans, cutf)
ship6 = SHIPBasis(SparseSHIP(6,  8; wL = 1.5), trans, cutf)
ships = [ship2, ship3, ship4, ship5, ship6]

@info("Test (de-)dictionisation of basis sets")
for ship in ships
   print_tf(@test (decode_dict(Dict(ship)) == ship))
   tmpf = tempname() * ".json"
   save_json(tmpf, Dict(ship))
   print_tf(@test decode_dict(load_json(tmpf)) == ship)
end

@info("Test isometry invariance for 3B-6B 🚢 s")
for ntest = 1:30
   Rs, Zs = randR(20)
   BB = [ eval_basis(🚢, Rs, Zs, 0) for 🚢 in ships ]
   RsX = randiso(Rs)
   BBX = [ eval_basis(🚢, RsX, Zs, 0) for 🚢 in ships ]
   for (B, BX) in zip(BB, BBX)
      print_tf(@test B ≈ BX)
   end
end
println()

##
@info("Test gradients for 3-6B 🚢-basis")
for 🚢 in ships
   @info("  body-order = $(SHIPs.bodyorder(🚢)):")
   Rs, Zs = randR(20)
   tmp = SHIPs.alloc_temp_d(🚢, Rs)
   # SHIPs.precompute_grads!(tmp, 🚢, Rs, Zs)
   B = eval_basis(🚢, Rs, Zs, 0)
   dB = SHIPs.alloc_dB(🚢, Rs)
   SHIPs.eval_basis_d!(dB, tmp, 🚢, Rs, Zs, 0)
   @info("      finite-difference test into random directions")
   for ndirections = 1:20
      Us, Zs = randR(length(Rs))
      errs = Float64[]
      for p = 2:10
         h = 0.1^p
         Bh = eval_basis(🚢, Rs+h*Us, Zs, 0)
         dBh = (Bh - B) / h
         dBxU = sum( dot.(Ref(Us[n]), dB[n,:])  for n = 1:length(Rs) )
         push!(errs, norm(dBh - dBxU, Inf))
      end
      success = (/(extrema(errs)...) < 1e-3) || (minimum(errs) < 1e-10)
      print_tf(@test success)
   end
   println()
end

##
verbose=false
@info("Test gradients for 3B with and R near the pole")
🚢 = ship2
@info("  body-order = $(SHIPs.bodyorder(🚢)):")
# Rs = [ randR(5); [SVector(1e-14*rand(), 1e-14*rand(), 1.1+1e-6*rand())] ]
Rs = [ randR(5)[1]; [SVector(0, 0, 1.1+0.5*rand())]; [SVector(1e-14*rand(), 1e-14*rand(), 0.9+0.5*rand())] ]
_, Zs = randR(length(Rs))
tmp = SHIPs.alloc_temp_d(🚢, Rs)
B = eval_basis(🚢, Rs, Zs, 0)
dB = SHIPs.alloc_dB(🚢, Rs)
SHIPs.eval_basis_d!(dB, tmp, 🚢, Rs, Zs, 0)
@info("      finite-difference test into random directions")
for ndirections = 1:30
   Us, _ = randR(length(Rs))
   errs = Float64[]
   for p = 2:10
      h = 0.1^p
      Bh = eval_basis(🚢, Rs+h*Us, Zs, 0)
      dBh = (Bh - B) / h
      dBxU = sum( dot.(Ref(Us[n]), dB[n,:])  for n = 1:length(Rs) )
      push!(errs, norm(dBh - dBxU, Inf))
      verbose && (@printf("  %2d | %.2e \n", p, errs[end]))
   end
   success = (/(extrema(errs)...) < 1e-3) || (minimum(errs) < 1e-10)
   print_tf(@test success)
end
println()

##
@info("Check Correctness of SHIPBasis calculators")

randcoeffs(B) = 2 * (rand(length(B)) .- 0.5) .* (1:length(B)).^(-2)

naive_energy(basis::SHIPBasis, at) =
      sum( eval_basis(basis, R, zeros(Int16, length(R)), 0)
            for (i, j, R) in sites(at, cutoff(basis)) )

for basis in ships
   @info("   body-order = $(SHIPs.bodyorder(basis))")
   at = bulk(:Si) * 3
   at.Z[:] .= 0   # this set of tests is species-agnostic!
   rattle!(at, 0.1)
   print("     energy: ")
   println(@test energy(basis, at) ≈ naive_energy(basis, at) )
   print("site-energy: ")
   println(@test energy(basis, at) ≈ sum( site_energy(basis, at, n)
                                         for n = 1:length(at) ) )
   # we can test consistency of forces, site energy etc by taking
   # random inner products with coefficients
   @info("     a few random combinations")
   for n = 1:10
      c = randcoeffs(basis)
      sh = JuLIP.MLIPs.combine(basis, c)
      print_tf(@test energy(sh, at) ≈ dot(c, energy(basis, at)))
      print_tf(@test forces(sh, at) ≈ sum(c*f for (c, f) in zip(c, forces(basis, at))) )
      print_tf(@test site_energy(sh, at, 5) ≈ dot(c, site_energy(basis, at, 5)))
      print_tf(@test site_energy_d(sh, at, 5) ≈ sum(c*f for (c, f) in zip(c, site_energy_d(basis, at, 5))) )
   end
   println()
end

##
end

max_stars_repo_name: cortner/SHIPs

max_issues_repo_name: cortner/SHIPs

max_forks_repo_name: cortner/SHIPs


text:
using Test
using Revise
using PolynomialBasis

function interpolating_gradient(poly)
    order = PolynomialBasis.order(poly)
    interpgrad = InterpolatingPolynomial(2, 2, order)

    points =
        PolynomialBasis.interpolation_points(PolynomialBasis.basis(interpgrad))
    interpvals = mapslices(x -> vec(gradient(poly, x)), points, dims = 1)

    PolynomialBasis.update!(interpgrad, interpvals)
    return interpgrad
end

function allapprox(u, v)
    return all(u .≈ v)
end

function linear(coeffs, v)
    x, y = v
    a, b, c, d = coeffs

    return a * x * y + b * x + c * y + d
end

function gradlinear(coeffs, v)
    x, y = v
    a, b, c, d = coeffs

    fx = a * y + b
    fy = a * x + c
    return [fx, fy]
end

coeffs = [1.0, 6.0, 2.0, 7.0]
poly = InterpolatingPolynomial(1, 2, 1)
points = PolynomialBasis.interpolation_points(PolynomialBasis.basis(poly))
interpcoeffs = vec(mapslices(x -> linear(coeffs, x), points, dims = 1))
update!(poly, interpcoeffs)

testp = rand(2, 10)
vals = mapslices(poly, testp, dims = 1)
testvals = mapslices(x -> linear(coeffs, x), testp, dims = 1)
@test allapprox(vals, testvals)

interpgrad = interpolating_gradient(poly)
gvals = mapslices(interpgrad, testp, dims = 1)
testgvals = mapslices(x -> gradlinear(coeffs, x), testp, dims = 1)
@test allapprox(gvals, testgvals)

function quadratic(coeffs, v)
    x, y = v
    a, b, c, d, e, f, g, h, i = coeffs

    return a * x^2 * y^2 +
           b * x^2 * y +
           c * x^2 +
           d * x * y^2 +
           e * x * y +
           f * x +
           g * y^2 +
           h * y +
           i
end

function gradquadratic(coeffs, v)
    x, y = v
    a, b, c, d, e, f, g, h, i = coeffs

    fx = 2 * a * x * y^2 + 2 * b * x * y + 2 * c * x + d * y^2 + e * y + f
    fy = 2 * a * x^2 * y + b * x^2 + 2 * d * x * y + e * x + 2 * g * y + h

    return [fx, fy]
end

coeffs = [3.0, 4.0, 1.0, 11.0, 7.0, 3.0, 5.0, 13.0, 15.0]
poly = InterpolatingPolynomial(1, 2, 2)
points = PolynomialBasis.interpolation_points(PolynomialBasis.basis(poly))
interpcoeffs = vec(mapslices(x -> quadratic(coeffs, x), points, dims = 1))
update!(poly, interpcoeffs)

testp = rand(2, 10)
vals = mapslices(poly, testp, dims = 1)
testvals = mapslices(x -> quadratic(coeffs, x), testp, dims = 1)
@test allapprox(vals, testvals)

interpgrad = interpolating_gradient(poly)
gvals = mapslices(interpgrad, testp, dims = 1)
testgvals = mapslices(x -> gradquadratic(coeffs, x), testp, dims = 1)
@test allapprox(gvals, testgvals)

function cubic(coeffs, v)
    x, y = v
    a, b, c, d, e, f, g, h, i, j = coeffs

    return a * x^3 +
           b * x^2 * y +
           c * x * y^2 +
           d * y^3 +
           e * x^2 +
           f * x * y +
           g * y^2 +
           h * x +
           i * y +
           j
end

function gradcubic(coeffs, v)
    x, y = v
    a, b, c, d, e, f, g, h, i, j = coeffs

    fx = 3 * a * x^2 + 2 * b * x * y + c * y^2 + 2 * e * x + f * y + h
    fy = b * x^2 + 2 * c * x * y + 3 * d * y^2 + f * x + 2 * g * y + i

    return [fx,fy]
end

# coeffs = [7.,2,4,9,6,3,5,7,1,11]
coeffs = [10,12,6,17,9,5,2,1,8,11.0]
poly = InterpolatingPolynomial(1, 2, 3)
points = PolynomialBasis.interpolation_points(PolynomialBasis.basis(poly))
interpcoeffs = vec(mapslices(x -> cubic(coeffs, x), points, dims = 1))
update!(poly, interpcoeffs)

testp = rand(2, 10)
vals = mapslices(poly, testp, dims = 1)
testvals = mapslices(x -> cubic(coeffs, x), testp, dims = 1)
@test allapprox(vals, testvals)

interpgrad = interpolating_gradient(poly)
gvals = mapslices(interpgrad, testp, dims = 1)
testgvals = mapslices(x -> gradcubic(coeffs, x), testp, dims = 1)
@test allapprox(gvals, testgvals)


text:
using Test
using JSExpr
using JSExpr: @crawl, deparse, JSTerminal, JSString

@testset "interpolation" begin
    @testset "interpolation of simple datatypes" begin
        foostr = "foo"
        @test @js(foo = $foostr) == JSString("foo = \"foo\"")

        answer = 42
        @test @js(theAnswer = $answer) == JSString("theAnswer = 42")
    end

    @testset "interpolation of compound datatypes" begin
        options = Dict("foo" => "bar")
        @test @js(options = $options) == JSString("""options = {"foo":"bar"}""")

        numbers = [1, 2, 3, 4]
        @test @js(numbers = $numbers) == JSString("numbers = [1,2,3,4]")
    end

    @testset "interpolation of javascript" begin
        callback = JSString("""() => console.log("done!")""")
        @test @js(callback = $callback) == JSString("callback = $(callback.s)")
    end
end


text:
S = [:healthy, :in_situ, :invasive, :death]
A = [:wait, :test, :treat]
O = [:pos, :neg]
γ = 0.99
s0 = Deterministic(:healthy)
term = Set([:death])

function T(s, a, sp)
    if s == :healthy
        if sp == :in_situ
            return 0.02
        elseif sp == s
            return 0.98
        else
            return 0.0
        end
    elseif s == :in_situ
        if a == :treat
            if sp == :healthy
                return 0.60
            elseif sp == s
                return 0.40
            else
                return 0.0
            end
        else #a == :test || a == :wait
            if sp == :invasive
                return 0.10
            elseif sp == s
                return 0.90
            else
                return 0.0
            end
        end
    elseif s == :invasive
        if a == :treat
            if sp == :healthy
                return 0.20
            elseif sp == :death
                return 0.20
            elseif sp == s
                return 0.60
            else
                return 0.0
            end
        else
            if sp == :death
                return 0.60
            elseif sp == s
                return 0.40
            else
                return 0.0
            end
        end
    else # s == :death
        return 0.25
    end
end

function Z(a, sp, o)
    if a == :test
        if sp == :healthy
            if o == :pos
                return 0.05
            else
                return 0.95
            end
        elseif sp == :in_situ
            if o == :pos
                return 0.80
            else
                return 0.20
            end
        else
            if o == :pos
                return 1.0
            else
                return 0.0
            end
        end
    elseif a == :treat
        if sp == :in_situ || sp == :invasive
            if o == :pos
                return 1.0
            else
                return 0.0
            end
        else
            if o == :pos
                return 0.0
            else
                return 1.0
            end
        end
    else #a == :wait
        if o == :pos
            return 0.0
        else
            return 1.0
        end
    end
end

function R(s, a)
    if s == :death
        return 0.0
    elseif a == :wait
        return 1.0
    elseif a == :test
        return 0.80
    else
        return 0.10
    end
end

c = DiscreteExplicitPOMDP(S, A, O, T, Z, R, γ, s0, terminals=term);

function evalSolver(m)
    results = Dict{String, Float64}()
    for (key, solver) in ["SARSOP"=>SARSOPSolver()]
        policy = solve(solver, m)
        N = 10000
        rsum = 0.0
        for i in 1:N
            rsum += simulate(RolloutSimulator(max_steps=500), m, policy)
        end
        results[key] = rsum/N
    end;
    return results
end

results_c = evalSolver(c);
@show results_c["SARSOP"]
@test results_c["SARSOP"] >= 62.5


text:
module RocketAsyncOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: async()" begin

    println("Testing: operator async()")

    run_testset([
        (
            source = from(1:5),
            values = @ts([1, 2, 3, 4, 5, c])
        ),
        (
            source = from(1:5) |> async(),
            values = @ts(0 ~ [ 1, 2, 3, 4, 5, c ])
        ),
        (
            source = from(1:5) |> async(0),
            values = @ts([1] ~ [2] ~ [3] ~ [4] ~ [5] ~ c)
        ),
        (
            source = from(1:5) |> async(1),
            values = @ts([1] ~ [2] ~ [3] ~ [4] ~ [5] ~ c)
        ),
        (
            source = from(1:5) |> async(2),
            values = @ts([ 1, 2 ] ~ [ 3, 4 ] ~ [ 5, c ])
        ),
        (
            source = from(1:5) |> async(3),
            values = @ts([ 1, 2, 3 ] ~ [ 4, 5, c ])
        ),
        (
            source = completed() |> async(0),
            values = @ts(c)
        ),
        (
            source = faulted("e") |> async(0),
            values = @ts(e("e")),
        ),
        (
            source = never() |> async(0),
            values = @ts()
        )
    ])

end

end


text:
const n = 10

@testset "vectors" begin
    @testset "BipolarHDV" begin

        hdv = BipolarHDV(n)

        @test length(hdv) == n
        @test eltype(hdv) <: Int
        @test hdv[2] isa Int
        @test all(-1 .≤ hdv .≤ 1)
        @test hdv == BipolarHDV(hdv.v)
        @test similar(hdv) isa BipolarHDV
        @test sum(hdv) ≈ sum(hdv.v)

        hdv_offset = BipolarHDV(hdv.v, 2)
        @test hdv_offset.offset == 2
        @test [hdv[i-2] for i in 1:5] ≈ [hdv_offset[i] for i in 1:5]
        hdv_offset[5] = 0
        @test hdv[3] == 0
    end

    @testset "BinaryHDV" begin

        hdv = BinaryHDV(n)

        @test length(hdv) == n
        @test eltype(hdv) <: Bool
        @test hdv[2] isa Bool
        @test hdv == BinaryHDV(hdv.v)
        @test similar(hdv) isa BinaryHDV
        @test sum(hdv) ≈ sum(hdv.v)

        hdv_offset = BinaryHDV(hdv.v, 2)
        @test hdv_offset.offset == 2
        @test [hdv[i-2] for i in 1:5] ≈ [hdv_offset[i] for i in 1:5]
        hdv_offset[5] = 0
        @test hdv[3] == 0
    end

    @testset "GradedBipolarHDV" begin

        hdv = GradedBipolarHDV(n)

        @test length(hdv) == n
        @test eltype(hdv) <: Real
        @test hdv[2] isa Real
        @test all(-1 .≤ hdv .≤ 1)
        @test hdv == GradedBipolarHDV(hdv.v)
        @test similar(hdv) isa GradedBipolarHDV
        @test sum(hdv) ≈ sum(hdv.v)
        @test eltype(GradedBipolarHDV(Float32, n)) <: Float32

        @test all(-0.4 .≤ GradedBipolarHDV(l=-0.4, u=0.6) .≤ 0.6)

        hdv_offset = GradedBipolarHDV(hdv.v, 2)
        @test hdv_offset.offset == 2
        @test [hdv[i-2] for i in 1:5] ≈ [hdv_offset[i] for i in 1:5]
        hdv_offset[5] = 0
        @test hdv[3] == 0
    end

    @testset "GradedHDV" begin

        hdv = GradedHDV(n)

        @test length(hdv) == n
        @test eltype(hdv) <: Real
        @test hdv[2] isa Real
        @test all(0 .≤ hdv .≤ 1)
        @test hdv == GradedHDV(hdv.v)
        @test similar(hdv) isa GradedHDV
        @test sum(hdv) ≈ sum(hdv.v)
        @test eltype(GradedHDV(Float32, n)) <: Float32

        @test all(0.4 .≤ GradedHDV(l=0.4, u=0.6) .≤ 0.6)

        hdv_offset = GradedHDV(hdv.v, 2)
        @test hdv_offset.offset == 2
        @test [hdv[i-2] for i in 1:5] ≈ [hdv_offset[i] for i in 1:5]
        hdv_offset[5] = 0
        @test hdv[3] == 0
    end

    @testset "RealHDV" begin

        hdv = RealHDV(n)

        @test length(hdv) == n
        @test eltype(hdv) <: Real
        @test hdv[2] isa Real

        @test hdv == RealHDV(hdv.v)
        @test similar(hdv) isa RealHDV
        @test sum(hdv) ≈ sum(hdv.v)
        @test eltype(RealHDV(Float32, n)) <: Float32

        hdv_offset = RealHDV(hdv.v, 2)
        @test hdv_offset.offset == 2
        @test [hdv[i-2] for i in 1:5] ≈ [hdv_offset[i] for i in 1:5]
        hdv_offset[5] = 0
        @test hdv[3] == 0
    end
end


text:
Base.@kwdef struct Path
    directory::String = ""
    filename::String = ""
end


text:
# module Bukdu.Plug

immutable Static
end

import HttpServer: Response, mimetypes
import ..Octo: Assoc

function FileResponse(filename)
    if isfile(filename)
        s = open(Base.read, filename)
        (_, ext) = splitext(filename)
        mime = length(ext)>1 && haskey(mimetypes,ext[2:end]) ? mimetypes[ext[2:end]] : "application/octet-stream"
        Response(200, Dict{AbstractString,AbstractString}([("Content-Type",mime)]), s)
    else
        Response(404, "Not Found - file $filename could not be found")
    end
end

import ..ApplicationController, ..Routing, ..Conn

immutable StaticController <: ApplicationController
end

function read(c::StaticController)
    filepath = c[:assigns][:filepath]
    resp = FileResponse(filepath)
    params = Assoc()
    query_params = Assoc()
    private = Dict{Symbol,String}()
    assigns = Dict{Symbol,String}()
    Conn(resp.status, Dict{String,String}(resp.headers), resp.data, params, query_params, private, assigns)
end

"""
plug `Plug.Static` to serve the static files.

```julia
Endpoint() do
    plug(Plug.Static, at= "/", from= "public")
end
```
"""
function plug(::Type{Plug.Static}; kw...)
    # at::String, from::String, only::Vector{String})
    opts = Dict(kw)
    at = opts[:at]
    from = opts[:from]
    try_index_html = true
    if haskey(opts, :only)
        only = opts[:only]
        has_only = !isempty(only)
    else
        has_only = false
    end
    if haskey(opts, :try_index_html)
        try_index_html = opts[:try_index_html]
    end
    for (root, dirs, files) in walkdir(from)
        for file in files
            filepath = joinpath(root, file)
            opts = Dict(:assigns => Dict{Symbol,String}(:filepath=>filepath))
            if try_index_html && "index.html" == file
                Routing.match(get, "/", StaticController, read, opts)
            end
            if has_only
                !any(x->startswith(filepath, x), only) && continue
            end
            reqpath = filepath[length(from)+1:end]
            Routing.match(get, reqpath, StaticController, read, opts)
        end
    end
end