text:
# Author: Ritchie Lee, ritchie.lee@sv.cmu.edu
# Date: 12/15/2014

# Correlated Encounter Model for Cooperative Aircraft in the National Airspace
# exposed as DBN.  Samples are generated at runtime at each step.

module CorrAEMDBNImpl

export
    AddObserver,

    getInitialState,
    initialize,
    update,
    get,

    CorrAEMDBN

import Compat.ASCIIString

using AbstractEncounterDBNImpl
using AbstractEncounterDBNInterfaces
using CommonInterfaces
using Util
using Encounter
using CorrAEMImpl
using RLESUtils, Observers

import CommonInterfaces.addObserver
import CommonInterfaces.initialize
import CommonInterfaces.update
import AbstractEncounterDBNInterfaces.get
import AbstractEncounterDBNInterfaces.getInitialState

import Base.convert

include(Pkg.dir("SISLES/src/Encounter/CorrAEMImpl/corr_aem_sample.jl"))

type CorrAEMDBN <: AbstractEncounterDBN

  number_of_aircraft::Int64

  encounter_file::ASCIIString
  initial_sample_file::ASCIIString
  transition_sample_file::ASCIIString

  encounter_number::Int64
  command_method::Symbol #:ENC = read from encounter file, :DBN = generate from DBN on-the-fly

  aem::CorrAEM
  dirichlet_transition

  #initial states
  init_aem_dstate::Vector{Int64} #discrete current state
  init_aem_dyn_cstate::Vector{Float64} #continuous of current dynamic variables

  #current states
  t::Int64
  aem_dstate::Vector{Int64} #discrete current state
  aem_dyn_cstate::Vector{Float64} #continuous of current dynamic variables

  #caching and reuse
  dynamic_variables0::Vector{Int64}
  dynamic_variables1::Vector{Int64}
  parents_cache::Dict{Int64, Vector{Bool}}
  weights_cache::Dict{Tuple{Int64, Int64}, Vector{Float64}}
  cumweights_cache::Dict{Tuple{Int64, Int64}, Vector{Float64}}

  #pre-allocated output to avoid repeating reallocations
  output_commands::Vector{CorrAEMCommand}

  logProb::Float64 #log probability of output

  function CorrAEMDBN(number_of_aircraft::Int, encounter_file::AbstractString, initial_sample_file::AbstractString,
                      transition_sample_file::AbstractString,
                      encounter_number::Int, encounter_seed::UInt64, command_method::Symbol)
    dbn = new()

    @assert number_of_aircraft == 2 #need to revisit the code if this is not true
    dbn.number_of_aircraft     = number_of_aircraft

    dbn.encounter_file         = encounter_file
    dbn.initial_sample_file    = initial_sample_file
    dbn.transition_sample_file = transition_sample_file
    dbn.encounter_number       = encounter_number
    dbn.command_method         = command_method

    dbn.aem = CorrAEM(encounter_file, initial_sample_file, transition_sample_file)

    dbn.t = 0

    srand(encounter_seed) #There's a rand inside generateEncounter
    generateEncounter(dbn.aem, sample_number=encounter_number) #To optimize: This allocates a lot of memory

    #compute initial states of variables
    p = dbn.aem.parameters
    dbn.dynamic_variables0 = p.temporal_map[:,1]
    dbn.dynamic_variables1 = p.temporal_map[:,2]

    aem_initial_unconverted = unconvertUnitsAemState(dbn.aem.initial)
    aem_initial_dstate = Int64[ val2ind(p.boundaries[i], p.r_transition[i], val)
                               for (i, val) in enumerate(aem_initial_unconverted)]
    dbn.init_aem_dstate = [aem_initial_dstate; aem_initial_dstate[dbn.dynamic_variables0]] #bins, [11:14] are updated with time, append space for t+1 variables
    dbn.init_aem_dyn_cstate = dbn.aem.initial[dbn.dynamic_variables0] #continuous variables.
    dbn.dirichlet_transition = bn_dirichlet_prior(p.N_transition)

    dbn.aem_dstate = deepcopy(dbn.init_aem_dstate)
    dbn.aem_dyn_cstate = deepcopy(dbn.init_aem_dyn_cstate)

    #precompute and cache these quantities
    dbn.parents_cache = Dict{Int64,Vector{Bool}}()
    dbn.weights_cache = Dict{Tuple{Int64,Int64}, Vector{Float64}}()
    dbn.cumweights_cache = Dict{Tuple{Int64,Int64}, Vector{Float64}}()

    for i = 1:length(p.N_transition)
      dbn.parents_cache[i] = p.G_transition[:, i]

      for j = 1:1:size(dbn.dirichlet_transition[i], 2)
        dbn.weights_cache[(i, j)] = p.N_transition[i][:, j] + dbn.dirichlet_transition[i][:, j]
        dbn.weights_cache[(i, j)] /= sum(dbn.weights_cache[(i, j)])
        dbn.cumweights_cache[(i, j)] = cumsum(dbn.weights_cache[(i, j)])
      end
    end

    dbn.output_commands = CorrAEMCommand[ CorrAEMCommand(0.0, 0.0, 0.0, 0.0) for i = 1:number_of_aircraft ]
    dbn.logProb = 0.0

    return dbn
  end
end

addObserver(dbn::CorrAEMDBN, f::Function) = add_observer(aem.observer, f)
addObserver(dbn::CorrAEMDBN, tag::AbstractString, f::Function) = add_observer(aem.observer, tag, f)

function initialize(dbn::CorrAEMDBN)
  #reset to initial state
  copy!(dbn.aem_dstate, dbn.init_aem_dstate)
  copy!(dbn.aem_dyn_cstate, dbn.init_aem_dyn_cstate)
  dbn.t = 0

  #reset aem indices
  initialize(dbn.aem)
end

function getInitialState(dbn::CorrAEMDBN, index::Int)
  return Encounter.getInitialState(dbn.aem, index)
end

function update(dbn::CorrAEMDBN)

  if dbn.command_method == :DBN
    logProb = step_dbn(dbn)
  elseif dbn.command_method == :ENC
    logProb = step_enc(dbn)
  else
    error("CorrAEMDBNImpl::Step: No such command method")
  end

  dbn.t += 1

  return logProb
end

function step_dbn(dbn::CorrAEMDBN)
  p = dbn.aem.parameters
  aem_dstate = dbn.aem_dstate #entire state, discrete bins
  aem_dyn_cstate = dbn.aem_dyn_cstate #dynamic states, continuous

  logProb = 0.0

  for (o,i) in enumerate(dbn.dynamic_variables1)
    if !isempty(find(dbn.parents_cache[i]))
      dims = tuple(p.r_transition[dbn.parents_cache[i]]...)
      j = sub2ind(dims, aem_dstate[dbn.parents_cache[i]]...)
      aem_dstate[i] = select_random_cumweights(dbn.cumweights_cache[(i,j)])
      logProb += log(dbn.weights_cache[(i,j)][aem_dstate[i]])
      #Resampling and dediscretizing process
      i_t = dbn.dynamic_variables0[o]
      if (aem_dstate[i] != aem_dstate[i_t]) || #compare to state at last time step, #Different bin, do resample
        (aem_dstate[i] == aem_dstate[i_t] && rand() < p.resample_rates[i_t]) #Same bin but meets resample rate
        aem_dyn_cstate[o],_ = dediscretize(aem_dstate[i],p.boundaries[i_t],p.zero_bins[i_t])
        if in(i,[17,18]) #these need unit conversion
          aem_dyn_cstate[o] /= 60 #convert units
        end
      end
      #Else same bin and does not meet rate, just set equal to previous (no update)
    end
  end

  # copy over x(t+1) to x(t)
  aem_dstate[dbn.dynamic_variables0] = aem_dstate[dbn.dynamic_variables1]

  #Just a reminder, this will break if number_of_aircraft != 2
  @assert dbn.number_of_aircraft == 2

  dbn.output_commands[1].t = dbn.t
  dbn.output_commands[1].v_d = getInitialSample(dbn.aem, :v1d)
  dbn.output_commands[1].h_d = dbn.aem_dyn_cstate[1]
  dbn.output_commands[1].psi_d = dbn.aem_dyn_cstate[3]

  dbn.output_commands[2].t = dbn.t
  dbn.output_commands[2].v_d = getInitialSample(dbn.aem, :v2d)
  dbn.output_commands[2].h_d = dbn.aem_dyn_cstate[2]
  dbn.output_commands[2].psi_d = dbn.aem_dyn_cstate[4]

  return dbn.logProb = logProb
end

#TODO: remove hardcoding
convert_units(x::AbstractFloat, i::Int) = in(i, [17,18]) ? x / 60 : x

convert(::Type{Vector{Float64}}, command_1::CorrAEMCommand, command_2::CorrAEMCommand) =
  [ command_1.h_d, command_2.h_d, command_1.psi_d, command_2.psi_d ]

function step_enc(dbn::CorrAEMDBN)
    aem = dbn.aem
    p = aem.parameters

    for i = 1:dbn.number_of_aircraft
        cmd = Encounter.update(aem, i)
        if cmd != nothing
            dbn.output_commands[i] = cmd
        end
    end

  #Just a reminder, this will break if number_of_aircraft != 2
  #@assert dbn.number_of_aircraft == 2

  #= FIXME: skip probability calc for now...
  #prepare t+1 from encounter commands
  aem_dyn_cstate = convert(Vector{Float64}, dbn.output_commands[1], dbn.output_commands[2])
  aem_dstate = dbn.aem_dstate #only copies pointer
  #load into the (t+1) slots
  #boundaries are specified at t
  aem_dyn_cstate_unconverted = unconvertUnitsDynVars(aem_dyn_cstate) #need to unconvert units
  aem_dstate[dbn.dynamic_variables1] = Int64[ val2ind(p.boundaries[i], p.r_transition[i], aem_dyn_cstate_unconverted[o])
                                                for (o, i) in enumerate(dbn.dynamic_variables0) ]

  # compute the probability of this transition
  logProb = 0.0
  for (o, i) in enumerate(dbn.dynamic_variables1)
    j = 1
    parents = dbn.parents_cache[i]
    if !isempty(find(parents))
      dims = tuple(p.r_transition[parents]...)
      indices = aem_dstate[parents]
      j = sub2ind(dims, indices...)
    end
    weights = dbn.weights_cache[(i, j)]
    logProb += log(weights[aem_dstate[i]])

    #probability from continuous sampling process
    #two components: resample prob, dediscretize prob
    i_prev = dbn.dynamic_variables0[o]
    if aem_dstate[i] != aem_dstate[i_prev] #not the same bin, resample wp 1
      logProb += dediscretize_prob(aem_dyn_cstate[o], aem_dstate[i], p.boundaries[i_prev], p.zero_bins[i_prev])
    elseif isapprox(aem_dyn_cstate[o], dbn.aem_dyn_cstate[o], atol=0.0001) #same bin same value, did not resample
      logProb += log(1.0 - p.resample_rates[i_prev])
    else #Same bin different value, got resampled
      logProb += log(p.resample_rates[i_prev])
      logProb += dediscretize_prob(aem_dyn_cstate[o], aem_dstate[i], p.boundaries[i_prev], p.zero_bins[i_prev])
    end
  end

  # copy over x(t+1) to x(t)
  aem_dstate[dbn.dynamic_variables0] = aem_dstate[dbn.dynamic_variables1]

  #push to sim
  dbn.aem_dstate = aem_dstate
  dbn.aem_dyn_cstate = aem_dyn_cstate

  #return
  dbn.logProb = logProb
  =#
  dbn.logProb = 0.0 #for now... FIXME
end

function get(dbn::CorrAEMDBN, aircraft_number::Int)
  return dbn.output_commands[aircraft_number]
end

function val2ind(boundariesi, ri, value)
  if !isempty(boundariesi)
    index = findfirst(x -> (x > value), boundariesi) - 1
    if index == -1
      index = ri
    end
  else
    index = value
  end
  return index
end

function dediscretize(dval::Int64, boundaries::Vector{Float64}, zero_bin::Int64)
  val_min = boundaries[dval]
  val_max = boundaries[dval + 1]

  if dval == zero_bin
    val = 0.0
    prob = 1.0
  elseif val_max == val_min
    val = val_min
    prob = 1.0
  else
    val = val_min +  rand() * (val_max - val_min)
    prob = 1.0 / (val_max - val_min)
    #this is a density so it won't be normalized to [0,1]
  end

  return (val, prob)
end

function dediscretize_prob(val::Float64, dval::Int64, boundaries::Vector{Float64}, zero_bin::Int64)
  val_min = boundaries[dval]
  val_max = boundaries[dval + 1]

  if dval == zero_bin
    @assert val == 0.0
    prob = 1.0
  elseif val_max == val_min
    @assert val == val_min
    prob = 1.0
  else
    @assert val_min <= val <= val_max
    prob = 1.0 / (val_max - val_min)
    #this is a density so it won't be normalized to [0,1]
  end
  return prob
end

function unconvertUnitsDynVars(v)
  return [v[1:2] * 60.0, v[3:end]]
end

function unconvertUnitsAemState(state_) #from CorrAEM
  state = deepcopy(state_)

  state[7] /= 1.68780
  state[8] /= 1.68780
  state[9] /= 1.68780
  state[10] /= 1.68780
  state[11] *= 60
  state[12] *= 60
  state[15] /= 6076.12

  return state
end

function select_random_cumweights(cweights::Vector{Float64})
#select randomly according to cumulative weights vector

  r = cweights[end] * rand()

  return findfirst(x -> (x >= r), cweights)
end

end #module

ext: jl

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


text:
function isprerelease(version::VersionNumber)::Bool
    prerelease = version.prerelease
    isempty(prerelease) && return false
    for x in prerelease
        if !isempty(strip(x))
            return true
        end
    end
    return false
end

function _calculate_increment(before::VersionNumber,
                              after::VersionNumber)::VersionNumber
    before_first3 = VersionNumber(before.major,
                                  before.minor,
                                  before.patch)
    after_first3 = VersionNumber(after.major,
                                 after.minor,
                                 after.patch)
    # @debug("before: ", before)
    # @debug("before_first3: ", before_first3)
    # @debug("after:", after)
    # @debug("after_first3:", after_first3)
    always_assert(after > before, "after > before")
    always_assert(after_first3 >= before_first3,
                  "after_first3 >= before_first3")
    if before.major == after.major
        if before.minor == after.minor
            always_assert(after.patch >= before.patch,
                          "after.patch >= before.patch")
            return VersionNumber(0, 0, after.patch - before.patch)
        else
            always_assert(after.minor >= before.minor,
                          "after.minor >= before.minor")
            return VersionNumber(0,
                                 after.minor - before.minor,
                                 after.patch - 0)
        end
    else
        always_assert(after.major >= before.major,
                      "after.major >= before.major")
        return VersionNumber(after.major - before.major,
                             after.minor - 0,
                             after.patch - 0)
    end
end

function check_version_increment(master_version::VersionNumber,
                                 head_version::VersionNumber;
                                 allow_skipping_versions::Bool,
                                 gh_set_output::Bool = get(ENV, "GITHUB_ACTIONS", "") == "true",
                                 gh_set_output_io::IO = stdout)::Nothing
    always_assert(head_version > master_version, "head_version > master_version")
    _gh_set_output_println(gh_set_output, gh_set_output_io, "compare_versions", "success")
    @info("Version number has increased")
    increment = _calculate_increment(master_version, head_version)
    if increment in (v"0.0.0", v"0.0.1", v"0.1.0", v"1.0.0")
        @info("Increment is good",
              master_version,
              head_version,
              increment)
    else
        if allow_skipping_versions
            @warn("Increment is bad, but `allow_skipping_versions` is true, so we will allow it",
                  master_version,
                  head_version,
                  increment,
                  allow_skipping_versions)
        else
            @error("Increment is bad",
                  master_version,
                  head_version,
                  increment,
                  allow_skipping_versions)
            error("Bad increment")
        end
    end
    return nothing
end

function compare_versions(master_version::VersionNumber,
                          head_version::VersionNumber;
                          allow_unchanged_prerelease::Bool,
                          allow_skipping_versions::Bool,
                          gh_set_output::Bool = get(ENV, "GITHUB_ACTIONS", "") == "true",
                          gh_set_output_io::IO = stdout)::Nothing
    if head_version > master_version
        check_version_increment(master_version,
                                head_version;
                                allow_skipping_versions = allow_skipping_versions,
                                gh_set_output = gh_set_output,
                                gh_set_output_io = gh_set_output_io)
    elseif head_version == master_version
        if isprerelease(head_version) && isprerelease(master_version) && allow_unchanged_prerelease
            _gh_set_output_println(gh_set_output, gh_set_output_io, "compare_versions", "success")
            @info("Version number did not change, but it is a prerelease so this is allowed")
        else
            _gh_set_output_println(gh_set_output, gh_set_output_io, "compare_versions", "failure")
            throw(ErrorException("Version number is unchanged, which is not allowed"))
        end
    else
        _gh_set_output_println(gh_set_output, gh_set_output_io, "compare_versions", "failure")
        throw(ErrorException("Version number decreased, which is not allowed"))
    end
    return nothing
end


text:
##Operators
# TODO: REMOVE!
for op in (:Derivative,:Integral)
    @eval begin
        function ($op)(d::AbstractVector{T}) where T<:IntervalOrSegment
            n=length(d)
            R=zeros(Operator{mapreduce(eltype,promote_type,d)},n,n)
            for k=1:n
                R[k,k]=$op(d[k])
            end

            R
        end
    end
end

function Evaluation(d::AbstractVector{T},x...) where T<:IntervalOrSegment
    n=length(d)
    R=zeros(Operator{mapreduce(eltype,promote_type,d)},n,n)
    for k=1:n
        R[k,k]=Evaluation(d[k],x...)
    end

    R
end

## Construction
function diagm_container(kv::Pair{<:Integer,<:AbstractVector{O}}...) where O<:Operator
    T = mapreduce(x -> mapreduce(eltype,promote_type,x.second),
                  promote_type, kv)
    n = mapreduce(x -> length(x.second) + abs(x.first), max, kv)
    zeros(Operator{T}, n, n)
end

##TODO: unify with other blockdiag
function blockdiag(d1::AbstractVector{T},d2::AbstractVector{T}) where T<:Operator
  if isempty(d1)&&isempty(d2)
    error("Empty blockdiag")
  end
  if isempty(d1)
    TT=mapreduce(eltype,promote_type,d2)
  elseif isempty(d2)
    TT=mapreduce(eltype,promote_type,d1)
  else
     TT=promote_type(mapreduce(eltype,promote_type,d1),
                    mapreduce(eltype,promote_type,d2))
  end

  D=zeros(Operator{TT},length(d1)+length(d2),2)
  D[1:length(d1),1]=d1
  D[length(d1)+1:end,2]=d2
  D
end

blockdiag(a::Operator,b::Operator) = blockdiag(Operator{promote_type(eltype(a),eltype(b))}[a],
                                                Operator{promote_type(eltype(a),eltype(b))}[b])

## broadcase

broadcast(::typeof(*),A::AbstractArray{N},D::Operator) where {N<:Number} =
    Operator{promote_type(N,eltype(D))}[A[k,j]*D for k=1:size(A,1),j=1:size(A,2)]
broadcast(::typeof(*),D::Operator,A::AbstractArray{N}) where {N<:Number}=A.*D


text:
using Test

using SpinMonteCarlo

const SEED = 137
const SEED2 = 19937
const MCS = 10000
const Therm = MCS
const alpha = 0.001

@testset begin
    filenames = [
                 "classical.jl",
                 "quantum.jl",
                 "checkpoint.jl",
                ]
    for filename in filenames
        t = @elapsed include(filename)
        println("$(filename): $t sec")
    end
end


text:
"""
    ListBasedNonProjective()

Transition system for list-based non-projective dependency parsing.

Described in Nivre 2008, "Algorithms for Deterministic Incremental Dependency Parsing."
"""
struct ListBasedNonProjective <: AbstractTransitionSystem end

initconfig(s::ListBasedNonProjective, graph::DependencyTree) =
    ListBasedNonProjectiveConfig(graph)
initconfig(s::ListBasedNonProjective, deptype, words) =
    ListBasedNonProjectiveConfig{deptype}(words)

projective_only(::ListBasedNonProjective) = false

transition_space(::ListBasedNonProjective, labels=[]) =
    isempty(labels) ? [LeftArc(), RightArc(), NoArc(), Shift()] :
    [LeftArc.(labels)..., RightArc.(labels)..., NoArc(), Shift()]

struct ListBasedNonProjectiveConfig{T} <: AbstractParserConfiguration{T}
    λ1::Vector{Int} # right-headed
    λ2::Vector{Int} # left-headed
    β::Vector{Int}
    A::Vector{T}
end

function ListBasedNonProjectiveConfig{T}(words::Vector{String}) where {T}
    λ1 = [0]
    λ2 = Int[]
    β = 1:length(words)
    A = [unk(T, id, w) for (id,w) in enumerate(words)]
    ListBasedNonProjectiveConfig{T}(λ1, λ2, β, A)
end

function ListBasedNonProjectiveConfig{T}(gold::DependencyTree) where {T}
    λ1 = [0]
    λ2 = Int[]
    β = 1:length(gold)
    A = [dep(token, head=-1) for token in gold]
    ListBasedNonProjectiveConfig{T}(λ1, λ2, β, A)
end
ListBasedNonProjectiveConfig(gold::DependencyTree) =
    ListBasedNonProjectiveConfig{eltype(gold)}(gold)

buffer(cfg::ListBasedNonProjectiveConfig) = cfg.β

token(cfg::ListBasedNonProjectiveConfig, i) = iszero(i) ? root(deptype(cfg)) :
                                              i == -1   ? noval(deptype(cfg)) :
                                              cfg.A[i]
tokens(cfg::ListBasedNonProjectiveConfig) = cfg.A
tokens(cfg::ListBasedNonProjectiveConfig, is) = [token(cfg, i) for i in is]

function leftarc(cfg::ListBasedNonProjectiveConfig, args...; kwargs...)
    λ1, i = cfg.λ1[1:end-1], cfg.λ1[end]
    j, β = cfg.β[1], cfg.β[2:end]
    A = copy(cfg.A)
    i != 0 && (A[i] = dep(A[i], args...; head=j, kwargs...))
    ListBasedNonProjectiveConfig(λ1, [i ; cfg.λ2], [j ; β], A)
end

function rightarc(cfg::ListBasedNonProjectiveConfig, args...; kwargs...)
    λ1, i = cfg.λ1[1:end-1], cfg.λ1[end]
    j, β = cfg.β[1], cfg.β[2:end]
    A = copy(cfg.A)
    A[j] = dep(A[j], args...; head=i, kwargs...)
    ListBasedNonProjectiveConfig(λ1, [i ; cfg.λ2], [j ; β], A)
end

function noarc(cfg::ListBasedNonProjectiveConfig)
    λ1, i = cfg.λ1[1:end-1], cfg.λ1[end]
    λ2, β, A = cfg.λ2, cfg.β, cfg.A
    ListBasedNonProjectiveConfig(λ1, [i ; λ2], β, A)
end

function shift(cfg::ListBasedNonProjectiveConfig)
    λ1, λ2 = cfg.λ1, cfg.λ2
    i, β = cfg.β[1], cfg.β[2:end]
    ListBasedNonProjectiveConfig([λ1 ; λ2 ; i], Int[], β, cfg.A)
end

function isfinal(cfg::ListBasedNonProjectiveConfig)
    return all(a -> head(a) != -1, tokens(cfg)) && length(cfg.λ1) == length(cfg.A) + 1 &&
        length(cfg.λ2) == 0 && length(cfg.β) == 0
end

"""
    static_oracle(::ListBasedNonProjectiveConfig, tree)

Return a training oracle function which returns gold transition
operations from a parser configuration with reference to `graph`.
"""
function static_oracle(cfg::ListBasedNonProjectiveConfig, tree, arc=untyped)
    l = i -> arc(tree[i])
    if length(cfg.λ1) >= 1 && length(cfg.β) >= 1
        i, λ1 = cfg.λ1[end], cfg.λ1[1:end-1]
        j, β = cfg.β[1], cfg.β[2:end]
        if !iszero(i) && head(tree, i) == j
            return LeftArc(l(i)...)
        elseif head(tree, j) == i
            return RightArc(l(j)...)
        end
        j_deps = dependents(tree, j)
        if (!(any(x -> x < j, j_deps) && j_deps[1] < i)) && !(head(tree, j) < i)
            return Shift()
        end
    end
    if length(cfg.λ1) == 0
        return Shift()
    end
    return NoArc()
end

# todo?
possible_transitions(cfg::ListBasedNonProjectiveConfig, graph::DependencyTree, arc=untyped) =
    TransitionOperator[static_oracle(cfg, graph, arc)]

==(cfg1::ListBasedNonProjectiveConfig, cfg2::ListBasedNonProjectiveConfig) =
    cfg1.λ1 == cfg2.λ1 && cfg1.λ2 == cfg2.λ2 && cfg1.β == cfg2.β && cfg1.A == cfg2.A

function Base.show(io::IO, c::ListBasedNonProjectiveConfig)
    λ1 = join(c.λ1, ",")
    λ2 = join(c.λ2, ",")
    β = join(c.β, ",")
    print(io, "ListBasedNonProjectiveConfig([$λ1],[$λ2],[$β])\n$(join([join([id(t),form(t),head(t)],'\t') for t in tokens(c)],'\n'))")
end


text:
# Basic test script - generates the file sample_uam_traj.csv containing one uam trajectory

# Include the files
include("UAMTrajectoryGenerator.jl")
# Set the random seed
Random.seed!(1)
# Generate the trajectory file
generate_trajectory_file("output/uam_traj.csv")

max_stars_repo_name: Airspace-Encounter-Models/em-model-uam

max_issues_repo_name: Airspace-Encounter-Models/em-model-uam

max_forks_repo_name: Airspace-Encounter-Models/em-model-uam


text:
using FunctionalCollections: PersistentSet, PersistentHashMap, dissoc, assoc, conj, disj
using Gen

#########################################
# compiling a trace into a factor graph #
#########################################

struct Latent{T,U}
    domain::Vector{T}
    parent_addrs::Vector{U}
end

struct Observation{U}
    parent_addrs::Vector{U} # only include addrs that are selected
end

parent_addrs(info::Union{Latent,Observation}) = info.parent_addrs

function get_domain_to_idx(domain::Vector{T}) where {T}
    domain_to_idx = Dict{T,Int}()
    for (i, value) in enumerate(domain)
        domain_to_idx[value] = i
    end
    return domain_to_idx
end

function cartesian_product(value_lists)
    tuples = Vector{Tuple}()
    for value in value_lists[1]
        if length(value_lists) > 1
            append!(tuples,
                [(value, rest...) for rest in cartesian_product(value_lists[2:end])])
        else
            append!(tuples, [(value,)])
        end
    end
    return tuples
end

# addr could be latent or obserrved...
# latent_addrs is the set of latent variables that are involved, which may or
# may not include the actual addr itself..

function create_factor(
        trace, addr,
        latents::Dict{Any,Latent}, observations::Dict{Any,Observation},
        all_latent_addrs::Vector{Any})
    N = length(all_latent_addrs)
    in_factor = Vector{Bool}(undef, N)
    if haskey(latents, addr)
        for (i, a) in enumerate(all_latent_addrs)
            in_factor[i] = (a == addr || a in parent_addrs(latents[addr]))
        end
        num_vars = length(parent_addrs(latents[addr]))+1
    elseif haskey(observations, addr)
        for (i, a) in enumerate(all_latent_addrs)
            in_factor[i] = (a in parent_addrs(observations[addr]))
        end
        num_vars = length(parent_addrs(observations[addr]))
    end
    dims = map(i -> in_factor[i] ? length(latents[all_latent_addrs[i]].domain) : 1, 1:N)
    log_factor = Array{Float64,N}(undef, dims...)
    view_inds = map(i -> in_factor[i] ? Colon() : 1, 1:N)
    log_factor_view = view(log_factor, view_inds...)
    var_addrs = Vector{Any}(undef, num_vars)
    value_idx_lists = Vector{Any}(undef, num_vars)
    j = 1
    for i in 1:N
        if in_factor[i]
            a = all_latent_addrs[i]
            var_addrs[j] = a
            value_idx_lists[j] = collect(1:length(latents[a].domain))
            j += 1
        end
    end
    @assert j == num_vars + 1

    # populate factor with values by probing trace with update
    # the key idea is that this scales exponentially in maximum number of
    # parents of a variable, not the total number of variables

    cprod = cartesian_product(value_idx_lists)
    for value_idx_tuple in cprod
        choices = choicemap()
        for (a, value_idx) in zip(var_addrs, value_idx_tuple)
            choices[a] = latents[a].domain[value_idx]
        end
        (tmp_trace, _, _, _) = update(trace, get_args(trace), map((_)->NoChange(),get_args(trace)), choices)
        # NOTE: technically, the generative function of trace can use any
        # internal proposal, not only forward sampling, but this code is only
        # correct if the internal proposal uses forward sampling. enhancing the
        # trace interface with some more methods that specifically assume a
        # dependency graph would resolve this
        weight = project(tmp_trace, select(addr))
        log_factor_view[value_idx_tuple...] = weight
    end

    log_factor = log_factor .- logsumexp(log_factor[:])

    return (log_factor, var_addrs)
end

function compile_trace_to_factor_graph(
        trace, latents::Dict{Any,Latent}, observations::Dict{Any,Observation})

    # choose order of addresses (note, this is NOT the elimination order)
    # TODO does the order in which the addresses are indexed matter for e.g. cache performance? maybe?
    all_latent_addrs = collect(keys(latents))
    latent_addr_to_idx = Dict{Any,Int}()
    for (idx, addr) in enumerate(all_latent_addrs)
        latent_addr_to_idx[addr] = idx
    end

    # construct factor nodes, one for each latent and downstream variable
    N = length(all_latent_addrs)
    addr_to_factor_node = Dict{Any,FactorNode{N}}()
    factor_id = 1
    for addr in Iterators.flatten((keys(latents), keys(observations)))
        (log_factor, var_addrs) = create_factor(
            trace, addr, latents, observations, all_latent_addrs)
        addr_to_factor_node[addr] = FactorNode{N}(factor_id, Int[latent_addr_to_idx[a] for a in var_addrs], log_factor)
        factor_id += 1
    end
    num_factors = factor_id - 1

    # for each latent address, the set of addresses for factors that it is involved in
    children_and_self = [Set{Any}([all_latent_addrs[i]]) for i in 1:N]
    for (addr, addr_info) in Iterators.flatten((latents, observations))
        for parent_addr in parent_addrs(addr_info)
            push!(children_and_self[latent_addr_to_idx[parent_addr]], addr)
        end
    end

    # construct factor graph
    var_nodes = PersistentHashMap{Int,VarNode}()
    for (addr, addr_info) in latents
        i = latent_addr_to_idx[addr]
        factor_nodes = PersistentSet{FactorNode{N}}(
            [addr_to_factor_node[addr] for addr in children_and_self[i]])
        var_node = VarNode(addr, factor_nodes, addr_info.domain, get_domain_to_idx(addr_info.domain))
        var_nodes = assoc(var_nodes, latent_addr_to_idx[addr], var_node)
    end
    return FactorGraph{N}(num_factors, var_nodes, latent_addr_to_idx)
end


text:
#MAIN ENTRY POINT
#final columns for the decompressed file

#merges crsp and compustat
# note refresh merge will be automatically run if any of the other two are run
function loadccm(;
  datapath::String=DATA_PATH,
  ccmname::String = CCM_NAME,
  incsvstream::Function = IN_CSV_STREAM,
  csvextension::String = CSV_EXTENSION)

  local ccm::DataFrame

  #makes a serialization object so we don't have to keep parsing the CSV
  ccm = incsvstream("$datapath\\$ccmname.$csvextension") |> CSV.File |> DataFrame

  return ccm

end

function mergecompccm!(comp::DataFrame;
    testoutput::Bool = TEST_OUTPUT,
    months2stale::Int =  MONTHS_2_STALE_LONG)

  #first
  local ccm::DataFrame = loadccm()
  ccm = preprocessccm!(ccm)

  testoutput && CSV.write("output\\ccm.csv", ccm)

  select!(ccm, Not([:cusip]))
  #println("Comp rows before merge: $(size(comp, 1))")
  comp = innerjoin(comp, ccm, on=[:gvkey])
  comp.keep = trues(size(comp,1))

  comp.keep = (r::DataFrameRow->
    r.adate ∈ r.linkeffdate:Day(1):r.linkenddate).(eachrow(comp))

  comp = comp[comp.keep,:]
  select!(comp, Not(:keep))
  println("Comp rows after merge, post-filter: $(size(comp, 1))")

  #Note: we still will have overlapping lineffdates and linkenddates, which will be corrected
  #with the below method
  comp = reconcilecompccmintervals!(comp)

  #  lagwithin!(comp, :fdate, :gvkey, :fyear)
  sort!(comp, [:lpermno, :fdate])
  issorted(comp, [:lpermno, :adate]) || error("I thought comp was sorted by [:lpermno, :adate]")
  period2stale::Month = Month(months2stale)
  lagwithin2!(comp, [:fdate, :fyrmonth], :lpermno, date=:fdate, maxnotstale = period2stale)
  lagwithin2!(comp, [:Lfdate, :Lfyrmonth], :lpermno, date=:fdate, maxnotstale = period2stale)
  for r ∈ eachrow(comp)
    ismissing(r.Lfdate) && continue
    (r.Lfdate < r.fdate -  period2stale) && error("failed date check $r")
  end
  @info "past date validatation check"

  #=
  NOTE: we accept dup fyears due to changing fiscal years
  function testfordups(df::AbstractDataFrame, groupfields::Vector{Symbol})
    for sdf ∈ groupby(df, groupfields)
      if size(sdf,1) > 1
        println(sdf[!, [:gvkey, :lpermno, :fdate,
          :Lfdate, :adate, :Nadate, :linkeffdate, :linkenddate, :begindate, :enddate, :fyear]])
        error("duplicate fyear found")
      end
    end
  @info "passed dedup fyear test"
  end
  testfordups(comp, [:lpermno, :fyear])=#

  return comp
end

#formats the ccm table in a reasonable way
function preprocessccm!(ccm::DataFrame;
  retainedcolsccm::Vector{Symbol} = RETAINED_COLS_CCM,
  lastdatestring::String = LAST_DATE_STRING,
  testoutput::Bool = TEST_OUTPUT)

  #names ot lower case
  cleannames::Vector{String} = (s::String->lowercase(s)).(names(ccm))
  rename!(ccm, ((oldn, newn)->oldn=>newn).(names(ccm), cleannames))

  #fix the end date
  ccm.linkenddt .= (s::MString->
    coalesce(s,"")=="E" ? lastdatestring : s).(ccm.linkenddt)

  #below didn't work for some reason
  #filter!(r::DataFrameRow->((!ismissing(r.gvkey)) && (!ismissing(r.lpermno))), ccm)
  ccm = ccm[(r::DataFrameRow->(
    (!ismissing(r.gvkey)) && (!ismissing(r.lpermno)))).(eachrow(ccm)),:]
  ccm.linkprim = (s::String->length(s)==1 ? s[1] : error("invalid linkprim in ccm")).(ccm.linkprim)
  filter!(r::DataFrameRow->r.linkprim =='P' || r.linkprim =='C', ccm)

  #parse the dates
  ccm.linkeffdate = (s->parseccm(Date, s)).(ccm.linkdt)
  ccm.linkenddate = (s->parseccm(Date, s)).(ccm.linkenddt)

  #dedup if possible
  ccm.tokeep = trues(size(ccm,1))
  ccm.dateranges = ((linkeffdate::Date,linkenddate::Date)->
    linkeffdate:Day(1):linkenddate).(ccm.linkeffdate,ccm.linkenddate)

  gccm::GroupedDataFrame = groupby(ccm, [:gvkey, :lpermno])
  @mpar for i ∈ 1:length(gccm)
    sccm::SubDataFrame = gccm[i]
    Nsccm::Int = size(sccm,1)

    #construct the union of all of the date sets
    if Nsccm >  1
      local periodrange::StepRange{Date,Day} = minimum(
        sccm.linkeffdate):Day(1):maximum(sccm.linkenddate)
      local dateunionlength::Int = length(unique([sccm.dateranges...;]))

      if dateunionlength == length(periodrange)
        sccm.linkeffdate[1] = minimum(periodrange)
        sccm.linkenddate[1] = maximum(periodrange)
        sccm.tokeep[2:end] .= false
      end

    end
  end

  println("ccm size pre-dedup $(size(ccm))")
  ccm = ccm[ccm.tokeep,:]
  println("ccm size post-dedup $(size(ccm))")

  if testoutput
    testoutput && CSV.write("output\\ccm_preproc.csv", ccm)
  end
  #rename!(ccm, :indfmt=>:indfmtold)
  #ccm.indfmt = (s::MString->parsecomp(Symbol, s)).(ccm.indfmtold)

  #only need some of the fields
  ccm = select!(ccm, retainedcolsccm)

  return ccm
end

#this function reconciles the link intervals with the announcement intervals
#NOTE: A weird situation can arise if gvkey changes. The the linkenddate finishes too soon,
#we end up with a gap . Considering the laternatives, it seems like the lesser evil
#to extending the record beyond the linkenddate
function reconcilecompccmintervals!(comp::DataFrame)
  #subtract one day since the effective interval is inclusive of the
  #first date while the begin-end date is exclusve of the first date
  comp.begindate .= (max).(comp.begindate, comp.linkeffdate .- Day(1))
  comp.enddate .= (min).(comp.enddate, comp.linkenddate)

  #this gets rid of one day records
  #(this assumes the previous effective record is the current record)
  #not many of these ~140
  #println("records before interval reconciliation: $(size(comp,1))")
  #println(comp[comp.begindate .≥ comp.enddate,
  #  [:gvkey, :adateq, :begindate, :enddate, :linkeffdate, :linkenddate]])
  comp = comp[comp.begindate .< comp.enddate, :]
  println("records after interval reconciliation: $(size(comp,1))")

  return comp
end

function parseccm(::Type{T}, s::String;
    #=parsedmissings = CCM_PARSED_MISSINGS=#) where T

  #if it doesn't parse to the right type, set to missing
  v::Union{T,Missing} = something(tryparse(T, s), missing)

  #check against the list of missing codes
  #(!ismissing(v)) && (v ∈ parsedmissings) && (v=missing)

  return v
end

#the date case
parseccm(::Type{Date}, s::String;
    ccmdateformat::DateFormat = CCM_DATE_FORMAT) = Dates.Date(s, ccmdateformat)

#helper methods and special cases
parseccm(::Type{<:Any}, v::Missing) = missing
parseccm(::Type{Date}, i::Int) = parseccm(Date, "$i")
parseccm(::Type{Symbol}, s::String) = Symbol(s)

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
module ComradeSoss
#Turn off precompilations because of GG bug https://github.com/cscherrer/Soss.jl/issues/267
__precompile__(false)
using HypercubeTransform

using Reexport
@reexport using Soss
@reexport using Comrade

import Distributions
const Dists = Distributions
using MeasureTheory
using NamedTupleTools
using NestedTuples
using MacroTools

using PyCall
using Random
using Requires
using ParameterHandling
using StatsBase: median
using StatsBase
using StructArrays
using TupleVectors

# This is a hack for MeasureTheory since it wants the type to output
Base.rand(rng::AbstractRNG, ::Type{Float64}, d::Dists.ContinuousDistribution) = rand(rng, d)

const rad2μas = 180.0/π*3600*1e6

#These hold the pointers to the PyObjects that will store dynesty and ehtim
#and are loaded at initialization
const dynesty = PyNULL()

export ObsChar, scandata,
       create_joint,
       DynestyStatic,
       sample, optimize,
       threaded_optimize,
       chi2, ehtim

#include("ehtim.jl")
include("utility.jl")
include("hypercube.jl")
include("inference.jl")
include("dists.jl")
include("models.jl")
#include("grads.jl")
#include("nlopt.jl")
#include("hmc.jl")

function __init__()
    copy!(dynesty, pyimport("dynesty"))
    @require UltraNest="6822f173-b0be-4018-9ee2-28bf56348d09" include("ultranest.jl")
    @require NestedSamplers="41ceaf6f-1696-4a54-9b49-2e7a9ec3782e" include("nested.jl")
    @require BlackBoxOptim="a134a8b2-14d6-55f6-9291-3336d3ab0209" include("bboptim.jl")
    @require Metaheuristics="bcdb8e00-2c21-11e9-3065-2b553b22f898" include("metaheuristics.jl")
end

end #end ComradeSoss


text:
#############################################################################
# Joulia
# A Large-Scale Spatial Power System Model for Julia
# See https://github.com/JuliaEnergy/Joulia.jl
#############################################################################
# This file contains functions and structs for nodes

mutable struct Nodes
	id
	load
	exchange

	function Nodes(nodes_df::DataFrame, load_df::DataFrame, exchange_df::DataFrame)
		N = Symbol.(nodes_df[1])
		load_dict = df_to_dict(load_df)
		insert_default!(load_dict, fill(0.0, 8760), N)

		exchange_dict = df_to_dict(exchange_df)
		insert_default!(exchange_dict, fill(0.0, 8760), N)

		return new(N, load_dict, exchange_dict)
	end
end


text:
@with_kw type Body2D <: Body
name::String = "None"
parent::Body = ground2D()
connection::connection2D = free2D("None")
mass::SymFloat = SymFloat(symbols("m_None",nonnegative=true,real=true))
inertia::SymFloat = SymFloat(symbols("I_None",nonnegative=true,real=true))
end
function ==(b1::Body2D,b2::Body2D)
 f = fieldnames(b1)
 for i in f
   if getfield(b1,i) != getfield(b2,i)
     return false
   end
 end
 return true
end


text:
"""
    child = replace(parent, query => replacement)

Generates a `child` crystal structure by (i) searches the `parent` crystal structure for subgraphs that match the `query` then
(ii) replaces the substructures of the `parent` matching the `query` fragment with the `replacement` fragment.

Equivalent to calling `substructure_replace(query ∈ parent, replacement)`.

Accepts the same keyword arguments as [`substructure_replace`](@ref).
"""
replace(p::Crystal, pair::Pair; kwargs...) = substructure_replace(pair[1] ∈ p, pair[2]; kwargs...)

"""
    alignment = Alignment(rot::Matrix{Float64}, shift_1::Vector{Float64}, shift_2::Vector{Float64}, err::Float64)

Data structure for tracking alignment in substructure find/replace operations.
"""
struct Alignment
    rot::Matrix{Float64}
    # before rotation
    shift_1::Vector{Float64}
    # after rotation
    shift_2::Vector{Float64}
    # error
    err::Float64
end

struct Installation
    aligned_replacement::Crystal
    q2p::Dict{Int, Int}
    r2p::Dict{Int, Int}
end

function get_r2p_alignment(replacement::Crystal, parent::Crystal, r2p::Dict{Int, Int}, q2p::Dict{Int, Int})
    center = (X::Matrix{Float64}) -> sum(X, dims=2)[:] / size(X, 2)
    # when both centered to origin
    @assert replacement.atoms.n ≥ 3 && parent.atoms.n ≥ 3 "Parent and replacement must each be at least 3 atoms for SVD alignment."
    ###
    #   compute centered Cartesian coords of the atoms of
    #       replacement fragment involved in alignment
    ###
    atoms_r = Cart(replacement.atoms[[r for (r, p) in r2p]], replacement.box)
    X_r = atoms_r.coords.x
    x_r_center = center(X_r)
    X_r = X_r .- x_r_center

    ###
    #   compute centered Cartesian coords of the atoms of
    #       parent involved in alignment
    ###
    # handle fragments cut across the PB using the parent subset isomorphic to query
    parent_substructure = deepcopy(parent[[p for (q, p) in q2p]])
    conglomerate!(parent_substructure) # must do this for when replacement fragment is disconnected
    # prepare parent substructure having correspondence with replacement
    p2ps = Dict([p => i for (i, p) in enumerate([p for (q, p) in q2p])]) # parent to parent subset map

    parent_substructure_to_align_to = parent_substructure[[p2ps[p] for p in [p for (r, p) in r2p]]]
    atoms_p = Cart(parent_substructure_to_align_to.atoms, parent.box)
    X_p = atoms_p.coords.x
    x_p_center = center(X_p)
    X_p = X_p .- x_p_center

    # solve the orthogonal procrustes probelm via SVD
    F = svd(X_r * X_p')
    # optimal rotation matrix
    rot =  F.V * F.U'

    err = norm(rot * X_r - X_p)

    return Alignment(rot, - x_r_center, x_p_center, err)
end

function conglomerate!(parent_substructure::Crystal)
    # snip the cross-PB bonds
    bonds = deepcopy(parent_substructure.bonds)
    if length(connected_components(bonds)) > 1
        @warn "# connected components in parent substructure > 1. assuming the substructure does not cross the periodic boundary..."
        return
    end
    drop_cross_pb_bonds!(bonds)

    # find connected components of bonding graph without cross-PB bonds
    #    these are the components split across the boundary
    conn_comps = connected_components(bonds)

    # if substructure is entireline in the unit cell, it's already conglomerated :)
    if length(conn_comps) == 1
        return
    end

    # we wish to shift all connected components to a reference component,
    #   defined to be the largest component for speed.
    conn_comps_shifted = [false for c = 1:length(conn_comps)]
    ref_comp_id = argmax(length.(conn_comps))
    conn_comps_shifted[ref_comp_id] = true  # consider it shifted.

    # has atom p been shifted?
    function shifted_atom(p::Int)
        # loop over all connected components that have been shifted
        for conn_comp in conn_comps[conn_comps_shifted]
            # if parent substructure atom in this, yes!
            if p in conn_comp
                return true
            end
        end
        # reached this far, atom p is not in component that has been shifted.
        return false
    end

    # to which component does atom p belong?
    function find_component(p::Int)
        for c = 1:length(conn_comps)
            if p in conn_comps[c]
                return c
            end
        end
    end

    # until all components have been shifted to the reference component...
    while ! all(conn_comps_shifted)
        # loop over cross-PB edges in the parent substructure
        for ed in edges(parent_substructure.bonds)
            if get_prop(parent_substructure.bonds, ed, :cross_boundary)
                # if one edge belongs to unshifted component and another belogs to any component that has been shifted...
                if     shifted_atom(ed.src) && ! shifted_atom(ed.dst)
                    p_ref, p = ed.src, ed.dst
                elseif shifted_atom(ed.dst) && ! shifted_atom(ed.src)
                    p_ref, p = ed.dst, ed.src
                else
                    continue # both are shifted or both are unshifted. ignore this cross-PB edge
                end
                # here's the unshifted component we will shift next, to be next to the shifted components.
                comp_id = find_component(p)
                # find displacement vector for this cross-PB edge.
                dx = parent_substructure.atoms.coords.xf[:, p_ref] - parent_substructure.atoms.coords.xf[:, p]
                # get distance to nearest image
                n_dx = copy(dx)
                nearest_image!(n_dx)
                # shift all atoms in this component by this vector.
                for atom_idx in conn_comps[comp_id]
                    parent_substructure.atoms.coords.xf[:, atom_idx] .+= dx - n_dx
                end
                # mark that we've shifted this component.
                conn_comps_shifted[comp_id] = true
            end
        end
    end
    return
end

function aligned_replacement(replacement::Crystal, parent::Crystal, r2p_alignment::Alignment)
    # put replacement into cartesian space
    atoms_r = Cart(replacement.atoms, replacement.box)
    # rotate replacement to align with parent_subset
    atoms_r.coords.x[:, :] = r2p_alignment.rot * (atoms_r.coords.x .+ r2p_alignment.shift_1) .+ r2p_alignment.shift_2
    # cast atoms back to Frac
    return Crystal(replacement.name, parent.box, Frac(atoms_r, parent.box), Charges{Frac}(0), replacement.bonds, replacement.symmetry)
end

function effect_replacements(search::Search, replacement::Crystal, configs::Vector{Tuple{Int, Int}}, name::String)::Crystal
    nb_not_masked = sum(.! occursin.(rc[:r_tag], String.(search.query.atoms.species)))
    if replacement.atoms.n > 0
        q_unmasked_in_r = substructure_search(search.query[1:nb_not_masked], replacement)
        q2r = Dict([q => q_unmasked_in_r.isomorphisms[1][1][q] for q in 1:nb_not_masked])
    else
        q2r = Dict{Int, Int}()
    end

    installations = [optimal_replacement(search, replacement, q2r, loc_id, [ori_id]) for (loc_id, ori_id) in configs]

    child = install_replacements(search.parent, installations, name)

    # handle `missing` values in edge :cross_boundary attribute
    for edge in edges(child.bonds) # loop over edges
        # check if cross-boundary info is missing
        if ismissing(get_prop(child.bonds, edge, :cross_boundary))
            # check if bond crosses boundary
            distance_e = get_prop(child.bonds, edge, :distance) # distance in edge property
            dxa = Cart(Frac(child.atoms.coords.xf[:, src(edge)] - child.atoms.coords.xf[:, dst(edge)]), child.box) # Cartesian displacement
            distance_a = norm(dxa.x) # current euclidean distance by atom coords
            set_prop!(child.bonds, edge, :cross_boundary, !isapprox(distance_e, distance_a, atol=0.1))
        end
    end

    return child
end

function install_replacements(parent::Crystal, replacements::Vector{Installation}, name::String)::Crystal
    # create child w/o symmetry rules for sake of crystal addition
    child = Crystal(name, parent.box, parent.atoms, parent.charges, parent.bonds, Xtals.SymmetryInfo())

    obsolete_atoms = Int[] # to delete at the end

    # loop over replacements to install
    for installation in replacements
        replacement, q2p, r2p = installation.aligned_replacement, installation.q2p, installation.r2p
        #add into parent
        if replacement.atoms.n > 0
            child = +(child, replacement, check_overlap=false)
        end

        # reconstruct bonds
        for (r, p) in r2p # p is in parent_subst
            p_nbrs = neighbors(parent.bonds, p)
            for p_nbr in p_nbrs
                if ! (p_nbr in values(q2p)) # p_nbr not in parent_subst
                    # need bond nbr => r in child, where r is in replacement
                    e = (p_nbr, child.atoms.n - replacement.atoms.n + r)
                    # create edge
                    add_edge!(child.bonds, e)
                    # copy edge attributes from parent (:cross_boundary will need to be reassessed later)
                    set_props!(child.bonds, e[1], e[2], props(parent.bonds, p, p_nbr))
                    set_prop!(child.bonds, e[1], e[2], :cross_boundary, missing)
                end
            end
        end

        # accumulate atoms to delete
        obsolete_atoms = vcat(obsolete_atoms, values(q2p)...)
    end

    # delete obsolete atoms
    obsolete_atoms = unique(obsolete_atoms)
    keep_atoms = [p for p = 1:child.atoms.n if ! (p in obsolete_atoms)]
    child = child[keep_atoms]

    # restore symmetry rules
    child = Crystal(name, child.box, child.atoms, child.charges, child.bonds, parent.symmetry)

    # return result
    return child
end

function optimal_replacement(search::Search, replacement::Crystal, q2r::Dict{Int,Int}, loc_id::Int, ori_ids::Vector{Int})
    # unpack search arg
    isomorphisms, parent = search.isomorphisms, search.parent

    if q2r == Dict{Int, Int}() # "replace-with-nothing" operation
        q2p = isomorphisms[loc_id][1]
        r2p = Dict([0 => p for p in values(q2p)])
        return Installation(replacement, q2p, r2p)
    end

    if ori_ids == [0]
        ori_ids = [1:nb_ori_at_loc(search)[loc_id]...]
    end

    # loop over ori_ids to find best r2p_alignment
    r2p_alignment = Alignment(zeros(1,1), [0.], [0.], Inf)
    best_ori = 0
    best_r2p = Dict{Int, Int}()
    for ori_id in ori_ids
        # find r2p isom
        q2p = isomorphisms[loc_id][ori_id]
        r2p = Dict([r => q2p[q] for (q, r) in q2r])
        # calculate alignment
        test_alignment = get_r2p_alignment(replacement, parent, r2p, q2p)
        # keep best alignment and generating ori_id
        if test_alignment.err < r2p_alignment.err
            r2p_alignment = test_alignment
            best_ori = ori_id
            best_r2p = r2p
        end
    end

    opt_aligned_replacement = aligned_replacement(replacement, parent, r2p_alignment)

    # return the replacement modified according to r2p_alignment
    @assert ne(opt_aligned_replacement.bonds) == ne(replacement.bonds)
    return Installation(opt_aligned_replacement, isomorphisms[loc_id][best_ori], best_r2p)
end

@doc raw"""
    child = substructure_replace(search, replacement; random=false, nb_loc=0, loc=Int[], ori=Int[], name="new_xtal", verbose=false, remove_duplicates=false, periodic_boundaries=true)

Replace the substructures of `search.parent` matching the `search.query` fragment with the `replacement` fragment,
at locations and orientations specified by the keyword arguments `random`, `nb_loc`, `loc`, and `ori`.
Default behavior is to effect replacements at all "hit" locations in the parent structure and, at each location,
choose the orientation giving the optimal (lowest error) spatial aligment.

Returns a new `Crystal` with the specified modifications (returns `search.parent` if no replacements are made).

# Arguments
- `search::Search` the `Search` for a substructure moiety in the parent crystal
- `replacement::Crystal` the moiety to use for replacement of the searched substructure
- `random::Bool` set `true` to select random replacement orientations
- `nb_loc::Int` assign a value to select random replacement at `nb_loc` random locations
- `loc::Array{Int}` assign value(s) to select specific locations for replacement.  If `ori` is not specified, replacement orientation is random.
- `ori::Array{Int}` assign value(s) when `loc` is assigned to specify exact configurations for replacement. `0` values mean the configuration at that location should be selected for optimal alignment with the parent.
- `name::String` assign to give the generated `Crystal` a name ("new_xtal" by default)
- `verbose::Bool` set `true` to print console messages about the replacement(s) being performed
- `remove_duplicates::Bool` set `true` to automatically combine overlapping atoms of the same species in generated structures.
- `reinfer_bonds::Bool` set `true` to re-infer bonds after producing a structure
- `periodic_boundaries::Bool` set `false` to disable periodic boundary conditions when checking for atom duplication or re-inferring bonds
"""
function substructure_replace(search::Search, replacement::Crystal; random::Bool=false,
    nb_loc::Int=0, loc::Array{Int}=Int[], ori::Array{Int}=Int[], name::String="new_xtal", verbose::Bool=false,
    remove_duplicates::Bool=false, periodic_boundaries::Bool=true, reinfer_bonds::Bool=false, wrap::Bool=true)::Crystal
    # replacement at all locations (default)
    if nb_loc == 0 && loc == Int[] && ori == Int[]
        nb_loc = nb_locations(search)
        loc = [1:nb_loc...]
        if random
            ori = [rand(1:nb_ori_at_loc(search)[i]) for i in loc]
            if verbose
                @info "Replacing" q_in_p=search r=replacement.name mode="random ori @ all loc"
            end
        else
            ori = zeros(Int, nb_loc)
            if verbose
                @info "Replacing" q_in_p=search r=replacement.name mode="optimal ori @ all loc"
            end
        end
    # replacement at nb_loc random locations
    elseif nb_loc > 0 && ori == Int[] && loc == Int[]
        loc = sample([1:nb_locations(search)...], nb_loc, replace=false)
        if random
            ori = [rand(1:nb_ori_at_loc(search)[i]) for i in loc]
            if verbose
                @info "Replacing" q_in_p=search r=replacement.name mode="random ori @ $nb_loc loc"
            end
        else
            ori = zeros(Int, nb_loc)
            if verbose
                @info "Replacing" q_in_p=search r=replacement.name mode="optimal ori @ $nb_loc loc"
            end
        end
    # specific replacements
    elseif ori ≠ Int[] && loc ≠ Int[]
        @assert length(loc) == length(ori) "one orientation per location"
        nb_loc = length(ori)
        if verbose
            @info "Replacing" q_in_p=search r=replacement.name mode="loc: $loc\tori: $ori"
        end
    # replacement at specific locations
    elseif loc ≠ Int[]
        nb_loc = length(loc)
        if random
            ori = [rand(1:nb_ori_at_loc(search)[i]) for i in loc]
            if verbose
                @info "Replacing" q_in_p=search r=replacement.name mode="random ori @ loc: $loc"
            end
        else
            ori = zeros(Int, nb_loc)
            if verbose
                @info "Replacing" q_in_p=search r=replacement.name mode="optimal ori @ loc: $loc"
            end
        end
    end

    # remove charges from parent
    if search.parent.charges.n > 0
        @warn "Dropping charges from parent."
        p = search.parent
        search = Search(
            Crystal(p.name, p.box, p.atoms, Charges{Frac}(0), p.bonds, p.symmetry),
            search.query,
            search.isomorphisms
        )
    end

    # generate configuration tuples (location, orientation)
    configs = Tuple{Int,Int}[(loc[i], ori[i]) for i in 1:nb_loc]
    # process replacements
    child = effect_replacements(search, replacement, configs, name)

    if remove_duplicates
        child = Crystal(child.name, child.box,
            Xtals.remove_duplicates(child.atoms, child.box, periodic_boundaries),
            Xtals.remove_duplicates(child.charges, child.box, periodic_boundaries)
        )
    end

    if wrap
        # throw error if installed replacement fragment spans the unit cell
        if any(abs.(child.atoms.coords.xf) .> 2.0)
            error("installed replacement fragment too large for the unit cell; replicate the parent and try again.")
        end

        # wrap coordinates
        wrap!(child.atoms.coords)

        # check :cross_boundary edge attributes
        for edge in edges(child.bonds) # loop over edges
            distance_e = get_prop(child.bonds, edge, :distance) # distance in edge property
            dxa = Cart(Frac(child.atoms.coords.xf[:, src(edge)] - child.atoms.coords.xf[:, dst(edge)]), child.box) # Cartesian displacement
            distance_a = norm(dxa.x) # current euclidean distance by atom coords
            set_prop!(child.bonds, edge, :cross_boundary, !isapprox(distance_e, distance_a, atol=0.1))
        end
    end

    if reinfer_bonds
        remove_bonds!(child)
        infer_bonds!(child, periodic_boundaries)
    end

    return child
end

substructure_replace(search::Search, replacement::Nothing; kwargs...) = substructure_replace(search, moiety(nothing); kwargs...)


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antdalert

"""
    'feffery'_antdalert(;kwargs...)

An AntdAlert component.

Keyword arguments:
- `id` (String; optional)
- `className` (String; optional)
- `closable` (Bool; optional)
- `description` (String; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `message` (String | Array of Strings; optional)
- `messageRenderMode` (a value equal to: 'default', 'loop-text', 'marquee'; optional)
- `showIcon` (Bool; optional)
- `style` (Dict; optional)
- `type` (a value equal to: 'success', 'info', 'warning', 'error'; optional)
"""
function 'feffery'_antdalert(; kwargs...)
        available_props = Symbol[:id, :className, :closable, :description, :loading_state, :message, :messageRenderMode, :showIcon, :style, :type]
        wild_props = Symbol[]
        return Component("'feffery'_antdalert", "AntdAlert", "feffery_antd_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: ReilYoung/feffery-antd-components

max_issues_repo_name: ReilYoung/feffery-antd-components

max_forks_repo_name: ReilYoung/feffery-antd-components


text:
module ATI

import GetC.@getCFun

typealias GLenum Cuint
typealias GLboolean Cuchar
typealias GLbitfield Cuint
typealias GLvoid Void
typealias GLbyte Cuchar
typealias GLshort Cshort
typealias GLint Cint
typealias GLubyte Cuchar
typealias GLushort Cushort
typealias GLuint Cuint
typealias GLsizei Csize_t
typealias GLfloat Cfloat
typealias GLclampf Cfloat
typealias GLdouble Cdouble
typealias GLclampd Cdouble
typealias GLchar Cchar
typealias GLint64 Clonglong
typealias GLuint64 Culonglong
typealias GLhalf Cushort
typealias GLhalfARB Cushort
typealias GLhalfNV Cushort
typealias GLsync Ptr{Void}
typealias Pointer Ptr{Void}
typealias GLsizeiptr Cint
typealias GLintptr Cptrdiff_t
const DRAW_BUFFER0_ATI = 0x8825
export DRAW_BUFFER0_ATI
const DRAW_BUFFER10_ATI = 0x882F
export DRAW_BUFFER10_ATI
const DRAW_BUFFER11_ATI = 0x8830
export DRAW_BUFFER11_ATI
const DRAW_BUFFER12_ATI = 0x8831
export DRAW_BUFFER12_ATI
const DRAW_BUFFER13_ATI = 0x8832
export DRAW_BUFFER13_ATI
const DRAW_BUFFER14_ATI = 0x8833
export DRAW_BUFFER14_ATI
const DRAW_BUFFER15_ATI = 0x8834
export DRAW_BUFFER15_ATI
const DRAW_BUFFER1_ATI = 0x8826
export DRAW_BUFFER1_ATI
const DRAW_BUFFER2_ATI = 0x8827
export DRAW_BUFFER2_ATI
const DRAW_BUFFER3_ATI = 0x8828
export DRAW_BUFFER3_ATI
const DRAW_BUFFER4_ATI = 0x8829
export DRAW_BUFFER4_ATI
const DRAW_BUFFER5_ATI = 0x882A
export DRAW_BUFFER5_ATI
const DRAW_BUFFER6_ATI = 0x882B
export DRAW_BUFFER6_ATI
const DRAW_BUFFER7_ATI = 0x882C
export DRAW_BUFFER7_ATI
const DRAW_BUFFER8_ATI = 0x882D
export DRAW_BUFFER8_ATI
const DRAW_BUFFER9_ATI = 0x882E
export DRAW_BUFFER9_ATI
const MAX_DRAW_BUFFERS_ATI = 0x8824
export MAX_DRAW_BUFFERS_ATI
const ELEMENT_ARRAY_ATI = 0x8768
export ELEMENT_ARRAY_ATI
const ELEMENT_ARRAY_POINTER_ATI = 0x876A
export ELEMENT_ARRAY_POINTER_ATI
const ELEMENT_ARRAY_TYPE_ATI = 0x8769
export ELEMENT_ARRAY_TYPE_ATI
const BUMP_ENVMAP_ATI = 0x877B
export BUMP_ENVMAP_ATI
const BUMP_NUM_TEX_UNITS_ATI = 0x8777
export BUMP_NUM_TEX_UNITS_ATI
const BUMP_ROT_MATRIX_ATI = 0x8775
export BUMP_ROT_MATRIX_ATI
const BUMP_ROT_MATRIX_SIZE_ATI = 0x8776
export BUMP_ROT_MATRIX_SIZE_ATI
const BUMP_TARGET_ATI = 0x877C
export BUMP_TARGET_ATI
const BUMP_TEX_UNITS_ATI = 0x8778
export BUMP_TEX_UNITS_ATI
const DU8DV8_ATI = 0x877A
export DU8DV8_ATI
const DUDV_ATI = 0x8779
export DUDV_ATI
const X2X_BIT_ATI = 0x00000001
export X2X_BIT_ATI
const X4X_BIT_ATI = 0x00000002
export X4X_BIT_ATI
const X8X_BIT_ATI = 0x00000004
export X8X_BIT_ATI
const ADD_ATI = 0x8963
export ADD_ATI
const BIAS_BIT_ATI = 0x00000008
export BIAS_BIT_ATI
const BLUE_BIT_ATI = 0x00000004
export BLUE_BIT_ATI
const CND0_ATI = 0x896B
export CND0_ATI
const CND_ATI = 0x896A
export CND_ATI
const COLOR_ALPHA_PAIRING_ATI = 0x8975
export COLOR_ALPHA_PAIRING_ATI
const COMP_BIT_ATI = 0x00000002
export COMP_BIT_ATI
const CON_0_ATI = 0x8941
export CON_0_ATI
const CON_10_ATI = 0x894B
export CON_10_ATI
const CON_11_ATI = 0x894C
export CON_11_ATI
const CON_12_ATI = 0x894D
export CON_12_ATI
const CON_13_ATI = 0x894E
export CON_13_ATI
const CON_14_ATI = 0x894F
export CON_14_ATI
const CON_15_ATI = 0x8950
export CON_15_ATI
const CON_16_ATI = 0x8951
export CON_16_ATI
const CON_17_ATI = 0x8952
export CON_17_ATI
const CON_18_ATI = 0x8953
export CON_18_ATI
const CON_19_ATI = 0x8954
export CON_19_ATI
const CON_1_ATI = 0x8942
export CON_1_ATI
const CON_20_ATI = 0x8955
export CON_20_ATI
const CON_21_ATI = 0x8956
export CON_21_ATI
const CON_22_ATI = 0x8957
export CON_22_ATI
const CON_23_ATI = 0x8958
export CON_23_ATI
const CON_24_ATI = 0x8959
export CON_24_ATI
const CON_25_ATI = 0x895A
export CON_25_ATI
const CON_26_ATI = 0x895B
export CON_26_ATI
const CON_27_ATI = 0x895C
export CON_27_ATI
const CON_28_ATI = 0x895D
export CON_28_ATI
const CON_29_ATI = 0x895E
export CON_29_ATI
const CON_2_ATI = 0x8943
export CON_2_ATI
const CON_30_ATI = 0x895F
export CON_30_ATI
const CON_31_ATI = 0x8960
export CON_31_ATI
const CON_3_ATI = 0x8944
export CON_3_ATI
const CON_4_ATI = 0x8945
export CON_4_ATI
const CON_5_ATI = 0x8946
export CON_5_ATI
const CON_6_ATI = 0x8947
export CON_6_ATI
const CON_7_ATI = 0x8948
export CON_7_ATI
const CON_8_ATI = 0x8949
export CON_8_ATI
const CON_9_ATI = 0x894A
export CON_9_ATI
const DOT2_ADD_ATI = 0x896C
export DOT2_ADD_ATI
const DOT3_ATI = 0x8966
export DOT3_ATI
const DOT4_ATI = 0x8967
export DOT4_ATI
const EIGHTH_BIT_ATI = 0x00000020
export EIGHTH_BIT_ATI
const FRAGMENT_SHADER_ATI = 0x8920
export FRAGMENT_SHADER_ATI
const GREEN_BIT_ATI = 0x00000002
export GREEN_BIT_ATI
const HALF_BIT_ATI = 0x00000008
export HALF_BIT_ATI
const LERP_ATI = 0x8969
export LERP_ATI
const MAD_ATI = 0x8968
export MAD_ATI
const MOV_ATI = 0x8961
export MOV_ATI
const MUL_ATI = 0x8964
export MUL_ATI
const NEGATE_BIT_ATI = 0x00000004
export NEGATE_BIT_ATI
const NUM_FRAGMENT_CONSTANTS_ATI = 0x896F
export NUM_FRAGMENT_CONSTANTS_ATI
const NUM_FRAGMENT_REGISTERS_ATI = 0x896E
export NUM_FRAGMENT_REGISTERS_ATI
const NUM_INPUT_INTERPOLATOR_COMPONENTS_ATI = 0x8973
export NUM_INPUT_INTERPOLATOR_COMPONENTS_ATI
const NUM_INSTRUCTIONS_PER_PASS_ATI = 0x8971
export NUM_INSTRUCTIONS_PER_PASS_ATI
const NUM_INSTRUCTIONS_TOTAL_ATI = 0x8972
export NUM_INSTRUCTIONS_TOTAL_ATI
const NUM_LOOPBACK_COMPONENTS_ATI = 0x8974
export NUM_LOOPBACK_COMPONENTS_ATI
const NUM_PASSES_ATI = 0x8970
export NUM_PASSES_ATI
const QUARTER_BIT_ATI = 0x00000010
export QUARTER_BIT_ATI
const RED_BIT_ATI = 0x00000001
export RED_BIT_ATI
const REG_0_ATI = 0x8921
export REG_0_ATI
const REG_10_ATI = 0x892B
export REG_10_ATI
const REG_11_ATI = 0x892C
export REG_11_ATI
const REG_12_ATI = 0x892D
export REG_12_ATI
const REG_13_ATI = 0x892E
export REG_13_ATI
const REG_14_ATI = 0x892F
export REG_14_ATI
const REG_15_ATI = 0x8930
export REG_15_ATI
const REG_16_ATI = 0x8931
export REG_16_ATI
const REG_17_ATI = 0x8932
export REG_17_ATI
const REG_18_ATI = 0x8933
export REG_18_ATI
const REG_19_ATI = 0x8934
export REG_19_ATI
const REG_1_ATI = 0x8922
export REG_1_ATI
const REG_20_ATI = 0x8935
export REG_20_ATI
const REG_21_ATI = 0x8936
export REG_21_ATI
const REG_22_ATI = 0x8937
export REG_22_ATI
const REG_23_ATI = 0x8938
export REG_23_ATI
const REG_24_ATI = 0x8939
export REG_24_ATI
const REG_25_ATI = 0x893A
export REG_25_ATI
const REG_26_ATI = 0x893B
export REG_26_ATI
const REG_27_ATI = 0x893C
export REG_27_ATI
const REG_28_ATI = 0x893D
export REG_28_ATI
const REG_29_ATI = 0x893E
export REG_29_ATI
const REG_2_ATI = 0x8923
export REG_2_ATI
const REG_30_ATI = 0x893F
export REG_30_ATI
const REG_31_ATI = 0x8940
export REG_31_ATI
const REG_3_ATI = 0x8924
export REG_3_ATI
const REG_4_ATI = 0x8925
export REG_4_ATI
const REG_5_ATI = 0x8926
export REG_5_ATI
const REG_6_ATI = 0x8927
export REG_6_ATI
const REG_7_ATI = 0x8928
export REG_7_ATI
const REG_8_ATI = 0x8929
export REG_8_ATI
const REG_9_ATI = 0x892A
export REG_9_ATI
const SATURATE_BIT_ATI = 0x00000040
export SATURATE_BIT_ATI
const SECONDARY_INTERPOLATOR_ATI = 0x896D
export SECONDARY_INTERPOLATOR_ATI
const SUB_ATI = 0x8965
export SUB_ATI
const SWIZZLE_STQ_ATI = 0x8977
export SWIZZLE_STQ_ATI
const SWIZZLE_STQ_DQ_ATI = 0x8979
export SWIZZLE_STQ_DQ_ATI
const SWIZZLE_STRQ_ATI = 0x897A
export SWIZZLE_STRQ_ATI
const SWIZZLE_STRQ_DQ_ATI = 0x897B
export SWIZZLE_STRQ_DQ_ATI
const SWIZZLE_STR_ATI = 0x8976
export SWIZZLE_STR_ATI
const SWIZZLE_STR_DR_ATI = 0x8978
export SWIZZLE_STR_DR_ATI
const RENDERBUFFER_FREE_MEMORY_ATI = 0x87FD
export RENDERBUFFER_FREE_MEMORY_ATI
const TEXTURE_FREE_MEMORY_ATI = 0x87FC
export TEXTURE_FREE_MEMORY_ATI
const VBO_FREE_MEMORY_ATI = 0x87FB
export VBO_FREE_MEMORY_ATI
const COLOR_CLEAR_UNCLAMPED_VALUE_ATI = 0x8835
export COLOR_CLEAR_UNCLAMPED_VALUE_ATI
const RGBA_FLOAT_MODE_ATI = 0x8820
export RGBA_FLOAT_MODE_ATI
const MAX_PN_TRIANGLES_TESSELATION_LEVEL_ATI = 0x87F1
export MAX_PN_TRIANGLES_TESSELATION_LEVEL_ATI
const PN_TRIANGLES_ATI = 0x87F0
export PN_TRIANGLES_ATI
const PN_TRIANGLES_NORMAL_MODE_ATI = 0x87F3
export PN_TRIANGLES_NORMAL_MODE_ATI
const PN_TRIANGLES_NORMAL_MODE_LINEAR_ATI = 0x87F7
export PN_TRIANGLES_NORMAL_MODE_LINEAR_ATI
const PN_TRIANGLES_NORMAL_MODE_QUADRATIC_ATI = 0x87F8
export PN_TRIANGLES_NORMAL_MODE_QUADRATIC_ATI
const PN_TRIANGLES_POINT_MODE_ATI = 0x87F2
export PN_TRIANGLES_POINT_MODE_ATI
const PN_TRIANGLES_POINT_MODE_CUBIC_ATI = 0x87F6
export PN_TRIANGLES_POINT_MODE_CUBIC_ATI
const PN_TRIANGLES_POINT_MODE_LINEAR_ATI = 0x87F5
export PN_TRIANGLES_POINT_MODE_LINEAR_ATI
const PN_TRIANGLES_TESSELATION_LEVEL_ATI = 0x87F4
export PN_TRIANGLES_TESSELATION_LEVEL_ATI
const STENCIL_BACK_FAIL_ATI = 0x8801
export STENCIL_BACK_FAIL_ATI
const STENCIL_BACK_FUNC_ATI = 0x8800
export STENCIL_BACK_FUNC_ATI
const STENCIL_BACK_PASS_DEPTH_FAIL_ATI = 0x8802
export STENCIL_BACK_PASS_DEPTH_FAIL_ATI
const STENCIL_BACK_PASS_DEPTH_PASS_ATI = 0x8803
export STENCIL_BACK_PASS_DEPTH_PASS_ATI
const TEXT_FRAGMENT_SHADER_ATI = 0x8200
export TEXT_FRAGMENT_SHADER_ATI
const MODULATE_ADD_ATI = 0x8744
export MODULATE_ADD_ATI
const MODULATE_SIGNED_ADD_ATI = 0x8745
export MODULATE_SIGNED_ADD_ATI
const MODULATE_SUBTRACT_ATI = 0x8746
export MODULATE_SUBTRACT_ATI
const ALPHA_FLOAT16_ATI = 0x881C
export ALPHA_FLOAT16_ATI
const ALPHA_FLOAT32_ATI = 0x8816
export ALPHA_FLOAT32_ATI
const INTENSITY_FLOAT16_ATI = 0x881D
export INTENSITY_FLOAT16_ATI
const INTENSITY_FLOAT32_ATI = 0x8817
export INTENSITY_FLOAT32_ATI
const LUMINANCE_ALPHA_FLOAT16_ATI = 0x881F
export LUMINANCE_ALPHA_FLOAT16_ATI
const LUMINANCE_ALPHA_FLOAT32_ATI = 0x8819
export LUMINANCE_ALPHA_FLOAT32_ATI
const LUMINANCE_FLOAT16_ATI = 0x881E
export LUMINANCE_FLOAT16_ATI
const LUMINANCE_FLOAT32_ATI = 0x8818
export LUMINANCE_FLOAT32_ATI
const RGBA_FLOAT16_ATI = 0x881A
export RGBA_FLOAT16_ATI
const RGBA_FLOAT32_ATI = 0x8814
export RGBA_FLOAT32_ATI
const RGB_FLOAT16_ATI = 0x881B
export RGB_FLOAT16_ATI
const RGB_FLOAT32_ATI = 0x8815
export RGB_FLOAT32_ATI
const MIRROR_CLAMP_ATI = 0x8742
export MIRROR_CLAMP_ATI
const MIRROR_CLAMP_TO_EDGE_ATI = 0x8743
export MIRROR_CLAMP_TO_EDGE_ATI
const ARRAY_OBJECT_BUFFER_ATI = 0x8766
export ARRAY_OBJECT_BUFFER_ATI
const ARRAY_OBJECT_OFFSET_ATI = 0x8767
export ARRAY_OBJECT_OFFSET_ATI
const DISCARD_ATI = 0x8763
export DISCARD_ATI
const DYNAMIC_ATI = 0x8761
export DYNAMIC_ATI
const OBJECT_BUFFER_SIZE_ATI = 0x8764
export OBJECT_BUFFER_SIZE_ATI
const OBJECT_BUFFER_USAGE_ATI = 0x8765
export OBJECT_BUFFER_USAGE_ATI
const PRESERVE_ATI = 0x8762
export PRESERVE_ATI
const STATIC_ATI = 0x8760
export STATIC_ATI
const MAX_VERTEX_STREAMS_ATI = 0x876B
export MAX_VERTEX_STREAMS_ATI
const VERTEX_SOURCE_ATI = 0x8774
export VERTEX_SOURCE_ATI
const VERTEX_STREAM0_ATI = 0x876C
export VERTEX_STREAM0_ATI
const VERTEX_STREAM1_ATI = 0x876D
export VERTEX_STREAM1_ATI
const VERTEX_STREAM2_ATI = 0x876E
export VERTEX_STREAM2_ATI
const VERTEX_STREAM3_ATI = 0x876F
export VERTEX_STREAM3_ATI
const VERTEX_STREAM4_ATI = 0x8770
export VERTEX_STREAM4_ATI
const VERTEX_STREAM5_ATI = 0x8771
export VERTEX_STREAM5_ATI
const VERTEX_STREAM6_ATI = 0x8772
export VERTEX_STREAM6_ATI
const VERTEX_STREAM7_ATI = 0x8773
export VERTEX_STREAM7_ATI
@getCFun "libGL" glDrawBuffersATI glDrawBuffersATI(n::GLsizei, bufs::Ptr{GLenum})::Void
export glDrawBuffersATI
@getCFun "libGL" glElementPointerATI glElementPointerATI(type_::GLenum, pointer::Ptr{Void})::Void
export glElementPointerATI
@getCFun "libGL" glDrawElementArrayATI glDrawElementArrayATI(mode::GLenum, count::GLsizei)::Void
export glDrawElementArrayATI
@getCFun "libGL" glDrawRangeElementArrayATI glDrawRangeElementArrayATI(mode::GLenum, start::GLuint, END::GLuint, count::GLsizei)::Void
export glDrawRangeElementArrayATI
@getCFun "libGL" glTexBumpParameterivATI glTexBumpParameterivATI(pname::GLenum, param::Ptr{GLint})::Void
export glTexBumpParameterivATI
@getCFun "libGL" glTexBumpParameterfvATI glTexBumpParameterfvATI(pname::GLenum, param::Ptr{GLfloat})::Void
export glTexBumpParameterfvATI
@getCFun "libGL" glGetTexBumpParameterivATI glGetTexBumpParameterivATI(pname::GLenum, param::Ptr{GLint})::Void
export glGetTexBumpParameterivATI
@getCFun "libGL" glGetTexBumpParameterfvATI glGetTexBumpParameterfvATI(pname::GLenum, param::Ptr{GLfloat})::Void
export glGetTexBumpParameterfvATI
@getCFun "libGL" glGenFragmentShadersATI glGenFragmentShadersATI(range_::GLuint)::Cuint
export glGenFragmentShadersATI
@getCFun "libGL" glBindFragmentShaderATI glBindFragmentShaderATI(id::GLuint)::Void
export glBindFragmentShaderATI
@getCFun "libGL" glDeleteFragmentShaderATI glDeleteFragmentShaderATI(id::GLuint)::Void
export glDeleteFragmentShaderATI
@getCFun "libGL" glBeginFragmentShaderATI glBeginFragmentShaderATI()::Void
export glBeginFragmentShaderATI
@getCFun "libGL" glEndFragmentShaderATI glEndFragmentShaderATI()::Void
export glEndFragmentShaderATI
@getCFun "libGL" glPassTexCoordATI glPassTexCoordATI(dst::GLuint, coord::GLuint, swizzle::GLenum)::Void
export glPassTexCoordATI
@getCFun "libGL" glSampleMapATI glSampleMapATI(dst::GLuint, interp::GLuint, swizzle::GLenum)::Void
export glSampleMapATI
@getCFun "libGL" glColorFragmentOp1ATI glColorFragmentOp1ATI(op::GLenum, dst::GLuint, dstMask::GLuint, dstMod::GLuint, arg1::GLuint, arg1Rep::GLuint, arg1Mod::GLuint)::Void
export glColorFragmentOp1ATI
@getCFun "libGL" glColorFragmentOp2ATI glColorFragmentOp2ATI(op::GLenum, dst::GLuint, dstMask::GLuint, dstMod::GLuint, arg1::GLuint, arg1Rep::GLuint, arg1Mod::GLuint, arg2::GLuint, arg2Rep::GLuint, arg2Mod::GLuint)::Void
export glColorFragmentOp2ATI
@getCFun "libGL" glColorFragmentOp3ATI glColorFragmentOp3ATI(op::GLenum, dst::GLuint, dstMask::GLuint, dstMod::GLuint, arg1::GLuint, arg1Rep::GLuint, arg1Mod::GLuint, arg2::GLuint, arg2Rep::GLuint, arg2Mod::GLuint, arg3::GLuint, arg3Rep::GLuint, arg3Mod::GLuint)::Void
export glColorFragmentOp3ATI
@getCFun "libGL" glAlphaFragmentOp1ATI glAlphaFragmentOp1ATI(op::GLenum, dst::GLuint, dstMod::GLuint, arg1::GLuint, arg1Rep::GLuint, arg1Mod::GLuint)::Void
export glAlphaFragmentOp1ATI
@getCFun "libGL" glAlphaFragmentOp2ATI glAlphaFragmentOp2ATI(op::GLenum, dst::GLuint, dstMod::GLuint, arg1::GLuint, arg1Rep::GLuint, arg1Mod::GLuint, arg2::GLuint, arg2Rep::GLuint, arg2Mod::GLuint)::Void
export glAlphaFragmentOp2ATI
@getCFun "libGL" glAlphaFragmentOp3ATI glAlphaFragmentOp3ATI(op::GLenum, dst::GLuint, dstMod::GLuint, arg1::GLuint, arg1Rep::GLuint, arg1Mod::GLuint, arg2::GLuint, arg2Rep::GLuint, arg2Mod::GLuint, arg3::GLuint, arg3Rep::GLuint, arg3Mod::GLuint)::Void
export glAlphaFragmentOp3ATI
@getCFun "libGL" glSetFragmentShaderConstantATI glSetFragmentShaderConstantATI(dst::GLuint, value::Ptr{GLfloat})::Void
export glSetFragmentShaderConstantATI
@getCFun "libGL" glMapObjectBufferATI glMapObjectBufferATI(buffer::GLuint)::Ptr{Void}
export glMapObjectBufferATI
@getCFun "libGL" glUnmapObjectBufferATI glUnmapObjectBufferATI(buffer::GLuint)::Void
export glUnmapObjectBufferATI
@getCFun "libGL" glPNTrianglesiATI glPNTrianglesiATI(pname::GLenum, param::GLint)::Void
export glPNTrianglesiATI
@getCFun "libGL" glPNTrianglesfATI glPNTrianglesfATI(pname::GLenum, param::GLfloat)::Void
export glPNTrianglesfATI
@getCFun "libGL" glStencilOpSeparateATI glStencilOpSeparateATI(face::GLenum, sfail::GLenum, dpfail::GLenum, dppass::GLenum)::Void
export glStencilOpSeparateATI
@getCFun "libGL" glStencilFuncSeparateATI glStencilFuncSeparateATI(frontfunc::GLenum, backfunc::GLenum, ref::GLint, mask::GLuint)::Void
export glStencilFuncSeparateATI
@getCFun "libGL" glNewObjectBufferATI glNewObjectBufferATI(size::GLsizei, pointer::Ptr{Void}, usage::GLenum)::Cuint
export glNewObjectBufferATI
@getCFun "libGL" glIsObjectBufferATI glIsObjectBufferATI(buffer::GLuint)::Bool
export glIsObjectBufferATI
@getCFun "libGL" glUpdateObjectBufferATI glUpdateObjectBufferATI(buffer::GLuint, offset::GLuint, size::GLsizei, pointer::Ptr{Void}, preserve::GLenum)::Void
export glUpdateObjectBufferATI
@getCFun "libGL" glGetObjectBufferfvATI glGetObjectBufferfvATI(buffer::GLuint, pname::GLenum, params::Ptr{GLfloat})::Void
export glGetObjectBufferfvATI
@getCFun "libGL" glGetObjectBufferivATI glGetObjectBufferivATI(buffer::GLuint, pname::GLenum, params::Ptr{GLint})::Void
export glGetObjectBufferivATI
@getCFun "libGL" glFreeObjectBufferATI glFreeObjectBufferATI(buffer::GLuint)::Void
export glFreeObjectBufferATI
@getCFun "libGL" glArrayObjectATI glArrayObjectATI(array::GLenum, size::GLint, type_::GLenum, stride::GLsizei, buffer::GLuint, offset::GLuint)::Void
export glArrayObjectATI
@getCFun "libGL" glGetArrayObjectfvATI glGetArrayObjectfvATI(array::GLenum, pname::GLenum, params::Ptr{GLfloat})::Void
export glGetArrayObjectfvATI
@getCFun "libGL" glGetArrayObjectivATI glGetArrayObjectivATI(array::GLenum, pname::GLenum, params::Ptr{GLint})::Void
export glGetArrayObjectivATI
@getCFun "libGL" glVariantArrayObjectATI glVariantArrayObjectATI(id::GLuint, type_::GLenum, stride::GLsizei, buffer::GLuint, offset::GLuint)::Void
export glVariantArrayObjectATI
@getCFun "libGL" glGetVariantArrayObjectfvATI glGetVariantArrayObjectfvATI(id::GLuint, pname::GLenum, params::Ptr{GLfloat})::Void
export glGetVariantArrayObjectfvATI
@getCFun "libGL" glGetVariantArrayObjectivATI glGetVariantArrayObjectivATI(id::GLuint, pname::GLenum, params::Ptr{GLint})::Void
export glGetVariantArrayObjectivATI
@getCFun "libGL" glVertexAttribArrayObjectATI glVertexAttribArrayObjectATI(index::GLuint, size::GLint, type_::GLenum, normalized::GLboolean, stride::GLsizei, buffer::GLuint, offset::GLuint)::Void
export glVertexAttribArrayObjectATI
@getCFun "libGL" glGetVertexAttribArrayObjectfvATI glGetVertexAttribArrayObjectfvATI(index::GLuint, pname::GLenum, params::Ptr{GLfloat})::Void
export glGetVertexAttribArrayObjectfvATI
@getCFun "libGL" glGetVertexAttribArrayObjectivATI glGetVertexAttribArrayObjectivATI(index::GLuint, pname::GLenum, params::Ptr{GLint})::Void
export glGetVertexAttribArrayObjectivATI
@getCFun "libGL" glVertexStream1sATI glVertexStream1sATI(stream::GLenum, x::GLshort)::Void
export glVertexStream1sATI
@getCFun "libGL" glVertexStream1svATI glVertexStream1svATI(stream::GLenum, coords::Ptr{GLshort})::Void
export glVertexStream1svATI
@getCFun "libGL" glVertexStream1iATI glVertexStream1iATI(stream::GLenum, x::GLint)::Void
export glVertexStream1iATI
@getCFun "libGL" glVertexStream1ivATI glVertexStream1ivATI(stream::GLenum, coords::Ptr{GLint})::Void
export glVertexStream1ivATI
@getCFun "libGL" glVertexStream1fATI glVertexStream1fATI(stream::GLenum, x::GLfloat)::Void
export glVertexStream1fATI
@getCFun "libGL" glVertexStream1fvATI glVertexStream1fvATI(stream::GLenum, coords::Ptr{GLfloat})::Void
export glVertexStream1fvATI
@getCFun "libGL" glVertexStream1dATI glVertexStream1dATI(stream::GLenum, x::GLdouble)::Void
export glVertexStream1dATI
@getCFun "libGL" glVertexStream1dvATI glVertexStream1dvATI(stream::GLenum, coords::Ptr{GLdouble})::Void
export glVertexStream1dvATI
@getCFun "libGL" glVertexStream2sATI glVertexStream2sATI(stream::GLenum, x::GLshort, y::GLshort)::Void
export glVertexStream2sATI
@getCFun "libGL" glVertexStream2svATI glVertexStream2svATI(stream::GLenum, coords::Ptr{GLshort})::Void
export glVertexStream2svATI
@getCFun "libGL" glVertexStream2iATI glVertexStream2iATI(stream::GLenum, x::GLint, y::GLint)::Void
export glVertexStream2iATI
@getCFun "libGL" glVertexStream2ivATI glVertexStream2ivATI(stream::GLenum, coords::Ptr{GLint})::Void
export glVertexStream2ivATI
@getCFun "libGL" glVertexStream2fATI glVertexStream2fATI(stream::GLenum, x::GLfloat, y::GLfloat)::Void
export glVertexStream2fATI
@getCFun "libGL" glVertexStream2fvATI glVertexStream2fvATI(stream::GLenum, coords::Ptr{GLfloat})::Void
export glVertexStream2fvATI
@getCFun "libGL" glVertexStream2dATI glVertexStream2dATI(stream::GLenum, x::GLdouble, y::GLdouble)::Void
export glVertexStream2dATI
@getCFun "libGL" glVertexStream2dvATI glVertexStream2dvATI(stream::GLenum, coords::Ptr{GLdouble})::Void
export glVertexStream2dvATI
@getCFun "libGL" glVertexStream3sATI glVertexStream3sATI(stream::GLenum, x::GLshort, y::GLshort, z::GLshort)::Void
export glVertexStream3sATI
@getCFun "libGL" glVertexStream3svATI glVertexStream3svATI(stream::GLenum, coords::Ptr{GLshort})::Void
export glVertexStream3svATI
@getCFun "libGL" glVertexStream3iATI glVertexStream3iATI(stream::GLenum, x::GLint, y::GLint, z::GLint)::Void
export glVertexStream3iATI
@getCFun "libGL" glVertexStream3ivATI glVertexStream3ivATI(stream::GLenum, coords::Ptr{GLint})::Void
export glVertexStream3ivATI
@getCFun "libGL" glVertexStream3fATI glVertexStream3fATI(stream::GLenum, x::GLfloat, y::GLfloat, z::GLfloat)::Void
export glVertexStream3fATI
@getCFun "libGL" glVertexStream3fvATI glVertexStream3fvATI(stream::GLenum, coords::Ptr{GLfloat})::Void
export glVertexStream3fvATI
@getCFun "libGL" glVertexStream3dATI glVertexStream3dATI(stream::GLenum, x::GLdouble, y::GLdouble, z::GLdouble)::Void
export glVertexStream3dATI
@getCFun "libGL" glVertexStream3dvATI glVertexStream3dvATI(stream::GLenum, coords::Ptr{GLdouble})::Void
export glVertexStream3dvATI
@getCFun "libGL" glVertexStream4sATI glVertexStream4sATI(stream::GLenum, x::GLshort, y::GLshort, z::GLshort, w::GLshort)::Void
export glVertexStream4sATI
@getCFun "libGL" glVertexStream4svATI glVertexStream4svATI(stream::GLenum, coords::Ptr{GLshort})::Void
export glVertexStream4svATI
@getCFun "libGL" glVertexStream4iATI glVertexStream4iATI(stream::GLenum, x::GLint, y::GLint, z::GLint, w::GLint)::Void
export glVertexStream4iATI
@getCFun "libGL" glVertexStream4ivATI glVertexStream4ivATI(stream::GLenum, coords::Ptr{GLint})::Void
export glVertexStream4ivATI
@getCFun "libGL" glVertexStream4fATI glVertexStream4fATI(stream::GLenum, x::GLfloat, y::GLfloat, z::GLfloat, w::GLfloat)::Void
export glVertexStream4fATI
@getCFun "libGL" glVertexStream4fvATI glVertexStream4fvATI(stream::GLenum, coords::Ptr{GLfloat})::Void
export glVertexStream4fvATI
@getCFun "libGL" glVertexStream4dATI glVertexStream4dATI(stream::GLenum, x::GLdouble, y::GLdouble, z::GLdouble, w::GLdouble)::Void
export glVertexStream4dATI
@getCFun "libGL" glVertexStream4dvATI glVertexStream4dvATI(stream::GLenum, coords::Ptr{GLdouble})::Void
export glVertexStream4dvATI
@getCFun "libGL" glNormalStream3bATI glNormalStream3bATI(stream::GLenum, nx::GLbyte, ny::GLbyte, nz::GLbyte)::Void
export glNormalStream3bATI
@getCFun "libGL" glNormalStream3bvATI glNormalStream3bvATI(stream::GLenum, coords::Ptr{GLbyte})::Void
export glNormalStream3bvATI
@getCFun "libGL" glNormalStream3sATI glNormalStream3sATI(stream::GLenum, nx::GLshort, ny::GLshort, nz::GLshort)::Void
export glNormalStream3sATI
@getCFun "libGL" glNormalStream3svATI glNormalStream3svATI(stream::GLenum, coords::Ptr{GLshort})::Void
export glNormalStream3svATI
@getCFun "libGL" glNormalStream3iATI glNormalStream3iATI(stream::GLenum, nx::GLint, ny::GLint, nz::GLint)::Void
export glNormalStream3iATI
@getCFun "libGL" glNormalStream3ivATI glNormalStream3ivATI(stream::GLenum, coords::Ptr{GLint})::Void
export glNormalStream3ivATI
@getCFun "libGL" glNormalStream3fATI glNormalStream3fATI(stream::GLenum, nx::GLfloat, ny::GLfloat, nz::GLfloat)::Void
export glNormalStream3fATI
@getCFun "libGL" glNormalStream3fvATI glNormalStream3fvATI(stream::GLenum, coords::Ptr{GLfloat})::Void
export glNormalStream3fvATI
@getCFun "libGL" glNormalStream3dATI glNormalStream3dATI(stream::GLenum, nx::GLdouble, ny::GLdouble, nz::GLdouble)::Void
export glNormalStream3dATI
@getCFun "libGL" glNormalStream3dvATI glNormalStream3dvATI(stream::GLenum, coords::Ptr{GLdouble})::Void
export glNormalStream3dvATI
@getCFun "libGL" glClientActiveVertexStreamATI glClientActiveVertexStreamATI(stream::GLenum)::Void
export glClientActiveVertexStreamATI
@getCFun "libGL" glVertexBlendEnviATI glVertexBlendEnviATI(pname::GLenum, param::GLint)::Void
export glVertexBlendEnviATI
@getCFun "libGL" glVertexBlendEnvfATI glVertexBlendEnvfATI(pname::GLenum, param::GLfloat)::Void
export glVertexBlendEnvfATI

end


text:
using AWSCore
using AWSS3
using AWSLambda

JL_VERSION_BASE="0.6"
JL_VERSION_PATCH="4"
JL_VERSION="$JL_VERSION_BASE.$JL_VERSION_PATCH"

image_name = "octech/$(replace(basename(pwd()), "_", "")):$JL_VERSION"

lambda_name = basename(pwd())

source_bucket = "octech.com.au.ap-southeast-2.awslambda.jl.deploy"

base_zip = "$(lambda_name)_$(VERSION)_$(AWSLambda.aws_lamabda_jl_version).zip"

if length(ARGS) == 0 || ARGS[1] == "build"
    cp("../../src/AWSLambda.jl", "AWSLambda.jl"; remove_destination=true)
    run(`docker build
            --build-arg JL_VERSION_BASE=$JL_VERSION_BASE
            --build-arg JL_VERSION_PATCH=$JL_VERSION_PATCH
             -t $image_name .`)
end

if length(ARGS) > 0 && ARGS[1] == "shell"
    run(`docker run --rm -it -v $(pwd()):/var/host $image_name bash`)
end

if length(ARGS) > 0 && ARGS[1] == "zip"
    rm(base_zip; force=true)
    cmd = `zip --symlinks -r -9 /var/host/$base_zip .`
    run(`docker run --rm -it -v $(pwd()):/var/host $image_name $cmd`)
end

if length(ARGS) > 0 && ARGS[1] == "deploy"
    AWSCore.Services.s3("PUT", "/$source_bucket/$base_zip",
                       headers=Dict("x-amz-acl" => "public-read"),
                       Body=read(base_zip))
end

if length(ARGS) > 0 && ARGS[1] == "deploy_regions"
    lambda_regions = ["us-east-1",
                      "us-east-2",
                      "us-west-1",
                      "us-west-2",
                      "ap-northeast-2",
                      "ap-south-1",
                      "ap-southeast-1",
                      "ap-northeast-1",
                      "eu-central-1",
                      "eu-west-1",
                      "eu-west-2"]

    @sync for r in lambda_regions

        raws = merge(default_aws_config(), Dict(:region => r))
        bucket = "octech.com.au.$r.awslambda.jl.deploy"

        @async begin

            s3_create_bucket(raws, bucket)

            AWSS3.s3(default_aws_config(), "PUT", bucket;
                     path = base_zip,
                     headers = Dict(
                         "x-amz-copy-source" => "$source_bucket/$base_zip",
                         "x-amz-acl" => "public-read"))
        end
    end
end


text:
using BetaVQE
using Yao, Yao.EasyBuild, Yao.BitBasis
using Test, Random, StatsBase
using Zygote
using BetaVQE.VAN

@testset "VAN" begin
    include("VAN/VAN.jl")
end

@testset "sample" begin
    Random.seed!(3)
    nbits = 4
    nhiddens = [10, 20]
    nsamples = 1000
    β = 1.0

    h = heisenberg(nbits)
    c = dispatch!(variational_circuit(nbits), :random)
    model = AutoRegressiveModel(nbits, nhiddens)

    configs = bitarray(collect(0:(1<<nbits-1)), nbits)
    logp = get_logp(model, configs)
    @test isapprox(sum(exp.(logp)), 1.0, rtol=1e-2)

    f = sum(exp.(logp) .* free_energy_local(β, h, model, c, configs))

    samples = gen_samples(model, nsamples)
    f_sample = free_energy(β, h, model, c, samples)
    @test isapprox(f, f_sample, rtol=1e-1)
end

@testset "circuit diff" begin
    Random.seed!(3)
    nbits = 4
    nhiddens = [10]
    nsamples = 2000

    h = heisenberg(nbits)
    c = dispatch!(variational_circuit(nbits), :random)
    model = AutoRegressiveModel(nbits, nhiddens)
    samples = gen_samples(model, nsamples)

    # check the gradient of circuit parameters
    params = parameters(c)
    ϵ = 1e-5
    for k in 1:length(params)
        params[k] -= ϵ
        dispatch!(c, params)
        l1 = free_energy(2.0, h, model, c, samples)
        params[k] += 2ϵ
        dispatch!(c, params)
        l2 = free_energy(2.0, h, model, c, samples)
        params[k] -= ϵ
        dispatch!(c, params)
        g = (l2-l1)/2ϵ
        g2 = gradient(c->free_energy(2.0, h, model, c, samples), c)[1]
        @test isapprox(g, g2[k], rtol=1e-2)
    end
end

@testset "network diff" begin
    Random.seed!(7)
    nbits = 2
    nhiddens = [10]
    nsamples = 1000

    h = heisenberg(nbits)
    c = dispatch!(variational_circuit(nbits), :random)
    model = AutoRegressiveModel(nbits, nhiddens)

    ϵ = 1e-5
    # check the gradient of model
    configs = bitarray(collect(0:(1<<nbits-1)), nbits)
    function f(model)
        logp = get_logp(model, configs)
        sum(exp.(logp) .* free_energy_local(2.0, h, model, c, configs))
    end

    m, n = 1, 1
    params = model_parameters(model)
    params[m][n] -= ϵ
    model_dispatch!(model, params)
    l1 = f(model)
    params[m][n] += 2ϵ
    model_dispatch!(model, params)
    l2 = f(model)
    g = (l2-l1)/2ϵ

    #tune it back
    params[m][n] -= ϵ
    model_dispatch!(model, params)
    samples = gen_samples(model, nsamples)
    g2 = gradient(model->free_energy(2.0, h, model, c, samples), model)[1]
    @test isapprox(g, g2.W[1][n], rtol=1e-1)
    @show g, g2.W[1][n]
end

@testset "tns circuit diff" begin
    Random.seed!(3)
    nx = 2
    ny = 2
    nbits = nx*ny
    depth = 3
    nhiddens = [10]
    nsamples = 2000

    h = heisenberg(nbits)
    c = tns_circuit(nbits, depth, EasyBuild.pair_square(nx, ny; periodic=false); entangler=(n,i,j)->put(n,(i,j)=>general_U4(rand(15)*2π)))
    model = AutoRegressiveModel(nbits, nhiddens)
    samples = gen_samples(model, nsamples)

    # check the gradient of circuit parameters
    params = parameters(c)
    ϵ = 1e-5
    for k in 1:length(params)
        params[k] -= ϵ
        dispatch!(c, params)
        l1 = free_energy(2.0, h, model, c, samples)
        params[k] += 2ϵ
        dispatch!(c, params)
        l2 = free_energy(2.0, h, model, c, samples)
        params[k] -= ϵ
        dispatch!(c, params)
        g = (l2-l1)/2ϵ
        g2 = gradient(c->free_energy(2.0, h, model, c, samples), c)[1]
        @test isapprox(g, g2[k], atol=1e-2)
    end
end

@testset "train" begin
    nbits = 4
    h = hamiltonian(TFIM(nbits, 1; Γ=0.0, periodic=false))
    network = PSAModel(nbits)
    circuit = tns_circuit(nbits, 2, EasyBuild.pair_square(nbits, 1; periodic=false); entangler=(n,i,j)->put(n,(i,j)=>general_U4()))
    network_params, circuit_params = train(1.0, h, network, circuit; nbatch=1000, niter=100)
    samples = gen_samples(network, 1000)
    @test free_energy(1.0, h, network, circuit, samples) <= -3.2
end


text:
# constants
const a = 6.378137e6 # the semi-major axis of the earth
const F = 298.257222101 # the inverse 1st flattening
const n = 1/(2F-1)
const m₀ = 0.9999 # central meridian scale factor
A₀ = 1 + n^2/4 + n^4/64
const A̅ = A₀/(1+n) *m₀*a

## longitude, latitude of origin for zones in Japan
const Origin_LatLon_Japan = [ 33.0  129.5    ; # 1
                              33.0  131.0    ; # 2
                              36.0  132.0+1/6; # 3
                              33.0  133.5    ; # 4
                              36.0  134.0+1/3; # 5
                              36.0  136.0    ; # 6
                              36.0  137.0+1/6; # 7
                              36.0  138.5    ; # 8
                              36.0  139.0+5/6; # 9
                              40.0  140.0+5/6; # 10
                              44.0  140.25   ; # 11
                              44.0  142.25   ; # 12
                              44.0  144.25   ; # 13
                              26.0  142.0    ; # 14
                              26.0  127.5    ; # 15
                              26.0  124.0    ; # 16
                              26.0  131.0    ; # 17
                              20.0  136.0    ; # 18
                              26.0  154.0    ; # 19
                             ]