text:
using ProgressMeter
using Dates

function main()

    # Start of time measurement
    start_t = now()

    out_dir  = "result"         # output dir
    PARAMDAT = "PARAMDAT.json"  # directory
    fwrite   = "write"          # write file

    nval = 5         # number of conserved variables
    Rd   = 287.0     # gas constant of air, J/(kg K)
    R    = 8.314     # gas constant, J/(K mol)
    icell = 2        # imaginary cell

    # read grids and parameter
    xmax, ymax, zmax, nodes, vecAx, vecAy, vecAz = read_allgrid()
    out_file_front, out_ext, restartnum, restart_file, init_small, norm_ok,
    time_integ, nt, dt, every_outnum, in_nt, dtau, cfl, ad_scheme,
    init_rho, init_u, init_v, init_w, init_p, init_T, specific_heat_ratio, Rd, bdcon = input_para(PARAMDAT)

    # number of cells
    cellxmax = xmax - 1
    cellymax = ymax - 1
    cellzmax = zmax - 1

    # allocation
    Qbase, Qbase_ave, volume, cellcenter, wally, yplus, Qcon, Qcon_hat, mu, mut, mut_bd, lambda,
    E_adv_hat, F_adv_hat, G_adv_hat, E_vis_hat, F_vis_hat, G_vis_hat, RHS, QbaseU, QbaseD, QbaseL, QbaseR, QbaseF, QbaseB,
    QconU, QconD, QconL, QconR, QconF, QconB = common_allocation(cellxmax, cellymax, cellzmax, nval)

    # set initial condition
    Qbase, restartnum = set_initQbase(Qbase, cellxmax, cellymax, cellzmax, restart_file,
                                        init_rho, init_u, init_v, init_w, init_p, init_T,
                                        specific_heat_ratio, out_file_front, out_ext, out_dir, restartnum, Rd, nval, icell)

    # set initial condition for imaginary cell
    Qbase    = set_boundary(Qbase, cellxmax, cellymax, cellzmax, vecAx, vecAy, vecAz, bdcon, Rd, specific_heat_ratio, nval, icell)

    # set volume
    volume = set_volume(nodes, cellxmax, cellymax, cellzmax, vecAx, vecAy, vecAz, volume)
    cellcenter = set_cellcenter(cellcenter, nodes, cellxmax, cellymax, cellzmax)
    reset_write(fwrite)

    #println(volume[3,3,3])

    # wally
    wally, swith_wall = set_wally(nodes, bdcon, wally, cellcenter, cellxmax, cellymax, cellzmax, icell)

    # AUSM+up
    Minf = 0.0
    if ad_scheme == 2
        Minf = set_Minf(bdcon, specific_heat_ratio, Rd, nval)
    end

    # write number of threads
    print("threads num : ")
    println(Threads.nthreads())

    # check boundary condition
    check_bd(bdcon)

    # main loop
    loop_ite = 0
    if time_integ == "1"
        # exlicit scheme
        prog = Progress(nt,1)
        @time for t in 1:nt
            next!(prog)

            loop_ite += 1

            # step number
            stepnum = t + restartnum

            # set conserved variables in the general coordinate system
            Qbase    = set_boundary(Qbase, cellxmax, cellymax, cellzmax, vecAx, vecAy, vecAz,
                                    bdcon, Rd, specific_heat_ratio, nval, icell)
            Qcon     = base_to_conservative(Qbase, Qcon, cellxmax, cellymax, cellzmax, specific_heat_ratio)
            Qcon_hat = setup_Qcon_hat(Qcon, Qcon_hat, cellxmax, cellymax, cellzmax, volume, nval)

            # muscl
            QbaseU, QbaseD, QbaseL, QbaseR, QbaseF, QbaseB,
            QconU, QconD, QconL, QconR, QconF, QconB = muscl(Qbase, QbaseU, QbaseD, QbaseL, QbaseR, QbaseF, QbaseB,
                                                            QconU, QconD, QconL, QconR, QconF, QconB,
                                                            cellxmax, cellymax, cellzmax, nval, icell, specific_heat_ratio, volume, nodes)

            # set viscosity and thermal Conductivity
            mu     = set_mu(mu, Qbase, cellxmax, cellymax, cellzmax, specific_heat_ratio, Rd)
            lambda = set_lambda(lambda, Qbase, cellxmax, cellymax, cellzmax, mu, specific_heat_ratio, Rd)

            # yplus
            yplus = cal_yplus(yplus, Qbase, wally, swith_wall, mu, cellxmax, cellymax, cellzmax, vecAx, vecAy, vecAz, volume, icell)

            # advection_term
            E_adv_hat, F_adv_hat, G_adv_hat = AUSM(E_adv_hat, F_adv_hat, G_adv_hat, QbaseF, QbaseB, QbaseU, QbaseD, QbaseL, QbaseR,
                                                QconF, QconB, QconU, QconD, QconL, QconR, cellxmax, cellymax, cellzmax,
                                                vecAx, vecAy, vecAz, specific_heat_ratio, volume, nval, Minf, ad_scheme, icell)

            # viscos_term
            E_vis_hat, F_vis_hat, G_vis_hat, mut = central_diff(E_vis_hat, F_vis_hat, G_vis_hat, QbaseU, QbaseD, QbaseL, QbaseR, QbaseF, QbaseB,
                                                            QconU, QconD, QconL, QconR, QconF, QconB, cellxmax, cellymax, cellzmax, mu, mut, mut_bd, lambda,
                                                            vecAx, vecAy, vecAz, specific_heat_ratio, volume, Rd, nval, yplus, swith_wall, icell)

            # RHS
            RHS = setup_RHS(RHS, cellxmax, cellymax, cellzmax, E_adv_hat, F_adv_hat, G_adv_hat,
                            E_vis_hat, F_vis_hat, G_vis_hat, nval, volume, icell)

            # time integral
            Qcon_hat = time_integration_explicit(dt, Qcon_hat, RHS, cellxmax, cellymax, cellzmax, nval, icell)

            # calculate primitive variables
            Qcon  = Qhat_to_Q(Qcon, Qcon_hat, cellxmax, cellymax, cellzmax, volume, nval)
            Qbase = conservative_to_base(Qbase, Qcon, cellxmax, cellymax, cellzmax, specific_heat_ratio)
            Qbase_ave = cal_Qave(Qbase, Qbase_ave, cellxmax, cellymax, cellzmax, nval, stepnum)

            # output
            if round(stepnum) % every_outnum == 0
                println("\n")
                println("nt_______________________________"*string(round(stepnum)))
                output_result(stepnum, Qbase, cellxmax, cellymax, cellzmax, specific_heat_ratio,
                                out_file_front, out_ext, out_dir, Rd, nval, icell)
                output_ave(Qbase_ave, cellxmax, cellymax, cellzmax, out_file_front, out_ext, out_dir, Rd, nval, icell)
            end

            # Find out if the results were divergent
            check_divrege(Qbase, cellxmax, cellymax, cellzmax, Rd, fwrite, icell)
        end
    elseif time_integ == "2"
        # implicit

        # allocation for implicit
        Qbasen, Qconn, Qconn_hat, Qbasem, dtau, lambda_facex, lambda_facey, lambda_facez, dx, dy, dz,
        A_adv_hat_m, A_adv_hat_p, B_adv_hat_m, B_adv_hat_p, C_adv_hat_p, C_adv_hat_m, A_beta_shig, B_beta_shig, C_beta_shig,
        jalphaP, jbetaP, jgammaP, delta_Q, delta_Q_temp, D, Lx, Ly, Lz, Ux, Uy, Uz, LdQ, UdQ, RHS_temp, res,
        norm2, I = allocation_implicit(cellxmax, cellymax, cellzmax, nval)

        dx, dy, dz = set_dx_lts(dx, dy, dz, nodes, cellxmax, cellymax, cellzmax, icell)

        prog = Progress(nt,1)
        @time for t in 1:nt
            next!(prog)

            # step number
            loop_ite += 1
            stepnum = t + restartnum

            # write physical time
            output_physicaltime(fwrite, t, dt)

            # copy
            for l in 1:nval
                for k in 1:cellzmax
                    for j in 1:cellymax
                        for i in 1:cellxmax
                            Qbasen[i,j,k,l] = Qbase[i,j,k,l]
                            Qbasem[i,j,k,l] = Qbase[i,j,k,l]
                        end
                    end
                end
            end

            # set conserved variables in the general coordinate system for inner iteration
            Qconn     = base_to_conservative(Qbasen, Qconn, cellxmax, cellymax, cellzmax, specific_heat_ratio)
            Qconn_hat = setup_Qcon_hat(Qconn, Qconn_hat, cellxmax, cellymax, cellzmax, volume, nval)

            # start inner iteration
            for tau in 1:in_nt
                # set conserved variables in the general coordinate system
                Qbasem    = set_boundary(Qbasem, cellxmax, cellymax, cellzmax, vecAx, vecAy, vecAz,
                bdcon, Rd, specific_heat_ratio, nval, icell)
                Qcon     = base_to_conservative(Qbasem, Qcon, cellxmax, cellymax, cellzmax, specific_heat_ratio)
                Qcon_hat = setup_Qcon_hat(Qcon, Qcon_hat, cellxmax, cellymax, cellzmax, volume, nval)

                #println("Qbasem")
                #println(Qbasem[4,4,:,4])

                # muscl
                QbaseU, QbaseD, QbaseL, QbaseR, QbaseF, QbaseB,
                QconU, QconD, QconL, QconR, QconF, QconB = muscl(Qbasem, QbaseU, QbaseD, QbaseL, QbaseR, QbaseF, QbaseB,
                                                        QconU, QconD, QconL, QconR, QconF, QconB,
                                                        cellxmax, cellymax, cellzmax, nval, icell, specific_heat_ratio, volume, nodes)

                # set viscosity and thermal Conductivity
                mu     = set_mu(mu, Qbasem, cellxmax, cellymax, cellzmax, specific_heat_ratio, Rd)
                lambda = set_lambda(lambda, Qbasem, cellxmax, cellymax, cellzmax, mu, specific_heat_ratio, Rd)

                # yplus
                yplus = cal_yplus(yplus, Qbasem, wally, swith_wall, mu, cellxmax, cellymax, cellzmax, vecAx, vecAy, vecAz, volume, icell)

                # advection_term
                E_adv_hat, F_adv_hat, G_adv_hat = AUSM(E_adv_hat, F_adv_hat, G_adv_hat, QbaseF, QbaseB, QbaseU, QbaseD, QbaseL, QbaseR,
                                                    QconF, QconB, QconU, QconD, QconL, QconR, cellxmax, cellymax, cellzmax,
                                                    vecAx, vecAy, vecAz, specific_heat_ratio, volume, nval, Minf, ad_scheme, icell)

                # viscos_term
                E_vis_hat, F_vis_hat, G_vis_hat, mut = central_diff(E_vis_hat, F_vis_hat, G_vis_hat, QbaseU, QbaseD, QbaseL, QbaseR, QbaseF, QbaseB,
                                                                QconU, QconD, QconL, QconR, QconF, QconB, cellxmax, cellymax, cellzmax, mu, mut, mut_bd, lambda,
                                                                vecAx, vecAy, vecAz, specific_heat_ratio, volume, Rd, nval, yplus, swith_wall, icell)

                # RHS
                RHS = setup_RHS(RHS, cellxmax, cellymax, cellzmax, E_adv_hat, F_adv_hat, G_adv_hat,
                                E_vis_hat, F_vis_hat, G_vis_hat, nval, volume, icell)

                #println(RHS[10,10,:,:])

                # 粘性の更新
                for k in 1:cellzmax
                    for j in 1:cellymax
                        for i in 1:cellxmax
                            mu[i,j,k] = mu[i,j,k] + mut[i,j,k]
                        end
                    end
                end

                # set inner time step by local time stepping
                dtau   = set_lts(dtau, lambda_facex, lambda_facey, lambda_facez, Qbasem, cellxmax, cellymax, cellzmax,
                                    mu, dx, dy, dz, vecAx, vecAy, vecAz, volume, specific_heat_ratio, cfl, icell)

                # lusgs_advection_term
                A_adv_hat_p, A_adv_hat_m, B_adv_hat_p, B_adv_hat_m, C_adv_hat_p, C_adv_hat_m,
                A_beta_shig, B_beta_shig, C_beta_shig = one_wave(A_adv_hat_m, A_adv_hat_p, B_adv_hat_m, B_adv_hat_p, C_adv_hat_m, C_adv_hat_p,
                                                                A_beta_shig, B_beta_shig, C_beta_shig, I, Qbase, Qcon, cellxmax, cellymax, cellzmax,
                                                                vecAx, vecAy, vecAz, specific_heat_ratio, volume, nval)
                # lusgs_viscos_term
                jalphaP, jbetaP, jgammaP = central_diff_jacobian(jalphaP, jbetaP, jgammaP, Qbase, Qcon, cellxmax, cellymax, cellzmax, mu, lambda,
                                                        vecAx, vecAy, vecAz, specific_heat_ratio, volume, nval)

                #println(RHS[:,:,6,:])
                #println(A_adv_hat_m[4,2,:,:,:])
                #println(dtau[4,4,:])
                #=
                println(E_adv_hat[10,10,5,:])
                println(F_adv_hat[10,10,5,:])
                println(G_adv_hat[10,10,5,:])
                =#

                #=

                for l in 1:nval
                    for k in 1+icell:cellzmax-icell
                        for j in 1+icell:cellymax-icell
                            for i in 1+icell:cellxmax-icell
                                #if isequal(RHS[i,j,k,l], NaN) == true
                                if isequal(F_adv_hat[i,j,k,l], NaN) == true
                                    println("  ")
                                    println(i)
                                    println(j)
                                    println(k)
                                end
                            end
                        end
                    end
                end
                =#

                # LUSGS
                ite = 0
                while true
                    # copy delta_Q
                    for l in 1:nval
                        for k in 1:cellzmax
                            for j in 1:cellymax
                                for i in 1:cellxmax
                                    delta_Q_temp[i,j,k,l] = delta_Q[i,j,k,l]
                                end
                            end
                        end
                    end
                    #println(delta_Q[4,4,:,:])

                    # Reversing the left-hand side by lusgs
                    delta_Q = lusgs(D, Lx, Ly, Lz, Ux, Uy, Uz, LdQ, UdQ, RHS_temp, I, dt, dtau, Qcon_hat, Qconn_hat, delta_Q,
                    A_adv_hat_p, A_adv_hat_m, B_adv_hat_p, B_adv_hat_m, C_adv_hat_p, C_adv_hat_m, A_beta_shig, B_beta_shig, C_beta_shig,
                    jalphaP, jbetaP, jgammaP, RHS, cellxmax, cellymax, cellzmax, volume, nval, icell)

                    #println(delta_Q[:,:,:,:])

                    # cal Residuals by norm-2
                    res   = set_res(res, delta_Q, delta_Q_temp, cellxmax, cellymax, cellzmax, nval, icell)
                    norm2 = check_converge(res, RHS, cellxmax, cellymax, cellzmax, init_small, nval, icell)

                    # Find out if the results converged
                    if norm2[1] < norm_ok && norm2[2] < norm_ok && norm2[3] < norm_ok && norm2[4] < norm_ok && norm2[5] < norm_ok
                        break
                    end
                    if ite % 100 ==0
                        println(" now cal norm2 ")
                        println(norm2)
                        if isequal(norm2[1], NaN) == true
                            println("  ")
                            println(" norm2 = NaN ")
                            println(" stop cal ")
                            throw(UndefVarError(:x))
                        end
                    end

                    ite += 1
                end

                # output inner time
                output_innertime(fwrite, tau, norm2, nval)

                # Updating
                for l in 1:nval
                    for k in 1+icell:cellzmax-icell
                        for j in 1+icell:cellymax-icell
                            for i in 1+icell:cellxmax-icell
                                Qcon_hat[i,j,k,l] = Qcon_hat[i,j,k,l] + delta_Q[i,j,k,l]
                            end
                        end
                    end
                end
                #println(Qbasem[4,4,:,4])

                # calculate primitive variables
                Qcon  = Qhat_to_Q(Qcon, Qcon_hat, cellxmax, cellymax, cellzmax, volume, nval)
                Qbasem = conservative_to_base(Qbasem, Qcon, cellxmax, cellymax, cellzmax, specific_heat_ratio)

                # println(delta_Q[4,4,:,4])
                #println(Qbasem[4,4,:,4])

                # Find out if the results were divergent
                check_divrege(Qbasem, cellxmax, cellymax, cellzmax, Rd, fwrite, icell)
            end
            # End of the inner iteration

            # Updating in physical time
            for l in 1:nval
                for k in 1:cellzmax
                    for j in 1:cellymax
                        for i in 1:cellxmax
                            Qbase[i,j,k,l] = Qbasem[i,j,k,l]
                        end
                    end
                end
            end

            Qbase_ave = cal_Qave(Qbase_ave, Qbase, cellxmax, cellymax, cellzmax, nval, stepnum)

            # output
            if round(stepnum) % every_outnum == 0
                println("\n")
                println("nt_______________________________"*string(round(stepnum)))
                output_result(stepnum, Qbase, cellxmax, cellymax, cellzmax, specific_heat_ratio,
                                out_file_front, out_ext, out_dir, Rd, nval, icell)
                output_ave(Qbase_ave, cellxmax, cellymax, cellzmax, out_file_front, out_ext, out_dir, Rd, nval, icell)
            end

            # Find out if the results were divergent
            check_divrege(Qbase, cellxmax, cellymax, cellzmax, Rd, fwrite, icell)
        end
    end

    # end of time measurement
    end_t = now()

    # output of calculation time
    output_fin(fwrite, start_t, end_t, nt, dt, in_nt, cellxmax, cellymax, cellzmax)
end

# -- main --
main()
#throw(UndefVarError(:x))


text:
###########################  FinancialTimeSeries & TickData

# single row
function getindex(ft::Union(FinancialTimeSeries, TickData), n::Int)
    isa(ft, FinancialTimeSeries) ?
    FinancialTimeSeries(ft.timestamp[n], ft.values[n,:], ft.colnames, ft.instrument) :
    TickData(ft.timestamp[n], ft.values[n,:], ft.colnames, ft.instrument)
end

# range of rows
function getindex(ft::Union(FinancialTimeSeries, TickData), r::Range1{Int})
    isa(ft, FinancialTimeSeries) ?
    FinancialTimeSeries(ft.timestamp[r], ft.values[r,:], ft.colnames, ft.instrument) :
    TickData(ft.timestamp[r], ft.values[r,:], ft.colnames, ft.instrument)
end

# array of rows
function getindex(ft::Union(FinancialTimeSeries, TickData), a::Array{Int})
    isa(ft, FinancialTimeSeries) ?
    FinancialTimeSeries(ft.timestamp[a], ft.values[a,:], ft.colnames, ft.instrument) :
    TickData(ft.timestamp[a], ft.values[a,:], ft.colnames, ft.instrument)
end

# single column fty name
function getindex(ft::Union(FinancialTimeSeries, TickData), s::ASCIIString)
    n = findfirst(ft.colnames, s)
    isa(ft, FinancialTimeSeries) ?
    FinancialTimeSeries(ft.timestamp, ft.values[:, n], ASCIIString[s], ft.instrument) :
    TickData(ft.timestamp, ft.values[:, n], ASCIIString[s], ft.instrument)
end

# array of columns fty name
function getindex(ft::Union(FinancialTimeSeries, TickData), args::ASCIIString...)
    ns = [findfirst(ft.colnames, a) for a in args]
    isa(ft, FinancialTimeSeries) ?
    FinancialTimeSeries(ft.timestamp, ft.values[:,ns], ASCIIString[a for a in args], ft.instrument) :
    TickData(ft.timestamp, ft.values[:,ns], ASCIIString[a for a in args], ft.instrument)
end

# single date
function getindex(ft::Union(FinancialTimeSeries, TickData), d::DateTime{ISOCalendar,UTC})
   for i in 1:length(ft)
     if [d] == ft[i].timestamp
       return ft[i]
     else
       nothing
     end
   end
 end

# range of dates
function getindex(ft::Union(FinancialTimeSeries, TickData), dates::Array{DateTime{ISOCalendar,UTC},1})
  counter = Int[]
  for i in 1:length(dates)
    if findfirst(ft.timestamp, dates[i]) != 0
      push!(counter, findfirst(ft.timestamp, dates[i]))
    end
  end
  ft[counter]
end

#### DOESN"T WORK #########

function getindex(ft::Union(FinancialTimeSeries, TickData), r::DateTimeRange{ISOCalendar,UTC})
    ft[[r]]
end

# day of week
# getindex{T,N}(ft::Union(FinancialTimeSeries, TickData){T,N}, d::DAYOFWEEK) = ft[dayofweek(ft.timestamp) .== d]

ext: jl


text:
#### piecewise linear warping function.

mutable struct Piecewise2DLineType{T}
    xs::Vector{T}
    ys::Vector{T}
    ms::Vector{T}
    bs::Vector{T}
    len_s::Vector{T}
    len_z::Vector{T}
end

"""
    getpiecewiselines(z_st::Vector{T}, z_fin::Vector{T}, c::T;
        lb::T = -one(T), ub::T = one(T))

Computes the parameters of a piece-wise linear function that contain the intervals in z_st and z_fin in the range, over the range in c.

f = xx->MonotoneMaps.evalpiecewise2Dlinearfunc(xx, xs, ys, ms, bs)*scale
maps [lb, ub] to [lb, ub]*scale.

# Examples
```jldoctest
julia> c = 1.8
1.8

julia> z_st = [ -0.12; ]
1-element Vector{Float64}:
 -0.12

julia> z_fin = [ 0.76; ]
1-element Vector{Float64}:
 0.76

julia> xs, ys, ms, bs, len_s, len_z = MonotoneMaps.getpiecewiselines(z_st, z_fin, c)
([-1.0, -0.842857142857143, 0.9571428571428571, 1.0], [-1.0, -0.12, 0.76, 1.0], [5.600000000000001, 0.4888888888888889, 5.600000000000001], [4.600000000000001, 0.29206349206349214, -4.600000000000001], [1.8], [0.88])

See piece-wise_linear.jl in the /examples folder for an example.
```
"""
function getpiecewiselines(z_st::Vector{T}, z_fin::Vector{T}, c::T;
    lb::T = -one(T), ub::T = one(T)) where T

    N_zones = length(z_st)

    @assert zero(T) < c < (ub-lb)

    ## set up s.
    @assert N_zones == length(z_fin)
    len_z = collect( z_fin[i]-z_st[i] for i = 1:N_zones )

    # sanity check.
    @assert all( len_z .> zero(T) )

    Z = sum(len_z)
    len_s = collect( c * len_z[i]/Z for i = 1:N_zones )

    ## get slope.
    m_zones, m0 = getslope(len_z, len_s, ub-lb)

    xs, ys, ms, bs = buildtransitionpts(lb, ub, m0, m_zones, z_st, z_fin)

    # numerical sensitive when z_st or z_fin contains lb or ub as an element.
    clamp!(xs, lb, ub)
    clamp!(ys, lb, ub)

    return xs, ys, ms, bs, len_s, len_z
end

function evalpiecewise2Dlinearfunc(x, A::Piecewise2DLineType{T})::T where T

    return evalpiecewise2Dlinearfunc(x, A.xs, A.ys, A.ms, A.bs)
end

function evalpiecewise2Dlinearfunc(x_inp::T, A::Piecewise2DLineType{T}, scale::T)::T where T
    return evalpiecewise2Dlinearfunc(x, A.xs, A.ys, A.ms, A.bs, scale)
end

"""
    evalpiecewise2Dlinearfunc(x_inp::T,
        xs::Vector{T}, ys::Vector{T}, ms::Vector{T}, bs::Vector{T}, scale::T)::T where T

Evaluates a piece-wise linear function.

# Select arguments
Two-element 1-D array for the following. The first entry relates to the a parameter, and the second relates to the b parameter.
- `x_inp::T`: input argument to function.
- `xs::Vector{T}`: x-coordinates for the line endpoints. Length L + 1.
- `ys::Vector{T}`: y-coordinates for the line endpoints. Length L + 1.
- `ms::Vector{T}`: slope for the lines. Length L.
- `bs::Vector{T}`: y-intercepts for the lines. Length L.
- `scale::T`: Re-scaling parameter to `x` so that it is in [-1,1], before applying the piece-wise linear function on the domain [-1,1].

***
evalpiecewise2Dlinearfunc(x_inp::T,
    xs::Vector{T}, ys::Vector{T}, ms::Vector{T}, bs::Vector{T})::T where T

In the case when the `scale` argument is absent, it is taken to be 1.

See piece-wise_linear.jl in the /examples folder for an example.
"""
function evalpiecewise2Dlinearfunc(x_inp::T,
    xs::Vector{T}, ys::Vector{T}, ms::Vector{T}, bs::Vector{T}, scale::T)::T where T

    x = clamp(x_inp/scale, -one(T), one(T))
    return evalpiecewise2Dlinearfunc(x, xs, ys, ms, bs)*scale
end

function evalpiecewise2Dlinearfunc(x::T,
    xs::Vector{T}, ys::Vector{T}, ms::Vector{T}, bs::Vector{T})::T where T

    @assert length(xs)-1 == length(ys)-1 == length(ms) == length(bs)

    ind = findfirst(xx->xx>x, xs)

    if typeof(ind) == Int

        return ms[ind-1]*x + bs[ind-1]
    end

    if x < xs[1]
        return ys[1]
    end

    return ys[end]
end

function evalinversepiecewise2Dlinearfunc(y, A::Piecewise2DLineType{T})::T where T

    return evalinversepiecewise2Dlinearfunc(y, A.xs, A.ys, A.ms, A.bs)
end

function evalinversepiecewise2Dlinearfunc(x_inp::T, A::Piecewise2DLineType{T}, scale::T)::T where T
    return evalinversepiecewise2Dlinearfunc(x, A.xs, A.ys, A.ms, A.bs, scale)
end

"""
    evalinversepiecewise2Dlinearfunc(x_inp::T,
        xs::Vector{T}, ys::Vector{T}, ms::Vector{T}, bs::Vector{T}, scale::T)::T where T

Evaluates the inverse of a (forward) piece-wise linear function, given the parameters of the forward function.

# Select arguments
Two-element 1-D array for the following. The first entry relates to the a parameter, and the second relates to the b parameter.
- `x_inp::T`: input argument to function.
- `xs::Vector{T}`: x-coordinates for the line endpoints of the forward function. Length L + 1.
- `ys::Vector{T}`: y-coordinates for the line endpoints of the forward function. Length L + 1.
- `ms::Vector{T}`: slope for the lines of the forward function. Length L.
- `bs::Vector{T}`: y-intercepts for the lines of the forward function. Length L.
- `scale::T`: Re-scaling parameter to `x` so that it is in [-1,1], before applying the piece-wise linear inverse function on the domain [-1,1].

***
evalinversepiecewise2Dlinearfunc(x_inp::T,
    xs::Vector{T}, ys::Vector{T}, ms::Vector{T}, bs::Vector{T})::T where T

In the case when the `scale` argument is absent, it is taken to be 1.

See piece-wise_linear.jl in the /examples folder for an example.
"""
function evalinversepiecewise2Dlinearfunc(x_inp::T,
    xs::Vector{T}, ys::Vector{T}, ms::Vector{T}, bs::Vector{T}, scale::T)::T where T

    x = clamp(x_inp/scale, -one(T), one(T))
    return evalinversepiecewise2Dlinearfunc(x, xs, ys, ms, bs)*scale
end

function evalinversepiecewise2Dlinearfunc(y::T,
    xs::Vector{T}, ys::Vector{T}, ms::Vector{T}, bs::Vector{T})::T where T

    @assert length(xs)-1 == length(ys)-1 == length(ms) == length(bs)

    ind = findfirst(yy->yy>y, ys)

    if typeof(ind) == Int && ind > 1

        return findx2Dline(ms[ind-1], bs[ind-1], y)
    end

    if y < ys[1]
        return xs[1]
    end

    return xs[end]
end

## build transition points.
function buildtransitionpts(lb::T, ub::T, m0::T, m_zones::Vector{T}, z_st::Vector{T}, z_fin::Vector{T}) where T
    N_zones = length(m_zones)
    @assert length(z_st) == length(z_fin) == N_zones

    # set up y coodinates of the transition points.
    ys::Vector{T} = [lb;]
    xs::Vector{T} = [lb;]

    ms = Vector{T}(undef, 0)
    bs = Vector{T}(undef, 0)

    for i = 1:N_zones
        push!(ms, m0)
        push!(bs, findyintercept(ms[end], xs[end], ys[end]))

        push!(ys, z_st[i])
        push!(ms, m_zones[i])
        push!(xs, findx2Dline(m0, bs[end], ys[end]))

        push!(bs, findyintercept(ms[end], xs[end], ys[end]))
        push!(ys, z_fin[i])
        push!(xs, findx2Dline(m_zones[i], bs[end], ys[end]))
    end
    push!(ys, ub)
    push!(xs, ub)
    push!(ms, m0)
    push!(bs, findyintercept(ms[end], xs[end], ys[end]))

    return xs, ys, ms, bs
end

function getslope(len_z::Vector{T}, len_s::Vector{T}, total_len::T) where T

    ms = collect( len_z[i]/len_s[i] for i = 1:length(len_s) )

    m0 = (total_len-sum(len_z))/(total_len-sum(len_s))

    return ms, m0
end

function findyintercept(m::T, x::T, y::T)::T where T
    return y -m*x
end
## test.
# m = randn()
# b = randn()
# f = xx->(m*xx+b)
#
# x = 3.2
# y = f(x)
# b_rec = findyintercept(m,x,y)
# x_rec = findx2Dline(m,b_rec,y)

function findx2Dline(m::T, b::T, y::T)::T where T
    return (y-b)/m
end

"""
example usage:
D = 5 #2 # number of regions.
lb = -1
ub = 1
scale = 1.0 #0.5

# amount of input region used to map to z_lens.
input_range_percentage = 0.95
c = input_range_percentage*(ub-lb)

# generate random boundary points that define the regions.
rand_gen_func = xx->MonotoneMaps.convertcompactdomain(xx, 0.0, 1.0, scale*lb, scale*ub)
z = sort(collect(rand_gen_func(rand()) for d = 1:D))

# sums to c.
z_lens = ones(D) .* 2.0 #0.2
z_lens = z_lens .* (c/sum(z_lens))

# process z points.
z_pts = collect( [z[i]-z_lens[i]/2; z[i]+z_lens[i]/2] for i = 1:D )
z_st0 = collect( z_pts[i][1] for i = 1:length(z_pts) )
z_fin0 = collect( z_pts[i][2] for i = 1:length(z_pts) )

# combine intervals such that they are not overlapping.
z_st, z_fin = MonotoneMaps.mergesuchthatnooverlap(z_st0, z_fin0)
"""
function mergesuchthatnooverlap(z_st_in::Vector{T}, z_fin_in::Vector{T};
    lb::T = -one(T), ub::T = one(T),
    portion_for_shift::T = 0.95) where T

    inds = sortperm(z_st_in)
    z_st = z_st_in[inds]
    z_fin = z_fin_in[inds]

    @assert length(z_st) == length(z_fin)

    out_st = Vector{T}(undef, 0)
    out_fin = Vector{T}(undef, 0)

    fin = -Inf
    for j = 1:length(z_st)

        if fin < z_st[j]
            st = z_st[j]
            fin = findendnooverlap(z_fin[j], z_st[j:end], z_fin[j:end])

            push!(out_st, st)
            push!(out_fin, fin)
        end
    end

    # # sanity check.
    # c = (ub-lb)*portion_for_shift
    # C = sum( out_fin[i]-out_st[i] for i = 1:length(out_fin) )
    # @assert norm(C - c) < 1e-12

    ϵ = 1e-9
    clamp!(out_st, lb+ϵ, ub-ϵ)
    clamp!(out_fin, lb+ϵ, ub-ϵ)

    return out_st, out_fin
end

function findendnooverlap(fin::T, z_st::Vector{T}, z_fin::Vector{T})::T where T
    for i = 2:length(z_st)
        if z_st[i] < fin
            fin = z_fin[i]
        else
            return fin
        end
    end

    return fin
end

"""
Given a few candidate shift zones (a center and window per zone), return a monotonic map that
emphasizes the candidate shift zones. The amount of emphasis is controlled by portion_for_shift.
"""
function setupshiftwarpmap(shift_centers::Vector{T},
    shift_windows::Vector{T};
    portion_for_shift::T = 0.95) where T

    portion_for_shift = clamp(portion_for_shift, 1e-5, one(T)-1e-5)

    # set up.
    N_zones = length(shift_centers)
    @assert length(shift_windows) == N_zones

    c = portion_for_shift*2

    z_pts = collect( [shift_centers[i]-shift_windows[i]/2;
    shift_centers[i]+shift_windows[i]/2] for i = 1:N_zones )

    z_st0 = collect( z_pts[i][1] for i = 1:N_zones )
    z_fin0 = collect( z_pts[i][2] for i = 1:N_zones )

    # combine intervals such that they are not overlapping.
    z_st, z_fin = mergesuchthatnooverlap(z_st0, z_fin0)

    xs, ys, ms, bs, len_s, len_z = getpiecewiselines(z_st, z_fin, c)
    #f = xx->evalpiecewise2Dlinearfunc(xx, xs, ys, ms, bs)

    return xs, ys, ms, bs, len_s, len_z
end

####### front end.

function getunityobj(dummy_val::T, ::Val{:PiecewiseLinear}) where T <: Real

    c::T = 0.5 # amount of input region used to map to z_lens.
    z::Vector{T} = [-0.5; 0.5]
    z_lens::Vector{T} = ones(2) .* 0.5

    return Piecewise2DLineType(setupshiftwarpmap(z, z_lens; portion_for_shift = c)...)
end

"""
allocates γ percentage of range to be around p0.
"""
function getproximityobj(p0_inp::T,
    input_range_percentage::T, ::Val{:PiecewiseLinear};
    guard_lb = -0.99, guard_ub = 0.99) where T <: Real

    if p0_lb < p0 < p0_ub
        # TODO I am here.
    end
    c = input_range_percentage*2
    z::Vector{T} = [-0.5; 0.5]
    z_lens::Vector{T} = ones(2) .* 0.5

    return Piecewise2DLineType(setupshiftwarpmap(z, z_lens; portion_for_shift = c)...)
end

# function fetchwarpparms(shift_centers_set::Vector{Vector{T}},
#     shift_windows_set::Vector{Vector{T}};
#     portion_for_shift::T = 0.95) where T <: Real
#
#     N_shifts = length(shift_centers_set)
#     @assert length(shift_windows_set) == N_shifts
#
#     warp_param_set = Vector{Piecewise2DLineType{T}}(undef, N_shifts)
#
#     for i = 1:N_shifts
#         shift_centers = shift_centers_set[i]
#         shift_windows = shift_windows_set[i]
#         @assert length(shift_centers) == length(shift_windows)
#
#         tmp = setupshiftwarpmap(shift_centers, shift_windows; portion_for_shift = portion_for_shift)
#         warp_param_set[i] = Piecewise2DLineType(tmp...)
#     end
#
#     return warp_param_set
# end

# piecve-wise lienar version.
function setupsimilarmonotonemap(shift_centers_set::Vector{Vector{T}},
    shift_windows_set::Vector{Vector{T}};
    portion_for_shift::T = 0.95) where T <: Real

    N_shifts = length(shift_centers_set)
    @assert length(shift_windows_set) == N_shifts

    warp_param_set = Vector{Piecewise2DLineType{T}}(undef, N_shifts)

    for i = 1:N_shifts
        shift_centers = shift_centers_set[i]
        shift_windows = shift_windows_set[i]
        @assert length(shift_centers) == length(shift_windows)

        tmp = setupshiftwarpmap(shift_centers, shift_windows; portion_for_shift = portion_for_shift)
        warp_param_set[i] = Piecewise2DLineType(tmp...)
    end

    return warp_param_set
end

################# application as a proximity-encouraging map.

# uses lb = -1, ub = 1.
function setupproximitymap(focus_center::T, focus_radius::T;
    input_range_percentage = 0.9,
    scale::T = one(T)) where T <: Real

    @assert focus_radius < 1
    @assert -1 < focus_center < 1

    c::T = input_range_percentage*2

    focus_lower = focus_center - focus_radius
    focus_upper = focus_center + focus_radius

    xs, ys, ms, bs, len_s, len_z = MonotoneMaps.getpiecewiselines([ focus_lower; ], [ focus_upper; ], c)
    f = xx->MonotoneMaps.evalpiecewise2Dlinearfunc(xx, xs, ys, ms, bs)*scale
    finv = xx->MonotoneMaps.evalinversepiecewise2Dlinearfunc(xx/scale, xs, ys, ms, bs)

    return f, finv
end

function setupproximitymapsitp(focus_lb::T, focus_ub::T, focus_radius::T;
    input_range_percentage::T = 0.9,
    scale::T = one(T)) where T <: Real

    @assert (focus_ub-focus_lb) > focus_radius

    focus_center = focus_lb + focus_radius
    f1, finv1 = setupproximitymap(focus_center, focus_radius;
        input_range_percentage = input_range_percentage,
        scale = scale)

    focus_center = focus_ub - focus_radius
    f2, finv2 = setupproximitymap(focus_center, focus_radius;
        input_range_percentage = input_range_percentage,
        scale = scale)

    #
    f = itpproximitymap(f1, f2, focus_center)

    return f
end

# function evallinearitp(x::T, x1::T, y1::T, x2::T, y2::T)::T where T
#     return (y1*(x2-x) + y2*(x-x1))/(x2-x1)
# end

function itpproximitymap(f1, f2, y::T) where T
    f2_t = f2(t)
    f1_t = f1(t)

    a = (f2_t - y) / (f2_t - f1_t)
    if f1_t > y
        a = zero(T)
    else
        f2_t < y
        a = one(T)
    end

    f = tt->(a*f1(tt) + (1-a)*f2(t))

    return f
end

#####################

"""
    getendomorphismpiecewiselinear(p_lb::T, p_ub::T, range_percentage::T;
        N_itp_samples::Int = 10,
        domain_percentage::T = 0.9) where T

Get the parameters for a set of `N_itp_samples` piece-wise linear functions that has a focus interval in its range between `p_lb + range_percentage` and `p_ub - range_percentage`.

f = xx->MonotoneMaps.evalpiecewise2Dlinearfunc(xx, xs, ys, ms, bs)*scale
maps [lb, ub] to [lb, ub]*scale.

# Arguments
Constraint: -1 <= `p_lb`` < `p_ub`` <= 1
- `p_lb::T`: lower bound for the domain and range for each function.
- `p_ub::T`: upper bound for the domain and range for each function.
- `range_percentage::T`: The percentage of the rain to for the focus interval.
- `N_itp_samples::Int`: The number of functions to fit, each with a focus interval that uniformly covers `p_lb + range_percentage` and `p_ub - range_percentage`.
- `domain_percentage::T`: The percentage of the domain to map to the focus interval in the range.

# Example
See /examples/logistic-logit_fit.jl and https://royccwang.github.io/MonotoneMaps.jl/

"""
function getendomorphismpiecewiselinear(p_lb::T, p_ub::T, range_percentage::T;
    N_itp_samples::Int = 10,
    domain_percentage::T = 0.9) where T

    # set up.
    c = domain_percentage*(p_ub-p_lb)
    window = range_percentage*(p_ub-p_lb)

    #p_range = LinRange(p_lb + window + ϵ, p_ub - window - ϵ, N_itp_samples)
    p_range = LinRange(p_lb + window, p_ub - window, N_itp_samples)

    infos = Vector{Piecewise2DLineType{T}}(undef, length(p_range))
    zs = Vector{Vector{T}}(undef, length(p_range))

    for (i,p) in enumerate(p_range)

        z_st = [p - window;]
        z_fin = [p + window;]

        xs, ys, ms, bs, len_s, len_z = getpiecewiselines(z_st, z_fin, c; lb = p_lb, ub = p_ub)
        infos[i] = Piecewise2DLineType(xs, ys, ms, bs, len_s, len_z)

        zs[i] = [z_st[1]; z_fin[1]]
    end

    return infos, zs, p_range
end


text:
mutable struct GFUpdateState
    prev_trace::DynamicDSLTrace
    trace::DynamicDSLTrace
    constraints::Any
    weight::Float64
    visitor::AddressVisitor
    params::Dict{Symbol,Any}
    discard::DynamicChoiceMap
end

function GFUpdateState(gen_fn, args, prev_trace, constraints, params)
    visitor = AddressVisitor()
    discard = choicemap()
    trace = DynamicDSLTrace(gen_fn, args)
    GFUpdateState(prev_trace, trace, constraints,
        0., visitor, params, discard)
end

function traceat(state::GFUpdateState, dist::Distribution{T},
                 args::Tuple, key) where {T}

    local prev_retval::T
    local retval::T

    # check that key was not already visited, and mark it as visited
    visit!(state.visitor, key)

    # check for previous choice at this key
    has_previous = has_choice(state.prev_trace, key)
    if has_previous
        prev_choice = get_choice(state.prev_trace, key)
        prev_retval = prev_choice.retval
        prev_score = prev_choice.score
    end

    # check for constraints at this key
    constrained = has_value(state.constraints, key)
    !constrained && check_no_submap(state.constraints, key)

    # record the previous value as discarded if it is replaced
    if constrained && has_previous
        set_value!(state.discard, key, prev_retval)
    end

    # get return value
    if constrained
        retval = get_value(state.constraints, key)
    elseif has_previous
        retval = prev_retval
    else
        retval = random(dist, args...)
    end

    # compute logpdf
    score = logpdf(dist, retval, args...)

    # update the weight
    if has_previous
        state.weight += score - prev_score
    elseif constrained
        state.weight += score
    end

    # add to the trace
    add_choice!(state.trace, key, retval, score)

    retval
end

function traceat(state::GFUpdateState, gen_fn::GenerativeFunction{T,U},
                 args::Tuple, key) where {T,U}

    local prev_subtrace::U
    local subtrace::U
    local retval::T

    # check key was not already visited, and mark it as visited
    visit!(state.visitor, key)

    # check for constraints at this key
    check_no_value(state.constraints, key)
    constraints = get_submap(state.constraints, key)

    # get subtrace
    has_previous = has_call(state.prev_trace, key)
    if has_previous
        prev_call = get_call(state.prev_trace, key)
        prev_subtrace = prev_call.subtrace
        get_gen_fn(prev_subtrace) === gen_fn || gen_fn_changed_error(key)
        (subtrace, weight, _, discard) = update(prev_subtrace,
            args, map((_) -> UnknownChange(), args), constraints)
    else
        (subtrace, weight) = generate(gen_fn, args, constraints)
    end

    # update the weight
    state.weight += weight

    # update discard
    if has_previous
        set_submap!(state.discard, key, discard)
    end

    # add to the trace
    add_call!(state.trace, key, subtrace)

    # get return value
    retval = get_retval(subtrace)

    retval
end

function splice(state::GFUpdateState, gen_fn::DynamicDSLFunction,
                args::Tuple)
    prev_params = state.params
    state.params = gen_fn.params
    retval = exec(gen_fn, state, args)
    state.params = prev_params
    retval
end

function update_delete_recurse(prev_trie::Trie{Any,ChoiceOrCallRecord},
                               visited::EmptySelection)
    score = 0.
    for (key, choice_or_call) in get_leaf_nodes(prev_trie)
        score += choice_or_call.score
    end
    for (key, subtrie) in get_internal_nodes(prev_trie)
        score += update_delete_recurse(subtrie, EmptySelection())
    end
    score
end

function update_delete_recurse(prev_trie::Trie{Any,ChoiceOrCallRecord},
                               visited::DynamicSelection)
    score = 0.
    for (key, choice_or_call) in get_leaf_nodes(prev_trie)
        if !(key in visited)
            score += choice_or_call.score
        end
    end
    for (key, subtrie) in get_internal_nodes(prev_trie)
        subvisited = visited[key]
        score += update_delete_recurse(subtrie, subvisited)
    end
    score
end

function update_delete_recurse(prev_trie::Trie{Any,ChoiceOrCallRecord},
                               visited::AllSelection)
    0.
end

function add_unvisited_to_discard!(discard::DynamicChoiceMap,
                                   visited::DynamicSelection,
                                   prev_choices::ChoiceMap)
    for (key, value) in get_values_shallow(prev_choices)
        if !(key in visited)
            @assert !has_value(discard, key)
            @assert isempty(get_submap(discard, key))
            set_value!(discard, key, value)
        end
    end
    for (key, submap) in get_submaps_shallow(prev_choices)
        @assert !has_value(discard, key)
        if key in visited
            # the recursive call to update already handled the discard
            # for this entire submap
            continue
        else
            subvisited = visited[key]
            if isempty(subvisited)
                # none of this submap was visited, so we discard the whole thing
                @assert isempty(get_submap(discard, key))
                set_submap!(discard, key, submap)
            else
                subdiscard = get_submap(discard, key)
                add_unvisited_to_discard!(
                    isempty(subdiscard) ? choicemap() : subdiscard,
                    subvisited, submap)
                set_submap!(discard, key, subdiscard)
            end
        end
    end
end

function update(trace::DynamicDSLTrace, arg_values::Tuple, arg_diffs::Tuple,
                constraints::ChoiceMap)
    gen_fn = trace.gen_fn
    state = GFUpdateState(gen_fn, arg_values, trace, constraints, gen_fn.params)
    retval = exec(gen_fn, state, arg_values)
    set_retval!(state.trace, retval)
    visited = get_visited(state.visitor)
    state.weight -= update_delete_recurse(trace.trie, visited)
    add_unvisited_to_discard!(state.discard, visited, get_choices(trace))
    if !all_visited(visited, constraints)
        error("Did not visit all constraints")
    end
    (state.trace, state.weight, UnknownChange(), state.discard)
end


text:
using Plots

@userplot PlotOscillator
@recipe function f(p::PlotOscillator)
    if length(p.args) != 2 || !(typeof(p.args[1]) <: AbstractVector) || !(typeof(p.args[2]) <: AbstractMatrix)
        error("Oscillator plots should be called with two arrays. Got: $(typeof(p.args))")
    end
    local t = p.args[1]
    local q = p.args[2]
    local h = zero(t)

    for i in eachindex(t)
        h[i] = hamiltonian(t[i], q[:,i], k)
    end

    size   := (800,600)
    layout := @layout [xPlot
                       yPlot
                       hPlot]
    legend := :none

    @series begin
        subplot := 1
        xlabel := "t"
        ylabel := "x"
        t, q[1,:]
    end

    @series begin
        subplot := 2
        xlabel := "t"
        ylabel := "y"
        t, q[2,:]
    end

    @series begin
        subplot := 3
        xlabel := "t"
        ylabel := "h"
        t, h
    end
end

@userplot PlotOscillatorOverview
@recipe function f(p::PlotOscillatorOverview)
    if length(p.args) != 1 || !(typeof(p.args[1]) <: AbstractArray)
        error("Oscillator plots should be called with an array. Got: $(typeof(p.args))")
    end

    local q = p.args[1]

    size   := (800, div(size(q,2), 2)*300)

    layout := (div(size(q,2), 2), 2)
    legend := :none

    for k in axes(q,2)
        @series begin
            subplot := k
            xlabel  := "x"
            ylabel  := "y"
            q[1,k,:], q[2,k,:]
        end
    end
end


text:
@testset "Evolution" begin
    @testset "Exponent function" begin
    @test evolve_operator(zeros(4, 4), 1.) ≈ Matrix{Float64}(I,4,4)
    A = [1. im 2;
    3 4 5im;
    6 7 8]
    @test evolve_operator(A,4.2) ≈ exp(4.2*A)

    # type check
    @test_throws MethodError evolve_operator(spzeros(4, 4), 1.)
  end

  @testset "Evolution functions" begin
    #size 4x4
    A = zeros(ComplexF64,4,4)
    A[2,2] = 1/2
    A[4,4] = 1/2
    S = [0.0+0.0im 0.5+0.0im 0.0+0.0im 0.0+0.0im 0.5+0.0im 0.5+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im;
        -0.5+0.0im -0.25+0.0im 0.0+0.5im 0.0+0.0im 0.0+0.0im 0.5+0.0im 0.5+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im;
        0.0+0.0im 0.0+0.5im -0.25+0.0im 0.0+0.5im 0.0+0.0im 0.0+0.0im 0.5+0.0im 0.5+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im;
        0.0+0.0im 0.0+0.0im 0.0+0.5im -0.25+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.5+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im;
        -0.5+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im -0.25+0.0im 0.5+0.0im 0.0+0.0im 0.0+0.0im 0.0-0.5im 0.5+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im;
        0.0+0.0im -0.5+0.0im 0.0+0.0im 0.0+0.0im -0.5+0.0im -0.5+0.0im 0.0+0.5im 0.0+0.0im 0.0+0.0im 0.0-0.5im 0.5+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im;
        0.0+0.0im 0.0+0.0im -0.5+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.5im -0.5+0.0im 0.0+0.5im 0.0+0.0im 0.0+0.0im 0.0-0.5im 0.5+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im;
        0.0+0.0im 0.0+0.0im 0.0+0.0im -0.5+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.5im -0.5+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0-0.5im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im;
        0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0-0.5im 0.0+0.0im 0.0+0.0im 0.0+0.0im -0.25+0.0im 0.5+0.0im 0.0+0.0im 0.0+0.0im 0.0-0.5im 0.5+0.0im 0.0+0.0im 0.0+0.0im;
        0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0-0.5im 0.0+0.0im 0.0+0.0im -0.5+0.0im -0.5+0.0im 0.0+0.5im 0.0+0.0im 0.0+0.0im 0.0-0.5im 0.5+0.0im 0.0+0.0im;
        0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0-0.5im 0.0+0.0im 0.0+0.0im 0.0+0.5im -0.5+0.0im 0.0+0.5im 0.0+0.0im 0.0+0.0im 0.0-0.5im 0.5+0.0im;
        0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0-0.5im 0.0+0.0im 0.0+0.0im 0.0+0.5im -0.5+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0-0.5im;
        0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0-0.5im 0.0+0.0im 0.0+0.0im 0.0+0.0im -0.25+0.0im 0.5+0.0im 0.0+0.0im 0.0+0.0im;
        0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0-0.5im 0.0+0.0im 0.0+0.0im -0.5+0.0im -0.5+0.0im 0.0+0.5im 0.0+0.0im;
         0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0-0.5im 0.0+0.0im 0.0+0.0im 0.0+0.5im -0.5+0.0im 0.0+0.5im;
         0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0-0.5im 0.0+0.0im 0.0+0.0im 0.0+0.5im -0.5+0.0im]

    result_matrix42 = [0.495478+1.2483e-17im -0.162938+0.0285877im 0.0104949-0.055023im -0.000950804+0.00066496im;
     -0.162938-0.0285877im 0.271419-1.3821e-17im -0.000400553+0.107483im 0.00474894+0.00823693im;
     0.0104949+0.055023im -0.000400553-0.107483im 0.155099-2.82874e-18im 0.00737685+0.0355666im;
    -0.000950804-0.00066496im 0.00474894-0.00823693im 0.00737685-0.0355666im 0.0780051-2.51094e-18im]
    result_matrix10 = [0.34111+5.57946e-18im -0.121394+0.0872953im -0.0693429+0.0173517im -0.0338042-0.0543639im;
     -0.121394-0.0872953im 0.248919-4.2406e-17im 0.0565208+0.0998314im -0.0158367+0.0694652im
     -0.0693429-0.0173517im 0.0565208-0.0998314im 0.223367-1.11563e-17im 0.0486222+0.0838723im;
     -0.0338042+0.0543639im -0.0158367-0.0694652im 0.0486222-0.0838723im 0.186604-2.90591e-18im]
    state = zeros(4,4)
    state[3,3] = 1
    #trivial evolutions
    @test evolve(zeros(16, 16), A, 0.) ≈ A
    @test evolve(zeros(16, 16), sparse(A), 0.) ≈ A
    @test evolve(spzeros(16, 16), A, 0.) ≈ A
    @test evolve(spzeros(16, 16), A, [0., 5., 10.])[1] ≈ A
    @test evolve(spzeros(16, 16), A, [0., 5., 10.])[2] ≈ A
    @test evolve(spzeros(16, 16), A, [0., 5., 10.])[3] ≈ A

    @test evolve(S, A, 0.) ≈ A
    @test evolve(Matrix{Float64}(I, 16, 16), A) ≈ A
    @test isapprox(evolve(S, A, 4.2), result_matrix42, atol=1e-5)
    @test isapprox(evolve(S, A, [0., 4.2, 10.]), [A, result_matrix42, result_matrix10], atol=1e-5)

  end
end


text:
function SparseFieldVector(arg0::Field)
    return SparseFieldVector((Field,), arg0)
end

function SparseFieldVector(arg0::Field, arg1::Vector{FieldElement})
    return SparseFieldVector((Field, Vector{FieldElement}), arg0, arg1)
end

function SparseFieldVector(arg0::Field, arg1::jint)
    return SparseFieldVector((Field, jint), arg0, arg1)
end

function SparseFieldVector(arg0::Field, arg1::jint, arg2::jint)
    return SparseFieldVector((Field, jint, jint), arg0, arg1, arg2)
end

function SparseFieldVector(arg0::SparseFieldVector)
    return SparseFieldVector((SparseFieldVector,), arg0)
end

function add(obj::SparseFieldVector, arg0::FieldVector)
    return jcall(obj, "add", FieldVector, (FieldVector,), arg0)
end

function add(obj::SparseFieldVector, arg0::SparseFieldVector)
    return jcall(obj, "add", FieldVector, (SparseFieldVector,), arg0)
end

function append(obj::SparseFieldVector, arg0::FieldElement)
    return jcall(obj, "append", FieldVector, (FieldElement,), arg0)
end

function append(obj::SparseFieldVector, arg0::FieldVector)
    return jcall(obj, "append", FieldVector, (FieldVector,), arg0)
end

function append(obj::SparseFieldVector, arg0::SparseFieldVector)
    return jcall(obj, "append", FieldVector, (SparseFieldVector,), arg0)
end

function copy(obj::SparseFieldVector)
    return jcall(obj, "copy", FieldVector, ())
end

function dot_product(obj::SparseFieldVector, arg0::FieldVector)
    return jcall(obj, "dotProduct", FieldElement, (FieldVector,), arg0)
end

function ebe_divide(obj::SparseFieldVector, arg0::FieldVector)
    return jcall(obj, "ebeDivide", FieldVector, (FieldVector,), arg0)
end

function ebe_multiply(obj::SparseFieldVector, arg0::FieldVector)
    return jcall(obj, "ebeMultiply", FieldVector, (FieldVector,), arg0)
end

function equals(obj::SparseFieldVector, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_dimension(obj::SparseFieldVector)
    return jcall(obj, "getDimension", jint, ())
end

function get_entry(obj::SparseFieldVector, arg0::jint)
    return jcall(obj, "getEntry", FieldElement, (jint,), arg0)
end

function get_field(obj::SparseFieldVector)
    return jcall(obj, "getField", Field, ())
end

function get_sub_vector(obj::SparseFieldVector, arg0::jint, arg1::jint)
    return jcall(obj, "getSubVector", FieldVector, (jint, jint), arg0, arg1)
end

function hash_code(obj::SparseFieldVector)
    return jcall(obj, "hashCode", jint, ())
end

function map_add(obj::SparseFieldVector, arg0::FieldElement)
    return jcall(obj, "mapAdd", FieldVector, (FieldElement,), arg0)
end

function map_add_to_self(obj::SparseFieldVector, arg0::FieldElement)
    return jcall(obj, "mapAddToSelf", FieldVector, (FieldElement,), arg0)
end

function map_divide(obj::SparseFieldVector, arg0::FieldElement)
    return jcall(obj, "mapDivide", FieldVector, (FieldElement,), arg0)
end

function map_divide_to_self(obj::SparseFieldVector, arg0::FieldElement)
    return jcall(obj, "mapDivideToSelf", FieldVector, (FieldElement,), arg0)
end

function map_inv(obj::SparseFieldVector)
    return jcall(obj, "mapInv", FieldVector, ())
end

function map_inv_to_self(obj::SparseFieldVector)
    return jcall(obj, "mapInvToSelf", FieldVector, ())
end

function map_multiply(obj::SparseFieldVector, arg0::FieldElement)
    return jcall(obj, "mapMultiply", FieldVector, (FieldElement,), arg0)
end

function map_multiply_to_self(obj::SparseFieldVector, arg0::FieldElement)
    return jcall(obj, "mapMultiplyToSelf", FieldVector, (FieldElement,), arg0)
end

function map_subtract(obj::SparseFieldVector, arg0::FieldElement)
    return jcall(obj, "mapSubtract", FieldVector, (FieldElement,), arg0)
end

function map_subtract_to_self(obj::SparseFieldVector, arg0::FieldElement)
    return jcall(obj, "mapSubtractToSelf", FieldVector, (FieldElement,), arg0)
end

function outer_product(obj::SparseFieldVector, arg0::FieldVector)
    return jcall(obj, "outerProduct", FieldMatrix, (FieldVector,), arg0)
end

function outer_product(obj::SparseFieldVector, arg0::SparseFieldVector)
    return jcall(obj, "outerProduct", FieldMatrix, (SparseFieldVector,), arg0)
end

function projection(obj::SparseFieldVector, arg0::FieldVector)
    return jcall(obj, "projection", FieldVector, (FieldVector,), arg0)
end

function set(obj::SparseFieldVector, arg0::FieldElement)
    return jcall(obj, "set", void, (FieldElement,), arg0)
end

function set_entry(obj::SparseFieldVector, arg0::jint, arg1::FieldElement)
    return jcall(obj, "setEntry", void, (jint, FieldElement), arg0, arg1)
end

function set_sub_vector(obj::SparseFieldVector, arg0::jint, arg1::FieldVector)
    return jcall(obj, "setSubVector", void, (jint, FieldVector), arg0, arg1)
end

function subtract(obj::SparseFieldVector, arg0::FieldVector)
    return jcall(obj, "subtract", FieldVector, (FieldVector,), arg0)
end

function subtract(obj::SparseFieldVector, arg0::SparseFieldVector)
    return jcall(obj, "subtract", SparseFieldVector, (SparseFieldVector,), arg0)
end

function to_array(obj::SparseFieldVector)
    return jcall(obj, "toArray", Vector{FieldElement}, ())
end

function walk_in_default_order(obj::SparseFieldVector, arg0::FieldVectorChangingVisitor)
    return jcall(obj, "walkInDefaultOrder", FieldElement, (FieldVectorChangingVisitor,), arg0)
end

function walk_in_default_order(obj::SparseFieldVector, arg0::FieldVectorChangingVisitor, arg1::jint, arg2::jint)
    return jcall(obj, "walkInDefaultOrder", FieldElement, (FieldVectorChangingVisitor, jint, jint), arg0, arg1, arg2)
end

function walk_in_default_order(obj::SparseFieldVector, arg0::FieldVectorPreservingVisitor)
    return jcall(obj, "walkInDefaultOrder", FieldElement, (FieldVectorPreservingVisitor,), arg0)
end

function walk_in_default_order(obj::SparseFieldVector, arg0::FieldVectorPreservingVisitor, arg1::jint, arg2::jint)
    return jcall(obj, "walkInDefaultOrder", FieldElement, (FieldVectorPreservingVisitor, jint, jint), arg0, arg1, arg2)
end

function walk_in_optimized_order(obj::SparseFieldVector, arg0::FieldVectorChangingVisitor)
    return jcall(obj, "walkInOptimizedOrder", FieldElement, (FieldVectorChangingVisitor,), arg0)
end

function walk_in_optimized_order(obj::SparseFieldVector, arg0::FieldVectorChangingVisitor, arg1::jint, arg2::jint)
    return jcall(obj, "walkInOptimizedOrder", FieldElement, (FieldVectorChangingVisitor, jint, jint), arg0, arg1, arg2)
end

function walk_in_optimized_order(obj::SparseFieldVector, arg0::FieldVectorPreservingVisitor)
    return jcall(obj, "walkInOptimizedOrder", FieldElement, (FieldVectorPreservingVisitor,), arg0)
end

function walk_in_optimized_order(obj::SparseFieldVector, arg0::FieldVectorPreservingVisitor, arg1::jint, arg2::jint)
    return jcall(obj, "walkInOptimizedOrder", FieldElement, (FieldVectorPreservingVisitor, jint, jint), arg0, arg1, arg2)
end


text:
struct Uninomial <: AbstractMonomial
    v::IDType
    d::UInt
end

Base.isless(x::Uninomial, y::Uninomial) = isless(degree(x), degree(y))
degree(m::Uninomial) = Int(m.d)
var(m::Uninomial) = m.v

coeff(m::Uninomial) = 1
Base.one(m::Uninomial) = Uninomial(m.v, 0)
Base.show(io::IO, m::Uninomial) = print_single_monomial(io, var(m), degree(m))

struct Uniterm{T,M<:AbstractMonomial} <: AbstractTerm
    coeff::T
    uninomial::M
end
Base.convert(::Type{<:Uniterm}, t::Uninomial) = Uniterm(t)
#Base.convert(::Type{<:Uniterm}, t::CoeffType) = Uniterm(t, Uninomial())
Uniterm(t::Uninomial) = Uniterm(coeff(t), t)

coeff(t::Uniterm) = t.coeff
monomial(t::Uniterm) = t.uninomial
var(t::Uniterm) = var(monomial(t))

struct SparsePoly{T} <: AbstractPolynomial
    coeffs::Vector{T}
    exps::Vector{UInt}
    v::IDType

    SparsePoly(coeffs::Vector{T}, exps::Vector{UInt}, v::IDType) where T = SparsePoly{T}(coeffs, exps, v)
    function SparsePoly{T}(coeffs::Vector{T}, exps::Vector{UInt}, v::IDType) where T
        length(coeffs) == length(exps) || error("coeffs and exps' length must match!")
        new{T}(coeffs, exps, v)
    end
end
Base.convert(::Type{<:SparsePoly}, m::Uninomial) = SparsePoly(Uniterm(m))
Base.convert(::Type{<:SparsePoly}, t::Uniterm) = SparsePoly(t)
#Base.convert(::Type{<:SparsePoly}, t::CoeffType) = SparsePoly(convert(Uniterm, t))
SparsePoly(t::Uniterm) = SparsePoly([coeff(t)], UInt[degree(t)], var(t))
SparsePoly(c::Number, v) = SparsePoly([c], [zero(UInt)], v)
const EMPTY_EXPS = UInt[]
#SparsePoly(t::Uniterm, id::IDType) = SparsePoly(typeof(coeff(t))[], EMPTY_EXPS, id)
Base.similar(p::SparsePoly) = SparsePoly(similar(coeffs(p)), similar(p.exps), var(p))

var(p::SparsePoly) = p.v
coeffs(p::SparsePoly) = p.coeffs
coeff(p::SparsePoly, i) = coeffs(p)[i]
term(p::SparsePoly, i) = Uniterm(p.coeffs[i], Uninomial(var(p), p.exps[i]))
terms(p::SparsePoly) = (term(p, i) for i in eachindex(coeffs(p)))
lt(p::SparsePoly) = term(p, 1)
lc(p::SparsePoly) = coeff(p, 1)
degree(p::SparsePoly) = iszero(p) ? -1 : Int(p.exps[1])
Base.iszero(p::SparsePoly) = isempty(p.exps) || iszero(lt(p))

Base.zero(t::SparsePoly) = SparsePoly(empty(t.coeffs), empty(t.exps), var(t))
# TODO: this fails when t is zero
Base.one(t::SparsePoly) = one(lt(t))

Base.deleteat!(p::SparsePoly, i::Int) = (deleteat!(p.coeffs, i); deleteat!(p.exps, i); nothing)
Base.copy(p::SparsePoly) = SparsePoly(copy(coeffs(p)), copy(p.exps), var(p))
Base.copy(p::SparsePoly{<:AbstractPolynomial}) = SparsePoly(map(copy, coeffs(p)), copy(p.exps), var(p))
Base.:(==)(p::SparsePoly, q::SparsePoly) = all(x->x[1] == x[2], zip(terms(p), terms(q)))

function check_poly(x, y)
    v = var(x)
    v == var(y) || error("$x and $y contain different variables!")
    nothing
end

function print_poly_term(io, t)
    need_par = !isone(monomial(t))
    need_par && print(io, '(')
    show(io, coeff(t))
    need_par && print(io, ')')
    need_par && show(io, monomial(t))
    nothing
end

function Base.show(io::IO, p::SparsePoly{<:AbstractPolynomial})
    ts = terms(p)
    if isempty(ts)
        print(io, 0)
        return
    end
    n = length(ts)
    t1, tr = Iterators.peel(ts)
    print_poly_term(io, t1)
    for t in tr
        print(io, " + ")
        print_poly_term(io, t)
    end
end

#########################
# Arithmetic/Algorithms #
#########################

Base.:(^)(x::Uninomial, n::Integer) = Uninomial(var(x), degree(x)+n-1)

function Base.:(+)(x::Uninomial, y::Uninomial)
    check_poly(x, y)
    if x < y
        x, y = y, x
    end
    SparsePoly([1, 1], [degree(x), degree(y)], var(x))
end
function Base.:(+)(x::Uniterm, y::Uniterm)
    check_poly(x, y)
    if ismatch(x, y)
        c = x.coeff + y.coeff
        return iszero(c) ? SparsePoly(typeof(x)[], EMPTY_EXPS, var(x)) : SparsePoly(Uniterm(c, monomial(x)))
    else
        if iszero(x) && iszero(y)
            c = coeff(y)
            return SparsePoly(typeof(c)[], EMPTY_EXPS, var(x))
        elseif iszero(x)
            exp = degree(y)
            c = coeff(y)
            return SparsePoly([c], UInt[exp], var(x))
        elseif iszero(y)
            exp = degree(x)
            c = coeff(x)
            return SparsePoly([c], UInt[exp], var(x))
        else
            if x < y
                x, y = y, x
            end
            return SparsePoly([coeff(x), coeff(y)], UInt[degree(x), degree(y)], var(x))
        end
    end
end
Base.:(*)(x::CoeffType, y::Uninomial) = Uniterm(x, y)
Base.:(*)(x::Uninomial, y::CoeffType) = y * x

function Base.:(*)(x::Uninomial, y::Uninomial)
    check_poly(x, y)
    Uninomial(var(x), degree(x) + degree(y))
end
Base.:(*)(x::Uniterm, y::Uniterm) = Uniterm(coeff(x) * coeff(y), monomial(x) * monomial(y))

function Base.:*(p::SparsePoly, x::Uniterm)
    check_poly(p, x)
    if iszero(x)
        return SparsePoly(x)
    elseif isone(x)
        return p
    else
        cfs = similar(coeffs(p))
        exps = similar(p.exps)
        for (i, t) in enumerate(terms(p))
            nt = t * x
            cfs[i] = coeff(nt)
            exps[i] = degree(nt)
        end
        return SparsePoly(cfs, exps, var(x))
    end
end
Base.:*(x::CoeffType, y::Uniterm) = Uniterm(x * coeff(y), monomial(y))
Base.:*(x::Uniterm, y::CoeffType) = y * x

function Base.:*(x::SparsePoly, y::SparsePoly)
    if iszero(x)
        return x
    elseif iszero(y)
        return y
    else
        return sum(t->x * t, terms(y))
    end
end
smul(x, y::SparsePoly) = SparsePoly(x * coeffs(y), y.exps, var(y))
Base.:*(x::CoeffType, y::SparsePoly) = smul(x, y)
Base.:*(x::SparsePoly, y::CoeffType) = y * x
Base.:*(x::T, y::SparsePoly{T}) where {T<:AbstractPolynomial} = smul(x, y)
Base.:*(x::SparsePoly{T}, y::T) where {T<:AbstractPolynomial} = y * x

addcoef(x::Uniterm, c) = (c += coeff(x); return iszero(c), Uniterm(c, monomial(x)))
addcoef(x::Uniterm, c::Uniterm) = addcoef(x, coeff(c))
Base.:(-)(x::Uniterm) = Uniterm(-coeff(x), monomial(x))
sub!(p::SparsePoly, x::Uniterm) = add!(p, -x)
function add!(p::SparsePoly, x::Uniterm)
    iszero(x) && return p
    for (i, t) in enumerate(terms(p))
        if ismatch(t, x)
            iz, t = addcoef(t, x)
            if iz
                deleteat!(p, i)
            else
                p.coeffs[i] = coeff(t)
            end
            return p
        elseif t < x
            insert!(p.coeffs, i, coeff(x))
            insert!(p.exps, i, degree(x))
            return p
        end
    end
    push!(p.coeffs, coeff(x))
    push!(p.exps, degree(x))
    return p
end

Base.div(x::Uninomial, y::Uninomial) = (x/y)[1]
function Base.:(/)(x::Uninomial, y::Uninomial)
    check_poly(x, y)
    xd, yd = degree(x), degree(y)
    xd >= yd ? (Uninomial(var(x), xd - yd), false) : (x, true)
end

# dest = (p * a).^n
function mulpow!(dest, p::SparsePoly, a, n::Integer)
    @assert n >= 0
    for i in eachindex(p.exps)
        dest.coeffs[i] = p.coeffs[i] * a
        dest.exps[i] = p.exps[i] + n
    end
    return dest
end

function mul!(p::SparsePoly, l)
    cs = coeffs(p)
    for i in eachindex(cs)
        cs[i] *= l
    end
end

function pseudorem(p::SparsePoly, d::SparsePoly)
    check_poly(p, d)
    degree(p) < degree(d) && return p
    k = degree(p) - degree(d) + 1
    l = lc(d)
    dd = similar(d)
    p = copy(p)
    while !iszero(p) && degree(p) >= degree(d)
        s = mulpow!(dd, d, lc(p), degree(p) - degree(d))
        mul!(p, l)
        sub!(p, s)
        k -= 1
    end
    return l^k * p
end

function termwise_content(a::MPoly)
    ts = terms(a)
    length(ts) == 1 && return ts[1]
    g = gcd(ts[1], ts[2])
    isone(g) || for i in 3:length(ts)
        g = gcd(g, ts[i])
        isone(g) && break
    end
    g
end

function content(a::SparsePoly)
    cfs = coeffs(a)
    length(cfs) == 1 && return first(cfs)
    # If any coeff here is a multivariate polynomial with only one term, then we
    # just need to compute the termwize content.
    MP = eltype(cfs)
    if MP <: AbstractPolynomial
        for i in eachindex(cfs)
            ts = terms(cfs[i])
            if length(ts) == 1
                g = gcd(termwise_content(cfs[1]), termwise_content(cfs[2]))
                isone(g) || for i in 3:length(cfs)
                    g = gcd(g, termwise_content(cfs[i]))
                    isone(g) && break
                end
                return MP(g)
            end
        end
    end

    g = gcd(cfs[1], cfs[2])
    isone(g) || for i in 3:length(cfs)
        g = gcd(g, cfs[i])
        isone(g) && break
    end
    g
end

#function univariate_gcd(x::AbstractPolynomial, y::AbstractPolynomial)
#    while !iszero(y)
#        x = pseudorem(x, y)
#        x, y = y, x
#    end
#    return x#, a / x, b / x
#end

function Base.divrem(p::SparsePoly{T}, d::SparsePoly{T}) where T<:CoeffType
    p = copy(p)
    q = zero(p)
    r = zero(p)
    while !iszero(p)
        nx, fail = lt(p) / lt(d)
        if fail
            break
        else
            p -= d * nx
            q += nx
        end
    end
    return q, r
end

function divexact(a::SparsePoly{T}, b::T) where {T<:AbstractPolynomial}
    cfs = [divexact(c, b) for c in coeffs(a)]
    SparsePoly(cfs, a.exps, var(a))
end

divexact(c, b) = ((d, r) = divrem(c, b); @assert iszero(r); d;)
function divexact(a::SparsePoly{T}, b::T) where T
    cfs = [divexact(c, b) for c in coeffs(a)]
    SparsePoly(cfs, a.exps, var(a))
end

function Base.gcd(x::SparsePoly, y::SparsePoly)
    if degree(y) > degree(x)
        x, y = y, x
    end
    if iszero(y)
        return x
    elseif isone(y)
        return y
    end
    c1, x = contprim(x)
    c2, y = contprim(y)
    c = gcd(c1, c2)
    g = one(eltype(coeffs(x)))
    h = one(eltype(coeffs(x)))
    while true
        r = pseudorem(x, y)
        iszero(r) && break
        degree(r) == 0 && return SparsePoly(c, var(x))

        d = degree(x) - degree(y)
        x, y = y, divexact(r, g*h^d)
        g = lc(x)
        h = d > 1 ? divexact(g^d, h^(d - 1)) : h^(1 - d)*g^d
    end
    return c*primpart(y)
end

primpart(p) = divexact(p, content(p))
contprim(p) = (c = content(p); (c, divexact(p, c)))

##############
# Conversion #
##############

# mpoly2poly(5x⁴ + 2x³ + x² + xy + y^2, 0x00001 [y])
# ->
#              ---------------- constants
#              |
# y^2 + xy + (5x⁴ + 2x³ + x²)
#
# M2P(a, v)
function emplace_back!(poly, ts, perm, chunk_start_idx, idx, olddegree)
    v = var(poly)
    if olddegree > 0
        t = ts[perm[chunk_start_idx]]
        coeff = MPoly(term2polycoeff(t, v))
        for i in chunk_start_idx+1:idx-1
            t = ts[perm[i]]
            add!(coeff, term2polycoeff(t, v))
        end
    else
        coeff = MPoly(ts[perm[chunk_start_idx]])
        for i in chunk_start_idx+1:idx-1
            add!(coeff, ts[perm[i]])
        end
    end
    push!(poly.exps, olddegree)
    push!(poly.coeffs, coeff)
    nothing
end
term2polycoeff(t::Term, v::IDType) = Term(coeff(t), rmid(monomial(t), v))
function SparsePoly(p::MPoly, v::IDType)
    ts = terms(p)
    pows = map(ts) do t
        degree(monomial(t), v)
    end
    perm = sortperm(pows, rev=true)
    coeffs = typeof(p)[]
    exps = UInt[]
    poly = SparsePoly(coeffs, exps, v)

    olddegree = pows[perm[1]]
    chunk_start_idx = idx = 1
    while idx <= length(perm)
        degree = pows[perm[idx]]
        if olddegree != degree # new chunk
            emplace_back!(poly, ts, perm, chunk_start_idx, idx, olddegree)
            chunk_start_idx = idx
            olddegree = degree
        end
        idx += 1
    end
    if !isempty(chunk_start_idx:idx-1) # remember to handle the last one
        emplace_back!(poly, ts, perm, chunk_start_idx, idx, olddegree)
    end
    return poly
end

univariate_to_multivariate(g::SparsePoly{<:AbstractPolynomial}) = univariate_to_multivariate(g, monomialtype(lc(g)))
function univariate_to_multivariate(g, ::Type{<:Monomial})
    cfs = coeffs(g)
    eps = g.exps
    v = var(g)
    @assert !isempty(eps)
    #sum(zip(cfs, eps)) do (c, e)
    #    c * Monomial(fill(v, e))
    #end
    s = cfs[1] * Monomial(fill(v, eps[1]))
    for i in 2:length(cfs)
        c = cfs[i]
        e = eps[i]
        add!(s, c * Monomial(fill(v, e)))
    end
    return s
end

function univariate_to_multivariate(g, P::Type{<:PackedMonomial})
    cfs = coeffs(g)
    eps = g.exps
    v = var(g)
    @assert !isempty(eps)
    s = cfs[1] * P(v)^eps[1]
    for i in 2:length(cfs)
        c = cfs[i]
        e = eps[i]
        add!(s, c * P(v)^e)
    end
    return s
end


text:
# Convert json Data form Dimbreath/GenshinData into Matrix/Array

module JSONPhraser
export ReliAffixECDConverter,ReliLevelECDConverter

using JSON

if isdefined(@__MODULE__, :LookUpTable)
    VectorFunctions isa Module || error("LookUpTable is present and it is not a Module")
else
    include("LookUpTable.jl")
end
using .LookUpTable

# Convert ReliquaryAffixExcelConfigData.json into Matrix
# Substat of Artifacts
function ReliAffixECDConverter()

    # Read form JSON File
    dict = Dict()
    open("ReliquaryAffixExcelConfigData.json", "r") do f
        dict=JSON.parse(f)  # parse and transform data
    end

    # Initialize with empty Array
    # Rank, Substats, Sequence
    dataSet=fill(0f0,5,length(LookUpTable.substatProject),4);

    # fill in the data
    for i in 1:length(dict)

        # read ID for identification
        this_id=string(dict[i]["Id"]);
        this_rank=parse(Int64,this_id[1]);
        this_seq=parse(Int64,this_id[6]);
        # ID e.g. 501244 first digit is Rank, last digit is The Sequence Number in same substat

        # Skip the 9xxxxxx ID
        if !(1<=this_rank<=5)
            continue;
        end

        # Skip the Substat not covered in LookUpTable.substatProject
        if !SubstatHasKey(dict[i]["PropType"])
            continue;
        end

        # Sort Data
        dataSet[ this_rank, SubstatLookUp(dict[i]["PropType"]), this_seq ] = dict[i]["PropValue"];
        end
    return dataSet
end

# Convert ReliquaryLevelExcelConfigData.json into Matrix
# Mainstat of Artifacts
function ReliLevelECDConverter()
    dict = Dict()

    # Read form JSON File
    open("ReliquaryLevelExcelConfigData.json", "r") do f
        dict=JSON.parse(f)  # parse and transform data
    end

    # Rank, Substats, Level
    dataSet=fill(0f0,5,length(LookUpTable.mainstatProject),21);

    for i in 1:length(dict)

        # Skip the Data with no rank
        if !(haskey(dict[i],"Rank"))
            continue;
        end

        # Sort Data
        for j in 1:length(dict[i]["AddProps"])
            # Skip the Mainstat not covered in LookUpTable.mainstatProject
            if !MainstatHasKey(dict[i]["AddProps"][j]["PropType"])
                continue;
            end
            dataSet[ dict[i]["Rank"], MainstatLookUp(dict[i]["AddProps"][j]["PropType"]), dict[i]["Level"] ] = dict[i]["AddProps"][j]["Value"];
        end
    end
    return dataSet
end

end


text:
"""
    wilson_rst(g, distmx=weights(g); seed=-1, rng=GLOBAL_RNG)

Randomly sample a spanning tree from an undirected connected graph g using
[Wilson's algorithm](https://en.wikipedia.org/wiki/Loop-erased_random_walk).
The tree will be sampled with probability measure proportional to the product
of the edge weights (given by distmx).

The exact probability of the tree may be determined by computing the
normalization constant of the distribution via [Kirchoff's](https://en.wikipedia.org/wiki/Kirchhoff%27s_theorem)
or [Tutte's Matrix Tree](https://personalpages.manchester.ac.uk/staff/mark.muldoon/Teaching/DiscreteMaths/LectureNotes/MatrixTreeProof.pdf)
theorem, which is equivalent to the determinant of any minor of the Laplacian matrix.

e.g.

tree = wilson_rst(g)
probability_of_tree = prod([weights(g)[src(e), dst(e)] for e in tree])
probability_of_tree /= det(laplacian_matrix(g)[2:nv(g), 2:nv(g)]))
"""
function wilson_rst end
@traitfn function wilson_rst(g::AG::(!IsDirected),
    distmx::AbstractMatrix{T}=weights(g);
    seed::Int=-1,
    rng::AbstractRNG=GLOBAL_RNG
) where {T <: Real, U, AG <: AbstractGraph{U}}

    if seed >= 0
        rng = getRNG(seed)
    end

    start1 = rand(rng, 1:nv(g))
    start2 = rand(rng, 1:nv(g))

    walk = loop_erased_randomwalk(g, start1, distmx=distmx, f=[start2], rng=rng)

    tree = SimpleGraph(nv(g))
    for i = 1:length(walk)-1
        add_edge!(tree, walk[i], walk[i+1])
    end

    visited_vertices = Set(walk)
    unvisited_vertices = setdiff(Set([i for i = 1:nv(g)]), visited_vertices)

    while length(unvisited_vertices) > 0
        v = rand(rng, unvisited_vertices)
        walk = loop_erased_randomwalk(g, v, distmx=distmx, f=visited_vertices,
                                      rng=rng)

        for i = 1:length(walk)-1
            add_edge!(tree, walk[i], walk[i+1])
        end
        walk_set = Set(walk)
        union!(visited_vertices, walk_set)
        unvisited_vertices = setdiff(unvisited_vertices, walk_set)
    end
    return [e for e in edges(tree)]
end


text:
using Oceananigans: AbstractGrid
using Oceananigans.Grids: RegularCartesianGrid

@inline δx_c2f(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i,   j, k] - f[i-1, j, k]
@inline δx_f2c(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i+1, j, k] - f[i,   j, k]
@inline δx_e2f(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i+1, j, k] - f[i,   j, k]
@inline δx_f2e(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i,   j, k] - f[i-1, j, k]

@inline δy_c2f(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j,   k] - f[i, j-1, k]
@inline δy_f2c(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j+1, k] - f[i, j,   k]
@inline δy_e2f(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j+1, k] - f[i, j,   k]
@inline δy_f2e(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j,   k] - f[i, j-1, k]

@inline δz_c2f(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j,   k] - f[i, j, k-1]
@inline δz_f2c(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j, k+1] - f[i, j,   k]
@inline δz_e2f(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j, k+1] - f[i, j,   k]
@inline δz_f2e(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j,   k] - f[i, j, k-1]

@inline avgx_c2f(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i,   j, k] + f[i-1,  j, k])
@inline avgx_f2c(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i+1, j, k] + f[i,    j, k])
@inline avgx_f2e(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i,   j, k] + f[i-1,  j, k])

@inline avgy_c2f(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i,   j, k] + f[i,  j-1, k])
@inline avgy_f2c(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i, j+1, k] + f[i,    j, k])
@inline avgy_f2e(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i,   j, k] + f[i,  j-1, k])

@inline avgz_c2f(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i, j,   k] + f[i, j, k-1])
@inline avgz_f2c(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i, j, k+1] + f[i, j,   k])
@inline avgz_f2e(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i, j,   k] + f[i, j, k-1])

@inline ∂x_p(i, j, k, grid::RegularCartesianGrid, p) = δx_c2f(grid, p, i, j, k) / grid.Δx
@inline ∂y_p(i, j, k, grid::RegularCartesianGrid, p) = δy_c2f(grid, p, i, j, k) / grid.Δy
@inline ∂z_p(i, j, k, grid::RegularCartesianGrid, p) = δz_c2f(grid, p, i, j, k) / grid.Δz

@inline div_f2c(g::RegularCartesianGrid, fx, fy, fz, i, j, k) =
    (δx_f2c(g, fx, i, j, k) / g.Δx) + (δy_f2c(g, fy, i, j, k) / g.Δy) + (δz_f2c(g, fz, i, j, k) / g.Δz)

@inline div_c2f(g::RegularCartesianGrid, fx, fy, fz, i, j, k) =
    (δx_c2f(g, fx, i, j, k) / g.Δx) + (δy_c2f(g, fy, i, j, k) / g.Δy) + (δz_c2f(g, fz, i, j, k) / g.Δz)

@inline function δx_f2c_ab̄ˣ(g::RegularCartesianGrid, a, b, i, j, k)
    @inbounds (a[i+1, j, k] * avgx_c2f(g, b, i+1, j, k) -
               a[i,   j, k] * avgx_c2f(g, b, i,   j, k))
end

@inline function δy_f2c_ab̄ʸ(g::RegularCartesianGrid, a, b, i, j, k)
    @inbounds (a[i, j+1, k] * avgy_c2f(g, b, i, j+1, k) -
               a[i,   j, k] * avgy_c2f(g, b, i, j,   k))
end

@inline function δz_f2c_ab̄ᶻ(g::RegularCartesianGrid, a, b, i, j, k)
    @inbounds (a[i, j, k+1] * avgz_c2f(g, b, i, j, k+1) -
               a[i, j,   k] * avgz_c2f(g, b, i, j,   k))
end

@inline div_flux(g::RegularCartesianGrid, u, v, w, Q, i, j, k) =
    (δx_f2c_ab̄ˣ(g, u, Q, i, j, k) / g.Δx) + (δy_f2c_ab̄ʸ(g, v, Q, i, j, k) / g.Δy) + (δz_f2c_ab̄ᶻ(g, w, Q, i, j, k) / g.Δz)

@inline δx_c2f_ūˣūˣ(g::RegularCartesianGrid, u, i, j, k) =
    avgx_f2c(g, u, i, j, k)^2 - avgx_f2c(g, u, i-1, j, k)^2

@inline function δy_e2f_v̄ˣūʸ(g::RegularCartesianGrid, u, v, i, j, k)
    avgx_f2e(g, v, i, j+1, k) * avgy_f2e(g, u, i, j+1, k) -
    avgx_f2e(g, v, i,   j, k) * avgy_f2e(g, u, i,   j, k)
end

@inline function δz_e2f_w̄ˣūᶻ(g::RegularCartesianGrid, u, w, i, j, k)
    @inbounds avgx_f2e(g, w, i, j, k+1) * avgz_f2e(g, u, i, j, k+1) -
              avgx_f2e(g, w, i, j,   k) * avgz_f2e(g, u, i, j,   k)
end

@inline u∇u(g::RegularCartesianGrid, u, v, w, i, j, k) =
    (δx_c2f_ūˣūˣ(g, u, i, j, k) / g.Δx) + (δy_e2f_v̄ˣūʸ(g, u, v, i, j, k) / g.Δy) + (δz_e2f_w̄ˣūᶻ(g, u, w, i, j, k) / g.Δz)

@inline function δx_e2f_ūʸv̄ˣ(g::RegularCartesianGrid, u, v, i, j, k)
    avgy_f2e(g, u, i+1, j, k) * avgx_f2e(g, v, i+1, j, k) -
    avgy_f2e(g, u, i,   j, k) * avgx_f2e(g, v, i,   j, k)
end

@inline δy_c2f_v̄ʸv̄ʸ(g::RegularCartesianGrid, v, i, j, k) =
    avgy_f2c(g, v, i, j, k)^2 - avgy_f2c(g, v, i, j-1, k)^2

@inline function δz_e2f_w̄ʸv̄ᶻ(g::RegularCartesianGrid, v, w, i, j, k)
    @inbounds avgy_f2e(g, w, i, j, k+1) * avgz_f2e(g, v, i, j, k+1) -
              avgy_f2e(g, w, i, j,   k) * avgz_f2e(g, v, i, j,   k)
end

@inline u∇v(g::RegularCartesianGrid, u, v, w, i, j, k) =
    (δx_e2f_ūʸv̄ˣ(g, u, v, i, j, k) / g.Δx) + (δy_c2f_v̄ʸv̄ʸ(g, v, i, j, k) / g.Δy) + (δz_e2f_w̄ʸv̄ᶻ(g, v, w, i, j, k) / g.Δz)

@inline function δx_e2f_ūᶻw̄ˣ(g::RegularCartesianGrid, u, w, i, j, k)
    avgz_f2e(g, u, i+1, j, k) * avgx_f2e(g, w, i+1, j, k) -
    avgz_f2e(g, u, i,   j, k) * avgx_f2e(g, w, i,   j, k)
end

@inline function δy_e2f_v̄ᶻw̄ʸ(g::RegularCartesianGrid, v, w, i, j, k)
    avgz_f2e(g, v, i, j+1, k) * avgy_f2e(g, w, i, j+1, k) -
    avgz_f2e(g, v, i,   j, k) * avgy_f2e(g, w, i,   j, k)
end

@inline δz_c2f_w̄ᶻw̄ᶻ(g::RegularCartesianGrid, w, i, j, k) =
    avgz_f2c(g, w, i, j, k)^2 - avgz_f2c(g, w, i, j, k-1)^2

@inline u∇w(g::RegularCartesianGrid, u, v, w, i, j, k) =
    (δx_e2f_ūᶻw̄ˣ(g, u, w, i, j, k) / g.Δx) + (δy_e2f_v̄ᶻw̄ʸ(g, v, w, i, j, k) / g.Δy) + (δz_c2f_w̄ᶻw̄ᶻ(g, w, i, j, k) / g.Δz)

@inline δx²_c2f2c(g::RegularCartesianGrid, f, i, j, k) = δx_c2f(g, f, i+1, j, k) - δx_c2f(g, f, i, j, k)
@inline δy²_c2f2c(g::RegularCartesianGrid, f, i, j, k) = δy_c2f(g, f, i, j+1, k) - δy_c2f(g, f, i, j, k)
@inline δz²_c2f2c(g::RegularCartesianGrid, f, i, j, k) = δz_c2f(g, f, i, j, k+1) - δz_c2f(g, f, i, j, k)

@inline κ∇²(g::RegularCartesianGrid, Q, κh, κv, i, j, k) =
    ((κh/g.Δx^2) * δx²_c2f2c(g, Q, i, j, k)) + ((κh/g.Δy^2) * δy²_c2f2c(g, Q, i, j, k)) + ((κv/g.Δz^2) * δz²_c2f2c(g, Q, i, j, k))

@inline δx²_f2c2f(g::RegularCartesianGrid, f, i, j, k) = δx_f2c(g, f, i, j, k) - δx_f2c(g, f, i-1, j, k)
@inline δy²_f2c2f(g::RegularCartesianGrid, f, i, j, k) = δy_f2c(g, f, i, j, k) - δy_f2c(g, f, i, j-1, k)
@inline δz²_f2c2f(g::RegularCartesianGrid, f, i, j, k) = δz_f2c(g, f, i, j, k) - δz_f2c(g, f, i, j, k-1)

@inline δx²_f2e2f(g::RegularCartesianGrid, f, i, j, k) = δx_f2e(g, f, i+1, j, k) - δx_f2e(g, f, i, j, k)
@inline δy²_f2e2f(g::RegularCartesianGrid, f, i, j, k) = δy_f2e(g, f, i, j+1, k) - δy_f2e(g, f, i, j, k)
@inline δz²_f2e2f(g::RegularCartesianGrid, f, i, j, k) = δz_f2e(g, f, i, j, k+1) - δz_f2e(g, f, i, j, k)

@inline 𝜈∇²u(g::RegularCartesianGrid, u, 𝜈h, 𝜈v, i, j, k) =
    ((𝜈h/g.Δx^2) * δx²_f2c2f(g, u, i, j, k)) + ((𝜈h/g.Δy^2) * δy²_f2e2f(g, u, i, j, k)) + ((𝜈v/g.Δz^2) * δz²_f2e2f(g, u, i, j, k))

@inline 𝜈∇²v(g::RegularCartesianGrid, v, 𝜈h, 𝜈v, i, j, k) =
    ((𝜈h/g.Δx^2) * δx²_f2e2f(g, v, i, j, k)) + ((𝜈h/g.Δy^2) * δy²_f2c2f(g, v, i, j, k)) + ((𝜈v/g.Δz^2) * δz²_f2e2f(g, v, i, j, k))

@inline 𝜈∇²w(g::RegularCartesianGrid, w, 𝜈h, 𝜈v, i, j, k) =
    ((𝜈h/g.Δx^2) * δx²_f2e2f(g, w, i, j, k)) + ((𝜈h/g.Δy^2) * δy²_f2e2f(g, w, i, j, k)) + ((𝜈v/g.Δz^2) * δz²_f2c2f(g, w, i, j, k))

@inline ∇²(g::RegularCartesianGrid, f, i, j, k) =
    (δx²_c2f2c(g, f, i, j, k) / g.Δx^2) + (δy²_c2f2c(g, f, i, j, k) / g.Δy^2) + (δz²_c2f2c(g, f, i, j, k) / g.Δz^2)

@inline ∇h_u(i, j, k, grid, u, v) = δx_f2c(grid, u, i, j, k) / grid.Δx + δy_f2c(grid, v, i, j, k) / grid.Δy


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE

import Base: length, size

function length{T<:AbstractBasis}(B::T)
    return length(T)
end

function size{T<:AbstractBasis}(B::T)
    return size(T)
end

function eval_basis!{T<:AbstractBasis}(B::T, N, xi)
    return eval_basis!(T, N, xi)
end

function eval_dbasis!{T<:AbstractBasis}(B::T, dN, xi)
    return eval_dbasis!(T, dN, xi)
end

function eval_basis(B, xi)
    N = zeros(1, length(B))
    eval_basis!(B, N, xi)
    return N
end

function eval_dbasis(B, xi)
    dN = zeros(size(B)...)
    eval_dbasis!(B, dN, xi)
    return dN
end

"""
    interpolate(B, T, xi)

Given basis B, interpolate T at xi.

# Example
```jldoctest
B = Quad4()
X = ((0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0))
T = (1.0, 2.0, 3.0, 4.0)
interpolate(B, T, (0.0, 0.0))

# output

2.5

```
"""
function interpolate(B, T, xi)
    Ti = sum(b*t for (b, t) in zip(eval_basis(B, xi), T))
    return Ti
end

"""
    jacobian(B, X, xi)

Given basis B, calculate jacobian at xi.

# Example
```jldoctest
B = Quad4()
X = ([0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0])
jacobian(B, X, (0.0, 0.0))

# output

2×2 Array{Float64,2}:
 0.5  0.0
 0.0  0.5

```
"""
function jacobian(B, X, xi)
    dB = eval_dbasis(B, xi)
    dim1, nbasis = size(B)
    dim2 = length(first(X))
    J = zeros(dim1, dim2)
    for i=1:dim1
        for j=1:dim2
            for k=1:nbasis
                J[i,j] += dB[i,k]*X[k][j]
            end
        end
    end
    return J
end

"""
    grad(B, X, xi)

Given basis B, calculate gradient dB/dX at xi.

# Example
```jldoctest
B = Quad4()
X = ([0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0])
grad(B, X, (0.0, 0.0))

# output

2×4 Array{Float64,2}:
 -0.5   0.5  0.5  -0.5
 -0.5  -0.5  0.5   0.5

```
"""
function grad(B, X, xi)
    J = jacobian(B, X, xi)
    dB = eval_dbasis(B, xi)
    G = inv(J) * dB
    return G
end

"""
    grad(B, T, X, xi)

Calculate gradient of `T` with respect to `X` in point `xi` using basis `B`.

# Example
```jldoctest
B = Quad4()
X = ([0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0])
u = ([0.0, 0.0], [1.0, -1.0], [2.0, 3.0], [0.0, 0.0])
grad(B, u, X, (0.0, 0.0))

# output

2×2 Array{Float64,2}:
 1.5  0.5
 1.0  2.0

```
"""
function grad(B, T, X, xi)
    G = grad(B, X, xi)
    nbasis = length(B)
    dTdX = sum(kron(G[:,i], T[i]') for i=1:nbasis)
    return dTdX'
end

"""
Data type for fast FEM.
"""
type BasisInfo{B<:AbstractBasis,T}
    N::Matrix{T}
    dN::Matrix{T}
    grad::Matrix{T}
    J::Matrix{T}
    invJ::Matrix{T}
    detJ::T
end

function length{T<:AbstractBasis}(B::BasisInfo{T})
    return length(T)
end

function size{T<:AbstractBasis}(B::BasisInfo{T})
    return size(T)
end

"""
Initialization of data type `BasisInfo`.

# Examples

```jldoctest

BasisInfo(Tri3)

# output

FEMBasis.BasisInfo{FEMBasis.Tri3,Float64}([0.0 0.0 0.0], [0.0 0.0 0.0; 0.0 0.0 0.0], [0.0 0.0 0.0; 0.0 0.0 0.0], [0.0 0.0; 0.0 0.0], [0.0 0.0; 0.0 0.0], 0.0)

```

"""
function BasisInfo{B<:AbstractBasis}(::Type{B}, T=Float64)
    dim, nbasis = size(B)
    N = zeros(T, 1, nbasis)
    dN = zeros(T, dim, nbasis)
    grad = zeros(T, dim, nbasis)
    J = zeros(T, dim, dim)
    invJ = zeros(T, dim, dim)
    detJ = zero(T)
    return BasisInfo{B,T}(N, dN, grad, J, invJ, detJ)
end

"""
Evaluate basis, gradient and so on for some point `xi`.

# Examples

```jldoctest

b = BasisInfo(Quad4)
X = ((0.0,0.0), (1.0,0.0), (1.0,1.0), (0.0,1.0))
xi = (0.0, 0.0)
eval_basis!(b, X, xi)

# output

FEMBasis.BasisInfo{FEMBasis.Quad4,Float64}([0.25 0.25 0.25 0.25], [-0.25 0.25 0.25 -0.25; -0.25 -0.25 0.25 0.25], [-0.5 0.5 0.5 -0.5; -0.5 -0.5 0.5 0.5], [0.5 0.0; 0.0 0.5], [2.0 -0.0; -0.0 2.0], 0.25)

```
"""
function eval_basis!{B}(bi::BasisInfo{B}, X, xi)

    # evaluate basis and derivatives
    eval_basis!(B, bi.N, xi)
    eval_dbasis!(B, bi.dN, xi)

    dim1, nbasis = size(B)
    dim2 = length(first(X))
    dims = (dim1, dim2)
    if size(bi.J) != dims
        bi.J = zeros(dim1, dim2)
    else
        fill!(bi.J, 0.0)
    end

    # calculate Jacobian

    for i=1:dim1
        for j=1:dim2
            for k=1:nbasis
                bi.J[i,j] += bi.dN[i,k]*X[k][j]
            end
        end
    end

    # calculate determinant of Jacobian + gradient operator

    if dims == (3, 3)
        a, b, c, d, e, f, g, h, i = bi.J
        bi.detJ = a*(e*i-f*h) + b*(f*g-d*i) + c*(d*h-e*g)
        bi.invJ[1] = 1.0 / bi.detJ * (e*i - f*h)
        bi.invJ[2] = 1.0 / bi.detJ * (c*h - b*i)
        bi.invJ[3] = 1.0 / bi.detJ * (b*f - c*e)
        bi.invJ[4] = 1.0 / bi.detJ * (f*g - d*i)
        bi.invJ[5] = 1.0 / bi.detJ * (a*i - c*g)
        bi.invJ[6] = 1.0 / bi.detJ * (c*d - a*f)
        bi.invJ[7] = 1.0 / bi.detJ * (d*h - e*g)
        bi.invJ[8] = 1.0 / bi.detJ * (b*g - a*h)
        bi.invJ[9] = 1.0 / bi.detJ * (a*e - b*d)
        A_mul_B!(bi.grad, bi.invJ, bi.dN)
    elseif dims == (2, 2)
        a, b, c, d = bi.J
        bi.detJ = a*d - b*c
        bi.invJ[1] = 1.0 / bi.detJ * d
        bi.invJ[2] = 1.0 / bi.detJ * -b
        bi.invJ[3] = 1.0 / bi.detJ * -c
        bi.invJ[4] = 1.0 / bi.detJ * a
        A_mul_B!(bi.grad, bi.invJ, bi.dN)
    elseif dims == (1, 1)
        bi.detJ = bi.J[1]
        bi.invJ[1] = 1.0 / bi.detJ
        bi.grad = bi.invJ * bi.dN
    elseif dim1 == 1 # curve
        bi.detJ = norm(bi.J)
    elseif dim1 == 2 # manifold
        bi.detJ = norm(cross(bi.J[1,:], bi.J[2,:]))
    end

    return bi
end

"""
    grad!(bi, gradu, u)

Evalute gradient ∂u/∂X and store result to matrix `gradu`. It is assumed
that `eval_basis!` has been already run to `bi` so it already contains
all necessary matrices evaluated with some `X` and `xi`.

# Example

First setup and evaluate basis using `eval_basis!`:
```jldoctest ex1
B = BasisInfo(Quad4)
X = ((0.0,0.0), (1.0,0.0), (1.0,1.0), (0.0,1.0))
xi = (0.0, 0.0)
eval_basis!(B, X, xi)

# output

FEMBasis.BasisInfo{FEMBasis.Quad4,Float64}([0.25 0.25 0.25 0.25], [-0.25 0.25 0.25 -0.25; -0.25 -0.25 0.25 0.25], [-0.5 0.5 0.5 -0.5; -0.5 -0.5 0.5 0.5], [0.5 0.0; 0.0 0.5], [2.0 -0.0; -0.0 2.0], 0.25)

```

Next, calculate gradient of `u`:
```jldoctest ex1
u = ((0.0, 0.0), (1.0, -1.0), (2.0, 3.0), (0.0, 0.0))
gradu = zeros(2, 2)
grad!(B, gradu, u)

# output

2×2 Array{Float64,2}:
 1.5  0.5
 1.0  2.0
```
"""
function grad!{B}(bi::BasisInfo{B}, gradu, u)
    dim, nbasis = size(B)
    fill!(gradu, 0.0)
    for i=1:dim
        for j=1:dim
            for k=1:nbasis
                gradu[i,j] += bi.grad[j,k]*u[k][i]
            end
        end
    end
    return gradu
end


text:
@testset "961.n-repeated-element-in-size-2n-array.jl" begin
    @test repeated_n_times([5,1,5,2,5,3,5,4]) == 5
    @test repeated_n_times([2,1,2,5,3,2]) == 2
    @test repeated_n_times([1,2,3,3]) == 3
end


text:
using PowerSystems
using NLsolve
const PSY = PowerSystems

include(joinpath(dirname(@__FILE__), "dynamic_test_data.jl"))
include(joinpath(dirname(@__FILE__), "data_utils.jl"))
############### Data Network ########################
omib_file_dir = joinpath(dirname(@__FILE__), "OMIB_DARCO_PSR.raw")
omib_sys = System(omib_file_dir, runchecks = false)
add_source_to_ref(omib_sys)

############### Data Dynamic devices ########################
function inv_gfoll(static_device)
    return PSY.DynamicInverter(
        get_name(static_device),
        1.0, #ω_ref
        converter_low_power(), #converter
        outer_control_gfoll(), #outercontrol
        current_mode_inner(), #inner_control
        dc_source_lv(),
        reduced_pll(),
        filt_gfoll(),
    ) #pss
end

#Attach dynamic generator. Currently use PSS/e format based on bus #.
device = [g for g in get_components(Generator, omib_sys)][1]
case_inv = inv_gfoll(device)
add_component!(omib_sys, case_inv, device)


text:
using DynamicMovementPrimitives
using Test
import DynamicMovementPrimitives: get_sched_sig, get_centers_linear, kernel_matrix, _1

# Setup ============================================
Nbasis = 20
αzC = 25.
αzPT = 10.
αx = 1.#αz/3
# bz = 1
# close("all")

optsC = DMPopts(Nbasis,αx,αzC, :canonical)
optsT = DMPopts(Nbasis,αx,αzPT, :time)
optsP = DMPopts(Nbasis,αx,αzPT, :position)
@test optsC.βz == αzC/4

y = [LinRange(0,1,100); ones(10)]
T = length(y)
t = LinRange(0,T,T)
τ = t[end]/3
h = t[2]-t[1]
y = [y 0.5y]
ẏ = centraldiff(y) ./ h
ÿ = centraldiff(ẏ) ./ h
g = y[end,:][:]

# Tests canonical ============================================

x = get_sched_sig(:canonical,αx,τ,t,y,g)
c,σ2    = get_centers_linear(Nbasis,x)
Ψ = kernel_matrix(x,c,σ2)
@test size(x) == (T,)
@test size(c) == (Nbasis,1)
@test size(σ2) == (Nbasis,1)
@test size(Ψ) == (T,Nbasis)

dmp = fit(y,ẏ,ÿ,t,optsC)
f = force(dmp,x)
@test size(f) == (T,2)
@test all(isfinite.(f))
f = force(dmp,x[1],1)
@test isa(f,Number)
@test isfinite(f)

a = acceleration(dmp,y[1],ẏ[1],x[1],g[1],1)
@test isa(a,Number)
@test isfinite(a)
a = acceleration(dmp,_1(y),_1(ẏ),x[1],g)
@test size(a) == (2,)
a = acceleration(dmp,y,ẏ,x,g)
@test size(a) == (T,2)
@test all(isfinite.(a))

# Tests time     ============================================
# x = get_sched_sig(:time,αx,τ,t,y,g)
# c,σ2    = get_centers_linear(Nbasis,x)
# Ψ = kernel_matrix(x,c,σ2)
# @test size(x) == (T,)
# @test size(c) == (Nbasis,1)
# @test size(σ2) == (Nbasis,1)
# @test size(Ψ) == (T,Nbasis)
#
# dmp = fit(y,ẏ,ÿ,t,optsT)
# f = force(dmp,x)
# @test size(f) == (T,2)
# @test all(isfinite(f))
# f = force(dmp,x[1],1)
# @test isa(f,Number)
# @test isfinite(f)

# Tests position ============================================
# x = get_sched_sig(:position,αx,τ,t,y,g)
# c,σ2    = get_centers_linear(Nbasis,x)
# Ψ = kernel_matrix(x[:,1],c[:,1],σ2[:,1])
# @test size(x) == (T,2)
# @test size(c) == (Nbasis,2)
# @test size(σ2) == (Nbasis,2)
# @test size(Ψ) == (T,Nbasis)
#
# dmp = fit(y,ẏ,ÿ,t,optsP)
# f = force(dmp,x)
# @test size(f) == (T,2)
# @test all(isfinite(f))
# f = force(dmp,x[1],1)
# @test isa(f,Number)
# @test isfinite(f)
#
# a = acceleration(dmp,y[1],ẏ[1],x[1],g[1],1)
# @test isa(a,Number)
# @test isfinite(a)
# a = acceleration(dmp,_1(y),_1(ẏ),x[1,:][:],g)
# @test size(a) == (2,)
# a = acceleration(dmp,y,ẏ,x,g)
# @test size(a) == (T,2)
# @test all(isfinite(a))

dmp = fit(y,ẏ,ÿ,t,optsC)
tout,youtC,ẏout,xout = solve(dmp)
@test tout == t
@test abs.(youtC .- y) |> sum < 2
@test abs.(ẏout .- ẏ) |> sum < 0.3
# plotdmp(dmp)

# dmp = fit(y,ẏ,ÿ,t,optsT)
# tout,youtT,ẏout,xout = solve(dmp)
# @test abs(youtT .- y) |> sum < 2
# @test abs(ẏout .- ẏ) |> sum < 0.3
# # plotdmp(dmp)

# dmp = fit(y,ẏ,ÿ,t,optsP)
# tout,youtP,ẏout,xout = solve(dmp)
# @test abs(youtP .- y) |> sum < 6
# @test abs(ẏout .- ẏ) |> sum < 0.3
# # plotdmp(dmp)
# # plotdmpphase(dmp)

# Test twolink
include("../src/two_link.jl")
using Main.TwoLink
p,pd,pdd = traj(0.,1.,0:100)
@test p[1] == 0
@test p[end] == 1
@test pd[1] == 0
@test pd[end] == 0
@test length(p) == length(pd) == length(pdd) == 101

cpoints = [0.5 -0.5;
        0.5 0;
        1.5 0;
        1.5 -0.5]

ctraj = connect_points(cpoints,20)[1]
jtraj = inverse_kin(ctraj,:up)
ctraj2 = forward_kin(jtraj)

@test ctraj ≈ ctraj2

# TEST 2 dof controller

using DynamicMovementPrimitives
Nbasis  = 15
αz      = 25.
αx      = 1.
opts    = DMPopts(Nbasis,αx,αz)

y       = [zeros(10);LinRange(0,2,1000); 2ones(100)]
T       = length(y)
t       = LinRange(0,10,T)
h       = t[2]-t[1] # Sample interval
y       = [y 0.5y]
ẏ       = centraldiff(y) ./h # Differentiate position to get velocity
ÿ       = centraldiff(ẏ) ./h
dmp     = fit(y,ẏ,ÿ,t,opts)
dmp2opts = DMP2dofopts(kp = 25,kv = 10,kc = 10_000,αe = 5)
dmp2    = DMP2dof(dmp, dmp2opts) # Upgrade dmp to 2DOF version

t,yc,ẏc,x,ya,ẏa,e = solve(dmp2,t)

import OrdinaryDiffEq

condition(u,t,integrator) = 2.5 <= t < 4
affect!(integrator) = (integrator.u[3] = integrator.uprev[3])
cb = OrdinaryDiffEq.DiscreteCallback(condition,affect!)
t,yc,ẏc,x,ya,ẏa,e = solve(dmp2,t, solver=OrdinaryDiffEq.Euler(), callback=cb)

# plot(t,ẏc, lab="\$ẏ_c\$", c=:red, l=(:dash, 3), layout=(2,2), subplot=1)
# plot!(t,yc, lab="\$y_c\$", c=:red, l=(:dash, 3), subplot=2)
# plot!(t,ẏa, lab="\$ẏ_a\$", c=:blue, subplot=1)
# plot!(t,ya, lab="\$y_a\$", c=:blue, subplot=2)
# plot!(t,e, lab="\$e\$", c=:green, subplot=3)
# plot!(t,2.5 .<= t .< 4, lab="Disturbance", c=:green, subplot=4, fillrange=0)
#
# t,yc,ẏc,x,ya,ẏa,e = solve(dmp2,t)
# plot!(t,ẏc, lab="\$ẏ_u\$", c=:black, l=(:dashdot, 3), subplot=1)
# plot!(t,yc, lab="\$y_u\$", c=:black, l=(:dashdot, 3), subplot=2)


text:
# carts moving on tracks.  move until they hit a crossroad, then their behavior
# depends on some state (turn left, then straight, then right, then repeat)

test_input = """
/->-\\
|   |  /----\\
| /-+--+-\\  |
| | |  | v  |
\\-+-/  \\-+--/
  \\------/
"""

mutable struct Cart
    pos::CartesianIndex{2}
    vel::CartesianIndex{2}
    turn::Int
end

Base.isless(a::Cart, b::Cart) = isless(a.pos, b.pos)

mutable struct Board
    map::Matrix{Char}
    tmp::Matrix{Char}
    carts::Vector{Cart}
    crashes::BitVector
end

Board(map::String) = Board(reduce(hcat, collect.(split(strip(map, '\n'), '\n'))))
function Board(map::Matrix{Char})
    board = Board(map, similar(map), Cart[], falses(0))
    for (char, dir, rep) in (('>', (1,0), '-'), ('<', (-1,0), '-'),
                             ('^', (0,-1), '|'), ('v', (0,1), '|'))
        for i in findall(isequal(char), board.map)
            push!(board.carts, Cart(i, CartesianIndex(dir), 0))
            board.map[i] = rep
        end
    end
    sort!(board.carts)
    return board
end

board = Board(read("13.input", String))

# right (CW) turn: up -> right, right -> down, down -> left, left -> up
# (0,-1) -> (1,0)
# (1,0)  -> (0,1)
# (0,1)  -> (-1,0)
# (-1,0) -> (0,-1)

# x1 = -1 * y0
# y1 = x0
right(x::CartesianIndex{2}) = CartesianIndex(-1*x.I[2], x.I[1])
left(x::CartesianIndex{2}) = CartesianIndex(x.I[2], -1*x.I[1])

function turn!(cart::Cart, location::Char)
    if location == '+'
        # turn depends on state
        cart.vel =
            cart.turn == 0 ? left(cart.vel) :
            cart.turn == 2 ? right(cart.vel) :
            cart.vel
        cart.turn = (cart.turn + 1) % 3
    elseif location == '/'
        cart.vel = -1 * CartesianIndex(reverse(cart.vel.I))
    elseif location == '\\'
        cart.vel = CartesianIndex(reverse(cart.vel.I))
    end
    cart
end

function move!(board::Board)
    resize!(board.crashes, length(board.carts))
    fill!(board.crashes, false)

    for (i,cart) in enumerate(board.carts)
        turn!(cart, board.map[cart.pos])
        cart.pos += cart.vel
        # check for collisions
        for (other_i, other_cart) in enumerate(board.carts)
            board.crashes[other_i] && continue
            if i != other_i && cart.pos == other_cart.pos
                board.crashes[i] = board.crashes[other_i] = true
                println("Collision at $(cart.pos.I .- 1)")
                @show board.crashes
            end
        end
    end

    deleteat!(board.carts, board.crashes)
    sort!(board.carts)
    return board
end

function Base.show(io::IO, board::Board)
    copyto!(board.tmp, board.map)
    for cart in board.carts
        board.tmp[cart.pos] =
            cart.vel.I == (1,0) ? '>' :
            cart.vel.I == (-1,0) ? '<' :
            cart.vel.I == (0,1) ? 'v' :
            '^'
    end
    xs, ys = axes(board.tmp)
    for y in ys
        println(String(board.tmp[:, y]))
    end
end

function star1(input)
    board = Board(input)
    n_carts = length(board.carts)
    while length(board.carts) == n_carts
        move!(board)
    end
end

star1(test_input)
star1(read("13.input", String))

test_input_2 = """
/>-<\\
|   |
| /<+-\\
| | | v
\\>+</ |
  |   ^
  \\<->/
"""

function star2(input)
    board = Board(input)
    while length(board.carts) > 1
        move!(board)
    end
    return first(board.carts).pos.I .- 1
end

star2(test_input_2)
star2(read("13.input", String))


text:
const JULIALANG_REPOS = MirrorUpdater.Types.SrcDestPair[
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/JuliaLang/METADATA.jl",
        destination_repo_name = "JuliaLang-METADATA.jl",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/JuliaLang/Pkg.jl",
        destination_repo_name =
            "JuliaLang-Pkg.jl-44cfe95a-1eb2-52ea-b672-e2afdf69b78f",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/JuliaLang/julia",
        destination_repo_name = "JuliaLang-julia",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/JuliaLang/libosxunwind",
        destination_repo_name = "JuliaLang-libosxunwind",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/JuliaLang/libuv",
        destination_repo_name = "JuliaLang-libuv",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/JuliaLang/readline",
        destination_repo_name = "JuliaLang-readline",
        ),
    ]


text:
import Pkg
Pkg.instantiate() # downloads all dependencies for the current project

Pkg.add("Plots")

using Plots, DataFrames, CSV

df = CSV.read("../externals/Core.Math.Data/data/Pejcic_318.csv", copycols = true)

a = df.ATT # use something to plot;
b = df.ATV
c = df.AOP
d = df.ANN

# display(plot(histogram(a)));

p1 = (histogram(a))
p2 = (histogram(b))
p3 = (histogram(c))
p4 = (histogram(d))

plot(p1, p2, p3, p4, layout = (2, 2), legend = false)


text:
using IMPES
using Test

@testset "IMPES.jl" begin
    # Write your tests here.
end