text:
module TestMetadataUtils

using MLJBase, Test

@mlj_model mutable struct FooRegressor <: Deterministic
    a::Int = 0::(_ ≥ 0)
    b
end
metadata_pkg(FooRegressor,
    name="FooRegressor",
    uuid="10745b16-79ce-11e8-11f9-7d13ad32a3b2",
    url="http://existentialcomics.com/",
    julia=true,
    license="MIT",
    is_wrapper=false
    )
metadata_model(FooRegressor,
    input=MLJBase.Table(Continuous),
    target=AbstractVector{Continuous},
    descr="La di da")

infos = info_dict(FooRegressor)

@test infos[:input_scitype] == MLJBase.Table(Continuous)
@test infos[:target_scitype] == AbstractVector{Continuous}
@test infos[:is_pure_julia]
@test !infos[:is_wrapper]
@test infos[:docstring] == raw"""La di da
    → based on [FooRegressor](http://existentialcomics.com/).
    → do `@load FooRegressor pkg="FooRegressor"` to use the model.
    → do `?FooRegressor` for documentation."""
@test infos[:name] == "FooRegressor"

@test infos[:is_supervised]
@test infos[:prediction_type] == :deterministic

end
true


text:
using Test

tests = [
    "timearray",
    "split",
    "apply",
    "broadcast",
    "combine",
    "modify",
    "meta",
    "readwrite",
    "timeseriesrc",
    "basemisc",
    "tables",
    "plotrecipes",
]

@testset "TimeSeries" begin
    @info("Running tests:")

    for test ∈ tests
        @info("\t* $test ...")
        include("$test.jl")
    end
end

ext: jl


text:
# This file is part of the TaylorIntegration.jl package; MIT licensed

using TaylorIntegration
using Test
using LinearAlgebra: Transpose, norm
using Logging
import Logging: Warn

@testset "Testing `many_ode.jl`" begin

    local _order = 28
    local _abstol = 1.0E-20
    local tT = Taylor1(_order)

    max_iters_reached() = "Maximum number of integration steps reached; exiting.\n"

    @testset "Tests: dot{x}=x.^2, x(0) = [3, 1]" begin
        function eqs_mov!(Dx, x, p, t)
            for i in eachindex(x)
                Dx[i] = x[i]^2
            end
            nothing
        end
        exactsol(t, x0) = x0/(1.0-x0*t)
        t0 = 0.0
        q0 = [3, 1]
        q0T = Taylor1.(Float64.(q0), _order)
        xdotT = Array{Taylor1{Float64}}(undef, length(q0))
        xaux = Array{Taylor1{Float64}}(undef, length(q0))
        tT[1] = t0
        TaylorIntegration.jetcoeffs!(eqs_mov!, tT, q0T, xdotT, xaux, nothing)
        @test q0T[1].coeffs[end] == 3.0^(_order+1)
        @test q0T[2].coeffs[end] == 1.0
        δt = (_abstol/q0T[1].coeffs[end-1])^inv(_order-1)
        @test TaylorIntegration.stepsize(q0T, _abstol) == δt

        tv, xv = (@test_logs (Warn, max_iters_reached()) taylorinteg(
            eqs_mov!, q0, 0.0, 0.5, _order, _abstol, nothing))
        @test length(tv) == 501
        @test isa(xv, SubArray)
        @test xv[1,:] == q0
        @test tv[end] < 1/3

        trange = 0.0:1/8:1.0
        xv = (@test_logs (Warn, max_iters_reached()) taylorinteg(
            eqs_mov!, q0, trange, _order, _abstol))
        @test size(xv) == (9,2)
        @test q0 == [3.0, 1.0]
        @test typeof(xv) == Transpose{Float64, Array{Float64,2}}
        @test xv[1,1:end] == q0
        @test (isnan(xv[4,1]) && isnan(xv[4,2]))
        @test (isnan(xv[end,1]) && isnan(xv[end,2]))
        @test abs(xv[3,2] - 4/3) ≤ eps(4/3)
        @test abs(xv[2,1] - 4.8) ≤ eps(4.8)

        tarray = vec(trange)
        xv2 = (@test_logs (Warn, max_iters_reached()) taylorinteg(
            eqs_mov!, q0, tarray, _order, _abstol, nothing))
        @test xv[1:3,:] == xv2[1:3,:]
        @test xv2[1:3,:] ≈ xv[1:3,:] atol=eps() rtol=0.0
        @test size(xv2) == (9,2)
        @test q0 == [3.0, 1.0]
        @test typeof(xv2) == Transpose{Float64, Array{Float64,2}}
        @test xv2[1,1:end] == q0
        @test (isnan(xv2[4,1]) && isnan(xv2[4,2]))
        @test (isnan(xv2[end,1]) && isnan(xv2[end,2]))
        @test abs(xv2[3,2] - 4/3) ≤ eps(4/3)
        @test abs(xv2[2,1] - 4.8) ≤ eps(4.8)

        # Output includes Taylor polynomial solution
        tv, xv, polynV = (@test_logs (Warn, max_iters_reached()) taylorinteg(
            eqs_mov!, q0, 0, 0.5, _order, _abstol, Val(true), nothing, maxsteps=2))
        @test size(polynV) == (3, 2)
        @test xv[1,:] == q0
        @test polynV[1,:] == Taylor1.(q0, _order)
        @test xv[2,:] == evaluate.(polynV[2, :], tv[2]-tv[1])
        @test polynV[2,2] == Taylor1(ones(_order+1))
    end

    @testset "Tests: dot{x}=x.^2, x(0) = [3, 3]" begin
        function eqs_mov!(Dx, x, p, t)
            for i in eachindex(x)
                Dx[i] = x[i]^2
            end
            nothing
        end
        exactsol(t, x0) = x0/(1.0-x0*t)

        q0 = [3.0, 3.0]
        tmax = 0.3
        tv, xv = (@test_logs min_level=Logging.Warn taylorinteg(
            eqs_mov!, q0, 0, tmax, _order, _abstol, nothing))
        @test length(tv) < 501
        @test length(xv[:,1]) < 501
        @test length(xv[:,2]) < 501
        @test length(tv) == 14
        @test length(xv[:,1]) == 14
        @test length(xv[:,2]) == 14
        @test xv[1,1:end] == q0
        @test tv[end] < 1/3
        @test tv[end] == tmax
        @test xv[end,1] == xv[end,2]
        @test abs(xv[end,1]-exactsol(tv[end], xv[1,1])) < 5e-14
        @test abs(xv[end,2]-exactsol(tv[end], xv[1,2])) < 5e-14

        tmax = 0.33
        tv, xv = (@test_logs min_level=Logging.Warn taylorinteg(
            eqs_mov!, [3, 3], 0.0, tmax, _order, _abstol))
        @test length(tv) < 501
        @test length(xv[:,1]) < 501
        @test length(xv[:,2]) < 501
        @test length(tv) == 28
        @test length(xv[:,1]) == 28
        @test length(xv[:,2]) == 28
        @test xv[1,1:end] == q0
        @test tv[end] < 1/3
        @test tv[end] == tmax
        @test xv[end,1] == xv[end,2]
        @test abs(xv[end,1]-exactsol(tv[end], xv[1,1])) < 1.0e-11
        @test abs(xv[end,2]-exactsol(tv[end], xv[1,2])) < 1.0e-11

        # Output includes Taylor polynomial solution
        tv, xv, polynV = taylorinteg(eqs_mov!, q0, 0, tmax, _order, _abstol, Val(true), nothing)
        @test size(polynV) == size(xv)
        @test polynV[1,:] == Taylor1.(q0, _order)
        @test xv[end,1] == evaluate.(polynV[end, 1], tv[end]-tv[end-1])
        @test polynV[:,1] == polynV[:,2]
    end

    @testset "Test non-autonomous ODE (2): dot{x}=cos(t)" begin
        function f!(Dx, x, p, t)
            Dx[1] = one(t)
            Dx[2] = cos(t)
            nothing
        end
        t0 = 0//1
        tmax = 10.25*(2pi)
        abstol = 1e-20
        order = 25
        x0 = [t0, 0.0] #initial conditions such that x(t)=sin(t)

        tv, xv = (@test_logs min_level=Logging.Warn taylorinteg(
            f!, x0, t0, tmax, order, abstol))
        @test length(tv) < 501
        @test length(xv[:,1]) < 501
        @test length(xv[:,2]) < 501
        @test xv[1,1:end] == x0
        @test tv[1] < tv[end]
        @test tv[end] == xv[end,1]
        @test abs(sin(tmax)-xv[end,2]) < 1e-14

        # Backward integration
        tb, xb = (@test_logs min_level=Logging.Warn taylorinteg(
            f!, [tmax, sin(tmax)], tmax, t0, order, abstol))
        @test length(tb) < 501
        @test length(xb[:,1]) < 501
        @test length(xb[:,2]) < 501
        @test tb[1] > tb[end]
        @test xb[1,1:end] == [tmax, sin(tmax)]
        @test tb[end] == xb[end,1]
        @test abs(sin(t0)-xb[end,2]) < 5e-14

        # Tests with a range, for comparison with backward integration
        tmax = 15*(2pi)
        Δt = (tmax-t0)/1024
        tspan = t0:Δt:tmax
        xv = (@test_logs min_level=Logging.Warn taylorinteg(
            f!, x0, tspan, order, abstol, nothing))
        @test xv[1,1:end] == x0
        @test tmax == xv[end,1]
        @test abs(sin(tmax)-xv[end,2]) < 1e-14

        # Backward integration
        xback = (@test_logs min_level=Logging.Warn taylorinteg(
            f!, xv[end, :], reverse(tspan), order, abstol, nothing))
        @test xback[1,:] == xv[end, :]
        @test abs(xback[end,1]-x0[1]) < 5.0e-14
        @test abs(xback[end,2]-x0[2]) < 5.0e-14
        @test abs(sin(t0)-xback[end,2]) < 5.0e-14
        @test norm((xv[:,:]-xback[end:-1:1,:]), Inf) < 5.0e-14

        # Tests if trange is properly sorted
        @test_throws AssertionError taylorinteg(f!, x0, rand(t0:Δt:tmax, 100), order, abstol)
    end

    @testset "Falling ball (stepsize)" begin
        function fallball!(dx, x, p, t)
            dx[1] = x[2]
            dx[2] = -one(x[1])
            nothing
        end
        exactsol(t, x0, v0) = x0 + v0*t - 0.5*t^2
        t0 = 0.0
        tmax = 5.0
        abstol = 1e-20
        order = 10
        x0 = [10.0, 0.0] #initial conditions such that x(t)=sin(t)
        tv, xv = (@test_logs min_level=Logging.Warn taylorinteg(
            fallball!, x0, t0, tmax, order, abstol))
        @test length(tv) < 501
        @test length(xv[:,1]) < 501
        @test exactsol.(tv, x0[1], x0[2]) ≈ xv[:,1]
    end

end


text:
module ParseTorrent

using Bencode
using SHA

export parsetorrent

Bytes = Vector{UInt8}

function parsetorrent(b::Bytes)::Dict
    torrent = bdecode(b)
    infohash = sha1(bencode(torrent["info"]))
    torrent["infohash"] = bytes2hex(infohash)
    pieces = torrent["info"]["pieces"]
    torrent["info"]["pieces"] = pieces2hashes(Vector{UInt8}(pieces))
    torrent
end

function pieces2hashes(pieces::Bytes)::Vector{String}
    len = length(pieces)
    pieces = reshape(pieces, (20, length(pieces) ÷ 20))
    map(bytes2hex, eachcol(pieces))
end

end # module


text:
# Description: startup function for solving an equation

import PDESolver.solvePDE
#=
"""
  This function invokes the solver for the Euler equations, using the
  specified input file

  Inputs:
    input_file: a string containing the path to an input file, or just a file
                name.  If it is just a file name, it is taken to be in the
                users pwd.

  Outputs:
    mesh: an AbstractMesh object
    sbp: the SBP operator used in solving the equation
    eqn: the AbstractSolutionData object used to solve the equation
    opts: the options dictonary

"""
function run_euler(input_file::AbstractString)

  mesh, sbp, eqn, opts, pmesh = createObjects(input_file)
  solve_euler(mesh, sbp, eqn, opts, pmesh)

  return mesh, sbp, eqn, opts
end
=#

"""
  This function creates and initializes the mesh, sbp, eqn, and opts objects

  Inputs:
    opts: options dictionary

  Outputs:
    mesh: an AbstractMesh.  The concrete type is determined by the options
          dictionary
    sbp: an AbstractOperator.  The concrete type is determined by the options
         dictionary
    eqn: an EulerData object
    opts: the options dictionary
    pmesh: mesh used for preconditioning, can be same object as mesh
"""
function createObjects(opts::Dict)

  Tdim = opts["dimensions"]
  dofpernode = Tdim + 2

  sbp, mesh, pmesh, Tsol, Tres, Tmsh, mesh_time = createMeshAndOperator(opts, dofpernode)

  # TODO: input argument for dofpernode

  # create euler equation
  var_type = opts["variable_type"]
  eqn = EulerData_{Tsol, Tres, Tdim, Tmsh, var_type}(mesh, sbp, opts)

  # initialize physics module and populate any fields in mesh and eqn that
  # depend on the physics module
  init(mesh, sbp, eqn, opts, pmesh)

#  if opts["need_adjoint"]
#    init_revm(mesh, sbp, eqn, opts, pmesh)
#  end

  return mesh, sbp, eqn, opts, pmesh
end

"""
  Constructs a the EulerData object given a mesh, sbp, and options dictionary.

  Used for submesh solves.

  **Inputs**

   * mesh
   * sbp
   * opts

  **Outputs**

   * mesh
   * sbp
   * eqn
   * opts
   * pmesh: currently, always the same as mesh
"""
function createObjects(mesh::AbstractMesh, sbp::AbstractOperator, opts::Dict)

  var_type = opts["variable_type"]

  Tdim = mesh.dim
  Tmsh, Tsbp, Tsol, Tres = getDataTypes(opts)

  eqn = EulerData_{Tsol, Tres, Tdim, Tmsh, var_type}(mesh, sbp, opts)

  init(mesh, sbp, eqn, opts, mesh)

  return mesh, sbp, eqn, opts, mesh
end

"""
  Given fully initialized mesh, sbp, eqn, opts, this function solves
  the Euler equations.  The 4 object should be obtained from createObjects().

  Specifically, it applies an initial condition and invokes a nonlinear
  solver according to the options dictionary.

  Inputs:
    mesh: an AbstractMesh
    sbp: an AbstractOperator
    eqn: an AbstractEulerData
    opts: the options dictionary.  This must be the options dictionary returned
          by createObjects().  Changing values in the options dictionary after
          calling createObjects() results in undefined behavior.
    pmesh: mesh used for preconditioning, can be same object as mesh.
           default value of mesh

"""
function solvePDE(mesh::AbstractMesh, sbp::AbstractOperator, eqn::AbstractEulerData, opts::Dict, pmesh::AbstractMesh=mesh)
  #delta_t = opts["delta_t"]   # delta_t: timestep for RK

  myrank = mesh.myrank
  # flag determines whether to calculate u, dR/du, or dR/dx (1, 2, or 3)
  t_max = opts["t_max"]
  flag = opts["run_type"]
  var_type = opts["variable_type"]

  # timestepping parameters
  order = opts["order"]       # order of accuracy

  # zero some arrays out in case the user forgot

  fill!(eqn.res, 0.0)
  fill!(eqn.res_vec, 0.0)
  res_vec = eqn.res_vec
  q_vec = eqn.q_vec       # solution at current timestep

  # calculate residual of some other function for res_reltol0
  # TODO: add a boolean options here?
  Relfunc_name = opts["Relfunc_name"]
  if haskey(ICDict, Relfunc_name)
    @mpi_master println("\ncalculating residual for relative residual tolerance")
    Relfunc = ICDict[Relfunc_name]
    @mpi_master println("Relfunc = ", Relfunc)
    Relfunc(mesh, sbp, eqn, opts, q_vec)

    if var_type == :entropy
      @mpi_master println("converting to entropy variables")
      for i=1:mesh.numDofPerNode:mesh.numDof
        q_view = sview(q_vec, i:(i+mesh.numDofPerNode-1))
        convertFromNaturalToWorkingVars(eqn.params, q_view, q_view)
      end
    end
    tmp = physicsRhs(mesh, sbp, eqn, opts, eqn.res_vec, (evalResidual,))
    res_real = real(eqn.res_vec)
    opts["res_reltol0"] = tmp
    println("res_reltol0 = ", tmp)

    saveSolutionToMesh(mesh, res_real)
    writeVisFiles(mesh, "solution_relfunc")
  end

  # populate u0 with initial condition
  @mpi_master println("\nEvaluating initial condition")
  ICfunc_name = opts["IC_name"]
  ICfunc = ICDict[ICfunc_name]
  @mpi_master println("ICfunc = ", ICfunc)
  ICfunc(mesh, sbp, eqn, opts, q_vec)

  if var_type == :entropy
    for i=1:mesh.numDofPerNode:mesh.numDof
      q_view = sview(q_vec, i:(i+mesh.numDofPerNode-1))
      convertFromNaturalToWorkingVars(eqn.params, q_view, q_view)
    end
  end
  println("\n\nFinished evaluating initial condition")

  #TODO: TESTING
  if opts["freeze_viscosity"]
    freezeViscosity(mesh, sbp, eqn, opts)
  end

  if opts["calc_error"]
    @mpi_master println("\ncalculating error of file ",
                       opts["calc_error_infname"],
                      " compared to initial condition")

    # read in this processors portion of the solution
    vals = copy(eqn.q_vec)
    readSolutionFiles(mesh, sbp, eqn, opts, opts["calc_error_infname"], vals)

    err_vec = abs.(vals - eqn.q_vec)
    err = calcNorm(eqn, err_vec)

    # calculate avg mesh size
    h_avg = calcMeshH(mesh, sbp, eqn, opts)

    @mpi_master begin
      outname = opts["calc_error_outfname"]
      println("printint err = ", err, " to file ", outname)
      f = open(outname, "w")
      println(f, err, " ", h_avg)
      close(f)
    end

    # write visualization
    saveSolutionToMesh(mesh, vec(err_vec))
    writeVisFiles(mesh, "solution_error")
  end

  if opts["calc_trunc_error"]  # calculate truncation error
    @mpi_master println("\nCalculating residual for truncation error")
    tmp = physicsRhs(mesh, sbp, eqn, opts, eqn.res_vec, (evalResidual,))

    @mpi_master begin
      f = open("error_trunc.dat", "w")
      println(f, tmp)
      close(f)
    end
  end

  if opts["perturb_ic"]
    @mpi_master println("\nPerturbing initial condition")
    perturb_mag = opts["perturb_mag"]
    for i=1:mesh.numDof
      q_vec[i] += perturb_mag*rand()
    end
  end

  res_vec_exact = deepcopy(q_vec)

  #writeSolutionFiles(mesh, sbp, eqn, opts, "IC")
  saveSolutionToMesh(mesh, q_vec)

  writeVisFiles(mesh, "solution_ic")
  if opts["calc_dt"]
    wave_speed = EulerEquationMod.calcMaxWaveSpeed(mesh, sbp, eqn, opts)
    @mpi_master println("max wave speed = ", wave_speed)
    @mpi_master println("min element size = ", mesh.min_el_size)
    delta_t = opts["CFL"]*mesh.min_el_size/wave_speed
    @mpi_master println("for a CFL of ", opts["CFL"], " delta_t = ", delta_t)
    opts["delta_t"] = delta_t
  end

  # this must be last, because all previous calculations need to use the
  # original initial condition
  if opts["is_restart"]
    loadRestartState(mesh, sbp, eqn, opts, pmesh)
  end

  call_nlsolver(mesh, sbp, eqn, opts, pmesh)
  postproc(mesh, sbp, eqn, opts)

  cleanup(eqn, opts)

  MPI.Barrier(mesh.comm)

  return mesh, sbp, eqn, opts
end  # end function
#runtest(1)

"""
  This function does post processing, if requested by the input options.
  Typical post processing includes calculation of errors, norms of important
  quantities, writing of files. etc.

  Inputs:
    mesh
    sbp
    eqn
    opts
"""
function postproc(mesh, sbp, eqn, opts)

  ##### Do postprocessing ######
  myrank = mesh.myrank
  @mpi_master println("\nDoing postprocessing")

  Tsol = eltype(eqn.q)
  Tres = eltype(eqn.res)
  Tmsh = eltype(mesh.dxidx)
  myrank = mesh.myrank

  if opts["do_postproc"] && opts["solve"]
    exfname = opts["exact_soln_func"]
    if haskey(ICDict, exfname)
      exfunc = ICDict[exfname]
      q_exact = zeros(Tsol, mesh.numDof)
      exfunc(mesh, sbp, eqn, opts, q_exact)
#    if var_type == :entropy
#      println("converting to entropy variables")
#      for i=1:mesh.numDofPerNode:mesh.numDof
#        q_view = sview(q_vec, i:(i+mesh.numDofPerNode-1))
#        convertFromNaturalToWorkingVars(eqn.params, q_view, q_view)
#      end
#    end

      q_diff = eqn.q_vec - q_exact
      saveSolutionToMesh(mesh, abs.(real(q_diff)))
      writeVisFiles(mesh, "solution_error")

      diff_norm = calcNorm(eqn, q_diff)
      @mpi_master println("solution error norm = ", diff_norm)
      # TODO: make this mesh.min_el_size?
      h_avg = calcMeshH(mesh, sbp, eqn, opts)

      # print to file
      @mpi_master begin
        outname = opts["calc_error_outfname"]
        f = open(outname, "w")
        println(f, diff_norm, " ", h_avg)
        close(f)
      end
    end  # end if haskey(ICname)
  end  # end if do_postproc

  return nothing
end


text:
# Chemfiles.jl, a modern library for chemistry file reading and writing
# Copyright (C) Guillaume Fraux and contributors -- BSD license

using Compat: contains
export id, residue_for_atom, atoms, contains

__ptr(residue::Residue) = __ptr(residue.__handle)
__const_ptr(residue::Residue) = __const_ptr(residue.__handle)

"""
Create a new residue with the given `name` and optional residue identifier `id`.
"""
function Residue(name::String, id=nothing)
    if id === nothing
        ptr = @__check_ptr(lib.chfl_residue(pointer(name)))
    else
        ptr = @__check_ptr(lib.chfl_residue_with_id(pointer(name), Int64(id)))
    end
    return Residue(CxxPointer(ptr, is_const=false))
end

"""
Get a copy of the residue at `index` from a `topology`.

The residue index in the topology is not always the same as the residue
identifier.
"""
function Residue(topology::Topology, index::Integer)
    ptr = @__check_ptr(lib.chfl_residue_from_topology(__ptr(topology), UInt64(index)))
    residue = Residue(CxxPointer(ptr, is_const=true))
    copy = deepcopy(residue)
    finalize(residue)
    return copy
end

"""
Get a copy of the residue containing the atom at `index` in the `topology`.

This function will return `nothing` if the atom is not in a residue, or if the
`index` is bigger than the number of atoms in the topology.
"""
function residue_for_atom(topology::Topology, index::Integer)
    ptr = lib.chfl_residue_for_atom(__ptr(topology), UInt64(index))
    if Int(ptr) == 0
        return nothing
    else
        residue = Residue(CxxPointer(ptr, is_const=true))
        copy = deepcopy(residue)
        finalize(residue)
        return copy
    end
end

"""
Get the name of a `residue`.
"""
function name(residue::Residue)
    return __call_with_growing_buffer(
        (buffer, size) -> __check(lib.chfl_residue_name(
            __const_ptr(residue), buffer, size
        ))
    )
end

"""
Get the identifier of a `residue` in the initial topology.
"""
function id(residue::Residue)
    resid = Ref{Int64}(0)
    __check(lib.chfl_residue_id(__const_ptr(residue), resid))
    return resid[]
end

"""
Get the number of atoms in a `residue`.
"""
function Base.size(residue::Residue)
    count = Ref{UInt64}(0)
    __check(lib.chfl_residue_atoms_count(__const_ptr(residue), count))
    return Int(count[])
end

"""
Add the atom at the given `index` in the `residue`.
"""
function add_atom!(residue::Residue, index::Integer)
    __check(lib.chfl_residue_add_atom(__ptr(residue), UInt64(index)))
    return nothing
end

"""
Check if the atom at the given `index` is in the `residue`.
"""
function Compat.contains(residue::Residue, index::Integer)
    result = Ref{UInt8}(0)
    __check(lib.chfl_residue_contains(__const_ptr(residue), UInt64(index), result))
    return convert(Bool, result[])
end

"""
Get the atoms in a ``residue``. This function returns a list of indexes.
"""
function atoms(residue::Residue)
    count = UInt64(size(residue))
    result = Array{UInt64}(undef, count)
    __check(lib.chfl_residue_atoms(__const_ptr(residue), pointer(result), count))
    return result
end

"""
Set a named property for the given residue.
"""
function set_property!(residue::Residue, name::String, value::PropertyValue)
    property = Property(value)
    __check(lib.chfl_residue_set_property(
        __ptr(residue), pointer(name), __const_ptr(property)
    ))
    return nothing
end

"""
Get a named property for the given residue.
"""
function property(residue::Residue, name::String)::PropertyValue
    ptr = lib.chfl_residue_get_property(__const_ptr(residue), pointer(name))
    return extract(Property(CxxPointer(ptr, is_const=false)))
end

"""
Get the number of properties associated with a residue.
"""
function properties_count(residue::Residue)
    count = Ref{UInt64}(0)
    __check(lib.chfl_residue_properties_count(__const_ptr(residue), count))
    return Int(count[])
end

"""
Get the names of all properties associated with a residue.
"""
function list_properties(residue::Residue)
    count = UInt64(properties_count(residue))
    names = Array{Ptr{UInt8}}(undef, count)
    __check(lib.chfl_residue_list_properties(__const_ptr(residue), pointer(names), count))
    return map(unsafe_string, names)
end

"""
Make a deep copy of a `residue`.
"""
function Base.deepcopy(residue::Residue)
    ptr = lib.chfl_residue_copy(__const_ptr(residue))
    return Residue(CxxPointer(ptr, is_const=false))
end


text:
using DataFrames
using DataFramesMeta
using LinearAlgebra
using NearestNeighbors
using Statistics
using StatsBase: countmap

mutable struct BmmData{L}

    # Static data
    ## Segmentation data
    x::DataFrame;
    position_data::Matrix{Float64};
    composition_data::Vector{<:Union{Int, Missing}};
    confidence::Vector{Float64}

    cluster_per_molecule::Vector{Int}
    segment_per_molecule::Vector{Int}
    nuclei_prob_per_molecule::Vector{Float64}

    ## MRF
    adjacent_points::Array{Vector{Int}, 1};
    adjacent_weights::Array{Vector{Float64}, 1};
    real_edge_weight::Float64;
    position_knn_tree::KDTree;
    knn_neighbors::Array{Vector{Int}, 1};

    # Distribution-related
    components::Array{Component{L}, 1};
    distribution_sampler::Component{L};
    assignment::Vector{Int};
    max_component_guid::Int;

    noise_density::Float64;

    cluster_per_cell::Vector{Int};

    # Prior segmentation

    n_molecules_per_segment::Vector{Int}
    main_segment_per_cell::Vector{Int}

    # Utils
    tracer::Dict{Symbol, Any}
    misc::Dict{Symbol, Any}
    center_sample_cache::Vector{Int}

    # Parameters
    prior_seg_confidence::Float64
    cluster_penalty_mult::Float64
    use_gene_smoothing::Bool
    min_nuclei_frac::Float64
    mrf_strength::Float64

    """
    ...
    # Arguments
    - `components::Array{Component, 1}`:
    - `x::DataFrame`:
    - `adjacent_points::Array{Array{Int, 1}, 1}`:
    - `adjacent_weights::Array{Array{Float64, 1}, 1}`: edge weights, used for smoothness penalty
    - `real_edge_weight::Float64`: weight of an edge for "average" real point
    - `distribution_sampler::Component`:
    - `assignment::Array{Int, 1}`:
    """
    function BmmData(components::Array{Component{N}, 1}, x::DataFrame, adjacent_points::Array{Vector{Int}, 1}, adjacent_weights::Array{Vector{Float64}, 1},
                     assignment::Vector{Int}, distribution_sampler::Component{N}; real_edge_weight::Float64=1.0, k_neighbors::Int=20,
                     cluster_penalty_mult::Float64=0.25, use_gene_smoothing::Bool=true, prior_seg_confidence::Float64=0.5,
                     min_nuclei_frac::Float64=0.1, mrf_strength::Float64=0.1, na_genes::Vector{Int}=Int[]) where N
        @assert maximum(assignment) <= length(components)
        @assert minimum(assignment) >= 0
        @assert length(assignment) == size(x, 1)

        if !all(s in propertynames(x) for s in [:x, :y, :gene])
            error("`x` data frame must have columns 'x', 'y' and 'gene'")
        end

        p_data = position_data(x)
        @assert size(p_data, 1) == N

        position_knn_tree = KDTree(p_data)
        knn_neighbors = knn(position_knn_tree, p_data, k_neighbors)[1]

        comp_data = [ifelse(g in na_genes, missing, g) for g in composition_data(x)]
        n_genes = maximum(skipmissing(composition_data(x)))

        x = deepcopy(x)
        if !(:confidence in propertynames(x))
            x[!, :confidence] .= 0.95
        end

        nuclei_probs = :nuclei_probs in propertynames(x) ? x[!, :nuclei_probs] : Float64[]
        cluster_per_molecule = :cluster in propertynames(x) ? x.cluster : Int[]
        self = new{N}(
            x, p_data, comp_data, confidence(x), cluster_per_molecule, Int[], nuclei_probs,
            adjacent_points, adjacent_weights, real_edge_weight, position_knn_tree, knn_neighbors,
            components, deepcopy(distribution_sampler), assignment, length(components), 0.0,
            Int[],
            Int[], Int[], # prior segmentation info
            Dict{Symbol, Any}(), Dict{Symbol, Any}(), Int[],
            prior_seg_confidence, cluster_penalty_mult, use_gene_smoothing, min_nuclei_frac, mrf_strength
        )

        for c in self.components
            c.n_samples = 0
        end

        for c_id in assignment[assignment .> 0]
            self.components[c_id].n_samples += 1
        end

        # Resolve component guids
        guids = [c.guid for c in self.components]

        if maximum(guids) <= 0
            self.max_component_guid = length(self.components)
            for (i,c) in enumerate(self.components)
                c.guid = i
            end
        else
            if minimum(guids) <= 0
                error("Either all or no guids can be <= 0")
            end

            self.max_component_guid = maximum(guids)
        end

        if :prior_segmentation in propertynames(x)
            self.segment_per_molecule = deepcopy(x.prior_segmentation);
            self.n_molecules_per_segment = count_array(self.segment_per_molecule, drop_zero=true);
            update_n_mols_per_segment!(self);
        end

        return self
    end
end

function position_data(df::AbstractDataFrame)::Matrix{Float64}
    if (:z in propertynames(df))
        return copy(Matrix{Float64}(df[:, [:x, :y, :z]])')
    end

    return copy(Matrix{Float64}(df[:, [:x, :y]])')
end
position_data(data::BmmData)::Matrix{Float64} = data.position_data
composition_data(df::AbstractDataFrame)::Union{Vector{Int}, Vector{Union{Missing, Int}}} = df.gene
composition_data(data::BmmData) = data.composition_data
confidence(df::AbstractDataFrame)::Vector{Float64} = df.confidence
confidence(data::BmmData)::Vector{Float64} = data.confidence

num_of_molecules_per_cell(data::BmmData) = count_array(data.assignment, max_value=length(data.components), drop_zero=true)

function assign!(data::BmmData, point_ind::Int, component_id::Int)
    old_id = data.assignment[point_ind]
    if old_id == component_id
        return
    end

    @assert component_id <= length(data.components) "Too large component id: $component_id, maximum available: $(length(data.components))"

    segment_id = isempty(data.segment_per_molecule) ? 0 : data.segment_per_molecule[point_ind]
    if segment_id > 0
        if component_id > 0
            data.components[component_id].n_molecules_per_segment[segment_id] = get(data.components[component_id].n_molecules_per_segment, segment_id, 0) + 1
        end

        if old_id > 0
            data.components[old_id].n_molecules_per_segment[segment_id] -= 1
        end
    end

    data.assignment[point_ind] = component_id
end

function estimate_assignment_by_history(data::BmmData)
    # TODO: it doesn't guarantee connectivity. Can try to run deterministic EM, or use some better estimate here
    if !(:assignment_history in keys(data.tracer)) || (length(data.tracer[:assignment_history]) == 0)
        @warn "Data has no saved history of assignments. Fall back to the basic assignment"
        return data.assignment, ones(length(data.assignment)) / 2
    end

    guid_map = Dict(c.guid => i for (i,c) in enumerate(data.components))
    current_guids = Set(vcat(collect(keys(guid_map)), [0]));
    assignment_mat = hcat(data.tracer[:assignment_history]...);

    reassignment = mapslices(assignment_mat, dims=2) do row
        c_row = row[in.(row, Ref(current_guids))]
        if length(c_row) == 0
            return 0
        end

        c_counts = countmap(c_row);
        count_vals = collect(values(c_counts))
        return maximum(collect(keys(c_counts))[count_vals .== maximum(count_vals)])
    end

    return get.(Ref(guid_map), vec(reassignment), 0), vec(mean(assignment_mat .== reassignment, dims=2))
end

function get_cell_qc_df(segmented_df::DataFrame, cell_assignment::Vector{Int}=segmented_df.cell; sigdigits::Int=4, max_cell::Int=maximum(cell_assignment), dapi_arr::Union{Matrix{<:Real}, Nothing}=nothing)
    seg_df_per_cell = split(segmented_df, cell_assignment .+ 1; max_factor=max_cell+1)[2:end];
    pos_data_per_cell = [position_data(df)[1:2,:] for df in seg_df_per_cell];

    df = DataFrame(:n_transcripts => size.(pos_data_per_cell, 2));
    large_cell_mask = (df.n_transcripts .> 2)

    df[!,:density] = fill(NaN, size(df, 1))
    df[!,:elongation] = fill(NaN, size(df, 1))
    df[!,:area] = fill(NaN, size(df, 1))

    df.area[large_cell_mask] = round.(area.(convex_hull.(pos_data_per_cell[large_cell_mask])), sigdigits=sigdigits);
    df.density[large_cell_mask] = round.(df.n_transcripts[large_cell_mask] ./ df.area[large_cell_mask], sigdigits=sigdigits);
    df.elongation[large_cell_mask] = [round(x[2] / x[1], sigdigits=sigdigits) for x in eigvals.(cov.(transpose.(pos_data_per_cell[large_cell_mask])))];
    if :confidence in propertynames(segmented_df)
        df[!,:avg_confidence] = round.([mean(df.confidence) for df in seg_df_per_cell], sigdigits=sigdigits)
    end

    if dapi_arr !== nothing # TODO: forward DAPI from CLI
        brightness_per_mol = staining_value_per_transcript(segmented_df, dapi_arr);
        df[!, :mean_dapi_brightness] = mean.(trim.(split(brightness_per_mol, cell_assignment .+ 1)[2:end]; prop=0.2))
    end

    return df
end

function get_cell_stat_df(data::BmmData, segmented_df::Union{DataFrame, Nothing}=nothing; add_qc::Bool=true, sigdigits::Int=4)
    df = DataFrame(:cell => 1:length(data.components))

    centers = hcat([Vector(c.position_params.μ) for c in data.components]...)

    for s in [:x, :y]
        df[!,s] = mean.(split(data.x[!,s], data.assignment .+ 1, max_factor=length(data.components) + 1)[2:end])
    end

    if !isempty(data.cluster_per_cell)
        df[!, :cluster] = data.cluster_per_cell
    end

    if add_qc
        if segmented_df === nothing
            segmented_df = get_segmentation_df(data);
        end

        df = hcat(df, get_cell_qc_df(segmented_df; sigdigits=sigdigits, max_cell=length(data.components)))
    end

    return df[num_of_molecules_per_cell(data) .> 0,:]
end

function get_segmentation_df(data::BmmData, gene_names::Union{Nothing, Array{String, 1}}=nothing; use_assignment_history::Bool=true)
    df = deepcopy(data.x)
    df[!,:cell] = deepcopy(data.assignment);

    if use_assignment_history && (:assignment_history in keys(data.tracer)) && (length(data.tracer[:assignment_history]) > 1)
        df[!,:cell], df[!,:assignment_confidence] = estimate_assignment_by_history(data)
        df.assignment_confidence .= round.(df.assignment_confidence, digits=5)
    end

    if :confidence in propertynames(df)
        df.confidence = round.(df.confidence, digits=5)
    end

    df[!,:is_noise] = (df.cell .== 0);

    if gene_names !== nothing
        df[!,:gene] = gene_names[df[!,:gene]]
    end

    if !isempty(data.cluster_per_molecule)
        df[!, :cluster] = data.cluster_per_molecule
    end

    return df
end

function global_assignment_ids(data::BmmData)::Vector{Int}
    cur_guids = [c.guid for c in data.components]
    res = deepcopy(data.assignment)
    non_noise_mask = (res .> 0)
    res[non_noise_mask] .= cur_guids[res[non_noise_mask]]

    return res
end

function update_n_mols_per_segment!(bm_data::BmmData)
    if isempty(bm_data.segment_per_molecule)
        return
    end

    # Estimate number of molecules per segment per component
    for comp in bm_data.components
        empty!(comp.n_molecules_per_segment)
    end

    for i in 1:length(bm_data.assignment)
        c_cell = bm_data.assignment[i]
        c_seg = bm_data.segment_per_molecule[i]
        if (c_cell == 0) || (c_seg == 0)
            continue
        end

        bm_data.components[c_cell].n_molecules_per_segment[c_seg] = get(bm_data.components[c_cell].n_molecules_per_segment, c_seg, 0) + 1
    end

    # Estimate the main segment per cell
    resize!(bm_data.main_segment_per_cell, length(bm_data.components))
    bm_data.main_segment_per_cell .= 0

    for ci in 1:length(bm_data.components)
        if isempty(bm_data.components[ci].n_molecules_per_segment)
            continue
        end

        i_max = 0
        f_max = 0.0
        for (si,nms) in bm_data.components[ci].n_molecules_per_segment
            seg_size = bm_data.n_molecules_per_segment[si]
            if ((nms / seg_size) > (f_max + 1e-10)) || ((nms == seg_size) && (seg_size > bm_data.n_molecules_per_segment[i_max]))
                f_max = nms / seg_size
                i_max = si
            end
        end
        bm_data.main_segment_per_cell[ci] = i_max
    end
end

max_stars_repo_name: repushko/Baysor

max_issues_repo_name: repushko/Baysor

max_forks_repo_name: repushko/Baysor


text:
type TextWidget
    content::AbstractString
end

get_comm(w::TextWidget) = Comm(symbol("widget",object_id(w)))

function Base.writemime(io, ::MIME"text/html", w::TextWidget)
    id = object_id(w)
    content = w.content
    html = """

$content

    """
    write(io, html)
end

function update(w::TextWidget, content::AbstractString)
    comm = get_comm(w)
    w.content = content
    msg = Dict("content" => content)
    send_comm(comm, msg)
end


text:
abstract type AbstractKernel end

radial(k, θ::Tuple, τ) = radial(k, θ..., τ)

struct Gaussian <: AbstractKernel end

radial(::Gaussian, α::Number, τ) = α * exp(-τ / 2)

struct Spherical <: AbstractKernel end

function radial(::Spherical, α::Number, r::Number, τ)
    τ′ = τ / r
    return τ′ > 1 ? zero(τ) : α * (1 - (3 / 2) * τ′ + (1 / 2) * τ′ ^ 3)
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
"""
    centroid([geojson::T])::Point where {T <: AbstractGeometry}

Takes one or more features and calculates the centroid using the mean of all vertices.
This lessens the effect of small islands and artifacts when calculating the centroid of a set of polygons.

# Examples
```jldoctest
julia> line = LineString([[1, 2], [4, 6], [8, 9.5], [12, 13.4]])
LineString(Array{Float64,1}[[1.0, 2.0], [4.0, 6.0], [8.0, 9.5], [12.0, 13.4]])

julia> centroid(line)
Point([6.25, 7.725])
```
"""
function centroid(geojson::T)::Point where {T<:AbstractGeometry}
    x = 0.
    y = 0.
    len = 0.
    # TODO: support MultiPolygon

    data = geojson.coordinates
    type = geotype(geojson)

    if type === :Point
        x = data[1]
        y = data[2]
        len = 1.
    elseif type === :LineString || type === :MultiPoint
        for i in eachindex(data)
            x += data[i][1]
            y += data[i][2]
            len += 1
        end
    elseif type === :Polygon || type === :MultiLineString
        for i in eachindex(data[1])
            x += data[1][i][1]
            y += data[1][i][2]
            len += 1
        end
    end

    return Point([x / len, y / len])
end

"""
    centroid([geojson::T])::Point where {T <: AbstractFeatureCollection}

Takes one or more features and calculates the centroid using the mean of all vertices.
This lessens the effect of small islands and artifacts when calculating the centroid of a set of polygons.

# Examples
```jldoctest
julia> line = LineString([[1, 2], [4, 6], [8, 9.5], [12, 13.4]])
LineString(Array{Float64,1}[[1.0, 2.0], [4.0, 6.0], [8.0, 9.5], [12.0, 13.4]])

julia> centroid(line)
Point([6.25, 7.725])
```
"""
function centroid(geojson::T)::Point where {T <: AbstractFeatureCollection}
    x = 0.
    y = 0.
    len = 0.

    for feat in geojson.features
        geometry = feat.geometry

        data = geometry.coordinates
        type = geotype(geometry)

        if type === :Point
            x = data[1]
            y = data[2]
            len = 1.
        elseif type === :LineString || type === :MultiPoint
            for i in eachindex(data)
                x += data[i][1]
                y += data[i][2]
                len += 1.
            end
        elseif type === :Polygon || type === :MultiLineString
            for i in eachindex(data[1])
                x += data[1][i][1]
                y += data[1][i][2]
                len += 1.
            end
        end
    end

    return Point([x / len, y / len])
end

"""
    center(geojson::T) where {T <: AbstractGeometry}

Take a GeoJson Geometry and return the absolute center point.

# Examples
```jldoctest
julia> line = LineString([[1, 2], [4, 6], [8, 9.5], [12, 13.4]])
LineString(Array{Float64,1}[[1.0, 2.0], [4.0, 6.0], [8.0, 9.5], [12.0, 13.4]])

julia> center(line)
Point([6.5, 7.7])
```
"""
function center(geojson::T)::Point where {T <: Union{AbstractGeometry, AbstractFeatureCollection}}
    box = bbox(geojson)

    x = (box[1] + box[3]) / 2
    y = (box[2] + box[4]) / 2

    return Point([x, y])
end

"""
    masscenter(geojson::T) where {T <: AbstractGeometry}

Take any GeoJson Geometry and return its [center of mass](https://en.wikipedia.org/wiki/Center_of_mass) using this formula:
[Centroid of Polygon](https://en.wikipedia.org/wiki/Centroid#Centroid_of_polygon).

# Examples
```jldoctest
julia> line = LineString([[1, 2], [4, 6], [8, 9.5], [12, 13.4]])
LineString(Array{Float64,1}[[1.0, 2.0], [4.0, 6.0], [8.0, 9.5], [12.0, 13.4]])

julia> masscenter(line)
Point([6.25, 7.725])
```
"""
function masscenter(geojson::T) where {T <: AbstractGeometry}
    type = geotype(geojson)

    type === :Point && return Point(geojson.coordinates)
    if type === :Polygon
        coords = deepcopy(geojson.coordinates)

        center::Point = centroid(geojson)

        trans = center.coordinates
        sx = 0.
        sy = 0.
        sArea = 0.

        for i in eachindex(coords[1])
            coords[1][i][1] = coords[1][i][1] - trans[1]
            coords[1][i][2] = coords[1][i][2] - trans[2]
        end

        for i in 1:length(coords[1]) - 1
            pi = coords[1][i]
            xi, yi = pi

            pj = coords[1][i + 1]
            xj, yj = pj

            a = xi * yj - xj * yi
            sArea += a

            sx += (xi + xj) * a
            sy += (yi + yj) * a
        end

        sArea == 0 && return center

        area = sArea * 0.5
        factor = 1 / (6 * area)

        return Point([trans[1] + factor * sx, trans[2] + factor * sy])

    else
        # TODO: add convex hull
        return centroid(geojson)

    end
end

"""
    meancenter(geojson::T, weight::Real=1.) where {T <: Union{AbstractGeometry, AbstractFeatureCollection}}

Take a GeoJson Geometry and return the mean center. Can be weighted.

# Examples
```jldoctest
julia> line = LineString([[1, 2], [4, 6], [8, 9.5], [12, 13.4]])
LineString(Array{Float64,1}[[1.0, 2.0], [4.0, 6.0], [8.0, 9.5], [12.0, 13.4]])

julia> meancenter(line)
Point([6.25, 7.725])
```
"""
function meancenter(geojson::T, weight::Real=1.) where {T <: Union{AbstractGeometry, AbstractFeatureCollection}}
    # TODO: so ugly! Refactor?
    xs = 0.
    ys = 0.
    ns = 0.

    type = geotype(geojson)
    coors = []

    weight <= 0 && throw(error("Invalid weight. Weight must be > 0."))

    if type === :FeatureCollection

        for feat in geojson.features
            geom = feat.geometry
            coors = geom.coordinates

            if geotype(geom) === :Point
                xs += coors[1] * weight
                ys += coors[2] * weight
                ns += weight

            elseif geotype(geom) === :Polygon || geotype(geom) === :MultiLineString
                for i in eachindex(coors[1])
                    xs += coors[1][i][1] * weight
                    ys += coors[1][i][2] * weight
                    ns += weight
                end
            else
                for i in eachindex(coors)
                    xs += coors[i][1] * weight
                    ys += coors[i][2] * weight
                    ns += weight
                end

            end
        end

    else
        coors = geojson.coordinates

        if type === :Point
            xs += coors[1] * weight
            ys += coors[2] * weight
            ns += weight

        elseif type === :Polygon || type === :MultiLineString
            for i in eachindex(coors[1])
                xs += coors[1][i][1] * weight
                ys += coors[1][i][2] * weight
                ns += weight
            end
        else
            for i in eachindex(coors)
                xs += coors[i][1] * weight
                ys += coors[i][2] * weight
                ns += weight
            end
        end
    end

    return Point([xs / ns, ys / ns])
end

"""
    mediancenter(geojson::T, weight::Real=1., tol::Real=0.001, count::Integer=10) where {T <: AbstractFeatureCollection}

Take a FeatureCollection of points and calculate the median center,
algorithimically. The median center is understood as the point that is
requires the least total travel from all other points.

# Examples
```jldoctest
julia> fc = FeatureCollection([Feature(Point([0, 0])),Feature(Point([9, 9])), Feature(Point([9.25, 9.25])),
             Feature(Point([9.5, 9.5])), Feature(Point([9.75, 9.75])), Feature(Point([10, 10]))])
FeatureCollection{Feature}(Feature[Feature(Point([0.0, 0.0]), Dict{String,Any}()), Feature(Point([9.0, 9.0]), Dict{String,Any}()), Feature(Point([9.25, 9.25]), Dict{String,Any}()), Feature(Point([9.5, 9.5]), Dict{String,Any}()), Feature(Point([9.75, 9.75]), Dict{String,Any}()), Feature(Point([10.0, 10.0]), Dict{String,Any}())], nothing, nothing)

julia> mediancenter(fc)
Point([9.25425, 9.25425])
```
"""
function mediancenter(geojson::T, weight::Real=1., tol::Real=0.001, count::Integer=10) where {T <: AbstractFeatureCollection}

    mean = meancenter(geojson)

    feat1 = Feature(nothing)

    feat1.properties = Dict(
        "tolerance" => tol,
        "weight" => weight,
        "candidates" => [])

    centroids = FeatureCollection([feat1])

    for feat in geojson.features
        push!(centroids.features, Feature(centroid(feat.geometry)))
    end

    return findMedian(mean.coordinates, [0., 0.], centroids, count)

end

function findMedian(candidate::Position, previous::Position, centroids::FeatureCollection, counter::Integer)

    i = findfirst(x -> x.properties["candidates"] != nothing, centroids.features)
    feat = centroids.features[i]

    tol = feat.properties["tolerance"]
    weight = feat.properties["weight"]
    xSum = 0.
    ySum = 0.
    kSum = 0.
    count = 0

    weight <= 0 && throw(error("Invalid weight. Weight must be > 0."))

    for feats in centroids.features

        feats.geometry == nothing && continue

        count += 1
        dist = weight * distance(feats.geometry, Point(candidate))

        dist === 0 && (dist = 1.)

        k = weight / dist

        xSum += feats.geometry.coordinates[1] * k
        ySum += feats.geometry.coordinates[2] * k

        kSum += k
    end

    count < 1 && throw(error("No features to measure."))

    candidateX = xSum / kSum
    candidateY = ySum / kSum

    (count === 1 || counter === 0 || (abs(candidateX - previous[1])) < tol && (abs(candidateY - previous[2])) < tol) &&
        return Point([candidateX, candidateY])

    push!(feat.properties["candidates"], [candidateX, candidateY])

    return findMedian([candidateX, candidateY], candidate, centroids, counter - 1)
end


text:
module CodecBzip2

export
    Bzip2Compressor,
    Bzip2CompressorStream,
    Bzip2Decompressor,
    Bzip2DecompressorStream

import TranscodingStreams:
    TranscodingStreams,
    TranscodingStream,
    Memory,
    Error,
    initialize,
    finalize,
    splitkwargs
using Libdl

const libbz2path = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if !isfile(libbz2path)
    error("CodecBzip2.jl is not installed properly, run Pkg.build(\"CodecBzip2\") and restart Julia.")
end
include(libbz2path)
check_deps()

include("libbz2.jl")
include("compression.jl")
include("decompression.jl")

end # module


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later
include("compile.jl")

# env var to force compilation from source, for testing purposes
const forcecompile = get(ENV, "FORCE_COMPILE_NLOPT", "no") == "yes"

# Parse some basic command-line arguments
const verbose = ("--verbose" in ARGS) || forcecompile
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, String["libnlopt_cxx"], :libnlopt),
]
verbose && forcecompile && @info("Forcing compilation from source.")

# Download binaries from hosted location
bin_prefix = "https://github.com/stevengj/NLoptBuilder/releases/download/v2.5.0"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, :glibc) => ("$bin_prefix/NLopt.v2.5.0.aarch64-linux-gnu.tar.gz", "2347e5e240b9b95cfe6717d356e1f1948758e983574bc97378bee0fc90e21936"),
    Linux(:aarch64, :musl) => ("$bin_prefix/NLopt.v2.5.0.aarch64-linux-musl.tar.gz", "0bf5cf5676a8df605cae199cd2a07b541f54a2c6c78f98109a863811d77cc268"),
    Linux(:armv7l, :glibc, :eabihf) => ("$bin_prefix/NLopt.v2.5.0.arm-linux-gnueabihf.tar.gz", "2436981433c7bd7bab8853d26341a01c82a553528232ae950bf390b399678f6f"),
    Linux(:armv7l, :musl, :eabihf) => ("$bin_prefix/NLopt.v2.5.0.arm-linux-musleabihf.tar.gz", "65cba09664827d6f7baf2324594b1e652912ef38dc0cb15adab2f99cfde4f1ac"),
    Linux(:i686, :glibc) => ("$bin_prefix/NLopt.v2.5.0.i686-linux-gnu.tar.gz", "61602ad8e8924badc2836e5be8279b9f85afb8608bfbb6c1bebe13eef755c2f4"),
    Linux(:i686, :musl) => ("$bin_prefix/NLopt.v2.5.0.i686-linux-musl.tar.gz", "cebfc02353c3cc16778789d4654dff193d186cd3c453b42712cc888ca387a0f6"),
    Windows(:i686) => ("$bin_prefix/NLopt.v2.5.0.i686-w64-mingw32.tar.gz", "b1d3f963866f383195bf9a0cf1c280e4817f8d5238911c830b8a6cdda949e817"),
    Linux(:powerpc64le, :glibc) => ("$bin_prefix/NLopt.v2.5.0.powerpc64le-linux-gnu.tar.gz", "a3515af1d464d420a0e78e61dabf6735a1320acfc898f8a5284093c7bb5fa57d"),
    MacOS(:x86_64) => ("$bin_prefix/NLopt.v2.5.0.x86_64-apple-darwin14.tar.gz", "623b8be18dba99021befbebaa50f73f70b9949849e8e916253c657038f014f14"),
    Linux(:x86_64, :glibc) => ("$bin_prefix/NLopt.v2.5.0.x86_64-linux-gnu.tar.gz", "ea0ee00881fa1517f2640a1e05ff9ea88757c1c148fde1c42dd124791684e0b4"),
    Linux(:x86_64, :musl) => ("$bin_prefix/NLopt.v2.5.0.x86_64-linux-musl.tar.gz", "d20ab1cfb9b5561ffaa70da919df4e409ca7ed7fef8e4b049d3ef6e24344a8a8"),
    FreeBSD(:x86_64) => ("$bin_prefix/NLopt.v2.5.0.x86_64-unknown-freebsd11.1.tar.gz", "078e9a08283d471473da5c6f465d8f907f068da8ef15ad36305a231cf776a309"),
    Windows(:x86_64) => ("$bin_prefix/NLopt.v2.5.0.x86_64-w64-mingw32.tar.gz", "776ea5d6777242f86355de2fdd6d5be8203e60b38ea5794c53f45c36551fde46"),
)

# source code tarball and hash for fallback compilation
source_url = "https://github.com/stevengj/nlopt/archive/v2.5.0.tar.gz"
source_hash = "c6dd7a5701fff8ad5ebb45a3dc8e757e61d52658de3918e38bab233e7fd3b4ae"

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
if haskey(download_info, platform_key()) && !forcecompile
    url, tarball_hash = download_info[platform_key()]
    if !isinstalled(url, tarball_hash; prefix=prefix)
        # Download and install binaries
        install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose)

        # check again whether the dependency is satisfied, which
        # may not be true if dlopen fails due to a libc++ incompatibility (#50)
        unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
    end
end

if unsatisfied || forcecompile
    # Fall back to building from source, giving the library a different name
    # so that it is not overwritten by BinaryBuilder downloads or vice-versa.
    libname = "libnlopt_from_source"
    products = [ LibraryProduct(prefix, [libname], :libnlopt) ]
    source_path = joinpath(prefix, "downloads", "src.tar.gz")
    if !isfile(source_path) || !verify(source_path, source_hash; verbose=verbose) || !satisfied(products[1]; verbose=verbose)
        compile(libname, source_url, source_hash, prefix=prefix, verbose=verbose)
    end
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)


text:
module MyUtils

export flatten, not_nan, myColors, cartesian2matrix, cartesian2polar, polar2cartesian, append_sessions, convexHullCircum, foo

using Plots
using Polyhedra
import GLPK
lib = DefaultLibrary{Float64}(GLPK.Optimizer)

function flatten(x)
    return collect(Iterators.flatten(x))
end

function myColors()
    return distinguishable_colors(25, [parse(Colorant,c) for c in  ["green", "yellow"]])[[1, 3, 5, 8, 12, 13, 14, 15, 16, 17, 18, 19, 21, 22, 23, 24]]
end

function cartesian2matrix(ci::AbstractArray)
    ### array of cartesian indices to matrix N x DIM
    dims = collect(1:length(ci[1]))
    sparse_idx = reduce(vcat, [[t[d] for d in dims]' for t in ci])
end

function cartesian2polar(x, y)
    ### cartesian coordinates to polar coordinates
    if size(x) == size(y)
        rphi = zeros(length(x), 2)
        rphi[:,1] = sqrt.(x.^2 .+ y.^2)
        rphi[:,2] = atan.(x ./ y)
        return rphi
    else
        @warn "Dimensions not matching"
    end
end

function polar2cartesian(r, phi)
    ###
    xy = zeros(length(r), 2)
    xy[:,1] = r.*cos.(phi)
    xy[:,2] = r.*sin.(phi)
    return xy
end

function append_sessions( rd_a, rd; which_sessions)

    #
    # INPUT
    # rd_a: empty DataFrame() !!!(required since recursive implementation)!!!
    # rd: Array{Newbe.Experiment.data_of_run,1}
    # which_sessions: array with indicies of sessions in rd to concatenate
    #
    # OUTPUT
    # rd_a: concatenated DataFrame()
    #

    if !isempty(which_sessions)
        rd_a = [rd_a; getfield( rd[which_sessions[1]], :trial_list)]
        which_sessions = which_sessions[2:end]
        append_sessions( rd_a, rd; which_sessions=which_sessions)
    else
        return rd_a
    end

end

function convexHullCircum(xy::Array{Array{Float64,1},1})

    xy_ = copy(xy)
    [push!(xy_[ii], xy_[ii][1]) for ii in [1,2]]

    P1 = polyhedron(vrep(
            hcat(xy_[1], xy_[2])
            ), lib)
    removevredundancy!(P1)

    ext = [P1.vrep.V; [P1.vrep.V[1,1] P1.vrep.V[1,2]]]

    diffs = diff(ext, dims=1).^2 # compute the differences along x and y dimensions
    return [P1, sum(sqrt.(sum(diffs, dims=2)))] #return of the length of the hull circumference

end

function non_nan(x_::Array{T,N}) where {T<:Number,N}
    return filter( el -> !isnan(el), x_ )
end

end

max_stars_repo_name: maalaria/Newbe

max_issues_repo_name: maalaria/Newbe

max_forks_repo_name: maalaria/Newbe


text:
# make your scripts automatically re-activate your project
cd(@__DIR__)
using Pkg; Pkg.activate("."); Pkg.instantiate()

# load packages
using Flux, DiffEqFlux, DiffEqSensitivity
using DifferentialEquations
using Plots
using BSON
using Statistics
using ReverseDiff, Zygote, ForwardDiff
using LinearAlgebra
using Random
using BenchmarkTools, Test

#################################################
# read parameters from command line
lr = 0.0015f0 #parse(Float64,ARGS[1]) #0.05
epochs = 400 #parse(Int,ARGS[2]) #100
slurmidx = 1 #parse(Int,ARGS[3]) #1

numtraj = 256 # number of trajectories in parallel simulations for training
numtrajplot = 256 # .. for plotting

# time range for the solver
dt = 0.01f0 #0.001f0
tinterval = 0.02f0
tstart = 0.0f0
Nintervals = 150 # total number of intervals, total time = t_interval*Nintervals
tspan = (tstart,tinterval*Nintervals)

# Hamiltonian parameters
Δ = 20.0f0
Ωmax = 10.0f0 # control parameter (maximum amplitude)

# loss hyperparameters
C1 = Float32(1.0)  # evolution state fidelity
C2 = Float32(0.0) # action amplitudes
C3 = Float32(0.0) # evolution state fidelity for last few steps!

struct Parameters{flType,intType,tType}
	lr::flType
	epochs::intType
	numtraj::intType
	numtrajplot::intType
	dt::flType
	tinterval::flType
	tspan::tType
	Nintervals::intType
	Δ::flType
	Ωmax::flType
	C1::flType
	C2::flType
	C3::flType
end

myparameters = Parameters{typeof(dt),typeof(numtraj), typeof(tspan)}(
  lr, epochs, numtraj, numtrajplot, dt, tinterval, tspan, Nintervals,
  Δ, Ωmax, C1, C2, C3)

################################################
# Define Neural Network

# state-aware
nn = FastChain(
	FastDense(4, 256, relu, initW = Flux.glorot_uniform, initb = Flux.glorot_uniform),
	FastDense(256, 64, relu, initW = Flux.glorot_uniform, initb = Flux.glorot_uniform),
	FastDense(64, 1, softsign, initW = Flux.glorot_uniform, initb = Flux.glorot_uniform))
p_nn = initial_params(nn)

###############################################
# initial state anywhere on the Bloch sphere
function prepare_initial(dt, n_par)
  # shape 4 x n_par
  # input number of parallel realizations and dt for type inference
  # random position on the Bloch sphere
  theta = acos.(2*rand(typeof(dt),n_par).-1)  # uniform sampling for cos(theta) between -1 and 1
  phi = rand(typeof(dt),n_par)*2*pi  # uniform sampling for phi between 0 and 2pi
  # real and imaginary parts ceR, cdR, ceI, cdI
  u0 = [cos.(theta/2), sin.(theta/2).*cos.(phi), false*theta, sin.(theta/2).*sin.(phi)]
  return vcat(transpose.(u0)...) # build matrix
end

# target state
# ψtar = |up>

u0 = prepare_initial(myparameters.dt, myparameters.numtraj)

###############################################
# Define ODE

function qubit(u,p,t)
  # expansion coefficients |Ψ> = ce |e> + cd |d>
  ceR, cdR, ceI, cdI = u # real and imaginary parts

  # Δ: atomic frequency
  # Ω: Rabi frequency for field in x direction
  Δ = p[end-1]
  Ωmax = p[end]
  nn_weights = p[1:end-2]

  Ω = (nn(u, nn_weights).*Ωmax)[1]

  dceR =  1//2*(ceI*Δ+cdI*Ω)
  dcdR = -cdI*Δ/2 +ceI*Ω/2
  dceI = 1//2*(-ceR*Δ-cdR*Ω)
  dcdI = cdR*Δ/2-ceR*Ω/2

  return [dceR, dcdR, dceI, dcdI]
end

# normalization callback

condition(u,t,integrator) = true
function affect!(integrator)
  integrator.u=integrator.u/norm(integrator.u)
end
cb = DiscreteCallback(condition,affect!,save_positions=(false,false))

# get control pulses
p_all = [p_nn; myparameters.Δ; myparameters.Ωmax]
# define ODE problem
prob = ODEProblem{false}(qubit, vec(u0[:,1]), myparameters.tspan, p_all,
   callback=cb
   )

#########################################
# compute loss
function g(u,p,t)
  ceR = @view u[1,:,:]
  cdR = @view u[2,:,:]
  ceI = @view u[3,:,:]
  cdI = @view u[4,:,:]
  p[1]*mean((cdR.^2 + cdI.^2) ./ (ceR.^2 + cdR.^2 + ceI.^2 + cdI.^2))
end

function loss(p, u0, myparameters::Parameters; sensealg = ForwardDiffSensitivity())
  pars = [p; myparameters.Δ; myparameters.Ωmax]

  function prob_func(prob, i, repeat)
    # prepare initial state and applied control pulse
	u0tmp = deepcopy(vec(u0[:,i]))
    remake(prob,
		p = pars,
		u0 = u0tmp,
		callback = cb
		)
  end

  ensembleprob = EnsembleProblem(prob,
   prob_func = prob_func,
   safetycopy = true
   )

  _sol = solve(ensembleprob, Tsit5(), EnsembleThreads(),
    sensealg=sensealg,
    saveat=myparameters.tinterval,
    dt=myparameters.dt,
    adaptive=true, abstol=1e-6, reltol=1e-6,
    trajectories=myparameters.numtraj, batch_size=myparameters.numtraj)

  A = convert(Array,_sol)
  loss = g(A,[myparameters.C1,myparameters.C2,myparameters.C3],nothing)

  return loss
end

#########################################
# visualization -- run for new batch
function visualize(p, u0, myparameters::Parameters; all_traj = true)
  # # initialization
  u = deepcopy(u0)

  pars = [p; myparameters.Δ; myparameters.Ωmax]

  function prob_func(prob, i, repeat)
    # prepare initial state and applied control pulse
	remake(prob,
		p = pars,
		u0 = vec(u0[:,i]),
		callback = cb
		)
  end

  ensembleprob = EnsembleProblem(prob,
   prob_func = prob_func,
   safetycopy = true
   )

  u = solve(ensembleprob, Tsit5(), ensemblealg=EnsembleThreads(),
	saveat=myparameters.tinterval,
	dt=myparameters.dt,
	adaptive=true, abstol=1e-6, reltol=1e-6,
	trajectories=myparameters.numtrajplot, batch_size=myparameters.numtrajplot)

  ceR = @view u[1,:,:]
  cdR = @view u[2,:,:]
  ceI = @view u[3,:,:]
  cdI = @view u[4,:,:]
  infidelity = @. cdR^2 + cdI^2 / (ceR^2 + cdR^2 + ceI^2 + cdI^2)
  meaninfidelity = mean(infidelity)
  loss = myparameters.C1*meaninfidelity

  @info "Loss: " loss

  fidelity = @. ceR^2 + ceI^2 / (ceR^2 + cdR^2 + ceI^2 + cdI^2)

  mf = mean(fidelity, dims=2)[:]
  sf = std(fidelity, dims=2)[:]

  # re-compute actions
  arrayu = Array(u)
  Ωlist = []
  for i = 1:(size(arrayu)[2])
    Ω = vec(nn(arrayu[:,i,:],p).*myparameters.Ωmax)
    push!(Ωlist, Ω)
  end
  Ωlist = hcat(Ωlist...)

  ma = mean(Ωlist, dims=1)[:]
  sa = std(Ωlist, dims=1)[:]

  pl1 = plot(0:myparameters.Nintervals, mf,
		ribbon = sf,
		ylim = (0,1), xlim = (0,myparameters.Nintervals),
		c=1, lw = 1.5, xlabel = "steps", ylabel="Fidelity", legend=false)
  pl2 = plot(0:myparameters.Nintervals, ma,
		ribbon = sa,
		ylim=(-myparameters.Ωmax,myparameters.Ωmax), xlim = (0,myparameters.Nintervals),
		c=2, lw = 1.5, xlabel = "steps", ylabel="Ω(t)", legend=false)
  if all_traj
     plot!(pl1, fidelity, legend=false, c=:gray, alpha=0.1)
     plot!(pl2, Ωlist', legend=false, c=:gray, alpha=0.1)
  else
     plot!(pl1, 0:myparameters.Nintervals, fidelity[:,end],  c=:gray, lw = 1.5, legend=false)
     plot!(pl2, 0:myparameters.Nintervals, Ωlist[end,:], c=:gray, lw = 1.5, legend=false)
  end

  pl = plot(pl1, pl2, layout = (1, 2), legend = false, size=(800,360))
  return pl, loss
end

###################################
# training loop

# optimize the parameters for a few epochs with ADAM on time span Nint
opt = ADAM(myparameters.lr)
list_plots = []
losses = []
for epoch in 1:myparameters.epochs
  println("epoch: $epoch / $(myparameters.epochs)")
  local u0 = prepare_initial(myparameters.dt, myparameters.numtraj)
  _dy, back = @time Zygote.pullback(p -> loss(p, u0, myparameters,
  	sensealg=InterpolatingAdjoint()), p_nn)
  gs = @time back(one(_dy))[1]
  push!(losses, _dy)
  if epoch % myparameters.epochs == 0
    # plot every xth epoch
    local u0 = prepare_initial(myparameters.dt, myparameters.numtrajplot)
    pl, test_loss = visualize(p_nn, u0, myparameters)
    println("Loss (epoch: $epoch): $test_loss")
    display(pl)
    push!(list_plots, pl)
  end
  Flux.Optimise.update!(opt, p_nn, gs)
  println("")
end

# plot training loss
pl = plot(losses, lw = 1.5, xlabel = "some epochs", ylabel="Loss", legend=false)

###################################
# Serialization

bson("Data/ODEControlCont-epochs="*string(myparameters.epochs)*"_numtraj="*join(string.([myparameters.numtraj,myparameters.numtrajplot]), '_')*"_"*string(slurmidx)*".bson",
  Dict(:losses => losses, :popt => p_nn, :lr =>myparameters.lr, :epoch =>myparameters.epochs))

savefig(pl,"./Figures/ODELossCont_epochs="*string(myparameters.epochs)*"_numtraj="*join(string.([myparameters.numtraj,myparameters.numtrajplot]), '_')*"_"*string(slurmidx)*".png")

for (i,plt) in enumerate(list_plots)
  savefig(plt,"./Figures/ODEControlCont"*string(i)*"_epochs="*string(myparameters.epochs)*"_numtraj="*join(string.([myparameters.numtraj,myparameters.numtrajplot]), '_')*"_"*string(slurmidx)*".png")
end

exit()


text:
@testset "event" begin

    event = T.KeyPressedEvent('c', T.CTRL)
    @test event.key === 'c'
    @test event.adjoint_key === T.CTRL
    @test event === T.KeyPressedEvent('c', T.CTRL)
    @test repr(event) === "KeyPressedEvent(CTRL+'c')"

    event = T.KeyPressedEvent(T.ESC)
    @test event.key === T.ESC
    @test event.adjoint_key === T.NO_ADJOINT
    @test repr(event) === "KeyPressedEvent('ESC')"

    event = T.PasteEvent("Paste Event")
    @test event.content === "Paste Event"
    @test event === T.PasteEvent("Paste Event")
    @test repr(event) === "PasteEvent(\"Paste Event\")"

end


text:
module EBayes

using Reexport

@reexport using StatsBase
import StatsBase:predict,
                 fit,
                 dof,
                 leverage,
				 response

@reexport using Distributions
import Distributions:location, RealInterval, support

import Base: eltype,
             firstindex,
             getindex,
             getproperty,
             iterate,
             lastindex,
             length,
             rand,
             size,
             zero,
             zeros
import Base.Broadcast: broadcastable

import MLJBase # to avoid conflicts with StatsBase
import MLJBase:predict_mean
import MLDataPattern:kfolds,FoldsView,shuffleobs

using LaTeXStrings
using NLSolversBase
using Optim
using RecipesBase
using Statistics
using StructArrays

import Statistics:var

include("ebayes_types.jl")
include("ebayes_samples.jl")
include("predict.jl")
include("shrinkage_locations.jl")
include("sure.jl")
include("ebcf.jl")
include("example_priors.jl")
include("simulations.jl")

export EBayesSample,
       NormalSample,
       StandardNormalSample,
       EBayesSamples,
       AbstractNormalSamples,
       NormalSamples,
       NormalEBayesSimulationResult,
       var, #reexport from Statistics
       predict_mean, #reexport from MLJBase
       FixedLocation,
       GrandMeanLocation,
       SURE,
       EBayesCrossFit

end # module


text:
using LibGit2

struct LoginBuffer <: IO
    secret::Base.SecretBuffer
end

function Base.readavailable(data::LoginBuffer)
    size = data.secret.size
    data.secret.ptr = size + 1
    return @view data.secret.data[1:size]
end
Base.eof(data::LoginBuffer) = Base.eof(data.secret)

function withbuf(f::Function, data::LoginBuffer)
    try
        return f(data)
    finally
        Base.shred!(data.secret)
    end
end

function LoginBuffer(cred::LibGit2.UserPasswordCredential)
    input = Base.SecretBuffer()
    write(input, "{\"username\" : \"")
    write(input, cred.user)
    write(input, "\", \"password\" : \"")
    write(input, cred.pass)
    write(input, "\"}")
    seekstart(input)
    return LoginBuffer(input)
end

function login_prompt(user::AbstractString = Base.prompt("username"))
    return LibGit2.UserPasswordCredential(user, Base.getpass("password"))
end

function erase_cred(user::AbstractString)
    helper = LibGit2.GitCredentialHelper(`git credential-store`)
    userpass = LibGit2.UserPasswordCredential(user)
    cred = LibGit2.GitCredential(userpass, get_endpoint())
    LibGit2.reject(helper, cred)
    Base.shred!(cred)
    Base.shred!(userpass)
    return nothing
end

function write_cred(userpass::LibGit2.UserPasswordCredential)
    helper = LibGit2.GitCredentialHelper(`git credential-store`)
    cred = LibGit2.GitCredential(userpass, get_endpoint())
    LibGit2.approve(helper, cred)
    Base.shred!(cred)
    return nothing
end

function read_cred(user::AbstractString)
    helper = LibGit2.GitCredentialHelper(`git credential-store`)
    userpass = LibGit2.UserPasswordCredential(user)
    cred = LibGit2.GitCredential(userpass, get_endpoint())
    LibGit2.fill!(helper, cred)
    Base.shred!(userpass)
    return cred
end

function maybe_login_prompt(user::AbstractString = Base.prompt("username"))
    cred = read_cred(user)
    if isempty(cred.password)
        Base.shred!(cred)
        return login_prompt(user)
    else
        return LibGit2.UserPasswordCredential(cred.username, cred.password)
    end
end

auth_pair(auth_token) = "authorization"=>"Bearer $auth_token"
auth_header(auth_token) = [auth_pair(auth_token)]
auth_header!(headers, auth_token) = push!(headers, auth_pair(auth_token))
auth_pair(::Nothing) = nothing
auth_header(::Nothing) = Pair{String,String}[]
auth_header!(headers, ::Nothing) = headers

ensure_token(::Nothing) = error("No token provided and saved token not found")
ensure_token(token) = token

function get_token()
    path = get_token_path()
    return isfile(path) ? read(path, String) : nothing
end

function save_token(token)
    path = get_token_path()
    mkpath(dirname(path))
    write(path, token)
    return nothing
end

delete_token() = rm(get_token_path(); force=true)


text:
include("../utils.jl")
include("group_utils.jl")

using OrdinaryDiffEq
import Manifolds: invariant_metric_dispatch, default_metric_dispatch, local_metric

struct TestInvariantMetricBase <: Metric end

function local_metric(
    ::MetricManifold{𝔽,<:Manifold,TestInvariantMetricBase},
    ::Identity,
) where {𝔽}
    return Diagonal([1.0, 2.0, 3.0])
end
function local_metric(
    ::MetricManifold{𝔽,<:Manifold,<:InvariantMetric{TestInvariantMetricBase}},
    p,
) where {𝔽}
    return Diagonal([1.0, 2.0, 3.0])
end

struct TestBiInvariantMetricBase <: Metric end

function invariant_metric_dispatch(
    ::MetricManifold{𝔽,<:Manifold,<:InvariantMetric{TestBiInvariantMetricBase}},
    ::ActionDirection,
) where {𝔽}
    return Val(true)
end

function local_metric(
    ::MetricManifold{𝔽,<:Manifold,<:TestBiInvariantMetricBase},
    ::Identity,
) where {𝔽}
    return Diagonal(0.4I, 3)
end

struct TestInvariantMetricManifold <: Manifold{ℝ} end

struct TestDefaultInvariantMetricManifold <: Manifold{ℝ} end

function default_metric_dispatch(
    ::MetricManifold{
        ℝ,
        TestDefaultInvariantMetricManifold,
        RightInvariantMetric{TestInvariantMetricBase},
    },
)
    return Val(true)
end

invariant_metric_dispatch(::TestDefaultInvariantMetricManifold, ::RightAction) = Val(true)

@testset "Invariant metrics" begin
    base_metric = TestInvariantMetricBase()
    metric = InvariantMetric(base_metric)
    lmetric = LeftInvariantMetric(base_metric)
    rmetric = RightInvariantMetric(base_metric)

    @test InvariantMetric(base_metric) === InvariantMetric(base_metric, LeftAction())
    @test lmetric === InvariantMetric(base_metric, LeftAction())
    @test rmetric === InvariantMetric(base_metric, RightAction())
    @test sprint(show, lmetric) == "LeftInvariantMetric(TestInvariantMetricBase())"
    @test sprint(show, rmetric) == "RightInvariantMetric(TestInvariantMetricBase())"

    @test direction(lmetric) === LeftAction()
    @test direction(rmetric) === RightAction()

    G = MetricManifold(TestInvariantMetricManifold(), lmetric)
    @test (@inferred invariant_metric_dispatch(G, LeftAction())) === Val(true)
    @test (@inferred invariant_metric_dispatch(G, RightAction())) === Val(false)

    G = MetricManifold(TestInvariantMetricManifold(), rmetric)
    @test (@inferred invariant_metric_dispatch(G, RightAction())) === Val(true)
    @test (@inferred invariant_metric_dispatch(G, LeftAction())) === Val(false)

    @test Manifolds.invariant_metric_dispatch(
        TestInvariantMetricManifold(),
        RightAction(),
    ) === Val{false}()
    @test Manifolds.invariant_metric_dispatch(
        TestInvariantMetricManifold(),
        LeftAction(),
    ) === Val{false}()

    G = MetricManifold(
        TestDefaultInvariantMetricManifold(),
        LeftInvariantMetric(TestInvariantMetricBase()),
    )
    @test !is_default_metric(G)
    G = MetricManifold(
        TestDefaultInvariantMetricManifold(),
        RightInvariantMetric(TestInvariantMetricBase()),
    )
    @test is_default_metric(G)

    e = Matrix{Float64}(I, 3, 3)
    @testset "inner/norm" begin
        SO3 = SpecialOrthogonal(3)
        p = exp(hat(SO3, Identity(SO3, e), [1.0, 2.0, 3.0]))
        X = hat(SO3, Identity(SO3, e), [2.0, 3.0, 4.0])
        Y = hat(SO3, Identity(SO3, e), [3.0, 4.0, 1.0])

        G = MetricManifold(SO3, lmetric)
        @test inner(G, p, X, Y) ≈ dot(X, Diagonal([1.0, 2.0, 3.0]) * Y)
        @test norm(G, p, X) ≈ sqrt(inner(G, p, X, X))

        G = MetricManifold(SO3, rmetric)
        @test inner(G, p, X, Y) ≈ dot(p * X * p', Diagonal([1.0, 2.0, 3.0]) * p * Y * p')
        @test norm(G, p, X) ≈ sqrt(inner(G, p, X, X))
    end

    @testset "log/exp bi-invariant" begin
        SO3 = SpecialOrthogonal(3)
        p = exp(hat(SO3, Identity(SO3, e), [1.0, 2.0, 3.0]))
        q = exp(hat(SO3, Identity(SO3, e), [3.0, 4.0, 1.0]))
        X = hat(SO3, Identity(SO3, e), [2.0, 3.0, 4.0])

        G = MetricManifold(SO3, InvariantMetric(TestBiInvariantMetricBase(), LeftAction()))
        @test isapprox(SO3, exp(G, p, X), exp(SO3, p, X))
        @test isapprox(SO3, p, log(G, p, q), log(SO3, p, q); atol = 1e-6)

        G = MetricManifold(SO3, InvariantMetric(TestBiInvariantMetricBase(), RightAction()))
        @test isapprox(SO3, exp(G, p, X), exp(SO3, p, X))
        @test isapprox(SO3, p, log(G, p, q), log(SO3, p, q); atol = 1e-6)
    end

    @testset "exp τ-invariant" begin
        T3 = TranslationGroup(3)
        p = [1.0, 2.0, 3.0]
        X = [3.0, 5.0, 6.0]
        @test isapprox(T3, exp(MetricManifold(T3, lmetric), p, X), p .+ X)
        @test isapprox(T3, exp(MetricManifold(T3, rmetric), p, X), p .+ X)
    end
end


text:
export plot_atom
export plot_atom!
export plot_graphene

@recipe function f(g::Vector{T}) where T<:AbstractGPrimitive
         get_x.(g), get_y.(g)
         end

function plot_atom!(plt, g::Vector{T}; kw...) where T<:AbstractGPrimitive
    atoms = filter(isatom, g)
    bonds = filter(isbond, g)
    polygons = filter(ispolygon, g)
    scatter!(plt, get_x.(atoms), get_y.(atoms); kw...)
end
plot_atom!(plt, g::T; kw...) where T<:AbstractGPrimitive = plot_atom!(plt, [g], kw...)

function plot_atom!(g::Vector{T}; kw...) where T<:AbstractGPrimitive
    atoms = filter(isatom, g)
    bonds = filter(isbond, g)
    polygons = filter(ispolygon, g)
    scatter!(get_x.(atoms), get_y.(atoms); kw...)
end
plot_atom!(g::T; kw...) where T<:AbstractGPrimitive = plot_atom!([g], kw...)

function plot_atom(g::Vector{T}; kw...) where T<:AbstractGPrimitive
    atoms = filter(isatom, g)
    bonds = filter(isbond, g)
    polygons = filter(ispolygon, g)
    scatter(get_x.(atoms), get_y.(atoms); kw...)
end
plot_atom(g::T; kw...) where T<:AbstractGPrimitive = plot_atom([g], kw...)

function plot_graphene(graphene, g::Vector{T}; kw...) where T<:AbstractGPrimitive
    plt = plot([], aspect_ratio = 1, xlims=(0,256), ylims=(0,256), framestyle=:box, legend=false)

    g_palette = [
    colorant"#6E6E6E", # 1
    colorant"#DF0101", # 2
    colorant"#FBD606", # 3
    colorant"#BF2626", # 4
    colorant"#FEE081", # 5
    colorant"#B3BED2", # 6
    colorant"#353D52", # 7
    colorant"#0D621E", # 8
    colorant"#934BB2", # 9
    colorant"#3EAE7E", # 10
    colorant"#EFBAA0", # 11
    colorant"#602E56", # 12
    colorant"#3f9778", # 13
    colorant"#FF8000", # 14
    colorant"#5FB404", # 15
    colorant"#182854", # 16
    colorant"#ACAEB5", # 17
    colorant"#5A7460", # 18
    colorant"#B93F01", # 19
    colorant"#1897AB", # 20
    colorant"#FFFFFF"] # >20

    gatoms = filter(isatom, g)
    # gbonds = filter(isbond, g)
    gpolygons = sort(filter(ispolygon, g), rev=true)
    for n in gpolygons
        points = [[get_x(p), get_y(p)] for p in filter_relatives_by_type(graphene, n, "Atom")]
        hull = convex_hull(points)
        noa = get_noa(n)
        plot!(plt, VPolygon(hull), alpha=1, color=g_palette[min(noa,21)])
    end
    if !isempty(gatoms)
        plot_atom!(plt, gatoms, color=g_palette[1])
    end
    return plt
end
plot_graphene(graphene, g::T; kw...) where T<:AbstractGEntry = plot_polygon(graphene, [g]; kw...)

function plot_graphene(graphene, g_vector::Vector{GDefect}; kw...)
    temp = AbstractGPrimitive[]
    for g in g_vector
        temp = [graphene[collect(get_members(g))]; temp]
    end
    plot_graphene(graphene, temp; kw...)
end


text:
function get_ylabel(name::String)
    if name == "Phosphorylated_MEKc"
        return "Phosphorylated MEK\n(cytoplasm)"
    elseif name == "Phosphorylated_ERKc"
        return "Phosphorylated ERK\n(cytoplasm)"
    elseif name == "Phosphorylated_RSKw"
        return "Phosphorylated RSK\n(whole cell)"
    elseif name == "Phosphorylated_CREBw"
        return "Phosphorylated CREB\n(whole cell)"
    elseif name == "dusp_mRNA"
        return L"$\it{dusp}$"*" mRNA\nexpression"
    elseif name == "cfos_mRNA"
        return L"$\it{c}$"*"-"*L"$\it{fos}$"*" mRNA\nexpression"
    elseif name == "cFos_Protein"
        return "c-Fos Protein\nexpression"
    elseif name == "Phosphorylated_cFos"
        return "Phosphorylated c-Fos\nProtein expression"
    end
end

function plotFunc_timecourse(Sim::Module, n_file::Vector{Int}, viz_type::String,
                                show_all::Bool, stdev::Bool, simulations_all::Array{Float64, 4})
    if !isdir("./figure/simulation/$viz_type")
        mkdir("./figure/simulation/$viz_type")
    end

    rc("figure", figsize = (20, 8))
    rc("font", family = "Arial")
    rc("mathtext", fontset = "custom")
    rc("mathtext", it = "Arial:italic")
    rc("font", size = 18)
    rc("axes", linewidth = 2)
    rc("xtick.major", width = 2)
    rc("ytick.major", width = 2)
    rc("xtick", direction = "in")
    rc("ytick", direction = "in")
    rc("lines", linewidth = 2.5)
    rc("lines", markersize = 12)

    cmap = ["mediumblue", "red"]
    shape = ["D", "s"]

    subplots_adjust(wspace=0.5, hspace=0.5)

    for (i, name) in enumerate(observables)
        subplot(2, 4, i)
        if show_all
            for j in eachindex(n_file)
                for l in eachindex(Sim.conditions)
                    plot(
                        Sim.t, simulations_all[i, j, :, l]./maximum(simulations_all[i, j, :, :]),
                        color=cmap[l], alpha=0.05
                    )
                end
            end
        end
        if viz_type != "average"
            for l in eachindex(Sim.conditions)
                plot(
                    Sim.t, Sim.simulations[i, :, l]/(maximum(Sim.simulations[i, :, :])),
                    color=cmap[l]
                )
            end
        else
            normalized = Array{Float64, 4}(
                undef, length(observables), length(n_file), length(Sim.t), length(Sim.conditions)
            )
            for j in eachindex(n_file)
                for l in eachindex(Sim.conditions)
                    normalized[i, j, :, l] = simulations_all[i, j, :, l]./maximum(simulations_all[i, j, :, :])
                end
            end
            for l in eachindex(Sim.conditions)
                plot(
                    Sim.t, [mean(filter(!isnan, normalized[i, :, k, l])) for k in eachindex(Sim.t)],
                    color=cmap[l]
                )
            end
            if stdev
                for l in eachindex(Sim.conditions)
                    ymean = [
                        mean(
                            filter(
                                !isnan, normalized[i, :, k, l]
                            )
                        ) for k in eachindex(Sim.t)
                    ]
                    yerr = [
                        std(
                            filter(
                                !isnan, normalized[i, :, k, l]
                            )
                        ) for k in eachindex(Sim.t)
                    ]
                    fill_between(
                        Sim.t, ymean - yerr, ymean + yerr,
                        lw=0, color=cmap[l], alpha=0.1
                    )
                end
            end
        end

        if isassigned(Exp.experiments, i)
            exp_t = Exp.get_timepoint(name)
            if isassigned(Exp.standard_error, i)
                for (l, condition) in enumerate(Sim.conditions)
                    if condition in keys(Exp.experiments[i])
                        exp_data = errorbar(
                            exp_t./60., Exp.experiments[i][condition],
                            yerr=Exp.standard_error[i][condition],
                            ecolor=cmap[l], elinewidth=1, capsize=8,
                            markerfacecolor="None", markeredgecolor=cmap[l],
                            fmt=shape[l], clip_on=false
                        )
                        for capline in exp_data[2]
                            capline.set_clip_on(false)
                        end
                        for barlinecol in exp_data[3]
                            barlinecol.set_clip_on(false)
                        end
                    end
                end
            else
                for (l, condition) in enumerate(Sim.conditions)
                    if condition in keys(Exp.experiments[i])
                        plot(
                            exp_t./60., Exp.experiments[i][condition], shape[l],
                            markerfacecolor="None", markeredgecolor=cmap[l],
                            clip_on=false
                        )
                    end
                end
            end
        end

        xlim(0, 90)
        xticks([0, 30, 60, 90])
        yticks([0, 0.3, 0.6, 0.9, 1.2])
        ylim(0, 1.2)
        xlabel("Time (min)")
        ylabel(get_ylabel(name))
    end

    savefig(
        "./figure/simulation/$viz_type.png", dpi=300, bbox_inches="tight"
    )
end

max_stars_repo_name: himoto/paramEstim

max_issues_repo_name: himoto/paramEstim

max_forks_repo_name: himoto/paramEstim