text:
"""
SymbolicVector
"""

mutable struct SymbolicVector{Tv} <: AbstractVectorSymbolic{Tv}
    params::Set{Symbol}
    vec::Vector{<:AbstractNumberSymbolic{Tv}}
    dim::Int
end

mutable struct SymbolicVectorExpression{Tv} <: AbstractVectorSymbolic{Tv}
    params::Set{Symbol}
    op::Symbol
    args::Vector{<:AbstractSymbolic{Tv}}
    dim::Int
    others
end

"""
convert
"""

function Base.convert(::Type{<:AbstractSymbolic{T}}, vec::Vector{<:AbstractNumberSymbolic{T}}) where {T<:Number}
    s = union([u.params for u = vec]...)
    SymbolicVector{T}(s, vec, length(vec))
end

function Base.convert(::Type{<:AbstractSymbolic{T}}, s::SymbolicVector{S}) where {T<:Number,S<:Number}
    v = [convert(AbstractNumberSymbolic{T}, x) for x = s.vec]
    SymbolicVector{T}(s.params, v, length(v))
end

"""
promotion
"""

function Base.promote_rule(::Type{<:AbstractVectorSymbolic{T}}, ::Type{S}) where {T<:Number,S<:Number}
    AbstractVectorSymbolic{promote_type(T,S)}
end

function Base.promote_rule(::Type{<:AbstractVectorSymbolic{T}}, ::Type{<:AbstractVectorSymbolic{S}}) where {T<:Number,S<:Number}
    AbstractVectorSymbolic{promote_type(T,S)}
end

"""
IO for structures
"""

function Base.show(io::IO, x::AbstractVectorSymbolic)
    expr = _toexpr(x)
    Base.show(io, expr)
end

function _toexpr(x::SymbolicVector)
    args = [_toexpr(e) for e = x.vec]
    Expr(:vect, args...)
end

function _toexpr(x::SymbolicVectorExpression)
    args = [_toexpr(e) for e = x.args]
    Expr(:call, x.op, args...)
end

"""
operators
"""

function Base.getindex(x::SymbolicVector{T}, i::Number) where {T<:Number}
    x.vec[i]
end

function Base.getindex(x::SymbolicVector{T}, i::Any) where {T<:Number}
    v = getindex(x.vec, i)
    s = union([u.params for u = v]...)
    SymbolicVector{T}(s, v, length(v))
end

function Base.getindex(x::AbstractVectorSymbolic{T}, i::Number) where {T<:Number}
    SymbolicExpression{T}(x.params, :getindex, [x], [i])
end

function Base.getindex(x::AbstractVectorSymbolic{T}, i::Any) where {T<:Number}
    SymbolicVectorExpression{T}(x.params, :getindex, [x], length(i), [i])
end

# function Base.setindex!(x::SymbolicVector{T}, v::AbstractNumberSymbolic{T}, i::Number) where {T<:Number}
#     x.vec[i] = v
# end

# """
# for macro
# """

# function _toexpr(x::Vector{<:AbstractSymbolic{Tv}}) where Tv
#     args = [_toexpr(e) for e = x]
#     Expr(:vect, args...)
# end

"""
seval(f, env, cache)
Return the value for expr f
"""

function seval(s::SymbolicVector{Tv}, env::SymbolicEnv, cache::SymbolicCache)::Vector{Tv} where {Tv <: Number}
    Tv[seval(x, env, cache) for x = s.vec]
end

function seval(f::AbstractVectorSymbolic{Tv}, env::SymbolicEnv, cache::SymbolicCache)::Vector{Tv} where Tv
    get(cache, f) do
        retval = _eval(Val(f.op), f, env, cache)
        cache[f] = retval
    end
end

"""
seval(f, dvar, env, cache)
Return the first derivative of expr f with respect to dvar
"""

function seval(s::SymbolicVector{Tv}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::Vector{Tv} where {Tv <: Number}
    Tv[seval(x, dvar, env, cache) for x = s.vec]
end

function seval(f::AbstractVectorSymbolic{Tv}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::Vector{Tv} where Tv
    (dvar in f.params) || return zeros(Tv, f.dim)
    get(cache, (f,dvar)) do
        retval = _eval(Val(f.op), f, dvar, env, cache)
        cache[(f,dvar)] = retval
    end
end

"""
seval(f, dvar, env, cache)
Return the second derivative of expr f with respect to dvar1 and dvar2
"""

function seval(s::SymbolicVector{Tv}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::Vector{Tv} where {Tv <: Number}
    Tv[seval(x, dvar, env, cache) for x = s.vec]
end

function seval(f::AbstractVectorSymbolic{Tv}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::Vector{Tv} where Tv
    (dvar[1] in f.params) || (dvar[2] in f.params) || return zeros(Tv, f.dim)
    get(cache, (f,dvar)) do
        retval = _eval(Val(f.op), f, dvar, env, cache)
        cache[(f,dvar)] = retval
    end
end

"""
seval for vec
"""

function seval(vec::Vector{<:AbstractNumberSymbolic{Tv}}, env::SymbolicEnv, cache::SymbolicCache)::Vector{Tv} where {Tv <: Number}
    seval(convert(SymbolicVector{Tv}, vec), env, cache)
end

function seval(vec::Vector{<:AbstractNumberSymbolic{Tv}}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::Vector{Tv} where {Tv <: Number}
    seval(convert(SymbolicVector{Tv}, vec), dvar, env, cache)
end

function seval(vec::Vector{<:AbstractNumberSymbolic{Tv}}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::Vector{Tv} where {Tv <: Number}
    seval(convert(SymbolicVector{Tv}, vec), dvar, env, cache)
end

"""
getindex1
"""

function _eval(::Val{:getindex}, f::SymbolicExpression{Tv}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    x = seval(f.args[1], env, cache)
    getindex(x, f.others[1])
end

function _eval(::Val{:getindex}, f::SymbolicExpression{Tv}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    x = seval(f.args[1], dvar, env, cache)
    getindex(x, f.others[1])
end

function _eval(::Val{:getindex}, f::SymbolicExpression{Tv}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    x = seval(f.args[1], dvar, env, cache)
    getindex(x, f.others[1])
end

"""
getindex
"""

function _eval(::Val{:getindex}, f::SymbolicVectorExpression{Tv}, env::SymbolicEnv, cache::SymbolicCache)::Vector{Tv} where Tv
    x = seval(f.args[1], env, cache)
    getindex(x, f.others[1])
end

function _eval(::Val{:getindex}, f::SymbolicVectorExpression{Tv}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::Vector{Tv} where Tv
    x = seval(f.args[1], dvar, env, cache)
    getindex(x, f.others[1])
end

function _eval(::Val{:getindex}, f::SymbolicVectorExpression{Tv}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::Vector{Tv} where Tv
    x = seval(f.args[1], dvar, env, cache)
    getindex(x, f.others[1])
end

"""
plus
"""

function _eval(::Val{:+}, f::SymbolicVectorExpression{Tv}, env::SymbolicEnv, cache::SymbolicCache)::Vector{Tv} where Tv
    args = [seval(x, env, cache) for x = f.args]
    +(args...)
end

function _eval(::Val{:+}, f::SymbolicVectorExpression{Tv}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::Vector{Tv} where Tv
    args = [seval(x, dvar, env, cache) for x = f.args]
    +(args...)
end

function _eval(::Val{:+}, f::SymbolicVectorExpression{Tv}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::Vector{Tv} where Tv
    args = [seval(x, dvar, env, cache) for x = f.args]
    +(args...)
end

"""
minus
"""

function _eval(::Val{:-}, f::SymbolicVectorExpression{Tv}, env::SymbolicEnv, cache::SymbolicCache)::Vector{Tv} where Tv
    args = [seval(x, env, cache) for x = f.args]
    -(args...)
end

function _eval(::Val{:-}, f::SymbolicVectorExpression{Tv}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::Vector{Tv} where Tv
    args = [seval(x, dvar, env, cache) for x = f.args]
    -(args...)
end

function _eval(::Val{:-}, f::SymbolicVectorExpression{Tv}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::Vector{Tv} where Tv
    args = [seval(x, dvar, env, cache) for x = f.args]
    -(args...)
end


text:
## Auto-correlation
plot(
    autocor(chains[1].planets[1].e, 1:100:1000)
)

## Trace plots
plot()
plot!(chains[1].planets[1].a[1:end], title="d")
plot!(chains[1].planets[2].a[1:end], title="e")
## Trace plots
plot()
plot!(chains[1].planets[1].e[1:end], title="d")
plot!(chains[1].planets[2].e[1:end], title="e")

##
mean(getproperty.(stats[1], :acceptance_rate))
mean(getproperty.(stats[1], :tree_depth))
maximum(getproperty.(stats[1], :tree_depth))
count(getproperty.(stats[1], :numerical_error))
mean(getproperty.(stats[1], :n_steps))
maximum(getproperty.(stats[1], :log_density))

##
cd(@__DIR__)

ext: jl


text:
function getEnergyList(hit)
        nParticles = icxx"$(hit)->getNMCParticles();"
        energyList = Float64[]
        for i in 1:nParticles
                push!(energyList, icxx"$(hit)->getEnergyCont($(i)-1);")
        end
        return energyList
end

function getTimeList(hit)
        nParticles = icxx"$(hit)->getNMCParticles();"
        timeList = Float64[]
        for i in 1:nParticles
                push!(timeList, icxx"$(hit)->getTimeCont($(i)-1);")
        end
        return timeList
end

function getPDGList(hit)
        nParticles = icxx"$(hit)->getNMCParticles();"
        PDGList = Float64[]
        for i in 1:nParticles
                push!(PDGList, icxx"$(hit)->getPDGCont($(i)-1);")
        end
        return PDGList
end


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni.
#

function print_JSON_render(
     png_output::String,
     pfm_output::Union{String, Nothing},
     scenefile::String,
     time_of_start::String,
     algorithm::Renderer,
     camera::Camera,
     α::Float64,
     width::Int64,
	height::Int64,
     samples_per_side::Int64,
     declare_float::Union{Dict, Nothing},
     rendering_time_s::Float64,
     a::Float64,
     γ::Float64,
     lum::Union{Number, Nothing},
     )

     dict_camera = if typeof(camera) == OrthogonalCamera
               Dict(
                    "projection" => "orthogonal",
                    "aspect_ratio" => camera.a,
                    "transformation" => camera.T.M,
               )

          elseif typeof(camera) == PerspectiveCamera
               Dict(
                    "projection" => "orthogonal",
                    "distance" => camera.d,
                    "aspect ratio" => camera.a,
                    "transformation" => camera.T.M,
               )
          end

     dict_renderer = if typeof(algorithm) == OnOffRenderer
              Dict(
                    "algorithm" => "On-Off Renderer",
                    "background color" => algorithm.background_color,
                    "color" => algorithm.color,
               )

          elseif typeof(algorithm) == FlatRenderer
              Dict(
                    "algorithm" => "Flat Renderer",
                    "background color" => algorithm.background_color,
               )
          elseif typeof(algorithm) == PathTracer
              Dict(
                    "algorithm" => "Path-Tracing Renderer",
                    "background color" => algorithm.background_color,
                    "PCG" => Dict(
                         "initial state" => get_state(algorithm.pcg),
                         "initial sequence" => get_inc(algorithm.pcg),
                         ),
                    "number of rays" => algorithm.num_of_rays,
                    "max depth" => algorithm.max_depth,
                    "russian roulette limit" => algorithm.russian_roulette_limit,
               )
          elseif typeof(algorithm) == PointLightRenderer
               Dict(
                    "algorithm" => "Point-Light Renderer",
                    "background color" => algorithm.background_color,
                    "color" => algorithm.ambient_color,
                    "dark parameter" => algorithm.dark_parameter,
               )
          end

     tonemapping = Dict(
          "normalization" => a,
          "gamma"=> γ,
          "input average luminosity" => lum,
     )

     data = Dict(
          "scene file" => scenefile,
          "time of start" => time_of_start,
          "png output" => png_output,
          "pfm output" => pfm_output,
          "camera" => dict_camera,
          "renderer" => dict_renderer,
          "alpha" => α,
          "width" => width,
          "height" => height,
          "tonemapping" => tonemapping,
          "samples per pixel (0 means no antialiasing)" => samples_per_side^2,
          "declared float from Command Line" => declare_float,
          "rendering time (in s)"=> @sprintf("%.3f", rendering_time_s),
     )

     open( join(map(x->x*".", split(png_output,".")[1:end-1])) * "json","w") do f
          JSON.print(f, data, 4)
     end
end

##########################################################################################92

function render(x::(Pair{T1,T2} where {T1,T2})...)
	render( parse_render_settings(  Dict( pair for pair in [x...]) )... )
end

function render(
          scenefile::String,
          renderer::Renderer = FlatRenderer(),
     	camera_type::Union{String, Nothing} = nothing,
		camera_position::Union{Point, Vec, Nothing} = nothing,
     	α::Float64 = 0.,
     	width::Int64 = 640,
     	height::Int64 = 480,
          a::Float64 = 0.18,
          γ::Float64 = 1.27,
          lum::Union{Number, Nothing} = nothing,
     	pfm_output::String = "scene.pfm",
        	png_output::String = "scene.png",
          samples_per_pixel::Int64 = 0,
		bool_print::Bool = true,
		bool_savepfm::Bool = true,
          declare_float::Union{Dict{String,Float64}, Nothing} = nothing,
          ONLY_FOR_TESTS::Bool = false,
     )

     (bool_print==true) && println("\n\nStarting the image rendering of \"$(scenefile)\"...")

     (ONLY_FOR_TESTS==false) || (return nothing)

     time_of_start = Dates.format(now(), DateFormat("Y-m-d : H:M:S"))
     time_1 = time()

     scene = open(scenefile, "r") do stream
          if isnothing(declare_float)
               inputstream = InputStream(stream, scenefile)
               parse_scene(inputstream)
          else
               inputstream = InputStream(stream, scenefile)
               parse_scene(inputstream, declare_float)
          end
     end

     (bool_print==true) && println("\nReaded and parsed \"$(scenefile)\", now initialize camera and renderer...\n")

     renderer.world = scene.world

     samples_per_side = string2rootint64(string(samples_per_pixel))

     observer_vec = isnothing(camera_position) ?
          nothing :
          typeof(camera_position) == Point ?
		camera_position - Point(0., 0., 0.) :
		camera_position

     aspect_ratio = width / height

     if isnothing(camera_type) && isnothing(observer_vec) && isnothing(scene.camera)
          camera = PerspectiveCamera(-1.0, aspect_ratio, rotation_z(deg2rad(α)))

     elseif isnothing(camera_type) && isnothing(observer_vec)
          camera = scene.camera
          camera.a = aspect_ratio

     elseif isnothing(camera_type) && isnothing(scene.camera)
          camera_tr = rotation_z(deg2rad(α)) * translation(observer_vec)
          camera = PerspectiveCamera(-1.0, aspect_ratio, camera_tr)

     elseif isnothing(observer_vec) && isnothing(scene.camera)
          if camera_type == "per"
		     (bool_print==true) && (println("Choosen perspective camera..."))
		     camera = PerspectiveCamera(1.0, aspect_ratio, rotation_z(deg2rad(α)))
	     elseif camera_type == "ort"
		     (bool_print==true) && (println("Choosen orthogonal camera..."))
		     camera = OrthogonalCamera(aspect_ratio, rotation_z(deg2rad(α)))
	     else
		     throw(ArgumentError("Unknown camera: \"$(camera_type)\""))
	     end

     elseif isnothing(camera_type)
          camera_tr = rotation_z(deg2rad(α)) * translation(observer_vec)
          if typeof(scene.camera) == OrthogonalCamera
               (bool_print==true) && (println("Choosen perspective camera..."))
               camera = OrthogonalCamera(aspect_ratio, camera_tr)
          elseif typeof(scene.camera) == PerspectiveCamera
               (bool_print==true) && (println("Choosen orthogonal camera..."))
               camera = PerspectiveCamera(scene.camera.d, aspect_ratio, camera_tr)
          else
		     throw(ArgumentError("Unknown camera: \"$(camera_type)\""))
	     end

     elseif isnothing(observer_vec)
          if camera_type == "per"
		     (bool_print==true) && (println("Choosen perspective camera..."))
		     camera = PerspectiveCamera(scene.camera.d, aspect_ratio, rotation_z(deg2rad(α)) * scene.camera.T)
	     elseif camera_type == "ort"
		     (bool_print==true) && (println("Choosen orthogonal camera..."))
		     camera = OrthogonalCamera(aspect_ratio, rotation_z(deg2rad(α)) * scene.camera.T)
	     else
		     throw(ArgumentError("Unknown camera: \"$(camera_type)\""))
	     end

     elseif isnothing(scene.camera)
          camera_tr = rotation_z(deg2rad(α)) * translation(observer_vec)
          if camera_type == "per"
		     (bool_print==true) && (println("Choosen perspective camera..."))
		     camera = PerspectiveCamera(1.0, aspect_ratio, camera_tr)
	     elseif camera_type == "ort"
		     (bool_print==true) && (println("Choosen orthogonal camera..."))
		     camera = OrthogonalCamera(aspect_ratio, camera_tr)
	     else
		     throw(ArgumentError("Unknown camera: \"$(camera_type)\""))
	     end

     else
          camera_tr = rotation_z(deg2rad(α)) * translation(observer_vec)
          if camera_type == "per"
		     (bool_print==true) && (println("Choosen perspective camera..."))
		     camera = PerspectiveCamera(scene.camera.d, aspect_ratio, camera_tr)
	     elseif camera_type == "ort"
		     (bool_print==true) && (println("Choosen orthogonal camera..."))
		     camera = OrthogonalCamera(aspect_ratio, camera_tr)
	     else
		     throw(ArgumentError("Unknown camera: \"$(camera_type)\""))
	     end

     end

     if typeof(renderer) == OnOffRenderer
		(bool_print==true) && (println("Choosen on-off renderer..."))
	elseif typeof(renderer) == FlatRenderer
		(bool_print==true) && (println("Choosen flat renderer..."))
	elseif typeof(renderer) == PathTracer
		(bool_print==true) && (println("Choosen path-tracing renderer..."))
	elseif typeof(renderer) == PointLightRenderer
          (bool_print==true) && (println("Choosen point-light renderer..."))
	else
		throw(ArgumentError("Unknown renderer: \"$(typeof(renderer))\""))
	end

	image = HDRimage(width, height)
	tracer = ImageTracer(image, camera, samples_per_side)

     algorithm = copy(renderer)

     (bool_print==true) && println("\nNow starts the rendering!\n")

     if (bool_print==true)
	     fire_all_rays!(tracer, renderer, (r,c) -> print_progress(r,c,image.height, image.width))
     else
          fire_all_rays!(tracer, renderer)
     end

     (bool_print==true) && println("\nRendering completed! Now saving the files...")

     img = tracer.img

	(bool_savepfm==true) && (open(pfm_output, "w") do outf; write(outf, img); end)
	(bool_print==true) && (println("\nHDR demo image written to \"$(pfm_output)\"!"))

     normalize_image!(img, a, lum)
	clamp_image!(img)
	γ_correction!(img, γ)

	if (typeof(query(png_output)) == File{DataFormat{:UNKNOWN}, String})
		(bool_print==true) && (
			println(
				"File{DataFormat{:UNKNOWN}, String} for \"$(png_output)\"n"*
				"Written as a .png file.\n"
			)
		)
     	Images.save(File{format"PNG"}(png_output), get_matrix(img))
	else
		Images.save(png_output, get_matrix(img))
	end

	(bool_print==true) && (println("\nHDR demo image written to \"$(png_output)\"!"))

     time_2 = time()
     rendering_time_s = time_2 - time_1

     (bool_savepfm==true) && print_JSON_render(
          png_output,
          pfm_output,
          scenefile,
          time_of_start,
          algorithm,
          camera,
          α,
     	width,
     	height,
          samples_per_side,
          declare_float,
          rendering_time_s,
          a, γ, lum,
     )

     name_json = join(map(x->x*".", split(png_output,".")[1:end-1])) * "json"
     (bool_print==true) && println("\nJSON file \"$(name_json)\" correctly created.")
     (bool_print==true) && println("\nEND OF RENDERING\n")
end

"""
	render(
          scenefile::String,
          renderer::Renderer = FlatRenderer(),
     	camera_type::Union{String, Nothing} = nothing,
		camera_position::Union{Point, Vec, Nothing} = nothing,
     	α::Float64 = 0.,
     	width::Int64 = 640,
     	height::Int64 = 480,
          a::Float64 = 0.18,
          γ::Float64 = 1.27,
          lum::Union{Number, Nothing} = nothing,
     	pfm_output::String = "scene.pfm",
        	png_output::String = "scene.png",
          samples_per_pixel::Int64 = 0,
		bool_print::Bool = true,
		bool_savepfm::Bool = true,
          declare_float::Union{Dict{String,Float64}, Nothing} = nothing,
          ONLY_FOR_TESTS::Bool = false,
          )

	render(x::(Pair{T1,T2} where {T1,T2})...) =
	     render( parse_render_settings(  Dict( pair for pair in [x...]) )... )

Render the input `scenefile` with the specified options, and creates the following
three files:
- the PFM image (`scene.pfm` is the default name, if none is specified from the command line)
- the LDR image (`scene.png` is the default name, if none is specified from the command line)
- the JSON file (which has the same name of the LDR image and `.json` estention, so
  `scene.json` is the default name, if none LDR image name is specified from the command line),
  that saves some datas about input commands, rendering time etc.

## Arguments

The following input arguments refers to the first method presented in the signature;
it's obviously very uncomfortable to use that method, so it's recommended to take
advantage of the second one, which allows to write the input values in a dictionary
like syntax with arbitrary order and comfort. See the documentation of
[`parse_render_settings`](@ref) to learn how to use the keys:

- `renderer::Renderer = FlatRenderer()` : renderer to be used in the rendering, with all
  the settings already setted (exception made for the `world`, that will be overridden
  and created here)

- `camera_type::String = "per"` : set the perspective projection view:
  - `camera_type=="per"` -> set [`PerspectiveCamera`](@ref)  (default value)
  - `camera_type=="ort"`  -> set [`OrthogonalCamera`](@ref)

- `camera_position::Union{Point, Vec} = Point(-1.,0.,0.)` : set the point of observation
  in (`X`,`Y,`Z`) coordinates, as a `Point`or a `Vec` input object.

- `α::Float64 = 0.` : angle of rotation _*IN RADIANTS*_, relative to the vertical
  (i.e. z) axis with a right-handed rule convention (clockwise rotation for entering (x,y,z)-axis
  corresponds to a positive input rotation angle)

- `width::Int64 = 640` and `height::Int64 = 480` : pixel dimensions of the demo image;
  they must be both even positive integers.

- `a::Float64 = 0.18` : normalization scale factor for the tone mapping.

- `γ::Float64 = 1.27` : gamma factor for the tone mapping.

- `lum::Union{Number, Nothing} = nothing ` : average luminosity of the image; iIf not specified or equal to 0,
  it's calculated through [`avg_lum`](@ref)

- `pfm_output::String = "demo.pfm"` : name of the output pfm file

- `png_output::String = "demo.png"` : name of the output LDR file

- `samples_per_pixel::Int64 = 0` : number of rays per pixel to be used (antialiasing);
  it must be a perfect square positive integer (0, 1, 4, 9, ...) and if is set to
  0 (default value) is choosen, no anti-aliasing occurs, and only one pixel-centered
  ray is fired for each pixel.

- `bool_print::Bool = true` : specifies if the WIP messages of the demo
  function should be printed or not (useful option for [`render_animation`](@ref))

- `bool_savepfm::Bool = true` : bool that specifies if the pfm file should be saved
  or not (useful option for [`render_animation`](@ref))

- `declare_float::Union{Dict{String,Float64}, Nothing} = nothing` : an option (for the
  command line in particularly) to manually override the values of the float variables in
  the scene file; each overriden variable name (the key) is associated with its float value
  (i.e. `declare_float = Dict("var1"=>0.1, "var2"=>2.5)`)

- `ONLY_FOR_TESTS::Bool = false` : it's a bool variable conceived only to
  test the correct behaviour of the renderer for the input arguments; if set to `true`,
  no rendering is made!

See also: [`Point`](@ref), [`Vec`](@ref), [`Renderer`](@ref) [`OnOffRenderer`](@ref),
[`FlatRenderer`](@ref), [`PathTracer`](@ref), [`PointLightRenderer`](@ref),
[`parse_render_settings`](@ref)
"""
render


text:
using Oceananigans.Operators: index_and_interp_dependencies
using Oceananigans.Utils: tupleit, user_function_arguments

import Oceananigans: location

"""
    struct ContinuousBoundaryFunction{X, Y, Z, I, F, P, D, N, ℑ} <: Function

A wrapper for the user-defined boundary condition function `func` at location
`X, Y, Z`. `I` denotes the boundary-normal index (`I=1` at western boundaries,
`I=grid.Nx` at eastern boundaries, etc). `F, P, D, N, ℑ` are, respectively, the
user-defined function, parameters, field dependencies, indices of the field dependencies
in `model_fields`, and interpolation operators for interpolating `model_fields` to the
location at which the boundary condition is applied.
"""
struct ContinuousBoundaryFunction{X, Y, Z, I, F, P, D, N, ℑ} <: Function
                          func :: F
                    parameters :: P
            field_dependencies :: D
    field_dependencies_indices :: N
     field_dependencies_interp :: ℑ

    """ Returns a location-less wrapper for `func`, `parameters`, and `field_dependencies`."""
    function ContinuousBoundaryFunction(func::F, parameters::P, field_dependencies) where {F, P}
        field_dependencies = tupleit(field_dependencies)
        D = typeof(field_dependencies)
        return new{Nothing, Nothing, Nothing, Nothing, F, P, D, Nothing, Nothing}(func, parameters, field_dependencies, nothing, nothing)
    end

    function ContinuousBoundaryFunction{X, Y, Z, I}(func::F,
                                                    parameters::P,
                                                    field_dependencies::D,
                                                    field_dependencies_indices::N,
                                                    field_dependencies_interp::ℑ) where {X, Y, Z, I, F, P, D, ℑ, N}

        return new{X, Y, Z, I, F, P, D, N, ℑ}(func, parameters, field_dependencies, field_dependencies_indices, field_dependencies_interp)
    end
end

location(::ContinuousBoundaryFunction{X, Y, Z}) where {X, Y, Z} = X, Y, Z

#####
##### "Regularization" for IncompressibleModel setup
#####

"""
    regularize_boundary_condition(bc::BoundaryCondition{C, <:ContinuousBoundaryFunction},
                                  topo, loc, dim, I, prognostic_field_names) where C

Regularizes `bc.condition` for location `loc`, boundary index `I`, and `prognostic_field_names`,
returning `BoundaryCondition(C, regularized_condition)`.

The regularization of `bc.condition::ContinuousBoundaryFunction` requries

1. Setting the boundary location to `LX, LY, LZ`.
   The location in the boundary-normal direction is `Nothing`.

2. Setting the boundary-normal index `I` for indexing into `field_dependencies`.
   `I` is either `1` (for left boundaries) or
   `size(grid, n)` for a boundary in the `n`th direction where `n ∈ (1, 2, 3)` corresponds
   to `x, y, z`.

3. Determining the `indices` that map `model_fields` to `field_dependencies`.

4. Determining the `interps` functions that interpolate field_dependencies to the location
   of the boundary.
"""
function regularize_boundary_condition(bc::BoundaryCondition{C, <:ContinuousBoundaryFunction},
                                       topo, loc, dim, I, prognostic_field_names) where C

    boundary_func = bc.condition

    # Set boundary-normal location to Nothing:
    LX, LY, LZ = Tuple(i == dim ? Nothing : loc[i] for i = 1:3)

    indices, interps = index_and_interp_dependencies(LX, LY, LZ,
                                                     boundary_func.field_dependencies,
                                                     prognostic_field_names)

    regularized_boundary_func = ContinuousBoundaryFunction{LX, LY, LZ, I}(boundary_func.func,
                                                                          boundary_func.parameters,
                                                                          boundary_func.field_dependencies,
                                                                          indices, interps)

    return BoundaryCondition(C, regularized_boundary_func)
end

@inline boundary_index(i) = ifelse(i == 1, 1, i + 1) # convert near-boundary Center() index into boundary Face() index

#####
##### Kernel functions
#####

# Return ContinuousBoundaryFunction on east or west boundaries.
@inline function (bc::ContinuousBoundaryFunction{Nothing, LY, LZ, i})(j, k, grid, clock, model_fields) where {LY, LZ, i}
    args = user_function_arguments(i, j, k, grid, model_fields, bc.parameters, bc)

    i′ = boundary_index(i)

    return bc.func(ynode(Face(), LY(), LZ(), i′, j, k, grid),
                   znode(Face(), LY(), LZ(), i′, j, k, grid),
                   clock.time, args...)
end

# Return ContinuousBoundaryFunction on south or north boundaries.
@inline function (bc::ContinuousBoundaryFunction{LX, Nothing, LZ, j})(i, k, grid, clock, model_fields) where {LX, LZ, j}
    args = user_function_arguments(i, j, k, grid, model_fields, bc.parameters, bc)

    j′ = boundary_index(j)

    return bc.func(xnode(LX(), Face(), LZ(), i, j′, k, grid),
                   znode(LX(), Face(), LZ(), i, j′, k, grid),
                   clock.time, args...)
end

# Return ContinuousBoundaryFunction on bottom or top boundaries.
@inline function (bc::ContinuousBoundaryFunction{LX, LY, Nothing, k})(i, j, grid, clock, model_fields) where {LX, LY, k}
    args = user_function_arguments(i, j, k, grid, model_fields, bc.parameters, bc)

    k′ = boundary_index(k)

    return bc.func(xnode(LX(), LY(), Face(), i, j, k′, grid),
                   ynode(LY(), LY(), Face(), i, j, k′, grid),
                   clock.time, args...)
end

# Don't re-convert ContinuousBoundaryFunctions passed to BoundaryCondition constructor
BoundaryCondition(Classification::DataType, condition::ContinuousBoundaryFunction) = BoundaryCondition(Classification(), condition)

Adapt.adapt_structure(to, bf::ContinuousBoundaryFunction{LX, LY, LZ, I}) where {LX, LY, LZ, I} =
    ContinuousBoundaryFunction{LX, LY, LZ, I}(Adapt.adapt(to, bf.func),
                                              Adapt.adapt(to, bf.parameters),
                                              nothing,
                                              Adapt.adapt(to, bf.field_dependencies_indices),
                                              Adapt.adapt(to, bf.field_dependencies_interp))


text:
using LinearAlgebra, LaTeXStrings, Plots; pyplot()

L = 10
p0, p1 = 1/2, 3/4
beta = 0.75
epsilon = 0.001

function valueIteration(kappa)
    P0 = diagm(1=>fill(p0,L-1)) + diagm(-1=>fill(1-p0,L-1))
    P0[1,1], P0[L,L] = 1 - p0, p0

    P1 = diagm(1=>fill(p1,L-1)) + diagm(-1=>fill(1-p1,L-1))
    P1[1,1], P1[L,L] = 1 - p1, p1

    R0 = collect(1:L)
    R1 = R0 .- kappa

    bellmanOperator(Vprev)=
        max.(R0 + beta*P0*Vprev, R1 + beta*P1*Vprev)
    optimalPolicy(V,state)=
        (R0+beta*P0*V)[state] >= (R1+beta*P1*V)[state] ? 0 : 1

    V, Vprev = fill(0,L), fill(1,L)
    while norm(V-Vprev) > epsilon
        Vprev = V
        V = bellmanOperator(Vprev)
    end

    return [optimalPolicy(V,s) for s in 1:L]
end

kappaGrid = 0:0.1:2.0
policyMap = zeros(L,length(kappaGrid))

for (i,kappa) in enumerate(kappaGrid)
    policyMap[:,i] = valueIteration(kappa)
end
heatmap(policyMap, fill=cgrad([:blue, :red]),
        xticks=(0:1:21, -0.1:0.1:2), yticks=(0:L, 0:L),
        xlabel=L"\kappa", ylabel="State", colorbar_entry=false)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
"""
Tests for NBP dynamics.
"""
module NBPTests

    using AstrodynamicalModels, ModelingToolkit, Test

    @testset "NBP Model Constructors" begin
        NBP(10; stm=false, structural_simplify=false)
        NBP(2; stm=true, structural_simplify=true)
        @test true
    end

    @testset "NBP Function Constructors" begin
        NBPFunction(10; stm=false, structural_simplify=false)
        NBPFunction(2; stm=true, structural_simplify=true)
        @test true
    end

end


text:
#################### Phylogenetic No-U-Turn Sampler ####################

#################### Types and Constructors ####################

mutable struct PNUTSTune <: SamplerTune
  logfgrad::Union{Function, Missing}
  adapt::Bool
  alpha::Float64
  epsilon::Float64
  epsilonbar::Float64
  gamma::Float64
  Hbar::Float64
  kappa::Float64
  m::Int
  mu::Float64
  nalpha::Int
  t0::Float64
  delta::Float64
  target::Float64
  moves::Int
  tree_depth::Int
  nniprime::Int
  targetNNI::Int

  PNUTSTune() = new()

  function PNUTSTune(x::Vector{T}, epsilon::Float64, logfgrad::Union{Function, Missing};
                    target::Real=0.6, tree_depth::Int=10, targetNNI::Int=5) where T<:Node
    new(logfgrad, false, 0.0, epsilon, 1.0, 0.05, 0.0, 0.75, 0, NaN, 0, 10.0,0.003,
        target,0, tree_depth,0, targetNNI)
  end
end

PNUTSTune(x::Vector{T}, logfgrad::Function, ::NullFunction, delta::Float64=0.003; args...) where T<:Node =
  PNUTSTune(x, nutsepsilon(x[1], logfgrad, delta), logfgrad; args...)

PNUTSTune(x::Vector{T}, logfgrad::Function, delta::Float64; args...) where T<:Node =
  PNUTSTune(x, nutsepsilon(x[1], logfgrad, delta), logfgrad; args...)

PNUTSTune(x::Vector; epsilon::Real, args...) =
    NUTSTune(x, epsilon, missing, args...)

const PNUTSVariate = SamplerVariate{PNUTSTune}

#################### Sampler Constructor ####################

"""
    PNUTS(params::ElementOrVector{Symbol}; args...)

Construct a `Sampler` object for PNUTS sampling. The Parameter is assumed to be
a tree.

Returns a `Sampler{PNUTSTune}` type object.

* params: stochastic node to be updated with the sampler.

* args...: additional keyword arguments to be passed to the PNUTSVariate constructor.
"""
function PNUTS(params::ElementOrVector{Symbol}; args...)
  samplerfx = function(model::Model, block::Integer)
    block = SamplingBlock(model, block, true)
    f = (x, sz, ll, gr) -> mlogpdfgrad!(block, x, sz, ll, gr)
    v = SamplerVariate(block, f, NullFunction(); args...)

    sample!(v::PNUTSVariate, f, adapt=model.iter <= model.burnin)

    relist(block, v)
  end
  Sampler(params, samplerfx, PNUTSTune())
end

#################### Sampling Functions ####################

function mlogpdfgrad!(block::SamplingBlock, x::T, sz::Int64, ll::Bool=false, gr::Bool=false)::Tuple{Float64, Vector{Float64}}  where T<:Node
  grad = Vector{Float64}(undef, sz)
  lp = zero(Float64)

  if gr
    lp, grad = gradlogpdf!(block, x)::Tuple{Float64, Vector{Float64}}
  elseif ll
    lp = logpdf!(block, x)::Float64
  end
  lp, grad
end
sample!(v::PNUTSVariate; args...) = sample!(v, v.tune.logfgrad; args...)

function sample!(v::PNUTSVariate, logfgrad::Function; adapt::Bool=false)

  tune = v.tune
  setadapt!(v, adapt)
  if tune.adapt
     tune.m += 1
     tune.nniprime = 0

     nuts_sub!(v, tune.epsilon, logfgrad)
     Ht = (tune.target - tune.alpha / tune.nalpha)
     avgnni = tune.targetNNI - tune.nniprime/tune.nalpha
     HT2 = - avgnni / (1 + abs(avgnni))
     HT = (Ht + HT2)/2
     p = 1.0 / (tune.m + tune.t0)
     tune.Hbar = (1.0 - p) * tune.Hbar +
                 p * HT
     tune.epsilon = exp(tune.mu - sqrt(tune.m) * tune.Hbar / tune.gamma)
     p = tune.m^-tune.kappa
     tune.epsilonbar = exp(p * log(tune.epsilon) +
                           (1.0 - p) * log(tune.epsilonbar))
  else
     if (tune.m > 0) tune.epsilon = tune.epsilonbar end

     nuts_sub!(v, tune.epsilon, logfgrad)
  end
  v
end

function setadapt!(v::PNUTSVariate, adapt::Bool)
  tune = v.tune
  if adapt && !tune.adapt
    tune.m = 0
    tune.mu = log(10.0 * tune.epsilon)
  end
  tune.adapt = adapt
  v
end

function nuts_sub!(v::PNUTSVariate, epsilon::Float64, logfgrad::Function)

  mt = v.value[1]
  nl = size(mt)[1]-1
  delta = v.tune.delta

  r = randn(nl)
  g = zeros(nl)

  x, r , logf, grad, nni = refraction(mt, r, 1, g, epsilon, logfgrad, delta, nl)

  lu = log(rand())

  logp0 = logf - 0.5 * dot(r)
  logu0 = logp0 + lu
  rminus = rplus = r
  gradminus = gradplus = grad

  xminus = xplus = x

  j = 0
  n = 1
  s = true

  while s && j < v.tune.tree_depth

    pm =2 * (rand() > 0.5) - 1

    if pm == -1

      xminus, rminus, gradminus, _, _, _, xprime, nprime, sprime, alpha,
        nalpha, nni1, lpp, nniprime = buildtree(xminus, rminus, gradminus, pm, j, epsilon, logfgrad,
                           logp0, logu0, delta, nl,lu)

    else

      _, _, _, xplus, rplus, gradplus, xprime, nprime, sprime, alpha, nalpha, nni1, lpp, nniprime =
        buildtree(xplus, rplus, gradplus, pm, j, epsilon, logfgrad, logp0,
                  logu0, delta, nl, lu)

    end#if pm

    if sprime && rand() < nprime / n
        v.value[1]= xprime
    end
    j += 1
    n += nprime
    s = sprime && nouturn(xminus, xplus, rminus, rplus, gradminus, gradplus, epsilon, logfgrad, delta, nl, j)
    v.tune.alpha, v.tune.nalpha, v.tune.nniprime = alpha, nalpha, nniprime
    nni += nni1
  end
  #println(j)
  v.tune.moves += nni
  v
end

function refraction(v::T, r::Vector{Float64}, pm::Int64,
                    grad::Vector{Float64}, epsilon::Float64, logfgrad::Function,
                    delta::Float64, sz::Int64)  where T<:Node

    v1 = deepcopy(v)

    ref_r = pm*r

    blenvec = get_branchlength_vector(v1)
    fac = scale_fac.(blenvec, delta)

    ref_r = @. ref_r + (epsilon * 0.5) * (grad * fac)
    tmpB = @. blenvec + (epsilon * ref_r)

    nni = 0

    if minimum(tmpB) <= 0
        v1, tmpB, ref_r, nni = ref_NNI(v1, tmpB, ref_r, epsilon, blenvec, delta, logfgrad, sz)

    end

    blenvec = molifier.(tmpB, delta)

    set_branchlength_vector!(v1, blenvec)

    logf, grad = logfgrad(v1, sz, true, true)

    fac = scale_fac.(blenvec, delta)
    ref_r = @. ref_r + (epsilon * 0.5) * (grad * fac)

    r = pm*ref_r

    return v1, r, logf, grad, nni
end

function ref_NNI(v::T, tmpB::Vector{Float64}, r::Vector{Float64}, epsilon::Float64, blv::Vector{Float64},
                 delta::Float64, logfgrad::Function, sz::Int64)  where T<:Node

  intext = internal_external(v)
  t = 0.0
  nni = 0

  while minimum(tmpB)<=0.0

     timelist = tmpB./abs.(r)
     ref_index = argmin(timelist)

     temp = epsilon-t+timelist[ref_index]
     blv = @. blv + (temp * r)

     r[ref_index] *= -1.0

     if intext[ref_index] == 1

       blv1 =molifier.(blv, delta)
       set_branchlength_vector!(v, blv1)

       # use thread parallelism
       res_before = @spawn logfgrad(v, sz, true, false) # still with molified branch length
       #U_before_nni = logfgrad(v, sz, true, false) # still with molified branch length
       v_copy = deepcopy(v)
       tmp_NNI_made = NNI!(v_copy, ref_index)

       # fetch the results from the parallel part
       U_before_nni, _ = fetch(res_before)
       U_before_nni *= -1

       if tmp_NNI_made != 0

            U_after_nni, _ = logfgrad(v_copy, sz, true, false)
            U_after_nni *= -1
            delta_U::Float64 = 2.0*(U_after_nni - U_before_nni)
            my_v::Float64 = r[ref_index]^2
            if my_v > delta_U
              nni += tmp_NNI_made
              r[ref_index] = sqrt(my_v - delta_U)
              v = v_copy
            end # if my_v
       end #if NNI
      end #non leave
      t = epsilon+timelist[ref_index]

      tmpB = @. blv + (epsilon-t) * r

  end #while

  v, tmpB, r, nni
end

function buildtree(x::T, r::Vector{Float64},
                   grad::Vector{Float64}, pm::Int64, j::Integer,
                   epsilon::Float64, logfgrad::Function, logp0::Real, logu0::Real,
                   delta::Float64, sz::Int64, lu::Float64)  where T<:Node

  if j == 0

    xprime, rprime, logfprime, gradprime, nni = refraction(x, r, pm, grad, epsilon,
                                          logfgrad, delta, sz)

    logpprime = logfprime - 0.5 * dot(rprime)

    nprime = Int(logu0 < logpprime)

    sprime = logu0 < logpprime + 1000.0
    xminus = xplus = xprime
    rminus = rplus = rprime
    gradminus = gradplus = gradprime
    alphaprime = min(1.0, exp(logpprime - logp0))
    nniprime = nni
    nalphaprime = 1

  else
    xminus, rminus, gradminus, xplus, rplus, gradplus, xprime, nprime, sprime,
      alphaprime, nalphaprime, nni, logpprime, nniprime = buildtree(x, r, grad, pm, j - 1, epsilon,
                                          logfgrad, logp0, logu0, delta,  sz, lu)
    if sprime

      if pm == -1

        xminus, rminus, gradminus, _, _, _, xprime2, nprime2, sprime2,
          alphaprime2, nalphaprime2 , nni, logpprime, nniprime2 = buildtree(xminus, rminus, gradminus, pm,
                                                j - 1, epsilon, logfgrad, logp0,
                                                logu0, delta,  sz, lu)
      else
        _, _, _, xplus, rplus, gradplus, xprime2, nprime2, sprime2,
          alphaprime2, nalphaprime2, nni, logpprime, nniprime2 = buildtree(xplus, rplus, gradplus, pm,
                                                j - 1, epsilon, logfgrad, logp0,
                                                logu0, delta, sz,lu)
      end # if pm

      if rand() < nprime2 / (nprime + nprime2)
        xprime = xprime2
      end
      nprime += nprime2
      sprime = sprime2 && nouturn(xminus, xplus, rminus, rplus, gradminus, gradplus, epsilon, logfgrad, delta, sz, j)
      alphaprime += alphaprime2
      nalphaprime += nalphaprime2
      nniprime += nniprime2
    end #if sprime
  end #if j

  xminus, rminus, gradminus, xplus, rplus, gradplus, xprime, nprime, sprime,
    alphaprime, nalphaprime, nni, logpprime, nniprime
end

function nouturn(xminus::T, xplus::T,
                rminus::Vector{Float64}, rplus::Vector{Float64}, gradminus::Vector{Float64},gradplus::Vector{Float64},
                epsilon::Float64, logfgrad::Function, delta::Float64, sz::Int64, j::Int64)  where T<:Node

        curr_l, curr_h = BHV_bounds(xminus, xplus)

        # use thread parallelism to calculuate both directions at once
        res_minus = Base.Threads.@spawn refraction(deepcopy(xminus), deepcopy(rminus), -1, gradminus, epsilon, logfgrad, delta, sz)
        xplus_bar,_,_,_,_ = refraction(deepcopy(xplus), deepcopy(rplus),1, gradplus, epsilon, logfgrad, delta,sz)

        # fetch the results
        xminus_bar,_,_,_,_ = fetch(res_minus)

        curr_t_l, curr_t_h = BHV_bounds(xminus_bar, xplus_bar)
        return curr_h < curr_t_l
end

#################### Auxilliary Functions ####################

function nutsepsilon(x::Node, logfgrad::Function, delta::Float64)

  x0 = deepcopy(x)
  n = size(x)[1] - 1

  _, r0, logf0, grad0,_ = refraction(x0, randn(n), 1, zeros(n), 0.0, logfgrad, delta, n)

  x0 = deepcopy(x)
  epsilon = 1.0
  _, rprime, logfprime, gradprime,_ = refraction(x0, r0, 1, grad0,  epsilon, logfgrad, delta, n)

  prob = exp(logfprime - logf0 - 0.5 * (dot(rprime) - dot(r0)))

  pm = 2 * (prob > 0.5) - 1
  while prob^pm > 0.5^pm
    epsilon *= 2.0^pm
    x0 = deepcopy(x)
    _, rprime, logfprime, _ ,_ = refraction(x0, r0, 1, grad0, epsilon, logfgrad, delta, n)

    prob = exp(logfprime - logf0 - 0.5 * (dot(rprime) - dot(r0)))

  end

  epsilon
end

@inline function scale_fac(x::T, delta::T) where T<:Float64
  x < delta ? x/delta : 1.0
end

@inline function molifier(x::Float64, delta::Float64)::Float64
    x >= delta ? x : (x^2+delta^2)/(2.0*delta)
end # function


text:
######## Part 1 Useful Information Retrieval ###################

#=
F: extracting BioSym & Verb of each sentence;
I: each sentence in Array of tuples, sentenceType dictionary;
O: sentence type, & array of biosymbol groups;
BA
=#
function extractBioSymGroups(tokenTagPairs::Array, keyVerbs::Set)
  BioSymGroupsArray = []
  typeMarker = ""
  parameterSetting = []
  i = 1
  while i <= length(tokenTagPairs)
    if tokenTagPairs[i][2] == "BioSym"  # find a BioSym group
      tmpBioSymGroup = [tokenTagPairs[i][1]]
      # repeat until find the whole group
      while (i+1 <= length(tokenTagPairs)) && (tokenTagPairs[i+1][2] == "BioSym")
        i += 1  # go for next one
        push!(tmpBioSymGroup, tokenTagPairs[i][1])
      end
      push!(BioSymGroupsArray, tmpBioSymGroup)
      # i += 1  # go for next one
    elseif in(tokenTagPairs[i][2], keyVerbs)
      typeMarker = tokenTagPairs[i][2]
      # i += 1  # go for next one
    elseif tokenTagPairs[i][2] == "reversible"
      push!(parameterSetting, "reversible")
      # i += 1  # go for next one
    elseif tokenTagPairs[i][2] == "irreversible"
      push!(parameterSetting, "irreversible")
    end
    i += 1  # go for next one
  end
  if typeMarker == ""
    println("ERROR: SENTENCE TYPE UNDEFINED for: $(tokenTagPairs)")
  end
  return typeMarker, BioSymGroupsArray, parameterSetting
end

#=
F: extracting BioSym & Key Verb information from tagged sentences, call
   extractBioSymGroups() to process each sentence;
I: Array of sentence, with each sentence in form of array of tuples (token, tag);
O: Array of Tuples, with each tuple in form of (sentence type, array of biosymbol groups);
BA
=#
function extractBioSymVerbInformation(taggedSenteces::Array,
                                      sentenceTypeKeyVerbs::Set, error::Dict)
  BioSymGroupsKeyVerbArray = []
  for i in 1:length(taggedSenteces)
    senType, BioSym, paraSet = extractBioSymGroups(taggedSenteces[i][1], sentenceTypeKeyVerbs)
    lineNum = taggedSenteces[i][2]
    if length(senType) == 0 # catch SENTENCE TYPE UNDEFINED error
        if haskey(error, lineNum)
            push!(error[lineNum], "ERROR: SENTENCE TYPE UNDEFINED")
        else
            error[lineNum] = ["ERROR: SENTENCE TYPE UNDEFINED"]
        end
    else  # successfully extract useful information
        if length(paraSet) == 0
          push!(BioSymGroupsKeyVerbArray, ((senType, lineNum), BioSym))
        else
          push!(BioSymGroupsKeyVerbArray, ((senType, lineNum), BioSym, paraSet))
        end
    end
  end
  return BioSymGroupsKeyVerbArray
end

########## Part 2 Symbol Decoding ###############################

#=
F: decoding biosymbol groups by enumerating all possible combinations
I: Array of tuples (Verb, biosymbol groups), "biosymbol groups" is array of
   biosymbol arrays.
O: biosymbol groups are replaced with corresponding array of biosymbol array of
   biosymbol arrays. (2-layer array for simple case, 3-layer array for symbol
   groups with parenthesis)
BA
=#
function decodingBioSymGroups(VerbBioSymArray::Array, error::Dict)
  # println(VerbBioSymArray)
  isID = Set{Int64}()
  for (id, vbsa) in enumerate(VerbBioSymArray)
    lineNum = vbsa[1][2]
    for i in 1:length(vbsa[2])
      errSym = join(vbsa[2][i], " ")
      vbsa[2][i], err = decodingABioSymGroupByLogicalRules(vbsa[2][i])
      if length(err) != 0
          push!(isID, id)
          if !haskey(error, lineNum)
              error[lineNum] = []
          end
          for e in err
              push!(error[lineNum], "in $(errSym), " * e)
          end
      end
    end
  end
  # delete this from array if has errors?
  deleteat!(VerbBioSymArray, isID)
end

#= ambiguity problem
F: just handling simple cases (one-layer of parenthesis),
   e.g. (A or B) & C & D & (E or F), (A & B) or C or D.
   --> currently, try to do the right thing, but does not consider all possible cases.
I: array of biological tokens (including logcial tokens & bio-tokens)
O: Array of Strings for simple cases; or array of arrays of strings for cases
   with parenthesis.
BA
=#
function decodingABioSymGroupByLogicalRules(tokens::Array)
  println("\n\ninside decodingABioSymGroupByLogicalRules")
  println(tokens)
  err_mes = []  # collection of error messages
  andSet = Set(["and", "&", ","])
  orSet = Set(["or", "|"])

  finalBioSymGroups = []  # elements in parallel relation
  BioSymGroup = []  # for insiders, each element is in the form of [logicalRelation, token1, token2, ...]

  tokensOutParenthesis = []  # for the outsider
  i = 1
  while i <= length(tokens)
    if tokens[i] == "("  # find a ()
      tokensInsideParenthesis = []  # for tokens inside a pair of ()
      i += 1
      while (i<=length(tokens)) && (tokens[i] != ")") # assume one layer of ()
        push!(tokensInsideParenthesis, tokens[i])
        i += 1
      end
      # out of while loop, i is at ")" if i <= length(tokens)
      if i <= length(tokens) && tokens[i] == ")" && length(tokensInsideParenthesis) != 0
        #  has tokens insider () -> handle it
        logicInParenthesis, tmpBioSymGroup, errMes = decodingASequenceOfTokensWithOneLogicalRelation(
        tokensInsideParenthesis, andSet, orSet)
        if errMes != ""  # error occurs
            push!(err_mes, errMes * "inside ()")
        else
          if logicInParenthesis == ""
              logicInParenthesis = "and"
          end
          prepend!(tmpBioSymGroup, [logicInParenthesis])
          push!(BioSymGroup, tmpBioSymGroup)
        end
      else  # what kind of errors?
          push!(err_mes, "ERRORS: incomplete or empty group inside \'()\'")
          println("ERRORS: incomplete or empty group inside \'()\'")
      end
    else  # the outsiders
      push!(tokensOutParenthesis, tokens[i])
    end
    i += 1  # move forward
  end
  # handling outsiders
  logicOutParen, outsiderGroup, errMes = decodingASequenceOfTokensWithOneLogicalRelation(
  tokensOutParenthesis, andSet, orSet)
  if errMes != ""  # error occurs
      push!(err_mes, errMes * "outside ()")
  else
    println("outsiderGroup: $outsiderGroup")
    println("logicOutParen: $logicOutParen")
    if logicOutParen == ""  # no relation for outsiders
      if length(outsiderGroup) == 0  # nothing outside ()
        if length(BioSymGroup) == 1 # can only have 1 group
          finalBioSymGroups = reshapeBioSymbolGroups([], BioSymGroup[1][2:end], BioSymGroup[1][1])
          # println("Simple Bio-sym group")
        else
          println("ERROR: unspecified/unclear relationship")
          push!(err_mes, "ERROR: unspecified/unclear relationship")
        end
      elseif length(outsiderGroup) == 1 && length(BioSymGroup) == 0 # only one token
        finalBioSymGroups = deepcopy(outsiderGroup)
        println("ONE token biosymbol")
      else
        println("ERROR: unspecified/unclear relationship")
        push!(err_mes, "ERROR: unspecified/unclear relationship")
      end
    else  # check logics in/outside () and reshape BioSymGroup for return
      goOnAllow =  true
      for (id, biosym) in enumerate(BioSymGroup)  # check logics in/outside ()
        if biosym[1] == logicOutParen  # logic check
          push!(err_mes, "ERROR: logic inside and outside \'()\' should be different")
          println("ERROR: logic inside and outside \'()\' should be different")
          goOnAllow =  false
        end
      end
      if goOnAllow  # reshape
        BioSymGroup = [subg[2:end] for subg in BioSymGroup]  # reshape
        finalBioSymGroups = reshapeBioSymbolGroups(BioSymGroup, outsiderGroup, logicOutParen)
        println(finalBioSymGroups)
      end
    end
    # println("\nfinalBioSymGroups: ")
    # println(s for s in finalBioSymGroups)
  end

  return finalBioSymGroups, err_mes
end

#=
F: take in a sequence of tokens, seperate bio tokens and logical tokens, return
   logical relation and bio symbol array. Assume only one relation in this sequence.
IOBA
=#
function decodingASequenceOfTokensWithOneLogicalRelation(tokens::Array, andSet::Set, orSet::Set)
  # println("inside decodingASequenceOfTokensWithOneLogicalRelation")
  biosymGroup = []
  logicalRelation = [""]
  errorMessage = ""
  for i = 1:length(tokens)
    if in(tokens[i], andSet)  # find "and"
      if logicalRelation[1] == ""
        logicalRelation[1] = "and"
      elseif logicalRelation[1] == "or"
        errorMessage = "ERROR: LOGICAL INCONSISTENCY"
        println("LOGICAL INCONSISTENCY")
        return logicalRelation[1], BioSymGroup, errorMessage
      else
        continue
      end
    elseif in(tokens[i], orSet)  # find "or"
      if logicalRelation[1] == ""
        logicalRelation[1] = "or"
      elseif logicalRelation[1] == "and"
          errorMessage = "ERROR: LOGICAL INCONSISTENCY"
          println("LOGICAL INCONSISTENCY")
          return logicalRelation[1], BioSymGroup, errorMessage
      else
        continue
      end
    else
      push!(biosymGroup, tokens[i])
    end
  end
  # println(logicalRelation[1])
  # println("jump out decodingASequenceOfTokensWithOneLogicalRelation")
  return logicalRelation[1], biosymGroup, errorMessage
end

#=
F: get all permutations from "OR"/"AND" relations
IOBA
=#
function reshapeBioSymbolGroups(inParenGroups, outParenGroup, logicOutParen)
  finalGroups = []
  if logicOutParen == "and"
    println("inParenGroups: $inParenGroups")
    println("outParenGroup: $outParenGroup")
    if length(inParenGroups) == 0  # simple case: Array of strings
      finalGroups = deepcopy(outParenGroup)
      println("finalGroups: $finalGroups")
      println("Logic AND simple bio-symbol group")
    else  # complicated case: array of arrays of strings
      finalGroups = [outParenGroup]
      # println("finalGroups: $finalGroups")
      for i = 1:length(inParenGroups)
        tmp2 = []
        for biosym in inParenGroups[i]  # append each token to the previous finalGroups
          tmp1 = deepcopy(finalGroups)
          # println("tmp1: $tmp1")
          for j = 1:length(tmp1)
            push!(tmp1[j], biosym)
          end
          # println("tmp1: $tmp1")
          append!(tmp2, tmp1)
        end
        # println("tmp2: $tmp2")
        finalGroups = deepcopy(tmp2)  # update finalGroups after processing each outParenGroup element
        # println("finalGroups: $finalGroups")
      end
      println("Logic AND complicated bio-symbol group")
    end
  else  # "or" --> split outsiders into arrays, add insiders
    for biosym in outParenGroup
      push!(finalGroups, [biosym])
    end
    if length(inParenGroups) == 0 # simple case
      # finalGroups = outParenGroup
      println("Logic OR simple bio-symbol group")
    else
      for biosymgroup in inParenGroups
        push!(finalGroups, biosymgroup)
      end
      println("Logic OR complicated bio-symbol group")
    end
  end
  return finalGroups
end

#=
F: print out an array of tuples, each tuple is a pair of a token and a multiple-
   layer arrays.
I: an array of tuples, each tuple is a pair of a token and a multiple-
   layer arrays. Array of arrays of strings (2-layers) for simple cases,
   array of arrays of arrays of strings (3-layers) for complicated cases.
O: see terminal.
BA
=#
function printArrayOfTripleToken(tokenPlusMultiLayerArray::Array, targetDataType::Type)
  println("\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
  for tuple in tokenPlusMultiLayerArray
    println(tuple[1])
    println()
    if length(tuple[2]) == 0  # empty biosymbols
      println("THIS sentence is incorrect")
    else  # non-empty
      for symArray1 in tuple[2]  # 1st layer
        if length(symArray1) == 0
          println("THIS sentence is incorrect")
        else  # FIXME: handle default return data structure '[]'
          if typeof(symArray1[1]) <: targetDataType || length(symArray1[1]) == 0
            # println("print 1: simple case --> Array of strings")
            println(symArray1)
          else
            for symArray2 in symArray1  # 2nd layer
              if typeof(symArray2[1]) <: targetDataType
                # println("print 2: complicated case --> Array of arrays")
                println(symArray2)
              else
                for symArray3 in symArray2  # 3rd layer
                  if typeof(symArray3[1]) <: targetDataType
                    println("print 3: too ambiguous to handle")
                    println(symArray3)
                  else
                    for symArray4 in symArray3  # 4th layer
                      if typeof(symArray4[1]) <: targetDataType
                        println("print 4: too ambiguous to handle")
                        println(symArray4)
                      else
                        println(symArray4)
                        println("NEED PRINT 5: world of no man")
                      end
                    end
                  end
                end
              end
            end
          end
        end
        println("")
      end
    end
    if isdefined(tuple, 3)
      println(tuple[3])
    end
    println("-----------------")
  end
end

########### Part 3 Symbol Conversion Dictionary & general biosym replacement ####

#=
F: set up symbol conversion table, delete symbol conversion sentences
I: verb-SymbolGroup pairs
O: conversion dictionary, also modify the input array;
BA
=#
function setUpSymbolConvertionDict(verbSymArray::Array, error::Dict)
  # verbSymArray = Array of pairs of ((senType, line#), BioSym, paraSet)
  conversionTable = Dict{AbstractString, AbstractString}()
  isID = Array{Int64, 1}()
  for (id, vs) in enumerate(verbSymArray)
    if vs[1][1] == "is"  #  find conversion sentence --> check legality
      if ((length(vs[2]) == 2) && (typeof(vs[2][1]) <: Array) && (typeof(vs[2][2]) <: Array)
        && (typeof(vs[2][1][1]) <: AbstractString) && (typeof(vs[2][2][1]) <: AbstractString)
        && (length(vs[2][1]) == length(vs[2][2])))
        for i = 1:length(vs[2][1])
          conversionTable[vs[2][1][i]] = vs[2][2][i]
        end
      else  # incorrect
        lineNum = vs[1][2]
        if !haskey(error, lineNum)
          error[lineNum] = []
        end
        push!(error[lineNum], "incorrect symbol conversion sentence")
        print("incorrect symbol conversion sentence found")
      end
      push!(isID, id)
    end
  end
  deleteat!(verbSymArray, isID)
  conversionTable["[]"] = "[]"
  conversionTable["SITE"] = "SITE"
  return conversionTable
end

mutable struct generalBioSym
  oriBioName::AbstractString
  bioType::AbstractString
  coeff::Float64  # Float64 has default value 0.0?
  function generalBioSym()
    new()
  end
end

#=
F: decompose a token into generalBioSym, check coefficient and type
I: a token, conversion table
O: a generalBioSym if success; otherwise, an empty array []
BA
=#
function BioSym2GeneralBioSym(token::AbstractString, convDict::Dict, err::Array)
  bioSym = generalBioSym()
  if token == "[]"  # special symbol '[]'
    bioSym.oriBioName = "[]"
    bioSym.bioType = "[]"
    bioSym.coeff = 0.0
  elseif (!occursin("*", token) && !occursin("_", token))  # site
    bioSym.oriBioName = token
    bioSym.bioType = "SITE"
    bioSym.coeff = 1.0
  else  # normal symbol
    tokenPiece1 = split(token, "*")
    if length(tokenPiece1) == 0 || length(tokenPiece1) > 2  # "*" is reserved for specifying coefficient
      println("error while parsing coefficient for \'$token\'")
      push!(err, "error while parsing coefficient for \'$token\'")
    else  # maybe 1) correct 2) default 3) error
      if length(tokenPiece1) == 2
        if isnull(tryparse(Float64, tokenPiece1[1]))  # parse() failed
          println("error while parsing coefficient for \'$token\'")
          push!(err, "error while parsing coefficient for \'$token\'")
        else  # got name and coefficient
          bioSym.coeff = parse(Float64, tokenPiece1[1])
          bioSym.oriBioName = tokenPiece1[2]
        end
      else  # default coefficient
        bioSym.coeff = 1.0
        bioSym.oriBioName = tokenPiece1[1]
      end
      # how to continue?
      if isdefined(bioSym, :oriBioName)  # means success in finding coefficient
        tokenPiece2 = split(bioSym.oriBioName, "_")
        if length(tokenPiece2) < 3  # no type or no name specified
          println("no type/compartment specified for \'$token\'")
          push!(err, "no type/compartment specified for \'$token\'")
        else
          tmpType = tokenPiece2[1]*"_"
          if tmpType in keys(convDict)  # want to change into system-defined type? -> YES
            bioSym.bioType = tmpType
          else
            println("incorrect type for \'$token\'")
            push!(err, "incorrect type for \'$token\'")
          end
        end
      end
    end
  end
  # only return when all fields are parsed correctly, otherwise, return empty array []
  if (isdefined(bioSym, :oriBioName) && isdefined(bioSym, :bioType)) # && bioSym.coeff != 0
    return bioSym
  else
    return []
  end
end

#=
F: replace bio-symbol in input array with general biosymbol
I: array of biosymbols or array of arrays of biosymbols, conversion table
O: array of generlBiosymbols or array of arrays of generalBiosymbols
BA
=#
function BioSymArray2generalBioSymArray(inputArr::Array, convDict::Dict, errM::Array)
  postInputArr = []
  if typeof(inputArr[1]) <: AbstractString  # one-layer (simple case)
    postInputArr = [BioSym2GeneralBioSym(s, convDict, errM) for s in inputArr]
  else
    for arr in inputArr
      if typeof(arr[1]) <: AbstractString  # two-layer (complicated case)
        push!(postInputArr, [BioSym2GeneralBioSym(s, convDict, errM) for s in arr])
      else  # three-layer or more, no consider
        println("Currently in land of no man")
        push!(errM, "sorry, currently in land of no man")
      end
    end
  end
  return postInputArr
end

#=
F: replace bio-symbol in input array with general biosymbol
I: array of verb-biosymbolArray pairs, conversion table
O: array of verb-generalBioSymArray pairs
BA
=#
function replaceEachBioSymWithGeneralBioSym(biosymArray::Array,
                     conversionDict::Dict, error_report::Dict)
  # biosymArray = Array of pairs of ((senType, line#), BioSym, paraSet)
  newVerbSymArray = []
  for bs in biosymArray
    err_mes = []
    lineNum = bs[1][2]
    if ((bs[1][1] == "is") || (bs[1][1] == ""))  # fail to delete it
      println("failed to modify the input array?")
      push!(err_mes, "failed to modify the input array in previous steps")
    else
      tmpArr = [BioSymArray2generalBioSymArray(s, conversionDict, err_mes) for s in bs[2]]
      if length(err_mes) == 0
        if isdefined(bs, 3)  # has parameter setting
          push!(newVerbSymArray, (bs[1], tmpArr, bs[3]))
        else
          push!(newVerbSymArray, (bs[1], tmpArr))
        end
      end
    end
    # error occurs
    if length(err_mes) != 0
      if !haskey(error_report, lineNum)
        error_report[lineNum] = []
      end
      append!(error_report[lineNum], err_mes)
    end
  end
  return newVerbSymArray
end

# # NOTE test
# include("preprocessor3.jl")
# println("\n---------loading sentences------------")
# inputSentencesArray = getRidOfNewlineHashcommentBlankline("../test/fbacase.txt")
# tokenizedInputSentencesArray = sentenceTokenization2(inputSentencesArray)
# println("\n------------normalizing and tagging------------")
# reservedWordsPath = "reservedWords.jl"
# taggedSentencesArray = tokenClassification2(tokenizedInputSentencesArray, reservedWordsPath)
# # reshape for output observation
# printTagSen = [["$(y[1])/$(y[2])" for y in x[1]] for x in taggedSentencesArray]
# println(typeof(printTagSen))
# foreach(println, [join(ts, "  ") for ts in printTagSen])
#
# println("\n------------information extraction------------")
# include(reservedWordsPath)
# error_report = Dict{Int64, Array}()
# BioSymVerbInfo = extractBioSymVerbInformation(taggedSentencesArray, reservedWords["SentenceType"], error_report)
# # println("print by foreach:")
# foreach(println, BioSymVerbInfo)
# println("\n------------decoding bio symbols------------")
# decodingBioSymGroups(BioSymVerbInfo, error_report)
# println("\n------------decoding bio symbols results------------")
# printArrayOfTripleToken(BioSymVerbInfo, AbstractString)
# println("\n---------type conversion dictionary---------------")
# typeConversionDict = setUpSymbolConvertionDict(BioSymVerbInfo, error_report)
# foreach(println, typeConversionDict)
# println("\n------replace each biosymbol string with generalBioSym composite-----------------")
# newVerbBiosymInfo = replaceEachBioSymWithGeneralBioSym(BioSymVerbInfo, typeConversionDict, error_report)
# printArrayOfTripleToken(newVerbBiosymInfo, generalBioSym)