text:
#File which defined all the algorithms. Each algorithm takes in an ABCtype

"""
    runabc(ABCsetup::ABCtype, targetdata; progress = false, verbose = false, parallel = true)

Run ABC with ABCsetup defining the algorithm and inputs to algorithm, targetdata is the data we wish to fit the model to and will be used as an input for the simulation function defined in ABCsetup. If progress is set to `true` a progress meter will be shown. Inference will be run in parallel via multithreading if `parallel = true`. The environmental variable JULIA_NUM_THREADS needs to be set prior to launching a julia session.
"""
function runabc(ABCsetup::ABCRejection, targetdata; progress = false, verbose = false, parallel = false)

  #initalize array of particles
  particlesall = Array{ParticleRejection}(undef, ABCsetup.maxiterations)

  if progress
    p = Progress(ABCsetup.nparticles, 1, "Running ABC rejection algorithm...", 30)
  end

  if parallel
    Printf.@printf("Preparing to run in parallel on %i processors\n", nthreads())

    particles = Array{ParticleRejection}(undef, ABCsetup.maxiterations)
    distvec = zeros(Float64, ABCsetup.maxiterations) #store distances in an array
    i = Atomic{Int64}(0)
    cntr = Atomic{Int64}(0)
    @threads for its = 1:ABCsetup.maxiterations

      if i[] > ABCsetup.nparticles
        break
      end

      #get new proposal parameters
      newparams = getproposal(ABCsetup.prior, ABCsetup.nparams)
      #simulate with new parameters
      dist, out = ABCsetup.simfunc(newparams, ABCsetup.constants, targetdata)
      #keep track of all particles incase we don't reach nparticles with dist < ϵ
      particlesall[its] = ParticleRejection(newparams, dist, out)

      #if simulated data is less than target tolerance accept particle
      if dist < ABCsetup.ϵ
        particles[its] = ParticleRejection(newparams, dist, out)
        distvec[its] = dist
        atomic_add!(i, 1)
      end
      atomic_add!(cntr,1)

    end
    # Remove particles that are still #undef and corresponding distances
    idx = [isassigned(particles,ii) for ii in eachindex(particles)]
    particles = particles[idx]
    distvec = distvec[idx]
    i = length(particles)    # Number of accepted particles
    its = cntr[]    # Total number of simulations

  else
    Printf.@printf("Preparing to run in serial on %i processor\n", 1)

    particles = Array{ParticleRejection}(undef, ABCsetup.nparticles)
    distvec = zeros(Float64, ABCsetup.nparticles) #store distances in an array
    i = 1 #set particle indicator to 1
    its = 0 #keep track of number of iterations
    while (i < (ABCsetup.nparticles + 1)) & (its < ABCsetup.maxiterations)

      its += 1
      #get new proposal parameters
      newparams = getproposal(ABCsetup.prior, ABCsetup.nparams)
      #simulate with new parameters
      dist, out = ABCsetup.simfunc(newparams, ABCsetup.constants, targetdata)
      #keep track of all particles incase we don't reach nparticles with dist < ϵ
      particlesall[its] = ParticleRejection(newparams, dist, out)

      #if simulated data is less than target tolerance accept particle
      if dist < ABCsetup.ϵ
        particles[i] = ParticleRejection(newparams, dist, out)
        distvec[i] = dist
        i +=1
        if progress
          next!(p)
        end
      end
    end
    i -= 1    # Correct to total number of particels
  end

  if i < ABCsetup.nparticles
    @warn "Only accepted $(i) particles with ϵ < $(ABCsetup.ϵ). \n\tYou may want to increase ϵ or increase maxiterations. \n\t Resorting to taking the $(ABCsetup.nparticles) particles with smallest distance"
    d = map(p -> p.distance, particlesall)
    particles = particlesall[sortperm(d)[1:ABCsetup.nparticles]]
    distvec = map(p -> p.distance, particles)
  elseif i > ABCsetup.nparticles
    particles = particles[1:ABCsetup.nparticles]
    distvec = distvec[1:ABCsetup.nparticles]
  end

  out = ABCrejectionresults(particles, its, ABCsetup, distvec)
  return out
end

function runabc(ABCsetup::ABCRejectionModel, targetdata; progress = false, verbose = false)

  ABCsetup.nmodels > 1 || error("Only 1 model specified, use ABCRejection method to estimate parameters for a single model")

  #initalize array of particles
  particles = Array{ParticleRejectionModel}(undef,
                ABCsetup.Models[1].nparticles)

  i = 1 #set particle indicator to 1
  its = 0 #keep track of number of iterations
  distvec = zeros(Float64, ABCsetup.Models[1].nparticles) #store distances in an array

  if progress
    p = Progress(ABCsetup.Models[1].nparticles, 1, "Running ABC rejection algorithm...", 30)
  end

  while (i < (ABCsetup.Models[1].nparticles + 1)) & (its < ABCsetup.Models[1].maxiterations)

    its += 1
    #sample uniformly from models
    model = rand(1:ABCsetup.nmodels)
    #get new proposal parameters
    newparams = getproposal(ABCsetup.Models[model].prior, ABCsetup.Models[model].nparams)
    #simulate with new parameters
    dist, out = ABCsetup.Models[model].simfunc(newparams, ABCsetup.Models[model].constants, targetdata)

    #if simulated data is less than target tolerance accept particle
    if dist < ABCsetup.Models[1].ϵ
      particles[i] = ParticleRejectionModel(newparams, model, dist, out)
      distvec[i] = dist
      i +=1
      if progress
        next!(p)
      end
    end
  end

  i > ABCsetup.Models[1].nparticles || error("Only accepted $(i) particles with ϵ < $(ABCsetup.Models[1].ϵ). \n\tDecrease ϵ or increase maxiterations ")

  out = ABCrejectionmodelresults(particles, its, ABCsetup, distvec)
  return out
end

function runabc(ABCsetup::ABCSMC, targetdata; verbose = false, progress = false, parallel = false)

  #run first population with parameters sampled from prior
  if verbose
    println("##################################################")
    println("Use ABC rejection to get first population")
  end
  ABCrejresults = runabc(ABCRejection(ABCsetup.simfunc, ABCsetup.nparams,
                  ABCsetup.ϵ1, ABCsetup.prior; nparticles = ABCsetup.nparticles,
                  maxiterations = ABCsetup.maxiterations,
                  constants = ABCsetup.constants), targetdata,
                  progress = progress, parallel = parallel);

  oldparticles, weights = setupSMCparticles(ABCrejresults, ABCsetup)
  ABCsetup.kernel.kernel_parameters = (maximum(ABCrejresults.parameters, dims = 1) - minimum(ABCrejresults.parameters, dims = 1) ./2)[:]
  ϵ = quantile(ABCrejresults.dist, ABCsetup.α) # set new ϵ to αth quantile
  ϵvec = [ϵ] #store epsilon values
  numsims = [ABCrejresults.numsims] #keep track of number of simualtions
  particles = Array{ParticleSMC}(undef, ABCsetup.nparticles) #define particles array

  if verbose
    println("Running ABC SMC... \n")
  end

  if parallel
    Printf.@printf("Preparing to run in parallel on %i processors\n", nthreads())
  else
    Printf.@printf("Preparing to run in serial on %i processor\n", 1)
  end

  popnum = 1
  finalpop = false

  if sum(ABCrejresults.dist .< ABCsetup.ϵT) == ABCsetup.nparticles
      @warn "Target ϵ reached with ABCRejection algorithm, no need to use ABC SMC algorithm, returning ABCRejection output..."
      return ABCrejresults
  end

  while (ϵ > ABCsetup.ϵT) & (sum(numsims) < ABCsetup.maxiterations)

    if progress
      p = Progress(ABCsetup.nparticles, 1, "ABC SMC population $(popnum), new ϵ: $(round(ϵ, digits = 2))...", 30)
    end

    if parallel

      # Arrays initialised with length maxiterations to enable use of unique index ii
      particles = Array{ParticleSMC}(undef, ABCsetup.maxiterations)
      distvec = zeros(Float64, ABCsetup.maxiterations)
      i = Atomic{Int64}(0)
      its = Atomic{Int64}(0)

      @threads for ii = 1:ABCsetup.maxiterations

        if i[] > ABCsetup.nparticles
          break
        end

        j = wsample(1:ABCsetup.nparticles, weights)
        particle = oldparticles[j]
        newparticle = perturbparticle(particle, ABCsetup.kernel)
        priorp = priorprob(newparticle.params, ABCsetup.prior)
        if priorp == 0.0 #return to beginning of loop if prior probability is 0
          continue
        end

        #simulate with new parameters
        dist, out = ABCsetup.simfunc(newparticle.params, ABCsetup.constants, targetdata)

        #if simulated data is less than target tolerance accept particle
        if dist < ϵ
          particles[ii] = newparticle
          distvec[ii] = dist
          particles[ii].other = out
          particles[ii].distance = dist
          atomic_add!(i, 1)
        end

        atomic_add!(its,1)
      end

      # Remove particles that are still #undef and corresponding distances
      idx = [isassigned(particles,ii) for ii in eachindex(particles)]
      particles = particles[idx][1:ABCsetup.nparticles]
      distvec = distvec[idx][1:ABCsetup.nparticles]
      its = its[]

    else
      particles = Array{ParticleSMC}(undef, ABCsetup.nparticles)
      distvec = zeros(Float64, ABCsetup.nparticles)
      i = 1 #set particle indicator to 1
      its = 1

      while i < ABCsetup.nparticles + 1

        j = wsample(1:ABCsetup.nparticles, weights)
        particle = oldparticles[j]
        newparticle = perturbparticle(particle, ABCsetup.kernel)
        priorp = priorprob(newparticle.params, ABCsetup.prior)
        if priorp == 0.0 #return to beginning of loop if prior probability is 0
          continue
        end

        #simulate with new parameters
        dist, out = ABCsetup.simfunc(newparticle.params, ABCsetup.constants, targetdata)

        #if simulated data is less than target tolerance accept particle
        if dist < ϵ
          particles[i] = newparticle
          particles[i].other = out
          particles[i].distance = dist
          distvec[i] = dist
          i += 1
          if progress
            next!(p)
          end
        end
        its += 1
      end
    end

    particles, weights = smcweights(particles, oldparticles, ABCsetup.prior, ABCsetup.kernel)
    ABCsetup.kernel.kernel_parameters = ABCsetup.kernel.calculate_kernel_parameters(particles)
    oldparticles = particles

    if finalpop == true
      break
    end

    ϵ = quantile(distvec, ABCsetup.α)

    if ϵ < ABCsetup.ϵT
      ϵ = ABCsetup.ϵT
      push!(ϵvec, ϵ)
      push!(numsims, its)
      popnum = popnum + 1
      finalpop = true
      continue
    end

    push!(ϵvec, ϵ)
    push!(numsims, its)

    if ((( abs(ϵvec[end - 1] - ϵ )) / ϵvec[end - 1]) < ABCsetup.convergence) == true
      if verbose
        println("\nNew ϵ is within $(round(ABCsetup.convergence * 100, digits=2))% of previous population, stop ABC SMC\n")
      end
      break
    end

    popnum = popnum + 1

    if verbose
      println("##################################################")
      show(ABCSMCresults(particles, numsims, ABCsetup, ϵvec))
      println("##################################################\n")
    end
  end
  if sum(numsims) >= ABCsetup.maxiterations
    if verbose
      println("\nPassed maxiterations=$(ABCsetup.maxiterations), stop ABC SMC\n")
    end
  end

  out = ABCSMCresults(particles, numsims, ABCsetup, ϵvec)
  return out
end

"""
    runabc(ABCsetup::ABCtype, targetdata; progress = false, verbose = false)

When the SMC algorithms are used, a print out at the end of each population will be made if verbose = true.
"""
function runabc(ABCsetup::ABCSMCModel, targetdata; verbose = false, progress = false)

  ABCsetup.nmodels > 1 || error("Only 1 model specified, use ABCSMC method to estimate parameters for a single model")

  #run first population with parameters sampled from prior
  if verbose
    println("##################################################")
    println("Use ABC rejection to get first population")
  end
  ABCrejresults = runabc(ABCRejectionModel(
            map(x -> x.simfunc, ABCsetup.Models),
            map(x -> x.nparams, ABCsetup.Models),
            ABCsetup.Models[1].ϵ1,
            map(x -> x.prior, ABCsetup.Models),
            constants = map(x -> x.constants, ABCsetup.Models),
            nparticles = ABCsetup.Models[1].nparticles,
            maxiterations = ABCsetup.Models[1].maxiterations),
            targetdata);

  oldparticles, weights = setupSMCparticles(ABCrejresults, ABCsetup)
  ϵ = quantile(ABCrejresults.dist, ABCsetup.α) # set new ϵ to αth quantile
  ϵvec = [ϵ] #store epsilon values
  numsims = [ABCrejresults.numsims] #keep track of number of simulations
  particles = Array{ParticleSMCModel}(undef, ABCsetup.nparticles) #define particles array
  weights, modelprob = getparticleweights(oldparticles, ABCsetup)
  ABCsetup = modelselection_kernel(ABCsetup, oldparticles)

  modelprob = ABCrejresults.modelfreq

  if verbose
    println("Run ABC SMC \n")
  end

  popnum = 1
  finalpop = false

  if verbose
    println(ABCSMCmodelresults(oldparticles, numsims, ABCsetup, ϵvec))
  end

  if sum(ABCrejresults.dist .< ABCsetup.ϵT) == ABCsetup.nparticles
      @warn "Target ϵ reached with ABCRejection algorithm, no need to use ABC SMC algorithm, returning ABCRejection output..."
      return ABCrejresults
  end

  while (ϵ >= ABCsetup.ϵT) & (sum(numsims) <= ABCsetup.maxiterations)

    i = 1 #set particle indicator to 1
    particles = Array{ParticleSMCModel}(undef, ABCsetup.nparticles)
    distvec = zeros(Float64, ABCsetup.nparticles)
    its = 1

    if progress
      p = Progress(ABCsetup.nparticles, 1, "ABC SMC population $(popnum), new ϵ: $(round(ϵ, digits = 2))...", 30)
    end
    while i < ABCsetup.nparticles + 1

      #draw model from previous model probabilities
      mstar = wsample(1:ABCsetup.nmodels, modelprob)
      #perturb model
      mdoublestar = perturbmodel(ABCsetup, mstar, modelprob)
      # sample particle with correct model
      j = wsample(1:ABCsetup.nparticles, weights[mdoublestar, :])
      particletemp = oldparticles[j]
      #perturb particle
      newparticle = perturbparticle(particletemp, ABCsetup.Models[mdoublestar].kernel)
      #calculate priorprob
      priorp = priorprob(newparticle.params, ABCsetup.Models[mdoublestar].prior)

      if priorp == 0.0 #return to beginning of loop if prior probability is 0
        continue
      end

      #simulate with new parameters
      dist, out = ABCsetup.Models[mdoublestar].simfunc(newparticle.params, ABCsetup.Models[mdoublestar].constants, targetdata)

      #if simulated data is less than target tolerance accept particle
      if dist < ϵ
        particles[i] = newparticle
        particles[i].other = out
        particles[i].distance = dist
        distvec[i] = dist
        i += 1
        if progress
          next!(p)
        end
      end

      its += 1
    end

    particles, weights = smcweightsmodel(particles, oldparticles, ABCsetup, modelprob)
    weights, modelprob = getparticleweights(particles, ABCsetup)
    ABCsetup = modelselection_kernel(ABCsetup, particles)
    oldparticles = particles

    if finalpop == true
      break
    end

    ϵ = quantile(distvec, ABCsetup.α)

    if ϵ < ABCsetup.ϵT
      ϵ = ABCsetup.ϵT
      push!(ϵvec, ϵ)
      push!(numsims, its)
      popnum = popnum + 1
      finalpop = true
      continue
    end

    push!(ϵvec, ϵ)
    push!(numsims, its)

    if verbose
      println("##################################################")
      println(ABCSMCmodelresults(particles, numsims, ABCsetup, ϵvec))
      println("##################################################\n")
    end

    if ((( abs(ϵvec[end - 1] - ϵ )) / ϵvec[end - 1]) < ABCsetup.convergence) == true
      println("New ϵ is within $(round(ABCsetup.convergence * 100, digits = 2))% of previous population, stop ABC SMC")
      break
    end

    popnum = popnum + 1
  end

  out = ABCSMCmodelresults(particles, numsims, ABCsetup, ϵvec)
  return out
end


text:
module TestSpecificityFunctions
  using Test
  using EcologicalNetworks
  using LinearAlgebra

  # specificity

  @test EcologicalNetworks.pdi(vec([1.0 0.0 0.0])) ≈ 1.0
  @test EcologicalNetworks.pdi(vec([0.0 1.0 0.0])) ≈ 1.0
  @test EcologicalNetworks.pdi(vec([0.0 0.2 0.0])) ≈ 1.0

  N = BipartiteNetwork(Matrix(I, (10,10)))
  @test specificity(N)[species(N; dims=1)[1]] ≈ 1.0
  N = BipartiteQuantitativeNetwork(Matrix{Float64}(I, (10,10)))
  @test specificity(N)[species(N; dims=1)[1]] ≈ 1.0

end

ext: jl


text:
export create_agent_template

template_model = mt"""
using Agents
using AgentsX
using Random
using InteractiveDynamics
using CairoMakie

include("agent.jl")

\"\"\"
    initialise()
Initialise function for the ABM.
This function will return the ABM.
Arguments necessary for declaring the space will need to be passed manually.
\"\"\"

function initialise(step_function; numagents = {{:num_agents}}, seed = {{:random_seed}} )
    space = {{:space}}
    properties = Dict(:step_function => step_function)
    rng = Random.MersenneTwister(seed)

    model = ABM(
        {{:agent_class}}, space;
        properties, rng
    )

    for n in 1:numagents
        agent = {{:agent_class}}(
            #TODO Complete Agent Initialisation
        )
        add_agent_single!(agent, model)
    end

    return model

end

#TODO reorder agent actions in expected execution order if required.
#TODO insert keyword arguments if required
model = initialise([{{#:actions}}{{.}}{{^.[end]}}, {{/.[end]}}{{/:actions}}])
"""

template_agent = mt"""
using Agents
using AgentsX

include("actions.jl")
{{#:layers}}

\"\"\"
    {{.}} <:ParamLayer
The {{.}} layer for {{:agent_class}}.
\"\"\"
mutable struct {{.}} <: ParamLayer
    #TODO enter parameters for the {{.}} Layer
end

{{/:layers}}

\"\"\"
    {{:agent_class}} <:AbstractAgent
A {{:agent_class}} in a {{:space}}.
The agents by default possess the following layers:
{{#:layers}}
    - {{.}}
{{/:layers}}
Spatial functions are handled by the respective space of Agents.jl
\"\"\"

mutable struct {{:agent_class}} <:AbstractAgent
    id::Int
    #TODO customise 'pos' as necessary
    #pos::NTuple{2, Int}
    {{#:layers}}
    {{.}}::{{.}}
    {{/:layers}}
end
"""

template_action = mt"""
{{#:actions}}
\"\"\"
{{.}}(agent, model)
TODO <Enter action description>
\"\"\"
function {{.}}(agent, model)

    #INPUT <Layer>

    #TODO complete {{.}} action code
    #TODO enter input/output layer information

    #OUTPUT <Layer>
end
{{/:actions}}
"""

function create_agent_template(;agent_class, num_agents, space, random_seed = 250, layers = Nothing, actions = Nothing, path = "src")
    if layers == Nothing
        layers = ["SocioEconomic", "Cognitive", "Physiological"]
    end
    dict = Dict(:agent_class => agent_class, :num_agents => num_agents, :space => space, :random_seed => random_seed, :layers => layers, :actions => actions)
    open(path * "/model.jl", "w") do io
        write(io, render(template_model, dict))
    end

    open(path * "/agent.jl", "w") do io
        write(io, render(template_agent, dict))
    end

    open(path * "/actions.jl", "w") do io
        write(io, render(template_action, dict))
    end

end


text:
# functions related to student's T distribution

# R implementations
# For pdf and logpdf we use the Julia implementation
using .RFunctions:
    tdistcdf,
    tdistccdf,
    tdistlogcdf,
    tdistlogccdf,
    tdistinvcdf,
    tdistinvccdf,
    tdistinvlogcdf,
    tdistinvlogccdf

# Julia implementations
tdistpdf(ν::Real, x::Real) = exp(tdistlogpdf(ν, x))

tdistlogpdf(ν::Real, x::Real) = tdistlogpdf(promote(ν, x)...)
function tdistlogpdf(ν::T, x::T) where {T<:Real}
    isinf(ν) && return normlogpdf(x)
    νp12 = (ν + 1) / 2
    return loggamma(νp12) - (logπ + log(ν)) / 2 - loggamma(ν / 2) - νp12 * log1p(x^2 / ν)
end


text:
## Entry point of the executable

using ArgParse

function parse_error(msg)
    println(stderr, msg)
    println(stderr, "See --help for usage.")
    return 5
end

function invalid_input_error(msg, e, stacktrace)
    println(stderr, msg)
    Base.display_error(stderr, e, stacktrace)
    return 4
end

function internal_error(msg, e, stacktrace)
    println(stderr, msg)
    Base.display_error(stderr, e, stacktrace)
    return 3
end

function unhandled_error(msg, e, stacktrace)
    println(stderr, msg)
    Base.display_error(stderr, e, stacktrace)
    return 2
end

function parse_commandline()
    s = ArgParseSettings(prog = "CrystalNets" * (Sys.iswindows() ? ".exe" : ""),
                         description = "Automatic reckognition of crystal net topologies.",
                         epilog = """\n\n\n\nCLUSTERING_MODE options:\n\n
                         \ua0\ua0\ua0* input: use the input residues as clusters. Fail if some atom does not belong to a residue.\n\n
                         \ua0\ua0\ua0* atom: each atom is its own cluster.\n\n
                         \ua0\ua0\ua0* mof: discard the input residues and consider the input as a MOF. Identify organic and inorganic clusters using a simple heuristic based on the atom types.\n\n
                         \ua0\ua0\ua0* guess: discard the input residues and try clustering mode "mof". If it fails, fall back to "atom".\n\n
                         \ua0\ua0\ua0* auto: if the input assigns each atom to a residue, equivalent to "input". Otherwise, fall back to "guess". Default option.\n\n
                         \n\n
                         BONDING_MODE options:\n\n
                         \ua0\ua0\ua0* input: use the bonds explicitly given by the input file. Fail if bonds are not provided by the input.\n\n
                         \ua0\ua0\ua0* chemfiles: guess bonds using chemfiles built-in bond detection mechanism.
                         \ua0\ua0\ua0* auto: if the input possesses explicit bonds, use them unless they are suspicious. Otherwise, fall back to "chemfiles". Default option.\n\n
                         \n\n
                         CREATE_MODE options:\n\n
                         \ua0\ua0\ua0* empty: empty archive, unable to reckognize any topological structure.\n\n
                         \ua0\ua0\ua0* rcsr: RCSR Systre archive.\n\n
                         \ua0\ua0\ua0* zeolites: zeolites topologies from the database of zeolite structures.\n\n
                         \ua0\ua0\ua0* full: combined rcsr and zeolites archives. Default option.\n\n
                         """,
                         preformatted_epilog = false,
                         autofix_names = true,
                         add_help = false,
                         usage = """
                   usage: CrystalNets [-a ARCHIVE_PATH [-u NAME [-f] | -r [-f]]]
                                      [[-c CLUSTERING_MODE] [-b BONDING_MODE]] | -g]
                                      [--no-export | -e DIR_PATH] CRYSTAL_FILE   (Form A)
                          CrystalNets -a ARCHIVE_PATH -n [CREATE_MODE] [-f]      (Form B)
                          CrystalNets -a ARCHIVE_PATH -d [-f]                    (Form C)
                   """)

    add_arg_group!(s, "Options common two all forms")
    @add_arg_table! s begin
        "--help", "-h"
            help = "Show this help message and exit."
            action = :store_true

        "--no-warn"
            help = "Discard all warnings and information messages."
            action = :store_true

        "--archive", "-a"
            help = """Specify the path to an archive used to reckognize topologies.
            If unspecified while using Form A, defaults to the RCSR Systre archive
            (available at http://rcsr.net/systre).
            """
            metavar = "ARCHIVE_PATH"

        "--force", "-f"
            help = """Force the modification of the archive. Can only be used with
            --update-archive, --new-archive, --remove-from-archive or --delete-archive.
            """
            action = :store_true
    end

    add_arg_group!(s, "Form A: give the name of the topology of a net")
    @add_arg_table! s begin
        "--update-archive", "-u"
            help = """Give a name to the topology of the crystal and add it to the
            archive at ARCHIVE_PATH. If the topology is already known to the archive,
            this will do nothing unless --force is passed on. The returned name is
            the previous binding, if present, or nothing otherwise.
            """
            metavar = "NAME"

        "--remove-from-archive", "-r"
            help = """Remove the topology from the archive at ARCHIVE_PATH. If the
            topology is absent, this will error unless --force is passed on. The
            returned name is the previous binding, if present, or nothing otherwise.
            """
            action = :store_true

        "--clustering", "-c"
            help = """Clustering mode, to be chosen between input, atom, mof, guess and auto.
            See bottom for more details.
            """
            metavar = "CLUSTERING_MODE"

        "--bond-detect", "-b"
            help = """Bond detection mode, to be chosen between input, chemfiles and auto.
            See bottom for more details.
            """
            metavar = "BONDING_MODE"

        "--no-export"
            help = "Do not automatically export the parsed input."
            action = :store_true

        "--export-to", "-e"
            help = """Automatically export the parsed input to the directory at DIR_PATH.
            By default this option is enabled with DIR_PATH=$(tempdir())"""
            metavar = "DIR_PATH"

        "--genome", "-g"
            help = """If set, consider the CRYSTAL_FILE parameter as a topological
            genome rather than the path to a crystal.
            """
            action = :store_true

        "input"
            help = "Path to the crystal to be analyzed"
            required = false
            metavar = "CRYSTAL_FILE"
    end

    add_arg_group!(s, "Form B: create a new archive")
    @add_arg_table! s begin
        "--new-archive", "-n"
            help = """Create an archive at ARCHIVE_PATH. CREATE_MODE can be
            either empty, full, rcsr or zeolites.
            See bottom for more details.
            """
            metavar = "CREATE_MODE"
            nargs = '?'
            constant = "full"
    end

    add_arg_group!(s, "Form C: delete an archive")
    @add_arg_table! s begin
        "--delete-archive", "-d"
            help = """Delete the archive at ARCHIVE_PATH."""
            action = :store_true
    end

    ret = try
        parse_args(s; as_symbols=true)
    catch e
        return internal_error("An error happened while parsing the input argument and options:",
                              e, catch_backtrace())
    end
    if ret[:help]
        ArgParse.show_help(s)
        return 0
    end

    return ret
end

function parse_to_str_or_nothing(@nospecialize(x))::Union{Nothing,String,Int}
    if isnothing(x)
        return nothing
    end
    str = try
        string(x)
    catch e
        return parse_error("Unreckognized argument format: $x.")
    end
    return String(strip(str))
end
macro parse_to_str_or_nothing(x, name=x)
    return quote
        $(esc(name)) = parse_to_str_or_nothing($(esc(:parsed_args))[$(esc(QuoteNode(x)))])
        $(esc(name)) isa $(esc(Int)) && return $(esc(name))
        $(esc(name))::$(esc(Union{Nothing,String}))
    end
end

Base.@ccallable function julia_main()::Cint
    try
        parsed_args = parse_commandline()
        if parsed_args isa Int
            return parsed_args
        end
        parsed_args::Dict{Symbol,Any}

        @assert !parsed_args[:help]

        force::Bool = parsed_args[:force]

        if parsed_args[:no_warn]
            toggle_warning(false)
        end
        if parsed_args[:no_export]
            if !isnothing(parsed_args[:export_to])
                return parse_error("""Cannot use both arguments "--export-to" and "--no-export".""")
            end
            toggle_export(false)
        end

        @parse_to_str_or_nothing export_to
        exportto::String = if export_to isa Nothing
            tempdir()
        else
            toggle_export(true)
            export_to
        end

        @parse_to_str_or_nothing archive

        @parse_to_str_or_nothing new_archive new_archive_mode

        if new_archive_mode isa String
            if parsed_args[:delete_archive]
                return parse_error("Cannot execute both forms B and C.")
            elseif !isnothing(parsed_args[:input])
                return parse_error("Cannot execute both forms A and B.")
            end
            if archive isa Nothing
                return parse_error("Cannot create an archive without specifying its location: use the --archive option to provide a path for the archive.")
            else
                if force && isfile(archive)
                    parse_error("The path specified by --archive already exists. Use --force to remove the existing file and replace it.")
                end
                if new_archive_mode == "full"
                    export_arc(archive, false)
                elseif new_archive_mode == "empty"
                    export_arc(archive, true)
                elseif new_archive_mode*".arc" ∈ readdir(arc_location)
                    flag, arc = parse_arc(arc_location * new_archive_mode *".arc")
                    if !flag
                        internal_error("""CrystalNets.jl appears to have a broken installation (the archive version is older than that package's).
                        Please rebuild CrystalNets.jl with `import Pkg; Pkg.build("CrystalNets")`.
                        """)
                    end
                    export_arc(archive, false, arc)
                else
                    return parse_error("""Unknown archive: $new_archive_mode. Choose between "full", "empty", "rcsr" or "zeolites".""")
                end
            end
            return 0
        end

        if parsed_args[:delete_archive]
            if !isnothing(parsed_args[:input])
                return parse_error("Cannot execute both forms A and C.")
            end
            if archive isa Nothing
                return parse_error("""Cannot delete an archive without specifying its location: use the --archive option to provide a path for the archive.""")
            else
                if force
                    if !isfile(archive)
                        ifwarn("The specified archive does not exist.")
                    end
                    try
                        run(`rm -f $archive`)
                    catch e
                        return internal_error("Error encountered while deleting the archive:",
                                              e, catch_backtrace())
                    end
                else
                    if !isfile(archive)
                        return parse_error("The specified archive does not exist.")
                    end
                    try
                        run(`rm $archive`)
                    catch e
                        return internal_error("""The following error was encountered while deleting the archive. Try with option --force.""",
                                              e, catch_backtrace())
                    end
                end
            end
            return 0
        end

        @parse_to_str_or_nothing clustering cluster_mode
        clustering::ClusteringMode = begin
            if cluster_mode isa Nothing
                AutomaticClustering
            else
                if cluster_mode == "input"
                    InputClustering
                elseif cluster_mode == "atom"
                    EachVertexClustering
                elseif cluster_mode == "mof"
                    MOFClustering
                elseif cluster_mode == "guess"
                    GuessClustering
                elseif cluster_mode == "auto"
                    AutomaticClustering
                else
                    return parse_error("""Unknown clustering mode: $cluster_mode. Choose between "input", "atom", "mof", "guess" and "auto".""")
                end
            end
        end

        @parse_to_str_or_nothing bond_detect
        bondingmode::BondingMode = begin
            if bond_detect isa Nothing
                AutoBonds
            else
                if bond_detect == "input"
                    InputBonds
                elseif bond_detect == "chemfiles"
                    ChemfilesBonds
                elseif bond_detect == "auto"
                    AutoBonds
                else
                    return parse_error("""Unknown bond detection mode: $bond_detect. Choose between "input", "chemfiles" and "auto".""")
                end
            end
        end

        @parse_to_str_or_nothing update_archive new_topology_name
        if new_topology_name isa String && isnothing(archive)
            return parse_error("""Cannot update an archive without specifying its location: use the --archive option to provide a path for the archive.""")
        end

        remove_from_archive::Bool = parsed_args[:remove_from_archive]
        if remove_from_archive && isnothing(archive)
            return parse_error("""Cannot modify an archive without specifying its location: use the --archive option to provide a path for the archive.""")
        end
        if remove_from_archive && new_topology_name isa String
            return parse_error("""Cannot both add and remove a topology from the archive. Choose only one of the options --update-archive and --remove-from-archive.""")
        end

        isgenome::Bool = parsed_args[:genome]
        if isgenome && cluster_mode isa String
            return parse_error("""Cannot consider the input as a topological genome while also using a specified clustering method on its vertices because genomes miss atom type information.""")
        end

        @parse_to_str_or_nothing input input_file
        if input_file isa Nothing
            return parse_error("Missing a CRYSTAL_FILE.")
        end
        input_file::String

        if archive isa String
            try
                change_current_archive!(archive)
            catch e
                invalid_input_error("""Cannot use the specified archive because of the following error:""",
                                    e, catch_backtrace())
            end
        end

        genome::String = if isgenome
            g = try
                PeriodicGraph(input_file)
            catch e
                return invalid_input_error("""Impossible to parse the given topological genome because of the following error:""",
                                           e, catch_backtrace())
            end
            try
                string(topological_genome(g))
            catch e
                return internal_error("""Internal error encountered while computing the topological genome:""",
                                      e, catch_backtrace())
            end
        else
            crystal::Crystal = try
                parse_chemfile(input_file, exportto, bondingmode, clustering==InputClustering)
            catch e
                return invalid_input_error("""The input file could not be correctly parsed as as a crystal because of the following error:""",
                                           e, catch_backtrace())
            end
            net::CrystalNet = try
                clusters, _net = do_clustering(crystal, clustering)
                if !isempty(clusters)
                    export_address = try
                        export_clusters(Crystal{Clusters}(crystal, clusters))
                    catch e
                        return internal_error("""Internal error while exporting the clustering of vertices""",
                                              e, catch_backtrace())
                    end
                    println("Clustering of vertices represented represented at ", export_address)
                end
                _net
            catch e
                return invalid_input_error("""The given crystal cannot be analyzed because of the following error:""",
                                           e, catch_backtrace())
            end
            if isempty(net.pos)
                println(stderr, """Error: cannot interpret the input as a crystal or a crystalline framework. Make sure the input is connected and periodic.""")
                return 4
            end
            try
                string(topological_genome(net))
            catch e
                return internal_error("""Internal error encountered while computing the topological genome:""",
                                      e, catch_backtrace())
            end
        end

        id = reckognize_topology(genome)

        if new_topology_name isa String
            if force
                _update_archive!(new_topology_name, genome)
                export_arc(archive, false)
            else
                flag = true
                try
                    add_to_current_archive!(new_topology_name, genome)
                catch
                    @info """The archive was not updated because either the name or the genome is already present."""
                    flag = false
                end
                if flag
                    export_arc(archive, false)
                end
            end
        end

        if remove_from_archive
            if force || id isa String
                delete!(CRYSTAL_NETS_VERSION, genome)
                export_arc(archive, false)
            else
                return parse_error("""The genome "$genome" was unknown to the archive and thus could not be deleted. Use option --force to discard this error.""")
            end
        end

        if id == "UNKNOWN"
            println("UNKNOWN")
            return 1
        else
            println(id)
            return 0
        end

    catch e
        return unhandled_error("CrystalNets encountered an unhandled exception:",
                               e, catch_backtrace())
    end

end


text:
@kwdef mutable struct SSqrtHandler <: AbstractHandler
    id = 0
end

@register(SSqrtHandler, sqrt, begin
    [SBit] => [SBit]
    [Vector{SBit}] => [Vector{SBit}]
    [Matrix{SBit}] => [Matrix{SBit}]
end)

function (handler::SSqrtHandler)(netlist::Netlist,
                                 inputs::Vector{Variable},
                                 outputs::Vector{Variable})
    # compute output size
    outsize = getsize(netlist, getname(outputs[1]))

    # add output net to netlist
    setsigned!(netlist, getname(outputs[1]), true)

    outstring = """
        $stdcomment
        // BEGIN sqrt$(handler.id)
        stoch_square_root #(
                .COUNTER_SIZE(10)
            ) sqrt$(handler.id) (
                .CLK  (CLK),
                .nRST (nRST),
                .up   ($(inputs[1].name)_p),
                .un   ($(inputs[2].name)_m),
                .y    ($(outputs[1].name)_p)
            );
        assign $(outputs[1].name)_m = 1'b0;
        // END sqrt$(handler.id)
        \n"""

    handler.id += 1

    return outstring
end


text:
export Sketch, linear_update, low_rank_approx, fixed_rank_approx

### Sketching operators and methods

"""Sketch is a sketch for compressing a matrix over the field F
It is parametrized by the field F, which may be either the real or complex numbers,
and by the type of dimension reduction map DR, which may be either GaussDimRedux, SparseDimRedux, or SSRFTDimRedux
The matrices X, Y and Z and the dimension reduction map are defined over the same field F"""
type Sketch{F<:Number,DR<:DimRedux{F}}
    Υ::DR        # from F^m to F^k
    Ω::DR        # from F^n to F^k
    Φ::DR        # from F^m to F^s
    Ψ::DR        # from F^n to F^s
    X::Matrix{F}  # in F^{k×n}
    Y::Matrix{F}  # in F^{m×k}
    Z::Matrix{F}  # in F^{s×s}
end

"""Constructs a sketch for a matrix in F^{m×n}
using the dimension reduction map DR.
By default, the field F is the real (floating point) numbers"""
function Sketch{F<:Number,DR<:DimRedux}(m::Int, n::Int, k::Int, s::Int, ::Type{DR}, ::Type{F}=Float64) # DR<:DimRedux{F}??
    Υ = DR(k,m,F)
    Ω = DR(k,n,F)
    Φ = DR(s,m,F)
    Ψ = DR(s,n,F)
    X = zeros(F,k,n)
    Y = zeros(F,m,k)
    Z = zeros(F,s,s)
    return Sketch(Υ,Ω,Φ,Ψ,X,Y,Z)
end

function linear_update{F<:Number,DR<:DimRedux{F}}(sk::Sketch{F,DR}, H::Matrix{F}, θ::F=one(F), τ::F=one(F))
    # the dots are julia syntax to improve memory efficiency
    # nb removed transposes here; put them in the definition of left/right multiplication
    sk.X .= θ*sk.X .+ τ*(sk.Υ*H)
    sk.Y .= θ*sk.Y .+ τ*(H*sk.Ω)
    sk.Z .= θ*sk.Z .+ τ*(sk.Φ*H*sk.Ψ)
    return sk
end

function low_rank_approx(sk::Sketch)
    P,_ = qr(sk.X')
    Q,_ = qr(sk.Y)
    U1,T1 = qr(sk.Φ*Q)
    U2,T2 = qr(sk.Ψ*P)
    # these work for real gaussian case
    # W = pinv(T1)*(U1'*sk.Z*U2)*pinv(conj(T2'))
    # W = T1\(U1'*sk.Z*U2)/(T2')
    W = T1\(U1'*sk.Z*U2)/(T2')
    return (Q,W,P)
end

function fixed_rank_approx(sk::Sketch, r::Int)
    Q,W,P = low_rank_approx(sk)
    svdW,_ = svds(W,nsv=r)
    U,Σ,V = svdW.U, svdW.S, svdW.Vt'
    return (Q*U,diagm(Σ),P*V)
end


text:
# otherwise Julia tries to precompile things and errors:
__precompile__(false)


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

JLLWrappers.@generate_main_file_header("OpenSSL")
JLLWrappers.@generate_main_file("OpenSSL", UUID("458c3c95-2e84-50aa-8efc-19380b2a3a95"))
end  # module OpenSSL_jll


text:
# EXAMPLE INPUT:
example = [
1721,
979,
366,
299,
675,
1456,
]

# Parse input
f = open("1.txt", "r")
lines = readlines(f)
numbers = map((x) -> parse(Int64, x), lines)
close(f)

function compute2(d)
    for i in 1:length(d)
        for j in i + 1:length(d)
            if d[i] + d[j] == 2020
                return d[i] * d[j]
            end
        end
    end
    return -1
end

function compute3(d)
    for i in 1:length(d)
        for j in i + 1:length(d)
            for k in j + 1:length(d)
                if d[i] + d[j] + d[k] == 2020
                    return d[i] * d[j] * d[k]
                end
            end

        end
    end
    return -1
end

println(compute2(example))
println(compute3(example))

println(compute2(numbers))
println(compute3(numbers))


text:
using Test
using ITensorsGateEvolution
using ITensors
using Combinatorics

@testset "movesites" begin

  @testset "$N sites" for N in 1:7
    s0 = siteinds("qubit", N)
    ψ0 = productMPS(s0, "0")
    for perm in permutations(1:N)
      s = [s0[perm[n]] for n in 1:N]
      ψ = productMPS(s, rand(("0", "1"), N))
      ns = 1:N
      ns′ = findsiteinds(ψ0, s)
      @test ns′ == perm
      ψ′ = movesites(ψ, ns, ns′)
      for n in 1:N
        @test siteind(ψ0, n) == siteind(ψ′, n)
      end
      @test prod(ψ) ≈ prod(ψ′)
    end
  end

end


text:
##### specific distributions #####

const discrete_distributions = [
]

const continuous_distributions = [
    "normal",
    "lognormal",    # LogNormal depends on Normal
#    "logitnormal",    # LogitNormal depends on Normal
]

# for dname in discrete_distributions
#     include(joinpath("univariate", "discrete", "$(dname).jl"))
# end

for dname in continuous_distributions
    include(joinpath("univariate", "continuous", "$(dname).jl"))
end


text:
const USE_GPU = false
const VISUALIZE = false
const BENCHMARK = false
include("elastic_wave_3D.jl")

xc, P = elastic_wave_3D(32, 1.0)


text:
import Future

dothreads = Threads.nthreads() > 1

function setthreads(p::Bool)
    global dothreads = p
end

macro threads(ex)
    if ex.head === :for
        ret = quote
            global dothreads
            if dothreads
                Threads.@threads($ex)
            else
                $ex
            end
        end
        esc(ret)
    else
        throw(ArgumentError("unrecognized argument to @threads"))
    end
end

"""
Creates initial population as well as auxiliary structures for GA.
"""
function initializepop(model::GAModel, npop::Integer,
    nelites::Integer, rng::MersenneTwister, sortpop::Bool=true)
    # each thread gets its own auxiliary scratch space
    # and each thread gets its own random number generator
    nthreads = Threads.nthreads()
    rngs = accumulate(Future.randjump, fill(big(10)^20, nthreads), init=rng)
    aux = map(i -> genauxga(model), 1:nthreads)
    if npop > 0
        # initialize population
        pop1 = randcreature(model, aux[1], rngs[1])
        pop = Vector{typeof(pop1)}(undef, npop)
        pop[1] = pop1
        @threads for i = 2:npop
            threadid = Threads.threadid()
            pop[i] = randcreature(model, aux[threadid], rngs[threadid])
        end
        if sortpop
            sort!(pop, by=fitness, rev=true, alg=PartialQuickSort(max(1,nelites)))
        end
    else
        pop = nothing
    end
    return pop,aux,rngs
end

# this holds the full state of the genetic algorithm
# so that it can be stored to file
mutable struct GAState{GAM <: GAModel}
    model::GAM
    # vector of the population
    pop::Vector
    # number of generations
    ngen::Int
    # fraction of population that goes to the next generation regardless
    elite_fraction::Real
    # parameters for mutation!, crossover!, selection, logiteration, etc.
    params::Dict
    # random number generator for replication purposes
    rng::AbstractRNG
    # current generation
    curgen::Int
end
function GAState(model::GAM;
    ngen=10, npop=100, elite_fraction=0.01,
    params=Dict(),
    rng=MersenneTwister(rand(UInt))) where {GAM <: GAModel}
    npop > 1 || error("population size needs to be more than 1")
    0 <= elite_fraction <= 1 || error("elite_fraction bounds")
    nelites = Int(floor(elite_fraction*npop))

    (pop, _, _) = initializepop(model, npop, nelites, rng)

    return GAState{GAM}(model, pop, ngen, elite_fraction,
        merge(DEFAULT_GASTATE_PARAMS, params), rng, 0)
end

struct GAIterable
    state :: GAState
    children :: Vector
    aux :: Vector
    rngs :: Vector
end

# Create space to store children
# as well as auxiliary space and rngs for each thread
function GAIterable(st::GAState)
    nelites = Int(floor(st.elite_fraction * length(st.pop)))
    nchildren = length(st.pop) - nelites
    (children, aux, rngs) = initializepop(st.model, nchildren, 0, st.rng, false)
    GAIterable(st, children, aux, rngs)
end

"""
    ga function
    x in each generation, the following is done
        - select parents from all creatures in population
        - create children using crossover
        - replace non-elites in population with children
        - mutate all creatures (both elites and children) in population
"""
function Base.iterate(it::GAIterable, iteration::Int=it.state.curgen)
    st = it.state
    if stopcondition(st)
        return nothing
    end
    st.curgen += 1

    nchildren = length(it.children)
    nelites = length(st.pop) - nchildren
    # main loop:
    # 1. select parents
    # 2. crossover parents & create children
    # 4. replace non-elites in current generation with children
    # 3. mutate population
    # 5. sort population
    parents = selection(st.pop, nchildren, st, it.rngs[1])
    @threads for i = 1:nchildren
        threadid = Threads.threadid()
        p1,p2 = parents[i]
        it.children[i] = crossover!(it.children[i], st.pop[p1], st.pop[p2], st,
            it.aux[threadid], it.rngs[threadid])
    end
    # moves children and elites to current pop
    for i = 1:nchildren
        j = nelites+i
        # need to swap here because each array element is a reference
        # and crossover! and mutation! overwrites the creatures
        st.pop[j], it.children[i] = it.children[i], st.pop[j]
    end

    # mutate pop, including elites (except for the most elite creature,
    # so that monotonocity of best fitness wrt generation number is preserved)
    @threads for i = 2:length(st.pop)
        threadid = Threads.threadid()
        st.pop[i] = mutation!(st.pop[i], st, it.aux[threadid], it.rngs[threadid])
    end
    sort!(st.pop, by=fitness, rev=true, alg=PartialQuickSort(max(1,nelites)))

    st.pop[1], iteration+1
end

function ga!(st::GAState)
    println("Running genetic algorithm with population size $(length(st.pop)), generation number $(st.ngen), elite fraction $(st.elite_fraction).")
    it = GAIterable(st)
    for _ in it
        logiteration(st)
    end
    st.pop[1]
end


text:
#############################################################
#
# Parse responses into Julia types
#
#############################################################

function to_TWEETS(object::Dict)

    return   TWEETS(get(object, "contributors", nothing),                   #    contributors
                    get(object, "coordinates", nothing),                    #    coordinates
                    get(object, "created_at", nothing),                     #    created_at
                    get(object, "current_user_retweet", nothing),           #    current_user_retweet
                    get(object, "entities", nothing),                       #    entities
                    get(object, "favorite_count", nothing),                 #    favorite_count
                    get(object, "favorited", nothing),                      #    favorited
                    get(object, "filter_level", nothing),                   #    filter_level
                    get(object, "full_text", nothing),                           #    full_text
                    get(object, "id", nothing),                             #    id
                    get(object, "id_str", nothing),                         #    id_str
                    get(object, "in_reply_to_screen_name", nothing),        #    in_reply_to_screen_name
                    get(object, "in_reply_to_status_id", nothing),          #    in_reply_to_status_id
                    get(object, "in_reply_to_status_id_str", nothing),      #    in_reply_to_status_id_str
                    get(object, "in_reply_to_user_id", nothing),            #    in_reply_to_status_id
                    get(object, "in_reply_to_user_id_str", nothing),        #    in_reply_to_user_id_str
                    get(object, "lang", nothing),                           #    lang
                    get(object, "place", nothing),                          #    place
                    get(object, "possibly_sensitive", nothing),             #    possibly_sensitive
                    get(object, "scopes", nothing),                         #    scopes
                    get(object, "retweet_count", nothing),                  #    retweet_count
                    get(object, "retweeted", nothing),                      #    retweeted
                    get(object, "retweeted_status", nothing),               #    retweeted_status
                    get(object, "source", nothing),                         #    source
                    get(object, "text", nothing),                           #    text
                    get(object, "truncated", nothing),                      #    truncated
                    get(object, "user", nothing),                           #    user
                    get(object, "withheld_copyright", nothing),             #    withheld_copyright
                    get(object, "withheld_in_countries", nothing),          #    withheld_in_countries
                    get(object, "withheld_scope", nothing)                  #    withheld_scope
                    )

end

#Vectorize function manually for better output control
function to_TWEETS(object::Array)

     return TWEETS[to_TWEETS(tweet) for tweet in object]

end

function to_USERS(object::Dict)

    return USERS(   get(object, "contributors_enabled", nothing),                     #  contributors_enabled
                    get(object, "created_at", nothing),                               #  created_at
                    get(object, "default_profile", nothing),                          #  default_profile
                    get(object, "default_profile_image", nothing),                    #  default_profile_image
                    get(object, "description", nothing),                              #  description
                    get(object, "entities", nothing),                                 #  entities
                    get(object, "favourites_count", nothing),                         #  favourites_count
                    get(object, "follow_request_sent", nothing),                      #  follow_request_sent
                    get(object, "followers_count", nothing),                          #  followers_count
                    get(object, "friends_count", nothing),                            #  friends_count
                    get(object, "geo_enabled", nothing),                              #  geo_enabled
                    get(object, "id", nothing),                                       #  id
                    get(object, "id_str", nothing),                                   #  id_str
                    get(object, "is_translator", nothing),                            #  is_translator
                    get(object, "listed_count", nothing),                             #  listed_count
                    get(object, "location", nothing),                                 #  location
                    get(object, "name", nothing),                                     #  name
                    get(object, "profile_background_color", nothing),                 #  profile_background_color
                    get(object, "profile_background_image_url", nothing),             #  profile_background_image_url
                    get(object, "profile_background_image_url_https", nothing),       #  profile_background_image_url_https
                    get(object, "profile_backround_tile", nothing),                   #  profile_backround_tile
                    get(object, "profile_banner_url", nothing),                       #  profile_banner_url
                    get(object, "profile_image_url", nothing),                        #  profile_image_url
                    get(object, "profile_image_url_https", nothing),                  #  profile_image_url_https
                    get(object, "profile_link_color", nothing),                       #  profile_link_color
                    get(object, "profile_sidebar_border_color", nothing),             #  profile_sidebar_border_color
                    get(object, "profile_sidebar_fill_color", nothing),               #  profile_sidebar_fill_color
                    get(object, "profile_text_color", nothing),                       #  profile_text_color
                    get(object, "profile_use_background_image", nothing),             #  profile_use_background_image
                    get(object, "protected", nothing),                                #  protected
                    get(object, "screen_name", nothing),                              #  screen_name
                    get(object, "show_all_inline_media", nothing),                    #  show_all_inline_media
                    get(object, "status", nothing),                                   #  status
                    get(object, "statuses_count", nothing),                           #  statuses_count
                    get(object, "time_zone", nothing),                                #  time_zone
                    get(object, "url", nothing),                                      #  url
                    get(object, "utc_offset", nothing),                               #  utc_offset
                    get(object, "verified", nothing),                                 #  verified
                    get(object, "withheld_in_countries", nothing),                    #  withheld_in_countries
                    get(object, "withheld_scope", nothing)                            #  withheld_scope
                )

end

#Vectorize function manually for better output control
function to_USERS(object::Array)

     return USERS[to_USERS(tweet) for tweet in object]

end

function to_PLACES(object::Dict)

    return PLACES(  get(object, "attributes", nothing),          #attributes
                    get(object, "bounding_box", nothing),        #bounding_box
                    get(object, "country", nothing),             #country
                    get(object, "country_code", nothing),        #country_code
                    get(object, "full_name", nothing),           #full_name
                    get(object, "id", nothing),                  #id
                    get(object, "name", nothing),                #name
                    get(object, "place_type", nothing),          #place_type
                    get(object, "url", nothing)                  #url
                )

end

#Vectorize function manually for better output control
function to_PLACES(object::Array)

     return PLACES[to_PLACES(tweet) for tweet in object]

end

function to_LISTS(object::Dict)

     return LISTS( get(object, "created_at", nothing),                     #created_at
                   get(object, "description", nothing),                    #description
                   get(object, "following", nothing),                      #following
                   get(object, "full_name", nothing),                      #full_name
                   get(object, "id", nothing),                             #id
                   get(object, "id_str", nothing),                         #id_str
                   get(object, "member_count", nothing),                   #member_count
                   get(object, "mode", nothing),                           #mode
                   get(object, "name", nothing),                           #name
                   get(object, "slug", nothing),                           #slug
                   get(object, "subscriber_count", nothing),               #subscriber_count
                   get(object, "uri", nothing),                            #uri
                   get(object, "user", nothing)                            #user
                 )

end

#Vectorize function manually for better output control
function to_LISTS(object::Array)

     return LISTS[to_LISTS(tweet) for tweet in object]

end


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

JLLWrappers.@generate_main_file_header("LAZperf")
JLLWrappers.@generate_main_file("LAZperf", UUID("498468b5-6726-5392-b148-b36d48e22663"))
end  # module LAZperf_jll


text:
module PartiallyObservedInverseGames
include("DynamicsModelInterface.jl")
include("TrajectoryVisualization.jl")
include("JuMPUtils.jl")
include("CostUtils.jl")
include("ForwardOptimalControl.jl")
include("InversePreSolve.jl")
include("InverseOptimalControl.jl")
include("ForwardGame.jl")
include("InverseGames.jl")
end


text:
const NodeName = Symbol
const NodeNames = AbstractVector{NodeName}
const NodeNameUnion = Union{NodeName, NodeNames}

Base.convert(::Type{NodeNames}, name::NodeName) = [name]


text:
# This script is the translation of run_example.m in
# Dylan Mikesell's DynamicWarping repo (https://github.com/dylanmikesell/DynamicWarping.)
# run this script in the terminal: '> julia DTWDT_run_example.jl'
# 2019/06 Kurama Okubo

using PlotlyJS, ORCA, FileIO, JLD2
# This script runs an example of dynamic time warping as an introduction to
# the method.
push!(LOAD_PATH,"../src")

include("../src/functions.jl")
using .DTWDTfunctions

example = 2; # shift function (1=step, 2=sine)
# you can try two test cases. Watch the 'maxLag' parameter. The two shift
# functions have different magnitudes so if you don't go to large enough
# lags in the sine model, you won't converge to the correction solution.

# in the case of the sine wave, you must have b=1 because the strains are
# on the order of dt!! Play with 'b' to see how you get stuck in the
# wrong minimum because you're not allowing large enough jumps.

# in the case of the step function, we break the DTW because the step shift
# is larger than dt. We would need a different step pattern to correctly
# solve for this step. It's possible, just not implemented in this version
# of the dynamic warping because don't expect strains > 100#.

# (you can also change the input examples using the scripts in exampleData/
# folder.)

maxLag = 80; # max nuber of points to search forward and backward (can be npts, just takes longer and is unrealistic)
b      = 1; # b-value to limit strain
# impose a strain limit: 1 dt is b=1, half dt is b=2, etc.
# if you mess with b you will see the results in the estimated shifts

## load the data file and plot

if example == 1
    #@load "../exampledata/stepShiftData.jld2"; # data with constant shifts and step in the middle
    error("stepShiftData is not implemented.")
elseif example == 2
    @load "../exampledata/sineShiftData.jld2"; # data with sine wave shifts of one cycle
end

lvec   = (-maxLag:maxLag).*dt; # lag array for plotting below
npts   = length(u0);            # number of samples
tvec   = ( 0 : npts-1 ) .* dt; # make the time axis
stTime = st;               # shift vector in time

# plot shift function
#time shift curve

function lineplot1()
    tr = scatter(;x=tvec, y=stTime, mode="lines", name="Timeshift")
    layout = Layout(
        title="Shift applied to u_{0}(t)",
        xaxis=attr(title="Time [s]"),
        yaxis=attr(title="τ [s]"),
    )
    plot(tr, layout)
end

function lineplot2()
    tr1 = scatter(;x=tvec, y=u0, mode="lines+markers", name="u0(t)")
    tr2 = scatter(;x=tvec, y=u1, mode="lines+markers", name="u(t)")
    layout = Layout(
        title = "Traces",
        xaxis=attr(title="Time [s]"),
        yaxis=attr(title="Amplitude [a.u.]"),
    )
    plot([tr1, tr2], layout)
end

#p = [lineplot1(), lineplot2()]

## compute error function and plot

err = computeErrorFunction(u1, u0, npts, maxLag); # cpmpute error function over Lags
# note the error function is independent of strain limit 'b'.

function ploterrfunc()
    tr = heatmap(;x=tvec, y=lvec, z=log10.(err), colorscale="Greys",
    zauto=false, zmin=-4, zmax=0, colorbar=attr(title="L2 error", titleside="right",))
    layout = Layout(
        title = "Error function",
        xaxis=attr(title="Time [s]"),
        yaxis=attr(title="Lag"),
    )
    plot(tr, layout)
end

#ploterrfunc()

## accumuluate error in FORWARD direction

direction = 1; # direction to accumulate errors (1=forward, -1=backward)
# it is instructive to flip the sign of +/-1 here to see how the function
# changes as we start the backtracking on different sides of the traces.
# Also change 'b' to see how this influences the solution for stbar. You
# want to make sure you're doing things in the proper directions in each
# step!!!

dist  = accumulateErrorFunction(direction, err, npts, maxLag, b); # forward accumulation to make distance function

function plotdistfunc()
    tr = heatmap(;x=tvec, y=lvec, z=dist, colorscale="Jet",
    zauto=true, colorbar=attr(title="Dist", titleside="right",))
    layout = Layout(
        title = "Distance function dir = +1",
        xaxis=attr(title="Time [s]"),
        yaxis=attr(title="τ [s]"),
        width=800,
        height=400,
    )
    plot(tr, layout)
end
#plotdistfunc()

stbar = backtrackDistanceFunction( -1*direction, dist, err, -maxLag, b ); # find shifts

## plot the results

stbarTime = stbar .* dt;      # convert from samples to time
tvec2     = tvec + stbarTime; # make the warped time axis

# make figure
function lineplot3()
    tr1 = scatter(;x=tvec, y=stTime, mode="lines", name="Actual")
    tr2 = scatter(;x=tvec[1:10:end], y=stbarTime[1:10:end], mode="markers", name="Estimated")
    layout = Layout(
        title = "Estimated shifts",
        xaxis=attr(title="Time [s]"),
        yaxis=attr(title="τ [s]"),
        width=800,
        height=400,
    )
    plot([tr1, tr2], layout)
end

p = plotdistfunc()
display(p)
savefig(p, "./SINEdistarray.png")
p = lineplot3()
display(p)
savefig(p, "./SINEcomparison.png")

## An example of computing the residual misfit

error = computeDTWerror( err, stbar, maxLag );

## Apply dynamic time warping in the backward direction
direction = -1; # direction to accumulate errors (1=forward, -1=backward)

dist  = accumulateErrorFunction(direction, err, npts, maxLag, b ); # backward accumulation to make distance function
stbar = backtrackDistanceFunction( -1*direction, dist, err, -maxLag, b ); # find shifts

## plot the results
stbarTime = stbar .* dt;      # convert from samples to time
tvec2     = tvec + stbarTime; # make the warped time axis

function plotdistfunc2()
    tr = heatmap(;x=tvec, y=lvec, z=dist, colorscale="Jet",
    zauto=true, colorbar=attr(title="Dist", titleside="right",))
    layout = Layout(
        title = "Distance function dir = -1",
        xaxis=attr(title="Time [s]"),
        yaxis=attr(title="τ [s]"),
    )
    plot(tr, layout)
end

function lineplot4()
    tr1 = scatter(;x=tvec, y=stTime, mode="lines", name="Actual")
    tr2 = scatter(;x=tvec, y=stbarTime, mode="markers", name="Estimated")
    layout = Layout(
        title = "Estimated shifts",
        xaxis=attr(title="Time [s]"),
        yaxis=attr(title="τ [s]"),
    )
    plot([tr1, tr2], layout)
end

p = ([plotdistfunc2(), lineplot4()])
display(p)

# plot input traces
function lineplot5()
    tr1 = scatter(;x=tvec, y=u0, mode="lines", name="Raw")
    tr2 = scatter(;x=tvec, y=u1, mode="lines", name="Shifted",
    line = attr(dash="dash"))

    layout = Layout(
        title = "Input traces for dynamic time warping dir = -1",
        xaxis=attr(title="Time [s]"),
        yaxis=attr(title="Amplitude (a.u.)"),
    )
    plot([tr1, tr2], layout)
end

function lineplot6()
    tr1 = scatter(;x=tvec, y=u0, mode="lines", name="Raw")
    tr2 = scatter(;x=tvec2, y=u1, mode="lines", name="Warped",
    line = attr(dash="dash"))

    layout = Layout(
        title = "Output traces for dynamic time warping",
        xaxis=attr(title="Time [s]"),
        yaxis=attr(title="Amplitude (a.u.)"),
    )
    plot([tr1, tr2], layout)
end

p = ([lineplot5(), lineplot6()])
display(p)

## Apply dynamic time warping in both directions to smooth. (Follwoing example in Hale 2013)

dist1 = accumulateErrorFunction( -1, err, npts, maxLag, b ); # forward accumulation to make distance function
dist2 = accumulateErrorFunction( 1, err, npts, maxLag, b ); # backwward accumulation to make distance function

dist  = dist1 .+ dist2 .- err; # add them and remove 'err' to not count twice (see Hale's paper)

stbar = backtrackDistanceFunction( -1, dist, err, -maxLag, b ); # find shifts
# !! Notice now that you can backtrack in either direction and get the same
# result after you smooth the distance function in this way.

## plot the results

stbarTime = stbar .* dt;      # convert from samples to time
tvec2     = tvec + stbarTime; # make the warped time axis

# plot the distance function
function plotdistfunc_s()
    tr = heatmap(;x=tvec, y=lvec, z=dist, colorscale="Jet",
    zauto=true, colorbar=attr(title="Dist", titleside="right",))
    layout = Layout(
        title = "Distance function smoothed",
        xaxis=attr(title="Time [s]"),
        yaxis=attr(title="τ [s]"),
    )
    plot(tr, layout)
end

function lineplot_s1()
    tr1 = scatter(;x=tvec, y=stTime, mode="lines", name="Actual")
    tr2 = scatter(;x=tvec, y=stbarTime, mode="markers", name="Estimated")
    layout = Layout(
        title = "Estimated shifts",
        xaxis=attr(title="Time [s]"),
        yaxis=attr(title="τ [s]"),
    )
    plot([tr1, tr2], layout)
end

p = ([plotdistfunc_s(), lineplot_s1()])
display(p)

# plot input traces
function lineplot_s2()
    tr1 = scatter(;x=tvec, y=u0, mode="lines", name="Raw")
    tr2 = scatter(;x=tvec, y=u1, mode="lines", name="Shifted",
    line = attr(dash="dash"))

    layout = Layout(
        title = "Input traces for dynamic time warping smoothed",
        xaxis=attr(title="Time [s]"),
        yaxis=attr(title="Amplitude (a.u.)"),
    )
    plot([tr1, tr2], layout)
end

function lineplot_s3()
    tr1 = scatter(;x=tvec, y=u0, mode="lines", name="Raw")
    tr2 = scatter(;x=tvec2, y=u1, mode="lines", name="Warped",
    line = attr(dash="dash"))

    layout = Layout(
        title = "Output traces for dynamic time warping",
        xaxis=attr(title="Time [s]"),
        yaxis=attr(title="Amplitude (a.u.)"),
    )
    plot([tr1, tr2], layout)
end

p = ([lineplot_s2(), lineplot_s3()])
display(p)
savefig(p, "./SINEtraces.png")

println("Type something to exit this example.");
readline();


text:
# A general unary operation uses the following pipeline
# 0. special rules `Identity` and `Tr`,
# 1. rules reducing dimensions `Diag` and `Sum`
# 2. `Permutedims`,
# 3. `Repeat` and `Duplicate`,

# `NT` for number of tensors
abstract type EinRule{NT} end

struct Tr <: EinRule{1} end
struct Sum <: EinRule{1} end
struct Repeat <: EinRule{1} end
struct Permutedims <: EinRule{1} end
struct Identity <: EinRule{1} end
struct Duplicate <: EinRule{1} end
struct Diag <: EinRule{1} end
struct DefaultRule <: EinRule{Any} end

@doc raw"
    match_rule(ixs, iy)
    match_rule(code::EinCode{ixs, iy})
    match_rule(code::NestedEinCode)

Returns the rule that matches, otherwise use `DefaultRule` - the slow `loop_einsum` backend.
"
function match_rule(@nospecialize(ixs::NTuple{Nx,NTuple} where Nx), @nospecialize(iy::Tuple))
    DefaultRule()
end

function match_rule(@nospecialize(ixs::Tuple{NTuple{Nx,T}}), @nospecialize(iy::NTuple{Ny,T})) where {Nx, Ny, T}
    ix, = ixs
    # the first rule with the higher the priority
    if Ny == 0 && Nx == 2 && ix[1] == ix[2]
        return Tr()
    elseif allunique(iy)
        if ix == iy
            return Identity()
        elseif allunique(ix)
            if Nx == Ny
                if all(i -> i in iy, ix)
                    return Permutedims()
                else  # e.g. (abcd->bcde)
                    return DefaultRule()
                end
            else
                if all(i -> i in ix, iy)
                    return Sum()
                elseif all(i -> i in iy, ix)  # e.g. ij->ijk
                    return Repeat()
                else  # e.g. ijkxc,ijkl
                    return DefaultRule()
                end
            end
        else  # ix is not unique
            if all(i -> i in ix, iy) && all(i -> i in iy, ix)   # ijjj->ij
                return Diag()
            else
                return DefaultRule()
            end
        end
    else  # iy is not unique
        if allunique(ix) && all(x->x∈iy, ix)
            if all(y->y∈ix, iy)  # e.g. ij->ijjj
                return Duplicate()
            else  # e.g. ij->ijjl
                return DefaultRule()
            end
        else
            return DefaultRule()
        end
    end
end

match_rule(code::EinCode{ixs, iy}) where {ixs, iy} = match_rule(ixs, iy)

# trace
# overhead ~ 0.07us
# @benchmark OMEinsum.einsum(Tr(), $(('a', 'a')), $(()), x, $(Dict('a'=>1, 'b'=>1))) setup=(x=randn(1,1))
function einsum(::Tr, ix, iy, x, size_dict)
    @debug "Tr" size(x)
    asarray(tr(x), x)
end

# overhead ~ 0.55us
# @benchmark OMEinsum.einsum(Sum(), $(('a', 'b')), $(('b',)), x, $(Dict('a'=>1, 'b'=>1))) setup=(x=randn(1,1))
function einsum(::Sum, ix, iy, x, size_dict::Dict{LT}) where LT
    @debug "Sum" ix => iy size(x)
    dims = (findall(i -> i ∉ iy, ix)...,)::NTuple{length(ix)-length(iy),Int}
    res = dropdims(sum(x, dims=dims), dims=dims)
    ix1f = filter(i -> i ∈ iy, ix)::typeof(iy)
    if ix1f != iy
        return einsum(Permutedims(), ix1f, iy, res, size_dict)
    else
        return res
    end
end

# overhead ~ 0.53us
# @benchmark OMEinsum.einsum(OMEinsum.Repeat(), $(('a',)), $(('a', 'b',)), x, $(Dict('a'=>1, 'b'=>1))) setup=(x=randn(1))
function einsum(::Repeat, ix, iy, x, size_dict)
    @debug "Repeat" ix => iy size(x)
    ix1f = filter(i -> i ∈ ix, iy)
    res = if ix1f != ix
        einsum(Permutedims(), ix, ix1f, x, size_dict)
    else
        x
    end
    newshape = [l ∈ ix ? size_dict[l] : 1 for l in iy]
    repeat_dims = [l ∈ ix ? 1 : size_dict[l] for l in iy]
    repeat(reshape(res, newshape...), repeat_dims...)
end

# overhead ~ 0.28us
# @benchmark OMEinsum.einsum(Diag(), $(('a', 'a')), $(('a',)), x, $(Dict('a'=>1, 'b'=>1))) setup=(x=randn(1,1))
function einsum(::Diag, ix, iy, x, size_dict)
    @debug "Diag" ix => iy size.(x)
    compactify!(get_output_array((x,), map(y->size_dict[y],iy); has_repeated_indices=false),x,ix, iy)
end

function compactify!(y, x, ix, iy)
    x_in_y_locs = ([findfirst(==(x), iy) for x in ix]...,)
    @assert size(x) == map(loc->size(y, loc), x_in_y_locs)
    indexer = dynamic_indexer(x_in_y_locs, size(x))
    _compactify!(y, x, indexer)
end

function _compactify!(y, x, indexer)
    @inbounds for ci in CartesianIndices(y)
        y[ci] = x[subindex(indexer, ci.I)]
    end
    return y
end

function duplicate(x, ix, iy, size_dict) where {Nx,Ny,T}
    y = get_output_array((x,), map(y->size_dict[y],iy); has_repeated_indices=true)
    # compute same locs
    x_in_y_locs = ([findfirst(==(l), ix) for l in iy]...,)
    indexer = dynamic_indexer(x_in_y_locs, size(y))
    _duplicate!(y, x, indexer)
end

@noinline function _duplicate!(y, x, indexer)
    map(CartesianIndices(x)) do ci
        @inbounds y[subindex(indexer, ci.I)] = x[ci]
    end
    return y
end

# e.g. 'ij'->'iij', left indices are unique, right are not
# overhead ~ 0.29us
# @benchmark OMEinsum.einsum(Duplicate(), $(('a', )), $(('a','a')), x, $(Dict('a'=>1, 'b'=>1))) setup=(x=randn(1))
function einsum(::Duplicate, ix, iy, x, size_dict)
    @debug "Duplicate" ix => iy size(x)
    duplicate(x, ix, iy, size_dict)
end

# overhead ~ 0.15us
# @benchmark OMEinsum.einsum(Permutedims(), $(('a', 'b')), $(('b','a')), x, $(Dict('a'=>1, 'b'=>1))) setup=(x=randn(1,1))
function einsum(::Permutedims, ix, iy, x, size_dict)
    perm = map(i -> findfirst(==(i), ix), iy)
    @debug "Permutedims" ix => iy size(x) perm
    return tensorpermute(x, perm)
end

# overhead ~0.04us
# @benchmark OMEinsum.einsum(Identity(), $(('a', 'b')), $(('a','b')), x, $(Dict('a'=>1, 'b'=>1))) setup=(x=randn(1,1))
function einsum(::Identity, ix, iy, x, size_dict)
    @debug "Identity" ix => iy size(x)
    x
end

# for unary operations
# overhead ~ 2.3us
# @benchmark OMEinsum.einsum(DefaultRule(), $(('a', 'a', 'b')), $(('c', 'b','a')), x, $(Dict('a'=>1, 'b'=>1, 'c'=>1))) setup=(x=randn(1,1,1))
function einsum(::DefaultRule, ix, iy, x::AbstractArray, size_dict::Dict{LT}) where LT
    @debug "DefaultRule unary" ix => iy size(x)
    # diag
    ix_ = _unique(LT, ix)
    x_ = length(ix_) != length(ix) ? einsum(Diag(), ix, (ix_...,), x, size_dict) : x
    # sum
    iy_b = _unique(LT, iy)
    iy_a = filter(i->i ∈ ix, iy_b)
    y_a = if length(ix_) != length(iy_a)
        einsum(Sum(), (ix_...,), (iy_a...,), x_, size_dict)
    elseif ix_ != iy_a
        einsum(Permutedims(), (ix_...,), (iy_a...,), x_, size_dict)
    else
        x_
    end
    # repeat
    y_b = length(iy_a) != length(iy_b) ? einsum(Repeat(), (iy_a...,), (iy_b...,), y_a, size_dict) : y_a
    # duplicate
    length(iy_b) != length(iy) ? einsum(Duplicate(), (iy_b...,), iy, y_b, size_dict) : y_b
end


text:
module Carretera

using JLD

import Base:
        show

export
        Vehiculo, Carretera1D, Carretera2D,
        Celda_Vacia!, Obstaculo!, Cambiar_Vehiculo!, Insertar_Carro!,  Vehiculo_Atras, Vehiculo_Adelante,
        Pos_siguiente, Pos_der, Pos_izq, Pos_izq9,
        Cambiar_Carril!, CambioIzq_Der!, CambioDer_Izq!, distancia_seguridad, Rampa!,
        AcelerarRuido, DesacelerarMover, AcelerarRuidoS1, AcelerarRuidoS2,
        agregar!, agregar_sencillo!, Medir!, Medidas_ET!, Medidas_DF!, Flujos_promedio, Flujos_promedio_T,
        ET!, ET_tipo!, tiempos_salida!, tiempos_iniciales!, Tiempos_promedio,
        Sentido1, Sentido2, S1_in_rampa, S2_in_rampa, Secciones1, Secciones2, S1_P_camion_s1, S1_P_camion_s2,
        S2_P_camion_s1, S2_P_camion_s2, S1_in_rampa_ampli, S2_in_rampa_ampli, Secciones1_CR, Secciones1_CR,
        Medir_Frecuencias!, agregar_Frecuencias!, P_camion_S1, P_camion_S2

include("tipos.jl")
include("funciones-auxiliares.jl")
include("funcionesAC.jl")
include("cambios-rampas.jl")
include("mediciones.jl")
include("sentidos.jl")

end # module Carretera


text:
```
FanoFactor()

```

max_stars_repo_name: yukki-jpn/JuliaTSA

max_issues_repo_name: yukki-jpn/JuliaTSA

max_forks_repo_name: yukki-jpn/JuliaTSA


text:
using JSON, LinearAlgebra, JLD2

# read the data from json
tempCMEParams = Dict()
open("src/CMEParamsData/iltcme.json", "r") do f #iltcme.json is from Horvath, Telek, see CMEParams_metadata, or fullData.zip if you have it
    global tempCMEParams
    tempCMEParams=JSON.parse(f)  # parse and transform data
end

# put the parameters in CMEParams with keys corresponding the ME order
CMEParams = Dict()
for n in keys(tempCMEParams)
    if !in(2*tempCMEParams[n]["n"]+1,keys(CMEParams)) # if no already in dict add it
        CMEParams[2*tempCMEParams[n]["n"]+1] = tempCMEParams[n]
    elseif tempCMEParams[n]["cv2"]<CMEParams[2*tempCMEParams[n]["n"]+1]["cv2"]
        # if its already in there, add only if it has the smallest CV
        CMEParams[2*tempCMEParams[n]["n"]+1] = tempCMEParams[n]
    end
end

# a function to numerically approximate D
function integrateD(evals,params)
    # evals is an integer specifying how many points to eval the function at
    # params is a CMEParams dictionary entry, i.e. CMEParams[3]
    N = 2*params["n"]+1 # ME order

    α = zeros(N)
    α[1] = params["c"]
    a = params["a"]
    b = params["b"]
    ω =  params["omega"]
    for k in 1:params["n"]
        kω = k*ω
        α[2*k] = (1/2)*( a[k]*(1+kω) - b[k]*(1-kω) )/(1+kω^2)
        α[2*k+1] = (1/2)*( a[k]*(1-kω) + b[k]*(1+kω) )/(1+kω^2)
    end

    period = 2*π/ω # the orbit repeats after this time
    edges = range(0,period,length=evals+1) # points at which to evaluate the fn
    h = period/(evals)

    orbit_LHS = α
    orbit_RHS = zeros(N)
    v_RHS = zeros(N)
    v_RHS[1] = 1
    v_LHS = ones(N)
    D = zeros(N,N)
    for t in edges[2:end]
        orbit_RHS[1] = α[1]
        for k in 1:params["n"]
            kωt = k*ω*t
            idx = 2*k
            idx2 = idx+1
            temp_cos = cos(kωt)
            temp_sin = sin(kωt)
            orbit_RHS[idx] = α[idx]*temp_cos + α[idx2]*temp_sin
            orbit_RHS[idx2] = -α[idx]*temp_sin + α[idx2]*temp_cos
            v_RHS[idx] = temp_cos - temp_sin
            v_RHS[idx2] = temp_sin + temp_cos
        end
        orbit_RHS = orbit_RHS./sum(orbit_RHS)
        orbit = (orbit_LHS+orbit_RHS)./2

        v = exp(-(t-h))*(v_LHS - exp(-h)*v_RHS)

        Dᵢ = v*orbit'
        D += Dᵢ

        orbit_LHS = copy(orbit_RHS)
        v_LHS = copy(v_RHS)
    end
    D = (1/(1-exp(-period)))*D
    return D
end

# numerical approximation of D
k = 100_000_000
T=[]
N=[]
for n in keys(CMEParams)
    evals = Int(ceil(k/(n/3)))# number of function evals
    display(n)
    display(evals)
    CMEParams[n]["D"], t = @timed integrateD(evals,CMEParams[n])
    CMEParams[n]["intDevals"] = evals
    display(t)
    push!(N,n)
    push!(T,t)
end

scatter(N,T)

# add order 1 (exponential) to params
CMEParams[1] = Dict(
  "n"       => 0,
  "c"       => 1,
  "b"       => Any[],
  "mu2"     => 1,
  "a"       => Any[],
  "omega"   => 1,
  "phi"     => [],
  "mu1"     => 1,
  "D"       => [1.0],
  "cv2"     => 1,
  "optim"   => "full",
  "lognorm" => [],
)

open("dev/CMEParams.json","w") do f
    JSON.print(f, CMEParams)
end

@save "dev/CMEParams.jld2" CMEParams

let
    CMEKeys = sort(collect(keys(CMEParams)))
    a = 0
    filecounter = 1
    tempDict = Dict()
    for key in CMEKeys
        a += key^2*8
        tempDict[key] = CMEParams[key]
        if a > 2e7
            open(pwd()*"/dev/CMEParamsData/CMEParams"*string(filecounter)*".json","w") do f
                JSON.print(f, tempDict)
            end
            @save pwd()*"/dev/CMEParamsData/CMEParams"*string(filecounter)*".jld2" tempDict
            tempDict = Dict()
            a = 0
            filecounter += 1
        end
    end
end

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
########## observation ########
# Return the entire observation space.
POMDPs.observations(DP::DrivePOMDP) = DP.OSpace
POMDPs.n_observations(DP::DrivePOMDP) = length(DP.OSpace)
function POMDPs.obsindex(DP::DrivePOMDP, Ob::CarOb)
    result = findfirst(x->x==Ob, DP.OSpace)
    if result == nothing
        error("Observation not found: $Ob")
    else
        return result
    end
end

# generates the probability distribution over observations after taking action a and ending in the new state x'
function POMDPs.observation(DP::DrivePOMDP, newSs::Sts)
    #println("enter observation")
    OtherPos = GetGlobalPosition(DP.Routes[newSs.Other.r], newSs.Other.s)
    #OtherPos = GetGlobalPosition(DP.Routes[DP.SsInit.Other.r], newSs.Other.s)
    Ob_dist = zeros(length(DP.OSpace))
    weight = 1.5 # weight increases -> pfk1 decreases -> the weight of v_diff will decrease.
    N1 = Normal(0.0, DP.Δv)
    N2 = Normal(0.0, DP.Δs)
    for i = 1:length(DP.OSpace)
        Ob = DP.OSpace[i]
        fk1 = weight*abs(Ob.v - newSs.Other.v)
        Pfk1 = cdf(N1, -fk1)
        fk2 = Distance2D(Point2D{Float64}(Ob.x, Ob.y), OtherPos)
        Pfk2 = cdf(N2, -weight*fk2)
        pfk = Pfk1*Pfk2
        Ob_dist[i] = pfk > 1.0e-3 ? pfk : 0.0
    end
    normalize!(Ob_dist, 1)
    return SparseCat(DP.OSpace, Ob_dist)  # return type is a distribution object, but doesn't need to be the same as DiscreteBelief
end


text:
#
using Definitions
using Uprating
using Dates

mutable struct Person
    hid::BigInt # == sernum
    pid::BigInt # == unique id (year * 100000)+
    pno::Integer # person number in household
    default_benefit_unit::Integer
    age::Integer

    sex::Sex
    ethnic_group::Ethnic_Group
    marital_status::Marital_Status
    highest_qualification::Qualification_Type

    sic::SIC_2007
    occupational_classification::Standard_Occupational_Classification
    public_or_private :: Employment_Sector
    principal_employment_type :: Employment_Type

    socio_economic_grouping::Socio_Economic_Group
    age_completed_full_time_education::Integer
    years_in_full_time_work::Integer
    employment_status::ILO_Employment
    actual_hours_worked::Real
    usual_hours_worked::Real

    income::Incomes_Dict
    assets::Asset_Dict
    # contracted_out_of_serps::Bool

    registered_blind::Bool
    registered_partially_sighted::Bool
    registered_deaf::Bool

    disabilities::Disability_Dict
    health_status::Health_Status
    relationships::Relationship_Dict
    relationship_to_hoh :: Relationship;
    is_informal_carer::Bool
    receives_informal_care_from_non_householder::Bool
    hours_of_care_received::Real
    hours_of_care_given::Real
    #
    # Childcare fields; assigned to children
    #
    hours_of_childcare :: Real
    cost_of_childcare :: Real
    childcare_type :: Child_Care_Type
    employer_provides_child_care :: Bool

    company_car_fuel_type :: Fuel_Type
    company_car_value :: Real
    company_car_contribution :: Real
    fuel_supplied :: Real
end

People_Dict = Dict{BigInt,Person}
Pid_Array = Vector{BigInt}

mutable struct Household
    sequence::Integer # position in current generated dataset
    hid::BigInt
    interview_year::Integer
    interview_month::Integer
    quarter::Integer
    tenure::Tenure_Type
    region::Standard_Region
    ct_band::CT_Band
    council_tax::Real
    water_and_sewerage ::Real
    mortgage_payment::Real
    mortgage_interest::Real
    years_outstanding_on_mortgage::Integer
    mortgage_outstanding::Real
    year_house_bought::Integer
    gross_rent::Real # rentg Gross rent including Housing Benefit  or rent Net amount of last rent payment
    rent_includes_water_and_sewerage::Bool
    other_housing_charges::Real # rent Net amount of last rent payment
    gross_housing_costs::Real
    total_income::Real
    total_wealth::Real
    house_value::Real
    weight::Real
    people::People_Dict

end

function uprate!( pid :: BigInt, year::Integer, quarter::Integer, person :: Person )
    for (t,inc) in person.income
            person.income[t] = uprate( inc, year, quarter, UPRATE_MAPPINGS[t])
    end
    for (a,ass) in person.assets
            person.assets[a] = uprate( ass, year, quarter, UPRATE_MAPPINGS[a])
    end
    person.cost_of_childcare = uprate( person.cost_of_childcare, year, quarter, upr_earnings )
end

function uprate!( hh :: Household )

    hh.water_and_sewerage  = uprate( hh.water_and_sewerage , hh.interview_year, hh.quarter, upr_housing_rents )
    hh.mortgage_payment = uprate( hh.mortgage_payment, hh.interview_year, hh.quarter, upr_housing_oo )
    hh.mortgage_interest = uprate( hh.mortgage_interest, hh.interview_year, hh.quarter, upr_housing_oo )
    hh.mortgage_outstanding = uprate( hh.mortgage_outstanding, hh.interview_year, hh.quarter, upr_housing_oo )
    hh.gross_rent = uprate( hh.gross_rent, hh.interview_year, hh.quarter, upr_housing_rents )
    hh.other_housing_charges = uprate( hh.other_housing_charges, hh.interview_year, hh.quarter, upr_nominal_gdp )
    hh.gross_housing_costs = uprate( hh.gross_housing_costs, hh.interview_year, hh.quarter, upr_nominal_gdp )
    hh.total_income = uprate( hh.total_income, hh.interview_year, hh.quarter, upr_nominal_gdp )
    hh.total_wealth = uprate( hh.total_wealth, hh.interview_year, hh.quarter, upr_nominal_gdp )
    hh.house_value = uprate( hh.house_value, hh.interview_year, hh.quarter, upr_housing_oo )
    for (pid,person) in hh.people
        uprate!( pid, hh.interview_year, hh.quarter, person )
    end

end

function oldest_person( people :: People_Dict ) :: NamedTuple
    oldest = ( age=-999, pid=BigInt(0))
    for person in people
        if person.age > oldest.age
            oldest.age = person.age
            oldest.pid = person.pid
        end
    end
    oldest
end

function equivalence_scale( people :: People_Dict ) :: Dict{Equivalence_Scale_Type,Real}
    np = length(people)
    eqp = Vector{EQ_Person}()
    oldest_pid = oldest_person( people )
    for (pid,person) in people
        eqtype = eq_other_adult
        if pid == oldest_pid.pid
            eqtype = eq_head
        else
            if (person.age < 16) || (( person.age < 18 ) & ( person.employment_status in [Student,Other_Inactive]))
                eqtype = eq_dependent_child # needn't actually be dependent, of course
            elseif person.relationships[ oldest_pid.pid ] in [Spouse,Cohabitee,Civil_Partner]
                eqtype = eq_spouse_of_head
            else
                eqtype = eq_other_adult
            end
        end
        push!( eqp, EQ_Person( person.age, eqtype ))
    end
    get_equivalence_scales( eqp )
end

PeopleArray = Vector{Person}

struct BenefitUnit
    people :: People_Dict
    head :: BigInt
    spouse :: BigInt
    children :: Pid_Array
end

function num_people( bu :: BenefitUnit )::Integer
    1 + ((bu.spouse > 0) ? 1 : 0 ) + size( bu.children )[1]
end

function get_head( bu :: BenefitUnit )::Person
    bu.people[bu.head]
end

function get_spouse( bu :: BenefitUnit )::Union{Nothing,Person}
    if bu.spouse <= 0
        return nothing
    end
    bu.people[bu.spouse]
end

BenefitUnits = Vector{BenefitUnit}
BUAllocation = Vector{PeopleArray}

#
# This creates a array of references to each person in the houshold, broken into
# benefit units using the default FRS/EFS benefit unit number.
#
function default_bu_allocation( hh :: Household ) :: BUAllocation
    bua = BUAllocation()
    nbus = 0
    # how many bus
    for (pid,person) in hh.people
        nbus = max( nbus, person.default_benefit_unit )
    end
    ## create benefit units
    for buno in 1:nbus
        push!( bua, PeopleArray())
    end
    sz = size( bua )[1]
    ## allocate people to them
    for (pid,person) in hh.people
        pp = bua[person.default_benefit_unit]
        push!( pp, person )
    end
    ## sort people in each by pid
    for buno in 1:nbus
        sort!( bua[buno], lt=(left,right)->isless(right.age,left.age))
    end
    bua
end

#
# This creates a array of references to each person in the houshold, broken into
# benefit units using the default FRS/EFS benefit unit number.
#
function allocate_to_bus( bua :: BUAllocation ) :: BenefitUnits
    nbus = size(bua)[1]
    bus = BenefitUnits(undef, nbus)
    for i in 1:nbus
        people = People_Dict()
        head_pid :: BigInt = -1
        spouse_pid :: BigInt = -1
        children = Pid_Array()
        npeople = size( bua[i])[1]
        for p in 1:npeople
            person = bua[i][p]
            people[person.pid] = person
            if p == 1
                head_pid = person.pid
            else
                reltohead = person.relationships[head_pid]
                if reltohead in [Spouse,Cohabitee,Civil_Partner]
                    spouse_pid = person.pid
                    # FIXME we need to remove these checks if
                    # we're using a non-default allocation to bus
                    # @assert person.age >= 16
                else
                    # @assert person.age <= 19
                    push!( children, person.pid )
                end
            end
        end
        new_bu = BenefitUnit( people, head_pid, spouse_pid, children )
        bus[i] = new_bu
    end
    bus
end

function get_benefit_units(
    hh :: Household,
    allocator :: Function=default_bu_allocation ) :: BenefitUnits
    allocate_to_bus( allocator(hh))
end

# simple diagnostic prints for testing allocation

function printpids( pers :: Person )
    println( "pid: $(pers.pid) age $(pers.age) sex $(pers.sex) default_benefit_unit $(pers.default_benefit_unit) ")
end

function printpids( bu::BenefitUnit)
      head = get_head( bu )
      print( "HEAD: ")
      printpids( head )
      spouse = get_spouse( bu )
      if spouse != nothing
            print( "SPOUSE: ")
            printpids( spouse )

      end
      for chno in bu.children
            print( "CHILD_$(chno):")
            child = bu.people[chno]
            printpids( child )
      end
end

function printpids( people :: People_Dict )
    ks = sort(collect(keys(people)))
    for k in ks
        print( "PERSON_$k")
        printpids( people[k] )
    end
end

function printpids( buas::BUAllocation )
    nbus = size( buas )[1]
    for i in 1:nbus
        npeople = size( buas[i] )[1]
        println( "bu[$i]; num_people = $npeople:")
        for p in 1:npeople
            print( "PERS_$(p): ")
            printpids( buas[i][p])
        end
    end
end


text:
"""
# Temporal

Time series data analysis package for Julia.

Type and methods for time series data structures with the `TS` object.

Inspired by the pandas DataFrame in python and the xts package in R.

"""
module Temporal
    using Dates
    export
        # foundational
        TS, ts,
        # alteration
        rename!, rename,
        # combining
        ojoin, ijoin, ljoin, rjoin, merge, hcat, vcat, head, tail, subset,
        # missingness
        nanrows, nancols, dropnan, fillnan, fillnan!, ffill!, bfill!, linterp!,
        # operations
        TemporalBroadcastStyle,
        numfun, arrfun, operation,
        ones, zeros, trues, falses, isnan, countnz, sign, round,
        sum, mean, maximum, minimum, prod, cumsum, cumprod, diff, lag,
        shift, pct_change,
        # aggregation
        mondays, tuesdays, wednesdays, thursdays, fridays, saturdays, sundays,
        bow, eow, bom, eom, boq, eoq, boy, eoy, collapse, apply,
        # financial
        has_open, has_high, has_low, has_close, has_volume, is_ohlc, is_ohlcv,
        op, hi, lo, cl, vo, ohlc, ohlcv, hlc, hl, hl2, hlc3, ohlc4,
        # models
        acf,
        # data
        tsread, tswrite, quandl, quandl_auth, quandl_meta, quandl_search, yahoo, google
    include("data/ts.jl")
    include("data/indexing.jl")
    include("data/stringrange.jl")
    include("data/combine.jl")
    include("data/collapse.jl")
    include("data/filter.jl")
    include("util/alter.jl")
    include("util/utils.jl")
    include("util/ohlc.jl")
    include("util/convert.jl")
    include("util/show.jl")
    include("util/viz.jl")
    include("feed/utils.jl")
    include("feed/yahoo.jl")
    include("feed/google.jl")
    include("feed/quandl.jl")
    include("feed/text.jl")
    include("calc/operations.jl")
    include("calc/shift.jl")
    include("calc/models.jl")
end


text:
module BaseFilters

using Random
using Dates
using DataFrames
using Statistics

using AutoMLPipeline.Utils
using AutoMLPipeline.AbsTypes: Machine, Transformer, Learner, Workflow, Computer

import AutoMLPipeline.AbsTypes: fit!, transform!
export fit!,transform!
export OneHotEncoder

"""
    OneHotEncoder(Dict(
       # Nominal columns
       :nominal_columns => Int[],

       # Nominal column values map. Key is column index, value is list of
       # possible values for that column.
       :nominal_column_values_map => Dict{Int,Any}()
    ))

Transforms instances with nominal features into one-hot form
and coerces the instance matrix to be of element type Float64.
"""
mutable struct OneHotEncoder <: Transformer
  name::String
  model::Dict
  args::Dict

  function OneHotEncoder(args::Dict=Dict())
    default_args = Dict(
                        :name => "ohe",
                        # Nominal columns
                        :nominal_columns => Int[],
                        # Nominal column values map. Key is column index, value is list of
                        # possible values for that column.
                        :nominal_column_values_map => Dict{Int,Any}()
                       )
    cargs=nested_dict_merge(default_args,args)
    cargs[:name] = cargs[:name]*"_"*randstring(3)
    new(cargs[:name],Dict(),cargs)
  end
end

function fit!(ohe::OneHotEncoder, instances::DataFrame, labels::Vector=[])
  # Obtain nominal columns
  nominal_columns = ohe.args[:nominal_columns]
  if nominal_columns == Int[]
    nominal_columns,_ = find_catnum_columns(instances)
  end

  # Obtain unique values for each nominal column
  nominal_column_values_map = ohe.args[:nominal_column_values_map]
  if nominal_column_values_map == Dict{Int,Any}()
    for column in nominal_columns
      nominal_column_values_map[column] = unique(instances[:, column])
    end
  end

  # Create model
  ohe.model = Dict(
    :nominal_columns => nominal_columns,
    :nominal_column_values_map => nominal_column_values_map
  )
end

function transform!(ohe::OneHotEncoder, instances::DataFrame)
  nominal_columns = ohe.model[:nominal_columns]
  nominal_column_values_map = ohe.model[:nominal_column_values_map]

  # Create new transformed instance matrix of type Float64
  num_rows = size(instances, 1)
  num_columns = (size(instances, 2) - length(nominal_columns))
  if !isempty(nominal_column_values_map)
    num_columns += sum(map(x -> length(x), values(nominal_column_values_map)))
  end
  transformed_instances = zeros(Float64, num_rows, num_columns)

  # Fill transformed instance matrix
  col_start_index = 1
  for column in 1:size(instances, 2)
    if !in(column, nominal_columns)
      transformed_instances[:, col_start_index] = instances[:, column]
      col_start_index += 1
    else
      col_values = nominal_column_values_map[column]
      for row in 1:size(instances, 1)
        entry_value = instances[row, column]
        entry_value_index = findfirst(isequal(entry_value),col_values)
        if entry_value_index == 0
          warn("Unseen value found in OneHotEncoder,
                for entry ($row, $column) = $(entry_value).
                Patching value to $(col_values[1]).")
          entry_value_index = 1
        end
        entry_column = (col_start_index - 1) + entry_value_index
        transformed_instances[row, entry_column] = 1
      end
      col_start_index += length(nominal_column_values_map[column])
    end
  end

  return transformed_instances |> DataFrame
end

"""
    Imputer(
       Dict(
          # Imputation strategy.
          # Statistic that takes a vector such as mean or median.
          :strategy => mean
       )
    )

#Imputes NaN values from Float64 features.
"""
mutable struct Imputer <: Transformer
  name::String
  model::Dict
  args::Dict

  function Imputer(args=Dict())
    default_args = Dict(
      # Imputation strategy.
      # Statistic that takes a vector such as mean or median.
      :strategy => mean
    )
    cargs=nested_dict_merge(default_args,args)
    cargs[:name] = cargs[:name]*"_"*randstring(3)
    new(cargs[:name],Dict(),cargs)
  end
end

function fit!(imp::Imputer, instances::DataFrame, labels::Vector=[])
  imp.model = imp.args
end

function transform!(imp::Imputer, instances::DataFrame)
  new_instances = deepcopy(instances)
  strategy = imp.model[:strategy]

  for column in 1:size(instances, 2)
    column_values = instances[:, column]
    col_eltype = infer_eltype(column_values)

    if <:(col_eltype, Real)
      na_rows = map(x -> isnan(x), column_values)
      if any(na_rows)
        fill_value = strategy(column_values[.!na_rows])
        new_instances[na_rows, column] .= fill_value
      end
    end
  end
  return new_instances |> DataFrame
end

"""
    Wrapper(
       default_args = Dict(
          # Transformer to call.
          :transformer => OneHotEncoder(),
          # Transformer args.
          :transformer_args => nothing
       )
    )

Wraps around a TSML transformer.
"""
mutable struct Wrapper <: Transformer
  name::String
  model::Dict
  args::Dict

  function Wrapper(args=Dict())
    default_args = Dict(
      # Transformer to call.
      :transformer => OneHotEncoder(),
      # Transformer args.
      :transformer_args => nothing
    )
    cargs=nested_dict_merge(default_args,args)
    cargs[:name] = cargs[:name]*"_"*randstring(3)
    new(cargs[:name],Dict(),cargs)
  end
end

function fit!(wrapper::Wrapper, instances::DataFrame, labels::Vector)
  transformer_args = wrapper.args[:transformer_args]
  transformer = createtransformer(
    wrapper.args[:transformer],
    transformer_args
  )

  if transformer_args != nothing
    transformer_args = mergedict(transformer.args, transformer_args)
  end
  fit!(transformer, instances, labels)

  wrapper.model = Dict(
    :transformer => transformer,
    :transformer_args => transformer_args
  )
end

function transform!(wrapper::Wrapper, instances::DataFrame)
  transformer = wrapper.model[:transformer]
  return transform!(transformer, instances)
end

"""
    createtransformer(prototype::Transformer, args=nothing)

Create transformer

- `prototype`: prototype transformer to base new transformer on
- `options`: additional options to override prototype's options

Returns: new transformer.
"""
function createtransformer(prototype::Transformer, args=nothing)
  new_args = copy(prototype.args)
  if args != nothing
    new_args = mergedict(new_args, args)
  end

  prototype_type = typeof(prototype)
  return prototype_type(new_args)
end

end


text:
# Autogenerated wrapper script for NetCDF_jll for x86_64-linux-gnu
export libnetcdf

using HDF5_jll
using Zlib_jll
using XML2_jll
using LibCURL_jll
JLLWrappers.@generate_wrapper_header("NetCDF")
JLLWrappers.@declare_library_product(libnetcdf, "libnetcdf.so.19")
function __init__()
    JLLWrappers.@generate_init_header(HDF5_jll, Zlib_jll, XML2_jll, LibCURL_jll)
    JLLWrappers.@init_library_product(
        libnetcdf,
        "lib/libnetcdf.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()