text:
# JuliaCon 2019 Demonstration
if (isdir("Z:"))
    include("Z:\\Ultimate Datetime Datatype\\Ultimate Datetime Git\\Ultimate-Datetime\\Julia Integration\\UltimateDT.jl")
elseif (isdir("C:"))
    include("C:\\Users\\jay\\Documents\\Ultimate Datetime Datatype\\Ultimate Datetime\\Julia Integration\\UltimateDT.jl")
end
# Create a struct to contain code fragments
struct CodeFragment
    description::String
    expression::Expr
end
# Create a function to describe then demonstrate a code fragment
function demo_code(codeFragment::CodeFragment)
    if (codeFragment.description != "")
        printstyled(stdout, "\n *** ", codeFragment.description, " ***\n", bold=true, color=:light_green)
    end
    # The replace functions removes comment lines inserted by the system
    print(stdout, replace(string(codeFragment.expression), r"#.*#\n" => ""), "  ")
    if eval(codeFragment.expression) != nothing
        printstyled(stdout, eval(codeFragment.expression), color=:light_cyan)
    end
    readline() # Wait for a carraige return
    nothing
end
# Use the Unicode plotting package
using UnicodePlots
# Create a function to create the vector of vectors required to generate a boxplot
#   comparing two UncertainFloat64s.
function boxCompareUF(uf1::UncertainFloat64, uf2::UncertainFloat64)
    val1 = uf1.value - uf1.precision*uf1.uncertainty
    val2 = uf1.value - uf1.precision*uf1.uncertainty
    val3 = uf1.value
    val4 = uf1.value + uf1.precision*uf1.uncertainty
    val5 = uf1.value + uf1.precision*uf1.uncertainty
    val6 = uf2.value - uf2.precision*uf2.uncertainty
    val7 = uf2.value - uf2.precision*uf2.uncertainty
    val8 = uf2.value
    val9 = uf2.value + uf2.precision*uf2.uncertainty
    val10 = uf2.value + uf2.precision*uf2.uncertainty
    return [[val1,val2,val3,val4,val5],[val6,val7,val8,val9,val10]]
end
# Create an array with strings for all of the comparison operators
compops = ["==","!=","==′","!=′","==″","!=″","==‴","!=‴",">","<=″",
           ">′","<=′",">″","<=",">=","<″",">=′","<′",">=″","<"]
# Create an array with descriptions of the operators and print them in boxAndCompare
compopsDesc = ["Must equal", "May not equal", "Has same midpoint as", "Does not have same midpoint as",
               "May equal", "Cannot equal", "Has same range as", "Does not have same range as",
               "Must be greater than", "Can be less than or equal to", "Has larger midpoint than",
               "Has midpoint less than or equal to that of", "Can be greater than", "Must be less than or equal to",
               "Must be greater than or equal to", "Can be less than", "Has midpoint greater than or equal to that of",
               "Has smaller midpoint than", "Can be greater than or equal to", "Must be less than"
               ]
# Create a function to draw a box plot and display the results of the comparisons
# uf1 and uf2 need to be in the global scope so that the eval in boxAndCompare receives the proper values
uf1 = UncertainFloat64(1.,p=1,u=1)
uf2 = uf1
function boxAndCompare(uf1t::UncertainFloat64, uf2t::UncertainFloat64)
    # Need to pass the global variables to eval
    global uf1 = uf1t
    global uf2 = uf2t
    # When used in an expression to be evaluated, boxplot works if it is the last line in
    #  the expression.  Otherwise, must use show(boxplot).
    show(boxplot([string(uf1), string(uf2)], boxCompareUF(uf1,uf2), border=:none))
    print(stdout,"\n")
    for i = 1:20
        code = Meta.parse("uf1" * compops[i] * "uf2")
        print(stdout, "\n", rpad(string(code),11), "   ")
        printstyled(stdout, rpad(compopsDesc[i],45), "  ", color=:light_green)
        if (eval(code) == true)
            printstyled(stdout, eval(code), color=:light_cyan)
        else
            printstyled(stdout, eval(code), color=:light_red)
        end
    end
    return nothing
end
# Create an array to contain the code fragments
code = Array{CodeFragment,1}()
# Uncertain Floats
# Create an Uncertain Float with default precision and uncertainty
push!(code,CodeFragment("Uncertain float with default precision and uncertainty",:(UncertainFloat64(1005.25))))
# Show change in default precision and uncertainty with magnitude
push!(code,CodeFragment("Default precision and uncertainty change with magnitude",:(UncertainFloat64(5005.25))))
# Create 2 uncertain floats to demonstrate arithmetic functions
push!(code,CodeFragment("Create uf1 with specified precision and uncertainty",:(uf1 = UncertainFloat64(10,p=1,u=2))))
push!(code,CodeFragment("Create uf2 with specified precision and uncertainty",:(uf2 = UncertainFloat64(8,p=1,u=1))))
push!(code,CodeFragment("Addition - add absolute uncertainties, use minimum absolute precision",:(uf1 + uf2)))
push!(code,CodeFragment("Subtraction - add absolute uncertainties, use minimum absolute precision",:(uf1 - uf2)))
push!(code,CodeFragment("Create uf3 with specified precision and uncertainty",:(uf3 = UncertainFloat64(100,p=1,u=2))))
push!(code,CodeFragment("Multiplication - add fractional uncertainties, use minimum fractional precision",:(uf3 * uf2)))
push!(code,CodeFragment("Division - add fractional uncertainties, use minimum fractional precision",:(uf3 / uf2)))
# Exact values have a precision of 1 and uncertainy of 0 (i.e., as Integers)
push!(code,CodeFragment("Create uf1 with specified precision and uncertainty",:(uf1 = UncertainFloat64(1005.25,p=.01,u=2))))
push!(code,CodeFragment("Create uf2 as an exact value",:(uf2 = UncertainFloat64(4,p=1,u=0))))
push!(code,CodeFragment("Multiply the 2 values.  Uf1 is only contributor to p and u.",:(uf1 * uf2)))
# Compare with integber result
push!(code,CodeFragment("Multiply an uncertain value by an integer.  Uf1 is only contributor to p and u.",:(4 * uf1)))
# Zero values require special handling
push!(code,CodeFragment("Zero value with an uncertainty of 1.",:(uf1 = UncertainFloat64(0,p=1,u=1))))
push!(code,CodeFragment("Zero value with an uncertainty of 1 billion",:(uf2 = UncertainFloat64(0,p=1e9,u=1))))
push!(code,CodeFragment("Non-zero value with a 90% uncertainty",:(uf3 = UncertainFloat64(10.,p=1,u=9))))
push!(code,CodeFragment("Multiplication multiplies the uncertainty by the largest possible value of the other factor.",:(uf3 * uf1)))
push!(code,CodeFragment("Division increases the uncertainty by (1 + fractional uncertainty of the denominator).",:(uf2 / uf3)))
# Demonstrate range comparisons
push!(code,CodeFragment("Range based comparison operators with overlapping values",
    :(boxAndCompare(UncertainFloat64(8.,p=1,u=2),UncertainFloat64(9.,p=1,u=2)))))
push!(code,CodeFragment("Range based comparison operators with exact values",
    :(boxAndCompare(UncertainFloat64(8,p=1,u=0),UncertainFloat64(8,p=1,u=0)))))
push!(code,CodeFragment("Range based comparison operators with disjoint values",
    :(boxAndCompare(UncertainFloat64(5,p=1,u=2),UncertainFloat64(10,p=1,u=2)))))
push!(code,CodeFragment("Range based comparison operators with a single overlapping value",
    :(boxAndCompare(UncertainFloat64(5,p=1,u=2),UncertainFloat64(9,p=1,u=2)))))
push!(code,CodeFragment("Range based comparison operators with the same range",
    :(boxAndCompare(UncertainFloat64(9,p=1,u=2),UncertainFloat64(9,p=1,u=2)))))
# Demonstrate that relative time is different from UTCDatetime
push!(code,CodeFragment("Create a UTCDatetime with default precision and uncertainty.",
    :(utc1 = UTCDatetime(y=2019,m=7,d=23,h=11,min=45))))
push!(code,CodeFragment("Create a 2nd UTCDatetime with default precision and uncertainty.",
    :(utc2 = UTCDatetime(y=2019,m=7,d=23,h=10,min=45))))
push!(code,CodeFragment("Subtract the 2 UTCDatetimes to produce a relative datetime.",:(rel1 = utc1 - utc2)))
push!(code,CodeFragment("Relative datetimes can be positive or negative.",:(rel2 = utc2 - utc1)))
# Show allowed operations on relative datetimes
push!(code,CodeFragment("Create a UTCDatetime with nanosecond precision and uncertainty.",
    :(utc1 = UTCDatetime(y=2019,m=7,d=23,h=11,min=45,p=-9,u=1))))
push!(code,CodeFragment("Create a 2nd UTCDatetime with nanosecond precision and uncertainty.",
    :(utc2 = UTCDatetime(y=2019,m=7,d=23,h=10,min=45,p=-9,u=1))))
push!(code,CodeFragment("Subtract the 2 UTCDatetimes to produce a relative datetime.",:(rel1 = utc1 - utc2)))
push!(code,CodeFragment("Relative datetimes can be multiplied by integers and floating point numbers.",:(2*rel1 + 1.5*rel1)))
push!(code,CodeFragment("A relative datetime can be divided by another relative datetime.",:( x = rel1 / (2*rel1 + 1.5*rel1))))
push!(code,CodeFragment("The result of this division is an UncertianFloat64 (which is the original reason the type was created.)!", :(typeof(x))))
push!(code,CodeFragment("When the UncertainFloat64 is multiplied by the denominator, the precision is restored with the propogated uncertainty.",
    :(x * (2*rel1 + 1.5*rel1))))
# Demonstrate Julian to Gregorian transition
push!(code,CodeFragment("Create a UTCDatetime just before the earliest transtion to the Gregorian calendar.",
    :(utcj = UTCDatetime(y=1582,m=10,d=4,h=23,min=0,p=0,u=1))))
push!(code,CodeFragment("Attempt to create a UTCDatetime during the skipped period.",
    :(utcs = UTCDatetime(y=1582,m=10,d=5,h=1,min=0,p=0,u=1))))
push!(code,CodeFragment("Create a UTCDatetime just after the earliest transition to the Gregorian calendar.",
    :(utcg = UTCDatetime(y=1582,m=10,d=15,h=1,min=0,p=0,u=1))))
push!(code,CodeFragment("Very little time elapsed between those two dates!",
    :(utcg - utcj)))
# Demonstrate leap seconds
push!(code,CodeFragment("Leap seconds occur at midnight UTC, either on June 30 or December 31.",
    :(lsminus1 = UTCDatetime(y=2015,m=6,d=30,h=23,min=59,s=59))))
push!(code,CodeFragment("This is the actual leap second.  Notice the seconds field.",
    :(ls = UTCDatetime(y=2015,m=6,d=30,h=23,min=59,s=60))))
push!(code,CodeFragment("This is the next second.",
    :(lsplus1 = UTCDatetime(y=2015,m=7,d=1,h=0,min=0,s=0))))
push!(code,CodeFragment("The arithmetic works.",
    :(lsplus1 - lsminus1)))
push!(code,CodeFragment("Only minutes containing leap seconds have a 60th second.",
    :(ls = UTCDatetime(y=2015,m=6,d=30,h=23,min=58,s=60))))
# Move to LocalCalCoords
push!(code,CodeFragment("Construct a Gregorian date in the proleptic America/New_York time zone.",
    :(lccj = LocalCalCoordsDT(10, 15, 1582, "America/New_York", 2; cal=0, h=10))))
push!(code,CodeFragment("The Julian calendar continued on - the last country to adopt the Gregorian calendar was Greece in 1923.",
    :(lccg = LocalCalCoordsDT(10, 5, 1582, "America/New_York", 2; cal=1, h=10))))
push!(code,CodeFragment("For every Gregorian datetime, there is a corresponding Julian datetime.",
    :(lccg == lccj)))
# UTC Datetime Formatting
push!(code,CodeFragment("Formatting is automatic based on precision.",
    :(for i = -18:15;
        utc = UTCDatetime(gy=99,y=999999999,m=12,d=31,h=23,min=59,s=59,ns=999999999,as=999999999,p=i,u=1);
        printstyled(stdout,"\n", utc,color=:light_cyan);
      end)))
map(demo_code, code)
nothing

max_stars_repo_path: Julia Integration/JuliaCon 2019 Demonstration - Saved.jl

max_stars_repo_name: JayDweck/Ultimate-Datetime

max_issues_repo_name: JayDweck/Ultimate-Datetime

max_forks_repo_name: JayDweck/Ultimate-Datetime


text:
# utility/z-list.jl

include("downsample.jl")
include("eql_root.jl")
include("interp1.jl")
include("jinc.jl")
include("map_many.jl")
include("mask.jl")
include("max_percent_diff.jl")
include("ndgrid.jl")
include("rect.jl")
include("reverser.jl")
include("rmsd100.jl")

ext: jl


text:
#=============================#
#== put needed modules here ==#
#=============================#
import JuliaChem
import Statistics
#import HypothesisTests
import MPI
using BenchmarkTools

#================================#
#== JuliaChem execution script ==#
#================================#
function script(input_file)
  #== initialize JuliaChem runtime ==#
  JuliaChem.initialize()

  try
    #== read in input file ==#
    input_time1_t1 = time_ns()/1e9
    molecule, driver, model, keywords = JuliaChem.JCInput.run(input_file;
      output=0)
    input_time1_t2 = time_ns()/1e9
    input_time1 = input_time1_t2 - input_time1_t1

    input_time2_t1 = time_ns()/1e9
    molecule, driver, model, keywords = JuliaChem.JCInput.run(input_file;
      output=0)
    input_time2_t2 = time_ns()/1e9
    input_time2 = input_time2_t2 - input_time2_t1

    input_jit = input_time1 - input_time2
    molecule, driver, model, keywords = JuliaChem.JCInput.run(input_file;
      output=0)

    #== generate basis set ==#
    basis_time1_t1 = time_ns()/1e9
    basis = JuliaChem.JCBasis.run(molecule, model; output=0)
    basis_time1_t2 = time_ns()/1e9
    basis_time1 = basis_time1_t2 - basis_time1_t1

    basis_time2_t1 = time_ns()/1e9
    basis = JuliaChem.JCBasis.run(molecule, model; output=0)
    basis_time2_t2 = time_ns()/1e9
    basis_time2 = basis_time2_t2 - basis_time2_t1

    basis_jit = basis_time1 - basis_time2
    mol, basis = JuliaChem.JCBasis.run(molecule, model; output=0)

    #== perform scf benchmark ==#
    timeof = Vector{Float64}(undef,0)
    scf_time1 = 0.0
    if (driver == "energy")
      if (model["method"] == "RHF")
        #scf_time1_t1 = time_ns()/1e9
        #scf = JuliaChem.JCRHF.run(mol, basis, keywords["scf"];
        #  output=2)
        #scf_time1_t2 = time_ns()/1e9
        #scf_time1 = scf_time1_t2 - scf_time1_t1

        molecule, driver, model, keywords = JuliaChem.JCInput.run(input_file;
          output=0)
        mol, basis = JuliaChem.JCBasis.run(molecule, model; output=0)

        scf_timeof_t1 = time_ns()/1e9
        #scf = JuliaChem.JCRHF.Energy.run(mol, basis, keywords["scf"];
        #  output=2)
        scf = BenchmarkTools.@benchmark begin
          JuliaChem.JCRHF.Energy.run($mol, $basis, $(keywords["scf"]);
            output=2) #initial run
        end
        display(scf)
        scf_timeof_t2 = time_ns()/1e9
        push!(timeof, scf_timeof_t2 - scf_timeof_t1)
      end
    end
    scf_jit = scf_time1 - Statistics.mean(timeof)

    MPI.Barrier(MPI.COMM_WORLD)
    if (MPI.Comm_rank(MPI.COMM_WORLD) == 0) display(timeof) end

    #if (MPI.Comm_rank(MPI.COMM_WORLD) == 0)
      #== output relevant information ==#
    #  println("Input JIT: ", input_jit)
    #  println("Basis JIT: ", basis_jit)
    #  println("SCF JIT: ", scf_jit)
    #  println("Total JIT: ", input_jit + basis_jit + scf_jit)
    #  println("")

      #== perform t-test to compare to GAMESS ==#
    #  p = HypothesisTests.OneSampleTTest(timeof,2.94+0.19)
    #  println(p)
    #end
    MPI.Barrier(MPI.COMM_WORLD)

    #== finalize JuliaChem runtime ==#
    JuliaChem.finalize()

  catch e
    bt = catch_backtrace()
    msg = sprint(showerror, e, bt)
    println(msg)

    JuliaChem.finalize()
    exit()
  end
end

script(ARGS[1])


text:
using Mimi, CSVFiles, DataFrames, Query, Interpolations

@defcomp SSPs begin

    model   = Parameter{String}() # can be one of IIASA GDP, OECD Env-Growth, and PIK GDP_32
    ssp     = Parameter{String}() # can be one of SSP1, SSP2, SSP3, SSP5
    rcp     = Parameter{String}() # can be one of RCP1.9, RCP2.6, RCP3.7, RCP4.5, or RCP8.5

    country_names = Parameter{String}(index=[countries]) # need the names of the countries from the dimension

    # TODO double check units on gases, do we want any other gases or parameters?
    population      = Variable(index=[time, countries], unit="million")
    gdp             = Variable(index=[time, countries], unit="billion US\$2005/yr")

    co2_emissions   = Variable(index=[time], unit="GtC")
    ch4_emissions   = Variable(index=[time], unit="MtCH4")
    n2o_emissions   = Variable(index=[time], unit="MtN")
    sf6_emissions   = Variable(index=[time], unit="MtSF6")

    function init(p,v,d)

        # ----------------------------------------------------------------------
        # Checks

        model_options = ["IIASA GDP", "OECD Env-Growth", "PIK GDP_23"]
        !(p.model in model_options) && error("Model $(p.model) provided to SSPs component model parameter not found in available list: $(model_options)")

        ssp_options = ["SSP1", "SSP2", "SSP3", "SSP5"]
        !(p.ssp in ssp_options) && error("SSP $(p.ssp) provided to SSPs component ssp parameter not found in available list: $(ssp_options)")

        rcp_options = ["RCP1.9", "RCP2.6", "RCP3.7", "RCP4.5", "RCP8.5"]
        !(p.rcp in rcp_options) && error("RCP $(p.rcp) provided to SSPs component rcp parameter not found in available list: $(rcp_options)")

        # ----------------------------------------------------------------------
        # Settings

        emissions_path_keys = Dict(
            "RCP1.9" => "ssp119",  # generally paired with SSP1
            "RCP2.6" => "ssp126",  # generally paired with SSP1
            "RCP4.5"=> "ssp245",  # generally paired with SSP2
            "RCP3.7" => "ssp370",  # generally paired with SSP3
            "RCP8.5" => "ssp585"   # generally paired with SSP4
        )

        socioeconomic_path = joinpath(@__DIR__, "..", "..", "data", "socioeconomic", "ssp_projections_$(p.model)_$(p.ssp).csv")
        emissions_path = joinpath(@__DIR__, "..", "..", "data", "emissions", "rcmip_$(emissions_path_keys[p.ssp])_emissions_1750_to_2500.csv")

        # ----------------------------------------------------------------------
        # Load Data as Needed

        dict_key = Symbol(p.model, "-", p.ssp)
        if !haskey(g_datasets, dict_key)

            # interpolate socioeconomic data to annual

            socioeconomic_data = load(socioeconomic_path) |>
                DataFrame|>
                @select(:year, :region, :pop, :gdp) |>
                DataFrame

            socioeconomic_data_interp = DataFrame(:year => [], :region => [], :pop => [], :gdp => [])
            all_years = collect(minimum(socioeconomic_data.year):maximum(socioeconomic_data.year))
            for region in unique(socioeconomic_data.region)

                region_data = socioeconomic_data |> @filter(_.region == region) |> @orderby(:year) |> DataFrame
                pop_itp = LinearInterpolation(region_data.year, region_data.pop)
                gdp_itp = LinearInterpolation(region_data.year, region_data.gdp)

                append!(socioeconomic_data_interp, DataFrame(
                    :year => all_years,
                    :region => fill(region, length(all_years)),
                    :pop => pop_itp[all_years],
                    :gdp => gdp_itp[all_years]
                ))
            end
            for (i, type) in enumerate([Int64, String, Float64, Float64])
                socioeconomic_data_interp[:,i] = convert.(type, socioeconomic_data_interp[:,i])
            end

            # add interpolated socioeconomic data, and emissions data which is already
            # annual, to the datasets global variable
            g_datasets[dict_key] = Dict(
                :socioeconomic => socioeconomic_data_interp,
                :emissions => load(emissions_path, skiplines_begin = 6)|>
                    DataFrame |>
                    i -> rename!(i, Symbol.([:year, names(i)[2:end]...])) |>
                    DataFrame |>
                    @select(:year, :carbon_dioxide, :nitrous_oxide, :methane, :sf6) |>
                    DataFrame
                )
        end

        # ----------------------------------------------------------------------
        # Check Countries - each country found in the model countries parameter
        # must exist in the ssp socioeconomics dataframe

        missing_countries = []
        for country in p.country_names
            !(country in unique(g_datasets[dict_key][:socioeconomic].region)) && push!(missing_countries, country)
        end
        !isempty(missing_countries) && error("All countries in countries parameter must be found in SSPs component Socioeconomic Dataframe, the following were not found: $(missing_countries)")
    end

    function run_timestep(p,v,d,t)

        dict_key = Symbol(p.model, "-", p.ssp)
        year_label = gettime(t)

        # check that we only run the component where we have data
        if !(year_label in unique(g_datasets[dict_key][:socioeconomic].year))
            error("Cannot run SSP component in year $(year_label), SSP socioeconomic variables only available for 2010 through 2300.")
        end
        if !(year_label in unique(g_datasets[dict_key][:emissions].year))
            error("Cannot run SSP component in year $(year_label), SSP emissions variables only available for 1750 through 2500.")
        end

        # ----------------------------------------------------------------------
        # Socioeconomic

        # filter the dataframe for values with the year matching timestep
        # t and only the ssp countries found in the model countries list,
        # already checked that all model countries are in ssp countries list
        subset = g_datasets[dict_key][:socioeconomic] |>
            @filter(_.year == gettime(t) && _.region in p.country_names) |>
            DataFrame

        # get the ordered indices of the ssp countries within the parameter
        # of the model countries, already checked that all model countries
        # are in ssp countries list
        order = indexin(p.country_names, subset.region)

        v.population[t,:] = subset.pop[order]
        v.gdp[t,:] = subset.gdp[order]

        # ----------------------------------------------------------------------
        # Emissions

        subset = g_datasets[dict_key][:emissions] |>
                    @filter(_.year == gettime(t)) |>
                    DataFrame

        v.co2_emissions[t] = subset.carbon_dioxide[1]
        v.ch4_emissions[t] = subset.methane[1]
        v.n2o_emissions[t] = subset.nitrous_oxide[1]
        v.sf6_emissions[t] = subset.sf6[1]

    end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule SHTOOLS_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("SHTOOLS")
JLLWrappers.@generate_main_file("SHTOOLS", UUID("ff580fc6-ede1-5b7e-942f-684e38381daf"))
end  # module SHTOOLS_jll


text:
# Autogenerated wrapper script for WhiteboxTools_jll for x86_64-linux-musl
export conditional_evaluation, edge_contamination, exposure_towards_wind_flux, qin_flow_accumulation, quinn_flow_accumulation, raster_calculator, rho8_flow_accumulation, split_vector_lines, whitebox_tools

JLLWrappers.@generate_wrapper_header("WhiteboxTools")
JLLWrappers.@declare_executable_product(conditional_evaluation)
JLLWrappers.@declare_executable_product(edge_contamination)
JLLWrappers.@declare_executable_product(exposure_towards_wind_flux)
JLLWrappers.@declare_executable_product(qin_flow_accumulation)
JLLWrappers.@declare_executable_product(quinn_flow_accumulation)
JLLWrappers.@declare_executable_product(raster_calculator)
JLLWrappers.@declare_executable_product(rho8_flow_accumulation)
JLLWrappers.@declare_executable_product(split_vector_lines)
JLLWrappers.@declare_executable_product(whitebox_tools)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        conditional_evaluation,
        "bin/conditional_evaluation",
    )

    JLLWrappers.@init_executable_product(
        edge_contamination,
        "bin/edge_contamination",
    )

    JLLWrappers.@init_executable_product(
        exposure_towards_wind_flux,
        "bin/exposure_towards_wind_flux",
    )

    JLLWrappers.@init_executable_product(
        qin_flow_accumulation,
        "bin/qin_flow_accumulation",
    )

    JLLWrappers.@init_executable_product(
        quinn_flow_accumulation,
        "bin/quinn_flow_accumulation",
    )

    JLLWrappers.@init_executable_product(
        raster_calculator,
        "bin/raster_calculator",
    )

    JLLWrappers.@init_executable_product(
        rho8_flow_accumulation,
        "bin/rho8_flow_accumulation",
    )

    JLLWrappers.@init_executable_product(
        split_vector_lines,
        "bin/split_vector_lines",
    )

    JLLWrappers.@init_executable_product(
        whitebox_tools,
        "bin/whitebox_tools",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Documenter, gazel

makedocs(;
    modules=[gazel],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/dnabanita7/gazel.jl/blob/{commit}{path}#L{line}",
    sitename="gazel.jl",
    authors="Nabanita Dash, Biswajit Ghosh, Gautam Mishra",
    assets=String[],
)

deploydocs(;
    repo="github.com/dnabanita7/gazel.jl",
)


text:
using FeatureRegistries
using ReTest

FeatureRegistries.runtests()


text:
using CGP
using Logging
using ArgParse

pmin = -100
pmax = 100

function regression(c::Chromosome, data::Array{Float64}, nin::Int64, nout::Int64)
    error = 0
    nsamples = size(data, 2)
    for d in 1:nsamples
        outputs = process(c, data[1:nin, d])
        for p in eachindex(outputs)
            error += (outputs[p] - data[nin+p, d])^2
        end
    end
    error /= nsamples
    -error
end

function izhikevich(u::Float64, v::Float64, I::Float64; a::Float64=0.02,
                    b::Float64=0.2, c::Float64=-65.0, d::Float64=2.0)
    if v >= 30
        v = c
        u = u + d
    end

    v = v+0.5*(0.04*v^2+5*v+140-u+I)
    v = v+0.5*(0.04*v^2+5*v+140-u+I)
    u = u+a*(b*v-u)

    v, u
end

function get_args()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--seed"
        arg_type = Int
        default = 0
        "--log"
        arg_type = String
        required = true
        "--function"
        arg_type = String
        default = "izhikevich"
        "--cfg"
        arg_type = String
        default = "cfg/base.yaml"
    end

    CGP.Config.add_arg_settings!(s)
end

args = parse_args(get_args())
println(args)
CGP.Config.init(Dict([k=>args[k] for k in setdiff(
    keys(args), ["seed", "log", "function", "cfg"])]...))

CGP.Config.init("cfg/base.yaml")
CGP.Config.init("cfg/classic.yaml")
# CGP.Config.init(args["cfg"])

srand(args["seed"])
Logging.configure(filename=args["log"], level=INFO)
ea = oneplus
ctype = CGP.PCGPChromo

nin = 3
nout = 2
# f = eval(parse(args["function"]))
fit = x->fmatch(x, izhikevich, nin, nout)
maxfit, best = ea(ctype, nin, nout, fit; seed=args["seed"])

best_ind = ctype(best, nin, nout)
test_fit = fitness(best_ind, test, nin, nout)
Logging.info(@sprintf("T: %d %0.8f %0.8f %d %d %s %s",
                      args["seed"], maxfit, test_fit,
                      sum([n.active for n in best_ind.nodes]),
                      length(best_ind.nodes), args["ea"], args["chromosome"]))

Logging.info(@sprintf("E%0.8f", -maxfit))


text:
const SymAnyDict = Dict{Symbol, Any}
const Vec{N, T} = SVector{N, T}
const Point = Vec
const Point3 = Point{3}
const Vec3   = Vec{3}
const Mat3{T} = SMatrix{3, 3, T}
const Mat4{T} = SMatrix{4, 4, T}

Point{N, T}(x::T) where {N, T} = Point{N, T}(x, x, x)

Base.convert(::Type{Point3{T}}, x::Vector{T}) where T<:AbstractFloat = Point3{T}(x[1], x[2], x[3])

abstract type Band end

"""
Energy band from DFT calculation.
"""
mutable struct DFBand{T<:AbstractFloat,K} <: Band
    k_points_cart  ::Vector{Vec3{K}}
    k_points_cryst ::Vector{Vec3{T}}
    eigvals        ::Vector{T}
    extra          ::Dict{Symbol, Any}
end
DFBand(k_points_cart, k_points_cryst, eigvals) = DFBand(k_points_cart, k_points_cryst, eigvals, Dict{Symbol,Any}())
DFBand(::Type{T}, vlength::Int) where T = DFBand(Vector{Vec3{T}}(undef, vlength), Vector{Vec3{T}}(undef, vlength), Vector{T}(undef, vlength), SymAnyDict())
DFBand(vlength::Int) = DFBand(Float64, vlength)

kpoints(band::DFBand, kind=:cryst) = kind == :cart ? band.k_points_cart : band.k_points_cryst
eigvals(band::DFBand) = band.eigvals

"""
	bandgap(bands::AbstractVector{DFBand}, fermi=0.0)

Calculates the bandgap (possibly indirect) around the fermi level.
"""
function bandgap(bands::AbstractVector{<:Band}, fermi=0.0)
	max_valence = -Inf
	min_conduction = Inf
	for b in bands
		max = maximum(eigvals(b).-fermi)
		min = minimum(eigvals(b).-fermi)
		if max_valence <= max <= 0.0
			max_valence = max
		end
		if 0.0 <= min <= min_conduction
			min_conduction = min
		end
	end
	return min_conduction - max_valence
end

mutable struct ExecFlag
    symbol     ::Symbol
    name       ::String
    typ        ::Type
    description::String
    value
end

ExecFlag(e::ExecFlag, value) = ExecFlag(e.symbol, e.name, e.typ, e.description, value)
ExecFlag(p::Pair{Symbol, T}) where T = ExecFlag(first(p), String(first(p)), T, "", last(p))

const QE_EXECFLAGS = ExecFlag[
    ExecFlag(:nk, "kpoint-pools", Int, "groups k-point parallelization into nk processor pools", 0),
    ExecFlag(:ntg, "task-groups", Int, "FFT task groups", 0),
    ExecFlag(:ndiag, "diag", Int, "Number of processes for linear algebra", 0),
    ExecFlag(:ni, "images", Int, "Number of processes used for the images", 0)
]

qeexecflag(flag::AbstractString) = getfirst(x -> x.name==flag, QE_EXECFLAGS)
qeexecflag(flag::Symbol) = getfirst(x -> x.symbol==flag, QE_EXECFLAGS)

function parse_qeexecflags(line::Vector{<:AbstractString})
    flags = ExecFlag[]
    i=1
    while i<=length(line)
        s = strip(line[i], '-')
        push!(flags, ExecFlag(qeexecflag(Symbol(s)), parse(Int, line[i+1])))
        i += 2
    end
    flags
end

const WAN_EXECFLAGS = ExecFlag[
    ExecFlag(:pp, "preprocess", Nothing, "Whether or not to preprocess the wannier input", nothing),
]

wan_execflag(flag::AbstractString) = getfirst(x -> x.name==flag, WAN_EXECFLAGS)
wan_execflag(flag::Symbol) = getfirst(x -> x.symbol==flag, WAN_EXECFLAGS)
function parse_wan_execflags(line::Vector{<:AbstractString})
    flags = ExecFlag[]
    i=1
    while i<=length(line)
        s = strip(line[i], '-')
        push!(flags, ExecFlag(wan_execflag(Symbol(s)), nothing))
        i += 1
    end
    flags
end

include(joinpath(depsdir, "mpirunflags.jl"))
const MPI_FLAGS = _MPI_FLAGS()

mpi_flag(flag::AbstractString) = getfirst(x -> x.name==flag, MPI_FLAGS)
mpi_flag(flag::Symbol) = getfirst(x -> x.symbol==flag, MPI_FLAGS)

function mpi_flag_val(::Type{String}, line, i)
    v = line[i+1]
    return v, i + 2
end

function mpi_flag_val(::Type{Vector{String}}, line, i)
    tval = String[]
    while i+1 <= length(line) && !occursin('-', line[i+1])
        push!(tval, line[i+1])
        i += 1
    end
    return tval, i + 1
end

function mpi_flag_val(::Type{Int}, line, i)
    if line[i+1][1] == '\$'
        tval = line[i+1]
    else
        tval = parse(Int, line[i+1])
    end
    return tval, i + 2
end

function mpi_flag_val(::Type{Vector{Int}}, line, i)
    tval = Union{Int, String}[]
    while i+1 <= length(line) && !occursin('-', line[i+1])
        if line[i+1][1] == '\$'
            push!(tval, line[i+1])
        else
            push!(tval, parse(Int, line[i+1]))
        end
        i += 1
    end
    return tval, i + 1
end

function mpi_flag_val(::Type{Pair{Int, String}}, line, i)
    tval = Union{Int, String}[]
    while i+1<=length(line) && !occursin('-', line[i+1])
        if line[i+1][1] == '\$'
            push!(tval, line[i+1])
        else
            tparse = tryparse(Int, line[i+1])
            if tparse != nothing
                push!(tval, tparse)
            else
                push!(tval, string(line[i+1]))
            end
        end
        i += 1
    end
    return tval, i + 1
end
function parse_mpi_flags(line::Vector{<:SubString})
    eflags = ExecFlag[]
    i = 1
    while i <= length(line)
        s = line[i]
        if s[1] != '-'
            break
        end
        if s[2] == '-'
            mflag = mpi_flag(strip(s, '-'))
        else
            mflag = mpi_flag(Symbol(strip(s, '-')))
        end

        @assert mflag != nothing "$(strip(s, '-')) is not a recognized mpi_flag"
        val, i = mpi_flag_val(mflag.typ, line, i)
        push!(eflags, ExecFlag(mflag, val))
    end
    eflags
end

const RUN_EXECS = ["mpirun","mpi","mpiexec"]
allexecs() = vcat(RUN_EXECS, QE_EXECS, WAN_EXECS, ELK_EXECS)
parseable_execs() = vcat(QE_EXECS, WAN_EXECS, ELK_EXECS)
has_parseable_exec(l::String) = occursin(">", l) && any(occursin.(parseable_execs(), (l,)))

mutable struct Exec
    exec ::String
    dir  ::String
    flags::Vector{ExecFlag}
end

Exec() = Exec("")
Exec(exec::String) = Exec(exec, "")
Exec(exec::String, dir::String) = Exec(exec, dir, ExecFlag[])
Exec(exec::String, dir::String, flags...) = Exec(exec, dir, SymAnyDict(flags))
function Exec(exec::String, dir::String, flags::SymAnyDict)
    _flags = ExecFlag[]
    for (f, v) in flags
        push!(_flags, ExecFlag(f => v))
    end
    return Exec(exec, dir, _flags)
end

isparseable(exec::Exec) = exec.exec ∈ parseable_execs()

is_qe_exec(exec::Exec)      = exec.exec ∈ QE_EXECS
is_wannier_exec(exec::Exec) = exec.exec ∈ WAN_EXECS
is_elk_exec(exec::Exec)     = exec.exec ∈ ELK_EXECS
is_mpi_exec(exec::Exec)     = occursin("mpi", exec.exec)

function inputparser(exec::Exec)
    if is_qe_exec(exec)
        qe_read_input
    elseif is_wannier_exec(exec)
        wan_read_input
    elseif is_elk_exec(exec)
        elk_read_input
    end
end

function setflags!(exec::Exec, flags...)
    for (f, val) in flags
        flag = isa(f, String) ? getfirst(x -> x.name == f, exec.flags) : getfirst(x -> x.symbol == f, exec.flags)
        if flag != nothing
            flag.value = convert(flag.typ, val)
        else
	        for (f1, f2) in zip((is_qe_exec, is_wannier_exec, is_mpi_exec), (qeexecflag, wan_execflag, mpi_flag))
	        	if f1(exec)
		        	def_flag = f2(f)
			        if def_flag != nothing
				        push!(exec.flags, ExecFlag(def_flag, val))
				        break
			        end
		        end
	        end
        end
    end
    exec.flags
end

function rmflags!(exec::Exec, flags...)
    for f in flags
        if isa(f, String)
            filter!(x -> x.name != f, exec.flags)
        elseif isa(f, Symbol)
            filter!(x -> x.symbol != f, exec.flags)
        end
    end
    exec.flags
end
hasflag(exec::Exec, s::Symbol) = findfirst(x->x.symbol == s, exec.flags) != nothing
setexecdir!(exec::Exec, dir) = exec.dir = dir


text:
export Propensity, TimeInvariantPropensity, TimeVaryingPropensity, propensity, propensity, istimeseparable, istimevarying
export StandardTimeInvariantPropensity, SeparableTimeVaryingPropensity, JointTimeVaryingPropensity
export PropensityGradient, TimeVaryingPropensityGradient, JointTimeVaryingPropensityGradient, SeparableTimeVaryingPropensityGradient, propensitygrad, propensitygrad_timevarying

"""
Base type for all propensities.
"""
abstract type Propensity end

abstract type TimeInvariantPropensity <: Propensity end
abstract type TimeVaryingPropensity <: Propensity end

istimevarying(α::Propensity) = false
istimeseparable(α::Propensity) = false

"""
    StandardTimeInvariantPropensity

Propensity function for a reaction with rate not dependent on time.

# Examples

A standard time-invariant propensity can be created by wrapping around a Julia function. The function must accept two inputs `x`, `p` where `x` is the CME state vector and `p` is the vector of parameters.
```@repl
julia> f(x, p) = p[1]*x[1]
f (generic function with 1 method)
julia> α = StandardTimeInvariantPropensity(f)
StandardTimeInvariantPropensity(f)
julia> α([1,2], [0.1,0.2])
0.1
```

One can also input `f` into the `propensity()` function. Provided that `f` is a function of two arguments with the CME state argument followed by the parameter argument, `propensity()` will output a standard time-invariant propensity. Continuing the example above,
```@repl
julia> β = propensity(f)
StandardTimeInvariantPropensity(f)
julia> β([1,2], [0.1,0.2]) == α([1,2], [0.1,0.2])
true
```

One can also take advantage of Julia's `do` block to make the code more concise.
```@repl
julia> γ = propensity() do x,p
        x[1]*p[1]
    end
julia> γ([1,2], [0.1,0.2])
0.1
```
"""
struct StandardTimeInvariantPropensity <: TimeInvariantPropensity
    f::Any
end

"""
    SeparableTimeVaryingPropensity

Reaction propensity with a time-varying rate that can be factorized into a function that depends only on the time variable (and parameters) and another function that does not depend on time. Subsequent types that build upon this type (such as `FspMatrixSparse`) takes advantage of time-space separability to avoid repeating expensive computations.

# Examples

A separable time-varying propensity can be created by specifying two Julia functions. The first function represents the time-varying factor and accepts two input arguments `t` (for time) and `p` (for parameters), the second function represents the time-invariant factor and accepts two input arguments `x` (for CME state) and `p` (for parameters).
```@repl
julia> c(t,p) = (1.0+cos(π*t/p[2]))
c (generic function with 1 method)
julia> f(x, p) = p[1]*x[1]
f (generic function with 1 method)
julia> α = SeparableTimeVaryingPropensity(c,f)
SeparableTimeVaryingPropensity(c,f)
julia> α(20.0,[1,2], [0.1,0.2])
0.2
```

One can also use a method of the `propensity()` function to create a `SeparableTimeVaryingPropensity` instance. Provided that `f` is a function of two arguments with the CME state argument followed by the parameter argument, `propensity(f,c)` (note that the time-varying factor comes _after_ the time-invariant factor) will output a standard time-invariant propensity. Continuing the example above,
```@repl
julia> β = propensity(f,c)
SeparableTimeVaryingPropensity(f, c)
julia> β(20.0, [1,2], [0.1,0.2]) == α(20.0, [1,2], [0.1,0.2])
true
```

One can also take advantage of Julia's `do` block to make the code more concise.
```@repl
julia> γ = propensity((1.0+cos(π*t/p[2]))) do x,p
        x[1]*p[1]
    end
julia> SeparableTimeVaryingPropensity(var"#19#21"(), var"#18#20"())
julia> γ(20.0, [1,2], [0.1,0.2])
0.2
```
"""
struct SeparableTimeVaryingPropensity <: TimeVaryingPropensity
    tfactor::Union{Any,Nothing}
    statefactor::Any
end

"""
    JointTimeVaryingPropensity

Reaction propensity with a time-varying rate. For speed, we recommend using `SeparableTimeVaryingPropensity` if the propensity function has a separable structure.

# Examples

A time-varying propensity can be created by wrapping a function with three input arguments, `t` (for time) `x` (for CME state) and `p` (for parameters).

```@repl
julia> f(t, x, p) = (1.0+cos(π*t/p[2]))*p[1]*x[1]
f (generic function with 1 method)
julia> α = JointTimeVaryingPropensity(c,f)
JointTimeVaryingPropensity(c,f)
julia> α(20.0,[1,2], [0.1,0.2])
0.2
```

One can also use a method of the `propensity()` function to create a `JointTimeVaryingPropensity` instance. Provided that `f` is a function of three arguments with time as the first argument, the CME state as the second argument followed by the parameter argument.

```@repl
julia> β = propensity(f)
julia> β(20.0, [1,2], [0.1,0.2]) == α(20.0, [1,2], [0.1,0.2])
true
```
"""
struct JointTimeVaryingPropensity <: TimeVaryingPropensity
    f::Any
end

istimevarying(α::TimeVaryingPropensity) = true
istimeseparable(α::SeparableTimeVaryingPropensity) = true

function propensity(f::Any)
    m = methods(f)[1]
    if m.nargs == 3
        return StandardTimeInvariantPropensity(f)
    elseif m.nargs == 4
        return JointTimeVaryingPropensity(f)
    else
        throw(ArgumentError("The callable passed to `propensity()` must have either two arguments (x,p) or three arguments (t,x,p)."))
    end
end

function propensity(xfactor::Any, tfactor::Any)
    return SeparableTimeVaryingPropensity(tfactor, xfactor)
end

# Call methods to make Propensity objects callable
function (a::StandardTimeInvariantPropensity)(x::AbstractVector, p=[])
    a.f(x,p)
end
function (a::JointTimeVaryingPropensity)(t::Real, x::AbstractVector, p=[])
    a.f(t,x,p)
end
function (a::SeparableTimeVaryingPropensity)(t::Real, x::AbstractVector, p=[])
    a.tfactor(t,p)*a.statefactor(x,p)
end


text:
export WishartGaussian, GaussianDiagonal, DirichletMultinomial, GammaNormal, NormalNormal, BetaBernoulli, ConjugatePostDistribution

abstract ConjugatePostDistribution

abstract UnivariateConjugatePostDistribution <: ConjugatePostDistribution
abstract DiscreteUnivariateConjugatePostDistribution <: UnivariateConjugatePostDistribution
abstract ContinuousUnivariateConjugatePostDistribution <: UnivariateConjugatePostDistribution

abstract MultivariateConjugatePostDistribution <: ConjugatePostDistribution
abstract DiscreteMultivariateConjugatePostDistribution <: MultivariateConjugatePostDistribution
abstract ContinuousMultivariateConjugatePostDistribution <: MultivariateConjugatePostDistribution

# Gaussian with Normal Inverse Wishart Prior
type WishartGaussian <: ContinuousMultivariateConjugatePostDistribution

    D::Int

    # sufficient statistics
    n::Int
    sums::Vector{Float64}
    ssums::Array{Float64}

    # base model parameters
    mu0::Vector{Float64}
    kappa0::Float64
    nu0::Float64
    Sigma0::Array{Float64}

    function WishartGaussian(mu::Vector{Float64}, kappa::Float64,
            nu::Float64, Sigma::Array{Float64})

        d = length(mu)

        new(d, 0, zeros(d), zeros(d, d), mu, kappa, nu, Sigma)
    end

    function WishartGaussian(D::Int, mu::Vector{Float64}, kappa::Float64,
            nu::Float64, Sigma::Array{Float64})

        new(D, 0, zeros(d), zeros(d, d), mu, kappa, nu, Sigma)
    end

    function WishartGaussian(D::Int,
            n::Int, sums::Vector{Float64}, ssums::Array{Float64},
            mu::Vector{Float64}, kappa::Float64,
            nu::Float64, Sigma::Array{Float64})

        new(D, n, sums, ssums, mu, kappa, nu, Sigma)
    end

end

Base.show(io::IO, d::WishartGaussian) =
    show_multline(io, d, [(:dim, d.D), (:μ0, d.mu0), (:Σ0, d.Sigma0), (:κ0, d.kappa0), (:ν0, d.nu0)])

# Normal with Gamma prior
type GammaNormal <: ContinuousUnivariateConjugatePostDistribution

	# sufficient statistics
	n::Int
	sums::Float64
  ssums::Float64

	# model parameters
	μ0::Float64
	λ0::Float64
	α0::Float64
	β0::Float64

	function GammaNormal(;μ = 0.0, λ = 1.0, α = 1.0, β = 1.0)
		new(0, 0.0, 0.0, μ, λ, α, β)
	end

end

Base.show(io::IO, d::GammaNormal) =
    show_multline(io, d, [(:μ0, d.μ0), (:λ0, d.λ0), (:α0, d.α0), (:β0, d.β0)])

# Normal with Normal prior
type NormalNormal <: ContinuousUnivariateConjugatePostDistribution

	# sufficient statistics
	n::Int
	sums::Float64
  ssums::Float64

	# model parameters
	μ0::Float64
	σ0::Float64

	function NormalNormal(;μ0 = 0.0, σ0 = 1.0)
		new(0, 0.0, 0.0, μ0, σ0)
	end

end

Base.show(io::IO, d::NormalNormal) =
    show_multline(io, d, [(:μ0, d.μ0), (:σ0, d.σ0)])

# Gaussian with Diagonal Covariance
type GaussianDiagonal{T <: ContinuousUnivariateConjugatePostDistribution} <: ContinuousMultivariateConjugatePostDistribution

    # sufficient statistics
    dists::Vector{T}

    function GaussianDiagonal(dists::Vector{T})
        new(dists)
    end

end

# Multinomial with Dirichlet Prior
type DirichletMultinomial <: DiscreteMultivariateConjugatePostDistribution

    D::Int

    # sufficient statistics
    n::Int
    counts::SparseMatrixCSC{Int,Int}

    # base model parameters
    alpha0::Float64

  	# cache
  	dirty::Bool
  	Z2::Float64
  	Z3::Array{Float64}

    function DirichletMultinomial(D::Int, alpha::Float64)
        new(D, 0, sparsevec(zeros(D)), alpha, true, 0.0, Array{Float64}(0))
    end

    function DirichletMultinomial(N::Int, counts::Vector{Int},
                            D::Int, alpha::Float64)
        new(D, N, sparsevec(counts), alpha, true, 0.0, Array{Float64}(0))
    end

    function DirichletMultinomial(N::Int, counts::SparseMatrixCSC{Int,Int},
                            D::Int, alpha::Float64)
        new(D, N, counts, alpha, true, 0.0, Array{Float64}(0))
    end

end

# Bernoulli with Beta Prior
type BetaBernoulli <: DiscreteUnivariateConjugatePostDistribution

	# sufficient statistics
	successes::Int
  n::Int

	# beta distribution parameters
	α::Float64
	β::Float64

	function BetaBernoulli(;α = 1.0, β = 1.0)
		new(0, 0, α, β)
	end

end


text:
"""
Acquisition geometry for each supersource.
"""
mutable struct AGeomss
    s::NamedStack{Vector{Float64}} # x, y, z coordinates
    r::NamedStack{Vector{Float64}} # x, y, z coordinates
    # sx::Vector{Float64}
    # sz::Vector{Float64}
    # rx::Vector{Float64}
    # rz::Vector{Float64}
    ns::Int64 # change to +ve only integer later
    nr::Int64 # change to +ve only integer later
    "adding conditions that are to be false while construction"
    AGeomss(s, r, ns, nr) =
        (
            any(names(s)[1] != names(r)[1]) |
            any([length(ss) != ns for ss in s]) |
            any([length(rr) != nr for rr in r])
        ) ? error("AGeomss construct") : new(s, r, ns, nr)
end # type

include("gallery.jl")

"""
Randomly place a given number of source and receivers in a `mgrid`.
"""
function AGeomss(mgrid::AbstractArray{T}, s::Srcs, r::Recs) where {T<:StepRangeLen}
    nd = length(mgrid)
    names = dim_names(nd)
    mins = [m[1] for m in mgrid]
    maxs = [m[end] for m in mgrid]
    return AGeomss(
        NamedArray([Random.rand(Uniform(mins[i], maxs[i]), s.n) for i = 1:nd], (names,)),
        NamedArray([Random.rand(Uniform(mins[i], maxs[i]), r.n) for i = 1:nd], (names,)),
        s.n,
        r.n,
    )
end

"""
```julia
ageom=AGeom(mgrid, SSrcs(2), [Srcs(3), Srcs(2)], [Recs(10), Recs(20)])
```
Initialize an acquisition, with 2 supersources, by randomly placing a
sources and receivers in `mgrid`.
"""
function Array{AGeomss,1}(
    mgrid::AbstractVector{T},
    ss::SSrcs,
    s::Vector{Srcs},
    r::Vector{Recs},
) where {T<:StepRangeLen}
    return [AGeomss(mgrid, s[i], r[i]) for i = 1:ss.n]
end

"""
```julia
ageom=AGeom(mgrid, SSrcs(10), Srcs(10), Recs(20))
```
Initialize an acquisition, with 10 supersources, by randomly placing a
sources and receivers in `mgrid`. In this case, all the supersources have same source and receiver positions.
"""
function Array{AGeomss,1}(
    mgrid::AbstractVector{T},
    ss::SSrcs,
    s::Srcs,
    r::Recs,
) where {T<:StepRangeLen}
    return Array{AGeomss,1}(mgrid, ss, fill(s, ss.n), fill(r, ss.n))
end

function Base.isequal(ageom1::AGeomss, ageom2::AGeomss, attrib = :all)
    @assert dims(ageom1) == dims(ageom2) "attempt to compare acquisitions of different dimensions"
    names = [:sx, :sz, :ns]
    recfnames = [:rx, :rz, :nr]
    src =
        all([(isequal(ageom1.s[name], ageom2.s[name])) for name in names(ageom1.s)[1]]) & (ageom1.ns == ageom2.ns)
    rec =
        all([
            (isequal(getfield(ageom1, name), getfield(ageom2, name))) for name in recfnames
        ]) & (ageom1.nr == ageom2.nr)
    if (attrib == :sources)
        return src
    elseif (attrib == :receivers)
        return rec
    else
        return src & rec
    end
end

"""
```julia
isequal(ageom1, ageom2)
```
Assert if `ageom1` equals `ageom2`.
"""
function Base.isequal(ageom1::AGeom, ageom2::AGeom, attrib = :all)
    return (length(ageom1) == length(ageom2)) ?
           all([isequal(ageom1[iss], ageom2[iss], attrib) for iss = 1:length(ageom1)]) :
           false
end

function Base.show(io::Base.IO, ageomss::AGeomss)
    println(
        io,
        string(
            "> ",
            ndims(ageomss),
            "-D acquisition with ",
            ageomss.ns,
            " source(s) and ",
            ageomss.nr,
            " receiver(s)",
        ),
    )
end

function Base.show(io::Base.IO, ageom::AGeom)
    println(
        io,
        string(
            "> ",
            ndims(ageom),
            "-D acquisition with ",
            length(ageom),
            " supersource(s), ",
            getfield.(ageom, :ns),
            " source(s) and ",
            getfield.(ageom, :nr),
            " receiver(s)",
        ),
    )
end

function Base.ndims(ageomss::AGeomss)
    return length(ageomss.s)
end

function Base.ndims(ageom::Array{AGeomss,1})
    nd = ndims.(ageom)
    @assert all(y -> y == nd[1], nd)
    return nd[1]
end

"""
Output `n` interpolated sources/ receivers between positions `p1` and `p2`.
"""
function spread(
    n,
    p1::Vector{T1},
    p2::Vector{T2},
) where {T1<:Real,T2<:Real}
    @assert (length(p1) == length(p2)) # check if 2D or 3D
    nd = length(p1)
    if (n == 1)
        return Tuple([[p1[i]] for i = 1:nd])
    else
        knots = ([1, n],)
        itp = [
            Interpolations.interpolate(
                knots,
                [Float64(p1[i]), Float64(p2[i])],
                Gridded(Linear()),
            ) for i = 1:nd
        ]
        return Tuple([[it(i) for i = 1:n] for it in itp])
    end
end

"""
Spread using difference sources/ receivers in each dimension given by `n_per_dim`.
"""
function spread(
    n,
    p1::Vector{T1},
    p2::Vector{T2},
    n_per_dim::Vector{T3}
) where {T1<:Real,T2<:Real,T3<:Int}
    @assert length(p1)==length(p2)==length(n_per_dim)
    nd = length(p1)
    @assert prod(n_per_dim) == n
    I=collect(Iterators.product([spread(n_per_dim[i], [p1[i]], [p2[i]])[1] for i in 1:length(n_per_dim)]...))
    return Tuple([[I[i][id] for i = 1:n] for id=1:nd])
end

"""
theta=0 means in x-z plane
"""
function spread(
    n,
    p0::Vector{T1},
    rad::T2,
    phi::Vector{T4},
    theta::Vector{T3} = [0, 0],
    rand_flag::Bool = false,
) where {T1<:Real,T2<:Real,T3<:Real,T4<:Real}
    phi = spread(n, [phi[1]], [phi[2]], rand_flag)[1]
    theta = spread(n, [theta[1]], [theta[2]], rand_flag)[1]
    nd = length(p0)
    p = [zeros(n) for i = 1:nd]

    for i = 1:n
        p[1][i] = p0[1] + rad * cos(phi[i])
        if (nd > 1)
            p[2][i] = p0[2] + rad * sin(phi[i]) * sin(theta[i])
        end
        if (nd > 2)
            p[3][i] = p0[3] + rad * sin(phi[i]) * cos(theta[i])
        end
    end
    return Tuple(p)
end

"""
```julia
update!(ageom[1], Srcs(10), p1, p2)
```
Update source geometry for the first supersource, with 10 sources placed regularly between `p1=[z1,x1]` and `p2=[z2,x2]`

```julia
update!(ageom[2], Srcs(10), p0, rad, angles)
```
Update source geometry for the second supersource, with 10 sources placed equidistant from `p0=[z0,x0]` given radius `rad`
and angles `angles=[0,2pi]`.
"""
function update!(ageomss::AGeomss, ::Srcs, args...)
    p = spread(ageomss.ns, args...)
    for id = 1:ndims(ageomss)
        copyto!(ageomss.s[id], p[id])
    end
    return ageomss
end

"""
```julia
update!(ageom[1], Recs(10), args...)
```
Similar to updating source positions, but for receivers.
"""
function update!(ageomss::AGeomss, ::Recs, args...)
    p = spread(ageomss.nr, args...)
    for id = 1:ndims(ageomss)
        copyto!(ageomss.r[id], p[id])
    end
    return ageomss
end

"""
```julia
update!(ageom, SSrcs(), p1, p2)
update!(ageom, SSrcs(), p0, rad, angles)
```
Update all source positions in each supersources.
"""
function update!(ageom::AGeom, ::SSrcs, args...)
    p = spread(length(ageom), args...)
    pz = collect(zip(p...))
    for iss = 1:length(ageom)
        update!(ageom[iss], Srcs(), collect(pz[iss]), collect(pz[iss]))
    end
    return ageom
end

"""
```julia
update!(ageom, Recs(), p1, p2)
update!(ageom, Recs(), p0, rad, angles)
```
Update receiver positions for all supersources.
"""
function update!(ageom::AGeom, ::Recs, args...)
    for iss = 1:length(ageom)
        update!(ageom[iss], Recs(), args...)
    end
    return ageom
end

"""
```julia
ageom ∈ mgrid
in(ageom, mgrid)
```
Assert if all the sources and receivers in `AGeom` are within bounds of `mgrid`.
"""
function Base.in(ageom::AGeom, mgrid::AbstractVector{T}) where {T<:StepRangeLen}
    # xmin, zmin, xmax, zmax = mgrid[2][1], mgrid[1][1], mgrid[2][end], mgrid[1][end]
    checkvec = fill(false, length(ageom))
    for iss = 1:length(ageom)
        nd = ndims(ageom[iss])
        @assert nd == length(mgrid)
        names = dim_names(nd)
        checkvec[iss] = any(
            hcat(
                [
                    vcat(
                        (
                            (
                                (ageom[iss].s[dn] .- mgrid[id][1]) .*
                                (mgrid[id][end] .- ageom[iss].s[dn])
                            ) .< 0.0
                        ),
                        (
                            (
                                (ageom[iss].r[dn] .- mgrid[id][1]) .*
                                (mgrid[id][end] .- ageom[iss].r[dn])
                            ) .< 0.0
                        ),
                        isnan.(ageom[iss].s[dn]),
                        isnan.(ageom[iss].r[dn]),
                    ) for (id, dn) in enumerate(names)
                ]...,
            ),
        )
    end
    return !(any(checkvec))
end

"""
```julia
A=SparseArrays.SparseMatrixCSC(ageom[1],mgrid)
D=A*P
```
Return a sparse matrix that restricts the field `P` on `mgrid` to receiver positions, to give `D`.
"""
function SparseArrays.SparseMatrixCSC(
    ageomss::AGeomss,
    mgrid::AbstractVector{T},
) where {T<:StepRangeLen}
    @assert length(mgrid) == 2 #  need to update this routine for 3D acq
    nz, nx = length.(mgrid)
    ACQ = spzeros(ageomss.nr, prod(length.(mgrid)))
    for ir = 1:ageomss.nr
        irx = argmin(abs.(mgrid[2] .- ageomss.r[:x][ir]))
        irz = argmin(abs.(mgrid[1] .- ageomss.r[:z][ir]))
        ACQ[ir, irz+(irx-1)*nz] = 1.0
    end
    return ACQ
end

#=

"""
Advance either source or receiver array in an acquisition ageometry in horizontal or vertical directions.

# Arguments
* `ageom::AGeom` : acquisition ageometry that is updated
* `advances::Vector{Float64}=[[0.,0.], [0.,0.,]]` : source and receiver advancements
"""
function AGeom_advance(ageom::AGeom, advances::Vector{Vector{Float64}}=[[0.,0.], [0.,0.]])
	ageom_out=deepcopy(ageom)
	for iss=1:ageom.nss
		ageom_out.s[:x][iss][:] .+= advances[1][2]
		ageom_out.s[:z][iss][:] .+= advances[1][1]
		ageom_out.r[:x][iss][:] .+= advances[2][2]
		ageom_out.r[:z][iss][:] .+= advances[2][1]
	end
	return ageom_out
end

"""
return a vector of the order
"""
function AGeom_getvec(ageom::Vector{AGeom}, field::Symbol)
	np = length(ageom);
	vect = [getfield(ageom[ip],field) for ip=1:np]
	return vec(vcat(vcat(vect...)...))
end

function save(ageom, folder)
	!(isdir(folder)) && error("invalid directory")
	for iss in 1:ageom.nss
		file=joinpath(folder, string("s", iss, ".csv"))
		CSV.write(file,DataFrame(hcat(ageom.s[:x][iss], ageom.s[:z][iss])))
		file=joinpath(folder, string("r", iss, ".csv"))
		CSV.write(file,DataFrame(hcat(ageom.r[:x][iss], ageom.r[:z][iss])))
	end
end

"""
Appends the input a vector of acquisition ageometries.
"Merge" `AGeom` objects in `ageomv` to `ageom1`.
Typically you have one acquisition,
you change it and want to create the macro acquisition
where the supersources of both acquisitions are appended.
"""
function AGeom_add!(ageom1::AGeom, ageomv::Vector{AGeom})
	for iageom=1:length(ageomv)
		append!(ageom1.s[:x], ageomv[iageom].s[:x])
		append!(ageom1.s[:z], ageomv[iageom].s[:z])
		append!(ageom1.r[:x], ageomv[iageom].r[:x])
		append!(ageom1.r[:z], ageomv[iageom].r[:z])
		append!(ageom1.ns, ageomv[iageom].ns)
		append!(ageom1.nr, ageomv[iageom].nr)
		ageom1.nss += ageomv[iageom].nss
	end
	return ageom1
end

"""
Adds input positions as either sources or receivers of every supershot.
"""
function AGeom_add!(ageom::AGeom, zpos::Vector{Float64}, xpos::Vector{Float64}, attrib::Symbol)
	length(xpos) == length(zpos) ? nothing : error("same length vectors needed")
	for iss=1:ageom.nss
		if(attrib == :receivers)
			ageom.r[:x][iss] = vcat(ageom.r[:x][iss][:], xpos)
			ageom.r[:z][iss] = vcat(ageom.r[:z][iss][:], zpos)
			ageom.nr[iss] += length(xpos)
		elseif(attrib == :sources)
			ageom.s[:x][iss] = vcat(ageom.s[:x][iss][:], xpos)
			ageom.s[:z][iss] = vcat(ageom.s[:z][iss][:], zpos)
			ageom.ns[iss] += length(xpos)
		end
	end
end

"""
A fixed spread acquisition has same set of sources and
receivers for each supersource.
This method constructs a
fixed spread acquisition ageometry using either a
horizontal or vertical array of supersources/ receivers.
Current implementation has only one source for every supersource.

# Arguments

* `ssmin::Float64` : minimum coordinate for sources
* `ssmax::Float64` : maximum coordinate for sources
* `ss0::Float64` : consant coordinate for sources
* `rmin::Float64` : minimum coordinate for receivers
* `rmax::Float64` : maximum coordinate for receivers
* `r0::Float64` : consant coordinate for receivers
* `nss::Int64` : number of supersources
* `nr::Int64` : number of receivers
* `ssattrib::Symbol=:horizontal` : supersource array kind
  `=:vertical` : vertical array of supersources
  `=:horizontal` horizontal array of supersources
* `rattrib::Symbol=:horizontal` : receiver array kind
  `=:vertical` : vertical array of receivers
  `=:horizontal` horizontal array of receivers
* `rand_flags::Vector{Bool}=[false, false]` : decide placement of supersources and receivers
  `=[true, false]` : randomly place supersources for regularly spaced receivers
  `=[true, true]` : randomly place supersources and receivers
  `=[false, false]` : regularly spaced supersources and receivers
  `=[false, true]` : randomly place receivers for regularly spaced supersources

# Return
* a fixed spread acquisition ageometry `AGeom`
"""
function AGeom_fixed(
	      ssmin::Real,
	      ssmax::Real,
	      ss0::Real,
	      rmin::Real,
	      rmax::Real,
	      r0::Real,
	      nss::Int64,
	      nr::Int64,
	      ssattrib::Symbol=:horizontal,
	      rattrib::Symbol=:horizontal,
 	      rand_flags::Vector{Bool}=[false, false];
	      ssα::Real=0.0,
	      rα::Real=0.0,
	      ns::Vector{Int64}=ones(Int,nss),
	      srad::Real=0.0
	     )

	ssmin=Float64(ssmin); rmin=Float64(rmin)
	ssmax=Float64(ssmax); rmax=Float64(rmax)
	ss0=Float64(ss0); r0=Float64(r0)
	ssα=Float64(ssα)*pi/180.
	rα=Float64(rα)*pi/180.

	ssarray = isequal(nss,1) ? fill(ssmin,1) : (rand_flags[1] ?
					     Random.rand(Uniform(minimum([ssmin,ssmax]),maximum([ssmin,ssmax])),nss) : range(ssmin,stop=ssmax,length=nss))
	if(ssattrib==:horizontal)
		ssz = ss0.+(ssarray.-minimum(ssarray)).*sin(ssα)/cos(ssα); ssx=ssarray
	elseif(ssattrib==:vertical)
		ssx = ss0.+(ssarray.-minimum(ssarray)).*sin(ssα)/cos(ssα); ssz=ssarray
	else
		error("invalid ssattrib")
	end

	rarray = isequal(nr,1) ? fill(rmin,1) : (rand_flags[2] ?
					  Random.rand(Uniform(minimum([rmin,rmax]),maximum([rmin,rmax])),nr) : range(rmin,stop=rmax,length=nr))
	if(rattrib==:horizontal)
		rz = r0.+(rarray.-minimum(rarray)).*sin(rα)/cos(rα); rx = rarray
	elseif(rattrib==:vertical)
		rx = r0.+(rarray.-minimum(rarray)).*sin(rα)/cos(rα); rz = rarray
	else
		error("invalid rattrib")
	end
	rxall = [rx for iss=1:nss];
	rzall = [rz for iss=1:nss];
	ssxall=[zeros(ns[iss]) for iss=1:nss];
	sszall=[zeros(ns[iss]) for iss=1:nss];
	for iss in 1:nss
		for is in 1:ns[iss]
			θ=Random.rand(Uniform(-Float64(pi),Float64(pi)))
			r = iszero(srad) ? 0.0 : Random.rand(Uniform(0,srad))
			x=r*cos(θ)
			z=r*sin(θ)
			ssxall[iss][is]=ssx[iss]+x
			sszall[iss][is]=ssz[iss]+z
		end
	end
	return AGeom(ssxall, sszall, rxall, rzall, nss, ns, fill(nr,nss))
end

=#


text:
@testset "Reaction" begin
    m1 = Metabolite("m1")
    m1.formula = "C2H3"
    m2 = Metabolite("m2")
    m2.formula = "H3C2"
    m3 = Metabolite("m3")
    m4 = Metabolite("m4")
    m5 = Metabolite("m5")
    m6 = Metabolite("m6")
    m7 = Metabolite("m7")
    m8 = Metabolite("m8")
    m9 = Metabolite("m9")
    m10 = Metabolite("m10")
    m11 = Metabolite("m11")
    m12 = Metabolite("m12")

    g1 = Gene("g1")
    g2 = Gene("g2")
    g3 = Gene("g3")

    r1 = Reaction()
    r1.id = "r1"
    r1.name = "reaction 1"
    r1.metabolites = Dict(m1.id => -1.0, m2.id => 1.0)
    r1.lb = -100.0
    r1.ub = 100.0
    r1.grr = [["g1", "g2"], ["g3"]]
    r1.subsystem = "glycolysis"
    r1.notes = Dict("notes" => ["blah", "blah"])
    r1.annotations = Dict("sboterm" => ["sbo"], "biocyc" => ["ads", "asds"])
    r1.objective_coefficient = 1.0

    @test all(
        contains.(
            sprint(show, MIME("text/plain"), r1),
            ["r1", "reaction 1", "100.0", "g1 and g2", "glycolysis", "blah", "biocyc"],
        ),
    )

    rlongfor = Reaction(
        "rlongfor",
        Dict(
            m1.id => -1.0,
            m2.id => -1.0,
            m3.id => -1.0,
            m4.id => -1.0,
            m5.id => -1.0,
            m6.id => -1.0,
            m7.id => 1.0,
            m8.id => 1.0,
            m9.id => 1.0,
            m10.id => 1.0,
            m11.id => 1.0,
            m12.id => 1.0,
        ),
        :forward,
    )
    @test contains(sprint(show, MIME("text/plain"), rlongfor), "...")

    rlongrev = Reaction(
        "rlongrev",
        Dict(
            m1.id => -1.0,
            m2.id => -1.0,
            m3.id => -1.0,
            m4.id => -1.0,
            m5.id => -1.0,
            m6.id => -1.0,
            m7.id => 1.0,
            m8.id => 1.0,
            m9.id => 1.0,
            m10.id => 1.0,
            m11.id => 1.0,
            m12.id => 1.0,
        ),
        :reverse,
    )
    @test occursin("...", sprint(show, MIME("text/plain"), rlongrev))

    r2 = Reaction("r2", Dict(m1.id => -2.0, m4.id => 1.0), :reverse)
    @test r2.lb == -1000.0 && r2.ub == 0.0

    r3 = Reaction("r3", Dict(m3.id => -1.0, m4.id => 1.0), :forward)
    @test r3.lb == 0.0 && r3.ub == 1000.0

    r4 = Reaction("r4", Dict(m3.id => -1.0, m4.id => 1.0), :bidirectional)
    r4.annotations = Dict("sboterm" => ["sbo"], "biocyc" => ["ads", "asds"])
    @test r4.lb == -1000.0 && r4.ub == 1000.0

    rd = OrderedDict(r.id => r for r in [r1, r2, r3, r4])
    @test issetequal(["r1", "r4"], ambiguously_identified_items(annotation_index(rd)))

    id = check_duplicate_reaction(r4, rd)
    @test id == "r3"

    r5 = Reaction("r5", Dict(m3.id => -11.0, m4.id => 1.0), :bidirectional)
    id = check_duplicate_reaction(r5, rd)
    @test id == "r3"

    r5 = Reaction("r5", Dict(m3.id => -11.0, m4.id => 1.0), :bidirectional)
    id = check_duplicate_reaction(r5, rd; only_metabolites = false)
    @test isnothing(id)

    r5 = Reaction("r5", Dict(m3.id => -1.0, m4.id => 1.0), :bidirectional)
    id = check_duplicate_reaction(r5, rd; only_metabolites = false)
    @test id == "r3"
end


text:
"""
    modelDimension(q::Int)

Returns the minimal dimension of the model at parameter q.
"""
function modelDimension(q::Int)
    q/2 + 1 |> floor |> Int64
end

"""
    modelInteraction(p::Int,q::Int)

Returns the interaction matrix J of the model at parameter (p,q).
Notice:
- some terms appear with a factor 1/2 with respect to Equation (S8). This is because Equation (S8) is not completely symmetric in the magnatizations;
- ```dot(ms,J,ms)``` is the energy of Equation (S8) divided by ell.
"""
function modelInteraction(p::Int,q::Int)
    q2 = modelDimension(q)
    J = zeros(q2,q2)
    J[1,1] = 1
    for i in 2:q2
        J[1,i] = 2
        J[i,1] = 2

        for j in i:q2
            J[i,j] = 4 * cos( 2 * pi * p / q * (i-1) * (j-1) )
            J[j,i] = 4 * cos( 2 * pi * p / q * (i-1) * (j-1) )
        end
    end

    if q%2==0
        J[q2,q2] = cos(pi*p*q/2)
        J[1,q2] = 1
        J[q2,1] = 1

        for j in 2:q2-1
            J[j,q2] = 2 * cos( pi * p * (j-1) )
            J[q2,j] = 2 * cos( pi * p * (j-1) )
        end
    end

    J /= q
end

"""
    magnetizationEntropy(m::Real)

Returns the entropy per degree of freedom associated to a magnetization m; see Equation (S9)
"""
function magnetizationEntropy(m::Real)
    if m < -1 || m > 1
        throw(ArgumentError(string("m = ", m, " outside of [-1,1] domain!")))
    else
        log(2) - 1/2*xlog(1-m) - 1/2*xlog(1+m)
    end
end

"""
    gradMagnetizationEntropy(m::Real)
    gradMagnetizationEntropy(m::Real; regularizer::Float64)

Returns the derivative of the entropy per degree of freedom associated to a magnetization m; see Equation (S9).
Optional arguments:
-   ```regularizer```: small constant to deal with the logarithmic divergences at m=1 and m=-1. Default value ```1e-3```.
"""
function gradMagnetizationEntropy(m::Real;
                                  regularizer::Real=1e-3)
    if m < -1 || m > 1
        throw(ArgumentError("m outside of [-1,1] domain!"))
    else
        1/2*log(1-m+regularizer) - 1/2*log(1+m+regularizer)
    end
end

"""
    modelEntropy(ms::Vector{T}, q::Int) where T <: Real

Returns the entropy of the model (divided by ell) at magnetizations ms, and at parameter q; see Equation (S11).
"""
function modelEntropy(ms::Vector{T},q::Int) where T <: Real
    q2 = length(ms)
    if q2 == modelDimension(q)
        entropy = 2*magnetizationEntropy(ms[1])
        if q2 - 1 >= 2
            entropy += 4*sum( magnetizationEntropy(ms[i]) for i in 2:q2-1)
        end
        return entropy += ( q%2==0 ? 2 : 4 )*magnetizationEntropy(ms[q2])
    else
        throw(ArgumentError("Length of ms incompatible with q!"))
    end
end

"""
    modelFreeEnergy(ms::Vector{T}, p::Int, q::Int, t::Real) where T <: Real

Returns the free energy of the model (divided by ell) at magnetizations ms, at parameters (p,q) and at temperature t; see Equation (S8) and (S11).
"""
function modelFreeEnergy(ms::Vector{T},p::Int,q::Int,t::Real) where T <: Real
    q2 = length(ms)
    if q2 == modelDimension(q)
        J = modelInteraction(p,q)
        return dot(ms,J,ms) - t * modelEntropy(ms,q)
    else
        throw(ArgumentError("Length of ms incompatible with q!"))
    end
end

"""
    gradModelFreeEnergy(ms::Vector{T}, p::Int, q::Int, t::Real) where T <: Real

Returns the gradient of the free energy of the model (divided by ell) at magnetizations ms, at parameters (p,q) and at temperature t; see Equation (S8) and Equation (S11).
"""
function gradModelFreeEnergy(ms::Vector{T},p::Int,q::Int,t::Real) where T <: Real
    q2 = length(ms)
    if q2 == modelDimension(q)
        J = modelInteraction(p,q)

        gradF = 2*(J*ms)
        gradF[1] += -t * 2 * gradMagnetizationEntropy(ms[1])
        for i in 2:q2-1
            gradF[i] += -t * 4 * gradMagnetizationEntropy(ms[i])
        end
        gradF[q2] += -t * ( q%2==0 ? 2 : 4 ) * gradMagnetizationEntropy(ms[q2])
        return gradF
    else
        throw(ArgumentError("Length of ms incompatible with q!"))
    end
end


text:
@testset "equal.jl" begin
    @testset "EqualConstant()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(2, 6, "x", trailer)

        constraint = SeaPearl.EqualConstant(x, 3, trailer)

        @test constraint in x.onDomainChange
        @test constraint.active.value
    end
    @testset "propagate!(::EqualConstant)" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(2, 6, "x", trailer)
        ax = SeaPearl.IntVarViewMul(x, 3, "3x")

        constraint = SeaPearl.EqualConstant(ax, 6, trailer)

        toPropagate = Set{SeaPearl.Constraint}()
        prunedDomains = SeaPearl.CPModification()

        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test length(ax.domain) == 1
        @test 6 in ax.domain
        @test !(9 in ax.domain)
        @test !(10 in ax.domain)
        @test prunedDomains == SeaPearl.CPModification(
            "3x" => [9, 12, 15, 18],
            "x"  => [3, 4, 5, 6]
        )

        cons2 = SeaPearl.EqualConstant(ax, 9, trailer)

        @test !SeaPearl.propagate!(cons2, toPropagate, prunedDomains)

        @test isempty(ax.domain)

        y = SeaPearl.IntVar(2, 6, "y", trailer)
        constraint1 = SeaPearl.EqualConstant(y, 3, trailer)
        constraint2 = SeaPearl.EqualConstant(y, 4, trailer)

        toPropagate2 = Set{SeaPearl.Constraint}()
        SeaPearl.propagate!(constraint1, toPropagate2, prunedDomains)

        @test constraint2 in toPropagate2

    end

    @testset "pruneEqual!()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(2, 6, "x", trailer)
        y = SeaPearl.IntVar(5, 8, "y", trailer)

        SeaPearl.pruneEqual!(y, x)

        @test length(y.domain) == 2
        @test !(8 in y.domain) && 5 in y.domain && 6 in y.domain

    end

    @testset "propagate!(::Equal)" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(2, 6, "x", trailer)
        y = SeaPearl.IntVar(5, 8, "y", trailer)

        constraint = SeaPearl.Equal(x, y, trailer)
        toPropagate = Set{SeaPearl.Constraint}()
        prunedDomains = SeaPearl.CPModification()

        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test length(x.domain) == 2
        @test length(y.domain) == 2
        @test !(2 in x.domain) && 5 in x.domain && 6 in x.domain
        @test !(8 in y.domain) && 5 in y.domain && 6 in y.domain
        @test prunedDomains == SeaPearl.CPModification("x" => [2, 3, 4],"y" => [7, 8])

        # Propagation test
        z = SeaPearl.IntVar(5, 15, "z", trailer)
        constraint2 = SeaPearl.Equal(y, z, trailer)
        SeaPearl.propagate!(constraint2, toPropagate, prunedDomains)

        # Domain not reduced => not propagation
        @test !(constraint in toPropagate)
        @test !(constraint2 in toPropagate)

        # Domain reduced => propagation
        SeaPearl.remove!(z.domain, 5)
        SeaPearl.propagate!(constraint2, toPropagate, prunedDomains)
        @test constraint in toPropagate
        @test !(constraint2 in toPropagate)

        #Unfeasible test
        t = SeaPearl.IntVar(15, 30, "t", trailer)
        constraint3 = SeaPearl.Equal(z, t, trailer)
        @test !SeaPearl.propagate!(constraint3, toPropagate, prunedDomains)
    end
end


text:
module test_jive_runtests_get_all_files

using Test
using Jive

dir = normpath(@__DIR__, "..") # Jive/test/jive/
skip = String[]
targets = ["//./pipe", "USE_REVISE", "USE_PLOTPANE", "runtests/get_all_files"]
(all_files, start_idx) = Jive.get_all_files(dir, skip, targets)
@test all_files == ["runtests/get_all_files.jl"]

dir = normpath(@__DIR__, ".") # Jive/test/jive/runtests/
skip = String[]
targets = ["//./pipe", "USE_REVISE", "USE_PLOTPANE", "get_all_files"]
(all_files, start_idx) = Jive.get_all_files(dir, skip, targets)
@test all_files == ["get_all_files.jl"]

dir = normpath(@__DIR__, "../../Example/test/example") # Jive/test/Example/test/example
skip = String[]

targets = ["test"]
(all_files, start_idx) = Jive.get_all_files(dir, skip, targets)
@test all_files == ["test1.jl", "test2.jl"]

targets = ["test1"]
(all_files, start_idx) = Jive.get_all_files(dir, skip, targets)
@test all_files == ["test1.jl"]

targets = ["abc"]
(all_files, start_idx) = Jive.get_all_files(dir, skip, targets)
@test isempty(all_files)

targets = ["abc", "test1"]
(all_files, start_idx) = Jive.get_all_files(dir, skip, targets)
@test all_files == ["test1.jl"]

targets = ["abc", "test"]
(all_files, start_idx) = Jive.get_all_files(dir, skip, targets)
@test all_files == ["test1.jl", "test2.jl"]

targets = ["abc", "test1", "test2"]
(all_files, start_idx) = Jive.get_all_files(dir, skip, targets)
@test all_files == ["test1.jl", "test2.jl"]

targets = String[]
(all_files, start_idx) = Jive.get_all_files(dir, skip, targets)
@test all_files == ["test1.jl", "test2.jl"]

targets = [""]
(all_files, start_idx) = Jive.get_all_files(dir, skip, targets)
@test all_files == ["test1.jl", "test2.jl"]

targets = ["."]
(all_files, start_idx) = Jive.get_all_files(dir, skip, targets)
@test all_files == ["test1.jl", "test2.jl"]
@test start_idx == 1

targets = [".", "start=2"]
(all_files, start_idx) = Jive.get_all_files(dir, skip, targets)
@test all_files == ["test1.jl", "test2.jl"]
@test start_idx == 2

end # module test_jive_runtests_get_all_files


text:
function binary_op(op, W, @nospecialize(T))
    ty = LLVM_TYPES[T]
    if isone(W)
        V = T
    else
        ty = "<$W x $ty>"
        V = NTuple{W,VecElement{T}}
    end
    instrs = "%res = $op $ty %0, %1\nret $ty %res"
    call = :($LLVMCALL($instrs, $V, Tuple{$V,$V}, data(v1), data(v2)))
    W > 1 && (call = Expr(:call, :Vec, call))
    Expr(:block, Expr(:meta, :inline), call)
end

# Integer
for (op,f) ∈ [("add",:+),("sub",:-),("mul",:*),("shl",:<<)]
    ff = Symbol('v', op); ff_fast = Symbol(ff, :_fast)
    _ff = Symbol('_', ff)
    _ff_fast = Symbol('_', ff_fast)
    @eval begin
        @generated $ff_fast(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:Integer} = binary_op($op * (T <: Signed ? " nsw" : " nuw"), W, T)
        @generated $ff(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:Integer} = binary_op($op, W, T)

        @generated $_ff_fast(v1::T, v2::T) where {T<:Integer} = binary_op($op * (T <: Signed ? " nsw" : " nuw"), 1, T)
        @inline $ff_fast(v1::T, v2::T) where {T<:Vec} = $_ff_fast(v1, v2)
        @inline $ff_fast(v1::T, v2::T) where {T} = $f(v1, v2)#fallback
        @inline $ff(x::T,y::T) where {T<:IntegerTypesHW} = $_ff_fast(x,y)
    end
end
for (op,f) ∈ [("div",:÷),("rem",:%)]
    ff = Symbol('v', op); _ff = Symbol(:_, ff)
    @eval begin
        @generated $ff(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:Integer} = binary_op((T <: Signed ? 's' : 'u') * $op, W, T)
        @generated $_ff(v1::T, v2::T) where {T<:Integer} = binary_op((T <: Signed ? 's' : 'u') * $op, 1, T)
        @inline $ff(v1::T, v2::T) where {T<:IntegerTypesHW} = $_ff(v1, v2)
    end
end
@inline vcld(x, y) = vadd(vdiv(vsub(x, one(x)), y), one(x))
@inline function vdivrem(x, y)
    d = vdiv(x, y)
    r = vsub(x, vmul(d, y))
    d, r
end
for (op,sub) ∈ [
    ("ashr",:SignedHW),
    ("lshr",:UnsignedHW),
    ("lshr",:IntegerTypesHW),
    ("and",:IntegerTypesHW),
    ("or",:IntegerTypesHW),
    ("xor",:IntegerTypesHW)
]
    ff = sub === :UnsignedHW ? :vashr : Symbol('v', op)
    @eval begin
        @generated $ff(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:$sub}  = binary_op($op, W, T)
        @generated $ff(v1::T, v2::T) where {T<:$sub}  = binary_op($op, 1, T)
    end
end

for (op,f) ∈ [("fadd",:vadd),("fsub",:vsub),("fmul",:vmul),("fdiv",:vfdiv),("frem",:vrem)]
    ff = Symbol(f, :_fast)
    @eval begin
        @generated  $f(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:Union{Float32,Float64}} = binary_op($(op * ' ' * fast_flags(false)), W, T)
        @generated $ff(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:Union{Float32,Float64}} = binary_op($(op * ' ' * fast_flags( true)), W, T)
    end
end
@inline vsub(a::T,b::T) where {T<:Union{Float32,Float64}} = Base.sub_float(a,b)
@inline vadd(a::T,b::T) where {T<:Union{Float32,Float64}} = Base.add_float(a,b)
@inline vsub_fast(a::T,b::T) where {T<:Union{Float32,Float64}} = Base.FastMath.sub_float_fast(a,b)
@inline vadd_fast(a::T,b::T) where {T<:Union{Float32,Float64}} = Base.FastMath.add_float_fast(a,b)

@inline vdiv(v1::AbstractSIMD{W,T}, v2::AbstractSIMD{W,T}) where {W,T<:FloatingTypes} = vfdiv(vsub(v1, vrem(v1, v2)), v2)
@inline vdiv_fast(v1::AbstractSIMD{W,T}, v2::AbstractSIMD{W,T}) where {W,T<:FloatingTypes} = vfdiv_fast(vsub_fast(v1, vrem_fast(v1, v2)), v2)
@inline vdiv_fast(v1::AbstractSIMD{W,T}, v2::AbstractSIMD{W,T}) where {W,T<:IntegerTypesHW} = trunc(T, vfloat_fast(a) / vfloat_fast(b))
@inline function vdiv_fast(v1, v2)
    v3, v4 = promote_div(v1, v2)
    vdiv_fast(v3, v4)
end

@inline vfdiv(a::AbstractSIMDVector{W}, b::AbstractSIMDVector{W}) where {W} = vfdiv(vfloat(a), vfloat(b))
@inline vfdiv_fast(a::AbstractSIMDVector{W}, b::AbstractSIMDVector{W}) where {W} = vfdiv_fast(vfloat_fast(a), vfloat_fast(b))
@inline vfdiv(a, b) = a / b
@inline vfdiv_fast(a, b) = Base.FastMath.div_fast(a, b)

for f ∈ [:vadd,:vadd_fast,:vsub,:vsub_fast,:vmul,:vmul_fast]
    @eval begin
        @inline function $f(a, b)
            c, d = promote(a, b)
            $f(c, d)
        end
    end
end


text:
# An simple example of how Memento can be plugged in as a backend to MicroLogging

# ------------------------------------------------------------------------------
module MementoShim

using MicroLogging
using Memento

# Confusingly, this is the thing which goes by the name "logger" in
# MicroLogging terminology.  Here it's just used as a placeholder to dispatch
# to the global Memento logger hierarchy.
immutable LogShim <: AbstractLogger
end

function MicroLogging.handle_message(::LogShim, ml_level, msg, _module, group, id, filepath, line; kwargs...)
    # TODO: Map the `group` keyword (or something equivalent) into the right logger.
    # For now, assume a module-based logger.
    logger = get_logger(_module)

    if     ml_level == MicroLogging.Debug ; level = "debug"
    elseif ml_level == MicroLogging.Info  ; level = "info"
    elseif ml_level == MicroLogging.Warn  ; level = "warn"
    elseif ml_level == MicroLogging.Error ; level = "error"
    end

    # TODO: Capture all the keyword arguments, and other call site metadata
    rec = logger.record(logger.name, level, logger.levels[level], msg)
    @sync log(logger, rec)
end

end

# ------------------------------------------------------------------------------
# A test module containing some logging
module A

using MicroLogging

function foo()
    @debug "message from micrologging"
    @info  "message from micrologging"
    @warn  "message from micrologging"
    @error "message from micrologging"
end

end

# ------------------------------------------------------------------------------
# Log configuration
using MicroLogging
using Memento

# MicroLogging setup - set Memento as backend
global_logger(MementoShim.LogShim())

# Configure Memento
Memento.config("info")
# Do something special for messages coming from module A
A_logger = get_logger(A)
add_handler(A_logger, DefaultHandler(Syslog(:local0, "julia"), DefaultFormatter("{level} {name}: {msg}")))
set_level(A_logger, "warn")

# Now, do the actual logging
@info "Global log message"
A.foo()


text:
function rectilinear_grid(filename::AbstractString, x::AbstractVector,
                          y::AbstractVector, z::AbstractVector;
                          compress=true, append::Bool=true,
                          extent=nothing)
    Ni, Nj, Nk = length(x), length(y), length(z)
    Npts = Ni*Nj*Nk
    Ncls = num_cells_structured(Ni, Nj, Nk)
    ext = extent_attribute(Ni, Nj, Nk, extent)

    xvtk = XMLDocument()
    vtk = DatasetFile(xvtk, add_extension(filename, ".vtr"), "RectilinearGrid",
                      Npts, Ncls, compress, append)

    # VTKFile node
    xroot = vtk_xml_write_header(vtk)

    # RectilinearGrid node
    xGrid = new_child(xroot, vtk.grid_type)
    set_attribute(xGrid, "WholeExtent", ext)

    # Piece node
    xPiece = new_child(xGrid, "Piece")
    set_attribute(xPiece, "Extent", ext)

    # Coordinates node
    xPoints = new_child(xPiece, "Coordinates")

    # DataArray node
    data_to_xml(vtk, xPoints, x, "x")
    data_to_xml(vtk, xPoints, y, "y")
    data_to_xml(vtk, xPoints, z, "z")

    return vtk::DatasetFile
end

"""
    vtk_grid(filename::AbstractString,
             x::AbstractVector{T}, y::AbstractVector{T}, [z::AbstractVector{T}];
             kwargs...)

Create 2D or 3D rectilinear grid (`.vtr`) file.

Coordinates are specified by separate vectors `x`, `y`, `z`.

# Examples

```jldoctest
julia> vtk = vtk_grid("abc", [0., 0.2, 0.5], collect(-2.:0.2:3), [1., 2.1, 2.3])
VTK file 'abc.vtr' (RectilinearGrid file, open)
```

"""
vtk_grid(filename::AbstractString, x::AbstractVector{T},
         y::AbstractVector{T}, z::AbstractVector{T}; kwargs...) where T =
    rectilinear_grid(filename, x, y, z; kwargs...)

# 2D variant
vtk_grid(filename::AbstractString, x::AbstractVector{T},
         y::AbstractVector{T}; kwargs...) where T =
    rectilinear_grid(filename, x, y, zeros(T, 1); kwargs...)


text:
function oid2string(path::String)
    id = LibGit2.head_oid(LibGit2.GitRepo(path))
    return prod(i -> string(i; base=16, pad=2), id.val)
end

function oid2string(paths::Vector{String})
    function get_oid(path)
        oid = ""
        try
            oid = "-" * oid2string(path)
        catch _
            @warn "The target $path is not a git repository. Commit id for this target will be ignored to name the output."
        end
        return oid
    end

    return prod(get_oid, paths)
end

version2string() = string(Pkg.project().version)

function version2string(paths::Vector{String})
    function get_version(path)
        Pkg.activate(path)
        v = version2string()
        return if isnothing(v)
            @warn "The target $path is not a project folder. Versioning for this target will be ignored to name the output."
            ""
        else
            "-" * v
        end
    end

    save_path = dirname(Pkg.project().path)
    output = prod(get_version, paths)
    Pkg.activate(save_path)

    return output
end

function smart_paths(paths)
    splitted_paths = map(splitpath ∘ normpath, paths)

    @info "debug" paths splitted_paths
    common = paths |> first |> dirname |> splitpath
    for path in splitted_paths
        to_pop = length(common)
        for name in Iterators.zip(common, path)
            name[1] == name[2] || break
            to_pop -= 1
        end
        foreach(_ -> pop!(common), 1:to_pop)
    end

    for path in splitted_paths
        foreach(_ -> popfirst!(path), 1:length(common))
    end

    return joinpath(common...), map(joinpath, splitted_paths)
end


text:
import Base: convert, ==, getproperty, setproperty!, getindex, setindex!, haskey,
             length, iterate, eltype,
             get, get!, getkey, delete!, pop!, keys, values, merge, merge!,
             sizehint!, keytype, valtype

# AbstractSynapse is used for all PyObjects that are wrapped in a Julia type
abstract type AbstractSynapse end
abstract type AbstractEntity <: AbstractSynapse end
abstract type AbstractSynapseDict <: AbstractDict{Any,Any} end # Would have preferred multiple inheritance. Use traits?

# utility function for making sure AbstractSynapse objects are passed as PyObjects to python
unwrap(a::Any) = a # fallback
unwrap(a::AbstractSynapse) = getfield(a,:po) #a.po
unwrap(a::AbstractSynapseDict) = getfield(a,:po) #a.po

# # unwrap for keyword arguments
# unwrappair(t::Pair{Symbol,Any}) = t # fallback
# unwrappair(t::Pair{Symbol,AbstractSynapse}) = (t[1],t[2].po)
# unwrappair(t::Pair{Symbol,AbstractSynapseDict}) = (t[1],t[2].po)

# wrap is more or less the inverse of unwrap()
wrap(x::Any) = x # fallback
function wrap(po::PyObject)
	if pyisinstance(po, synapseclient.Entity)
		pyisinstance(po, synapseclient.Project) && return Project(po)
		pyisinstance(po, synapseclient.Folder)  && return Folder(po)
		pyisinstance(po, synapseclient.File)    && return File(po)
		pyisinstance(po, synapseclient.Link)    && return Link(po)
		pyisinstance(po, synapseclient.Schema)  && return Schema(po)
		return Entity(po)
	elseif pyisinstance(po, pydict_object.DictObject)
		pyisinstance(po, synapseclient.Evaluation) && return Evaluation(po)
		pyisinstance(po, synapseclient.Wiki)       && return Wiki(po)
		pyisinstance(po, synapseclient.Team)       && return Team(po)
		return DictObject(po) # IMPORTANT - this will make sure the internal dictionaries are modified in-place and not copied
	end
	# Activity and Annotations inherit from dict, not DictObject
	pyisinstance(po, synapseclient.Activity)    && return Activity(po)
	pyisinstance(po, pyannotations.Annotations) && return annotations.Annotations(po)
	pyisinstance(po, pyclient.Synapse)          && return Synapse(po)
	pyisinstance(po, pycache.Cache)             && return cache.Cache(po)

	# IMPORTANT: do not convert python dicts because that will make a copy
	pyisinstance(po, pybuiltin(:dict)) && return PyDict(po)

	convert(pytype_query(po),po) # fallback to default python -> julia conversion
end

# return type conversion is done by wrap, since PyCall conversion is a bit too eager to convert classes to pure dicts
# function synapsecall(obj::PyObject,method::Symbol,args...;kwargs...)
# 	wrap(pycall(obj[method], PyObject,
# 	                        map(unwrap,args)...;
# 	                        map(unwrappair,kwargs)...))
# end
function synapsecall(func::PyObject,args...;kwargs...)
	# wrap(pycall(func, PyObject,
	#                         map(unwrap,args)...;
	#                         map(unwrappair,kwargs)...))
	wrap(pycall(func, PyObject,
	                        unwrap.(args)...;
	                        (k=>unwrap(v) for (k,v) in kwargs)... ))
end
synapsecall(obj::PyObject,method::Symbol,args...;kwargs...) = synapsecall(getproperty(obj,method),args...;kwargs...)
synapsecall(obj::AbstractSynapse,method::Symbol,args...;kwargs...) = synapsecall(unwrap(obj),method,args...;kwargs...)
# synapsecall(obj::AbstractSynapseDict,method::Symbol,args...;kwargs...) = synapsecall(unwrap(obj),method,args...;kwargs...)
synapsecall(obj::AbstractSynapseDict,method::Symbol,args...;kwargs...) = synapsecall(convert(PyObject,obj),method,args...;kwargs...)
#utilcall(method::Symbol,args...;kwargs...) = synapsecall(synapseclient.utils,method,args...;kwargs...)

lowercasesymbol(s::Symbol) = Symbol(lowercase(string(s)))

macro synapsefunction(T, name)
	x = Expr(:quote,name) # x is now :(:val) where val is $name
	name = lowercasesymbol(name)
	esc(:(  $name(t::$T,args...;kwargs...) = synapsecall(t,$x,args...;kwargs...)  ))
end

macro standalonefunction(parent, name)
	x = Expr(:quote,name) # x is now :(:val) where val is $name
	name = lowercasesymbol(name)
	esc(:(  $name(args...;kwargs...) = synapsecall($parent,$x,args...;kwargs...)  ))

end
macro utilfunction(name)
	esc(:( @standalonefunction(synapseclient.utils, $name) ))
end
macro entityfunction(name)
	esc(:( @standalonefunction(synapseclient.entity, $name) ))
end
macro annotationsfunction(name)
	esc(:( @standalonefunction(synapseclient.annotations, $name) ))
end
macro cachefunction(name)
	esc(:( @standalonefunction(synapseclient.cache, $name) ))
end

macro createtype(name, super, wrappedClass, storageClass, doAssert)
	# assert = doAssert ? :(@assert pytypeof(po)==$wrappedClass) : (:;) # require exact match
	assert = doAssert ? :(@assert pyisinstance(po,$wrappedClass)) : (:;) # require only that it is a subtype

	esc(quote
		struct $name <: $super
			po::$storageClass
			function $name(po::Union{PyObject,$storageClass})
				$assert
				new(po)
			end
		end
		$name(x::$name) = $name(x.po)
		#$name(args...;kwargs...) = $name(pycall($wrappedClass,PyObject,args...;kwargs...))
		$name(args...;kwargs...) = synapsecall($wrappedClass,args...;kwargs...)
	end)
end

macro entitytype(name)
	esc(:( @createtype($name,AbstractEntity,synapseclient.$name,PyObject,true) ))
end
macro dicttype(name, args...)
	pyType = isempty(args) ? (:(synapseclient.$name)) : args[1]
	esc(:( @createtype($name,AbstractSynapseDict,$pyType,PyDict,true) ))
end

convert(::Type{PyObject}, a::AbstractSynapse) = unwrap(a)
convert(::Type{PyObject}, a::AbstractSynapseDict) = convert(PyObject,unwrap(a))
==(x::AbstractSynapse, y::AbstractSynapse) = unwrap(x)==unwrap(y)
==(x::AbstractSynapseDict, y::AbstractSynapseDict) = unwrap(x)==unwrap(y)

# old
# getindex(e::AbstractSynapse, key::AbstractString)              = wrap( unwrap(e)[key] )
# setindex!(e::AbstractSynapse, value, key::AbstractString)      = unwrap(e)[key] = value
# # getindex(a::AbstractSynapseDict, key::AbstractString)         = wrap( unwrap(a)[key] )
# getindex(a::AbstractSynapseDict, key::AbstractString)         = wrap( get(unwrap(a).o,Any,key) )
# setindex!(a::AbstractSynapseDict, value, key::AbstractString) = unwrap(a)[key] = value

getproperty(e::AbstractSynapse, key::Symbol)         = wrap( getproperty(unwrap(e),key) )
setproperty!(e::AbstractSynapse, key::Symbol, value) = setproperty!(unwrap(e),key,value)
getproperty(a::AbstractSynapseDict, key::Symbol)         = wrap( getproperty(unwrap(a).o,key) )
setproperty!(a::AbstractSynapseDict, key::Symbol, value) = setproperty!(unwrap(a).o,key,value)

getindex(a::AbstractSynapseDict, key::AbstractString)         = wrap(unwrap(a)[key])
setindex!(a::AbstractSynapseDict, value, key::AbstractString) = unwrap(a)[key] = value

pyhasattr = PyNULL() # pybuiltin(:hasattr)
hasattr(e::AbstractSynapse, key)      = synapsecall(pyhasattr, e, key)
# hasattr(d::AbstractSynapseDict, key) = haskey(unwrap(d), key)
hasattr(d::AbstractSynapseDict, key) = synapsecall(pyhasattr, d, key)
haskey(d::AbstractSynapseDict, key) = hasattr(d,key)

length(a::AbstractSynapseDict) = length(unwrap(a))
iterate(a::AbstractSynapseDict) = iterate(unwrap(a))
iterate(a::AbstractSynapseDict,state) = iterate(unwrap(a),state)
eltype(::Type{AbstractSynapseDict}) = Pair{Any,Any}#eltype(PyDict) # PyDict bug?

get(a::AbstractSynapseDict, key, default) = get(unwrap(a),key,default)
get(f::Function, a::AbstractSynapseDict, key) = get(f,unwrap(a),key)
get!(a::AbstractSynapseDict, key, default) = get!(unwrap(a),key,default)
get!(f::Function, a::AbstractSynapseDict, key) = get!(f,unwrap(a),key)
getkey(a::AbstractSynapseDict, key, default) = getkey(unwrap(a),key,default)
delete!(a::AbstractSynapseDict, key) = delete!(unwrap(a),key)
pop!(a::AbstractSynapseDict, key) = pop!(unwrap(a),key)
pop!(a::AbstractSynapseDict, key, default) = pop!(unwrap(a),key,default)
keys(a::AbstractSynapseDict) = keys(unwrap(a))
values(a::AbstractSynapseDict) = values(unwrap(a))
merge(a::T, others...) where {T<:AbstractSynapseDict} = T(merge(unwrap(a),others))
merge!(a::T, others...) where {T<:AbstractSynapseDict} = merge!(unwrap(a),others)
sizehint!(a::AbstractSynapseDict, n) = sizehint!(unwrap(a),n)
keytype(a::AbstractSynapseDict) = keytype(unwrap(a))
valtype(a::AbstractSynapseDict) = valtype(unwrap(a))


text:
using Test
using Polyhedra
using JuMP
const MOIT = MOI.Test
const MOIB = MOI.Bridges

@testset "Continuous Linear problems with VRepOptimizer" begin
    optimizer = VRepOptimizer{Float64}()
    @testset "SolverName" begin
        @test MOI.get(optimizer, MOI.SolverName()) == "VRep"
    end
    cache = MOIU.UniversalFallback(Polyhedra._MOIModel{Float64}())
    cached = MOIU.CachingOptimizer(cache, optimizer)
    bridged = MOIB.full_bridge_optimizer(cached, Float64)
    config = MOIT.TestConfig(duals=false)
    MOIT.contlineartest(bridged, config,
                        # linear8a and linear12 will be solved by https://github.com/JuliaOpt/MathOptInterface.jl/pull/702
                        ["linear8a", "linear12", "partial_start"])
end


text:
module ExaTron

using Libdl
using LinearAlgebra
using CUDA
using CUDA.CUBLAS

using MPI

using Printf
using DelimitedFiles

export dtron, solveProblem, createProblem, setOption, getOption, ExaTronProblem

const BLAS_LIBRARY = :Tron
const EXATRON_LIBRARY = "libtron"

has_c_library() = !isnothing(Libdl.dlopen(EXATRON_LIBRARY; throw_error=false))
tron_zeros(S, n) = fill!(S(undef, Int64(n)), zero(eltype(S)))
tron_zeros(S, dims::Tuple) = fill!(S(undef, Int64(dims[1]), Int64(dims[1])), zero(eltype(S)))

include("daxpy.jl")
include("dcopy.jl")
include("ddot.jl")
include("dmid.jl")
include("dnrm2.jl")
include("dgpnorm.jl")
include("dscal.jl")
include("dsel2.jl")
include("dssyax.jl")
include("ihsort.jl")
include("insort.jl")
include("dnsol.jl")
include("dtsol.jl")
include("dtrqsol.jl")
include("dbreakpt.jl")
include("dgpstep.jl")
include("dicf.jl")
include("dicfs.jl")
include("dprsrch.jl")
include("dcauchy.jl")
include("dtrpcg.jl")
include("dspcg.jl")
include("dtron.jl")
include("driver.jl")

include("admm/opfdata.jl")
include("admm/environment.jl")
include("admm/proxal_wrapper.jl")
include("admm/bus_kernel.jl")
include("admm/generator_kernel.jl")
include("admm/generator_kernel_proxal.jl")
include("admm/polar_kernel.jl")
include("admm/tron_kernel.jl")
include("admm/eval_kernel.jl")
include("admm/auglag_kernel.jl")
include("admm/acopf_admm_gpu.jl")
include("admm/acopf_admm_gpu_two_level.jl")

end # module