text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _mphase{T<:AbstractFloat}(jd::T)
    ram, decm, dism = moonpos(jd, radians=true)
    ras, decs = sunpos(jd, radians=true)
    # phi: geocentric elongation of the Moon from the Sun
    # inc: selenocentric (Moon centered) elongation of the Earth from the Sun
    phi = acos(sin(decs)*sin(decm) + cos(decs)*cos(decm)*cos(ras - ram))
    # "dism" is in kilometers, AU in meters
    inc = atan2(AU*sin(phi), dism*1e3 - AU*cos(phi))
    return (1.0 + cos(inc))/2.0
end

"""
    mphase(jd) -> k

### Purpose ###

Return the illuminated fraction of the Moon at given Julian date(s).

### Arguments ###

* `jd`: the Julian ephemeris date.  It can be either a scalar or an array.

### Output ###

The illuminated fraction \$k\$ of Moon's disk, with \$0 \\leq k \\leq 1\$. \$k
= 0\$ indicates a new moon, while \$k = 1\$ stands for a full moon.  If `jd` is
given as an array, an array of the same number of elements as `jd` is returned.

### Method ###

Algorithm from Chapter 46 of "Astronomical Algorithms" by Jean Meeus
(Willmann-Bell, Richmond) 1991.  `sunpos` and `moonpos` are used to get
positions of the Sun and the Moon, and the Moon distance.  The selenocentric
elongation of the Earth from the Sun (phase angle) is then computed, and used to
determine the illuminated fraction.

### Example ###

Plot the illuminated fraction of the Moon for every day in January 2018 with a
hourly sampling.  Use [PyPlot.jl](https://github.com/stevengj/PyPlot.jl) for
plotting

``` julia
using PyPlot
points = DateTime(2018,01,01):Dates.Hour(1):DateTime(2018,01,31,23,59,59);
plot(points, mphase(jdcnv(points)))
```

Note that in this calendar month there are two full moons, this event is called
[blue moon](https://en.wikipedia.org/wiki/Blue_moon).

### Notes ###

Code of this function is based on IDL Astronomy User's Library.
"""
mphase(jd::Real) = _mphase(float(jd))

function mphase{J<:Real}(jd::AbstractArray{J})
    k = similar(jd, float(J))
    for i in eachindex(jd)
        k[i] = mphase(jd[i])
    end
    return k
end


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

function resubplot(infile::AbstractString, newcols::Int, newrows::Int, subplots::Vector{Int64})

  @assert length(subplots) <= newcols * newrows

  fileroot, fileext = splitext(infile)
  outfile = string(fileroot, "_resubplot", fileext)

  s = readall(infile)

  #find group plot command
  m = match(r"\\begin\{groupplot\}\[group style=\{.*group size=\s*(\d)\s*by\s*(\d)\s*\}\]", s)

  if m == nothing
    error("resubplot::Group Plot not found!")
  end

  oldcols, oldrows = Int(m.captures)
  coloffset, rowoffset = Int(m.offsets)

  part1 = string(s[1:(coloffset - 1)], newcols,
                 s[(coloffset + 1):rowoffset-1], newrows,
                 s[(rowoffset + 1): m.offset + length(m.match)]) #manual regex replace

  transitions = Int64[]
  idx = 1
  while true

    # try matching between groupplots
    m = match(r"(\\nextgroupplot|\\end\{groupplot\})", s, idx)

    if m == nothing
      break;
    end

    push!(transitions, m.offsets[1])
    idx = m.offsets[1] + length(m.captures[1])
  end

  groupplots =ASCIIString[]
  for (r1, r2) in zip(transitions[1:end - 1], transitions[2:end] - 1)
    push!(groupplots, s[r1:r2])
  end

  # only include the ones in subplots
  part2 = string(groupplots[subplots]...)

  part3 = s[transitions[end]:end]

  f = open(outfile, "w")
  write(f, part1, part2, part3)
  close(f)

  return outfile
end

function tex2embedded(infile::AbstractString; scaletext::Bool = true)
  #To use the embedded tikz tex, input it inside a figure
  #\begin{figure}
  #\centering
  #\resizebox{\columnwidth}{!}{%
  #                            \input{fig1.tikz}
  #                            }%
  #\caption{figure1}
  #\label{fig:figure1}
  #\end{figure}

  s = readall(infile)

  begin_matches = eachmatch(r"\\begin\{tikzpicture\}", s)
  end_matches = eachmatch(r"\\end\{tikzpicture\}\n", s)

  if begin_matches == nothing || end_matches == nothing
      error("tex2embedded:tikzpicture not found!")
    end

  #FIXME: this is pretty inefficient
  @test length(collect(begin_matches)) == length(collect(end_matches)) #begin and ends should correspond

  i = 1
  for (b, e) in zip(begin_matches, end_matches)

    b_ind = b.offset
    e_ind = e.offset + length(e.match)

    #FIXME: this is pretty inefficient. does s[b:e] create new substrings?
    m = match(r"(\\begin\{tikzpicture\}[\s\S]*\\end\{tikzpicture\}\n)", s[b_ind:e_ind])

    if m == nothing
      error("tex2embedded::not matched!")
    end

    tikz = m.captures[1]

    if scaletext
      tikz = replace(tikz, "\\begin{tikzpicture}[]", "\\begin{tikzpicture}[transform shape]")
    end

    outfile = string(splitext(infile)[1], "_$i.tikz")
    f = open(outfile, "w")
    write(f, tikz)
    close(f)

    i += 1
  end

end

ext: jl

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS


text:
include("euler/euler.jl")
using .BigIntegers: BigInteger

function compute()::BigInteger
    number = BigInteger(138_902_663)
    j = 3
    while(!all(parse(Int, (number * number)[i * 2 + 1]) == i + 1 for i ∈ 0:8))
        if j == 3
            number -= 6
            j = 7
        elseif j == 7
            number -= 4
            j = 3
        end
    end
    return number * 10
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
ITensors.itensor(A::ITensor) = A

# Insert missing diagonal blocks
function insert_diag_blocks!(T::Tensor)
  for b in eachdiagblock(T)
    blockT = blockview(T, b)
    if isnothing(blockT)
      # Block was not found in the list, insert it
      insertblock!(T, b)
    end
  end
end
insert_diag_blocks!(T::ITensor) = insert_diag_blocks!(tensor(T))

# Reshape into an order-2 ITensor
matricize(T::ITensor, inds::Index...) = matricize(T, inds)

function matricize(T::ITensor, inds)
  left_inds = commoninds(T, inds)
  right_inds = uniqueinds(T, inds)
  return matricize(T, left_inds, right_inds)
end

function matricize(T::ITensor, left_inds, right_inds)
  CL = combiner(left_inds; dir=ITensors.Out, tags="CL")
  CR = combiner(right_inds; dir=ITensors.In, tags="CR")
  M = (T * CL) * CR
  return M, CL, CR
end

function setdims(t::NTuple{N,Pair{QN,Int}}, dims::NTuple{N,Int}) where {N}
  return first.(t) .=> dims
end

# XXX: generalize this function
function _getindex(T::DenseTensor{ElT,N}, I1::Colon, I2::UnitRange{Int64}) where {ElT,N}
  A = array(T)[I1, I2]
  return tensor(Dense(vec(A)), setdims(inds(T), size(A)))
end

function getblock(i::Index, n::Integer)
  return ITensors.space(i)[n]
end

# Make `Pair{QN,Int}` act like a regular `dim`
NDTensors.dim(qnv::Pair{QN,Int}) = last(qnv)

Base.:*(qnv::Pair{QN,Int}, d::ITensors.Arrow) = qn(qnv) * d => dim(qnv)

function getblock_preserve_qns(T::Tensor, b::Block)
  # TODO: make `T[b]` preserve QNs
  Tb = T[b]
  indsTb = getblock.(inds(T), Tuple(b)) .* dir.(inds(T))
  return ITensors.setinds(Tb, indsTb)
end

function blocksparsetensor(blocks::Dict{B,TB}) where {B,TB}
  b1, Tb1 = first(pairs(blocks))
  N = length(b1)
  indstypes = typeof.(inds(Tb1))
  blocktype = eltype(Tb1)
  indsT = getindex.(indstypes)
  # Determine the indices from the blocks
  for (b, Tb) in pairs(blocks)
    indsTb = inds(Tb)
    for n in 1:N
      bn = b[n]
      indsTn = indsT[n]
      if bn > length(indsTn)
        resize!(indsTn, bn)
      end
      indsTn[bn] = indsTb[n]
    end
  end
  T = BlockSparseTensor(blocktype, indsT)
  for (b, Tb) in pairs(blocks)
    if !isempty(Tb)
      T[b] = Tb
    end
  end
  return T
end

function _nullspace_hermitian(M::Tensor; atol::Real=0.0)
  tol = atol
  # Insert any missing diagonal blocks
  insert_diag_blocks!(M)
  #D, U = eigen(Hermitian(M))
  Dᵢₜ, Uᵢₜ = eigen(itensor(M); ishermitian=true)
  D = tensor(Dᵢₜ)
  U = tensor(Uᵢₜ)
  nullspace_blocks = Dict()
  for bU in nzblocks(U)
    bM = Block(bU[1], bU[1])
    bD = Block(bU[2], bU[2])
    # Assume sorted from largest to smallest
    indstart = sum(d -> abs(d) .> tol, storage(D[bD])) + 1
    Ub = getblock_preserve_qns(U, bU)
    indstop = lastindex(Ub, 2)
    Nb = _getindex(Ub, :, indstart:indstop)
    nullspace_blocks[bU] = Nb
  end
  return blocksparsetensor(nullspace_blocks)
end

function LinearAlgebra.nullspace(M::Hermitian{<:Number,<:Tensor}; kwargs...)
  return _nullspace_hermitian(parent(M); kwargs...)
end

function LinearAlgebra.nullspace(::Order{2}, M::ITensor, left_inds, right_inds; kwargs...)
  @assert order(M) == 2
  M² = prime(dag(M), right_inds) * M
  M² = permute(M², right_inds'..., right_inds...)
  M²ₜ = tensor(M²)
  Nₜ = nullspace(Hermitian(M²ₜ); kwargs...)
  indsN = (Index(ind(Nₜ, 1); dir=ITensors.In), Index(ind(Nₜ, 2); dir=ITensors.In))
  N = dag(itensor(ITensors.setinds(Nₜ, indsN)))
  # Make the index match the input index
  Ñ = replaceinds(N, (ind(N, 1),) => right_inds)
  return Ñ
end

function LinearAlgebra.nullspace(T::ITensor, is...; kwargs...)
  M, CL, CR = matricize(T, is...)
  @assert order(M) == 2
  cL = commoninds(M, CL)
  cR = commoninds(M, CR)
  N₂ = nullspace(Order(2), M, cL, cR; kwargs...)
  return N₂ * CR
end


text:
function benchmark(dataset::D,truth::Array{Array{Int64,1},1},s::Sorting,cal_length=15.0, sample_rate=20000) where D

    (buf,nums)=output_buffer(1)

    #Calibrate

    cal_samples=round(Int,sample_rate*cal_length)
    v=zeros(eltype(dataset),div(cal_samples,2),1)

    #Threshold calibration for first half of calibration time
    v[:,1]=dataset[1:size(v,1),1]
    cal!(s,v,buf,nums,0)

    #Clustering / DM calibration for second half of calibration time
    v[:,1]=dataset[(size(v,1)+1):(2*size(v,1)),1]
    cal!(s,v,buf,nums,2)

    iter=div(sample_rate,20)
    spikes=[Array{Int64,1}() for i=0:0]
    v=zeros(eltype(dataset),iter,1)
    num_clusters=0

    for i=cal_samples+1:iter:(round(Int64,length(dataset)/sample_rate)*sample_rate-sample_rate)

        v[:,1]=dataset[i:(i+iter-1),1]
        onlinesort!(s,v,buf,nums)

        for k=1:nums[1]
            if buf[k,1].id>num_clusters
                num_clusters=buf[k,1].id
                push!(spikes,zeros(Int64,0))
            end
            push!(spikes[buf[k,1].id],buf[k,1].inds[1]+i+25)
            buf[k,1]=Spike()
        end
        nums[1]=0
    end

    #ISI violations

    #Accuracy due to overlap, clustering, and detection phases
    (spike_totals,TP,myFP,myFN)=accuracy_bench(dataset,spikes,truth,cal_length,sample_rate)

    #speed calculations

    (spikes,spike_totals,TP,myFP,myFN)
end

#=
function benchmark_all(dataset::Array{Int64,2},newstep::Algorithm)

    masterlist=Array(Any,5)

    steps=subtypes(Algorithm)

    for i=1:length(masterlist)
        if steps[i]!=super(newstep)
            masterlist[i]=subtypes(mystep)
        else
            masterlist[i]=newstep
        end
    end
end
=#

function accuracy_bench(dataset::D,spikes::Array{Array{Int64,1},1},truth::Array{Array{Int64,1},1},cal_length::Float64, sample_rate::Int64) where D

    #first go through detected spikes to classify as false positives or true positives
    #then go through ground truth to look for false negatives

    #need to figure out what cluster corresponds to what neuron
    corrmat=zeros(Float64,length(truth),length(spikes))

    cal_samples=round(Int,sample_rate*cal_length)
    data_s=zeros(eltype(dataset),length(dataset))
    data_t=zeros(eltype(dataset),length(dataset))

    for i=1:length(spikes)

        data_s[:]=0
        for k in spikes[i]
            if k+25<length(dataset)
                if k>cal_samples
                    data_s[k-25:k+25]=1
                end
            end
        end

        for j=1:length(truth)
            data_t[:]=0
            spike_num=0
            for k in truth[j]
                if k+25<length(dataset)
                    if k>cal_samples
                        data_t[k-25:k+25]=1
                        spike_num+=1
                    end
                end
            end

            test=data_s+data_t

            l=1
            while l<length(test)-1
                if test[l]>1
                    l+=1
                    forward=true
                    while forward
                        if test[l]>1
                            test[l]=0
                            l+=1
                        else
                            forward=false
                        end
                    end
                end
                l+=1
            end
            corrmat[j,i]=sum(test.>1)/spike_num*100.0
        end
    end

    assignedd=zeros(Int64,length(spikes)) #zero indicates it doesn't correspond to anything
    assigneds=falses(length(truth))
    count=1

    overlaps=zeros(Int64,length(dataset))
    spike_totals=zeros(Int64,length(truth))

    for i=1:length(truth)
        for k in truth[i]
            if k+25<length(dataset)
                if k>cal_samples
                    overlaps[k-25:k+25]+=1
                    spike_totals[i]+=1
                end
            end
        end
    end

    mytotals=zeros(Int64,length(dataset))
    for i=1:length(truth)
        for k in truth[i]
            if k+25<length(dataset)
                if k>cal_samples
                    mytotals[k]+=2^i
                end
            end
        end
    end

    if length(truth)>=length(spikes) #fewer neurons (or equal) detected than actually exist
        while count <= length(spikes)
            myinds=ind2sub(size(corrmat),indmax(corrmat))

            assignedd[myinds[2]]=myinds[1]
            assigneds[myinds[1]]=true
            count+=1

            corrmat[myinds[1],:]=0.0
            corrmat[:,myinds[2]]=0.0
        end

    else #more neurons detected than actually exist
        while count <= length(truth)
            myinds=ind2sub(size(corrmat),indmax(corrmat))

            assignedd[myinds[2]]=myinds[1]
            assigneds[myinds[1]]=true
            count+=1

            corrmat[myinds[1],:]=0.0
            corrmat[:,myinds[2]]=0.0
        end
    end

    TP=0
    FP_C=0
    FP_O=0
    FP_N=0
    FN_T=0
    FN_O=0
    #1 yet undetected
    #-1 detected
    tol=10

    missed=[zeros(Int64,0) for i=1:length(truth)]

    for i=1:length(spikes)
        if assignedd[i]==0
            FP_C+=sum(spikes[i].>cal_samples)
        else
            data_t[:]=0
            for k in truth[assignedd[i]]
                if k+tol<length(dataset)
                    if k>cal_samples
                        data_t[k]=k
                    end
                end
            end
            for j=1:length(spikes[i])
                found=false
                for k=(spikes[i][j]-tol):(spikes[i][j]+tol)
                    if data_t[k]>0 #found
                        TP+=1
                        found=true
                        mytotals[k] -= 2^assignedd[i]
                        break
                    end
                end
                if found==false
                    ov=false
                    clus=false
                    fn=false
                    for k=(spikes[i][j]-tol):(spikes[i][j]+tol)
                        if overlaps[k]>1
                            ov=true
                        elseif mytotals[k]>0
                            mytotals[k]=0
                            clus=true
                        end
                    end
                    if ov==true
                        FP_O+=1
                    elseif clus==true
                        FP_C+=1
                    else
                        FP_N+=1
                    end
                end
            end
        end
    end

    #find missed spikes
    for i=1:length(truth)
        for j in truth[i]
            if j>cal_samples
                md=digits(mytotals[j],2)
                if length(md)>i
                    if md[i+1]==1
                        push!(missed[i],j)
                    end
                end
            end
        end
    end

    for i=1:length(missed)
        for j=1:length(missed[i])
            ov=false
            for k=(missed[i][j]-tol):(missed[i][j]+tol)
                if overlaps[k]>1
                    ov=true
                    break
                end
            end
            if ov==true
                FN_O+=1
            else
                FN_T+=1
            end
        end
    end

    println("Spike totals: ", spike_totals)
    println("False Positive Total: ", FP_C+FP_O+FP_N)
    println("False Positive Clustering: ", FP_C)
    println("False Positive Overlap: ", FP_O)
    println("False Positive Noise: ", FP_N)
    println("True Positive: ", TP)
    println("Total False Negative: ", FN_O+FN_T)
    println("False Negative Overlap: ", FN_O)
    println("False Negative Threshold: ", FN_T)

myFP=[FP_C,FP_O,FP_N]
myFN=[FN_O,FN_T]

    (spike_totals,TP,myFP,myFN)
end

#ISI violation calculation
#=
function ISI_violations()

end
=#

#Speed / neuron

#parallel vs not parallel


text: Genie.secret_token!("73059b46d5f07e8bbf19b78b4e4e05a38feea1737ae66cf363083a62c3d4a75b")


text:
mutable struct Engine
	renderer::GRU.Renderer
	defs::Dict{Symbol, Any}
	dataPath::String
	assets::Dict{Symbol, Any}
	events::EventHandlers
	shouldClose::Bool
	timePrev::Float64
	timeNow::Float64
	api::Symbol
	apiVersion::Tuple{Int, Int}
	window::GLFW.Window

	function Engine(dataPath::String)
		now = time()
		new(GRU.Renderer(), Dict{Symbol, Any}(), dataPath, Dict{Symbol, Any}(), Dict{Symbol, Vector{Function}}(), false, now, now, :none, (0, 0))
	end
end

init(engine::Engine, defname::String) = init(engine, json_load(asset_path(engine, defname*".json")))

function init(engine::Engine, def::Dict{Symbol, Any} = Dict{Symbol, Any}())
	init_api(engine, def)
	init_window(engine, def)
	GRU.init(engine.renderer)
	preamble = get(def, :shader_preamble) do; get_shader_preamble(engine) end
	GRU.set_shader_preamble(engine.renderer, preamble)
	init(engine.renderer, def)
	FTFont.init()
	GLHelper.gl_info()
end

function done(engine::Engine)
	while !isempty(engine.assets)
		id, res = first(engine.assets)
		remove_asset(engine, id)
	end

	FTFont.done()
	GRU.done(engine.renderer)
	done_window(engine)
end

function init_api(engine::Engine, def::Dict{Symbol, Any})
	api = (Sys.islinux() && Sys.ARCH == :arm) ? :gles : :gl
	engine.api = get(def, :api, api)
	local ver
	if api == :gl
		ver = (3, 3)
	elseif api == :gles
		ver = (3, 0)
	else
		ver = (0, 0)
	end
	engine.apiVersion = map(Int, get(def, :api_version, ver))
end

function get_shader_preamble(engine::Engine)
	preamble =
		"""
		#version $(engine.apiVersion[1])$(engine.apiVersion[2])0 $(engine.api == :gles ? "es" : "")
		precision mediump float;

		"""
	Vector{UInt8}(preamble)
end

should_close(engine::Engine) = engine.shouldClose || GLFW.WindowShouldClose(engine.window)

function render(engine::Engine)
	call_event(engine, :render)
	GRU.render_frame(engine.renderer)
end

deltatime(engine::Engine) = engine.timeNow - engine.timePrev

function update(engine::Engine)
	engine.timePrev = engine.timeNow
	engine.timeNow = time()
	call_event(engine, :update)
end

function input(engine::Engine)
	GLFW.PollEvents()
	call_event(engine, :input)
end

function run(engine::Engine)
	init_viewport(engine)

	while !should_close(engine)
		render(engine)
		update(engine)

		GLFW.SwapBuffers(engine.window)

		input(engine)

		yield()
	end
end

asset_path(engine::Engine, path::String) = joinpath(engine.dataPath, path)

asset_id(filename::String, args...) = filename * reduce((v1, v2)->"$(v1)_$v2", args; init = "")
id_count = let count = 0
	()->count += 1
end
next_id(prefix::String) = asset_id(prefix, id_count())

function add_asset(engine::Engine, id::Symbol, val)
	@assert !haskey(engine.assets, id)
	engine.assets[id] = val
end

function remove_asset(engine, id)
	delete!(engine.assets, id)
end

function set_typed(val, def::Dict{Symbol, Any}, key::Symbol, dataType::DataType)
	if !isa(val, dataType)
		if isa(val, String) && !(dataType <: AbstractString || dataType == Symbol)
			val = eval(Meta.parse(val))
		end
		if !isa(val, dataType)
			val = dataType(val)
		end
		def[Symbol(string(key)*"#org")] = def[key] # remember the original value in case we need to resave the json
		def[key] = val
	end
	val
end

get_typed(def::Dict{Symbol, Any}, key::Symbol, dataType::DataType) = set_typed(def[key], def, key, dataType)
get_typed!(default::Function, def::Dict{Symbol, Any}, key::Symbol, dataType::DataType) = set_typed(get!(default, def, key), def, key, dataType)
get_typed!(def::Dict{Symbol, Any}, key::Symbol, defVal, dataType::DataType = typeof(defVal)) = set_typed(get!(()->defVal, def, key), def, key, dataType)

function get_transform(def::Dict{Symbol, Any})
	mat = Matrix{Float32}(I, 4, 4)
	if haskey(def, :transform)
		mat[:] = get_typed(def, :transform, Vector{Float32})
	else
		if haskey(def, :scale)
			mat = Math3D.scale(get_typed(def, :scale, Vector{Float32}))
		end
		if haskey(def, :rot_axis_angle)
			axis_angle = get_typed(def, :rot_axis_angle, Vector{Float32})
			mat = Math3D.rot(axis_angle[1:3], axis_angle[4]) * mat
		end
		if haskey(def, :position)
			mat = Math3D.trans(get_typed(def, :position, Vector{Float32})) * mat
		end
	end
	mat
end

get_id!(def::Dict{Symbol, Any}, prefix::String = string(def[:type])) = get_typed!(def, :id, Symbol) do; next_id(prefix) end

function resolve_def(engine::Engine, defpath::String, def::Dict{Symbol, Any})
	def[:defpath] = defpath
end

function merge_rec!(dst::AbstractDict, src::AbstractDict)
	for (k, v) in src
		if !haskey(dst, k)
			dst[k] = v
		elseif isa(dst[k], Associative) && isa(v, Associative)
			merge_rec!(dst[k], v)
		else
			@assert dst[k] == v
		end
	end
end

function def_container(engine::Engine, path::Vector{Symbol})
	container = engine.defs
	for i = 1:length(path)
		container = get!(container, path[i]) do; Dict{Symbol, Any}() end
	end
	container
end

function try_load(engine::Engine, path::Vector{Symbol})
	defname = join(map(string, path), '/')
	filename = asset_path(engine, defname*".json")
	if !isfile(filename)
		return nothing
	end
	json = json_load(filename)
	resolve_def(engine, defname, json)
	container = def_container(engine, path)
	merge_rec!(container, json)
	container
end

json_load(path::String) = JSON.parsefile(path, dicttype=Dict{Symbol, Any})

function get_def(engine::Engine, defname::String)
	path = map(Symbol, split(defname, '/'))
	container = def_container(engine, path[1:end-1])
	if haskey(container, path[end])
		return container[path[end]]
	end
	# if the def is not present, we start trying to load defs from prefixes of its path, and check if a loaded prefix contains the required key
	for i = length(path):-1:1
		loaded = try_load(engine, path[1:i])
		if loaded != nothing # a def was loaded, now we check if the required key was merged by this def into the container
			if haskey(container, path[end])
				return container[path[end]]
			end
		end
	end
	nothing
end

load_def(engine::Engine, defname::String) = load_def(engine, get_def(engine, defname))

function load_def(engine::Engine, def::Dict{Symbol, Any})
	instance = get(def, :instance, false)
	!isa(instance, Bool) && return instance
	objType = get_typed(def, :type, DataType)
	obj = init(engine, objType, def)
	if instance == true
		def[:instance] = obj
	end
	obj
end

function original_def(def::Dict{Symbol, Any})
	dst = Dict{Symbol, Any}()
	for (k,v) in def
		val = v
		if k == :instance
			val = isa(val, Bool) ? v : true
		elseif k == :defpath || endswith(string(k), "#org")
			continue
		else
			orgKey = Symbol(string(k)*"#org")
			if haskey(def, orgKey)
				val = def[orgKey]
			end
		end
		if isa(val, Dict{Symbol, Any})
			val = original_def(val)
		end
		dst[k] = val
	end
	dst
end

function save_def(engine::Engine, def::Dict{Symbol, Any}, target::String = def[:defpath])
	write(asset_path(engine, target*".json"), JSON.json(original_def(def), 2))
end

function init(engine::Engine, ::Type{T}, def::Dict{Symbol, Any}) where T
	obj = T()
	for i = 1:nfields(T)
		field = fieldname(T, i)
		if haskey(def, field)
			val = def[field]
			if isa(val, String)
				reference = get_def(engine, val)
				if isa(reference, Dict{Symbol, Any}) && haskey(val, :type) && get_typed(val, :type, DataType) <: fieldtype(T, i)
					val = reference
					def[field] = reference
				end
			end
			dstType = fieldtype(T, i)
			if isa(val, Dict{Symbol, Any}) && haskey(val, :type) && vget_typed(val, :type, DataType) <: dstType
				val = load_def(engine, val)
			end
			if dstType <: AbstractArray
				copy!(getfield(obj, field), val)
			else
				setfield!(obj, field, val)
			end
		end
	end
	obj
end


text:
using Random, Distributions

function normal_approximation_to_binomial(n, p)
    mu = n * p
    sigma = sqrt(n * (1-p) * p)
    return mu, sigma
end

function normal_cdf(lo, mu=0, sigma=1)
   return cdf(Normal(mu, sigma), lo)
end

# function normal_probability_below(lo, mu=0, sigma=1)
#   return normal_cdf(lo, mu, sigma)
# end

function normal_probability_above(lo, mu=0, sigma=1)
   return 1 - normal_cdf(lo, mu, sigma)
end

function normal_probability_between(lo, hi, mu=0, sigma=1)
   return normal_cdf(hi,mu, sigma) - normal_cdf(lo, mu, sigma)
end

function normal_probability_outside(lo, hi, mu=0, sigma=1)
    return 1 - normal_probability_between(lo, hi, mu, sigma)
end

@assert normal_cdf(0, 0 ,1)==0.5

function normal_upper_bound(probability, mu=0, sigma=1)
   return quantile(Normal(mu, sigma), probability)
end

function normal_lower_bound(probability, mu=0, sigma=1)
   return quantile(Normal(mu, sigma), 1 - probability)
end

function normal_two_sided_bounds(probability, mu=0, sigma=1)
   tail_probability = (1 - probability)/2

   upper_bound = normal_lower_bound(tail_probability, mu, sigma)

   lower_bound = normal_upper_bound(tail_probability, mu, sigma)

   return  lower_bound, upper_bound
end

mu_0, sigma_0 = normal_approximation_to_binomial(1000, 0.5)

@assert mu_0 == 500
@assert 15.8 < sigma_0 < 15.9
# 95% bounds based on assumption p is 0.5
lo, hi = normal_two_sided_bounds(0.95, mu_0, sigma_0)

# actual mu and sigma based on p = 0.55
mu_1, sigma_1 = normal_approximation_to_binomial(1000, 0.55)

# a type 2 error means we fail to reject the null hypothesis
# which will happen when X is still in our original interval
type_2_probability = normal_probability_between(lo, hi, mu_1, sigma_1)
power = 1 - type_2_probability      # 0.887
@assert 0.886 < power < 0.888

function two_sided_p_value(x, mu = 0, sigma = 1)
    """
    How likely are we to see a value at least as extreme as x (in either
    direction) if our values are from a N(mu, sigma)?
    """
   if (x >= mu)
        # x is greater than the mean, so the tail is everything greater than x
        return 2 * normal_probability_above(x, mu, sigma)
   else
        # x is less than the mean, so the tail is everything less than x
        return 2 * normal_probability_below(x, mu, sigma)
   end
end

two_sided_p_value(529.5, mu_0, sigma_0)   # 0.062

# TODO: Add julia code Monte-Carlo simulation

two_sided_p_value(531.5, mu_0, sigma_0)   # 0.0463

# TODO: Two sided tests

p_hat = 525 / 1000
mu = p_hat
sigma = sqrt(p_hat * (1 - p_hat) / 1000)   # 0.0158

normal_two_sided_bounds(0.95, mu, sigma)        # [0.4940, 0.5560]

p_hat = 540 / 1000
mu = p_hat
sigma = sqrt(p_hat * (1 - p_hat) / 1000) # 0.0158
normal_two_sided_bounds(0.95, mu, sigma) # [0.5091, 0.5709]


text:
using Distributions
using Turing
using Turing: gradient_forward, invlink, link, realpart
using ForwardDiff
using ForwardDiff: Dual
using Test

# Define model
@model ad_test_3() = begin
  v ~ Wishart(7, [1 0.5; 0.5 1])
  v
end
# Turing.TURING[:modelex]
# Call Turing's AD
# The result out is the gradient information on R
ad_test_3_f = ad_test_3()
vi = ad_test_3_f()
vvn = collect(Iterators.filter(vn -> vn.sym == :v, keys(vi)))[1]
_v = map(d -> realpart(d), vi[vvn])
_, grad_Turing = gradient_forward(vi[nothing], vi, ad_test_3_f)

dist_v = Wishart(7, [1 0.5; 0.5 1])

# Hand-written logp
function logp3(x::Vector)
  v = [x[1] x[3]; x[2] x[4]]
  lp = Turing.logpdf_with_trans(dist_v, v, false)
  lp
end

# Call ForwardDiff's AD
g = x -> ForwardDiff.gradient(logp3, vec(x));
# _s = link(dist_v, _s)
grad_FWAD = -g(_v)

# Compare result
@test grad_Turing ≈ grad_FWAD atol=1e-9


text:
using Plots, Statistics

function sigm1(x) 1/(1+exp(-x)) end
function isigm1(x) -log((x-1)/(x)) end

x=range(0,1,length=100)
y=sigm1.(10x.-7) # vari apo 0 ws 100, b apo 0 ws -100
plot(x,y)
y=sigm1.(100x.-70) # vari apo 0 ws 100, b apo 0 ws -100
plot!(x,y)

y=sigm1.(100x.-30)
plot(x,y)
y=sigm1.(100x.-50)
plot!(x,y)
y=sigm1.(100x.-30)+sigm1.(100x.-50)
plot!(x,y)
y=sigm1.(100x.-30)-sigm1.(100x.-50)
plot!(x,y)

y=sigm1.(100x.-30)-2sigm1.(100x.-50)
plot(x,y)
y=sigm1.(100x.-30)-2sigm1.(100x.-50)+3sigm1.(100x.-70)
plot(x,y)
y=sigm1.(10x.-3)-2sigm1.(10x.-5)+3sigm1.(10x.-7)
plot!(x,y)

x1=rand(1000)
x2=rand(1000)
y=0sigm1.(100x1.-0)+sigm1.(100x2.-80)
scatter3d(x1,x2,y,zcolor=y,camera=(50,40))
y=sigm1.(100x1.-40)+sigm1.(100x2.-80)
scatter3d(x1,x2,y,zcolor=y,camera=(5,40))
y=3sigm1.(100x1.-40)-20sigm1.(100x2.-80)
scatter3d(x1,x2,y,zcolor=y,camera=(7,40))
y=3sigm1.(10x1.-4)-20sigm1.(10x2.-8)
scatter3d(x1,x2,y,zcolor=y,camera=(7,40))

# counter example, why regression cannot make smooth curves (+... higher order polynomial)
x=range(0,1,length=100)
y=(2x.-0.8).^2
plot(x,y)
y=(2x.-1.2).^2
plot!(x,y)
y=(2x.-0.8).^2-0.2(2x.-1.2).^2
plot!(x,y)

obs=1000
vars=2
x=rand(obs,vars)
y=sin.(exp.(x[:,1])) .+cos.(5x[:,2])
y=(y.-minimum(y))./(maximum(y)-minimum(y))
y.*=0.8;y.+=0.1
scatter3d(x[:,1],x[:,2],y)
pred=1.3sigm1.(10x.-8)
scatter!(x,pred)
pred=1.2sigm1.(10x.-9).+0.07
scatter!(x,pred)
pred=0.7sigm1.(10x.-8)+0.3sigm1.(10x.-9).+0.01
scatter!(x,pred)
mean(abs.(y-pred))

opti_ww=rand(7)
opti_ww=[-9.91842, -1.59284, 4.58076, 1.2246, 8.79918, 8.4979, 0.185864]
pred=opti_ww[1]sigm1.(opti_ww[2]x[:,1].-opti_ww[3])+opti_ww[4]sigm1.(opti_ww[5]x[:,2].-opti_ww[6]).+opti_ww[7]
mean(abs.(y-pred))
scatter(x,y);scatter!(x,pred)

# func_evals=Int64(10000)
# include("run_its_4.jl")
# lb=-10ones(7)
# ub=-lb
# opti_ww,lb,ub=run_its_4(7,x,y,obs,func_evals,lb,ub,vars,opti_ww)
# pred=opti_ww[1]sigm1.(opti_ww[2]x.-opti_ww[3])+opti_ww[4]sigm1.(opti_ww[5]x.-opti_ww[6]).+opti_ww[7]
# mean(abs.(y-pred))
# scatter!(x,pred)


text:
function MatrixDFT(n::Int, ::Type{T} = Float64) where T
    out = fill(one(Complex{T}), n, n)
    for j ∈ 2:n
        out[j, 2] = exp( -2π*im*(j-1) / n )
    end
    for i ∈ 3:n, j ∈ i:n
        out[j,i] = out[j,i-1] * out[j,2]
    end
    Symmetric(out, :L) #./= √n
end

### Broadly, the strategy used here is to try and move as much of the work to compile time as possible.
### In particular, calculating the values of cos and sin at different points on the unit circle.
### This is enabled by aggressive loop unrolling in the functions, so all terms are known.
### This approach probably does not scale well to larger problems.

ωr(N, i) = cospi(-2(i % N)/N)
ωi(N, i) = sinpi(-2(i % N)/N)
ω(N, i) = ωr(N, i) + im*ωi(N, i)

ωr( ::Val{N}, ::Val{0} ) where N = 1.0
@generated ωr( ::Val{N}, ::Val{i} ) where {N,i} = ωr(N, i)
ωi( ::Val{N}, ::Val{0} ) where N = 0.0
@generated ωi( ::Val{N}, ::Val{i} ) where {N,i} = ωi(N, i)
@generated ωr( ::Val{N}, ::Val{i} ) where {N,i} = cos(-2π*i/N)
@generated ωi( ::Val{N}, ::Val{i} ) where {N,i} = sin(-2π*i/N)
ω( ::Val{N}, ::Val{0} ) where N = 1.0 + 0.0im
@generated ω( ::Val{N}, ::Val{i} ) where {N,i} = ωr(Val{N}(), Val{i}()) + im*ωi(Val{N}(), Val{i}())

function dft_sum_expr(N, ::Type{T} = Float64) where T
    d = N ÷ 2
    out = :( (ffte[1] + ffto[1], ) )
    for i ∈ 2:N
        m = mod1(i,d)
        g = Complex{T}(ω(Val{N}(), Val{i-1}()))
        expr = :( ffte[$m] + $g * ffto[$m])
        push!(out.args, expr )
    end
    out
end

function dft_even_odd_split_expr(N)
    d = N ÷ 2
    e_arg = :( (x[1], ) )
    o_arg = :( (x[2], ) )
    for i ∈ 2:d
        push!(e_arg.args, :(x[$(2i-1)]))
        push!(o_arg.args, :(x[$(2i)]))
    end
    e_arg, o_arg
end

@inline dft(x::NTuple{1,T}) where T = x
@generated function dft(x::NTuple{n,T}) where {n,T}
    out = dft_sum_expr(n, T)
    e_arg, o_arg = dft_even_odd_split_expr(n)
    quote
        $(Expr(:meta, :inline))
        ffte = dft($e_arg)
        ffto = dft($o_arg)
        $out
    end
end

function recursion_level(d_array, recursion_levels)
    end_recursion_size = d_array[recursion_levels + 1] #Smallest component in which we can break this down.

    step = 1 << recursion_levels
    initial_eval = quote end
    for i ∈ 1:step
        base_step!(initial_eval.args, end_recursion_size, i, step, d_array[1])
    end
    initial_eval
end

function base_step(blength, offset, step, N, x = :x_, r = :r_)
    base = :( $(Symbol(r, offset)) = $(Symbol(x, offset)) + $(Symbol(x, offset + step)) * $(ω(N, step)))
    for i ∈ 2:blength-1
        push!(base.args[2].args, :( $(Symbol(x, offset + i*step)) * $(ω(N, i*step)) ))
    end
    base
end

function base_step!(initial, blength, offset, step, N, x = :x_, r = :r_)
    for j ∈ 1:blength
        push!(initial.args, base_step(blength, offset, step, N, x, r))
    end
    initial
end

function compose( level, step )

end

qe = quote
        x = a + b
        y = c * d
        z = x * y + 4
    end

function recursion_meta(n)
    recursion_levels = 0
    d, r = divrem(n, 2)
    d_array = [n,d]
    while r != 1
        recursion_levels += 1
        d, r = divrem(d, 2)
        push!(d_array, d)
        d == 1 && break
    end
    d_array, recursion_levels
end

@generated function dfte(x::NTuple{n,T}) where {n,T}

    quote
        @nextract $n x x

    end
end

@nexprs $N i -> begin
    @nexprs $P-i j -> begin
        i * j
    end
end

@generated function dft(x::NTuple{n,T}, ::Val{N} = Val{n}()) where {n,T,N}
    d, r = divrem(N, 2)

    quote
        $(Expr(:meta, :inline))
        @nexprs $d r -> begin

        end
    end
end
@inline function dft(x::NTuple{1,T}, ::Val{N} = Val{1}()) where {T,N} = x
@inline function dft(x::NTuple{2,T}, ::Val{N} = Val{2}()) where {T,N}
    ωr(Val{1}())
end

function dft_call_expr(N)
    d = N ÷ 2
    e_arg = :( dftva(x[1]) )
    o_arg = :( dftva(x[2]) )
    for i ∈ 2:d
        push!(e_arg.args, :(x[$(2i-1)]))
        push!(o_arg.args, :(x[$(2i)]))
    end
    e_arg, o_arg
end
@inline dftva(x::T) where T <: Number = x
@generated function dftva(x::Vararg{T,n}) where {n,T}
    out = dft_sum_expr(n, T)
    e_arg, o_arg = dft_call_expr(n)
    quote
        $(Expr(:meta, :inline))
        ffte = $e_arg
        ffto = $o_arg
        $out
    end
end


text:
function create_compute_context()
    ctx    = create_some_context()
    device = first(devices(ctx))
    queue  = CmdQueue(ctx)
    return (device, ctx, queue)
end

opencl_version(obj :: CLObject) = api.parse_version(obj[:version])
opencl_version(c :: Context)  = opencl_version(first(devices(c)))
opencl_version(q :: CmdQueue) = opencl_version(q[:device])

const _versionDict = Dict{Ptr{Nothing}, VersionNumber}()

_deletecached!(obj :: CLObject) = delete!(_versionDict, pointer(obj))

function check_version(obj :: CLObject, version :: VersionNumber)
    version <= get!(_versionDict, pointer(obj)) do
        opencl_version(obj)
    end
end

"""
Format string using dict-like variables, replacing all accurancies of
`%(key)` with `value`.

Example:
    s = "Hello, %(name)"
    format(s, name="Tom")  ==> "Hello, Tom"
"""
function format(s::String; vars...)
    for (k, v) in vars
        s = replace(s, "%($k)" => v)
    end
    s
end

function build_kernel(ctx::Context, program::String,
                      kernel_name::String; vars...)
    src = format(program; vars...)
    p = Program(ctx, source=src)
    build!(p)
    return Kernel(p, kernel_name)
end

# cache for kernels; dict of form `(program_file, kernel_name, vars) -> kernel`
const CACHED_KERNELS = Dict{Tuple{String, String, Dict}, Kernel}()

function get_kernel(ctx::Context, program_file::String,
                    kernel_name::String; vars...)
    key = (program_file, kernel_name, Dict(vars))
    if in(key, keys(CACHED_KERNELS))
        return CACHED_KERNELS[key]
    else
        kernel = build_kernel(ctx, Base.read(program_file, String), kernel_name; vars...)
        CACHED_KERNELS[key] = kernel
        return kernel
    end
end

min_v11(elem) = check_version(elem, v"1.1")
min_v12(elem) = check_version(elem, v"1.2")
min_v20(elem) = check_version(elem, v"2.0")
min_v21(elem) = check_version(elem, v"2.1")
min_v22(elem) = check_version(elem, v"2.2")


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct StaticReserve{T <: ReserveDirection} <: Reserve{T}
        name::String
        available::Bool
        time_frame::Float64
        requirement::Float64
        ext::Dict{String, Any}
        operation_cost::Union{Nothing, TwoPartCost}
        internal::InfrastructureSystemsInternal
    end

Data Structure for a proportional reserve product for system simulations.

# Arguments
- `name::String`
- `available::Bool`
- `time_frame::Float64`: the relative saturation time_frame, validation range: `(0, nothing)`, action if invalid: `error`
- `requirement::Float64`: the static value of required reserves in system p.u., validation range: `(0, nothing)`, action if invalid: `error`
- `ext::Dict{String, Any}`
- `operation_cost::Union{Nothing, TwoPartCost}`: Cost for providing reserves  [`TwoPartCost`](@ref)
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct StaticReserve{T <: ReserveDirection} <: Reserve{T}
    name::String
    available::Bool
    "the relative saturation time_frame"
    time_frame::Float64
    "the static value of required reserves in system p.u."
    requirement::Float64
    ext::Dict{String, Any}
    "Cost for providing reserves  [`TwoPartCost`](@ref)"
    operation_cost::Union{Nothing, TwoPartCost}
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function StaticReserve{T}(name, available, time_frame, requirement, ext=Dict{String, Any}(), operation_cost=nothing, ) where T <: ReserveDirection
    StaticReserve{T}(name, available, time_frame, requirement, ext, operation_cost, InfrastructureSystemsInternal(), )
end

function StaticReserve{T}(; name, available, time_frame, requirement, ext=Dict{String, Any}(), operation_cost=nothing, internal=InfrastructureSystemsInternal(), ) where T <: ReserveDirection
    StaticReserve{T}(name, available, time_frame, requirement, ext, operation_cost, internal, )
end

# Constructor for demo purposes; non-functional.
function StaticReserve{T}(::Nothing) where T <: ReserveDirection
    StaticReserve{T}(;
        name="init",
        available=false,
        time_frame=0.0,
        requirement=0.0,
        ext=Dict{String, Any}(),
        operation_cost=TwoPartCost(nothing),
    )
end

InfrastructureSystems.get_name(value::StaticReserve) = value.name
"""Get [`StaticReserve`](@ref) `available`."""
get_available(value::StaticReserve) = value.available
"""Get [`StaticReserve`](@ref) `time_frame`."""
get_time_frame(value::StaticReserve) = value.time_frame
"""Get [`StaticReserve`](@ref) `requirement`."""
get_requirement(value::StaticReserve) = value.requirement
"""Get [`StaticReserve`](@ref) `ext`."""
get_ext(value::StaticReserve) = value.ext
"""Get [`StaticReserve`](@ref) `operation_cost`."""
get_operation_cost(value::StaticReserve) = value.operation_cost
"""Get [`StaticReserve`](@ref) `internal`."""
get_internal(value::StaticReserve) = value.internal

InfrastructureSystems.set_name!(value::StaticReserve, val) = value.name = val
"""Set [`StaticReserve`](@ref) `available`."""
set_available!(value::StaticReserve, val) = value.available = val
"""Set [`StaticReserve`](@ref) `time_frame`."""
set_time_frame!(value::StaticReserve, val) = value.time_frame = val
"""Set [`StaticReserve`](@ref) `requirement`."""
set_requirement!(value::StaticReserve, val) = value.requirement = val
"""Set [`StaticReserve`](@ref) `ext`."""
set_ext!(value::StaticReserve, val) = value.ext = val
"""Set [`StaticReserve`](@ref) `operation_cost`."""
set_operation_cost!(value::StaticReserve, val) = value.operation_cost = val


text:
"""
Weight Stabilizing Parallel Tempering

See: https://link.springer.com/article/10.1007/s11222-019-09863-3
"""
module WSPT

function gentempers(maxtemp, ntemps; degree=1)
  i = collect(1:ntemps)
  powers = (i.^ degree .- 1) / (ntemps ^ degree - 1)
  return maxtemp .^ powers
end

function compute_log_accept_ratio(loglike::Function, states, temps; verbose=0)
  s1, s2 = states
  t1, t2 = temps

  ll_s1_t2 = loglike(s1, t2)
  ll_s2_t1 = loglike(s2, t1)
  ll_s1_t1 = loglike(s1, t1)
  ll_s2_t2 = loglike(s2, t2)

  if verbose > 0
    println("ll_s1_t2: $(ll_s1_t2)")
    println("ll_s2_t1: $(ll_s2_t1)")
    println("ll_s1_t1: $(ll_s1_t1)")
    println("ll_s2_t2: $(ll_s2_t2)")
  end

  loglike_cand = ll_s1_t2 + ll_s2_t1
  loglike_curr = ll_s1_t1 + ll_s2_t2

  return loglike_cand - loglike_curr
end

function compute_accept_ratio(loglike::Function, states, temps; verbose=0)
  return exp(compute_log_accept_ratio(loglike, states, temps, verbose=verbose))
end

end


text:
function Shape(g::GeometricObject)
    Shape(g, nothing, nothing, Attributes(), Clip(nothing), Opacity(1), Operator())
end

function Shape(::Type{T}, args...) where T <: GeometricObject
    Shape(T(args...), nothing, nothing, Attributes(), Clip(nothing), Opacity(1), Operator())
end

function Shape(f::Function, ::Type{T}, deps...) where T
    Shape{T}((f, [deps...]), nothing, nothing, Attributes(), Clip(nothing), Opacity(1), Operator())
end

function Shapes(arr::Array{T}) where T <: GeometricObject
    Shapes(arr, nothing, nothing, Attributes(), Clip(nothing), Opacity(1), Operator())
end

function Shapes(::Type{T}, args...) where T <: GeometricObject
    Shapes(T.(args...), nothing, nothing, Attributes(), Clip(nothing), Opacity(1), Operator())
end

function Shapes(f::Function, ::Type{T}, deps...) where T
    Shapes{T}((f, [deps...]), nothing, nothing, Attributes(), Clip(nothing), Opacity(1), Operator())
end

const Shapelike = Union{Shape, Shapes}

Base.Broadcast.broadcastable(s::Shapelike) = Ref(s)
Base.Broadcast.broadcastable(s::Shapes) = Ref(s)

function Base.copy!(l::Layer, s::Shape{T}) where T
    news = Shape{T}(s.content, nothing, nothing, Attributes(), s.clip, s.opacity, s.operator)
    for (Ta, attr) in s.attrs.attrs
        news + attr
    end
    push!(l, news)
    news
end

function upward_transform(s::Shapelike)
    return upward_transform(s.parent)
end

function solve!(s::Shape{T}, cc::Union{Nothing, Cairo.CairoContext}=nothing) where T

    # the shape was already solved
    if !isnothing(s.solved)
        return s.solved
    end

    # the shape's content was given directly without a closure function
    if typeof(s.content) <: GeometricObject
        s.solved = s.content

        # extra logic for text, get the extents
        if T <: Txt
            t = s.solved
            if isnothing(t.extent)
                t_with_extent = Txt(t, TextExtent(cc, t))
                s.solved = t_with_extent
            end
        end

        return s.solved

    # the shape's content was given as a function with dependencies
    else
        closure = s.content[1]
        dependencies = s.content[2]

        solved_deps = []
        for d in dependencies
            if typeof(d) <: Shape
                solved = solve!(d, cc)

                t = transform_from_to(d, s, cc)
                transformed = t * solved
                push!(solved_deps, transformed)
            else
                # for other arguments that are not shapes
                push!(solved_deps, d)
            end
        end

        function transform_function(lay, x)
            transform_from_to(lay, s, cc) * x
        end
        function transform_function(x)
            inverse(upward_transform(s.parent, cc)) * x
        end

        return_values = if applicable(closure, solved_deps...)
            closure(solved_deps...)
        elseif applicable(closure, solved_deps..., transform_function)
            closure(solved_deps..., transform_function)
        else
            error("Arguments don't fit closure")
        end

        if typeof(return_values) <: Tuple
            s.solved = return_values[1]
            for a in return_values[2:end]
                if !(typeof(a) <: Attribute)
                    error("You can't return a non-attribute with type $(typeof(a)) in a tuple from a closure.")
                end
                s.attrs[typeof(a)] = a
            end
        else
            s.solved = return_values
        end

        if typeof(s.solved) <: Txt
            t = s.solved
            if isnothing(t.extent)
                t_with_extent = Txt(t, TextExtent(cc, t))
                s.solved = t_with_extent
            end
        end

        return s.solved
    end
    # the dependencies should be converted into this shapes's transform
    # then the solution is also in this shape's transform
end

function solve!(s::Shapes{T}, cc::Union{Nothing, Cairo.CairoContext}=nothing) where T <: GeometricObject

    # the shape was already solved
    if !isnothing(s.solved)
        return s.solved
    end

    # the shape's content was given directly without a closure function
    if typeof(s.content) <: Array{T}
        s.solved = s.content

        # extra logic for text, get the extents
        if T <: Txt
            ts = s.solved
            for (i, t) in enumerate(ts)
                if isnothing(t.extent)
                    t_with_extent = Txt(t, TextExtent(cc, ts))
                    ts[i] = t_with_extent
                end
            end
        end

        return s.solved

    # the shape's content was given as a function with dependencies
    else
        closure = s.content[1]
        dependencies = s.content[2]

        solved_deps = []
        for d in dependencies
            if typeof(d) <: Shape
                solved = solve!(d, cc)

                t = transform_from_to(d, s, cc)
                transformed = t * solved
                push!(solved_deps, transformed)
            elseif typeof(d) <: Shapes
                solved = solve!(d, cc)

                t = transform_from_to(d, s, cc)
                transformed = t .* solved
                push!(solved_deps, transformed)
            else
                # for other arguments that are not shapes
                push!(solved_deps, d)
            end
        end

        function transform_function(lay, x)
            transform_from_to(lay, s, cc) * x
        end
        function transform_function(x)
            inverse(upward_transform(s.parent, cc)) * x
        end

        return_values = if applicable(closure, solved_deps...)
            closure(solved_deps...)
        elseif applicable(closure, solved_deps..., transform_function)
            closure(solved_deps..., transform_function)
        else
            error("Arguments don't fit closure")
        end

        if typeof(return_values) <: Tuple
            s.solved = return_values[1]
            for a in return_values[2:end]
                if !(typeof(a) <: Attribute)
                    error("You can't return a non-attribute with type $(typeof(a)) in a tuple from a closure.")
                end
                s.attrs[typeof(a)] = a
            end
        else
            s.solved = return_values
        end

        if typeof(s.solved) <: Array{Txt}
            ts = s.solved
            for (i, t) in enumerate(ts)
                if isnothing(t.extent)
                    t_with_extent = Txt(t, TextExtent(cc, t))
                    ts[i] = t_with_extent
                end
            end
        end

        return s.solved
    end
    # the dependencies should be converted into this shapes's transform
    # then the solution is also in this shape's transform
end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

@given("put attribute type: is-alive, with value type: boolean") do context
    @fail "Implement me"
end

@given("put attribute type: age, with value type: long") do context
    @fail "Implement me"
end

@given("put attribute type: score, with value type: double") do context
    @fail "Implement me"
end

@given("put attribute type: birth-date, with value type: datetime") do context
    @fail "Implement me"
end

@given("put attribute type: name, with value type: string") do context
    @fail "Implement me"
end

@given("put attribute type: email, with value type: string") do context
    @fail "Implement me"
end

@given("attribute(email) as(string) set regex: \\S+@\\S+\\.\\S+") do context
    @fail "Implement me"
end
#=

from behave import *
from hamcrest import *

from typedb.api.concept.type.attribute_type import AttributeType
from tests.behaviour.config.parameters import parse_value_type, parse_list, parse_label
from tests.behaviour.context import Context

@step("put attribute type: {type_label}, with value type: {value_type}")
def step_impl(context: Context, type_label: str, value_type: str):
    context.tx().concepts().put_attribute_type(type_label, parse_value_type(value_type))

@step("attribute({type_label}) get value type: {value_type}")
def step_impl(context: Context, type_label: str, value_type: str):
    assert_that(context.tx().concepts().get_attribute_type(type_label).get_value_type(), is_(parse_value_type(value_type)))

@step("attribute({type_label}) get supertype value type: {value_type}")
def step_impl(context: Context, type_label: str, value_type: str):
    supertype = context.tx().concepts().get_attribute_type(type_label).as_remote(context.tx()).get_supertype()
    assert_that(supertype.get_value_type(), is_(parse_value_type(value_type)))

def attribute_type_as_value_type(context: Context, type_label: str, value_type: AttributeType.ValueType):
    attribute_type = context.tx().concepts().get_attribute_type(type_label)
    if value_type is AttributeType.ValueType.OBJECT:
        return attribute_type
    elif value_type is AttributeType.ValueType.BOOLEAN:
        return attribute_type.as_boolean()
    elif value_type is AttributeType.ValueType.LONG:
        return attribute_type.as_long()
    elif value_type is AttributeType.ValueType.DOUBLE:
        return attribute_type.as_double()
    elif value_type is AttributeType.ValueType.STRING:
        return attribute_type.as_string()
    elif value_type is AttributeType.ValueType.DATETIME:
        return attribute_type.as_datetime()
    else:
        raise ValueError("Unrecognised value type: " + str(value_type))

@step("attribute({type_label}) as({value_type}) get subtypes contain")
def step_impl(context: Context, type_label: str, value_type: str):
    sub_labels = [parse_label(s) for s in parse_list(context.table)]
    attribute_type = attribute_type_as_value_type(context, type_label, parse_value_type(value_type))
    actuals = list(map(lambda tt: tt.get_label(), attribute_type.as_remote(context.tx()).get_subtypes()))
    for sub_label in sub_labels:
        assert_that(sub_label, is_in(actuals))

@step("attribute({type_label}) as({value_type}) get subtypes do not contain")
def step_impl(context: Context, type_label: str, value_type: str):
    sub_labels = [parse_label(s) for s in parse_list(context.table)]
    attribute_type = attribute_type_as_value_type(context, type_label, parse_value_type(value_type))
    actuals = list(map(lambda tt: tt.get_label(), attribute_type.as_remote(context.tx()).get_subtypes()))
    for sub_label in sub_labels:
        assert_that(sub_label, not_(is_in(actuals)))

@step("attribute({type_label}) as({value_type}) set regex: {regex}")
def step_impl(context: Context, type_label: str, value_type, regex: str):
    value_type = parse_value_type(value_type)
    assert_that(value_type, is_(AttributeType.ValueType.STRING))
    attribute_type = attribute_type_as_value_type(context, type_label, value_type)
    attribute_type.as_remote(context.tx()).set_regex(regex)

@step("attribute({type_label}) as({value_type}) unset regex")
def step_impl(context: Context, type_label: str, value_type):
    value_type = parse_value_type(value_type)
    assert_that(value_type, is_(AttributeType.ValueType.STRING))
    attribute_type = attribute_type_as_value_type(context, type_label, value_type)
    attribute_type.as_remote(context.tx()).set_regex(None)

@step("attribute({type_label}) as({value_type}) get regex: {regex}")
def step_impl(context: Context, type_label: str, value_type, regex: str):
    value_type = parse_value_type(value_type)
    assert_that(value_type, is_(AttributeType.ValueType.STRING))
    attribute_type = attribute_type_as_value_type(context, type_label, value_type)
    assert_that(attribute_type.as_remote(context.tx()).get_regex(), is_(regex))

@step("attribute({type_label}) as({value_type}) does not have any regex")
def step_impl(context: Context, type_label: str, value_type):
    value_type = parse_value_type(value_type)
    assert_that(value_type, is_(AttributeType.ValueType.STRING))
    attribute_type = attribute_type_as_value_type(context, type_label, value_type)
    assert_that(attribute_type.as_remote(context.tx()).get_regex(), is_(None))

@step("attribute({type_label}) get key owners contain")
def step_impl(context: Context, type_label: str):
    owner_labels = [parse_label(s) for s in parse_list(context.table)]
    attribute_type = context.tx().concepts().get_attribute_type(type_label)
    actuals = list(map(lambda tt: tt.get_label(), attribute_type.as_remote(context.tx()).get_owners(only_key=True)))
    for owner_label in owner_labels:
        assert_that(actuals, has_item(owner_label))

@step("attribute({type_label}) get key owners do not contain")
def step_impl(context: Context, type_label: str):
    owner_labels = [parse_label(s) for s in parse_list(context.table)]
    attribute_type = context.tx().concepts().get_attribute_type(type_label)
    actuals = list(map(lambda tt: tt.get_label(), attribute_type.as_remote(context.tx()).get_owners(only_key=True)))
    for owner_label in owner_labels:
        assert_that(actuals, not_(has_item(owner_label)))

@step("attribute({type_label}) get attribute owners contain")
def step_impl(context: Context, type_label: str):
    owner_labels = [parse_label(s) for s in parse_list(context.table)]
    attribute_type = context.tx().concepts().get_attribute_type(type_label)
    actuals = list(map(lambda tt: tt.get_label(), attribute_type.as_remote(context.tx()).get_owners(only_key=False)))
    for owner_label in owner_labels:
        assert_that(actuals, has_item(owner_label))

@step("attribute({type_label}) get attribute owners do not contain")
def step_impl(context: Context, type_label: str):
    owner_labels = [parse_label(s) for s in parse_list(context.table)]
    attribute_type = context.tx().concepts().get_attribute_type(type_label)
    actuals = list(map(lambda tt: tt.get_label(), attribute_type.as_remote(context.tx()).get_owners(only_key=False)))
    for owner_label in owner_labels:
        assert_that(actuals, not_(has_item(owner_label)))

=#


text:
# ---
# title: 780. Reaching Points
# id: problem780
# author: Indigo
# date: 2021-01-28
# difficulty: Hard
# categories: Math
# link: <https://leetcode.com/problems/reaching-points/description/>
# hidden: true
# ---
#
# A move consists of taking a point `(x, y)` and transforming it to either `(x,
# x+y)` or `(x+y, y)`.
#
# Given a starting point `(sx, sy)` and a target point `(tx, ty)`, return `True`
# if and only if a sequence of moves exists to transform the point `(sx, sy)` to
# `(tx, ty)`. Otherwise, return `False`.
#
#
#
#     **Examples:**
#     Input: sx = 1, sy = 1, tx = 3, ty = 5
#     Output: True
#     Explanation:
#     One series of moves that transforms the starting point to the target is:
#     (1, 1) -> (1, 2)
#     (1, 2) -> (3, 2)
#     (3, 2) -> (3, 5)
#
#     Input: sx = 1, sy = 1, tx = 2, ty = 2
#     Output: False
#
#     Input: sx = 1, sy = 1, tx = 1, ty = 1
#     Output: True
#
#
#
# **Note:**
#
#   * `sx, sy, tx, ty` will all be integers in the range `[1, 10^9]`.
#
#
## @lc code=start
using LeetCode

function reaching_points(sx, sy, tx, ty)
    while tx >= sx && ty >= sy
        if tx == ty
            break
        end
        if tx > ty
            if ty > sy
                tx %= ty;
            else
                return (tx - sx) % ty == 0;
            end
        else
            if tx > sx
                ty %= tx;
            else
                return (ty - sy) % tx == 0;
            end
        end
    end
    tx == sx && ty == sy
end
## @lc code=end


text:
### A Pluto.jl notebook ###
# v0.16.4

using Markdown
using InteractiveUtils

# ╔═╡ 0fbc9f16-3669-11ec-155a-bbbd20bd8062
using Chess

# ╔═╡ 94c467ae-bb18-4386-af49-340bf5620d95
startboard()

# ╔═╡ 857ddcf5-e61f-4e40-a806-fee8f2e9d75f

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Chess = "717200cc-f167-4fd3-b4bf-b5e480529844"

[compat]
Chess = "~0.7.2"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[BinaryProvider]]
deps = ["Libdl", "Logging", "SHA"]
git-tree-sha1 = "ecdec412a9abc8db54c0efc5548c64dfce072058"
uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
version = "0.5.10"

[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "3533f5a691e60601fe60c90d8bc47a27aa2907ec"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.11.0"

[[Chess]]
deps = ["Artifacts", "Crayons", "Dates", "DefaultApplication", "Formatting", "Hiccup", "HypertextLiteral", "InteractiveUtils", "JSON", "Printf", "Random", "SQLite", "StaticArrays", "StatsBase", "UUIDs"]
git-tree-sha1 = "efad67292aa88bbf81cfb7331846e4899228da4e"
uuid = "717200cc-f167-4fd3-b4bf-b5e480529844"
version = "0.7.2"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "31d0151f5716b655421d9d75b7fa74cc4e744df2"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.39.0"

[[Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[DBInterface]]
git-tree-sha1 = "9b0dc525a052b9269ccc5f7f04d5b3639c65bca5"
uuid = "a10d1c49-ce27-4219-8d33-6db1a4562965"
version = "2.5.0"

[[DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DefaultApplication]]
deps = ["InteractiveUtils"]
git-tree-sha1 = "fc2b7122761b22c87fec8bf2ea4dc4563d9f8c24"
uuid = "3f0dd361-4fe0-5fc6-8523-80b14ec94d85"
version = "1.0.0"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[Hiccup]]
deps = ["MacroTools", "Test"]
git-tree-sha1 = "6187bb2d5fcbb2007c39e7ac53308b0d371124bd"
uuid = "9fb69e20-1954-56bb-a84f-559cc56a8ff7"
version = "0.2.2"

[[HypertextLiteral]]
git-tree-sha1 = "bc09b8d183505e128545c2a3b00cae5e5ea89c52"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.6.0"

[[InlineStrings]]
deps = ["Parsers"]
git-tree-sha1 = "19cb49649f8c41de7fea32d089d37de917b553da"
uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48"
version = "1.0.1"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "f0c6489b12d28fb4c2103073ec7452f3423bd308"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.1"

[[IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[LogExpFunctions]]
deps = ["ChainRulesCore", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "6193c3815f13ba1b78a51ce391db8be016ae9214"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.4"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "5a5bc6bf062f0f95e62d0fe0a2d99699fed82dd9"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.8"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "f19e978f81eca5fd7620650d7dbea58f825802ee"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.1.0"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[SQLite]]
deps = ["BinaryProvider", "DBInterface", "Dates", "Libdl", "Random", "SQLite_jll", "Serialization", "Tables", "Test", "WeakRefStrings"]
git-tree-sha1 = "ff81f9a47cbfaad52693ce0c74bbd86f2166f64d"
uuid = "0aa819cd-b072-5ff4-a722-6bc24af294d9"
version = "1.3.0"

[[SQLite_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "e601efdf2ce1df61147a1d1fed83cc428af1b53a"
uuid = "76ed43ae-9a5d-5a62-8c75-30186b810ce8"
version = "3.36.1+0"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "da4cf579416c81994afd6322365d00916c79b8ae"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "0.12.5"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[StatsAPI]]
git-tree-sha1 = "1958272568dc176a1d881acb797beb909c785510"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.0.0"

[[StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "eb35dcc66558b2dda84079b9a1be17557d32091a"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.12"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "fed34d0e71b91734bf0a7e10eb1bb05296ddbcd0"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.6.0"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[WeakRefStrings]]
deps = ["DataAPI", "InlineStrings", "Parsers"]
git-tree-sha1 = "c69f9da3ff2f4f02e811c3323c22e5dfcb584cfa"
uuid = "ea10d353-3f73-51f8-a26c-33c1cb351aa5"
version = "1.4.1"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╠═0fbc9f16-3669-11ec-155a-bbbd20bd8062
# ╠═94c467ae-bb18-4386-af49-340bf5620d95
# ╠═857ddcf5-e61f-4e40-a806-fee8f2e9d75f
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
export cg

function cg(A::Union(AbstractMatrix, Function),b::Vector;
	tol::Real=1e-6,maxIter::Int64=1000,
	M::Union(AbstractMatrix, Function)=x->x ,x0::Vector=[],out::Int64=0)

	n = length(b)
	nb = norm(b)

	# allocate memory once
	Ap = zeros(eltype(b),n)

	Af =  isa(A,Function) ? A : x->A_mul_B!(Ap, A, x)
	Mf =  isa(M,Function) ? M : x->M\x

	if isempty(x0)
		x0 = zeros(eltype(b),n)
		r = copy(b)
	else
		r = b - Af(x0)
	end

	z = Mf(r)
	p = copy(z)

	resvec = zeros(maxIter)

	iter   = 1
	flag   = 1

	@inbounds for iter=1:maxIter
		Ap = Af(p)
		gamma = dot(r,z)
		alpha = gamma/dot(p,Ap)

		# if A is ill-conditioned
		if alpha==Inf || alpha<0
			flag = 2; break
		end

		BLAS.axpy!(n,alpha,p,1,x0,1) # x = alpha*p+x
		BLAS.axpy!(n,-alpha,Ap,1,r,1) # r -= alpha*Ap

		resvec[iter] = norm(r)/nb
		# converge
		if  resvec[iter] <= tol
			flag = 0; break
		end

		z    = Mf(r)
		beta = dot(z,r)/gamma

		# p = z + beta*p
		p = BLAS.scal!(n,beta,p,1)
		p = BLAS.axpy!(n,1.0,z,1,p,1)
	end
	# control output
	if out>=0
		if flag == 1
			println("Does not converge within maximum iterations.")
		elseif flag == 2
			println("Matrix A in cg has to be positive definite.")
		end
	end
	return x0,flag, resvec[iter],  resvec, iter
end


text:
type FiniteDifferenceNewtonSafe <: Solver1D
  solverInfo::SolverInfo
end
FiniteDifferenceNewtonSafe(maxEvals::Int = 100, lowerBoundEnforced::Bool = false, upperBoundEnforced::Bool = false, lowerBound::Float64 = 0.0, upperBound::Float64 = 0.0) =
  FiniteDifferenceNewtonSafe(SolverInfo(maxEvals, lowerBoundEnforced, upperBoundEnforced, lowerBound, upperBound))

# Finite Differences Solver
function _solve(solver::FiniteDifferenceNewtonSafe, f::Function, accuracy::Float64, xMin::Float64, xMax::Float64, fxMin::Float64, fxMax::Float64,
                root::Float64, eval_num::Int)

  max_evals = solver.solverInfo.maxEvals
  # orienting search such that f(xl) < 0
  if fxMin < 0.0
    xl = xMin
    xh = xMax
  else
    xh = xMin
    xl = xMax
  end

  froot = f(root)
  eval_num += 1

  # first order finite difference derivative
  dfroot = xMax - root < root - xMin ? (fxMax - froot) / (xMax - root) : (fxMin - froot) / (xMin - root)

  dx = xMax - xMin
  while eval_num <= max_evals
    frootold = froot
    rootold = root
    dxold = dx
    # bisect if out of range or not decreasing fast enough
    if (((root - xh) * dfroot - froot) * ((root - xl) * dfroot - froot)) > 0.0 || abs(2.0 * froot) > abs(dxold * dfroot)
      dx = (xh - xl) / 2.0
      root = xl + dx
      # if root estimate computed is close to previous one, calculate dfroot at root ane xh rather than
      # root and rootold
      if is_close(root, rootold, 2500)
        rootold = xh
        frootold = f(xh)
      end
    else
      # Newton
      dx = froot / dfroot
      root -= dx
    end

    # Convergence check
    if abs(dx) < accuracy
      return root
    end

    froot = f(root)
    eval_num += 1
    dfroot = (frootold - froot) / (rootold - root)

    if (froot < 0.0)
      xl = root
    else
      xh = root
    end
  end

  error("Maximum number of function evals exceeded!")
end


text:
testfiles = [
    "misc.jl",
    "infer_geometry_based_bonds!.jl",
    "primitive_cell.jl"]

@assert VERSION.major == 1
@assert VERSION.minor ≥ 4

using XtalsPyTools, Test
XtalsPyTools.banner()

for testfile ∈ testfiles
    @info "Running test/$testfile"
    include(testfile)
end

@info "Done."

max_stars_repo_name: SimonEnsemble/XtalsPyTools

max_issues_repo_name: SimonEnsemble/XtalsPyTools