text:
#
# Tests for Organism related functions like lookAt(), toMem(),...
#
module TestCodeOrganism
  using FactCheck
  import Code
  import Creature
  import Helper
  import Mutator
  import Config
  import World
  include("Helper.jl")
  #
  # lookAt
  #
  facts("Testing lookAt()") do
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}())
    local org = d.orgs[1]

    code(d, :lookAt, org)
    @fact Code.eval(org.code)(d.cfg, org) --> 0.0

    val(org, Float16(100.0))
    @fact Code.eval(org.code)(d.cfg, org) --> 0.0

    val(org, Float16(0.0))
    @fact Code.eval(org.code)(d.cfg, org) --> 0.0

    val(org, Float16(5.0))
    World.setEnergy(d.man.world, Helper.Point(5,5), UInt16(10))
    @fact Code.eval(org.code)(d.cfg, org) --> 10.0

    Manager.destroy(d.man)
  end
  #
  # Eating
  #
  facts("Checking if left organism can eat right one") do
    local d     = create([Helper.Point(1,1), Helper.Point(2,1)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100,:orgStartEnergy=>5))
    local left  = d.orgs[1]
    local right = d.orgs[2]

    code(d, :eatRight, left)
    @fact left.energy === right.energy === 5 --> true
    consume(d.task)
    consume(d.task)
    @fact left.energy  --> 10
    @fact right.energy --> 0

    Manager.destroy(d.man)
  end
  facts("Checking if left and right organisms can eat eath other and step") do
    local d     = create([Helper.Point(1,1), Helper.Point(2,1)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100,:orgStartEnergy=>10))
    local left  = d.orgs[1]
    local right = d.orgs[2]

    code(d, :eatRight, left)
    code(d, :stepRight, left)
    code(d, :eatRight, left)
    code(d, :eatRight, left)

    code(d, :eatLeft, right)
    code(d, :stepLeft, right)
    code(d, :eatLeft, right)
    code(d, :eatLeft, right)

    @fact left.energy === right.energy === 10 --> true
    consume(d.task)
    consume(d.task)
    consume(d.task)
    @fact left.energy  --> 10
    @fact right.energy --> 10

    Manager.destroy(d.man)
  end
  facts("Checking if left organism can eat right one until death") do
    local d     = create([Helper.Point(1,1), Helper.Point(2,1)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100,:orgStartEnergy=>5))
    local left  = d.orgs[1]
    local right = d.orgs[2]

    code(d, :eatRight, left)
    code(d, :eatRight, left)
    code(d, :eatRight, left)
    @fact left.energy === right.energy === 5 --> true
    consume(d.task)
    consume(d.task)
    @fact left.energy  --> 10
    @fact right.energy --> 0

    Manager.destroy(d.man)
  end
  facts("Checking if left organism can donate right one") do
    local d     = create([Helper.Point(1,1), Helper.Point(2,1)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100))
    local left  = d.orgs[1]
    local right = d.orgs[2]

    code(d, :eatRight, left)
    val(left, Float16(-5.0))
    updateCode(left)

    @fact left.energy === right.energy === 100 --> true
    consume(d.task)
    consume(d.task)
    @fact left.energy --> 95
    @fact right.energy --> 105

    Manager.destroy(d.man)
  end
  facts("Checking if left organism can eat zero energy from right one") do
    local d = create([Helper.Point(1,1), Helper.Point(2,1)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100))
    local left = d.orgs[1]
    local right = d.orgs[2]

    code(d, :eatRight, left)
    val(left, Float16(0.0))

    @fact left.energy === right.energy === 100 --> true
    consume(d.task)
    consume(d.task)
    @fact left.energy --> 100
    @fact right.energy --> 100

    Manager.destroy(d.man)
  end
  facts("Checking if left organism eats right and right eats left") do
    local d = create([Helper.Point(1,1), Helper.Point(2,1)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100))
    local left = d.orgs[1]
    local right = d.orgs[2]

    code(d, :eatRight, left)
    code(d, :eatLeft, right)

    @fact left.energy === right.energy === 100 --> true
    consume(d.task)
    consume(d.task)
    consume(d.task)
    @fact left.energy --> 100
    @fact right.energy --> 100

    Manager.destroy(d.man)
  end

  facts("Checking if up organism can eat down one") do
    local d = create([Helper.Point(1,1), Helper.Point(1,2)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100))
    local up = d.orgs[1]
    local down = d.orgs[2]

    code(d, :eatDown, up)

    @fact up.energy === down.energy === 100 --> true
    consume(d.task)
    consume(d.task)
    @fact up.energy --> 105
    @fact down.energy --> 95

    Manager.destroy(d.man)
  end
  facts("Checking if up organism can donate down one") do
    local d = create([Helper.Point(1,1), Helper.Point(1,2)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100))
    local up = d.orgs[1]
    local down = d.orgs[2]

    code(d, :eatDown, up)
    val(up, Float16(-5.0))

    @fact up.energy === down.energy === 100 --> true
    consume(d.task)
    consume(d.task)
    @fact up.energy --> 95
    @fact down.energy --> 105

    Manager.destroy(d.man)
  end
  facts("Checking if up organism can eat zero energy from down one") do
    local d = create([Helper.Point(1,1), Helper.Point(1,2)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100))
    local up = d.orgs[1]
    local down = d.orgs[2]

    code(d, :eatDown, up)
    val(up, Float16(0.0))

    @fact up.energy === up.energy === 100 --> true
    consume(d.task)
    consume(d.task)
    @fact up.energy --> 100
    @fact down.energy --> 100

    Manager.destroy(d.man)
  end
  facts("Checking if up organism eats down and down eats up") do
    local d = create([Helper.Point(1,1), Helper.Point(1,2)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100))
    local up = d.orgs[1]
    local down = d.orgs[2]

    code(d, :eatDown, up)
    code(d, :eatUp, down)

    @fact up.energy === down.energy === 100 --> true
    consume(d.task)
    consume(d.task)
    consume(d.task)
    @fact up.energy --> 100
    @fact down.energy --> 100

    Manager.destroy(d.man)
  end

  facts("Checking if organism can eat empty space") do
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100))
    local up = d.orgs[1]

    code(d, :eatDown, up)

    @fact up.energy === 100 --> true
    consume(d.task)
    consume(d.task)
    @fact up.energy --> 100

    Manager.destroy(d.man)
  end
  facts("Checking if organism can eat out of borders") do
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100))
    local up = d.orgs[1]

    code(d, :eatUp, up)
    updateCode(up)

    @fact up.energy === 100 --> true
    consume(d.task)
    consume(d.task)
    @fact up.energy --> 100

    Manager.destroy(d.man)
  end
  facts("Checking if organism can eat energy") do
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100))
    local up = d.orgs[1]

    code(d, :eatRight, up)
    World.setEnergy(d.man.world, Helper.Point(2,1), UInt16(30))

    @fact up.energy === 100 --> true
    @fact World.getEnergy(d.man.world, Helper.Point(2,1)) --> 30
    consume(d.task)
    consume(d.task)
    @fact up.energy --> 105
    @fact World.getEnergy(d.man.world, Helper.Point(2,1)) --> 25
    consume(d.task)
    consume(d.task)
    @fact up.energy --> 115
    @fact World.getEnergy(d.man.world, Helper.Point(2,1)) --> 15
	  consume(d.task)
    @fact up.energy --> 120
    @fact World.getEnergy(d.man.world, Helper.Point(2,1)) --> 10

    Manager.destroy(d.man)
  end
  #
  # Moving
  #
  facts("Checking if organism can move left") do
    local d = create([Helper.Point(3,3)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100))
    local org = d.orgs[1]

    code(d, :stepLeft, org)
    consume(d.task)
    @fact World.getEnergy(d.man.world, Helper.Point(3,3)) --> UInt16(100)
    @fact World.getEnergy(d.man.world, Helper.Point(2,3)) --> UInt16(0)
    consume(d.task)
    @fact World.getEnergy(d.man.world, Helper.Point(3,3)) --> UInt16(0)
    @fact World.getEnergy(d.man.world, Helper.Point(2,3)) --> UInt16(100)

    code(d, :stepLeft, org)
    @fact World.getEnergy(d.man.world, Helper.Point(2,3)) --> UInt16(100)
    @fact World.getEnergy(d.man.world, Helper.Point(3,3)) --> UInt16(0)
    consume(d.task)
    consume(d.task)
    @fact World.getEnergy(d.man.world, Helper.Point(1,3)) --> UInt16(100)
    @fact World.getEnergy(d.man.world, Helper.Point(2,3)) --> UInt16(0)
    @fact World.getEnergy(d.man.world, Helper.Point(3,3)) --> UInt16(0)

    Manager.destroy(d.man)
  end
  facts("Checking if organism can move right") do
    local d = create([Helper.Point(3,3)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100))
    local org = d.orgs[1]

    code(d, :stepRight, org)
    consume(d.task)
    @fact World.getEnergy(d.man.world, Helper.Point(3,3)) --> UInt16(100)
    @fact World.getEnergy(d.man.world, Helper.Point(4,3)) --> UInt16(0)
    consume(d.task)
    @fact World.getEnergy(d.man.world, Helper.Point(3,3)) --> UInt16(0)
    @fact World.getEnergy(d.man.world, Helper.Point(4,3)) --> UInt16(100)

    code(d, :stepRight, org)
    updateCode(org)
    @fact World.getEnergy(d.man.world, Helper.Point(4,3)) --> UInt16(100)
    @fact World.getEnergy(d.man.world, Helper.Point(3,3)) --> UInt16(0)
    consume(d.task)
    consume(d.task)
    @fact World.getEnergy(d.man.world, Helper.Point(5,3)) --> UInt16(100)
    @fact World.getEnergy(d.man.world, Helper.Point(4,3)) --> UInt16(0)
    @fact World.getEnergy(d.man.world, Helper.Point(3,3)) --> UInt16(0)

    Manager.destroy(d.man)
  end
  facts("Checking if organism can move up") do
    local d = create([Helper.Point(3,3)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100))
    local org = d.orgs[1]

    code(d, :stepUp, org)
    consume(d.task)
    @fact World.getEnergy(d.man.world, Helper.Point(3,3)) --> UInt16(100)
    @fact World.getEnergy(d.man.world, Helper.Point(3,2)) --> UInt16(0)
    consume(d.task)
    @fact World.getEnergy(d.man.world, Helper.Point(3,3)) --> UInt16(0)
    @fact World.getEnergy(d.man.world, Helper.Point(3,2)) --> UInt16(100)

    code(d, :stepUp, org)
    updateCode(org)
    @fact World.getEnergy(d.man.world, Helper.Point(3,2)) --> UInt16(100)
    @fact World.getEnergy(d.man.world, Helper.Point(3,3)) --> UInt16(0)
    consume(d.task)
    consume(d.task)
    @fact World.getEnergy(d.man.world, Helper.Point(3,1)) --> UInt16(100)
    @fact World.getEnergy(d.man.world, Helper.Point(3,2)) --> UInt16(0)
    @fact World.getEnergy(d.man.world, Helper.Point(3,3)) --> UInt16(0)

    Manager.destroy(d.man)
  end
  facts("Checking if organism can move down") do
    local d = create([Helper.Point(3,3)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100))
    local org = d.orgs[1]

    code(d, :stepDown, org)
    consume(d.task)
    @fact World.getEnergy(d.man.world, Helper.Point(3,3)) --> UInt16(100)
    @fact World.getEnergy(d.man.world, Helper.Point(3,4)) --> UInt16(0)
    consume(d.task)
    @fact World.getEnergy(d.man.world, Helper.Point(3,3)) --> UInt16(0)
    @fact World.getEnergy(d.man.world, Helper.Point(3,4)) --> UInt16(100)

    code(d, :stepDown, org)
    @fact World.getEnergy(d.man.world, Helper.Point(3,4)) --> UInt16(100)
    @fact World.getEnergy(d.man.world, Helper.Point(3,3)) --> UInt16(0)
    consume(d.task)
    consume(d.task)
    @fact World.getEnergy(d.man.world, Helper.Point(3,5)) --> UInt16(100)
    @fact World.getEnergy(d.man.world, Helper.Point(3,4)) --> UInt16(0)
    @fact World.getEnergy(d.man.world, Helper.Point(3,3)) --> UInt16(0)

    Manager.destroy(d.man)
  end
  # TODO: same tests for right, up, down
  facts("Checking if organism can move left out of the border in cyclic mode") do
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100, :worldCyclical=>true))
    local org = d.orgs[1]

    code(d, :stepLeft, org)
    consume(d.task)
    @fact World.getEnergy(d.man.world, Helper.Point(1,1)) --> UInt16(100)
    @fact World.getEnergy(d.man.world, Helper.Point(10,1)) --> UInt16(0)
    consume(d.task)
    @fact World.getEnergy(d.man.world, Helper.Point(1,1)) --> UInt16(0)
    @fact World.getEnergy(d.man.world, Helper.Point(10,1)) --> UInt16(100)

    Manager.destroy(d.man)
  end
  facts("Checking if organism can't move left, because of energy block") do
    local d = create([Helper.Point(2,1)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100))
    local org = d.orgs[1]

    code(d, :stepLeft, org)
    World.setEnergy(d.man.world, Helper.Point(1,1), UInt16(123))
    consume(d.task)
    @fact World.getEnergy(d.man.world, Helper.Point(2,1)) --> UInt16(100)
    @fact World.getEnergy(d.man.world, Helper.Point(1,1)) --> UInt16(123)
    consume(d.task)
    @fact World.getEnergy(d.man.world, Helper.Point(2,1)) --> UInt16(100)
    @fact World.getEnergy(d.man.world, Helper.Point(1,1)) --> UInt16(123)

    Manager.destroy(d.man)
  end
  facts("Checking if organism can move left out of the border in cyclic mode and can't, because of energy block") do
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100, :worldCyclical=>true))
    local org = d.orgs[1]

    code(d, :stepLeft, org)
    #
    # Sets energy on the right side of the world
    #
    World.setEnergy(d.man.world, Helper.Point(10,1), UInt16(123))
    consume(d.task)
    @fact World.getEnergy(d.man.world, Helper.Point(1,1)) --> UInt16(100)
    @fact World.getEnergy(d.man.world, Helper.Point(10,1)) --> UInt16(123)
    consume(d.task)
    @fact World.getEnergy(d.man.world, Helper.Point(1,1)) --> UInt16(100)
    @fact World.getEnergy(d.man.world, Helper.Point(10,1)) --> UInt16(123)
    #
    # Removes energy on the right side of the world
    #
    World.setEnergy(d.man.world, Helper.Point(10,1), UInt16(0))
    consume(d.task)
    consume(d.task)
    @fact World.getEnergy(d.man.world, Helper.Point(1,1)) --> UInt16(0)
    @fact World.getEnergy(d.man.world, Helper.Point(10,1)) --> UInt16(100)

    Manager.destroy(d.man)
  end
  facts("Checking if organism can move left out of the border in non cyclic mode") do
    local d = create([Helper.Point(1,1)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100, :worldCyclical=>false))
    local org = d.orgs[1]

    code(d, :stepLeft, org)

    consume(d.task)
    @fact World.getEnergy(d.man.world, Helper.Point(1,1)) --> UInt16(100)
    @fact World.getEnergy(d.man.world, Helper.Point(10,1)) --> UInt16(0)
    consume(d.task)
    consume(d.task)
    @fact World.getEnergy(d.man.world, Helper.Point(1,1)) --> UInt16(100)
    @fact World.getEnergy(d.man.world, Helper.Point(10,1)) --> UInt16(0)

    Manager.destroy(d.man)
  end
  #
  # Memory
  #
  facts("Checking if organism can save something to internal memory") do
    local d   = create([Helper.Point(1,1)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100))
    local org = d.orgs[1]

    code(d, :toMem, org)
    consume(d.task)
    @fact haskey(org.mem, val(org)) --> false
    consume(d.task)
    @fact haskey(org.mem, val(org)) --> true

    Manager.destroy(d.man)
  end
  facts("Checking if organism can load something from internal memory") do
    local d   = create([Helper.Point(1,1)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100))
    local org = d.orgs[1]

    code(d, :fromMem, org)
    org.mem[val(org)] = 12.0
    consume(d.task)
    consume(d.task)
    @fact Code.eval(org.code)(d.cfg, org) --> 12.0

    Manager.destroy(d.man)
  end
  # #
  # # id
  # #
  # facts("Checking if organism can check id of other organism") do
  #   local d     = create([Helper.Point(1,1), Helper.Point(2,1)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100))
  #   local left  = d.orgs[1]
  #   local right = d.orgs[2]
  #
  #   code(d, :idLeft, right)
  #   consume(d.task)
  #   consume(d.task)
  #   @fact Code.eval(right.code)(d.cfg, right) --> 1.0
  #
  #   Manager.destroy(d.man)
  # end
  # facts("Checking if organism can check id of nothing") do
  #   local d   = create([Helper.Point(1,1)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100))
  #   local org = d.orgs[1]
  #
  #   code(d, :idRight, org)
  #   consume(d.task)
  #   consume(d.task)
  #   @fact Code.eval(org.code)(d.cfg, org) --> 0.0
  #
  #   Manager.destroy(d.man)
  # end
  # facts("Checking if organism can check id out of borders") do
  #   local d   = create([Helper.Point(1,1)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100))
  #   local org = d.orgs[1]
  #
  #   code(d, :idLeft, org)
  #   consume(d.task)
  #   consume(d.task)
  #   @fact Code.eval(org.code)(d.cfg, org) --> 0.0
  #
  #   Manager.destroy(d.man)
  # end
  # facts("Checking if organism can check id of energy block") do
  #   local d   = create([Helper.Point(1,1)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100))
  #   local org = d.orgs[1]
  #
  #   World.setEnergy(d.man.world, Helper.Point(2,1), UInt16(10))
  #   code(d, :idRight, org)
  #   consume(d.task)
  #   consume(d.task)
  #   @fact Code.eval(org.code)(d.cfg, org) --> 0.0
  #
  #   Manager.destroy(d.man)
  # end
  #
  # Check energy
  #
  facts("Checking if organism can get energy of other organism") do
    local d     = create([Helper.Point(1,1), Helper.Point(2,1)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100))
    local left  = d.orgs[1]
    local right = d.orgs[2]

    code(d, :energyLeft, right)
    consume(d.task)
    consume(d.task)
    @fact Code.eval(right.code)(d.cfg, right) --> 100.0

    Manager.destroy(d.man)
  end
  facts("Checking if organism can get energy of nothing") do
    local d   = create([Helper.Point(1,1)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100))
    local org = d.orgs[1]

    code(d, :energyRight, org)
    consume(d.task)
    consume(d.task)
    @fact Code.eval(org.code)(d.cfg, org) --> 0.0

    Manager.destroy(d.man)
  end
  facts("Checking if organism can get energy out of borders") do
    local d   = create([Helper.Point(1,1)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100))
    local org = d.orgs[1]

    code(d, :energyLeft, org)
    consume(d.task)
    consume(d.task)
    @fact Code.eval(org.code)(d.cfg, org) --> 0.0

    Manager.destroy(d.man)
  end
  facts("Checking if organism can get energy from real energy block") do
    local d   = create([Helper.Point(1,1)], Dict{Symbol, Any}(:orgEnergySpendPeriod=>100))
    local org = d.orgs[1]

    World.setEnergy(d.man.world, Helper.Point(2,1), UInt16(10))
    code(d, :energyRight, org)
    consume(d.task)
    consume(d.task)
    @fact Code.eval(org.code)(d.cfg, org) --> 0.0

    Manager.destroy(d.man)
  end
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
module PhotonicCrystal

using ...Bravais,
...Shapes,
...DielectricFunctions,
...DefineSimulation,
...BoundaryConditions,
..ConstructionToolsBase

export build_pc_domain,
line_defect

"""
    sim = Simulation(domain; bnd=Boundary(bc=:p), dis, sct, tls, Na=1, Nb=1)

Simulation from a single domain.

Intended for case where `domain` has a finite lattice structure, in which case
    it defines the simulation to be `Na`×`Nb` unit cells.

Care must be taken with boundary layer specifications. It is recommended that
`bnd` be left unspecified.
"""
function DefineSimulation.Simulation(dom::Domain; bnd::Boundary=Boundary(bc=PeriodicBC), dis::Discretization, sct::Scattering=Scattering(), tls::TwoLevelSystem=TwoLevelSystem(), Na::Int=1, Nb::Int=1)
    lattice = dom.lattice

    # bnd = deepcopy(bnd)
    if !isinf(lattice.a)
        if typeof(bnd.bc[1][1])<:PeriodicBC
            ∂Ω1 = (lattice.x0,
                    lattice.x0+Na*lattice.a*(cos(lattice.α)-sin(lattice.α)*cot(lattice.β)))
        else
            ∂Ω1 = bnd.∂Ω[1]
        end
    end
    if !isinf(lattice.b)
        if typeof(bnd.bc[2][1])<:PeriodicBC
            ∂Ω2 = (lattice.y0,
                    lattice.y0+Nb*lattice.b*sin(lattice.β))
        else
            ∂Ω2 = bnd.∂Ω[2]
        end
    end
    return deepcopy(Simulation(System(dom), Boundary(bnd; ∂Ω=(∂Ω1,∂Ω2)), dis, tls, BravaisLattice(lattice; :a=>Na*lattice.a, :b=>Nb*lattice.b)))
end
"""
    sim = Simulation(domain, sim; Na=1, Nb=1)

Build simulation out of a single `domain`, with fields other than `bnd` the same
as in `sim`, except it is `Na`×`Nb` unit cells.

Intended for case where `domain` has a finite lattice structure.
"""
DefineSimulation.Simulation(dom::Domain, sim::Simulation; Na::Int=1, Nb::Int=1) = Simulation(dom; dis=sim.dis, tls=sim.tls, Na=Na, Nb=Nb)
"""
    sim = Simulation(domain_num, sim; Na=1, Nb=1)

Build simulation out of domain number `domain_num`, with fields other than `bnd`
the same as in `sim`, except it is `Na`×`Nb` unit cells.

Intended for case where specified domain has a finite lattice structure.
"""
Simulation(domain_num::Int, sim::Simulation; Na::Int=1, Nb::Int=1) = deepcopy(Simulation(sim.sys.domains[domain_num]; dis=sim.dis, tls=sim.tls, Na=Na, Nb=Nb))
"""
    domain = build_pc_domain
"""
function build_pc_domain(
            lattice::BravaisLattice,
            subdomains::Tsd = Subdomains();
            is_in_domain::Tsh = Universe(),
            domain_params::Td = Dict(:n₁ => 1, :n₂ => 0, :F => 0),
            domain_type::Symbol = :pc,
            domain_ε::Function = piecewise_constant_ε,
            domain_F::Function = piecewise_constant_F,
            which_asymptote::Symbol = :none,
            which_waveguide::Int = 0
            ) where Tsd<:Subdomains where Tsh<:AbstractShape where Td<:Dict

    return build_domain(
                subdomains;
                is_in_domain=is_in_domain,
                domain_params=domain_params,
                domain_type=domain_type,
                domain_ε=domain_ε,
                domain_F=domain_F,
                lattice=lattice,
                which_asymptote=which_asymptote,
                which_waveguide=which_waveguide
                )
end

##########################################################################################
### DEFECT CONSTRUCTION
##########################################################################################
"""
    line_defect(pc_domain, direction, N_start, N_stop, N_transverse, params, ε, F, subdomains; N_width=1)
"""
function line_defect(
    pc_domain::Tdom,
    direction::Symbol,
    N_start::Int,
    N_stop::Int,
    N_transverse::Int,
    params=pc_domain.domain_params,
    ε::Function=pc_domain.domain_ε,
    F::Function=pc_domain.domain_F,
    subdomains::Tsd=Subdomains();
    N_width::Int=1
    ) where Tdom<:Domain where Tsd<:Subdomains

    @assert direction ∈ [:a, :b] "invalid direction $(direction)."
    if direction == :a
        Na_start = N_start
        Na_stop = N_stop
        Nb_start = N_transverse
        Nb_stop = N_transverse+N_width-1
    else direction == :b
        Nb_start = N_start
        Nb_stop = N_stop
        Na_start = N_transverse
        Na_stop = N_transverse+N_width-1
    end
    lattice = pc_domain.lattice

    if Na_start==Na_stop && Nb_start==Nb_stop
        defect_type = :site_defect
    else
        defect_type = :line_defect
    end

    a, b = lattice.a*(Na_stop-Na_start+1), lattice.a*(Nb_stop-Nb_start+1)
    α, θ = lattice.β-lattice.α, lattice.α
    x0 = lattice.x0 + Na_start*lattice.a*lattice.v1[1] + lattice.x0 + Nb_start*lattice.b*lattice.v2[1]
    y0 = lattice.y0 + Na_start*lattice.a*lattice.v1[2] + lattice.x0 + Nb_start*lattice.b*lattice.v2[2]

    return build_domain(subdomains;
                is_in_domain = Parallelogram(a,b,α,x0,y0,θ),
                domain_params = params,
                domain_type = defect_type,
                domain_ε = ε,
                domain_F = F,
                lattice = lattice)
end

"""
    site_defect
"""
function site_defect(pc_domain::Domain, Na::Int, Nb::Int,
        params=pc_domain.domain_params, ε=pc_domain.domain_ε, F=pc_domain.domain_F,
        sudomains::Subdomains=Subdomains() )

    defect_domain = line_defect_domain(pc_domain, :a, Na, Na, Nb, params, ε, F, sudomains)
    return defect_domain
end
function site_defect(pc_domain::Domain, Na::Int, Nb::Array{Int},
    params=pc_domain.domain_params, ε=pc_domain.domain_ε, F=pc_domain.domain_F,
    sudomains::Subdomains=Subdomains() )

    return defect_domain.(Ref(pc_domain), Na, Nb, Ref(params), Ref(ε), Ref(F), Ref(sudomains))
end
function site_defect(pc_domain::Domain, Na::Array{Int}, Nb,
    params=pc_domain.domain_params, ε=pc_domain.domain_ε, F=pc_domain.domain_F,
    sudomains::Subdomains=Subdomains())

    return defect_domain.(Ref(pc_domain), Na, Nb, Ref(params), Ref(ε), Ref(F), Ref(sudomains))
end

#################################################################################################################
### SHAPE OVERLOADING
#################################################################################################################
"""
    Circle(R, lattice, referece, params[, subdomains; x0, y0, type, ε, F])
"""
function Shapes.Circle(R::Real, reference::Symbol,
            params::Dict,
            domain::Domain;
            x0::Real=0,
            y0::Real=0,
            type::Symbol=:background,
            ε::Function = piecewise_constant_ε,
            F::Function = piecewise_constant_F)

    x0_new, y0_new = get_reference(reference, x0, y0, domain.lattice)
    for i ∈ eachindex(x0_new)
        domain = Circle(R, x0_new[i], y0_new[i], params, domain; type=type, ε=ε, F=F)
    end
    return domain
end

"""
    Ellipse(a, b, θ, lattice, referece, params[, subdomains; x0, y0, type, ε, F])
"""
function Shapes.Ellipse(a::Real, b::Real, θ::Real, reference::Symbol,
            params::Dict,
            domain::Domain;
            x0::Number=0,
            y0::Number=0,
            type::Symbol=:background,
            ε::Function = piecewise_constant_ε,
            F::Function = piecewise_constant_F)

    x0_new, y0_new = get_reference(reference, x0, y0, domain.lattice)
    for i ∈ eachindex(x0_new)
        domain = Ellipse(a, b, x0_new[i], y0_new[i], θ, params, domain; type=type, ε=ε, F=F)
    end
    return domain
end

"""
    Square(a, θ, lattice, referece, params[, subdomains; x0, y0, type, ε, F])
"""
function Shapes.Square(a::Real, θ::Real, reference::Symbol,
            params::Dict,
            domain::Domain;
            x0::Number=0,
            y0::Number=0,
            type::Symbol=:background,
            ε::Function = piecewise_constant_ε,
            F::Function = piecewise_constant_F)

    x0_new, y0_new = get_reference(reference, x0, y0, domain.lattice)
    for i ∈ eachindex(x0_new)
        domain = Square(a, x0_new[i], y0_new[i], θ, params, domain; type=type, ε=ε, F=F)
    end
    return domain
end

"""
    Rectangle(a, b, θ, lattice, referece, params[, subdomains; x0, y0, type, ε, F])
"""
function Shapes.Rectangle(a::Real, b::Real, θ::Real, reference::Symbol,
            params::Dict,
            domain::Domain;
            x0::Number=0,
            y0::Number=0,
            type::Symbol=:background,
            ε::Function = piecewise_constant_ε,
            F::Function = piecewise_constant_F)

    x0_new, y0_new = get_reference(reference, x0, y0, domain.lattice)
    for i ∈ eachindex(x0_new)
        domain = Rectangle(a, b, x0_new[i], y0_new[i], θ, params, domain; type=type, ε=ε, F=F)
    end
    return domain
end

"""
    Parallelogram(a, b, α, θ, lattice, referece, params[, subdomains; x0, y0, type, ε, F])
"""
function Shapes.Parallelogram(a::Real, b::Real, α::Real, θ::Real, reference::Symbol,
            params::Dict,
            domain::Domain;
            x0::Number=0,
            y0::Number=0,
            type::Symbol=:background,
            ε::Function = piecewise_constant_ε,
            F::Function = piecewise_constant_F)

    x0_new, y0_new = get_reference(reference, x0, y0, domain.lattice)
    for i ∈ eachindex(x0_new)
        domain = Parallelogram(a, b, α, x0_new[i], y0_new[i], θ, params, domain; type=type, ε=ε, F=F)
    end
    return domain
end

"""
    get_reference(reference::Symbol, x0, y0, lattice::BravaisLattice)
"""
function get_reference(reference::Symbol, x0::Real, y0::Real, lattice::BravaisLattice)
    x0_new, y0_new = [], []
    v1, v2 = lattice.v1, lattice.v2
    a, b = lattice.a, lattice.b
    if reference == :center
        push!(x0_new, x0 - v1[1]*a/2 + v2[1]*b/2)
        push!(y0_new, y0 - v1[2]*a/2 + v2[2]*b/2)

        push!(x0_new, x0 + v1[1]*a/2 - v2[1]*b/2)
        push!(y0_new, y0 + v1[2]*a/2 - v2[2]*b/2)

        push!(x0_new, x0 + 3v1[1]*a/2 + v2[1]*b/2)
        push!(y0_new, y0 + 3v1[2]*a/2 + v2[2]*b/2)

        push!(x0_new, x0 + v1[1]*a/2 + 3v2[1]*b/2)
        push!(y0_new, y0 + v1[2]*a/2 + 3v2[2]*b/2)

        push!(x0_new, x0 + v1[1]*a/2 + v2[1]*b/2)
        push!(y0_new, y0 + v1[2]*a/2 + v2[2]*b/2)
    elseif reference ∈ [:corner, :corners]
        push!(x0_new, v1[1]*a)
        push!(y0_new, v1[2]*a)

        push!(x0_new, v2[1]*b)
        push!(y0_new, v2[2]*b)

        push!(x0_new, v1[1]*a + v2[1]*b)
        push!(y0_new, v1[2]*a + v2[2]*b)

        push!(x0_new, 0)
        push!(y0_new, 0)
    elseif reference ∈ [:bottom_edge, :top_edge, :bottom, :top]
        push!(x0_new, v1[1]*a/2 + v2[1]*b)
        push!(y0_new, v1[2]*a/2 + v2[2]*b)

        push!(x0_new, v1[1]*a/2)
        push!(y0_new, v1[2]*a/2)
    elseif reference ∈ [:left_edge, :right_edge, :left, :right]
        push!(x0_new, v1[1]*a + v2[1]*b/2)
        push!(y0_new, v1[2]*a + v2[2]*b/2)

        push!(x0_new, v2[1]*b/2)
        push!(y0_new, v2[2]*b/2)
    elseif reference ∈ [:southwest, :sw]
        push!(x0_new, x0)
        push!(y0_new, y0)
    elseif reference ∈ [:south, :s]
        push!(x0_new, x0 + v1[1]*a/2)
        push!(y0_new, y0 + v1[1]*a/2)
    elseif reference ∈ [:southeast, :se]
        push!(x0_new, x0 + v1[1]*a)
        push!(y0_new, y0 + v1[2]*a)
    elseif reference ∈ [:east, :e]
        push!(x0_new, x0 + v1[1]*a + v2[1]*b/2)
        push!(y0_new, y0 + v1[2]*a + v2[2]*b/2)
    elseif reference ∈ [:northeast, :ne]
        push!(x0_new, x0 + v1[1]*a + v2[1]*b)
        push!(y0_new, y0 + v1[2]*a + v2[2]*b)
    elseif reference ∈ [:north, :n]
        push!(x0_new, x0 + v1[1]*a/2 + v2[1]*b)
        push!(y0_new, y0 + v1[2]*a/2 + v2[2]*b)
    elseif reference ∈ [:northwest, :nw]
        push!(x0_new, x0 + v2[1]*b)
        push!(y0_new, y0 + v2[2]*b)
    elseif reference ∈ [:west, :w]
        push!(x0_new, x0 + v2[1]*b/2)
        push!(y0_new, y0 + v2[2]*b/2)
    else
        throw(ArgumentError("invalid unit cell position reference $(reference).
        Should be one of :center, :corner or a compass direction"))
    end
    return x0_new, y0_new
end

end #module

ext: jl


text:
using DashBootstrapComponents, DashHtmlComponents

buttons = html_div([
    dbc_button("Large button", size="lg", className="mr-1"),
    dbc_button("Regular button", className="mr-1"),
    dbc_button("Small button", size="sm"),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
# Compare benchmarks
using GPUBenchmarks, BenchmarkTools
db = GPUBenchmarks.get_database()
using GPUBenchmarks: name, timestamp

ts = last(sort(unique(timestamp.(db))))

newbenchs = filter(x-> x.timestamp == ts, db)

unique(GPUBenchmarks.device.(db))

oldbenchs = filter(x-> x.device == "opencl", db)
old_newest_ts = last(sort(unique(timestamp.(oldbenchs))))
oldbenchs = filter(x-> x.timestamp == old_newest_ts, oldbenchs)

sort(unique(GPUArrays.name.(db))

for elem in unique(GPUBenchmarks.name.(newbenchs))
    A = filter(x-> x.name == elem, oldbenchs)
    B = filter(x-> x.name == elem, oldbenchs)
    for (a, b) in zip(A, B)
        x = judge(minimum(a.benchmark), minimum(b.benchmark))
        if x.time != :invariant
            println(name)
            @show x
        end
    end
end


text:
include("header.jl")
using Statistics

@testset "core" begin

    # Issue #71.1: AutoGrad error in backprop when iterating dict
    u = [rand(2,3), rand(2)]
    v = [rand(1,2), rand(1)]
    m = Dict(:u=>u, :v=>v)
    x,y = rand(3,4),rand(1,4)
    pred(m,x) = foldl((x,w)->w[1]*x .+ w[2], [m[:u],m[:v]], init=x)
    loss(m,x,y) = mean(abs2, pred(m,x)-y)
    ∇ = grad(loss)
    @test isa(loss(m,x,y), Real)
    @test isa(∇(m,x,y), Dict)
    l2(ws) = mean(mean.(abs2, ws))
    loss(m,x,y) = mean(abs2, pred(m,x)-y) + mean(l2.(collect(values(m))))
    ∇ = grad(loss)
    @test isa(loss(m,x,y), Real)
    @test isa(∇(m,x,y), Dict)

    # Issue #71.2
    m = [rand(1,3), rand(1)]
    x,y = rand(3,4),rand(1,4)
    pred(m,x) = m[1]*x .+ m[2]
    loss(m,x,y) = mean(abs2, pred(m,x) .- y)
    ∇ = grad(loss)
    @test isa(∇(m,x,y), Array)
    loss(m,x,y) = mean(abs2, pred(m,x) .- y) + mean(mean.(abs2, m))
    ∇ = grad(loss)
    @test isa(∇(m,x,y), Array)

    # differentiate vs @diff with expression arguments
    x = Param(1)
    @test grad((@diff sin(sqrt(x))),x) == cos(x)/2
    @test grad(AutoGrad.differentiate(sin,sqrt(x)),x) != cos(x)/2

    # Issue #101.1
    x = Param(1.0); f1(x)=x
    @test grad((@diff f1(x)), x) == 1

    # Issue #101.2
    x = Param([1.,2.]); f2(x)=1x
    @test grad((@diff sum(f2.(x))), x) == [1.0,1.0]

    # Issue #103
    f(x) = exp(x)
    @test grad(x->f(x))(1.) == exp(1.)
    @test grad(x -> sum(f.(x)))([1.]) == [ exp(1.) ]

    # Double broadcasting
    x = Param([1.,2.]); f3(x)=sin(x); f4(x)=sin.(x)
    @test grad((@diff sum(f3.(x))), x) == grad((@diff sum(f4.(x))), x) == grad((@diff sum(f4(x))), x)

    # array-scalar mul
    a = Param(rand(2,3)); s = Param(rand())
    @test @gcheck sum(a .* s)
    @test @gcheck sum(a * s)

    # @zerograd needs to handle Bcasted
    x = Param(rand(2,3))
    f(x) = sign(x)
    @test @gcheck sum(f.(x))

    # result may not always be last on tape
    x = Param(rand(2,3))
    f(x) = (x1=sum(x); x2=2x; x1)
    @test @gcheck f(x)

    # Issue #106:
    h(x) = exp(-x); h′(x,y) = -y
    𝓁(x,y) = sum(abs2,x-y)/2
    function neural_net(mparams, input; h=h, h′=h′, N=length(mparams))
        δ = [];
        X = Any[input];
        for i=1:N
            x = sum(mparams[i] .* [X[i],1])
            y = h.(x)
            push!(δ, h′.(x,y))
            push!(X,y)
        end
        return X,δ
    end
    mparams =[[randn(),randn()] for i=1:3]
    P = Param(mparams)
    loss(P,x,y)= 𝓁(neural_net(P,x)[1][end],y)
    x,y=randn(),randn()
    J = @diff loss(P,x,y)
    @test isa(J, AutoGrad.Tape)
    @test_broken @gcheck loss(P,x,y)

end


text:
"""
    UnnormalizedCategorical

**NOTE:** This assumes (and does not verify) that `p` is positive.
"""
struct UnnormalizedCategorical{T} <: Distribution{Univariate,Discrete}
    p::T
    s::Float64
end

function UnnormalizedCategorical(p::Array{Float64})
    UnnormalizedCategorical(p, sum(p))
end

function rand(d::UnnormalizedCategorical)
    draw = Base.rand() * d.s
    cp = 0.0
    i = 0
    while cp < draw
        cp += d.p[i += 1]
    end
    Tuple(CartesianIndices(size(d.p))[max(1,i)])
end

rand(d::UnnormalizedCategorical, n::Int) = map(_ -> rand(d), 1:n)

max_stars_repo_name: maxmouchet/Statistikk

max_issues_repo_name: maxmouchet/Statistikk

max_forks_repo_name: maxmouchet/Statistikk


text:
import Base.Threads
import LinearAlgebra

"""
    get_frame(frame_idx::Int64, trajectory::Trajectory)
    ::Frame

Return Frame object which correspond to `framd_idx` frame from `trajectory`.
"""
function get_frame(frame_idx::IntT, trj::Trajectory{RealT}
    )::Frame where {IntT <: Integer, RealT <: Real}

    Frame(trj.coordinates[:,frame_idx], trj.attributes)
end

"""
    get_atom(query::Int64, trajectory::Trajectory)
    ::Vector{Atom}

Return Vector of `query`th Atom object from `trajectory`.
"""
function get_atom(query::IntT, trj::Trajectory)::Vector{Atom} where IntT <: Integer
    map(xyz -> Atom(xyz, trj.attributes[query]), trj.coordinates[query, :])
end

"""
    get_atom(query::Int64, frame::Frame)
    ::Atom

Return `query`th Atom object from `frame`.
"""
function get_atom(query::IntT, frame::Frame)::Atom where IntT <: Integer
    Atom(frame.coordinates[query], frame.attributes[query])
end

"""
    clip_trajectory(query::Union{Integer, Vector{IntT}, OrdinalRange}, trj::Trajectory;
                    query_key::Symbol = :frame)
    ::Trajecotory where IntT <: Integer

Clip a part of trajectory you specified.
Expected query key is one of the following

- `:frame` : in this case, query specify the frame index.
- `:atom`  : in this case, query specify the atom index.
"""
function clip_trajectory(query::IntT, trj::Trajectory{RealT};
    query_key = query_key::Symbol = :frame
    )::Trajectory{RealT} where {IntT <: Integer, RealT <: Real}

    if query_key == :frame
        Trajectory(reshape(trj.coordinates[:, query], (trj.natom, 1)), trj.attributes)
    elseif query_key == :atom
        Trajectory(reshape(trj.coordinates[query, :], (1, trj.nframe)), [trj.attributes[query]])
    else
        throw(ArgumentError("""
                            clip_trajectory: invalid query key :$(query_key) here
                            expected query key is one of the following.
                            - :frame
                            - :atom
                            """))
    end
end

function clip_trajectory(query::Union{Vector{IntT}, OrdinalRange}, trj::Trajectory{RealT};
    query_key::Symbol = :frame
    )::Trajectory{RealT} where {IntT <: Integer, RealT <: Real}

    if query_key == :frame
        Trajectory(trj.coordinates[:, query], trj.attributes)
    elseif query_key == :atom
        Trajectory(trj.coordinates[query, :], trj.attributes[query])
    else
        throw(ArgumentError("""
                            clip_trajectory: invalid query key :$(query_key) here
                            expected query key is one of the following.
                            - :frame
                            - :atom
                            """))
    end
end

"""
    geometric_center_of_mass(coordinates::Vector{Coordinate};
                             atom_indices::Union{Vector, OrdinalRange, Colon} = :)
    ::Coordinate

Calculate the geometric center of mass of Coordinate vector for specified atom indices.

> For trajectory case, geometric center of mass is calculated by `center_of_mass` method with geometric flag.
"""
function geometric_center_of_mass(coordinates::Vector{<:Coordinate{RealT}};
    atom_indices::Union{Vector, OrdinalRange, Colon} = :)::Coordinate{RealT} where RealT <: Real

    selected_coord = view(coordinates, atom_indices)
    sum(selected_coord) / length(selected_coord)
end

"""
    center_of_mass(coordinates::Vector{Coordiante}, mass_vec{Real};
                   atom_indices::Union{Vector, OrdinalRange, Colon} = :)
    ::Coordinate{RealT}

Calculate the (non-geometric) center of mass of Coordinate vector for specified atom indices.
You have to pass mass information vector which have same length to coordinates vector.
"""
function center_of_mass(coordinates::Vector{<:Coordinate{RealT}}, mass_vec::Vector{<:Real};
    atom_indices::Union{Vector, OrdinalRange, Colon} = :)::Coordinate{RealT} where RealT <: Real

    if length(coordinates) != length(mass_vec)
        throw(AssertionError("""
                             coordinate vector and mass vector should be same length.
                             """))
    end

    selected_coord = view(coordinates, atom_indices)
    selected_mass  = view(mass_vec,    atom_indices)
    sum(selected_coord .* selected_mass) / sum(selected_mass)
end

"""
    center_of_mass(query::Trajectory;
                   frame_indices::Union{Vector, OrdinalRange, Colon} = :,
                   atom_indices::Union{Vector, OrdinalRange, Colon} = :,
                   geometric::Bool = false)
    ::Vector{Coordinate{Real}}

Calculate the center of mass of trajectory for specified frame and atom indices.
If you set `geometric` is `true`, this function calculate geometric center of mass.
"""
function center_of_mass(trj::Trajectory{RealT};
    frame_indices::Union{Vector, OrdinalRange, Colon} = :,
    atom_indices::Union{Vector, OrdinalRange, Colon} = :,
    geometric::Bool = false)::Vector{Coordinate{RealT}} where RealT <: Real

    if !geometric
        mass_vec = map(attributes -> attributes.mass, trj.attributes)
        selected_mass_vec = view(mass_vec, atom_indices)
        if nothing ∈ selected_mass_vec || length(selected_mass_vec) == 0
            throw(ArgumentError("""
                                center_of_mass: geometric flag is false but attributes of trajectory have Nothing value.
                                you can not use mass to calculation.
                                """))
        end
        total_mass = sum(selected_mass_vec)
        weited_sum = reshape(selected_mass_vec, (1, length(selected_mass_vec))) * view(trj.coordinates, atom_indices, frame_indices)
        vec(weited_sum / total_mass)
    else
        selected_coord = view(trj.coordinates, atom_indices, frame_indices)
        vec(sum(selected_coord, dims = 1) / size(selected_coord, 1))
    end
end

"""
    pair_length_matrix(corrds1::Vector{Coordinate}, coords2::Vector{Coordinate})
    ::Matrix{Real}

Calculate distance matrix for all combination between `coords1` and `coord2`.
The row of returned matrix coresspond to `coords1`, and the column correspond to `coords2`.
"""
function pair_length_matrix(first_coords::ArrayT1, second_coords::ArrayT2
    )::Matrix{<:Real} where {ArrayT1 <: AbstractArray{<:Coordinate{<:Real}, 1},
                             ArrayT2 <: AbstractArray{<:Coordinate{<:Real}, 1}}

    distance.(first_coords, reshape(second_coords, (1, length(second_coords))))
end

"""
     pair_length_matrix_pbc(corrds1::Vector{Coordinate}, coords2::Vector{Coordinate},
                            box_size::Coordinate)
    ::Matrix{Real}

Calculate distance matrix for all combination between `coords1` and `coord2` considering periodic boundary condition.
The box size information used for calculation is specified by `box_size` argument.
The row of returned matrix coresspond to `coords1`, and the column correspond to `coords2`.
"""
function pair_length_matrix_pbc(first_coords::ArrayT1, second_coords::ArrayT2,
    box_size::Coordinate{<:Real}
    )::Matrix{<:Real} where {ArrayT1 <: AbstractArray{<:Coordinate{<:Real}, 1},
                             ArrayT2 <: AbstractArray{<:Coordinate{<:Real}, 1}}

    distance_pbc.(first_coords, reshape(second_coords, (1, length(second_coords))), box_size)
end

"""
    pair_length_matrix(trj::Trajectory;
                       frame_indices::Union{Vector, OrdinalRange, Colon}       = :,
                       first_atom_indices::Union{Vector, OrdinalRange, Colon}  = :,
                       second_atom_indices::Union{Vector, OrdinalRange, Colon} = atom_indices1)
    ::Vector{Matrix{Coordinate}}

Calculate distance matrix for all combination between `first_atom_indices` and `second_atom_indices`.
This cauculation apply to each frame of trajectory and the result matrices are stored in Vector.
The target frame can be restricted by pass indeces vector or range to `frame_indices`.
"""
function pair_length_matrix(trj::Trajectory{<:Real};
    frame_indices::Union{Vector, OrdinalRange, Colon} = :,
    first_atom_indices::Union{Vector, OrdinalRange, Colon} = :,
    second_atom_indices::Union{Vector, OrdinalRange, Colon} = first_atom_indices
    )::Vector{<:Matrix{<:Real}}

    zip_iterate4frame = zip(eachcol(view(trj.coordinates, first_atom_indices, frame_indices)),
                            eachcol(view(trj.coordinates, second_atom_indices, frame_indices)))
    map(coords_vec_pair -> pair_length_matrix(coords_vec_pair...), zip_iterate4frame)
end

"""
    pair_length_matrix_pbc(trj::Trajectory, box_size::Coordinate;
                           frame_indices::Union{Vector, OrdinalRange, Colon}       = :,
                           first_atom_indices::Union{Vector, OrdinalRange, Colon}  = :,
                           second_atom_indices::Union{Vector, OrdinalRange, Colon} = atom_indices1)
    ::Vector{Matrix{Coordinate}}

Calculate distance matrix for all combinations between `first_atom_indices` and `second_atom_indices` considering periodic boundary condition.
The box size information used for calculation is specified by `box_size` argument.
This cauculation apply to each frame of trajectory and the result matrices are stored in Vector.
The target frame can be restricted by pass indeces vector or range to `frame_indices`.
"""
function pair_length_matrix_pbc(trj::Trajectory{<:Real}, box_size::Coordinate{<:Real};
    frame_indices::Union{Vector, OrdinalRange, Colon} = :,
    first_atom_indices::Union{Vector, OrdinalRange, Colon} = :,
    second_atom_indices::Union{Vector, OrdinalRange, Colon} = first_atom_indices
    )::Vector{<:Matrix{<:Real}}

    zip_iterate4frame = zip(eachcol(view(trj.coordinates, first_atom_indices, frame_indices)),
                            eachcol(view(trj.coordinates, second_atom_indices, frame_indices)))
    map(coords_vec_pair -> pair_length_matrix_pbc(coords_vec_pair..., box_size),
        zip_iterate4frame)
end

"""
    pair_length_matrix_parallel(trj::Trajectory;
                                frame_indices::Union{Vector, OrdinalRange, Colon}       = :,
                                first_atom_indices::Union{Vector, OrdinalRange, Colon}  = :,
                                second_atom_indices::Union{Vector, OrdinalRange, Colon} = atom_indices1)
    ::Vector{Matrix{Coordinate}}

Multi-threads version of pair_length_matrix. If you set available threads number to `Threads.nthreads()`, this function would faster than non-parallel version.
"""
function pair_length_matrix_parallel(trj::Trajectory{RealT};
    frame_indices::Union{Vector, OrdinalRange, Colon} = :,
    first_atom_indices::Union{Vector, OrdinalRange, Colon} = :,
    second_atom_indices::Union{Vector, OrdinalRange, Colon} = first_atom_indices
    )::Vector{Matrix{RealT}} where RealT <: Real

    first_target_coords  = view(trj.coordinates, first_atom_indices,  frame_indices)
    second_target_coords = view(trj.coordinates, second_atom_indices, frame_indices)
    target_frame_num = size(first_target_coords, 2)
    return_vec = Vector{Matrix{Float32}}(undef, target_frame_num)
    Threads.@threads for frame_idx in 1:target_frame_num
        return_vec[frame_idx] =
            pair_length_matrix(view(first_target_coords, :, frame_idx), view(second_target_coords, :, frame_idx))
    end
    return_vec
end

"""
    contact_bool_matrix(threshold::Real, trj::Trajectory;
                        frame_indices::Union{Array, OrdinalRange, Colon}       = :,
                        first_atom_indices::Union{Array, OrdinalRange, Colon}  = :,
                        second_atom_indices::Union{Array, OrdinalRange, Colon} = first_atom_indices)
    ::Vector{Matrix{Bool}}

Judge contact is formed or not. If the distance between two coordinate is shorter than threshold, contact is considered to be formed. In returned vector of matrices, each matrix correspond to contact matrix of each frame.
You can specify the target frames or atoms by `frame_indices`, `first_atom_indices` or `second_atom_indices`. When you specify the target atoms, the row of matrices corresponds to first_atom_indices and column of matrices corresponds to second_atom_indices.
"""
function contact_bool_matrix(threshold::RealT,
    first_coords::AbstractArray{<:Coordinate, 1},
    second_coords::AbstractArray{<:Coordinate, 1}
    )::Matrix{Bool} where RealT <: Real

    length_matrix = pair_length_matrix(first_coords, second_coords)
    map(length -> length < threshold, length_matrix)
end

function contact_bool_matrix(threshold::RealT1, trj::Trajectory{RealT2};
    frame_indices::Union{Vector, OrdinalRange, Colon} = :,
    first_atom_indices::Union{Vector, OrdinalRange, Colon} = :,
    second_atom_indices::Union{Vector, OrdinalRange, Colon} = first_atom_indices
    )::Vector{Matrix{Bool}} where {RealT1 <: Real, RealT2 <: Real}

    length_mat_arr = pair_length_matrix(trj, frame_indices = frame_indices,
                                        first_atom_indices = first_atom_indices,
                                        second_atom_indices = second_atom_indices)
    map(length_mat_arr) do length_matrix
        map(length -> length < threshold, length_matrix)
    end
end

"""
    contact_bool_matrix_pbc(threshold::Real, trj::Trajectory, box_size::Coordinate;
                            frame_indices::Union{Array, OrdinalRange, Colon}       = :,
                            first_atom_indices::Union{Array, OrdinalRange, Colon}  = :,
                            second_atom_indices::Union{Array, OrdinalRange, Colon} = first_atom_indices)
    ::Vector{Matrix{Bool}}

Judge contact is formed or not considering periodic boundary condition. If the distance between two coordinate is shorter than threshold, contact is considered to be formed. In returned vector of matrices, each matrix correspond to contact matrix of each frame.
The box size information used for calculation is specified `box_size` argument.
You can specify the target frames or atoms by `frame_indices`, `first_atom_indices` or `second_atom_indices`. When you specify the target atoms, the row of matrices corresponds to first_atom_indices and column of matrices corresponds to second_atom_indices.
"""
function contact_bool_matrix_pbc(threshold::RealT, trj::Trajectory{<:Real},
    box_size::Coordinate{<:Real};
    frame_indices::Union{Vector, OrdinalRange, Colon} = :,
    first_atom_indices::Union{Vector, OrdinalRange, Colon} = :,
    second_atom_indices::Union{Vector, OrdinalRange, Colon} = first_atom_indices
    )::Vector{Matrix{Bool}} where RealT <: Real

    length_mat_arr = pair_length_matrix_pbc(trj, box_size,
                                            frame_indices = frame_indices,
                                            first_atom_indices = first_atom_indices,
                                            second_atom_indices = second_atom_indices)
    map(length_mat_arr) do length_matrix
        map(length -> length < threshold, length_matrix)
    end
end

"""
    contact_bool_matrix_parallel(threshold::Real, trj::Trajectory;
                                 frame_indices::Union{Vector, OrdinalRange, Colon}       = :,
                                 first_atom_indices::Union{Vector, OrdinalRange, Colon}  = :,
                                 second_atom_indices::Union{Vector, OrdinalRange, Colon} = first_atom_indices)
    ::Vector{Matrix{Bool}}

Multi-threads version of contact_bool_matrix. If you set available threads number to `Threads.nthreads()`, this function would faster than non-parallel version.
"""
function contact_bool_matrix_parallel(threshold::RealT1, trj::Trajectory{RealT2};
    frame_indices::Union{Vector, OrdinalRange, Colon} = :,
    first_atom_indices::Union{Vector, OrdinalRange, Colon} = :,
    second_atom_indices::Union{Vector, OrdinalRange, Colon} = first_atom_indices
    )::Vector{Matrix{Bool}} where {RealT1 <: Real, RealT2 <: Real}

    length_mat_arr = pair_length_matrix_parallel(trj, frame_indices = frame_indices,
                                                 first_atom_indices = first_atom_indices,
                                                 second_atom_indices = second_atom_indices)
    target_frame_num = length(length_mat_arr)
    return_vec = Vector{Matrix{Float32}}(undef, target_frame_num)
    threads_num = Threads.nthreads()
    iteration_in_thread = Int32(floor(target_frame_num / threads_num))
    Threads.@threads for frame_idx in 1:target_frame_num
        return_vec[frame_idx] = map(length -> length < threshold, length_mat_arr[frame_idx])
    end
    return_vec
end

"""
    contact_count_matrix(threshold::Real, trj::Trajectory;
                         frame_indices::Union{Vector, OrdinalRange, Colon}       = :,
                         first_atom_indices::Union{Vector, OrdinalRange, Colon}  = :,
                         second_atom_indices::Union{Vector, OrdinalRange, Colon} = first_atom_indices)
    ::Matrix{Real}

Calculate contact formation count over the trajectory. If the distance between two coordinate is shorter than threshold, contact is considered to be formed.
You can specify the target frames or atoms by `frame_indices`, `first_atom_indices` or `second_atom_indices`. When you specify the target atoms, the row of matrices corresponds to first_atom_indices and column of matrices corresponds to second_atom_indices.
"""
function contact_count_matrix(threshold::RealT1, trj::Trajectory{RealT2};
    frame_indices::Union{Vector, OrdinalRange, Colon} = :,
    first_atom_indices::Union{Vector, OrdinalRange, Colon} = :,
    second_atom_indices::Union{Vector, OrdinalRange, Colon} = first_atom_indices
    )::Matrix{RealT2} where {RealT1 <: Real, RealT2 <: Real}

    first_target_coords  = view(trj.coordinates, first_atom_indices,  frame_indices)
    second_target_coords = view(trj.coordinates, second_atom_indices, frame_indices)
    target_frame_num = size(first_target_coords, 2)
    count_mat = zeros(RealT2, (size(first_target_coords, 1), size(second_target_coords, 1)))
    for frame_idx in 1:target_frame_num
        count_mat +=
            contact_bool_matrix(threshold, view(first_target_coords, :, frame_idx), view(second_target_coords, :, frame_idx))
    end
    count_mat
end

"""
    contact_count_matrix_parallel(threshold::Real, trj::Trajectory;
                                  frame_indices::Union{Vector, OrdinalRange, Colon}       = :,
                                  first_atom_indices::Union{Vector, OrdinalRange, Colon}  = :,
                                  second_atom_indices::Union{Vector, OrdinalRange, Colon} = first_atom_indices)
    ::Matrix{Real}

Multi-threads version of contact_count_matrix. If you set available threads number to `Threads.nthreads()`, this function would faster than non-parallel version.
"""
function contact_count_matrix_parallel(threshold::RealT1, trj::Trajectory{RealT2};
    frame_indices::Union{Vector, OrdinalRange, Colon} = :,
    first_atom_indices::Union{Vector, OrdinalRange, Colon} = :,
    second_atom_indices::Union{Vector, OrdinalRange, Colon} = first_atom_indices
    )::Matrix{RealT2} where {RealT1 <: Real, RealT2 <: Real}

    first_target_coords  = view(trj.coordinates, first_atom_indices,  frame_indices)
    second_target_coords = view(trj.coordinates, second_atom_indices, frame_indices)
    target_frame_num = size(first_target_coords, 2)
    first_atoms_len = size(first_target_coords, 1)
    second_atoms_len = size(second_target_coords, 1)
    count_mat_arr =
        [zeros(RealT2, first_atoms_len, second_atoms_len) for thread in 1:Threads.nthreads()]
    Threads.@threads for frame_idx in 1:target_frame_num
        count_mat_arr[Threads.threadid()] +=
            contact_bool_matrix(threshold,
                                view(first_target_coords, :, frame_idx), view(second_target_coords, :, frame_idx))
    end
    sum(count_mat_arr)
end

"""
    contact_probability_matrix(threshold::Real, trj::Trajectory;
                               frame_indices::Union{Vector, OrdinalRange, Colon}       = :,
                               first_atom_indices::Union{Vector, OrdinalRange, Colon}  = :,
                               second_atom_indices::Union{Vector, OrdinalRange, Colon} = first_atom_indices)
    ::Matrix{Real}

Calculate contact formation probability over the trajectory. If the distance between two coordinate is shorter than threshold, contact is considered to be formed.
You can specify the target frames or atoms by `frame_indices`, `first_atom_indices` or `second_atom_indices`. When you specify the target atoms, the row of matrices corresponds to first_atom_indices and column of matrices corresponds to second_atom_indices.
"""
function contact_probability_matrix(threshold::RealT1, trj::Trajectory{RealT2};
    frame_indices::Union{Vector, OrdinalRange, Colon} = :,
    first_atom_indices::Union{Vector, OrdinalRange, Colon} = :,
    second_atom_indices::Union{Vector, OrdinalRange, Colon} = first_atom_indices
    )::Matrix{RealT2} where {RealT1 <: Real, RealT2 <: Real}

    count_mat= contact_count_matrix(threshold, trj,
                                    frame_indices       = frame_indices,
                                    first_atom_indices  = first_atom_indices,
                                    second_atom_indices = second_atom_indices)
    count_mat / size(view(trj.coordinates, :, frame_indices), 2)
end

"""
    contact_probability_matrix_parallel(threshold::Real, trj::Trajectory;
                                        frame_indices::Union{Vector, OrdinalRange, Colon}       = :,
                                        first_atom_indices::Union{Vector, OrdinalRange, Colon}  = :,
                                        second_atom_indices::Union{Vector, OrdinalRange, Colon} = first_atom_indices)
    ::Matrix{Real}

Multi-threads version of contact_probability_matrix. If you set available threads number to `Threads.nthreads()`, this function would faster than non-parallel version.
"""
function contact_probability_matrix_parallel(threshold::RealT1, trj::Trajectory{RealT2};
    frame_indices::Union{Vector, OrdinalRange, Colon} = :,
    first_atom_indices::Union{Vector, OrdinalRange, Colon} = :,
    second_atom_indices::Union{Vector, OrdinalRange, Colon} = first_atom_indices
    )::Matrix{RealT2} where {RealT1 <: Real, RealT2 <: Real}

    count_mat_arr = contact_count_matrix_parallel(threshold, trj,
                                                  frame_indices       = frame_indices,
                                                  first_atom_indices  = first_atom_indices,
                                                  second_atom_indices = second_atom_indices)
    count_mat_arr / size(view(trj.coordinates, :,frame_indices), 2)
end

"""
    radius_of_gyration(trj::Trajectory;
                       frame_indices::Union{Vector, OrdinalRange, Colon} = :,
                       atom_indices ::Union{Vector, OrdinalRange, Colon} = :,
                       geometric::Bool = false)
    ::Vector{Real}

Calculate radius of gyration over the trajectory.
If you set `geometric` is `true`, this function calculate radius of gyration without mass info.
"""
function radius_of_gyration(trj::Trajectory{RealT};
             frame_indices::Union{Vector, OrdinalRange, Colon} = :,
             atom_indices ::Union{Vector, OrdinalRange, Colon} = :,
             geometric::Bool = false
             )::Vector{RealT} where RealT <: Real

    com = center_of_mass(trj,
                         frame_indices = frame_indices,
                         atom_indices  = atom_indices,
                         geometric     = geometric)

    target_coordinates = view(trj.coordinates, atom_indices, frame_indices)
    target_frame_num = length(com)
    return_vec = Vector{RealT}(undef, target_frame_num)

    dist_from_com_vec = norm.(target_coordinates .- permutedims(com))

    if !geometric
        mass_vec = map(attributes -> attributes.mass, view(trj.attributes, atom_indices))
        mass_sum = sum(mass_vec)
        for frame_idx in 1:target_frame_num
            return_vec[frame_idx] = sqrt(sum(view(dist_from_com_vec, :, frame_idx).^2 .* mass_vec) / mass_sum)
        end
    else target_atom_num = size(target_coordinates)[1]
        for frame_idx in 1:target_frame_num
            return_vec[frame_idx] = sqrt(sum(view(dist_from_com_vec, :, frame_idx).^2) / target_atom_num)
        end
    end
    return_vec
end

"""
    atom_mass(name::AbstractString)
    ::Float32

Return mass of the `name` atom.
The kind of each atom is judged by PDB format atom names rule.
https://cdn.rcsb.org/wwpdb/docs/documentation/file-format/PDB_format_1992.pdf
"""
function atom_mass(name::AbstractString)::Float32
    if occursin(r"^(H|[1-4]H)", name)
        ATOMNAME2MASS["H"]
    elseif occursin(r"^C($|[^L]+)", name)
        ATOMNAME2MASS["C"]
    elseif occursin(r"^O", name)
        ATOMNAME2MASS["O"]
    elseif occursin(r"^N($|[^A]+)", name)
        ATOMNAME2MASS["N"]
    elseif occursin(r"^S", name)
        ATOMNAME2MASS["S"]
    elseif occursin(r"^P", name)
        ATOMNAME2MASS["P"]
    else
        ATOMNAME2MASS[name]
    end
end

"""
    residue_mass(name::AbstractString)
    ::Float32

Return mass of the `name` residue.
The kind of each residue is judged by amino acid 3-letter abbreviation.
"""
function residue_mass(name::AbstractString)::Float32
    RESNAME2MASS[name]
end

"""
    distance_pdc(first_atom::Coordinate,  second_atom::Coordinate,
                 box_size::Coordinate)
    ::Real

Calculate distance `first_atom` and `second_atom` considering periodic boundary condition.
The box size information used calculatin is specified by `box_size` argument.
"""
function distance_pbc(first_coord::Coordinate{RealT},  second_coord::Coordinate{RealT},
    box_size::Coordinate{<:Real}
    )::RealT where RealT <: Real

    half_box_size = box_size * 0.5
    dist_vec = first_coord  - second_coord
    x = abs(dist_vec.x) < half_box_size.x ? dist_vec.x : box_size.x - abs(dist_vec.x)
    y = abs(dist_vec.y) < half_box_size.y ? dist_vec.y : box_size.y - abs(dist_vec.y)
    z = abs(dist_vec.z) < half_box_size.z ? dist_vec.z : box_size.z - abs(dist_vec.z)
    sqrt(x^2 + y^2 + z^2)
end

# This function get half_box_size, because if process will apply to all trajectory,
# half box size calculation will occur every time frame, this is overhead.
function fix_pbc(coordinates::Vector{<:Coordinate{RealT}}, groupid_vec::Vector{<:Integer},
    box_size::Coordinate{<:Real}, half_box_size::Coordinate{<:Real};
    x::Bool= true, y::Bool = true, z::Bool = true
    )::Vector{<:Coordinate{RealT}} where RealT <: Real

    if length(coordinates) != length(groupid_vec)
        throw(AssertionError("""
                             coordinate vector and groupid vector should be same length.
                             """))
    end

    new_coords = deepcopy(coordinates)
    unique_groupid_vec = unique(groupid_vec)

    for groupid in unique_groupid_vec
        same_group_ids = findall(id->id==groupid, groupid_vec)
        @views sbj_coords = new_coords[same_group_ids[2:end]]
        @views dist2first = sbj_coords .- coordinates[same_group_ids[1]]
        sbj_dist_zip      = zip(sbj_coords, dist2first)

        if x
            for (coord, dist) in sbj_dist_zip
                coord.x = abs(dist.x) < half_box_size.x ? coord.x : coord.x - sign(dist.x) * box_size.x
            end
        end

        if y
            for (coord, dist) in sbj_dist_zip
                coord.y = abs(dist.y) < half_box_size.y ? coord.y : coord.y - sign(dist.y) * box_size.y
            end
        end

        if z
            for (coord, dist) in sbj_dist_zip
                coord.z = abs(dist.z) < half_box_size.z ? coord.z : coord.z - sign(dist.z) * box_size.z
            end
        end
    end
    new_coords
end

"""
    fix_pbc(coordinates::Vector{Coordinate}, groupid_vec::Vector{Integer},
    box_size::Coordinate{<:Real};
    x::Bool = true, y::Bool = true, z::Bool = true)::Vector{Coordinate}

Fix the atom group splited by periodic boundary box. This group is defined by `groupid_vec` which contain each particles group id.
For example, if your system have 3 atom, and atom 1 and 2 are group 1, and atom 3 is group 2, this `groupid_vec` is [1, 1, 2].
If the group over the boundary of the box, the atoms in the group which separated from first atom of that move to the position where the atom locate without periodic boundary condition.
You can specify the axies which applied fixing by x,y and option. If you set `z = false`, only x and y coordinate fixed and z don't change.
"""
function fix_pbc(coordinates::Vector{<:Coordinate{RealT}}, groupid_vec::Vector{<:Integer},
    box_size::Coordinate{<:Real};
    x::Bool = true, y::Bool = true, z::Bool = true
    )::Vector{<:Coordinate{RealT}} where RealT <: Real
    half_box_size = box_size * 0.5
    fix_pbc(coordinates, groupid_vec, box_size, half_box_size, x=x, y=y, z=z)
end

"""
    fix_pbc(trj::Trajectory, groupid_vec::Vector{Integer},
    box_size::Coordinate;
    x::Bool = true, y::Bool = true, z::Bool = true)::Trajectory

Fix residues splited by periodic boundary condition. This is more specific version of `fix_pbc` function for trajectory handling.
For example, if your system have 3 atom, and atom 1 and 2 are group 1, and atom 3 is group 2, this `groupid_vec` is [1, 1, 2].
If the group over the boundary of the box, the atoms in the group which separated from first atom of that move to the position where the atom locate without periodic boundary condition.
You can specify the axies which applied fixing by x,y and option. If you set `z = false`, only x and y coordinate fixed and z don't change.
"""
function fix_pbc(trj::Trajectory{RealT}, groupid_vec::Vector{<:Integer},
    box_size::Coordinate{<:Real};
    x::Bool = true, y::Bool = true, z::Bool = true
    )::Trajectory{RealT} where RealT <: Real

    new_trj = deepcopy(trj)
    half_box_size = box_size * 0.5

    unique_groupid_vec = unique(groupid_vec)
    coordinates        = new_trj.coordinates
    for groupid in unique_groupid_vec
        same_mol_indices = findall(id->id==groupid, groupid_vec)
        @views sbj_coords = coordinates[same_mol_indices[2:end], :]
        @views dist2first_mat = sbj_coords .- permutedims(coordinates[same_mol_indices[1], :])
        if x
            for (coord, dist) in zip(sbj_coords, dist2first_mat)
                coord.x = abs(dist.x) < half_box_size.x ? coord.x : coord.x - sign(dist.x) * box_size.x
            end
        end

        if y
            for (coord, dist) in zip(sbj_coords, dist2first_mat)
                coord.y = abs(dist.y) < half_box_size.y ? coord.y : coord.y - sign(dist.y) * box_size.y
            end
        end

        if z
            for (coord, dist) in zip(sbj_coords, dist2first_mat)
                coord.z = abs(dist.z) < half_box_size.z ? coord.z : coord.z - sign(dist.z) * box_size.z
            end
        end
    end
    new_trj
end

function fix_pbc!(trj::Trajectory{RealT}, box_size::Coordinate{<:Real};
    x::Bool = true, y::Bool = true, z::Bool = true
    ) where RealT <: Real

    # fix atom in over pbc residue based on first atom of the residue
    half_box_size = box_size * 0.5

    resid_vec  = map(attr->attr.resid, trj.attributes)
    unique_resid_vec = unique(resid_vec)
    coordinates = trj.coordinates
    attributes  = trj.attributes
    for resid in unique_resid_vec
        same_mol_indices = findall(id->id==resid, resid_vec)
        @views sbj_coords = coordinates[same_mol_indices[2:end], :]
        @views dist2first_mat = sbj_coords .- permutedims(coordinates[same_mol_indices[1], :])
        if x
            for (coord, dist) in zip(sbj_coords, dist2first_mat)
                coord.x = abs(dist.x) < half_box_size.x ? coord.x : coord.x - sign(dist.x) * box_size.x
            end
        end
        if y
            for (coord, dist) in zip(sbj_coords, dist2first_mat)
                coord.y = abs(dist.y) < half_box_size.y ? coord.y : coord.y - sign(dist.y) * box_size.y
            end
        end
        if z
            for (coord, dist) in zip(sbj_coords, dist2first_mat)
                coord.z = abs(dist.z) < half_box_size.z ? coord.z : coord.z - sign(dist.z) * box_size.z
            end
        end
    end
end

"""
    move_pbc_center(coordinates::Vector{Coordinate},
                    new_center::Coordinate, box_size::Coordinate)
                    )::Vector{Coordinate}
Fix coordinate to consistent with new box center and box size.
Box size should match to original periodic boundary box size.
"""
function move_pbc_center(coordinates::Vector{<:Coordinate{RealT}},
    center::Coordinate{<:Real}, box_size::Coordinate{<:Real}
    )::Vector{Coordinate{RealT}} where RealT <: Real

    new_coords = deepcopy(coordinates)
    half_box = box_size * 0.5
    for coord in new_coords
        dist = coord - center
        coord.x = abs(dist.x) < half_box.x ? coord.x : coord.x - sign(dist.x) * box_size.x
        coord.y = abs(dist.y) < half_box.y ? coord.y : coord.y - sign(dist.y) * box_size.y
        coord.z = abs(dist.z) < half_box.z ? coord.z : coord.z - sign(dist.z) * box_size.z
    end
    new_coords
end


text:
nc_date(ds::NCdata) = ds["time"][:]

function nc_date(file::NCfiles)
    nc_open(file) do ds
        nc_date(ds)
    end
end

function nc_calendar(file::AbstractString)
    nc_open(file) do ds
        # typeof(ds["time"][1])
        ds["time"].attrib["calendar"]
    end
end


text:
module GammaTest

using Test
using ForneyLab
using ForneyLab: prod!, unsafeMean, unsafeVar, outboundType, isApplicable, dims, naturalParams, standardDistribution
using ForneyLab: SPGammaOutNPP, VBGammaOut, VBGammaA, VBGammaB

@testset "dims" begin
    @test dims(ProbabilityDistribution(Univariate, Gamma, a=1.0, b=1.0)) == ()
end

@testset "vague" begin
    @test vague(Gamma) == ProbabilityDistribution(Univariate, Gamma, a=1.0, b=tiny)
end

@testset "prod!" begin
    @test ProbabilityDistribution(Univariate, Gamma, a=1.0, b=2.0) * ProbabilityDistribution(Univariate, Gamma, a=3.0, b=4.0) == ProbabilityDistribution(Univariate, Gamma, a=3.0, b=6.0)
    @test ProbabilityDistribution(Univariate, Gamma, a=1.0, b=2.0) * ProbabilityDistribution(Univariate, PointMass, m=1.0) == ProbabilityDistribution(Univariate, PointMass, m=1.0)
    @test ProbabilityDistribution(Univariate, PointMass, m=1.0) * ProbabilityDistribution(Univariate, Gamma, a=1.0, b=2.0) == ProbabilityDistribution(Univariate, PointMass, m=1.0)
    @test_throws Exception ProbabilityDistribution(Univariate, PointMass, m=-1.0) * ProbabilityDistribution(Univariate, Gamma, a=1.0, b=2.0)
end

@testset "unsafe mean and variance" begin
    @test unsafeMean(ProbabilityDistribution(Univariate, Gamma, a=1.0, b=2.0)) == 0.5
    @test unsafeVar(ProbabilityDistribution(Univariate, Gamma, a=1.0, b=2.0)) == 0.25
end

@testset "natural parameters" begin
    d = ProbabilityDistribution(Univariate, Gamma, a=2.0, b=5.0)
    η = naturalParams(d)
    s = standardDistribution(Univariate, Gamma, η=η)
    @test d.params[:a] == s.params[:a] # Test conversion consistency
    @test d.params[:b] == s.params[:b]

    x = [0.1, 2.0, 15.0]
    d_x = logPdf.([d], x)
    η_x = logPdf.(Univariate, Gamma, x; η=η)
    @test isapprox(d_x, η_x) # Test pdf consistency
end

#-------------
# Update rules
#-------------

@testset "SPGammaOutNPP" begin
    @test SPGammaOutNPP <: SumProductRule{Gamma}
    @test outboundType(SPGammaOutNPP) == Message{Gamma}
    @test isApplicable(SPGammaOutNPP, [Nothing, Message{PointMass}, Message{PointMass}])

    @test ruleSPGammaOutNPP(nothing, Message(Univariate, PointMass, m=1.0), Message(Univariate, PointMass, m=2.0)) == Message(Univariate, Gamma, a=1.0, b=2.0)
end

@testset "VBGammaOut" begin
    @test VBGammaOut <: NaiveVariationalRule{Gamma}
    @test outboundType(VBGammaOut) == Message{Gamma}
    @test isApplicable(VBGammaOut, [Nothing, ProbabilityDistribution, ProbabilityDistribution])
    @test !isApplicable(VBGammaOut, [ProbabilityDistribution, ProbabilityDistribution, Nothing])

    @test ruleVBGammaOut(nothing, ProbabilityDistribution(Univariate, PointMass, m=1.5), ProbabilityDistribution(Univariate, PointMass, m=3.0)) == Message(Univariate, Gamma, a=1.5, b=3.0)
end

@testset "VBGammaA" begin
    @test VBGammaA <: NaiveVariationalRule{Gamma}
    @test outboundType(VBGammaA) == Message{Function}
    @test !isApplicable(VBGammaA, [Nothing, ProbabilityDistribution, ProbabilityDistribution])
    @test isApplicable(VBGammaA, [ProbabilityDistribution, Nothing, ProbabilityDistribution])
end

@testset "VBGammaB" begin
    @test VBGammaB <: NaiveVariationalRule{Gamma}
    @test outboundType(VBGammaB) == Message{Gamma}
    @test !isApplicable(VBGammaB, [Nothing, ProbabilityDistribution, ProbabilityDistribution])
    @test isApplicable(VBGammaB, [ProbabilityDistribution, ProbabilityDistribution, Nothing])

    @test ruleVBGammaB(ProbabilityDistribution(Univariate, PointMass, m=1.5), ProbabilityDistribution(Univariate, PointMass, m=3.0), nothing) == Message(Univariate, Gamma, a=4.0, b=1.5)
end

@testset "averageEnergy and differentialEntropy" begin
    @test differentialEntropy(ProbabilityDistribution(Univariate, Gamma, a=1.0, b=2.0)) == averageEnergy(Gamma, ProbabilityDistribution(Univariate, Gamma, a=1.0, b=2.0), ProbabilityDistribution(Univariate, PointMass, m=1.0), ProbabilityDistribution(Univariate, PointMass, m=2.0))
end

end #module


text:
# This code is based on https://github.com/JuliaLang/Pkg.jl/blob/53b3580707309d6170651cfc79ddb17b8d4ff8f9/src/Types.jl

using Pkg
using UUIDs

const STDLIB = Ref{Dict{UUID, String}}()
const STDLIB_LOCK = ReentrantLock()

function projectfile_path(env_path::String; strict=false)
    for name in Base.project_names
        maybe_file = joinpath(env_path, name)
        isfile(maybe_file) && return maybe_file
    end
    return strict ? nothing : joinpath(env_path, "Project.toml")
end

stdlib_dir() = normpath(joinpath(Sys.BINDIR, "..", "share", "julia", "stdlib", "v$(VERSION.major).$(VERSION.minor)"))

stdlib_path(stdlib::String) = joinpath(stdlib_dir(), stdlib)

function load_stdlib()
    stdlib = Dict{UUID,String}()
    for name in readdir(stdlib_dir())
        projfile = projectfile_path(stdlib_path(name); strict=true)
        nothing === projfile && continue
        project = Pkg.TOML.parsefile(projfile)
        uuid = get(project, "uuid", nothing)
        nothing === uuid && continue
        stdlib[UUID(uuid)] = name
    end
    return stdlib
end

function stdlib()
    lock(STDLIB_LOCK) do
        if !isassigned(STDLIB)
            STDLIB[] = load_stdlib()
        end
        stdlib = deepcopy(STDLIB[])
    end
end

function stdlibmodules(m::Module)::Vector{Module}
    result = Vector{Module}(undef, 0)
    for x in values(stdlib())
        s = Symbol(x)
        if isdefined(m, s)
            m_dot_s = getproperty(m, s)
            if m_dot_s isa Module
                push!(result, m_dot_s)
            end
        end
    end
    return result
end


text:
using Plots
pyplot();
using Images

using FileIO: load, save, loadstreaming, savestreaming
import LibSndFile
using SampledSignals
#SampledSignals.embed_javascript();

using DSP
using FFTW
using Statistics

using Distributed, SharedArrays

using Printf

function make_animation(filename, framesize, zero_pad_to, FPS)
  snd = load(filename);
  fs = snd.samplerate;
  x = vec(snd.data);
  xx = zeros((length(x) + framesize));
  xx[1:length(x)] = x;
  w = hanning(framesize);

  framestep = convert(Int, fs/FPS);
  zz = SharedArray{Float64}(floor(Int, length(xx)/framestep), round(Int, max(nextpow(2, framesize), zero_pad_to)/2+1));
  @sync @distributed for n = 1:floor(Int, length(x)/framestep)
    xxx = xx[(n-1)*framestep+1 : (n-1)*framestep+framesize];
    z = zeros(max(nextpow(2, framesize), zero_pad_to));
    z[1:framesize] = xxx .* w;
    y = FFTW.fft(z) / length(z);
    y[2:round(Int, length(y)/2)] *= 2.0;
    y = y[1:round(Int, length(y)/2+1)];
    zz[n, :] = abs.(y);
  end
  zmax = 20*log10.(maximum(zz));
  t = range(0, stop=length(xx)-1, length=size(zz, 1)) / fs;
  f = range(0, stop=fs/2, length=size(zz, 2));

  anim = @animate for n=1:size(zz, 1)
    zframe = 20*log10.(zz[n, :]);
    plot(f, zframe .- zmax,
         xlim=(0, 1500), ylim=(-75, 3),
         xlabel="Frequency (Hz)", ylabel="Power (dB)",
         #title=sprintf("%s t=%.2fsec", filename, round(t[n]*100)/100),
         legend=false)
  end
  gif(anim, "$(filename)_$(round(Int, fs/framestep))fps.gif", fps=round(Int, fs/framestep));
end

filename = "electric_bass.flac"
framesize = 2^14;
zero_pad_to = 2^18;
FPS = 15;
make_animation(filename, framesize, zero_pad_to, FPS);


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate(x^-1*arccsc(a*x^5), x) == :(-1//5*arccsc(a*x^5)*log(1+-1*exp(2*im*arccsc(a*x^5)))+1//10*im*arccsc(a*x^5)^2+1//10*im*Polylogarithms.polylog(2,exp(2*im*arccsc(a*x^5))))
@test integrate(x^3*arccsc(x^(1/2)), x) == :(1//4*(-1+x)^(1/2)+1//4*(-1+x)^3//2+1//28*(-1+x)^7//2+3//20*(-1+x)^5//2+1//4*x^4*arccsc(x^(1/2)))
@test integrate(x^2*arccsc(x^(1/2)), x) == :(1//3*(-1+x)^(1/2)+1//15*(-1+x)^5//2+2//9*(-1+x)^3//2+1//3*x^3*arccsc(x^(1/2)))
@test integrate(x*arccsc(x^(1/2)), x) == :((1/2)*(-1+x)^(1/2)+1//6*(-1+x)^3//2+(1/2)*x^2*arccsc(x^(1/2)))
@test integrate(arccsc(x^(1/2)), x) == :((-1+x)^(1/2)+x*arccsc(x^(1/2)))
@test integrate(x^-1*arccsc(x^(1/2)), x) == :(im*arccsc(x^(1/2))^2+im*Polylogarithms.polylog(2,exp(2*im*arccsc(x^(1/2))))+-2*arccsc(x^(1/2))*log(1+-1*exp(2*im*arccsc(x^(1/2)))))
@test integrate(x^-2*arccsc(x^(1/2)), x) == :(-1//2*arctan((-1+x)^(1/2))+-1*x^-1*arccsc(x^(1/2))+-1//2*x^-1*(-1+x)^(1/2))
@test integrate(x^-3*arccsc(x^(1/2)), x) == :(-3//16*arctan((-1+x)^(1/2))+-3//16*x^-1*(-1+x)^(1/2)+-1//2*x^-2*arccsc(x^(1/2))+-1//8*x^-2*(-1+x)^(1/2))
@test integrate(x^2*arccsc(a*x^-1), x) == :(-1//9*a^3*(1+-1*a^-2*x^2)^3//2+1//3*a^3*(1+-1*a^-2*x^2)^(1/2)+1//3*x^3*arcsin(x*a^-1))
@test integrate(x*arccsc(a*x^-1), x) == :((1/2)*x^2*arcsin(x*a^-1)+-1//4*a^2*arcsin(x*a^-1)+1//4*a*x*(1+-1*a^-2*x^2)^(1/2))
@test integrate(arccsc(a*x^-1), x) == :(a*(1+-1*a^-2*x^2)^(1/2)+x*arcsin(x*a^-1))
@test integrate(x^-1*arccsc(a*x^-1), x) == :(arcsin(x*a^-1)*log(1+-1*exp(2*im*arcsin(x*a^-1)))+-1//2*im*arcsin(x*a^-1)^2+-1//2*im*Polylogarithms.polylog(2,exp(2*im*arcsin(x*a^-1))))
@test integrate(x^-2*arccsc(a*x^-1), x) == :(-1*a^-1*arctanh((1+-1*a^-2*x^2)^(1/2))+-1*x^-1*arcsin(x*a^-1))
@test integrate(x^-3*arccsc(a*x^-1), x) == :(-1//2*x^-2*arcsin(x*a^-1)+-1//2*a^-1*x^-1*(1+-1*a^-2*x^2)^(1/2))
@test integrate(x^-4*arccsc(a*x^-1), x) == :(-1//3*x^-3*arcsin(x*a^-1)+-1//6*a^-3*arctanh((1+-1*a^-2*x^2)^(1/2))+-1//6*a^-1*x^-2*(1+-1*a^-2*x^2)^(1/2))
@test integrate(x^-1*arccsc(a*x^n), x) == :((1/2)*im*n^-1*arccsc(a*x^n)^2+(1/2)*im*n^-1*Polylogarithms.polylog(2,exp(2*im*arccsc(a*x^n)))+-1*n^-1*arccsc(a*x^n)*log(1+-1*exp(2*im*arccsc(a*x^n))))
@test integrate(x^4*arccsc(a+b*x), x) == :(1//5*x^5*arccsc(a+b*x)+1//5*a^5*b^-5*arccsc(a+b*x)+1//40*b^-5*(3+40*a^2+40*a^4)*arctanh((1+-1*(a+b*x)^-2)^(1/2))+1//20*b^-2*x^3*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)+1//120*b^-5*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)^2*(9+58*a^2)+-11//60*a*b^-3*x^2*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)+-1//30*a*b^-5*(1+-1*(a+b*x)^-2)^(1/2)*(20+53*a^2)*(a+b*x))
@test integrate(x^3*arccsc(a+b*x), x) == :(1//4*x^4*arccsc(a+b*x)+-1//4*a^4*b^-4*arccsc(a+b*x)+-1//2*a*b^-4*(1+2*a^2)*arctanh((1+-1*(a+b*x)^-2)^(1/2))+-1//3*a*b^-4*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)^2+1//12*b^-4*(1+-1*(a+b*x)^-2)^(1/2)*(2+17*a^2)*(a+b*x)+1//12*b^-2*x^2*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x))
@test integrate(x^2*arccsc(a+b*x), x) == :(1//3*x^3*arccsc(a+b*x)+1//3*a^3*b^-3*arccsc(a+b*x)+1//6*b^-3*(1+6*a^2)*arctanh((1+-1*(a+b*x)^-2)^(1/2))+-5//6*a*b^-3*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)+1//6*x*b^-2*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x))
@test integrate(x*arccsc(a+b*x), x) == :((1/2)*x^2*arccsc(a+b*x)+(1/2)*b^-2*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)+-1*a*b^-2*arctanh((1+-1*(a+b*x)^-2)^(1/2))+-1//2*a^2*b^-2*arccsc(a+b*x))
@test integrate(arccsc(a+b*x), x) == :(b^-1*arctanh((1+-1*(a+b*x)^-2)^(1/2))+b^-1*(a+b*x)*arccsc(a+b*x))
@test integrate(x^-1*arccsc(a+b*x), x) == :(arccsc(a+b*x)*log(1+im*a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x)))+arccsc(a+b*x)*log(1+im*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x)))+(1/2)*im*Polylogarithms.polylog(2,exp(2*im*arccsc(a+b*x)))+-1*im*Polylogarithms.polylog(2,-1*im*a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x)))+-1*im*Polylogarithms.polylog(2,-1*im*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x)))+-1*arccsc(a+b*x)*log(1+-1*exp(2*im*arccsc(a+b*x))))
@test integrate(x^-2*arccsc(a+b*x), x) == :(-1*x^-1*arccsc(a+b*x)+-1*b*a^-1*arccsc(a+b*x)+-2*b*a^-1*(1+-1*a^2)^-1//2*arctan((1+-1*a^2)^-1//2*(a+-1*tan((1/2)*arccsc(a+b*x)))))
@test integrate(x^-3*arccsc(a+b*x), x) == :(-1//2*x^-2*arccsc(a+b*x)+(1/2)*a^-2*b^2*arccsc(a+b*x)+a^-2*b^2*(1+-1*a^2)^-3//2*(1+-2*a^2)*arctan((1+-1*a^2)^-1//2*(a+-1*tan((1/2)*arccsc(a+b*x))))+-1//2*b*a^-1*x^-1*(1+-1*a^2)^-1*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x))
@test integrate(x^-4*arccsc(a+b*x), x) == :(-1//3*x^-3*arccsc(a+b*x)+-1//3*a^-3*b^3*arccsc(a+b*x)+-1//3*a^-3*b^3*(1+-1*a^2)^-5//2*(2+-5*a^2+6*a^4)*arctan((1+-1*a^2)^-1//2*(a+-1*tan((1/2)*arccsc(a+b*x))))+-1//6*b*a^-1*x^-2*(1+-1*a^2)^-1*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)+1//6*a^-2*b^2*x^-1*(1+-1*a^2)^-2*(1+-1*(a+b*x)^-2)^(1/2)*(2+-5*a^2)*(a+b*x))
@test integrate(x^-5*arccsc(a+b*x), x) == :(-1//4*x^-4*arccsc(a+b*x)+1//4*a^-4*b^4*arccsc(a+b*x)+1//4*a^-4*b^4*(1+-1*a^2)^-7//2*(2+-8*a^6+-7*a^2+8*a^4)*arctan((1+-1*a^2)^-1//2*(a+-1*tan((1/2)*arccsc(a+b*x))))+-1//12*b*a^-1*x^-3*(1+-1*a^2)^-1*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)+-1//24*a^-3*b^3*x^-1*(1+-1*a^2)^-3*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)*(6+-17*a^2+26*a^4)+1//24*a^-2*b^2*x^-2*(1+-1*a^2)^-2*(1+-1*(a+b*x)^-2)^(1/2)*(3+-8*a^2)*(a+b*x))
@test integrate(x^3*arccsc(a+b*x)^2, x) == :(1//3*b^-4*log(a+b*x)+1//4*x^4*arccsc(a+b*x)^2+1//12*b^-4*(a+b*x)^2+-1*a*x*b^-3+3*a^2*b^-4*log(a+b*x)+-1//4*a^4*b^-4*arccsc(a+b*x)^2+im*a*b^-4*Polylogarithms.polylog(2,-1*exp(im*arccsc(a+b*x)))+-1*im*a*b^-4*Polylogarithms.polylog(2,exp(im*arccsc(a+b*x)))+-4*a^3*b^-4*arccsc(a+b*x)*arctanh(exp(im*arccsc(a+b*x)))+-2*im*a^3*b^-4*Polylogarithms.polylog(2,exp(im*arccsc(a+b*x)))+-2*a*b^-4*arccsc(a+b*x)*arctanh(exp(im*arccsc(a+b*x)))+2*im*a^3*b^-4*Polylogarithms.polylog(2,-1*exp(im*arccsc(a+b*x)))+1//3*b^-4*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)*arccsc(a+b*x)+1//6*b^-4*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)^3*arccsc(a+b*x)+-1*a*b^-4*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)^2*arccsc(a+b*x)+3*a^2*b^-4*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)*arccsc(a+b*x))
@test integrate(x^2*arccsc(a+b*x)^2, x) == :(1//3*x*b^-2+1//3*x^3*arccsc(a+b*x)^2+-2*a*b^-3*log(a+b*x)+-1//3*im*b^-3*Polylogarithms.polylog(2,-1*exp(im*arccsc(a+b*x)))+1//3*im*b^-3*Polylogarithms.polylog(2,exp(im*arccsc(a+b*x)))+1//3*a^3*b^-3*arccsc(a+b*x)^2+2//3*b^-3*arccsc(a+b*x)*arctanh(exp(im*arccsc(a+b*x)))+-2*im*a^2*b^-3*Polylogarithms.polylog(2,-1*exp(im*arccsc(a+b*x)))+2*im*a^2*b^-3*Polylogarithms.polylog(2,exp(im*arccsc(a+b*x)))+4*a^2*b^-3*arccsc(a+b*x)*arctanh(exp(im*arccsc(a+b*x)))+1//3*b^-3*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)^2*arccsc(a+b*x)+-2*a*b^-3*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)*arccsc(a+b*x))
@test integrate(x*arccsc(a+b*x)^2, x) == :(b^-2*log(a+b*x)+(1/2)*x^2*arccsc(a+b*x)^2+-1//2*a^2*b^-2*arccsc(a+b*x)^2+b^-2*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)*arccsc(a+b*x)+-4*a*b^-2*arccsc(a+b*x)*arctanh(exp(im*arccsc(a+b*x)))+-2*im*a*b^-2*Polylogarithms.polylog(2,exp(im*arccsc(a+b*x)))+2*im*a*b^-2*Polylogarithms.polylog(2,-1*exp(im*arccsc(a+b*x))))
@test integrate(arccsc(a+b*x)^2, x) == :(b^-1*arccsc(a+b*x)^2*(a+b*x)+-2*im*b^-1*Polylogarithms.polylog(2,-1*exp(im*arccsc(a+b*x)))+2*im*b^-1*Polylogarithms.polylog(2,exp(im*arccsc(a+b*x)))+4*b^-1*arccsc(a+b*x)*arctanh(exp(im*arccsc(a+b*x))))
@test integrate(x^-1*arccsc(a+b*x)^2, x) == :(2*Polylogarithms.polylog(3,-1*im*a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x)))+2*Polylogarithms.polylog(3,-1*im*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x)))+-1//2*Polylogarithms.polylog(3,exp(2*im*arccsc(a+b*x)))+arccsc(a+b*x)^2*log(1+im*a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x)))+arccsc(a+b*x)^2*log(1+im*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x)))+-1*arccsc(a+b*x)^2*log(1+-1*exp(2*im*arccsc(a+b*x)))+im*Polylogarithms.polylog(2,exp(2*im*arccsc(a+b*x)))*arccsc(a+b*x)+-2*im*Polylogarithms.polylog(2,-1*im*a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x)))*arccsc(a+b*x)+-2*im*Polylogarithms.polylog(2,-1*im*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x)))*arccsc(a+b*x))
@test integrate(x^-2*arccsc(a+b*x)^2, x) == :(-1*x^-1*arccsc(a+b*x)^2+-1*b*a^-1*arccsc(a+b*x)^2+-2*b*a^-1*(1+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*im*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x)))+2*b*a^-1*(1+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*im*a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x)))+-2*im*b*a^-1*(1+-1*a^2)^-1//2*arccsc(a+b*x)*log(1+im*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x)))+2*im*b*a^-1*(1+-1*a^2)^-1//2*arccsc(a+b*x)*log(1+im*a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x))))
@test integrate(x^2*arccsc(a+b*x)^3, x) == :(b^-3*Polylogarithms.polylog(3,-1*exp(im*arccsc(a+b*x)))+b^-3*arctanh((1+-1*(a+b*x)^-2)^(1/2))+-1*b^-3*Polylogarithms.polylog(3,exp(im*arccsc(a+b*x)))+1//3*x^3*arccsc(a+b*x)^3+b^-3*arccsc(a+b*x)^2*arctanh(exp(im*arccsc(a+b*x)))+b^-3*(a+b*x)*arccsc(a+b*x)+-6*a^2*b^-3*Polylogarithms.polylog(3,exp(im*arccsc(a+b*x)))+6*a^2*b^-3*Polylogarithms.polylog(3,-1*exp(im*arccsc(a+b*x)))+1//3*a^3*b^-3*arccsc(a+b*x)^3+im*b^-3*Polylogarithms.polylog(2,exp(im*arccsc(a+b*x)))*arccsc(a+b*x)+(1/2)*b^-3*(1+-1*(a+b*x)^-2)^(1/2)*(a+b*x)^2*arccsc(a+b*x)^2+-1*im*b^-3*Polylogarithms.polylog(2,-1*exp(im*arccsc(a+b*x)))*arccsc(a+b*x)+-3*im*a*b^-3*arccsc(a+b*x)^2+-3*im*a*b^-3*Polylogarithms.polylog(2,exp(2*im*arccsc(a+b*x)))+6*a*b^-3*arccsc(a+b*x)*log(1+-1*exp(2*im*arccsc(a+b*x)))+6*a^2*b^-3*arccsc(a+b*x)^2*arctanh(exp(im*arccsc(a+b*x)))+-6*im*a^2*b^-3*Polylogarithms.polylog(2,-1*exp(im*arccsc(a+b*x)))*arccsc(a+b*x)+-3*a*b^-3*(1+-1*(a+b*x)^-2)^(1/2)*arccsc(a+b*x)^2*(a+b*x)+6*im*a^2*b^-3*Polylogarithms.polylog(2,exp(im*arccsc(a+b*x)))*arccsc(a+b*x))
@test integrate(x*arccsc(a+b*x)^3, x) == :((1/2)*x^2*arccsc(a+b*x)^3+-6*a*b^-2*Polylogarithms.polylog(3,-1*exp(im*arccsc(a+b*x)))+-3*b^-2*arccsc(a+b*x)*log(1+-1*exp(2*im*arccsc(a+b*x)))+6*a*b^-2*Polylogarithms.polylog(3,exp(im*arccsc(a+b*x)))+-1//2*a^2*b^-2*arccsc(a+b*x)^3+3//2*im*b^-2*arccsc(a+b*x)^2+3//2*im*b^-2*Polylogarithms.polylog(2,exp(2*im*arccsc(a+b*x)))+(1/2)*b^-2*(1+-1*(a+b*x)^-2)^(1/2)*arccsc(a+b*x)^2*(3a+3*b*x)+-6*a*b^-2*arccsc(a+b*x)^2*arctanh(exp(im*arccsc(a+b*x)))+-6*im*a*b^-2*Polylogarithms.polylog(2,exp(im*arccsc(a+b*x)))*arccsc(a+b*x)+6*im*a*b^-2*Polylogarithms.polylog(2,-1*exp(im*arccsc(a+b*x)))*arccsc(a+b*x))
@test integrate(arccsc(a+b*x)^3, x) == :(-6*b^-1*Polylogarithms.polylog(3,exp(im*arccsc(a+b*x)))+6*b^-1*Polylogarithms.polylog(3,-1*exp(im*arccsc(a+b*x)))+b^-1*arccsc(a+b*x)^3*(a+b*x)+6*b^-1*arccsc(a+b*x)^2*arctanh(exp(im*arccsc(a+b*x)))+-6*im*b^-1*Polylogarithms.polylog(2,-1*exp(im*arccsc(a+b*x)))*arccsc(a+b*x)+6*im*b^-1*Polylogarithms.polylog(2,exp(im*arccsc(a+b*x)))*arccsc(a+b*x))
@test integrate(x^-1*arccsc(a+b*x)^3, x) == :(arccsc(a+b*x)^3*log(1+im*a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x)))+arccsc(a+b*x)^3*log(1+im*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x)))+-1*arccsc(a+b*x)^3*log(1+-1*exp(2*im*arccsc(a+b*x)))+6*im*Polylogarithms.polylog(4,-1*im*a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x)))+6*im*Polylogarithms.polylog(4,-1*im*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x)))+6*Polylogarithms.polylog(3,-1*im*a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x)))*arccsc(a+b*x)+6*Polylogarithms.polylog(3,-1*im*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x)))*arccsc(a+b*x)+-3//2*Polylogarithms.polylog(3,exp(2*im*arccsc(a+b*x)))*arccsc(a+b*x)+-3//4*im*Polylogarithms.polylog(4,exp(2*im*arccsc(a+b*x)))+-3*im*arccsc(a+b*x)^2*Polylogarithms.polylog(2,-1*im*a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x)))+-3*im*arccsc(a+b*x)^2*Polylogarithms.polylog(2,-1*im*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x)))+3//2*im*arccsc(a+b*x)^2*Polylogarithms.polylog(2,exp(2*im*arccsc(a+b*x))))
@test integrate(x^-2*arccsc(a+b*x)^3, x) == :(-1*x^-1*arccsc(a+b*x)^3+-1*b*a^-1*arccsc(a+b*x)^3+-6*im*b*a^-1*(1+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*im*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x)))+-6*b*a^-1*(1+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*im*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x)))*arccsc(a+b*x)+6*im*b*a^-1*(1+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*im*a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x)))+6*b*a^-1*(1+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*im*a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x)))*arccsc(a+b*x)+-3*im*b*a^-1*(1+-1*a^2)^-1//2*arccsc(a+b*x)^2*log(1+im*a*(1+-1*(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x)))+3*im*b*a^-1*(1+-1*a^2)^-1//2*arccsc(a+b*x)^2*log(1+im*a*(1+(1+-1*a^2)^(1/2))^-1*exp(im*arccsc(a+b*x))))
@test integrate(x^3*arccsc(a+b*x^4), x) == :(1//4*b^-1*arctanh((1+-1*(a+b*x^4)^-2)^(1/2))+1//4*b^-1*(a+b*x^4)*arccsc(a+b*x^4))
@test integrate(x^(-1+n)*arccsc(a+b*x^n), x) == :(b^-1*n^-1*arctanh((1+-1*(a+b*x^n)^-2)^(1/2))+b^-1*n^-1*(a+b*x^n)*arccsc(a+b*x^n))
@test integrate(arccsc(c*exp(a+b*x)), x) == :((1/2)*im*b^-1*arccsc(c*exp(a+b*x))^2+(1/2)*im*b^-1*Polylogarithms.polylog(2,exp(2*im*arccsc(c*exp(a+b*x))))+-1*b^-1*arccsc(c*exp(a+b*x))*log(1+-1*exp(2*im*arccsc(c*exp(a+b*x)))))
@test integrate(x^-2*exp(arccsc(a*x)), x) == :(-1//2*x^-1*exp(arccsc(a*x))+-1//2*a*(1+-1*a^-2*x^-2)^(1/2)*exp(arccsc(a*x)))
@test integrate(x^-3*exp(arccsc(a*x)), x) == :(-1//10*a^2*exp(arccsc(a*x))*sin(2*arccsc(a*x))+1//5*a^2*cos(2*arccsc(a*x))*exp(arccsc(a*x)))
@test integrate(x^-4*exp(arccsc(a*x)), x) == :(-1//8*a^2*x^-1*exp(arccsc(a*x))+-1//8*a^3*(1+-1*a^-2*x^-2)^(1/2)*exp(arccsc(a*x))+1//40*a^3*cos(3*arccsc(a*x))*exp(arccsc(a*x))+3//40*a^3*exp(arccsc(a*x))*sin(3*arccsc(a*x)))
@test integrate(x^-5*exp(arccsc(a*x)), x) == :(-1//20*a^4*exp(arccsc(a*x))*sin(2*arccsc(a*x))+-1//34*a^4*cos(4*arccsc(a*x))*exp(arccsc(a*x))+1//10*a^4*cos(2*arccsc(a*x))*exp(arccsc(a*x))+1//136*a^4*exp(arccsc(a*x))*sin(4*arccsc(a*x)))
@test integrate((d*x+a*d*b^-1)^-1*arccsc(a+b*x), x) == :((1/2)*im*d^-1*arccsc(a+b*x)^2+(1/2)*im*d^-1*Polylogarithms.polylog(2,exp(2*im*arccsc(a+b*x)))+-1*d^-1*arccsc(a+b*x)*log(1+-1*exp(2*im*arccsc(a+b*x))))

max_stars_repo_path: test/5 Inverse trig functions/5.6 Inverse cosecant/5.6.2 Inverse cosecant functions.jl


text:
using CCBlade
using Test

@testset "normal operation" begin

# --- verification using example from "Wind Turbine Blade Analysis using the Blade Element Momentum Method"
# --- by Grant Ingram: https://community.dur.ac.uk/g.l.ingram/download/wind_turbine_design.pdf

# Note: There were various problems with the pdf data. Fortunately the excel
# spreadsheet was provided: http://community.dur.ac.uk/g.l.ingram/download.php

# - They didn't actually use the NACA0012 data.  According to the spreadsheet they just used cl = 0.084*alpha with alpha in degrees.
# - There is an error in the spreadsheet where CL at the root is always 0.7.  This is because the classical approach doesn't converge properly at that station.
# - the values for gamma and chord were rounded in the pdf.  The spreadsheet has more precise values.
# - the tip is not actually converged (see significant error for Delta A).  I converged it further using their method.

# --- rotor definition ---
turbine = true
Rhub = 0.01
Rtip = 5.0
Rtip_eff = 5.0*100  # to eliminate tip effects as consistent with their study.
B = 3  # number of blades

rotor = Rotor(Rhub, Rtip_eff, B, turbine)

# --- section definitions ---

r = [0.2, 1, 2, 3, 4, 5]
gamma = [61.0, 74.31002131, 84.89805553, 89.07195504, 91.25038415, 92.58003871]
theta = (90.0 .- gamma)*pi/180
chord = [0.7, 0.706025153, 0.436187551, 0.304517933, 0.232257636, 0.187279622]

function affunc(alpha, Re, M)

    cl = 0.084*alpha*180/pi

    return cl, 0.0
end

sections = Section.(r, chord, theta, affunc)

# --- inflow definitions ---

Vinf = 7.0
tsr = 8
Omega = tsr*Vinf/Rtip
rho = 1.0

inflows = simpleinflow.(Vinf, Omega, r, rho)

# --- evaluate ---

outputs = solve.(rotor, sections, inflows)

Np, Tp = loads(outputs)
avec = getfield.(outputs, :a)
apvec = getfield.(outputs, :ap)
phivec = getfield.(outputs, :phi)

ivec = phivec*180/pi .- theta*180/pi
betavec = 90 .- phivec*180/pi

# outputs[1] is uncomparable because the classical method fails at the root so they fixed cl

@test isapprox(avec[2], 0.2443, atol=1e-4)
@test isapprox(apvec[2], 0.0676, atol=1e-4)
@test isapprox(avec[3], 0.2497, atol=1e-4)
@test isapprox(apvec[3], 0.0180, atol=1e-4)
@test isapprox(avec[4], 0.2533, atol=1e-4)
@test isapprox(apvec[4], 0.0081, atol=1e-4)
@test isapprox(avec[5], 0.2556, atol=1e-4)
@test isapprox(apvec[5], 0.0046, atol=1e-4)
@test isapprox(avec[6], 0.25725, atol=1e-4)  # note that their spreadsheet is not converged so I ran their method longer.
@test isapprox(apvec[6], 0.0030, atol=1e-4)

@test isapprox(betavec[2], 66.1354, atol=1e-3)
@test isapprox(betavec[3], 77.0298, atol=1e-3)
@test isapprox(betavec[4], 81.2283, atol=1e-3)
@test isapprox(betavec[5], 83.3961, atol=1e-3)
@test isapprox(betavec[6], 84.7113, atol=1e-3)  # using my more converged solution

#
#

# idx = 6

# aguess = avec[idx]
# apguess = apvec[idx]

# # aguess = 0.2557
# # apguess = 0.0046

# for i = 1:100

#     sigmap = B*chord[idx]/(2*pi*r[idx])
#     lambdar = inflows[idx].Vy/inflows[idx].Vx
#     beta2 = atan(lambdar*(1 + apguess)/(1 - aguess))
#     inc2 = gamma[idx]*pi/180 - beta2
#     cl2, cd2 = affunc(inc2, 1.0, 1.0)
#     global aguess = 1.0/(1 + 4*cos(beta2)^2/(sigmap*cl2*sin(beta2)))
#     global apguess = sigmap*cl2/(4*lambdar*cos(beta2))*(1 - a2)
#     # global apguess = 1.0 / (4*sin(beta2)/(sigmap*cl2) - 1)

# end

# aguess
# apguess

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

# -------- verification: propellers.  using script at http://www.aerodynamics4students.com/propulsion/blade-element-propeller-theory.php ------
# I increased their tolerance to 1e-6

# inputs
chord = 0.10
D = 1.6
RPM = 2100
rho = 1.225
pitch = 1.0  # pitch distance in meters.

# --- rotor definition ---
turbine = false
Rhub = 0.0
Rtip = D/2
Rhub_eff = 1e-6  # something small to eliminate hub effects
Rtip_eff = 100.0  # something large to eliminate tip effects
B = 2  # number of blades

rotor_no_F = Rotor(Rhub_eff, Rtip_eff, B, turbine)
rotor = Rotor(Rhub, Rtip, B, turbine)

# --- section definitions ---

R = D/2.0
r = range(R/10, stop=R, length=11)
theta = atan.(pitch./(2*pi*r))

function affunc(alpha, Re, M)

    cl = 6.2*alpha
    cd = 0.008 - 0.003*cl + 0.01*cl*cl

    return cl, cd
end

sections = Section.(r, chord, theta, affunc)

# --- inflow definitions ---

tsim = 1e3*[1.045361193032356, 1.025630300048415, 1.005234466788998, 0.984163367036026, 0.962407923825140, 0.939960208707079, 0.916813564966455, 0.892962691000145, 0.868403981825492, 0.843134981103815, 0.817154838249790, 0.790463442573673, 0.763063053839278, 0.734956576558370, 0.706148261507327, 0.676643975451150, 0.646450304160057, 0.615575090105131, 0.584027074365864, 0.551815917391907, 0.518952127358381, 0.485446691671386, 0.451311288662196, 0.416557935286392, 0.381199277009438, 0.345247916141561, 0.308716772800348, 0.271618894441869, 0.233967425339051, 0.195775319296371, 0.157055230270717, 0.117820154495231, 0.078082266879117, 0.037854059197644, -0.002852754149850, -0.044026182837742, -0.085655305814570, -0.127728999394140, -0.170237722799272, -0.213169213043848, -0.256515079286031, -0.300266519551194, -0.344414094748869, -0.388949215983616, -0.433863576642539, -0.479150401337354, -0.524801553114807, -0.570810405128802, -0.617169893200684, -0.663873474163182, -0.710915862524620, -0.758291877949762, -0.805995685105502, -0.854022273120508, -0.902366919041604, -0.951025170820984, -0.999992624287163, -1.049265666456123, -1.098840222937414, -1.148712509929845]
qsim = 1e2*[0.803638686218187, 0.806984572453978, 0.809709290183008, 0.811743686838315, 0.813015017103876, 0.813446921530685, 0.812959654049620, 0.811470393912576, 0.808893852696513, 0.805141916379142, 0.800124489784850, 0.793748780791057, 0.785921727832179, 0.776548246109426, 0.765532528164390, 0.752778882688809, 0.738190986274448, 0.721673076180745, 0.703129918771009, 0.682467282681955, 0.659592296506578, 0.634413303042323, 0.606840565246423, 0.576786093366321, 0.544164450503912, 0.508891967461804, 0.470887571011192, 0.430072787279711, 0.386371788290446, 0.339711042057184, 0.290019539402947, 0.237229503458026, 0.181274942660876, 0.122093307308376, 0.059623821454727, -0.006190834182631, -0.075406684829235, -0.148076528546541, -0.224253047813501, -0.303980950928302, -0.387309291734422, -0.474283793689904, -0.564946107631716, -0.659336973911858, -0.757495165410553, -0.859460291551374, -0.965266648683888, -1.074949504731187, -1.188540970723477, -1.306072104649531, -1.427575034895290, -1.553080300508925, -1.682614871422754, -1.816205997296014, -1.953879956474228, -2.095662107769925, -2.241576439746701, -2.391647474158875, -2.545897099743367, -2.704346566395035]

for i = 1:60

  Vinf = float(i)
  Omega = RPM * pi/30

  inflows = simpleinflow.(Vinf, Omega, r, rho)

  # --- evaluate ---

  outputs = solve.(rotor_no_F, sections, inflows)

  Np, Tp = loads(outputs)

  # T, Q = thrusttorque(r[1], r, r[end], Np, Tp, B)
  # their spreadsheet did not use trapzezoidal rule, so just do a rect sum.
  T = sum(Np*(r[2]-r[1]))*B
  Q = sum(r.*Tp*(r[2]-r[1]))*B

  @test isapprox(T, tsim[i], atol=1e-2)  # 2 decimal places
  @test isapprox(Q, qsim[i], atol=2e-3)

end

Vinf = 20.0
Omega = RPM * pi/30

inflows = simpleinflow.(Vinf, Omega, r, rho)

outputs = solve.(rotor_no_F, sections, inflows)

Np, Tp = loads(outputs)

T = sum(Np*(r[2]-r[1]))*B
Q = sum(r.*Tp*(r[2]-r[1]))*B
eff, CT, CQ = nondim(T, Q, Vinf, Omega, rho, rotor)

@test isapprox(CT, 0.056110238632657, atol=1e-7)
@test isapprox(CQ, 0.004337202960642, atol=1e-8)
@test isapprox(eff, 0.735350632777002, atol=1e-6)

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

# ------ hover verification --------

end


text:
module CancerSeqSim

using Distributions
using DataFrames
using GLM
using HypothesisTests
using StatsBase
using Printf
using Statistics
using Random
using RecipesBase

import Base.show

export
  # types
  Simulation,
  AllMetrics,
  MetricObj,
  RsqObj,
  AnalysedData,
  SampledData,
  InputParameters,
  SimResult,
  RawOutput,
  cancercell,
  bdprocess,
  SimulationSamples,

  #functions
  simulate,
  getmetrics,
  getsummary,
  vafhistogram,
  cumulativeplot,
  cumulativedist,
  birthdeathprocess,
  getsamples,
  tumourgrow_stemcell,
  run1simulationstem,
  simulatestemcells,
  simulatedifferentmutations,
  simulatedifferentmutationsmoran

### source files
include("readdata.jl")
include("runsims.jl")
include("runsimssample.jl")
include("runstemcellsim.jl")
include("runsimsdifferentmu.jl")
include("samplesim.jl")
include("util.jl")
include("multiplesamples.jl")
include("plots.jl")

end


text:
# ------------------------------------------------------------------- #
#                          FlicsFlow.jl                               #
#                       -----------------                             #
#     A General Three-dimensional Unstructured Finite Volume Solver   #
#                                                                     #
#                    Copyright(c) 2021 NUMFLICS                       #
# ------------------------------------------------------------------- #
using Pkg
Pkg.activate("C:\\Users\\arifa\\.julia\\dev\\FlicsFlow")

using FlicsFlow

# Get the current case directory
casedir = @__DIR__

# Read the OpenFOAM Mesh
mesh = readopenfoammesh(casedir)