text:
export days

#backset() = Set([:b, :back, :backward])
#fwdset() = Set([:f, :fwd, :forward])
#isforward(direction) = in(direction, fwdset()) ? true  :
#             in(direction, backset())? false :
#             throw(ArgumentError("bad argument direction=$direction"))

"""
using MT: days

days( Date(2016,04,04), 3) |> first # 2016-04-04

days( Date(2016,04,04), 3, direction=:b) |> first # 2016-04-06

days( Date(2016,04,04), 3, shift=1) |> first # 2016-04-05

days( Date(2016,04,04), 3, shift=-1) |> first # 2016-04-03

Keyword argument \"direction\" in days() and Hours() may be one of:

    :f, :fwd, :forward

    :b, :back, :backward

    f.e.:

    days( Date(2016,04,04), 3, direction=:b) |> first # 2016-04-06

"""
function days{T<:Date}(d::T, n::Int; direction::Symbol=:forward, shift::Int=0)
 forward = isforward(direction)
 nextfun = forward ? d->d+Dates.Day(1) : d->d-Dates.Day(1)
 jumpstart = n>0 && forward ? d-Dates.Day(1-shift) :
             n>0 && !forward ? d+Dates.Day(n+shift) :
             n<=0 && forward ? d+Dates.Day(n-1+shift) :
             d+Dates.Day(shift) # <-- n<=0 && !forward
 ManyTimesIter(jumpstart, abs(n), nextfun, forward, Dates.Day(1))
end

"days( Date(2016,04,04), Dates.Day(3)) |> first # 2016-04-04"
days{T<:Date}(d::T, nd::Base.Dates.Day; direction::Symbol=:forward, shift::Int=0) =
    days(d, nd.value, direction=direction, shift=shift)

"""
days( DateTime(2016,04,04), 3) |> first # 2016-04-04

days( DateTime(2016,04,04,6,5,4), 3) |> first # 2016-04-04
"""
days{T<:DateTime}(t::T, n::Int; direction::Symbol=:forward, shift::Int=0) =
    days(Date(t), n, direction=direction, shift=shift)

"days( DateTime(2016,04,04), Dates.Day(3)) |> first # 2016-04-04"
days{T<:DateTime}(t::T, nd::Base.Dates.Day; direction::Symbol=:forward, shift::Int=0) =
    days(Date(t), nd.value, direction=direction, shift=shift)

"""
days(Date(2016,04,04), Date(2016,04,06)) |> first # 2016-04-04

days(Date(2016,04,04), Date(2016,04,06), shift=1) |> first # 2016-04-05

days(Date(2016,04,06), Date(2016,04,04)) |> first # 2016-04-05

days(Date(2016,04,06), Date(2016,04,04), shift=1) |> first # 2016-04-06
"""
function days{T<:Date}(d1::T, d2::T; shift::Int=0)
 if d2 > d1
    dstart = d1
    direction = :forward
 else
    dstart = d2
    direction = :backward
 end
 n = abs( Int( Dates.Day( d2 - d1 )))
 days(dstart, n, direction=direction, shift=shift)
end

"days(Date(2016,04,04), DateTime(2016,04,06,12,20,0)) |> first # 2016-04-04"
days{ T1<:Date, T2<:DateTime }(d1::T1, t2::T2; shift::Int=0) =
    days(d1, Date(t2), shift=shift)

"days(DateTime(2016,04,04,12,30,0), Date(2016,04,06)) |> first # 2016-04-04"
days{ T1<:DateTime, T2<:Date }(t1::T1, d2::T2; shift::Int=0) =
    days(Date(t1), d2, shift=shift)

"""
days(DateTime(2016,04,04,12,30,0), DateTime(2016,04,06,12,20,0)) |>
    first # 2016-04-04
"""
days{ T1<:DateTime, T2<:DateTime }(t1::T1, t2::T2; shift::Int=0) =
    days(Date(t1), Date(t2), shift=shift)

"""
Days from tuples:

days((2016,08,15),3)

days((2016,08,15),(2016-08-17))
"""
days(d1::Tuple{Int64,Int64,Int64}, d2::Tuple{Int64,Int64,Int64}, args...) =
 days(Date(d1...), Date(d2...), args...)

days(d1::Tuple{Int64,Int64,Int64}, args...) =
 days(Date(d1...), args...)

max_stars_repo_name: closescreen/Many

max_issues_repo_name: closescreen/Many

max_forks_repo_name: closescreen/Many


text:
using RefractiveIndexDatabase
using Test

@testset "RefractiveIndexDatabase.jl" begin
    @testset "Dispersion formulas" begin
        # Sellmeier
        @test isapprox(RefractiveIndex(get_material("main", "Ar", "Grace-liquid-90K"))(0.3885), 1.22809258, rtol=1e-6)

        # Sellmeier-2
        @test isapprox(RefractiveIndex(get_material("main", "CdTe", "Marple"))(1.68), 2.72726378, rtol=1e-6)

        # Polynomial
        @test isapprox(RefractiveIndex(get_material("other", "Mg-LiTaO3", "Moutzouris-o"))(1.0), 2.13370331, rtol=1e-6)

        # RefractiveIndex.INFO
        @test isapprox(RefractiveIndex(get_material("main", "ZnTe", "Li"))(15.275), 2.66048474, rtol=1e-6)

        # Cauchy
        @test isapprox(RefractiveIndex(get_material("main", "SF6", "Vukovic"))(0.95), 1.00072070, rtol=1e-6)

        # Gases
        @test isapprox(RefractiveIndex(get_material("main", "He", "Mansfield"))(1.26935), 1.00003472, rtol=1e-6)

        # Herzberger
        @test isapprox(RefractiveIndex(get_material("main", "Si", "Edwards"))(13.71865), 3.42084598, rtol=1e-6)

        # Retro
        @test isapprox(RefractiveIndex(get_material("main", "AgBr", "Schröter"))(0.5825), 2.26004419, rtol=1e-6)

        # Exotic
        @test isapprox(RefractiveIndex(get_material("organic", "urea", "Rosker-e"))(0.68), 1.60004980, rtol=1e-6)
    end

    @testset "Tabular data" begin
        # FormulaNK
        m = RefractiveIndex(get_material("main","Si","Chandler-Horowitz"))
        @test isapprox(m(10.), 3.41807041 + 7.4e-5im,  rtol=1e-6)

        # TabulatedNK
        m = RefractiveIndex(get_material("main", "ZnO", "Stelling"))
        @test isapprox(m(0.866), 1.59478654 + 0.01177447im,  rtol=1e-6)

    end

    @testset "Load data from url" begin

        url = "https://refractiveindex.info/database/data/main/Ag/Johnson.yml"
        m = RefractiveIndex(load_url(url))
        @test isapprox(m(0.49), 0.05087409 + 3.04074578im,  rtol=1e-6)

    end
    @testset "Load data from file" begin

        path = joinpath(@__DIR__, "data", "BAF2.yml")
        m = RefractiveIndex(load_file(path))
        @test isapprox(m(0.82), 1.56104714 + 1.28383192e-8im,  rtol=1e-6)

    end

end

ext: jl


text:
struct IsBitsStruct
    a::Int32
    b::Int64
    c::Float32
end

mutable struct MutableBitsStruct
    a::Int32
    b::Int32
    c::Int64
end

@testset "Testing Persist" begin
    import Persistence.Lib
    using Persistence.Persist

    poolfile = "test.pool"
    ispath(poolfile) && rm(poolfile)

    pool = Lib.create(poolfile, 2^24)

    #####
    ##### isbits testing
    #####

    # int64
    A = 1
    p = persist(pool, A)
    @test isa(p, Persistent{Int})
    @test Lib.isnull(Persist.getoid(p)) == false
    B = retrieve(p)
    @test A === B

    # float32
    A = Float32(10)
    p = persist(pool, A)
    @test isa(p, Persistent{Float32})
    @test Lib.isnull(Persist.getoid(p)) == false
    B = retrieve(p)
    @test A === B

    # isbits tuple
    A = (1, 1.0, UInt8(5), Float32(-1.0), Int16(-100))
    p = persist(pool, A)
    @test Lib.isnull(Persist.getoid(p)) == false
    B = retrieve(p)
    @test A === B

    # custom `isbits` struct
    @test isbitstype(IsBitsStruct)
    A = IsBitsStruct(typemin(Int32), typemax(Int64), -0.0)
    p = persist(pool, A)
    B = retrieve(p)
    @test A === B

    #####
    ##### string testing
    #####

    # String
    A = "hello world"
    p = persist(pool, A)
    B = retrieve(p)
    @test A == B

    # Symbol
    A = :symbol
    p = persist(pool, A)
    B = retrieve(p)
    @test typeof(A) == typeof(B)
    @test A === B

    #####
    ##### array testing
    #####

end


text:
module Systems

using SymEngine
using ..HopfNormalForm: calculate_normal_form
const I=1im
const i=1im

function example1()
    norder = 3; # Choose the order of power series representation of Centermanifold
    cm = 3; # Dimension of Centermanifold of the dynamical system
    m1 = 1; # Number of the zero eigenvalue corresponding to the center space of the Jacobian
    m2 = 1; # Number of pair of complex conjugate eigenvalue corresponding to the center space of the Jacobian
    m3 = 4; # Dimension of the stable manifold of the dynmical system
    n = m1+2*m2+m3; # Total dimension of the dynamical system
    ω_h = 15.4793; # Eigenvalue that crosses unit circle at Hopf bifurcation point ± ω_h
    zero_tol = 1e-13; # zero tolerance to chop out very small coeffs of symbolic computation

    v = [symbols("v_$i") for i in 1:n]; # Define state-space variables
    f = Array{Basic}(undef, n);

    @vars ka2, ka3, ν  # Define symbolic variables which is needed to define the equation of motion
    # In this example the variables are nonlinear stifness parameters

    # Define the equation of motion at the Hopf bifurcation point where the dimension of centermanifold is 2
    # Let the dimension of control parameter vector be 1, and express it as a one of state-space variable v[1]
    # d/dt v[1]=0 and let rest of the state-space variable be v[2],v[3],…,v[n].
    f[1]=0

    f[2]=v[2]*15.4793i - 1.0*v[5]*(- v[1]^3*(0.0000344811 - 0.0000411838i) + v[1]^2*(0.00239738 + 0.00406001i) + v[1]*(0.128199 + 0.108645i)) - 1.0*v[7]*(v[1]^3*(0.00028883 - 0.000720132i) + v[1]^2*(0.0132244 - 0.0199877i) + v[1]*(0.167761 - 0.219131i)) + ka2*(- v[2]*v[3]*(0.00329 + 0.00786i) + v[2]*v[4]*(0.00457 + 0.00683i) - v[2]*v[5]*(0.00165 + 0.00805i) - v[3]*v[4]*(0.00457 + 0.00683i) - v[2]*v[6]*(0.00018 - 0.00007i) + v[3]*v[5]*(0.00165 + 0.00805i) - v[2]*v[7]*(0.00006 - 0.00002i) + v[3]*v[6]*(0.00018 - 0.00007i) - v[4]*v[5]*(0.00306 + 0.00731i) + v[3]*v[7]*(0.00006 - 0.00002i) - v[4]*v[6]*(0.00015 - 0.0001i) - v[4]*v[7]*(0.00005 - 0.00003i) + v[5]*v[6]*(0.00018 - 0.00004i) + v[5]*v[7]*(0.00006 - 0.00001i) + v[2]^2*(0.00164 + 0.00393i) + v[3]^2*(0.00164 + 0.00393i) + v[4]^2*(0.00279 + 0.00281i) + v[5]^2*(0.00004 + 0.00396i)) - 1.0*v[6]*(- v[1]^3*(0.0000447835 - 0.000111657i) + v[1]^2*(0.00116393 - 0.00491519i) + v[1]*(0.113807 - 0.0877385i)) + v[4]*(v[1]^3*(0.00000353325 - 0.0000535964i) + v[1]^2*(0.00520619 - 0.00123418i) + v[1]*(0.167966 + 0.0119843i)) - 1.0*v[3]*(- v[1]^3*(6.45125e-7 - 0.000040738i) + v[1]^2*(0.00381858 + 0.00281874i) + v[1]*(0.150471 + 0.0573303i)) + v[2]*(v[1]^3*(0.0000276837 - 0.0000298934i) + v[1]^2*(0.00508047 + 0.000946023i) + v[1]*(0.136146 + 0.057017i)) + ka3*(v[2]*v[3]^2*(0.0007 - 0.00029i) - v[2]^2*v[3]*(0.0007 - 0.00029i) + v[2]*v[4]^2*(0.0005 - 0.0005i) + v[2]^2*v[4]*(0.00061 - 0.00041i) + v[2]*v[5]^2*(0.00071 - 0.00001i) - v[3]*v[4]^2*(0.0005 - 0.0005i) - v[2]^2*v[5]*(0.00072 - 0.00015i) + v[3]^2*v[4]*(0.00061 - 0.00041i) - v[3]*v[5]^2*(0.00071 - 0.00001i) + v[2]^2*v[6]*(0.00001 + 0.00002i) - v[3]^2*v[5]*(0.00072 - 0.00015i) + v[4]*v[5]^2*(0.00067 - 0.00014i) + v[3]^2*v[6]*(0.00001 + 0.00002i) - v[4]^2*v[5]*(0.00057 - 0.00038i) + v[4]^2*v[6]*(0.00001 + 0.00001i) + v[5]^2*v[6]*0.00002i + v[2]^3*(0.00023 - 0.0001i) - v[3]^3*(0.00023 - 0.0001i) + v[4]^3*(0.00013 - 0.00019i) - v[5]^3*(0.00022 + 0.00004i) - v[2]*v[3]*v[4]*(0.00122 - 0.00082i) + v[2]*v[3]*v[5]*(0.00144 - 0.0003i) - v[2]*v[3]*v[6]*(0.00001 + 0.00003i) - v[2]*v[4]*v[5]*(0.00131 - 0.00055i) - v[2]*v[3]*v[7]*0.00001i + v[2]*v[4]*v[6]*(0.00002 + 0.00003i) + v[3]*v[4]*v[5]*(0.00131 - 0.00055i) + v[2]*v[4]*v[7]*(0.00001 + 0.00001i) - v[2]*v[5]*v[6]*(0.00001 + 0.00003i) - v[3]*v[4]*v[6]*(0.00002 + 0.00003i) - v[2]*v[5]*v[7]*0.00001i - v[3]*v[4]*v[7]*(0.00001 + 0.00001i) + v[3]*v[5]*v[6]*(0.00001 + 0.00003i) + v[3]*v[5]*v[7]*0.00001i - v[4]*v[5]*v[6]*(0.00001 + 0.00003i) - v[4]*v[5]*v[7]*0.00001i)

    f[3]=- v[3]*15.4793i - 1.0*v[4]*(- v[1]^3*(0.0000344811 + 0.0000411838i) + v[1]^2*(0.00239738 - 0.00406001i) + v[1]*(0.128199 - 0.108645i)) - 1.0*v[7]*(v[1]^3*(0.00028883 + 0.000720132i) + v[1]^2*(0.0132244 + 0.0199877i) + v[1]*(0.167761 + 0.219131i)) + ka2*(- v[2]*v[3]*(0.00329 - 0.00786i) + v[2]*v[4]*(0.00165 - 0.00805i) - v[2]*v[5]*(0.00457 - 0.00683i) - v[3]*v[4]*(0.00165 - 0.00805i) + v[2]*v[6]*(0.00018 + 0.00007i) + v[3]*v[5]*(0.00457 - 0.00683i) + v[2]*v[7]*(0.00006 + 0.00002i) - v[3]*v[6]*(0.00018 + 0.00007i) - v[4]*v[5]*(0.00306 - 0.00731i) - v[3]*v[7]*(0.00006 + 0.00002i) + v[4]*v[6]*(0.00018 + 0.00004i) + v[4]*v[7]*(0.00006 + 0.00001i) - v[5]*v[6]*(0.00015 + 0.0001i) - v[5]*v[7]*(0.00005 + 0.00003i) + v[2]^2*(0.00164 - 0.00393i) + v[3]^2*(0.00164 - 0.00393i) + v[4]^2*(0.00004 - 0.00396i) + v[5]^2*(0.00279 - 0.00281i)) - 1.0*v[6]*(- v[1]^3*(0.0000447835 + 0.000111657i) + v[1]^2*(0.00116393 + 0.00491519i) + v[1]*(0.113807 + 0.0877385i)) + v[5]*(v[1]^3*(0.00000353325 + 0.0000535964i) + v[1]^2*(0.00520619 + 0.00123418i) + v[1]*(0.167966 - 0.0119843i)) - 1.0*v[2]*(- v[1]^3*(6.45125e-7 + 0.000040738i) + v[1]^2*(0.00381858 - 0.00281874i) + v[1]*(0.150471 - 0.0573303i)) + v[3]*(v[1]^3*(0.0000276837 + 0.0000298934i) + v[1]^2*(0.00508047 - 0.000946023i) + v[1]*(0.136146 - 0.057017i)) + ka3*(- v[2]*v[3]^2*(0.0007 + 0.00029i) + v[2]^2*v[3]*(0.0007 + 0.00029i) - v[2]*v[4]^2*(0.00071 + 0.00001i) - v[2]^2*v[4]*(0.00072 + 0.00015i) - v[2]*v[5]^2*(0.0005 + 0.0005i) + v[3]*v[4]^2*(0.00071 + 0.00001i) + v[2]^2*v[5]*(0.00061 + 0.00041i) - v[3]^2*v[4]*(0.00072 + 0.00015i) + v[3]*v[5]^2*(0.0005 + 0.0005i) + v[2]^2*v[6]*(0.00001 - 0.00002i) + v[3]^2*v[5]*(0.00061 + 0.00041i) - v[4]*v[5]^2*(0.00057 + 0.00038i) + v[3]^2*v[6]*(0.00001 - 0.00002i) + v[4]^2*v[5]*(0.00067 + 0.00014i) - v[4]^2*v[6]*0.00002i + v[5]^2*v[6]*(0.00001 - 0.00001i) - v[2]^3*(0.00023 + 0.0001i) + v[3]^3*(0.00023 + 0.0001i) - v[4]^3*(0.00022 - 0.00004i) + v[5]^3*(0.00013 + 0.00019i) + v[2]*v[3]*v[4]*(0.00144 + 0.0003i) - v[2]*v[3]*v[5]*(0.00122 + 0.00082i) - v[2]*v[3]*v[6]*(0.00001 - 0.00003i) + v[2]*v[4]*v[5]*(0.00131 + 0.00055i) + v[2]*v[3]*v[7]*0.00001i + v[2]*v[4]*v[6]*(0.00001 - 0.00003i) - v[3]*v[4]*v[5]*(0.00131 + 0.00055i) - v[2]*v[4]*v[7]*0.00001i - v[2]*v[5]*v[6]*(0.00002 - 0.00003i) - v[3]*v[4]*v[6]*(0.00001 - 0.00003i) - v[2]*v[5]*v[7]*(0.00001 - 0.00001i) + v[3]*v[4]*v[7]*0.00001i + v[3]*v[5]*v[6]*(0.00002 - 0.00003i) + v[3]*v[5]*v[7]*(0.00001 - 0.00001i) - v[4]*v[5]*v[6]*(0.00001 - 0.00003i) + v[4]*v[5]*v[7]*0.00001i)

    f[4]=- v[4]*(3.22105 - 16.4334i) + ka3*(0.00001*v[2]^2*v[7] + v[2]^2*v[3]*(0.00042 + 0.00167i) - v[2]*v[4]^2*(0.00095 + 0.00129i) - v[2]^2*v[4]*(0.00071 + 0.0015i) + v[2]*v[5]^2*(0.00022 - 0.00158i) + v[3]*v[4]^2*(0.00095 + 0.00129i) + v[2]^2*v[5]*(0.00009 + 0.00165i) - v[3]^2*v[4]*(0.00071 + 0.0015i) - v[3]*v[5]^2*(0.00022 - 0.00158i) + v[2]^2*v[6]*(0.00004 - 0.00001i) + v[3]^2*v[5]*(0.00009 + 0.00165i) - v[4]*v[5]^2*(0.00008 + 0.00154i) - v[2]*v[3]^2*(0.00042 + 0.00167i) + v[3]^2*v[6]*(0.00004 - 0.00001i) + v[4]^2*v[5]*(0.00066 + 0.00139i) + 0.00001*v[3]^2*v[7] + v[4]^2*v[6]*(0.00003 - 0.00002i) + v[4]^2*v[7]*(0.00001 - 0.00001i) + 0.00004*v[5]^2*v[6] + 0.00001*v[5]^2*v[7] - v[2]^3*(0.00014 + 0.00056i) + v[3]^3*(0.00014 + 0.00056i) - v[4]^3*(0.00038 + 0.00035i) - v[5]^3*(0.00017 - 0.00049i) + v[2]*v[3]*v[4]*(0.00142 + 0.003i) - v[2]*v[3]*v[5]*(0.00018 + 0.00331i) - v[2]*v[3]*v[6]*(0.00007 - 0.00002i) + v[2]*v[4]*v[5]*(0.00078 + 0.0031i) - v[2]*v[3]*v[7]*(0.00002 - 0.00001i) + v[2]*v[4]*v[6]*(0.00007 - 0.00003i) - v[3]*v[4]*v[5]*(0.00078 + 0.0031i) + v[2]*v[4]*v[7]*(0.00002 - 0.00001i) - 0.00007*v[2]*v[5]*v[6] - v[3]*v[4]*v[6]*(0.00007 - 0.00003i) - 0.00002*v[2]*v[5]*v[7] - v[3]*v[4]*v[7]*(0.00002 - 0.00001i) + 0.00007*v[3]*v[5]*v[6] + 0.00002*v[3]*v[5]*v[7] - v[4]*v[5]*v[6]*(0.00007 - 0.00002i) - v[4]*v[5]*v[7]*(0.00002 - 0.00001i)) - 1.0*v[2]*(v[1]^3*(0.0000420472 - 0.0000200759i) + v[1]^2*(0.00501926 + 0.00319374i) + v[1]*(0.152099 + 0.126793i)) + v[7]*(v[1]^3*(0.000611125 - 0.000643319i) + v[1]^2*(0.0231344 - 0.0150345i) + v[1]*(0.280706 - 0.152463i)) + v[6]*(- v[1]^3*(0.0000947556 - 0.0000997473i) + v[1]^2*(0.00321416 - 0.00482834i) + v[1]*(0.169289 - 0.0352994i)) + v[3]*(- v[1]^3*(0.0000178927 - 0.0000430216i) + v[1]^2*(0.00282464 + 0.00459i) + v[1]*(0.104275 + 0.115717i)) + ka2*(- v[2]*v[3]*(0.01859 - 0.00471i) + v[2]*v[4]*(0.01672 - 0.00791i) - v[2]*v[5]*(0.01846 - 0.00101i) - v[3]*v[4]*(0.01672 - 0.00791i) + v[2]*v[6]*(0.00011 + 0.00042i) + v[3]*v[5]*(0.01846 - 0.00101i) + v[2]*v[7]*(0.00003 + 0.00013i) - v[3]*v[6]*(0.00011 + 0.00042i) - v[4]*v[5]*(0.01728 - 0.00438i) - v[3]*v[7]*(0.00003 + 0.00013i) + v[4]*v[6]*(0.00018 + 0.00037i) + v[4]*v[7]*(0.00006 + 0.00012i) - v[5]*v[6]*(0.00002 + 0.00041i) - v[5]*v[7]*(0.00001 + 0.00013i) + v[2]^2*(0.0093 - 0.00236i) + v[3]^2*(0.0093 - 0.00236i) + v[4]^2*(0.00717 - 0.00529i) + v[5]^2*(0.00883 + 0.00123i)) + v[5]*(- v[1]^3*(0.00005404 - 0.0000292036i) + v[1]^2*(0.000804306 + 0.00527276i) + v[1]*(0.0564505 + 0.15966i)) - 1.0*v[4]*(v[1]^3*(0.0000263932 - 0.0000554668i) + v[1]^2*(0.00610947 + 0.000914594i) + v[1]*(0.207125 + 0.0938587i))

    f[5]=- v[5]*(3.22105 + 16.4334i) + ka3*(v[2]*v[3]^2*(0.00042 - 0.00167i) - v[2]^2*v[3]*(0.00042 - 0.00167i) - v[2]*v[4]^2*(0.00022 + 0.00158i) + v[2]^2*v[4]*(0.00009 - 0.00165i) + v[2]*v[5]^2*(0.00095 - 0.00129i) + v[3]*v[4]^2*(0.00022 + 0.00158i) - v[2]^2*v[5]*(0.00071 - 0.0015i) + v[3]^2*v[4]*(0.00009 - 0.00165i) - v[3]*v[5]^2*(0.00095 - 0.00129i) + v[2]^2*v[6]*(0.00004 + 0.00001i) - v[3]^2*v[5]*(0.00071 - 0.0015i) + v[4]*v[5]^2*(0.00066 - 0.00139i) + 0.00001*v[2]^2*v[7] + v[3]^2*v[6]*(0.00004 + 0.00001i) - v[4]^2*v[5]*(0.00008 - 0.00154i) + 0.00001*v[3]^2*v[7] + 0.00004*v[4]^2*v[6] + 0.00001*v[4]^2*v[7] + v[5]^2*v[6]*(0.00003 + 0.00002i) + v[5]^2*v[7]*(0.00001 + 0.00001i) + v[2]^3*(0.00014 - 0.00056i) - v[3]^3*(0.00014 - 0.00056i) - v[4]^3*(0.00017 + 0.00049i) - v[5]^3*(0.00038 - 0.00035i) - v[2]*v[3]*v[4]*(0.00018 - 0.00331i) + v[2]*v[3]*v[5]*(0.00142 - 0.003i) - v[2]*v[3]*v[6]*(0.00007 + 0.00002i) - v[2]*v[4]*v[5]*(0.00078 - 0.0031i) - v[2]*v[3]*v[7]*(0.00002 + 0.00001i) + 0.00007*v[2]*v[4]*v[6] + v[3]*v[4]*v[5]*(0.00078 - 0.0031i) + 0.00002*v[2]*v[4]*v[7] - v[2]*v[5]*v[6]*(0.00007 + 0.00003i) - 0.00007*v[3]*v[4]*v[6] - v[2]*v[5]*v[7]*(0.00002 + 0.00001i) - 0.00002*v[3]*v[4]*v[7] + v[3]*v[5]*v[6]*(0.00007 + 0.00003i) + v[3]*v[5]*v[7]*(0.00002 + 0.00001i) - v[4]*v[5]*v[6]*(0.00007 + 0.00002i) - v[4]*v[5]*v[7]*(0.00002 + 0.00001i)) - 1.0*v[3]*(v[1]^3*(0.0000420472 + 0.0000200759i) + v[1]^2*(0.00501926 - 0.00319374i) + v[1]*(0.152099 - 0.126793i)) + v[7]*(v[1]^3*(0.000611125 + 0.000643319i) + v[1]^2*(0.0231344 + 0.0150345i) + v[1]*(0.280706 + 0.152463i)) + v[6]*(- v[1]^3*(0.0000947556 + 0.0000997473i) + v[1]^2*(0.00321416 + 0.00482834i) + v[1]*(0.169289 + 0.0352994i)) + v[2]*(- v[1]^3*(0.0000178927 + 0.0000430216i) + v[1]^2*(0.00282464 - 0.00459i) + v[1]*(0.104275 - 0.115717i)) + ka2*(- v[2]*v[3]*(0.01859 + 0.00471i) + v[2]*v[4]*(0.01846 + 0.00101i) - v[2]*v[5]*(0.01672 + 0.00791i) - v[3]*v[4]*(0.01846 + 0.00101i) - v[2]*v[6]*(0.00011 - 0.00042i) + v[3]*v[5]*(0.01672 + 0.00791i) - v[2]*v[7]*(0.00003 - 0.00013i) + v[3]*v[6]*(0.00011 - 0.00042i) - v[4]*v[5]*(0.01728 + 0.00438i) + v[3]*v[7]*(0.00003 - 0.00013i) - v[4]*v[6]*(0.00002 - 0.00041i) - v[4]*v[7]*(0.00001 - 0.00013i) + v[5]*v[6]*(0.00018 - 0.00037i) + v[5]*v[7]*(0.00006 - 0.00012i) + v[2]^2*(0.0093 + 0.00236i) + v[3]^2*(0.0093 + 0.00236i) + v[4]^2*(0.00883 - 0.00123i) + v[5]^2*(0.00717 + 0.00529i)) + v[4]*(- v[1]^3*(0.00005404 + 0.0000292036i) + v[1]^2*(0.000804306 - 0.00527276i) + v[1]*(0.0564505 - 0.15966i)) - 1.0*v[5]*(v[1]^3*(0.0000263932 + 0.0000554668i) + v[1]^2*(0.00610947 - 0.000914594i) + v[1]*(0.207125 - 0.0938587i))

    f[6]=- 35.42*v[6] + ka2*(0.00156*v[2]^2 + v[2]*v[4]*(0.00296 - 0.00058i) - v[2]*v[5]*(0.00296 + 0.00058i) - v[3]*v[4]*(0.00296 - 0.00058i) + v[2]*v[6]*0.00007i + v[3]*v[5]*(0.00296 + 0.00058i) + v[2]*v[7]*0.00002i - v[3]*v[6]*0.00007i - 0.00291*v[4]*v[5] - v[3]*v[7]*0.00002i + v[4]*v[6]*(0.00001 + 0.00007i) + v[4]*v[7]*0.00002i + v[5]*v[6]*(0.00001 - 0.00007i) - v[5]*v[7]*0.00002i - 0.00313*v[2]*v[3] + 0.00156*v[3]^2 + v[4]^2*(0.00135 - 0.00055i) + v[5]^2*(0.00135 + 0.00055i)) + v[4]*(- v[1]^3*(0.0000100928 - 0.00000329552i) + v[1]^2*(0.00800912 - 0.00360892i) + v[1]*(0.277632 - 0.00848122i)) + v[5]*(- v[1]^3*(0.0000100928 + 0.00000329552i) + v[1]^2*(0.00800912 + 0.00360892i) + v[1]*(0.277632 + 0.00848122i)) - 1.0*ka3*(v[2]^3*0.00009i - v[2]^2*v[3]*0.00028i + v[2]^2*v[4]*(0.00005 + 0.00027i) + v[2]^2*v[5]*(0.00005 - 0.00027i) - 0.00001*v[6]*v[2]^2 + v[2]*v[3]^2*0.00028i - v[2]*v[3]*v[4]*(0.0001 + 0.00053i) - v[2]*v[3]*v[5]*(0.0001 - 0.00053i) + 0.00001*v[6]*v[2]*v[3] + v[2]*v[4]^2*(0.0001 + 0.00024i) - v[2]*v[4]*v[5]*0.00052i - 0.00001*v[6]*v[2]*v[4] - v[2]*v[5]^2*(0.0001 - 0.00024i) + 0.00001*v[6]*v[2]*v[5] - v[3]^3*0.00009i + v[3]^2*v[4]*(0.00005 + 0.00027i) + v[3]^2*v[5]*(0.00005 - 0.00027i) - 0.00001*v[6]*v[3]^2 - v[3]*v[4]^2*(0.0001 + 0.00024i) + v[3]*v[4]*v[5]*0.00052i + 0.00001*v[6]*v[3]*v[4] + v[3]*v[5]^2*(0.0001 - 0.00024i) - 0.00001*v[6]*v[3]*v[5] + v[4]^3*(0.00005 + 0.00007i) - v[4]^2*v[5]*(0.00005 + 0.00025i) - 0.00001*v[6]*v[4]^2 - v[4]*v[5]^2*(0.00005 - 0.00025i) + 0.00001*v[6]*v[4]*v[5] + v[5]^3*(0.00005 - 0.00007i) - 0.00001*v[6]*v[5]^2) - 1.0*v[7]*(- v[1]^3*(0.000153368 - 3.3098e-20i) + v[1]^2*(0.126015 + 1.35887e-18i) + v[1]*(1.95201 + 1.68276e-17i)) + v[2]*(- v[1]^3*(0.00000752126 + 0.00000287923i) + v[1]^2*(0.0061892 + 0.00149802i) + v[1]*(0.0833156 - 0.0381373i)) + v[3]*(- v[1]^3*(0.00000752126 - 0.00000287923i) + v[1]^2*(0.0061892 - 0.00149802i) + v[1]*(0.0833156 + 0.0381373i)) - 1.0*v[6]*(v[1]^3*(0.0000237799 - 5.13187e-21i) - v[1]^2*(0.0212456 - 1.57651e-19i) + v[1]*(1.96145 + 1.07687e-17i))

    f[7]=- 5.41329*v[7] - 1.0*v[6]*(v[1]^3*(0.0000032869 - 4.60936e-22i) - v[1]^2*(0.00321098 + 2.22918e-19i) + v[1]*(0.297817 + 2.42034e-17i)) + v[4]*(- v[1]^3*(0.00000139504 - 4.55513e-7i) + v[1]^2*(0.00122348 - 0.000536855i) + v[1]*(0.0383104 - 0.000859108i)) + v[5]*(- v[1]^3*(0.00000139504 + 4.55513e-7i) + v[1]^2*(0.00122348 + 0.000536855i) + v[1]*(0.0383104 + 0.000859108i)) + v[2]*(- v[1]^3*(0.0000010396 + 3.97973e-7i) + v[1]^2*(0.000942262 + 0.000240279i) + v[1]*(0.00883636 - 0.00537487i)) + v[3]*(- v[1]^3*(0.0000010396 - 3.97973e-7i) + v[1]^2*(0.000942262 - 0.000240279i) + v[1]*(0.00883636 + 0.00537487i)) - 1.0*v[7]*(- v[1]^3*(0.0000211988 - 2.9728e-21i) + v[1]^2*(0.0191876 + 1.65108e-18i) + v[1]*(0.29695 + 2.49701e-17i)) + ka2*(0.00138*v[2]^2 + v[2]*v[4]*(0.00261 - 0.00051i) - v[2]*v[5]*(0.00261 + 0.00051i) - v[3]*v[4]*(0.00261 - 0.00051i) + v[2]*v[6]*0.00006i + v[3]*v[5]*(0.00261 + 0.00051i) + v[2]*v[7]*0.00002i - v[3]*v[6]*0.00006i - 0.00256*v[4]*v[5] - v[3]*v[7]*0.00002i + v[4]*v[6]*(0.00001 + 0.00006i) + v[4]*v[7]*0.00002i + v[5]*v[6]*(0.00001 - 0.00006i) - v[5]*v[7]*0.00002i - 0.00276*v[2]*v[3] + 0.00138*v[3]^2 + v[4]^2*(0.00119 - 0.00048i) + v[5]^2*(0.00119 + 0.00048i)) - 1.0*ka3*(v[2]^3*0.00008i - v[2]^2*v[3]*0.00025i + v[2]^2*v[4]*(0.00005 + 0.00023i) + v[2]^2*v[5]*(0.00005 - 0.00023i) - 0.00001*v[6]*v[2]^2 + v[2]*v[3]^2*0.00025i - v[2]*v[3]*v[4]*(0.00009 + 0.00047i) - v[2]*v[3]*v[5]*(0.00009 - 0.00047i) + 0.00001*v[6]*v[2]*v[3] + v[2]*v[4]^2*(0.00009 + 0.00021i) - v[2]*v[4]*v[5]*0.00046i - 0.00001*v[6]*v[2]*v[4] - v[2]*v[5]^2*(0.00009 - 0.00021i) + 0.00001*v[6]*v[2]*v[5] - v[3]^3*0.00008i + v[3]^2*v[4]*(0.00005 + 0.00023i) + v[3]^2*v[5]*(0.00005 - 0.00023i) - 0.00001*v[6]*v[3]^2 - v[3]*v[4]^2*(0.00009 + 0.00021i) + v[3]*v[4]*v[5]*0.00046i + 0.00001*v[6]*v[3]*v[4] + v[3]*v[5]^2*(0.00009 - 0.00021i) - 0.00001*v[6]*v[3]*v[5] + v[4]^3*(0.00004 + 0.00006i) - v[4]^2*v[5]*(0.00004 + 0.00022i) - v[4]*v[5]^2*(0.00004 - 0.00022i) + 0.00001*v[6]*v[4]*v[5] + v[5]^3*(0.00004 - 0.00006i))

    f[2] = f[2]/ω_h; f[3] = f[3]/ω_h; f[4] = f[4]/ω_h; f[5] = f[5]/ω_h; f[6] = f[6]/ω_h; f[7] = f[7]/ω_h

    return calculate_normal_form(f, v, norder, cm, m1, m2, m3, n, ω_h, zero_tol)
end

function example2()
    norder = 3; # Choose the order of power series representation of Centermanifold
    cm = 3; # Dimension of Centermanifold of the dynamical system
    m1 = 1; # Number of the zero eigenvalue corresponding to the center space of the Jacobian
    m2 = 1; # Number of pair of complex conjugate eigenvalue corresponding to the center space of the Jacobian
    m3 = 4; # Dimension of the stable manifold of the dynmical system
    n = m1+2*m2+m3; # Total dimension of the dynamical system
    ω_h = 15.1829; # Eigenvalue that crosses unit circle at Hopf bifurcation point ± ω_h
    zero_tol = 1e-13; # zero tolerance to chop out very small coeffs of symbolic computation

    v = [symbols("v_$i") for i in 1:n]; # Define state-space variables
    f = Array{Basic}(undef, n);

    @vars ka2, ka3, ν  # Define symbolic variables which is needed to define the equation of motion
    # In this example the variables are nonlinear stifness parameters

    # Define the equation of motion at the Hopf bifurcation point where the dimension of centermanifold is 2
    # Let the dimension of control parameter vector be 1, and express it as a one of state-space variable v[1]
    # d/dt v[1]=0 and let rest of the state-space variable be v[2],v[3],…,v[n].

    f[1]=0

    f[1]=0

    f[2]=v[2]*15.1829i + v[4]*(v[1]^3*(0.00000311958 - 0.000034438i) + v[1]^2*(0.00343708 - 0.00109031i) + v[1]*(0.143011 + 0.00225607i)) - 1.0*v[7]*(v[1]^3*(0.000169667 - 0.000348283i) + v[1]^2*(0.009276 - 0.0122307i) + v[1]*(0.154139 - 0.176441i)) - 1.0*ka2*(v[2]*v[3]*(0.00211 + 0.00537i) - v[2]*v[4]*(0.00357 + 0.00452i) + v[2]*v[5]*(0.00047 + 0.00574i) + v[3]*v[4]*(0.00357 + 0.00452i) - v[2]*v[6]*(0.00015 - 0.00006i) - v[3]*v[5]*(0.00047 + 0.00574i) + v[2]*v[7]*(0.00012 - 0.00005i) + v[3]*v[6]*(0.00015 - 0.00006i) + v[4]*v[5]*(0.00211 + 0.00535i) - v[3]*v[7]*(0.00012 - 0.00005i) - v[4]*v[6]*(0.00013 - 0.0001i) + v[4]*v[7]*(0.0001 - 0.00008i) + v[5]*v[6]*(0.00016 - 0.00001i) - v[5]*v[7]*(0.00013 - 0.00001i) - v[2]^2*(0.00106 + 0.00268i) - v[3]^2*(0.00106 + 0.00268i) - v[4]^2*(0.00235 + 0.00165i) + v[5]^2*(0.0006 - 0.00281i)) + v[2]*(v[1]^3*(0.0000207736 - 0.0000108138i) + v[1]^2*(0.00314984 + 0.00104467i) + v[1]*(0.103133 + 0.0572328i)) - 1.0*v[3]*(v[1]^3*(0.00000428964 + 0.0000230234i) + v[1]^2*(0.00240965 + 0.00190241i) + v[1]*(0.120388 + 0.0527565i)) + v[6]*(- v[1]^3*(0.0000259479 - 0.0000532643i) + v[1]^2*(0.00106845 - 0.00323481i) + v[1]*(0.0949347 - 0.0606467i)) - 1.0*v[5]*(- v[1]^3*(0.0000251947 - 0.000023684i) + v[1]^2*(0.000895386 + 0.00299059i) + v[1]*(0.0879635 + 0.11187i)) + ka3*(v[2]*v[3]^2*(0.00044 - 0.00017i) - v[2]^2*v[3]*(0.00044 - 0.00017i) + v[2]*v[4]^2*(0.00027 - 0.00039i) + v[2]^2*v[4]*(0.00037 - 0.00029i) + v[2]*v[5]^2*(0.00046 + 0.0001i) - v[3]*v[4]^2*(0.00027 - 0.00039i) - v[2]^2*v[5]*(0.00047 - 0.00004i) + v[3]^2*v[4]*(0.00037 - 0.00029i) - v[3]*v[5]^2*(0.00046 + 0.0001i) - v[2]^2*v[6]*0.00001i - v[3]^2*v[5]*(0.00047 - 0.00004i) + v[4]*v[5]^2*(0.00047 - 0.00004i) + v[2]^2*v[7]*0.00001i - v[3]^2*v[6]*0.00001i - v[4]^2*v[5]*(0.00037 - 0.00029i) + v[3]^2*v[7]*0.00001i - v[4]^2*v[6]*(0.00001 + 0.00001i) + v[4]^2*v[7]*(0.00001 + 0.00001i) - v[5]^2*v[6]*0.00001i + v[5]^2*v[7]*0.00001i + v[2]^3*(0.00015 - 0.00006i) - v[3]^3*(0.00015 - 0.00006i) + v[4]^3*(0.00005 - 0.00015i) - v[5]^3*(0.00014 + 0.00008i) - v[2]*v[3]*v[4]*(0.00074 - 0.00059i) + v[2]*v[3]*v[5]*(0.00094 - 0.00008i) + v[2]*v[3]*v[6]*(0.00001 + 0.00003i) - v[2]*v[4]*v[5]*(0.00088 - 0.00035i) - v[2]*v[3]*v[7]*(0.00001 + 0.00002i) - v[2]*v[4]*v[6]*(0.00002 + 0.00002i) + v[3]*v[4]*v[5]*(0.00088 - 0.00035i) + v[2]*v[4]*v[7]*(0.00001 + 0.00002i) + v[2]*v[5]*v[6]*0.00003i + v[3]*v[4]*v[6]*(0.00002 + 0.00002i) - v[2]*v[5]*v[7]*0.00002i - v[3]*v[4]*v[7]*(0.00001 + 0.00002i) - v[3]*v[5]*v[6]*0.00003i + v[3]*v[5]*v[7]*0.00002i + v[4]*v[5]*v[6]*(0.00001 + 0.00003i) - v[4]*v[5]*v[7]*(0.00001 + 0.00002i))

    f[3]=- v[3]*15.1829i + v[5]*(v[1]^3*(0.00000311958 + 0.000034438i) + v[1]^2*(0.00343708 + 0.00109031i) + v[1]*(0.143011 - 0.00225607i)) - 1.0*v[7]*(v[1]^3*(0.000169667 + 0.000348283i) + v[1]^2*(0.009276 + 0.0122307i) + v[1]*(0.154139 + 0.176441i)) - 1.0*ka2*(v[2]*v[3]*(0.00211 - 0.00537i) - v[2]*v[4]*(0.00047 - 0.00574i) + v[2]*v[5]*(0.00357 - 0.00452i) + v[3]*v[4]*(0.00047 - 0.00574i) + v[2]*v[6]*(0.00015 + 0.00006i) - v[3]*v[5]*(0.00357 - 0.00452i) - v[2]*v[7]*(0.00012 + 0.00005i) - v[3]*v[6]*(0.00015 + 0.00006i) + v[4]*v[5]*(0.00211 - 0.00535i) + v[3]*v[7]*(0.00012 + 0.00005i) + v[4]*v[6]*(0.00016 + 0.00001i) - v[4]*v[7]*(0.00013 + 0.00001i) - v[5]*v[6]*(0.00013 + 0.0001i) + v[5]*v[7]*(0.0001 + 0.00008i) - v[2]^2*(0.00106 - 0.00268i) - v[3]^2*(0.00106 - 0.00268i) + v[4]^2*(0.0006 + 0.00281i) - v[5]^2*(0.00235 - 0.00165i)) + v[3]*(v[1]^3*(0.0000207736 + 0.0000108138i) + v[1]^2*(0.00314984 - 0.00104467i) + v[1]*(0.103133 - 0.0572328i)) - 1.0*v[2]*(v[1]^3*(0.00000428964 - 0.0000230234i) + v[1]^2*(0.00240965 - 0.00190241i) + v[1]*(0.120388 - 0.0527565i)) + v[6]*(- v[1]^3*(0.0000259479 + 0.0000532643i) + v[1]^2*(0.00106845 + 0.00323481i) + v[1]*(0.0949347 + 0.0606467i)) - 1.0*v[4]*(- v[1]^3*(0.0000251947 + 0.000023684i) + v[1]^2*(0.000895386 - 0.00299059i) + v[1]*(0.0879635 - 0.11187i)) - 1.0*ka3*(v[2]*v[3]^2*(0.00044 + 0.00017i) - v[2]^2*v[3]*(0.00044 + 0.00017i) + v[2]*v[4]^2*(0.00046 - 0.0001i) + v[2]^2*v[4]*(0.00047 + 0.00004i) + v[2]*v[5]^2*(0.00027 + 0.00039i) - v[3]*v[4]^2*(0.00046 - 0.0001i) - v[2]^2*v[5]*(0.00037 + 0.00029i) + v[3]^2*v[4]*(0.00047 + 0.00004i) - v[3]*v[5]^2*(0.00027 + 0.00039i) - v[2]^2*v[6]*0.00001i - v[3]^2*v[5]*(0.00037 + 0.00029i) + v[4]*v[5]^2*(0.00037 + 0.00029i) + v[2]^2*v[7]*0.00001i - v[3]^2*v[6]*0.00001i - v[4]^2*v[5]*(0.00047 + 0.00004i) + v[3]^2*v[7]*0.00001i - v[4]^2*v[6]*0.00001i + v[4]^2*v[7]*0.00001i + v[5]^2*v[6]*(0.00001 - 0.00001i) - v[5]^2*v[7]*(0.00001 - 0.00001i) + v[2]^3*(0.00015 + 0.00006i) - v[3]^3*(0.00015 + 0.00006i) + v[4]^3*(0.00014 - 0.00008i) - v[5]^3*(0.00005 + 0.00015i) - v[2]*v[3]*v[4]*(0.00094 + 0.00008i) + v[2]*v[3]*v[5]*(0.00074 + 0.00059i) - v[2]*v[3]*v[6]*(0.00001 - 0.00003i) - v[2]*v[4]*v[5]*(0.00088 + 0.00035i) + v[2]*v[3]*v[7]*(0.00001 - 0.00002i) - v[2]*v[4]*v[6]*0.00003i + v[3]*v[4]*v[5]*(0.00088 + 0.00035i) + v[2]*v[4]*v[7]*0.00002i - v[2]*v[5]*v[6]*(0.00002 - 0.00002i) + v[3]*v[4]*v[6]*0.00003i + v[2]*v[5]*v[7]*(0.00001 - 0.00002i) - v[3]*v[4]*v[7]*0.00002i + v[3]*v[5]*v[6]*(0.00002 - 0.00002i) - v[3]*v[5]*v[7]*(0.00001 - 0.00002i) - v[4]*v[5]*v[6]*(0.00001 - 0.00003i) + v[4]*v[5]*v[7]*(0.00001 - 0.00002i))

    f[4]=- v[4]*(5.05926 - 16.7804i) + v[7]*(v[1]^3*(0.000397449 - 0.000246334i) + v[1]^2*(0.0177999 - 0.00583005i) + v[1]*(0.280611 - 0.0666978i)) - 1.0*v[4]*(v[1]^3*(0.0000253568 - 0.00003315i) + v[1]^2*(0.00427065 + 0.00114243i) + v[1]*(0.176935 + 0.10919i)) - 1.0*v[5]*(v[1]^3*(0.0000409701 - 0.00000795858i) + v[1]^2*(0.0010276 - 0.00355808i) + v[1]*(0.0149366 - 0.156251i)) - 1.0*ka3*(v[2]*v[3]^2*(0.00006 + 0.00136i) - v[2]^2*v[3]*(0.00006 + 0.00136i) + v[2]*v[4]^2*(0.00079 + 0.0011i) + v[2]^2*v[4]*(0.00045 + 0.00128i) - v[2]*v[5]^2*(0.0007 - 0.00116i) - v[3]*v[4]^2*(0.00079 + 0.0011i) + v[2]^2*v[5]*(0.00034 - 0.00131i) + v[3]^2*v[4]*(0.00045 + 0.00128i) + v[3]*v[5]^2*(0.0007 - 0.00116i) + 0.00004*v[2]^2*v[6] + v[3]^2*v[5]*(0.00034 - 0.00131i) - v[4]*v[5]^2*(0.00034 - 0.00131i) - 0.00003*v[2]^2*v[7] + 0.00004*v[3]^2*v[6] - v[4]^2*v[5]*(0.00044 + 0.00128i) - 0.00003*v[3]^2*v[7] + v[4]^2*v[6]*(0.00003 - 0.00002i) - v[4]^2*v[7]*(0.00002 - 0.00002i) + v[5]^2*v[6]*(0.00003 + 0.00002i) - v[5]^2*v[7]*(0.00003 + 0.00002i) + v[2]^3*(0.00002 + 0.00045i) - v[3]^3*(0.00002 + 0.00045i) + v[4]^3*(0.00036 + 0.00027i) + v[5]^3*(0.00033 - 0.0003i) - v[2]*v[3]*v[4]*(0.00089 + 0.00256i) - v[2]*v[3]*v[5]*(0.00067 - 0.00263i) - 0.00008*v[2]*v[3]*v[6] - v[2]*v[4]*v[5]*(0.00011 + 0.0027i) + 0.00006*v[2]*v[3]*v[7] + v[2]*v[4]*v[6]*(0.00007 - 0.00003i) + v[3]*v[4]*v[5]*(0.00011 + 0.0027i) - v[2]*v[4]*v[7]*(0.00006 - 0.00002i) - v[2]*v[5]*v[6]*(0.00008 + 0.00002i) - v[3]*v[4]*v[6]*(0.00007 - 0.00003i) + v[2]*v[5]*v[7]*(0.00006 + 0.00001i) + v[3]*v[4]*v[7]*(0.00006 - 0.00002i) + v[3]*v[5]*v[6]*(0.00008 + 0.00002i) - v[3]*v[5]*v[7]*(0.00006 + 0.00001i) - 0.00008*v[4]*v[5]*v[6] + 0.00006*v[4]*v[5]*v[7]) - 1.0*v[2]*(v[1]^3*(0.0000281704 + 0.00000233443i) + v[1]^2*(0.00254065 + 0.0031441i) + v[1]*(0.0954459 + 0.135962i)) - 1.0*v[6]*(- v[1]^3*(0.0000607835 - 0.0000376728i) + v[1]^2*(0.00310379 - 0.00271927i) + v[1]*(0.147457 + 0.0138217i)) + ka2*(v[2]^2*(0.00827 - 0.00035i) - v[2]*v[3]*(0.01654 - 0.0007i) + v[2]*v[4]*(0.0156 - 0.00543i) - v[2]*v[5]*(0.016 + 0.0041i) - v[2]*v[6]*(0.00002 + 0.00048i) + v[7]*v[2]*(0.00002 + 0.00036i) + v[3]^2*(0.00827 - 0.00035i) - v[3]*v[4]*(0.0156 - 0.00543i) + v[3]*v[5]*(0.016 + 0.0041i) + v[3]*v[6]*(0.00002 + 0.00048i) - v[7]*v[3]*(0.00002 + 0.00036i) + v[4]^2*(0.00667 - 0.00484i) - v[4]*v[5]*(0.01648 - 0.00069i) - v[4]*v[6]*(0.00016 + 0.00045i) + v[7]*v[4]*(0.00012 + 0.00034i) + v[5]^2*(0.00706 + 0.00426i) - v[5]*v[6]*(0.00012 - 0.00046i) + v[7]*v[5]*(0.00009 - 0.00035i) - 0.00001*v[6]^2 + 0.00001*v[7]*v[6]) + v[3]*(- v[1]^3*(0.0000104436 - 0.0000262669i) + v[1]^2*(0.0011885 + 0.00347674i) + v[1]*(0.0617113 + 0.120387i))

    f[5]=- v[5]*(5.05926 + 16.7804i) + v[7]*(v[1]^3*(0.000397449 + 0.000246334i) + v[1]^2*(0.0177999 + 0.00583005i) + v[1]*(0.280611 + 0.0666978i)) - 1.0*v[5]*(v[1]^3*(0.0000253568 + 0.00003315i) + v[1]^2*(0.00427065 - 0.00114243i) + v[1]*(0.176935 - 0.10919i)) - 1.0*v[4]*(v[1]^3*(0.0000409701 + 0.00000795858i) + v[1]^2*(0.0010276 + 0.00355808i) + v[1]*(0.0149366 + 0.156251i)) - 1.0*ka3*(0.00004*v[2]^2*v[6] + v[2]^2*v[3]*(0.00006 - 0.00136i) + v[2]*v[4]^2*(0.0007 + 0.00116i) + v[2]^2*v[4]*(0.00034 + 0.00131i) - v[2]*v[5]^2*(0.00079 - 0.0011i) - v[3]*v[4]^2*(0.0007 + 0.00116i) + v[2]^2*v[5]*(0.00045 - 0.00128i) + v[3]^2*v[4]*(0.00034 + 0.00131i) + v[3]*v[5]^2*(0.00079 - 0.0011i) - v[2]*v[3]^2*(0.00006 - 0.00136i) + v[3]^2*v[5]*(0.00045 - 0.00128i) - v[4]*v[5]^2*(0.00044 - 0.00128i) - 0.00003*v[2]^2*v[7] + 0.00004*v[3]^2*v[6] - v[4]^2*v[5]*(0.00034 + 0.00131i) - 0.00003*v[3]^2*v[7] + v[4]^2*v[6]*(0.00003 - 0.00002i) - v[4]^2*v[7]*(0.00003 - 0.00002i) + v[5]^2*v[6]*(0.00003 + 0.00002i) - v[5]^2*v[7]*(0.00002 + 0.00002i) - v[2]^3*(0.00002 - 0.00045i) + v[3]^3*(0.00002 - 0.00045i) + v[4]^3*(0.00033 + 0.0003i) + v[5]^3*(0.00036 - 0.00027i) - v[2]*v[3]*v[4]*(0.00067 + 0.00263i) - v[2]*v[3]*v[5]*(0.00089 - 0.00256i) - 0.00008*v[2]*v[3]*v[6] + v[2]*v[4]*v[5]*(0.00011 - 0.0027i) + 0.00006*v[2]*v[3]*v[7] + v[2]*v[4]*v[6]*(0.00008 - 0.00002i) - v[3]*v[4]*v[5]*(0.00011 - 0.0027i) - v[2]*v[4]*v[7]*(0.00006 - 0.00001i) - v[2]*v[5]*v[6]*(0.00007 + 0.00003i) - v[3]*v[4]*v[6]*(0.00008 - 0.00002i) + v[2]*v[5]*v[7]*(0.00006 + 0.00002i) + v[3]*v[4]*v[7]*(0.00006 - 0.00001i) + v[3]*v[5]*v[6]*(0.00007 + 0.00003i) - v[3]*v[5]*v[7]*(0.00006 + 0.00002i) - 0.00008*v[4]*v[5]*v[6] + 0.00006*v[4]*v[5]*v[7]) - 1.0*v[3]*(v[1]^3*(0.0000281704 - 0.00000233443i) + v[1]^2*(0.00254065 - 0.0031441i) + v[1]*(0.0954459 - 0.135962i)) - 1.0*v[6]*(- v[1]^3*(0.0000607835 + 0.0000376728i) + v[1]^2*(0.00310379 + 0.00271927i) + v[1]*(0.147457 - 0.0138217i)) + ka2*(v[2]^2*(0.00827 + 0.00035i) - v[2]*v[3]*(0.01654 + 0.0007i) + v[2]*v[4]*(0.016 - 0.0041i) - v[2]*v[5]*(0.0156 + 0.00543i) + v[2]*v[6]*(0.00002 - 0.00048i) - v[7]*v[2]*(0.00002 - 0.00036i) + v[3]^2*(0.00827 + 0.00035i) - v[3]*v[4]*(0.016 - 0.0041i) + v[3]*v[5]*(0.0156 + 0.00543i) - v[3]*v[6]*(0.00002 - 0.00048i) + v[7]*v[3]*(0.00002 - 0.00036i) + v[4]^2*(0.00706 - 0.00426i) - v[4]*v[5]*(0.01648 + 0.00069i) - v[4]*v[6]*(0.00012 + 0.00046i) + v[7]*v[4]*(0.00009 + 0.00035i) + v[5]^2*(0.00667 + 0.00484i) - v[5]*v[6]*(0.00016 - 0.00045i) + v[7]*v[5]*(0.00012 - 0.00034i) - 0.00001*v[6]^2 + 0.00001*v[7]*v[6]) + v[2]*(- v[1]^3*(0.0000104436 + 0.0000262669i) + v[1]^2*(0.0011885 - 0.00347674i) + v[1]*(0.0617113 - 0.120387i))

    f[6]=- 47.1636*v[6] - 1.0*ka2*(0.00112*v[2]^2 + v[2]*v[4]*(0.00215 - 0.00065i) - v[2]*v[5]*(0.00215 + 0.00065i) - v[3]*v[4]*(0.00215 - 0.00065i) - v[2]*v[6]*0.00006i + v[3]*v[5]*(0.00215 + 0.00065i) + v[2]*v[7]*0.00005i + v[3]*v[6]*0.00006i - 0.00224*v[4]*v[5] - v[3]*v[7]*0.00005i - v[4]*v[6]*(0.00002 + 0.00006i) + v[4]*v[7]*(0.00001 + 0.00005i) - v[5]*v[6]*(0.00002 - 0.00006i) + v[5]*v[7]*(0.00001 - 0.00005i) - 0.00225*v[2]*v[3] + 0.00112*v[3]^2 + v[4]^2*(0.00093 - 0.00062i) + v[5]^2*(0.00093 + 0.00062i)) - 1.0*v[4]*(- v[1]^3*(0.00000809455 - 0.00000307436i) + v[1]^2*(0.00771859 - 0.00375081i) + v[1]*(0.296768 - 0.00813346i)) - 1.0*v[5]*(- v[1]^3*(0.00000809455 + 0.00000307436i) + v[1]^2*(0.00771859 + 0.00375081i) + v[1]*(0.296768 + 0.00813346i)) + v[7]*(- v[1]^3*(0.0000970094 + 4.23489e-20i) + v[1]^2*(0.0959083 + 2.31888e-18i) + v[1]*(1.99451 + 5.6881e-17i)) - 1.0*v[2]*(- v[1]^3*(0.00000471246 + 0.00000349051i) + v[1]^2*(0.004705 + 0.0027992i) + v[1]*(0.0243829 - 0.0146621i)) - 1.0*v[3]*(- v[1]^3*(0.00000471246 - 0.00000349051i) + v[1]^2*(0.004705 - 0.0027992i) + v[1]*(0.0243829 + 0.0146621i)) - 1.0*v[6]*(v[1]^3*(0.000014836 + 6.47659e-21i) - v[1]^2*(0.0160897 + 9.75407e-19i) + v[1]*(1.96279 + 6.86328e-17i)) + ka3*(v[2]*v[3]^2*0.00018i - v[2]^2*v[3]*0.00018i + v[2]*v[4]^2*(0.0001 + 0.00015i) + v[2]^2*v[4]*(0.00005 + 0.00018i) - v[2]*v[5]^2*(0.0001 - 0.00015i) - v[3]*v[4]^2*(0.0001 + 0.00015i) + v[2]^2*v[5]*(0.00005 - 0.00018i) + v[3]^2*v[4]*(0.00005 + 0.00018i) + v[3]*v[5]^2*(0.0001 - 0.00015i) + 0.00001*v[2]^2*v[6] + v[3]^2*v[5]*(0.00005 - 0.00018i) - v[4]*v[5]^2*(0.00005 - 0.00018i) + 0.00001*v[3]^2*v[6] - v[4]^2*v[5]*(0.00005 + 0.00018i) + v[2]^3*0.00006i - v[3]^3*0.00006i + v[4]^3*(0.00005 + 0.00004i) + v[5]^3*(0.00005 - 0.00004i) - v[2]*v[3]*v[4]*(0.00011 + 0.00035i) - v[2]*v[3]*v[5]*(0.00011 - 0.00035i) - 0.00001*v[2]*v[3]*v[6] - v[2]*v[4]*v[5]*0.00037i + 0.00001*v[2]*v[3]*v[7] + 0.00001*v[2]*v[4]*v[6] + v[3]*v[4]*v[5]*0.00037i - 0.00001*v[2]*v[4]*v[7] - 0.00001*v[2]*v[5]*v[6] - 0.00001*v[3]*v[4]*v[6] + 0.00001*v[2]*v[5]*v[7] + 0.00001*v[3]*v[4]*v[7] + 0.00001*v[3]*v[5]*v[6] - 0.00001*v[3]*v[5]*v[7] - 0.00001*v[4]*v[5]*v[6] + 0.00001*v[4]*v[5]*v[7])

    f[7]=- 7.17103*v[7] + v[4]*(- v[1]^3*(0.00000200279 - 7.60675e-7i) + v[1]^2*(0.0010852 - 0.000614867i) + v[1]*(0.0374856 - 0.00440386i)) + v[5]*(- v[1]^3*(0.00000200279 + 7.60675e-7i) + v[1]^2*(0.0010852 + 0.000614867i) + v[1]*(0.0374856 + 0.00440386i)) - 1.0*v[7]*(- v[1]^3*(0.0000240026 + 1.83924e-20i) + v[1]^2*(0.0138481 - 1.01014e-18i) + v[1]*(0.290099 - 1.68169e-17i)) - 1.0*v[2]*(v[1]^3*(0.00000116598 + 8.63639e-7i) - v[1]^2*(0.000684091 + 0.000337026i) + v[1]*(0.000578732 + 0.00526632i)) - 1.0*v[3]*(v[1]^3*(0.00000116598 - 8.63639e-7i) - v[1]^2*(0.000684091 - 0.000337026i) + v[1]*(0.000578732 - 0.00526632i)) + v[6]*(v[1]^3*(0.00000367081 + 2.81283e-21i) - v[1]^2*(0.00246968 + 1.15122e-19i) + v[1]*(0.288282 - 1.04975e-17i)) - 1.0*ka3*(v[2]*v[3]^2*0.00024i - v[2]^2*v[3]*0.00024i + v[2]*v[4]^2*(0.00013 + 0.0002i) + v[2]^2*v[4]*(0.00007 + 0.00023i) - v[2]*v[5]^2*(0.00013 - 0.0002i) - v[3]*v[4]^2*(0.00013 + 0.0002i) + v[2]^2*v[5]*(0.00007 - 0.00023i) + v[3]^2*v[4]*(0.00007 + 0.00023i) + v[3]*v[5]^2*(0.00013 - 0.0002i) + 0.00001*v[2]^2*v[6] + v[3]^2*v[5]*(0.00007 - 0.00023i) - v[4]*v[5]^2*(0.00007 - 0.00022i) - 0.00001*v[2]^2*v[7] + 0.00001*v[3]^2*v[6] - v[4]^2*v[5]*(0.00007 + 0.00022i) - 0.00001*v[3]^2*v[7] + 0.00001*v[4]^2*v[6] + 0.00001*v[5]^2*v[6] + v[2]^3*0.00008i - v[3]^3*0.00008i + v[4]^3*(0.00006 + 0.00005i) + v[5]^3*(0.00006 - 0.00005i) - v[2]*v[3]*v[4]*(0.00014 + 0.00045i) - v[2]*v[3]*v[5]*(0.00014 - 0.00045i) - 0.00001*v[2]*v[3]*v[6] - v[2]*v[4]*v[5]*0.00047i + 0.00001*v[2]*v[3]*v[7] + 0.00001*v[2]*v[4]*v[6] + v[3]*v[4]*v[5]*0.00047i - 0.00001*v[2]*v[4]*v[7] - 0.00001*v[2]*v[5]*v[6] - 0.00001*v[3]*v[4]*v[6] + 0.00001*v[2]*v[5]*v[7] + 0.00001*v[3]*v[4]*v[7] + 0.00001*v[3]*v[5]*v[6] - 0.00001*v[3]*v[5]*v[7] - 0.00001*v[4]*v[5]*v[6] + 0.00001*v[4]*v[5]*v[7]) + ka2*(0.00144*v[2]^2 + v[2]*v[4]*(0.00275 - 0.00083i) - v[2]*v[5]*(0.00275 + 0.00083i) - v[3]*v[4]*(0.00275 - 0.00083i) - v[2]*v[6]*0.00008i + v[3]*v[5]*(0.00275 + 0.00083i) + v[2]*v[7]*0.00006i + v[3]*v[6]*0.00008i - 0.00287*v[4]*v[5] - v[3]*v[7]*0.00006i - v[4]*v[6]*(0.00002 + 0.00008i) + v[4]*v[7]*(0.00002 + 0.00006i) - v[5]*v[6]*(0.00002 - 0.00008i) + v[5]*v[7]*(0.00002 - 0.00006i) - 0.00288*v[2]*v[3] + 0.00144*v[3]^2 + v[4]^2*(0.0012 - 0.00079i) + v[5]^2*(0.0012 + 0.00079i))

    f[2] = f[2]/ω_h; f[3] = f[3]/ω_h; f[4] = f[4]/ω_h; f[5] = f[5]/ω_h; f[6] = f[6]/ω_h; f[7] = f[7]/ω_h;

    return calculate_normal_form(f, v, norder, cm, m1, m2, m3, n, ω_h, zero_tol)
end

end # module


text:
using Test, SBML, EzXML, MathML, Symbolics
# TODO: ANAND delete all of this or make into readme for copasi feature comparison
fn = "data/01193-sbml-l3v1.xml"
d = sbml_to_sysinfo(fn)
@test isa(d, Dict{String,Union{Nothing,Vector{EzXML.Node}}})
@test !all(isnothing, values(d))

fn = "data/BIOMD0000000346_url.xml"
xml = readxml(fn)
d = sbml_to_sysinfo(fn)
td = make_nice(d)
@test isa(td, Dict)

# Namespacing issue: fbc:ChemicalFormula
@test isempty(collect(skipmissing(td["listOfSpecies"][end].second)))

# start of ReactionNetwork bind

# https://github.com/SciML/CellMLToolkit.jl/blob/dev/src/cellml.jl
# * make find_iv work
fn = "data/BIOMD0000000704_url.xml"
fn = "data/BIOMD0000000469_url.xml" # big boi
d = sbml_to_sysinfo(fn)
td = make_nice(d)

rs = d["listOfRules"]
vars = map(x -> Symbol(x["variable"]), rs)

# todo find iv, make vars => vars(t)
dvs = @. Num(Variable(vars))
D = Differential(t)

eqs = @. parse_node(getproperty(rs, :firstelement))
eqs = reduce(vcat, eqs)

# iv problem, cant ODESystem(deqs)
deqs = D.(vars) .~ eqs

pm = build_map(d["listOfParameters"], "id", "value")
sm = build_map(d["listOfSpecies"], "id", "initialConcentration")
cm = build_map(d["listOfCompartments"], "id", "size")

# do locally
dfs = DataFrame.(values(td))
ddd = Dict(keys(td) .=> dfs)


text:
using DifferentialEquations, Plots

#Reference: https://diego.assencio.com/?index=e5ac36fcb129ce95a61f8e8ce0572dbf

#Solving the double pendulum with a traditional ODE method
#==========================================================#
function doublependulum(du, u, params, t)
    l1 = params[1]
    l2 = params[2]
    m1 = params[3]
    m2 = params[4]
    g  = params[5]

    p1 = u[1]
    p2 = u[2]
    θ1 = u[3]
    θ2 = u[4]

    h1 = p1*p2*sin(θ1-θ2)/(l1*l2*(m1+m2*sin(θ1-θ2)^2))
    h2 = (m2*l2^2*p1^2 + (m1 + m2) * l1^2 * p2^2 - 2*m2*l1*l2*p1*p2*cos(θ1-θ2))/(2 * l1^2 * l2^2 *(m1 + m2*sin(θ1-θ2)^2)^2)

    d_p1 = -(m1 + m2)* g*l1*sin(θ1) - h1 + h2*sin(2*(θ1-θ2))
    d_θ1 = (l2*p1 - l1*p2*cos(θ1-θ2))/(l1^2 * l2*(m1 + m2*sin(θ1 - θ2)^2))
    d_p2 = -m2*g*l2*sin(θ2) + h1 - h2*sin(2*(θ1-θ2))
    d_θ2 = (-m2*l2*p1*cos(θ1 - θ2) + (m1 + m2)*l1*p2) / (m2*l1*l2^2*(m1 + m2*sin(θ1-θ2)^2))

    du .= [d_p1, d_p2, d_θ1, d_θ2]
    return nothing
end

l1 = 1. #length of pendulum1
l2 = 2. #length of pendulum2
m1 = 1. #mass of pendulum1
m2 = 1. #mass of pendulum2
g = 9.81 #gravity

params = [l1, l2, m1, m2, g]
times = (0., 25.)
u0 = [1.,1.,1.,1.]
prob = ODEProblem(doublependulum, u0, times, params)
sol1 = solve(prob, AutoVern7(Rodas5()), dt = .005)

#plot solution
plot(sol1, vars=1, xlim=(0,20), label="Momentum1")
plot!(sol1, vars=2, xlim=(0,20), label="Momentum2")
plot!(sol1, vars=3, xlim=(0,20), label="theta1")
plot!(sol1, vars=4, xlim=(0,20), label="theta2")
#==========================================================#

#Now with HamiltonianProblem()
#==========================================================#
function H(p, θ, params)
    l1 = params[1]
    l2 = params[2]
    m1 = params[3]
    m2 = params[4]
    g  = params[5]

    return  (m2*l2^2*p[1]^2 + (m1+m2)*l1^2*p[2]^2 - 2*m2*l1*l2*p[1]*p[2]*cos(θ[1]-θ[2]) ) /
        (2*m2*l1^2*l2^2*(m1+m2*sin(θ[1]-θ[2])^2)) -
        (m1+m2)*g*l1*cos(θ[1]) - m2*g*l2*cos(θ[2])
end

l1 = 1. #length of pendulum1
l2 = 2. #length of pendulum2
m1 = 1. #mass of pendulum1
m2 = 1. #mass of pendulum2
g = 9.81 #gravity

params = [l1, l2, m1, m2, g]
q0 = [1.0,1.0]
p0 = [1.0,1.0]
times = (0.,25.)
prob = HamiltonianProblem(H, q0, p0, times, params)
sol2 = solve(prob, SofSpa10(), dt = .05)

plot(sol2, vars=1, xlim=(0,20), label="Momentum1")
plot!(sol2, vars=2, xlim=(0,20), label="Momentum2")
plot!(sol2, vars=3, xlim=(0,20), label="theta1")
plot!(sol2, vars=4, xlim=(0,20), label="theta2")
#==========================================================#

#Animation
#==========================================================#
function make_pretty_gif(sol)
    timepoints = sol.t

    x1 = l1*sin.(sol[3,:])
    y1 = -l1*cos.(sol[3,:])
    x2 = x1 + l2*sin.(sol[4,:])
    y2 = y1 - l2*cos.(sol[4,:])

    axis_lim = (l1+l2)*1.2

    anim = Animation()
    for i =1:length(timepoints)
        str = string("Time = ", round(timepoints[i],digits=1), " sec")
        plot([0,x1[i]], [0,y1[i]], size=(400,300), xlim=(-axis_lim,axis_lim), ylim=(-axis_lim,1), markersize = 10, markershape = :circle,label ="",axis = [])
        plot!([x1[i],x2[i]], [y1[i],y2[i]], markersize = 10, markershape = :circle,label ="",title = str, title_location = :left)

        if i > 8 #rainbow trail
            plot!([x2[i-2:i]],   [y2[i-2:i]],  alpha = 0.15, linewidth = 2, color = :red, label=nothing)
            plot!([x2[i-3:i-2]], [y2[i-3:i-2]],alpha = 0.15, linewidth = 2, color = :orange, label=nothing)
            plot!([x2[i-4:i-3]], [y2[i-4:i-3]],alpha = 0.15, linewidth = 2, color = :yellow, label=nothing)
            plot!([x2[i-6:i-4]], [y2[i-6:i-4]],alpha = 0.15, linewidth = 2, color = :green, label=nothing)
            plot!([x2[i-7:i-6]], [y2[i-7:i-6]],alpha = 0.15, linewidth = 2, color = :blue, label=nothing)
            plot!([x2[i-8:i-7]], [y2[i-8:i-7]],alpha = 0.15, linewidth = 2, color = :purple, label=nothing)
        end
        frame(anim)
    end
    gif(anim, fps = 30)
end

make_pretty_gif(sol1)
make_pretty_gif(sol2)


text:
module Rotations3D

using StaticArrays
using LinearAlgebra: norm, rank, svd, Diagonal, tr
using Combinatorics: permutations

export ClebschGordan, Rot3DCoeffs, ri_basis, rpi_basis, R3DC, Rot3DCoeffsEquiv

"""
`ClebschGordan: ` storing precomputed Clebsch-Gordan coefficients; see
`?clebschgordan` for the convention that is use.
"""
struct ClebschGordan{T}
	vals::Dict{Tuple{Int, Int, Int, Int, Int, Int}, T}
end

"""
`Rot3DCoeffs: ` storing recursively precomputed coefficients for a
rotation-invariant basis.
"""
abstract type R3DC{T} end

struct Rot3DCoeffs{T} <: R3DC{T}
   vals::Vector{Dict}
   cg::ClebschGordan{T}
end

# -----------------------------------
# iterating over an m collection
# -----------------------------------

_mvec(::CartesianIndex{0}) = SVector(Int(0))

_mvec(mpre::CartesianIndex) = SVector(Tuple(mpre)..., - sum(Tuple(mpre)))

struct MRange{N, T2}
   ll::SVector{N, Int}
   cartrg::T2
end

Base.length(mr::MRange) = sum(_->1, _mrange(mr.ll))

"""
Given an l-vector `ll` iterate over all combinations of `mm` vectors  of
the same length such that `sum(mm) == 0`
"""
_mrange(ll) = MRange(ll, Iterators.Stateful(
                     CartesianIndices(ntuple(i -> -ll[i]:ll[i], length(ll)-1))))

function Base.iterate(mr::MRange{1}, args...)
   if isempty(mr.cartrg)
      return nothing
   end
   while !isempty(mr.cartrg)
      popfirst!(mr.cartrg)
   end
   return SVector{1, Int}(0), nothing
end

function Base.iterate(mr::MRange, args...)
   while true
      if isempty(mr.cartrg)
         return nothing
      end
      mpre = popfirst!(mr.cartrg)
      if abs(sum(mpre.I)) <= mr.ll[end]
         return _mvec(mpre), nothing
      end
   end
   error("we should never be here")
end

# ----------------------------------------------------------------------
#     ClebschGordan code
# ----------------------------------------------------------------------

cg_conditions(j1,m1, j2,m2, J,M) =
	cg_l_condition(j1, j2, J)   &&
	cg_m_condition(m1, m2, M)   &&
	(abs(m1) <= j1) && (abs(m2) <= j2) && (abs(M) <= J)

cg_l_condition(j1, j2, J) = (abs(j1-j2) <= J <= j1 + j2)

cg_m_condition(m1, m2, M) = (M == m1 + m2)

"""
`clebschgordan(j1, m1, j2, m2, J, M, T=Float64)` :

A reference implementation of Clebsch-Gordon coefficients based on

https://hal.inria.fr/hal-01851097/document
Equation (4-6)

This heavily uses BigInt and BigFloat and should therefore not be employed
for performance critical tasks, but only precomputation.

The ordering of parameters corresponds to the following convention:
```
clebschgordan(j1, m1, j2, m2, J, M) = C_{j1m1j2m2}^{JM}
```
where
```
   D_{m1k1}^{l1} D_{m2k2}^{l2}}
	=
	∑_j  C_{l1m1l2m2}^{j(m1+m2)} C_{l1k1l2k2}^{j2(k1+k2)} D_{(m1+m2)(k1+k2)}^{j}
```
"""
function clebschgordan(j1, m1, j2, m2, J, M, T=Float64)
	if !cg_conditions(j1, m1, j2, m2, J, M)
		return zero(T)
	end

   N = (2*J+1) *
       factorial(big(j1+m1)) * factorial(big(j1-m1)) *
       factorial(big(j2+m2)) * factorial(big(j2-m2)) *
       factorial(big(J+M)) * factorial(big(J-M)) /
       factorial(big( j1+j2-J)) /
       factorial(big( j1-j2+J)) /
       factorial(big(-j1+j2+J)) /
       factorial(big(j1+j2+J+1))

   G = big(0)
   # 0 ≦ k ≦ j1+j2-J
   # 0 ≤ j1-m1-k ≤ j1-j2+J   <=>   j2-J-m1 ≤ k ≤ j1-m1
   # 0 ≤ j2+m2-k ≤ -j1+j2+J  <=>   j1-J+m2 ≤ k ≤ j2+m2
   lb = (0, j2-J-m1, j1-J+m2)
   ub = (j1+j2-J, j1-m1, j2+m2)
   for k in maximum(lb):minimum(ub)
      bk = big(k)
      G += (-1)^k *
           binomial(big( j1+j2-J), big(k)) *
           binomial(big( j1-j2+J), big(j1-m1-k)) *
           binomial(big(-j1+j2+J), big(j2+m2-k))
   end

   return T(sqrt(N) * G)
end

ClebschGordan(T=Float64) =
	ClebschGordan{T}(Dict{Tuple{Int,Int,Int,Int,Int,Int}, T}())

_cg_key(j1, m1, j2, m2, J, M) = (j1, m1, j2, m2, J, M)

function (cg::ClebschGordan{T})(j1, m1, j2, m2, J, M) where {T}
	if !cg_conditions(j1,m1, j2,m2, J,M)
		return zero(T)
	end
	key = _cg_key(j1, m1, j2, m2, J, M)
	if haskey(cg.vals, key)
		return cg.vals[key]
	end
	val = clebschgordan(j1, m1, j2, m2, J, M, T)
	cg.vals[key] = val
	return val
end

# ----------------------------------------------------------------------
#     Rot3DCoeffs code: generalized cg coefficients
#
#  Note: in this section kk is a tuple of m-values, it is not
#        related to the k index in the 1-p basis (or radial basis)
# ----------------------------------------------------------------------

dicttype(N::Integer) = dicttype(Val(N))

dicttype(::Val{N}) where {N} =
   Dict{Tuple{SVector{N,Int}, SVector{N,Int}, SVector{N,Int}}, Float64}

Rot3DCoeffs(T=Float64) = Rot3DCoeffs(Dict[], ClebschGordan(T))

function get_vals(A::R3DC, valN::Val{N}) where {N}
	if length(A.vals) < N
		for n = length(A.vals)+1:N
			push!(A.vals, dicttype(n)())
		end
	end
   return A.vals[N]::dicttype(valN)
end

_key(ll::StaticVector{N}, mm::StaticVector{N}, kk::StaticVector{N}) where {N} =
      (SVector{N, Int}(ll), SVector{N, Int}(mm), SVector{N, Int}(kk))

function (A::Rot3DCoeffs{T})(ll::StaticVector{N},
                             mm::StaticVector{N},
                             kk::StaticVector{N}) where {T, N}
   if       sum(mm) != 0 ||
            sum(kk) != 0 ||
            !all(abs.(mm) .<= ll) ||
            !all(abs.(kk) .<= ll)
      return T(0)
   end
   vals = get_vals(A, Val(N))  # this should infer the type!
   key = _key(ll, mm, kk)
   if haskey(vals, key)
      val  = vals[key]
   else
      val = _compute_val(A, key...)
      vals[key] = val
   end
   return val
end

# the recursion has two steps so we need to define the
# coupling coefficients for N = 1, 2
# TODO: actually this seems false; it is only one recursion step, and a bit
#       or reshuffling should allow us to get rid of the {N = 2} case.

function (A::Rot3DCoeffs{T})(ll::StaticVector{1},
                            mm::StaticVector{1},
                            kk::StaticVector{1}) where {T}
   if ll[1] == mm[1] == kk[1] == 0
      return T(1)
   else
      return T(0)
   end
end

function (A::Rot3DCoeffs{T})(ll::StaticVector{2},
                            mm::StaticVector{2},
                            kk::StaticVector{2}) where {T}
   if ll[1] != ll[2] || sum(mm) != 0 || sum(kk) != 0
      return T(0)
   else
      return T( 8 * pi^2 / (2*ll[1]+1) * (-1)^(mm[1]-kk[1]) )
   end
end

# next comes the recursion step for N ≧ 3

function _compute_val(A::Rot3DCoeffs{T}, ll::StaticVector{N},
                                        mm::StaticVector{N},
                                        kk::StaticVector{N}) where {T, N}
	val = T(0)
   llp = ll[1:N-2]
   mmp = mm[1:N-2]
   kkp = kk[1:N-2]
   for j = abs(ll[N-1]-ll[N]):(ll[N-1]+ll[N])
      if abs(kk[N-1]+kk[N]) > j || abs(mm[N-1]+mm[N]) > j
         continue
      end
		cgk = try
			A.cg(ll[N-1], kk[N-1], ll[N], kk[N], j, kk[N-1]+kk[N])
		catch
			@show (ll[N-1], kk[N-1], ll[N], kk[N], j, kk[N-1]+kk[N])
			T(0)
		end
		cgm = A.cg(ll[N-1], mm[N-1], ll[N], mm[N], j, mm[N-1]+mm[N])
		if cgk * cgm  != 0
			val += cgk * cgm * A( SVector(llp..., j),
								       SVector(mmp..., mm[N-1]+mm[N]),
								       SVector(kkp..., kk[N-1]+kk[N]) )
		end
   end
   return val
end

# ----------------------------------------------------------------------
#   construction of a possible set of generalised CG coefficient;
#   numerically via SVD
# ----------------------------------------------------------------------

function ri_basis(A::R3DC{T}, ll::SVector; ordered=false) where {T}
	CC = compute_Al(A, ll, Val(ordered))
	svdC = svd(CC)
	rk = rank(Diagonal(svdC.S))
	return svdC.U[:, 1:rk]'
end

# unordered
function compute_Al(A::R3DC{T}, ll::SVector, ::Val{false}) where {T}
	len = length(_mrange(ll))
	CC = zeros(T, len, len)
	for (im, mm) in enumerate(_mrange(ll)), (ik, kk) in enumerate(_mrange(ll))
		CC[ik, im] = A(ll, mm, kk)
	end
	return CC
end

# TODO: this could use some documentation

rpi_basis(A::R3DC, zz, nn, ll) =
			rpi_basis(A, SVector(zz...), SVector(nn...), SVector(ll...))

function rpi_basis(A::R3DC,
						 nn::SVector{N, TN},
						 ll::SVector{N, Int}) where {N, TN}
	Uri = ri_basis(A, ll)
	Mri = collect( _mrange(ll) )   # rows...
	G = _gramian(nn, ll, Uri, Mri)
    S = svd(G)
    rk = rank(G; rtol =  1e-7)
	Urpi = S.U[:, 1:rk]'
	return Diagonal(sqrt.(S.S[1:rk])) * Urpi * Uri, Mri
end

function _gramian(nn, ll, Uri, Mri)
   N = length(nn)
   nri = size(Uri, 1)
   @assert size(Uri, 1) == nri
   G = zeros(Complex{Float64}, nri, nri)
   for σ in permutations(1:N)
      if (nn[σ] != nn) || (ll[σ] != ll); continue; end
      for (iU1, mm1) in enumerate(Mri), (iU2, mm2) in enumerate(Mri)
         if mm1[σ] == mm2
            for i1 = 1:nri, i2 = 1:nri
               G[i1, i2] += conj(Uri[i1, iU1]) * Uri[i2, iU2]
            end
         end
      end
   end
   return G
end

## Matthias' code
include("rotations3d-equiv.jl")

## Covariant construction for SphericalVector/Matrix - Liwei
include("rotations3d-spher-vec.jl")

include("rotations3d-spher-mat.jl")

end


text:
type OPCSPActionSpace <: AbstractSpace
    coll::Array{Int} # this may be larger than the actual action space to avoid reallocations - make sure to use len
    len::Int
end
length(as::OPCSPActionSpace) = as.len
OPCSPActionSpace() = OPCSPActionSpace(Array(Int,0),0)
## iteration ##
Base.start(as::OPCSPActionSpace) = 1
Base.done(as::OPCSPActionSpace, state) = state > as.len
Base.next(as::OPCSPActionSpace, state) = (OPCSPAction(as.coll[state]), state+1)
iterator(as::OPCSPActionSpace) = as

function rand!(rng::AbstractRNG, a::OPCSPAction, as::OPCSPActionSpace)
    a.next = rand(rng, as.coll[1:as.len])
    return a
end
function rand(rng::AbstractRNG, as::OPCSPActionSpace, a::OPCSPAction=OPCSPAction())
    a.next = rand(rng, as.coll[1:as.len])
    return a
end

actions(op::Union{OPCSP,OPCSPBeliefMDP}) = actions(op, op.start, delete!(IntSet(1:length(op)), op.start), op.distance_limit)

actions(op::Union{OPCSP,OPCSPBeliefMDP}, b::Union{OPCSPDistribution, OPCSPState, OPCSPBelief}, as::OPCSPActionSpace=OPCSPActionSpace()) = actions(op, b.i, b.open, b.remaining, as)

actions(op::Union{OPCSP,OPCSPBeliefMDP}, b::Tuple{Int,IntSet,Float64,Vector{Int}}) = actions(op, b[1], b[2], b[3])

function actions(op::Union{OPCSP,OPCSPBeliefMDP}, i::Int, open::IntSet, remaining::Float64, as::OPCSPActionSpace=OPCSPActionSpace())
    if length(as.coll) < length(op)
        as.coll = Array(Int, length(op))
        # resize!(as.coll, length(op)) # this can run into a problem in parallel
    end
    f = filter(j->within_range(op,
                               i,
                               op.stop,
                               j,
                               remaining),
               open)
    k = 0
    for a in f
        k += 1
        @inbounds as.coll[k] = a
    end
    as.len = k
    return as
end

### performance may be better with a BitArray of the right length instead of an IntSet


text:
s1. 100x14400 - strength array
2. reshape each layer from distance array
3. radar_counter = 0
  while radar_counter < 100:
    strength[radar_counter] = np.interp(distance, range_bin, data_array[radar_counter])

4. strength array is updated
5. flatten strength array down along the pulse side

max_stars_repo_name: jjimmykang/bwsi-backprojection

max_issues_repo_name: jjimmykang/bwsi-backprojection

max_forks_repo_name: jjimmykang/bwsi-backprojection


text:
using Test
@testset "SummarizedExperiments.jl" begin
using SummarizedExperiments, Documenter
doctest(SummarizedExperiments)
include("class.jl")
end


text:
# coding:utf-8

doc=
"""
  Response Generation Module Clustering for Infomation-Seeking Dialogue System
  2017/02/15 ~
  ・同一話者による連続発話を1発話に統合
  ・1ファイルに必要な情報をすべてマージする
  ・各発話に，IDを振り素性を参照できるようにする
"""

include("switch_depwarn.jl")
switch_depwarn!(false)

using DataFrames
using DataStructures
using PyCall
using ScikitLearn
using JSON
using ScikitLearn: fit!, predict, transform

@sk_import cluster: KMeans
@sk_import metrics: classification_report

# コーパス読み込み
d = open(deserialize, "./data/corpus.dic")

# 対話ID
session_ids = collect(keys(d))

# unigram & bigram 辞書
function word_dict(d, session_ids)
  word_d = Dict(:unigram => [], :bigram => [], :pos => [])
  for ids in session_ids
     for i in  collect(length(d[ids][:unigram]))
       for x in d[ids][:unigram][i]; push!(word_d[:unigram], x) ; end
       for x in d[ids][:bigram][i]; push!(word_d[:bigram], x) ; end
       for x in d[ids][:pos][i]; push!(word_d[:pos], x) ; end
    end
  end
  return Dict(:unigram => counter(word_d[:unigram]),
    :bigram => counter(word_d[:bigram]), :pos => counter(word_d[:pos]))
end

word_d = word_dict(d, session_ids)
d[:unigram] = word_d[:unigram]
d[:bigram] = word_d[:bigram]

# 対話行為ラベル
df = readtable("./data/dialogueactlabel.dat", separator = '\t', header = true)
da_dic = Dict()
da_dic[2000] = ["NULL","NULL"]
for i in collect(1:length(df[1]))
  da = df[:dialogueActLabelId][i]
  label = df[:label][i]
  label_nosub = df[:label][i]
  if length(split(label_nosub, ':')) == 3
    label_nosub = join(split(label_nosub, ':')[1:2],':')
  end
  da_dic[da] = [label_nosub, label]
end

# 追加情報
for sid in session_ids
  normarized_sequence, message_length = [], []
  da_function_label, da_domain_label = [], []
  for (i,cont) in enumerate(d[sid][:unigram])
    push!(normarized_sequence, d[sid][:sequence][i]/length(cont))
    push!(message_length, length(cont))
    push!(da_function_label, da_dic[d[sid][:function_da][i]])
    push!(da_domain_label, da_dic[d[sid][:domain_da][i]])
    d[sid][:speaker][i] = d[sid][:speaker][i] - 1
  end
  d[sid][:norm_sequence] = normarized_sequence
  d[sid][:message_length] = message_length
  d[sid][:function_da_label] = da_function_label
  d[sid][:domain_da_label] = da_domain_label
end

# corpus dictionary
open(io -> serialize(io, d), "./data/annotated_corpus.dic", "w")
JSON.print(open("./data/annotated_corpus.json", "w"), JSON.json(d))


text:
@testset "WeightedSum fit!" begin
    test_fit(WeightedSum{Float64}, x, w, sum, (x, w) -> sum(x .* w))
    test_fit(WeightedSum{Float32}, x, w, sum, (x, w) -> sum(x .* w))

    test_fit(WeightedSum{Float64}, xmis, wmis, sum, (x, w) -> sum(skipmissing(x .* w)))

    s = sum(broadcast(*, x, w))

    szip = sum(fit!(WeightedSum(), zip(x, w)))
    sval = sum(fit!(WeightedSum(), x, w))

    s2val = sum(fit!(WeightedSum(Float32), x, w))
    s2zip = sum(fit!(WeightedSum(Float32), zip(x, w)))

    @test s ≈ szip
    @test s ≈ sval

    @test s ≈ s2zip
    @test s ≈ s2val
    @test typeof(s2val) == Float32
    @test typeof(s2zip) == Float32
end

@testset "WeightedSum merge!" begin
    s = sum(broadcast(*, x, w))

    smap = map(x, w) do xi, wi
        fit!(WeightedSum(), xi, wi)
    end;

    rs = reduce(merge!, smap) |> sum
    rs2 = merge!(
        fit!(WeightedSum(), x[1:end ÷ 2], w[1:end ÷ 2]),
        fit!(WeightedSum(), x[end ÷ 2 + 1:end], w[end ÷ 2 + 1:end])
    ) |> sum

    # Float32
    smap2 = map(x, w) do xi, wi
        fit!(WeightedSum(), xi, wi)
    end;

    rs_32 = reduce(merge!, smap2, init = WeightedSum(Float32)) |> sum
    rs2_32 = merge!(
        fit!(WeightedSum(Float32), x[1:end ÷ 2], w[1:end ÷ 2]),
        fit!(WeightedSum(), x[end ÷ 2 + 1:end], w[end ÷ 2 + 1:end])
    ) |> sum

    @test rs ≈ s
    @test rs2 ≈ s

    @test rs_32 ≈ s
    @test rs2_32 ≈ s

    @test typeof(rs) == Float64
    @test typeof(rs2) == Float64
    @test typeof(rs_32) == Float32
    @test typeof(rs2_32) == Float32
end

@testset "WeightedSum constructor" begin
    @test WeightedSum{Float64}() == WeightedSum()
    @test WeightedSum{Float32}() == WeightedSum(Float32)
    @test WeightedSum() == WeightedSum(0.0, 0.0, 0)
end