text:
# # include("quaternions.jl")
# # using Quaternions
# export
#       trim_controls
#
# struct Quadrotor{Q,T} <: AbstractModel #RigidBody{Quat{VectorPart}}
#       n::Int
#       m::Int
#       mass::T
#       J::SMatrix{3,3,T,9}
#       Jinv::SMatrix{3,3,T,9}
#       gravity::SVector{3,T}
#       motor_dist::T
#       kf::T
#       km::T
#       info::Dict{Symbol,Any}
# end
#
# function Quadrotor(;mass=0.5,
#                    J=SMatrix{3,3}(Diagonal([0.0023, 0.0023, 0.004])),
#                    gravity=SVector(0,0,-9.81),
#                    motor_dist=0.1750,
#                    kf=1.0,
#                    km=0.0245,
#                    info=Dict{Symbol,Any}(),
#                    use_quat=false)
#
#     Quadrotor{use_quat,Float64}(13,4,mass,J,inv(J),gravity,motor_dist,kf,km,info)
# end
#
# function dynamics(quad::Quadrotor, x::AbstractVector, u::AbstractVector)
#       # Quaternion representation
#       # Modified from D. Mellinger, N. Michael, and V. Kumar,
#       # "Trajectory generation and control for precise aggressive maneuvers with quadrotors",
#       # In Proceedings of the 12th International Symposium on Experimental Robotics (ISER 2010), 2010.
#
#       ## States: X ∈ R^13; q = [s;v]
#       # x
#       # y
#       # z
#       # q0
#       # q1
#       # q2
#       # q3
#       # xdot
#       # ydot
#       # zdot
#       # omega1
#       # omega2
#       # omega3
#
#       # x = X[1:3]
#       q = normalize(UnitQuaternion(x[4],x[5],x[6],x[7]))
#       # q = Rotations.Quat(x[4],x[5],x[6],x[7])
#       # q = Quaternions.Quaternion(x[4],x[5],x[6],x[7])
#       # q = normalize(q)
#       # q = normalize(Quaternion(view(x,4:7)))
#       # q = view(x,4:7)
#       # normalize!(q)
#       v = view(x,8:10)
#       omega = x[@SVector [11,12,13]]
#
#       # Parameters
#       m = quad.mass # mass
#       J = quad.J # inertia matrix
#       Jinv = quad.Jinv # inverted inertia matrix
#       g = quad.gravity # gravity
#       L = quad.motor_dist # distance between motors
#       kf = quad.kf # 6.11*10^-8;
#       km = quad.km
#
#       w1 = u[1]
#       w2 = u[2]
#       w3 = u[3]
#       w4 = u[4]
#
#       F1 = kf*w1;
#       F2 = kf*w2;
#       F3 = kf*w3;
#       F4 = kf*w4;
#       F = @SVector [0., 0., F1+F2+F3+F4] #total rotor force in body frame
#
#       M1 = km*w1;
#       M2 = km*w2;
#       M3 = km*w3;
#       M4 = km*w4;
#       tau = @SVector [L*(F2-F4), L*(F3-F1), (M1-M2+M3-M4)] #total rotor torque in body frame
#
#       # ẋ[4:7] = 0.5*qmult(q,[0;omega]) #quaternion derivative
#       # qdot = SVector(0.5*q*Quaternion(zero(x[1]), omega))
#       qdot = kinematics(q, ω)
#       ωhat = Quaternion(zero(eltype(omega)), omega[1], omega[2], omega[3])
#       qdot = SVector(0.5*q*ωhat)
#       vdot = g + (1/m)*(q*F) #acceleration in world frame
#       omdot = Jinv*(tau - cross(omega,J*omega)) #Euler's equation: I*ω + ω x I*ω = constraint_decrease_ratio
#       @SVector [v[1], v[2], v[3], qdot[1], qdot[2], qdot[3], qdot[4], vdot[1], vdot[2], vdot[3], omdot[1], omdot[2], omdot[3]]
# end
#
#
# function state_diff(quad::Quadrotor, x, x0)
#       inds = @SVector [4,5,6,7]
#       q = x[inds]
#       q0 = x0[inds]
#       δx = x - x0
#       δq = quat_diff(q, q0)
#       δx = @SVector [δx[1], δx[2], δx[3], δq[1], δq[2], δq[3], δx[8], δx[9], δx[10], δx[11], δx[12], δx[13]]
# end
#
# state_diff(::Quadrotor{false}, x, x0) = x - x0
#
# function quat_diff(q2::SVector{4,T1}, q1::SVector{4,T2}) where {T1,T2}
#     # equivalent to q2 - q1
#     # same as inv(q1)*q2
#     vec = @SVector [2,3,4]
#     s1,v1 = q1[1],-q1[vec]
#     s2,v2 = q2[1], q2[vec]  # needs an inverse
#     # this is q1*q2
#     s1*v2 + s2*v1 + v1 × v2
# end
#
# function quat_diff_jacobian(q::SVector{4,T}) where T
#     w,x,y,z = q
#     x,y,z = -x,-y,-z  # invert q
#     @SMatrix [x  w -z  y;
#               y  z  w -x;
#               z -y  x  w];
# end
#
# function state_diff_jacobian(quad::Quadrotor, x0::SVector{N,T}) where {N,T}
#     inds = @SVector [4,5,6,7]
#     q0 = x0[inds]
#     G = quat_diff_jacobian(q0)
#     I1 = @SMatrix [1 0 0 0 0 0 0 0 0 0 0 0 0;
#                    0 1 0 0 0 0 0 0 0 0 0 0 0;
#                    0 0 1 0 0 0 0 0 0 0 0 0 0;
#                    0 0 0 G[1] G[4] G[7] G[10] 0 0 0 0 0 0;
#                    0 0 0 G[2] G[5] G[8] G[11] 0 0 0 0 0 0;
#                    0 0 0 G[3] G[6] G[9] G[12] 0 0 0 0 0 0;
#                    0 0 0 0 0 0 0 1 0 0 0 0 0;
#                    0 0 0 0 0 0 0 0 1 0 0 0 0;
#                    0 0 0 0 0 0 0 0 0 1 0 0 0;
#                    0 0 0 0 0 0 0 0 0 0 1 0 0;
#                    0 0 0 0 0 0 0 0 0 0 0 1 0;
#                    0 0 0 0 0 0 0 0 0 0 0 0 1.]
#    return I1'
# end
#
# state_diff_jacobian(::Quadrotor{false}, x0::SVector) = I
#
# function TrajectoryOptimization.state_diff_jacobian!(G, model::Quadrotor{true}, Z::Traj)
#     for k in eachindex(Z)
#         G[k] = state_diff_jacobian(model, state(Z[k]))
#     end
# end
#
# @inline state_diff_size(::Quadrotor{false}) = 13
# @inline state_diff_size(::Quadrotor{true}) = 12

struct Quadrotor2{B,R} <: TestRigidBody{B,R}
    n::Int
    m::Int
    mass::Float64
    J::Diagonal{Float64,SVector{3,Float64}}
    Jinv::Diagonal{Float64,SVector{3,Float64}}
    gravity::SVector{3,Float64}
    motor_dist::Float64
    kf::Float64
    km::Float64
    info::Dict{Symbol,Any}
end

Base.size(::Quadrotor2{B,<:UnitQuaternion}) where B = 13,4
Base.size(::Quadrotor2) = 12,4

function Quadrotor2{R}(;use_rot=R<:UnitQuaternion,
        mass=0.5,
        J=Diagonal(@SVector [0.0023, 0.0023, 0.004]),
        gravity=SVector(0,0,-9.81),
        motor_dist=0.1750,
        kf=1.0,
        km=0.0245,
        info=Dict{Symbol,Any}()) where R
    if R <: UnitQuaternion{T,IdentityMap} where T
        if use_rot == true
            use_rot = false
        end
    end
    Quadrotor2{use_rot,R}(13,4,mass,J,inv(J),gravity,motor_dist,kf,km,info)
end

(::Type{Quadrotor2})(;kwargs...) = Quadrotor2{UnitQuaternion{Float64,CayleyMap}}(;kwargs...)

function trim_controls(model::Quadrotor2)
    @SVector fill(-model.gravity[3]*model.mass/4.0, size(model)[2])
end

function forces(model::Quadrotor2, x, u)
    q = orientation(model, x)
    kf = model.kf
    g = model.gravity
    m = model.mass

    w1 = u[1]
    w2 = u[2]
    w3 = u[3]
    w4 = u[4]

    F1 = max(0,kf*w1);
    F2 = max(0,kf*w2);
    F3 = max(0,kf*w3);
    F4 = max(0,kf*w4);
    F = @SVector [0., 0., F1+F2+F3+F4] #total rotor force in body frame

    m*g + q*F # forces in world frame
end

function moments(model::Quadrotor2, x, u)

    kf, km = model.kf, model.km
    L = model.motor_dist

    w1 = u[1]
    w2 = u[2]
    w3 = u[3]
    w4 = u[4]

    F1 = max(0,kf*w1);
    F2 = max(0,kf*w2);
    F3 = max(0,kf*w3);
    F4 = max(0,kf*w4);

    M1 = km*w1;
    M2 = km*w2;
    M3 = km*w3;
    M4 = km*w4;
    tau = @SVector [L*(F2-F4), L*(F3-F1), (M1-M2+M3-M4)] #total rotor torque in body frame
end

function wrenches(model::Quadrotor2, x::SVector, u::SVector)
    F = forces(model, x, u)
    M = moments(model, x, u)
    return F,M

    q = orientation(model, x)
    C = forceMatrix(model)
    mass, g = model.mass, model.gravity

    # Calculate force and moments
    w = max.(u, 0)  # keep forces positive
    fM = forceMatrix(model)*w
    f = fM[1]
    M = @SVector [fM[2], fM[3], fM[4]]
    e3 = @SVector [0,0,1]
    F = mass*g - q*(f*e3)
    return F,M
end

function forceMatrix(model::Quadrotor2)
    kf, km = model.kf, model.km
    L = model.motor_dist
    @SMatrix [
       -kf  -kf  -kf  -kf;
        0    L*kf 0   -L*kf;
       -L*kf 0    L*kf 0;
        km  -km   km  -km;
    ]
end

inertia(model::Quadrotor2, x, u) = model.J
inertia_inv(model::Quadrotor2, x, u) = model.Jinv
mass_matrix(model::Quadrotor2, x, u) = Diagonal(@SVector fill(model.mass,3))

# TrajectoryOptimization.state_diff_size(::Quadrotor{<:UnitQuaternion,false}) = 13
#
# function state_diff(::Quadrotor2{R,false}, x::SVector, x0::SVector) where R
#       return x-x0
# end
#
# function state_diff_jacobian(::Quadrotor2{R,false}, x0::SVector) where R
#       return I
# end
#
# function TrajectoryOptimization.∇²differential(model::Quadrotor2{R,false}, x::SVector,
#             dx::SVector) where R
#       return I*0
# end

############################################################################################
#                                   QUATERNION SLACK
############################################################################################
Base.size(::Quadrotor2{:slack,<:UnitQuaternion}) = 13,5

function dynamics(model::Quadrotor2{:slack,<:UnitQuaternion}, x, u) where D

    r,q,v,ω = parse_state(model, x, false)
    q = q*u[5]

    F = forces(model, x, u)
    τ = moments(model, x, u)
    M = mass_matrix(model, x, u)
    J = inertia(model, x, u)
    Jinv = inertia_inv(model, x, u)

    xdot = v
    qdot = kinematics(q,ω)
    vdot = M\F
    ωdot = Jinv*(τ - ω × (J*ω))

    build_state(model, xdot, qdot, vdot, ωdot)
end

function forces(model::Quadrotor2{:slack,<:UnitQuaternion}, x, u)
    q = orientation(model, x, false)
    q = q*u[5]
    kf = model.kf
    g = model.gravity
    m = model.mass

    w1 = u[1]
    w2 = u[2]
    w3 = u[3]
    w4 = u[4]

    F1 = max(0,kf*w1);
    F2 = max(0,kf*w2);
    F3 = max(0,kf*w3);
    F4 = max(0,kf*w4);
    F = @SVector [0., 0., F1+F2+F3+F4] #total rotor force in body frame

    m*g + q*F # forces in world frame
end

function state_diff(::Quadrotor2{:slack,<:UnitQuaternion}, x::SVector, x0::SVector)
    return x-x0
end

function state_diff_jacobian(model::Quadrotor2{:slack,<:UnitQuaternion}, s::SVector)
    return I
end

function TrajectoryOptimization.∇²differential(model::Quadrotor2{:slack,<:UnitQuaternion},
        x::SVector, dx::SVector) where R
    return I*0
end

function TrajectoryOptimization.inverse_map_jacobian(
        model::Quadrotor2{:slack,<:UnitQuaternion}, x::SVector)
    return I
end

function TrajectoryOptimization.inverse_map_∇jacobian(
        model::Quadrotor2{:slack,<:UnitQuaternion}, x::SVector, b::SVector)
    return I*0
end

############################################################################################
#                                       Lee Quad
############################################################################################

@with_kw struct LeeQuad{B,R,T} <: TestRigidBody{B,R}
    mass::T = 4.34
    J::Diagonal{T,SVector{3,T}} = Diagonal(@SVector [0.0820, 0.0845, 0.1377])
    d::T = 0.315  # motor distance
    cf::T = 8.004e-4 # torque constant
    g::T = -9.81
end

Base.size(::LeeQuad{B,<:UnitQuaternion}) where B = 13,4
Base.size(::LeeQuad) = 12,4

(::Type{<:LeeQuad})(Rot; use_rot=true, kwargs...) = LeeQuad{use_rot,Rot,Float64}(;kwargs...)

inertia(model::LeeQuad, x, u) = model.J
inertia_inv(model::LeeQuad, x, u) = inv(model.J)
mass_matrix(model::LeeQuad, x, u) = Diagonal(@SVector fill(model.mass,3))

function Dynamics.Quadrotor2(model::LeeQuad)
    Dynamics.Quadrotor2{UnitQuaternion{Float64,CayleyMap}}(
        mass=model.mass, J=model.J, motor_dist=model.d, km=model.cf)
end

function forceMatrix(model::LeeQuad)
    d,c = model.d, model.cf
    @SMatrix [
        1  1  1  1;
        0 -d  0  d;
        d  0 -d  0;
       -c  c -c  c;
    ]
end

function wrenches(model::LeeQuad, x::SVector, u::SVector)
    mass,g = model.mass, model.g
    q = orientation(model, x)
    R = rotmat(q)

    # Calculate force and moments
    u = min.(u, 0)  # keep forces negative
    fM = forceMatrix(model)*u
    f = fM[1]
    M = @SVector [fM[2], fM[3], fM[4]]
    e3 = @SVector [0,0,1]
    F = mass*g*e3 - f*R*e3
    return F,M
end

trim_controls(model::LeeQuad) = @SVector fill(model.g*model.mass/4.0, size(model)[2])

############################################################################################
#                                  BODY FRAME LEE QUAD
############################################################################################

@with_kw struct BodyFrameLeeQuad{B,R,T} <: TestRigidBody{B,R}
    mass::T = 4.34
    J::Diagonal{T,SVector{3,T}} = Diagonal(@SVector [0.0820, 0.0845, 0.1377])
    d::T = 0.315  # motor distance
    cf::T = 8.004e-4 # torque constant
    g::T = -9.81
end

Base.size(::BodyFrameLeeQuad{B,<:UnitQuaternion}) where B = 13,4
Base.size(::BodyFrameLeeQuad) = 12,4

(::Type{<:BodyFrameLeeQuad})(Rot; use_rot=true, kwargs...) =
    BodyFrameLeeQuad{use_rot,Rot,Float64}(;kwargs...)

inertia(model::BodyFrameLeeQuad, x, u) = model.J
inertia_inv(model::BodyFrameLeeQuad, x, u) = inv(model.J)
mass_matrix(model::BodyFrameLeeQuad, x, u) = Diagonal(@SVector fill(model.mass,3))

function forceMatrix(model::BodyFrameLeeQuad)
    d,c = model.d, model.cf
    @SMatrix [
        1  1  1  1;
        0 -d  0  d;
        d  0 -d  0;
       -c  c -c  c;
    ]
end

function dynamics(model::BodyFrameLeeQuad, x::SVector, u::SVector)
    # r position in inertial frame
    # q rotation from body to inertial frame
    # v velocity in body frame
    # ω angular velocity in body frame
    r,q,v,ω = parse_state(model, x)

    F,τ = wrenches(model, x, u)
    M = mass_matrix(model, x, u)
    J = inertia(model, x, u)
    Jinv = inertia_inv(model, x, u)

    xdot = q*v
    qdot = kinematics(q,ω)
    vdot = (q\F)/model.mass - ω × v
    ωdot = Jinv*(τ - ω × (J*ω))

    build_state(model, xdot, qdot, vdot, ωdot)
end

function wrenches(model::BodyFrameLeeQuad, x::SVector, u::SVector)
    mass,g = model.mass, model.g
    q = orientation(model, x)
    R = rotmat(q)

    # Calculate force and moments
    u = min.(u, 0)  # keep forces negative
    fM = forceMatrix(model)*u
    f = fM[1]
    M = @SVector [fM[2], fM[3], fM[4]]
    e3 = @SVector [0,0,1]
    F = mass*g*e3 - f*R*e3
    return F,M
end

trim_controls(model::BodyFrameLeeQuad) = @SVector fill(model.g*model.mass/4.0, size(model)[2])


text:
export pfork, pfork_shm

function pfork(n::Integer, f::Function, args...)
    if nprocs() > 1
        error("pfork is incompatible with worker processes. Execute a `rmprocs(workers())` and try again.")
    end

    pipes = cell(n)
    cpids = cell(n)
    responses=cell(n)

    for i in 1:n
        pipe_fds = Array(Cint,2)
        assert( 0 == ccall(:pipe, Cint, (Ptr{Cint},), pipe_fds))
        pipes[i] = pipe_fds

        cpids[i] = ccall(:fork, Cint, ())

        if cpids[i] == 0
            try
                # child
                # close the read end in the child
                ccall(:close, Cint, (Cint, ),  pipes[i][1])

                # seed the random number generator differently for each child...
                srand([uint32(time()), uint32(i)])

                rv = f(i, args...)

                iob=IOBuffer()
                if rv == nothing
                    serialize(iob, ())
                else
                    serialize(iob, rv)
                end

                arr=takebuf_array(iob)
                nb = length(arr)
                nb_pending = nb
                while nb_pending > 0
                    #ssize_t write(int fd, const void *buf, size_t count);
                    bw = ccall(:write, Csize_t, (Cint, Ptr{Uint8}, Csize_t), pipes[i][2], pointer(arr, (nb - nb_pending) + 1), nb_pending)
                    if bw <= 0
                        exit(1)
                    end
                    nb_pending -= bw
                end

                #println("child process exiting")
                ccall(:close, Cint, (Cint, ),  pipes[i][2])
            catch e
                println("forked child exiting with exception $e")
            finally
                exit()
            end
        else
            # close the write end in the parent
            ccall(:close, Cint, (Cint, ),  pipes[i][2])
        end
    end

    # Read all responses and wait for children to exit
    @sync begin
        for i in 1:n
            #read all data in the pipe...
            @async begin
                buf = zeros(Uint8,65536)
                iob = IOBuffer()
                while true
                    br = ccall(:read, Csize_t, (Cint, Ptr{Uint8}, Csize_t), pipes[i][1], buf, 65536)
                    if br > 0
                        write(iob, buf[1:br])
                    else
                        break;
                    end
                end

                #close the read end
                ccall(:close, Cint, (Cint, ),  pipes[i][1])

                #retcode = Cint[0]
                #cpid_exited = ccall(:waitpid, Cint, (Cint, Ptr{Cint}, Cint), cpids[i], retcode, 0)
                # Sometimes cpid_exited is -1 and errno is ECHILD, removing check on return value
                # of waitpid for now.
#                assert(cpid_exited == cpids[i])
                #println("child process $(cpids[i]) exited")

                responses[i] = iob
            end
        end
    end

    for i in 1:n
        seekstart(responses[i])
        responses[i] = deserialize(responses[i])
    end

    responses
end

pfork_shm(n::Int, return_type::Type, return_size::Tuple, f::Function, args...) =
    pfork_shm(n::Int, {(return_type, return_size)}, f::Function, args...)[1]
# forks the current process n times, invoking f,
# and providing a shared memory segment for the return value
#
# n                  ... number of forks. passed to pfork
# sharedmems         ... type of the return value. passed to ShmCfg
# retreturn_sizesize ... size of the return value. passed to ShmCfg
# f                  ... function with the signature is f(ind, result, args...)
#                        where result is a shared memory segment of type return_type
#                        and size return_size

function pfork_shm(n::Int, shmems::Array, f::Function, args...)
  # f                  ... function with the signature is f(ind, result1, .., resultN, args...)

  map(x->assert(isa(x,Tuple)),shmems)

  # set up the shared memory segment
  return_shms = map(x->ShmCfg(gensym("pfork_shm"),x[1],x[2]), shmems)
  unlink_shm(return_shms)
  setup_shm(return_shms)

  # wrap f, passing the return memory segment
  function g(ind, a...)
    rets = map(x->task_local_storage(x.sname), return_shms)
    #@show rets tuple(ind, rets..., a...)
    f(ind, rets..., a...)
    nothing
  end

  # invoke the function
  #@show n g args
  pfork(n, g, args...)

  # release the shared memory segment and return
  unlink_shm(return_shms)
  r = map(x->task_local_storage(x.sname), return_shms)
  map(x->task_local_storage(x.sname, nothing), return_shms)
  r
end

ext: jl


text:
#  Marek Lipert 2015

mean(sequence) = reduce(+,sequence)/length(sequence)
sigma(sequence,mean) = sqrt(reduce(+,map((x) -> (x-mean)^2,sequence))/(length(sequence)-1))

# filters

function kalman(sequence)
  x0 = mean(sequence)
  s = 3.0                # measurement error
  e = 1000
  r = 0.0
  for i in sequence
    p = e/(e+s)
    r = r + p * (i - r)
    e = (1 - p)*e
  end
  r
end

function motion_kalman(sequence)
  dt = 1.0
  sigma_process = 0.3 # process error
  r = 3.0    # measurement error
  L = 1000
  F = [1 dt; 0 1]
  Q = [dt^3 / 3.0 dt^2 / 2; dt^2 / 2 dt ] .* sigma_process
  H = [1 0]
  R = [r]

  x0 = [0.0 ; 0.0]
  P0 = [L 0; 0 L]
  x1 = x0
  P1 = P0
  U = [1 0 ; 0 1]
  for z in sequence[2:]
    # predict
    x1 = F * x1
    P1 = F * P1 * F' + Q
    # update
    y = [z] - H * x1
    S = H * P1 * H' + R
    K = P1 * H'  ./ S[1,1]
    x1 = x1 + K .* y[1]
    P1 = (U - K*H) * P1
  end
  (H * x1)[1]
end

s0 = [25 7 14 8 13 20 ]
s1 = [3.2 3.1 3.2 3.3 3.5]
s2 = [3.2 4.2 2.8 6.1 2.7]

measures = (s0,s1,s2)

for l in measures
  print(l)
  @printf("mean %.2f sd %.2f kalman %.2f motion_kalman %.2f\n",mean(l),sigma(l,mean(l)),kalman(l),motion_kalman(l))
  println("------------------")
end

max_stars_repo_name: mll/iBeacon-supporting-calculations

max_issues_repo_name: mll/iBeacon-supporting-calculations

max_forks_repo_name: mll/iBeacon-supporting-calculations


text:
export runmodel

"Returns a tuple whose entries are the mean fitness of each species."
function mean_fitness_per_species(model::ABM)
  mean_fitness = Array{Float64}(undef, model.nspecies)
  for species in 1:model.nspecies
    fitness = mean([i.W for i in values(model.agents) if i.species == species])
    mean_fitness[species] = fitness
  end

  return mean_fitness
end

function species_N(model::ABM)
  allagents = model.agents
  if length(allagents) == 0
    return fill(0, model.nspecies)
  else
    counts = countmap([a.species for a in values(model.agents)])
    output = fill(0, model.nspecies)
    for (k, v) in counts
      output[k] = v
    end
    return output
  end
end

"""
    runmodel(param_file::AbstractString; kwargs)

Creates and runs a model given `parameters`. Returns a `DataFrame` of collected data, which are specified by `kwargs`.

# Keywords

* adata=[] agent data to be collected. Either agent fields or functions that accept an agent as input can be put in the array. To aggregate collected data, provide tuples inside the array. For example, to collect mean and median fitness of individuals which is in field `W`, your array will be [(:W,mean), (:W,median)].
* mdata=[mean_fitness_per_species] model data to be collected. By default, collects mean population fitness per species. Each row of the output DataFrame corresponds to all agents and each column is the value function applied to a field. The functions in the array are applied to the model object.
* when=nothing The generations from which data are collected. By default collect at all steps.
* replicates::Int = 0 Number of replicates per simulation.
* parallel::Bool = false Whether to run replicates in parallel. If `true`, you should add processors to your julia session (e.g. by `addprocs(n)`) and define your parameters and `EvoDynamics` on all workers. To do that, add `@everywhere` before them. For example, `@everywhere EvoDynamics`.
* seeds = optionally, provide an array of integers as seeds for each replicate.
"""
function runmodel(param_file::AbstractString;
  adata=nothing, mdata=[mean_fitness_per_species, species_N],
  when=nothing,
  replicates::Int = 0,
  parallel::Bool = false,
  seeds = nothing
  )

  if !isnothing(seeds) && replicates > 0
    @assert length(seeds) == replicates "If seeds is not nothing, it has to be an array of ints as long as replicates."
  end

  dd = load_parameters(param_file)
  if isnothing(when)
    whenn = 0:dd[:model]["generations"]
  else
    whenn = when
  end
  # create model
  if replicates == 0
    model = model_initiation(dd)

    # run model and collect data
    agdata, modata = run!(model, agent_step!, model_step!, model.generations, adata=adata, mdata=mdata, when=whenn, agents_first=false)
    return agdata, modata, [model]
  else
    models = [model_generator(i, seeds, param_file) for i in 1:replicates]

    agdata, modata, models = ensemblerun!(models, agent_step!, model_step!, dd[:model]["generations"], adata=adata, mdata=mdata, when=whenn, parallel=parallel, agents_first=false)
    return agdata, modata, models
  end
end

function model_generator(counter, seeds, param_file)
  dd = load_parameters(param_file)
  if isnothing(seeds)
    dd[:model]["seed"] = rand(1:1000_000)
  else
    dd[:model]["seed"] = seeds[counter]
  end
  model = model_initiation(dd)
  return model
end

"""
    expand_iterable_columns(agdata, modata)

Expand columns that have arrays (multiple values). Inputs are agent data and model data outputs.
"""
function expand_iterable_columns!(agdata, modata)
  allnames = Agents.names(agdata)
  tuple_colsa = String[]
  for nam in allnames
    if eltype(agdata[!, nam]) <: AbstractArray
      push!(tuple_colsa, nam)
    end
  end
  allnames = Agents.names(modata)
  tuple_colsm = String[]
  for nam in allnames
    if eltype(modata[!, nam]) <: AbstractArray
      push!(tuple_colsm, nam)
    end
  end

  for nam in tuple_colsa
    nitems = length(agdata[1, nam])
    for item in 1:nitems
      agdata[!, Symbol(nam * "_$item")] = getindex.(agdata[!, Symbol(nam)], item)
    end
    Agents.select!(agdata, Agents.Not(nam))  # remove the column with tuples
  end

  for nam in tuple_colsm
    nitems = length(modata[1, nam])
    for item in 1:nitems
      modata[!, Symbol(nam * "_$item")] = getindex.(modata[!, Symbol(nam)], item)
    end
    Agents.select!(modata, Agents.Not(Symbol(nam)))  # remove the column with tuples
  end
  return agdata, modata
end


text:
function x0_PT_flash(model::EoSModel,T,P,x0)
    #check each T with T_scale, if treshold is over, replace Pi with inf

 _0 = zero(P+T+first(x0))
 pure = split_model(model)
 crit = crit_pure.(pure)

 T_c = [tup[1] for tup in crit]
 V_c = [tup[3] for tup in crit]
 nan = _0/_0
 sat_nan = (nan,nan,nan)
 replaceP = ifelse.(T_c .< T,true,false)

 sat = [if !replaceP[i] sat_pure(pure[i],T) else sat_nan end for i in 1:length(pure)]

 P_sat = [tup[1] for tup in sat]
 Ki = zeros(eltype(x0),length(x0))
 for i in 1:length(x0)
     if !replaceP[i]
         Ki[i] = P_sat[i][1]/P
     else
         Ki[i] = pressure(pure[i],V_c[i],T)/P
     end
 end

return Ki
end

function PT_flash(model::EoSModel,P,T,z,K0=nothing)
    _0 = zero(P+T+first(z))
    _1 = one(_0)
    if K0 === nothing
        K0 = x0_PT_flash(model,T,P,z)
    end
    g_0 = rr_flash_eval(K0,z,_0)
    g_1 = rr_flash_eval(K0,z,_1)

    if g_0 <= 0  #bubble point assumption
        println("bubble point")
        β = _0
        β0 = _0
        xil = copy(z)
        xiv = rr_flash_vapor(K0,z,_0)
    elseif g_1 >= 0 #dew point assumption
        println("dew point")
        β = _1
        β0 = _1
        xil = rr_flash_liquid(K0,z,_1)
        xiv = copy(z)
    else #two phase assumption
        println("flash point")
        β = rr_vle_vapor_fraction(K0,z)

        xil = rr_flash_liquid(K0,z,β)
        xiv = rr_flash_vapor(K0,z,β)
    end
    logϕl = similar(xil)
    logϕv = similar(xiv)

    vv = volume(model,P,T,xiv,phase=:v)
    vl = volume(model,P,T,xil,phase=:l)
    @show vl
    @show vv
    @show β

    return xil, xiv
end


text:
include("divide_singularity.jl")
include("JacobiWeight.jl")
include("JacobiWeightOperators.jl")
include("JacobiWeightChebyshev.jl")
include("LogWeight.jl")
include("ExpWeight.jl")


text:
using Test
using OpenQuantumSystems
using Random, SparseArrays, LinearAlgebra, StableRNGs

import OrdinaryDiffEq

@testset "liouville" begin

    Random.seed!(StableRNG(0), 1)

    D(op1::Array, op2::Array) = abs(norm(op1 - op2))
    D(x1::StateVector, x2::StateVector) = norm(x2 - x1)
    D(op1::AbstractOperator, op2::AbstractOperator) =
        abs(tracedistance_nh(dense(op1), dense(op2)))
    D(op1::AbstractSuperOperator, op2::AbstractSuperOperator) =
        abs(tracedistance_nh(dense(op1), dense(op2)))

    mode1 = Mode(0.2, 1.0)
    Energy = [0.0, 200.0]
    mol1 = Molecule([mode1], 2, Energy)
    mol2 = Molecule([mode1], 2, Energy)
    agg = Aggregate([mol1, mol2])
    aggInds = getIndices(agg)
    aggIndsLen = length(aggInds)
    basis = GenericBasis([aggIndsLen])
    FCFact = getFranckCondonFactors(agg, aggInds)
    Ham = getAggHamiltonian(agg, aggInds, FCFact)

    ket0 = randstate(basis)
    rho0 = dm(ket0)
    # tests have to be quick enough
    tspan = [0.0:0.1:1.0;]

    T, rho_t = liouvilleVonNeumann(
        rho0,
        tspan,
        Ham;
        reltol = 1e-10,
        abstol = 1e-10,
        alg = OrdinaryDiffEq.Tsit5(),
    )
    for t_i = 1:length(tspan)
        U_op = evolutionOperator(Ham, tspan[t_i])
        rho = U_op * rho0 * U_op'
        @test 1e-8 > D(rho, rho_t[t_i])
        # println(t_i, " ", D(rho.data, rho_t[t_i].data))
    end

    # diabled for now for, it takes about 20 s to run this simulation
    #=
    T, rho_t = liouvilleVonNeumann(rho0, Ham, tspan; reltol=1e-12, abstol=1e-12, alg=OrdinaryDiffEq.Vern7())
    for t_i in 1:length(tspan)
        U_op = evolutionOperator(Ham, tspan[t_i])
        rho = U_op * rho0 * U_op'
        @test 1e-14 > D(rho, rho_t[t_i])
        # println(t_i, " ", D(rho.data, rho_t[t_i].data))
    end
    =#

end # testset


text:
# ---
# title: 406. Queue Reconstruction by Height
# id: problem406
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Greedy
# link: <https://leetcode.com/problems/queue-reconstruction-by-height/description/>
# hidden: true
# ---
#
# You are given an array of people, `people`, which are the attributes of some
# people in a queue (not necessarily in order). Each `people[i] = [hi, ki]`
# represents the `ith` person of height `hi` with **exactly** `ki` other people
# in front who have a height greater than or equal to `hi`.
#
# Reconstruct and return _the queue that is represented by the input
# array_`people`. The returned queue should be formatted as an array `queue`,
# where `queue[j] = [hj, kj]` is the attributes of the `jth` person in the queue
# (`queue[0]` is the person at the front of the queue).
#
#
#
# **Example 1:**
#
#
#
#     Input: people = [[7,0],[4,4],[7,1],[5,0],[6,1],[5,2]]
#     Output: [[5,0],[7,0],[5,2],[6,1],[4,4],[7,1]]
#     Explanation:
#     Person 0 has height 5 with no other people taller or the same height in front.
#     Person 1 has height 7 with no other people taller or the same height in front.
#     Person 2 has height 5 with two persons taller or the same height in front, which is person 0 and 1.
#     Person 3 has height 6 with one person taller or the same height in front, which is person 1.
#     Person 4 has height 4 with four people taller or the same height in front, which are people 0, 1, 2, and 3.
#     Person 5 has height 7 with one person taller or the same height in front, which is person 1.
#     Hence [[5,0],[7,0],[5,2],[6,1],[4,4],[7,1]] is the reconstructed queue.
#
#
# **Example 2:**
#
#
#
#     Input: people = [[6,0],[5,0],[4,0],[3,2],[2,2],[1,4]]
#     Output: [[4,0],[5,0],[2,2],[3,2],[1,4],[6,0]]
#
#
#
#
# **Constraints:**
#
#   * `1 <= people.length <= 2000`
#   * `0 <= hi <= 106`
#   * `0 <= ki < people.length`
#   * It is guaranteed that the queue can be reconstructed.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using ProgressMeter
using Printf
using MATLAB
using BSON: @save, @load
using Dates
include("Dimer_Lagrangian_Descriptor.jl")

@everywhere begin
    include("Dimer_Lagrangian_Descriptor.jl")
    push!(LOAD_PATH, "/home/brandon_behring/Desktop/Leap_Frog_Code")
    using ProgressMeter
    using Printf
    using MATLAB
    using .Dimer_Lagrangian_Descriptor
    using BSON: @save, @load
    using Dates
    # H=range(.11,stop=0.145,length=107)
    # Energy=0.112012
    # Energy=.125789
    # Energy=.12501
    # Energy=.12501
    # Energy=.1249
    # Energy=.125
    Yfind(h)=sqrt(h/(2h+1));

    H = range(.128, stop = .132, length = 5)
    # H = range(.11, stop = .21, length = 15)

    count = 1
    right_now = replace(replace(replace(string(Dates.now()),"."=>"_"),":"=>"_"),"-"=>"_")

end
while count <= length(H)
    # @everywhere t_end=15
      @everywhere t_end=150

    @everywhere P=Yfind(H[count])

    @everywhere     N=500;#50
    # Q_start=.20
    # Q_end=2.5
    # @everywhere  Q_end=5e-2
    # @everywhere  Q_end=2.25
    @everywhere  Q_end=4.5e-2
    @everywhere  Q_start=-Q_end
    # @everywhere  Q_end=.00002
    @everywhere  Q_start=0.0

    # @everywhere  n_iter_P=2001
    # @everywhere  P_start=-.5
    @everywhere  n_iter_P=N
    @everywhere  n_iter_Q=N
    @everywhere  P_start=P-5e-3
    @everywhere  P_end=P+4e-3

    # @everywhere  P_start=0.0
    # @everywhere  P_end=.45
    # # P_end=1
    @everywhere  t_end_mesh = t_end * ones(n_iter_Q,n_iter_P)
    @everywhere  ArrP=range(P_start,stop=P_end,length=n_iter_P)
    @everywhere  ArrQ=range(Q_start,stop=Q_end,length=n_iter_Q)
    @everywhere location_fig="/media/brandon_behring/Extra_Space/MATLAB_FIGURES/"
    @everywhere location_bson="/media/brandon_behring/Extra_Space/BSON_FILES/"
    @everywhere  mesh = [(Q, P) for P in ArrP, Q in ArrQ]
    @everywhere  mesh_list = reshape(mesh, 1, :)
@everywhere Energy_mesh = H[count] * ones(n_iter_Q,n_iter_P)
@everywhere h=replace(@sprintf("%.15f",H[count]),"."=>"_")
@everywhere nQ = @sprintf("_%5d",n_iter_Q)
@everywhere Q_win = @sprintf("_%.3f",Q_end)
@everywhere P_win = @sprintf("_%.3f",P_end)
@everywhere tend= @sprintf("_%d",t_end)
@everywhere tend_T= @sprintf("%4f",t_end)
@everywhere h_BSON= location_bson * "BandF_Dimer_lagrangian_descriptors_data"*h*nQ*tend*Q_win*P_win*right_now*".bson"
@everywhere file_name1=location_fig* "Lagrangian_Descriptor_Dimer_1_"*h*nQ*tend*Q_win*P_win*right_now*".fig"
@everywhere file_name2=location_fig* "Lagrangian_Descriptor_Dimer_2_"*h*nQ*tend*Q_win*P_win*right_now*".fig"
@everywhere file_name3=location_fig* "Lagrangian_Descriptor_Dimer_3_"*h*nQ*tend*Q_win*P_win*right_now*".fig"

@everywhere h_title = @sprintf("h= %.6f",H[count])*" with tend="*tend_T*" and "*nQ
println(h_title)

LD= @showprogress pmap(Dimer_Lagrangian_Descriptor.Dimer_Lagrangian_Descriptor_Function, mesh, Energy_mesh, t_end_mesh)

dx=(Q_end-Q_start)/N
dy=(P_end-P_start)/N

gM=Dimer_Lagrangian_Descriptor.gradient_matrix_4(LD,dx,dy)
gradM_X,gradM_Y=Dimer_Lagrangian_Descriptor.gradient_matrix(LD,dx,dy)

Log_gradM_X=log10.(gradM_X)
Log_gradM_Y=log10.(gradM_Y)

biggest_X=maximum(filter(!isnan,Log_gradM_X))
biggest_Y=maximum(filter(!isnan,Log_gradM_Y))
smallest_X=minimum(filter(!isnan,Log_gradM_X))
smallest_Y=minimum(filter(!isnan,Log_gradM_Y))
Nl,Nl=size(gradM_X)
gradM=zeros(Nl,Nl)
for i=1:Nl
    for j=1:Nl
        cutoff_X=2
        cutoff_Y=2
        X=maximum([Log_gradM_X[i,j], cutoff_X])
        Y=maximum([Log_gradM_Y[i,j], cutoff_Y])
        gradM[i,j]=log10(10^(2*X)+10^(2*Y))
    end
end

# gradM=10 .^gradM
SAVE_DATA=Dict("LD"=>LD,"gradM"=>gradM,"Q_end"=>Q_end,"Q_start"=>Q_start, "P_start"=>P_start, "P_end"=>P_end,"N"=>N,"Energy"=>H[count],"ArrP"=>ArrP, "ArrQ"=>ArrQ,"t_end"=>t_end)
@save h_BSON SAVE_DATA

mat"figure();set(gcf, 'Position',  [0, 0, 1500, 1000]); hold on;"
mat"title($h_title)"
mat"axis([ -$Q_end,$Q_end,$P_start,$P_end ])"

# mat"imagesc([$Q_start , $Q_end ],[$P_start, $P_end],$gradM, $clims)"
# mat"imagesc([-$Q_start, -$Q_end ],[-$P_start, -$P_end], $gradM, $clims)"
mat"imagesc([$Q_start , $Q_end ],[$P_start, $P_end],$LD)"
mat"imagesc([-$Q_start, -$Q_end ],[$P_start, $P_end], $LD)"
# mat"imagesc([$Q_start, $Q_end ],[-$P_start, -$P_end], $LD)"
# mat"imagesc([-$Q_start, -$Q_end ],[-$P_start, -$P_end], $LD)"
# mat"axis([ -$Q_end,$Q_end,-$P_end,$P_end ])"

# mat"imagesc([-$Q_start, -$Q_end ],[-$P_start, -$P_end], $gradM)"
# mat"imagesc([0,-$Q_end ],[$P_start, $P_end],$gradM, $clims)"
# mat"colormap winter"
# mat" cmocean('balance')"
mat"colormap(twilight)"

mat"colorbar"
Q_fix=sqrt(6)/3
mat"plot($Q_fix,0,'b.','MarkerSize',30)"
mat"plot(-$Q_fix,0,'b.','MarkerSize',30)"

mat"plot(0,$P,'r.','MarkerSize',10)"
mat"plot(0,-$P,'r.','MarkerSize',30)"
# mat"axis([ -$Q_end,$Q_end,$P_start,$P_end ])"
# mat"axis([ $Q_start,$Q_end,-$P_end,$P_end ])"
# mat"axis([ $Q_start, $Q_end,$P_start,$P_end ])"

mat"savefig($file_name1)"

SAVE_DATA=Dict("LD"=>LD,"gradM"=>gradM,"Q_end"=>Q_end,"Q_start"=>Q_start, "P_start"=>P_start, "P_end"=>P_end,"N"=>N,"Energy"=>H[count],"ArrP"=>ArrP, "ArrQ"=>ArrQ,"t_end"=>t_end)
@save h_BSON SAVE_DATA

mat"figure();set(gcf, 'Position',  [0, 0, 1500, 1000]); hold on;"
mat"title($h_title)"
mat"axis([ -$Q_end,$Q_end,$P_start,$P_end ])"

# mat"imagesc([$Q_start , $Q_end ],[$P_start, $P_end],$gradM, $clims)"
# mat"imagesc([-$Q_start, -$Q_end ],[-$P_start, -$P_end], $gradM, $clims)"
mat"imagesc([$Q_start , $Q_end ],[$P_start, $P_end],$gradM)"
mat"imagesc([-$Q_start, -$Q_end ],[$P_start, $P_end], $gradM)"
# mat"imagesc([$Q_start, $Q_end ],[-$P_start, -$P_end], $gradM)"
# mat"imagesc([-$Q_start, -$Q_end ],[-$P_start, -$P_end], $gradM)"
# mat"axis([ -$Q_end,$Q_end,-$P_end,$P_end ])"

# mat"imagesc([-$Q_start, -$Q_end ],[-$P_start, -$P_end], $gradM)"
# mat"imagesc([0,-$Q_end ],[$P_start, $P_end],$gradM, $clims)"
# mat"colormap winter"
mat" cmocean('balance')"

mat"colorbar"
Q_fix=sqrt(6)/3
mat"plot($Q_fix,0,'b.','MarkerSize',30)"
mat"plot(-$Q_fix,0,'b.','MarkerSize',30)"

mat"plot(0,$P,'r.','MarkerSize',10)"
mat"plot(0,-$P,'r.','MarkerSize',30)"
# mat"axis([ -$Q_end,$Q_end,$P_start,$P_end ])"
# mat"axis([ $Q_start,$Q_end,-$P_end,$P_end ])"
# mat"axis([ $Q_start, $Q_end,$P_start,$P_end ])"

mat"savefig($file_name2)"

@everywhere global count += 1

end


text:
#=

ESSE equations

Ignacio Quintero Mächler

t(-_-t)

September 16 2017

=#

"""
    make_esse_s(k::Int64, x::Array{Float64,1}, y::Array{Float64,1})

Speciation ESSE equation function for integration, one z(t).
"""
function make_esse_s(k ::Int64,
                     af::Function)

  # make Qf function
  Qf! = makeQ(k)

  # preallocate matrices & vectors
  bk::Int64 = k + (k*k)
  rt::Int64 = bk + k + 1
  Q ::Array{Float64,2} = zeros(k, k)
  Qs::Array{Float64,1} = zeros(k)
  Qe::Array{Float64,1} = zeros(k)

  # estimate indices previously
  qsi::UnitRange{Int64} = 1:k
  qei::UnitRange{Int64} = (k+1):2k

  function f(du::Array{Float64,1},
             u::Array{Float64,1},
             p::Array{Float64,1},
             t::Float64)

    Qf!(Q, p)
    @views BLAS.gemv!('N', 1.0, Q, u[qsi], 0.0, Qs)
    @views BLAS.gemv!('N', 1.0, Q, u[qei], 0.0, Qe)

    aft = af(t)::Float64

    for i in Base.OneTo(k)
      st      = p[i]*exp(p[bk+i]*aft)
      du[i]   =        - (st + p[k+i])*u[i]   + 2.0*st*u[i+k]*u[i]   + Qs[i]
      du[i+k] = p[k+i] - (st + p[k+i])*u[i+k] +     st*u[i+k]*u[i+k] + Qe[i]
    end

    return nothing
  end

  return f
end

"""
    make_esse_s(t,u,p,du)

Speciation ESSE equation function for integration, for `k` `z(t)` functions.
"""
function make_esse_s(k  ::Int64,
                     af!::Function,
                     md ::Bool)

  # make Qf function
  Qf! = makeQ(k)

  # preallocate matrices & vectors
  bk::Int64 = k+k*k
  Q ::Array{Float64,2} = zeros(k, k)
  Qs::Array{Float64,1} = zeros(k)
  Qe::Array{Float64,1} = zeros(k)
  r ::Array{Float64,1} = Array{Float64}(undef,k)

  # estimate indices previously
  qsi::UnitRange{Int64} = 1:k
  qei::UnitRange{Int64} = (k+1):2k

  function f(du, u, p, t)

    Qf!(Q, p)
    @views BLAS.gemv!('N', 1.0, Q, u[qsi], 0.0, Qs)
    @views BLAS.gemv!('N', 1.0, Q, u[qei], 0.0, Qe)

    af!(t, r)

    for i in Base.OneTo(k)
      st      = p[i]*exp(p[bk+i]*r[i])::Float64
      du[i]   =        - (st + p[k+i])*u[i]   + 2.0*st*u[i+k]*u[i]   + Qs[i]
      du[i+k] = p[k+i] - (st + p[k+i])*u[i+k] +     st*u[i+k]*u[i+k] + Qe[i]
    end

    return nothing
  end

  return f
end

"""
    make_esse_e(t,u,p,du)

Extinction ESSE equation function for integration, one `z(t)`.
"""
function make_esse_e(k ::Int64,
                     af::Function)

  # make Qf function
  Qf! = makeQ(k)

  # preallocate matrices & vectors
  bk::Int64 = k+k*k
  Q ::Array{Float64,2} = zeros(k, k)
  Qs::Array{Float64,1} = zeros(k)
  Qe::Array{Float64,1} = zeros(k)

  # estimate indices previously
  qsi::UnitRange{Int64} = 1:k
  qei::UnitRange{Int64} = (k+1):2k

  function f(du, u, p, t)

    Qf!(Q, p)
    @views BLAS.gemv!('N', 1.0, Q, u[qsi], 0.0, Qs)
    @views BLAS.gemv!('N', 1.0, Q, u[qei], 0.0, Qe)

    aft = af(t)::Float64

    for i in Base.OneTo(k)
      et      = p[k+i]*exp(p[bk+i]*aft)::Float64
      du[i]   =    - (p[i] + et)*u[i]   + 2.0*p[i]*u[i+k]*u[i]   + Qs[i]
      du[i+k] = et - (p[i] + et)*u[i+k] +     p[i]*u[i+k]*u[i+k] + Qe[i]
    end

    return nothing
  end

  return f
end

"""
    make_esse_e(t,u,p,du)

Extinction ESSE equation function for integration, for `k` `z(t)` functions.

"""
function make_esse_e(k  ::Int64,
                     af!::Function,
                     md ::Bool)

  # make Qf function
  Qf! = makeQ(k)

  # preallocate matrices & vectors
  bk::Int64 = k+k*k
  Q ::Array{Float64,2} = zeros(k, k)
  Qs::Array{Float64,1} = zeros(k)
  Qe::Array{Float64,1} = zeros(k)
  r ::Array{Float64,1} = Array{Float64}(undef, k)

  # estimate indices previously
  qsi::UnitRange{Int64} = 1:k
  qei::UnitRange{Int64} = (k+1):2k

  function f(du, u, p, t)

    Qf!(Q, p)
    @views BLAS.gemv!('N', 1.0, Q, u[qsi], 0.0, Qs)
    @views BLAS.gemv!('N', 1.0, Q, u[qei], 0.0, Qe)

    af!(t, r)

    for i in Base.OneTo(k)
      et      = p[k+i]*exp(p[bk+i]*r[i])::Float64
      du[i]   =    - (p[i] + et)*u[i]   + 2.0*p[i]*u[i+k]*u[i]   + Qs[i]
      du[i+k] = et - (p[i] + et)*u[i+k] +     p[i]*u[i+k]*u[i+k] + Qe[i]
    end

    return nothing
  end

  return f
end

"""
    make_esse_q(t,u,p,du)

Rates ESSE equation function for integration for one `z(t)` function.
"""
function make_esse_q(k ::Int64,
                     af::Function)

  # make Qt function
  Qft! = makeQt(k)

  Q ::Array{Float64,2} = zeros(k, k)
  Qs::Array{Float64,1} = zeros(k)
  Qe::Array{Float64,1} = zeros(k)

  # estimate indices previously
  qsi::UnitRange{Int64} = 1:k
  qei::UnitRange{Int64} = (k+1):2k

  function f(du, u, p, t)

    Qft!(Q, p, af(t))
    @views BLAS.gemv!('N', 1.0, Q, u[qsi], 0.0, Qs)
    @views BLAS.gemv!('N', 1.0, Q, u[qei], 0.0, Qe)

    for i in Base.OneTo(k)
      du[i]   =        - (p[i] + p[k+i])*u[i]   + 2.0*p[i]*u[i+k]*u[i]   + Qs[i]
      du[i+k] = p[k+i] - (p[i] + p[k+i])*u[i+k] +     p[i]*u[i+k]*u[i+k] + Qe[i]
    end

    return nothing
  end

  return f
end

"""
    make_esse_q(t,u,p,du)

Rates ESSE equation function for integration for multiple time functions.
"""
function make_esse_q(k  ::Int64,
                     af!::Function,
                     md ::Bool)

  # make Qt function
  Qft! = makeQt(k, af!)

  Q ::Array{Float64,2} = zeros(k, k)
  Qs::Array{Float64,1} = zeros(k)
  Qe::Array{Float64,1} = zeros(k)

  # estimate indices previously
  qsi::UnitRange{Int64} = 1:k
  qei::UnitRange{Int64} = (k+1):2k

  function f(du, u, p, t)

    Qft!(Q, p, t)
    @views BLAS.gemv!('N', 1.0, Q, u[qsi], 0.0, Qs)
    @views BLAS.gemv!('N', 1.0, Q, u[qei], 0.0, Qe)

    for i in Base.OneTo(k)
      du[i]   =        - (p[i] + p[k+i])*u[i]   + 2.0*p[i]*u[i+k]*u[i]   + Qs[i]
      du[i+k] = p[k+i] - (p[i] + p[k+i])*u[i+k] +     p[i]*u[i+k]*u[i+k] + Qe[i]
    end

    return nothing
  end

  return f
end

"""
    makeQt(k)

Transform instantaneous transition rates into
the instantaneous transition matrix `Q` for **one**
`z(t)` functions.
"""
function makeQt(k::Int64)

  pidx::UnitRange{Int64} = (2k+1):(k+k*k)
  Didx::StepRange{Int64} = 1:(k+1):(k*k)
  lidx::Int64            = lastindex(pidx)
  bk  ::Int64            = k+k*k
  Qidx::Array{Int64,1}   = setdiff(1:k*k, Didx)

  function f(Q  ::Array{Float64,2},
             p  ::Array{Float64,1},
             aft::Float64)

    @inbounds begin
      # assign parameters
      for i in Base.OneTo(lidx)
        Q[Qidx[i]] = p[pidx[i]]*exp(p[bk + i]*aft)::Float64
      end

      #assign 0s
      for i in Didx
        Q[i] = 0.0
      end

      #row sums
      for j = Base.OneTo(k), i = Base.OneTo(k)
        if i == j
          continue
        end
        Q[Didx[i]] -= Q[i,j]::Float64
      end

    end

    return nothing
  end
end

"""
    makeQt(k::Int64, af!::Function)

Transform instantaneous transition rates into
the instantaneous transition matrix `Q` for **multiple**
`z(t)` functions.
"""
function makeQt(k::Int64, af!::Function)

  pidx::UnitRange{Int64} = (2k+1):(k+k*k)
  Didx::StepRange{Int64} = 1:(k+1):(k*k)
  lidx::Int64            = lastindex(pidx)
  bk  ::Int64            = k+k*k
  Qidx::Array{Int64,1}   = setdiff(1:k*k, Didx)
  r   ::Array{Float64,1} = Array{Float64}(undef, k*k - k)

  function f(Q  ::Array{Float64,2},
             p  ::Array{Float64,1},
             t  ::Float64)

    @inbounds begin
      # estimate z(t)s
      af!(t, r)

      for i in Base.OneTo(lidx)
        Q[Qidx[i]] = p[pidx[i]]*exp(p[bk+i]*r[i])::Float64
      end

      #assign 0s
      for i in Didx
        Q[i] = 0.0
      end

      #row sums
      for j = Base.OneTo(k), i = Base.OneTo(k)
        if i == j
          continue
        end
        Q[Didx[i]] -= Q[i,j]::Float64
      end
    end

    return nothing
  end

end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
# The Fibonacci sequence is defined by the recurrence relation:
#
# F_n = F_n−1 + F_n−2, where F_1 = 1 and F_2 = 1.
# Hence the first 12 terms will be:
#
# F_1 = 1
# F_2 = 1
# F_3 = 2
# F_4 = 3
# F_5 = 5
# F_6 = 8
# F_7 = 13
# F_8 = 21
# F_9 = 34
# F_10 = 55
# F_11 = 89
# F_12 = 144
# The 12th term, F_12, is the first term to contain three digits.
#
# What is the index of the first term in the Fibonacci sequence to contain 1000
# digits?

using ProjectEulerSolutions

# There's a pattern in the number of digits in Fibonacci numbers.  It appears
# to be a new digit every 5 numbers or so.  We use that to estimate the upper
# bound for how far to generate numbers.
#
# map(x -> ndigits(x), fibonacci_count(48))
# [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 5,
#  6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10]
function p025solution(n_digits::Integer=1)::Integer

    fibs = fibonacci_bound(BigInt(10)^(n_digits-1))
    n = length(fibs) + 1

    if ndigits(fibs[end]) < n_digits
        return n + 1
    else
        return n
    end
end

p025 = Problems.Problem(p025solution)

Problems.benchmark(p025, 1000)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
module JuDoc

using JuDocTemplates

using Markdown
using Markdown: htmlesc
using Dates # see jd_vars
using DelimitedFiles: readdlm
using OrderedCollections
using Pkg
using DocStringExtensions: SIGNATURES, TYPEDEF

import Logging
import LiveServer
import Base.push!
import NodeJS
import Literate
import HTTP

export serve, publish, cleanpull, newsite, optimize, jd2html, literate_folder, verify_links

# -----------------------------------------------------------------------------
#
# CONSTANTS
#

"""Big number when we want things to be far."""
const BIG_INT = typemax(Int)

"""Flag for debug mode."""
const DEBUG_MODE = Ref(false)

"""Flag for the initial pass over pages"""
const FULL_PASS = Ref(true)

"""Flag for re-evaluation of all code blocks"""
const FORCE_REEVAL = Ref(false)

"""Flag for error suppression mode (set and unset in optimize only)."""
const SUPPRESS_ERR = Ref(false)

"""Dict to keep track of languages and how comments are indicated and their extensions."""
const CODE_LANG = LittleDict{String,NTuple{2,String}}(
    "c"          => (".c",    "//"),
    "cpp"        => (".cpp",  "//"),
    "fortran"    => (".f90",  "!"),
    "go"         => (".go",   "//"),
    "javascript" => (".js",   "//"),
    "julia"      => (".jl",   "#"),
    "julia-repl" => (".jl",   "#"),
    "lua"        => (".lua",  "--"),
    "matlab"     => (".m",    "%"),
    "python"     => (".py",   "#"),
    "r"          => (".r",    "#"),
    "toml"       => (".toml", "#"),
    )

"""Flag to suppress all statements including eval ones."""
const SILENT_MODE = Ref(false)

"""Extreme negative offset for ordering of newcommands defined in config."""
const OFFSET_GLOB_LXDEFS = Ref(-BIG_INT)

# copied from Base/path.jl
if Sys.isunix()
    """Indicator for directory separation on the OS."""
    const PATH_SEP = "/"
elseif Sys.iswindows()
    const PATH_SEP = "\\"
else
    error("Unhandled OS")
end

"""Type of the containers for page variables (local and global)."""
const PageVars = LittleDict{String,Pair{K,NTuple{N, DataType}} where {K, N}}

"""Relative path to the current file being processed by JuDoc."""
const CUR_PATH = Ref("")

"""Relative path to the last with eval'd code blocks (to know what's in scope)."""
const CUR_PATH_WITH_EVAL = Ref("")

"""Shorter name for a type that we use everywhere"""
const AS = Union{String,SubString{String}}

"""Convenience constants for an automatic message to add to code files."""
const MESSAGE_FILE_GEN     = "This file was generated, do not modify it."
const MESSAGE_FILE_GEN_LIT = "# $MESSAGE_FILE_GEN\n\n"
const MESSAGE_FILE_GEN_JMD = "# $MESSAGE_FILE_GEN # hide\n"

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

include("build.jl") # check if user has Node/minify

# PARSING
include("parser/tokens.jl")
include("parser/ocblocks.jl")
# > latex
include("parser/lx_tokens.jl")
include("parser/lx_blocks.jl")
# > markdown
include("parser/md_tokens.jl")
include("parser/md_validate.jl")
# > html
include("parser/html_tokens.jl")
include("parser/html_blocks.jl")

# CONVERSION
# > markdown
include("converter/md_blocks.jl")
include("converter/md_utils.jl")
include("converter/md.jl")
# > latex
include("converter/lx.jl")
include("converter/lx_simple.jl")
# > html
include("converter/html_functions.jl")
include("converter/html.jl")
# > fighting Julia's markdown parser
include("converter/link_fixer.jl")
# > javascript
include("converter/js_prerender.jl")

# FILE PROCESSING
include("jd_paths.jl")
include("jd_vars.jl")

# FILE AND DIR MANAGEMENT
include("manager/rss_generator.jl")
include("manager/dir_utils.jl")
include("manager/file_utils.jl")
include("manager/judoc.jl")
include("manager/post_processing.jl")

# MISC UTILS
include("misc_utils.jl")
include("misc_html.jl")

# ERROR TYPES
include("error_types.jl")

# INTEGRATION
include("integration/literate.jl")

end # module


text:
# Code snippets from https://github.com/JuliaLang/Compat.jl

function _compat(ex::Expr)
    if VERSION < v"0.7.0-DEV.2562"
        if ex.head == :call && ex.args[1] == :finalizer
            ex.args[2], ex.args[3] = ex.args[3], ex.args[2]
        end
    end
    return Expr(ex.head, map(_compat, ex.args)...)
end

_compat(ex) = ex

macro compat(ex)
    esc(_compat(ex))
end

@static if !isdefined(Base, :pairs)
    pairs(collection) = Base.Generator(=>, keys(collection), values(collection))
    pairs(a::Associative) = a
end


text:
mutable struct RelaxedSequencingNode{T<:Real,U<:Int,V<:Int}<:Node
    uuid                    ::UUID
    length                  ::U
    state                   ::V
    allUp                   ::BitVector
    allDown                 ::BitVector
    someUp                  ::BitVector
    someDown                ::BitVector
    lengthToRoot            ::T
    lengthToTerminal        ::T
    bestSeq                 ::Vector{V}
    exact                   ::Bool
    arcsIn                  ::Vector{ExactArc}
    arcsOut                 ::Vector{ExactArc}
end

#change hashing of nodes to use uuids
Base.hash(node::RelaxedSequencingNode) = hash(node.uuid)
Base.:(==)(node1::RelaxedSequencingNode, node2::RelaxedSequencingNode) = node1.uuid == node2.uuid
Base.length(node::RelaxedSequencingNode) = node.length
#do UUIDs automatically
function RelaxedSequencingNode(pathLength::U,state::V,allUp::BitVector,allDown::BitVector,someUp::BitVector,someDown::BitVector,lengthToRoot::T,lengthToTerminal::T,bestSeq::Vector{V}, exact::Bool,arcsIn::Vector{ExactArc},arcsOut::Vector{ExactArc})where{T<:Real,U<:Int,V<:Int}
    return RelaxedSequencingNode(uuid4(), pathLength,state,allUp,allDown,someUp,someDown,lengthToRoot,lengthToTerminal,bestSeq, exact,arcsIn,arcsOut)
end

"""
    deepCopy(node::RelaxedSequencingNode)

Deep copy a node

# Arguments
- `node::RelaxedSequencingNode`: the node
"""
function deepCopy(node::RelaxedSequencingNode)
    return RelaxedSequencingNode(
        node.uuid, #uuid
        node.length, #path length
        node.state, #current state
        copy(node.allUp), #allUp
        copy(node.allDown), #allDown
        copy(node.someUp), #someUp
        copy(node.someDown), #someDown
        node.lengthToRoot, #to root
        node.lengthToTerminal, #to terminal
        copy(node.bestSeq), #best Seq
        node.exact, #exact
        deepCopy(node.arcsIn), #arcsin
        deepCopy(node.arcsOut) #arcsout
    )
end

Base.show(io::IO,node::RelaxedSequencingNode) = Base.print(
    io,
    "\n   pathLength: ", node.length, "\n",
    "   state: ", node.state, "\n",
    "   allUp: ", node.allUp, "\n",
    "   allDown: ", node.allDown, "\n",
    "   someUp: ", node.someUp, "\n",
    "   someDown: ", node.someDown, "\n",
    "   lengthToRoot: ", node.lengthToRoot, "\n",
    "   lengthToTerminal: ", node.lengthToTerminal, "\n",
    "   best Sequence: ", node.bestSeq, "\n",
    "   isExact: ", node.exact, "\n",
    "   arcsIn: ", showInArcOrigins(node.arcsIn), "\n",
    "   arcsOut: ", showOutArcAssignments(node.arcsOut), "\n",
)

function showOutArcAssignments(arcs::Vector{ExactArc})
    str = ""
    for arc in arcs
        str = string(str, arc.assignment, " ")
    end
    return str
end

function showInArcOrigins(arcs::Vector{ExactArc})
    str = ""
    for arc in arcs
        str = string(str, arc.origin.state, " ")
    end
    return str
end

"""
    createRootRelaxedSequencingNodeFromModel(model::SequencingModel)

Create a root node for a relaxed DD from the sequencing model

# Arguments
- `model::SequencingModel`: the sequencing problem to be evaluated
"""
function createRootRelaxedSequencingNodeFromModel(model::SequencingModel)
    restrictedNode = createRestrictedRootNode(model)
    relaxedNode = RelaxedSequencingNode(
        length(restrictedNode.path), #length
        last(restrictedNode.path), #current state
        falses(length(model)), #allUp
        copy(restrictedNode.visited), #allDown
        falses(length(model)), #someUp
        copy(restrictedNode.visited), #someDown
        0, #to root
        0, #to terminal
        copy(restrictedNode.path),
        true, #exact
        Vector{ExactArc}(), #arcsin
        Vector{ExactArc}() #arcsout
    )
    return relaxedNode, restrictedNode.domain
end

"""
    createPartialRelaxedNode(;length::T=0, state::U=-1, allUp::BitVector=BitVector(), allDown::BitVector=BitVector(), someUp::BitVector=BitVector(), someDown::BitVector=BitVector(), bestSeq::Vector{U}=Vector{Int}(), lengthToRoot::V=0, lengthToTerminal::V=0, exact::Bool=false, arcsIn::Vector{ExactArc}= Vector{ExactArc}(), arcsOut::Vector{ExactArc}= Vector{ExactArc}())where{T<:Int, U<:Int, V<:Real}

Create a node using default values for any combination of fields

# Arguments
- `length::T`: the length field
- `state::U`: the state field
- `allUp::BitVector`: the allUp field
- `allDown::BitVector`: the allDown field
- `someUp::BitVector`: the someUp field
- `someDown::BitVector`: the someDown field
- `bestSeq::Vector{Int}`: the bestSeq field
- `lengthToRoot::V`: the lengthToRoot field
- `lengthToTerminal::V`: the lengthToTerminal field
- `exact::Bool`: the exact field
- `arcsIn::Vector{ExactArc}`: the arcsIn field
- `arcsOut::Vector{ExactArc}`: the arcsOut field
"""
function createPartialRelaxedNode(;length::T=0, state::U=-1, allUp::BitVector=BitVector(), allDown::BitVector=BitVector(), someUp::BitVector=BitVector(), someDown::BitVector=BitVector(), bestSeq::Vector{U}=Vector{Int}(), lengthToRoot::V=0, lengthToTerminal::V=0, exact::Bool=false, arcsIn::Vector{ExactArc}= Vector{ExactArc}(), arcsOut::Vector{ExactArc}= Vector{ExactArc}())where{T<:Int, U<:Int, V<:Real}
    return RelaxedSequencingNode(
        length, #length
        state, #current state
        allUp, #allUp
        allDown, #allDown
        someUp, #someUp
        someDown, #someDown
        lengthToRoot, #to root
        lengthToTerminal, #to terminal
        bestSeq,
        exact, #exact
        arcsIn, #arcsin
        arcsOut #arcsout
    )
end

"""
    getParents(node::RelaxedSequencingNode)

Find and return the parents of a node

# Arguments
- `node::RelaxedSequencingNode`: the node whose parents to find
"""
function getParents(node::RelaxedSequencingNode)
    parents = Set{RelaxedSequencingNode}()
    for arc in node.arcsIn
        push!(parents, arc.origin)
    end
    return parents
end

"""
    getChildren(node::RelaxedSequencingNode)

Find and return the children
# Arguments
- `node::RelaxedSequencingNode`: the node whose children to find
"""
function getChildren(node::RelaxedSequencingNode)
    children = Set{RelaxedSequencingNode}()
    for arc in node.arcsOut
        push!(children, arc.destination)
    end
    return children
end

"""
    updateNodeDownVariables!(node::RelaxedSequencingNode,model::SequencingModel)

Update the state variable for allDown and someDown and lengthToRoot and exact

# Arguments
- `node::RelaxedSequencingNode`: the node to update
- `model::SequencingModel`: the sequencing problem being evaluated
"""
function updateNodeDownVariables!(node::RelaxedSequencingNode,model::SequencingModel)
    node.allDown = falses(length(model))
    node.someDown = falses(length(model))
    node.exact = true

    node.someDown[node.state] = true
    node.allDown[node.state] = true
    #update the best path variables
    bestArc, node.lengthToRoot = findBestArc(model, node.arcsIn)
    node.bestSeq = push!(copy(bestArc.origin.bestSeq), node.state)
    #union domains with parents and check for exactness
    for parent in getParents(node)
        node.someDown = node.someDown .| parent.someDown
        node.allDown = node.allDown .| parent.allDown
        if !parent.exact
            node.exact = false
        end
    end
    #check each all down to make sure its valid
    for i in 1:length(node.allDown)
        if node.allDown[i]
            for arc in node.arcsIn
                if arc.assignment != i
                    if !arc.origin.allDown[i]
                        node.allDown[i]=false
                        break
                    end
                end
            end
        end
    end
    #final check for exactness (first part is just a check to see if the second part should run)
    if node.exact && (node.someDown != node.allDown || node.length!=sum(node.allDown))
        node.exact = false
    end
end

"""
    updateNodeUpVariables!(node::RelaxedSequencingNode,model::SequencingModel)

Update the state variable for someUp, allUp, and lengthToTerminal

# Arguments
- `node::RelaxedSequencingNode`: the node to update
- `model::SequencingModel`: the sequencing problem being evaluated
"""
function updateNodeUpVariables!(node::RelaxedSequencingNode,model::SequencingModel)
    node.allUp = falses(length(model))
    node.someUp = falses(length(model))
    node.lengthToTerminal = node.arcsOut[1].value + node.arcsOut[1].destination.lengthToTerminal
    for arc in node.arcsOut
        node.someUp[arc.assignment] = true
        node.allUp[arc.assignment] = true
        #update the length to terminal
        if isBetterSolutionValue(model, arc.value + arc.destination.lengthToTerminal, node.lengthToTerminal)
            node.lengthToTerminal = arc.value + arc.destination.lengthToTerminal
        end
    end

    #union domains with children
    for child in getChildren(node)
        node.someUp = node.someUp .| child.someUp
        node.allUp = node.allUp .| child.allUp
    end
    #check each all up to make sure its valid
    for i in 1:length(node.allUp)
        if node.allUp[i]
            for arc in node.arcsOut
                if arc.assignment != i
                    if !arc.destination.allUp[i]
                        node.allUp[i]=false
                        break
                    end
                end
            end
        end
    end
end

"""
    createSquareStartDD(model::SequencingModel;rootNode::Union{Nothing, RelaxedSequencingNode}=nothing, masterDomain::BitVector=BitVector())

Create an initial dd from the model where all the arcs have exact values

# Arguments
- `model::SequencingModel`: the sequencing problem being evaluated
- `rootNode::Union{Nothing, RelaxedSequencingNode}`: a root node if not starting from the model
- `masterDomain::BitVector`: the domain of the root node if there is a root node, must be included if root node is included
"""
function createSquareStartDD(model::SequencingModel;rootNode::Union{Nothing, RelaxedSequencingNode}=nothing, masterDomain::BitVector=BitVector())
    if isnothing(rootNode)
        rootNode, masterDomain = createRootRelaxedSequencingNodeFromModel(model)
    end
    relaxedDD = Vector{Vector{RelaxedSequencingNode}}()
    push!(relaxedDD, [rootNode])

    #for each layer going down
    terminalNode = createPartialRelaxedNode(
        length = length(model),
        state=0,
        allUp=falses(length(masterDomain)),
        allDown=falses(length(masterDomain)),
        someUp=falses(length(masterDomain)),
        someDown=falses(length(masterDomain))
    )
    if hasEnd(model)
        lastLayerLoopIndex = length(model)-length(rootNode)
        terminalNode.state = model.seqEnd
    else
        lastLayerLoopIndex = length(model)-length(rootNode)+1
    end
    for layer in 2:lastLayerLoopIndex
        push!(relaxedDD,Vector{RelaxedSequencingNode}())
        #new node starts as a state and length
        for (index, inDomain) in enumerate(masterDomain)
            if inDomain
                newNode = createPartialRelaxedNode(length = layer, state=index)
                push!(relaxedDD[layer],newNode)
            end
        end
        #update domain based on precedence
        if model.precedenceConstraints
            for (i, inDomain) in enumerate(masterDomain)
                if !inDomain && (!hasStart(model) || i!=model.seqStart)
                    #add satisfied precedence constraints to the domain
                    if layer >= get(model.precedenceNumbers,i,0)
                        masterDomain[i] = true
                    end
                end
            end
        end
    end#layer iteration
    push!(relaxedDD, [terminalNode])

    #add all arcs except terminal arcs and update down values
    for layer in 1:length(relaxedDD)-2
        for parent in relaxedDD[layer]
            if layer != 1
                updateNodeDownVariables!(parent,model)
            end
            for child in relaxedDD[layer+1]
                if hasKey(model.objective, parent.state, child.state) && !parent.allDown[child.state]
                    newArc = ExactArc(parent, child, child.state, evaluateDecision(model.objective, parent.state, child.state))
                    push!(parent.arcsOut, newArc)
                    push!(child.arcsIn, newArc)
                end
            end
        end
        toRemove = Vector{RelaxedSequencingNode}()
        for child in relaxedDD[layer+1]
            if isempty(child.arcsIn)
                push!(toRemove, child)
            end
        end
        filter!(x->!(x in toRemove), relaxedDD[layer+1])
    end
    #add arcs to terminal
    if hasEnd(model)
        for parent in relaxedDD[length(relaxedDD)-1]
            if hasKey(model.objective, parent.state, model.seqEnd)
                newArc = ExactArc(parent, terminalNode, model.seqEnd, evaluateDecision(model.objective, parent.state, model.seqEnd))
                push!(parent.arcsOut, newArc)
                push!(terminalNode.arcsIn, newArc)
            end
            updateNodeDownVariables!(parent,model)
        end
    else
        for parent in relaxedDD[length(relaxedDD)-1]
            newArc = ExactArc(parent, terminalNode, -1, 0)
            push!(parent.arcsOut, newArc)
            push!(terminalNode.arcsIn, newArc)
            updateNodeDownVariables!(parent,model)
        end
    end
    updateNodeDownVariables!(terminalNode,model)

    #update all the up variables
    for layer in Iterators.reverse(1:length(relaxedDD)-1)
        for node in relaxedDD[layer]
            updateNodeUpVariables!(node,model)
        end
    end

    return relaxedDD
end


text:
@testset "The Gr(4,6)" begin
  import Base: zeros, one
  import Random: seed!
  seed!(42); #set seed -> at least always the same random numbers.
  M = Grassmannian(4, 6)
  x = GrPoint([1. 0. 0. 0.; 0. 1. 0. 0.; 0. 0. 1. 0.; 0. 0. 0. 1.; 0. 0. 0. 0.; 0. 0. 0. 0.])
  yanti = GrPoint([1. 0. 0. 0.; 0. 0. 0. 0.; 0. 1. 0. 0.; 0. 0. 1. 0.; 0. 0. 0. 1.; 0. 0. 0. 0.])
  y = randomMPoint(M)
  z = randomMPoint(M)
  η = randomTVector(M,x)
  ξ = log(M, x ,y)
  μ = log(M, x, z)
  y2 = exp(M, x, ξ)
  z2 = exp(M, x, μ)
  x2 = addNoise(M,x,0.5)
  x3 = addNoise(M,x,500.)
  w = randomMPoint(M)
  s = rand(Float64)
  μ2 = log(M,x,w)
  w2 = exp(M,x,μ2)

  # Test Grassmannian
  @test_throws ErrorException Grassmannian(6, 4)
  # Test Dimension
  @test manifoldDimension(x)==8
  @test manifoldDimension(M)==8
  # Test Minus operator
  @test getValue( (-η) ) == getValue(-η)
  # Test Cauchy-Schwarz-Inequation
  @test dot(M,x,ξ,η)^2 <= norm(M,x,ξ)^2*norm(M,x,η)^2
  @test dot(M,x,ξ,η) <= norm(M,x,ξ)*norm(M,x,η)
  # Test dot
  @test dot(M,x,ξ,η) ≈ dot(M,x,η,ξ) atol = 10. ^(-16)
  @test dot(M,x,μ,η+ξ) ≈ dot(M,x,μ,η) + dot(M,x,μ,ξ) atol = 10. ^(-15)
  @test dot(M,x,μ,s*η) ≈ s*dot(M,x,μ,η) atol = 10. ^(-15)
  # Test norm
  @test norm(M,x,μ)^2 ≈ dot(M,x,μ,μ) atol = 10. ^(-12)
  @test norm(M,x,μ) - sqrt(dot(M,x,μ,μ)) ≈ 0 atol = eps()
  # Test distance
  @test distance(M,x,x) ≈ 0 atol = 10. ^(-7)
  @test distance(M,x,z) - norm(M,x,μ) ≈ 0 atol = 10. ^(-7)
  @test distance(M,x,yanti) ≈ pi/2 atol = 10. ^(-16)
  # Test randomMPoint and randomTVector
  @test norm( getValue(x)'*getValue(η) )≈ 0 atol = 10. ^(-6)
  # Test addNoise
  @test norm( getValue(x2)'*getValue(x2) - one(getValue(x2)'*getValue(x2)) ) ≈ 0 atol = 10. ^(-15)
  @test norm( getValue(x3)'*getValue(x3) - one(getValue(x3)'*getValue(x3)) ) ≈ 0 atol = 10. ^(-15)
  # Test Exponential and logarithm
  @test norm( getValue(log(M, x, y2)) - getValue(ξ)) ≈ 0 atol = 10. ^(-5)
  @test norm( getValue(log(M, x, w2)) - getValue(μ2) ) ≈ 0 atol =  10. ^(-13)
  @test norm( getValue(log(M, x, y)) - getValue(ξ) ) ≈ 0 atol = 10.0^(-16)
  @test norm( getValue(log(M, x, z)) - getValue(μ) ) ≈ 0 atol = 10.0^(-16)
  @test distance(M, exp(M,x,log(M,x,y)) , y ) ≈ 0 atol = 5*10. ^(-8)
  @test norm( getValue(log(M,x,exp(M,x,ξ))) - getValue(ξ) ) ≈ 0 atol = 10. ^(-15)
  @test_throws ErrorException log(M,x,yanti)
  #Test retraction
  @test norm( getValue(inverseRetraction(M,x,retraction(M,x,ξ))) - getValue(ξ)) ≈ 0 atol = 10.0^(-14)
  @test norm( transpose(getValue(retraction(M,x,ξ))) * getValue(retraction(M,x,ξ)) - one(transpose(getValue(x))*getValue(x)) ) ≈ 0 atol = 10.0^(-14)
  # Test parallelTransport
  @test norm(getValue(parallelTransport(M,x,z,η)) - projection(M,z,getValue(η))) ≈ 0 atol = 10.0^(-16)
  # Test zeroTVector
  @test norm(M,x,zeroTVector(M,x)) ≈ 0 atol = 10.0^(-16)
  # Test validateMPoint and validateTVector
  ynot1 = GrPoint([1. 0. 0. 0.; 0. 1. 0. 0.; 0. 0. 1. 0.; 0. 0. 0. 1.; 0. 0. 0. 0.; 0. 0. 0. 0.;  0. 0. 0. 0.])
  ynot2 = GrPoint([1. 0. 0. 0. 0.; 0. 1. 0. 0. 0.; 0. 0. 1. 0. 0.; 0. 0. 0. 1. 0.; 0. 0. 0. 0. 0.; 0. 0. 0. 0. 0.])
  ynot3 = GrPoint([1. 0. 0. 0.; 0. 1. 0. 0.; 0. 0. 1. 0.; 0. 0. 0. 0.; 0. 0. 0. 0.; 0. 0. 0. 0.])
  ynot4 = GrPoint([1. 0. 0. 0.; 0. 2. 0. 0.; 0. 0. 4. 0.; 0. 0. 0. 8.; 0. 0. 0. 0.; 0. 0. 0. 0.])
  ξnot1 = GrTVector([1. 0. 0. 0.; 0. 1. 0. 0.; 0. 0. 1. 0.; 0. 0. 0. 1.; 0. 0. 0. 0.; 0. 0. 0. 0.;  0. 0. 0. 0.])
  ξnot2 = GrTVector([1. 0. 0. 0. 0.; 0. 1. 0. 0. 0.; 0. 0. 1. 0. 0.; 0. 0. 0. 1. 0.; 0. 0. 0. 0. 0.; 0. 0. 0. 0. 0.])
  ξnot3 = GrTVector([ 1. 2. 3. 0.; 4. 5. 6. 0.; 7. 8. 9. 0.; 0. 0. 0. 0.; 1. 0. 0. 0.; 0. 1. 0. 0.])
  @test_throws ErrorException validateMPoint(M,ynot1)
  @test_throws ErrorException validateMPoint(M,ynot2)
  @test_throws ErrorException validateMPoint(M,ynot3)
  @test_throws ErrorException validateMPoint(M,ynot4)
  @test_throws ErrorException validateTVector(M,x,ξnot1)
  @test_throws ErrorException validateTVector(M,x,ξnot2)
  @test_throws ErrorException validateTVector(M,x,ξnot3)

  #
  # Complex manifold
  #
  N = Grassmannian{Complex{Float64}}(4,6)
  xcompl = randomMPoint(N)
  ycompl = randomMPoint(N)
  zcompl = randomMPoint(N)
  ηcompl = randomTVector(N, xcompl)
  ωcompl = randomTVector(N, xcompl)
  ξcompl = randomTVector(N, xcompl)
  νcompl = randomTVector(N,xcompl)
  wcompl = randomMPoint(N)
  x2compl = addNoise(N,xcompl,0.5)
  x3compl = addNoise(N,xcompl,500.0)

  # Test Minus operator
  @test getValue( (-ηcompl) ) == getValue(-ηcompl)
  # Test Cauchy-Schwarz-Inequation
  @test abs(dot(N,xcompl,ωcompl,ηcompl)^2) <= abs(norm(N,xcompl,ωcompl)^2*norm(N,xcompl,ηcompl)^2)
  @test dot(N,xcompl,ωcompl,ηcompl) <= norm(N,xcompl,ωcompl)*norm(N,xcompl,ηcompl)
  # Test dot
  @test dot(N,xcompl,ωcompl,ηcompl) ≈ dot(N,xcompl,ηcompl,ωcompl) atol = 10.0^(-16)
  @test dot(N,xcompl,ωcompl,ηcompl+ξcompl) ≈ dot(N,xcompl,ωcompl,ηcompl) + dot(N,xcompl,ωcompl,ξcompl) atol = 10.0^(-15)
  @test dot(N,xcompl,ωcompl,s*ηcompl) ≈ s*dot(N,xcompl,ωcompl,ηcompl) atol = 10.0^(-14)
  # Test norm
  @test norm(N,xcompl,ωcompl)^2 ≈ dot(N,xcompl,ωcompl,ωcompl) atol = 10.0^(-14)
  @test norm(N,xcompl,ωcompl) ≈ sqrt(dot(N,xcompl,ωcompl,ωcompl)) atol = 10.0^(-14)
  # Test randomMPoint and randomTVector
  @test norm( getValue(wcompl)'*getValue(wcompl) - one(getValue(wcompl)'*getValue(wcompl)) ) ≈ 0 atol = 10.0^(-14)
  @test norm( getValue(xcompl)'*getValue(νcompl) + (getValue(xcompl)'*getValue(νcompl))' ) ≈ 0 atol = 10.0^(-5)
  # Test addNoise
  @test norm( getValue(x2compl)'*getValue(x2compl) - one(getValue(x2compl)'*getValue(x2compl)) ) ≈ 0 atol=10.0^(-5)
  @test norm( getValue(x3compl)'*getValue(x3compl) - one(getValue(x3compl)'*getValue(x3compl)) ) ≈ 0 atol=10.0^(-6)

  #Test manifoldDimension
  @test manifoldDimension(N) == manifoldDimension(wcompl)
  @test manifoldDimension(N) == 16
  @test manifoldDimension(xcompl) == 16
end


text:
cd(@__DIR__); allNGWPs = false; include("setups/toronto.jl")
using Plots.PlotMeasures
gr(dpi = 200)

## (a) pedestrian volume graph signal
fp = load("../datasets/new_toronto.jld", "fp")
G_Sig.f = reshape(fp, (N, 1))
gplot(A, X; width=1); plot!(size = (600, 600), right_margin = 5mm)
signal_plt = scatter_gplot!(X; marker = fp, plotOrder = :s2l, ms = 3)
savefig(signal_plt, "../figs/Toronto_fp.png")

## (b) pedestrian signal relative l2 approximation error by various methods
# uncomment the following if `allNGWPs` = true
# DVEC = getall_expansioncoeffs2(G_Sig, GP_dual, GP_dual_Lsym, VM_NGWP, PC_NGWP, LP_NGWP,
#                                VM_NGWP_Lsym, PC_NGWP_Lsym, LP_NGWP_Lsym, 𝚽, 𝚽sym)
# use precomputed results
DVEC = load("../datasets/Toronto_fp_DVEC.jld", "DVEC")
approx_error_plot2(DVEC);
plt = plot!(xguidefontsize = 14, yguidefontsize = 14, legendfontsize = 11, size = (600, 600))
savefig(plt, "../figs/Toronto_fp_DAG_approx.png")


text:
using Pkg;
using DelimitedFiles, DataFrames
using LinearAlgebra, Statistics, JuMP, Ipopt, Juniper, Optim
using Gadfly, Cairo, Compose;
cd(@__DIR__) #src

include("func_synth.jl")
dat, ~ = readdlm("MLAB_data.txt", header = true);

# Index to group the data: Control vs. treated ============================
idx_control = 1:38;
idx_treat   = 39;

# Group data ==============================================================
# Endogenous time-invariant predictor -------------------------------------
idx_predictor = 1:7
xₒ = dat[idx_predictor, idx_control];         # Control group
x₁ = dat[idx_predictor, idx_treat];           # Treatment group

# Rescale predictors
global xₒ, x₁ = normalz(xₒ, x₁, idx_control, idx_treat);

# Outcome data to evaluate ------------------------------------------------
idx_outcome = 8:38;
yₒ = dat[idx_outcome, idx_control];           # Control group
y₁ = dat[idx_outcome, idx_treat];             # Treatment group

idx_pretreat = 1:19;                          # Rows of pre-treat period
global yₒpre = yₒ[idx_pretreat, :];
global y₁pre = y₁[idx_pretreat, :];

# Optimization
s = std(hcat(xₒ, x₁)', dims = 1)';
s2 = s[2:end];
s1 = s[1];
ν₀ = ((s1./s2).^2);

# Optimization to estimate weights ========================================
#https://julianlsolvers.github.io/Optim.jl/stable/#user/minimization/
#https://julianlsolvers.github.io/Optim.jl/stable/#algo/samin/
q = (length(ν₀))
opt = optimize(loss_ssr, fill(0.0, q), fill(5e8, q), ν₀,
                Fminbox(BFGS()), Optim.Options(iterations=100) # or 100_000
            );
# NOTE: Requires further improvement, as the optimization takes a long time.
# Can try with different optimization methods
opt = optimize(loss_ssr, fill(0.0, q), fill(5e8, q), ν₀,
                  SAMIN(), Optim.Options(iterations=1_000_000) # or 100_000
            );

# Save results
ν = vcat(1, Optim.minimizer(opt));

#MATLAB optimal solutions - to check the accuracy of solveQadProg(H, f)
#=
ν = [1,
      780.438610055690,
      24909.6348724115,
      37806.0921140682,
      124.042538397502,
      60413.0041399961,
      63630.9996910117];
=#
# Recover weights: w-weights
D = (eye(length(ν)) .* ν);

H = xₒ' * D * xₒ;
f = -1 .* (x₁' * D * xₒ);

ω = abs.(solveQadProg(H, f)); # we know this works now.

# Build counterfactual output.
ŷₒ = yₒ * ω;
ȳₒ = mean(yₒ, dims = 2); # Res of US, for plotting purpose

# plot ====================================================================
yr = 1970:1:2000;
plt_yavg = layer(x = yr, y = y₁, Geom.line,
                  Theme(default_color=colorant"deepskyblue",
                        line_style = [:solid]),
                  xintercept = [1988],
                        Geom.vline(color = "black", style=[[0.5mm,0.5mm]])
            );
plt_y    = layer(x = yr, y = ȳₒ, Geom.line,
                  Theme(default_color=colorant"red", line_style = [:solid]));
plt_yhat = layer(x = yr, y = ŷₒ, Geom.line, Theme(line_style = [:dash]));
plt_fin = plot(plt_yavg, plt_y, plt_yhat,
            Guide.xlabel("Year"),
            Guide.ylabel("Somking per capita (in packs)"),
            Guide.annotation(compose(context(),
                        Compose.text(1988.5, 25, "← Proposition 99"))),
            Guide.annotation(compose(context(),
                        Compose.text(1993, 110, "Rest of U.S."))),
            Guide.annotation(compose(context(),
                        Compose.text(1971, 95, "Solid: California (actual) \n Dash: California (synthetic)")))
      )
pngout = SVG("plt_fin.svg", 5inch, 4inch)
draw(pngout, plt_fin)


text:
# This package is a work in progress for the estimation and identification of Vector Autoregressive (VAR) models.
# Credits:
# Kilian and Kim 2011, Cremfi codes

module VectorAutoregressions
using Parameters

type Intercept end

type VAR
    Y::Array
    X::Array
    β::Array
    ϵ::Array
    Σ::Array
    p::Int64
    inter::Intercept
    VAR(Y,X,β,ϵ,Σ,p,inter) = p <= 0 ? error("Lag-length error: 'p' must be strictly positive") : new(Y,X,β,ϵ,Σ,p,inter)
end

function VAR(y::Array,p::Int64,i::Bool)
    i == false ? ((Y,X,β,ϵ,Σ,p) = fit(y,p)) : ((Y,X,β,ϵ,Σ,p) = fit(y,p,Intercept()))
    return VAR(Y,X,β,ϵ,Σ,p,Intercept())
end

abstract type CI end

type IRFs
    IRF::Array
    CI::CI
end

type CI_asy <: CI
    CIl::Array
    CIh::Array
end

type CI_boot <: CI
    CIl::Array
    CIh::Array
end

function IRFs_a(V::VAR,H::Int64,i::Bool)
    if i == true
        mVar1 = get_VAR1_rep(V,V.inter)
        mIRF = irf_chol(V, mVar1, H)
        mStd = irf_ci_asymptotic(V, H, V.inter)
    else
        mVar1 = get_VAR1_rep(V)
        mIRF = irf_chol(V, mVar1, H)
        mStd = irf_ci_asymptotic(V, H)
    end
    mCIl = mIRF - 1.96.*mStd
    mCIh = mIRF + 1.96.*mStd
    return IRFs(mIRF,CI_asy(mCIl,mCIh))
end

function IRFs_b(V::VAR,H::Int64,nrep::Int64,i::Bool)
    if i == true
        mVar1 = get_VAR1_rep(V,V.inter)
        mIRF = irf_chol(V, mVar1, H)
        CI = irf_ci_bootstrap(V, H, nrep,V.inter)
    else
        mVar1 = get_VAR1_rep(V)
        mIRF = irf_chol(V, mVar1, H)
        CI = irf_ci_bootstrap(V, H, nrep)
    end
    return IRFs(mIRF,CI)
end

function fit(y::Array,p::Int64)
    (T,K) = size(y)
    T < K && error("error: there are more covariates than observation")
    X = y
    y = transpose(y)
    Y = y[:,p+1:T]
    X = lagmatrix(X,p)'
    β = (Y*X')/(X*X')
    ϵ = Y - β*X
    Σ = ϵ*ϵ'/(T-p-p*K-1)
    return Y,X,β,ϵ,Σ,p
end

function fit(y::Array,p::Int64,inter::Intercept)
    (T,K) = size(y)
    T < K && error("error: there are more covariates than observation")
    X = y
    y = transpose(y)
    Y = y[:,p+1:T]
    X = lagmatrix(X,p,inter)'
    β = (Y*X')/(X*X')
    ϵ = Y - β*X
    Σ = ϵ*ϵ'/(T-p-p*K-1)
    return Y,X,β,ϵ,Σ,p
end

function lagmatrix{F}(x::Array{F},p::Int64,inter::Intercept)
    sk = 1
    T, K = size(x)
    k    = K*p+1
    idx  = repmat(1:K, p)
    X    = Array{F}(T-p, k)
    # building X (t-1:t-p) allocating data from D matrix - avoid checking bounds
    for j = 1+sk:(sk+K*p)
        for i = 1:(T-p)
            lg = round(Int, ceil((j-sk)/K)) - 1 # create index [0 0 1 1 2 2 ...etc]
            @inbounds X[i, j] = x[i+p-1-lg, idx[j-sk]]
        end
    end
    for j=1:T-p
        @inbounds X[j,1] = 1.0
    end
    return X
end

function lagmatrix{F}(x::Array{F},p::Int64)
    sk = 1
    T, K = size(x)
    k    = K*p+1
    idx  = repmat(1:K, p)
    X    = Array{F}(T-p, k)
    # building X (t-1:t-p) allocating data from D matrix - avoid checking bounds
    for j = 1+sk:(sk+K*p)
        for i = 1:(T-p)
            lg = round(Int, ceil((j-sk)/K)) - 1 # create index [0 0 1 1 2 2 ...etc]
            @inbounds X[i, j] = x[i+p-1-lg, idx[j-sk]]
        end
    end
    return X[:,2:end]
end

function lagmatrix{F}(x::Vector{F},p::Int64)
    sk = 1
    T = length(x)
    K = 1
    k    = K*p+1
    idx  = repmat(1:K, p)
    X    = Array{F}(T-p, k)
    # building X (t-1:t-p) allocating data from D matrix - avoid checking bounds
    for j = 1+sk:(sk+K*p)
        for i = 1:(T-p)
            lg = round(Int, ceil((j-sk)/K)) - 1 # create index [0 0 1 1 2 2 ...etc]
            @inbounds X[i, j] = x[i+p-1-lg, idx[j-sk]]
        end
    end
    return X[:,2:end]
end

function commutation(n::Int64, m::Int64)
    # returns Magnus and Neudecker's commutation matrix of dimensions n by m
    k = reshape(kron(vec(eye(n)), eye(m)), n*m, n*m)
    return k
end

function duplication(n::Int64)
    # Returns Magnus and Neudecker's duplication matrix of size n
    # VERY AMBIGUOUS FUNC
    a = tril(ones(n,n))::Array{Float64}
    i = find(a)::Vector{Int64}
    a[i] = 1:size(i,1)
    a = a + tril(a,-1)'
    j = convert(Vector{Int64}, vec(a))::Vector{Int64}
    m = trunc.(Int,(n*(n+1)/2))::Int64
    d = zeros(n*n,m)
    for r = 1:size(d,1)
        d[r, j[r]] = 1.0
    end
    return d
end

function elimat(m::Int64)
    # elimat(m) returns the elimination matrix Lm
    # The elimination matrix Lm is for any matrix F, Vech(F)=Lm Vec(F)
    A = eye(m^2)
    L = A[1:m,:]
    for n in 2:m
        S = A[m*(n-1)+1:n*m,:]
        S = S[n:end,:]
        L = [L;S]
    end
    return L
end

function get_lag_length(D::Array, pbar::Integer)
    IC   = zeros(pbar,1)
    for p = 1:pbar
        V = VAR(D,p,false)
        n,m = size(V.X)
        t = n-p
        sig = V.Σ/t
        if ic == "aic"
            IC[p] = log(det(sig))+2*p*m^2/t
        elseif ic == "bic"
            IC[p] = log(det(sig))+(m^2*p)*log(t)/t
        elseif ic == "hqc"
            IC[p] = log(det(sig))+2*log(log(t))*m^2*p/t
        elseif ic == "aicc"
            b = t/(t-(p*m+m+1))
            IC[p] = t*(log(det(sig))+m)+2*b*(m^2*p+m*(m+1)/2)
        elseif error("'ic' must be aic, bic, hqc or aicc")
        end
    end
    length_ic = indmin(IC)
    println("Using $ic the best lag-length is $length_ic")
    return length_ic
end

function get_lag_length_aic(D::Array, pbar::Integer, inter::Intercept)
    IC   = zeros(pbar,1)
    for p = 1:pbar
        V = VAR(D,p,true)
        n,m = size(V.X)::Tuple{Int64,Int64}
        t = convert(Float64,n-p)
        sig = (V.Σ./t)::Array{Float64}
        IC[p] = log(det(sig))+2*p*m^2/t
    end
    length_ic = indmin(IC)
    println("The best lag-length is $length_ic")
    return length_ic
end

function get_lag_length_aic(D::Array, pbar::Integer)
    IC   = zeros(pbar,1)
    for p = 1:pbar
        V = VAR(D,p,false)
        n,m = size(V.X)::Tuple{Int64,Int64}
        t = convert(Float64,n-p)
        sig = (V.Σ./t)::Array{Float64}
        IC[p] = log(det(sig))+2*p*m^2/t
    end
    length_ic = indmin(IC)
    println("The best lag-length is $length_ic")
    return length_ic
end

function get_lag_length_bic(D::Array, pbar::Integer, inter::Intercept)
    IC   = zeros(pbar,1)
    for p = 1:pbar
        V = VAR(D,p,true)
        n,m = size(V.X)::Tuple{Int64,Int64}
        t = convert(Float64,n-p)
        sig = (V.Σ./t)::Array{Float64}
        IC[p] = log(det(sig))+(m^2*p)*log(t)/t
    end
    length_ic = indmin(IC)
    println("The best lag-length is $length_ic")
    return length_ic
end

function get_lag_length_bic(D::Array, pbar::Integer)
    IC   = zeros(pbar,1)
    for p = 1:pbar
        V = VAR(D,p,false)
        n,m = size(V.X)::Tuple{Int64,Int64}
        t = convert(Float64,n-p)
        sig = (V.Σ./t)::Array{Float64}
        IC[p] = log(det(sig))+(m^2*p)*log(t)/t
    end
    length_ic = indmin(IC)
    println("The best lag-length is $length_ic")
    return length_ic
end

function get_lag_length_hqc(D::Array, pbar::Integer, inter::Intercept)
    IC   = zeros(pbar,1)
    for p = 1:pbar
        V = VAR(D,p,true)
        n,m = size(V.X)::Tuple{Int64,Int64}
        t = convert(Float64,n-p)
        sig = (V.Σ./t)::Array{Float64}
        IC[p] = log(det(sig))+2*log(log(t))*m^2*p/t
    end
    length_ic = indmin(IC)
    println("The best lag-length is $length_ic")
    return length_ic
end

function get_lag_length_hqc(D::Array, pbar::Integer)
    IC   = zeros(pbar,1)
    for p = 1:pbar
        V = VAR(D,p,false)
        n,m = size(V.X)::Tuple{Int64,Int64}
        t = convert(Float64,n-p)
        sig = (V.Σ./t)::Array{Float64}
        IC[p] = log(det(sig))+2*log(log(t))*m^2*p/t
    end
    length_ic = indmin(IC)
    println("The best lag-length is $length_ic")
    return length_ic
end

function get_lag_length_aicc(D::Array, pbar::Integer, inter::Intercept)
    IC   = zeros(pbar,1)
    for p = 1:pbar
        V = VAR(D,p,true)
        n,m = size(V.X)::Tuple{Int64,Int64}
        t = convert(Float64,n-p)
        sig = (V.Σ./t)::Array{Float64}
        b = t/(t-(p*m+m+1))
        IC[p] = t*(log(det(sig))+m)+2*b*(m^2*p+m*(m+1)/2)
    end
    length_ic = indmin(IC)
    println("The best lag-length is $length_ic")
    return length_ic
end

function get_lag_length_aicc(D::Array, pbar::Integer)
    IC   = zeros(pbar,1)
    for p = 1:pbar
        V = VAR(D,p,false)
        n,m = size(V.X)::Tuple{Int64,Int64}
        t = convert(Float64,n-p)
        sig = (V.Σ./t)::Array{Float64}
        b = t/(t-(p*m+m+1))
        IC[p] = t*(log(det(sig))+m)+2*b*(m^2*p+m*(m+1)/2)
    end
    length_ic = indmin(IC)
    println("The best lag-length is $length_ic")
    return length_ic
end

function get_VAR1_rep(V::VAR)
    K = size(V.Σ,1)
    B = vcat(V.β, hcat(eye(K*(V.p-1)), zeros(K*(V.p-1),K))::Array{Float64,2})::Array{Float64,2}
    B = convert(Array{Float64,2},B)
end

function get_VAR1_rep(V::VAR,inter::Intercept)
    K = size(V.Σ,1)
    # v = [V.β[:,1]; zeros(K*(V.p-1),1)]
    B = vcat(V.β[:,2:end], hcat(eye(K*(V.p-1)), zeros(K*(V.p-1),K)))::Array{Float64,2}
    B = convert(Array{Float64,2},B)
end

function irf_ci_asymptotic(V::VAR, H::Int64)
    K,T = size(V.Y)::Tuple{Int64,Int64}
    SIGa = kron(inv(V.X*V.X'/(T-V.p)),V.Σ)
    # Calculation of stdev follows Lutkepohl(2005) p.111,93
    A = get_VAR1_rep(V)
    A0inv = cholfact(V.Σ)[:L]
    STD   = zeros(K^2,H+1)
    COV2   = zeros(K^2,H+1)
    J = [eye(K) zeros(K,K*(V.p-1))]
    L = elimat(K)
    Kk = commutation(K,K)
    Hk = L'/(L*(eye(K^2)+Kk)*kron(A0inv,eye(K))*L')
    D = duplication(K)
    Dplus = (D'*D)\D'
    SIGsig = 2*Dplus*kron(V.Σ,V.Σ)*Dplus'
    Cbar0 = kron(eye(K),J*eye(K*V.p)*J')*Hk
    STD[:,1] = vec((reshape(diag(real(sqrt.(complex(Cbar0*SIGsig*Cbar0'/(T-V.p))))),K,K))')
    COV2[:,1] = vec((reshape(diag((Cbar0*SIGsig*Cbar0'/(T-V.p))),K,K))')
    for h=1:H
        Gi = zeros(K^2,K^2*V.p)
        for m=0:(h-1)
            Gi += kron(J*(A')^(h-1-m),J*(A^m)*J')
        end
        C = kron(A0inv',eye(K))*Gi
        Cbar = kron(eye(K),J*A^h*J')*Hk
        STD[:,h+1] = vec((reshape(diag(real(sqrt.(complex(C*SIGa*C'+Cbar*SIGsig*Cbar')/(T-V.p)))),K,K))')
        COV2[:,h+1] = vec((reshape(diag(((Cbar*SIGsig*Cbar')/(T-V.p))),K,K))')
    end
    return STD
end

function irf_ci_asymptotic(V::VAR, H::Int64, inter::Intercept)
    K,T = size(V.Y)::Tuple{Int64,Int64}
    SIGa = kron(inv(V.X*V.X'/(T-V.p)),V.Σ)
    SIGa = SIGa[K+1:end,K+1:end]
    # Calculation of stdev follows Lutkepohl(2005) p.111,93
    A = get_VAR1_rep(V,V.inter)
    A0inv = cholfact(V.Σ)[:L]
    STD   = zeros(K^2,H+1)
    COV2   = zeros(K^2,H+1)
    J = [eye(K) zeros(K,K*(V.p-1))]
    L = elimat(K)
    Kk = commutation(K,K)
    Hk = L'/(L*(eye(K^2)+Kk)*kron(A0inv,eye(K))*L')
    D = duplication(K)
    Dplus = (D'*D)\D'
    SIGsig = 2*Dplus*kron(V.Σ,V.Σ)*Dplus'
    Cbar0 = kron(eye(K),J*eye(K*V.p)*J')*Hk
    STD[:,1] = vec((reshape(diag(real(sqrt.(complex(Cbar0*SIGsig*Cbar0'/(T-V.p))))),K,K))')
    COV2[:,1] = vec((reshape(diag((Cbar0*SIGsig*Cbar0'/(T-V.p))),K,K))')
    for h=1:H
        Gi = zeros(K^2,K^2*V.p)
        for m=0:(h-1)
            Gi += kron(J*(A')^(h-1-m),J*(A^m)*J')
        end
        C = kron(A0inv',eye(K))*Gi
        Cbar = kron(eye(K),J*A^h*J')*Hk
        STD[:,h+1] = vec((reshape(diag(real(sqrt.(complex(C*SIGa*C'+Cbar*SIGsig*Cbar')/(T-V.p)))),K,K))')
        COV2[:,h+1] = vec((reshape(diag(((Cbar*SIGsig*Cbar')/(T-V.p))),K,K))')
    end
    return STD
end

get_boot_init_int_draw(T::Int64,p::Int64) = Int64(trunc.((T-p+1)*rand()+1))
get_boot_init_vector_draw(T::Int64,p::Int64) = Array{Int64}(trunc.((T-p)*rand(T-p)+1))

function build_sample(V::VAR)
    K,T = size(V.Y)::Tuple{Int64,Int64}
    # Draw block of initial pre-sample values
    y = zeros(K,T)                                       # bootstrap data
    u = zeros(K,T)                             # bootstrap innovations
    iDraw = get_boot_init_int_draw(T,V.p)            # position of initial draw
    y[:,1:V.p] = V.Y[:,iDraw:iDraw+V.p-1]                   # drawing pre-sample obs
    # Draw innovations
    vDraw = get_boot_init_vector_draw(T,V.p)    # index for innovation draws
    u[:, V.p+1:T] = V.ϵ[:,vDraw]                  # drawing innovations
    @inbounds for i = V.p+1:T
        y[:,i] = u[:,i]
        for j =  1:V.p
            y[:,i] += V.β[:,(j-1)*K + 1:j*K]*y[:,i-j]
        end
    end
    return y
end

function build_sample(V::VAR,inter::Intercept)
    K,T = size(V.Y)::Tuple{Int64,Int64}
    # Draw block of initial pre-sample values
    y = zeros(K,T)                                       # bootstrap data
    u = zeros(K,T)                             # bootstrap innovations
    iDraw = get_boot_init_int_draw(T,V.p)            # position of initial draw
    y[:,1:V.p] = V.Y[:,iDraw:iDraw+V.p-1]                   # drawing pre-sample obs
    # Draw innovations
    vDraw = get_boot_init_vector_draw(T,V.p)    # index for innovation draws
    u[:, V.p+1:T] = V.ϵ[:,vDraw]                  # drawing innovations
    @inbounds for i = V.p+1:T
        y[:,i] = V.β[:,1] + u[:,i]
        for j =  1:V.p
            y[:,i] += V.β[:,(j-1)*K + 2:j*K+1]*y[:,i-j]
        end
    end
    return y
end

col_mean(x::Array) = mean(x,2)
test_bias_correction(x::Array) =  any(abs.(eigvals(x)).<1)

function get_boot_ci(V::VAR,H::Int64,nrep::Int64, bDo_bias_corr::Bool,inter::Intercept)
    K,T = size(V.Y)::Tuple{Int64,Int64}
    mIRFbc = zeros(nrep, K^2*(H+1))
    @inbounds for j = 1:nrep
        # Recursively construct sample
        yr = build_sample(V,V.inter)
        yr = (yr .- col_mean(yr))'  # demean yr bootstrap data
        #pr = V.p # also using lag length selection
        Vr = VAR(yr,V.p,true)
        # Bias correction: if the largest root of the companion matrix
        # is less than 1, do BIAS correction
        mVar1 = get_VAR1_rep(Vr,V.inter)::Array{Float64,2}
        bBias_corr_test = test_bias_correction(mVar1)
        if all([bDo_bias_corr, bBias_corr_test])
            mVar1 = bias_correction(Vr,mVar1)
        end
        mIRF = irf_chol(Vr,mVar1,H)
        mIRFbc[j,:] = vec(mIRF')'
    end                     # end bootstrap
    return mIRFbc
end

function get_boot_ci(V::VAR,H::Int64,nrep::Int64, bDo_bias_corr::Bool)
    K,T = size(V.Y)::Tuple{Int64,Int64}
    mIRFbc = zeros(nrep, K^2*(H+1))
    @inbounds for j = 1:nrep
        # Recursively construct sample
        yr = build_sample(V)
        yr = (yr .- col_mean(yr))'  # demean yr bootstrap data
        #pr = V.p # also using lag length selection
        Vr = VAR(yr,V.p,false)
        # Bias correction: if the largest root of the companion matrix
        # is less than 1, do BIAS correction
        mVar1 = get_VAR1_rep(Vr)::Array{Float64,2}
        bBias_corr_test = test_bias_correction(mVar1)
        if all([bDo_bias_corr, bBias_corr_test])
            mVar1 = bias_correction(Vr,mVar1)
        end
        mIRF = irf_chol(Vr,mVar1,H)
        mIRFbc[j,:] = vec(mIRF')'
    end                     # end bootstrap
    return mIRFbc
end

function get_companion_vcv(V::VAR)
    K,T = size(V.Y)::Tuple{Int64,Int64}
    mSigma = vcat(hcat(V.Σ, zeros(K,K*V.p-K)), zeros(K*V.p-K,K*V.p))::Array{Float64,2}
    mSigma = convert(Array{Float64,2},mSigma)
end

function get_sigma_y(V::VAR, mVar1::Array, mSigma::Array)
    K,T = size(V.Y)::Tuple{Int64,Int64}
    vSigma = (eye((K*V.p)^2)::Array{Float64,2}-kron(mVar1,mVar1)::Array{Float64,2})\vec(mSigma)::Vector{Float64}    # Lutkepohl p.29 (2.1.39)
    mSigma_y = reshape(vSigma, K*V.p::Int64, K*V.p::Int64)::Array{Float64,2}
    return convert(Array{Float64,2},mSigma_y)
end

function get_bias(mSigma::Array,mSigma_y::Array,B::Array,I::Array,mSum_eigen::Array)
    return mBias = mSigma*(inv(I - B) + B/(I-B^2) + mSum_eigen)/(mSigma_y)
end

function bias_correction(V::VAR,mVar1::Array)
    # Bias-correction Pope (1990)
    K,T = size(V.Y)::Tuple{Int64,Int64}
    mSigma = get_companion_vcv(V)
    mSigma_y = get_sigma_y(V,mVar1,mSigma)
    I = eye(K*V.p, K*V.p)
    B = mVar1'
    vEigen = eigvals(mVar1)
    mSum_eigen = zeros(K*V.p,K*V.p)
    for h = 1:K*V.p
        mSum_eigen += vEigen[h].\(I - vEigen[h]*B)
    end
    mBias = get_bias(mSigma,mSigma_y,B,I,mSum_eigen)
    mAbias = -mBias/T
    return mBcA = get_proportional_bias_corr(mVar1,mAbias)
end

function get_proportional_bias_corr(mVar1::Array,Abias::Array;iBc_stab::Int64 = 9, iδ::Int64 = 1)
    mBcA = zeros(size(mVar1))
    while iBc_stab >= 1
        # Adjust bias-correction proportionately
        mBcA = real(mVar1-iδ*Abias)
        bcmod = abs.(eigvals(mBcA))
        any(bcmod.>= 1) ?  iBc_stab = 1 : iBc_stab = 0
        iδ += -0.01
        iδ <= 0 && (iBc_stab = 0)
    end
    return mBcA
end

function get_boot_conf_interval(mIRFbc::Array,H::Int64,K::Int64)
    N = size(mIRFbc,2)
    mCILv = zeros(1,N)
    mCIHv = zeros(1,N)
    for i = 1:N
        mCILv[:,i] = quantile(vec(mIRFbc[:,i]),0.025)
        mCIHv[:,i] = quantile(vec(mIRFbc[:,i]),0.975)
    end
    mCIL  = reshape(mCILv',H+1,K^2)'
    mCIH  = reshape(mCIHv',H+1,K^2)'
    return mCIL, mCIH
end

function irf_ci_bootstrap(V::VAR, H::Int64, nrep::Int64; bDo_bias_corr::Bool=true)
    K,T = size(V.Y)::Tuple{Int64,Int64}
    iScale_ϵ = sqrt((T-V.p)/(T-V.p-K*V.p-1))
    u = V.ϵ*iScale_ϵ   # rescaling residual (Stine, JASA 1987)
    mIRFbc = get_boot_ci(V,H,nrep,bDo_bias_corr)
    # Calculate 95 perccent interval endpoints
    mCIL, mCIH = get_boot_conf_interval(mIRFbc::Array,H::Int64,K::Int64)
    return CI_boot(mCIL,mCIH)
end

function irf_ci_bootstrap(V::VAR, H::Int64, nrep::Int64, inter::Intercept; bDo_bias_corr::Bool=true)
    K,T = size(V.Y)::Tuple{Int64,Int64}
    iScale_ϵ = sqrt((T-V.p)/(T-V.p-K*V.p-1))
    u = V.ϵ*iScale_ϵ   # rescaling residual (Stine, JASA 1987)
    mIRFbc = get_boot_ci(V,H,nrep,bDo_bias_corr,V.inter)
    # Calculate 95 perccent interval endpoints
    mCIL, mCIH = get_boot_conf_interval(mIRFbc::Array,H::Int64,K::Int64)
    return CI_boot(mCIL,mCIH)
end

function irf_chol(V::VAR, mVar1::Array, H::Int64)
    K = size(V.Σ,1)
    mSigma = cholfact(V.Σ)[:L]::LowerTriangular  # Cholesky or reduced form
    J = [eye(K,K) zeros(K,K*(V.p-1))]
    mIRF = zeros(K^2,H+1)
    mIRF[:,1] = reshape((J*mVar1^0*J'*mSigma)',K^2,1)::Array
    for i = 1:H
        mIRF[:,i+1] = reshape((J*mVar1^i*J'*mSigma)',K^2,1)::Array
    end
    return mIRF
end

function irf_reduce_form(V::VAR, mVar1::Array, H::Int64)
    K = size(V.Σ,1)
    mSigma = eye(K,K)
    J = [eye(K,K) zeros(K,K*(V.p-1))]
    mIRF = zeros(K^2,H+1)
    mIRF[:,1] = reshape((J*mVar1^0*J'*mSigma)',K^2,1)::Array
    for i = 1:H
        mIRF[:,i+1] = reshape((J*mVar1^i*J'*mSigma)',K^2,1)::Array
    end
    return mIRF
end

function t_test(V::VAR)
    K = size(V.Σ,1)
    Kp = size(V.X,1)
    H = kron(inv(V.X*V.X'),V.Σ)
    SE = sqrt.(diag(H))
    T = reshape(vec(V.β)./SE,K,Kp)
    return T
end

function gen_var1_data!(y::Array,mR::Array,mP,burnin::Int64)
    T,K = size(y)
     for j = 2:T
         y[j,:] = mR*y[j-1,:] + mP*randn(K,1)
     end
     y = y[burnin+1:end,:]
     return y .- mean(y,1)
 end

export VAR, IRFs_a, IRFs_b, gen_var1_data!
end # end of the module


text:
function mark_field(r)
    side=Ost
    mark_row!(r,side)

    while isborder(r,side)==false
        move!(r,Nord)
        side=inverse(side)
        mark_row!(r,side)
    end
end

function mark_row(r::Robot,side::HorizonSide)
    putmarker!(r)
    while move_if_posible!(r,side) == true
        putmarker!(r)
    end
end
function move_if_possible!(r::Robot, direct_side::HorizonSide)::Bool
    orthogonal_side = next(direct_side)
    reverse_side = reverse(orthogonal_side)
    num_steps=0
    while isborder(direct_side) == true
        if isborder(r, orthogonal_side) == false
            move(r, orthogonal_side)
            num_steps += 1
        else
            break
        end
    end

    if isborder(r,direct_side) == false
        while isborder(r,reverse_side) == true
            move!(r,direct_side)
        end
        result = true
    else
        result = false
    end
    move!(r,reverse_side)
    return result
end

movements!(r::Robot, side::HorizonSide, num_steps::Int) =
for _ in 1:num_steps
    move_if_posible!(r,side)
end


text:
using LightGraphs
using MetaGraphs
using DataFrames
using GraphDataFrameBridge
using Test

import LightGraphs.SimpleGraphs: SimpleGraph, SimpleDiGraph
testdir = dirname(@__FILE__)

testgraphs(g) = [g, SimpleGraph{UInt8}(g), SimpleGraph{Int16}(g)]
testdigraphs(g) = [g, SimpleDiGraph{UInt8}(g), SimpleDiGraph{Int16}(g)]

tests = [
    "graphdataframebridge"
]

@testset "GraphDataFrameBridge" begin
    for t in tests
        tp = joinpath(testdir, "$(t).jl")
        include(tp)
    end
end


text:
# This is the simplest example of full waveform inversion. It is
# built in the 3D poststack domain using a simple convolution, which
# highlights the theory of FWI, but is inaccurate on real data.
#
# by GRAM | 23 Apr 2016

print("\n3D Poststack full waveform inversion algorithm\n\n
In this program FWI updates a velocity model one sample at a time,
one trace at a time. The velocity sample is perturbed positively,
negatively, and not at all. Then the algorithm models new seismic
from all three of the perturbed velocity models, and compares them to
the input (real) seismic data. Whichever perturbation generates the
least amount of error is chosen to be correct, and input into the
updated velocity model.
This algorithm accepts six inputs in the following order:\n
	1. Initial velocity model
	2. Poststack 3D seismic volume in Seis format
	3. One dimensional wavelet in Seis format
	4. Output file name (updated vel model in Seis format)
	5. Velocity update increment percentage in decimal
	6. Maximum number of velocity update iterations
	7. Verbosity (see note)\n
If verbose is 0 operation is silent. If verbose is 1 updates will
print to stdout. If verbose is 2 debugging info will print to stdout.
Here is an example to get you started on the syntax:\n\n
fwi(\"../dat/vel_incorrect\", \"../dat/image_correct\", \"../dat/wav\", \"../dat/updated_vel\", .1, 5, 0)\n\n
Please note this algorithm is written for clarity not speed, so
use forgivingly small velocity models.\n\n")

using Seismic
include("compare_wavelet.jl")	# wavefield comparison algorithm
include("model_conv.jl")	# seismic modeling algorithm

function fwi(vel="",sei="",wav="",out="",inc=.1,itr_max=5,verbose=0)

	# read initial vel model, seismic image, and wavelet
	vel, vel_h, ex = SeisRead(vel)
	sei, sei_h = SeisRead(sei)
	wav, wav_h = SeisRead(wav)

	# initialize
	nz = size(vel)[1]
	nx = size(vel)[2]
	ny = size(vel)[3]
	dz = vel_h[1].d1
	dx = 20
	dy = 20
	thrs = 0.01		# change in vel threshold
	wvln = size(wav)[1]	# wavelength of the wavelet

	#### start FWI ##########################
	for i = 1 : ny, j = 1 : nx

		# tell the user where we are in the volume
		if verbose > 0
			print("\nNow updating trace x = ", j, ", y = ", i, " \n")
		end

		# calc velocity intervals
		layer = [0]
		for kk = 2 : nz
			current = vel[kk - 1, j, i]
			if vel[kk, j, i] > current + thrs
				append!(layer, [kk - 1])
			end
		end

		if verbose == 2 # debugging information
			print("\nlayer = ", layer,"\n")
		end

		# work over all vel intervals
		for k = 1 : size(layer)[1] - 1

			# iteratively update the velocity model
			for itr = 1 : itr_max

				# build perturbed velocity models
				v_p = vel[:,j,i]
				v_p[layer[k] + 1 : layer[k + 1]] = v_p[layer[k] + 1] * (1 + inc)

				v_n = vel[:,j,i]
				v_n[layer[k] + 1 : layer[k + 1]] = v_n[layer[k] + 1] * (1 - inc)

				v_o = vel[:,j,i]

				# model perturbed seismic data
				s_p, s_n, s_o = model(v_p, v_n, v_o, nx, ny, nz, wav)

				# compare images
				comp = compare(s_p, s_n, s_o, sei, i, j, k, wvln, layer, verbose)

				# move in the direction of improvement
				dir = 0
				dir = find(comp .== minimum(comp))[1]

				if verbose == 2 # debugging information
					print("\n",comp,"\n")
					print("\n",dir,"\n")
				end

				if dir == 1
					vel[:,j,i] = v_p
				elseif dir == 2
					vel[:,j,i] = v_n
				elseif dir == 3
					break
				else
					print("\n diff symmetric @", k,"\n")
				end

			end # end of vel iterations

		end # end of vel intervals

	end #### end of FWI ###########################

	# write out the FWI updated velocity model
	SeisWrite(out,vel,vel_h,ex)

end # end of function


text:
module SalsaTest

using Test
using Salsa

# The Salsa.jl test is configurable to run with arbitrary Runtime instances, by
# providing your own implementation of this function before running `"test/Salsa.jl"`.
new_test_rt() = Runtime()
new_test_rt(ctx::Context) where Context = Runtime{Context}(ctx)

@testset "Salsa.jl" begin
    include("Salsa.jl")
end

end # SalsaTest


text:
# TODO
@testset "size" begin
    for S in (LinSRange(1,3,3),
              StepSRangeLen(1,1,3),
              StepSRange(1,1,3),
              UnitSRange(1,3),
              OneToSRange(3)
             )
        @test Size(typeof(S)) === Size{(3,)}()
    end
end