text:
module LavaCore

include("common/Common.jl")

using VulkanCore
include("common/VkExt.jl")
using .VkExt
export VkExt
include("common/Utils.jl")
using .Utils
export Utils
include("camera/Camera.jl")
using .Camera
export Camera

using GeometryTypes
include("common/FormatInfo.jl")
include("common/ShaderExtension.jl")
include("common/FeatureBase.jl")

include("objects/ShaderModule.jl")

include("createinfos/DescriptorSetLayoutCreateInfo.jl")
include("createinfos/DescriptorPoolCreateInfo.jl")
include("createinfos/PipelineLayoutCreateInfo.jl")
include("createinfos/RenderPassCreateInfo.jl")
include("createinfos/PipelineColorBlendStateCreateInfo.jl")
include("createinfos/PipelineDepthStencilStateCreateInfo.jl")
include("createinfos/PipelineDynamicStateCreateInfo.jl")
include("createinfos/PipelineInputAssemblyStateCreateInfo.jl")
include("createinfos/PipelineMultisampleStateCreateInfo.jl")
include("createinfos/PipelineRasterizationStateCreateInfo.jl")
include("createinfos/PipelineShaderStageCreateInfo.jl")
include("createinfos/PipelineTessellationStateCreateInfo.jl")
include("createinfos/PipelineVertexInputStateCreateInfo.jl")
include("createinfos/PipelineViewportStateCreateInfo.jl")
include("createinfos/GraphicsPipelineCreateInfo.jl")
include("createinfos/ComputePipelineCreateInfo.jl")
include("createinfos/ImageCreateInfo.jl")
include("createinfos/BufferCreateInfo.jl")
include("createinfos/SamplerCreateInfo.jl")

include("gpuselection/SelectionStrategy.jl")

include("objects/Queue.jl")
include("objects/DescriptorPool.jl")
include("objects/DescriptorSetLayout.jl")
include("objects/PipelineLayout.jl")
include("objects/RenderPass.jl")
include("objects/GraphicsPipeline.jl")
include("objects/ComputePipeline.jl")

include("objects/MemoryChunk.jl")
include("objects/Suballocator.jl")

include("objects/Device.jl")
include("objects/Image.jl")
include("objects/ImageView.jl")
include("objects/Framebuffer.jl")
include("objects/CommandBuffer.jl")
include("objects/Instance.jl")
include("objects/Buffer.jl")
include("objects/Sampler.jl")
include("objects/DescriptorSetWriter.jl")
include("objects/DescriptorSet.jl")

include("raii/ActiveRenderPass.jl")
include("raii/Barriers.jl")

end


text:
plane = ComplexPlane(Quaternion(rand(), ℝ³(rand(3)))) # Construct a unit Quaternion
w, z = plane.z₁, plane.z₂
B = Complex(10rand(), 10rand())
clifford = Clifford(w, z, B)

@test typeof(clifford.total) <: S³
@test typeof(clifford.base) <: S²
@test isapprox(clifford.total, plane)

ext: jl


text:
using OED
using Base.Test
using MAT
using jInv.Mesh

r = 3;
l = 3;
nrhs = 10
Mimg  = getRegularMesh([-1 1 -1 1.],[4,4])
Mgrid = getRegularMesh([-1.2 2.4 -1.6 1.4],[8,8])

Ftrue = randn(Mimg.nc,nrhs).^2
for k=1:nrhs
	Ftrue[k,k] = 10;
end
OEDparam = getOEDTypeBParam(r,l,Ftrue,Mimg,Mgrid)
display(OEDparam)
L = speye(Mimg.nc)

IPmodes = (:unconstr,:nonnegconstr,:nonnegconstr,:boxconstr)

Ce = ones(1,Mimg.nc); ce = Ce*Ftrue;
Ci = [speye(Mimg.nc);-speye(Mimg.nc)];
ci = [minimum(Ftrue)*ones(Mimg.nc,nrhs); -maximum(Ftrue)*ones(Mimg.nc,nrhs)]
m0 = rand(l)
mOpt = []
for ip=IPmodes
mO = runOEDTypeB([OEDparam],Ftrue,L,ip,copy(m0);Ce=Ce,ce=ce,Ci=Ci,ci=ci,pcgMaxIter=10,maxIter=3,maxStep=10.0,exName="/tmp/OEDb")
push!(mOpt,mO)
end


text:
module Test3D

using Test
using Trixi2Vtk

include("test_trixi2vtk.jl")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

# Create artifacts directory where all files that should be preserved will be stored
artifacts_dir = joinpath(pathof(Trixi2Vtk) |> dirname |> dirname, "artifacts")
if !isdir(artifacts_dir)
  mkdir(artifacts_dir)
end

@testset "3D" begin
  @testset "TreeMesh" begin
    isdir(outdir) && rm(outdir, recursive=true)
    run_trixi(joinpath("3d", "elixir_advection_extended.jl"), maxiters=1)

    @testset "uniform mesh" begin
      test_trixi2vtk("solution_000000.h5", outdir,
          hashes=[("solution_000000.vtu", "6ab3aa525851187ee0839e1d670a254a66be4ad7"),
                  ("solution_000000_celldata.vtu", "99c782d732a4d1f6764013c3fba2cdaddf3927ab")])

      # Store output files as artifacts to facilitate debugging of failing tests
      outfiles = ("solution_000000.vtu", "solution_000000_celldata.vtu")
      testname = "3d-tree-mesh-uniform-mesh"
      for outfile in outfiles
        println("Copying '", abspath(joinpath(outdir, outfile)),
                "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)),
                "'...")
        cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile), force=true)
      end
    end
  end

  @testset "CurvedMesh" begin
    isdir(outdir) && rm(outdir, recursive=true)
    run_trixi(joinpath("3d", "elixir_advection_basic_curved.jl"), maxiters=1)

    @testset "basic" begin
      test_trixi2vtk("solution_000000.h5", outdir,
          hashes=[("solution_000000.vtu", "57c58480d8f99b7e9d365cb3fa71790db42de750")])

      # Store output files as artifacts to facilitate debugging of failing tests
      outfiles = ("solution_000000.vtu",)
      testname = "3d-curved-mesh-basic"
      for outfile in outfiles
        println("Copying '", abspath(joinpath(outdir, outfile)),
                "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)),
                "'...")
        cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile), force=true)
      end
    end
  end
end

# Clean up afterwards: delete Trixi output directory
@test_skip rm(outdir, recursive=true)

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "libsharp2"
 # ↓ This version number is a lie to be able to build for experimental
 # ↓ platforms, but it shouldn't be a problem as the library is abandoned.
version = v"1.0.2"

# Collection of sources required to complete build
sources = [
    GitSource("https://gitlab.mpcdf.mpg.de/mtr/libsharp.git", "54856313a5fcfb6a33817b7dfa28c4b1965ffbd1"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd libsharp/
sed -i 's/LT_INIT.*/LT_INIT\(\[win32-dll\]\)/g' configure.ac
sed -i 's/libsharp2_la_LDFLAGS = -version-info 0:0:0/libsharp2_la_LDFLAGS = -version-info 0:0:0 -no-undefined/g' Makefile.am
autoreconf -i
if [[ "${target}" == x86_64-* ]] && [[ "${target}" != *-apple-* ]] && [[ "${target}" != *-w64-* ]]; then
    # Enable runtime multiarch, but not for macOS nor Windows. For Mac OS X, see
    # https://gitlab.mpcdf.mpg.de/mtr/libsharp/-/blob/6374a3a1ffb935443c56f09b371b11cf982a7e28/COMPILE
    # For Windows, see https://github.com/ziotom78/Libsharp.jl/issues/3
    export CFLAGS="-DMULTIARCH"
fi
export CFLAGS="${CFLAGS} -O3"
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               # GCC 7 needed to enable multiarch.
               preferred_gcc_version=v"7", julia_compat="1.6")


text:
"""
    _foldl_blockarray(rf, acc, coll::BlockArrays.BlockArray)
"""
_foldl_blockarray(rf, acc, coll) =
    complete(rf, @return_if_reduced _foldl_blockarray(
        rf, acc, coll,
        (),  # (outer) block indices
        (),  # (outer) offset indices
        Val(ndims(coll) - 1),  # stop recursion when `== Val(length(block))`
    ))

# Looping over outer/right dimensions; i.e., one of `b`, `c`, `d` in
# `coll[a, b, c, d]` but not `a`.
@inline function _foldl_blockarray(rf, acc, coll, block, offset, ValN₋₁)
    nblocks = BlockArrays.nblocks
    blocksize = BlockArrays.blocksize

    i = ndims(coll) - length(block)
    for j in 1:nblocks(coll, i)
        for k in 1:blocksize(coll, i, j)
            acc = @return_if_reduced _foldl_blockarray(
                rf, acc, coll, (j, block...), (k, offset...), ValN₋₁,
            )
        end
    end
    return acc
end

# Finally, looping over the inner/left -most dimension:
@inline function _foldl_blockarray(
    rf,
    acc,
    coll,
    block::NTuple{N₋₁, Int},
    offset::NTuple{N₋₁, Int},
    ::Val{N₋₁},
) where {N₋₁}

    nblocks = BlockArrays.nblocks
    blocksize = BlockArrays.blocksize
    Block = BlockArrays.Block

    @inbounds for j in 1:nblocks(coll, 1)
        array = coll[Block(j, block...)]
        @simd_if rf for k in 1:blocksize(coll, 1, j)
            acc = @next(rf, acc, array[k, offset...])
        end
    end
    return acc
end

# TODO: write reduce for BlockArrays which can be done in the "natural" order


text:
# Julia script to run set of models in tuple modeltypes on region
# User sets total number of mcmc steps, folders for input and output and temperature
# run using:  julia -p #nprocs SICRPoissonscrptBW.jl

t0 = time()

using CSV
using DataFrames
using Dates
using DelimitedFiles
using DifferentialEquations
using Distributed
using Distributions
using LinearAlgebra
using StatsBase
using Statistics
using NLsolve
using StochasticDiffEq
@everywhere using CSV
@everywhere using DataFrames
@everywhere using Dates
@everywhere using DelimitedFiles
@everywhere using DifferentialEquations
@everywhere using Distributed
@everywhere using Distributions
@everywhere using LinearAlgebra
@everywhere using StatsBase
@everywhere using Statistics
@everywhere using NLsolve
@everywhere using StochasticDiffEq

@everywhere include("SICRPoisson.jl")

const txtstr = ".txt"

function getdatafiles(datafolder)
    datafiles = readdir("data")
    datafiles = datafiles[occursin.("covid_timeseries",datafiles)]
    regions = Array{String,1}(undef,0)
    for datafile in datafiles
        push!(regions,replace(replace(datafile,".csv"=>""),"covid_timeseries_"=>""))
    end
    return datafiles,regions
end

# beta,sigmac,sigmar,sigmad,q,f,mbase,mlocation,mrate,cmax,c50 = p
function getpfit(modelname)
    if modelname == "r"
        return  [1;2;3;4;12]
    end
    if modelname == "rq"
        return [1;2;3;4;5;12]
    end
    if modelname == "rf"
        return [1;2;3;4;6;12]
    end
    if modelname == "rqf"
        return [1;2;3;4;5;6;12]
    end
    if modelname == "rm"
        return [1;2;3;4;7;8;9;12]
    end
    if modelname == "rqm"
        return [1;2;3;4;5;7;8;9;12]
    end
    if modelname == "rfm"
        return [1;2;3;4;6;7;8;9;12]
    end
    if modelname == "rqfm"
        return [1;2;3;4;5;6;7;8;9;12]
    end
    if modelname == "rmt"
        return [1;2;3;4;5;6;7;8;9;10;11;12]
    end
    if modelname == "rqmt" || modelname == "rfm"
        return [1;2;3;4;5;6;7;8;9;10;11;12]
    end
    if modelname == "rfmt"
        return [1;2;3;4;5;6;7;8;9;10;11;12]
    end
    if modelname == "rqfmt"
        return [1;2;3;4;5;6;7;8;9;10;11;12]
    end
    if modelname == "u"
        return [.25,.1,.1,.1]
    end
    if modelname == "uq" || modelname == "uf"
        return [.25,.1,.1,.1,.1]
    end
    if modelname == "um"
        return [.25,.1,.1,.1,5.,.1]
    end
    if modelname == "uqm" || modelname == "ufm"
        return [.25,.1,.1,.1,.1,5.,.1]
    end
    if modelname == "uqfm"
        return [.25,.1,.1,.1,.1,.1,5.,.1]
    end
end

# name = Dict(sicr! => "r", sicrq! => "rq", sicrf! => "rf", sicrqf! => "rqf", sicrm! => "rm", sicrqm! => "rqm", sicrfm! => "rfm", sicrqfm! => "rqfm",sicu! => "u", sicuq! => "uq", sicuf! => "uf", sicuq! => "uq", sicuqf! => "uqf", sicum! => "um", sicuqm! => "uqm", sicufm! => "ufm", sicuqfm! => "uqfm")

datestrin = "2020-04-08linear"
datestrout = "2020-04-08linear"
println(datestrin,"->",datestrout)
if ~isdir("Results")
    mkdir("Results")
end
pathin = "Results/$datestrin/"
pathout = "Results/$datestrout/"
if ~isdir(pathout)
    mkdir(pathout)
end

datafiles,regions = getdatafiles("data")

colsr = [5,6,7]
colsu = [5,6]

N = 1e7
total = Int(2.e3)
temp = 1.
const pinstr = "pl"

println(total,' ',temp)
println(pinstr)

# modeltype = (sicr!,sicrq!,sicrf!,sicrm!,sicu!,sicuq!,sicuf!,sicum!)
# modeltype = (sicr!,sicrq!,sicrf!,sicrfq!,sicrfqm!,sicurfqmt!)
modeltype = ["rqfmt"]

pstatstr = pathout * "pstats" * txtstr
fpstat = open(pstatstr,"w")

out = Array{Future,2}(undef,length(modeltype),length(datafiles))

for j in eachindex(datafiles)
    datafile = datafiles[j]
    df = DataFrame(CSV.File("data/" * datafile))
    region = regions[j]
    println(region)

    for i in eachindex(modeltype)

        name = modeltype[i]

        if occursin("u",modeltype[i])
            data = sicuprep(df)
            cols = colsu
        else
            data = df
            cols = colsr
        end

        infile = pathin * pinstr * "_" * region * "_" * name * txtstr

        if isfile(infile) && ~isempty(read(infile))
            p = readdlm(infile)[2:end]
            pfit = getpfit(name)
        else
            pfit = getpfit(name)
            p = [.25,.1,.01,.01,.01,1.,.1,1.,1.,.1,10.,1.]
            println(name,": No prior")
        end

        out[i,j] = @spawn mcmc(data,cols,p,pfit,N,total,linearrd!,temp)
        # println(name,": maxll: ",outtuple[4],", Accept: ",outtuple[3])
    end
end

for j in eachindex(datafiles), i in eachindex(modeltype)

    outtuple = fetch(out[i,j])

    name = modeltype[i]
    region = regions[j]

    mlstr = pathout * "maxll" * "_" * region * "_" * name * txtstr
    llstr = pathout * "llsamples" * "_" * region * "_" * name * txtstr
    pstr = pathout * "pl" * "_" * region * "_" * name * txtstr

    pstats,R0 = measures(outtuple[2],outtuple[5])
    fpstat = open(pstatstr,"a")
    writedlm(fpstat,[region * "_" * name])
    writedlm(fpstat,[outtuple[4] R0[1] R0[2] R0[3] R0[4]'])
    for ptup in pstats
        writedlm(fpstat, [ptup[1] ptup[2] ptup[3] ptup[4]'])
    end
    close(fpstat)

    fml=open(mlstr,"w")
    fml=open(mlstr,"a")
    writedlm(fml,[outtuple[4];outtuple[5]])
    close(fml)

    fchil=open(llstr,"w")
    fchil=open(llstr,"a")
    writedlm(fchil,outtuple[1])
    close(fchil)

    fpl=open(pstr,"w")
    fpll=open(pstr,"a")
    writedlm(fpl,[outtuple[4];outtuple[2][end,:]])
    close(fpl)

end

println(" Time: ",time()-t0)

max_stars_repo_name: cannin/covid-sicr-test

max_issues_repo_name: cannin/covid-sicr-test

max_forks_repo_name: cannin/covid-sicr-test


text:
"""
    $(SIGNATURES)
Save a DFG to a folder. Will create/overwrite folder if it exists.

DevNotes:
- TODO remove `compress` kwarg.

# Example
```julia
using DistributedFactorGraphs, IncrementalInference
# Create a DFG - can make one directly, e.g. LightDFG{NoSolverParams}() or use IIF:
dfg = initfg()
# ... Add stuff to graph using either IIF or DFG:
v1 = addVariable!(dfg, :a, ContinuousScalar, tags = [:POSE], solvable=0)
# Now save it:
saveDFG(dfg, "/tmp/saveDFG.tar.gz")
```
"""
function saveDFG(folder::AbstractString, dfg::AbstractDFG)

    # TODO: Deprecate the folder functionality

    # Clean up save path if a file is specified
    savepath = folder[end] == '/' ? folder[1:end-1] : folder
    savepath = splitext(splitext(savepath)[1])[1] # In case of .tar.gz

    variables = getVariables(dfg)
    factors = getFactors(dfg)
    varFolder = "$savepath/variables"
    factorFolder = "$savepath/factors"
    # Folder preparations
    if !isdir(savepath)
        @debug "Folder '$savepath' doesn't exist, creating..."
        mkpath(savepath)
    end
    !isdir(varFolder) && mkpath(varFolder)
    !isdir(factorFolder) && mkpath(factorFolder)
    # Clearing out the folders
    map(f -> rm("$varFolder/$f"), readdir(varFolder))
    map(f -> rm("$factorFolder/$f"), readdir(factorFolder))
    # Variables
    for v in variables
        vPacked = packVariable(dfg, v)
        io = open("$varFolder/$(v.label).json", "w")
        JSON2.write(io, vPacked)
        close(io)
    end
    # Factors
    for f in factors
        fPacked = packFactor(dfg, f)
        io = open("$factorFolder/$(f.label).json", "w")
        JSON2.write(io, fPacked)
        close(io)
    end

    savedir = dirname(savepath) # is this a path of just local name? #344 -- workaround with unique names
    savename = basename(string(savepath))
    @assert savename != ""
    destfile = joinpath(savedir, savename*".tar.gz")
    if length(savedir) != 0
      run( pipeline(`tar -zcf - -C $savedir $savename`, stdout="$destfile"))
    else
      run( pipeline(`tar -zcf - $savename`, stdout="$destfile"))
    end
    Base.rm(joinpath(savedir,savename), recursive=true)
end
# support both argument orders, #581
saveDFG(dfg::AbstractDFG, folder::AbstractString) = saveDFG(folder, dfg)

"""
    $(SIGNATURES)
Load a DFG from a saved folder. Always provide the IIF module as the second
parameter.

# Example
```julia
using DistributedFactorGraphs, IncrementalInference
# Create a DFG - can make one directly, e.g. LightDFG{NoSolverParams}() or use IIF:
dfg = initfg()
# Load the graph
loadDFG!(dfg, "/tmp/savedgraph.tar.gz")
# Use the DFG as you do normally.
ls(dfg)
```
"""
function loadDFG!(dfgLoadInto::AbstractDFG, dst::AbstractString)

    #
    # loaddir gets deleted so needs to be unique
    loaddir=split(joinpath("/","tmp","caesar","random", string(uuid1())), '-')[1]
    # Check if zipped destination (dst) by first doing fuzzy search from user supplied dst
    folder = dst  # working directory for fileDFG variable and factor operations
    dstname = dst # path name could either be legacy FileDFG dir or .tar.gz file of FileDFG files.
    unzip = false
    # add if doesn't have .tar.gz extension
    lastdirname = splitpath(dstname)[end]
    if !isdir(dst)
        unzip = true
        sdst = split(lastdirname, '.')
        if sdst[end] != "gz" #  length(sdst) == 1 &&
            dstname *= ".tar.gz"
            lastdirname *= ".tar.gz"
        end
    end
    # check the file actually exists
    @assert isfile(dstname) "cannot find file $dstname"
    # TODO -- what if it is not a tar.gz but classic folder instead?
    # do actual unzipping
    filename = lastdirname[1:(end-length(".tar.gz"))] |> string
    if unzip
        @show sfolder = split(dstname, '.')
        Base.mkpath(loaddir)
        folder = joinpath(loaddir, filename) #splitpath(string(sfolder[end-2]))[end]
        @info "loadDFG! detected a gzip $dstname -- unpacking via $loaddir now..."
        Base.rm(folder, recursive=true, force=true)
        # unzip the tar file
        run(`tar -zxf $dstname -C $loaddir`)
    end
    # extract the factor graph from fileDFG folder
    variables = DFGVariable[]
    factors = DFGFactor[]
    varFolder = "$folder/variables"
    factorFolder = "$folder/factors"
    # Folder preparations
    !isdir(folder) && error("Can't load DFG graph - folder '$folder' doesn't exist")
    !isdir(varFolder) && error("Can't load DFG graph - folder '$varFolder' doesn't exist")
    !isdir(factorFolder) && error("Can't load DFG graph - folder '$factorFolder' doesn't exist")

    varFiles = readdir(varFolder)
    factorFiles = readdir(factorFolder)
    for varFile in varFiles
        open("$varFolder/$varFile") do io
            packedData = JSON2.read(io, Dict{String, Any})
            push!(variables, unpackVariable(dfgLoadInto, packedData))
        end
    end
    @info "Loaded $(length(variables)) variables - $(map(v->v.label, variables))"
    @info "Inserting variables into graph..."
    # Adding variables
    map(v->addVariable!(dfgLoadInto, v), variables)

    for factorFile in factorFiles
        open("$factorFolder/$factorFile") do io
            packedData = JSON2.read(io, Dict{String, Any})
            push!(factors, unpackFactor(dfgLoadInto, packedData))
        end
    end
    @info "Loaded $(length(variables)) factors - $(map(f->f.label, factors))"
    @info "Inserting factors into graph..."
    # # Adding factors
    map(f->addFactor!(dfgLoadInto, f), factors)

    # Finally, rebuild the CCW's for the factors to completely reinflate them
    # NOTE CREATES A NEW DFGFactor IF  CCW TYPE CHANGES
    @info "Rebuilding CCW's for the factors..."
    for factor in factors
        rebuildFactorMetadata!(dfgLoadInto, factor)
    end

    # remove the temporary unzipped file
    if unzip
      @info "DFG.loadDFG! is deleting a temp folder created during unzip, $loaddir"
      # need this because the number of files created in /tmp/caesar/random is becoming redonkulous.
      Base.rm(loaddir, recursive=true, force=true)
    end

    return dfgLoadInto
end

# to be extended by users with particular choices in dispatch.
function loadDFG end


text:
"""
    eom_logistic2(dx, x, p, n) -> function

Equations of motions for a system consisting of two coupled logistic maps where
X unidirectionally influences Y [1].

## Equations of motion

The equations of motion are

```math
\\begin{aligned}
x(t+1) &= r_1 x(t)(1 - x(t)) \\\\
y(t+1) &= r_2 f(x,y)(1 - f(x,y)),
\\end{aligned}
```

with

```math
\\begin{aligned}
f(x,y) = \\dfrac{y + \\frac{c_{xy}(x \\xi )}{2}}{1 + \\frac{c_{xy}}{2}(1+ \\sigma )}
\\end{aligned}
```

The parameter `c_xy` controls how strong the dynamical forcing is. If `σ > 0`,
dynamical noise masking the influence of  `x` on `y` equivalent to
``\\sigma \\cdot \\xi`` is added at each iteration. Here,``\\xi`` is a draw from a
flat distribution on ``[0, 1]``. Thus, setting `σ = 0.05` is equivalent to
add dynamical noise corresponding to a maximum of ``5 \\%`` of the possible
range of values of the logistic map.

1. Diego, David, Kristian Agasøster Haaga, and Bjarte Hannisdal. "Transfer entropy computation
    using the Perron-Frobenius operator." Physical Review E 99.4 (2019): 042212.
"""

function eom_logistic2_unidir(dx, x, p, n)
    c_xy, r₁, r₂, σ = (p...,)
    ξ = rand() # random number from flat distribution on [0, 1]
    x, y = x[1], x[2]
    f_xy = (y +  (c_xy*(x + σ*ξ)/2) ) / (1 + (c_xy/2)*(1+σ))

    dx[1] = r₁ * x * (1 - x)
    dx[2] = r₂ * (f_xy) * (1 - f_xy)
    return
end

function logistic2_unidir(u₀, c_xy, r₁, r₂, σ)
    p = [c_xy, r₁, r₂, σ]
    DiscreteDynamicalSystem(eom_logistic2_unidir, u₀, p)
end

"""
    logistic2(;u₀ = rand(2), c_xy = 0.1, σ = 0.05,
        r₁ = 3.78, r₂ = 3.66) -> DiscreteDynamicalSystem

Initialise a system consisting of two coupled logistic maps where X
unidirectionally influences Y. By default, the parameters `r₁` and `r₂` are set
to values yielding chaotic behaviour.

## Equations of motion

The equations of motion are

```math
\\begin{aligned}
x(t+1) &= r_1 x(t)(1 - x(t)) \\\\
y(t+1) &= r_2 f(x,y)(1 - f(x,y)),
\\end{aligned}
```

with

```math
\\begin{aligned}
f(x,y) = \\dfrac{y + \\frac{c_{xy}(x \\xi )}{2}}{1 + \\frac{c_{xy}}{2}(1+ \\sigma )}
\\end{aligned}
```

The parameter `c_xy` controls how strong the dynamical forcing is. If `σ > 0`,
dynamical noise masking the influence of  `x` on `y` equivalent to
``\\sigma \\cdot \\xi`` is added at each iteration. Here,``\\xi`` is a draw from a
flat distribution on ``[0, 1]``. Thus, setting `σ = 0.05` is equivalent to
add dynamical noise corresponding to a maximum of ``5 \\%`` of the possible
range of values of the logistic map.

## References

1. Diego, David, Kristian Agasøster Haaga, and Bjarte Hannisdal. "Transfer entropy computation
    using the Perron-Frobenius operator." Physical Review E 99.4 (2019): 042212.
"""
logistic2_unidir(;u₀ = rand(2), c_xy = 0.1, r₁ = 3.78, r₂ = 3.66, σ = 0.05) =
    logistic2_unidir(u₀, c_xy, r₁, r₂, σ)

#To get chaotic realisation, check that the orbit doesn't settle to a few unique values
function good_logistic_unidir_trajectory(npts::Int;
        Ttr = 1000, dt = 1,
        c_xy = 0.5,
        Dr₁ = Uniform(3.6, 4.0),
        Dr₂ = Uniform(3.6, 4.0),
        σ = 0.0,
        n_maxtries = 300)

    n_tries = 0
    while n_tries <= n_maxtries
        s = logistic2_unidir(u₀ = rand(2),
            c_xy = c_xy,
            σ = σ,
            r₁ = rand(Dr₁),
            r₂ = rand(Dr₂))

        o = trajectory(s, npts * dt - 1, Ttr = Ttr, dt = dt)

        # Ensure there are not too many repeated values, so we don't have trivial behaviour

        if length(unique(o[:, 1])) > npts * 0.9 && length(unique(o[:, 2])) > npts * 0.9
            return o
        end

        n_tries += 1
    end
end


text:
########################################################################################################
## Details
# This file implements the sampled-based estimation proceedure for UQ inference problems.
#       We make heavy use of Turing.jl and Distributions.jl, many of the functions below are wrappers
#########################################################################################################
function sample(pdist :: SNAPDistribution{T}; num_samples = 600, return_type = "NamedTuple", verbose = true) where T<: AbstractFloat
    # Setup optimization utilities
    trans_x, x0 = namedtp_to_vec(pdist.x)
    n = length(x0)
    trans_p, params = namedtp_to_vec(pdist.p)
    f = distribution_wrapper(pdist, trans_x, trans_p)

    Turing.@model function regression(params)
        x ~ MvNormal(zeros(n), I(n))
        Turing.@addlogprob! -1.0*f(x, params)
    end

    model = regression(params)
    chain = Turing.sample(model, Turing.NUTS(0.65), num_samples)

    chain_mat = zeros(num_samples, n)
    for i = 1:n
        chain_mat[:, i] = chain[:, i, 1]
    end

    mean_chain = mean(chain_mat, dims = 1)
    pdist.x = transform(trans_x, vec(mean_chain))

    if return_type == "NamedTuple"
        posterior_samples = [ transform(trans_x, chain_mat[i, :]) for i = 1:num_samples]
        if verbose println("Saving   posterior samples as NamedTuple") end
        return posterior_samples
    elseif return_type == "Matrix"
        if verbose println("Returning samples as matrix, along with transform to NamedTuple") end
        return chain_mat, trans_x
    else
        if verbose println("Returning chains along with transform to NamedTuple") end
        return chain, trans_x
    end
end


text:
module TestUtils

import Wordle
import Test: @test, @test_throws

@test Wordle.to_index('a') === 1
@test Wordle.to_index('z') === 26

counts = Wordle.LetterCounts(Wordle.Candidate("alone"))
@test counts['a'] === 1
@test counts['l'] === 1
@test counts['o'] === 1
@test counts['n'] === 1
@test counts['e'] === 1
@test sum(counts.counts) === 5

end # module


text:
using TextAnalysis

#function to clean article
function cleanText(rawText)
    sd = StringDocument(rawText) # convert to SD
  try
    #remove cases
    prepare!(sd, strip_case)
    # remove punctuations
    prepare!(sd, strip_punctuation)
    # remove numbers
    prepare!(sd, strip_numbers )
    # remove preposition
    prepare!(sd, strip_prepositions)
    # remove prepos
    prepare!(sd, strip_pronouns)
    # remove stop words
    prepare!(sd, strip_stopwords)
    # remove articles
    prepare!(sd, strip_articles)
    # remove white space
    #  prepare!(sd, strip_whitespace)
    # remove corrupt
    remove_corrupt_utf8!(sd)
    # remove non letters
    # remove_nonletters!(sd)
    # remove custom stop words
     remove_words!(sd, ["et","al", "journal", "https", "org","pmid", "fig", "doi",
                        "nin", "na", "nand", "nof", "pone", "nn", "nthe", "ns",
                        "pmc", "nih", "author", "pa", "manuscript", "nanuscript",
                        "nu", "nscript", "nauthor", "nuthor", "nan", "page", "nmed",
                        "available","nint", "nm", #unit
                        "nam", "pubmed", "wwww", "http","nhttp", "www", "biomedcentral",
                        "com", "npage", "citation", "purposes",
                        "ip", "address", "national" ,"library",
                        "support", "center", "web", "site", "ncbi", "web","terms",
                        "copyright","copyright", "notice", "search", "databasepmc",
                        "content", "skip", "toolkit", "search",  "datasets", "geo",

                        "nplos", "based", "using", "uthor", "script", "uthor"])

  catch err
      # supress error
        warn(err)
  end
  sx = TextAnalysis.text(sd)
  remove_nonletters!(sd) # this brings error
  return(sd, sx)
end

max_stars_repo_name: jmulloka/Breast-Cancer-Text-Mining

max_issues_repo_name: jmulloka/Breast-Cancer-Text-Mining

max_forks_repo_name: jmulloka/Breast-Cancer-Text-Mining


text:
import DataFrames.skipmissing

sim_param = setup_sim_param_model()

include("misc_functions.jl")

##### To load the DR25 KOI catalog and apply the necessary cuts:

Q1Q17_DR25 = CSV.read("q1_q17_dr25_koi.tab_selectcols_new.csv", header=19, allowmissing=:all)

#If using our own cuts applied to the stellar and koi catalogs:

#=
Q1Q17_DR25_stellar = CSV.read("q1_q17_dr25_stellar_koi.tab_selectcols.csv", header=30, allowmissing=:all)
Q1Q17_DR25_stellar_all = CSV.read("q1_q17_dr25_stellar_koi.tab_all.csv", header=1, allowmissing=:all)

N_Kepler_targets = sum((Q1Q17_DR25_stellar_all[:teff] .> 4000.) .& (Q1Q17_DR25_stellar_all[:teff] .< 7000.) .& (Q1Q17_DR25_stellar_all[:logg] .> 4.))
println("Total number of Kepler targets satisfying our cuts: ", N_Kepler_targets)

table_confirmed = Q1Q17_DR25[(Q1Q17_DR25[:koi_disposition] .== "CONFIRMED") .| (Q1Q17_DR25[:koi_disposition] .== "CANDIDATE"), :] #Table containing only the confirmed and candidate objects
table_stellar = Q1Q17_DR25_stellar

table_confirmed = table_confirmed[(table_confirmed[:koi_period] .> get_real(sim_param,"min_period")) .& (table_confirmed[:koi_period] .< get_real(sim_param,"max_period")), :] #to make additional cuts in period P to be comparable to our simulated sample
table_confirmed = table_confirmed[(table_confirmed[:koi_prad] .> get_real(sim_param,"min_radius")/ExoplanetsSysSim.earth_radius) .& (table_confirmed[:koi_prad] .< get_real(sim_param,"max_radius")/ExoplanetsSysSim.earth_radius) .& (.~ismissing.(table_confirmed[:koi_prad])), :] #to make additional cuts in planetary radii to be comparable to our simulated sample

#To make cuts based on stellar properties of T_eff and logg:
teff_confirmed = zeros(Int64, size(table_confirmed,1)) #list to be filled with the T_eff (K) for the objects
logg_confirmed = zeros(Float64, size(table_confirmed,1)) #list to be filled with the logg(cgs) for the objects
cdpp5_confirmed = zeros(Float64, size(table_confirmed,1)) #list to be filled with the RMS CDPP 5h values for the objects
cdpp6_confirmed = zeros(Float64, size(table_confirmed,1)) #list to be filled with the RMS CDPP 6h values for the objects
for (i,KepID) in enumerate(table_confirmed[:kepid])
    teff_confirmed[i] = table_stellar[:teff][table_stellar[:kepid] .== KepID, :][1]
    logg_confirmed[i] = table_stellar[:logg][table_stellar[:kepid] .== KepID, :][1]
    cdpp5_confirmed[i] = table_stellar[:rrmscdpp05p0][table_stellar[:kepid] .== KepID, :][1]
    cdpp6_confirmed[i] = table_stellar[:rrmscdpp06p0][table_stellar[:kepid] .== KepID, :][1]
end

cdpp_cut = 250.
println("Fraction of CONFIRMED and CANDIDATE planets after CDPP cut: ", size(table_confirmed[(teff_confirmed .> 4000.) .& (teff_confirmed .< 7000.) .& (logg_confirmed .> 4.) .& (cdpp5_confirmed .< cdpp_cut), :], 1), "/", size(table_confirmed[(teff_confirmed .> 4000.) .& (teff_confirmed .< 7000.) .& (logg_confirmed .> 4.), :], 1))

table_confirmed = table_confirmed[(teff_confirmed .> 4000.) .& (teff_confirmed .< 7000.) .& (logg_confirmed .> 4.) .& (cdpp5_confirmed .< cdpp_cut), :]
=#

#If using a stellar table with cuts already made and just matching kepids to get a koi catalog:

stellar_catalog = setup_star_table_christiansen(sim_param)
N_Kepler_targets = get_int(sim_param,"num_kepler_targets")

table_confirmed = Q1Q17_DR25[(Q1Q17_DR25[:koi_disposition] .== "CONFIRMED") .| (Q1Q17_DR25[:koi_disposition] .== "CANDIDATE"), :] #Table containing only the confirmed and candidate objects

in_stellar_catalog = [] #will be filled with booleans indicating whether each koi in 'table_confirmed' is found in the 'stellar_catalog' or not
for i in 1:length(table_confirmed[:kepid])
    if any(x->x==table_confirmed[i,:kepid], stellar_catalog[:kepid])
        push!(in_stellar_catalog, true)
    else
        push!(in_stellar_catalog, false)
    end
end
in_stellar_catalog_indices = findall(in_stellar_catalog)

table_confirmed = table_confirmed[in_stellar_catalog_indices, :]

for (i,kepid) in enumerate(table_confirmed[:kepid]) #to replace the stellar and planetary radii in 'table_confirmed' with the more reliable values as derived from the stellar properties in 'stellar_catalog'
    stellar_radii_new = stellar_catalog[stellar_catalog[:kepid] .== kepid, :radius][1]
    table_confirmed[i, :koi_srad] = stellar_radii_new
    table_confirmed[i, :koi_prad] = (1 ./ExoplanetsSysSim.earth_radius)*stellar_radii_new*sqrt(table_confirmed[i, :koi_depth]/(1e6))
end

table_confirmed = table_confirmed[(table_confirmed[:koi_period] .> get_real(sim_param,"min_period")) .& (table_confirmed[:koi_period] .< get_real(sim_param,"max_period")), :] #to make additional cuts in period P to be comparable to our simulated sample
table_confirmed = table_confirmed[(table_confirmed[:koi_prad] .> get_real(sim_param,"min_radius")/ExoplanetsSysSim.earth_radius) .& (table_confirmed[:koi_prad] .< get_real(sim_param,"max_radius")/ExoplanetsSysSim.earth_radius) .& (.~ismissing.(table_confirmed[:koi_prad])), :] #to make additional cuts in planetary radii to be comparable to our simulated sample

#=
#If we want to write the stellar radii to a file:
CSV.write("stellar_radii_q1_q17_dr25_gaia_fgk.txt", stellar_catalog[[:kepid, :radius]])
=#

#=
#If we want to write the kep_oi's of the remaining koi catalog to a file:
f = open("kepoi_names.txt", "w")
println(f, table_confirmed[:kepoi_name])
close(f)
=#

#####To compute arrays of the observables (multiplicities, periods, period ratios, transit durations, transit depths, period-normalized transit duration ratios (xi), and transit depth ratios) from the remaining sample of planets:

KOI_systems = [x[1:6] for x in table_confirmed[:kepoi_name]]
checked_bools = zeros(size(table_confirmed,1)) #0's denote KOI that were not checked yet; 1's denote already checked KOI

M_confirmed = Int64[] #list to be filled with the planet multiplicities of the systems
R_confirmed = Float64[] #list to be filled with period ratios of adjacent planet pairs
xi_confirmed = Float64[] #list to be filled with the period-normalized transit duration ratios of adjacent planet pairs
D_ratio_confirmed = Float64[] #list to be filled with the transit depth ratios of adjacent planet pairs
P_confirmed = table_confirmed[:koi_period] #array of the periods (days)
P_confirmed = collect(skipmissing(P_confirmed))
t_D_confirmed = table_confirmed[:koi_duration] #array of the transit durations (hrs)
t_D_confirmed = collect(skipmissing(t_D_confirmed))
D_confirmed = table_confirmed[:koi_depth]/(1e6) #array of the transit depths (fraction)
D_confirmed = D_confirmed[ismissing.(D_confirmed) .== false] #to get rid of NA values

D_above_confirmed = Float64[] #list to be filled with the transit depths of planets above the photoevaporation boundary in Carrera et al 2018
D_below_confirmed = Float64[] #list to be filled with the transit depths of planets below the boundary
D_ratio_above_confirmed = Float64[] #list to be filled with the transit depth ratios of adjacent planet pairs, both above the boundary
D_ratio_below_confirmed = Float64[] #list to be filled with the transit depth ratios of adjacent planet pairs, both below the boundary
D_ratio_across_confirmed = Float64[] #list to be filled with the transit depth ratios of adjacent planet pairs, across the boundary

for i in 1:length(KOI_systems)
    if checked_bools[i] == 0 #if the KOI has not been checked (included while looking at another planet in the same system)
        system_i = (1:length(KOI_systems))[KOI_systems .== KOI_systems[i]]
        checked_bools[system_i] .= 1

        #To get the periods and transit durations in this system:
        system_P = table_confirmed[:koi_period][system_i] #periods of all the planets in this system
        system_t_D = table_confirmed[:koi_duration][system_i] #transit durations of all the planets in this system
        system_D = table_confirmed[:koi_depth][system_i] #transit depths (in ppm) of all the planets in this system
        system_radii = table_confirmed[:koi_prad][system_i] #radii of all the planets in this system
        system_sort_i = sortperm(system_P) #indices that would sort the periods of the planets in this system
        system_P = system_P[system_sort_i] #periods of all the planets in this system, sorted
        system_t_D = system_t_D[system_sort_i] #transit durations of all the planets in this system, sorted by period
        system_D = system_D[system_sort_i] #transit depths of all the planets in this system, sorted by period
        system_radii = system_radii[system_sort_i] #radii of all the planets in this system, sorted by period

        #To count the total number of planets in this system:
        push!(M_confirmed, length(system_P))

        #To compute the period ratios, period-normalized transit duration ratios, and transit depth ratios in this system:
        system_R = system_P[2:end]./system_P[1:end-1] #period ratios of all the adjacent planet pairs in this system
        system_D_ratio = system_D[2:end]./system_D[1:end-1] #transit depth ratios of all the adjacent planet pairs in this system
        system_xi = (system_t_D[1:end-1]./system_t_D[2:end]).*(system_P[2:end]./system_P[1:end-1]).^(1//3) #period-normalized transit duration ratios of all the adjacent planet pairs in this system

        append!(R_confirmed, system_R)
        append!(D_ratio_confirmed, system_D_ratio)
        append!(xi_confirmed, system_xi)

        #To separate the planets in the system as above and below the boundary:
        system_above_bools = [photoevap_boundary_Carrera2018(system_radii[x], system_P[x]) for x in 1:length(system_P)]
        #if length(system_above_bools) > 1 println(system_above_bools) end

        #To record the transit depths of the planets above and below the boundary:
        for (j,D) in enumerate(system_D)
            if system_above_bools[j] == 1
                append!(D_above_confirmed, D)
            elseif system_above_bools[j] == 0
                append!(D_below_confirmed, D)
            end
        end

        #To record the transit depth ratios of the planets above, below, and across the boundary:
        for (j,D_ratio) in enumerate(system_D_ratio)
            if system_above_bools[j] + system_above_bools[j+1] == 2 #both planets are above the boundary
                append!(D_ratio_above_confirmed, D_ratio)
            elseif system_above_bools[j] + system_above_bools[j+1] == 1 #one planet is above, the other is below the boundary
                append!(D_ratio_across_confirmed, D_ratio)
            elseif system_above_bools[j] + system_above_bools[j+1] == 0 #both planets are below the boundary
                append!(D_ratio_below_confirmed, D_ratio)
            end
        end
    end
end


text:
######### StanVariational Bernoulli example  ###########

using StanVariational, Test

bernoulli_model = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

bernoulli_data = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

# Use tmpdir across multiple runs to prevent re-compilation
#tmpdir = joinpath(@__DIR__, "tmp")

sm = VariationalModel("variational", bernoulli_model)

(sample_file, log_file) = stan_variational(sm; data=bernoulli_data)

if sample_file !== Nothing

  (chns, cnames) = read_variational(sm)

  # Show the same output in DataFrame format
  sdf = read_summary(sm)

  # Retrieve mean value of theta from the summary
  sdf[:theta, :mean]
  @test sdf[:theta, :mean][1] ≈ 0.33 atol=0.1

end


text:
using PeaceVote
using PeaceFounder
using PeaceFounder: BraiderConfig, CertifierConfig, BraidChainConfig, SystemConfig, Port

#using Sockets: @ip_str

### First we need to make a DemeSpec file

name = "PeaceDeme"

crypto = quote
G = CryptoGroups.Scep256k1Group()
hash(x::AbstractString) = parse(BigInt,Nettle.hexdigest("sha256",x),base=16)

Signer() = CryptoSignatures.Signer(G)
Signature(x::AbstractString,signer) = CryptoSignatures.DSASignature(hash(x),signer)
verify(data,signature) = CryptoSignatures.verify(signature,G) && hash(data)==signature.hash ? hash("$(signature.pubkey)") : nothing

(Signer,Signature,verify,hash)
end

deps = Symbol[:Nettle,:CryptoGroups,:CryptoSignatures]

notary = Notary(crypto,deps)
demespec = DemeSpec(name,crypto,deps,:PeaceFounder,notary)
save(demespec)

### The next step is to configure the server

# Maintainer sends the demespec file to the server. The server generates a server key pair and returns ID. Also configures the server such that it automatically starts to serve demespec file. That makes the server to stay in the listening state which allows maintainer to upload the configuration file specifying how the machine should work.

# SSH = "pi@192.1.1.1"
# SERVER_ID = configure(SSH,demespec)

deme = Deme(demespec,nothing)
server = Signer(deme,"server")

SERVER_ID = server.invoke.id

# Now maintainer finishes the configuration

maintainer = PeaceVote.Signer(deme,"maintainer")
MAINTAINER_ID = maintainer.invoke.id

MIXER_ID = (deme.invoke.spec.uuid,SERVER_ID) # Self mixing
CA_ID = (deme.invoke.spec.uuid,SERVER_ID)
TOOKEN_CA = MAINTAINER_ID

N = 3 # Number of participants per ballot/braid

### For local tests

MIXER_PORT = 2001 # Self mixing
BRAIDER_PORT = 2000
REGISTRATOR_PORT = 2002
VOTING_PORT = 2003
PROPOSAL_PORT = 2004
SYNC_PORT = 2005
TOOKEN_PORT = 2006
CERTIFIER_PORT = 2007

### For some remote action

# MIXER_PORT = Port(SERVER_IP,2001) # Self mixing
# BRAIDER_PORT = Port(SERVER_IP,2000)
# REGISTRATOR_PORT = Port(SERVER_IP,2002)
# VOTING_PORT = Port(SERVER_IP,2003)
# PROPOSAL_PORT = Port(SERVER_IP,2004)
# SYNC_PORT = Port(SERVER_IP,2005)

certifier = CertifierConfig(TOOKEN_CA,SERVER_ID,TOOKEN_PORT,CERTIFIER_PORT)
braider = BraiderConfig(BRAIDER_PORT,MIXER_PORT,N,SERVER_ID,MIXER_ID)
braidchain = BraidChainConfig(MAINTAINER_ID,[CA_ID,],SERVER_ID,REGISTRATOR_PORT,VOTING_PORT,PROPOSAL_PORT)

systemconfig = SystemConfig(MIXER_PORT,SYNC_PORT,certifier,braider,braidchain)

PeaceFounder.save(systemconfig,maintainer)

### After that one may upload the config file to the server.
# upload(deme,SERVER_IP)

max_stars_repo_name: akels/DemocracyCredit

max_issues_repo_name: akels/DemocracyCredit

max_forks_repo_name: akels/DemocracyCredit


text:
module PushFilterTest1

using HPAT

#HPAT.CaptureAPI.set_debug_level(3)
#HPAT.DomainPass.set_debug_level(3)
#ParallelAccelerator.set_debug_level(3)
#using CompilerTools
#CompilerTools.OptFramework.set_debug_level(3)

@acc hpat function test1(file_name, file_name2)
    t1 = DataSource(DataTable{:userid=Int64, :val2=Float64}, HDF5, file_name)
    t2 = DataSource(DataTable{:userid=Int64, :val3=Float64}, HDF5, file_name2)
    t3 = join(t1,t2, :userid==:userid, :userid)
    t4 = t3[:val2>2.0]
    return t4[:userid], t4[:val3]
end

using HDF5
using MPI
using Base.Test

function main()

if MPI.Comm_rank(MPI.COMM_WORLD)==0
    h5write("test1_1.hdf5", "/userid", [1,2,3,1,2])
    h5write("test1_1.hdf5", "/val2", [1.1,2.1,3.1,3.2,1.9])
    h5write("test1_2.hdf5", "/userid", [1,3])
    h5write("test1_2.hdf5", "/val3", [7.1,8.3])
end
a,b = test1("test1_1.hdf5","test1_2.hdf5")
println(a)
println(b)

@test a==[1,3]
@test_approx_eq b [7.1,8.3]

if MPI.Comm_rank(MPI.COMM_WORLD)==0
    rm("test1_1.hdf5")
    rm("test1_2.hdf5")
end

end

end

PushFilterTest1.main()


text:
export @theory!, @term, @rules, @rewrite

const THEORY = Dict{Σ,AbstractTheory}()

macro theory!(theories)
    @assert isa(theories, Expr) && theories.head === :block
    Base.remove_linenums!(theories)

    for line ∈ theories.args
        @assert line.head === :call && line.args[1] === :(=>)
        THEORY[line.args[2]] = Base.eval(__module__, line.args[3])
    end

    nothing
end

function _to_theory(th, expr; strict=false)
    (root, args) = if isa(expr, Symbol)
        (expr, [])
    elseif isa(expr, Expr)
        expr.head === :call || throw(ArgumentError("invalid non-call Expr"))
        (expr.args[1], expr.args[2:end])
    elseif isa(expr, Union{Variable,AbstractTerm})
        return expr
    else
        throw(ArgumentError("invalid expression: $(repr(expr))"))
    end

    root_theory = haskey(th, root) ? th[root] :
        strict ? throw(ArgumentError("$root undefined in theory")) : FreeTheory()

    for i ∈ eachindex(args)
        args[i] = _to_theory(th, args[i]; strict=strict)
    end

    term(root_theory, root, args)
end

function Base.show(io::IO, t::Union{Variable, AbstractTerm})
    term_expr = convert(Expr, t)
    Base.show_call(io, :call, Symbol("@term"), [term_expr], 0)
end

macro term(expr)
    :(_to_theory(THEORY, $(Meta.quot(expr))))
end

function _rules_replace!(ex, vars)
    isa(ex, Symbol) && return get(vars, ex, ex)

    if isa(ex, Expr)
        for i ∈ eachindex(ex.args)
            ex.args[i] = _rules_replace!(ex.args[i], vars)
        end
    end
    return ex
end
function _clean(ex, vars)
    ex′ = _rules_replace!(ex, vars)
    _to_theory(THEORY, ex′; strict=true)
end
macro rules(name, varnames, body)
    @assert isa(name, Symbol)
    @assert isa(varnames, Expr) && varnames.head === :vect
    @assert all(isa(x, Symbol) for x ∈ varnames.args)

    @assert isa(body, Expr) && body.head === :block
    Base.remove_linenums!(body)

    vars = Dict(x => Variable() for x ∈ varnames.args)

    rw = Rewriter()
    for i ∈ eachindex(body.args)
        line = body.args[i]
        @assert line.head === :(:=)
        lhs, rhs = _clean.(line.args, Ref(vars))
        push!(rw, lhs => replace(rhs))
    end

    fn, expr = compile(rw)

    struct_name = gensym(name)
    return quote
        $expr
        $(esc(name))(t) = $fn(t)
    end
end

macro rewrite(rw, expr)
    :($(esc(rw))(_to_theory(THEORY, $(Meta.quot(expr)))))
end


text:
using LinearAlgebra
using SparseArrays

# using Flux

⊗ = kron

# LITTLE ENDIAN !
# 1st qubit = far right

# if IittleEndian then the last qubit is the control
function Controlled(gate::AbstractMatrix)
  sz = size(gate)
  n = sz[1]
  ret = spzeros(ComplexF64, 2*n, 2*n)
  ret[1:n,1:n] += I
  ret[n+1:2*n,n+1:2*n] = gate
  return ret
end

# returns sparse permutation matrix P of width 2^pad
# input indices of qubits are 1-based
# if G acts on qubits 1,2,... then P^(-1)*G*P acts on ind[1], ind[2],...
# LittleEndian!
function Permute(ind::Array{Int}; pad::Int=0)
  # assert maximum(ind) <= pad
  # assert ind[i] all different

  k = length(ind)
  pad = max(pad, maximum(ind))

  # reconstruct full permutation sigma
  sigma = append!(copy(ind), zeros(Int64, pad-k))
  # rest = setdiff(1:pad, ind)
  rest = filter(x -> !(x in ind), 1:pad)
  sigma[k+1:pad] = rest
  # todo: more stable sigma?

  # println(sigma)

  rows = zeros(Int64, 2^pad)
  cols = zeros(Int64, 2^pad)
  # vals = zeros(ComplexF64, 2^pad)
  vals = zeros(Int64, 2^pad)

  function permute_binary(b::Int64)
    ret = 0

    b -= 1

    for i=1:pad
      digit = b & 1
      ret += digit << (sigma[i]-1)
      b = b >> 1
    end

    return ret+1
  end

  for b in 1:2^pad
    rows[b] = b
    cols[b] = permute_binary(b)
    vals[b] = 1
  end

  # println("cols", cols)

  return sparse(rows, cols, vals)
end

function ApplyGateOn(gate::AbstractMatrix, ind::Array{Int}; pad::Int=0)
  # assert gate is acting on length(ind) qubits

  k = length(ind)
  pad = max(pad, maximum(ind))

  # w = size(gate)[1]
  P = Permute(ind, pad=pad)

  ret = gate
  for i in k+1:pad
    ret = [1 0; 0 1] ⊗ ret
  end

  return P'*ret*P
end

# https://docs.microsoft.com/en-us/qsharp/api/qsharp/microsoft.quantum.intrinsic.rx
function Rx(θ::Real)
  return [cos(θ/2) -im*sin(θ/2); -im*sin(θ/2) cos(θ/2)]
end

# https://docs.microsoft.com/en-us/qsharp/api/qsharp/microsoft.quantum.intrinsic.ry
function Ry(θ::Real)
  return [cos(θ/2) -sin(θ/2); sin(θ/2) cos(θ/2)]
end

# https://docs.microsoft.com/en-us/qsharp/api/qsharp/microsoft.quantum.intrinsic.rz
function Rz(θ::Real)
  return [exp(-im*θ/2) 0; 0 exp(im*θ/2)]
end

# https://docs.microsoft.com/en-us/qsharp/api/qsharp/microsoft.quantum.machinelearning.controlledrotation
# little endian!
# function ControlledRotation(ind::Tuple{Int, Array{}}, Pauli::Function,  param::Real)
function ControlledRotation(target::Int, control::Int, Pauli::Function, param::Real; pad::Int=0)
  # indices are 1 based
  # id = [1 0; 0 1]

  if target==control
    return ApplyGateOn( Pauli(param), [target], pad=pad)
    # return Rz(param)
  else
    C = Controlled(Pauli(param))
    return ApplyGateOn( C, [target, control], pad=pad)
  end

end

max_stars_repo_name: dandanua/QML

max_issues_repo_name: dandanua/QML

max_forks_repo_name: dandanua/QML


text:
import Term.layout: pad, vLine, hLine
import Term: Panel

println("pad")
@time pad("aaa", 20, :left)
@time pad("aaa", 20, :right)
@time pad("aaa", 20, :center)

@time pad("aaa", 10, 20)
p = Panel(; width=20, height=10)
@time pad(p.segments, 10, 10);

println("string stacking")
s1 = "123"
s2 = "12345"

@time s1 * s2
@time s1 / s2

println("renderables stacking")
p2 = Panel(; width=5, height=12)
@time p * p2
@time p / p2

println("vLine")
@time vLine(10; style="red")
@time vLine(p2)

println("hLine")
@time hLine(10; style="red")
@time hLine(50, "title")

# TODO test hLine and vLine sizes correct


text:
#TODO: example with large deviation simulation
#TODO: update example of contact process with the implementation from sonnet
using MonteCarloX
using StatsBase
using Random
using LightGraphs
using Plots

"""
main program that implements full example simulation of 2D Ising model and compares to the Beale solution
call this from within julia -O3
@time main(flag_plot=false) should run O(4s) comparable with C code from Martin?!
"""
function main(;flag_plot::Bool = true)
    println("Initialize system (8x8)")
    rng = MersenneTwister(1000)
    system = constructIsing([8,8], rng)

    if flag_plot
        p_E = plot(xaxis = ("beta"), yaxis = ("energy per spin"))
        p_C = plot(xaxis = ("beta"), yaxis = ("specific heat"))
        display(plot(p_E, p_C, layout = (2, 1)))
    end

    # eval Beale solution
    println("Evaluate exact solution")
    ana_beta = [beta for beta in 0.0:0.001:0.8]
    ana_E = []
    ana_C = []
    for beta in ana_beta
        E, C  = analytic(beta, log_dos_beale_8x8)
        push!(ana_E, E)
        push!(ana_C, C)
    end
    if flag_plot
        display(plot!(ana_beta, ana_E / 64.0, label = "Beale", subplot = 1))
        display(plot!(ana_beta, ana_C / 64.0, label = "Beale", subplot = 2))
    end

    println("Multicanonical initialization")
    alg = Multicanonical()
    E_max =  8*8
    E_min = -E_max
    dE = 4
    log_weight = Histogram(E_min:dE:E_max)
    histogram  = zero(log_weight)

    println("Multicanonical iteration")
    for i=1:100
      histogram.weights .= 0
      run_iteration(rng, system, log_weight, histogram, Int(1e5), 1000, rng)
      update_weights(alg, log_weight, histogram)
    end

    #run_production(rng, system, log_weight, histogram, Int(1e5), 1000, rng)
    #list_est_E = reweight()...
    ## plot energy compared to Beale
    #if flag_plot
    #    display(plot!(list_beta, list_est_E / 64.0, yerr = list_err_E / 64.0, seriestype = :scatter, label = "MC", subplot = 1))
    #    display(plot!(list_beta, list_est_C / 64.0, yerr = list_err_C / 64.0, seriestype = :scatter, label = "MC", subplot = 2))
    #end
end

"""
run simulation
"""
function run_iteration(rng::AbstractRNG, system::IsingSystem, log_weight::Histogram, histogram::Histogram, n_meas::Int64, n_therm::Int64)
    current_E = E(system)
    N = length(system.spins)
    for sweep in 1:n_meas + n_therm
        # sweep
        for step in 1:N
            dE = update_spin_flip(rng, system, log_weight, current_E)
            current_E += dE
            push!(histogram,current_E)
        end
        #current_E = E(system)
    end
end

function update_spin_flip(rng::AbstractRNG, system::IsingSystem, log_weight::Histogram, E_old::Int)::Int
    index = rand(rng, 1:length(system.spins))
    dE    = -2 * E_local(system, index)
    if accept(Multicanonical(), rng, log_weight, E_new+dE, E_old)
        system.spins[index] *= -1
    else
        dE = 0
    end
    return dE
end

function mean_err(X::Vector{Float64})::Tuple{Float64,Float64}
    mean  = 0.0
    mean2 = 0.0
    for x in X
        mean  += x
        mean2 += x * x
    end
    N = length(X)
    mean  /= N
    mean2 /= N
    return mean, sqrt((mean2 - mean) / (N - 1))
end

###########################################################################
struct IsingSystem{F}
    dims::Vector{Int}
    lattice::SimpleGraph
    spins::Vector{Int}
    nearest_neighbors::F
end

function constructIsing(dims::Vector{Int}, rng::AbstractRNG)::IsingSystem
    lattice = LightGraphs.SimpleGraphs.grid(dims, periodic = true)
    spins   = rand(rng, [-1,1], nv(lattice))
    nearest_neighbors(i) = outneighbors(lattice, i)
    return IsingSystem(dims, lattice, spins, nearest_neighbors)
end

function analytic(beta::Float64, log_dos::Vector{Tuple{Int,Float64}})::Tuple{Float64,Float64}
    mean_E  = 0.0
    mean_E2 = 0.0
    norm    = 0.0
    for (E, log_d) in log_dos
        mean_E  +=   E * exp(log_d - beta * E)
        mean_E2 += E^2 * exp(log_d - beta * E)
        norm    +=     exp(log_d - beta * E)
    end
    mean_E  /= norm
    mean_E2 /= norm
    return mean_E, beta^2 * (mean_E2 - mean_E^2)
end

# J=1
# E=-J*sum s_i s_j
function E(system::IsingSystem)::Int
    E = 0
    for i in 1:length(system.spins)
        E += E_local(system, i)
    end
    return E / 2
end

function E_local(system::IsingSystem, index::Int)::Int
    e = 0
    for j in system.nearest_neighbors(index)
        e -= system.spins[index] * system.spins[j]
    end
    return e
end

# logarithmic density of states for the 2D 8x8 Ising model from Bealse solution
log_dos_beale_8x8 = [(-128, 0.6931471805599453), (-120, 4.852030263919617), (-116, 5.545177444479562), (-112, 8.449342524508063), (-108, 9.793672686528922), (-104, 11.887298863200714), (-100, 13.477180596840947), (-96, 15.268195474147658), (-92, 16.912371686315282), (-88, 18.59085846191256), (-84, 20.230089202801466), (-80, 21.870810400320693), (-76, 23.498562234123614), (-72, 25.114602234581373), (-68, 26.70699035290573), (-64, 28.266152815389898), (-60, 29.780704423363996), (-56, 31.241053997806176), (-52, 32.63856452513369), (-48, 33.96613536105969), (-44, 35.217576663643314), (-40, 36.3873411250109), (-36, 37.47007844691906), (-32, 38.46041522581422), (-28, 39.35282710786369), (-24, 40.141667825183845), (-20, 40.82130289691285), (-16, 41.38631975325592), (-12, 41.831753810069756), (-8, 42.153328313883975), (-4, 42.34770636939425), (0, 42.41274640460084), (4, 42.34770636939425), (8, 42.153328313883975), (12, 41.831753810069756), (16, 41.38631975325592), (20, 40.82130289691285), (24, 40.141667825183845), (28, 39.35282710786369), (32, 38.46041522581422), (36, 37.47007844691906), (40, 36.3873411250109), (44, 35.217576663643314), (48, 33.96613536105969), (52, 32.63856452513369), (56, 31.241053997806176), (60, 29.780704423363996), (64, 28.266152815389898), (68, 26.70699035290573), (72, 25.114602234581373), (76, 23.498562234123614), (80, 21.870810400320693), (84, 20.230089202801466), (88, 18.59085846191256), (92, 16.912371686315282), (96, 15.268195474147658), (100, 13.477180596840947), (104, 11.887298863200714), (108, 9.793672686528922), (112, 8.449342524508063), (116, 5.545177444479562), (120, 4.852030263919617), (128, 0.6931471805599453)];


text:
@testset "JS_WG_HNE_DL" begin
    θ₀ = [0.7, 0.9, 3.3, 5.0, 0.3, 0.4]
    ω = 1.1
    @testset "Preallocation" begin
        model = JS_WG_HNE_DL{1,40}(θ₀)
        @test model.r_over4 == model.r/4
        @test model.ωₚ² == model.ωₚ^2
        @test model.ωₚ³ == model.ωₚ^3
        @test model.ωₚ⁴ == model.ωₚ^4
        @test model.logγ == log(model.γ)
        @test model.cosϕₘ == cos(model.ϕₘ)
        @test model.sinϕₘ == sin(model.ϕₘ)
        @test model.cos2ϕₘ == cos(2model.ϕₘ)
        @test model.sin2ϕₘ == sin(2model.ϕₘ)
        @test model.invexp2σ² == exp(-2.0 * model.σ^2)
        @test model.invexphalfσ² == exp(-0.5 * model.σ^2)
        @test model.r_over4 == model.r/4
        @test model.ωₚ² == model.ωₚ^2
        @test model.ωₚ³ == model.ωₚ^3
        @test model.ωₚ⁴ == model.ωₚ^4
        @test model.logγ == log(model.γ)
        @test model.cosϕₘ == cos(model.ϕₘ)
        @test model.sinϕₘ == sin(model.ϕₘ)
        @test model.cos2ϕₘ == cos(2model.ϕₘ)
        @test model.sin2ϕₘ == sin(2model.ϕₘ)
    end
    @testset "Gradient" begin
        @test approx_gradient(θ -> sdf(JS_WG_HNE_DL{1,40}(θ), ω), θ₀)      ≈ grad_sdf(JS_WG_HNE_DL{1,40}(θ₀), ω)
        @test approx_gradient(θ -> sdf(JS_WG_HNE_DL{1,40}(θ), 0), θ₀)      ≈ grad_sdf(JS_WG_HNE_DL{1,40}(θ₀), 0)
    end
    @testset "coherancy" begin
        S = sdf(JS_WG_HNE_DL{1,40}(θ₀), ω)
        @test coherancy(JS_WG_HNE_DL{1,40}(θ₀), ω) ≈ [S[i,j]/sqrt(S[i,i]*S[j,j]) for i in 1:size(S,1), j in 1:size(S,2)]
    end
    @testset "Error handling" begin
        @test_throws ArgumentError JS_WG_HNE_DL{1,40}(-0.7,0.9, 3.3, 5.0, 0, 0.4)
        @test_throws ArgumentError JS_WG_HNE_DL{1,40}(0.7,-0.9, 3.3, 5.0, 0, 0.4)
        @test_throws ArgumentError JS_WG_HNE_DL{1,40}(0.7, 0.9, 0.8, 5.0, 0, 0.4) # 0.8 < 1
        @test_throws ArgumentError JS_WG_HNE_DL{1,40}(0.7, 0.9, 3.3, 0.9, 0, 0.4) # 0.9 < 1
        @test_throws ArgumentError JS_WG_HNE_DL{1,40}(0.7, 0.9, 3.3, 5.0, 0,-0.4)
        @test_throws ArgumentError JS_WG_HNE_DL{1,40}(ones(5))
        @test_throws ArgumentError JS_WG_HNE_DL{1,40}(ones(7))
        @test_throws ErrorException JS_WG_HNE_DL(ones(6))
        @test_throws ErrorException JS_WG_HNE_DL(0.7,0.9, 3.3, 5.0, 0, 0.4)
    end
end


text:
# Interface to libLLVMCore, which implements the LLVM intermediate representation as well as
# other related types and utilities.

# NOTE: we don't completely stick to the C API's organization here

include("core/context.jl")
include("core/type.jl")
include("core/value.jl")
include("core/metadata.jl")
include("core/module.jl")
include("core/basicblock.jl")
include("core/attributes.jl")
include("core/function.jl")
include("core/instructions.jl")


text:
type Subplots <: PlotWidget
  widget::PyObject
  plots::Vector{Plot}
  n::Int
  nrowsOverride::Int
  ncolsOverride::Int
  nrows::Int
  ncols::Int

  Subplots() = new(PLOT.SubplotWidget(), [], 0, 0, 0)
end

function updateGrid(sp::Subplots)
  if sp.nrowsOverride == 0
    if sp.ncolsOverride == 0
      sp.nrows = round(Int, sqrt(sp.n))
      sp.ncols = ceil(Int, sp.n / sp.nrows)
    else
      sp.ncols = sp.ncolsOverride
      sp.nrows = ceil(Int, sp.n / sp.ncols)
    end
  else
    sp.nrows = sp.nrowsOverride
    sp.ncols = ceil(Int, sp.n / sp.nrows)
  end
end

function getplot(sp::Subplots, c::Int)
  while c > length(sp.plots)
    push!(sp.plots, Plot())
    sp.n += 1
    updateGrid(sp)
  end
  sp.plots[c]
end

function refresh(sp::Subplots)
  for plt in sp.plots
    refresh(plt)
  end
  widgets = reshape([[plt.widget for plt in sp.plots] ; fill(nothing, sp.nrows*sp.ncols - sp.n)], sp.ncols, sp.nrows)'
  sp.widget[:addFigures](widgets)
  nothing
end

subplot(y; kvs...) = subplot(; y = y, kvs...)
subplot(x, y; kvs...) = subplot(; x = x, y = y, kvs...)

function subplot(; kvs...)
  sp = Subplots()
  resizewidget(sp, 800, 600)
  moveToLastScreen(sp)

  d = Dict(kvs)
  sp.nrowsOverride = get(d, :nrows, 0)
  sp.ncolsOverride = get(d, :ncols, 0)

  if !((:show, false) in kvs)
    push!(kvs, (:show, true))
  end

  oplot(sp; kvs...)
  sp
end


text:
using ConvexTests, Hypatia

@info "Starting Hypatia tests"

do_tests("Hypatia", () -> Hypatia.Optimizer(verbose=false); exclude = [r"mip"])

max_stars_repo_name: ericphanson/ConvexTests

max_issues_repo_name: ericphanson/ConvexTests

max_forks_repo_name: ericphanson/ConvexTests


text:
# m5.6s.jl

using Pkg, DrWatson

using MonteCarloMeasurements
using StanSample
using StatisticalRethinking

# ### snippet 5.29

df = CSV.read(sr_datadir("milk.csv"), DataFrame; delim=';');
df = filter(row -> !(row[:neocortex_perc] == "NA"), df);
df[!, :neocortex_perc] = parse.(Float64, df[:, :neocortex_perc])
df[!, :lmass] = log.(df[:, :mass])
scale!(df, [:kcal_per_g, :neocortex_perc, :lmass])

stan5_6 = "
data {
     int < lower = 1 > N; // Sample size
     vector[N] K; // Outcome
     vector[N] M; // Predictor
}

parameters {
     real a; // Intercept
     real bM; // Slope (regression coefficients)
     real < lower = 0 > sigma;    // Error SD
}

model {
      vector[N] mu;               // mu is a vector
      a ~ normal(0, 0.2);           //Priors
      bM ~ normal(0, 0.5);
      sigma ~ exponential(1);
      mu = a + bM * M;
      K ~ normal(mu , sigma);     // Likelihood
}
";

# Define the SampleModel and set the output format to :mcmcchains.

m5_6s = SampleModel("m5.6", stan5_6);

# Input data for cmdstan

m5_6_data = Dict("N" => size(df, 1), "M" => df[!, :lmass_s],
    "K" => df[!, :kcal_per_g_s]);

# Sample using StanSample

rc5_6s = stan_sample(m5_6s, data=m5_6_data);

if success(rc5_6s)

      part5_6s = read_samples(m5_6s, :particles)
      part5_6s |> display

      rethinking = "
                 mean   sd  5.5% 94.5%
        a      0.05 0.15 -0.20  0.29
        bM    -0.28 0.19 -0.59  0.03
        sigma  0.95 0.16  0.70  1.20
      "

end


text:
struct PersistentDisjointSet{T}
    table::PersistentHashMap{T,T}
    heights::PersistentHashMap{T,UInt8}
end
PersistentDisjointSet{T}() where {T} = PersistentDisjointSet{T}(PersistentHashMap{T,T}(),PersistentHashMap{T,UInt8}())

height(t::PersistentDisjointSet{T},i::T) where {T} =  get(t.heights,i,zero(UInt8))

function Base.getindex(t::PersistentDisjointSet{T},i::T) where {T}
    while haskey(t.table,i) # In a persistent setting, writes cost memory, so path compression may not be worth the allocations.
        i = t.table[i]      # Height balancing already guarentees finding the root in floor(lg(N)) lookups.
    end
    i
end

function Base.union(t::PersistentDisjointSet{T},i::T,j::T) where {T}
    i = t[i];j=t[j]
    if i == j return t end
    if height(t,i) > height(t,j)
        PersistentDisjointSet(assoc(t.table,j,i),t.heights)
    elseif height(t,i) < height(t,j)
        PersistentDisjointSet(assoc(t.table,i,j),t.heights)
    else
        PersistentDisjointSet(assoc(t.table,j,i),assoc(t.heights,i,height(t,i) + 1))
    end
end

function Base.show(io::IO,t::PersistentDisjointSet{T}) where {T}
    Base.print(io,"PersistentDisjointSet{$T}")
    Base.print(io,sort!(map(x->first(x) => t[last(x)],collect(pairs(t.table)));by=x->(last(x),first(x))))
end


text:
module Agents

using Distributed
using LightGraphs
using DataFrames
using Random
import Base.Iterators.product
import Base.iterate
import Base.length

# Core structures of Agents.jl
include("core/model.jl")
include("core/space.jl")
include("core/agent_space_interaction.jl")
include("core/continuous_space.jl")

# Stepping and data collection functionality
include("simulations/step.jl")
include("simulations/collect.jl")
include("simulations/paramscan.jl")
include("simulations/sample.jl")

end # module


text:
# Copyright Wen-Chih Chen
# The example determine the efficiencies for all DMUs based on the CRS input-oriented, CCR, model (Charnes et al., 1978)

using JuMP
using Gurobi # Gurobi is used as the LP solver
using FrontierEfficiencyAnalysis

# read data from a csv file
data = readcsv("example.csv") # input User's (.csv) data path
# each row is for a DMU; each column is an input or output
# the first two columns are the inputs while the others are the outputs
# determine the size of the data
scale, dimension = size(data) # scale is the number of DMU, dimension is the total number of inputs and outputs

# compute efficiencies for all DMUs
for t = 1 : scale
    ### Modeling section
    # Here is the CRS input-oriented model (CCR model) to evaluate DMU t
    model = Model(solver = GurobiSolver()) # Gurobi is used as the LP solver here. Users can choose their favorite solver.
    @variable(model, Lambda[1:scale] >= 0)
    @variable(model, Theta)
    @objective(model, Min, Theta)
    @constraint(model, inputCon[i=1:2], sum(Lambda[r]*data[r,i] for r = 1:scale) <= Theta*data[t,i])
    @constraint(model, outputCon[j=3:dimension], sum(Lambda[r]*data[r,j] for r = 1:scale) >= data[t,j])
    ## uncomment the following line to add the convexity constraint for the VRS model
    # @constraint(model, sum(Lambda[r] for r = 1:scale) == 1)

    ### Problem solving
    solveDEA(model)

    ### Display
    println("Results for DMU $t")
    println("The efficicncy: $(getobjectivevalue(model))")
    println("The efficicncy: $(getvalue(Theta))")
    println("The lambdas: $(getvalue(Lambda)))")

    ## Use the following if returning dual values and slacks is needed

    ### Problem solving
    # duals, slacks = solveDEA(model)
    ### Display
    # println("Results for DMU $t")
    # println("The efficicncy: $(getobjectivevalue(model))")
    # println("The efficicncy: $(getvalue(Theta))")
    # println("The lambdas: $(getvalue(Lambda)))")

    # println("The dual values (weights): $duals") # a vector associated with the constraints you define from the top to the bottom
    # println("the slack values: $slacks") # a vector associated with the constraints you define from the top to the bottom

end


text:
module Gomen

export Game, sparam, tparam, play, Games, quadrant, name
export SimpleGraph, GraphGenerator, nv, params
export CycleGraphGenerator, WheelGraphGenerator, StarGraphGenerator, GridGraphGenerator
export BarabasiAlbertGenerator, ErdosRenyiGenerator
export AbstractRule, apply, Sigmoid, Heaviside
export AbstractScheme, RandomScheme, CounterFactual, MimicRandom, MimicBest, MimicBiased, decide
export AbstractArena, game, graph, scheme, Arena, payoffs
export Aggregator, Sum, HarmonicMean
export Scores, Rescorer, GammaRescorer, CLRRescorer
export EdgeEvidence, Scorer, infer
export SymmetricScorer
export ChisqScorer, LaggedChisqScorer, SymLaggedChisqScorer
export MIScorer, LaggedMIScorer, SymLaggedMIScorer
export TEScorer, SymTEScorer
export SignificanceScorer
export PerformanceCurve, xlabel, ylabel, xvalues, yvalues, roc, auc
export ROC, PRC, tpr, fpr, recall, precision
export logloss, brier
export SimulationError, ScoringError, PerformanceError
export gomen

using Base.Iterators, Base.Meta, LinearAlgebra, Random, Statistics
using Distributions
using DrWatson
using HypothesisTests, Imogen
using LightGraphs, LightGraphs.SimpleGraphs
using MLBase
using RecipesBase
using Requires
using StaticArrays

include("core.jl")
include("games.jl")
include("graphs.jl")
include("inference.jl")
include("performance.jl")

struct SimulationError
    rng::AbstractRNG
    arena::Arena
    rounds::Int
    replicates::Int
    err
end

struct ScoringError
    rng::AbstractRNG
    arena::Arena
    series::Array{Int,3}
    scorer::Union{Nothing,Scorer}
    significance::Bool
    rescorer::Union{Nothing,Rescorer}
    err
end

struct PerformanceError
    arena::Arena
    series::Array{Int,3}
    scorer::Union{Nothing,Scorer}
    significance::Bool
    rescorer::Union{Nothing,Rescorer}
    scores::Vector{EdgeEvidence}
    err
end

function gomen(arena, scorer, significance, rescorer, rounds, replicates; kwargs...)
    rng = deepcopy(Random.default_rng())
    series = try
        play(arena; rounds, replicates)
    catch err
        throw(SimulationError(rng, arena, rounds, replicates, err))
    end

    scorerrng = deepcopy(Random.default_rng())
    scores = try
        if isnothing(rescorer)
            infer(significance ? SignificanceScorer(scorer) : scorer, series; kwargs...)
        else
            infer(significance ? SignificanceScorer(scorer) : scorer, rescorer, series; kwargs...)
        end
    catch err
        throw(ScoringError(scorerrng, arena, series, scorer, significance, rescorer, err))
    end

    roc, rauc, prc, pauc, logloss, brier = try
        roc = Gomen.roc(ROC, arena, scores)
        prc = MLBase.roc(PRC, arena, scores)
        roc, auc(roc), prc, auc(prc), Gomen.logloss(arena, scores), Gomen.brier(arena, scores)
    catch err
        throw(PerformanceError(arena, series, scorer, significance, rescorer, scores, err))
    end

    @dict game rng scorerrng series scorer significance rescorer roc rauc prc pauc logloss brier
end

function __init__()
    @require DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" begin
        include("dataframes.jl")
    end
end

end


text:
e_x, _ = noisy_sin(50; noise = 0.)
e_X = expand_poly(e_x, degree = 5)
df = DataFrame(A=rand(10), B=collect(1:10), C=[string(x) for x in 1:10])
df_na = deepcopy(df)
df_na[!, :A] = allowmissing(df_na[!, :A])
df_na[1, :A] = missing

@testset "Test expand_poly" begin
    @test size(e_X) == (5, 50)
end

@testset "Test center! and rescale!" begin
    # Center Vectors
    xa = copy(e_x)
    @test center!(xa) ≈ mean(e_x)
    @test abs(mean(xa)) <= 10e-10

    xa = copy(e_x)
    mu = mean(xa)
    center!(xa, mu, obsdim=1)
    @test abs(mean(xa)) <= 10e-10

    xa = copy(e_x)
    mu = vec(fill!(similar(xa),1))
    center!(xa, mu, obsdim=1)
    @test sum(e_x .- mean(xa)) ≈ length(mu)

    # Center Matrix w/o mu
    Xa = copy(e_X)
    center!(Xa)
    @test abs(sum(mean(Xa, dims=2))) <= 10e-10

    Xa = copy(e_X)
    center!(Xa, obsdim=1)
    @test abs(sum(mean(Xa, dims=1))) <= 10e-10

    Xa = copy(e_X)
    center!(Xa, ObsDim.First())
    @test abs(sum(mean(Xa, dims=1))) <= 10e-10

    Xa = copy(e_X)
    center!(Xa, obsdim=2)
    @test abs(sum(mean(Xa, dims=2))) <= 10e-10

    Xa = copy(e_X)
    center!(Xa, ObsDim.Last())
    @test abs(sum(mean(Xa, dims=2))) <= 10e-10

    # Center Matrix with mu as input
    Xa = copy(e_X)
    mu = vec(mean(Xa, dims=1))
    center!(Xa, mu, obsdim=1)
    @test abs(sum(mean(Xa, dims=1))) <= 10e-10

    Xa = copy(e_X)
    mu = vec(mean(Xa, dims=2))
    center!(Xa, mu, obsdim=2)
    @test abs(sum(mean(Xa, dims=2))) <= 10e-10

    Xa = copy(e_X)
    mu = vec(mean(Xa, dims=2))
    center!(Xa, mu, ObsDim.Last())
    @test abs(sum(mean(Xa, dims=2))) <= 10e-10

    # Center DataFrame
    D = copy(df)
    mu_check = [mean(D[!, colname]) for colname in names(D)[1:2]]
    mu = center!(D)
    @test length(mu) == 2
    @test abs(sum(mu .- mu_check)) <= 10e-10

    D = copy(df)
    mu_check = [mean(D[!, colname]) for colname in names(D)[1:2]]
    mu = center!(D, [:A, :B])
    @test abs(sum(mu .- mu_check)) <= 10e-10

    D = copy(df)
    mu_check = [mean(D[!, colname]) for colname in names(D)[1:2]]
    mu = center!(D, [:A, :B], mu_check)
    @test abs(sum([mean(D[!, colname]) for colname in names(D)[1:2]])) <= 10e-10

    # skip columns that contain missing values
    D = copy(df_na)
    mu = center!(D, [:A, :B])
    @test ismissing(D[1, :A])
    @test all(D[2:end, :A] .== df_na[2:end, :A])
    @test abs(mean(D[!, :B])) < 10e-10

    D = copy(df_na)
    mu_check = [mean(D[!, colname]) for colname in names(D)[1:2]]
    mu = center!(D, [:A, :B], mu_check)
    @test ismissing(D[1, :A])
    @test all(D[2:end, :A] .== df_na[2:end, :A])
    @test abs(mean(D[!, :B])) < 10e-10

    # Rescale Vector
    xa = copy(e_x)
    mu, sigma = rescale!(xa)
    @test mu ≈ mean(e_x)
    @test sigma ≈ std(e_x)
    @test abs(mean(xa)) <= 10e-10
    @test std(xa) ≈ 1

    xa = copy(e_x)
    mu, sigma = rescale!(xa, mu, sigma)
    @test abs(mean(xa)) <= 10e-10
    @test std(xa) ≈ 1

    xa = copy(e_x)
    mu, sigma = rescale!(xa, mu, sigma, obsdim=1)
    @test abs(mean(xa)) <= 10e-10
    @test std(xa) ≈ 1

    xa = copy(e_x)
    mu = copy(e_x) .- 1
    sigma = fill!(similar(e_x),1)
    mu, sigma = rescale!(xa, mu, sigma, obsdim=1)
    @test mean(xa) ≈ 1

    Xa = copy(e_X)
    rescale!(Xa)
    @test abs(sum(mean(Xa, dims=2))) <= 10e-10
    @test std(Xa, dims=2) ≈ [1, 1, 1, 1, 1]

    Xa = copy(e_X)
    rescale!(Xa, obsdim=2)
    @test abs(sum(mean(Xa, dims=2))) <= 10e-10
    @test std(Xa, dims=2) ≈ [1, 1, 1, 1, 1]

    Xa = copy(e_X)
    rescale!(Xa, obsdim=1)
    @test abs(sum(mean(Xa, dims=1))) <= 10e-10

    Xa = copy(e_X)
    mu = vec(mean(Xa, dims=1))
    sigma = vec(std(Xa, dims=1))
    rescale!(Xa, mu, sigma, obsdim=1)
    @test abs(sum(mean(Xa, dims=1))) <= 10e-10

    Xa = copy(e_X)
    mu = vec(mean(Xa, dims=2))
    sigma = vec(std(Xa, dims=2))
    rescale!(Xa, mu, sigma, obsdim=2)
    @test abs(sum(mean(Xa, dims=2))) <= 10e-10

    D = copy(df)
    mu, sigma = rescale!(D)
    @test abs(sum([mean(D[!, colname]) for colname in names(D)[1:2]])) <= 10e-10
    @test mean([std(D[!, colname]) for colname in names(D)[1:2]]) - 1 <= 10e-10

    D = copy(df)
    mu, sigma = rescale!(D, [:A, :B])
    @test abs(sum([mean(D[!, colname]) for colname in names(D)[1:2]])) <= 10e-10
    @test mean([std(D[!, colname]) for colname in names(D)[1:2]]) - 1 <= 10e-10

    D = copy(df)
    mu_check = [mean(D[!, colname]) for colname in names(D)[1:2]]
    sigma_check = [std(D[!, colname]) for colname in names(D)[1:2]]
    mu, sigma = rescale!(D, [:A, :B], mu_check, sigma_check)
    @test abs(sum([mean(D[!, colname]) for colname in names(D)[1:2]])) <= 10e-10
    @test mean([std(D[!, colname]) for colname in names(D)[1:2]]) - 1 <= 10e-10

    # skip columns that contain missing values
    D = copy(df_na)
    mu, sigma = rescale!(D, [:A, :B])
    @test ismissing(D[1, :A])
    @test all(D[2:end, :A] .== df_na[2:end, :A])
    @test abs(mean(D[!, :B])) < 10e-10
    @test abs(std(D[!, :B])) - 1 < 10e-10

    D = copy(df_na)
    mu_check = [mean(D[!, colname]) for colname in names(D)[1:2]]
    if VERSION >= v"0.7.0-DEV.2035"
        sigma_check = [std(D[!, colname]) for colname in names(D)[1:2]]
        mu, sigma = rescale!(D, [:A, :B], mu_check, sigma_check)
    end
    #= @test ismissing(D[1, :A]) =#
    #= @test all(D[2:end, :A] .== df_na[2:end, :A]) =#
    #= @test abs(mean(D[:B])) < 10e-10 =#
    #= @test (abs(std(D[:B])) - 1) < 10e-10 =#
end

e_x = collect(-5:.1:5)
e_X = [e_x e_x.^2 e_x.^3]'

@testset "Test FeatureNormalizer model" begin
    cs = fit(FeatureNormalizer, e_X)
    @test vec(mean(e_X, dims=2)) ≈ cs.offset
    @test vec(std(e_X, dims=2)) ≈ cs.scale

    Xa = predict(cs, e_X)
    @test Xa != e_X
    @test abs(sum(mean(Xa, dims=2))) <= 10e-10
    @test std(Xa, dims=2) ≈ [1, 1, 1]
end


text:
module ArrayFire

using Cxx
using Base.Meta

import Base: rand, show, randn, ones, diag, eltype, size, elsize,
    sizeof, length, showarray, convert, ndims, lu, qr, svd, setindex!
import Cxx: CppEnum
export AFArray, chol!, constant, aftype, @gpu, @icxx_str, __current_compiler__

# If you have a crash, enable this
const AF_DEBUG = true

function init_library()
    if !haskey(ENV, "AFMODE")
        info("Setting Unified Backend by default")
        Libdl.dlopen("libaf", Libdl.RTLD_GLOBAL)
    else
        if ENV["AFMODE"] == "CPU"
            info("Setting CPU Backend")
            Libdl.dlopen("libafcpu", Libdl.RTLD_GLOBAL)
        elseif ENV["AFMODE"] == "OPENCL"
            info("Setting OPENCL Backend")
            Libdl.dlopen("libafopencl", Libdl.RTLD_GLOBAL)
        elseif ENV["AFMODE"] == "CUDA"
            info("Setting GPU Backend")
            Libdl.dlopen("libafcuda", Libdl.RTLD_GLOBAL)
        else
            info("Invalid value for environment variable AFMODE. Setting to default.")
            Libdl.dlopen("libaf", Libdl.RTLD_GLOBAL)
        end
    end
    if !haskey(ENV, "AFPATH")
       addHeaderDir("/usr/local/include/"; kind = C_System)
    else
        addHeaderDir(ENV["APATH"]; kind = C_System)
    end
end
init_library()

cxx"""
#include <arrayfire.h>
"""

setDevice(i) = icxx"af::setDevice($i);"

const f32 = CppEnum{:af_dtype,UInt32}(0) # 32-bit floating point values
const c32 = CppEnum{:af_dtype,UInt32}(1) # 32-bit complex floating point values
const f64 = CppEnum{:af_dtype,UInt32}(2) # 64-bit complex floating point values
const c64 = CppEnum{:af_dtype,UInt32}(3) # 64-bit complex floating point values
const b8  = CppEnum{:af_dtype,UInt32}(4) #  8-bit boolean values
const s32 = CppEnum{:af_dtype,UInt32}(5) # 32-bit signed integral values
const u32 = CppEnum{:af_dtype,UInt32}(6) # 32-bit unsigned integral values
const u8  = CppEnum{:af_dtype,UInt32}(7) #  8-bit unsigned integral values
const s64 = CppEnum{:af_dtype,UInt32}(8) # 64-bit signed integral values
const u64 = CppEnum{:af_dtype,UInt32}(9) # 64-bit unsigned integral values

aftype(::Type{Float32})          = f32
aftype(::Type{Complex{Float32}}) = c32
aftype(::Type{Float64})          = f64
aftype(::Type{Complex{Float64}}) = c64
aftype(::Type{Bool})             = b8
aftype(::Type{Int32})            = s32
aftype(::Type{UInt32})           = u32
aftype(::Type{UInt8})            = u8
aftype(::Type{Int64})            = s64
aftype(::Type{UInt64})           = u64

function jltype(af_dtype)
    if af_dtype == f32
        return Float32
    elseif af_dtype == c32
        return Complex{Float32}
    elseif af_dtype == f64
        return Float64
    elseif af_dtype == c64
        return Complex{Float64}
    elseif af_dtype == b8
        return Bool
    elseif af_dtype == s32
        return Int32
    elseif af_dtype == u32
        return UInt32
    elseif af_dtype == u8
        return UInt8
    elseif af_dtype == s64
        return Int64
    elseif af_dtype == u64
        return UInt64
    end
end

function af_promote{T,S}(::Type{T},::Type{S})
    if T == S
        return T
    elseif T == Complex{Float64} || S == Complex{Float64}
        return Complex{Float64}
    elseif T == Complex{Float32} || S == Complex{Float32}
        (T == Float64 || S == Float64) && return Complex{Float64}
        return Complex{Float32}
    elseif T == Float64 || S == Float64
        return Float64
    elseif T == Float32 || S == Float32
        return Float32
    elseif T == UInt64 || S == UInt64
        return UInt64
    elseif T == Int64 || S == Int64
        return Int64
    elseif T == UInt32 || S == UInt32
        return UInt32
    elseif T == Int32 || S == Int32
        return Int32
    elseif T == UInt8 || S == UInt8
        return UInt8
    elseif T == Bool || S == Bool
        return Bool
    else
        return Float32
    end
end

abstract AFAbstractArray{T,N} <: AbstractArray{T,4}

immutable AFArray{T,N} <: AFAbstractArray{T,N}
    array::vcpp"af::array"
    function AFArray(array::vcpp"af::array")
        ret = new(array)
        if AF_DEBUG && backend_eltype(ret) != T
            error("Tried to create AFArray{$T} with array of eltype $(backend_eltype(ret))")
        end
        ret
    end
end
immutable AFSubArray{T,N} <: AFAbstractArray{T,N}
    array::vcpp"af::array::array_proxy"
    function AFSubArray(array::vcpp"af::array::array_proxy")
        ret = new(array)
        if AF_DEBUG && backend_eltype(ret) != T
            error("Tried to create AFArray{$T} with array of eltype $(backend_eltype(ret))")
        end
        ret
    end
end

AFArray() = icxx"af::array();"

call{T,N}(::Type{AFAbstractArray{T,N}}, proxy::vcpp"af::array::array_proxy") =
    AFSubArray{T,N}(proxy)
call{T,N}(::Type{AFAbstractArray{T,N}}, arr::vcpp"af::array") =
    AFArray{T,N}(arr)
ndims(arr::vcpp"af::array") = icxx"$arr.numdims();"
ndims(arr::vcpp"af::array::array_proxy") = icxx"$arr.numdims();"
(::Type{AFArray{T}}){T}(arr::vcpp"af::array") = AFArray{T,Int(ndims(arr))}(arr)
(::Type{AFSubArray{T}}){T}(arr::vcpp"af::array::array_proxy") = AFSubArray{T,Int(ndims(arr))}(arr)

Cxx.cppconvert{T}(x::AFAbstractArray{T}) = x.array

convert{T,N}(::Type{AFArray{T,N}}, arr::AFSubArray{T,N}) =
    AFArray{T,N}(icxx"(af::array)$arr;")
(::Type{AFArray})(arr::AFSubArray) = AFArray{backend_eltype(arr), ndims(arr)}(arr)
convert(::Type{AFArray}, arr::AFSubArray) = AFArray(arr)

eltype{T}(x::AFAbstractArray{T}) = T
backend_eltype(x) = jltype(icxx"$x.type();")
sizeof{T}(a::AFAbstractArray{T}) = elsize(a) * length(a)

# GPU to Host
function convert{T,N}(::Type{Array{T,N}},x::AFAbstractArray{T,N})
    ret = Array(UInt8,sizeof(x))
    icxx"$x.host($(pointer(ret)));"
    ret = reinterpret(T, ret)
    ret = reshape(ret, size(x)...)
    ret
end
convert{T,N}(::Type{Array},x::AFAbstractArray{T,N}) = convert(Array{T,N},x)

# Host to GPU
function convert{T,N}(::Type{AFArray{T,N}}, x::Array{T,N})
    AFArray{T,N}(icxx"return af::array{$(dims_to_dim4(size(x))),$(pointer(x)),afHost};")
end
convert{T,N}(::Type{AFArray}, x::Array{T,N}) = convert(AFArray{T,N}, x)

# Show this array, by getting a data pointer to it and using julia's printing
# mechanism. ArrayFire copies internally anyway, so there is nothing to be gained
# by using it's printing
function showarray{T,N}(io::IO, X::AFAbstractArray{T,N};
                   header::Bool=true, kwargs...)
    header && print(io, summary(X))
    !isempty(X) && println(io,":")
    showarray(io, convert(Array{T,N},X); header = false, kwargs...)
end

Base.show(io::IO, arr::vcpp"af::array") = print(io, "ArrayFire Array")

af_print(X::AFAbstractArray) = icxx"""af::print("",$X);"""

# Translate dimensions. Note that we need to translate between 0 and 1 based
# indexing
ndims(dim4) = icxx"$dim4.ndims();"
ndims(A::AFAbstractArray) = ndims(icxx"$A.dims();")
function dims_to_dim4(dims)
    if length(dims) == 1
        icxx"af::dim4{$(dims[1])};"
    elseif length(dims) == 2
        icxx"af::dim4{$(dims[1]),$(dims[2])};"
    elseif length(dims) == 3
        icxx"af::dim4{$(dims[1]),$(dims[2]),$(dims[3])};"
    elseif length(dims) == 4
        icxx"af::dim4{$(dims[1]),$(dims[2]),$(dims[3]),$(dims[4])};"
    else
        error("Too many dimensions")
    end
end
function dim4_to_dims(dim4)
    d = ndims(dim4)
    t = (Int(icxx"(int)$dim4[0];"),)
    d == 1 && return t
    t = tuple(t...,Int(icxx"(int)$dim4[1];"))
    d == 2 && return t
    t = tuple(t...,Int(icxx"(int)$dim4[2];"))
    d == 3 && return t
    t = tuple(t...,Int(icxx"(int)$dim4[3];"))
    d == 4 && return t
    error("Bad dim4 object")
end

size(x::AFAbstractArray) = dim4_to_dims(icxx"($x).dims();")

#Specialize vectors and matrices
typealias AFMatrix{T} AFAbstractArray{T,2}
typealias AFVector{T} AFAbstractArray{T,1}

#Functions to generate array
import Base:fill, eye, diag
export iota
constant{T<:Real}(val::T, dims::Integer...) = AFArray{T}(af_constant(val, dims_to_dim4(dims), aftype(typeof(val))))
rand{T}(::Type{AFArray{T}}, dims::Integer...) = AFArray{T}(af_randu(dims_to_dim4(dims), aftype(T)))
randn{T}(::Type{AFArray{T}}, dims::Integer...) = AFArray{T}(af_randn(dims_to_dim4(dims), aftype(T)))
eye{T}(::Type{AFArray{T}}, dims::Integer...) = AFArray{T}(af_identity(dims_to_dim4(dims), aftype(T)))
diag{T}(x::AFArray{T}, k = 0) = AFArray{T}(af_diag(x, k))
getSeed() = af_getSeed()
setSeed(a::Integer) = af_setSeed(a)
range{T}(::Type{AFArray{T}}, dims::Integer...) = AFArray{T}(af_range(dims_to_dim4(dims), length(dims)-1, aftype(T)))
iota{T}(::Type{AFArray{T}}, dims::Integer...) = AFArray{T}(af_iota(dims_to_dim4(dims), dims_to_dim4(1), aftype(T)))
moddims{T}(a::AFAbstractArray{T}, dims::Integer...) = AFArray{T}(af_moddims(a, dims_to_dim4(dims)))

Base.rand{A<:AFArray}(::Type{A},sz::Tuple)=rand(A, sz...)

#TODO : make `tile` compatible with `repeat` in base

export tile

tile{T}(a::AFAbstractArray{T}, inds::Integer...) = AFArray{T}(af_tile(a, dims_to_dim4(inds)))

import Base: getindex

# Note that we need to translate between 0 and 1 based indexing
to_af_idx(x::Real) = convert(Cint, x-1)
to_af_idx(x::Range) = icxx"af::seq($(first(x))-1,$(last(x))-1,$(step(x)));"
to_af_idx(x::Colon) = icxx"af::span;"
to_af_idx(x::rcpp"af::seq") = x
const tis = to_af_idx

IS = Union{Real,Range,Colon,rcpp"af::seq"}

#ArrayFire needs to index with Int32s or Float32s. Also, array_proxy wrap for indexing with arrays crashes currently.
function _getindex(x::AFAbstractArray,y::AFAbstractArray)
    idx = y - 1
    out = AFArray()
    icxx"$out = $x($idx);"
    AFArray{backend_eltype(out)}(out)
end
# Avoid an ambiguity
_getindex{T}(x::AFAbstractArray{T},s0::Real)             = icxx"$x($(tis(s0)));"

_getindex{T}(x::AFAbstractArray{T},s0::IS)               = icxx"$x($(tis(s0)));"

_getindex{T}(x::AFAbstractArray{T},s0::IS,s1::IS)        = icxx"$x($(tis(s0)),
                                                                   $(tis(s1)));"

_getindex{T}(x::AFAbstractArray{T},s0::IS,s1::IS,s2::IS) = icxx"$x($(tis(s0)),
                                                                   $(tis(s1)),
                                                                   $(tis(s2)));"
_getindex{T}(x::AFAbstractArray{T},
                    s0::IS,s1::IS,s2::IS, s3::IS)        = icxx"$x($(tis(s0)),
                                                                   $(tis(s1)),
                                                                   $(tis(s2)),
                                                                   $(tis(s3)));"

#_getindex{T}(x::AFAbstractArray{T}, i) = icxx"$x($i)"

function getindex{T}(x::AFAbstractArray{T}, idxs...)
    proxy = _getindex(x,idxs...)
    if isa(proxy, AFArray)
        return proxy
    end
    num_elems = icxx"$proxy.elements();"
    if num_elems == 1
        return Array(AFSubArray{T}(proxy))[1]
    else
        return AFSubArray{T}(proxy)
    end
end

#TODO : return 0-element array when b is all falses
function getindex(x::AFAbstractArray, b::AFAbstractArray{Bool})
    out = AFArray()
    if any(b)
        icxx"$out = $x($b);"
        AFArray{backend_eltype(out)}(out)
    else
        return 0
    end
end

function setindex!{T}(x::AFAbstractArray{T}, val, idxs...)
    proxy = _getindex(x,idxs...)
    icxx"$proxy = $val;"
end

# Avoid crash when referencing with boolean arrays with all falses
function setindex!(x::AFAbstractArray, val, b::AFAbstractArray{Bool})
    if any(b)
        icxx"$x($b) = $val;"
    else
        return val
    end
end

# Exception handling
import Base: showerror
@exception function showerror(io::IO, e::rcpp"af::exception")
    println(io, bytestring(icxx"$e.what();"))
end

immutable AFFeatures
    feat::vcpp"af::features"
    function AFFeatures(f::vcpp"af::features")
        new(f)
    end
end
AFFeatures() = icxx"af::features();"

Cxx.cppconvert(f::AFFeatures) = f.feat

immutable AFWindow
    win::vcpp"af::Window"
    function AFWindow(w::vcpp"af::Window")
        new(w)
    end
end

Cxx.cppconvert(w::AFWindow) = w.win

#import other files
include("AFWrap.jl")
include("math.jl")
include("image.jl")
include("stats.jl")
include("vector.jl")
include("linalg.jl")
include("signal.jl")
include("graphics.jl")
include("unified.jl")

#Info
export af_info, af_isDoubleAvailable, af_sync, getDevice
getDevice() = af_getDevice()

#Storing and Loading Arrays
save(a::AFAbstractArray, path::AbstractString, key::AbstractString) = af_saveArray(key, a, path)
function load(path::AbstractString, key::AbstractString)
    out = af_readArray(path, key)
    AFArray{backend_eltype(out)}(out)
end

macro gpu(ex...)
    ex = ex[1]
    top = ex.args[1]
    bottom = ex.args[2]
    iterator = top.args[1]
    start = :($(top.args[2]).start)
    step = :(step($(top.args[2])))
    stop = :($(top.args[2]).stop)
    x = bottom.args[2]
    body = Expr(:->,:i,Expr(:block, x, :nothing))
    esc(Expr(:let, quote icxx"""
                    gfor (af::seq i, $start - 1, $stop - 1, $step)
                    {
                        $body(i);
                    }
                    """
                    end , :(body = $body), :(stop = $stop), :(start = $start), :(step = $step)))
end

end # module


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# =================================================================================================
"""
    mutable struct TessInst
        ptr::Ptr{Cvoid}
    end

A wrapper for the Api object in the Tesseract library.

__Values:__

| Name | Description
| :--- | :----------
| ptr  | The pointer to the Api object allocated by the C library.

__Details:__

Most method calls cannot use this object until `tess_init()` called on it for initialization.

When the garbage collector collects this object the associated pointer object will be freed in the
library.  The object can also be manually freed by calling `tess_delete!()` on it.

See also: [`TessInst(languages::AbstractString, dataPath::AbstractString)`](@ref).
"""
mutable struct TessInst
    ptr::Ptr{Cvoid}
end

# =================================================================================================
"""
    TessInst(
        languages::AbstractString = "eng",
        dataPath::AbstractString = TESS_DATA
    )

Construct an initialize a TessInst object.

__Arguments:__

| T | Name      | Default    | Description
|:--| :-------- | :--------- | :----------
| O | languages | `eng`      | The language(s) to load.
| O | dataPath  | `tessdata` | The directory to look for the language files in.

__Details:__

To change the langauges identified by this instance you can call `tess_init()`.  Multiple
langagues can be specified by seperating them with a plus.  So if you want english and spanish you
could specify "eng+spa".  The language codes are (usually) the
[ISO 639-3](https://en.wikipedia.org/wiki/ISO_639-3) code.

__Example:__

```jldoctest
julia> using Tesseract

julia> download_languages("eng+fra")
true

julia> instance = TessInst("eng+fra")
Allocated Tesseract instance.
```

See also: [`tess_init`](@ref).
"""
function TessInst(
            languages::AbstractString = "eng",
            dataPath::AbstractString = TESS_DATA
        )
    local ptr    = ccall((:TessBaseAPICreate, TESSERACT), Ptr{Cvoid}, ())
    local retval = TessInst(ptr)

    finalizer(retval) do obj
        tess_delete!(obj)
    end

    tess_init(retval, languages, dataPath)

    retval
end

# =================================================================================================
"""
    show(
        io::IO,
        inst::TessInst
    )::Nothing

Display summary information about the tesseract instance.

__Arguments:__

| T | Name  | Default | Description
|---| :---- | :------ | :----------
| R | io    |         | The stream to write the information to.
| R | inst  |         | The tesseract instance to display info about.
"""
function Base.show(
            io::IO,
            inst::TessInst
        )::Nothing
    if is_valid(inst)
        print(io, "Allocated Tesseract instance.")
    else
        print(io, "Freed Tesseract instance.")
    end
    nothing
end

# =================================================================================================
"""
    is_valid(
        inst::TessInst
    )::Bool

Check if the instance has been freed or if it's still valid.

__Arguments:__

| T | Name | Default | Description
|---| :--- | :------ | :----------
| R | inst |         | The instance to check.
"""
function is_valid(
            inst::TessInst
        )::Bool
    return inst.ptr != C_NULL
end

# =================================================================================================
"""
    unsafe_convert(
        ::Type{Ptr{Cvoid}},
        inst::TessInst
    )::Ptr{Cvoid}

Convert the instance into a the handle pointer used by the Tesseract library.

__Arguments:__

| T | Name               | Default | Description
|---| :----------------- | :------ | :----------
| R | ::Type{Ptr{Cvoid}} |         | The type to convert into.
| R | inst               |         | The instance to return the Teseract handle for.
"""
Base.unsafe_convert(::Type{Ptr{Cvoid}}, inst::TessInst) = inst.ptr

# =================================================================================================
"""
    tess_delete!(
        inst::TessInst
    )::Nothing

Destroy the Tesseract object and release any associated memory.

__Arguments:__

| T | Name | Default | Description
|---| :--- | :------ | :----------
| R | inst |         | The instance to free.

__Details:__

This method is called automatically by the garbage collector but can be called manually to release
the object early.  This method can be called multiple times without any negative effects.

Once `tess_delete!()` is called on an object passing the object to any other library call will
result in an error.

Note: This method is not thread safe.
"""
function tess_delete!(
            inst::TessInst
        )::Nothing

    if is_valid(inst) == true
        ccall(
            (:TessBaseAPIDelete, TESSERACT),
            Cvoid,
            (Ptr{Cvoid},),
            inst
        )
        inst.ptr = C_NULL
    end
    nothing
end

# =================================================================================================
"""
    tess_recognize(
        inst::TessInst
    )::Bool

Perform the OCR extraction.  This will be called automatically if you call one of the retrieval
functions so you don't need to call it directly.  Returns `false` if there was an error.

__Arguments:__

| T | Name | Default | Description
|---| :--- | :------ | :----------
| R | inst |         | The instance to perform the recognition with.

__Example:__

```jldoctest
julia> using Tesseract

julia> download_languages("eng+fra")
true

julia> instance = TessInst("eng+fra")
Allocated Tesseract instance.

julia> pix = sample_pix()
Image (500, 600) at 32ppi

julia> tess_image(instance, pix)

julia> tess_resolution(instance, 72)

julia> tess_recognize(instance)
true
```

See also: [`tess_text`](@ref), [`tess_hocr`](@ref), [`tess_alto`](@ref), [`tess_tsv`](@ref),
[`tess_parsed_tsv`](@ref).
"""
function tess_recognize(
            inst::TessInst
        )::Bool
    local result = -1

    if is_valid(inst) == true
        result = @threadcall(
            (:TessBaseAPIRecognize, TESSERACT),
            Cint,
            (Ptr{Cvoid}, Ptr{Cvoid}),
            inst,
            C_NULL
        )
    else
        @error "Instance has been freed."
    end

    return result == 0
end