text:
#*******************************************************/
#* Copyright(c) 2018 by Artelys                        */
#* This source code is subject to the terms of the     */
#* MIT Expat License (see LICENSE.md)                  */
#*******************************************************/

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This example demonstrates how to use Knitro to solve the following
# simple nonlinear optimization problem.  This model is test problem
# HS40 from the Hock & Schittkowski collection.
#
# max   x0*x1*x2*x3        (obj)
# s.t.  x0^3 + x1^2 = 1    (c0)
#       x0^2*x3 - x2 = 0   (c1)
#       x3^2 - x1 = 0      (c2)
#
# This is the same as exampleNLP2.c, but it demonstrates using
# multiple callbacks for the nonlinear evaluations and computing
# some gradients using finite-differencs, while others are provided
# in callback routines.
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

using KNITRO, Test

function example_multiple_cb(; verbose=true)
    # The signature of this function matches KNITRO.KN_eval_callback in knitro.h.
    # Only "obj" is set in the KNITRO.KN_eval_result structure.
    function callbackEvalObj(kc, cb, evalRequest, evalResult, userParams)
        xind = userParams

        if evalRequest.evalRequestCode != KNITRO.KN_RC_EVALFC
            println("*** callbackEvalObj incorrectly called with eval type ", evalRequest.evalRequestCode)
            return -1
        end
        x = evalRequest.x

        # Evaluate nonlinear term in objective
        evalResult.obj[1] = x[xind[1]] * x[xind[2]] * x[xind[3]] * x[xind[4]]

        return 0
    end

    # The signature of this function matches KNITRO.KN_eval_callback in knitro.h.
    # Only "c0" is set in the KNITRO.KN_eval_result structure.
    function callbackEvalC0(kc, cb, evalRequest, evalResult, userParams)
        xind = userParams

        if evalRequest.evalRequestCode != KNITRO.KN_RC_EVALFC
            println("*** callbackEvalC0 incorrectly called with eval type ", evalRequest.evalRequestCode)
            return -1
        end
        x = evalRequest.x

        # Evaluate nonlinear terms in constraint, c0
        evalResult.c[1] = x[xind[1]] * x[xind[1]] * x[xind[1]]

        return 0
    end

    # The signature of this function matches KNITRO.KN_eval_callback in knitro.h.
    # Only "c1" is set in the KNITRO.KN_eval_result structure.
    function callbackEvalC1(kc, cb, evalRequest, evalResult, userParams)
        xind = userParams

        if evalRequest.evalRequestCode != KNITRO.KN_RC_EVALFC
            println("*** callbackEvalC1 incorrectly called with eval type %d" % evalRequest.evalRequestCode)
            return -1
        end
        x = evalRequest.x

        # Evaluate nonlinear terms in constraint, c1
        evalResult.c[1] = x[xind[1]] * x[xind[1]] * x[xind[4]]

        return 0
    end

    #*------------------------------------------------------------------*
    #*     GRADIENT EVALUATION CALLBACKS                                *
    #*------------------------------------------------------------------*

    # The signature of this function matches KNITRO.KN_eval_callback in knitro.h.
    # Only "objGrad" is set in the KNITRO.KN_eval_result structure.
    function callbackEvalObjGrad(kc, cb, evalRequest, evalResult, userParams)
        xind = userParams

        if evalRequest.evalRequestCode != KNITRO.KN_RC_EVALGA
            println("*** callbackEvalObjGrad incorrectly called with eval type %d" % evalRequest.evalRequestCode)
            return -1
        end
        x = evalRequest.x

        # Evaluate nonlinear terms in objective gradient
        evalResult.objGrad[xind[1]] = x[xind[2]] * x[xind[3]] * x[xind[4]]
        evalResult.objGrad[xind[2]] = x[xind[1]] * x[xind[3]] * x[xind[4]]
        evalResult.objGrad[xind[3]] = x[xind[1]] * x[xind[2]] * x[xind[4]]
        evalResult.objGrad[xind[4]] = x[xind[1]] * x[xind[2]] * x[xind[3]]

        return 0
    end

    # The signature of this function matches KNITRO.KN_eval_callback in knitro.h.
    # Only gradient of c0 is set in the KNITRO.KN_eval_result structure.
    function callbackEvalC0Grad(kc, cb, evalRequest, evalResult, userParams)
        xind = userParams

        if evalRequest.evalRequestCode != KNITRO.KN_RC_EVALGA
            println("*** callbackEvalC0Grad incorrectly called with eval type ", evalRequest.evalRequestCode)
            return -1
        end
        x = evalRequest.x

        # Evaluate nonlinear terms in c0 constraint gradients
        evalResult.jac[1] = 3.0 * x[xind[1]] * x[xind[1]] # *  derivative of x0^3 term  wrt x0

        return 0
    end

    #*------------------------------------------------------------------*
    #*     main                                                         *
    #*------------------------------------------------------------------*

    # Create a new Knitro solver instance.
    kc = KNITRO.KN_new()

    # Initialize Knitro with the problem definition.

    # Add the variables and specify initial values for them.
    # Note: any unset lower bounds are assumed to be
    # unbounded below and any unset upper bounds are
    # assumed to be unbounded above.
    xIndices = KNITRO.KN_add_vars(kc, 4)
    for x in xIndices
        KNITRO.KN_set_var_primal_init_values(kc, x, 0.8)
    end

    # Add the constraints and set the rhs and coefficients
    cIndices = KNITRO.KN_add_cons(kc, 3)
    KNITRO.KN_set_con_eqbnds(kc, cIndices[1], 1.0)
    KNITRO.KN_set_con_eqbnds(kc, cIndices[2], 0.0)
    KNITRO.KN_set_con_eqbnds(kc, cIndices[3], 0.0)

    # Coefficients for 2 linear terms
    lconIndexCons = Int32[1, 2]
    lconIndexVars = Int32[2, 1]
    lconCoefs = [-1.0, -1.0]
    KNITRO.KN_add_con_linear_struct(kc, lconIndexCons, lconIndexVars, lconCoefs)

    # Coefficients for 2 quadratic terms

    #* x1^2 term in c0
    #* x3^2 term in c2
    qconIndexCons = Int32[0, 2]
    qconIndexVars1 = Int32[1, 3]
    qconIndexVars2 = Int32[1, 3]
    qconCoefs = [1.0, 1.0]

    KNITRO.KN_add_con_quadratic_struct(kc, qconIndexCons, qconIndexVars1, qconIndexVars2, qconCoefs)

    # Add separate callbacks.

    # Set callback data for nonlinear objective term.
    cbObj = KNITRO.KN_add_objective_callback(kc, callbackEvalObj)
    KNITRO.KN_set_cb_grad(kc, cbObj, callbackEvalObjGrad,
                        objGradIndexVars=xIndices, nV=length(xIndices))

    # Set callback data for nonlinear constraint 0 term.
    cbC0 = KNITRO.KN_add_eval_callback(kc, false, [cIndices[1]], callbackEvalC0)
    indexCons = cIndices[1]  # constraint c0
    indexVars = xIndices[1]  # variable x0
    KNITRO.KN_set_cb_grad(kc, cbC0, callbackEvalC0Grad, jacIndexCons=[indexCons],
                        jacIndexVars=[indexVars])

    # Set callback data for nonlinear constraint 1 term
    cbC1 = KNITRO.KN_add_eval_callback(kc, false, [cIndices[2]], callbackEvalC1)
    indexCons = [cIndices[2], cIndices[2]]  # constraint c1
    indexVars = [xIndices[1], xIndices[4]]  # variables x0 and x3
    KNITRO.KN_set_cb_grad(kc, cbC1, nothing, jacIndexCons=indexCons, jacIndexVars=indexVars)
    # This one will be approximated via forward finite differences.
    KNITRO.KN_set_cb_gradopt(kc, cbC1, KNITRO.KN_GRADOPT_FORWARD)

    # Demonstrate passing a userParams structure to the evaluation
    # callbacks.  Here we pass back the variable indices set from
    # KNITRO.KN_add_vars() for use in the callbacks. Here we pass the same
    # userParams structure to each callback but we could define different
    # userParams for different callbacks.  This could be useful, for
    # instance, if different callbacks operate on different sets of
    # variables.
    # We increment the indices as Julia is 1-indexed
    userEval = xIndices .+ 1
    KNITRO.KN_set_cb_user_params(kc, cbObj, userEval)
    KNITRO.KN_set_cb_user_params(kc, cbC0, userEval)
    KNITRO.KN_set_cb_user_params(kc, cbC1, userEval)

    # Set minimize or maximize(if not set, assumed minimize)
    KNITRO.KN_set_obj_goal(kc, KNITRO.KN_OBJGOAL_MAXIMIZE)

    # Approximate hessian using BFGS
    KNITRO.KN_set_param(kc, "hessopt", KNITRO.KN_HESSOPT_BFGS)

    kn_outlev = verbose ? KNITRO.KN_OUTLEV_ALL : KNITRO.KN_OUTLEV_NONE
    KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_OUTLEV, kn_outlev)

    # Solve the problem.
    #
    # Return status codes are defined in "knitro.h" and described
    # in the Knitro manual.
    nStatus = KNITRO.KN_solve(kc)
    nStatus, objSol, x, lambda_ =  KNITRO.KN_get_solution(kc)

    # An example of obtaining solution information.
    if verbose
        println("Knitro converged with final status = ", nStatus)
        println("  optimal objective value  = ", objSol)
        println("  optimal primal values x  = ", x)
        println("  feasibility violation    = ", KNITRO.KN_get_abs_feas_error(kc))
        println("  KKT optimality violation = ", KNITRO.KN_get_abs_opt_error(kc))
    end

    # Delete the Knitro solver instance.
    KNITRO.KN_free(kc)

    @testset "Example multipleCB" begin
        @test nStatus == 0
        @test objSol ≈ 0.25
        @test x ≈ [0.793701, 0.707107, 0.529732, 0.840896] atol=1e-5
    end
end

example_multiple_cb(; verbose=isdefined(Main, :KN_VERBOSE) ? KN_VERBOSE : true)


text:
const fp_version_clean = r"[^\w-]"
const fp_cache_regex = r"^v[\w-]+m[0-9a-fA-F]+$"

function build_fingerprint(path::AbstractString, version, hash_value)
    path_parts = split(path, '/')
    (filename, extension) = split(path_parts[end], '.', limit = 2)
    return string(
        join(vcat(path_parts[1:end-1], filename), '/'),
        ".v", replace(string(version), fp_version_clean=>"_"),
        'm', hash_value,
        '.', extension
    )
end

function parse_fingerprint_path(path::AbstractString)
    path_parts = split(path, '/')
    name_parts = split(path_parts[end], '.')
    if length(name_parts) > 2 && occursin(fp_cache_regex, name_parts[2])
        origin_path = string(
            join(vcat(path_parts[1:end-1], [name_parts[1]]), '/'),
            '.', join(name_parts[3:end], '.')
        )
        return (origin_path, true)
    end
    return (path, false)
end

ext: jl


text:
# Autogenerated wrapper script for MPIR_jll for x86_64-linux-gnu-cxx03
export libmpir

JLLWrappers.@generate_wrapper_header("MPIR")
JLLWrappers.@declare_library_product(libmpir, "libmpir.so.23")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libmpir,
        "lib/libmpir.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Documenter

push!(LOAD_PATH,"../src/")
using FindPeaks1D

makedocs(
    sitename = "FindPeaks1D.jl",
    format = Documenter.HTML(prettyurls = false),
    pages = [
        "index.md",
        "findpeaks1d.md",
    ],
)

deploydocs(
  repo = "github.com/ymtoo/FindPeaks1D.jl",
  versions = ["stable" => "v^", "v#.#", "dev" => "master"],
  branch = "gh-pages",
)


text:
using AbstractNumbers, Makie, LinearAlgebra, GeometryTypes
import AbstractNumbers: number
import AbstractPlotting: limits, to_screen, to_world, scene_limits
function to_screen(scene::Scene, mpos)
    return Point2f0(mpos) .- Point2f0(minimum(pixelarea(scene)[]))
end

abstract type Unit{T} <: AbstractNumber{T} end

number(x::Unit) = x.value

"""
Unit space of the scene it's displayed on.
Also referred to as data units
"""
struct SceneSpace{T} <: Unit{T}
    value::T
end
AbstractNumbers.basetype(::Type{<: SceneSpace}) = SceneSpace

"""
Unit is relative to bounding frame.
E.g. if the area is IRect(0, 0, 100, 100)
Point(0.5rel, 0.5rel) == Point(50, 50)
"""
struct Relative{T <: Number} <: Unit{T}
    value::T
end
AbstractNumbers.basetype(::Type{<: Relative}) = Relative
const rel = Relative(1)

"""
https://en.wikipedia.org/wiki/Device-independent_pixel
A device-independent pixel (also: density-independent pixel, dip, dp) is a
physical unit of measurement based on a coordinate system held by a
computer and represents an abstraction of a pixel for use by an
application that an underlying system then converts to physical pixels.
"""
struct DeviceIndependentPixel{T <: Number} <: Unit{T}
    value::T
end
const DIP = DeviceIndependentPixel
const dip = DIP(1)
const dip_in_millimeter = 0.15875
const dip_in_inch = 1/160

AbstractNumbers.basetype(::Type{<: DIP}) = DIP

"""
Unit in pixels on screen.
This one is a bit tricky, since it refers to a static attribute (pixels on screen don't change)
but since every visual is attached to a camera, the exact scale might change.
So in the end, this is just relative to some normed camera - the value on screen, depending on the camera,
will not actually sit on those pixels. Only camera that guarantees the correct mapping is the
`:pixel` camera type.
"""
struct Pixel{T} <: Unit{T}
    value::T
end
AbstractNumbers.basetype(::Type{<: Pixel}) = Pixel
(::Type{Pixel{T}})(x::Pixel{T}) where T = x
const px = Pixel(1)

"""
Millimeter on screen. This unit respects the dimension and pixel density of the screen
to represent millimeters on the screen. This is the must use unit for layouting,
that needs to look the same on all kind of screens. Similar as with the [`Pixel`](@ref) unit,
a camera can change the actually displayed dimensions of any object using the millimeter unit.
"""
struct Millimeter{T} <: Unit{T}
    value::T
end
AbstractNumbers.basetype(::Type{<: Millimeter}) = Millimeter
const mm = Millimeter(1)

Base.show(io::IO, x::DIP) = print(io, number(x), "dip")
Base.:(*)(a::Number, b::DIP) = DIP(a * number(b))

dpi(scene::Scene) = events(scene).window_dpi[]

function pixel_per_mm(scene)
    dpi(scene) ./ 25.4
end

function Base.convert(::Type{<: Millimeter}, scene::Scene, x::SceneSpace)
    pixel = convert(Pixel, scene, x)
    Millimeter(number(pixel_per_mm(scene) / pixel))
end

function Base.convert(::Type{<: SceneSpace}, scene::Scene, x::Relative{T}) where T
    rel = maximum(widths(scene_limits(scene)[])) .* number(x)
    SceneSpace(rel)
end
function Base.convert(::Type{<: SceneSpace}, scene::Scene, x::Point{2, Relative{T}}) where T
    idx = Vec(1, 2)
    lims = scene_limits(scene)
    rel = widths(lims)[idx] .* number.(x)
    SceneSpace(origin(lims)[idx] .+ rel)
end

function Base.convert(::Type{<: SceneSpace}, scene::Scene, x::DIP)
    mm = convert(Millimeter, scene, x)
    SceneSpace(number(mm * dip_in_millimeter))
end

function Base.convert(::Type{<: Millimeter}, scene::Scene, x::DIP)
    Millimeter(number(x * dip_in_millimeter))
end

function Base.convert(::Type{<: Pixel}, scene::Scene, x::Millimeter)
    px = pixel_per_mm(scene) * x
    Pixel(number(px))
end
function Base.convert(::Type{<: Pixel}, scene::Scene, x::DIP)
    inch = (x * dip_in_inch)
    dots = dpi(scene) * inch
    Pixel(number(dots))
end
function Base.convert(::Type{<: SceneSpace}, scene::Scene, x::DIP)
    px = convert(Pixel, scene, x)
    convert(SceneSpace, scene, px)
end

function Base.convert(::Type{<: SceneSpace}, scene::Scene, x::Point{2, <: Pixel})
    s = to_world(scene, to_screen(scene, number.(x)))
    SceneSpace.(s)
end

function Base.convert(::Type{<: SceneSpace}, scene::Scene, x::Vec{2, <: Pixel})
    zero = to_world(scene, to_screen(scene, Point2f0(0)))
    s = to_world(scene, to_screen(scene, number.(Point(x))))
    SceneSpace.(Vec(s .- zero))
end

function Base.convert(::Type{<: SceneSpace}, scene::Scene, x::Pixel)
    zero = to_world(scene, to_screen(scene, Point2f0(0)))
    s = to_world(scene, to_screen(scene, Point2f0(number(x), 0.0)))
    SceneSpace(norm(s .- zero))
end

function Base.convert(::Type{<: SceneSpace}, scene::Scene, x::Millimeter)
    pix = convert(Pixel, scene, x)
    (SceneSpace, mm)
end

# function unpack_plot(plot::AbstractPlot; unlift = false, convert = true)
#     result = Dict()
#
#
#
# function get_boundingbox(x)

x  = [Point2(0px), Point2(3px), Point2(4px)]
Pixel <: AbstractNumber
typemax(Pixel{Int})
bb = Rect(x)
scene = scatter(rand(4))
to_world(scene, Point2f0(500))
a = number.(convert(SceneSpace, scene, Point2(1rel)))
b = number.(convert(SceneSpace, scene, Point2(0rel)))

lines!(Rect(a, b.-a), raw = true)


text:
# Solve model using ensemble parameters
# The reference for the model is Horvath et al. 2020
# paper: https://doi.org/10.1016/j.mec.2019.e00113
# github: https://github.com/varnerlab/Kinetic-CFPS-Model-Publication-Code
# modified by N Laohakunakorn
# University of Edinburgh 2020

include("./model/scripts/DataFile.jl");
include("./model/scripts/MassBalances.jl");
using Sundials
using Plots
using DataFrames
using CSV
using ProgressBars

# Global settings: paths and plot
PATH_OUT = "./results/"
PATH_NET = pwd()*"/model/network/Network.dat"
PATH_PARAMS = pwd()*"/model/params"
FN = "CAT_timeseries.pdf" # Output filename

# Important: select which parameter sets to solve with
#PARAMS_TO_SOLVE = collect(1:100) # all parameter sets
PARAMS_TO_SOLVE = [101] # best fit parameter set

Plots.pyplot()
fntsm = Plots.font("sans-serif", pointsize=round(10.0))
fntlg = Plots.font("sans-serif", pointsize=round(12.0))
default(titlefont=fntlg, guidefont=fntlg, tickfont=fntlg, legendfont=fntsm)

# 1. Load model stoichiometry, initial
# conditions, and parameters.

TSTART = 0.0
TSTOP = 3.0
Ts = 0.1

# Create the data_dictionary and initialise it with network stoichiometric matrix
# and all params to 0 (rates), 1 (initial conds, saturation, order),
# and 0.1 (gain).

# Then populate with parameters from file.

PARAMDICT = []
folders = readdir(PATH_PARAMS*"/Ensemble/") # We have 100 parameter sets

for folder in folders

    data_dictionary = DataFile(TSTART,TSTOP,Ts,PATH_NET)

    # Import initial conditions
    initial_condition_vector = vec(broadcast(abs, float(open(readdlm,PATH_PARAMS*"/initial_condition.dat"))));
    data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_vector;
    # Import rate constants
    rate_constant_vector = vec(broadcast(abs, float(open(readdlm,PATH_PARAMS*"/Ensemble/"*folder*"/rate_constant.dat"))));
    data_dictionary["RATE_CONSTANT_ARRAY"] = rate_constant_vector;
    # Import saturation constants
    saturation_constant_array = broadcast(abs, float(open(readdlm,PATH_PARAMS*"/Ensemble/"*folder*"/saturation_constant.dat")));
    data_dictionary["SATURATION_CONSTANT_ARRAY"] = saturation_constant_array;
    # Import control parameters
    control_constant_array = broadcast(abs, float(open(readdlm,PATH_PARAMS*"/Ensemble/"*folder*"/control_constant.dat")));
    data_dictionary["CONTROL_PARAMETER_ARRAY"] = control_constant_array;

    push!(PARAMDICT, data_dictionary)
end

# Ready to solve
# 2. Solve ODEs, plot and save results
# selected parameter sets defined at top of script

p = plot(legend=:false,xaxis="Time (h)",yaxis="CAT (μM)")
for idx in PARAMS_TO_SOLVE
    @info "Solving $(idx)"
    # Set simulation timespan and initial conditions
    tspan = (TSTART,TSTOP);
    initial_condition_vector = PARAMDICT[idx]["INITIAL_CONDITION_ARRAY"];
    # Use sundials to solve fluxes.
    prob = ODEProblem(MassBalances,initial_condition_vector,tspan,PARAMDICT[idx]);
    sol = solve(prob,CVODE_BDF(),abstol=1e-12,reltol=1e-12); # very small tolerances required
    df = DataFrame(sol);  # df is 352 rows x 2508 columns
    #CSV.write(string(PATH_OUT,"solution.csv"), df) # uncomment if you want to save

    # make arrays and plot
    t=sol.t
    x=Array(df[98,:])
    plot!(t,x*1000,alpha=0.7,color="skyblue1")
end

savefig(PATH_OUT*FN) # Save plot
p # Show plot

max_stars_repo_name: Laohakunakorn-Group/nllab-CFPSmetabolic

max_issues_repo_name: Laohakunakorn-Group/nllab-CFPSmetabolic

max_forks_repo_name: Laohakunakorn-Group/nllab-CFPSmetabolic


text:
using Test
using ExtendableGrids

##########################################
# Correctness test
function test_adj_correctness()

    println("A0: Adjacency transpose from matrix")
    A0=[i+3*(j-1) for i=1:3, j=1:3]

    A0t=atranspose(A0)
    println("A0:\n", A0)
    println("A0t:\n", A0t)

    @test A0t == VariableTargetAdjacency{Int64}([1, 1, 1, 2, 2, 2, 3, 3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10])

    println("A1: Adjacency transpose from variable target adj")
    A1=VariableTargetAdjacency()
    append!(A1, (1,2,3))
    append!(A1, (4,5,6,7))
    append!(A1, (8,9))
    A1t=atranspose(A1)

    println("A1:\n", A1)
    println("A1t:\n", A1t)

    @test A1t == VariableTargetAdjacency{Int64}([1, 1, 1, 2, 2, 2, 2, 3, 3], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10])

    println("A2: Adjacency transpose from variable target adj")
    A2=VariableTargetAdjacency()
    append!(A2, (1,2,3))
    append!(A2, (2,3,4))
    append!(A2, (3,4,5))

    A2t=atranspose(A2)
    println("A2:\n", A2)
    println("A2t:\n", A2t)
    @test A2t == VariableTargetAdjacency{Int64}([1, 1, 2, 1, 2, 3, 2, 3, 3], [1, 2, 4, 7, 9, 10])
end

#############################################
# Performance test

# Function to be called with adjaceny for checking access
# performance
function sum_adj(a)
    sum=0
    for irun=1:100
        for isource::Int64=1:num_sources(a)
            for itarget=1:num_targets(a,isource)
                sum+=a[itarget,isource]
            end
        end
    end
    sum
end

struct TestStruct{T}
    adj::Adjacency{T}
end

struct TestStructAny
    adj
end

function test_adj_performance(n=100_000)
    m=3
    matrix=[i*j for i=1:m, j=1:n]

    ########################################################
    # Access matrix sas FixedTargeAdjacency, should be fast
    print("Direct matrix access:")
    @time begin
        for irun=1:100
            sum=0
            for isource=1:num_sources(matrix)
                for itarget=1:num_targets(matrix,isource)
                    sum+=matrix[itarget,isource]
                end
            end
        end
    end

    # Use this in call, should be fast, but
    # in fact it is slower
    print("    Call with matrix:")
    sum_adj(matrix)
    @time sum_adj(matrix)
    println()

    ########################################################
    # Reorganize adjacency as variable target adjacency
    vadj=VariableTargetAdjacency(matrix)

    # Direct acces should be slower than matrix access
    print("  Direct vadj access:")
    @time begin
        for irun=1:100
            sum=0
            for isource=1:num_sources(vadj)
                for itarget=1:num_targets(vadj,isource)
                    sum+=vadj[itarget,isource]
                end
            end
        end
    end

    # Direct acces should be slower than matrix access
    print("      Call with vadj:")
    @time sum_adj(vadj)
    println()

    ############################################
    # Assign matrix to adj::Adjacency in struct
    # As Adjacency is a union, this could lead to perfromance problems,
    # but we have "UnionSplitting": https://julialang.org/blog/2018/08/union-splitting/
    tstruct=TestStruct(matrix)
    tadj=tstruct.adj

    # Here we see a slowdown compared to the matrix case as
    # there must be some dispatch in the Union
    print("  Direct tadj access:")
    @time begin
        for irun=1:100
            sum=0
            for isource=1:num_sources(tadj)
                for itarget=1:num_targets(tadj,isource)
                    sum+=tadj[itarget,isource]
                end
            end
        end
    end

    # Here the dispatch due to Union splitting
    # allows for the same performance as for the matrix case
    print("      Call with tadj:")
    @time sum_adj(tadj)
    println()

    ############################################
    # Assign matrix to adj::Adjacency in struct
    # As Adjacency is a union, this could lead to perfromance problems,
    # but we have "UnionSplitting": https://julialang.org/blog/2018/08/union-splitting/
    tstruct_v=TestStruct(vadj)
    tvadj=tstruct.adj

    # Here we see a slowdown compared to the matrix case as
    # there must be some dispatch in the Union
    print("  Direct tvadj access:")
    @time begin
        for irun=1:100
            sum=0
            for isource=1:num_sources(tvadj)
                for itarget=1:num_targets(tvadj,isource)
                    sum+=tvadj[itarget,isource]
                end
            end
        end
    end

    # Here the dispatch due to Union splitting
    # allows for the same performance as for the matrix case
    print("      Call with tvadj:")
    @time sum_adj(tvadj)
    println()

    ############################################
    # Assign matrix to adj::Any in struct

    tstruct_a=TestStructAny(matrix)
    aadj=tstruct_a.adj

    # Here we have unboxing with every access so this is
    # really expensive
    print("  Direct aadj access:")
    @time begin
        for irun=1:100
            sum=0
            for isource=1:num_sources(aadj)
                for itarget=1:num_targets(aadj,isource)
                    sum+=aadj[itarget,isource]
                end
            end
        end
    end

    # After "unboxing" when passin the data, specialized
    # code for matrix is called
    print("      Call with aadj:")
    @time sum_adj(aadj)
    """
    Summary: if we want to have a flexible structure, we could either use
    only VariableTargetAdjacency or Union{Matrix,VariableTargetAdjacency}
    """
    nothing
end


text:
function den_invariant(x,sigx_inv);
  # Compute Density --- scale included
  m = size(x,1);
  xc = chol(Hermitian(sigx_inv))*x;
  #xc = factorize(sigx_inv).factors*x;
  c = sqrt(det(sigx_inv))*0.5*gamma(m/2)/(pi^(m/2));
  den = c*((sum(xc.^2,1)).^(-m/2))';
  return den
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2020, DLR Institute of System Dynamics and Control

"""
    module ModiaPlot

Organize and plot simulation result data (time series).
The [`ModiaPlot.plot`](@ref) function of this
module allows to plot the result data by giving the signal names.
The legends/labels of the plots are automatically constructed by the
signal names and their unit. Example:

```julia
ModiaPlot.plot(result, [ (:phi,:r)      (:phi,:phi2,:w);
                         (:w,:w2,:phi2) (:phi,:w)      ],
               heading="Matrix of plots")
```

generates the following plot:

![Matrix-of-Plots](../../resources/images/matrix-of-plots.svg)

# Main developer

[Martin Otter](https://rmc.dlr.de/sr/en/staff/martin.otter/),
[DLR - Institute of System Dynamics and Control](https://www.dlr.de/sr/en)
"""
module ModiaPlot

export plot, showFigure, saveFigure, closeFigure, closeAllFigures, printResultInfo

# using/imports
import DataFrames
import Colors
import Measurements
import MonteCarloMeasurements
using  Unitful

# Constants
const headingSize = 10

const path = dirname(dirname(@__FILE__))   # Absolute path of package directory
const Version = "0.7.1"
const Date = "2021-03-31"

println("Importing ModiaPlot Version $Version ($Date) - this takes some time due to GLMakie import")

include("plot_base.jl")
include("plot_with_GLMakie.jl")

end


text:
type Conv <: Op; padding; stride; upscale; mode; algorithm; workSpace; workSpaceSizeInBytes; alpha; beta;
    function Conv(; padding=0, stride=1, upscale=1, mode=CUDNN_CONVOLUTION,
                  algorithm=CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_GEMM,
                  workSpace=nothing, workSpaceSizeInBytes=0, alpha=1.0, beta=0.0,
                  o...)
        new(padding, stride, upscale, mode, algorithm, workSpace, workSpaceSizeInBytes, alpha, beta)
    end
end

ninputs(::Conv)=2
canoverwrite(::Conv)=false
back_reads_x(::Conv)=true
back_reads_y(::Conv)=false

function forw(c::Conv, w, x, y; o...)
    if w == nothing
        error("Uninitialized filter")
    elseif x == nothing
        return nothing
    end
    ws = (c.workSpace == nothing ? C_NULL : c.workSpace)
    cudnnConvolutionForward(x, w, y; padding=c.padding, stride=c.stride, upscale=c.upscale, mode=c.mode,
                            algorithm=c.algorithm, workSpace=ws, workSpaceSizeInBytes=c.workSpaceSizeInBytes,
                            alpha=c.alpha, beta=c.beta)
    gpusync()
    return y
end

function back(c::Conv, dy, dw, dx; x=nothing, o...)
    dw == nothing && dx == nothing && return
    dw != nothing && (x[2] != nothing ? cudnnConvolutionBackwardFilter(x[2], dy, dw; padding=c.padding, stride=c.stride, upscale=c.upscale, mode=c.mode) : fillsync!(dw,0))
    dx != nothing && (x[1] != nothing ? cudnnConvolutionBackwardData(x[1], dy, dx; padding=c.padding, stride=c.stride, upscale=c.upscale, mode=c.mode) : error("Uninitialized filter"))
    gpusync()
end

# x: (x1,x2...,C,N)
# w: (w1,w2...,C,K)
# y: (y1,y2...,K,N)
# If padding=0 and stride=1: yi=xi-wi+1
# In general we have: yi = 1 + (xi + 2*padding - wi) / stride

function infersize(c::Conv,w,x,y)
    if w==x==y==nothing
        nothing
    elseif w==nothing || x==nothing || y==nothing
        n = (w!=nothing ? length(w) : x!=nothing ? length(x) : y!=nothing ? length(y) : error())
        w == nothing && (w = ntuple(i->0, n))
        x == nothing && (x = ntuple(i->0, n))
        y == nothing && (y = ntuple(i->0, n))
        infersize(c,w,x,y)
    else
        s = (isa(c.stride,  Integer) ? ntuple(i->c.stride,  length(y)-2) : c.stride)
        p = (isa(c.padding, Integer) ? ntuple(i->c.padding, length(y)-2) : c.padding)
        length(w) == length(x) == length(y) == length(s)+2 == length(p)+2 || throw(DimensionMismatch())
        w = [w...]; x = [x...]; y = [y...]
        for i=1:length(s)
            if w[i] > 0 && x[i] > 0 && y[i] > 0
                y[i] == 1 + div(x[i] + 2*p[i] - w[i], s[i]) || throw(DimensionMismatch())
            elseif w[i] > 0 && x[i] > 0 && y[i] == 0
                y[i] = 1 + div(x[i] + 2*p[i] - w[i], s[i])
            elseif w[i] > 0 && x[i] == 0 && y[i] > 0 && s[i] == 1
                x[i] = y[i] - 1 - 2*p[i] + w[i]
            end
        end
        n = length(x)
        equate!(x,n-1,w,n-1)
        equate!(x,n,y,n)
        equate!(w,n,y,n-1)
        w = tuple(w...); x = tuple(x...); y = tuple(y...)
        (w,x,y)
    end
end

function equate!(a,i,b,j)
    a[i] == b[j] ? nothing :
    a[i] == 0 ? a[i] = b[j] :
    b[j] == 0 ? b[j] = a[i] :
    throw(DimensionMismatch())
end

### DEAD CODE

# TODO: generalize to N-D
# TODO: cpu implementation
# TODO: upgrade to new cudnn version
# TODO: upgrade to new knet interface

# type Conv <: Op; w; x; ybuf; dx; Conv(p::KUparam)=new(p); end

# Conv(d...; o...)=Conv(KUparam(d...; o...))
# Conv(nout::Integer, width::Integer; o...)=Conv(KUparam(width, 0, nout; o...))

# params(l::Conv)=Any[l.w]
# ninputs(::Conv)=1
# overwrites(::Conv)=false
# back_reads_x(::Conv)=true
# back_reads_y(::Conv)=false

# # TODO: this unnecessarily allocates w and y
# ysize(l::Conv, x)=(isempty(l.w) && initforw(l,x,nothing); cudnnGetConvolutionNdForwardOutputDim(x,l.w))

# function forw(l::Conv, x; y=nothing, o...)
#     l.x = x
#     y = initforw(l,x,y)
#     cudnnConvolutionForward(x, l.w, y)
# end

# function back(l::Conv, dy; dx=nothing, x=l.x, incr=false, returndx=true, o...)
#     initback(l, dy, x, incr)
#     if incr
#         cudnnConvolutionBackwardFilter(x, dy, l.w.inc)
#         axpy!(1, l.w.inc, l.w.diff)
#     else
#         cudnnConvolutionBackwardFilter(x, dy, l.w.diff)
#     end
#     if returndx
#         dx = initbackx(l,x,dx)
#         cudnnConvolutionBackwardData(l.w, dy, dx)
#     end
# end

# function initback(l::Conv, dy, x, incr)
#     atype(dy) == atype(x) || error("atype mismatch")
#     eltype(dy) == eltype(x) || error("eltype mismatch")
#     size(dy) == ysize(l,x) || error("ysize mismatch")
#     similar!(l.w, :diff, l.w.arr)
#     incr && similar!(l.w, :inc, l.w.arr)
# end

# function initbackx(l::Conv, x, dx)
#     dx == nothing && (dx = similar!(l, :dx, x))
#     issimilar(dx,x) || error("Gradient mismatch")
#     return dx
# end

# # TODO: We should split up the w and y parts and share with Mmul

# function initforw(l::Conv, x, y)
#     n = ndims(x)
#     c = size(x)[n-1]  # x dims are (x1, x2, ..., channels, images)
#     if isempty(l.w)
#         nz(l.w,:init,nothing) || (l.w.init = xavier!)
#         r = size(l.w, 1)
#         o = size(l.w, ndims(l.w))
#         wsize = ntuple(i->(i<n-1 ? r : i==n-1 ? c : o), n)
#         init(l.w, eltype(x), wsize)
#     end
#     eltype(x) == eltype(l.w) || "$(eltype(x)) != $(eltype(l.w))"
#     n == ndims(l.w) || error("ndims mismatch")
#     c == size(l.w)[n-1] || error("channel mismatch")
#     ys = ysize(l,x)
#     y == nothing && (y = similar!(l, :ybuf, x, ys))
#     typeof(y) == typeof(x) || error("Type mismatch")
#     size(y) == ys || error("Size mismatch")
#     return y
# end

# xavier!(a)=(fanin = length(a) / (size(a)[end]); scale = sqrt(3 / fanin); rand!(a, -scale, scale); a)

# # Make things work with KUdense

# import CUDNN: cudnnGetConvolutionNdForwardOutputDim, cudnnConvolutionForward, cudnnConvolutionBackwardFilter, cudnnConvolutionBackwardData

# cudnnGetConvolutionNdForwardOutputDim(x::KUdense, w::KUparam)=cudnnGetConvolutionNdForwardOutputDim(x.arr, w.arr)
# cudnnConvolutionForward(x::KUdense, w::KUparam, y::KUdense)=(cudnnConvolutionForward(x.arr, w.arr, y.arr);y)
# cudnnConvolutionBackwardFilter(x::KUdense, dy::KUdense, w::BaseArray)=(cudnnConvolutionBackwardFilter(x.arr, dy.arr, w);w)
# cudnnConvolutionBackwardData(w::KUparam, dy::KUdense, dx::KUdense)=(cudnnConvolutionBackwardData(w.arr, dy.arr, dx.arr);dx)

# Make things work with CPU (for now)

# cudnnGetConvolutionNdForwardOutputDim(x::Array, w::Array)=cudnnGetConvolutionNdForwardOutputDim(CudaArray(x),CudaArray(w))
# cudnnConvolutionForward(x::Array, w::Array, y::Array)=(y1=CudaArray(y);cudnnConvolutionForward(CudaArray(x), CudaArray(w), y1);copysync!(y,1,y1,1,length(y)))
# cudnnConvolutionBackwardFilter(x::Array, dy::Array, w::Array)=(w1=CudaArray(w);cudnnConvolutionBackwardFilter(CudaArray(x), CudaArray(dy), w1); copysync!(w,1,w1,1,length(w)))
# cudnnConvolutionBackwardData(w::Array, dy::Array, dx::Array)=(dx1=CudaArray(dx);cudnnConvolutionBackwardData(CudaArray(w), CudaArray(dy), dx1); copysync!(dx,1,dx1,1,length(dx)))

#     else

#     @assert length(w) == length(x) == length(y)
#     nd = length(x)
#     x = [x...]
#     w = [w...]
#     x[nd-1] == 0 && (x[nd-1] = w[nd-1])
#     w[nd-1] == 0 && (w[nd-1] = x[nd-1])
#     @assert x[nd-1] == w[nd-1]
#     y = zeros(x)
#     for i=1:nd-2
#         w[i] > 0 && x[i] > 0 && (y[i] = x[i]-w[i]+1)
#     end
#     y[nd-1] = w[nd]
#     y[nd] = x[nd]
#     return (tuple(w...), tuple(x...), tuple(y...))
# end

# # TODO: write doc
# """
# @knet function conv(w, x; padding=0, stride=1, upscale=1, mode=CUDNN_CONVOLUTION)
# """
# function conv(w, x, y; padding=0, stride=1, upscale=1, mode=CUDNN_CONVOLUTION)
#     @assert in(mode, (CUDNN_CONVOLUTION, CUDNN_CROSS_CORRELATION))
#     (Conv(padding, stride, upscale, mode), w, x, y)
# end

# using CUDNN: cudnnConvolutionDescriptor_t


text:
using Zygote, Test, ChainRules

@testset "ChainRules Integration" begin
    @testset "basic" begin
        cr_inner_demo_rrule_hitcount = Ref(0)
        cr_inner_demo_pullback_hitcount = Ref(0)
        cr_inner_demo(x) = 5x
        function ChainRules.rrule(::typeof(cr_inner_demo), x)
            cr_inner_demo_rrule_hitcount[] += 1
            function cr_inner_demo_pullback(Δx)
                cr_inner_demo_pullback_hitcount[] += 1
                return ChainRules.NO_FIELDS, 5.0*Δx
            end
            return cr_inner_demo(x), cr_inner_demo_pullback
        end

        function cr_outer_demo(x)
            2 + 10cr_inner_demo(x)
        end

        @testset "gradient inner" begin
            cr_inner_demo_rrule_hitcount[] = 0
            cr_inner_demo_pullback_hitcount[] = 0
            @test (5.0,) == gradient(cr_inner_demo, 11)
            @test cr_inner_demo_rrule_hitcount[] == 1
            @test cr_inner_demo_pullback_hitcount[] == 1
        end

        @testset "gradient outer" begin
            cr_inner_demo_rrule_hitcount[] = 0
            cr_inner_demo_pullback_hitcount[] = 0
            @test (50.0,) == gradient(cr_outer_demo, 11)
            @test cr_inner_demo_rrule_hitcount[] == 1
            @test cr_inner_demo_pullback_hitcount[] == 1
        end

        @testset "pullback inner" begin
            cr_inner_demo_rrule_hitcount[] = 0
            cr_inner_demo_pullback_hitcount[] = 0
            y, pb = pullback(cr_inner_demo, 11)
            @test y == 55
            @test cr_inner_demo_rrule_hitcount[] == 1
            @test cr_inner_demo_pullback_hitcount[] == 0
            @test pb(1)==(5.0,);
            @test pb(2)==(10.0,);
            @test pb(3)==(15.0,);
            @test cr_inner_demo_pullback_hitcount[] == 3
            @test cr_inner_demo_rrule_hitcount[] == 1
        end
    end

    @testset "Multiple output single input" begin
        simo_rrule_hitcount = Ref(0)
        simo_pullback_hitcount = Ref(0)
        simo(x) = (5x, 7x)
        function ChainRules.rrule(::typeof(simo), x)
            simo_rrule_hitcount[] += 1
            function simo_pullback((Δa, Δb))
                simo_pullback_hitcount[] += 1
                return ChainRules.NO_FIELDS, 5*Δa + 7*Δb
            end
            return simo(x), simo_pullback
        end

        simo_outer(x) = sum(simo(x))

        @assert simo_rrule_hitcount[] == 0
        @assert simo_pullback_hitcount[] == 0
        @test (12,) == Zygote.gradient(simo_outer, π)
        @test simo_rrule_hitcount[] == 1
        @test simo_pullback_hitcount[] == 1
    end

    @testset "multiple input, Single output" begin
        miso_rrule_hitcount = Ref(0)
        miso_pullback_hitcount = Ref(0)
        miso(a, b) = 5a + 7b
        function ChainRules.rrule(::typeof(miso), a, b)
            miso_rrule_hitcount[] += 1
            function miso_pullback(Δy)
                miso_pullback_hitcount[] += 1
                return ChainRules.NO_FIELDS, 5Δy , 7Δy
            end
            return miso(a, b), miso_pullback
        end

        miso_outer(x) = miso(100x, 10x)

        @assert miso_rrule_hitcount[] == 0
        @assert miso_pullback_hitcount[] == 0
        @test (570,) == Zygote.gradient(miso_outer, π)
        @test miso_rrule_hitcount[] == 1
        @test miso_pullback_hitcount[] == 1
    end

    @testset "multiple input multiple output" begin
        mimo_rrule_hitcount = Ref(0)
        mimo_pullback_hitcount = Ref(0)
        mimo(a, b) = (5a + 7b, 100a, 10b)
        function ChainRules.rrule(::typeof(mimo), a, b)
            mimo_rrule_hitcount[] += 1
            function mimo_pullback((Δx, Δy, Δz))
                mimo_pullback_hitcount[] += 1
                return ChainRules.NO_FIELDS, 5Δx + 100Δy , 7Δx + 10Δz
            end
            return mimo(a, b), mimo_pullback
        end

        @assert mimo_rrule_hitcount[] == 0
        @assert mimo_pullback_hitcount[] == 0
        _, pb = Zygote.pullback(mimo, π, 2π)
        @test (105, 17) == pb((1, 1, 1))
        @test mimo_rrule_hitcount[] == 1
        @test mimo_pullback_hitcount[] == 1

        mimo_outer(x) = sum(mimo(x, x))

        mimo_rrule_hitcount[] = 0
        mimo_pullback_hitcount[] = 0
        @test (122,) == gradient(mimo_outer, π)
        @test mimo_rrule_hitcount[] == 1
        @test mimo_pullback_hitcount[] == 1
    end

    @testset "all AbstractZero partials" begin
        # while ChainRules always has a partial for every input, Zygote combined them all
        # to a single `nothing` if they are all zero-like.

        not_diff_eg(x, i) = [10, 20][i]
        function ChainRules.rrule(::typeof(not_diff_eg), x, i)
            function not_diff_eg_pullback(Δ)
                return ChainRules.NO_FIELDS, ChainRules.Zero(), ChainRules.DoesNotExist()
            end
            return not_diff_eg(x, i), not_diff_eg_pullback
        end

        _, pb = Zygote.pullback(not_diff_eg, 10.4, 2)
        @test pb(1.2) === nothing
    end

    @testset "nested AD hitting identity(::Tuple) pullback" begin
        # This is is  a particularly fiddly case.
        # Its kind of a simplified version of `sin'''(0.5)` but different in some places.

        f(x) = tuple(x, 2x, 3x)

        function g(y)
            a1, pb1 = Zygote.pullback(f, π)
            pb1((y,0,0))
        end

        @test (1,) == g(1)

        function h(n)
            a2, pb2 = Zygote.pullback(g, 1)
            pb2(n)
        end

        @test (1,) == h(1)

        if VERSION >= v"1.3"
            # broken on Julia 1.0 because of https://github.com/FluxML/Zygote.jl/issues/638
            a3, pb3 = Zygote.pullback(h, 1)
            @test ((1,),) == pb3(1)
        end
    end

    @testset "kwargs" begin
        kwfoo_rrule_hitcount = Ref(0)
        kwfoo_pullback_hitcount = Ref(0)
        kwfoo(x; k=10) = x + k
        function ChainRules.rrule(::typeof(kwfoo), x; k=10)
            kwfoo_rrule_hitcount[] += 1
            function kwfoo_pullback(Δy)
                kwfoo_pullback_hitcount[] += 1
                return ChainRules.NO_FIELDS, Δy
            end
            return kwfoo(x; k=k), kwfoo_pullback
        end

        kwfoo_outer_unused(x) = kwfoo(x)
        kwfoo_outer_used(x) = kwfoo(x; k=-15)

        @testset "$outer" for outer in (kwfoo_outer_used, kwfoo_outer_unused)
            kwfoo_rrule_hitcount[] = 0
            kwfoo_pullback_hitcount[] = 0
            @test (1,) == Zygote.gradient(outer, π)
            @test kwfoo_rrule_hitcount[] == 1
            @test kwfoo_pullback_hitcount[] == 1
        end
    end

    @testset "kwarg, with all AbstractZero partials" begin
        # while ChainRules always has a partial for every input, Zygote combined them all
        # to a single `nothing` if they are all zero-like.

        not_diff_kw_eg(x, i; kw=1.0) = [10, 20][i]
        function ChainRules.rrule(::typeof(not_diff_kw_eg), x, i; kwargs...)
            function not_diff_kw_eg_pullback(Δ)
                return ChainRules.NO_FIELDS, ChainRules.Zero(), ChainRules.DoesNotExist()
            end
            return not_diff_kw_eg(x, i; kwargs...), not_diff_kw_eg_pullback
        end

        @test (nothing,) == Zygote.gradient(x->not_diff_kw_eg(x, 2), 10.4)
        @test (nothing,) == Zygote.gradient(x->not_diff_kw_eg(x, 2; kw=2.0), 10.4)
    end
end

@testset "FastMath support" begin
    @test gradient(2.0) do x
      @fastmath x^2.0
    end == (4.0,)
end


text:
rstd = 3.0
astd = 0.1

rrm = 5.0
rrstd = 15.0
aam = 0.0
aastd = 0.2

fg = initfg()

addVariable!(fg, :t0, Polar)
addFactor!(fg, [:t0], PriorPolar(Normal(ranges[1],rstd), Normal(angles[1],astd)) )

for tidx in 1:(nposes-1)
  psym = Symbol("t$(tidx-1)")
  sym = Symbol("t$tidx")
  addVariable!(fg, sym, Polar)
  addFactor!(fg, [sym], PriorPolar(Normal(ranges[tidx+1],rstd), Normal(angles[tidx+1],astd)) )
  addFactor!(fg, [psym; sym], PolarPolar(Normal(rrm, rrstd), Normal(aam, aastd)) )
end


text:
import Pkg;
Pkg.add("JSON")
Pkg.add("JSON3")
Pkg.add("LazyJSON")

import JSON;
import JSON3;
import LazyJSON;

open("array_02_numeric.json", "r") do file
    global dicttxt = read(file, String)  # file information to string
    global a1 = LazyJSON.value(dicttxt)  # parse and transform data
end

a2 = JSON3.read(dicttxt)

println("a2 = ", a2)
println("a1 = ", a1)

Pkg.add("Statistics")
Pkg.add("StatsBase")
using Statistics;
using StatsBase;

# https://rosettacode.org/wiki/Averages/Arithmetic_mean#Julia

#m1 = mean(a1)
mean_arithmetic_2 = mean(a2)
mode_2 = StatsBase.mode(a2)

#println(m1)
println("mean arithmetic (m2)   = ", mean_arithmetic_2)
println("mode (m2)              = ", mode_2)

dict_result = Dict(
                    "mean arithmetic" => mean_arithmetic_2,
                    "mode" => mode_2
                )
string_data_results = JSON.json(dict_result)

# write the file with the stringdata variable information
open("results.json", "w") do f
    write(f, string_data_results)
end


text:
module Main

using LearnConvection.Data
using LearnConvection.GaussianProcess

model_output(x, time_index, ℳ, 𝒟) = GaussianProcess.model_output(𝒟.modify_predictor_fn(x, time_index), ℳ)

include("errors.jl")
export  get_me_true_check, # evolving forward from an arbitrary initial timestep
        get_me_greedy_check, # how well does the mean GP prediction fit the training data?
        disparity_vector

include("plot_profile.jl")
export  plot_profile,
        plot_model_output,
        animate_profile,
        animate_profile_and_model_output

include("predict.jl")
export predict

# plot hyperparameter landscapes for analysis / optimization
include("hyperparameter_landscapes.jl")
export  plot_landscapes_compare_error_metrics,
        plot_landscapes_compare_files_me,
        plot_error_histogram,
        get_min_gamma,
        get_min_gamma_alpha,
        train_validate_test

include("optimize_smp.jl")
export  optimize_SMP_kernel

end #module

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
@doc raw"""
    ∇distance(M,y,x[, p=2])

compute the (sub)gradient of the distance (squared)

```math
f(x) = \frac{1}{2} d^p_{\mathcal M}(x,y)
```

to a fixed point `y` on the manifold `M` and `p` is an
integer. The gradient reads

```math
  \nabla f(x) = -d_{\mathcal M}^{p-2}(x,y)\log_xy
```

for $p\neq 1$ or $x\neq  y$. Note that for the remaining case $p=1$,
$x=y$ the function is not differentiable. In this case, the function returns the
corresponding zero tangent vector, since this is an element of the subdifferential.

# Optional

* `p` – (`2`) the exponent of the distance,  i.e. the default is the squared
  distance
"""
∇distance(M,y,x,p::Int=2) = (p==2) ? -log(M,x,y) : -distance(M,x,y)^(p-2)*log(M,x,y)

@doc raw"""
    ∇u,∇v = ∇intrinsic_infimal_convolution_TV12(M,f,u,v,α,β)

compute (sub)gradient of the intrinsic infimal convolution model using the mid point
model of second order differences, see [`costTV2`](@ref), i.e. for some $f ∈ \mathcal M$
on a `PowerManifold` manifold $\mathcal M$ this function computes the (sub)gradient of

```math
E(u,v) =
\frac{1}{2}\sum_{i ∈ \mathcal G} d_{\mathcal M}(g(\frac{1}{2},v_i,w_i),f_i)
+ \alpha
\bigl(
\beta\mathrm{TV}(v) + (1-\beta)\mathrm{TV}_2(w)
\bigr),
```
where both total variations refer to the intrinsic ones, [`∇TV`](@ref) and
[`∇TV2`](@ref), respectively.
"""
function ∇intrinsic_infimal_convolution_TV12(M::mT,f,u,v,α,β) where {mT <: Manifold}
  c = mid_point(M,u,v,f)
  iL = log(M,c,f)
  return adjoint_differential_geodesic_startpoint(M,u,v,1/2,iL) + α*β*∇TV(M,u), adjoint_differential_geodesic_endpoint(M,u,v,1/2,iL) + α * (1-β) * ∇TV2(M,v)
end
@doc raw"""
    ∇TV(M,(x,y),[p=1])

compute the (sub) gradient of $\frac{1}{p}d^p_{\mathcal M}(x,y)$ with respect
to both $x$ and $y$.
"""
function ∇TV(M::MT, xT::Tuple{T,T}, p=1)where {MT <: Manifold, T}
  x = xT[1];
  y = xT[2];
  if p==2
      return (-log(M,x,y), -log(M,y,x))
  else
    d = distance(M,x,y);
    if d==0 # subdifferential containing zero
      return (zero_tangent_vector(M,x),zero_tangent_vector(M,y))
    else
      return (-log(M,x,y)/(d^(2-p)), -log(M,y,x)/(d^(2-p)))
    end
  end
end
@doc raw"""
    ξ = ∇TV(M,λ,x,[p])
Compute the (sub)gradient $\partial F$ of all forward differences orrucirng,
in the power manifold array, i.e. of the function

$F(x) = \sum_{i}\sum_{j ∈ \mathcal I_i} d^p(x_i,x_j)$

where $i$ runs over all indices of the `PowerManifold` manifold `M`
and $\mathcal I_i$ denotes the forward neighbors of $i$.

# Input
* `M` – a `PowerManifold` manifold
* `x` – a point.

# Ouput
* ξ – resulting tangent vector in $T_x\mathcal M$.
"""
function ∇TV(M::PowerManifold,x,p::Int=1)
    power_size = power_dimensions(M)
    rep_size = representation_size(M.manifold)
    R = CartesianIndices(Tuple(power_size))
    d = length(power_size)
    maxInd = last(R)
    X = zero_tangent_vector(M,x)
    c = costTV(M,x,p,0)
    for i in R # iterate over all pixel
        di = 0.
        for k in 1:d # for all direction combinations
            ek = CartesianIndex(ntuple(i  ->  (i==k) ? 1 : 0, d) ) #k th unit vector
            j = i+ek # compute neighbor
            if all( map(<=, j.I, maxInd.I)) # is this neighbor in range?
                if p != 1
                    g = (c[i]==0 ? 1 : 1/c[i]) .* ∇TV(M.manifold,(x[i],x[j]),p) # Compute TV on these
                else
                    g = ∇TV(M.manifold,(x[i],x[j]),p) # Compute TV on these
                end
                X[i] += g[1]
                X[j] += g[2]
            end
        end # directions
    end # i in R
    return X
end

@doc raw"""
    ξ = forward_logs(M,x)

compute the forward logs $F$ (generalizing forward differences) orrucirng,
in the power manifold array, the function

```math
$F_i(x) = \sum_{j ∈ \mathcal I_i} \log_{x_i} x_j,\quad i  ∈  \mathcal G,
```

where $\mathcal G$ is the set of indices of the `PowerManifold` manifold `M`
and $\mathcal I_i$ denotes the forward neighbors of $i$.

# Input
* `M` – a `PowerManifold` manifold
* `x` – a point.

# Ouput
* `ξ` – resulting tangent vector in $T_x\mathcal M$ representing the logs, where
  $\mathcal N$ is thw power manifold with the number of dimensions added to `size(x)`.
"""
function forward_logs(M::PowerManifold{𝔽,TM,TSize,TPR}, p) where {𝔽,TM,TSize,TPR}
    power_size = power_dimensions(M)
    R = CartesianIndices(Tuple(power_size))
    d = length(power_size)
    sX = size(p)
    maxInd = [last(R).I...] # maxInd as Array
    if d > 1
        d2 = fill(1,d+1)
        d2[d+1] = d
    else
        d2 = 1
    end
    N = PowerManifold(M.manifold, TPR(), power_size..., d)
    xT = repeat(p,inner=d2)
    X = zero_tangent_vector(N,xT)
    for i in R # iterate over all pixel
        for k in 1:d # for all direction combinations
            I = [i.I...] # array of index
            J = I .+ 1 .* (1:d .== k) #i + e_k is j
            if all( J .<= maxInd ) # is this neighbor in range?
                j = CartesianIndex{d}(J...) # neigbbor index as Cartesian Index
                set_component!(
                    N,
                    X,
                    log(M.manifold, get_component(M,p,i), get_component(M,p,j)),
                    i,
                    k,
                ) # Compute log and store in kth entry
            end
        end # directions
    end # i in R
    return X
end

@doc raw"""
    ∇TV2(M,(x,y,z),p)

computes the (sub) gradient of $\frac{1}{p}d_2^p(x,y,z)$ with respect
to $x$, $y$, and $z$, where $d_2$ denotes the second order absolute difference
using the mid point model, i.e. let
```math
  \mathcal C = \bigl\{ c ∈  \mathcal M \ |\ g(\tfrac{1}{2};x_1,x_3) \text{ for some geodesic }g\bigr\}
```
denote the mid points between $x$ and $z$ on the manifold $\mathcal M$.
Then the absolute second order difference is defined as

```math
d_2(x,y,z) = \min_{c ∈ \mathcal C_{x,z}} d(c,y).
```

While the (sub)gradient with respect to $y$ is easy, the other two require
the evaluation of an [`adjoint_Jacobi_field`](@ref).
See [Illustration of the Gradient of a Second Order Difference](@ref secondOrderDifferenceGrad)
for its derivation.
"""
function ∇TV2(M::MT, xT, p::Number=1) where {MT <: Manifold}
  x = xT[1];
  y = xT[2];
  z = xT[3];
  c = mid_point(M,x,z,y) # nearest mid point of x and z to y
  d = distance(M,y,c)
  innerLog = -log(M,c,y)
  if p==2
      return ( adjoint_differential_geodesic_startpoint(M,x,z,1/2,innerLog), -log(M,y,c), adjoint_differential_geodesic_endpoint(M,x,z,1/2,innerLog))
  else
    if d==0 # subdifferential containing zero
      return (zero_tangent_vector(M,x),zero_tangent_vector(M,y),zero_tangent_vector(M,z))
    else
      return ( adjoint_differential_geodesic_startpoint(M,x,z,1/2,innerLog/(d^(2-p))), -log(M,y,c)/(d^(2-p)), adjoint_differential_geodesic_endpoint(M,x,z,1/2,innerLog/(d^(2-p))) )
    end
  end
end
@doc raw"""
    ∇TV2(M,q [,p=1])

computes the (sub) gradient of $\frac{1}{p}d_2^p(x_1,x_2,x_3)$
with respect to all $x_1,x_2,x_3$ occuring along any array dimension in the
point `x`, where `M` is the corresponding `PowerManifold`.
"""
function ∇TV2(M::PowerManifold, q, p::Int=1)
    power_size = power_dimensions(M)
    rep_size = representation_size(M.manifold)
    R = CartesianIndices(Tuple(power_size))
    d = length(power_size)
    minInd, maxInd = first(R), last(R)
    X = zero_tangent_vector(M,q)
    c = costTV2(M,q,p,false)
    for i in R # iterate over all pixel
        di = 0.
        for k in 1:d # for all direction combinations (TODO)
            ek = CartesianIndex(ntuple(i  ->  (i==k) ? 1 : 0, d) ) #k th unit vector
            jF = i+ek # compute forward neighbor
            jB = i-ek # compute backward neighbor
            if all( map(<=, jF.I, maxInd.I) ) && all( map(>=, jB.I, minInd.I)) # are neighbors in range?
                if p != 1
                    g = (c[i] == 0 ? 1 : 1/c[i]) .* ∇TV2(M.manifold,(q[jB],q[i],q[jF]),p) # Compute TV2 on these
                else
                    g = ∇TV2(M.manifold,(q[jB],q[i],q[jF]),p) # Compute TV2 on these
                end
                set_component!(M,X,g[1],jB)
                set_component!(M,X,g[2],i)
                set_component!(M,X,g[3],jF)
            end
        end # directions
    end # i in R
    return X
end


text:
#########################################################################################################
#
#  The following shader functions have inputs that are all in the form (vx, vy, vz, nx, ny, nz, args...)
#  and return boolean, new ray, a multiplication/scattering factor (scalar or array)
#
#########################################################################################################

##############################################################################################
#
# Photon shaders... (single waveband, includes a probability to be absorbed)
#
##############################################################################################

function lambertian_full_fate_shader(Ix::T, Iy::T, Iz::T, vx::T, vy::T, vz::T, nx::T, ny::T, nz::T, rho::T) where {T<:AbstractFloat}
    # note that vx, vy, and vz are not used but only exist so that input arguments are of the same pattern as in other shader functions...
    scattering = false
    newray = Ray(Inf, Inf, Inf, Inf, Inf, Inf)
    r = 0.0
    if rand() <= rho
        scattering = true
        t, p = sample_f() # could be replaced with a function RPVSampler if we knew how to implement importance sampling for RPV...
	fx, fy, fz = spherical2carthesian(t, p)
	sx, sy, sz = project2normal(fx, fy, fz, nx, ny, nz)
	newray = propagate(Ix, Iy, Iz, sx, sy, sz)
	r = rho * cos(t)
    end
    return scattering, newray, r
end

function bilambertian_full_fate_shader(Ix::T, Iy::T, Iz::T, vx::T, vy::T, vz::T, nx::T, ny::T, nz::T, rho::T, tau::T) where {T<:AbstractFloat}
    # note that vx, vy, and vz are not used but only exist so that input arguments are of the same pattern as in other shader functions...
    scattering = false
    newray = Ray(Inf, Inf, Inf, Inf, Inf, Inf)
    r = 0.0
    omega = rho + tau # total scatter...
    if rand() <= omega
        scattering = true
        # choose either reflection or transmission...
        if rand() > (rho/omega)
            nx, ny, nz = -1.0*nx, -1.0*ny, -1.0*nz
        end
        # now, treats as BSDF and the total scatter is assigned to the direction of propagation
        t, p = sample_f()
        fx, fy, fz = spherical2carthesian(t, p)
        sx, sy, sz = project2normal(fx, fy, fz, nx, ny, nz)
        newray = propagate(Ix, Iy, Iz, sx, sy, sz)
        r = omega * cos(t)
    end
    return scattering, newray, r
end

function rpv_full_fate_shader(vx::T, vy::T, vz::T, nx::T, ny::T, nz::T, rho_zero::T, rho_c::T, captheta::T, k::T) where {T<:AbstractFloat}
    # note that vx, vy, and vz are not used but only exist so that input arguments are of the same pattern as in other shader functions...
    scattering = false
    newray = Ray(Inf, Inf, Inf, Inf, Inf, Inf)
    rho_sfc = 0.0
    if rand() <= rho_zero # 7SEPT2021: I PICKED RHO_ZERO BUT MAYBE IT SHOULD BE RHO_C INSTEAD?
        scattering = true
        t, p = sample_f() # could be replaced with a function RPVSampler if we knew how to implement importance sampling for RPV...
        fx, fy, fz = spherical2carthesian(t, p)
        vx, vy, vz = normalize(vx, vy, vz)
        nx, ny, nz = normalize(nx, ny, nz)
        sx, sy, sz = project2world(fx, fy, fz, nx, ny, nz)
        phi = relazi(vn, vy, vz, nx, ny, nz, sx, sy, sz)
        theta = abs(acos(vx*nx + vy*ny + vz*nz))
        theta_zero = abs(acos(sx*nx + sy*ny + sz*nz))
        rho_sfc = rpv(theta_zero, theta, phi, rho_zero, rho_c, captheta, k)
        newray = propagate(Ix, Iy, Iz, sx, sy, sz)
    end
    return scattering, newray, rho_sfc
end

function rpv_full_fate_shader(Ix::T, Iy::T, Iz::T, vx::T, vy::T, vz::T, nx::T, ny::T, nz::T, rho_zero::Array{T}, rho_c::Array{T}, captheta::Array{T}, k::Array{T}) where {T<:AbstractFloat}
    # note that vx, vy, and vz are not used but only exist so that input arguments are of the same pattern as in other shader functions...
    scattering = false
    newray = Ray(Inf, Inf, Inf, Inf, Inf, Inf)
    rho_sfc = similar(rho_zero)
    if rand() <= sum(rho_zero)/length(rho_zero) # 7SEPT2021: I PICKED RHO_ZERO BUT MAYBE IT SHOULD BE RHO_C INSTEAD?
        scattering = true
        t, p = sample_f() # could be replaced with a function RPVSampler if we knew how to implement importance sampling for RPV...
        fx, fy, fz = spherical2carthesian(t, p)
        vx, vy, vz = normalize(vx, vy, vz)
        nx, ny, nz = normalize(nx, ny, nz)
        sx, sy, sz = project2normal(fx, fy, fz, nx, ny, nz)
        phi = relazi(vx, vy, vz, nx, ny, nz, sx, sy, sz)
        theta = abs(acos(vx*nx + vy*ny + vz*nz))
        theta_zero = abs(acos(sx*nx + sy*ny + sz*nz))
        rho_sfc = rpv(theta_zero, theta, phi, rho_zero, rho_c, captheta, k)
        newray = propagate(Ix, Iy, Iz, sx, sy, sz)
    end
    return scattering, newray, rho_sfc
end

##############################################################################################
#
# Photon shaders... (spectra, includes probability of absorption)
#
##############################################################################################

function lambertian_full_fate_shader(Ix::T, Iy::T, Iz::T, vx::T, vy::T, vz::T, nx::T, ny::T, nz::T, rho::AbstractArray{T,1}) where {T<:AbstractFloat}
    # note that vx, vy, and vz are not used but only exist so that input arguments are of the same pattern as in other shader functions...
    scattering = false
    newray = Ray(Inf, Inf, Inf, Inf, Inf, Inf)
    r = 0.0
    if rand() <= sum(rho)/length(rho)
        scattering = true
        t, p = sample_f() # could be replaced with a function RPVSampler if we knew how to implement importance sampling for RPV...
	fx, fy, fz = spherical2carthesian(t, p)
	sx, sy, sz = project2normal(fx, fy, fz, nx, ny, nz)
	newray = propagate(Ix, Iy, Iz, sx, sy, sz)
	r = rho * cos(t)
    end
    return scattering, newray, r
end

function bilambertian_full_fate_shader(Ix::T, Iy::T, Iz::T, vx::T, vy::T, vz::T, nx::T, ny::T, nz::T, rho::AbstractArray{T,1}, tau::AbstractArray{T,1}) where {T<:AbstractFloat}
    # note that vx, vy, and vz are not used but only exist so that input arguments are of the same pattern as in other shader functions...
    scattering = false
    newray = Ray(Inf, Inf, Inf, Inf, Inf, Inf)
    r = 0.0
    omega = rho + tau # total scatter...
    if rand() <= sum(omega)/length(omega)
        scattering = true
        # choose either reflection or transmission...
        if rand() > (sum(rho)/sum(omega))
            nx, ny, nz = -1.0*nx, -1.0*ny, -1.0*nz
        end
        # now, treat as BSDF: the total scatter is assigned to the direction of propagation
        t, p = sample_f()
        fx, fy, fz = spherical2carthesian(t, p)
        sx, sy, sz = project2normal(fx, fy, fz, nx, ny, nz)
        newray = propagate(Ix, Iy, Iz, sx, sy, sz)
        r = omega * cos(t)
    end
    return scattering, newray, r
end

##############################################################################################
#
# Single-waveband shaders...
#
##############################################################################################

function lambertian_shader(Ix::T, Iy::T, Iz::T, vx::T, vy::T, vz::T, nx::T, ny::T, nz::T, rho::T) where {T<:AbstractFloat}
    # note that vx, vy, and vz are not used but only exist so that input arguments are of the same pattern as in other shader functions...
    t, p = sample_f() # could be replaced with a function RPVSampler if we knew how to implement importance sampling for RPV...
    fx, fy, fz = spherical2carthesian(t, p)
    sx, sy, sz = project2normal(fx, fy, fz, nx, ny, nz)
    newray = propagate(Ix, Iy, Iz, sx, sy, sz)
    r = rho * cos(t)
    return true, newray, r
end

function bilambertian_shader(Ix::T, Iy::T, Iz::T, vx::T, vy::T, vz::T, nx::T, ny::T, nz::T, rho::T, tau::T) where {T<:AbstractFloat}
    # note that vx, vy, and vz are not used but only exist so that input arguments are of the same pattern as in other shader functions...
    omega = rho + tau # total scatter...
    if rand() > (rho/omega)
        nx, ny, nz = -1.0*nx, -1.0*ny, -1.0*nz
    end
    # now, treats as BSDF and the total scatter is assigned to the direction of propagation
    t, p = sample_f()
    fx, fy, fz = spherical2carthesian(t, p)
    sx, sy, sz = project2normal(fx, fy, fz, nx, ny, nz)
    newray = propagate(Ix, Iy, Iz, sx, sy, sz)
    r = omega * cos(t)
    return true, newray, r
end

function pure_reflection_shader(vx::T, vy::T, vz::T, nx::T, ny::T, nz::T, rho::T) where {T<:AbstractFloat}
    sx, sy, sz = pure_reflection(vx, vy, vz, nx, ny, nz)
    newray = propagate(Ix, Iy, Iz, sx, sy, sz)
    return true, newray, rho
end

function rpv_shader(vx::T, vy::T, vz::T, nx::T, ny::T, nz::T, rho_zero::T, rho_c::T, captheta::T, k::T) where {T<:AbstractFloat}
    t, p = sample_f() # could be replaced with a function RPVSampler if we knew how to implement importance sampling for RPV...
    fx, fy, fz = spherical2carthesian(t, p)
    vx, vy, vz = normalize(vx, vy, vz)
    nx, ny, nz = normalize(nx, ny, nz)
    sx, sy, sz = project2world(fx, fy, fz, nx, ny, nz)
    phi = relazi(vn, vy, vz, nx, ny, nz, sx, sy, sz)
    theta = abs(acos(vx*nx + vy*ny + vz*nz))
    theta_zero = abs(acos(sx*nx + sy*ny + sz*nz))
    rho_sfc = rpv(theta_zero, theta, phi, rho_zero, rho_c, captheta, k)
    newray = propagate(Ix, Iy, Iz, sx, sy, sz)
    return true, newray, rho_sfc
end

###############################################################
#
#  Spectral shaders... (i.e. acting on spectra)
#
###############################################################

function lambertian_shader(Ix::T, Iy::T, Iz::T, vx::T, vy::T, vz::T, nx::T, ny::T, nz::T, rho::Array{T}) where {T<:AbstractFloat}
    # note that vx, vy, and vz are not used but only exist so that input arguments are of the same pattern as in other shader functions...
    t, p = sample_f() # could be replaced with a function RPVSampler if we knew how to implement importance sampling for RPV...
    fx, fy, fz = spherical2carthesian(t, p)
    sx, sy, sz = project2normal(fx, fy, fz, nx, ny, nz)
    newray = propagate(Ix, Iy, Iz, sx, sy, sz)
    r = rho * cos(t)
    return true, newray, r
end

function bilambertian_shader(Ix::T, Iy::T, Iz::T, vx::T, vy::T, vz::T, nx::T, ny::T, nz::T, rho::Array{T}, tau::Array{T}) where {T<:AbstractFloat}
    # note that vx, vy, and vz are not used but only exist so that input arguments are of the same pattern as in other shader functions...
    # note also that the whole spectrum undergoes one treatment: either reflect or transmit (obviously.)
    omega = rho + tau
    if rand() > (sum(rho)/sum(omega))
        nx, ny, nz = -1.0*nx, -1.0*ny, -1.0*nz
    end
    # treat as BSDF. All radiance goes to new direction...
    t, p = sample_f() # could be replaced with a function RPVSampler if we knew how to implement importance sampling for RPV...
    fx, fy, fz = spherical2carthesian(t, p)
    sx, sy, sz = project2normal(fx, fy, fz, nx, ny, nz)
    newray = propagate(Ix, Iy, Iz, sx, sy, sz)
    r = omega * cos(t)
    return true, newray, r
end

function rpv_shader(Ix::T, Iy::T, Iz::T, vx::T, vy::T, vz::T, nx::T, ny::T, nz::T, rho_zero::Array{T}, rho_c::Array{T}, captheta::Array{T}, k::Array{T}) where {T<:AbstractFloat}
    t, p = sample_f() # could be replaced with a function RPVSampler if we knew how to implement importance sampling for RPV...
    fx, fy, fz = spherical2carthesian(t, p)
    vx, vy, vz = normalize(vx, vy, vz)
    nx, ny, nz = normalize(nx, ny, nz)
    sx, sy, sz = project2normal(fx, fy, fz, nx, ny, nz)
    phi = relazi(vx, vy, vz, nx, ny, nz, sx, sy, sz)
    theta = abs(acos(vx*nx + vy*ny + vz*nz))
    theta_zero = abs(acos(sx*nx + sy*ny + sz*nz))
    rho_sfc = rpv(theta_zero, theta, phi, rho_zero, rho_c, captheta, k)
    newray = propagate(Ix, Iy, Iz, sx, sy, sz)
    return true, newray, rho_sfc
end

function pure_reflection_shader(Ix::T, Iy::T, Iz::T, vx::T, vy::T, vz::T, nx::T, ny::T, nz::T, rho::Array{T}) where {T<:AbstractFloat}
    sx, sy, sz = pure_reflection(vx, vy, vz, nx, ny, nz)
    newray = propagate(Ix, Iy, Iz, sx, sy, sz)
    return true, newray, rho
end

function lightsource_shader(Ix::T, Iy::T, Iz::T, vx::T, vy::T, vz::T, nx::T, ny::T, nz::T, emission::Array{T}) where {T<:AbstractFloat}
    return false, Ray(Ix, Iy, Iz, vx, vy, vz), emission
end

max_stars_repo_name: martinvanleeuwen/RenderJay

max_issues_repo_name: martinvanleeuwen/RenderJay

max_forks_repo_name: martinvanleeuwen/RenderJay


text:
using QuantEcon
using Plots
plotlyjs(size = (900, 900), legend = false)

include("jv.jl")

wp = JvWorker(grid_size=25)
v_init = wp.x_grid .* 0.5

f(x) = bellman_operator(wp, x)
V = compute_fixed_point(f, v_init, max_iter=200)

s_policy, phi_policy = bellman_operator(wp, V, ret_policies=true)

# === plot solution === #

plot(wp.x_grid, [phi_policy s_policy V],
  title = ["phi policy" "s policy" "value function"],
  color = [:orange :blue :green], width = 3,
  xaxis = ("x", (0.0, maximum(wp.x_grid))),
  yaxis = ((-0.1, 1.1)),
  layout = (3,1), bottom_margin = 20mm, show = true)

max_stars_repo_name: mwhchen/quantecon

max_issues_repo_name: mwhchen/quantecon

max_forks_repo_name: mwhchen/quantecon


text:
export Episode, get_inputs, step!, run!

struct Episode
    robot::Robot
    grid::Grid
    reward::Function
    meta::Dict
end

no_reward(e::Episode) = 0

function Episode(g::Grid; reward::Function=no_reward, meta::Dict=Dict())
    robot = Robot(g)
    meta["step"] = get(meta, "step", 0)
    meta["max_step"] = get(meta, "max_step", 100)
    meta["time"] = get(meta, "time", 0.0)
    meta["max_time"] = get(meta, "max_time", 10.0)
    meta["memory"] = get(meta, "memory", 0)
    meta["max_memory"] = get(meta, "max_memory", 1073741824)
    meta["reward"] = get(meta, "reward", 0.0)
    Episode(robot, g, reward, meta)
end

function get_inputs(e::Episode)
    cids = Array{Int64}(undef, (3, 3))
    linds = LinearIndices(e.grid.cells)
    i = 1
    for dy in -1:1
        for dx in -1:1
            cids[i] = get(linds, (e.robot.y + dy, e.robot.x + dx), 0)
            i += 1
        end
    end
    if e.robot.d > 0
        cids = [rotr90, rot180, rotl90][e.robot.d](cids)
    end
    inputs = Array{Float64}(undef, 10)
    inputs[1] = e.meta["reward"]
    i = 2
    for cid in cids
        c = get(e.grid.cells, cid, WallCell)
        inputs[i] = c.obj
        # inputs[i+1] = c.color
        i += 1
    end
    inputs
end

function step!(e::Episode, action::Int64)
    if action == 1
        move!(e.robot, e.grid)
    elseif action == 2
        move!(e.robot, e.grid; forward=false)
    elseif action == 3
        turn_right!(e.robot)
    elseif action == 4
        turn_left!(e.robot)
    elseif action == 5
        e.grid.cells[e.robot.y, e.robot.x].color = color_to_float("red")
    elseif action == 6
        e.grid.cells[e.robot.y, e.robot.x].color = color_to_float("green")
    end
    e.meta["reward"] = e.reward(e)
    e.meta["reward"]
end

function step!(e::Episode, actions::AbstractArray)
    step!(e, argmax(actions))
end

function step!(e::Episode, controller::Function)
    inputs = get_inputs(e)
    output = controller(inputs)
    step!(e, output)
end

function terminate(e::Episode)
    (e.meta["step"] >= e.meta["max_step"] ||
     e.meta["time"] >= e.meta["max_time"] ||
     e.meta["memory"] >= e.meta["max_memory"])
end

function run!(e::Episode, controller::Function; terminate::Function=terminate)
    linds = LinearIndices(e.grid.cells)
    total_reward = 0.0
    while true
        rew, t, bytes, gctime, memallocs = @timed step!(e, controller)
        total_reward += rew
        e.meta["step"] += 1
        e.meta["time"] += t
        e.meta["memory"] += bytes
        if (terminate(e) || (linds[e.robot.y, e.robot.x] in e.grid.exits))
            break
        end
    end
    total_reward
end


text:
"""
    Report

Internal structure that holds the temporary reported values. See [`report!`](@ref).
"""
struct Report
    data::LittleDict{Symbol,Vector}

    Report() = new(Dict{Symbol,Vector}())
end

function Base.show(io::IO, report::Report)
    print(io, "Report")
    if !isempty(report.data)
        print(":")
        keywidth = maximum(length.(string.(keys(report.data))))
        for (k, v) in report.data
            print("\n  $(lpad(k, keywidth)) => $v")
        end
    end
end
function Base.empty!(report::Report)
    foreach(empty!, values(report.data))
end
function DataFrames.DataFrame(report::Report)
    DataFrame(report.data; copycols=false)
end

const SymbolOrString = Union{Symbol,AbstractString}

"""
    report!(report, keys, values, id="")
    report!(report, pairs, id="")

Write `keys`, `values` pairs to `report` that will be converted to a `DataFrame` later.
Alternatively, a named tuple or a collection of pairs can be passed instead of `keys` and
`values`.

The value of `id` is appended to the name of the column, e.g.
`report!(report, :key, value, :_1)` will report `value` to a column named `:key_1`.
"""
function report!(report::Report, key::SymbolOrString, value)
    data = report.data
    if haskey(data, key)
        column = data[key]::Vector{typeof(value)}
        push!(column, value)
    else
        data[key] = [value]
    end
    return report
end
function report!(report::Report, key::SymbolOrString, value, postfix::SymbolOrString)
    report!(report, Symbol(key, postfix), value)
end
function report!(report::Report, keys, vals, postfix::SymbolOrString="")
    @assert length(keys) == length(vals)
    for (k, v) in zip(keys, vals)
        report!(report, k, v, postfix)
    end
    return report
end
function report!(report::Report, nt::NamedTuple, postfix::SymbolOrString="")
    report!(report, pairs(nt), postfix)
    return report
end
function report!(report::Report, kvpairs, postfix::SymbolOrString="")
    for (k, v) in kvpairs
        report!(report, k, v, postfix)
    end
    return report
end

function ensure_correct_lengths(report::Report)
    lens = length.(values(report.data))
    for i in 1:length(report.data)
        lens[i] == lens[1] || error("duplicate keys reported to `DataFrame`")
    end
end

"""
    ReportingStrategy

Abstract type for strategies for reporting data in a DataFrame with [`report!()`](@ref).

# Implemented strategies:

* [`ReportDFAndInfo`](@ref)
* [`ReportToFile`](@ref)

# Interface:

A `ReportingStrategy` can define any of the following:

* [`refine_r_strat`](@ref)
* [`report!`](@ref)
* [`report_after_step`](@ref)
* [`finalize_report!`](@ref)

"""
abstract type ReportingStrategy end

"""
    refine_r_strat(r_strat::ReportingStrategy) -> r_strat

Initialize the reporting strategy. This can be used to set up filenames or other attributes
that need to be unique for a run of FCIQMC.
"""
refine_r_strat(r_strat::ReportingStrategy) = r_strat

"""
     report!(::ReportingStrategy, step, report::Report, keys, values, id="")
     report!(::ReportingStrategy, step, report::Report, nt, id="")

Report `keys` and `values` to `report`, which will be converted to a `DataFrame` before
[`lomc!`](@ref) exits. Alternatively, a `nt::NamedTuple` can be passed in place of `keys`
and `values`. If `id` is specified, it is appended to all `keys`. This is used to
differentiate between values reported by different replicas.

To overload this function for a new `ReportingStrategy`, overload
`report!(::ReportingStrategy, step, args...)` and apply the report by calling
`report!(args...)`.
"""
function report!(::ReportingStrategy, _, args...)
    report!(args...)
    return nothing
end

"""
    report_after_step(::ReportingStrategy, step, report, state)

This function is called exactly once at the very end of a step. For example, it can be used
to print some information to `stdout`.
"""
function report_after_step(::ReportingStrategy, args...)
    return nothing
end

"""
    finalize_report!(::ReportingStrategy, report)

Finalize the report. This function is called after all steps in [`lomc!`](@ref) have
finished.
"""
function finalize_report!(::ReportingStrategy, report)
    DataFrame(report)
end

function print_stats(io::IO, step, state)
    print(io, "[ ", lpad(step, 11), " | ")
    shift = lpad(round(state.replicas[1].params.shift, digits=4), 10)
    norm = lpad(round(state.replicas[1].pnorm, digits=4), 10)
    println(io, "shift: ", shift, " | norm: ", norm)
    flush(io)
end

"""
    ReportDFAndInfo(; k=1, i=100, io=stdout, writeinfo=false) <: ReportingStrategy

The default [`ReportingStrategy`](@ref). Report every `k`th step to a `DataFrame` and write
info message to `io` every `i`th step (unless `writeinfo == false`). The flag `writeinfo` is
useful for controlling info messages in MPI codes, e.g. by setting
`writeinfo =`[`is_mpi_root()`](@ref).
"""
@with_kw struct ReportDFAndInfo <: ReportingStrategy
    k::Int = 1
    i::Int = 100
    io::IO = stdout
    writeinfo::Bool = false
end
function report!(s::ReportDFAndInfo, step, args...)
    step % s.k == 0 && report!(args...)
    return nothing
end
function report_after_step(s::ReportDFAndInfo, step, _, state)
    if s.writeinfo && step % s.i == 0
        print_stats(s.io, step, state)
    end
end

"""
    ReportToFile(; kwargs...) <: ReportingStrategy

Reporting strategy that writes the report directly to a file. Useful when dealing with long
jobs or large numbers of replicas, when the report can incur a significant memory cost.

# Keyword arguments

* `filename`: the file to report to. If the file already exists, a new file is created.
* `chunk_size = 1000`: the size of each chunk that is written to the file. A `DataFrame` of
  this size is collected in memory and written to disk. When saving, an info message is also
  printed to `io`.
* `save_if =`[`is_mpi_root()`](@ref): if this value is true, save the report, otherwise
  ignore it.
* `return_df = false`: if this value is true, read the file and return the data frame at the
  end of computation. Otherwise, an empty `DataFrame` is returned.
* `io = stdout`: The `IO` to print messages to. Set to `devnull` if you don't want to see
  messages printed out.
"""
@with_kw struct ReportToFile <: ReportingStrategy
    filename::String
    chunk_size::Int = 1000
    save_if::Bool = RMPI.is_mpi_root()
    return_df::Bool = false
    io::IO = stdout
end
function refine_r_strat(s::ReportToFile)
    if s.save_if
        # If filename exists, add -1 to the end of it. If that exists as well,
        # increment the number after the dash
        new_filename = s.filename
        while isfile(new_filename)
            base, ext = splitext(new_filename)
            m = match(r"(.*)-([0-9]+)$", base)
            if isnothing(m)
                new_filename = string(base, "-1", ext)
            else
                new_filename = string(m[1], "-", parse(Int, m[2]) + 1, ext)
            end
        end
        if s.filename ≠ new_filename
            println(s.io, "File `$(s.filename)` exists.")
            s = @set s.filename = new_filename
        end
        println(s.io, "Saving report to `$(s.filename)`.")
    end
    return s
end
function report!(s::ReportToFile, _, args...)
    if s.save_if
        report!(args...)
    end
    return nothing
end
function report_after_step(s::ReportToFile, step, report, state)
    if s.save_if && step % s.chunk_size == 0
        # Report some stats:
        print_stats(s.io, step, state)

        if isfile(s.filename)
            Arrow.append(s.filename, report.data)
        else
            Arrow.write(s.filename, report.data; file=false)
        end
        empty!(report)
    end
end
function finalize_report!(s::ReportToFile, report)
    if s.save_if
        println(s.io, "Finalizing.")
        if isfile(s.filename)
            Arrow.append(s.filename, report.data)
        else
            Arrow.write(s.filename, report.data; file=false)
        end
        if s.return_df
            return DataFrame(Arrow.Table(s.filename))
        end
    end
    return DataFrame()
end


text:
# ideally this should be in the GenerativeAD pkg
function compute_ranks(rt)
	mask_nan_max = (x) -> (isnan(x) ? -Inf : x)
	rs = []
	for row in eachrow(rt)
		push!(rs, StatsBase.competerank(mask_nan_max.(Vector(row)), rev = true))
	end
	# each row represents ranks of a method
	reduce(hcat, rs)
end

const MODEL_RENAME = Dict(
    "aae_full" => "aae",
    "wae_full" => "wae",
    "vae_full" => "vae")

"""
	sorted_rank(d, agg; val_metric=:val_auc, tst_metric=:tst_auc, downsample=Dict{String, Int}())

Computes rank table from `d` and aggregation `agg`. Applies model and dataset aliases.
Sorts models based on model_type. Returns the rank table and also the aggregated dataframe.
"""
function sorted_rank(d, agg, val_metric=:val_auc, tst_metric=:tst_auc, downsample=Dict{String, Int}(); verbose=true)
	df_agg = agg(d, val_metric; downsample=downsample, verbose=verbose)

	df_agg["model_type"] = copy(df_agg["modelname"])
	apply_aliases!(df_agg, col="modelname", d=MODEL_RENAME)
	apply_aliases!(df_agg, col="modelname", d=MODEL_ALIAS)
	apply_aliases!(df_agg, col="dataset", d=DATASET_ALIAS)
	apply_aliases!(df_agg, col="model_type", d=MODEL_TYPE)

	sort!(df_agg, (:dataset, :model_type, :modelname))
	rt = rank_table(df_agg, tst_metric)
	names(rt)[2:end], rt
end

"""
automatic aggregation of combined dataframes for single/multi class datasets

single class datasets are marked with anomaly_class = -1 and will be aggregated usin aggregate_stats_mean_max
multi class datasets are marked with anomaly_class > 0 and will be aggregated usin aggregate_stats_max_mean

"""
function aggregate_stats_auto(df::DataFrame, criterion_col=:val_auc; kwargs...)
	df_single = select(filter(x -> x.anomaly_class == -1, df), Not(:anomaly_class))
	df_multi = filter(x -> x.anomaly_class > 0, df)

	if nrow(df_single) > 0 && nrow(df_multi) > 0
		df_agg_single = aggregate_stats_mean_max(df_single, criterion_col; kwargs...)
		df_agg_multi = aggregate_stats_max_mean(df_multi, criterion_col; kwargs...)
		return	vcat(df_agg_single, df_agg_multi, cols=:intersect)
	elseif nrow(df_single) > 0
		return aggregate_stats_mean_max(df_single, criterion_col; kwargs...)
	else
		return aggregate_stats_max_mean(df_multi, criterion_col; kwargs...)
	end
end

# top 10 mean-max
function aggregate_stats_mean_max_top_10(df::DataFrame,
										criterion_col=:val_auc;
										min_samples=("anomaly_class" in names(df)) ? 10 : 3,
										add_col=nothing, verbose=true)
	agg_cols = vcat(_prefix_symbol.("val", BASE_METRICS), _prefix_symbol.("tst", BASE_METRICS))
	agg_cols = vcat(agg_cols, _prefix_symbol.("val", PAT_METRICS), _prefix_symbol.("tst", PAT_METRICS))
	agg_cols = vcat(agg_cols, _prefix_symbol.("val", PATN_METRICS), _prefix_symbol.("tst", PATN_METRICS))
	agg_cols = vcat(agg_cols, _prefix_symbol.("val", PAC_METRICS), _prefix_symbol.("tst", PAC_METRICS))
	agg_cols = vcat(agg_cols, Symbol.(TRAIN_EVAL_TIMES))
	agg_cols = (add_col !== nothing) ? vcat(agg_cols, add_col) : agg_cols
	top10_std_cols = _prefix_symbol.(agg_cols, "top_10_std")
	std_cols = _prefix_symbol.(agg_cols, "std")

	# agregate by seed over given hyperparameter and then choose best
	results = []
	for (dkey, dg) in pairs(groupby(df, :dataset))
		for (mkey, mg) in pairs(groupby(dg, :modelname))
			pg = groupby(mg, :phash)
			mg_suff = reduce(vcat, [g for g in pg if nrow(g) >= min_samples])
			pg_agg = combine(groupby(mg_suff, :phash),
						nrow => :psamples,
						agg_cols .=> mean .=> agg_cols,
						agg_cols .=> std,
						:parameters => unique => :parameters)

			# sort by criterion_col
			sort!(pg_agg, order(criterion_col, rev=true))

			# take best 10 models
			best_10 = first(pg_agg, 10)
			best = combine(best_10,
							agg_cols .=> mean .=> agg_cols,
							std_cols .=> mean .=> std_cols,
							agg_cols .=> std  .=> top10_std_cols)

			# add grouping keys
			best[:dataset] = dkey.dataset
			best[:modelname] = mkey.modelname

			push!(results, best)
		end
	end
	vcat(results...)
end

# top 10 max-mean
function aggregate_stats_max_mean_top_10(df::DataFrame,
										criterion_col=:val_auc; add_col=nothing, verbose=true)
	agg_cols = vcat(_prefix_symbol.("val", BASE_METRICS), _prefix_symbol.("tst", BASE_METRICS))
	agg_cols = vcat(agg_cols, _prefix_symbol.("val", PAT_METRICS), _prefix_symbol.("tst", PAT_METRICS))
	agg_cols = vcat(agg_cols, _prefix_symbol.("val", PATN_METRICS), _prefix_symbol.("tst", PATN_METRICS))
	agg_cols = vcat(agg_cols, _prefix_symbol.("val", PAC_METRICS), _prefix_symbol.("tst", PAC_METRICS))
	agg_cols = vcat(agg_cols, Symbol.(TRAIN_EVAL_TIMES))
	agg_cols = (add_col !== nothing) ? vcat(agg_cols, add_col) : agg_cols
	top10_std_cols = _prefix_symbol.(agg_cols, "top_10_std")
	std_cols = _prefix_symbol.(agg_cols, "std")

	agg_keys = ("anomaly_class" in names(df)) ? [:seed, :anomaly_class] : [:seed]

	# choose best for each seed/anomaly_class cobination and average over them
	results = []
	for (dkey, dg) in pairs(groupby(df, :dataset))
		for (mkey, mg) in pairs(groupby(dg, :modelname))
			partial_results = []

			if ("anomaly_class" in names(df))
				classes = unique(mg.anomaly_class)
				dif = setdiff(collect(1:10), classes)
				if (length(classes) < 10) && verbose
					@warn "$(mkey.modelname) - $(dkey.dataset): missing runs on anomaly_class $(dif)."
				end
			else
				seeds = unique(mg.seed)
				dif = setdiff(collect(1:5), seeds)
				if (length(seeds) < 3) && verbose
					@warn "$(mkey.modelname) - $(dkey.dataset): missing runs on seed $(dif)."
				end
			end

			# iterate over seed-anomaly_class groups
			for (skey, sg) in pairs(groupby(mg, agg_keys))
				ssg = sort(sg, order(criterion_col, rev=true))
				# best hyperparameter after sorting by criterion_col
				best_10 = first(ssg, 10)
				best = combine(best_10,
							agg_cols .=> mean .=> agg_cols,
							agg_cols .=> std  .=> top10_std_cols)

				push!(partial_results, best)
			end

			best_per_seed = reduce(vcat, partial_results)
			# average over seed-anomaly_class groups
			best = combine(best_per_seed,
						agg_cols .=> mean .=> agg_cols,
						top10_std_cols .=> mean .=> top10_std_cols,
						agg_cols .=> std)

			# add grouping keys
			best[:dataset] = dkey.dataset
			best[:modelname] = mkey.modelname

			push!(results, best)
		end
	end
	vcat(results...)
end


text: Dict{String, Any}("milliseconds" => Dict{String, Any}("value" => "10:32:00.555", "type" => "time-local"), "besttimeever" => Dict{String, Any}("value" => "17:45:00", "type" => "time-local"))


text:
"""
    selectCrossPoint(genotype::CGPGenotype, rng::Random.AbstractRNG;
                     fixedDepth::Integer = -1, minDepth::Integer = -1,
                     maxDepth::Integer = -1)

Selects a point in the genotype where the crossover is going to be performed
by taking some considerations into account, such as the maximum depth allowed.

# Arguments
- `genotype::CGPGenotype`: genotype from which the cross point is going to be
    selected.
- `rng::Random.AbstractRNG`: random number generator for random number consistency
    along an experiment.

# Keyword Arguments
- `fixedDepth::Integer = -1`: specifies the exact depth at where the crosspoint
    must be. If this argument is set to -1, it will be ignored.
- `minDepth::Integer = -1`: specifies the minimum depth at where the crosspoint
    can be. If this argument is set to -1, it will be ignored.
- `maxDepth::Integer = -1`: specifies the maximum depth at where the crosspoint
    can be. If this argument is set to -1, it will be ignored.
"""
function selectCrossPoint(genotype::CGPGenotype, rng::Random.AbstractRNG;
                          fixedDepth::Integer = -1, minDepth::Integer = -1,
                          maxDepth::Integer = -1)

    if fixedDepth != -1 && minDepth != -1
        minDepth = -1
    end

    representation = getRepresentation(genotype)
    lenRep = length(representation)

    # Sets depths for every node
    depths = Array{Integer}(undef, lenRep)
    visitedFromRoot = zeros(Integer, 1, 2)
    nodeIndex = 1
    depths[nodeIndex] = 0
    parent = 0
    depth = 0

    while nodeIndex != 0

        if nodeIndex > size(visitedFromRoot)[1]
            visitedFromRoot = vcat(visitedFromRoot, zeros(Integer, 1, 2))
            depths[nodeIndex] = depth
            visitedFromRoot[end, 2] = parent
        end

        if visitedFromRoot[nodeIndex, 1] < getArity(representation[nodeIndex])
            parent = nodeIndex
            depth = depths[nodeIndex] + 1
            nodeIndex = size(visitedFromRoot)[1] + 1
        else
            nodeIndex = visitedFromRoot[nodeIndex, 2]
            if nodeIndex != 0
                visitedFromRoot[nodeIndex, 1] += 1
            end
        end
    end

    # Sets heights for every node
    heights = Array{Integer}(undef, lenRep)
    heights[end] = 0
    stack = [lenRep]
    currIndex = lenRep - 1

    while currIndex > 0

        if typeof(representation[currIndex]) <: TerminalNode
            heights[currIndex] = 0
        elseif typeof(representation[currIndex]) <: FunctionNode
            childrenIndexes = Array{Integer}(undef, 0)
            for i=1:getArity(representation[currIndex])
                push!(childrenIndexes, pop!(stack))
            end
            heights[currIndex] = maximum(heights[childrenIndexes]) + 1
        end

        push!(stack, currIndex)
        currIndex -= 1
    end

    # Chooses cross point
    indexes = Random.shuffle(rng, collect(1:lenRep))
    pointIndex = 0

    if maxDepth != -1
        if fixedDepth != -1
            indexes = union(filter(x->depths[x]==fixedDepth, indexes),
                            filter(x->heights[x]+fixedDepth<=maxDepth, indexes))
        elseif minDepth != -1
            indexes = union(filter(x->depths[x]>=minDepth, indexes),
                            filter(x->heights[x]+minDepth<=maxDepth, indexes))
        end
    end

    for i in indexes
        random = rand(rng)
        prob = typeof(representation[i]) <: FunctionNode ? 0.9 : 0.1

        if random < prob
            pointIndex = i
            break
        end
    end

    if pointIndex == 0
        pointIndex = rand(rng, indexes)
    end

    # Subtree obtention
    subtree = Array{Node}(undef, 0)
    visitedFromRoot = zeros(Integer, 1, 2)
    push!(subtree, representation[pointIndex])
    nodeIndex = pointIndex
    parent = 0

    while nodeIndex != 0

        if (nodeIndex-pointIndex+1) > size(visitedFromRoot)[1]
            visitedFromRoot = vcat(visitedFromRoot, zeros(Integer, 1, 2))
            visitedFromRoot[end, 2] = parent
            push!(subtree, representation[nodeIndex])
        end

        if visitedFromRoot[nodeIndex-pointIndex+1, 1] < getArity(representation[nodeIndex])
            parent = nodeIndex
            nodeIndex = size(visitedFromRoot)[1] + pointIndex
        else
            nodeIndex = visitedFromRoot[nodeIndex-pointIndex+1, 2]
            if nodeIndex != 0
                visitedFromRoot[nodeIndex-pointIndex+1, 1] += 1
            end
        end
    end

    return subtree, pointIndex, depths[pointIndex]
end # function

"""
    onePointCross(parent1::CGPGenotype, parent2::CGPGenotype,
                  gpExperimentInfo::CGPInfo, rng::Random.AbstractRNG)

Performs a crossover between two parent trees in which the crosspoints are at the
same depth.

`Self-provided Arguments` are provided by the library, so only `User Arguments` must be provided.

# Self-provided Arguments
- `parent1::CGPGenotype`: first parent for the crossover operation.
- `parent2::CGPGenotype`: first parent for the crossover operation.
- `gpExperimentInfo::CGPInfo`: information about the GP experiment.
- `rng::Random.AbstractRNG`: random number generator for random number consistency along an experiment.

# User Arguments
None

# Returns
An array that contains the two resulting children.
"""
function onePointCross(parent1::CGPGenotype, parent2::CGPGenotype,
                       gpExperimentInfo::CGPInfo, rng::Random.AbstractRNG)

    maxDepth = gpExperimentInfo._maxTreeDepth

    subtree1, crossPoint1, pointDepth1 = selectCrossPoint(parent1, rng)
    subtree2, crossPoint2, _ = selectCrossPoint(parent2, rng, fixedDepth=pointDepth1, maxDepth=maxDepth)
    child1Tree = mergeSubtrees(parent1, subtree2, crossPoint1, length(subtree1))
    child2Tree = mergeSubtrees(parent2, subtree1, crossPoint2, length(subtree2))

    return [CGPGenotype(child1Tree), CGPGenotype(child2Tree)]
end # function

"""
    oneChildSubtreeCross(parent1::CGPGenotype, parent2::CGPGenotype,
                         gpExperimentInfo::CGPInfo, rng::Random.AbstractRNG)

Performs the classical subtree cross between two parents but resulting in only
one child. Selects a random crosspoint in the first parent and then selects
another crosspoint in the second parent, knowing that it should not generate a
child that violates the maximum depth established.

`Self-provided Arguments` are provided by the library, so only `User Arguments` must be provided.

# Self-provided Arguments
- `parent1::CGPGenotype`: first parent for the crossover operation.
- `parent2::CGPGenotype`: first parent for the crossover operation.
- `gpExperimentInfo::CGPInfo`: information about the GP experiment.
- `rng::Random.AbstractRNG`: random number generator for random number consistency along an experiment.

# User Arguments
None

# Returns
An array that contains the resulting child.
"""
function oneChildSubtreeCross(parent1::CGPGenotype, parent2::CGPGenotype,
                              gpExperimentInfo::CGPInfo, rng::Random.AbstractRNG)

    maxDepth = gpExperimentInfo._maxTreeDepth

    subtree1, crossPoint1, pointDepth1 = selectCrossPoint(parent1, rng)
    subtree2, _, _ = selectCrossPoint(parent2, rng, minDepth=pointDepth1, maxDepth=maxDepth)
    child1Tree = mergeSubtrees(parent1, subtree2, crossPoint1, length(subtree1))

    return [CGPGenotype(child1Tree)]
end # function

"""
    subtreeCross(parent1::CGPGenotype, parent2::CGPGenotype,
                 gpExperimentInfo::CGPInfo, rng::Random.AbstractRNG)

Performs the classical subtree cross between two parents. Selects a random
crosspoint in the first parent and then selects another crosspoint in the second
parent, knowing that it should not generate a child that violates the maximum
depth established. Then, repeats the same operation but taking the second parent
as it was the first one.

`Self-provided Arguments` are provided by the library, so only `User Arguments` must be provided.

# Self-provided Arguments
- `parent1::CGPGenotype`: first parent for the crossover operation.
- `parent2::CGPGenotype`: first parent for the crossover operation.
- `gpExperimentInfo::CGPInfo`: information about the GP experiment.
- `rng::Random.AbstractRNG`: random number generator for random number consistency along an experiment.

# User Arguments
None

# Returns
An array that contains the two resulting children.
"""
function subtreeCross(parent1::CGPGenotype, parent2::CGPGenotype,
                      gpExperimentInfo::CGPInfo, rng::Random.AbstractRNG)

    child1 = oneChildSubtreeCross(parent1, parent2, gpExperimentInfo, rng)[1]
    child2 = oneChildSubtreeCross(parent2, parent1, gpExperimentInfo, rng)[1]

    return [child1, child2]
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__()

module RandomQuantum

using Compat, QuantumInfo

import Base.rand

export  # types
        GinibreEnsemble,
        FubiniStudyPureState,
        FubiniStudyMixedState,
        HilbertSchmidtMixedState,
        BuresMixedState,
        ClosedHaarEnsemble,
        OpenHaarEnsemble,
        RandomClosedEvolution,
        RandomOpenEvolution,
        # functions
        rand

"""
RandomQuantum.GinibreEnsemble(rows,cols)

Type corresponding to the Ginibre distribution
of complex matrices.
"""
type GinibreEnsemble
    rows::Int64
    cols::Int64
end

GinibreEnsemble(dim::Number) = GinibreEnsemble(dim,dim)

function rand(dist::GinibreEnsemble)
    return (randn(dist.rows,dist.cols)+im*randn(dist.rows,dist.cols))/sqrt(2)
end

"""
RandomQuantum.FubiniStudyPureState(dim)

Type corresponding to unitarily-invariant distribution
of pure states for some Hilbert space dimension.
"""
type FubiniStudyPureState
    dim::Int64
end

function rand(dist::FubiniStudyPureState)
    return normalize( vec(rand(GinibreEnsemble(dist.dim,1))) )
end

"""
RandomQuantum.FubiniStudyMixedState(dim, bath_dim)

Type corresponding to distribution of mixed states obtained by
tracing out part of a pure state obtained from a FubiniStudy
distribution. When `dim == bath_dim`, this is identical to the
Hilbert-Schmidt distribution of mixed states.
See, e.g., Zyczkowski et al., [J. Math. Phys. 52, 062201
(2011)](http://dx.doi.org/10.1063/1.3595693)
[arXiv:1010.3570](http://arxiv.org/abs/1010.3570).
"""
type FubiniStudyMixedState
    dim::Int64
    bath_dim::Int64
end

function rand(dist::FubiniStudyMixedState)
    ψ = rand(FubiniStudyPureState(dist.dim*dist.bath_dim))
    return partialtrace(projector(ψ),[dist.dim,dist.bath_dim],2)
end

"""
RandomQuantum.HilbertSchmidtState(dim)

Type corresponding to distribution of mixed states according to the
Hilbert-Schmidt measure (induced by the Frobenius distance). This is
identical to the mixed state distribution induced by tracing out half
of a pure state distributed according to the Fubini-Study metric.
See, e.g., Zyczkowski et al., [J. Math. Phys. 52, 062201
(2011)](http://dx.doi.org/10.1063/1.3595693)
[arXiv:1010.3570](http://arxiv.org/abs/1010.3570).
"""
type HilbertSchmidtMixedState
    dim::Int64
end

function rand(dist::HilbertSchmidtMixedState)
    X = rand(GinibreEnsemble(dist.dim))
    M = X*X'
    return M/trace(M)
end

"""
RandomQuantum.BuresMixedState(dim)

Type corresponding to distribution of mixed states according to the
Bures metric.
See, e.g., Zyczkowski et al., [J. Math. Phys. 52, 062201
(2011)](http://dx.doi.org/10.1063/1.3595693)
[arXiv:1010.3570](http://arxiv.org/abs/1010.3570).
"""
type BuresMixedState
    dim::Int64
end

function rand(dist::BuresMixedState)
    G = rand(GinibreEnsemble(dist.dim))
    U = rand(ClosedHaarEnsemble(dist.dim))
    H = (eye(dist.dim)+U)*G
    return H*H'/trace(H*H')
end

"""
RandomQuantum.ClosedHaarEnsemble(dim)

Type corresponding to the unitarily invariant distribution of unitary
transformations for some Hilbert space dimension.  See, e.g.,
Mezzadri, [Notices Amer Math Soc 54 4 592
(2007)](http://www.ams.org/notices/200705/fea-mezzadri-web.pdf).
"""
type ClosedHaarEnsemble
    dim::Int64
end

function rand(dist::ClosedHaarEnsemble)
    X = rand(GinibreEnsemble(dist.dim))
    Q,_ = qr(X,thin=false)
    d = diag(Q)
    d = d./abs(d)
    Q = Q./d
    return Q
end

"""
RandomQuantum.OpenHaarEnsemble(dim, bath_dim)

Type corresponding to the unitarily invariant distribution of
unitary transformations for system and bath, such that the bath
(initially in the ground state) is traced out.
See, e.g., Bruzda et al., [Phys. Lett. A 373, 320-324
(2009)](http://dx.doi.org/10.1016/j.physleta.2008.11.043),
[arXiv:0804.2361](http://arxiv.org/abs/0804.2361).
"""
type OpenHaarEnsemble
    dim::Int64
    bath_dim::Int64
end

function rand(dist::OpenHaarEnsemble)
    X = rand(GinibreEnsemble(dist.dim^2,dist.bath_dim))
    W = X*X'
    W = W/trace(W)
    #Y = sqrtm(trace(W,[dist.dim,dist.dim],1))
    #IY = kron(eye(dist.dim),Y)
    Y = sqrtm(partialtrace(W,[dist.dim,dist.dim],2))
    IY = kron(Y,eye(dist.dim))
    R = IY\W/IY
    return choi2liou(R)/dist.dim
end

"""
RandomQuantum.GUE(dim)

Type corresponding to the Gaussian Unitary Ensemble of complex
matrices.
"""
type GUE
    size::Int64
end

function rand(dist::GUE)
    X = rand(GinibreEnsemble(dist.size))
    return (X+X')/2
end

"""
RandomQuantum.RandomClosedEvolution(dim, α)

Type corresponding to the integrated evolution of random Hamiltonians
(with unitarily invariant distribution) on some Hilbert space.
"""
type RandomClosedEvolution
    dim::Int64
    α::Float64
end

# TODO: Instead of matrix exponentiation, could this be done by
#       sampling eigenvalues, exponetiating them, and then multiplying
#       by a Haar random unitary? If so, that should be cheaper
function rand(dist::RandomClosedEvolution)
    evals,U = eig(rand(GUE(dist.dim)))
    return U*diagm(exp(-1im*dist.α*evals))*U'
end

"""
RandomQuantum.RandomOpenEvolution(dim, bath_dim, α)

Type corresponding to the distribution of completely positive trace
preserving maps obtained from the integrated evolution of random
Hamiltonians (with unitarily invariant distribution) on a
dimensional Hilbert space, such that bath (initially in the ground
state) is traced out after the evolution.
"""
type RandomOpenEvolution
    dim::Int64
    bath_dim::Int64
    α::Float64
end

function rand(dist::RandomOpenEvolution)
    d = dist.dim
    de = dist.bath_dim
    ru = rand(RandomClosedEvolution(d*de,dist.α)) * kron(eye(d),ket(0,de))
    k = Matrix{Complex128}[]
    for ii=1:de
        push!(k, kron(eye(d),bra(ii-1,de)) * ru )
    end
    return kraus2liou(k)
end

end