text:
# import Pkg; Pkg.activate(@__DIR__);
# Pkg.instantiate()

# using RigidBodyDynamics
# using StaticArrays
# using RobotDynamics
# using Altro
# using TrajectoryOptimization

# function mpc_update(altro, prob_mpc, Z_track, t0, k_mpc)
#     TrajectoryOptimization.set_initial_time!(prob_mpc, t0)

#     # Propagate the system forward w/ noise
#     # x0 = discrete_dynamics(TrajectoryOptimization.integration(prob_mpc),
#     #                             prob_mpc.model, prob_mpc.Z[1])
#     x0 = rk4(prob_mpc.model, state(prob_mpc.Z[1]), control(prob_mpc.Z[1]), prob_mpc.Z[1].dt)

#     TrajectoryOptimization.set_initial_state!(prob_mpc, x0)

#     # Update tracking cost
#     TrajectoryOptimization.update_trajectory!(prob_mpc.obj, Z_track, k_mpc)

#     # Shift the initial trajectory
#     RobotDynamics.shift_fill!(prob_mpc.Z)

#     # Shift the multipliers and penalties
#     Altro.shift_fill!(TrajectoryOptimization.get_constraints(altro.solver_al))
# end

# function update_constraints!(conSet::Altro.ALConstraintSet, cons::Altro.ALConstraintSet, start=1)
#     for c = 1:length(conSet)
#         prob_cons = conSet.convals[c]
#         inds = (start-1) .+ (1:length(conSet.convals[c].λ))
#         if (start + length(conSet.convals[c]) - 1) <= length(cons)
#             for (i,k) in enumerate(inds)
#                 conSet.convals[c].λ[i] = cons.convals[k].λ
#                 # TrajectoryOptimization.set_LQR_goal!(obj[i], state(Z[k]), control(Z[k]))
#             end
#         else
#             for (i,k) in enumerate(inds)
#                 conSet.convals[c].λ[i] = cons[min(k,length(cons))]
#                 # TrajectoryOptimization.set_LQR_goal!(obj[i], state(Z[min(k,length(Z))]), control(Z[min(k,length(Z))]))
#             end
#         end
#     end
# end

# function update_trajectory_fill!(obj::Objective{<:TrajectoryOptimization.QuadraticCostFunction}, Z::TrajectoryOptimization.AbstractTrajectory, start=1)
#     inds = (start-1) .+ (1:length(obj))
#     if (start + length(obj) - 1) <= length(Z)
#         for (i,k) in enumerate(inds)
#             TrajectoryOptimization.set_LQR_goal!(obj[i], state(Z[k]), control(Z[k]))
#         end
#     else
#         for (i,k) in enumerate(inds)
#             TrajectoryOptimization.set_LQR_goal!(obj[i], state(Z[min(k,length(Z))]), control(Z[min(k,length(Z))]))
#         end
#     end
# end

# function update_guess_trajectory!(prob_Z::Traj, Z::TrajectoryOptimization.AbstractTrajectory, start=1)
#     inds = (start-1) .+ (1:length(prob_Z))
#     if (start + length(prob_Z) - 1) <= length(Z)
#         for (i,k) in enumerate(inds)
#             prob_Z[i] = Z[k]
#             # TrajectoryOptimization.set_LQR_goal!(obj[i], state(Z[k]), control(Z[k]))
#         end
#     else
#         for (i,k) in enumerate(inds)
#             prob_Z[i] = Z[min(k,length(Z))]
#             # TrajectoryOptimization.set_LQR_goal!(obj[i], state(Z[min(k,length(Z))]), control(Z[min(k,length(Z))]))
#         end
#     end
# end

# function mpc_update(altro, prob_mpc, Z_track, cons, x0, t0)
#     TrajectoryOptimization.set_initial_time!(prob_mpc, t0)

#     # Propagate the system forward w/ noise
#     # x0 = discrete_dynamics(TrajectoryOptimization.integration(prob_mpc),
#     #                             prob_mpc.model, prob_mpc.Z[1])
#     # x0 = rk4(prob_mpc.model, state(prob_mpc.Z[1]), control(prob_mpc.Z[1]), prob_mpc.Z[1].dt)
#     k_mpc = argmin(norm.([(states(Z_track)[k] - x0) for k=1:length(Z_track)]))
#     println(k_mpc)

#     TrajectoryOptimization.set_initial_state!(prob_mpc, x0)

#     # Update tracking cost
#     update_trajectory_fill!(prob_mpc.obj, Z_track, k_mpc)

#     # Shift the initial trajectory
#     # update_guess_trajectory!(prob_mpc.Z, Z_track, k_mpc)
#     RobotDynamics.shift_fill!(prob_mpc.Z)

#     # Shift the multipliers and penalties
#     # update_constraints!(TrajectoryOptimization.get_constraints(altro.solver_al), cons, k_mpc)
#     Altro.shift_fill!(TrajectoryOptimization.get_constraints(altro.solver_al), 0)
# end

function ArthurProblem(Xref; params=MPC_Params())
    # Create initial guess of U trajectory
    Qref = Vector{Float64}[]
    Uref = Vector{Float64}[]
    N = length(Xref)
    n = params.n
    m = params.m
    tf = round(params.dt*(N-1), digits = 3)
    Qref = [Xref[1] for k=1:length(Xref)]
    set_configuration!(params.state, Xref[1][1:7])
    Uref = [inverse_dynamics(params.state, params.v̇) for k=1:(length(Xref)-1)]

    # push!(Qref, copy(Xref[1][1:14]))
    # for k = 1:N-1
    #     set_configuration!(params.state, Qref[k][1:7])
    #     set_velocity!(params.state, Qref[k][8:14])
    #     push!(Uref, inverse_dynamics(params.state, params.v̇))
    #     push!(Qref, rk4(params.model, Qref[k], Uref[k], params.dt))
    # end

    # Gather trajectory information
    # Create Empty ConstraintList
    conSet = ConstraintList(n,m,N)

    # Control Bounds based on Robot Specs (Joint torque limits)
    u_bnd = [39.0, 39.0, 39.0, 39.0, 9.0, 9.0, 9.0]
    control_bnd = BoundConstraint(n,m, u_min=-u_bnd, u_max=u_bnd)
    add_constraint!(conSet, control_bnd, 1:N-1)

    # State Bounds based on Robot Specs (Joint velocity and speed limits)
    x_bnd = zeros(n)
    x_bnd[1:7] = [Inf, deg2rad(128.9), Inf, deg2rad(147.8), Inf, deg2rad(120.3), Inf] # rad
    x_bnd[8:14] = [1.39, 1.39, 1.39, 1.39, 1.22, 1.22, 1.22] # rad/sec
    # x_bnd[8:14] = [0.25, 0.05, 0.05, 0.05, 0.25, 0.15, 0.15] # rad/sec
    # x_bnd[15:end] = [Inf for k=1:(n-14)] # Constraints on force elsewhere
    state_bnd = BoundConstraint(n,m, x_min=-x_bnd, x_max=x_bnd)
    add_constraint!(conSet, state_bnd, 1:N-1)

    # # Cartesian Velocity Bound
    # ẋ_max = 0.0005 # m/s
    # vel_bnd = NormConstraint(n, m, ẋ_max, Inequality(), 21:23)
    # add_constraint!(conSet, vel_bnd, 1:N)

    # # Force Bound (Fx Fy Fz)
    # F_max = 20 # Newtons
    # F_bnd = NormConstraint(n, m, F_max, Inequality(), 27:29)
    # add_constraint!(conSet, F_bnd, 1:N)

    dtref = [params.dt for k=1:N]
    traj = RobotDynamics.Traj(Xref, Uref, dtref, cumsum(dtref) .- dtref[1])
    obj = TrajectoryOptimization.TrackingObjective(params.Q, params.R, traj, Qf=params.Qf)
    prob = Problem(params.model, obj, Xref[end], tf, x0=Xref[1], constraints=conSet, X0=Qref, U0=Uref, integration=RK3)
    return prob
end

# function ArthurHorizonProblem(prob::TrajectoryOptimization.Problem, x0, N; start=1, params=MPC_Params())
#     # H = N
#     while (start+N-1) > length(prob.Z)
#         N -= 1
#     end

#     n,m = size(prob)
#     dt = prob.Z[1].dt
#     tf = (N-1)*dt
#     # Get sub-trajectory
#     if N > 1
#     # if N == H
#         Z = Traj(prob.Z[start:start+N-1])
#     # else
#         # Z = Traj([prob.Z[start:start+N-1]; [prob.Z[end] for k = 1:(H-N+1)]])
#     # end

#     # Generate a cost that tracks the trajectory
#         obj = TrajectoryOptimization.TrackingObjective(params.Q, params.R, Z, Qf=params.Qf)

#     # Use the same constraints, except the Goal constraint
#         cons = ConstraintList(n,m,N)
#         for (inds, con) in zip(prob.constraints)
#             if !(con isa GoalConstraint)
#                 if inds.stop > N
#                     inds = inds.start:N-(prob.N - inds.stop)
#                 end
#                 length(inds) > 0 && TrajectoryOptimization.add_constraint!(cons, con, inds)
#             end
#         end

#         prob = TrajectoryOptimization.Problem(prob.model, obj, state(Z[end]), tf, x0=x0, constraints=cons,
#             integration=TrajectoryOptimization.integration(prob)
#         )
#         initial_trajectory!(prob, Z)
#     else
#         obj = TrajectoryOptimization.LQRObjective(params.Q, params.R, params.Qf, state(prob.Z[end]), params.LQRH)
#          # # Gather trajectory information
#         # Create Empty ConstraintList
#         cons = ConstraintList(params.n,params.m,params.LQRH)

#         # Control Bounds based on Robot Specs (Joint torque limits)
#         u_bnd = [39.0, 39.0, 39.0, 39.0, 9.0, 9.0, 9.0]
#         control_bnd = BoundConstraint(params.n,params.m, u_min=-u_bnd, u_max=u_bnd)
#         add_constraint!(cons, control_bnd, 1:params.LQRH-1)

#         # State Bounds based on Robot Specs (Joint velocity and speed limits)
#         x_bnd = zeros(params.n)
#         x_bnd[1:7] = [Inf, deg2rad(128.9), Inf, deg2rad(147.8), Inf, deg2rad(120.3), Inf] # rad
#         x_bnd[8:14] = [1.39, 1.39, 1.39, 1.39, 1.22, 1.22, 1.22] # rad/sec
#         # x_bnd[8:14] = [0.25, 0.05, 0.05, 0.05, 0.25, 0.15, 0.15] # rad/sec
#         # x_bnd[15:end] = [Inf for k=1:(n-14)] # Constraints on force elsewhere
#         state_bnd = BoundConstraint(params.n,params.m, x_min=-x_bnd, x_max=x_bnd)
#         add_constraint!(cons, state_bnd, 1:params.LQRH-1)

#         # # # Cartesian Velocity Bound
#         # # ẋ_max = 0.0005 # m/s
#         # # vel_bnd = NormConstraint(n, m, ẋ_max, Inequality(), 21:23)
#         # # add_constraint!(cons, vel_bnd, 1:N)

#         # # # Force Bound (Fx Fy Fz)
#         # # F_max = 20 # Newtons
#         # # F_bnd = NormConstraint(n, m, F_max, Inequality(), 27:29)
#         # # add_constraint!(cons, F_bnd, 1:N)
#         Z = Traj([prob.Z[end] for k = 1:params.LQRH])
#         tf = params.LQRH*dt
#         prob = TrajectoryOptimization.Problem(prob.model, obj, state(prob.Z[end]), tf, x0=x0, constraints=cons,
#             integration=TrajectoryOptimization.integration(prob)
#         )
#         initial_trajectory!(prob, Z)
#     end

#     return prob
# end

# function ArthurHorizonProblem(prob::TrajectoryOptimization.Problem, x0, N; start=1, params=MPC_Params())
#     # H = N
#     # while (start+N-1) > length(prob.Z)
#     #     N -= 1
#     # end

#     n,m = size(prob)
#     dt = prob.Z[1].dt
#     tf = params.LQRH*dt
#     # tf = (N-1)*dt
#     if start >= length(prob.Z)
#         N = length(prob.Z)
#     else
#         N = start + 1
#     end

#     x_current = zeros(n)
#     x_current[1:7] = x0[1:7]
#     xf = zeros(n)
#     xf[1:7] .= state(prob.Z[N])[1:7]

#     obj = TrajectoryOptimization.LQRObjective(params.Q, params.R, params.Qf, xf, params.LQRH)

#     # Create Empty ConstraintList
#     cons = ConstraintList(params.n,params.m,params.LQRH)

#     # Control Bounds based on Robot Specs (Joint torque limits)
#     u_bnd = [39.0, 39.0, 39.0, 39.0, 9.0, 9.0, 9.0]
#     control_bnd = BoundConstraint(params.n,params.m, u_min=-u_bnd, u_max=u_bnd)
#     add_constraint!(cons, control_bnd, 1:params.LQRH-1)

#     # State Bounds based on Robot Specs (Joint velocity and speed limits)
#     x_bnd = zeros(params.n)
#     x_bnd[1:7] = [Inf, deg2rad(128.9), Inf, deg2rad(147.8), Inf, deg2rad(120.3), Inf] # rad
#     x_bnd[8:14] = [1.39, 1.39, 1.39, 1.39, 1.22, 1.22, 1.22] # rad/sec
#     # x_bnd[15:end] = [Inf for k=1:(n-14)] # Constraints on force elsewhere
#     state_bnd = BoundConstraint(params.n,params.m, x_min=-x_bnd, x_max=x_bnd)
#     add_constraint!(cons, state_bnd, 1:params.LQRH-1)

#     # # # Cartesian Velocity Bound
#     # # ẋ_max = 0.0005 # m/s
#     # # vel_bnd = NormConstraint(n, m, ẋ_max, Inequality(), 21:23)
#     # # add_constraint!(cons, vel_bnd, 1:N)

#     # # # Force Bound (Fx Fy Fz)
#     # # F_max = 20 # Newtons
#     # # F_bnd = NormConstraint(n, m, F_max, Inequality(), 27:29)
#     # # add_constraint!(cons, F_bnd, 1:N)

#     Qref = [x_current for k=1:params.LQRH]
#     set_configuration!(params.state, x_current[1:7])
#     Uref = [inverse_dynamics(params.state, params.v̇) for k=1:params.LQRH-1]
#     dtref = [params.dt for k=1:params.LQRH]
#     Z = Traj(Qref, Uref, dtref, cumsum(dtref) .- dtref[1])
#     # Z = Traj([prob.Z[end] for k = 1:params.LQRH])
#     tf = params.LQRH*dt
#     prob = TrajectoryOptimization.Problem(prob.model, obj, xf, tf, x0=x0, constraints=cons,
#         integration=TrajectoryOptimization.integration(prob)
#     )
#     initial_trajectory!(prob, Z)

#     return prob
# end

function ArthurHorizonProblem(Xref, x0, N; start=1, params=MPC_Params())
    # H = N
    # while (start+N-1) > length(prob.Z)
    #     N -= 1
    # end

    n = params.n
    m = params.m
    dt = params.dt
    tf = params.LQRH*dt
    # tf = (N-1)*dt
    if start >= length(Xref)
        N = length(Xref)
    else
        N = start + 1
    end

    x_current = zeros(n)
    x_current[1:7] = x0[1:7]
    xf = zeros(n)
    xf[1:7] .= Xref[N][1:7]
    # xf[1:7] .= x0[1:7]

    obj = TrajectoryOptimization.LQRObjective(params.Q, params.R, params.Qf, xf, params.LQRH)

    # Create Empty ConstraintList
    cons = ConstraintList(params.n,params.m,params.LQRH)

    # Control Bounds based on Robot Specs (Joint torque limits)
    u_bnd = [39.0, 39.0, 39.0, 39.0, 9.0, 9.0, 9.0]
    control_bnd = BoundConstraint(params.n,params.m, u_min=-u_bnd, u_max=u_bnd)
    add_constraint!(cons, control_bnd, 1:params.LQRH-1)

    # State Bounds based on Robot Specs (Joint velocity and speed limits)
    x_bnd = zeros(params.n)
    x_bnd[1:7] = [Inf, deg2rad(128.9), Inf, deg2rad(147.8), Inf, deg2rad(120.3), Inf] # rad
    x_bnd[8:14] = [1.39, 1.39, 1.39, 1.39, 1.22, 1.22, 1.22] # rad/sec
    # x_bnd[15:end] = [Inf for k=1:(n-14)] # Constraints on force elsewhere
    state_bnd = BoundConstraint(params.n,params.m, x_min=-x_bnd, x_max=x_bnd)
    add_constraint!(cons, state_bnd, 1:params.LQRH-1)

    # # # Cartesian Velocity Bound
    # # ẋ_max = 0.0005 # m/s
    # # vel_bnd = NormConstraint(n, m, ẋ_max, Inequality(), 21:23)
    # # add_constraint!(cons, vel_bnd, 1:N)

    # # # Force Bound (Fx Fy Fz)
    # # F_max = 20 # Newtons
    # # F_bnd = NormConstraint(n, m, F_max, Inequality(), 27:29)
    # # add_constraint!(cons, F_bnd, 1:N)

    Qref = [x_current for k=1:params.LQRH]
    set_configuration!(params.state, x_current[1:7])
    Uref = [inverse_dynamics(params.state, params.v̇) for k=1:params.LQRH-1]
    dtref = [params.dt for k=1:params.LQRH]
    Z = Traj(Qref, Uref, dtref, cumsum(dtref) .- dtref[1])
    # Z = Traj([prob.Z[end] for k = 1:params.LQRH])
    tf = params.LQRH*dt
    prob = TrajectoryOptimization.Problem(params.model, obj, xf, tf, x0=x0, constraints=cons,
        integration=RK3
    )
    initial_trajectory!(prob, Z)

    return prob
end

struct MPC_Params
    model::Arthur
    n::Int64
    m::Int64
    H::Int64
    LQRH::Int64
    tf::Float64
    dt::Float64
    state::MechanismState
    v̇::SegmentedVector{JointID, Float64, Base.OneTo{JointID}, Vector{Float64}}
    opts::SolverOptions{Float64}
    Q::Diagonal{Float64, SVector{14, Float64}}
    Qf::Diagonal{Float64, SVector{14, Float64}}
    R::Diagonal{Float64, SVector{7, Float64}}

    function MPC_Params()
        model = Arthur()
        n,m = size(model)

        tf = 0.25 # Time Horizon (seconds)
        dt = 0.05 # Time step (seconds)
        H = Int(round(tf/dt) + 1) # Time Horizon (discrete steps)
        LQRH = 5
        state = MechanismState(model.mechanism)
        zero!(state)
        v̇ = similar(velocity(state))
        for joint in joints(state.mechanism)
            if joint == joints(state.mechanism)[2] || joint == joints(state.mechanism)[3] || joint == joints(state.mechanism)[4] || joint == joints(state.mechanism)[5] || joint == joints(state.mechanism)[6] || joint == joints(state.mechanism)[7] || joint == joints(state.mechanism)[8]
                v̇[joint][1] = 0.0
            end
        end

        opts = SolverOptions(
            cost_tolerance_intermediate=1e-2,
            penalty_scaling=10.,
            penalty_initial=1.0,
            show_summary=true,
            projected_newton=true,
        )

        Q = 100.0*Diagonal(@SVector ones(n))
        Qf = 1000.0*Diagonal(@SVector ones(n))
        R = 1.0e-1*Diagonal(@SVector ones(m))
        new(model,n,m,H,LQRH,tf,dt,state,v̇,opts,Q,Qf,R)
    end
end

# function initialize_solver(params::MPC_Params)
#     traj = RobotDynamics.Traj(params.n, params.m, params.dt, params.H)
#     obj = TrajectoryOptimization.TrackingObjective(params.Q, params.R, traj, Qf=params.Qf)
#     prob = Problem(params.model, obj, Xref[end], tf, x0=Xref[1], constraints=params.conSet, X0=Xref, U0=Uref)
#     altro = ALTROSolver(prob, params.opts)
# end


text:
include("manchester.jl")

# String (UTF-8) to Binary: Hello, World!
test_string = "01001000011001010110110001101100011011110010110000100000010101110110111101110010011011000110010000100001"
println("\nGround truth:\n$test_string")

# The argument "method" is optional and specifies if you
# wan't some other encoding method, and in this case "invert".
# Remember the method when you are decoding the data, or else
# you will get the wrong data back

method = "invert"
println("\nEncoding method:\n$method")
encoded = encode(test_string, method)
println("\nEncoded data:\n$(encoded.data)")

# Here is the decoded data
decoded = decode(string(encoded.data, "11"), encoded.method)
println("\nDecoded data:\n$(decoded.data)")

# For reference
println("\nGround truth (again for reference):\n$test_string")

ext: jl

max_stars_repo_name: Irreq/manchester-coding

max_issues_repo_name: Irreq/manchester-coding

max_forks_repo_name: Irreq/manchester-coding


text:
module TimeSeriesResampler

export resample, ohlc, mean, sum
export TimeFrame, Begin, End

include("resample.jl")

end # module


text:
# # [Time-Domain Simulation of a Wind Turbine Blade](@id wind-turbine-blade)
#
#-
#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`wind-turbine-blade.ipynb`](@__NBVIEWER_ROOT_URL__/examples/wind-turbine-blade.ipynb).
#-

using GXBeam, LinearAlgebra

L = 60 # m

## create points
nelem = 10
x = range(0, L, length=nelem+1)
y = zero(x)
z = zero(x)
points = [[x[i],y[i],z[i]] for i = 1:length(x)]

## index of endpoints of each beam element
start = 1:nelem
stop = 2:nelem+1

## stiffness matrix for each beam element
stiffness = fill(
    [2.389e9  1.524e6  6.734e6 -3.382e7 -2.627e7 -4.736e8
     1.524e6  4.334e8 -3.741e6 -2.935e5  1.527e7  3.835e5
     6.734e6 -3.741e6  2.743e7 -4.592e5 -6.869e5 -4.742e6
    -3.382e7 -2.935e5 -4.592e5  2.167e7 -6.279e5  1.430e6
    -2.627e7  1.527e7 -6.869e5 -6.279e5  1.970e7  1.209e7
    -4.736e8  3.835e5 -4.742e6  1.430e6  1.209e7  4.406e8],
    nelem)

## mass matrix for each beam element
mass = fill(
    [258.053      0.0        0.0      0.0      7.07839  -71.6871
       0.0      258.053      0.0     -7.07839  0.0        0.0
       0.0        0.0      258.053   71.6871   0.0        0.0
       0.0       -7.07839   71.6871  48.59     0.0        0.0
       7.07839    0.0        0.0      0.0      2.172      0.0
     -71.6871     0.0        0.0      0.0      0.0       46.418],
     nelem)

## create assembly of interconnected nonlinear beams
assembly = Assembly(points, start, stop; stiffness=stiffness, mass=mass)

## prescribed conditions
prescribed_conditions = (t) -> begin
    Dict(
        ## fixed left side
        1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0),
        ## force on right side
        nelem+1 => PrescribedConditions(Fz = 1e5*sin(20*t))
    )
end

## simulation time
t = 0:0.001:2.0

system, history, converged = time_domain_analysis(assembly, t;
    prescribed_conditions = prescribed_conditions)

#!jl nothing #hide

# We can visualize tip displacements and the resultant forces in the root by accessing
# the post-processed results for each time step contained in the variable `history`.
# Note that the root resultant forces for this case are equal to the external
# forces/moments, but with opposite sign.

#md using Suppressor #hide

#md @suppress_err begin #hide

using Plots
pyplot()

point = vcat(fill(nelem+1, 6), fill(1, 6))
field = [:u, :u, :u, :theta, :theta, :theta, :F, :F, :F, :M, :M, :M]
direction = [1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3]
ylabel = ["\$u_x\$ (\$m\$)", "\$u_y\$ (\$m\$)", "\$u_z\$ (\$m\$)",
    "Rodriguez Parameter \$\\theta_x\$ (degree)",
    "Rodriguez Parameter \$\\theta_y\$ (degree)",
    "Rodriguez Parameter \$\\theta_z\$ (degree)",
    "\$F_x\$ (\$N\$)", "\$F_y\$ (\$N\$)", "\$F_z\$ (\$N\$)",
    "\$M_x\$ (\$Nm\$)", "\$M_y\$ (\$Nm\$)", "\$M_z\$ (\$N\$)"]
#nb ph = Vector{Any}(undef, 12)

for i = 1:12

    local y

#nb    ph[i] = plot(
#!nb    plot(
        xlim = (0, 2.0),
        xticks = 0:0.5:2.0,
        xlabel = "Time (s)",
        ylabel = ylabel[i],
        grid = false,
        overwrite_figure=false
        )
    y = [getproperty(state.points[point[i]], field[i])[direction[i]]
        for state in history]

    if field[i] == :theta
        ## convert to Rodriguez parameter
        @. y = 4*atan(y/4)
        ## convert to degrees
        @. y = rad2deg(y)
    end

    if field[i] == :F || field[i] == :M
        y = -y
    end

    plot!(t, y, label="")
#!nb     plot!(show=true)
#md     savefig("../assets/wind-turbine-blade-"*string(field[i])*string(direction[i])*".svg"); #hide
#md     closeall() #hide
end

#md end #hide
#md nothing #hide

#nb ph[1]
#nb #-
#nb ph[2]
#nb #-
#nb ph[3]
#nb #-
#nb ph[4]
#nb #-
#nb ph[5]
#nb #-
#nb ph[6]
#nb #-
#nb ph[7]
#nb #-
#nb ph[8]
#nb #-
#nb ph[9]
#nb #-
#nb ph[10]
#nb #-
#nb ph[11]
#nb #-
#nb ph[12]
#nb #-

#md # ![](../assets/wind-turbine-blade-u1.svg)
#md # ![](../assets/wind-turbine-blade-u2.svg)
#md # ![](../assets/wind-turbine-blade-u3.svg)
#md # ![](../assets/wind-turbine-blade-theta1.svg)
#md # ![](../assets/wind-turbine-blade-theta2.svg)
#md # ![](../assets/wind-turbine-blade-theta3.svg)
#md # ![](../assets/wind-turbine-blade-F1.svg)
#md # ![](../assets/wind-turbine-blade-F2.svg)
#md # ![](../assets/wind-turbine-blade-F3.svg)
#md # ![](../assets/wind-turbine-blade-M1.svg)
#md # ![](../assets/wind-turbine-blade-M2.svg)
#md # ![](../assets/wind-turbine-blade-M3.svg)

# These plots are identical to those presented by Qi Wang, Wenbin Yu, and Michael A.
# Sprague in "Geometric Nonlinear Analysis of Composite Beams Using Wiener-Milenkovic
# Parameters".
#
#-
#
# We can also visualize the time history of the system using ParaView.
#
root_chord = 1.9000
tip_chord =  0.4540
airfoil = [ # MH-104
    1.00000000  0.00000000;
    0.99619582  0.00017047;
    0.98515158  0.00100213;
    0.96764209  0.00285474;
    0.94421447  0.00556001;
    0.91510964  0.00906779;
    0.88074158  0.01357364;
    0.84177999  0.01916802;
    0.79894110  0.02580144;
    0.75297076  0.03334313;
    0.70461763  0.04158593;
    0.65461515  0.05026338;
    0.60366461  0.05906756;
    0.55242353  0.06766426;
    0.50149950  0.07571157;
    0.45144530  0.08287416;
    0.40276150  0.08882939;
    0.35589801  0.09329359;
    0.31131449  0.09592864;
    0.26917194  0.09626763;
    0.22927064  0.09424396;
    0.19167283  0.09023579;
    0.15672257  0.08451656;
    0.12469599  0.07727756;
    0.09585870  0.06875796;
    0.07046974  0.05918984;
    0.04874337  0.04880096;
    0.03081405  0.03786904;
    0.01681379  0.02676332;
    0.00687971  0.01592385;
    0.00143518  0.00647946;
    0.00053606  0.00370956;
    0.00006572  0.00112514;
    0.00001249 -0.00046881;
    0.00023032 -0.00191488;
    0.00079945 -0.00329201;
    0.00170287 -0.00470585;
    0.00354717 -0.00688469;
    0.00592084 -0.00912202;
    0.01810144 -0.01720842;
    0.03471169 -0.02488211;
    0.05589286 -0.03226730;
    0.08132751 -0.03908459;
    0.11073805 -0.04503763;
    0.14391397 -0.04986836;
    0.18067874 -0.05338180;
    0.22089879 -0.05551392;
    0.26433734 -0.05636585;
    0.31062190 -0.05605816;
    0.35933893 -0.05472399;
    0.40999990 -0.05254383;
    0.46204424 -0.04969990;
    0.51483073 -0.04637175;
    0.56767889 -0.04264894;
    0.61998250 -0.03859653;
    0.67114514 -0.03433153;
    0.72054815 -0.02996944;
    0.76758733 -0.02560890;
    0.81168064 -0.02134397;
    0.85227225 -0.01726049;
    0.88883823 -0.01343567;
    0.92088961 -0.00993849;
    0.94797259 -0.00679919;
    0.96977487 -0.00402321;
    0.98607009 -0.00180118;
    0.99640466 -0.00044469;
    1.00000000  0.00000000;
]

sections = zeros(3, size(airfoil, 1), length(points))
for ip = 1:length(points)
    chord = root_chord * (1 - x[ip]/L) + tip_chord * x[ip]/L
    sections[1, :, ip] .= 0
    sections[2, :, ip] .= chord .* (airfoil[:,1] .- 0.5)
    sections[3, :, ip] .= chord .* airfoil[:,2]
end

write_vtk("wind-turbine-blade-simulation", assembly, history, t; sections = sections)

# ![](../assets/wind-turbine-blade-simulation.gif)


text:
# Structure representing an instance of a Master side problem.
# Contains the configuration of the execution parameters and other management parameters, the final solving set and the final solution

mutable struct MasterInstance
    # n parameter
    n::Int64
    # k parameter
    k::Int64
    # m parameter
    m::Int64
    # Distance function type, used by the Slave to calculate distance betweeen 2 points
    dist::String
    # Slave dataset cardinality info
    datasetSize::Int64
    # Point dimension
    pointDim::Int64
    # Set containing the Solving Set
    solvingSet::Set{Point}
    # Heap containing the final Solution
    solution::Heap

    # Default costructor
    MasterInstance() = new(0, 0, 0, "", 0, 0, Set{Point}(), MinHeap())

    function MasterInstance(n::Int64, k::Int64, m::Int64, dist::String, dataset_size::Int64, solvingSet::Set{Point}, solution::Heap)
        new(n, k, m, dist, datasetSize, pointDim, Set{Point}(), MinHeap())
    end
end

# Factory method that creates a master instance from a configuration
function createMasterInstance(cfg::MasterConfiguration)
    masterInstance = MasterInstance()
    masterInstance.n = cfg.n
    masterInstance.k = cfg.k
    masterInstance.m = cfg.m
    masterInstance.dist = cfg.dist
    return masterInstance
end

# Factory method that creates a master instance
function createMasterInstance()
    return MasterInstance()
end


text:
function opt_add(self::Tensor{T, N}, other::Tensor{T, N}) where {T, N}
    outputs__ = Int[0]
    __cret = ccall((:atg_add, libtorch_capi),
                 Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
                   outputs__, self.pointer, other.pointer)
    return Tensor{T, N}(Ptr{Cvoid}(outputs__[1]), nothing)
end

function opt_index_select(self::Tensor{T, N}, dim::Int64, index::Int64) where {T, N}
    ptr = ccall((:tensor_method_index_select_int64, libtorch_capi),
                Ptr{Cvoid}, (Ptr{Cvoid}, Clonglong, Clonglong),
                self.pointer, dim, index)
    return Tensor{T, N}(ptr, nothing)
end

function opt_index_select(self::Tensor{T, N}, dim::Int64, i::Array{Int64}) where {T, N}
    return opt_index_select(self, dim, Tensor(i))
end

function opt_index_select(self::Tensor{T, N}, dim::Int64, index::Tensor) where {T, N}
    outputs__ = Int[0]
    __cret = ccall((:atg_index_select, libtorch_capi),
                 Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}, Clonglong, Ptr{Cvoid}),
                 outputs__, self.pointer, dim, index.pointer)
    return Tensor{T, N}(Ptr{Cvoid}(outputs__[1]), nothing)
end

function opt_reshape(self::Tensor{T}, shape_data::Array{Int64}) where T
    outputs__ = Int[0]
    shape_len = length(shape_data)
    __cret = ccall((:atg_reshape, libtorch_capi),
                 Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Cint),
                 outputs__, self.pointer, shape_data, shape_len)
    return Tensor{T, shape_len}(Ptr{Cvoid}(outputs__[1]), nothing)
end


text:
module MRIReco

using Reexport
using FileIO
using ProgressMeter
@reexport using ImageUtils
using NFFT
using Distributions
using LinearOperators
@reexport using RegularizedLeastSquares
using StatsBase
using LinearAlgebra
using Random
@reexport using FFTW
using Distributed
using Graphics: @mustimplement
using Wavelets
using HDF5
using LightXML
using NIfTI
using ThreadPools
@reexport using Unitful
import LowRankApprox.psvd

include("Trajectories/Trajectories.jl")
include("Sequences/Sequence.jl")
include("Datatypes/Datatypes.jl")
include("Tools/Tools.jl")
include("Operators/Operators.jl")
include("Simulation/Simulation.jl")
include("Reconstruction/Reconstruction.jl")
include("IO/IO.jl")
include("Sampling/Sampling.jl")

function __init__()
  if Threads.nthreads() > 1
    BLAS.set_num_threads(1)
    FFTW.set_num_threads(1)
  end
end

end # module


text:
struct Component
    variable::Int64
    poweridx::Int64
    s::Float64
    c::Float64
end

struct Frequency
    coefficients_idx::Array{Int64, 1}
    coefficients::SVector{14, Float64}
    components::Array{Component, 1}
end

function assemble_frequencies(tables)

    frequencies = Array{Frequency, 1}(0)

    for k in 1:length(tables)

        for i in 1:size(tables[k], 1)
            row = tables[k][i, :]

            comp = Component(k, row[1]+1, row[3], row[4])

            j = findfirst(f -> f.coefficients == row[5:18], frequencies)

            if j !== 0
                push!(frequencies[j].components, comp)
            else
                freq = Frequency(find(!iszero, row[5:18]), row[5:18], [comp])
                push!(frequencies, freq)
            end
        end

    end

    return SVector{length(frequencies), Frequency}(frequencies)
end

data_X = sortrows(readdlm(joinpath(@__DIR__, "data/nutation_2000A_X.dlm")), by = r -> hypot(r[3], r[4]))
data_Y = sortrows(readdlm(joinpath(@__DIR__, "data/nutation_2000A_Y.dlm")), by = r -> hypot(r[3], r[4]))
data_s = sortrows(readdlm(joinpath(@__DIR__, "data/nutation_2000A_s.dlm")), by = r -> hypot(r[3], r[4]))

const frequencies = assemble_frequencies([data_X, data_Y, data_s])


text:
using TOML, CSV, DataFrames, ArgParse, JSON

# Used to search through command line for particular arguments
function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--input","-i"
            help = "Specifies an inputfile to replace default parameters"
        "--verbosity", "-v"
            help = "Prints information to command line for debugging"
            action = :store_true
		"--spatial", "-s"
			help = "Use the spatial ploidy model"
			action = :store_true
		"--outputfile","-o"
			help = "Specifies the filename for saving data from the simulation"
			default = "output"
    end

    return parse_args(s)
end

# Convert string keys to symbols (needed for splatting dictionary into struct fields)
string_to_symbol(d::Dict) = Dict(Symbol(k) => v for (k,v) in d)

# Modified get function to get struct from the dict file
function get_structs_from_file(d::Dict,key::String,default::DataType,verbosity::Bool)

	local s
	try
		s = default(string_to_symbol(d[key]))
		if verbosity
			println("Replaced $key with inputfile")
		end
	catch
		s = default()
		if verbosity
			println("Did not replace $key with inputfile")
		end
	end

	return s

end

struct Options

	debugging::Int
	compartmentMinimum::Bool
	replating::Bool
	maxPop::Number
	progress_check::Bool
	saveat::Number
	interpolation_order::Int

end

struct Files
	copy_number_file::String		# Contains numeric array of copy numbers and header of which chromosomes
	birth_rate_file::String			# Contains numeric vector of birth rates
	blood_vessel_file::String		# Blood vessel coordinate info
end

struct UniversalParameters
	stepsize::Number
	minChrom::Number
	maxChrom::Number
	deathRate::Number
	misRate::Number
	finalDay::Number
	startPop::Number
	starting_copy_number::Vector{Int}
	max_cell_cycle_duration::Number
end

struct SpatialParameters
	Γ::Number						# Random cell motion
	Γₑ::Number						# Nutrient diffusion
	ϕ::Number						# Nutrient level needed for half-maximal growth rate
	Ξ::Number						# Energy needed for half maximal directional motion
	χ::Number						# Directed motion magnitude
	δ::Number						# Energy consumption
	Lp::Vector{<:Number}
	Np::Vector{Int}
	k::Number
	E_vessel::Number

end

# Struct containing input data needed for simulation with default values
struct WellMixedInput

	Options
	Files
	UniversalParameters

end

struct SpatialInput

	Options
	Files
	UniversalParameters
	SpatialParameters

end

struct WellMixedInputParameters

	Options
	UniversalParameters

end

struct SpatialInputParameters

	Options
	UniversalParameters
	SpatialParameters
	blood_vessel_file::String

end

function Files(;
	copy_number_file ="birthLandscapeBrainCancer.txt",
	birth_rate_file = "GrowthRate_brain_cancer_CLs.txt",
	blood_vessel_file ="13496_2_Slides_and_Data_xy_test.txt")

	Files(
		copy_number_file,
		birth_rate_file,
		blood_vessel_file
		)
end

Files(d::Dict) = Files(;d...)

function UniversalParameters(;stepsize= 1.0,
	minChrom=1.0,
	maxChrom=5.0,
	deathRate=0.1,
	misRate=0.15,
	finalDay=30.0,
	startPop=1e3,
	starting_copy_number = [1,1,1,1,1],
	max_cell_cycle_duration = 100)

	UniversalParameters(
					stepsize,
					minChrom,
					maxChrom,
					deathRate,
					misRate,
					finalDay,
					startPop,
					starting_copy_number,
					max_cell_cycle_duration
					)

end

UniversalParameters(d::Dict) = UniversalParameters(;d...)

function Options(;debugging = 3,
	compartmentMinimum = false,
	replating = true,
	maxPop = 1e6,
	progress_check = true,
	saveat = 0.1,
	interpolation_order = 2)

	Options(
			debugging,
			compartmentMinimum,
			replating,
			maxPop,
			progress_check,
			saveat,
			interpolation_order,
			)

end

Options(d::Dict) = Options(;d...)

function SpatialParameters(;
	Γ = 0.1,
	Γₑ = 1.0,
	ϕ = 0.02,
	Ξ = 0.02,
	χ = 0.5,
	δ = 0.05,
	Lp = [1000.0,1000.0],
	Np = [21,21],
	k = 0.5,
	E_vessel = 1.0)

	SpatialParameters(
					Γ,
					Γₑ,
					ϕ,
					Ξ,
					χ,
					δ,
					Lp,
					Np,
					k,
					E_vessel
					)

end

SpatialParameters(d::Dict) = SpatialParameters(;d...)

function WellMixedInput(;Options=Options(),
	Files=Files(),
	UniversalParameters=UniversalParameters())

	WellMixedInput(
					Options,
					Files,
					UniversalParameters
					)

end

function WellMixedInput(inputfile::String,verbosity::Bool=false)

	data = TOML.tryparsefile(inputfile)

	WellMixedInput(
					get_structs_from_file(data,"Options",Options,verbosity),
					get_structs_from_file(data,"Files",Files,verbosity),
					get_structs_from_file(data,"UniversalParameters",UniversalParameters,verbosity)
					)

end

function SpatialInput(;Options=Options(),
	Files=Files(),
	UniversalParameters=UniversalParameters(),
	SpatialParameters=SpatialParameters())

	SpatialInput(
					Options,
					Files,
					UniversalParameters,
					SpatialParameters
					)
end

function SpatialInput(inputfile::String,verbosity::Bool=false)

	data = TOML.tryparsefile(inputfile)

	SpatialInput(
					get_structs_from_file(data,"Options",Options,verbosity),
					get_structs_from_file(data,"Files",Files,verbosity),
					get_structs_from_file(data,"UniversalParameters",UniversalParameters,verbosity),
					get_structs_from_file(data,"SpatialParameters",SpatialParameters,verbosity)
					)

end

function initialize()

	# Grab command line args
	parsed_args = parse_commandline()

	# Prints info if set to true
	verbosity = parsed_args["verbosity"]

	spatial = parsed_args["spatial"]

	# Grab output file name
	outputFile = parsed_args["outputfile"]

	# # Grab initial condition (it is a string so we parse and convert)
	# starting_copy_number = Int.(JSON.parse(parsed_args["starting_copy_number"]))

	# Default values to be fed into ploidyMovement
	if spatial
		input = SpatialInput()
		# Replaces default parameters if input file is given
		if !isnothing(parsed_args["input"])
			input = SpatialInput(parsed_args["input"],verbosity)
		end
	else
		input = WellMixedInput()
		# Replaces default parameters if input file is given
		if !isnothing(parsed_args["input"])
			input = WellMixedInput(parsed_args["input"],verbosity)
		end
	end

	# Check to see if input file is given
    if verbosity
		println("parsed_args:")
		for (arg,val) in parsed_args
			println("  $arg => $val")
		end
		println()
	end

	return input, outputFile, verbosity, spatial
end

function extract_XY(X_filename::String, Y_filename::String)

	# Read the birth rate file to be used for the interpolation
	X_df = CSV.read(X_filename,DataFrame)
	Y_df = CSV.read(Y_filename, DataFrame)

	# We make all strings upper case to ensure consistency
	[X_df[!,name] = uppercase.(x) for (name,x) in zip(names(X_df),eachcol(X_df))
					if eltype(x) === String]
	[Y_df[!,name] = uppercase.(y) for (name,y) in zip(names(Y_df),eachcol(Y_df))
					if eltype(y) === String]

	# Join the data frames by cell line name
	XY_df = outerjoin(X_df, Y_df,on = intersect(names(X_df), names(Y_df)))

	# Error if the number of elements in X_df or Y_df do not match
	@assert nrow(X_df) == nrow(Y_df) == nrow(XY_df) "Cell line names likely did not match."

	# Once we are sure the sizes are the same, we remove missings
	dropmissing!(XY_df)

	return XY_df

end

function main()

	# Grab input parameters and whether to print info to terminal
	input, outputFile, verbosity, spatial = initialize()

	# Printing stuff to terminal
	if verbosity
		for name in fieldnames(typeof(input))
			println("$(name):")
			for var in fieldnames(typeof(getfield(input,name)))
				println("  $var => $(getfield(getfield(input,name),var))")
			end
			println()
		end
	end

	# Error checking
	if !isfile(input.Files.copy_number_file)
		error("$(input.Files.copy_number_file) file cannot be found")
	end
	if !isfile(input.Files.birth_rate_file)
		error("$(input.Files.birth_rate_file) file cannot be found")
	end

	# Extract copy number (X) and birth rate (Y) to be used in the interpolation in ploidyMovement.jl
	CN_birthrate_df = extract_XY(input.Files.copy_number_file,input.Files.birth_rate_file)

	# Get copy number and birth rates
	local copy_number,birth_rate,chromosome_header
	try
		copy_number,birth_rate = Matrix(CN_birthrate_df[!,r"chr"]),Matrix(CN_birthrate_df[!,r"birth"])
		chromosome_header = names(CN_birthrate_df[!,r"chr"])
	catch
		@error("chromosome X should be the column names")
	end

	# Drop dims is required to convert to a vector (from mat, required for polyharmonic)
	birth_rate = dropdims(birth_rate,dims=2)

	# Run ploidy movement and grab params from the input
	### FIX ME: We want to replace input file with dictionaries that we can grab the param chunk without deleting

	if spatial

		params = SpatialInputParameters(input.Options,
										input.UniversalParameters,
										input.SpatialParameters,
										input.Files.blood_vessel_file)

		sol, cnArray = runPloidyMovement(params,copy_number,birth_rate,spatial)

		# Write results to multiple files by time points
		for (index,t) in enumerate(sol.t)
			open(outputFile*"_states_time_"*replace(string(t),"." => "_")*".csv","w") do io
				for (i,cnstate) in enumerate(eachrow(cnArray))
					z = sol[index].s[cnstate...,:,:]
					cnstate = "#"*join(string(cnArray[i,:]...),":")
					writedlm(io,[cnstate])
					writedlm(io,z,',')
				end
			end
			open(outputFile*"_Energy_time_"*replace(string(t),"." => "_")*".csv","w") do io
				writedlm(io,sol[index].E,',')
			end
		end

		# Save discretized domain
		writedlm(outputFile*"_x.csv",range(0,1,length=params.SpatialParameters.Np[1]))
		writedlm(outputFile*"_y.csv",range(0,1,length=params.SpatialParameters.Np[2]))

	else

		params = WellMixedInputParameters(input.Options,input.UniversalParameters)

		sol, cnArray,time = runPloidyMovement(params,copy_number,birth_rate,spatial)

		# convert to array for output
		results = hcat(cnArray,sol)
		outputHeader = permutedims(vcat(chromosome_header,time))

		# concatnate the header with the results array
		output = vcat(outputHeader,results)

		# save to file
		writedlm( outputFile*".csv",  output, ',')

	end

end

@time main()

max_stars_repo_name: MathOnco/ploidyEvolution

max_issues_repo_name: MathOnco/ploidyEvolution

max_forks_repo_name: MathOnco/ploidyEvolution


text:
module DriftDiffusion

using ForwardDiff

export adapt_clicks, compute_LL, VGHfun

function compute_LL(bup_times::Array{<:AbstractFloat},  bup_side::AbstractArray,
    stim_dur::Array{<:AbstractFloat}, poke_r::Array{Bool}, input_params::Any ;
    nan_times=Array{Bool}(0), use_param=fill(true,1,9), param_default=[0 1 1 0 1 0.1 0 0.01 1],
    use_prior=zeros(1,9), prior_mu=zeros(1,9), prior_var=zeros(1,9), window_dt=0.01,
    adaptation_scales_perclick="var", nt=[])

    # validate inputs
    @assert size(bup_times,1)==size(bup_side,1) "bup_times and bup_side must have same 1st dimension length (i.e. number of maxbups)"
    @assert size(bup_side,2)==size(stim_dur,2) && size(bup_side,2)==size(poke_r,2) "bup_side, stim_dur, and poke_r must have the same 2nd dimension length (i.e. number of trials)"
    @assert length(use_param)==9 &&  length(use_param)==length(param_default) &&
        length(use_param)==length(use_prior) && length(use_param)==length(prior_mu) &&
        length(use_param)==length(prior_var) "use_param, param_default, use_prior, prior_mu, and prior_var must all have length 9 (because this is currently a 9-parameter model) ";
    if ~all(prior_var[use_prior.==0] .== 0)
        warn("some unused priors have nonzero variance")
    end
    if ~all(use_prior[prior_var.==0] .== 0)
        error("cannot use a prior with 0 variance")
    end

    # initialize variables
    # (Adrian: initializing some variables to be the same type as input_params is crucial, because
    # when you are using compute_LL to compute the Hessian, these must all be duals.)
    params_full              = zeros(eltype(input_params), 1,9)
    params_full[.!use_param] = param_default[.!use_param]
    params_full[use_param]   = input_params
    if params_full[end]<1 & isempty(nt)
        sz=size(bup_times);
        bup_times_normalized = broadcast(/,bup_times,stim_dur);
        nt = Int(round((1-params_full[end])/window_dt))+1;
        bup_times = repmat(bup_times,1,nt);
        in_window = Array{Bool}(size(bup_times));
        for t=1:nt
            inds=(1:sz[2])+sz[2]*(t-1);
            time_window = [0 params_full[end]] + (t-1)*window_dt ;
            in_window[:,inds] = (bup_times_normalized.>time_window[1]) .& (bup_times_normalized.<=time_window[2]);
        end
        bup_times[.~in_window] = NaN;
        nan_times = isnan.(bup_times);
        stim_dur  = stim_dur.*params_full[end];
    elseif isempty(nt)
        nt=1;
    end
    nansum(x::Array{Float64}) = sum(x.*.!isnan.(x),1);
    nTrials = length(poke_r);
    prob_poked_r = zeros(eltype(input_params),nt,nTrials);
    NLL = zeros(eltype(input_params),1,nTrials);
    adapted = ones(eltype(bup_side),size(bup_side))
    init_var = input_params[4];
    a_var = input_params[2];
    c_var = input_params[3];
    lambda = input_params[1];
    bias    = params_full[7]
    lapse   = params_full[8]
    bups_lambda_applied=1;
    if lambda ==0
        s2 = init_var + a_var*stim_dur
    else
        s2 = init_var*exp.(2*lambda*stim_dur) + (a_var./(2*lambda))*(exp.(2*lambda*stim_dur)-1)
    end
    # calculate LL simultaneously for all trials, looping over time points
    for t=1:nt
        inds=(1:nTrials)+nTrials*(t-1);
        curr_buptimes = bup_times[:,inds];
       if isempty(nan_times)
           curr_nantimes = isnan.(curr_buptimes);
       else
           curr_nantimes = nan_times[:,inds];
       end
        # adapt those clicks
        adapted = adapt_clicks(curr_buptimes, curr_nantimes, params_full[5],params_full[6])
        # apply integration timescale
        if lambda!=0
            bups_lambda_applied = exp(lambda * time_from_end[:,inds])
        end
        # compute mean of distribution
        mean_a = nansum(adapted.*bups_lambda_applied.*bup_side,1);
        # add per click variance
        if adaptation_scales_perclick=="std"
            var_a      = s2 + nansum(c_var .* adapted.^2 .* bups_lambda_applied .^ 2,1);
        elseif adaptation_scales_perclick=="var"
            var_a      = s2 + nansum(c_var .* adapted .* bups_lambda_applied .^ 2 ,1);
        elseif adaptation_scales_perclick=="none"
            var_a      = s2 + nansum(c_var .* bups_lambda_applied .^ 2,1);
        end
        no_noise_trials = (var_a .==0);
        if any(no_noise_trials)
            warn( @sprintf("Total noise variance is zero for %g trials. ",sum(no_noise_trials)),
                "Could be due to a model with no noise terms or only stimulus noise but no bups on a trial. ",
                "Adding eps to avoid NaNs but you may want to do some sanity checks.");
            var_a[no_noise_trials]=eps;
        end
        erfTerm[:,inds] =  (mean_a-bias)./sqrt.(2*var_a);
    end
    erfTerm = erf(erfTerm);
    if any(abs(erfTerm).==1)
        erfTerm[erfTerm.==1]=1-eps();
        erfTerm[erfTerm.==-1]=eps()-1;
    end
    prob_poked_r=((1-lapse).*(1+erfTerm)+lapse)/2 ;
    prob_poked_r = mean(prob_poked_r,1);
    NLL[poke_r] = - ( log.( prob_poked_r[poke_r] ) );
    NLL[.!poke_r] = - ( log.(1- prob_poked_r[.!poke_r] ) );
    NLL_total = sum(NLL)
    # increment likelihood for priors
    for pp = find(use_prior)
        NLL_total += -(params_full[pp]-prior_mu[pp])^2/(2*prior_var[pp])
    end
    return NLL_total
end

function adapt_clicks(bup_times, nan_times, phi,tau_phi) #phi, tau_phi)
    adapt   = zeros(eltype(phi),size(bup_times)); # this must be the same type as the params because these could be duals
    if phi<0 | tau_phi<0
        error("adaptation parameters must be non-negative")
    end
    adapt[.!nan_times] = 1
    ici     = diff(bup_times)
    for i = 2:size(bup_times,1)
        adapt[i, :] = 1 + exp.(-ici[i-1,:]./tau_phi).*(adapt[i-1,:]*phi-1)
        adapt[i, ici[i-1, :] .<= 0] = 0
        adapt[i-1, ici[i-1,:] .<= 0] = 0
    end
    return adapt
end

function VGHfun(bup_times,  bup_side, stim_dur, poke_r, input_params; nan_times=fill(Bool,0), use_param=fill(true,1,9),
    param_default=[0 1 1 0 1 0.1 0 0.01 1],use_prior=zeros(1,9), prior_mu=zeros(1,9), prior_var=zeros(1,9),
    window_dt=0.01, adaptation_scales_perclick="var", nt=[])
    # compute hessian using autodiff
    # (Adrian: For reasons I don't understand, you need to call compute_LL with explicit keyword declaration for ForwardDiff to run
    optimFun(x) = DriftDiffusion.compute_LL(bup_times, bup_side, stim_dur, poke_r, x;nan_times=nan_times, use_param=use_param,
        param_default=param_default,use_prior=use_prior, prior_mu=prior_mu, prior_var=prior_var,
        window_dt=window_dt, adaptation_scales_perclick=adaptation_scales_perclick, nt=nt);
    out = DiffResults.HessianResult(zeros(size(input_params)));
    out = ForwardDiff.hessian!(out,optimFun,input_params);
    return DiffResults.value(out),DiffResults.gradient(out),DiffResults.hessian(out)
end

end # module


text:
using CuNFFT
import CuNFFT.CuArray

m = 5
σ = 2.0

# test CuNFFT
if CuNFFT.CUDA.functional()
  @testset "CuNFFT in multiple dimensions" begin
    for (u,N) in enumerate([(256,), (32,32), (12,12,12)])
        eps = [1e-7, 1e-3, 1e-6, 1e-4]
        for (l,window) in enumerate([:kaiser_bessel, :gauss, :kaiser_bessel_rev, :spline])
            D = length(N)
            @info "Testing CuNFFT in $D dimensions using $window window"

            M = prod(N)
            x = rand(Float64,D,M) .- 0.5
            p = plan_nfft(Array, x, N; m, σ, window, precompute = NFFT.FULL,
                         fftflags = FFTW.ESTIMATE)
            p_d = plan_nfft(CuArray, x, N; m, σ, window, precompute = NFFT.FULL)
            pNDFT = NDFTPlan(x, N)

            fHat = rand(Float64,M) + rand(Float64,M)*im
            f = adjoint(pNDFT) * fHat
            fHat_d = CuArray(fHat)
            fApprox_d = adjoint(p_d) * fHat_d
            fApprox = Array(fApprox_d)
            e = norm(f[:] - fApprox[:]) / norm(f[:])
            @debug "error adjoint nfft "  e
            @test e < eps[l]

            gHat = pNDFT * f
            gHatApprox = Array( p_d * CuArray(f) )
            e = norm(gHat[:] - gHatApprox[:]) / norm(gHat[:])
            @debug "error nfft "  e
            @test e < eps[l]
        end
    end
  end
end


text:
using NamedDims
using NamedDims: names
using SparseArrays
using Test

@testset "get the parent array that was wrapped" begin
    orig = [1 2; 3 4]
    @test parent(NamedDimsArray(orig, (:x, :y))) === orig
end

@testset "get the named array that was wrapped" begin
    @test names(NamedDimsArray([10 20; 30 40], (:x, :y))) === (:x, :y)
end

@testset "Name-asserting constructor" begin
    orig_full = NamedDimsArray(ones(3,4,5), (:a, :b, :c))
    @test names(NamedDimsArray(orig_full, (:a, :b, :c))) == (:a, :b, :c)
    @test names(NamedDimsArray(orig_full, (:a, :b, :_))) == (:a, :b, :c)
    @test_throws DimensionMismatch NamedDimsArray(orig_full, (:a, :b, :wrong))
    @test_throws DimensionMismatch NamedDimsArray(orig_full, (:c, :a, :b))

    orig_partial = NamedDimsArray(ones(3,4,5), (:a, :_, :c))
    @test names(NamedDimsArray(orig_partial, (:a, :b, :c))) == (:a, :b, :c)
    @test names(NamedDimsArray(orig_partial, (:a, :_, :c))) == (:a, :_, :c)
end

@testset "getindex" begin
    nda = NamedDimsArray([10 20; 30 40], (:x, :y))

    @test nda[x=1, y=1] == nda[y=1, x=1] == nda[1, 1] == nda[CartesianIndex(1,1)] == 10
    @test nda[y=end,x=end] == nda[end, end] == 40

    # Unspecified dims become slices
    @test nda[y=1] == nda[y=1, x=:] == nda[:, 1] == [10; 30]

    @testset "name preservation" begin
        @test names(nda[y=1]) == (:x,)
        @test names(nda[y=1:1]) == (:x, :y)
    end

    # https://github.com/invenia/NamedDims.jl/issues/8
    @testset "with multiple-wildcards" begin
        nda_mw = NamedDimsArray{(:_,:_,:c)}(ones(10,20,30));
        @test nda_mw[c=2] == ones(10,20)
        @test names(nda_mw[c=2]) == (:_, :_)
    end
end

@testset "views" begin
    nda = NamedDimsArray([10 20; 30 40], (:x, :y))

    @test @view(nda[y=1]) == @view(nda[y=1, x=:]) == @view(nda[:, 1]) == [10; 30]

    @testset "name preservation" begin
        @test names(nda[y=1]) == (:x,)
        @test names(nda[y=1:1]) == (:x, :y)
    end
end

@testset "setindex!" begin
    @testset "by name" begin
        nda = NamedDimsArray([10 20; 30 40], (:x, :y))

        nda[x=1, y=1] = 100
        @test nda == [100 20; 30 40]

        nda[x=1] .= 1000
        @test nda == [1000 1000; 30 40]
    end

    @testset "by position" begin
        nda = NamedDimsArray([10 20; 30 40], (:x, :y))

        nda[1, 1] = 100
        @test nda == [100 20; 30 40]

        nda[1,:] .= 1000
        @test nda == [1000 1000; 30 40]
    end
end

@testset "IndexStyle" begin
    nda = NamedDimsArray([10 20; 30 40], (:x, :y))
    @test IndexStyle(typeof(nda)) == IndexLinear()

    sparse_nda = NamedDimsArray(spzeros(4,2), (:x, :y))
    @test IndexStyle(typeof(sparse_nda)) == IndexCartesian()
end

@testset "length/size/axes" begin
    nda = NamedDimsArray([10 20; 30 40; 50 60], (:x, :y))

    @test length(nda) == 6

    @test axes(nda) == (1:3, 1:2)
    @test axes(nda,:x) == (1:3) == axes(nda, 1)

    @test size(nda) == (3,2)
    @test size(nda,:x) == 3 == size(nda, 1)
end

@testset "similar" begin
    nda = NamedDimsArray{(:a, :b, :c, :d)}(ones(10, 20, 30, 40))

    @testset "content" begin
        ndb = similar(nda)
        @test parent(ndb) !== parent(nda)
        @test eltype(ndb) == Float64
        @test size(ndb) == (10, 20 ,30, 40)
        @test names(ndb) == (:a, :b, :c, :d)
    end
    @testset "eltype" begin
        ndb = similar(nda, Char)
        @test parent(ndb) !== parent(nda)
        @test eltype(ndb) == Char
        @test size(ndb) == (10, 20, 30, 40)
        @test names(ndb) == (:a, :b, :c, :d)
    end
    @testset "size" begin
        ndb = similar(nda, Float64, (15, 25, 35, 45))
        @test parent(ndb) !== parent(nda)
        @test eltype(ndb) == Float64
        @test size(ndb) == (15, 25, 35, 45)
        @test names(ndb) == (:a, :b, :c, :d)
    end
    @testset "dim names" begin
        ndb = similar(nda, Float64, (:w, :x, :y, :z))
        @test parent(ndb) !== parent(nda)
        @test eltype(ndb) == Float64
        @test size(ndb) == (10, 20, 30, 40)
        @test names(ndb) == (:w, :x, :y, :z)
    end

    @testset "dimensions" begin
        ndb = similar(nda, Float64, (w=11, x=22))
        @test parent(ndb) !== parent(nda)
        @test eltype(ndb) == Float64
        @test size(ndb) == (11,22)
        @test names(ndb) == (:w, :x,)
    end
end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
using .DataFiles

# =========================================================================================
"""
    update_languages(
        languages::AbstractString = "eng";
        target::AbstractString = "tessdata",
        frequency::Integer = 7,
        source::GitHubProject = DATA_REPO
    )::Bool

Update the training data files for the specified languages.  Returns false if there was an
error updating the files.

__Arguments:__

| T | Name      | Default    | Description
|:--| :-------- | :----------| :----------
| O | languages | `eng`      | The languages to update separated with "+".
| K | target    | `tessdata` | The directory to save the data files in.
| K | frequency | `7`        | How often to check for updates on the server in days.
| K | source    | ...        | The GitHub project to download the training data from.

__Details:__

By default the data files are saved in a "tessdata" directory under the current directory.
The training data is normally downloaded from the
https://github.com/tesseract-ocr/tessdata_best GitHub project.

See also: [`download_languages`](@ref)
"""
function update_languages(
            languages::AbstractString = "eng";
            target::AbstractString = "tessdata",
            frequency::Integer = 7,
            source::GitHubProject = DATA_REPO
        )::Bool
    local result = true
    local index = load_languages(target)

    try
        local list  = split(languages, "+")
        local force = !has_all_languages(index, list)

        refresh_languages(index, source, frequency, force)

        for lang in list
            if has_language(index, lang) == true
                local file = get_language_file(index, lang)

                if file_out_of_date(file, target)
                    local fullpath = joinpath(target, file.filename)
                    file.localsha1 = download_file(file.url, fullpath)
                    file.timestamp = mtime(fullpath)
                end
            else
                @warn "Could not find training data for language: $(lang)"
                result = false
            end

        end

    catch ex
        if isa(ex, FileError) == false && isa(ex, NetworkError) == false
            rethrow(ex)
        end
        result = false
    end

    save_languages(index, target)

    return result
end

# =========================================================================================
"""
    update_pdf_font(;
        target::AbstractString = "tessdata",
        frequency::Integer = 7,
        source::GitHubProject = PDF_FONT_REPO
    )::Bool

Update the PDF font file.  This is only needed when generating PDF files.  Returns false if
there was an error updating the files.

__Arguments:__

| T | Name      | Default    | Description
|:--| :-------- | :----------| :----------
| K | target    | `tessdata` | The directory to save the font file in.
| K | frequency | `7`        | How often to check for updates on the server in days.
| K | source    | ...        | The GitHub project to download the font file from.

__Details:__

By default the data files are saved in a "tessdata" directory under the current directory.
The training data is normally downloaded from the
https://github.com/tesseract-ocr/tessdata_best GitHub project.

See also: [`download_pdf_font`](@ref)
"""
function update_pdf_font(;
            target::AbstractString = "tessdata",
            frequency::Integer = 7,
            source::GitHubProject = PDF_FONT_REPO
        )::Bool
    local result = true
    local info = load_pdf_font_info(target)

    try
        refresh_pdf_font_info(info, source, frequency)

        if file_out_of_date(info.font, target)
            local fullpath = joinpath(target, info.font.filename)
            info.font.localsha1 = download_file(info.font.url, fullpath)
            info.font.timestamp = mtime(fullpath)
        end

    catch ex
        if isa(ex, FileError) == false && isa(ex, NetworkError) == false
            rethrow(ex)
        end
        result = false
    end

    save_pdf_font_info(info, target)

    return result
end


text:
using Flux, CuArrays, Base.Test

info("Testing Flux/CUDNN")

@testset "RNN" begin
  @testset for R in [RNN, GRU, LSTM]
    rnn = R(10, 5)
    curnn = mapleaves(gpu, rnn)
    @testset for batch_size in (1, 5)
      Flux.reset!(rnn)
      Flux.reset!(curnn)
      x = batch_size == 1 ?
        param(rand(10)) :
        param(rand(10,batch_size))
      cux = gpu(x)
      y = (rnn(x); rnn(x))
      cuy = (curnn(cux); curnn(cux))

      @test y.data ≈ collect(cuy.data)
      @test haskey(Flux.CUDA.descs, curnn.cell)

      Δ = randn(size(y))

      Flux.back!(y, Δ)
      Flux.back!(cuy, gpu(Δ))

      @test x.grad ≈ collect(cux.grad)
      @test rnn.cell.Wi.grad ≈ collect(curnn.cell.Wi.grad)
      @test rnn.cell.Wh.grad ≈ collect(curnn.cell.Wh.grad)
      @test rnn.cell.b.grad ≈ collect(curnn.cell.b.grad)
      @test rnn.cell.h.grad ≈ collect(curnn.cell.h.grad)
      if isdefined(rnn.cell, :c)
        @test rnn.cell.c.grad ≈ collect(curnn.cell.c.grad)
      end

      Flux.reset!(rnn)
      Flux.reset!(curnn)
      ohx = batch_size == 1 ?
        Flux.onehot(rand(1:10), 1:10) :
        Flux.onehotbatch(rand(1:10, batch_size), 1:10)
      cuohx = gpu(ohx)
      y = (rnn(ohx); rnn(ohx))
      cuy = (curnn(cuohx); curnn(cuohx))

      @test y.data ≈ collect(cuy.data)
    end
  end
end


text:
mutable struct Laser <: GameObject
    pos::Vector{<:Real}
    vel::Vector{<:Real}
    @init function Laser(x::Real, y::Real, ϕ::Real)
        this.pos = [x, y]
        this.vel = 10 * [cos(ϕ), sin(ϕ)]
    end
end

function update!(this::Laser)
    this.pos += this.vel
end

function render(this::Laser)
    @layer begin
        sethue("white")
        circle(this.pos..., 4, :fill)
    end
end

function -->(this::Laser, asteroid::Asteroid)
    d = distance(Point(this.pos...), Point(asteroid.pos...))
    return d < asteroid.r
end

max_stars_repo_name: jorge-brito/LynxGallery

max_issues_repo_name: jorge-brito/LynxGallery

max_forks_repo_name: jorge-brito/LynxGallery


text:
"""
    struct PackedHessenberg{T,V<:AbstractVector{T}} <: AbstractMatrix{T}
        data::V
        n::Int
    end

A custom struct to store a Hessenberg matrix in a packed format (without zeros). Hereto, the
non-zero entries are stored sequentially in vector `data` of length `n(n+1)/2`.
"""
struct PackedHessenberg{T,V<:AbstractVector{T}} <: AbstractMatrix{T}
    data::V
    n::Int
    function PackedHessenberg{T,V}(data::V, n::Int) where {T,V<:AbstractVector{T}}
        @assert length(data) >= ((n * n + 3 * n - 2) >> 1)
        return new{T,V}(data, n)
    end
end
PackedHessenberg(data::AbstractVector, n::Int) =
    PackedHessenberg{eltype(data),typeof(data)}(data, n)
Base.size(A::PackedHessenberg) = (A.n, A.n)

function Base.replace_in_print_matrix(
    A::PackedHessenberg,
    i::Integer,
    j::Integer,
    s::AbstractString
)
    return i <= j + 1 ? s : Base.replace_with_centered_mark(s)
end

function Base.getindex(A::PackedHessenberg{T}, i::Integer, j::Integer) where {T}
    @boundscheck checkbounds(A, i, j)
    if i > j + 1
        return zero(T)
    else
        return A.data[((j*j+j-2)>>1)+i]
    end
end
function Base.setindex!(A::PackedHessenberg{T}, v, i::Integer, j::Integer) where {T}
    @boundscheck checkbounds(A, i, j)
    if i > j + 1 && !iszero(v)
        throw(ReadOnlyMemoryError())
    else
        A.data[((j*j+j-2)>>1)+i] = v
    end
    return v
end

Base.IndexStyle(::Type{<:PackedHessenberg}) = Base.IndexCartesian()

# TODO: add more methods from the AbstractArray interface


text:
const Py_NotImplemented__ref = Ref(PyNULL)
Py_NotImplemented() = pyglobal(Py_NotImplemented__ref, :_Py_NotImplementedStruct)

PyNotImplemented_Check(o) = Py_Is(o, Py_NotImplemented())

PyNotImplemented_New() = (o = Py_NotImplemented(); Py_IncRef(o); o)

const Py_None__ref = Ref(PyNULL)
Py_None() = pyglobal(Py_None__ref, :_Py_NoneStruct)

PyNone_Check(o) = Py_Is(o, Py_None())

PyNone_New() = (o = Py_None(); Py_IncRef(o); o)

PyNone_TryConvertRule_nothing(o, ::Type{Nothing}) = putresult(nothing)
PyNone_TryConvertRule_missing(o, ::Type{Missing}) = putresult(missing)


text:
@testset "PoissonDDP" begin
    N, G, K0 = 10, 4, 1
    rng = MersenneTwister(1)
    m = PoissonDDP(rng, N, G; K0)
    @test m.G  == 4
    @test m.a0 == 2
    @test m.b0 == 4
    @test m.a1 == [2 * ones(G)]
    @test m.b1 == [4 * ones(G)]
    @test m.gamma  == ones(Bool, G)
end

@testset "PoissonDDP inherited accessors" begin
    N, G, K0 = 10, 4, 1
    rng = MersenneTwister(1)
    m = PoissonDDP(rng, N, G; K0)
    @test dp_mass(m) > 0.0
    @test dp_mass(m) < Inf
    @test n_clusters(m) == K0
    @test sum(cluster_sizes(m) .== 0) == 1
    @test sum(cluster_sizes(m)) == N
    @test active_clusters(m) == Set(1:K0)
    @test passive_clusters(m) == Set(K0 + 1)
    @test cluster_capacity(m) == K0 + 1
end

@testset "add_cluster!" begin
    N, G, K0 = 10, 4, 1
    rng = MersenneTwister(1)
    m = PoissonDDP(rng, N, G; K0)
    ANOVADDPTest.add_cluster!(m)
    @test length(m.a1) == 2
    @test length(m.b1) == 2
    @test length(m.sumlogu) == 2
end

@testset "update_suffstats! (1)" begin
    N, G, K0 = 1, 1, 1
    rng = MersenneTwister(1)
    data = PoissonData([1], [1])
    m = PoissonDDP(rng, N, G; K0)
    ANOVADDPTest.update_suffstats!(m, data)
    # TODO: Add tests
end

@testset "update_suffstats! (2)" begin
    N, G, K0 = 1, 1, 1
    rng = MersenneTwister(1)
    data = PoissonData([1], [1])
    m = PoissonDDP(rng, N, G; K0)
    ANOVADDPTest.update_suffstats!(m, data)
    ANOVADDPTest.update_suffstats!(m, data, 1, 1, 2)
    # TODO: Add tests
end

@testset "logpredlik (empty clusters)" begin
    N, G, K0 = 1, 1, 1
    rng = MersenneTwister(1)
    data = PoissonData([1], [1])
    m = PoissonDDP(rng, N, G; K0)
    ANOVADDPTest.update_suffstats!(m, data)
    ANOVADDPTest.logpredlik(m, data, 1, first(passive_clusters(m)))
    # TODO: Add tests
end

@testset "logpredlik (non-empty clusters)" begin
    N, G, K0 = 2, 1, 1
    rng = MersenneTwister(1)
    data = PoissonData([1, 1], [1, 0])
    m = PoissonDDP(rng, N, G; K0)
    ANOVADDPTest.update_suffstats!(m, data)
    ANOVADDPTest.logpredlik(m, data, 2, 1)
    # TODO: Add some tests
end

@testset "update! (1)" begin
    N, G, K0 = 2, 1, 1
    rng = MersenneTwister(1)
    data = PoissonData([1, 1], [1, 0])
    m = PoissonDDP(rng, N, G; K0)
    update!(rng, m, data)
end

@testset "update! (2)" begin
    rng = MersenneTwister(1)
    N, G, K0 = 1000, 3, 1
    x = rand(rng, 1:G, N)
    y = rand(rng, Poisson(2), N)
    data = PoissonData(x, y)
    m = PoissonDDP(rng, N, G; K0)
    update!(rng, m, data)
end

@testset "train (1)" begin
    rng = MersenneTwister(1)
    N, G, K0 = 1000, 3, 1
    x = rand(rng, 1:G, N)
    y = rand(rng, Poisson(2), N)
    data = PoissonData(x, y)
    predict =
        expandgrid(1:G, 1:8) |>
        x -> PoissonData(x...)
    m = PoissonDDP(rng, N, G; K0)
    chain = train(rng, m, data, predict)
    @test all(mode(chain.gamma) .== [true, false, false])
end

@testset "train (2)" begin
    rng = MersenneTwister(1)
    N, G, K0 = 1000, 3, 1
    x = rand(rng, 1:G, N)
    y = rand(rng, Poisson(2), N)
    for i = 1:N
        if x[i] == 2
            y[i] += 2
        end
    end
    data = PoissonData(x, y)
    predict =
        expandgrid(1:G, 1:8) |>
        x -> PoissonData(x...)
    m = PoissonDDP(rng, N, G; K0)
    chain = train(rng, m, data, predict)
    @test all(mode(chain.gamma) .== [true, true, false])
end

@testset "train (3)" begin
    rng = MersenneTwister(1)
    N, G, K0 = 1000, 3, 1
    x = rand(rng, 1:G, N)
    y = rand(rng, Poisson(2), N)
    for i = 1:N
        if x[i] == 3
            y[i] += 2
        end
    end
    data = PoissonData(x, y)
    predict =
        expandgrid(1:G, 1:8) |>
        x -> PoissonData(x...)
    m = PoissonDDP(rng, N, G; K0)
    chain = train(rng, m, data, predict)
    @test all(mode(chain.gamma) .== [true, false, true])
end

@testset "train (4)" begin
    rng = MersenneTwister(1)
    N, G, K0 = 1000, 3, 1
    x = rand(rng, 1:G, N)
    y = rand(rng, Poisson(3), N)
    for i = 1:N
        if x[i] != 1
            y[i] += 2
        end
    end
    data = PoissonData(x, y)
    predict =
        expandgrid(1:G, 1:8) |>
        x -> PoissonData(x...)
    m = PoissonDDP(rng, N, G; K0)
    chain = train(rng, m, data, predict)
    @test all(mode(chain.gamma) .== [true, true, true])
end


text:
@testset "test_rhs_compose" begin
#test rhs_compose.jl

TF=Float64
TI=Int64

#with BLAS serial
rhs=randn(100000)
rhs_2 = deepcopy(rhs)
rhs_3 = deepcopy(rhs)
p=2
i=10

y=Vector{Vector{Float64}}(p)
y[1]=randn(51000)
y[2]=randn(100000)

l=Vector{Vector{Float64}}(p)
l[1]=randn(51000)
l[2]=randn(100000)

rho=[1.234;10.23432]

TD_OP = Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joLinearFunction{TF,TF}}}(p)
TD_OP[1] = speye(51000,100000)*2.0
TD_OP[2] = speye(100000)

parallel=false
Blas_active = false
rhs_compose!(rhs,l,y,rho,TD_OP,p,Blas_active,parallel)

#with explicit BLAS callse, serial Julia
Blas_active = true
rhs_compose!(rhs_2,l,y,rho,TD_OP,p,Blas_active,parallel)

@test isapprox(rhs,rhs_2,rtol=10*eps(TF))

#test parallel version
parallel=true

#distributed quantities
TD_OP_d  = distribute(TD_OP)
rho_d    = distribute(rho)
y_d      = distribute(y)
l_d      = distribute(l)

rhs_3=rhs_compose!(rhs_3,l_d,y_d,rho_d,TD_OP_d,p,Blas_active,parallel)

@test isapprox(rhs,rhs_3,rtol=10*eps(TF))

end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
function place_sell_order_command(k, v, id_str, config, open_orders, closed_orders, bot_path, client, usr, fees)
    t = try Pinguu.examine_order(closed_orders[id_str], id_str, client, config)
    catch
        return nothing
    end
    qty = floor(t["use_vol"], digits=config["precision"])
    price = round((abs(t["current_profit"]) / qty *
        (1 + fees["maker"] +
            parse(Float64, config["take_profit"])/100)
            * t["price_correction"]), digits=config["sell_precision"])
    price = Pinguu.price2string(price, config["sell_precision"])

    account = client.get_account()["balances"];
    new = 0
    for bal in account
        if bal["asset"] == config["farm_coin"]
            if parse(Float64, bal["free"]) <= qty
                new = parse(Float64, bal["free"])
            end
        end
    end
    if qty > new && new != 0
        reference = qty * parse(Float64, price)
        qty = floor(new, digits=config["precision"])
        price = Pinguu.price2string(reference / qty, config["sell_precision"])
    end
    @info("Try to place sell")
    @info("Qty  = $qty")
    @info("Price = " * price)

    order = try client.order_limit_sell(
        symbol=config["symbol"],
        quantity=qty,
        price=price);
        catch
        nothing
    end
    if order !== nothing
        order["time"] = string(now())
        open_orders["SELL"][k] = order
        @info("Success")
    else
        v["fail"] = 1
        @info("Failed.. Placing another trigger")
    end
    open(joinpath(bot_path, "open_orders.json"), "w") do f
        JSON.print(f, open_orders, 4)
    end
    open(joinpath(bot_path, "trigger.json"), "w") do f
        JSON.print(f, Dict(), 4)
    end
end
export place_sell_order_command

function sell_filled(tmp, k, v, id_str, config, open_orders, closed_orders, bot_path, client, usr, fees)
    #
    ## Check if id is in closed_orders
    !(id_str in keys(closed_orders)) ? closed_orders[id_str] = Dict() : ""
    tmp["time"] = string(now())
    #
    ## Delete Sell from open_orders
    delete!(open_orders["SELL"], k)
    #
    ## Cancel all related open BUYs
    for (k2, v2) in open_orders["BUY"]
        if split(k2, "_")[end] == id_str
            sleep(0.25)
            #
            ## Load order status
            tmp_buy = try client.get_order(
                symbol=config["symbol"],
                orderId=v2["orderId"])
            catch
                nothing
            end
            #
            ## Update open/closed orders
            if tmp_buy !== nothing
                if !(tmp_buy["status"] in [client.ORDER_STATUS_FILLED, client.ORDER_STATUS_CANCELED])
                    sleep(0.25)
                    result = try client.cancel_order(
                        symbol=config["symbol"],
                        orderId=v2["orderId"])
                    catch
                        nothing
                    end
                    if result !== nothing
                        delete!(open_orders["BUY"], k2)
                    else
                        open_orders["BUY"][k2]["status"] = "cancel"
                    end
                else
                    delete!(open_orders["BUY"], k2)
                end
                if parse(Float64, tmp_buy["executedQty"]) > 0
                    closed_orders[id_str][k2] = tmp_buy
                end
            else
                open_orders["BUY"][k2]["status"] = "cancel"
            end
        end
    end
    #
    ## Write to files
    open(joinpath(bot_path, "closed_orders.json"), "w") do f
        JSON.print(f, closed_orders, 4)
    end
    open(joinpath(bot_path, "open_orders.json"), "w") do f
        JSON.print(f, open_orders, 4)
    end
    qty = 0     # Totel volume
    net_qty = 0 # Total worth
    val = 0     # Value in base coin
    for (k2,v2) in closed_orders[id_str]
        v2["side"] == "BUY" ? qty += parse(Float64, v2["executedQty"]) : ""
        v2["side"] == "BUY" ? net_qty += parse(Float64, v2["executedQty"]) : net_qty -= parse(Float64, v2["executedQty"])
        v2["side"] == "BUY" ? val -= parse(Float64, v2["cummulativeQuoteQty"]) * (1.0 + fees["taker"]) : val += parse(Float64, v2["cummulativeQuoteQty"]) * (1.0 - fees["maker"])
    end
    #
    ## Check if this was the only SELL
    if val < 0 && val + parse(Float64, tmp["cummulativeQuoteQty"]) < 0
        #
        ## If there was an extra SELL executed
        rnd_id = "SELL-" * lpad(rand(1:1:1000), 6, "0")
        while rnd_id in keys(closed_orders[id_str])
            rnd_id = "SELL-" * lpad(rand(1:1:1000), 6, "0")
        end
        closed_orders[id_str][rnd_id] = deepcopy(tmp);
        open(joinpath(bot_path, "closed_orders.json"), "w") do f
            JSON.print(f, closed_orders, 4)
        end
        #
        ## Create new sell
        if abs(val) > config["min_qty"]
            #
            ## qty is high enough to be sold
            new_qty = floor(net_qty, digits=config["precision"])
            new_price = Pinguu.price2string(abs(val)/net_qty * (1 + parse(Float64, config["take_profit"])/100), config["sell_precision"], round_up=true)
            sleep(0.2)
            order = try client.order_limit_sell(
                symbol=config["symbol"],
                quantity=new_qty,
                price=new_price);
                catch
                nothing
            end
            new_id = "SELL_" * id_str
            if order !== nothing
                order["time"] = string(now())
                open_orders["SELL"][new_id] = order
                open(joinpath(bot_path, "open_orders.json"), "w") do f
                    JSON.print(f, open_orders, 4)
                end
                #
                ## Place new buy
                extra_iteration = 0
                for (k2,v2) in closed_orders[id_str]
                    if length(split(k2, "SELL")) < 2
                        if parse(Int64, split(split(k2, "_")[1], "-")[end]) > extra_iteration
                            extra_iteration = parse(Int64, split(split(k2, "_")[1], "-")[end])
                        end
                    end
                end
                extra_iteration += 1
                if extra_iteration <= parse(Float64, config["extra_order_count"])
                    new_price = Pinguu.price2string(parse(Float64, new_price) * (1 - config["extra_order_step"]), config["sell_precision"])
                    quant = floor(config["order_volume"] * (config["extra_order_martingale"])^(extra_iteration + 1) / parse(Float64, new_price), digits=config["precision"])

                    new_buy = Dict()
                    new_buy["status"] = "place order"
                    new_buy["origQty"] = quant
                    new_buy["price"] = new_price
                    new_buy_id = "EXTRA-" * string(extra_iteration) * "_" * id_str
                    open_orders["BUY"][new_buy_id] = new_buy
                    open(joinpath(bot_path, "open_orders.json"), "w") do f
                        JSON.print(f, open_orders, 4)
                    end
                end
            else
                #
                ## If new sell could not be placed
                open_orders["SELL"][new_id] = Dict()
                open_orders["SELL"][new_id]["origQty"] = new_quant
                open_orders["SELL"][new_id]["price"]   = new_price
                open_orders["SELL"][new_id]["status"]  = "place order"
                open_orders["SELL"][new_id]["time"]    = string(now())
                open_orders["SELL"][new_id]["fail"]    = 0
                open_orders["SELL"][new_id]["if_fail"] = new_quant - 1/10^config["precision"]

                open(joinpath(bot_path, "open_orders.json"), "w") do f
                    JSON.print(f, open_orders, 4)
                end
            end
        else
            #
            ## qty is NOT high enough to be sold
            new_price = Pinguu.price2string(abs(val)/net_qty * (1 + parse(Float64, config["take_profit"])/100), digits=config["sell_precision"])

            extra_iteration = 0
            for (k2,v2) in closed_orders[id_str]
                if length(split(k2, "SELL")) < 2
                    if parse(Int64, split(split(k2, "_")[1], "-")[end]) > extra_iteration
                        extra_iteration = parse(Int64, split(split(k2, "_")[1], "-")[end])
                    end
                end
            end
            extra_iteration += 1
            if extra_iteration <= parse(Float64, config["extra_order_count"])
                new_price = string(round(parse(Float64, new_price) * (1 - config["extra_order_step"]), digits=config["sell_precision"]))
                quant = floor(config["order_volume"] * (config["extra_order_martingale"])^(extra_iteration + 1) / parse(Float64, new_price), digits=config["precision"])

                new_buy = Dict()
                new_buy["status"] = "place order"
                new_buy["origQty"] = quant
                new_buy["price"] = new_price
                new_buy_id = "EXTRA-" * string(extra_iteration) * "_" * id_str
                open_orders["BUY"][new_buy_id] = new_buy
                open(joinpath(bot_path, "open_orders.json"), "w") do f
                    JSON.print(f, open_orders, 4)
                end
            end
        end
    else
        #
        ## If everything went well and the profit is positive
        sell_id = !haskey(closed_orders[id_str], "SELL") ? "SELL" : "SELL-" * lpad(rand(1:1:100000), 6, "0")
        closed_orders[id_str][sell_id] = deepcopy(tmp);

        open(joinpath(bot_path, "closed_order-" * id_str * ".json"), "w") do f
            JSON.print(f, closed_orders[id_str], 4)
        end

        try for (k5,v5) in closed_orders
                if DateTime(split(k5, "--")[end]) < now() - Day(42)
                    delete!(closed_orders, k5)
                end
            end
        catch;
        end

        open(joinpath(bot_path, "closed_orders.json"), "w") do f
            JSON.print(f, closed_orders, 4)
        end
    end
end
export sell_filled

function sell_canceled(tmp, k, v, id_str, config, open_orders, closed_orders, bot_path, client, usr, fees)
    delete!(open_orders["SELL"], k)
    open(joinpath(bot_path, "open_orders.json"), "w") do f
        JSON.print(f, open_orders, 4)
    end
    if parse(Float64, tmp["executedQty"]) > 0
        closed_orders[id_str]["SELL-" * lpad(rand(1:1:100000), 6, "0")] = tmp
        open(joinpath(bot_path, "closed_orders.json"), "w") do f
            JSON.print(f, closed_orders, 4)
        end
    end
end
export sell_canceled

function sell_status_changed(tmp, k, v, id_str, config, open_orders, closed_orders, bot_path, client, usr, fees)
    open_orders["SELL"][k] = tmp
    open_orders["SELL"][k]["time"] = string(now())
    open(joinpath(bot_path, "open_orders.json"), "w") do f
        JSON.print(f, open_orders, 4)
    end
end
export sell_status_changed

function cancel_sell(tmp, k, v, id_str, config, open_orders, closed_orders, bot_path, client, usr, fees)
    tmp_sell = try client.cancel_order(
        symbol=config["symbol"],
        orderId=v["orderId"])
    catch
        nothing
    end
    if tmp_sell !== nothing || tmp["status"] == client.ORDER_STATUS_FILLED || tmp["status"] == client.ORDER_STATUS_CANCELED
        if parse(Float64, v2["executedQty"]) > 0
            closed_orders[id_str]["SELL-" * lpad(rand(1:1:10000), 5, "0")] = result
            open(joinpath(bot_path, "closed_orders.json"), "w") do f
                JSON.print(f, closed_orders, 4)
            end
        end
        delete!(open_orders["SELL"], k2)
    else
        @info("Cancel SELL failed")
    end
end
export cancel_sell


text:
"""
    ReadStatMeta

A collection of metadata parsed from a data file.
"""
struct ReadStatMeta
    labels::Dict{Symbol, String}
    formats::Dict{Symbol, String}
    val_label_keys::Dict{Symbol, String}
    val_label_dict::Dict{String, Dict{Any,String}}
    filelabel::String
    timestamp::DateTime
    fileext::String
end

"""
    varlabels(m::ReadStatMeta)
    varlabels(tb::ReadStatTable)

Retrieve the variable labels as a dictionary indexed by variable names.
"""
varlabels(m::ReadStatMeta) = m.labels

"""
    varformats(m::ReadStatMeta)
    varformats(tb::ReadStatTable)

Retrieve the variable format strings as a dictionary indexed by variable names.
"""
varformats(m::ReadStatMeta) = m.formats

"""
    val_label_keys(m::ReadStatMeta)
    val_label_keys(tb::ReadStatTable)

Retrieve the names of the collections of value labels applied to each variable
as a dictionary indexed by variable names.
"""
val_label_keys(m::ReadStatMeta) = m.val_label_keys

"""
    val_label_dict(m::ReadStatMeta)
    val_label_dict(tb::ReadStatTable)

Retrieve all collections of value labels stored in the data file
as a dictionary indexed by the names of the collections.
Each collection of value labels itself is again a dictionary
that maps the values to the associated labels.
"""
val_label_dict(m::ReadStatMeta) = m.val_label_dict

"""
    filelabel(m::ReadStatMeta)
    filelabel(tb::ReadStatTable)

Retrieve the label of the data file.
"""
filelabel(m::ReadStatMeta) = m.filelabel

"""
    filetimestamp(m::ReadStatMeta)
    filetimestamp(tb::ReadStatTable)

Retrieve the time stamp of the data file. (Time zone is not handled.)
"""
filetimestamp(m::ReadStatMeta) = m.timestamp

"""
    fileext(m::ReadStatMeta)
    fileext(tb::ReadStatTable)

Retrieve the file extension of the data file.
"""
fileext(m::ReadStatMeta) = m.fileext

varlabels(::Nothing) = nothing
varformats(::Nothing) = nothing
val_label_keys(::Nothing) = nothing
val_label_dict(::Nothing) = nothing
filelabel(::Nothing) = nothing
filetimestamp(::Nothing) = nothing
fileext(::Nothing) = nothing

Base.show(io::IO, m::ReadStatMeta) = print(io, typeof(m).name.name)
function Base.show(io::IO, ::MIME"text/plain", m::ReadStatMeta)
    println(io, typeof(m).name.name, ':')
    io = IOContext(io, :limit=>true)
    println(io, "  variable labels:    ", m.labels)
    println(io, "  variable formats:   ", m.formats)
    println(io, "  value label names:  ", m.val_label_keys)
    println(io, "  value labels:       ", m.val_label_dict)
    println(io, "  file label:         ", m.filelabel)
    println(io, "  file timestamp:     ", m.timestamp)
    print(io,   "  file extension:     ", m.fileext)
end

"""
    ReadStatTable <: Tables.AbstractColumns

A `Tables.jl`-compatible column table that collects data
from a Stata, SAS or SPSS file.
"""
struct ReadStatTable <: Tables.AbstractColumns
    columns::Vector{AbstractVector}
    names::Vector{Symbol}
    lookup::Dict{Symbol, Int}
    meta::Union{ReadStatMeta, Nothing}
    function ReadStatTable(columns::Vector{AbstractVector}, names::Vector{Symbol}, meta=nothing)
        lookup = Dict{Symbol, Int}(n=>i for (i, n) in enumerate(names))
        return new(columns, names, lookup, meta)
    end
end

# Accessor functions that are not intended to be exported
_columns(tb::ReadStatTable) = getfield(tb, :columns)
_names(tb::ReadStatTable) = getfield(tb, :names)
_lookup(tb::ReadStatTable) = getfield(tb, :lookup)

"""
    getmeta(tb::ReadStatTable)

Retrieve the metadata parsed from a data file.
"""
getmeta(tb::ReadStatTable) = getfield(tb, :meta)

varlabels(tb::ReadStatTable) = varlabels(getmeta(tb))
varformats(tb::ReadStatTable) = varformats(getmeta(tb))
val_label_keys(tb::ReadStatTable) = val_label_keys(getmeta(tb))
val_label_dict(tb::ReadStatTable) = val_label_dict(getmeta(tb))
filelabel(tb::ReadStatTable) = filelabel(getmeta(tb))
filetimestamp(tb::ReadStatTable) = filetimestamp(getmeta(tb))
fileext(tb::ReadStatTable) = fileext(getmeta(tb))

Tables.getcolumn(tb::ReadStatTable, i::Int) = _columns(tb)[i]
Tables.getcolumn(tb::ReadStatTable, n::Symbol) = _columns(tb)[_lookup(tb)[n]]
Tables.columnnames(tb::ReadStatTable) = copy(_names(tb))

Tables.schema(tb::ReadStatTable) =
    Tables.Schema{(_names(tb)...,), Tuple{(eltype(col) for col in _columns(tb))...}}()

Tables.columnindex(tb::ReadStatTable, n::Symbol) = _lookup(tb)[n]
Tables.columntype(tb::ReadStatTable, n::Symbol) = eltype(tb[n])

ncol(tb::ReadStatTable) = length(_names(tb))
nrow(tb::ReadStatTable) = ncol(tb) > 0 ? length(_columns(tb)[1])::Int : 0
Base.size(tb::ReadStatTable) = (nrow(tb), ncol(tb))
function Base.size(tb::ReadStatTable, i::Integer)
    if i == 1
        nrow(tb)
    elseif i == 2
        ncol(tb)
    else
        throw(ArgumentError("ReadStatTable only have two dimensions"))
    end
end

Tables.rowcount(tb::ReadStatTable) = nrow(tb)

Base.length(tb::ReadStatTable) = ncol(tb)
Base.isempty(tb::ReadStatTable) = size(tb, 1) == 0 || size(tb, 2) == 0

Base.values(tb::ReadStatTable) = _columns(tb)
Base.haskey(tb::ReadStatTable, key::Symbol) = haskey(_lookup(tb), key)
Base.haskey(tb::ReadStatTable, i::Int) = 0 < i <= length(_names(tb))

Base.show(io::IO, tb::ReadStatTable) = print(io, nrow(tb), '×', ncol(tb), " ReadStatTable")

function Base.show(io::IO, ::MIME"text/plain", tb::ReadStatTable)
    show(io, tb)
    if ncol(tb) > 0 && nrow(tb) > 0
        println(io, ':')
        nr, nc = displaysize(io)
        pretty_table(IOContext(io, :limit=>true, :displaysize=>(nr-3, nc)),
            tb, vlines=1:1, hlines=1:1, show_row_number=true,
            newline_at_end=false, vcrop_mode=:middle)
    end
end


text:
type SDF_MAXWELL_VD <: SDFDynamicSystem
    omg::Float64
    zta::Float64
    xi::Float64
    kpa::Float64
    alpha::Float64
    u0::Float64
    v0::Float64
    ud0::Float64

    SDF_MAXWELL_VD(omg::Float64, zta::Float64, xi::Float64, kpa::Float64, alpha::Float64) = new(omg, zta, xi, kpa, alpha, 0.0, 0.0, 0.0)
    SDF_MAXWELL_VD(omg::Float64, zta::Float64, xi::Float64, kpa::Float64) = new(omg, zta, xi, kpa, 1.0, 0.0, 0.0, 0.0)
end

function systemEquation(s::SDF_MAXWELL_VD, p::Function)
    function dydt(y::AbstractArray{Float64,1}, t::Float64)
        u, v, ud = y
        omg2 = s.omg*s.omg
        dvdt = p(t)-2.0*s.zta*s.omg*v-omg2*u-s.kpa*omg2*(u-ud)
        if s.alpha == 1.0
            duddt = s.kpa*s.omg/2.0/s.xi*(u-ud)
        else
            duddt = sign(u-ud)*abs(s.kpa*s.omg/2.0/s.xi*(u-ud))^(1.0/s.alpha)
        end
        [v; dvdt; duddt]
    end
    dydt
end

function response(s::SDF_MAXWELL_VD, p::Function, t::AbstractArray{Float64,1})
    y0 = [s.u0; s.v0; s.ud0]

    f = systemEquation(s, p)

    y = ode4(f, y0, t)

    u = y[:,1]
    v = y[:,2]
    ud = y[:,3]
    fdz = 2.0*s.zta*s.omg*v
    fdx = s.kpa*s.omg*s.omg*(u-ud)
    a = map(p,t)-fdz-fdx-s.omg*s.omg*u

    [u v a fdz fdx]
end

max_stars_repo_name: Panchatantra/DSR

max_issues_repo_name: Panchatantra/DSR

max_forks_repo_name: Panchatantra/DSR