text:
# AbstractFunctional

# allow sending of arbitrary face/element-based data
global const PARALLEL_DATA_NONE = 1000
global const PARALLEL_DATA_FACE = 1001
global const PARALLEL_DATA_ELEMENT = 1002

"""
  Given the string describing what data is shared in parallel, returns the
  enum value.  An exception is thrown if an unrecognized value is supplied.

  **Inputs**

   * pdata: a string describing what data is shared in parallel

  **Outputs**

   * an Int enum
"""
function getParallelDataEnum(pdata::String)
  if pdata == "none"
    return PARALLEL_DATA_NONE
  elseif pdata == "face"
    return PARALLEL_DATA_FACE
  elseif pdata == "element"
    return PARALLEL_DATA_ELEMENT
  else
    error("unrecognized parallel data string: $pdata")
  end
end

"""
  Inverse of [`getParallelDataEnum`](@ref), takes the enum and returns the
  string.

  **Inputs**

   * pdata: Int enum

  **Outputs**

   * a String describing the parallel data
"""
function getParallelDataString(pdata::Int)
  if pdata == PARALLEL_DATA_NONE
    return "none"
  elseif pdata == PARALLEL_DATA_FACE
    return "face"
  elseif pdata == PARALLEL_DATA_ELEMENT
    return "element"
  else
    error("unrecognized parallel data enum: $pdata")
  end
end

@doc """
### ODLCommonTools.AbstractFunctional

Abstract datatype for optimization related data structures. All data types
corresponding to optimization problems should be a subtype of this.

  **Static Parameters**

   * Topt

"""->
abstract type AbstractFunctional{Topt} end

"""
  Abstract type for integral objective functions.  Any integral objective
  function should be a subtype of this, which is a subtype of
  [`AbstractFunctional`](@ref).

  **Static Parameters**

   * Topt
"""
abstract type AbstractIntegralFunctional{Topt} <: AbstractFunctional{Topt} end

"""
  Abstract type for functionals that are boundary integrals

  **Static Parameters**

   * Topt
"""
abstract type AbstractBoundaryFunctional{Topt} <: AbstractIntegralFunctional{Topt} end

#------------------------------------------------------------------------------
# API for functionals

# evalFunctional, evalFunctionalDeriv_q and evalFunctionalDeriv_m are defined
# in more detail in PDESolver

"""
  Evaluates a functional.  This function always returns a value of type
  `Topt`, the static parameter of the [`AbstractFunctional`](@ref) being
  evaluated.
  """
function evalFunctional
end

"""
  Implementation of [`evalFunctional`](@ref). Each new functional
  should implement this function, not `evalFunctional`.
  The final positional argument of this function *must* be
  the `AbstractFunctional` to be evaluated.
"""
function _evalFunctional
end

"""
  Performs reverse-mode differentiation of a functional with respect to the
  metrics
"""
function evalFunctionalDeriv_m
end

"""
  Implementation of [`evalFunctionalDeriv_m`](@ref).  Each new functional
  should implement this function, not `evalFunctionalDeriv_m`
  The final positional argument of this function *must* be
  the `AbstractFunctional` to be evaluated.

"""
function _evalFunctionalDeriv_m
end

"""
  Computes the derivative of a functional with respect to the solution `q`.
  For efficiency, this should use reverse-mode differentiation internally.
"""
function evalFunctionalDeriv_q
end

"""
  Implementation of [`evalFunctionalDeriv_q`](@ref).  Each new functional
  should implement this function, not `evalFunctionalDeriv_q`.
  The final two positional argument of this function *must* be
  the `AbstractFunctional` to be evaluated and the [`Abstract3DArray`](@ref)
  to be overwritten with the solution
"""
function _evalFunctionalDeriv_q
end

"""
  A setup function called immediately before the functional or its derivatives
  are calculated.  This provides the opportunity to cache any
  expensive-to-compute values inside the functional object that might be
  needed when evaluating the integrand.

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
"""
function setupFunctional(mesh::AbstractMesh, sbp,
                         eqn::AbstractSolutionData, opts::Dict,
                         functional::AbstractFunctional)

  # in general nothing to do here
  return nothing
end

"""
  Function that describes what kind of solution parallel data is required to
  evaluate this functional.  The function should return one of:

   * `PARALLEL_DATA_NONE`
   * `PARALLEL_DATA_FACE`
   * `PARALLEL_DATA_ELEMENT`

  Note that it is not necessary to implement this method for
  [`AbstractBoundaryFunctional`](@ref).

  **Inputs**

   * obj: the [`AbstractFunctional`](@ref)

  **Outputs**

   * one of the above strings
"""
function getParallelData(obj::AbstractFunctional)

  error("generic fallback for getParallelData(::AbstractFunctional) reached.  Did you forget to extend getParallelData with a new method for your functional?")

end

function getParallelData(obj::AbstractBoundaryFunctional)

  return PARALLEL_DATA_NONE
end

#------------------------------------------------------------------------------
# Composite Functional

"""
  Represents a functional that is the sum of other functionals.  All functional
  must have the same parallel data requirement.

  **Constructors**

  ```
  CompositeFunctional(func1, func2)
  ```

  where `func1` and `func2` are `AbstractFunctionals`](@ref), or

  ```
  CompositeFunctional([func1, func2])
  ```
"""
struct CompositeFunctional{Topt} <: AbstractFunctional{Topt}
  funcs::Vector{AbstractFunctional{Topt}}

  function CompositeFunctional{Topt}(funcs::AbstractVector) where {Topt}
    pdata = getParallelData(funcs[1])
    for i=2:length(funcs)
      if getParallelData(funcs[i]) != PARALLEL_DATA_NONE
        @assert getParallelData(funcs[i]) == pdata
      end
    end

    return new(funcs)
  end
end

function CompositeFunctional(funcs::Vararg{AbstractFunctional{Topt}}) where {Topt}

  return CompositeFunctional{Topt}(collect(funcs))
end

function CompositeFunctional(funcs::AbstractVector{AbstractFunctional{Topt}}) where {Topt}

  return CompositeFunctional{Topt}(funcs)
end

function getParallelData(func::CompositeFunctional)
  return getParallelData(func.funcs[1])
end

# In an ideal world, these would use varags to be fully generic, but
# Julia only allows varags on the final argument, so we have to specify the
# exact signature
#function _evalFunctional(args..., func::CompositeFunctional{Topt}; kwargs...)::Topt where {Topt}
function _evalFunctional(mesh::AbstractMesh{Tmsh},
                         sbp::AbstractOperator, eqn::AbstractSolutionData{Tsol}, opts,
                         func::CompositeFunctional{Topt})::Topt where {Tmsh, Tsol, Topt}

  J = 0
  for f in func.funcs
    J += _evalFunctional(mesh, sbp, eqn, opts, f)
  end

  return J
end

function _evalFunctionalDeriv_m(mesh::AbstractDGMesh{Tmsh},
                         sbp::AbstractOperator, eqn::AbstractSolutionData{Tsol}, opts,
                         func::CompositeFunctional{Topt},
                         val_bar::Number=1) where {Tmsh, Tsol, Topt}

  for f in func.funcs
    _evalFunctionalDeriv_m(mesh, sbp, eqn, opts, f, val_bar)
  end

  return nothing
end

function _evalFunctionalDeriv_q(mesh::AbstractDGMesh{Tmsh},
                         sbp::AbstractOperator, eqn::AbstractSolutionData{Tsol}, opts,
                         func::CompositeFunctional{Topt},
                         func_deriv_arr::Abstract3DArray) where {Tmsh, Tsol, Topt}

  func_deriv_arr2 = zeros(func_deriv_arr)
  for f in func.funcs
    _evalFunctionalDeriv_q(mesh, sbp, eqn, opts, f, func_deriv_arr2)
    for i=1:length(func_deriv_arr)
      func_deriv_arr[i] += func_deriv_arr2[i]
    end
    fill!(func_deriv_arr2, 0)
  end

  return nothing
end


text:
#-------------------------------------------------------------------------#
#                                                                         #
#                               Simulating_nets                           #
#                                                                         #
#-------------------------------------------------------------------------#

# V1 18 may 2020
#-------------------------------------------------------------------------#
# OBSERVATIONS:

# 1. Bla bla bla
# 2. Bla bla bla
# 3. Bla bla bla
#
#-------------------------------------------------------------------------#
# using PhyloNetworks

## FUNCTIONS

using Random

"""
    testing if a one_taxa string is valid
        This is a function in progress
"""
function is_valid_1taxa(one_taxa::String)
    test = true
    if(tryparse(Float64, SubString(one_taxa, 1,1)) !== nothing)
        error("taxa should not start with a number")
    end
    return test
end

"""
    Simulation of a phylogenetic tree topology
        We start with a two taxa tree and start to paste subtrees by choosing
        randomly taxa to paste.
        Currently the branch lengths are uniformly distributed (0,1)
        Returs a the topology in Newick format
"""
function random_tree_topology(taxa)
    shuffle!(taxa)
    # Starting the tree
    x = Any["(",taxa[1],",",taxa[2],")",";"]

    for i in 3:length(taxa)
        # taxa already in the tree and taxa not in the tree
        used_taxa = taxa[1:(i-1)]
        left_taxa = taxa[i:end]
        # choosing a taxa in the tree to add new branches
        taxa_to_change = rand(used_taxa)
        # possition to add branches
        position_of_change = findfirst(x.==taxa_to_change)
        # updating tree
        insert!(x,position_of_change+1,")")
        insert!(x,position_of_change,",")
        insert!(x,position_of_change,taxa[i])
        insert!(x,position_of_change,"(")
    end
    return x
end

"""
    Adds random branch lengths to a phylogenetic tree topology
        Currently the branch lengths are uniformly distributed (0,1)
        Returs a the topology in Newick format
"""
function add_random_bl(x,taxa)
    # We get the index of the taxa and right parenthesis (to add the branch lengths)
    idx_right_parenthesis = findall(x.==")")
    idx_taxa = zeros(Int64,length(taxa))
    for i in 1:length(taxa)
        idx_taxa[i] = findfirst(x.==taxa[i])
    end
    idx_taxa_and_right_p = vcat(idx_taxa, idx_right_parenthesis)
    # To paste the new branch length we need to sort in decreasing order
    sort!(idx_taxa_and_right_p, rev=true)
    # Random uniform branch length
    bl = rand(length(idx_taxa_and_right_p))
    # Adding the branch length
    for i in 1:length(idx_taxa_and_right_p)
        insert!(x,idx_taxa_and_right_p[i]+1,bl[i])
        insert!(x,idx_taxa_and_right_p[i]+1,":")
    end
    return x
end

"""
    Given a topology and branch lengths, this function change the branch lengths
    in a way that the tree is now (insert word to say that all the paths from a
    leaf to the root are of the same length)
        Only keeps the topology and with a stick break methodology gives the
        new branch lengths.
        Returs a the topology in Newick format
"""
function balance_a_tree(x,total_bl)
    # Case of taxa-colon-branch length-semicolon; eg A:0.32;
    if length(x)==4
        x[3] = total_bl[1]
        return x
    end

    # Case of two taxa
    if length(x)==12
      partition = rand(1)
      # bl = total_bl.*partition
      x[[4,8]] .= total_bl.*partition
      aux = total_bl.*(Float64[1]-partition)
      x[11] = aux[1]
      return x
    end

    # We need to find the place of the comma that divides the tree
    aux = (x.==")") - (x.=="(")
    pop!(aux)
    popfirst!(aux)
    central_comma_place = findfirst(cumsum(aux).==0) +4

    # Stuff of the left tree
    left_tree = vcat(x[2:(central_comma_place-1)],";")
    partition = rand(1)
    left_tree = balance_a_tree(left_tree,total_bl.*partition)

    # Stuff of the right tree
    right_tree = vcat(x[(central_comma_place+1):(length(x)-4)],";")
    right_tree[length(right_tree)-1] = left_tree[length(left_tree)-1]
    right_tree = balance_a_tree(right_tree,total_bl.*partition)

    # Pasting the left and right trees
    pop!(left_tree)
    pop!(right_tree)
    x = vcat("(",left_tree,",",right_tree,")",":",total_bl.*(Float64[1]-partition),";")

    return x
end

"""
    Addition of a random hybrid node to a phylogenetic tree
        We start a tree "x" with taxa "taxa" and select two taxa to add the hybrid
        node.
        Currently the branch length of the new hybrid and its gamma parameter are
        uniformly distributed (0,1). The first node that is selected we currently
        TRANSFORM it to a hybrid node. For the second node we only paste the
        other hybrid node.
        For instance, if we have the tree "((a,b),(c,d));" and select "b" and "c",
        we get "((a,#H1),((c)#H1,d));", without considering the branch lengths and
        the gamma parameter.
        Currently, sometimes the hybrid is in the same branch.
        Returs the tree with the added hybrid, in Newick format.

        Future improvements:
            - Extra variable with the name of the hybrid, now is H1
            - Correcting that in some simulations the hybrid nodes are in the same
            path an this probuces a weird hybrid that follows the same path as
            the original tree
"""
function add_random_hybrid(x,taxa)
   # sampling two taxa to insert the hybrid
   taxa_to_change1,taxa_to_change2 = taxa[randperm(length(taxa))[1:2]]
   # We get the branch length and the gamma parameter for the new hybrid
   bl_hybrid = rand(1)
   gamma = rand(1)
   # Finding the index of the taxa to change
   idx_node1 = findfirst(x.==taxa_to_change1)
   idx_node2 = findfirst(x.==taxa_to_change2)
   # We allways prefere the index of node 1 to be smaller
   if idx_node1 > idx_node2
      idx_node1,idx_node2 = idx_node2,idx_node1
   end
   # Now we can make the changes on the array x
   x[idx_node2] = "#H1"

   insert!(x,idx_node1+3,gamma[1])
   insert!(x,idx_node1+3,":")
   insert!(x,idx_node1+3,":")
   insert!(x,idx_node1+3,bl_hybrid[1])
   insert!(x,idx_node1+3,":")
   insert!(x,idx_node1+3,"#H1")
   insert!(x,idx_node1+3,")")
   insert!(x,idx_node1,"(")
   return x
end

"""
    Simulating a phylogenetic tree/network
        Returs a the topology in Newick format

        Future improvements:
            - To add the hybrid a node is eliminated, this could be solved by
            since the begining ask the name of the hybrid and past it as part of
            the tree, but remember that this change also affects the function
            add_random_hybrid; in particular, since this function chose randomly
            which node of the tree is transformed to hybrid.
"""
function Sim_tree(taxa,total_bl,is_network)
    # Getting the topology
    x = random_tree_topology(taxa)
    # Adding some branch lengths
    x = add_random_bl(x,taxa)
    # Fixing the branch lengths to be "equal"
    x = balance_a_tree(x,total_bl)
    # Adding a random hybrid
    if is_network
        x = add_random_hybrid(x,taxa)
    end
    return x
end

## Funcitons of njWillems2014

function Julia_willems2014(x,name_file_dist_txt,
                           name_file_Input_Matrix_Distances,
                           name_file_Parameters_Distances,
                           name_file_new_Network_Distances,
                           folder_output_distances_willems)
   # Folder to save the output of willems
   output_willems = folder_output_distances_willems*"/"*name_file_new_Network_Distances
   # Getting the netwokr from the Newick format
   net = readTopology(join(x))
   # Distance matrix
   Sigma = sharedPathMatrix(net)
   Vy = Sigma[:Tips]

   # From Cecil
   M = pairwiseTaxonDistanceMatrix(net, keepInternal=false)
   tipLabels(net)

   # Generating file to run njWillems2014
   ind = parse.(Int, tipLabels(net))
   # ind =  tipLabels(net)

   Mext = hcat(ind,M)
   n = length(tipLabels(net))

   df = DataFrame(Mext)
   df[!,:x1] = convert.(Int,df[!,:x1])

   CSV.write(name_file_dist_txt, df, writeheader=false, delim=' ')

   # the first row with the number of tips
   f=readlines(name_file_dist_txt)
   g = open(name_file_Input_Matrix_Distances,"w")
   write(g,"$n \n")

   for l in f
      write(g,l)
      write(g,"\n")
   end
   close(g)

   # This "run" generates the file "Network_Distances.txt"
   text_to_run = `./njwillems2014 $name_file_Input_Matrix_Distances $name_file_Parameters_Distances`
   run(text_to_run)

   # We copy the file "Network_Distances.txt", that generates njwillems2014
   cp("Network_Distances.txt", output_willems, force=true)

   net_dist = readdlm("Network_Distances.txt")
   # rm(name_file_new_Network_Distances)

   # txt_file = open("Network_Distances.txt","r")
   # readline(txt_file)
   # close(txt_file)

   # print(join(net_dist[1:10]))
   # println(net_dist)

   idx_Main_text = findfirst(net_dist .== "Main")
   println(net_dist[CartesianIndex(1, 1):(idx_Main_text)])

   # println((net_dist[1:20]))
   return net_dist
end
## End of functions

ext: jl

max_stars_repo_path: Backup/V01 10 aug 2020/Fn_PhyloSim.jl

max_stars_repo_name: ArrigoCoen/PhyloSim

max_issues_repo_path: Backup/V01 10 aug 2020/Fn_PhyloSim.jl

max_issues_repo_name: ArrigoCoen/PhyloSim

max_forks_repo_path: Backup/V01 10 aug 2020/Fn_PhyloSim.jl

max_forks_repo_name: ArrigoCoen/PhyloSim


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

using Markdown
using InteractiveUtils

# ╔═╡ 54fed622-0b1f-11ec-1987-2db84b3dfb4f
begin
    import Pkg
    # activate a temporary environment
    Pkg.activate(mktempdir())

	Pkg.add("ProfileSVG")
	Pkg.add(url="https://github.com/gdalle/PlutoProfile.jl")
    using PlutoProfile, ProfileSVG
end

# ╔═╡ 1f43857d-e3d2-4927-b2ba-03d983dee7cd
md"""
Define a dummy test function.
"""

# ╔═╡ 173cc0f5-5d57-4ce2-9d8f-72d5995db5c3
function profile_test(n)
    for i = 1:n
        A = randn(100,100,20)
        m = maximum(A)
        Am = mapslices(sum, A; dims=2)
        B = A[:,:,5]
        Bsort = mapslices(sort, B; dims=1)
        b = rand(100)
        C = B.*b
    end
end

# ╔═╡ 8b92bd93-9a5d-48bd-94e2-c57638571d67
md"""
Run it once to precompile.
"""

# ╔═╡ 943949fd-c12c-4f9b-b26a-48e4891431e9
profile_test(1)

# ╔═╡ 85ca4fea-386e-4d87-a7ba-b32263a6df0c
md"""
Profile and visualize the flame graph using `ProfileSVG.@profview`.
"""

# ╔═╡ f0afdc27-2d5d-418e-b0a6-c8960abfe58d
@profview profile_test(10)

# ╔═╡ ec802f21-bcff-4fc3-b391-fc12bb8265d2
md"""
Profile and visualize the flame graph using `PlutoProfile.@plutoprofview`.
"""

# ╔═╡ b57caa51-03fb-445c-9a36-ce5d8c912947
@plutoprofview profile_test(10)

# ╔═╡ Cell order:
# ╠═54fed622-0b1f-11ec-1987-2db84b3dfb4f
# ╟─1f43857d-e3d2-4927-b2ba-03d983dee7cd
# ╠═173cc0f5-5d57-4ce2-9d8f-72d5995db5c3
# ╟─8b92bd93-9a5d-48bd-94e2-c57638571d67
# ╠═943949fd-c12c-4f9b-b26a-48e4891431e9
# ╟─85ca4fea-386e-4d87-a7ba-b32263a6df0c
# ╠═f0afdc27-2d5d-418e-b0a6-c8960abfe58d
# ╟─ec802f21-bcff-4fc3-b391-fc12bb8265d2
# ╠═b57caa51-03fb-445c-9a36-ce5d8c912947


text:
L = 50
D = UncertainDataset([UncertainValue(Normal, rand(), rand()) for i = 1:L])
V = UncertainValueDataset([UncertainValue(Normal, rand(), rand()) for i = 1:L])
I = UncertainIndexDataset([UncertainValue(Normal, rand(), rand()) for i = 1:L])

@test randomshuffle(D) isa Vector{Float64}
@test randomshuffle(V) isa Vector{Float64}
@test randomshuffle(I) isa Vector{Float64}

@test randomphases(D) isa Vector{Float64}
@test randomphases(V) isa Vector{Float64}
@test randomphases(I) isa Vector{Float64}

@test randomamplitudes(D) isa Vector{Float64}
@test randomamplitudes(V) isa Vector{Float64}
@test randomamplitudes(I) isa Vector{Float64}

@test aaft(D) isa Vector{Float64}
@test aaft(V) isa Vector{Float64}
@test aaft(I) isa Vector{Float64}

@test iaaft(D) isa Vector{Float64}
@test iaaft(V) isa Vector{Float64}
@test iaaft(I) isa Vector{Float64}


text:
using Manifolds
@testset "Subgradient Plan" begin
    M = Euclidean(2)
    x = [4., 2.]
    x0 = [5.,2.]
    o = SubGradientMethodOptions(M, x0, StopAfterIteration(200), DecreasingStepsize(0.1))
    o.∂ = [1., 0.]
    f = y -> distance(M,y,x)
    ∂f = y -> distance(M,x,y) == 0 ? zero_tangent_vector(M,y) : -2*log(M,y,x)/distance(M,x,y)
    p = SubGradientProblem(M,f, ∂f)
    oR = solve(p,o)
    xHat = get_solver_result(oR)
    @test get_initial_stepsize(p,o) == 0.1
    @test get_stepsize(p,o,1) == 0.1
    @test get_last_stepsize(p,o,1) == 0.1
    # Check Fallbacks of Problen
    @test get_cost(p,x) == 0.
    @test norm(M,x,get_subgradient(p,x)) == 0
    @test_throws ErrorException get_gradient(p,o.x)
    @test_throws ErrorException getProximalMap(p,1.,o.x,1)
end


text:
# HITRAN constants

# reference temperature [K]
const c_T_ref = 296.

# reference pressure [atm]
const c_p_ref = 1.

# conversion factor atm to Pa
const c_atm_factor = 101325.

#=
Hartmann-Tran reference temperatures for line parameters

Reference:
    Wcisło P, Gordon IE, Tran H, Tan Y, Hu S-M, Campargue A, et al.
    The im- plementation of non-Voigt line profiles in the HITRAN
    database: H2 case study. J Quant Spectrosc Radiat Transf 2016;177:75–91.
    http://dx.doi.org/10.1016/j.jqsrt.2016.01.024.
=#
const c_HT_T_ref = Base.ImmutableDict(
    (0., 100.) => 50.,
    (100., 200.) => 150.,
    (200., 400.) => 296.,
    (400., Inf) => 700.)

## SI constants according to CODATA2018

# speed of light in m/s (SI)
const c_c_SI = 299792458.

# Avogrado constant (dimensionless) (SI)
const c_NA_SI = 6.02214076e23

# Boltzmann constant in J/K (SI)
const c_kB_SI = 1.380649e-23

# Planck-constant in J*s (SI)
const c_h_SI = 6.62607015e-34

# c2 = hc/k (cgs units)
const c_c2 = 1e2*c_h_SI * c_c_SI / c_kB_SI

# natural log of 2
const c_log2 = log(2.)

### profile constants
const c_default_zero = zero(Float64)

#=
Constants for saturation water vapor pressure taken from:

IAPWS R7-97(2012)
Revised Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam
August 2007
=#
const c_n1 = 0.11670521452767e4
const c_n2 = -0.72421316703206e6
const c_n3 = -0.17073846940092e2
const c_n4 = 0.12020824702470e5
const c_n5 = -0.32325550322333e7
const c_n6 = 0.14915108613530e2
const c_n7 = -0.48232657361591e4
const c_n8 = 0.40511340542057e6
const c_n9 = -0.23855557567849
const c_n10 = 0.65017534844798e3

#=

Constants for ice saturation vaper pressure taken from:

Huang, Jianhua. "A Simple Accurate Formula for Calculating
Saturation Vapor Pressure of Water and Ice", Journal of Applied Meteorology and
Climatology 57, 6 (2018): 1265-1272

=#

const c_a = 43.494
const c_b = 6545.8
const c_d1 = 278.
const c_d2 = 868.

const c_α = 1.00062
const c_β = 3.14e-8
const c_γ = 5.60e-7


text:
using Genie,Genie.Router
import Genie.Router: route
using Genie.Renderer
import Genie.Renderer.Json

Genie.config.run_as_server = true

route("/hello") do
    json("Hello World")
end

Genie.startup()


text:
include("src/base.jl")
include("src/envs/toytext_mdp.jl")
include("src/policies/evolution_policy.jl")
include("src/solvers/evolution_policy_solver.jl")

gym_env = ToyTextMDP(gym.make("FrozenLake-v0"))
solver = EvolutionPolicySolver(training_experiment_count = 500, experiment_repeats = 100)
policy = solve(solver, gym_env, verbose = true)

print(policy.action_map)

max_stars_repo_name: dmitrijsc/practical-rl

max_issues_repo_name: dmitrijsc/practical-rl

max_forks_repo_name: dmitrijsc/practical-rl


text:
# a basic create robot type node example
using Caesar, IncrementalInference

## Uncomment out for command line operation
cloudGraph, addrdict = standardcloudgraphsetup(nparticles=true)
session = addrdict["session"]
Nparticles = parse(Int,addrdict["num particles"])

# interactive operation
# session = "SESSCLOUDTEST"
# Nparticles = 100
# include(joinpath(dirname(@__FILE__),"blandauthremote.jl"))

N=Nparticles profroot
profroot

fg = Caesar.initfg(sessionname=session, cloudgraph=cloudGraph)

# Robot navigation and inference type stuff
doors = [-100.0;0.0;100.0;300.0]'
cov = [3.0]

# robot style, add first pose vertex
v1 = addNode!(fg,:x1,doors,N=N,labels=["POSE"])

# add a prior for the initial position of the robot
f0  = addFactor!(fg,[v1], Obsv2(doors, cov', [1.0]))

# add second pose vertex
tem = 2.0*randn(1,N)+getVal(v1)+50.0
v2 = addNode!(fg, :x2, tem, N=N,labels=["POSE"])

# now add the odometry factor between them
f1 = addFactor!(fg,[v1;v2],Odo([50.0]',[2.0]',[1.0]))

v3=addNode!(fg,:x3,4.0*randn(1,N)+getVal(v2)+50.0, N=N,labels=["POSE"])
addFactor!(fg,[v2;v3],Odo([50.0]',[4.0]',[1.0]))
f2 = addFactor!(fg,[v3], Obsv2(doors, cov', [1.0]))

inc =  fg.g.inclist
using Graphs
tt = fg.g.vertices[f0.index]
for vv in out_neighbors(tt, fg.g)
  @show vv.label
end

# set this part of graph to ready for solving
setDBAllReady!(conn, fg.sessionname)

println("Waiting for initial solve to occur")
sleep(10)

v4=addNode!(fg,:x4,2.0*randn(1,N)+getVal(v3)+50.0, N=N,labels=["POSE"])
addFactor!(fg,[v3;v4],Odo([50.0]',[2.0]',[1.0]))

if true
    l1=addNode!(fg, :l1, 0.5*randn(1,N)+getVal(v3)+64.0, N=N,labels=["LANDMARK"])
    addFactor!(fg, [v3,l1], Ranged([64.0],[0.5],[1.0]))
    addFactor!(fg, [v4,l1], Ranged([16.0],[0.5],[1.0]))
end

v5=addNode!(fg,:x5,2.0*randn(1,N)+getVal(v4)+50.0, N=N,labels=["POSE"])
addFactor!(fg,[v4;v5],Odo([50.0]',[2.0]',[1.0]))

v6=addNode!(fg,:x6,1.25*randn(1,N)+getVal(v5)+40.0, N=N,labels=["POSE"])
addFactor!(fg,[v5;v6],Odo([40.0]',[1.25]',[1.0]))

v7=addNode!(fg,:x7,2.0*randn(1,N)+getVal(v6) +60.0, N=N,labels=["POSE"])
addFactor!(fg,[v6;v7],Odo([60.0]',[2.0]',[1.0]))

f3 = addFactor!(fg,[v7], Obsv2(doors, cov', [1.0]))

# release the rest
setDBAllReady!(conn, fg.sessionname)

# Get neighbors
# neighs = CloudGraphs.get_neighbors(fg.cg, cv1r)

gt = Dict{String, Array{Float64,2}}()
gt[:x1]=([0.0;1.97304 ]')' # -0.0342366
gt[:x2]=([50.0; 2.83153 ]')' # 49.8797
gt[:x3]=([100.0; 1.65557 ]')' # 99.8351
gt[:x4]=([150.0; 1.64945 ]')' # 148.637
gt[:x5]=([200.0; 1.77992 ]')' # 198.62
gt[:x6]=([240.0; 2.20466 ]')' # 238.492
gt[:x7]=([300.0; 2.14353 ]')' # 298.467
gt[:l1]=([165.0; 1.17284 ]')' # 164.102


text:
euler_angles(p) = (n = norm(p); [0.0, acos(-p[3] / n), atan(-p[2], -p[1])])

# Find the actual crossing point
x0, y0 = TheoryOfCISS.find_crossing(
    1.23,
    2.39,
    14,
    threshold=1e-6,
    attenuation=1e-2
);

# Calculate the D and C vectors
vs = TheoryOfCISS.atpoint(
    x -> x[:v],
    x0,
    y0,
    N = 14,
    calc = TheoryOfCISS.calc_data1
);

D = 2real.(vs)
C = 2imag.(vs)

# Calculate the contribution to Λ12 from each atom
vxs, vys, vzs = [
    TheoryOfCISS.atpoint(
        x -> (xx -> xx[i]).(x[:vs]),
        x0 - 0.05,
        y0,
        N = 14,
        calc = TheoryOfCISS.calc_data2
    )
    for i in 1:3
];

ds = dxs, dys, dzs = real.(vxs), real.(vys), real.(vzs)
cs = cxs, cys, czs = imag.(vxs), imag.(vys), imag.(vzs)

λ = 8e-3

# Normalize using the total sum of contributions
N_d = sqrt(sum(dxs)^2 + sum(dys)^2 + sum(dzs)^2);
N_c = sqrt(sum(cxs)^2 + sum(cys)^2 + sum(czs)^2);
# Cumulatively add up contributions
cumulated_d = collect(zip(cumsum(dxs), cumsum(dys), cumsum(dzs)))
cumulated_c = collect(zip(cumsum(cxs), cumsum(cys), cumsum(czs)))

# Compute the scaling and Euler angles
scales = norm.(polarizations)
angles = euler_angles.(polarizations)

script = """
import bpy

electron_vector = bpy.context.scene.objects["ElectronVector"]

scales = [$(join(scales, ","))]
angles = [$(join(angles, ","))]

bpy.ops.screen.frame_jump(end=False)

frame = 10

for (scale, (x, y, z)) in zip(scales, angles):
    bpy.context.scene.frame_set(frame)
    electron_vector.scale[2] = scale
    electron_vector.rotation_euler[0] = x
    electron_vector.rotation_euler[1] = y
    electron_vector.rotation_euler[2] = z
    bpy.ops.anim.keyframe_insert()
    frame += 10

bpy.ops.screen.frame_jump(end=False)
"""

write("animate.py", script)


text:
# Readinfo types
# A container for storage of extracted fastq info
struct FastqInfo
	quality::Float64
	length::Int64
	gc::Float64
end

# A container for storage of extracted BAM info
struct BAMInfo
	quality::Float64
	length::Int64
	identity::Float64
	gc::Float64
end


text:
using Test
using Bonsai
using Bonsai: register!, Node, match_middleware, split_route
using URIs
using HTTP
using HTTP: Request
using AbstractTrees

@testset "split_route" begin
    @test split_route("/") == ["/"]
    @test split_route("/a/b") == split_route("/a/b/")
end

@testset "register!" begin
    r = Node("*")

    # 3 get routes
    register!(r, GET, "/fish/{id}", x -> 1)
    register!(r, GET, "/fish/super", x -> 2)
    register!(r, GET, "/fish/**", x -> 3)
    # 1 put
    register!(r, PUT, "/fish/**", x -> 4)
    # no fish route
    register!(r, PUT, "/turtle/super", x -> 5)

    p = "/fish/super"
    params = Dict()
    segments = split(p, "/"; keepempty=false)
    ms = match_middleware(r, params, "GET", segments, 1)
    @test sort(map(x -> x(nothing), ms)) == [1, 2, 3]

    match(r, "GET", p)

end

@testset "match" begin
    app = App()
    app.get("/{x}") do stream
    end
    req = Request()
    req.url = URI("http://locahost:4040/pet")
    req.method = "GET"
    # mutates the request storing the match in ctx
    match(app, req)
    @test haskey(req.context[:params], :x)
end

@testset "app" begin
    app = App()

    app.get("/") do stream

    end

    node = app.paths
    method = "GET"

    handler, params = Base.match(node,  method, "/")
    @test !isnothing(handler) && isempty(params)

    req = HTTP.Request()
    req.method = "GET"
    req.target = "/"
    handler, middleware = match(app, req)
    @test !isnothing(handler) && isempty(middleware)
end

@testset "app with middleware" begin
    app = App()

    app.get("**") do stream, next
    end

    app.get("/") do stream, next
    end

    app.get("**") do stream

    end

    req = HTTP.Request()
    req.method = "GET"
    req.target = "/"
    handler, middleware = match(app, req)

    @test !isnothing(handler)
    @test length(middleware) == 2
end

@testset "/" begin
    app = App()

    app.get("/") do stream
    end

    req = HTTP.Request()
    req.method = "GET"
    req.target = "/"
    handler, middleware = match(app, req)
    @test !isnothing(handler)
end


text:
function gethydrodata( dataobject::InfoType,
                        Nnvarh::Int,
                        nvarh_corr::Array{Int,1},
                        lmax::Real,
                        ranges::Array{Float64,1},
                        print_filenames::Bool,
                        read_cpu::Bool,
                        read_level::Bool)

    println("Reading data...")
# Narraysize::Int,

    kind = Float64 #data type

    xmin, xmax, ymin, ymax, zmin, zmax = ranges # 6 elements

    idom = zeros(Int32, 8)
    jdom = zeros(Int32, 8)
    kdom = zeros(Int32, 8)

    bounding_min = zeros(Float64, 8)
    bounding_max = zeros(Float64, 8)

    cpu_min = zeros(Int32, 8)
    cpu_max = zeros(Int32, 8)

    #vars_1D = zeros(kind, Narraysize, Nnvarh)
    #pos_1D = zeros(Int32, Narraysize, 4)

    vars_1D = ElasticArray{Float64}(undef, Nnvarh, 0) # to append multidim array
    if read_level # if AMR
        pos_1D = ElasticArray{Int}(undef, 4, 0) # to append multidim array
    else # if uniform grid
        pos_1D = ElasticArray{Int}(undef, 3, 0) # to append multidim array
    end

    #if read_cpu cpus_1D = zeros(Int, Narraysize) end
    cpus_1D = zeros(Int, 0) # zero size to use for append function

    path = dataobject.path
    overview = dataobject.grid_info
    nvarh = dataobject.nvarh
    cpu_overview = dataobject.grid_info

    twotondim=2^dataobject.ndim
    twotondim_float=2. ^dataobject.ndim

    xbound=[round( dataobject.grid_info.nx/2, RoundDown),
            round( dataobject.grid_info.ny/2, RoundDown),
            round( dataobject.grid_info.nz/2, RoundDown)]

    ngridfile = zeros(Int32, dataobject.ncpu+dataobject.grid_info.nboundary, dataobject.grid_info.nlevelmax)
    ngridlevel = zeros(Int32, dataobject.ncpu, dataobject.grid_info.nlevelmax)
    if dataobject.grid_info.nboundary > 0
        ngridbound = zeros(Int32, dataobject.grid_info.nboundary, dataobject.grid_info.nlevelmax)
    end

    dmax= maximum([xmax-xmin,ymax-ymin,zmax-zmin])

    ilevel = 1 # define ilevel before loop
    for il=1:lmax
        ilevel=il
        dx=0.5^ilevel
        if dx < dmax break end
    end

    bit_length=ilevel-1
    maxdom=2^bit_length

    #for positive valuse: Julia floor == Fotran int
    if bit_length > 0
        imin=floor(Int32, xmin * maxdom)
        imax=imin+1
        jmin=floor(Int32, ymin * maxdom)
        jmax=jmin+1
        kmin=floor(Int32, zmin * maxdom)
        kmax=kmin+1
    else
        imin=0
        imax=0
        jmin=0
        jmax=0
        kmin=0
        kmax=0
    end

    dkey=(2^(dataobject.grid_info.nlevelmax+1)/maxdom)^dataobject.ndim

    if bit_length>0 ndom=8 else ndom=1 end

    idom = [imin, imax, imin, imax, imin, imax, imin, imax]
    jdom = [jmin, jmin, jmax, jmax, jmin, jmin, jmax, jmax]
    kdom = [kmin, kmin, kmin, kmin, kmax, kmax, kmax, kmax]

    order_min=0.0e0 #todo: predeclare
    for i=1:ndom
        if bit_length > 0
            order_min = hilbert3d(idom[i],jdom[i],kdom[i],bit_length,1)
        else
            order_min=0.0e0
        end
        bounding_min[i]=(order_min)*dkey
        bounding_max[i]=(order_min+1.)*dkey
    end

    for impi=1:dataobject.ncpu
        for i=1:ndom
            if (dataobject.grid_info.bound_key[impi] <= bounding_min[i] &&
                dataobject.grid_info.bound_key[impi+1] > bounding_min[i])
                cpu_min[i]=impi
            end
            if (dataobject.grid_info.bound_key[impi] < bounding_max[i] &&
                dataobject.grid_info.bound_key[impi+1] >= bounding_max[i])
                cpu_max[i]=impi
            end
        end
    end

    cpu_read = copy(dataobject.grid_info.cpu_read)
    cpu_list = zeros(Int32, dataobject.ncpu)
    ncpu_read=Int32(0)
    for i=1:ndom #Todo
        for j=(cpu_min[i]):(cpu_max[i])
          if cpu_read[j]==false
              ncpu_read=ncpu_read+1
              cpu_list[ncpu_read]=j
              cpu_read[j]=true
          end
        end
    end

# -----------------------------------------------------------

    grid = fill(LevelType(0,0,0,0,0,0), lmax)
    # Compute hierarchy
    for ilevel=1:lmax
        nx_full=Int32(2^ilevel)
        ny_full=nx_full
        nz_full=nx_full

        imin=floor(Int32, xmin * nx_full) +1
        imax=floor(Int32, xmax * nx_full) +1
        jmin=floor(Int32, ymin * ny_full) +1
        jmax=floor(Int32, ymax * ny_full) +1
        kmin=floor(Int32, zmin * nz_full) +1
        kmax=floor(Int32, zmax * nz_full) +1

        grid[ilevel] = LevelType( imin, imax, jmin, jmax, kmin, kmax )
    end

    fnames = createpath(dataobject.output, path)

    xc = zeros(kind, 8,3)
    dummy1 = 0
    read_data= 0
    lmax2 =2^lmax
    var_firsttime =1
    var1 = 0.
    var_x = 0
    var_level = 0
    #get_var_totsize = 0

    @showprogress 1 "" for k=1:ncpu_read #Reading files...
        icpu=cpu_list[k]
        #println("icpu ",icpu)

        # Open AMR file and skip header
        ###if print_filenames println(Full_Path_AMR_cpufile) end
        amrpath = getproc2string(fnames.amr, icpu)
        f_amr = FortranFile(amrpath)

        # Open AMR file and skip header
        skiplines(f_amr, 21)

        # Read grid numbers
        read(f_amr, ngridlevel)
        ngridfile[1:dataobject.ncpu, 1:dataobject.grid_info.nlevelmax] = ngridlevel

        skiplines(f_amr, 1)

        if dataobject.grid_info.nboundary > 0
            skiplines(f_amr, 2)
            read(f_amr, ngridbound)
            ngridfile[(dataobject.ncpu+1):(dataobject.ncpu+overview.nboundary),1:dataobject.grid_info.nlevelmax]=ngridbound
        end

        skiplines(f_amr, 6)

        # Open HYDRO file and skip header
        hydropath = getproc2string(fnames.hydro, icpu)
        if print_filenames println(hydropath) end
        #println(Full_Path_Hydro_cpufile)
        f_hydro = FortranFile(hydropath)

        skiplines(f_hydro, 6)

        # Loop over levels
        for ilevel=1:lmax

            # Geometry
            dx=0.5^ilevel
            nx_full=Int32(2^ilevel)
            ny_full=nx_full
            nz_full=nx_full
            xc = geometry(twotondim_float, ilevel, xc)

            # Allocate work arrays
            ngrida=ngridfile[icpu,ilevel] # integer
            if ngrida>0
                xg  = zeros(kind, ngrida, dataobject.ndim)
                son = zeros(Int32, ngrida, twotondim)
                vara = zeros(kind, ngrida, twotondim, Nnvarh)
            end

            # Loop over domains
            for j=1:(overview.nboundary+dataobject.ncpu)
                # Read AMR data
                if ngridfile[j,ilevel]>0
                    skiplines(f_amr, 3)  # Skip grid index

                    # Read grid center
                    for idim=1:dataobject.ndim
                        if j == icpu
                            xg[:,idim] = read(f_amr, (kind, ngrida))
                        else
                            skiplines(f_amr, 1)
                        end
                    end

                    # Skip father index + Skip nbor index
                    skiplines(f_amr, 1 + (2*dataobject.ndim))

                    # Read son index
                    for ind=1:twotondim
                        if j == icpu
                            son[:,ind] = read(f_amr, (Int32, ngrida))
                        else
                            skiplines(f_amr, 1)
                        end
                    end

                    # Skip cpu map + refinement map
                    skiplines(f_amr, twotondim * 2)
                end

                # Read hydro data
                skiplines(f_hydro, 2)

                if ngridfile[j,ilevel]>0
                    # Read hydro variables
                    for ind=1:twotondim
                        for ivar=1:nvarh

                            if j == icpu
                                if nvarh_corr[ivar] != 0
                                    vara[:,ind,nvarh_corr[ivar]] = read(f_hydro,(kind, ngrida) )
                                else
                                    skiplines(f_hydro, 1)
                                end
                            else
                                skiplines(f_hydro, 1)
                            end

                        end
                    end

                end

            end

            # Compute map
            if ngrida>0
                vars_1D, pos_1D, cpus_1D = loopovercellshydro(twotondim,
                                                                ngrida, ilevel, lmax,
                                                                xg, xc, son, xbound,
                                                                nx_full, ny_full, nz_full,
                                                                grid, vara, vars_1D, pos_1D,
                                                                read_cpu, cpus_1D, k, read_level) # for kind = Float64
            end
        end # End loop over levels

        close(f_amr)
        close(f_hydro)

    end # End loop over cpu files

    if read_cpu
        return vars_1D, pos_1D, cpus_1D #arrays
    else
        return vars_1D, pos_1D #arrays
    end
end

function geometry(twotondim_float::Float64, ilevel::Int, xc::Array{Float64,2})
    dx=0.5^ilevel
    for (ind,iind) =enumerate(1.:twotondim_float)
        iiz=round( (iind-1)/4, RoundDown) #floor(Int32, (ind-1)/4)
        iiy=round( (iind-1-4*iiz)/2, RoundDown) #floor(Int32, (ind-1-4*iiz)/2)
        iix=round( (iind-1-2*iiy-4*iiz), RoundDown) #floor(Int32, (ind-1-2*iiy-4*iiz))

        xc[ind,1]=(iix-0.5)*dx
        xc[ind,2]=(iiy-0.5)*dx
        xc[ind,3]=(iiz-0.5)*dx
    end
    return xc
end

function loopovercellshydro(twotondim::Int,
                            ngrida::Int32,
                            ilevel::Int,
                            lmax::Int,
                            xg::Array{Float64,2},
                            xc::Array{Float64,2},
                            son::Array{Int32,2},
                            xbound::Array{Float64,1},
                            nx_full::Int32,
                            ny_full::Int32,
                            nz_full::Int32,
                            grid::Array{LevelType,1},
                            vara::Array{Float64,3},
                            vars_1D::ElasticArray{Float64,2,1},
                            pos_1D::ElasticArray{Int,2,1},
                            read_cpu::Bool,
                            cpus_1D::Array{Int,1},
                            k::Int,
                            read_level::Bool)
    for ind=1:twotondim # Loop over cells
        # Store data cube
        for i=1:ngrida
            if !(son[i,ind]>0 && ilevel<lmax) #if !ref[i]

                #=
                ix = round( (xg[i,1]+xc[ind,1]-xbound[1]) *nx_full, RoundDown)  +1.
                iy = round( (xg[i,2]+xc[ind,2]-xbound[2]) *ny_full, RoundDown)  +1.
                iz = round( (xg[i,3]+xc[ind,3]-xbound[3]) *nz_full, RoundDown)  +1.
                =#
                # todo: simplify "fortran int"

                ix=floor(Int, (xg[i,1]+xc[ind,1]-xbound[1]) *nx_full)+1
                iy=floor(Int, (xg[i,2]+xc[ind,2]-xbound[2]) *ny_full)+1
                iz=floor(Int, (xg[i,3]+xc[ind,3]-xbound[3]) *nz_full)+1

                #println(ilevel, ",    ", ix, " ", iy, " ", iz, "     ", grid[ilevel] )

                if      ix>=(grid[ilevel].imin) &&
                        iy>=(grid[ilevel].jmin) &&
                        iz>=(grid[ilevel].kmin) &&
                        ix<=(grid[ilevel].imax) &&
                        iy<=(grid[ilevel].jmax) &&
                        iz<=(grid[ilevel].kmax)
                    #println()
                    #println(ilevel, ",    ", ix, " ", iy, " ", iz, "     ", grid[ilevel] )

                    append!(vars_1D, vara[i,ind,:])
                    if read_level # if AMR
                        append!(pos_1D, [ix, iy, iz, ilevel])
                    else # if uniform grid
                        append!(pos_1D, [ix, iy, iz])
                    end
                    if read_cpu append!(cpus_1D, k) end
                end
            end
        end
    end # End loop over cell

    return vars_1D, pos_1D, cpus_1D
end


text:
#======================== PAWN ATTACK TABLES ==================================#
# Used internally to generate possible pawn advances
function pawnAdvance(pawns::Bitboard, targets::Bitboard, color::Color)
    targets & (color == WHITE ? (pawns << 8) : (pawns >> 8))
end

# Used internally to generate possible double pawn advances
function pawnDoubleAdvance(pawns::Bitboard, targets::Bitboard, color::Color)
    if color == WHITE
        pawnAdvance(pawnAdvance(pawns & RANK_2, targets, color), targets, color)
    else
        pawnAdvance(pawnAdvance(pawns & RANK_7, targets, color), targets, color)
    end
end

# Used internally to generate left captures
function pawnLeftCaptures(pawns::Bitboard, targets::Bitboard, color::Color)
    targets & (color == WHITE ? ((pawns & ~FILE_A) << 9) : ((pawns & ~FILE_H) >> 9))
end

# Used internally to generate right captures
function pawnRightCaptures(pawns::Bitboard, targets::Bitboard, color::Color)
    targets & (color == WHITE ? ((pawns & ~FILE_H) << 7) : ((pawns & ~FILE_A) >> 7))
end

# functions used to pregenerate possible pawn moves
function pawnAdvancesWhite(pawns::Bitboard, targets::Bitboard)
    pawnAdvance(pawns, targets, WHITE)
end

function pawnAdvancesBlack(pawns::Bitboard, targets::Bitboard)
    pawnAdvance(pawns, targets, BLACK)
end

function pawnCapturesWhite(pawns::Bitboard, targets::Bitboard)
    pawnLeftCaptures(pawns, targets, WHITE) | pawnRightCaptures(pawns, targets, WHITE)
end

function pawnCapturesBlack(pawns::Bitboard, targets::Bitboard)
    pawnLeftCaptures(pawns, targets, BLACK) | pawnRightCaptures(pawns, targets, BLACK)
end

# Precomputed pawn advance tables
const PAWN_ADVANCES_WHITE = @SVector [pawnAdvancesWhite(Bitboard(sqr), FULL) for sqr in 1:64]
const PAWN_ADVANCES_BLACK = @SVector [pawnAdvancesBlack(Bitboard(sqr), FULL) for sqr in 1:64]
const PAWN_CAPTURES_WHITE = @SVector [pawnCapturesWhite(Bitboard(sqr), FULL) for sqr in 1:64]
const PAWN_CAPTURES_BLACK = @SVector [pawnCapturesBlack(Bitboard(sqr), FULL) for sqr in 1:64]

"""
    pawnAttacks(c::Color, sqr::Integer)

Generate the potential pawn attacks from a square given by `sqr`, from the point of view of the color `c`. Returned as a `Bitboard`.
"""
function pawnAttacks(c::Color, sqr::Integer)
    if c == WHITE
        return @inbounds PAWN_CAPTURES_WHITE[sqr]
    else
        return @inbounds PAWN_CAPTURES_BLACK[sqr]
    end
end

"""
    pawnAdvances(c::Color, sqr::Integer)

Generate the potential pawn pushes from a square given by `sqr`, from the point of view of the color `c`. Returned as a `Bitboard`.
"""
function pawnAdvances(c::Color, sqr::Integer)
    if c == WHITE
        return @inbounds PAWN_ADVANCES_WHITE[sqr]
    else
        return @inbounds PAWN_ADVANCES_BLACK[sqr]
    end
end

#======================== KNIGHT ATTACK TABLES ================================#

# steps in all knight compass directions
knightMove_NNE(knight::Bitboard) = ~FILE_A & (knight << 15)
knightMove_ENE(knight::Bitboard) = ~FILE_AB & (knight << 6)
knightMove_ESE(knight::Bitboard) = ~FILE_AB & (knight >> 10)
knightMove_SSE(knight::Bitboard) = ~FILE_A & (knight >> 17)
knightMove_SSW(knight::Bitboard) = ~FILE_H & (knight >> 15)
knightMove_WSW(knight::Bitboard) = ~FILE_GH & (knight >> 6)
knightMove_WNW(knight::Bitboard) = ~FILE_GH & (knight << 10)
knightMove_NNW(knight::Bitboard) = ~FILE_H & (knight << 17)

# collate all moves
function knightMove_all(knight::Bitboard)
    knightMove_NNE(knight) | knightMove_ENE(knight) | knightMove_ESE(knight) |
    knightMove_SSE(knight) | knightMove_SSW(knight) | knightMove_WSW(knight) |
    knightMove_WNW(knight) | knightMove_NNW(knight)
end

# pre-compute
const KNIGHT_MOVES = @SVector [knightMove_all(Bitboard(one(UInt64) << i)) for i in 0:63]

"""
    knightMoves(sqr::Integer)

Generate the potential knight moves from a square given by `sqr`, returned as a `Bitboard`.
"""
knightMoves(sqr::Integer) = @inbounds KNIGHT_MOVES[sqr]

#======================== BISHOP ATTACK TABLES ================================#

bishopMove_NE(bishop::Bitboard) = ~FILE_A & (bishop << 7)
bishopMove_SE(bishop::Bitboard) = ~FILE_A & (bishop >> 9)
bishopMove_SW(bishop::Bitboard) = ~FILE_H & (bishop >> 7)
bishopMove_NW(bishop::Bitboard) = ~FILE_H & (bishop << 9)

const BISHOP_MOVE_FUNCTIONS = @SVector Function[bishopMove_NE, bishopMove_SE, bishopMove_SW, bishopMove_NW]

# initialise the tables on startup
"""
    BISHOP_TABLE

Constant containing precomputed bishop moves for all occupancy configurations.
"""
const BISHOP_TABLE = SVector{64}(initSlidingTable(Vector{Vector{Bitboard}}(undef, 64), BISHOP_MAGICS, BISHOP_MOVE_FUNCTIONS))

bishopMoves(sqr::Int, occupied::Bitboard) = @inbounds BISHOP_TABLE[sqr][subindex(occupied, BISHOP_MAGICS[sqr])]
bishopMoves(bb::Bitboard, occupied::Bitboard) = bishopMoves(square(bb), occupied)

const BISHOP_MOVES_EMPTY = @SVector [bishopMoves(i, EMPTY) for i in 1:64]

bishopMoves_empty(sqr::Int) = @inbounds BISHOP_MOVES_EMPTY[sqr]

#========================== ROOK ATTACK TABLES ================================#

rookMove_N(rook::Bitboard) = rook << 8
rookMove_S(rook::Bitboard) = rook >> 8
rookMove_E(rook::Bitboard) = ~FILE_A & (rook >> 1)
rookMove_W(rook::Bitboard) = ~FILE_H & (rook << 1)

const ROOK_MOVE_FUNCTIONS = @SVector Function[rookMove_N, rookMove_S, rookMove_E, rookMove_W]

# initialise the tables on startup
"""
    ROOK_TABLE

Constant containing precomputed rook moves for all occupancy configurations.
"""
const ROOK_TABLE = SVector{64}(initSlidingTable(Vector{Vector{Bitboard}}(undef, 64), ROOK_MAGICS, ROOK_MOVE_FUNCTIONS))

rookMoves(sqr::Integer, occupied::Bitboard) = @inbounds ROOK_TABLE[sqr][subindex(occupied, ROOK_MAGICS[sqr])]
rookMoves(bb::Bitboard, occupied::Bitboard) = rookMoves(square(bb), occupied)

const ROOK_MOVES_EMPTY = @SVector [rookMoves(i, EMPTY) for i in 1:64]

rookMoves_empty(sqr::Integer) = @inbounds ROOK_MOVES_EMPTY[sqr]

#========================== QUEEN ATTACK TABLES ===============================#

# Queen moves are just the union of all bishop moves, plus all rook moves.
queenMoves(sqr::T, occupied::Bitboard) where T <: Union{Integer, Bitboard} = rookMoves(sqr, occupied) | bishopMoves(sqr, occupied)

#=========================== KING ATTACK TABLES ===============================#

# steps in all compass directions
kingMove_N(king::Bitboard) = king << 8
kingMove_NE(king::Bitboard) = ~FILE_A & (king << 7)
kingMove_E(king::Bitboard) = ~FILE_A & (king >> 1)
kingMove_SE(king::Bitboard) = ~FILE_A & (king >> 9)
kingMove_S(king::Bitboard) = king >> 8
kingMove_SW(king::Bitboard) = ~FILE_H & (king >> 7)
kingMove_W(king::Bitboard) = ~FILE_H & (king << 1)
kingMove_NW(king::Bitboard) = ~FILE_H & (king << 9)

# collate all moves
function kingMove_all(king::Bitboard)
    kingMove_N(king) | kingMove_NE(king) | kingMove_E(king) | kingMove_SE(king) |
    kingMove_S(king) | kingMove_SW(king) | kingMove_W(king) | kingMove_NW(king)
end

# pre-compute table
const KING_MOVES = @SVector [kingMove_all(Bitboard(one(UInt64) << i)) for i in 0:63]

"""
    kingMoves(sqr::Integer)

Generate the potential king moves from a square given by `sqr`, returned as a `Bitboard`.
"""
kingMoves(sqr::Integer) = @inbounds KING_MOVES[sqr]

#=========================== Additional Utilities =============================#

"""
    squareAttackers(board::Board, sqr::Integer)

Find all the enemy units that attack a given square, `sqr`.
"""
function squareAttackers(board::Board, sqr::Integer)
    occ = occupied(board)
    attackers = (pawns(board) & pawnAttacks(board.turn, sqr)) |
        (knightMoves(sqr) & knights(board)) |
        (kingMoves(sqr) & kings(board))
    if !isempty(bishoplike(board))
        attackers |= (bishopMoves(sqr, occ) & bishoplike(board))
    end
    if !isempty(rooklike(board))
        attackers |= (rookMoves(sqr, occ) & rooklike(board))
    end
    attackers &= enemy(board)
end
# function squareAttackers(board::Board, sqr::Integer)
#     occ = occupied(board)
#     pvec = Vec{5, UInt64}((pawns(board).val, knights(board).val, kings(board).val, bishoplike(board).val, rooklike(board).val))
#     avec = Vec{5, UInt64}((pawnAttacks(board.turn, sqr).val, knightMoves(sqr).val, kingMoves(sqr).val, bishopMoves(sqr, occ).val, rookMoves(sqr, occ).val))
#     attackers = Bitboard(reduce(|, pvec & avec) & enemy(board).val)
# end
squareAttackers(board::Board, bb::Bitboard) = squareAttackers(board, square(bb))

function squaresquareAttackers_through_king(board::Board, sqr::Integer)
    occ = occupied(board)
    occ &= ~(kings(board) & friendly(board))
    attackers = (pawns(board) & pawnAttacks(board.turn, sqr)) |
        (knightMoves(sqr) & knights(board))
    if !isempty(bishoplike(board))
        attackers |= (bishopMoves(sqr, occ) & bishoplike(board))
    end
    if !isempty(rooklike(board))
        attackers |= (rookMoves(sqr, occ) & rooklike(board))
    end
    attackers &= enemy(board)
end
# function squaresquareAttackers_through_king(board::Board, sqr::Integer)
#     occ = occupied(board)
#     occ &= ~(kings(board) & friendly(board))
#     pvec = Vec{4, UInt64}((pawns(board).val, knights(board).val, bishoplike(board).val, rooklike(board).val))
#     avec = Vec{4, UInt64}((pawnAttacks(board.turn, sqr).val, knightMoves(sqr).val, bishopMoves(sqr, occ).val, rookMoves(sqr, occ).val))
#     attackers = Bitboard(reduce(|, pvec & avec) & enemy(board).val)
# end

"""
    isattacked(board::Board, sqr::Integer)

Is the given square attacked?
"""
function isattacked(board::Board, sqr::Integer)
    occ = occupied(board)
    enemies = enemy(board)
    isempty(pawns(board) & enemies & pawnAttacks(board.turn, sqr)) == false && return true
    isempty(knightMoves(sqr) & knights(board) & enemies) == false && return true
    isempty(bishopMoves(sqr, occ) & bishoplike(board) & enemies) == false && return true
    isempty(rookMoves(sqr, occ) & rooklike(board) & enemies) == false && return true
    isempty(kingMoves(sqr) & kings(board) & enemies) == false && return true
    return false
end

function isattacked_through_king(board::Board, sqr::Integer)
    occ = occupied(board)
    occ &= ~(kings(board) & friendly(board))
    enemies = enemy(board)
    isempty(pawns(board) & enemies & pawnAttacks(board.turn, sqr)) == false && return true
    isempty(knightMoves(sqr) & knights(board) & enemies) == false && return true
    isempty(bishopMoves(sqr, occ) & bishoplike(board) & enemies) == false && return true
    isempty(rookMoves(sqr, occ) & rooklike(board) & enemies) == false && return true
    isempty(kingMoves(sqr) & kings(board) & enemies) == false && return true
    return false
end

"""
    kingAttackers(board::Board, sqr::Integer)

Find all the attackers of a given square, `sqr`.
"""
function kingAttackers(board::Board, sqr::Integer)
    enemies = enemy(board)
    occ = occupied(board)
    attackers = (pawns(board) & pawnAttacks(board.turn, sqr)) |
        (knightMoves(sqr) & knights(board))
    if !isempty(bishoplike(board))
        attackers |= (bishopMoves(sqr, occ) & bishoplike(board))
    end
    if !isempty(rooklike(board))
        attackers |= (rookMoves(sqr, occ) & rooklike(board))
    end
    attackers &= enemy(board)
end
kingAttackers(board::Board, bb::Bitboard) = kingAttackers(board, square(bb))
kingAttackers(board::Board) = kingAttackers(board, square(kings(board) & friendly(board)))

# function to precompte the masks for blocking squares of a sliding attack
function initBlockerMasks(blockermasks::Array{Bitboard, 2})
    for sqr1 in 1:64
        for sqr2 in 1:64
            if isempty(rookMoves(sqr1, EMPTY) & sqr2) == false
                @inbounds blockermasks[sqr1, sqr2] = rookMoves(sqr1, Bitboard(sqr2)) & rookMoves(sqr2, Bitboard(sqr1))
            end
            if isempty(bishopMoves(sqr1, EMPTY) & sqr2) == false
                @inbounds blockermasks[sqr1, sqr2] = bishopMoves(sqr1, Bitboard(sqr2)) & bishopMoves(sqr2, Bitboard(sqr1))
            end
        end
    end
    return blockermasks
end

# precomputed blocker masks.
const BLOCKERMASKS = initBlockerMasks(fill(EMPTY, (64,64)))

"""
    blockers(sqr_1::Integer, sqr_2::Integer)

Retrieve all the blocking squares between `sqr_1` and `sqr_2`, as a `Bitboard`.
"""
blockers(sqr_1::Integer, sqr_2::Integer) = @inbounds BLOCKERMASKS[sqr_1, sqr_2]

"""
    pins(board::Board)

Retrieves all the pinned pieces on the board, as a `Bitboard`.
"""
function findpins(board::Board)
    king = square(kings(board) & friendly(board))
    occ = occupied(board)
    # Find all sliding pieces that xray through to the king.
    sliders = EMPTY
    if !isempty(bishoplike(board))
        sliders |= (bishopMoves_empty(king) & bishoplike(board))
    end
    if !isempty(rooklike(board))
        sliders |= (rookMoves_empty(king) & rooklike(board))
    end
    sliders &= enemy(board)
    # Find all times where there is only one blocker on the xray path.
    pinned = EMPTY
    for sqr in sliders
        blocking = blockers(sqr, king) & occ
        if isone(blocking) && (isempty(blocking & friendly(board)) === false)
            pinned |= blocking
        end
    end
    return pinned
end

max_stars_repo_name: J-Revell/Astellarn

max_issues_repo_name: J-Revell/Astellarn

max_forks_repo_name: J-Revell/Astellarn


text:
struct TransDateTimeCovariate{T<:AbstractFloat} <: AbstractCovariate{T}
    name::String
    basevariate::AbstractDateTimeVariate
    transform::Function
end

Base.length(var::TransDateTimeCovariate) = length(var.basevariate)

function TransDateTimeCovariate(name::String, basevariate::AbstractDateTimeVariate, transform::Function)
    TransDateTimeCovariate{Float32}(name, basevariate, transform)
end

function TransDateTimeCovariate(::Type{T}, name::String, basevariate::AbstractDateTimeVariate, transform::Function) where {T<:AbstractFloat}
    TransDateTimeCovariate{T}(name, basevariate, transform)
end

function slice(covariate::TransDateTimeCovariate{T}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {T<:AbstractFloat}
    basevar = covariate.basevariate
    f = covariate.transform
    g = (ms::Int64) -> begin
        ms == zero(Int64) ? T(NaN32) : begin
            dt = Dates.epochms2datetime(ms)
            convert(T, f(dt))
        end
    end
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    slices = slice(basevar, fromobs, toobs, slicelength)
    mapslice(g, slices, slicelength, T)
end

function Base.map(covariate::TransDateTimeCovariate, dataframe::AbstractDataFrame)
    TransDateTimeCovariate(covariate.name, map(covariate.basevariate, dataframe), covariate.transform)
end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
@testset "test_sum" begin

@test 0 == my_sum(0,0)

end


text:
module TypedIO
export toVec, fromVec, @typedIO, mkAltnertive, altnertiveToUnion

using MLStyle

function toVec end
function fromVec end

toVec(::Type{Any}, x::Symbol) = ["Symbol", string(x)]
toVec(::Type{Any}, x::Expr) = ["Expr", string(x)]
toVec(::Type{Any}, x) = string(x)

fromVec(::Type{Any}, x::String) = Some(x)
function fromVec(::Type{Any}, x::Vector)
    length(x) === 2 || return nothing
    if x[1] == "Expr"
        Some(Meta.parse(x[2]))
    elseif x[1] == "Symbol"
        Some(Symbol(x[2]))
    else
        nothing
    end
end

toVec(::Type{UInt64}, a::UInt64) where A <: Signed = ["uint", sizeof(a), string(a)]
function fromVec(::Type{UInt64}, a::Integer)
    Some(UInt64(a))
end

function fromVec(::Type{UInt64}, a::Vector)
    tag = a[1]
    bit = a[2]
    tag  == "uint" && bit == 8 || return nothing
    return Some(parse(UInt64, a[3]))
end

toVec(::Type{Integer}, a::A) where A <: Signed = ["int", sizeof(a), string(a)]
toVec(::Type{Integer}, a::A) where A <: Unsigned = ["uint", sizeof(a), string(a)]
toVec(::Type{AbstractFloat}, a::A) where A <: AbstractFloat = ["float", sizeof(a), string(a)]

function fromVec(::Type{Integer}, a::Integer)
    Some(a)
end

function fromVec(::Type{Integer}, a::Vector)
    tag = a[1]
    tag in ("int", "uint") || return nothing
    bit = a[2]
    s = a[3]
    if tag == "int"
        if bit === 1
            return Some(parse(Int8, s))
        elseif bit === 2
            return Some(parse(Int16, s))
        elseif bit === 4
            return Some(parse(Int32, s))
        elseif bit === 8
            return Some(parse(Int64, s))
        else
            return nothing
        end
    else
        if bit === 1
            return Some(parse(UInt8, s))
        elseif bit === 2
            return Some(parse(UInt16, s))
        elseif bit === 4
            return Some(parse(UInt32, s))
        elseif bit === 8
            return Some(parse(UInt64, s))
        else
            return nothing
        end
    end
end

function fromVec(::Type{AbstractFloat}, a::AbstractFloat)
    Some(a)
end
function fromVec(::Type{AbstractFloat}, a::Vector)
    tag = a[1]
    tag == "float" || return nothing
    bit = a[2]
    s = a[3]
    if bit === 2
        return Some(parse(Float16, s))
    elseif bit === 4
        return Some(parse(Float32, s))
    elseif bit === 8
        return Some(parse(Float64, s))
    else
        return nothing
    end
end

toVec(::Type{Nothing}, ::Nothing) = nothing
fromVec(::Type{Nothing}, ::Nothing) = Some(nothing)

toVec(::Type{I}, x::I) where I <: Union{Int8, Int16, Int32, UInt8, UInt16, UInt32, Int64} = I(x)
fromVec(::Type{I}, x::Integer) where I <: Union{Int8, Int16, Int32, UInt8, UInt16, UInt32,  Int64} = Some(I(x))

toVec(::Type{Float64}, x::Float64) = x
fromVec(::Type{Float64}, x::Float64) = Some(x)

toVec(::Type{F}, x::F) where F <: AbstractFloat = x
fromVec(::Type{F}, x::AbstractFloat) where F <: AbstractFloat = Some(F(x))

toVec(::Type{String}, x::String) = x
fromVec(::Type{String}, x::String) = Some(x)

toVec(::Type{Char}, x::Char) = String([x])
fromVec(::Type{Char}, x::String) = Some(x[1])
fromVec(::Type{Char}, x::Char) = Some(x)

toVec(::Type{Symbol}, x::Symbol) = string(x)
fromVec(::Type{Symbol}, x::String) = Some(Symbol(x))

toVec(::Type{Bool}, x::Bool) = x
fromVec(::Type{Bool}, x::Bool) = Some(x)

toVec(::Type{Pair{A, B}}, x::Pair) where {A, B} = [toVec(A, x.first), toVec(B, x.second)]
function fromVec(::Type{Pair{A, B}}, x::Vector) where {A, B}
    length(x) === 2 || return nothing
    l = fromVec(A, x[1])
    l === nothing && return nothing
    r = fromVec(B, x[2])
    r === nothing && return nothing
    Some(Pair{A, B}(l.value, r.value))
end

toVec(::Type{Vector{E}}, x::Vector) where E = [toVec(E, each) for each in x]
fromVec(::Type{Vector{E}}, x::Vector) where E =
    let res = altnertiveToUnion(E)[]
        for v in x
            e = fromVec(E, v)
            e === nothing && return nothing
            push!(res, e.value)
        end
       Some(res)
    end

abstract type Altnertive end
struct UnionSig{A, B} <: Altnertive end
struct UnionEnd <: Altnertive end

toVec(::Type{UnionSig{A, B}}, x::A) where {A, B<:Altnertive} = toVec(A, x)
toVec(::Type{UnionSig{A, B}}, x) where {A, B <: Altnertive} = toVec(B, x)
toVec(::Type{UnionSig{A, UnionEnd}}, x::A) where A = toVec(A, x)

fromVec(t, a) = begin
    nothing

end

function fromVec(::Type{UnionSig{A, B}}, x) where {A, B<:Altnertive}
    l = fromVec(A, x)
    l === nothing && return fromVec(B, x)
    l
end

fromVec(::Type{UnionSig{A, UnionEnd}}, x) where A = fromVec(A, x)

function compileVariant(m::Module, T, targs)
    subpatterns = []
    yields = [Symbol(:_, i) for i = targs]
    fs = fieldnames(m.eval(T))
    @assert length(targs) == length(fs)

    tag = string(T)
    contentsToDict =  Any[tag]
    contentsFromDict = []
    n_targs = length(targs)
    for i in 1:n_targs
        targ = targs[i]
        targ = @match targ begin
            :(Union{$(uargs...)}) && if m.eval(Union) == Union end =>
                foldr(uargs, init=UnionEnd) do each, last
                    :($UnionSig{$each, $last})
                end
            _ => targ
        end
        arg = :(subject.$(fs[i]))
        push!(subpatterns, arg)
        push!(contentsToDict, :($toVec($targ, $arg)))
        push!(contentsFromDict,
            quote
                $(yields[i]) = $fromVec($targ, subject[$(i + 1)])
                $(yields[i]) === nothing && return nothing
                $(yields[i]).value
            end)
    end
    to = T => :[$(contentsToDict...)]
    from = tag => :($length(subject) !== $(1 + n_targs) ? nothing : Some($T($(contentsFromDict...))))
    to, from

end

function asTypedVec(m::Module, absT::Any, fields::AbstractVector)
    variantsToDict = []
    variantsFromDict = []
    for field in fields
        @switch field begin
        @case :($T($(targs...)))
            to, from = compileVariant(m, T, targs)
            push!(variantsToDict,  to)
            push!(variantsFromDict, from)
            nothing
        @case ::LineNumberNode
            nothing
        end
    end

    variantsToDictExpr = foldr(variantsToDict, init = Expr(:call, error, "unknown $absT")) do (T, ret), last
        Expr(:if, :(subject isa $T), ret, last)
    end

    variantsFromDictExpr = foldr(variantsFromDict, init = Expr(:call, error, "cannot parse to $absT")) do (tag, ret), last
        Expr(:if, :(subject[1] == $tag), ret, last)
    end

    quote
        function $TypedIO.fromVec(::Type{$absT}, subject::Vector)
            $variantsFromDictExpr
        end

        function $TypedIO.toVec(::Type{$absT}, subject::$absT)
            $variantsToDictExpr
        end
    end
end

function typedIORecord(__module__, tname, targs)
    (_, to), (tag, from) = compileVariant(__module__, tname, targs)
    quote
        function $TypedIO.fromVec(::Type{$tname}, subject::Vector)
            $tag !== subject[1] && return nothing
            $from
        end

        function $TypedIO.toVec(::Type{$tname}, subject::$tname)
            $to
        end
    end
end

function typedIO(ex, __module__)
    @match ex begin
        :($absT = $(Expr(:block, stmts...))) => TypedIO.asTypedVec(__module__, absT, stmts)
        :($T($(targs...))) => typedIORecord(__module__, T, targs)
        :($absT => $(from::Symbol)) =>
            quote
                @inline function $TypedIO.toVec(::Type{$absT}, subject::$absT)
                    $TypedIO.toVec(::Type{$from}, subject)
                end
                @inline function $TypedIO.fromVec(::Type{$absT}, subject)
                    $TypedIO.fromVec(::Type{$from}, subject)
                end
            end
    end
end

macro typedIO(ex)
    esc(typedIO(ex, __module__))
end

@static if VERSION < v"1.5"
    @typedIO LineNumberNode(Int64, Union{Nothing, String})
else
    @typedIO LineNumberNode(Int64, Union{Nothing, Symbol})
end

function mkAltnertive(head::Type, ts::Type...)
    isempty(ts) && return head
    foldr(ts, init = UnionSig{head, UnionEnd}) do each, last
        UnionSig{each, last}
    end
end

altnertiveToUnion(::Type{UnionSig{A, UnionEnd}}) where A = A
altnertiveToUnion(::Type{UnionSig{A, Tl}}) where {A, Tl} = Union{A, altnertiveToUnion(Tl)}
altnertiveToUnion(::Type{UnionEnd}) = Union{}
altnertiveToUnion(a) = a

const exprCommon = mkAltnertive(Expr, QuoteNode, Float64, Int, LineNumberNode, Char, Bool, Symbol, Nothing)

@typedIO QuoteNode(exprCommon)

@typedIO Expr(Symbol, Vector{exprCommon})

end

# using .TypedIO
# using MLStyle

# @data A begin
#     A1(Int, Int)
#     A2()
# end

# @typedIO A = begin
#     A1(Int, Int)
#     A2()
# end

# struct C
#     a :: Int
# end
# @typedIO C(Int)

# @data B begin
#     B1(A, Union{Nothing, A})
#     B2(Vector{A}, C)
# end

# @typedIO B = begin
#     B1(A, Union{Nothing, A})
#     B2(Vector{A}, C)
# end

# println(toVec(A, A1(1, 2)))
# println(fromVec(A, toVec(A, A1(1, 2))))

# a1 = A1(1, 2)
# a2 = A1(3, 4)
# a3 = A2()

# b1 = B2(A[a1, a2, a3], C(1))
# b2 = B1(a2, a3)
# b3 = B1(a2, nothing)

# println(toVec(Vector{B}, [b1, b2, b3]))
# println(fromVec(Vector{B}, toVec(Vector{B}, [b1, b2, b3])))


text:
# Nastran BIGMAT column limit
const COL = 65536
const iCOL = 1/COL

const MAT_ALL = "ALL"
const MAT_NAMES = "NAMES"

# PCH constants
const TITLE = "TITLE"
const SUBTITLE = "SUBTITLE"
const LABEL = "LABEL"
const SUBCASE_ID = "SUBCASE ID"
const EIGENVALUE = "EIGENVALUE"
const MODE = "MODE"
const POINT = "POINT ID"

const INCLUDE = "INCLUDE"
const SINGLEP = 7
const DOUBLEP = 15

const PCHL = 18
const PCHLL = 72
const DOLLAR = '\$'
const CONT = "-CONT-"
const EQUAL = '='

const DOF = 6

max_stars_repo_name: mgcth/SDNastran

max_issues_repo_name: mgcth/SDNastran

max_forks_repo_name: mgcth/SDNastran


text:
"""
translate([-24,0,0]) {
    union() {
        cube(15, center=true);
        sphere(10);
    }
}

intersection() {
    cube(15, center=true);
    sphere(10);
}

translate([24,0,0]) {
    difference() {
        cube(15, center=true);
        sphere(10);
    }
}
"""

using Descartes
using GeometryBasics

a = translate([-24,0,0])union(
        Cuboid([15,15,15], center=true),
        Descartes.Sphere(10))

b = intersect(
        Cuboid([15,15,15], center=true),
        Descartes.Sphere(10))

c = translate([24,0,0])diff(
        Cuboid([15,15,15], center=true),
        Descartes.Sphere(10))

m = Mesh(a,b,c)
@show typeof(m)
using GLMakie

scene = GLMakie.mesh(coordinates(m), faces(m); shading=true, color=:yellow)


text:
using Pkg
#Pkg.activate("/home/nagj/.julia/environments/sparse_discrete")
include("convexRelaxations.jl")

# TODO: edit function description to describe zero_indices and one_indices params
function solve_sparse_problem(U_tilde, mu, k_sparse;
    zero_indices=[], one_indices=[])
    """
    This function exactly solves the problem of approximating (in squared
    frobenius norm) an arbitrary matrix U_tilde by a sparse matrix and under a
    regularization penalty.

    :param U_tilde: An arbitrary n-by-n matrix.
    :param mu: The regularization parameter for the sparse matrix penalty
               (Float64).
    :param k_sparse: The maximum sparsity of the sparse matrix (Int64).

    :return: The n-by-n sparse matrix that solves the problem.
    """
    S = construct_binary_matrix(U_tilde, k_sparse,
        zero_indices = zero_indices, one_indices = one_indices)

    Y = (S .* U_tilde) ./ (1 + mu)

    return Y

end;

function solve_rank_problem(U_tilde, lambda, k_rank; exact_svd=true)
    """
    This function exactly solves the problem of approximating (in squared
    frobenius norm) an arbitrary matrix U_tilde by a low rank matrix and under
    a regularization penalty.

    :param U_tilde: An arbitrary n-by-n matrix.
    :param lambda: The regularization parameter for the low rank matrix penalty
                   (Float64).
    :param k_rank: The maximum rank of the low rank matrix (Int64).

    :return: The n-by-n low rank matrix that solves the problem.
    """
    return project_matrix(U_tilde ./ (1 + lambda), k_rank, exact_svd=exact_svd)

end;

function solve_rank_relaxed_problem(U, mu, lambda, k_sparse, k_rank;
    check_feasibility=true)
    """
    This function exactly solves the problem of interest in the absence of the
    low rank constraint. It then rounds this lower bound to a feasible solution.

    :param U: An arbitrary n-by-n matrix.
    :param mu: The regularization parameter for the sparse matrix penalty
                   (Float64).
    :param lambda: The regularization parameter for the low rank matrix penalty
                   (Float64).
    :param k_sparse: The maximum sparsity of the sparse matrix (Int64).
    :param k_rank: The maximum rank of the low rank matrix (Int64).

    :return: This function returns two values. The first value is a tuple of two
             n-by-n arrays that correspond to the solution of the rank
             constraint relaxed problem (the first element in the tuple is the
             matrix X and the second element is the matrix Y). The second value
             is a tuple of two n-by-n arrays that correspond to the rounded
             feasible solution to the original problem (the first element in the
             tuple is the matrix X and the second element is the matrix Y).
    """
    S = construct_binary_matrix(U, k_sparse)

    Y = S .* U
    Y = lambda * Y / (lambda + mu + mu * lambda)

    X = (U - Y) / (1 + lambda)

    # Round to a feasible solution
    projected_X = project_matrix(X, k_rank)

    if check_feasibility
        @assert is_feasible(projected_X, Y, k_sparse, k_rank)
    end

    return (X, Y), (projected_X, Y)

end;

function solve_sparse_relaxed_problem(U, mu, lambda, k_sparse, k_rank;
    check_feasibility=true)
    """
    This function exactly solves the problem of interest in the absence of the
    sparsity constraint. It then rounds this lower bound to a feasible solution.

    :param U: An arbitrary n-by-n matrix.
    :param mu: The regularization parameter for the sparse matrix penalty
                   (Float64).
    :param lambda: The regularization parameter for the low rank matrix penalty
                   (Float64).
    :param k_sparse: The maximum sparsity of the sparse matrix (Int64).
    :param k_rank: The maximum rank of the low rank matrix (Int64).

    :return: This function returns two values. The first value is a tuple of two
             n-by-n arrays that correspond to the solution of the sparsity
             constraint relaxed problem (the first element in the tuple is the
             matrix X and the second element is the matrix Y). The second value
             is a tuple of two n-by-n arrays that correspond to the rounded
             feasible solution to the original problem (the first element in the
             tuple is the matrix X and the second element is the matrix Y).
    """
    projected_U = project_matrix(U, k_rank)

    X = mu * projected_U / (lambda + mu + mu * lambda)

    Y = (U - X) / (1 + mu)

    # Round to a feasible solution
    S = construct_binary_matrix(Y, k_sparse)
    rounded_Y = S .* Y

    if check_feasibility
        @assert is_feasible(X, rounded_Y, k_sparse, k_rank)
    end

    return (X, Y), (X, rounded_Y)

end;

function local_search(X_init, Y_init, U, mu, lambda, k_sparse, k_rank;
    update_sparsity_pattern = true, min_improvement = 0.001)
    """
    This function implements an alternating minimization local search method to
    arrive at a locally optimal solution. Starting from an initial feasible
    solution (X_init, Y_init), this function parametrizes X as X = V * W^T and
    parametrizes Y as Y = S .* (U - V * W^T) / (1 + mu). The function then
    iteratively optimizes over V, W and S until the objective value converges.

    :param X_init: An n-by-n matrix with rank(X_init) <= k_rank.
    Lparam Y_init: An n-by-n matrix with ||Y||_0 <= k_sparse.
    :param U: An arbitrary n-by-n matrix.
    :param mu: The regularization parameter for the sparse matrix penalty
               (Float64).
    :param lambda: The regularization parameter for the low rank matrix penalty
                   (Float64).
    :param k_sparse: The maximum sparsity of the sparse matrix (Int64).
    :param k_rank: The maximum rank of the low rank matrix (Int64).
    :param min_improvement: The minimal fractional decrease in the objective
                            value required for the minimization procedure to
                            continue iterating.

    :return: This function returns 4 values. The first return value is a tuple
             of two n-by-n arrays that correspond to the locally optimal
             feasible solution (the first element in the tuple is the matrix X
             and the second element is the matrix Y). The second return value
             is the objective value achieved by the locally optimal solution
             (Float64). The third return value is the number of iterations
             performed (Int64). The final return value is the total fractional
             improvement of the final objective value compared to the initial
             objective value (Float64).
    """
    # Verify that the initial solution (X_init, Y_init) is a feasible solution.
    @assert is_feasible(X_init, Y_init, k_sparse, k_rank)

    # Initialize V, Wt and S
    svd_out = svd(X_init)
    svd_U = svd_out.U
    svd_S = svd_out.S
    svd_Vt = svd_out.Vt

    initial_V = svd_U[:, 1:k_rank] * Diagonal(svd_S[1:k_rank].^0.5)
    initial_Wt = Diagonal(svd_S[1:k_rank].^0.5) * svd_Vt[1:k_rank, :]
    S = abs.(Y_init) .> 1e-10

    old_objective = compute_objective_value(X_init, Y_init, U, mu, lambda)
    init_objective = old_objective

    # Take an initial step in the local search procedure
    V = optimize_V(initial_Wt, S, U, mu, lambda, k_rank)
    Wt = optimize_Wt(V, S, U, mu, lambda, k_rank)
    X_new = V * Wt
    if update_sparsity_pattern
        S = construct_binary_matrix(U - X_new, k_sparse)
    end

    Y_new = S .* (U - X_new)
    Y_new = Y_new / (1 + mu)
    new_objective = compute_objective_value(X_new, Y_new, U, mu, lambda)
    step_count = 1

    # If the last step improved the solution more than min_provement,
    # take another step
    while((old_objective - new_objective) / old_objective > min_improvement)

        V = optimize_V(Wt, S, U, mu, lambda, k_rank)
        Wt = optimize_Wt(V, S, U, mu, lambda, k_rank)
        X_new = V * Wt
        if update_sparsity_pattern
            S = construct_binary_matrix(U - X_new, k_sparse)
        end

        Y_new = S .* (U - X_new)
        Y_new = Y_new / (1 + mu)

        old_objective = new_objective
        new_objective = compute_objective_value(X_new, Y_new, U, mu, lambda)
        step_count = step_count + 1

    end

    total_improvement = (init_objective - new_objective) / init_objective

    return (X_new, Y_new), new_objective, step_count, total_improvement

end;

function optimize_V(Wt, S, U, mu, lambda, k_rank; solver_output=0)
    """
    This function finds the n-by-k_rank matrix V that minimizes the following
    objective:
        ||U - V * Wt -  S .* (U - V * WT) / (1 + mu)||_F^2 +
        lambda * ||V * Wt||_F^2 + mu * ||S .* (U - V * Wt) / (1 + mu)||_F^2

    :param Wt: An arbitrary k_rank-by-n matrix.
    :param S: An arbitrary n-by-n binary matrix.
    :param U: An arbitrary n-by-n matrix.
    :param mu: The regularization parameter for the sparse matrix penalty
               (Float64).
    :param lambda: The regularization parameter for the low rank matrix penalty
                   (Float64).
    :param k_rank: The maximum rank of the low rank matrix (Int64).

    :return: The optimal n-by-k_rank matrix V.
    """
    n = size(U)[1]

    # Construct the optimization problem in JuMP
    m = Model(with_optimizer(Gurobi.Optimizer, GUROBI_ENV))
    set_optimizer_attribute(m, "OutputFlag", solver_output)

    @variable(m, V[i=1:n, j=1:k_rank])
    @variable(m, X[i=1:n, j=1:n])
    @variable(m, Y[i=1:n, j=1:n])

    @constraint(m, X .== V * Wt)
    @constraint(m, [i = 1:n, j = 1:n],
        Y[i, j] == S[i, j] * (U[i, j]-X[i, j]) / (1+mu))

    @objective(m, Min, sum((U[i, j]-X[i, j]-Y[i, j])^2 for i=1:n for j=1:n) +
        lambda * sum(X[i, j]^2 for i=1:n for j=1:n) +
        mu * sum(Y[i, j]^2 for i=1:n for j=1:n))

    optimize!(m)

    return value.(V)

end;

function optimize_Wt(V, S, U, mu, lambda, k_rank; solver_output=0)
    """
    This function finds the K-rank-by-n matrix Wt that minimizes the following
    objective:
        ||U - V * Wt -  S .* (U - V * WT) / (1 + mu)||_F^2 +
        lambda * ||V * Wt||_F^2 + mu * ||S .* (U - V * Wt) / (1 + mu)||_F^2

    :param V: An arbitrary n-by-k_rank matrix.
    :param S: An arbitrary n-by-n binary matrix.
    :param U: An arbitrary n-by-n matrix.
    :param mu: The regularization parameter for the sparse matrix penalty
               (Float64).
    :param lambda: The regularization parameter for the low rank matrix penalty
                   (Float64).
    :param k_rank: The maximum rank of the low rank matrix (Int64).

    :return: The optimal k_rank-by-n matrix Wt.
    """
    n = size(U)[1]

    # Construct the optimization problem in JuMP
    m = Model(with_optimizer(Gurobi.Optimizer, GUROBI_ENV))
    set_optimizer_attribute(m, "OutputFlag", solver_output)

    @variable(m, Wt[i=1:k_rank, j=1:n])
    @variable(m, X[i=1:n, j=1:n])
    @variable(m, Y[i=1:n, j=1:n])

    @constraint(m, X .== V * Wt)
    @constraint(m, [i = 1:n, j = 1:n],
        Y[i, j] == S[i, j] * (U[i, j]-X[i, j]) / (1+mu))

    @objective(m, Min, sum((U[i, j]-X[i, j]-Y[i, j])^2 for i=1:n for j=1:n) +
        lambda * sum(X[i, j]^2 for i=1:n for j=1:n) +
        mu * sum(Y[i, j]^2 for i=1:n for j=1:n))

    optimize!(m)

    return value.(Wt)

end;


text:
# Generate right hand coefficients
function mycoeffs(f, n, lam)
    F = Fun(f, Ultraspherical(lam))
    m = coefficients(F)
    length(m) < n ? (return [m; zeros(n-length(m))]) : (return m[1:n])
end


text:
include("primer74graph.jl")
include("primer74BFS.jl")

G = emptyGraph(8)
G.V[1].ime = "Mesa"
G.V[1].prezime = "Selimovic"
G.V[1].pol = 'm'
G.V[2].ime = "Laza"
G.V[2].prezime = "Lazarevic"
G.V[2].pol = 'm'
G.V[3].ime = "Jovan"
G.V[3].prezime = "Jovanovic"
G.V[3].pol = 'm'
G.V[4].ime = "Isidora"
G.V[4].prezime = "Sekulic"
G.V[4].pol = 'z'
G.V[5].ime = "Dragoljub"
G.V[5].prezime = "Mihailovic"
G.V[5].pol = 'm'
G.V[6].ime = "Borislav"
G.V[6].prezime = "Pekic"
G.V[6].pol = 'm'
G.V[7].ime = "Ivo"
G.V[7].prezime = "Andric"
G.V[7].pol = 'm'
G.V[8].ime = "Danica"
G.V[8].prezime = "Maksimovic"
G.V[8].pol = 'z'
G.AdjMatrix = [ 0 1 0 0 0 0 1 0;
                1 0 1 0 1 0 0 1;
                0 1 0 0 1 0 1 0;
                0 0 0 0 0 1 1 1;
                0 1 1 0 0 1 0 1;
                0 0 0 1 1 0 1 0;
                1 0 1 1 0 1 0 0;
                0 1 0 1 1 0 0 0]
ime = "Danica"
prezime = "Maksimovic"

sug = primer74BFS(G, ime, prezime)

println("Sugestije za korisnika $ime $prezime su:")
for i in 1:length(sug)
  println("$i. $(sug[i].ime) $(sug[i].prezime)")
end


text:
module EMLcondenser
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...stage_separators.EMLtank
	import EMLtypes.length
	include("condenser/condenserSteady.jl")
	include("condenser/condenserSteady_fakeH.jl")
	include("condenser/condenserReact.jl")
	include("condenser/condenser.jl")
	include("condenser/condenser2.jl")
	include("condenser/condenserSubcooled.jl")
end


text:
module Kwonly
export @add_kwonly

try
    using Core: UndefKeywordError
catch
    # UndefKeywordError added in: RFC: required keyword arguments
    # https://github.com/JuliaLang/julia/pull/25830/files
    struct UndefKeywordError <: Exception
        var::Symbol
    end
    Base.showerror(io::IO, ex::UndefKeywordError) =
        print(io, "UndefKeywordError: keyword argument $(ex.var) not assigned")
end

"""
    @add_kwonly function_definition

Define keyword-only version of the `function_definition`.

    @add_kwonly function f(a, b; c=1, d=2)
        ...
    end

expands to:

    function f(x; y=1)
        ...
    end
    function f(; x = throw(UndefKeywordError(:x)), y=1)
        ...
    end
"""
macro add_kwonly(ex)
    esc(add_kwonly(ex))
end

add_kwonly(ex::Expr) = add_kwonly(Val{ex.head}, ex)

function add_kwonly(::Type{<: Val}, ex)
    error("add_kwonly does not work with expression $(ex.head)")
end

function add_kwonly(::Union{Type{Val{:function}},
                            Type{Val{:(=)}}}, ex::Expr)
    body = ex.args[2:end]  # function body
    default_call = ex.args[1]  # e.g., :(f(a, b=2; c=3))
    kwonly_call = add_kwonly(default_call)
    if kwonly_call === nothing
        return ex
    end

    return quote
        begin
            $ex
            $(Expr(ex.head, kwonly_call, body...))
        end
    end
end

function add_kwonly(::Type{Val{:where}}, ex::Expr)
    default_call = ex.args[1]
    rest = ex.args[2:end]
    kwonly_call = add_kwonly(default_call)
    if kwonly_call === nothing
        return nothing
    end
    return Expr(:where, kwonly_call, rest...)
end

function add_kwonly(::Type{Val{:call}}, default_call::Expr)
    # default_call is, e.g., :(f(a, b=2; c=3))
    funcname = default_call.args[1]  # e.g., :f
    required = []  # required positional arguments; e.g., [:a]
    optional = []  # optional positional arguments; e.g., [:(b=2)]
    default_kwargs = []
    for arg in default_call.args[2:end]
        if isa(arg, Symbol)
            push!(required, arg)
        elseif arg.head == :(::)
            push!(required, arg)
        elseif arg.head == :kw
            push!(optional, arg)
        elseif arg.head == :parameters
            @assert default_kwargs == []  # can I have :parameters twice?
            default_kwargs = arg.args
        else
            error("Not expecting to see: $arg")
        end
    end
    if isempty(required) && isempty(optional)
        # If the function is already keyword-only, do nothing:
        return nothing
    end
    if isempty(required)
        # It's not clear what should be done.  Let's not support it at
        # the moment:
        error("At least one positional mandatory argument is required.")
    end

    kwonly_kwargs = Expr(:parameters, [
        Expr(:kw, pa, :(throw($UndefKeywordError($(QuoteNode(pa))))))
        for pa in required
    ]..., optional..., default_kwargs...)
    kwonly_call = Expr(:call, funcname, kwonly_kwargs)
    # e.g., :(f(; a=error(...), b=error(...), c=1, d=2))

    return kwonly_call
end

end  # module


text:
"Create a new feature definition object to hold the field definitions."
createfeaturedefn(name::AbstractString) = GDAL.fd_create(name)

"Destroy a feature definition object and release all memory associated with it"
destroy(featuredefn::FeatureDefn) = GDAL.destroy(featuredefn)

"Drop a reference, and destroy if unreferenced."
release(featuredefn::FeatureDefn) = GDAL.release(featuredefn)

"Get name of the OGRFeatureDefn passed as an argument."
getname(featuredefn::FeatureDefn) = GDAL.getname(featuredefn)

"Fetch number of fields on the passed feature definition."
nfield(featuredefn::FeatureDefn) = GDAL.getfieldcount(featuredefn)

"Fetch field definition of the passed feature definition."
fetchfielddefn(featuredefn::FeatureDefn, i::Integer) =
    GDAL.getfielddefn(featuredefn, i)

"""
Find field by name.

### Returns
the field index, or -1 if no match found.
"""
getfieldindex(featuredefn::FeatureDefn, name::AbstractString) =
    GDAL.getfieldindex(featuredefn, name)

"Add a new field definition to the passed feature definition."
addfielddefn(featuredefn::FeatureDefn, fielddefn::FieldDefn) =
    GDAL.addfielddefn(featuredefn, fielddefn)

"Delete an existing field definition."
function deletefielddefn(featuredefn::FeatureDefn, i::Integer)
    result = GDAL.deletefielddefn(featuredefn, i)
    if result != GDAL.OGRERR_NONE
        error("Failed to delete field $i in the feature definition")
    end
end

"""
Reorder the field definitions in the array of the feature definition.

### Parameters
* **featuredefn**: handle to the feature definition.
* **indices**: an array of GetFieldCount() elements which is a permutation of
    [0, GetFieldCount()-1]. panMap is such that, for each field definition at
    position i after reordering, its position before reordering was panMap[i].
"""
function reorderfielddefns(featuredefn::FeatureDefn, indices::Vector{Cint})
    result = GDAL.reorderfielddefns(featuredefn, indices)
    if result != GDAL.OGRERR_NONE
        error("Failed to delete field $i in the feature definition")
    end
end

"Fetch the geometry base type of the passed feature definition."
getgeomtype(featuredefn::FeatureDefn) = GDAL.getgeomtype(featuredefn)

"""Assign the base geometry type for the passed layer (the same as the
feature definition)."""
setgeomtype(featuredefn::FeatureDefn, etype::GDAL.OGRwkbGeometryType) =
    GDAL.setgeomtype(featuredefn, etype)

"Determine whether the geometry can be omitted when fetching features."
isgeomignored(featuredefn::FeatureDefn) =
    GDAL.isgeometryignored(featuredefn)

"Set whether the geometry can be omitted when fetching features."
setgeomignored(featuredefn::FeatureDefn, ignore::Bool) =
    GDAL.setgeometryignored(featuredefn, ignore)

"Determine whether the style can be omitted when fetching features."
isstyleignored(featuredefn::FeatureDefn) =
    Bool(GDAL.isstyleignored(featuredefn))

# """
#     OGR_FD_SetStyleIgnored(OGRFeatureDefnH hDefn,
#                            int bIgnore) -> void
# Set whether the style can be omitted when fetching features.
# ### Parameters
# * **hDefn**: handle to the feature definition on witch OGRFeature are based on.
# * **bIgnore**: ignore state
# """
# function setstyleignored(arg1::Ptr{OGRFeatureDefnH},arg2::Integer)
#     ccall((:OGR_FD_SetStyleIgnored,libgdal),Void,(Ptr{OGRFeatureDefnH},Cint),arg1,arg2)
# end

# """
#     OGR_FD_Reference(OGRFeatureDefnH hDefn) -> int
# Increments the reference count by one.
# ### Parameters
# * **hDefn**: handle to the feature definition on witch OGRFeature are based on.
# ### Returns
# the updated reference count.
# """
# function reference(arg1::Ptr{OGRFeatureDefnH})
#     ccall((:OGR_FD_Reference,libgdal),Cint,(Ptr{OGRFeatureDefnH},),arg1)
# end

# """
#     OGR_FD_Dereference(OGRFeatureDefnH hDefn) -> int
# Decrements the reference count by one.
# ### Parameters
# * **hDefn**: handle to the feature definition on witch OGRFeature are based on.
# ### Returns
# the updated reference count.
# """
# function dereference(arg1::Ptr{OGRFeatureDefnH})
#     ccall((:OGR_FD_Dereference,libgdal),Cint,(Ptr{OGRFeatureDefnH},),arg1)
# end

# """
#     OGR_FD_GetReferenceCount(OGRFeatureDefnH hDefn) -> int
# Fetch current reference count.
# ### Parameters
# * **hDefn**: handle to the feature definition on witch OGRFeature are based on.
# ### Returns
# the current reference count.
# """
# function getreferencecount(arg1::Ptr{OGRFeatureDefnH})
#     ccall((:OGR_FD_GetReferenceCount,libgdal),Cint,(Ptr{OGRFeatureDefnH},),arg1)
# end

"Fetch number of geometry fields on the passed feature definition."
ngeomfield(featuredefn::FeatureDefn) = GDAL.getgeomfieldcount(featuredefn)

"Fetch geometry field definition of the passed feature definition."
fetchgeomfielddefn(featuredefn::FeatureDefn, i::Integer) =
    GDAL.getgeomfielddefn(featuredefn, i)

"""
Find geometry field by name.

### Returns
the geometry field index, or -1 if no match found.
"""
getgeomfieldindex(featuredefn::FeatureDefn, name::AbstractString) =
    GDAL.getgeomfieldindex(featuredefn, name)

"Add a new field definition to the passed feature definition."
add(featuredefn::FeatureDefn, geomfielddefn::GeomFieldDefn) =
    GDAL.addgeomfielddefn(featuredefn, geomfielddefn)

"Delete an existing geometry field definition."
function deletegeomfielddefn(featuredefn::FeatureDefn, i::Integer)
    result = GDAL.deletegeomfielddefn(featuredefn, i)
    if result != GDAL.OGRERR_NONE
        error("Failed to delete geometry field $i in the feature definition")
    end
end

"Test if the feature definition is identical to the other one."
issame(featuredefn1::FeatureDefn, featuredefn2::FeatureDefn) =
    Bool(GDAL.issame(featuredefn1, featuredefn2))

"""Returns the new feature object with null fields and no geometry

Starting with GDAL 2.1, returns NULL in case out of memory situation.
"""
function createfeature(featuredefn::FeatureDefn)
    result = GDAL.f_create(featuredefn)
    (result == C_NULL) && error("out of memory when creating feature")
    result
end

"Destroy the feature passed in."
destroy(feature::Feature) = GDAL.destroy(feature)

"Fetch feature definition."
getfeaturedefn(feature::Feature) = GDAL.getdefnref(feature)