text:
@testset "REST API" begin
    @testset "Direct endpoint wrapper" begin
        # Direct endpoint wrappers should return a Future.
        f = get_channel_message(c, 123, 456)
        @test f isa Future
        # Since we don't have a valid token, we shouldn't get anything.
        @test fetchval(f) === nothing
        r = fetch(f)
        # But the type should still be sound.
        @test r isa Response{Message}
        @test !r.ok
        @test r.val === nothing
        @test r.http_response !== nothing
        @test r.exception === nothing
    end

    @testset "CRUD" begin
        # This API should behave just like the direct endpoint wrappers.
        f = create(c, Guild; name="foo")
        @test f isa Future
        r = fetch(f)
        @test r isa Response{Guild}
        @test r.val === nothing
        @test !r.ok
        @test r.http_response !== nothing
        @test r.exception === nothing
    end

    @testset "Simultaneous requests" begin
        # We should be able to make a bunch of requests without deadlocking.
        fs = map(i -> retrieve(c, Guild, i), 1:10)
        @test all(f -> f isa Future, fs)
        rs = fetch.(fs)
        @test all(r -> r isa Response{Guild}, rs)
    end

    @testset "Selective caching" begin
        f(c, e) = nothing

        # The regexes should be strict.
        @test all(
            p -> startswith(p, '^') && endswith(p, '$'),
            map(t -> t[1].pattern, vcat(values(EVENTS_FIRED)...)),
        )

        # We don't need to worry about this fast path.
        @eval Discord Base.isopen(::Client) = true

        # GETs should always be cached, because they never fire events.
        @test should_put(c, :GET, "foo")

        # If we have the default handler in place, we shouldn't put.
        add_handler!(c, MessageCreate, f; tag=DEFAULT_HANDLER_TAG)
        @test !should_put(c, :POST, "/channels/1/messages")

        # But if we delete it, then we should.
        delete_handler!(c, MessageCreate, DEFAULT_HANDLER_TAG)
        @test should_put(c, :POST, "/channels/1/messages")

        # For endpoints which fire multiple events, put unless they all have a default.
        add_handler!(c, GuildMemberRemove, f; tag=DEFAULT_HANDLER_TAG)
        add_handler!(c, PresenceUpdate, f; tag=DEFAULT_HANDLER_TAG)
        @test !should_put(c, :DELETE, "/users/@me/guilds/1")
        delete_handler!(c, PresenceUpdate, DEFAULT_HANDLER_TAG)
        @test should_put(c, :DELETE, "/users/@me/guilds/1")
    end
end


text:
# Contains the physical specification of the model

# A physical specification of a model.
# Contains the geometry information, but no discretization parameters.
# The exact model used is defined by the list of terms.
struct Model{T <: Real}
    # Lattice and reciprocal lattice vectors in columns
    lattice::Mat3{T}
    recip_lattice::Mat3{T}
    unit_cell_volume::T
    recip_cell_volume::T
    dim::Int # Dimension of the system; 3 unless `lattice` has zero columns

    # Electrons, occupation and smearing function
    n_electrons::Int # usually consistent with `atoms` field, but doesn't have to

    # spin_polarization values:
    #     :none       No spin polarization, αα and ββ density identical,
    #                 αβ and βα blocks zero
    #     :collinear  Spin is polarized, but everywhere in the same direction.
    #                 αα ̸= ββ, αβ = βα = 0
    #     :full       Generic magnetization, non-uniform direction.
    #                 αβ, βα, αα, ββ all nonzero, different
    #     :spinless   No spin at all ("spinless fermions", "mathematicians' electrons").
    #                 The difference with :none is that the occupations are 1 instead of 2
    spin_polarization::Symbol

    # If temperature==0, no fractional occupations are used.
    # If temperature is nonzero, the occupations are
    # `fn = max_occ*smearing((εn-εF) / temperature)`
    temperature::T
    smearing::Smearing.SmearingFunction # see Smearing.jl for choices

    # Vector of pairs Element => vector of vec3 (positions, fractional coordinates)
    # Possibly empty. `atoms` just contain the information on the atoms, it's up to
    # the `terms` to make use of it (or not)
    atoms::Vector{Pair{Any, Vector{Vec3{T}}}}

    # each element t must implement t(basis), which instantiates a
    # term in a given basis and gives back a term (<: Term)
    # see terms.jl for some default terms
    term_types::Vector

    # list of symmetries of the model
    symops::Vector{SymOp}
end

"""
    Model(lattice; n_electrons, atoms, terms, temperature,
                   smearing, spin_polarization, symmetry)

Creates the physical specification of a model (without any
discretization information).

`n_electrons` is taken from `atoms` if not specified

`spin_polarization` is :none by default (paired electrons)

`smearing` is Fermi-Dirac if `temperature` is non-zero, none otherwise

The `symmetry` kwarg can be:
- auto: determine from terms if they respect the symmetry of the lattice/atoms.
- off: no symmetries at all
- force: force all the symmetries of the lattice/atoms.
Careful that in this last case, wrong results can occur if the
external potential breaks symmetries (this is not checked).
"""
function Model(lattice::AbstractMatrix{T};
               n_electrons=nothing,
               atoms=[],
               terms=[],
               temperature=T(0.0),
               smearing=nothing,
               spin_polarization=:none,
               symmetry=:auto
               ) where {T <: Real}

    lattice = Mat3{T}(lattice)

    if n_electrons === nothing
        # get it from the atom list
        isempty(atoms) && error("Either n_electrons or a non-empty atoms list should be provided")
        n_electrons = sum(length(pos) * n_elec_valence(spec) for (spec, pos) in atoms)
    else
        @assert n_electrons isa Int
    end

    # Special handling of 1D and 2D systems, and sanity checks
    d = 3-count(c -> norm(c) == 0, eachcol(lattice))
    d > 0 || error("Check your lattice; we do not do 0D systems")
    for i = d+1:3
        norm(lattice[:, i]) == norm(lattice[i, :]) == 0 || error(
            "For 1D and 2D systems, the non-empty dimensions must come first")
    end
    cond(lattice[1:d, 1:d]) > 1e-5 || @warn "Your lattice is badly conditioned, the computation is likely to fail"

    # Compute reciprocal lattice and volumes.
    # recall that the reciprocal lattice is the set of G vectors such
    # that G.R ∈ 2π ℤ for all R in the lattice
    recip_lattice = zeros(T, 3, 3)
    recip_lattice[1:d, 1:d] = 2T(π)*inv(lattice[1:d, 1:d]')
    recip_lattice = Mat3{T}(recip_lattice)
    # in the 1D or 2D case, the volume is the length/surface
    unit_cell_volume = abs(det(lattice[1:d, 1:d]))
    recip_cell_volume = abs(det(recip_lattice[1:d, 1:d]))

    @assert spin_polarization in (:none, :collinear, :full, :spinless)

    if smearing === nothing
        @assert temperature >= 0
        # Default to Fermi-Dirac smearing when finite temperature
        smearing = temperature > 0.0 ? Smearing.FermiDirac() : Smearing.None()
    end

    if !allunique(string.(nameof.(typeof.(terms))))
        error("Having several terms of the same name is not supported.")
    end

    @assert symmetry in (:auto, :force, :off)
    if symmetry == :auto
        # check dimension : in 1D or 2D we might have symmetries but it's not
        # detected by spglib so we ignore them
        compute_symmetry = !(any(breaks_symmetries, terms)) && (d == 3)
    else
        compute_symmetry = (symmetry == :force)
    end

    if compute_symmetry
        symops = symmetry_operations(lattice, atoms)
    else
        symops = [identity_symop()]
    end

    Model{T}(lattice, recip_lattice, unit_cell_volume, recip_cell_volume, d, n_electrons,
             spin_polarization, T(temperature), smearing, atoms, terms, symops)
end

"""
Convenience constructor, which builds a standard atomic (kinetic + atomic potential) model.
Use `extra_terms` to add additional terms.
"""
function model_atomic(lattice::AbstractMatrix, atoms::Vector; extra_terms=[], kwargs...)
    @assert !(:terms in keys(kwargs))
    @assert !(:atoms in keys(kwargs))
    terms = [Kinetic(),
             AtomicLocal(),
             AtomicNonlocal(),
             Ewald(),
             PspCorrection(),
             extra_terms...]
    if :temperature in keys(kwargs) && kwargs[:temperature] != 0
        terms = [terms..., Entropy()]
    end
    Model(lattice; atoms=atoms, terms=terms, kwargs...)
end

"""
Build a DFT model from the specified atoms, with the specified functionals.
"""
function model_DFT(lattice::AbstractMatrix, atoms::Vector, functionals; extra_terms=[], kwargs...)
    model_atomic(lattice, atoms; extra_terms=[Hartree(), Xc(functionals), extra_terms...], kwargs...)
end

"""
Build an LDA model (Teter93 parametrization) from the specified atoms.
"""
function model_LDA(lattice::AbstractMatrix, atoms::Vector; kwargs...)
    model_DFT(lattice, atoms, :lda_xc_teter93; kwargs...)
end

"""
Build an PBE-GGA model from the specified atoms.
"""
function model_PBE(lattice::AbstractMatrix, atoms::Vector; kwargs...)
    model_DFT(lattice, atoms, [:gga_x_pbe, :gga_c_pbe]; kwargs...)
end

"""
Maximal occupation of a state (2 for non-spin-polarized electrons, 1 otherwise).
"""
function filled_occupation(model)
    @assert model.spin_polarization in (:none, :spinless)
    if model.spin_polarization == :none
        @assert model.n_electrons % 2 == 0
        filled_occ = 2
    else
        filled_occ = 1
    end
    filled_occ
end

ext: jl


text:
#using Iterators

#############################################
###### Structs ##############################
#############################################
struct Reaction ## Reaction Object
	index::Int64 # List index of reaction
	reactants::Array{Int64,1} # List of reactant indices
	react_coef::Array{Int64,1} # List of reactant coefficients
	products::Array{Int64,1} # List of product indices
	prod_coef::Array{Int64,1} # List of product coefficients
	catalysts::Array{Int64,1} # Catalyst indicies
	cat_coef::Array{Float64,1} # Catalysts coefficients
	propensity::String # Propensity identifier
	rate_constant::Float64 # Rate constant
end

struct CRS # Chemical reaction system object
	molecule_list::Array{String,1} # List of molcule strings
	molecule_dict::Dict{String, Int64} # Dict maps molecule strings to index
	reaction_list::Array{Reaction,1} # List of reaction objects
	parameters::Dict{Any,Any} # Every other parameter
end

#############################################
###### CRS functions ########################
#############################################

function generate_all_binary_seqs(max_L)
    ### Generate all possible sequences
    all_seqs = []
    for j in 1:max_L
        comps = Array{String}(undef,j + 1)
        for l in 0:j
            comps[l+1] = "B"^l * "A"^(j-l)
        end

        for c in comps
            p = unique(collect(permutations(c)))
            all_seqs = vcat(p,all_seqs)
        end
    end
    seqs = [join(s) for s in all_seqs]
    sort!(seqs, by = length, rev= false)
    return seqs
end

### Count occurance of t in s Function ##
function countall(s, t)
    n = length(collect(eachmatch(Regex(t), s)))
    return n
end

function all_splits(seq)
    l = length(seq)
    left = Array{String}(undef,(l-1))
    right =Array{String}(undef,(l-1))

    for i in 1:(l-1)
        left[i] = seq[1:i]
        right[i] = seq[(i+1):end]
    end
    return left,right
end

function calculate_reduced_mass_binary(m1_string, m2_string, mA,mB)
    m1_mass = countall(m1_string,"A")*mA + countall(m1_string,"B")*mB
    m2_mass = countall(m2_string,"A")*mA + countall(m2_string,"B")*mB
    reduce = reduced_mass(m1_mass,m2_mass)
    return reduce
end

function reduced_mass(m1,m2)
    return (m1*m2)/(m1+m2)
end

function bimolecular_coef(k, T, R, V,m1, m2)
    m = reduced_mass(m1,m2)
    return (k/V)*sqrt(R*T/m)
end

function copy_CRS(original_CRS, volume_multiplier)
    reaction_list = similar(original_CRS.reaction_list)
    for rID in 1:length(reaction_list)
        rxn = original_CRS.reaction_list[rID]
        if sum(rxn.react_coef) ==1
            reaction_list[rID] = deepcopy(original_CRS.reaction_list[rID])
        elseif sum(rxn.react_coef) == 2
            new_rxn = Reaction(rxn.index, rxn.reactants, rxn.react_coef, rxn.products, rxn.prod_coef, rxn.catalysts, rxn.cat_coef, rxn.propensity, (1.0/volume_multiplier)*rxn.rate_constant)
            reaction_list[rID] = new_rxn
        end

    end
    new_parameters = original_CRS.parameters
    new_parameters[:volume] *= volume_multiplier

    new_CRS = CRS(original_CRS.molecule_list, original_CRS.molecule_dict, reaction_list, new_parameters )
    return new_CRS
end

function generate_tuples(max, sum)
    pairs = collect(Iterators.product(1:max, 1:max))
    my_pairs = []
    for pa in pairs
        if sum(pa) <= sum
            if !(pa in my_pairs) && !(reverse(pa) in my_pairs)
                push!(my_pairs,pa)
            end
        end
    end
    return my_pairs
end

max_stars_repo_name: croningp/clusterautocatalyticsets

max_issues_repo_name: croningp/clusterautocatalyticsets

max_forks_repo_name: croningp/clusterautocatalyticsets


text:
using BinaryBuilder, Pkg

name = "QuantumEspresso"
version = v"7.0.0"

sources = [
    ArchiveSource("https://gitlab.com/QEF/q-e/-/archive/qe-7.0/q-e-qe-7.0.tar.gz",
                  "85beceb1aaa1678a49e774c085866d4612d9d64108e0ac49b23152c8622880ee"),
    DirectorySource("bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd q-e-qe-*
atomic_patch -p1 ../patches/0000-pass-host-to-configure.patch

export BLAS_LIBS="-L${libdir} -lopenblas"
export LAPACK_LIBS="-L${libdir} -lopenblas"
export FFTW_INCLUDE=${includedir}
export FFT_LIBS="-L${libdir} -lfftw3"
export FC=mpif90
export CC=mpicc
export LD=

flags=(--enable-parallel=yes)
if [ "${nbits}" == 64 ]; then
    # Enable Libxc support only on 64-bit platforms
    atomic_patch -p1 ../patches/0001-libxc-prefix.patch
    flags+=(--with-libxc=yes --with-libxc-prefix=${prefix})
fi

if [[    "${target}" == powerpc64le-linux-* \
      || "${bb_full_target}" == armv6l-linux-* \
      || "${target}" == aarch64-apple-darwin* ]]; then
    # No scalapack binary available on these platforms
    flags+=(--with-scalapack=no)
else
    export SCALAPACK_LIBS="-L${libdir} -lscalapack"
    flags+=(--with-scalapack=yes)
fi

./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} ${flags[@]}
make all "${make_args[@]}" -j $nproc
make install
# Manually make all binary executables...executable.  Sigh
chmod +x "${bindir}"/*
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_gfortran_versions(supported_platforms())
filter!(!Sys.iswindows, platforms)
# On aarch64-apple-darwin we get
#    f951: internal compiler error: in doloop_contained_procedure_code, at fortran/frontend-passes.c:2464
filter!(p -> !(Sys.isapple(p) && arch(p) == "aarch64"), platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("pw.x", :pwscf),
    ExecutableProduct("bands.x", :bands),
    ExecutableProduct("plotband.x", :plotband),
    ExecutableProduct("plotrho.x", :plotrho),
    ExecutableProduct("dos.x", :density_of_states),
    ExecutableProduct("ibrav2cell.x", :ibrav_to_cell),
    ExecutableProduct("kpoints.x", :kpoints),
    ExecutableProduct("cp.x", :carparinello),
    ExecutableProduct("ph.x", :phonon),
    ExecutableProduct("q2r.x", :reciprocal_to_real),
    ExecutableProduct("matdyn.x", :dynamical_matrix_generic),
    ExecutableProduct("dynmat.x", :dynamical_matrix_gamma),
    ExecutableProduct("hp.x", :hubbardparams),
    ExecutableProduct("neb.x", :nudged_elastic_band),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
    Dependency("FFTW_jll"),
    Dependency("Libxc_jll"),
    Dependency("MPICH_jll"),
    Dependency(PackageSpec(name="OpenBLAS32_jll", uuid="656ef2d0-ae68-5445-9ca0-591084a874a2")),
    Dependency("SCALAPACK_jll"),
]

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

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
using Test
using PackageCompatUI: compute_compat, get_julia_versions

@testset "compute_compat" begin
    versions = ["0.0.1", "0.0.2", "0.1.0", "0.1.1", "0.2.0",
                "1.0.0", "1.0.1", "1.0.2", "1.0.3",
                "1.1.0", "1.1.1", "1.2.0", "1.3.0", "1.3.1",
                "1.4.0", "1.4.1", "1.4.2",
                "1.5.0", "1.5.1", "1.5.2", "1.5.3"]

    @test compute_compat(versions, versions[1:1]) == "0.0.1"
    @test compute_compat(versions, versions[1:2]) == "0.0.1, 0.0.2"
    @test compute_compat(versions, versions[3:3]) == "=0.1.0"
    @test compute_compat(versions, versions[3:4]) == "0.1"
    @test compute_compat(versions, versions[3:5]) == "0.1, 0.2"
    @test compute_compat(versions, versions[6:6]) == "=1.0.0"
    @test compute_compat(versions, versions[6:9]) == "~1.0"
    @test compute_compat(versions, versions[6:10]) == "~1.0, =1.1.0"
    @test compute_compat(versions, versions[6:11]) == "~1.0, ~1.1"
    @test compute_compat(versions, versions[12:21]) == "1.2"
    @test compute_compat(versions, versions[14:21]) == "1.3.1"
    @test compute_compat(versions, versions[[1; 3; 4; 16:end]]) == "0.0.1, 0.1, 1.4.1"
end

@testset "get_julia_versions" begin
    # Julia versions as of 2020-11-15.
    julia_versions = ["1.0.0", "1.0.1", "1.0.2", "1.0.3", "1.0.4", "1.0.5",
                      "1.1.0", "1.1.1", "1.2.0", "1.3.0", "1.3.1",
                      "1.4.0", "1.4.1", "1.4.2",
                      "1.5.0", "1.5.1", "1.5.2", "1.5.3"]
    @test issubset(julia_versions, get_julia_versions())
end


text:
using MPI
using MPIStragglers

MPI.Init()
const comm = MPI.COMM_WORLD
const root = 0
const rank = MPI.Comm_rank(comm)
# const nworkers = MPI.Comm_size(comm) - 1
isroot() = MPI.Comm_rank(comm) == root
const data_tag = 0
const control_tag = 1
const nelems = 100 # number of 64-bit floats sent in each direction per iteration and worker

function shutdown()
    for i in 1:(MPI.Comm_size(comm)-1)
        MPI.Isend(zeros(1), i, control_tag, comm)
    end
end

function root_main()

    # parse command-line-arguments
    nworkers = MPI.Comm_size(comm) - 1
    println((nworkers, rank, root))
    print("starting with arguments $ARGS\n")
    nsamples = parse(Int, ARGS[1])
    0 < nsamples || throw(ArgumentError("The number of samples must be positive"))
    filename = ARGS[2]

    # setup
    sendbuf = Vector{Float64}(undef, nelems)
    recvbuf = Vector{Float64}(undef, (nworkers+1)*nelems)

    # warmup
    print("warming up\n")
    MPI.Bcast!(sendbuf, root, comm)
    MPI.Gather!(sendbuf, recvbuf, root, comm)

    # output file
    f = open(filename, "a")
    write(f, "nworkers, nMB_m2w, nMB_w2m, rtt\n")
    nMB_m2w = nelems * 64/8 / 1e6 # MB sent from master to workers
    nMB_w2m = nelems * 64/8 / 1e6 # MB sent from workers to master

    # collect rtt samples
    print("starting test\n")
    for _ in 1:nsamples
        rtt = @elapsed begin
            MPI.Bcast!(sendbuf, root, comm)
            MPI.Gather!(sendbuf, recvbuf, root, comm)
        end
        write(f, "$nworkers, $nMB_m2w, $nMB_w2m, $rtt\n")
    end

    print("tests done\n")
    close(f)
    shutdown()
end

function worker_main()
    recvbuf = Vector{Float64}(undef, nelems)
    sendbuf = Vector{Float64}(undef, nelems)
    crreq = MPI.Irecv!(zeros(1), root, control_tag, comm)
    while true

        # check for exit message on control channel
        flag, _ = MPI.Test!(crreq)
        if flag
            break
        end

        # send/receive
        MPI.Bcast!(recvbuf, root, comm)
        sendbuf .= recvbuf
        MPI.Gather!(sendbuf, nothing, root, comm)
    end
end

if isroot()
    root_main()
else
    worker_main()
end
MPI.Barrier(comm)


text:
## PARSE_SPREADSHEET
#
# Log:
# 11.05.2021 - Created function [1]
#
# Authors:
# [1] Alexander Smolka (alexander.smolka@tum.de)

"""
    parse_preadsheet(raw::CDataRaw)
    parse_spreadsheet(path::String)

**Comments**
- Exported spreadsheet data from COMSOL will write the field names as `<name> (<unit>)`.
  Since this format is assumed to determine the unit of the respective field,
  the respective file has to be in that format for this function to work.

!!! note "COMSOL Version"
    5.3.0.223

!!! note "Julia Version"
    1.6.0
"""
function parse_spreadsheet(raw::CDataRaw)
    # iniate n, u, & v
    n = String[]
    u = String[]
    v = Dict{String, Any}()

    # get keys (names of old dictionary)
    old_keys = collect(keys(raw.v))

    # update value dictionary with new names
    for i in 1:length(raw.field)
        # create new key (name) and unit of field by splotting the field string
        new_key, unit = _separatenameunit(old_keys[i])

        # add name and unit to respective vector
        push!(n, new_key)
        push!(u, unit)

        # add values from old dictionary to new dictionary
        v[new_key] = raw.v[old_keys[i]]
    end

    #
    return CDataParsedSpreadsheet(raw.path, raw.meta, n, u, v)
end
parse_spreadsheet(path::String) = parse_spreadsheet(read_spreadsheet(path))

## INTERNAL FUNCTIONS ======================================================= ##
function _separatenameunit(s::String)
    for i in length(s):-1:1
        if s[i] == '('
            return s[1:i-2], s[i+1:end-1]
        end
    end
end


text:
using RegionTrees
using StaticArrays
using PyPlot
using Statistics
using LinearAlgebra
using SparseArrays
using VectorizedRoutines
using ProgressBars
using MAT
using Images

struct ray2d
    origin::Union{SVector{2,Float64},Vector{Float64}}
    dir::Union{SVector{2,Float64},Vector{Float64}}
    inv_direction::Union{SVector{2,Float64},Vector{Float64}}
    sign::Union{Array{Int64,1},SVector{2,Int64}}
end

function ray2d(p1::T, p2::T) where {T}
    dir::T = p2 - p1
    invdir::T = 1.0 ./ dir
    sign = (invdir .< 0) .* 1
    return ray2d(p1, dir, invdir, sign)
end

function ray2ddir(p1::T, dir::T) where {T}
    invdir::T = 1.0 ./ dir
    sign = (invdir .< 0) .* 1
    return ray2d(p1, dir, invdir, sign)
end

@inline @inbounds function intersects(
    r::ray2d,
    b::HyperRectangle{2,Float64},
)::Float64

    if (r.inv_direction[1] >= 0.0)
        #tmin = (parameters[1][1] - r.origin[1]) * r.inv_direction[1]
        tmin = (b.origin[1] - r.origin[1]) * r.inv_direction[1]
        #tmax = (parameters[2][1] - r.origin[1]) * r.inv_direction[1]
        tmax = (prevfloat(b.origin[1] + b.widths[1]) - r.origin[1]) * r.inv_direction[1]
    else
        #tmin = (parameters[2][1] - r.origin[1]) * r.inv_direction[1]
        tmin = ((b.origin[1] + b.widths[1]) - r.origin[1]) * r.inv_direction[1]
        #tmax = (parameters[1][1] - r.origin[1]) * r.inv_direction[1]
        tmax = (b.origin[1] - r.origin[1]) * r.inv_direction[1]
    end

    if (r.inv_direction[2] >= 0.0)
        #tymin = (parameters[1][2] - r.origin[2]) * r.inv_direction[2]
        tymin = (b.origin[2]  - r.origin[2]) * r.inv_direction[2]
        #tymax = (parameters[2][2]- r.origin[2]) * r.inv_direction[2]
        tymax = (prevfloat(b.origin[2] + b.widths[2]) - r.origin[2]) * r.inv_direction[2]
    else
        #tymin = (parameters[2][2] - r.origin[2]) * r.inv_direction[2]
        tymin = ((b.origin[2] + b.widths[2]) - r.origin[2]) * r.inv_direction[2]
        #tymax = (parameters[1][2] - r.origin[2]) * r.inv_direction[2]
        tymax = (b.origin[2]  - r.origin[2]) * r.inv_direction[2]
    end
    if ((tmin > tymax) || (tymin > tmax))
        return -1.0
    end
     if (tymin > tmin)
         tmin = tymin
     end
    if (tymax < tmax)
         tmax = tymax
    end

    if (isnan(tmin))
        #error("NaN detected")
       return b.widths[1]
    end

    if (isnan(tymin))
        #error("NaN detected")
       return  b.widths[2]
    end

    return norm(r.dir*(tmax-tmin))

end

@inline function intersects(
    a::HyperRectangle{N,Float64},
    b::HyperRectangle{N,Float64},
) where {N}
    maxA = a.origin + a.widths
    minA = a.origin
    maxB = b.origin + b.widths
    minB = b.origin
    return (all(maxA .>= minB) && all(maxB .> minA))

end

function inittree2d(n;lowx= -1.0, lowy = -1.0, wx = 2.0, wy = 2.0)
    empty::Array{Int64,1} = []
    oct = Cell(SVector(lowx, lowy), SVector(wx, wx), empty)
    for _ = 1:n
        for leaf in allleaves(oct)
            split!(leaf)
        end
    end
    for node in allcells(oct)
        node.data = copy(node.data)
    end
    return oct
end

function inserttotree!(
    box::HyperRectangle{N,Float64},
    tree::Cell,
    index::Int64,
) where {N}
    if (isleaf(tree))
        if (intersects(box, tree.boundary))
            push!(tree.data, index)
        end
    elseif (intersects(box, tree.boundary))
        for i = 1:length(tree.children)
            inserttotree!(box, tree[i], index)
        end
    end
    return nothing
end

function possiblepixels(r::ray2d, tree::Cell)::Vector{Int64}
    if (isleaf(tree) && intersects(r, tree.boundary) > 0.0)
        return tree.data

    elseif (~(tree.parent === nothing) && intersects(r, tree.boundary) > 0.0)
        N = length(tree.children)
        v = Vector{Int64}()
        for i = 1:N
            append!(v, possiblepixels(r, tree[i]))
        end
        return v

    elseif ((tree.parent === nothing) && intersects(r, tree.boundary) > 0.0)
        N = length(tree.children)
        v = Vector{Int64}()
        for i = 1:N
            append!(v, possiblepixels(r, tree[i]))
        end
        return collect(Set(v))

    else
        return []
    end

end

function pixelray(
    r::ray2d,
    vv::Vector{HyperRectangle{2,Float64}},
    checklist::Vector{Int64},
)
    N = length(checklist)
    indices = Vector{Int64}()
    tlens = Vector{Float64}()
    for i = 1:N
        q = intersects(r, vv[checklist[i]])
        if (q>0.0)
            push!(tlens, q)
            push!(indices, checklist[i])
        end
    end
    return (indices, tlens)
end

function constructmatrix(tree::Cell, vv::Vector{HyperRectangle{2,Float64}}, Nrays::Int64, theta)
    r = sqrt(2) # Radius of the transform.
    Nproj = length(theta)
    if (~isa(theta, Union{SArray,Array}))
        theta = [theta]
    end

    span = range(r, stop = -r, length = Nrays)

    Ntotal = Nrays * Nproj
    rows = Vector{Int64}()
    cols = Vector{Int64}()
    vals = Vector{Float64}()
    Nofpixels = length(vv)

    Nth = Threads.nthreads()
    rows = Vector{Vector{Int64}}(undef, Nth)
    cols = Vector{Vector{Int64}}(undef, Nth)
    vals = Vector{Vector{Float64}}(undef, Nth)
    for p = 1:Nth
        rows[p] = []
        cols[p] = []
        vals[p] = []
    end

    t = time()
    if (Nproj > 1)
        pb = ProgressBar(1:Nproj)
    else
        pb = 1:1
    end
    for a in pb
        dir = [cos(theta[a]), sin(theta[a])]
        aux = [-sin(theta[a]), cos(theta[a])]
        Threads.@threads for i = 1:Nrays
            rayindex = (a - 1) * Nrays + i
            p1 = dir + aux * span[i]
            p2 = p1 - 2 * dir * r
            or = ray2d(p1, p2)
            checklist = possiblepixels(or, tree) #Possible pixels.
            indices, tlens = pixelray(or, vv, checklist)
            Nel = length(indices)

            append!(rows[Threads.threadid()],rayindex * ones(Int64, length(tlens)))
            append!(cols[Threads.threadid()], indices)
            append!(vals[Threads.threadid()], tlens)

        end
    end
    rows = vcat(rows...)
    cols = vcat(cols...)
    vals = vcat(vals...)
    M = sparse(rows, cols, vals, Ntotal, Nofpixels)
    println("Matrix constructed in ", time() - t, " seconds.")
    return M
end

function constructmatrix(tree::Cell, vv::Vector{HyperRectangle{2,Float64}}, rays::Array{Array{ray2d,1},1})
    Nproj = length(rays)
    for i = 1:Nproj-1
        @assert length(rays[i]) == length(rays[i+1])
    end
    Nrays = length(rays[1])
    Ntotal = Nrays * Nproj
    rows = Vector{Int64}()
    cols = Vector{Int64}()
    vals = Vector{Float64}()
    Nofpixels = length(vv)

    Nth = Threads.nthreads()
    rows = Vector{Vector{Int64}}(undef, Nth)
    cols = Vector{Vector{Int64}}(undef, Nth)
    vals = Vector{Vector{Float64}}(undef, Nth)
    for p = 1:Nth
        rows[p] = []
        cols[p] = []
        vals[p] = []
    end

    t = time()
    if (Nproj > 1)
        pb = ProgressBar(1:Nproj)
    else
        pb = 1:1
    end
    for a in pb
        Threads.@threads for i = 1:Nrays
            rayindex = (a - 1) * Nrays + i
            or = rays[a][i]
            checklist = possiblepixels(or, tree) #Possible pixels.
            indices, tlens = pixelray(or, vv, checklist)
            Nel = length(indices)

            append!(rows[Threads.threadid()],rayindex * ones(Int64, length(tlens)))
            append!(cols[Threads.threadid()], indices)
            append!(vals[Threads.threadid()], tlens)

        end
        # if (a==Nproj)
        #     cols = vcat(cols...)
        #     ix = zeros(256^2,)
        #     ix[cols] .= 1
        #     figure()
        #     imshow(reshape(ix,256,256))
        #     error("")
        # end
    end
    rows = vcat(rows...)
    cols = vcat(cols...)
    vals = vcat(vals...)
    M = sparse(rows, cols, vals, Ntotal, Nofpixels)
    println("Matrix constructed in ", time() - t, " seconds.")
    return M
end

function sub2ind(sizes, i::Int64, j::Int64, k::Int64)
    @assert i > 0 && i <= sizes[1]
    @assert j > 0 && j <= sizes[2]
    @assert k > 0 && k <= sizes[3]
    return (k - 1) * sizes[1] * sizes[2] + sizes[1] * (j - 1) + i
end

function setuppixels(sizex::Int64,sizey::Int64,octsize::Int64;lowy = -1.0, lowx = -1.0, widthy = 2, widthx = 2)
    @assert widthy == widthx
    oct = inittree2d(octsize;lowx= lowx, lowy = lowy, wx = widthx, wy = widthy)
    size = [sizex, sizey]
    cellsize = [widthx,widthy] ./ size
    wds = [widthx,widthy]
    width = @SVector[cellsize[1], cellsize[2]]
    #l = Array([ [0.0, 0] [1, 0] [1, 1] [0, 1.0]]')
    #pm = mean(l, dims = 1)
    # for i = 1:2
    #     l[:, i] = cellsize[i] .* (l[:, i] .- pm[i]) .- cellsize[i] / 2.0 * size[i] .+ 0.5 * cellsize[i]
    # end
    l = [[lowx, lowy] [lowx + cellsize[1], lowy ] [lowx , lowy + cellsize[2]] [lowx + cellsize[1], lowy + cellsize[2]] ]'
    # println(l)
    # error("")
    pixelvector = Vector{HyperRectangle{2,Float64}}(undef, sizex * sizey)
    t = time()

    for j = 1:sizey
        for i = 1:sizex
            ll = l .+ reshape(cellsize .* [i - 1,  sizey - j + 1], 1, 2) # Flip Y-axis so that we have canonical Euclidean coordinate system.
            origin = @SVector[
                minimum(ll[:, 1]),
                minimum(ll[:, 2]),
            ]
            pixel = HyperRectangle(origin, width)
            index = j + (i - 1) * sizey
            inserttotree!(pixel, oct, index)
            pixelvector[index] = pixel
        end
    end

    println("Pixels initialized in ", time() - t, " seconds.")
    return (oct, pixelvector)
end

function constructparallelrays(Nrays,Nproj;rotations=range(-pi/2,stop=3/4*2*pi,length=Nproj),dete_plate_span=range(-sqrt(2), stop = sqrt(2), length = Nrays))

    @assert length(dete_plate_span) == Nrays
    @assert length(rotations) == Nproj
    rays = Vector{Vector{ray2d}}(undef,Nproj)

    # Parallel beam.

    for i = 1:Nproj
        span = dete_plate_span#range(r, stop = -r, length = Nrays)
        rays[i] = Vector{ray2d}(undef,Nrays)
        for j = 1:Nrays
            dir = [cos(rotations[i]), sin(rotations[i])]
            #aux = [-sin(rotations[i]), cos(rotations[i])]
            aux = [-sin(rotations[i]), cos(rotations[i])]*span[j] + center
            p1 =  aux #* span[j]
            #p2 = p1 - 2 * dir * r
            #ray = ray2d(p1, p2)
            ray = ray2ddir(p1,dir)
            rays[i][j] = ray

            #plot([ray.origin[1], ray.origin[1] +  ray.dir[1]],[ray.origin[2], ray.origin[2] +  ray.dir[2]] )
           # plot([rays[i][j].origin[1], rays[i][j].origin[1] +  rays[i][j].dir[1]],[rays[i][j].origin[2], rays[i][j].origin[2] +  rays[i][j].dir[2]] )

        end
    end

    return rays

end

function constructfanrays(Nrays,Nproj;src_shift_func = a->[0,0.0],translation=[0,0.0],src_to_det_init=[0,1],det_axis_init = nothing,det_shift_func=a->[0,0.0],apart=range(0,stop=pi,length=Nproj),det_radius=sqrt(2),src_radius=sqrt(2),dpart=range(-sqrt(2), stop = sqrt(2), length = Nrays))

    @assert length(dpart) == Nrays
    @assert length(apart) == Nproj
    rays = Vector{Vector{ray2d}}(undef,Nproj)

    sdinit = 2*pi*1/4+atan(src_to_det_init[2],src_to_det_init[1])
    apart = apart .+ sdinit
    translation = [translation[2], -translation[1]]
    dextra_rot = 0.0

    if (det_axis_init === nothing)
        dextra_rot = 0.0
    else
        dinit = 2*pi*1/4+atan(det_axis_init[2],det_axis_init[1])
        dextra_rot = mod((2*pi*1/4-(sdinit-dinit)),2*pi)
    end

    span = dpart

    for i = 1:Nproj

        rays[i] = Vector{ray2d}(undef,Nrays)
        forplots = zeros(2,2)

        src_tangent = [-sin(apart[i]), cos(apart[i])]
        src_orth = [cos(apart[i]), sin(apart[i])]
        srcshift = src_shift_func(apart[i])
        source = [cos(apart[i]), sin(apart[i])]*src_radius + srcshift[1]*src_orth + srcshift[2]*src_tangent
        p2 = source + translation

        detcenter = -[cos(apart[i]), sin(apart[i])]*det_radius
        det_tangent = [-sin(apart[i]), cos(apart[i])]
        det_orth = [cos(apart[i]), sin(apart[i])]
        detshift = det_shift_func(apart[i])
        totcenter = detcenter - detshift[1]*det_orth - detshift[2]*det_tangent

        for j = 1:Nrays
            aux = [-sin(apart[i]+dextra_rot), cos(apart[i]+dextra_rot)]*span[j] + totcenter
            p1 = aux + translation
            ray = ray2d(p1, p2)
            rays[i][j] = ray

            # if(j==1)
            #     forplots[1,:] = p1
            # elseif(j==Nrays)
            #     forplots[2,:] = p1
            # end

            # plot([rays[i][j].origin[1], rays[i][j].origin[1] +  rays[i][j].dir[1]],[rays[i][j].origin[2], rays[i][j].origin[2] +  rays[i][j].dir[2]] )
            # scatter(aux[1],aux[2])

        end
        ## plot([rays[i][1].origin[1], rays[i][1].origin[1] +  rays[i][1].dir[1]],[rays[i][1].origin[2], rays[i][1].origin[2] +  rays[i][1].dir[2]] )
        ## plot([rays[i][end].origin[1], rays[i][end].origin[1] +  rays[i][end].dir[1]],[rays[i][end].origin[2], rays[i][end].origin[2] +  rays[i][end].dir[2]] )
        # plot(forplots[:,1],forplots[:,2])
        # println(forplots)
        # scatter(p2[1],p2[2])
        # xlim([-9,9])
        # ylim([-9,9])
        # axis(:equal)
     end
    #  error("")
    return rays
end

function test()

    close("all")

    #(qt,pixelvector)=setuppixels(Nx,Ny,Os)   # Default grid span is a rectangle with extents of [-1,-1] and [1,1].
    # fn = "r.mat"
    # file = matopen(fn)
    # zn = read(file, "S")
    # zn = zn[:,1]

    fn = "simplephantom.png"
    zn = load(fn)
    zn = channelview(zn)
    zn = zn[:]

    ## When compared to ODL, the logic of the detector span is the same. The span vector should be increasing, so
    ## detector_partition = odl.uniform_partition(-8, 4, 100) equals dete_plate = range(-8, stop = 4, length = 100)

    ## However, with the coordinate axis are flipped in ODL. Depending on the case, π/2 radians plus a
    ## possbile angle shift must be added to obtain the same
    ## rotational span due to the inverted coordinate system.
    ## This is very important if the code is modified, but the current fanbeam function takes care of the axis differences.

    ## The vector src_to_det_init defines the shift to the initial angle:
    ## angle_partition = odl.uniform_partition(0,2*np.pi, 360)
    ## equals rotations = range(0,stop=2*pi,length=360) .+ 2*pi*1/4+atan(src_to_det_init[2],src_to_det_init[1])
    ## Again, this extra angle is added due to the flipped coordinate axis between ODL and this code.

    ## ODL's parameter det_axis_init=[1,-0.5] refers to extra rotation of the detector plate. That is,
    ## det_axis_init=[3,-0.7] equals dextra_rot = atan(-0.7,3)
    ## if the src_to_det_init is left to its default value.

    ## Source and detector radii have the same logic.
    ## By default, the center of rotation in ODL is the origin. The vector translation
    ## moves it. This is just added to  initial points of the rays and the detector plate
    ## position vectors.

    ## The function det_shift_func defines a shift for the detector at each projection angle
    ## in the parallel direction of the ray that goes through the COR and in the tangent direction of the ray.
    ## Function src_shift_func works also. It work in the same manner as det_shift_func, it just determines a shift for the source.
    ## Finally, having non-constant shift functions do not seem to behave properly in ODL. Perhaps there is a bug in ODL?

    ## Curved detector plates are not implemented.

    ## The domain of the object is odl_space and its correspondence to ODL should be clear. Since Astra toolbox and ODL do not seem to support
    ## anisotropic pixels together, the  number of pixels of the object and the geometry dimensions should be the same.

    Nx = 256; Ny = 256;  # Number of X and Y pixels,
    Os = 6 # Splitting factor of the quadtree. Does not affect the theory matrix, only the performance of building it. 6 seems optimal for 256x256.
    Nproj = 360 # Number of projections i.e. no. angles
    Nrays = 400 # Number of rays in a projection

    odl_space = (min_pt=[-1.4, -2], max_pt=[3.6, 3])
    src_to_det_init  = [1,5.5] #[0,1] is the default
    det_axis_init =   [-1,1.3] #[1,0] is the default
    det_shift_func(angle) = [-1.5, 0.4]
    src_shift_func(angle) = [0.7,1.5]
    rotations=  range(pi/3,stop=2*pi/3,length=Nproj)
    translation = [0.5,0.9]

    dete_radius = 5.7
    source_radius = 5.0
    dete_plate = range(-6, stop = 10, length = Nrays)

    (qt,pixelvector)=setuppixels(Nx,Ny,Os;lowx=odl_space.min_pt[2],lowy=-odl_space.max_pt[1],widthx=abs(odl_space.max_pt[2]-odl_space.min_pt[2]),widthy=abs(odl_space.max_pt[1]-odl_space.min_pt[1]))
    rays = constructfanrays(Nrays,Nproj;translation=translation, src_shift_func=src_shift_func, src_to_det_init = src_to_det_init, det_axis_init=det_axis_init, det_shift_func=det_shift_func,apart=rotations,det_radius=dete_radius,src_radius=source_radius, dpart=dete_plate)

    #imshow(reshape(zn[:,1],256,256),extent=[-1,1,-1,1])
    #xlim([-5,5])
    #ylim([-5,5])

    M = constructmatrix(qt, pixelvector, rays)

    yp = M * zn

    # For simulating a rotating object, but stationary device
    # Nlog = size(zn)[2]
    # rad = zeros(Nrays,Nlog)
    # for i = 1:Nlog
    #     p = zn[:,i]
    #     yp = M * p
    #     rad[:,i] = yp
    # end

    figure()
    rad = reshape(yp,  Nrays,Nproj)
    println(size(rad))
    imshow(rad, aspect = "auto")
    title("Julia")
    # figure()
    # imshow(reshape(zn,256,256))
    return rad,M
end

sinogram,M=test()

# Extremely simple reconstruction.
# using Pkg
# Pkg.add("RandomizedLinAlg")
# using RandomizedLinAlg
# Q = rsvd(M, 500)
# y = sinogram[:]
# x = (sparse(Q.Vt')*(spdiagm(0=>1.0./Q.S))*(sparse(Q.U')*y));
# figure()
# imshow(reshape(x,256,256))

#################################################################################
## ODL code for comparison

# import numpy as np
# import odl
# import scipy.io
# from PIL import Image
# import matplotlib.pyplot as plt

# plt.close("all")

# reco_space = odl.uniform_discr(min_pt=[-1.4, -2], max_pt=[3.6, 3], shape=[256, 256], dtype='float32')

# angle_partition = odl.uniform_partition(np.pi/3,2*np.pi/3, 360)
# detector_partition = odl.uniform_partition(-6, 10, 400)
# # If the first param of the det_shift_func is constant, the effect is the same as increasing or decreasing the detector radius.
# # Positive values mean decreasing the radius, negative values decrease it.

# dsfun = lambda angle:  np.array([[-1.5, 0.4]])
# ssfun = lambda angle: np.array([[0.7, 1.5]])
# geometry = odl.tomo.geometry.conebeam.FanBeamGeometry(angle_partition, detector_partition, src_shift_func=ssfun,  det_shift_func = dsfun, translation = [0.5,0.9], src_to_det_init=[1,5.5] , det_axis_init=[-1,1.3], src_radius=5, det_radius=5.7)

# # Ray transform (= forward projection).
# ray_trafo = odl.tomo.RayTransform(reco_space, geometry, impl='astra_cpu')

# mat = Image.open("simplephantom.png").convert('L')
# phantom = np.array(mat)/255.0
# proj_data = ray_trafo(phantom)
# p = proj_data.data.T

# print(p.shape)
# plt.figure()
# plt.imshow(p,aspect = "auto")
# plt.title("ODL")
# #plt.plot(p)
# plt.show()

max_stars_repo_name: suurj/radon

max_issues_repo_name: suurj/radon

max_forks_repo_name: suurj/radon


text:
struct CoordMode{D} <: IndexMode
    dims::D
end

dims(m::CoordMode) = m.dims
order(m::CoordMode) = Unordered()

struct Coord{T,Mo<:IndexMode,Me<:AllMetadata} <: Dimension{T,Mo,Me}
    val::T
    mode::Mo
    metadata::Me
end
function Coord(val::T, dims::Tuple, metadata::Me=NoMetadata()) where {T<:AbstractVector,Me<:AllMetadata}
    length(dims) = length(first(val))
    mode = CoordMode(key2dim(dims))
    Coord{T,typeof(mode),Me}(val, mode, metadata)
end
const SelOrStandard = Union{Selector,StandardIndices}
Coord(s1::SelOrStandard, s2::SelOrStandard, sels::SelOrStandard...) = Coord((s1, s2, sels...))
Coord(sel::Selector) = Coord((sel,))
Coord(d1::Dimension, dims::Dimension...) = Coord((d1, dims...))
Coord(val::T=:) where T = Coord{T,AutoMode,NoMetadata}(val, AutoMode(), NoMetadata())

dims(d::Coord) = dims(mode(d))
bounds(d::Coord) = ntuple(i -> extrema((x[i] for x in val(d))), length(first(d)))

# Return a Vector{Bool} for matching coordinates
sel2indices(dim::Coord, sel::DimTuple) = sel2indices(dim, sortdims(sel, dims(dim)))
sel2indices(dim::Coord, sel::Tuple) = [all(map(_matches, sel, x)) for x in val(dim)]
sel2indices(dim::Coord, sel::StandardIndices) = sel

_matches(sel::Dimension, x) = _matches(val(sel), x)
_matches(sel::Between, x) = (x >= first(sel)) & (x < last(sel))
_matches(sel::At, x) = x == val(sel)
_matches(sel::Colon, x) = true
_matches(sel::Nothing, x) = true

_tozero(xs) = map(x -> zero(x), xs)
@inline _reducedims(mode::CoordMode, dim::Dimension) =
    rebuild(dim, [_tozero(dim.val[1])], dim.mode)


text:
module FakeServer

using HTTP
using Sockets

abstract type AbstractFakeServer end;
abstract type AbstractFakeResource end;
abstract type AbstractConnection end;

export EchoResource, InteractiveResource, ErrorResource
export waitnconnected, waitnwritten, waitnreceived

struct DefaultServer <: AbstractFakeServer
    tcpserver::Sockets.TCPServer
    task::Task
    default::AbstractFakeResource
    resources::Dict{AbstractString, AbstractFakeResource}
end
DefaultServer(tcpserver, task, default) = DefaultServer(tcpserver, task, default, Dict{AbstractString, AbstractFakeServer}())

Base.get(server::DefaultServer, key::AbstractString) = get(server.resources, key,server.default)
Base.length(server::DefaultServer) = length(server.resources)
Base.getindex(server::DefaultServer, key::AbstractString) = getindex(server.resources, key)
Base.setindex!(server::DefaultServer,resource::AbstractFakeResource, key::AbstractString) = setindex!(server.resources, resource, key)
abstract type AbstractResourceType end;
struct Echo <: AbstractResourceType end
struct Error <: AbstractResourceType end
struct Interactive <: AbstractResourceType end

struct FakeResource{T <: AbstractResourceType}  <: AbstractFakeResource
    connections::Array{AbstractConnection}
    option
end
FakeResource{T}(option=nothing) where T = FakeResource{T}(AbstractConnection[], option)
Base.length(resource::FakeResource) = length(resource.connections)
Base.getindex(resource::FakeResource,args...) = getindex(resource.connections, args...)
Base.push!(resource::FakeResource, connection::AbstractConnection) = push!(resource.connections, connection)

waitnconnected(resource::FakeResource, n::Integer; timeout::Real=3.0) = timedwait(timeout) do
    length(resource.connections)>=n
end === :ok

EchoResource(;greeting::AbstractString="")=FakeResource{Echo}(greeting)
ErrorResource(;statuscode::UInt16=UInt16(404))=FakeResource{Error}(statuscode)
InteractiveResource(;autoclose::Bool=true)=FakeResource{Interactive}(autoclose)

mutable struct Connection <: AbstractConnection
    io::IO
    received::Array{Any}
    nwritten::Integer
    channel::Channel{Any}
end

Connection(io::IO)=Connection(io,Any[],0, Channel(32))
waitnwritten(conn::Connection, n::Integer; timeout::Real=3.0) = timedwait(timeout) do
    conn.nwritten >= n
end === :ok

waitnreceived(conn::Connection, n::Integer; timeout::Real=3.0) = timedwait(timeout) do
    length(conn.received) >= n
end === :ok

function listen(port::UInt16, args...;
                listenfn::Function=HTTP.WebSockets.listen,
                defaultresource=EchoResource(), kwargs...)
    tcpserver = Sockets.listen(ip"127.0.0.1", port)
    task = @async try
                    s = wait()
                    listenfn("127.0.0.1", port, args...;server=tcpserver, kwargs...) do io
                            serve(s,io)
                        end
                catch e
                    @debug "caught error" e
                    if !isa(e, Base.IOError)
                        rethrow()
                    end
                end
    server =  DefaultServer(tcpserver, task, defaultresource)
    while !istaskstarted(task)
        yield()
    end
    schedule(task, server)
    return server
end

function serve(server::AbstractFakeServer, io::IO)
    #each serve() is in an independent coroutine.
    #we need to catch and show errors to make debuging easier.
    try
        target = gettarget(io)
        resource = get(server, target)
        serve(resource, io)
    catch e
        @error "caught error in serve" e
        rethrow()
    end
end

function serve(resource::AbstractFakeResource, io::IO)
    connection = Connection(io)
    push!(resource, connection)
    handle(resource, connection)
end

function handle(resource::FakeResource{Error},connection::Connection)
    closewrite(connection.io; statuscode=resource.option)
end

function handle(resource::FakeResource{Echo}, connection::Connection)
    if length(resource.option) > 0
        write(connection, resource.option)
    end
    while !eof(connection.io)
        data = readavailable(connection)
        if length(data) > 0
            write(connection, data)
        end
    end
end

function handle(resource::FakeResource{Interactive}, connection::Connection)
    write_task = handle_write(connection)
    handle_read(connection)
    if resource.option
        put!(connection, nothing)
    end
    wait(write_task)
end

function handle_write(connection::Connection)
    @async try
    #catch exception in async
        while isopen(connection.channel)
            data = take!(connection.channel)
            if data === nothing
                break
            end
            write(connection, data)
        end
    catch e
        @error "caught error in handle write" e
    finally
        close(connection.channel)
    end
end

function handle_read(connection::Connection)
    try
        while !eof(connection.io)
            readavailable(connection)
        end
    catch e
        @error "caught error in handle_read(::Connection)" e
    end
end

function Base.readavailable(conn::Connection)
    data = readavailable(conn.io)
    if length(data) > 0
        push!(conn.received, data)
    end
    return data
end

function Base.write(conn::Connection, data)
    write(conn.io, data)
    conn.nwritten += 1
end

Base.put!(conn::Connection, x) = put!(conn.channel, x)

gettarget(ws::HTTP.WebSockets.WebSocket) = ws.request.target
gettarget(http::HTTP.Stream) = http.message.target

function Base.close(server::DefaultServer)
    try
        close(server.tcpserver)
        wait(server.task)
    catch e
        @error "caught error in stop(::DefaultServer" e
        rethrow()
    end
end

end #module FakeServer


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test, Core.IR

function f22938(a, b, x...)
    nothing
    nothing
    nothing
    nothing
    nothing
    nothing
    nothing
    d = 1
    a = d
    for i in 1:b
        d += i
    end
    return i * a
end

msig = Tuple{typeof(f22938),Int,Int,Int,Int}
world = typemax(UInt)
_, msp, m = Base._methods_by_ftype(msig, -1, world)[]
mi = Core.Compiler.specialize_method(m, msig, msp, false)
c0 = Core.Compiler.retrieve_code_info(mi)

@test isempty(Core.Compiler.validate_code(mi))
@test isempty(Core.Compiler.validate_code(c0))

@testset "INVALID_EXPR_HEAD" begin
    c = copy(c0)
    c.code[1] = Expr(:invalid, 1)
    errors = Core.Compiler.validate_code(c)
    @test length(errors) == 1
    @test errors[1].kind === Core.Compiler.INVALID_EXPR_HEAD
end

@testset "INVALID_LVALUE" begin
    c = copy(c0)
    c.code[1] = Expr(:(=), GotoNode(1), 1)
    c.code[2] = Expr(:(=), :x, 1)
    c.code[3] = Expr(:(=), 3, 1)
    errors = Core.Compiler.validate_code(c)
    @test length(errors) == 3
    @test all(e.kind === Core.Compiler.INVALID_LVALUE for e in errors)
end

@testset "INVALID_RVALUE" begin
    c = copy(c0)
    c.code[1] = Expr(:(=), SlotNumber(2), GotoNode(1))
    c.code[2] = Expr(:(=), SlotNumber(2), LineNumberNode(2))
    i = 2
    for h in (:gotoifnot, :line, :const, :meta)
        c.code[i+=1] = Expr(:(=), SlotNumber(2), Expr(h))
    end
    errors = Core.Compiler.validate_code(c)
    @test length(errors) == 6
    @test count(e.kind === Core.Compiler.INVALID_RVALUE for e in errors) == 6
end

@testset "INVALID_CALL_ARG" begin
    c = copy(c0)
    c.code[1] = Expr(:(=), SlotNumber(2), Expr(:call, GlobalRef(Base,:+), SlotNumber(2), GotoNode(1)))
    c.code[2] = Expr(:call, GlobalRef(Base,:-), Expr(:call, GlobalRef(Base,:sin), GotoNode(2)), 3)
    c.code[3] = Expr(:call, LineNumberNode(2))
    i = 3
    for h in (:gotoifnot, :line, :const, :meta)
        c.code[i+=1] = Expr(:call, GlobalRef(@__MODULE__,:f), Expr(h))
    end
    errors = Core.Compiler.validate_code(c)
    @test length(errors) == 7
    @test count(e.kind === Core.Compiler.INVALID_CALL_ARG for e in errors) == 7
end

@testset "EMPTY_SLOTNAMES" begin
    c = copy(c0)
    empty!(c.slotnames)
    errors = Core.Compiler.validate_code(c)
    @test length(errors) == 2
    @test any(e.kind === Core.Compiler.EMPTY_SLOTNAMES for e in errors)
    @test any(e.kind === Core.Compiler.SLOTFLAGS_MISMATCH for e in errors)
end

@testset "SLOTFLAGS_MISMATCH" begin
    c = copy(c0)
    push!(c.slotflags, 0x00)
    errors = Core.Compiler.validate_code(c)
    @test length(errors) == 1
    @test errors[1].kind === Core.Compiler.SLOTFLAGS_MISMATCH
end

@testset "SSAVALUETYPES_MISMATCH" begin
    c = code_typed(f22938, (Int,Int,Int,Int))[1][1]
    empty!(c.ssavaluetypes)
    errors = Core.Compiler.validate_code(c)
    @test length(errors) == 1
    @test errors[1].kind === Core.Compiler.SSAVALUETYPES_MISMATCH
end

@testset "SSAVALUETYPES_MISMATCH_UNINFERRED" begin
    c = copy(c0)
    c.ssavaluetypes -= 1
    errors = Core.Compiler.validate_code(c)
    @test length(errors) == 1
    @test errors[1].kind === Core.Compiler.SSAVALUETYPES_MISMATCH_UNINFERRED
end

@testset "SIGNATURE_NARGS_MISMATCH" begin
    old_sig = mi.def.sig
    mi.def.sig = Tuple{1,2}
    errors = Core.Compiler.validate_code(mi)
    mi.def.sig = old_sig
    @test length(errors) == 1
    @test errors[1].kind === Core.Compiler.SIGNATURE_NARGS_MISMATCH
end

@testset "NON_TOP_LEVEL_METHOD" begin
    c = copy(c0)
    c.code[1] = Expr(:method, :dummy)
    errors = Core.Compiler.validate_code(c)
    @test length(errors) == 1
    @test errors[1].kind === Core.Compiler.NON_TOP_LEVEL_METHOD
end

@testset "SLOTNAMES_NARGS_MISMATCH" begin
    mi.def.nargs += 20
    errors = Core.Compiler.validate_code(mi)
    mi.def.nargs -= 20
    @test length(errors) == 2
    @test count(e.kind === Core.Compiler.SLOTNAMES_NARGS_MISMATCH for e in errors) == 1
    @test count(e.kind === Core.Compiler.SIGNATURE_NARGS_MISMATCH for e in errors) == 1
end

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    pairwise(γ, domain)

Evaluate variogram `γ` between all elements in the `domain`.
"""
function pairwise(γ::Variogram, domain)
  u = first(domain)
  n = length(domain)
  R = result_type(γ, u, u)
  Γ = Matrix{R}(undef, n, n)
  pairwise!(Γ, γ, domain)
end

function pairwise!(Γ, γ::Variogram, domain)
  n = length(domain)
  @inbounds for j in 1:n
    vⱼ = domain[j]
    sⱼ = _sample(γ, vⱼ)
    for i in j+1:n
      vᵢ = domain[i]
      sᵢ = _sample(γ, vᵢ)
      Γ[i,j] = mean(γ(pᵢ, pⱼ) for pᵢ in sᵢ, pⱼ in sⱼ)
    end
    Γ[j,j] = mean(γ(pⱼ, pⱼ) for pⱼ in sⱼ, pⱼ in sⱼ)
    for i in 1:j-1
      Γ[i,j] = Γ[j,i] # leverage the symmetry
    end
  end
  Γ
end

"""
    pairwise(γ, domain₁, domain₂)

Evaluate variogram `γ` between all elements of `domain₁` and `domain₂`.
"""
function pairwise(γ::Variogram, domain₁, domain₂)
  u = first(domain₁)
  v = first(domain₂)
  m = length(domain₁)
  n = length(domain₂)
  R = result_type(γ, u, v)
  Γ = Matrix{R}(undef, m, n)
  pairwise!(Γ, γ, domain₁, domain₂)
end

function pairwise!(Γ, γ::Variogram, domain₁, domain₂)
  m = length(domain₁)
  n = length(domain₂)
  @inbounds for j in 1:n
    vⱼ = domain₂[j]
    sⱼ = _sample(γ, vⱼ)
    for i in 1:m
      vᵢ = domain₁[i]
      sᵢ = _sample(γ, vᵢ)
      Γ[i,j] = mean(γ(pᵢ, pⱼ) for pᵢ in sᵢ, pⱼ in sⱼ)
    end
  end
  Γ
end


text:
# Julia-specific high-level API

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

export AdiosFile
"""
    struct AdiosFile

Context for the high-level API for an ADIOS file
"""
struct AdiosFile
    adios::Adios
    io::AIO
    engine::Engine
    AdiosFile(adios::Adios, io::AIO, engine::Engine) = new(adios, io, engine)
end

export adios_open_serial
"""
    adios = adios_open_serial(filename::AbstractString, mode::Mode)
    adios::AdiosFile

Open an ADIOS file. Use `mode = mode_write` for writing and `mode =
mode_read` for reading.
"""
function adios_open_serial(filename::AbstractString, mode::Mode)
    adios = adios_init_serial()
    io = declare_io(adios, "IO")
    engine = open(io, filename, mode)
    return AdiosFile(adios, io, engine)
end

export adios_open_mpi
"""
    adios = adios_open_mpi(comm::MPI.Comm, filename::AbstractString, mode::Mode)
    adios::AdiosFile

Open an ADIOS file for parallel I/O. Use `mode = mode_write` for
writing and `mode = mode_read` for reading.
"""
function adios_open_mpi(comm::MPI.Comm, filename::AbstractString, mode::Mode)
    adios = adios_init_mpi(comm)
    io = declare_io(adios, "IO")
    engine = open(io, filename, mode)
    return AdiosFile(adios, io, engine)
end

"""
    flush(file::AdiosFile)

Flush an ADIOS file. When writing, flushing or closing a file is
necesssary to ensure that data are actually written to the file.
"""
function Base.flush(file::AdiosFile)
    flush(file.engine)
    return nothing
end

"""
    close(file::AdiosFile)

Close an ADIOS file. When writing, flushing or closing a file is
necesssary to ensure that data are actually written to the file.
"""
function Base.close(file::AdiosFile)
    close(file.engine)
    adios_finalize(file.adios)
    return nothing
end

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

export adios_subgroup_names
"""
    groups = adios_subgroup_names(file::AdiosFile, groupname::AbstractString)
    vars::Vector{String}

List (non-recursively) all subgroups in the group `groupname` in the
file.
"""
function adios_subgroup_names(file::AdiosFile, groupname::AbstractString)
    vars = inquire_subgroups(file.io, groupname)
    vars ≡ nothing && return String[]
    return vars
end

export adios_define_attribute
"""
   adios_define_attribute(file::AdiosFile, name::AbstractString,
                          value::AdiosType)

Write a scalar attribute.
"""
function adios_define_attribute(file::AdiosFile, name::AbstractString,
                                value::AdiosType)
    define_attribute(file.io, name, value)
    return nothing
end

"""
    adios_define_attribute(file::AdiosFile, name::AbstractString,
                           value::AbstractArray{<:AdiosType})

Write an array-valued attribute.
"""
function adios_define_attribute(file::AdiosFile, name::AbstractString,
                                value::AbstractArray{<:AdiosType})
    define_attribute_array(file.io, name, value)
    return nothing
end

"""
    adios_define_attribute(file::AdiosFile, path::AbstractString,
                           name::AbstractString, value::AdiosType)

Write a scalar attribute into the path `path` in the file.
"""
function adios_define_attribute(file::AdiosFile, path::AbstractString,
                                name::AbstractString, value::AdiosType)
    define_variable_attribute(file.io, name, value, path)
    return nothing
end

"""
    adios_define_attribute(file::AdiosFile, path::AbstractString,
                           name::AbstractString,
                           value::AbstractArray{<:AdiosType})

Write an array-valued attribute into the path `path` in the file.
"""
function adios_define_attribute(file::AdiosFile, path::AbstractString,
                                name::AbstractString,
                                value::AbstractArray{<:AdiosType})
    define_variable_attribute_array(file.io, name, value, path)
    return nothing
end

export adios_all_attribute_names
"""
    attrs = adios_all_attribute_names(file::AdiosFile)
    attrs::Vector{String}

List (recursively) all attributes in the file.
"""
function adios_all_attribute_names(file::AdiosFile)
    attrs = inquire_all_attributes(file.io)
    attrs ≡ nothing && return String[]
    return name.(attrs)
end

export adios_group_attribute_names
"""
    vars = adios_group_attribute_names(file::AdiosFile, groupname::AbstractString)
    vars::Vector{String}

List (non-recursively) all attributes in the group `groupname` in the
file.
"""
function adios_group_attribute_names(file::AdiosFile, groupname::AbstractString)
    vars = inquire_group_attributes(file.io, groupname)
    vars ≡ nothing && return String[]
    return name.(vars)
end

export adios_attribute_data
"""
    attr_data = adios_attribute_data(file::AdiosFile, name::AbstractString)
    attr_data::Union{Nothing,AdiosType}

Read an attribute from a file. Return `nothing` if the attribute is
not found.
"""
function adios_attribute_data(file::AdiosFile, name::AbstractString)
    attr = inquire_attribute(file.io, name)
    attr ≡ nothing && return nothing
    attr_data = data(attr)
    return attr_data
end

"""
    attr_data =  adios_attribute_data(file::AdiosFile, path::AbstractString,
                                      name::AbstractString)
    attr_data::Union{Nothing,AdiosType}

Read an attribute from a file in path `path`. Return `nothing` if the
attribute is not found.
"""
function adios_attribute_data(file::AdiosFile, path::AbstractString,
                              name::AbstractString)
    attr = inquire_variable_attribute(file.io, name, path)
    attr ≡ nothing && return nothing
    attr_data = data(attr)
    return attr_data
end

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

export adios_put!
"""
   adios_put!(file::AdiosFile, name::AbstractString, scalar::AdiosType)

Schedule writing a scalar variable to a file.

The variable is not written until `adios_perform_puts!` is called and
the file is flushed or closed.
"""
function adios_put!(file::AdiosFile, name::AbstractString, scalar::AdiosType)
    var = define_variable(file.io, name, scalar)
    put!(file.engine, var, scalar)
    return var
end

"""
    adios_put!(file::AdiosFile, name::AbstractString,
               array::AbstractArray{<:AdiosType}; make_copy::Bool=false)

Schedule writing an array-valued variable to a file.

`make_copy` determines whether to make a copy of the array, which is
expensive for large arrays. When no copy is made, then the array must
not be modified before `adios_perform_puts!` is called.

The variable is not written until `adios_perform_puts!` is called and
the file is flushed or closed.
"""
function adios_put!(file::AdiosFile, name::AbstractString,
                    array::AbstractArray{<:AdiosType}; make_copy::Bool=false)
    # 0-dimensional arrays need to be passed as scalars
    ndims(array) == 0 &&
        return adios_put!(file, name, make_copy ? copy(array)[] : array[])
    var = define_variable(file.io, name, array)
    put!(file.engine, var, make_copy ? copy(array) : array)
    return var
end

export adios_perform_puts!
"""
    adios_perform_puts!(file::AdiosFile)

Execute all scheduled `adios_put!` operations.

The data might not be in the file yet; they might be buffered. Call
`adios_flush` or `adios_close` to ensure all data are written to file.
"""
function adios_perform_puts!(file::AdiosFile)
    perform_puts!(file.engine)
    return nothing
end

export adios_all_variable_names
"""
    vars = adios_all_variable_names(file::AdiosFile)
    vars::Vector{String}

List (recursively) all variables in the file.
"""
function adios_all_variable_names(file::AdiosFile)
    vars = inquire_all_variables(file.io)
    vars ≡ nothing && return String[]
    return name.(vars)
end

export adios_group_variable_names
"""
    vars = adios_group_variable_names(file::AdiosFile, groupname::AbstractString)
    vars::Vector{String}

List (non-recursively) all variables in the group `groupname` in the
file.
"""
function adios_group_variable_names(file::AdiosFile, groupname::AbstractString)
    vars = inquire_group_variables(file.io, groupname)
    vars ≡ nothing && return String[]
    return name.(vars)
end

export IORef
"""
    mutable struct IORef{T,D}

A reference to the value of a variable that has been scheduled to be
read from disk. This value cannot be accessed bofre the read
operations have actually been executed.

Use `fetch(ioref::IORef)` to access the value. `fetch` will trigger
the actual reading from file if necessary. It is most efficient to
schedule multiple read operations at once.

Use `adios_perform_gets` to trigger reading all currently scheduled
variables.
"""
mutable struct IORef{T,D}
    engine::Union{Nothing,Engine}
    array::Array{T,D}
    function IORef{T,D}(engine::Engine, array::Array{T,D}) where {T,D}
        return new{T,D}(engine, array)
    end
end

"""
    isready(ioref::IORef)::Bool

Check whether an `IORef` has already been read from file.
"""
Base.isready(ioref::IORef) = ioref.engine ≡ nothing

"""
    value = fetch(ioref::IORef{T,D}) where {T,D}
    value::Array{T,D}

Access an `IORef`. If necessary, the variable is read from file and
then cached. (Each `IORef` is read at most once.)

Scalars are handled as zero-dimensional arrays. To access the value of
a zero-dimensional array, write `array[]` (i.e. use array indexing,
but without any indices).
"""
function Base.fetch(ioref::IORef)
    isready(ioref) || perform_gets(ioref.engine)
    @assert isready(ioref)
    # return 0-D arrays as scalars
    # D == 0 && return ioref.array[]
    return ioref.array
end

export adios_get
"""
    ioref = adios_get(file::AdiosFile, name::AbstractString)
    ioref::Union{Nothing,IORef}

Schedule reading a variable from a file.

The variable is not read until `adios_perform_gets` is called. This
happens automatically when the `IORef` is accessed (via `fetch`). It
is most efficient to first schedule multiple variables for reading,
and then executing the reads together.
"""
function adios_get(file::AdiosFile, name::AbstractString)
    var = inquire_variable(file.io, name)
    var ≡ nothing && return nothing
    T = type(var)
    T ≡ nothing && return nothing
    D = ndims(var)
    D ≡ nothing && return nothing
    sh = count(var)
    sh ≡ nothing && return nothing
    ioref = IORef{T,D}(file.engine, Array{T,D}(undef, Tuple(sh)))
    get(file.engine, var, ioref.array)
    push!(file.engine.get_tasks, () -> (ioref.engine = nothing))
    return ioref
end

export adios_perform_gets
"""
    adios_perform_gets(file::AdiosFile)

Execute all currently scheduled read opertions. This makes all pending
`IORef`s ready.
"""
function adios_perform_gets(file::AdiosFile)
    perform_gets(file.engine)
    return nothing
end


text:
### A test for the server. Verifies if all services are accessable.

using DemeAssemblies: addtooken, SystemConfig

using PeaceVote
setnamespace(@__MODULE__)

import DemeAssemblies
uuid = PeaceVote.uuid("DemeAssemblies")

maintainer = PeaceVote.Signer(uuid,"maintainer")

#SERVER_PORT =
SERVER_ID = 21908404092254970037667366726608484539958839406300700723035320376059498464496

config = SystemConfig(SERVER_ID)

#ca = PeaceVote.Signer(uuid,"ca")

### Getting certificates by proving identity
# certificates = []
# for i in 1:3
#     account = "account$i"
#     member = PeaceVote.KeyChain(uuid,account)
#     identification = PeaceVote.ID("$i","today",member.member.id)
#     cert = PeaceVote.Certificate(identification,ca)
#     push!(certificates,cert)
# end

# tookenrecords = []

# for

# tookens = Channel(Inf)

tooken = 223434

addtooken(config.bcconfig,tooken,maintainer)

# # Each member then participates.
# @sync for (cert,i) in zip(certificates,1:3)
#     @async begin
#         account = "account$i"
#         keychain = PeaceVote.KeyChain(uuid,account)
#         PeaceVote.register(uuid,cert)
#     end
# end

# sleep(1)

# @sync for i in 1:3
#     @async begin
#         account = "account$i"
#         keychain = PeaceVote.KeyChain(uuid,account)
#         PeaceVote.braid!(keychain)
#     end
# end

# pmember = PeaceVote.KeyChain(uuid,"account2")
# PeaceVote.propose("Found peace for a change?",["yes","no","maybe"],pmember);

# sleep(1)

# ledger = PeaceVote.Ledger(uuid)
# PeaceVote.sync!(ledger, uuid)
# braidchain = PeaceVote.braidchain(ledger,uuid)

# proposal = PeaceVote.proposals(braidchain)[1]

# @sync for i in 1:3
#     @async begin
#         account = "account$i"
#         keychain = PeaceVote.KeyChain(uuid,account)
#         option = PeaceVote.Option(proposal,rand(1:3))
#         PeaceVote.vote(option,keychain,braidchain)
#     end
# end

# sleep(1)

# # Need to add synchronization here

# ledger = PeaceVote.Ledger(uuid)
# PeaceVote.sync!(ledger, uuid)
# braidchain = PeaceVote.braidchain(ledger,uuid)

# members = PeaceVote.members(braidchain)

# for proposal in PeaceVote.proposals(braidchain)
#     tally = PeaceVote.count(uuid,proposal,braidchain)
#     @show proposal,tally
# end

max_stars_repo_name: akels/DemocracyCredit

max_issues_repo_name: akels/DemocracyCredit

max_forks_repo_name: akels/DemocracyCredit


text:
@inline function checked_sqrt(x::V) where {W, T, V <: Union{SIMDPirates.Vec{W,T},SIMDPirates.SVec{W,T}}}
    y = SIMDPirates.vsqrt(x)
    SIMDPirates.vifelse(SIMDPirates.visfinite(y), y, SIMDPirates.vbroadcast(V, T(0.005)))
end
@inline function checked_sqrt(x::T) where {T <: Number}
    y = @fastmath sqrt(x)
    isfinite(T) ? y : T(0.005)
end

@generated function process_big_prop_points!(X::ResizableMatrix{T}, BPP::AbstractMatrix{T}) where T
    quote
        resize!(X, size(BPP,1))
        @vectorize $T for i ∈ 1:size(BPP,1)
            y1 = BPP[i,1]
            y2 = BPP[i,2]
            y3 = BPP[i,3]

            R11 = one(T) / checked_sqrt(BPP[i,5])
            # R11 = SIMDPirates.vifelse(SIMDPirates.visfinite(R11)
            x1 = R11 * y1
            L21 = R11 * BPP[i,6]
            L31 = R11 * BPP[i,8]
            R22 = one(T) / checked_sqrt(BPP[i,7] - L21*L21)
            x2 = R22 * (y2 - L21*x1)

            L32 = R22 * (BPP[i,9] - L21 * L31)
            R33 = one(T) / checked_sqrt(BPP[i,10] - L31*L31 - L32*L32)

            X[i,1] = x1
            X[i,2] = x2
            X[i,3] = R33 * (y3 - L31*x1 - L32*x2)
        end
    end
end

@generated function process_big_prop_points!(X::AbstractMatrix{T}, BPP::AbstractMatrix{T}) where T
    quote
        @vectorize $T for i ∈ 1:size(BPP,1)
            y1 = BPP[i,1]
            y2 = BPP[i,2]
            y3 = BPP[i,3]

            R11 = one(T) / checked_sqrt(BPP[i,5])
            # R11 = SIMDPirates.vifelse(SIMDPirates.visfinite(R11)
            x1 = R11 * y1
            L21 = R11 * BPP[i,6]
            L31 = R11 * BPP[i,8]
            R22 = one(T) / checked_sqrt(BPP[i,7] - L21*L21)
            x2 = R22 * (y2 - L21*x1)

            L32 = R22 * (BPP[i,9] - L21 * L31)
            R33 = one(T) / checked_sqrt(BPP[i,10] - L31*L31 - L32*L32)

            X[i,1] = x1
            X[i,2] = x2
            X[i,3] = R33 * (y3 - L31*x1 - L32*x2)
        end
    end
end

function generate_rank1covariances!(rank1covs::AbstractVector{InverseWishart{T}}, X::AbstractMatrix{T}) where T
    resize!(rank1covs, size(X,1))
    @inbounds for i ∈ 1:size(X,1)
        rank1covs[i] = InverseWishart(X[i,1], X[i,2], X[i,3])
    end
end

function process_BPP!(X, rank1covs, mahals, BPP)
    process_big_prop_points!(X, BPP)
    generate_rank1covariances!(rank1covs, X)
    MahalanobisDistances!(mahals, X)
    nothing
end

function process_BPP!(X, rank1covs, mahals, icc::InvCholCovar, BPP, t)
    process_big_prop_points!(X, BPP)
    decorrelate_data!(icc, X, t)
    generate_rank1covariances!(rank1covs, X)
    MahalanobisDistances!(mahals, X)
    nothing
end


text:
#=
	Testing out the sundials package

=#

max_stars_repo_name: jepayne/NYU

max_issues_repo_name: jepayne/NYU

max_forks_repo_name: jepayne/NYU


text:
using ConvNet

CN = ConvNet
using Base.Test

function test_fc_forward()
    dimensions = (1, 1, 3)
    input = CN.Vol(reshape(float64([1.0, 2.0, 0.0]), dimensions), zeros(dimensions))
    layer = CN.FullyConnectedLayer(2, 1.0, dimensions)
    CN.flatten(layer.filters[1].w)[:] = [-5.0, 2.0, 3.0]
    CN.flatten(layer.filters[2].w)[:] = [3.0, 2.0, 5.0]
    CN.forward(layer, input)
    @test layer.inputActivation == input
    @test size(layer.outputActivation.w) == (1, 1, 2)
    @test CN.flatten(layer.outputActivation.w) == [0.0, 8.0]
end

function test_sig_forward()
    dimensions = (1, 1, 2)
    input = CN.Vol(reshape(float64([1.0, 5.0]), dimensions), zeros(dimensions))
    layer = CN.SigmoidLayer(dimensions)
    CN.forward(layer, input)
    @test layer.inputActivation == input
    @test size(layer.outputActivation.w) == dimensions
    @test CN.flatten(layer.outputActivation.w) == [0.7310585786300049, 0.9933071490757153]
end

function test_softmax_forward()
    dimensions = (1, 1, 5)
    input = CN.Vol(reshape(float64([1.0, 2.0, 3.0, 4.0, 5.0]), dimensions), zeros(dimensions))
    layer = CN.SoftMaxLayer(5)
    output = CN.forward(layer, input)
    @test output == layer.outputActivation
    @test layer.inputActivation == input
    @test size(layer.outputActivation.w) == dimensions
    @test layer.outputActivation.w |> CN.flatten |> x -> round(x, 2) == [0.01, 0.03, 0.09, 0.23, 0.64]
end

function test_softmax_backward()
    dimensions = (1, 1, 2)
    input = CN.Vol(reshape(float64([1.0, 5.0]), dimensions), zeros(dimensions))
    layer = CN.SoftMaxLayer(2)
    output = CN.forward(layer, input)

    @test layer.inputActivation == input
    @test CN.backward(layer, 2) < CN.backward(layer, 1)

    CN.backward(layer, 2)
    oldDw = copy(layer.inputActivation.dw)
    CN.backward(layer, 2)
    @test oldDw == layer.inputActivation.dw
    @test layer.inputActivation == input
    println(layer.inputActivation)
    @test CN.flatten(layer.dw) == [0.017986209962091555, -0.01798620996209155]
end

function test_softmax_backward()
    numOutputs = 2
    dimensions = (1, 1, numOutputs)

    input = CN.Vol(reshape(float64([1.0, 5.0]), dimensions), zeros(dimensions))
    layer = CN.SoftMaxLayer(numOutputs)
    output = CN.forward(layer, input)

    @test layer.inputActivation == input
    @test CN.backward(layer, 2) < CN.backward(layer, 1)

    CN.backward(layer, 2)
    oldDw = copy(layer.inputActivation.dw)
    CN.backward(layer, 2)
    @test oldDw == layer.inputActivation.dw
    @test layer.inputActivation == input
    @test CN.flatten(layer.inputActivation.dw) == [0.017986209962091555, -0.01798620996209155]
end

function test_xor()
    numInputs = 2
    numHidden = 2
    numClasses = 2
    net = CN.Net([CN.FullyConnectedLayer(numHidden, 1.0, (1, 1, numInputs)),
                  CN.FullyConnectedLayer(numClasses, 1.0, (1, 1, numHidden)),
                  CN.SoftMaxLayer(numClasses)])
    function instance()
        samples = rand(0:1, (1, 1, 2)) |> float64
        return CN.Vol(samples, zeros(size(samples))), (sum(samples) == 1.0 ? 1 : 2)
    end
    trainer = CN.Trainer(net, 0.01)
    iterations = 5
    for i in 1:iterations
        inputActivation, class = instance()
        @printf("Iteration %s\n\n", i)
        @printf("Instance: %s\nclass: %s\nprediction: %s\n\n",
                inputActivation, class, CN.predict(trainer.net, inputActivation))
        println("Before training")
        for (i, s) in enumerate(CN.gradientStates(trainer.net))
            @printf("Layer: %s, Gradients: %s\n", i, s)
        end

        CN.train(trainer, inputActivation, class)
        println("After training")
        for (i, s) in enumerate(CN.gradientStates(trainer.net))
            @printf("Layer: %s, Gradients: %s\n", i, s)
        end
        CN.updateGradients(trainer, 1)

        println("After BackPropagation")
        for (i, s) in enumerate(CN.gradientStates(trainer.net))
            @printf("Layer: %s, Gradients: %s\n", i, s)
        end
    end

    evaluations = 10
    evaluationSamples = [instance() for _ in 1:evaluations]
    errors = 0
    for _ in 1:evaluations
        inputActivation, class = instance()
        prediction = CN.predict(trainer, inputActivation)
        errors += int(prediction != class)
        println("Prediction: ", prediction)
        println("Class: ", class)
    end
    @printf("Accuracy: %s", errors / evaluations)
end

test_fc_forward()
test_sig_forward()
test_softmax_forward()
test_softmax_backward()
test_xor()


text:
@testset "bp on erdos_renyi" begin
    Random.seed!(16)
    N, c = 100, 4
    net = erdos_renyi(N, c / N, Network)
    eprop!(net, "J", e -> randn())
    vprop!(net, "H", v -> 1 + randn())
    res = Ising.run_bp(net, T=2, verbose=false)

    @test mean(res.mags) > 0.1
    # writenetwork("rrg.gml", net)
end

@testset "bp on rrg" begin
    ## NO COUPLINGS
    T = 3
    N, z = 100, 4
    Random.seed!(17)
    net = random_regular_graph(N, z, Network)
    eprop!(net, "J", e -> 0)
    vprop!(net, "H", v -> randn())
    res = Ising.run_bp(net, T=T, verbose=false)

    @test res.mags ≈ tanh.(vprop(net, "H").data ./ T)

    ## NO FIELDS, J=1
    T = 2
    N, z = 100, 4
    net = random_regular_graph(N, z, Network)
    eprop!(net, "J", e -> 1)
    vprop!(net, "H", v -> 0)
    res = Ising.run_bp(net, T=T, verbose=false)

    # X = rrg_to_mcgraph(net)
    # Es, σ, mc_mags = run_monte_carlo(X, β=1/T, infotime=10, sweeps=10^7);
    # Δ = mean(abs, res.mags .- mc_mags)
    # @test Δ <= 0.01
end

# 2D lattice

Tc = 2 / log(1+√2) #2.269185314213022 from Onsager's solution
T = 2
L = 16
net = Grid([L, L], Network, periodic=true)
eprop!(net, "J", e -> 1)
vprop!(net, "H", v -> 0)
res = Ising.run_bp(net, T=T, verbose=false)


text:
# In principle, FiniteCMPSTangent has the same fields as FiniteCMPS.
# However, instances of FiniteCMPSTangent constitute a vector space, and we can define the # necessary methods to allow using them as such in KrylovKit.jl.
# Furthermore, we explicitly store the base point and the indices of the original tent
# functions, as we might actually be working on a finer grid than the optimization grid.
mutable struct FiniteCMPSTangent{T<:MatrixFunction, N, V<:AbstractVector, I}
    base::FiniteCMPS{T, N, V}
    dQ::T
    dRs::NTuple{N,T}
    dvL::V
    dvR::V
    indices::I
    function FiniteCMPSTangent(Ψ::FiniteCMPS{T, N, V}, dQ::T, dRs::NTuple{N,T}, dvL::V, dvR::V, indices = 1:length(nodes(dQ))) where {T,N,V}
        grid = nodes(Ψ.Q)
        grid == nodes(dQ) || throw(DomainMismatch())
        for dR in dRs
            grid == nodes(dR) || throw(DomainMismatch())
        end
        a = first(grid)
        Qa = Ψ.Q(a)
        length(dvL) == length(dvR) == size(Qa,1) || throw(DimensionMismatch())
        size(dQ(a)) == size(Qa) || throw(DimensionMismatch())
        for dR in dRs
            size(dR(a)) == size(Qa) || throw(DimensionMismatch())
        end

        first(indices) == 1
        @assert last(indices) == length(grid)
        @assert issorted(indices)

        return new{T,N,V,typeof(indices)}(Ψ, dQ, dRs, dvL, dvR, indices)
    end
end
domain(Φ::FiniteCMPSTangent) = domain(Φ.dQ)
base(Φ::FiniteCMPSTangent) = Φ.base

Base.iterate(Φ::FiniteCMPSTangent, args...) = iterate((Φ.dQ, Φ.dRs, Φ.dvL, Φ.dvR), args...)

for f in (:copy, :zero, :similar)
    @eval Base.$f(Φ::FiniteCMPSTangent) =
        FiniteCMPSTangent(Φ.base, $f(Φ.dQ), $f.(Φ.dRs), $f(Φ.dvL), $f(Φ.dvR), Φ.indices)
end

function Base.copy!(Φ₁::FiniteCMPSTangent, Φ₂::FiniteCMPSTangent)
    (base(Φ₁) === base(Φ₂) && Φ₁.indices === Φ₂.indices) || throw(DomainMismatch())

    copy!(Φ₁.dvL, Φ₂.dvL)
    copy!(Φ₁.dvR, Φ₂.dvR)
    copy!(Φ₁.dQ, Φ₂.dQ)
    copy!.(Φ₁.dRs, Φ₂.dRs)
    return Φ₁
end

# Basic out-of-place arithmitic
Base.:-(Φ::FiniteCMPSTangent) =
    FiniteCMPSTangent(base(Φ), -Φ.dQ, .- Φ.dRs, -Φ.dvL, -Φ.dvR, Φ.indices)
Base.:*(α, Φ::FiniteCMPSTangent) =
        FiniteCMPSTangent(base(Φ), α * Φ.dQ, α .* Φ.dRs, α * Φ.dvL, α * Φ.dvR, Φ.indices)
Base.:\(α, Φ::FiniteCMPSTangent) =
        FiniteCMPSTangent(base(Φ), α \ Φ.dQ, α .\ Φ.dRs, α \ Φ.dvL, α \ Φ.dvR, Φ.indices)
Base.:*(Φ::FiniteCMPSTangent, α) =
        FiniteCMPSTangent(base(Φ), Φ.dQ * α, Φ.dRs .* α, Φ.dvL * α, Φ.dvR * α, Φ.indices)
Base.:/(Φ::FiniteCMPSTangent, α) =
        FiniteCMPSTangent(base(Φ), Φ.dQ / α, Φ.dRs ./ α, Φ.dvL / α, Φ.dvR / α, Φ.indices)
function Base.:+(Φ₁::FiniteCMPSTangent, Φ₂::FiniteCMPSTangent)
    (base(Φ₁) == base(Φ₂) && Φ₁.indices == Φ₂.indices) || throw(DomainMismatch())

    return FiniteCMPSTangent(base(Φ₁),
                                Φ₁.dQ + Φ₂.dQ,
                                Φ₁.dRs .+ Φ₂.dRs,
                                Φ₁.dvL + Φ₂.dvL,
                                Φ₁.dvR + Φ₂.dvR,
                                Φ₁.indices)
end
function Base.:-(Φ₁::FiniteCMPSTangent, Φ₂::FiniteCMPSTangent)
    (base(Φ₁) == base(Φ₂) && Φ₁.indices == Φ₂.indices) || throw(DomainMismatch())

    return FiniteCMPSTangent(base(Φ₁),
                                Φ₁.dQ - Φ₂.dQ,
                                Φ₁.dRs .- Φ₂.dRs,
                                Φ₁.dvL - Φ₂.dvL,
                                Φ₁.dvR - Φ₂.dvR,
                                Φ₁.indices)
end

# In-place arithmitic
function LinearAlgebra.axpy!(α, Φ₁::FiniteCMPSTangent, Φ₂::FiniteCMPSTangent)
    (base(Φ₁) == base(Φ₂) && Φ₁.indices == Φ₂.indices) || throw(DomainMismatch())

    axpy!(α, Φ₁.dvL, Φ₂.dvL)
    axpy!(α, Φ₁.dvR, Φ₂.dvR)
    axpy!(α, Φ₁.dQ, Φ₂.dQ)
    axpy!.(α, Φ₁.dRs, Φ₂.dRs)
    return Φ₂
end
function LinearAlgebra.axpby!(α, Φ₁::FiniteCMPSTangent, β, Φ₂::FiniteCMPSTangent)
    (base(Φ₁) == base(Φ₂) && Φ₁.indices == Φ₂.indices) || throw(DomainMismatch())

    axpby!(α, Φ₁.dvL, β, Φ₂.dvL)
    axpby!(α, Φ₁.dvR, β, Φ₂.dvR)
    axpby!(α, Φ₁.dQ, β, Φ₂.dQ)
    axpby!.(α, Φ₁.dRs, β, Φ₂.dRs)
    return Φ₂
end
function LinearAlgebra.lmul!(α, Φ::FiniteCMPSTangent)
    lmul!(α, Φ.dvL)
    lmul!(α, Φ.dvR)
    lmul!(α, Φ.dQ)
    lmul!.(α, Φ.dRs)
    return Φ
end
function LinearAlgebra.rmul!(Φ::FiniteCMPSTangent, α)
    rmul!(Φ.dvL, α)
    rmul!(Φ.dvR, α)
    rmul!(Φ.dQ, α)
    rmul!.(Φ.dRs, α)
    return Φ
end
function LinearAlgebra.mul!(Φ₁::FiniteCMPSTangent, α, Φ₂::FiniteCMPSTangent)
    (base(Φ₁) == base(Φ₂) && Φ₁.indices == Φ₂.indices) || throw(DomainMismatch())

    mul!(Φ₁.dvL, α, Φ₂.dvL)
    mul!(Φ₁.dvR, α, Φ₂.dvR)
    mul!(Φ₁.dQ, α, Φ₂.dQ)
    mul!.(Φ₁.dRs, α, Φ₂.dRs)
    return Φ₁
end
function LinearAlgebra.mul!(Φ₁::FiniteCMPSTangent, Φ₂::FiniteCMPSTangent, α)
    (base(Φ₁) == base(Φ₂) && Φ₁.indices == Φ₂.indices) || throw(DomainMismatch())

    mul!(Φ₁.dvL, Φ₂.dvL, α)
    mul!(Φ₁.dvR, Φ₂.dvR, α)
    mul!(Φ₁.dQ, Φ₂.dQ, α)
    mul!.(Φ₁.dRs, Φ₂.dRs, α)
    return Φ₁
end

# We choose `dot` to represent a standard Euclidean norm.
# The actual metric of the manifold will be implemented as a preconditioner.
function LinearAlgebra.dot(Φ₁::FiniteCMPSTangent, Φ₂::FiniteCMPSTangent)
    (base(Φ₁) == base(Φ₂) && Φ₁.indices == Φ₂.indices) || throw(DomainMismatch())

    ind = Φ₁.indices
    s = dot(Φ₁.dvL, Φ₂.dvL) + dot(Φ₁.dvR , Φ₂.dvR)
    s += dot(view(nodevalues(Φ₁.dQ), ind), view(nodevalues(Φ₂.dQ), ind))
    for (dR₁, dR₂) in zip(Φ₁.dRs, Φ₂.dRs)
        s += dot(view(nodevalues(dR₁), ind), view(nodevalues(dR₂), ind))
    end
    return s
end
function LinearAlgebra.norm(Φ::FiniteCMPSTangent)
    ind = Φ.indices
    s = hypot(norm(Φ.dvL), norm(Φ.dvR))
    s = hypot(s, norm(view(nodevalues(Φ.dQ), ind)))
    for dR in Φ.dRs
        s = hypot(s, norm(view(nodevalues(dR), ind)))
    end
    return s
end

# Given the functional derivatives of some object with respect to Q, R and ∂R as instances
# of `AbstractPiecewise`, compute the corresponding `PiecewiseLinear` version resulting from
# applying the chain rule
function _project(𝒬̅, ℛ̅s, ∂ℛ̅s = nothing; gradindices = 1:length(nodes(𝒬̅)))
    (a, b) = domain(𝒬̅)
    grid = collect(nodes(𝒬̅))

    # Compute gradients with respect to PiecewiseLinear parameters
    Q̄ = [zero(𝒬̅(a)) for _ = 1:length(gradindices)]
    R̄s = map(ℛ̅->[zero(ℛ̅(a)) for _ = 1:length(gradindices)], ℛ̅s)

    k = gradindices[1] # == 1
    knext = gradindices[2]
    xc = grid[k]
    xb = grid[knext]
    t = TaylorSeries([1,-1/(xb-xc)], xc)
    Q̄i = Q̄[1]
    for l = k:(knext-1)
        t = shift!(t, offset(𝒬̅[l]))
        Q̄i .+= integrate(𝒬̅[l] * t, (grid[l], grid[l+1]))
    end
    Q̄[1] = Q̄i

    for i = 2:length(gradindices)-1
        k = gradindices[i]
        kprev = gradindices[i-1]
        knext = gradindices[i+1]
        xa = grid[kprev]
        xc = grid[k]
        xb = grid[knext]

        Q̄i = Q̄[i]
        R̄is = getindex.(R̄s, i)
        t = TaylorSeries([0,1/(xc-xa)], xa)
        for l = kprev:(k-1)
            t = shift!(t, offset(𝒬̅[l]))
            Q̄i .+= integrate(𝒬̅[l] * t, (grid[l], grid[l+1]))
            for (R̄i, ℛ̅) in zip(R̄is, ℛ̅s)
                R̄i .+= integrate(ℛ̅[l] * t, (grid[l], grid[l+1]))
            end
            if !isnothing(∂ℛ̅s)
                for (R̄i, ∂ℛ̅) in zip(R̄is, ∂ℛ̅s)
                    R̄i .+= integrate(∂ℛ̅[l]/(xc-xa), (grid[l], grid[l+1]))
                end
            end
        end
        t = TaylorSeries([1,-1/(xb-xc)], xc)
        for l = k:(knext-1)
            t = shift!(t, offset(𝒬̅[l]))
            Q̄i .+= integrate(𝒬̅[l] * t, (grid[l], grid[l+1]))
            for (R̄i, ℛ̅) in zip(R̄is, ℛ̅s)
                R̄i .+= integrate(ℛ̅[l] * t, (grid[l], grid[l+1]))
            end
            if !isnothing(∂ℛ̅s)
                for (R̄i, ∂ℛ̅) in zip(R̄is, ∂ℛ̅s)
                    R̄i .+= integrate(∂ℛ̅[l]/(xc-xb), (grid[l], grid[l+1]))
                end
            end
        end
        setindex!(Q̄, Q̄i, i)
        setindex!.(R̄s, R̄is, i)
    end

    k = gradindices[end]
    kprev = gradindices[end-1]
    xa = grid[kprev]
    xc = grid[k]

    Q̄i = Q̄[end]
    t = TaylorSeries([0,1/(xc-xa)], xa)
    for l = kprev:(k-1)
        t = shift!(t, offset(𝒬̅[l]))
        Q̄i .+= integrate(𝒬̅[l] * t, (grid[l], grid[l+1]))
    end
    Q̄[end] = Q̄i

    if gradindices == 1:length(grid)
        ∇Q = PiecewiseLinear(grid, Q̄)
        ∇Rs = PiecewiseLinear.((grid,), R̄s)
    else
        grid2 = grid[gradindices]
        ∇Q = PiecewiseLinear(grid2, Q̄)
        ∇Rs = PiecewiseLinear.((grid2,), R̄s)
        ∇Q = PiecewiseLinear(grid, ∇Q.(grid))
        ∇Rs = map(∇Rs) do ∇R
            PiecewiseLinear(grid, ∇R.(grid))
        end
    end

    return ∇Q, ∇Rs
end

# Actual metric acting on a given tangent vector
function metric(Φ::FiniteCMPSTangent, Ψρs::FiniteCMPSData;
                δ = 0, Kmax = 50, tol = defaulttol(base(Φ)),
                left_boundary = :free, right_boundary = :free,
                gradindices = 1:length(nodes(Φ.dQ)))

    Ψ, ρL, ρR = Ψρs

    base(Φ) == Ψ || throw(DomainMismatch())

    Q, Rs, vL, vR = Ψ
    (a, b) = domain(Ψ)
    T = scalartype(Ψ)

    fL = ρL * Φ.dQ
    fR = Φ.dQ * ρR
    temp = zero(ρL)
    for (R, dR) in zip(Rs, Φ.dRs)
        RdρL = mul!(temp, R', ρL)
        fL = mul!(fL, RdρL, dR, one(T), one(T))
        ρRRd = mul!(temp, ρR, R')
        fR = mul!(fR, dR, ρRRd, one(T), one(T))
    end
    FL, = lefttransfer(zero(fL(a)), fL, Ψ; Kmax = Kmax, tol = tol)
    FR, = righttransfer(zero(fR(b)), fR, Ψ; Kmax = Kmax, tol = tol)

    𝒬̅ = zero(FL)
    𝒬̅ = mul!(𝒬̅, FL, ρR, one(T), one(T))
    𝒬̅ = mul!(𝒬̅, ρL, FR, one(T), one(T))
    ℛ̅s = map(_->zero(FL), Rs)
    for (R, ℛ̅, dR) in zip(Rs, ℛ̅s, Φ.dRs)
        FLR = mul!(temp, FL, R)
        ℛ̅ = mul!(ℛ̅, FLR, ρR, one(T), one(T))
        RFR = mul!(temp, R, FR)
        ℛ̅ = mul!(ℛ̅, ρL, RFR, one(T), one(T))
        dRρR = mul!(temp, dR, ρR)
        ℛ̅ = mul!(ℛ̅, ρL, dRρR, one(T), one(T))
    end
    return _project(𝒬̅, ℛ̅s; gradindices = gradindices)
end


text:
module LibCURL

typealias time_t Int
typealias size_t Csize_t
typealias curl_off_t Int64

include("lC_exports_h.jl")
include("lC_common_h.jl")

libcurl = :libcurl

include("lC_curl_h.jl")

@c CURLcode curl_easy_setopt (Ptr{CURL}, CURLoption, Int...) libcurl
@c CURLcode curl_easy_setopt (Ptr{CURL}, CURLoption,Ptr{Void}...) libcurl
@c CURLcode curl_easy_setopt (Ptr{CURL}, CURLoption,Ptr{Uint8}...) libcurl

@c CURLcode curl_multi_setopt (Ptr{CURL}, CURLoption, Int...) libcurl
@c CURLcode curl_multi_setopt (Ptr{CURL}, CURLoption,Ptr{Void}...) libcurl
@c CURLcode curl_multi_setopt (Ptr{CURL}, CURLoption,Ptr{Uint8}...) libcurl

@c CURLcode curl_easy_getinfo (Ptr{CURL}, CURLoption,Ptr{Int}...) libcurl
@c CURLcode curl_easy_getinfo (Ptr{CURL}, CURLoption,Ptr{Float64}...) libcurl
@c CURLcode curl_easy_getinfo (Ptr{CURL}, CURLoption,Ptr{Uint8}...) libcurl

include("lC_defines_h.jl")

include("Mime_ext.jl")
export Mime_ext

end


text:
# type annotations and scope:
x = 1.0 # x is Float64
x = 1 # now x is Int
# y::Float64 = 1.0 # LoadError: "y is not defined"

function scopetst()
    y::Float64 = 1.0 # y must be Float64, not possible in global scope
    # y = 1
end

# compound expressions:
x = begin
  a = 5
  2 * a
end # after: x = 10
println(a) #=> a = 5
x = (a = 5; 2 * a)

const GC = 6.67e-11 # gravitational constant in m3/kg s2

GC = 3.14 #=> Warning: redefining constant GC
# GC = 10 #=> ERROR: invalid redefinition of constant GC

julia> const ARR = [4,7,1]
julia> ARR[1] = 9
julia> show(ARR) #=> [9,7,1]
julia> ARR = [1, 2, 3] #=> Warning: redefining constant ARR


text:
"
    _ftIdx(ft, grp)

Finds the index of grp (string) in ft.labels which corresponds to ft.mat.
For Index i for the grp  returned by this function ft[i, :, :] returns the
2D array [[TP, FP], [FN, TN]] for that group.
"
function _ftIdx(ft::FairTensor, grp)
    idx = findfirst(x->x==string.(grp), ft.labels)
    if idx==nothing throw(ArgumentError("$grp not found in the fairness tensor")) end
    return idx
end

# Helper function for calculating TruePositive, FalseNegative, etc.
# If grp is :, it calculates combined value for all groups, else the specified group
_calcmetric(ft::FairTensor, grp, inds...) = typeof(grp)==Colon ? sum(ft[:, inds...]) : sum(ft[_ftIdx(ft, grp), inds...])

(::TruePositive)(ft::FairTensor; grp=:) = _calcmetric(ft, grp, 1, 1)
(::FalsePositive)(ft::FairTensor; grp=:) = _calcmetric(ft, grp, 1, 2)
(::FalseNegative)(ft::FairTensor; grp=:) = _calcmetric(ft, grp, 2, 1)
(::TrueNegative)(ft::FairTensor; grp=:) = _calcmetric(ft, grp, 2, 2)

# The functions true_positive, false_negative, etc are instances of TruePositive, FalseNegative, etc.
# So on using true_positive(fair_tensor) will use the above defined functions
(::TPR)(ft::FairTensor; grp=:) = 1/(1+false_negative(ft; grp=grp)/(1e-15 + true_positive(ft; grp=grp)))
(::TNR)(ft::FairTensor; grp=:) = 1/(1+false_positive(ft; grp=grp)/(1e-15 + true_negative(ft; grp=grp)))
(::FPR)(ft::FairTensor; grp=:) = 1-true_negative_rate(ft; grp=grp)
(::FNR)(ft::FairTensor; grp=:) = 1-true_positive_rate(ft; grp=grp)

(::FDR)(ft::FairTensor; grp=:) = 1/(1+false_positive(ft; grp=grp)/(1e-15 + true_positive(ft; grp=grp)))
(::Precision)(ft::FairTensor; grp=:) = 1 - false_discovery_rate(ft; grp=grp)
(::NPV)(ft::FairTensor; grp=:) =  1/(1+false_negative(ft; grp=grp)/(1e-15 + true_negative(ft; grp=grp)))

(::MLJBase.Accuracy)(ft::FairTensor; grp=:) = (true_positive(ft; grp = grp) + true_negative(ft; grp = grp))/(
        true_positive(ft; grp = grp) + true_negative(ft; grp = grp) +
        false_positive(ft; grp = grp) + false_negative(ft; grp = grp))

struct PredictedPositiveRate <: Measure end
PP = PredictedPositiveRate
predicted_positive_rate = PP()
ppr = predicted_positive_rate
(::PP)(ft::FairTensor; grp=:) = _calcmetric(ft, grp, 1, :)/_calcmetric(ft, grp, :, :)

struct FalseOmissionRate <: Measure end
FOR = FalseOmissionRate
false_omission_rate = FOR()
foar = false_omission_rate
(::FOR)(ft::FairTensor; grp=:) = 1/(1+true_negative(ft; grp=grp)/(1e-15 + false_negative(ft; grp=grp)))

struct TruePositiveRateDifference <: Measure end
TPRD = TruePositiveRateDifference
true_positive_rate_difference = TPRD()
tprd = true_positive_rate_difference
(::TPRD)(ft::FairTensor; grp1=:,grp2=:)=abs(true_positive_rate(ft;grp=grp1)-true_positive_rate(ft;grp=grp2))

struct FalsePositiveRateDifference <: Measure end
FPRD = FalsePositiveRateDifference
false_positive_rate_difference = FPRD()
fprd = false_positive_rate_difference
(::FPRD)(ft::FairTensor; grp1=:,grp2=:)=abs(false_positive_rate(ft;grp=grp1)-false_positive_rate(ft;grp=grp2))

struct FalseNegativeRateDifference <: Measure end
FNRD = FalseNegativeRateDifference
false_negative_rate_difference = FNRD()
fnrd = false_negative_rate_difference
(::FNRD)(ft::FairTensor; grp1=:, grp2=:)=abs(false_negative_rate(ft;grp=grp1)-false_negative_rate(ft;grp=grp2))

struct FalseOmissionRateDifference <: Measure end
FORD = FalseOmissionRateDifference
false_omission_rate_difference = FORD()
ford = false_omission_rate_difference
(::FORD)(ft::FairTensor; grp1=:, grp2=:)=abs(false_omission_rate(ft;grp=grp1)-false_omission_rate(ft;grp=grp2))

struct FalseDiscoveryRateDifference <: Measure end
FDRD = FalseDiscoveryRateDifference
false_discovery_rate_difference = FDRD()
fdrd = false_discovery_rate_difference
(::FDRD)(ft::FairTensor; grp1=:, grp2=:)=abs(false_discovery_rate(ft;grp=grp1)-false_discovery_rate(ft;grp=grp2))

struct FalsePositiveRateRatio <: Measure end
FPRR = FalsePositiveRateRatio
false_positive_rate_ratio = FPRR()
fprr = false_positive_rate_ratio
(::FPRR)(ft::FairTensor; gpr1=:, grp2=:)=false_positve_rate(ft;grp=grp1)/false_positive_rate(ft;grp=grp2)

struct FalseNegativeRateRatio <: Measure end
FNRR = FalseNegativeRateRatio
false_negative_rate_ratio = FNRR()
fnrr = false_negative_rate_ratio
(::FNRR)(ft::FairTensor; gpr1=:, grp2=:)=false_negative_rate(ft;grp=grp1)/false_negative_rate(ft;grp=grp2)

struct FalseOmissionRateRatio <: Measure end
FORR = FalseOmissionRateRatio
false_omission_rate_ratio = FORR()
forr = false_omission_rate_ratio
(::FORR)(ft::FairTensor; grp1=:, grp2=:)=false_omission_rate(ft;grp=grp1)/false_omission_rate(ft;grp=grp2)

struct FalseDiscoveryRateRatio <: Measure end
FDRR = FalseDiscoveryRateRatio
false_discovery_rate_ratio = FDRR()
fdrr = false_discovery_rate_ratio
(::FDRR)(ft::FairTensor; grp1=:, grp2=:)=false_discovery_rate(ft;grp=grp1)/false_discovery_rate(ft;grp=grp2)

struct AverageOddsDifference <: Measure end
AOD = AverageOddsDifference
average_odds_difference = AOD()
aod = average_odds_difference
(::AOD)(ft::FairTensor; grp1=:, grp2=:)=0.5*(false_positive_rate_difference(ft;grp1=grp1,grp2=grp2)+true_positive_rate_difference(ft;grp1=grp1,grp2=grp2))