text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $0))" --color=yes --startup-file=no -e 'include(popfirst!(ARGS))' \
    "${BASH_SOURCE[0]}" "$@"
    =#

#---Turing Machine-------------------------------------------------------------------------
# Adapted form https://rosettacode.org/wiki/Universal_Turing_machine#Julia

# define moving
"""
```julia
@enum Move Left=1 Stay Right
```

An enumerated type `Move`.
"""
@enum Move Left=1 Stay Right

"""
```julia
mutable struct MachineState <: MachineComponent
```

A snapshot of a state of a Turing machine.

Fields:
  - `state::String`: The name of the state of the Turing machine at the current machine state.
  - `tape::Dict{Int, String}``: The current state of the tape of the machine.  The head will be at `Int` value, and the contents at `String` value.
  - `headpos::Int`: The position of the head at the current state of the Turing Machine.
"""
mutable struct MachineState <: MachineComponent
    state::String
    tape::Dict{Int, String}
    headpos::Int
end

"""
```julia
struct Rule <: MachineComponent
```

Fields:
  - `instate::String`: The state name (`String`) where the rule starts.
  - `outstate::String`: The state name (`String`) where the rule ends.
  - `s1::String`: What the Turing machine will read.
  - `s2::String`: What the Turing machine will write.
  - `move::String`: How the Turing machine will move (see `Move`).

Note: `instate` and `outstate` are connected nodes in the Turing Machine.
"""
struct Rule <: MachineComponent
    instate::String
    outstate::String
    s1::String
    s2::String
    move::Move
end

"""
```julia
struct TMProgramme <: TuringMachine
```

Fields:
  - `title::String`: A descriptive name for your programme.
  - `initial::String`: The state name where the Turing machine begins.
  - `final::String`: The state name where the Turing machine halts.
  - `blank::String`: The symbol used to denote a blank cell.
  - `Rules::Vector{Rule}``: A list of rules of the programme.
"""
struct TMProgramme <: TuringMachine
    title::String
    initial::String
    final::String
    blank::String
    rules::Vector{Rule}
end

"""
```julia
Base.show(io::IO, mstate::MachineState)
```

Prints the current state of a Turing machine.
"""
function Base.show(io::IO, mstate::MachineState)
    ibracket(i, curpos, val) = isequal(i, curpos) ? "[$val]" : " $val "
    print(io, rpad("($(mstate.state))", 12))
    for i in sort(collect(keys(mstate.tape)))
        print(io, "   $(ibracket(i, mstate.headpos, mstate.tape[i]))")
    end
end

"""
```julia
run_turing_machine(TMProgramme, tape, verbose)
```

Run a specified Turing programme on a specified tape.  `verbose` is a boolean.
"""
function run_turing_machine(TMProgramme, tape, verbose)
    println("\u001b[1;38m$(TMProgramme.title)\u001b[0;38m")
    verbose && println(" State\t\t\tTape [head]\n", "-"^displaysize(stdout)[2])

    mstate = MachineState(TMProgramme.initial, tape, 1)
    stepcount = 0
    while true
        if ! haskey(mstate.tape, mstate.headpos)
            mstate.tape[mstate.headpos] = TMProgramme.blank
        end

        verbose && println(mstate)

        for rule in TMProgramme.rules
            if isequal(rule.instate, mstate.state) && isequal(rule.s1, mstate.tape[mstate.headpos])
                mstate.tape[mstate.headpos] = rule.s2
                if isequal(rule.move, Left)
                    mstate.headpos -= 1
                elseif isequal(rule.move, Right)
                    mstate.headpos += 1
                end

                mstate.state = rule.outstate
                break
            end
        end

        stepcount += 1

        if isequal(mstate.state, TMProgramme.final)
            break
        end
    end

    verbose && println(mstate)
    println("Total number of steps taken: $stepcount")
end

#--Register Machines--------------------------------------------------------------------------

"""
```julia
mutable struct RegisterMachine <: Machine
```

A register machine.  That is, a machine with finitely many registers, whose registers contain a natural number n ∈ ℕ ∪ {0}.

Fields:
  - `contents::AbastractArray`: The contents of the registers.

---
```julia
RegisterMachine(contents::AbstractArray)
```

A constructor function for the `struct` `RegisterMachine`.  Ensures the contents of the machine are indeed all natural numbers

---

```julia
RegisterMachine(T::Union{Tuple, UnitRange})
RegisterMachine(a::Integer, b::Integer...)
```

Constructor functions for the `struct` `RegisterMachine`.  Your contents are allowed to be tuples, ranges, or a series of integers.
"""
mutable struct RegisterMachine <: Machine
    contents::AbstractArray#Vector{<:Integer}

    function RegisterMachine(contents::AbstractArray)
        anylessthan(0, contents) && throw(error("Registers must contain non-negative numbers."))

        new(contents)
    end

    RegisterMachine(T::Union{Tuple, UnitRange}) = RegisterMachine([T...])
    RegisterMachine(a::Integer...) = RegisterMachine([a...])
end

"""
```julia
run_goto_programme(P::GoToProgramme, R::RegisterMachine) -> RegisterMachine.contents
run_goto_programme(P::GoToProgramme, T::Tuple) -> RegisterMachine.contents
run_goto_programme(P::Integer, T::Tuple) -> RegisterMachine.contents
run_goto_programme(P::Integer, R::RegisterMachine) -> RegisterMachine.contents
```

Takes in a *coded* programme `P` (for coding a programme, see notes), and runs the programme using register machine R.  It a tuple is given, will convert to a register machine whose contents is the tuple.  If an integer is given, will convert to a goto programme.

Notes: There are five possible instructions for a GoTo programme:
  - Increment Rᵢ  ⟵ `pair_tuple(0, i)`
  - Decrement Rᵢ  ⟵ `pair_tuple(1, i)`
  - Goto line k  ⟵ `pair_tuple(2, k)`
  - If Rᵢ is zero, goto line k  ⟵ `pair_tuple(3, (i, k))`
  - Halt  ⟵ `pair_tuple(4, 0)`

Then, the programme P, consisting of instructions I₁, I₂, ..., Iᵢ, is coded by

```julia
pair_tuple(i, pair_tuple(I₁, I₂, ..., Iᵢ))
```

For utilities regarding these instructions, see `pair_tuple`, `Sequence`, `Instruction`, `increment`, `decrement`, `goto`, `ifzero_goto`, `halt`, and `GoToProgramme`.

---

### Examples

```julia
julia> run_goto_programme(363183787614755732766753446033240)
(1, 0)

julia> run_goto_programme(363183787614755732766753446033240, Register(0, 0, 0 ,0))
(1, 0, 0, 0)
```
"""
function run_goto_programme(P::GoToProgramme, R::RegisterMachine)
    line_number = 0

    for instruction in P.instructions
        primary_identifier, secondary_identifier = instruction

        if isequal(primary_identifier, __increment_identifier)
            n = secondary_identifier
            R.contents[n + 1] += 1
        elseif isequal(primary_identifier, __decrement_identifier)
            n = secondary_identifier
            R.contents[n + 1] = R.contents[n + 1] ∸ 1
        elseif isequal(primary_identifier, __goto_identifier)
            k = secondary_identifier
            line_number = k
            continue
        elseif isequal(primary_identifier, __ifzero_goto_identifier)
            n, k = secondary_identifier
            if iszero(R.contents[n + 1])
                line_number = k
                continue
            end
        elseif isequal(primary_identifier, halt().first) && isequal(secondary_identifier, halt().second)
            break
        end

        line_number += 1
    end

    return tuple(R.contents...)
end
function run_goto_programme(P::GoToProgramme)
    max_register = extrema_tuple(P.instructions)[2][2] # the maximum value in the list of instructions in the second position
    R = RegisterMachine(zeros(Integer, max_register+1)) # fill register with zeros

    return run_goto_programme(P, R)
end
run_goto_programme(P::GoToProgramme, T::Tuple) = run_goto_programme(P, RegisterMachine(T))
run_goto_programme(P::Integer, T::Tuple) = run_goto_programme(GoToProgramme(P), RegisterMachine(T))
run_goto_programme(P::Integer, R::RegisterMachine) = run_goto_programme(GoToProgramme(P), R)
run_goto_programme(P::Integer) = run_goto_programme(GoToProgramme(P))


text:
"""
    Hidden Markov Mixture Model.
"""
struct H3M{T <: AbstractHMM}
    M::Vector{T}
    ω::Vector{Float64}
    H3M{T}(M, ω) where T = assert_h3m(M, ω) && new(M, ω)
end

H3M(M::AbstractVector{T}, ω) where T  = H3M{T}(M, ω)
H3M(M::AbstractVector{T}) where T = H3M{T}(M, ones(length(M))/length(M))

length(m::H3M) = length(m.M)
iterate(m::H3M, args...) = iterate(zip(m.M, m.ω), args...)

function assert_h3m(M, ω)
    @argcheck length(M) == length(ω)
    @argcheck isprobvec(ω)
    return true
end

ext: jl


text:
# Autogenerated wrapper script for Clang_assert_jll for x86_64-w64-mingw32-cxx03
export clang, libclang, libclang_cpp

using Zlib_jll
using libLLVM_assert_jll
JLLWrappers.@generate_wrapper_header("Clang_assert")
JLLWrappers.@declare_executable_product(clang)
JLLWrappers.@declare_library_product(libclang, "libclang.dll")
JLLWrappers.@declare_library_product(libclang_cpp, "libclang-cpp.dll")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll, libLLVM_assert_jll)
    JLLWrappers.@init_executable_product(
        clang,
        "tools\\clang.exe",
    )

    JLLWrappers.@init_library_product(
        libclang,
        "bin\\libclang.dll",
        nothing,
    )

    JLLWrappers.@init_library_product(
        libclang_cpp,
        "bin\\libclang-cpp.dll",
        nothing,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using PatModules
@auxmodule! utils () begin
    # Import external modules here
    # (not that Base is ever external, this is
    # just for illustration)
    import Base

    export say

    function say(msg)
        Base.println(msg)
    end
end


text:
######
# can pass a Dict to kwarg "userdata" in the solve function, this can be accessed via integrator.opts.userdata
# may be useful for something
######

# integrator options
int_options = (abstol=1e-5, reltol=1e-5, verbose=false)

# function to run simulation
objective(x) =
    try
        cost(simulate(x, prob)...)
    catch
        return 10000.0
    end

# simulate a stride from prob
function simulate(prob)
    # integrate step 1 to touchdown of step 2
    sol1 = solve(prob, Tsit5(), callback=cb1, saveat=0.001; int_options...)

    # check if a second step will occur
    # pocmy(sol1, sol1.t[end]) ≤ com_drop && return 2e6 # or cost(sol1)
    vocmy(sol1, sol1.t[end]) ≤ -1.0 && return sol1, -1
    sol1.t[end] < 0.1 && return sol1, -1
    sol1.retcode ≠ :Terminated && return sol1, -1

    # fit splines to grf
    T = sol1.t[end] # get last time which is multiple of saveat in sol1
    t = (sol1.t[1]:0.001:T) .+ T # callback means sol.t =/= 0:0.001:sol.t[end]

    # but need evenly spaced t for spline
    splX1 = CubicSplineInterpolation(t, [rx1(sol1, _t) for _t in t .- T], extrapolation_bc=0.0)
    splY1 = CubicSplineInterpolation(t, [ry1(sol1, _t) for _t in t .- T], extrapolation_bc=0.0)
    splX2 = CubicSplineInterpolation(t, [rx2(sol1, _t) for _t in t .- T], extrapolation_bc=0.0)
    splY2 = CubicSplineInterpolation(t, [ry2(sol1, _t) for _t in t .- T], extrapolation_bc=0.0)

    _vrx1(t)::Float64 = splX1(t)
    _vry1(t)::Float64 = splY1(t)
    _vrx2(t)::Float64 = splX2(t)
    _vry2(t)::Float64 = splY2(t)

    _virtual_force = VirtualForce(_vrx1, _vry1, _vrx2, _vry2)

    # update parameters and initial conditions for step 2
    pnew = setproperties(prob.p, (virtual_force = _virtual_force))
    unew = sol1.u[end]
    tnew = (sol1.t[end], 0.484) # TODO: make the end of integration time a global variable?
    newprob = remake(prob, u0=unew, p=pnew, tspan=tnew)

    # integrate step 2
    saveats = sol1.t[end-1]+0.001:0.001:tnew[2] # start saving at next value after end of last integration
    sol2 = solve(newprob, Tsit5(), saveat=saveats; int_options...)
    sol2.retcode ≠ :Success && return sol1, -1

    return sol1, sol2
end

# run a simulate with new parameters
function simulate(x, prob)
    # new parameter struct with updated parameters from x and reset torque generators, updated intial CC angles
    pnew, unew = updateParameters(prob.p, x, prob.u0)

    # remake problem with updated parameters and initial conditions
    newprob = remake(prob, u0=unew, p=pnew)

    # integrate new problem
    return simulate(newprob)
end
simulate(x::Vector{Float64}) = simulate(x, prob)

# cost function
# cost(sol) = 10(abs(VCMX - step_velocity(sol))) + abs(TSW - swing_time(sol))
# mean squared error between two vectors
mse(x, y) = (x .- y) .^ 2 |> mean
# function cost(sol)::Float64
#     global matching_data

#     ## match hip and knee angles
#     N = size(sol, 2)

#     # calculate hip and knee angles
#     θhip = π .+ view(sol, 7, :) .|> rad2deg
#     θknee = π .- view(sol, 6, :) .|> rad2deg

#     # calculate mse
#     hip_mse = mse(θhip, view(matching_data[:lhip], 1:N))
#     knee_mse = mse(θknee, view(matching_data[:lknee], 1:N))

#     # scale by range
#     sf_hip = view(matching_data[:lhip], :) |> extrema |> collect |> diff
#     sf_knee = view(matching_data[:lknee], :) |> extrema |> collect |> diff
#     hip_mse /= sf_hip[1]
#     knee_mse /= sf_knee[1]

#     return hip_mse + knee_mse

#     # ## mse between start and end of simulation
#     # weights = [0.0, 0.0, 1.0, 0.1, 0.5, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.25, 0.7, 0.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # weights vector
#     # return ((sol.u[end] .- u₀).^2) .* weights |> mean
# end

cost(sol1, sol2) = cost([Array(sol1(0:0.001:sol1.t[end])) Array(sol2)])
cost(sol1, err::Integer) = 10000.0

# FORCING SIMULATION TO BE SAME AS DATA
function cost(sim_data)::Float64
    global matching_data
    # calculate hip and knee angles
    θhip = π .+ view(sim_data, 7, :) .|> rad2deg
    θknee = π .- view(sim_data, 6, :) .|> rad2deg
    θhat = view(sim_data, 3, :) .|> rad2deg
    θankle = π .+ view(sim_data, 5, :) .|> rad2deg

    # calculate mse
    hip_mse = mse(θhip, matching_data[:lhip])
    knee_mse = mse(θknee, matching_data[:lknee])
    hat_mse = mse(θhat, matching_data[:lhat])
    ankle_mse = mse(θankle, matching_data[:lankle])

    # scale by range
    sf_hip = view(matching_data[:lhip], :) |> extrema |> collect |> diff
    sf_knee = view(matching_data[:lknee], :) |> extrema |> collect |> diff
    sf_hat = view(matching_data[:lhat], :) |> extrema |> collect |> diff
    hip_mse /= sf_hip[1]
    knee_mse /= sf_knee[1]
    hat_mse /= sf_hat[1]

    return hip_mse + knee_mse + hat_mse + 0.1ankle_mse

    # time cost
    # time_cost = 1000abs(sol.t[end] - 0.484)

    # periodicity cost
    # weights = [0.0, 1.0, 0.8, 0.2, 0.5, 1.0, 1.0, 0.5, 0.5, 0.4, 0.1, 0.25, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # weighting of state components
    # periodicity_cost = (sol.prob.u0 .- sol.u[end]) .^ 2 .* weights |> mean

    # return angles_cost #+ 0.5time_cost #+ 0.8periodicity_cost

end

## vector to parameters to be optimised and resets torque generators to inital values
function updateParameters(p, x, u₀)
    @unpack he, ke, ae, hf, kf, af = p
    @inbounds q1, q2, q3, q4, q5, q6, q7, u1, u2, u3, u4, u5, u6, u7 = u₀

    θh = π + q7
    θk = π - q6
    θa = π + q5

    ωh = u7
    ωk = -u6
    ωa = u5

    # manually separate out
    # he_α_p = ActivationProfile(x[1:7])
    # ke_α_p = ActivationProfile(x[8:14])
    # ae_α_p = ActivationProfile(x[15:21])
    # hf_α_p = ActivationProfile(x[22:28])
    # kf_α_p = ActivationProfile(x[29:35])
    # af_α_p = ActivationProfile(x[36:42])
    he_α_p = ActivationProfile(x[1:10])
    ke_α_p = ActivationProfile(x[11:20])
    ae_α_p = ActivationProfile(x[21:30])
    hf_α_p = ActivationProfile(x[31:40])
    kf_α_p = ActivationProfile(x[41:50])
    af_α_p = ActivationProfile(x[51:60])
    # he_α_p = ActivationProfile(x[1:13])
    # ke_α_p = ActivationProfile(x[14:26])
    # ae_α_p = ActivationProfile(x[27:39])
    # hf_α_p = ActivationProfile(x[40:52])
    # kf_α_p = ActivationProfile(x[53:65])
    # af_α_p = ActivationProfile(x[66:78])
    _k1, _k2, _k3, _k4, _k5, _k6, _k7, _k8 = x[61:68]

    # torque generators with inital values and new activation parameters
    _he = TorqueGenerator(2π - θh, -ωh, he.cc.cc_p, he.sec.sec_p, he_α_p)
    _ke = TorqueGenerator(2π - θk, -ωk, ke.cc.cc_p, ke.sec.sec_p, ke_α_p)
    _ae = TorqueGenerator(2π - θa, -ωa, ae.cc.cc_p, ae.sec.sec_p, ae_α_p)
    _hf = TorqueGenerator(θh, ωh, hf.cc.cc_p, hf.sec.sec_p, hf_α_p)
    _kf = TorqueGenerator(θk, ωk, kf.cc.cc_p, kf.sec.sec_p, kf_α_p)
    _af = TorqueGenerator(θa, ωa, af.cc.cc_p, af.sec.sec_p, af_α_p)

    # get new initial CC angles
    _u₀ = SVector(u₀[1:14]..., map(x -> x.cc.θ, [_he, _ke, _ae, _hf, _kf, _af])...)

    # returns new Params struct with torque generators ready for next simulation
    _p = setproperties(p, (he=_he, ke=_ke, ae=_ae, hf=_hf, kf=_kf, af=_af, k1=_k1, k2=_k2, k3=_k3, k4=_k4, k5=_k5, k6=_k6, k7=_k7, k8=_k8))

    return _p, _u₀
end

# reutns parameters from an Params struct
get_parameters(p) = mapreduce(tqgen -> tqgen.cc.α_p |> Vector, vcat, [p.he, p.ke, p.ae, p.hf, p.kf, p.af])


text:
include("TypeDefine.jl")

#see (4.46) from Toro
function p_guess(left::HydroStatus, right::HydroStatus, ::TwoRarefaction)
		@assert left.gamma == right.gamma
		gamma = left.gamma
		gamma_power = (gamma-1) / 2. / gamma
		delta_u = right.u - left.u
		return ( (left.c + right.c - 0.5 * (gamma - 1) * delta_u) /
			   (left.c/left.p^gamma_power + right.c/right.p^gamma_power) ) ^
		       (1. / gamma_power)
end

#see (4.47) from Toro
function p_guess(left::HydroStatus, right::HydroStatus, ::PrimitiveValue)
	delta_u = right.u - left.u
	mean_pressure = 0.5 * (left.p + right.p)
	mean_density = 0.5 * (left.rho + right.rho)
	mean_speed_of_sound = 0.5 * (left.c + right.c)
	return mean_pressure - delta_u * mean_density * mean_speed_of_sound
end

#see (4.48) from Toro
g(p, status::HydroStatus) = sqrt(status.A / (p + status.B))
function p_guess(left::HydroStatus, right::HydroStatus, ::TwoShock)
	p_hat = p_guess(left, right, PrimitiveValue)
	g_left = g(p_hat, left)
	g_right = g(p_hat, right)
	delta_u = right.u - left.u
	return (g_right * left.p + g_left * right.p - delta_u)/(g_right + g_left)
end

#see (4.49) from Toro
p_guess(left::HydroStatus, right::HydroStatus, ::MeanPressure) = 0.5*(left.p + right.p)

max_stars_repo_name: ryarazi/ExactRiemannProblemSolver

max_issues_repo_name: ryarazi/ExactRiemannProblemSolver

max_forks_repo_name: ryarazi/ExactRiemannProblemSolver


text:
"""
Split strings by an specified symbol.

# Arguments
 - **text** (String) Text to split
 - **symbol** (String) Symbol to use for spliting

# Returns
Array{String}

# Example
```julia
> text =
"
a = b
+ 3
c = d
"
> splitLines(text,"=")
["a = b + 3","c = b"]
```
"""
function splitLines(text::String,symbol::String)
    list=[string(i) for i in split(text,"\n") if i != ""]
    listN=[]
    for e in list
        if !occursin(symbol, e) #Append if just jump in line
            listN[end] = string(listN[end]," ",e)
        else
            append!(listN,[e])
        end
    end
    return listN
end

function splitUpdating(text::String)
    list=[string(i) for i in split(text,"\n") if i != ""]
    listN=[]
    for e in list
        app = true
        for i in Array(["+=";"-=";"*=";"/=";"÷=";"%=";"^=";" = "])
            if occursin(i, e)
                app = false
                append!(listN,[e])
                break
            end
        end
        if app
            listN[end] = string(listN[end]," ",e)
        end
    end
    return listN
end

"""
Recurrent function auxiliar to the splitEqs for the block of code adaptation
"""
function splits(income,level,l)
    out = nothing
    if income == :dt && level == 0
        l[1] = 1
    elseif income == :(-dt) && level == 0
        l[1] = -1
    elseif income == :dW && level == 0
        l[2] = 1
    elseif income == :(-dW) && level == 0
        l[2] = -1
    else
        args = income.args
        if args[1] == :+ || args[1] == :-
            for i in 1:length(args)
                if args[i] == :dt
                    args[i] = 1
                    out = "t"
                elseif args[i] == :dW
                    args[i] = 1
                    out = "W"
                elseif typeof(args[i]) == Expr
                    args[i],out = splits(args[i],level+1,l)
                end
                if  out == "t"
                    if l[1] == 0
                        l[1] = args[i]
                    else
                        if args[1] == :+
                            l[1] = :($(l[1])+$(args[i]))
                        else
                            l[1] = :($(l[1])-$(args[i]))
                        end
                    end
                    out = nothing
                elseif  out == "W"
                    if l[2] == 0
                        l[2] = args[i]
                    else
                        if args[2] == :+
                            l[2] = :($(l[2])+$(args[i]))
                        else
                            l[2] = :($(l[2])-$(args[i]))
                        end
                    end
                    out = nothing
                end
            end
        elseif args[1] == :*
            if :dt in args
                if length(args) > 3
                    income.args = args[[i for i in 1:length(args) if args[i]!=:dt]]
                    out = "t"
                else
                    income = :($(args[[i for i in 2:length(args) if args[i]!=:dt]][1]))
                    out = "t"
                end
            elseif :dW in args
                if length(args) > 3
                    income.args = args[[i for i in 1:length(args) if args[i]!=:dW]]
                    out = "W"
                else
                    income = :($(args[[i for i in 2:length(args) if args[i]!=:dW]][1]))
                    out = "W"
                end
            end
            if level == 0
                if  out == "t"
                    if l[1] == 0
                        l[1] = income
                    else
                        if args[1] == :+
                            l[1] = :($(l[1])+$(income))
                        else
                            l[1] = :($(l[1])-$(income))
                        end
                    end
                    out = nothing
                elseif  out == "W"
                    if l[2] == 0
                        l[2] = income
                    else
                        if args[2] == :+
                            l[2] = :($(l[2])+$(income))
                        else
                            l[2] = :($(l[2])-$(income))
                        end
                    end
                    out = nothing
                end
            end
        end
    end

    return income, out
end

"""
Prepares a block of code with differential equations to be vectorized and introduced in the integration algorithms.
Returns the adapted block of code and the list of the lines where there has been declared a Wienner random variable.

# Arguments
 - **value** (Expr) Block of code with differential equations

# Returns
Expr, Array{Int}
"""
function splitEqs(value)
    list = quote end
    nRand = []
    count = 1
    for i in 1:length(value.args)
        try value.args[i].head
            eq = value.args[i]
            l = Any[0,0]
            args,out = splits(eq.args[2],0,l)
            eq.args[2] = args
            if l == [0,0]
                push!(list.args,eq)
                count += 1
            else
                var = Meta.parse(string(eq.args[1])[2:end])
                varf = Meta.parse(string(var,"f_"))
                varg = Meta.parse(string(var,"g_"))
                ineq = nothing
                if l[1] != 0
                    push!(list.args,:($varf = $(l[1])))
                    ineq = :($varf*dt)
                    count += 1
                end
                if l[2] != 0
                    push!(list.args,:($varg = $(l[2])))
                    if ineq != nothing
                        ineq = :($ineq+$varg*dW)
                    else
                        ineq = :($varg*dW)
                    end
                    count += 1
                end
                push!(list.args,:($(eq.args[1]) = $ineq))
                count += 1
                if nRand == []
                    push!(nRand,count)
                    count = 0
                else
                    push!(nRand,nRand[end]+count)
                    count = 0
                end
            end
        catch
        end
    end

    return list, nRand
end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
@testset "code generation" begin

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

@testset "LLVM IR" begin

@testset "basic reflection" begin
    @eval llvm_valid_kernel() = return
    @eval llvm_invalid_kernel() = 1

    ir = sprint(io->CUDAnative.code_llvm(io, llvm_valid_kernel, Tuple{}; optimize=false, dump_module=true))

    # module should contain our function + a generic call wrapper
    @test occursin("define void @julia_llvm_valid_kernel", ir)
    @test !occursin("define %jl_value_t* @jlcall_", ir)

    @test CUDAnative.code_llvm(devnull, llvm_invalid_kernel, Tuple{}) == nothing
    @test_throws CUDAnative.CompilerError CUDAnative.code_llvm(devnull, llvm_invalid_kernel, Tuple{}; kernel=true) == nothing
end

@testset "exceptions" begin
    @eval codegen_exception() = throw(DivideError())
    ir = sprint(io->CUDAnative.code_llvm(io, codegen_exception, Tuple{}))

    # exceptions should get lowered to a plain trap...
    @test occursin("llvm.trap", ir)
    # not a jl_throw referencing a jl_value_t representing the exception
    @test !occursin("jl_value_t", ir)
    @test !occursin("jl_throw", ir)
end

@testset "sysimg" begin
    # bug: use a system image function

    @eval function codegen_call_sysimg(a,i)
        Base.pointerset(a, 0, mod1(i,10), 8)
    end

    ir = sprint(io->CUDAnative.code_llvm(io, codegen_call_sysimg, Tuple{Ptr{Int},Int}))
    @test !occursin("jlsys_", ir)
end

@testset "child functions" begin
    # we often test using `@noinline sink` child functions, so test whether these survive
    @eval @noinline codegen_child(i) = sink(i)
    @eval codegen_parent(i) = codegen_child(i)

    ir = sprint(io->CUDAnative.code_llvm(io, codegen_parent, Tuple{Int}))
    @test_broken occursin(r"call .+ @julia_codegen_child_", ir)
end

@testset "JuliaLang/julia#21121" begin
    @eval function codegen_tuple_leak()
        weight_matrix = @cuStaticSharedMem(Float32, (16, 16))
        sync_threads()
        weight_matrix[1, 16] *= 2
        sync_threads()
    end

    ir = sprint(io->CUDAnative.code_llvm(io, codegen_tuple_leak, Tuple{}))
    @test !occursin("inttoptr", ir)
end

@testset "kernel functions" begin
@testset "wrapper function aggregate rewriting" begin
    @eval codegen_aggregates(x) = nothing

    @eval struct Aggregate
        x::Int
    end

    ir = sprint(io->CUDAnative.code_llvm(io, codegen_aggregates, Tuple{Aggregate}))
    @test occursin(r"@julia_codegen_aggregates_\d+\({ i64 }( addrspace\(\d+\))?\*", ir)

    ir = sprint(io->CUDAnative.code_llvm(io, codegen_aggregates, Tuple{Aggregate}; kernel=true))
    @test occursin(r"@ptxcall_codegen_aggregates_\d+\({ i64 }\)", ir)
end

@testset "property_annotations" begin
    ir = sprint(io->CUDAnative.code_llvm(io, llvm_valid_kernel, Tuple{}; dump_module=true))
    @test !occursin("nvvm.annotations", ir)

    ir = sprint(io->CUDAnative.code_llvm(io, llvm_valid_kernel, Tuple{};
                                         dump_module=true, kernel=true))
    @test occursin("nvvm.annotations", ir)
    @test !occursin("maxntid", ir)
    @test !occursin("reqntid", ir)
    @test !occursin("minctasm", ir)
    @test !occursin("maxnreg", ir)

    ir = sprint(io->CUDAnative.code_llvm(io, llvm_valid_kernel, Tuple{};
                                         dump_module=true, kernel=true, maxthreads=42))
    @test occursin("maxntidx\", i32 42", ir)
    @test occursin("maxntidy\", i32 1", ir)
    @test occursin("maxntidz\", i32 1", ir)

    ir = sprint(io->CUDAnative.code_llvm(io, llvm_valid_kernel, Tuple{};
                                         dump_module=true, kernel=true, minthreads=42))
    @test occursin("reqntidx\", i32 42", ir)
    @test occursin("reqntidy\", i32 1", ir)
    @test occursin("reqntidz\", i32 1", ir)

    ir = sprint(io->CUDAnative.code_llvm(io, llvm_valid_kernel, Tuple{};
                                         dump_module=true, kernel=true, blocks_per_sm=42))
    @test occursin("minctasm\", i32 42", ir)

    ir = sprint(io->CUDAnative.code_llvm(io, llvm_valid_kernel, Tuple{};
                                         dump_module=true, kernel=true, maxregs=42))
    @test occursin("maxnreg\", i32 42", ir)
end
end

@testset "LLVM D32593" begin
    @eval struct llvm_D32593_struct
        foo::Float32
        bar::Float32
    end

    @eval llvm_D32593(arr) = arr[1].foo

    CUDAnative.code_llvm(devnull, llvm_D32593, Tuple{CuDeviceVector{llvm_D32593_struct,AS.Global}})
end

@testset "kernel names" begin
    @eval codegen_regular() = nothing
    @eval codegen_closure = ()->nothing

    function test_name(f, name; kwargs...)
        code = sprint(io->CUDAnative.code_llvm(io, f, Tuple{}; kwargs...))
        @test occursin(name, code)
    end

    test_name(codegen_regular, "julia_codegen_regular")
    test_name(codegen_regular, "julia_codegen_renamed"; alias="codegen_renamed")

    test_name(codegen_regular, "ptxcall_codegen_regular"; kernel=true)
    test_name(codegen_regular, "ptxcall_codegen_renamed"; kernel=true, alias="codegen_renamed")

    test_name(codegen_closure, "julia_anonymous")
    test_name(codegen_closure, "julia_codegen_renamed"; alias="codegen_renamed")

    test_name(codegen_closure, "ptxcall_anonymous"; kernel=true)
    test_name(codegen_closure, "ptxcall_codegen_renamed"; kernel=true, alias="codegen_renamed")
end

end

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

@testset "PTX assembly" begin

@testset "basic reflection" begin
    @eval ptx_valid_kernel() = nothing
    @eval ptx_invalid_kernel() = 1

    @test CUDAnative.code_ptx(devnull, ptx_valid_kernel, Tuple{}) == nothing
    @test CUDAnative.code_ptx(devnull, ptx_invalid_kernel, Tuple{}) == nothing
    @test_throws CUDAnative.CompilerError CUDAnative.code_ptx(devnull, ptx_invalid_kernel, Tuple{}; kernel=true) == nothing
end

@testset "child functions" begin
    # we often test using @noinline child functions, so test whether these survive
    # (despite not having side-effects)
    @eval @noinline ptx_child(i) = sink(i)
    @eval function ptx_parent(i)
        ptx_child(i)
        return
    end

    asm = sprint(io->CUDAnative.code_ptx(io, ptx_parent, Tuple{Int64}))
    @test_broken occursin(r"call.uni\s+julia_ptx_child_"m, asm)
end

@testset "kernel functions" begin
    @eval @noinline ptx_nonentry(i) = sink(i)
    @eval ptx_entry(i) = (ptx_nonentry(i); nothing)

    asm = sprint(io->CUDAnative.code_ptx(io, ptx_entry, Tuple{Int64}; kernel=true))
    @test occursin(r"\.visible \.entry ptxcall_ptx_entry_", asm)
    @test !occursin(r"\.visible \.func julia_ptx_nonentry_", asm)
    @test_broken occursin(r"\.func julia_ptx_nonentry_", asm)

@testset "property_annotations" begin
    asm = sprint(io->CUDAnative.code_ptx(io, ptx_entry, Tuple{Int64}; kernel=true))
    @test !occursin("maxntid", asm)

    asm = sprint(io->CUDAnative.code_ptx(io, ptx_entry, Tuple{Int64};
                                         kernel=true, maxthreads=42))
    @test occursin(".maxntid 42, 1, 1", asm)

    asm = sprint(io->CUDAnative.code_ptx(io, ptx_entry, Tuple{Int64};
                                         kernel=true, minthreads=42))
    @test occursin(".reqntid 42, 1, 1", asm)

    asm = sprint(io->CUDAnative.code_ptx(io, ptx_entry, Tuple{Int64};
                                         kernel=true, blocks_per_sm=42))
    @test occursin(".minnctapersm 42", asm)

    if LLVM.version() >= v"4.0"
        asm = sprint(io->CUDAnative.code_ptx(io, ptx_entry, Tuple{Int64};
                                             kernel=true, maxregs=42))
        @test occursin(".maxnreg 42", asm)
    end
end
end

@testset "IR validation" begin
@testset "delayed lookup" begin
    @eval codegen_ref_nonexisting() = nonexisting
    try
         CUDAnative.code_ptx(codegen_ref_nonexisting, Tuple{})
    catch err
        @test isa(err, CUDAnative.CompilerError)
        msg = sprint(io->showerror(io, err))
        @test occursin("could not compile codegen_ref_nonexisting() for GPU; unsupported LLVM IR", msg)
        @test occursin("Reason: unsupported call to the Julia runtime", msg)
        if VERSION >= v"0.7.0-alpha.37"
            @test occursin(r"\[1\] codegen_ref_nonexisting at .*codegen.jl", msg)
        end
    end
end

@testset "generic call" begin
    @eval codegen_call_nonexisting() = nonexisting()
    try
         CUDAnative.code_ptx(codegen_call_nonexisting, Tuple{})
    catch err
        @test isa(err, CUDAnative.CompilerError)
        msg = sprint(io->showerror(io, err))
        @test occursin("could not compile codegen_call_nonexisting() for GPU; unsupported LLVM IR", msg)
        @test occursin("Reason: unsupported call to the Julia runtime", msg)
        if VERSION >= v"0.7.0-alpha.37"
            @test occursin(r"\[1\] codegen_call_nonexisting at .*codegen.jl", msg)
        end
    end
end

@testset "compile traces" begin
    @eval codegen_call_inner() = error()
    @eval codegen_call_outer() = codegen_call_inner()
    try
         CUDAnative.code_ptx(codegen_call_outer, Tuple{})
    catch err
        @test isa(err, CUDAnative.CompilerError)
        msg = sprint(io->showerror(io, err))
        if VERSION >= v"0.7.0-alpha.37"
            @test occursin(r"\[1\] codegen_call_inner at .*codegen.jl", msg)
            @test occursin(r"\[2\] codegen_call_outer at .*codegen.jl", msg)
        end
    end
end
end

@testset "idempotency" begin
    # bug: generate code twice for the same kernel (jl_to_ptx wasn't idempotent)

    @eval codegen_idempotency() = nothing
    CUDAnative.code_ptx(devnull, codegen_idempotency, Tuple{})
    CUDAnative.code_ptx(devnull, codegen_idempotency, Tuple{})
end

@testset "child function reuse" begin
    # bug: depending on a child function from multiple parents resulted in
    #      the child only being present once

    @eval @noinline codegen_child_reuse_child(i) = sink(i)
    @eval function codegen_child_reuse_parent1(i)
        codegen_child_reuse_child(i)
        return
    end

    asm = sprint(io->CUDAnative.code_ptx(io, codegen_child_reuse_parent1, Tuple{Int}))
    @test_broken occursin(r".func julia_codegen_child_reuse_child_", asm)

    @eval function codegen_child_reuse_parent2(i)
        codegen_child_reuse_child(i+1)
        return
    end

    asm = sprint(io->CUDAnative.code_ptx(io, codegen_child_reuse_parent2, Tuple{Int}))
    @test_broken occursin(r".func julia_codegen_child_reuse_child_", asm)
end

@testset "child function reuse bis" begin
    # bug: similar, but slightly different issue as above
    #      in the case of two child functions
    @eval @noinline codegen_child_reuse_bis_child1(i) = sink(i)
    @eval @noinline codegen_child_reuse_bis_child2(i) = sink(i+1)
    @eval function codegen_child_reuse_bis_parent1(i)
        codegen_child_reuse_bis_child1(i) + codegen_child_reuse_bis_child2(i)
        return
    end
    asm = sprint(io->CUDAnative.code_ptx(io, codegen_child_reuse_bis_parent1, Tuple{Int}))

    @eval function codegen_child_reuse_bis_parent2(i)
        codegen_child_reuse_bis_child1(i+1) + codegen_child_reuse_bis_child2(i+1)
        return
    end
    asm = sprint(io->CUDAnative.code_ptx(io, codegen_child_reuse_bis_parent2, Tuple{Int}))
end

@testset "indirect sysimg function use" begin
    # issue #9: re-using sysimg functions should force recompilation
    #           (host fldmod1->mod1 throws, so the PTX code shouldn't contain a throw)

    # FIXME: Int64 because of #49

    @eval function codegen_recompile(out)
        wid, lane = fldmod1(unsafe_load(out), Int64(32))
        unsafe_store!(out, wid)
        return
    end

    asm = sprint(io->CUDAnative.code_ptx(io, codegen_recompile, Tuple{Ptr{Int64}}))
    @test !occursin("jl_throw", asm)
    @test !occursin("jl_invoke", asm)   # forced recompilation should still not invoke
end

@testset "compile for host after PTX" begin
    # issue #11: re-using host functions after PTX compilation
    @eval @noinline codegen_recompile_bis_child(i) = sink(i+1)

    @eval function codegen_recompile_bis_fromhost()
        codegen_recompile_bis_child(10)
    end

    @eval function codegen_recompile_bis_fromptx()
        codegen_recompile_bis_child(10)
        return
    end

    CUDAnative.code_ptx(devnull, codegen_recompile_bis_fromptx, Tuple{})
    @test codegen_recompile_bis_fromhost() == 11
end

@testset "LLVM intrinsics" begin
    # issue #13 (a): cannot select trunc
    @eval function codegen_issue_13(x)
        unsafe_trunc(Int, x)
        return
    end
    CUDAnative.code_ptx(devnull, codegen_issue_13, Tuple{Float64})
end

@testset "kernel names" begin
    @eval codegen_regular() = nothing
    @eval codegen_closure = ()->nothing

    function test_name(f, name; kwargs...)
        code = sprint(io->CUDAnative.code_ptx(io, f, Tuple{}; kwargs...))
        @test occursin(name, code)
    end

    test_name(codegen_regular, "julia_codegen_regular")
    test_name(codegen_regular, "julia_codegen_renamed"; alias="codegen_renamed")

    test_name(codegen_regular, "ptxcall_codegen_regular"; kernel=true)
    test_name(codegen_regular, "ptxcall_codegen_renamed"; kernel=true, alias="codegen_renamed")

    test_name(codegen_closure, "julia_anonymous")
    test_name(codegen_closure, "julia_codegen_renamed"; alias="codegen_renamed")

    test_name(codegen_closure, "ptxcall_anonymous"; kernel=true)
    test_name(codegen_closure, "ptxcall_codegen_renamed"; kernel=true, alias="codegen_renamed")
end

end

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

end


text:
md"# BOLD\_055\_500kHZ\_732nm\_100us"

md"## Read waveform"

md" - Parameters:
 - laser frequency= $f500kHz
 - DAQ window= $w100mus
 - Laser pulses every $tlaser500kHz
 - sampling rate= $(round(sr73n5k1m, sigdigits=2)) GHz
 - sampling period = $(round(sp73n5k1m, sigdigits=2) )ns
"

fn73n5k1m = 1  # test file

begin
	fls73n5k1m = glob("../labdata/BOLD_055_500kHZ_750nm_100us/C1*.csv")
	n73n5k1m = length(fls73n5k1m)

	ffl73n5k1m = 200.0e+6
	wvfm73n5k1m = read_waveform(fls73n5k1m,fn73n5k1m; wmus=wmus100mus, wsym=true)
	first(wvfm73n5k1m,3)
end
begin
	sr73n5k1m = sampling_rate(wvfm73n5k1m,wmus100mus)/GHz
	sp73n5k1m = sampling_period(wvfm73n5k1m,wmus100mus)/ns
end
md"
Unfiltered wvfm:
- mean = $(round(mean73n5k1m, sigdigits=4))
- std = $(round(std73n5k1m, sigdigits=4))
"
begin
	mean73n5k1m, std73n5k1m, thrp73n5k1m, thrn73n5k1m = thrx(wvfm73n5k1m; nsigma=3.0, waveform="raw");
	pfw73n5k1m = plot_filtered_waveform(wvfm73n5k1m; window=false, waveform="raw", thrp=thrp73n5k1m, thrn=thrn73n5k1m, scatter=false)
	pfft73n5k1m =plot_fft(wvfm73n5k1m,wmus100mus)
	plot(pfw73n5k1m, pfft73n5k1m)
end

md"### Filter signal
- Filter frequency (in Hz) = $ffl73n5k1m

Filtered wvfm:
- mean = $fmean73n5k1m
- std = $fstd73n5k1m
"
begin
    filter_signal_lp!(wvfm73n5k1m, wmus100mus; filtertype="Butterworth", flhz=ffl73n5k1m, n=4)
    first(wvfm73n5k1m,3)
end

begin
    fmean73n5k1m, fstd73n5k1m, fthrp73n5k1m, fthr73n5k1m = thrx(wvfm73n5k1m; nsigma=2.0, waveform="filtered")
    plot_filtered_waveform(wvfm73n5k1m; window=false, thrp=fthrp73n5k1m, thrn=fthr73n5k1m)
end

md"### Select peaks"

sprom73n5k1m = 2.0  #prominence cut

md"
- Peak Prominence cut = $sprom73n5k1m
"
begin
	r73n5k1m = select_filtered_peaks(wvfm73n5k1m, fthrp73n5k1m; promsel = sprom73n5k1m, wsel=0.0);
	if r73n5k1m != Nothing
		wvflt73n5k1m, spk73n5k1m = r73n5k1m
		plot_filtered_peaks(wvflt73n5k1m, spk73n5k1m)
	end
end

md"### Optimizing prominence cut"

begin
	P73n5k1m, W73n5k1m = proms(fls73n5k1m, 1, n73n5k1m;
	                 wmus=wmus100mus, wsym=true, flhz=ffl73n5k1m, promsel=0.5)
	hP73n5k1m = lfi.LaserLab.hist1d(P73n5k1m, 20, 0.0, 10.0)
	hW73n5k1m = lfi.LaserLab.hist1d(W73n5k1m, 20, 5.0, 15.0)
end
lfi.LaserLab.hist1d(hP73n5k1m, "Prominence (promsel=0.5)")
lfi.LaserLab.hist1d(hW73n5k1m, "Width")

md"### Selection"
spks73n5k1m = select_peaks(fls73n5k1m, 1, n73n5k1m; wmus=wmus100mus, wsym=true, flhz=ffl73n5k1m, promsel=sprom73n5k1m)
t73n5k1m = reduce(vcat,[subtract_laser(spks73n5k1m[i].xs, tlus2) for i in 1:length(spks73n5k1m)])
hp73n5k1m = lfi.LaserLab.hist1d(t73n5k1m, 20, 0.0, 2.0)
lfi.LaserLab.hist1d(hp73n5k1m, "times")

md"### Fits"
begin
	tfit73n5k1m = lfi.LaserLab.centers(hp73n5k1m)
	vfit73n5k1m = hp73n5k1m.weights
	fit73n5k1m = curve_fit(mexp, tfit73n5k1m, vfit73n5k1m, pa0)
	cofe73n5k1m = coef(fit73n5k1m)
	stder73n5k1m = stderror(fit73n5k1m)
	tft73n5k1m = expo.(tfit73n5k1m, coef(fit73n5k1m)...)
	println("")
end
begin
    ps73n5k1m = scatter(tfit73n5k1m, vfit73n5k1m, yerr=sqrt.(vfit73n5k1m), markersize=2,
            color = :black,
            label="data",
            fmt = :png)
    pp73n5k1m = plot(ps73n5k1m, tfit73n5k1m, tft73n5k1m, lw=2, label="μ = $(round(cofe73n5k1m[2]*1000, sigdigits=2)) ns", fmt = :png)
    xlabel!("t (μs)")
    ylabel!("frequency")
    end

max_stars_repo_name: jjgomezcadenas/laserlab

max_issues_repo_name: jjgomezcadenas/laserlab

max_forks_repo_name: jjgomezcadenas/laserlab


text:
# This file is a part of TextSearch.jl

#####
using CategoricalArrays
export EntropyWeighting

"""
    EntropyWeighting(; smooth=0.0, lowerweight=0.0, weights=:balance)

Entropy weighting uses the empirical entropy of the vocabulary along classes to produce a notion of importance for each token
"""
struct EntropyWeighting <: GlobalWeighting end

function entropy(dist)
    e = 0.0
    ipop = 1/sum(dist)

    for x in dist
        p = x * ipop

        if p > 0.0
            e -= p * log2(p)
        end
    end

    e
end

"""
    VectorModel(ent::EntropyWeighting, lw::LocalWeighting, corpus::BOW, labels;
        mindocs::Integer=1,
        smooth::Float64=0.0,
        weights=:balance,
        lowerweight=0.0
    )

Creates a vector model using the input corpus.
"""
function VectorModel(ent::EntropyWeighting, lw::LocalWeighting, corpus::AbstractVector{BOW}, labels=nothing;
        mindocs=1,
        smooth::Float64=0.0,
        weights=:balance
    )
    labels === nothing && error("EntropyWeighting requires labels as a categorical array to work")
    nclasses = length(levels(labels))
    D = Dict{UInt64,Vector{Float64}}()
    V = Vocabulary()

    for i in eachindex(corpus)
        vec = corpus[i]
        code = labels.refs[i]
        #M = sum(x for x in values(vec))
        for (t, occs) in vec
            s = get(V, t, UnknownTokenStats)
            V[t] = TokenStats(s.occs + occs, s.ndocs + 1, 0f0)
            dist = get(D, t, nothing)
            if dist === nothing
                D[t] = dist = fill(smooth, nclasses)
            end

            dist[code] += 1  # occs/M # log2(1 + occs)
        end
    end

    weights = weights isa String ? Symbol(weights) : weights

    if weights === :balance
        weights = sum(Vector{Float64}, values(D))
        weights .= sum(weights) ./ weights
    elseif weights === :none
        weights = ones(Float64, nclasses)
    end

    tokens, maxfreq = _create_vocabulary_with_entropy(V, D, weights, nclasses, mindocs)

    VectorModel(ent, lw, tokens, maxfreq, length(tokens), length(corpus))
end

function _create_vocabulary_with_entropy(V, D, weights, nclasses, mindocs)
    tokens = Vocabulary()
    tokenID = 0
    maxfreq = 0
    maxent = log2(nclasses)

    for (t, s) in V
        s.ndocs < mindocs && continue
        dist = D[t]
        dist .= dist .* weights
        e = 1.0 - entropy(dist) / maxent
        tokens[t] = TokenStats(s.occs, s.ndocs, e)
        maxfreq = max(maxfreq, s.occs)
    end

    tokens, maxfreq
end

global_weighting(::VectorModel{EntropyWeighting}, s::TokenStats) = s.weight


text:
#md# ### Words
using WordCloud
stwords = ["us"];
words_weights = processtext(open(pkgdir(WordCloud) * "/res/Barack Obama's First Inaugural Address.txt"), stopwords=WordCloud.stopwords_en ∪ stwords)
words_weights = Dict(zip(words_weights...))
#md# ### Embedding
#md# The positions of words can be initialized with pre-trained word vectors so that similar words will appear near each other.
using Embeddings
using TSne
const embtable = load_embeddings(GloVe{:en})
const get_word_index = Dict(word => ii for (ii, word) in enumerate(embtable.vocab))
function get_embedding(word)
    ind = get_word_index[word]
    emb = embtable.embeddings[:,ind]
    return emb
end
wordvec = Dict()
for k in keys(words_weights)
    if k in keys(get_word_index)
        wordvec[k] = get_embedding(k)
    elseif lowercase(k) in keys(get_word_index)
        wordvec[k] = get_embedding(lowercase(k))
    else
        pop!(words_weights, k)
        println("remove ", k)
    end
end
words = keys(wordvec) |> collect
vectors = hcat(values(wordvec)...)
embedded = tsne(vectors', 2)
#md#
wc = wordcloud(
    words_weights,
    mask=box,
    masksize=(1000, 1000),
    cornerradius=0,
    density=0.3,
    colors=0.3,
    backgroundcolor=:maskcolor,
    state=initwords!,
    # angles = (0,45), font = "Helvetica thin", maskcolor=0.98,
)

pos = embedded
mean = sum(pos, dims=1) / size(pos, 1)
r = maximum(sqrt.(pos[:,1].^2 + pos[:,2].^2))
pos = (pos .- mean) ./ 2r
sz = collect(reverse(size(wc.mask)))'
sz0 = collect(getparameter(wc, :masksize)[1:2])'
pos = round.(Int, pos .* sz0 .+ sz ./ 2)

setpositions!(wc, words, eachrow(pos), type=setcenter!)
setstate!(wc, :placewords!)
generate!(wc, reposition=false)
paint(wc, "semantic_embedding.png")
#md# ![](semantic_embedding.png)
#md# ### Clustering
#md# Words can be further colored according to semantic clustering
using Clustering
V = vectors
G = V' * V
H = sum(V.^2, dims=1)
D = max.(0, (H .+ H' .- 2G))
D ./= sum(D) / length(D)
D .= .√D # the distance matrix
tree = hclust(D, linkage=:ward)
lb = cutree(tree, h=2, k=10)
println("$(length(lb)) words are divided into $(length(unique(lb))) groups")
#md#
colors = parsecolor(:seaborn_dark)
setcolors!(wc, words, colors[lb .% length(colors) .+ 1])
recolor!(wc, style=:reset)
paint(wc, "semantic_clustering.png")
#md# ![](semantic_clustering.png)
wc
#eval# runexample(:semantic)


text:
# FIXME show methods here are terrible
function axes_string(a)
    return join(string.(collect(axes(a))), "×")
end

Base.show(io::IO, x::IndicesArray) = show(io, parent(x))
function Base.show(io::IO, m::MIME{Symbol("text/plain")}, x::IndicesArray)
    print(io, "IndicesArray - $(axes_string(x))\n")
    show(io, m, parent(x))
end


text:
# ---------------------------------------------------------------------------- #
#
#   mesh.jl
#
#   Abstract mesh type
#   This allows for writing an n-dimensional version of solver methods
#
#   λυτέος
#   Fall 2017
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

"""
    Mesh

Mesh abstract type:
Overarching abstract type for mesh types (2D and 3D).
"""
abstract type Mesh

end

include("mesh2D.jl")
include("mesh3D.jl")


text:
module NiBundleAdjustment

using NiLang, NiLang.AD
using LinearAlgebra
using ForwardDiff
using StaticArrays
using Requires

export compute_ba_J, Camera, read_ba_instance

struct Camera{T,V3,V2}
    r::V3
    c::V3
    f::T
    x0::V2
    κ::V2
end

include("forwarddiff.jl")
include("reversible.jl")
include("load.jl")

function __init__()
    @require KernelAbstractions="63c18a36-062a-441e-b654-da1e3ab1ce7c" include("cuda.jl")
end

end # module


text:
@safetestset "Dynamics Diff Tests" begin
    include("dynamics_test_functions.jl")
    include("dynamics_test.jl")
end

@safetestset "Joint Diff Tests" begin
    include("joint_test_functions.jl")
    include("joint_test.jl")
end

@safetestset "Joint Diff Tensor Tests" begin
    include("tensor_test.jl")
end


text:
using Plots

include("helper_functions.jl")

begin
    λ1 = 1
    λ2 = 1
    α = 1
    μ(α) = 1 / (1.5 * (α))
    ν(α) = 1 / (0.7 * (α))
    qs = 0:0.01:1
    parameters = (λ1, λ2, α, μ, ν)
end

ps = numerical_pstar.(qs, parameters...)
a_ps = pstar.(qs, parameters...)
costs = cost.(ps, qs, parameters...)
a_costs = cost.(a_ps, qs, parameters...)

ν_assumption = ν.(α) .> λ1 .* ps .+ λ2 .* qs
μ_assumption = μ.(α) .> λ1 .* (1 .- ps)

a_ν_assumption = ν.(α) .> λ1 .* a_ps .+ λ2 .* qs
a_μ_assumption = μ.(α) .> λ1 .* (1 .- a_ps)

plot(qs, ps, label = "autodiff")
plot!(qs, a_ps, label = "analytic")

# plot!(qs, ν_assumption, colour = :red, label = "pu o ν")
# plot!(qs, μ_assumption, colour = :black, label = "pr o μ")
# plot!(qs, a_ν_assumption, colour = :red, label = "pu o ν")
# plot!(qs, a_μ_assumption, colour = :black, label = "pr o μ")

plot(qs, costs, label = "autodiff_costs")
plot!(qs, a_costs, label = "analytic_costs")
plot!(ylims = (-50, 50))

p = plot(qs, qs, (p, q) -> clamp(cost(p, q, parameters...), -50, 50), st = :surface, zlims =  (-50, 50))

# begin
#     spin_plot = @animate for i in vcat(0:0.5:90, reverse(0:0.5:90))
#         plot(p, camera = (i,i))
#     end
#     gif(spin_plot, "spin_plot.gif", fps = 60)
# end

# gif(spin_plot, "spin_plot.gif", fps = 60)

plot!(shapes_from_intervals(shaded_regions(μ_assumption, qs)), colour=:red, label = "μ violated", alpha=0.3)
plot!(shapes_from_intervals(shaded_regions(ν_assumption, qs)), colour=:blue, label = "ν violated", alpha=0.3)
plot!(shapes_from_intervals(shaded_regions(a_μ_assumption, qs)), colour=:yellow, label = "a_μ violated", alpha=0.3)
plot!(shapes_from_intervals(shaded_regions(a_ν_assumption, qs)), colour=:gray, label = "a_ν violated", alpha=0.3)

plot!(shapes_from_intervals(push!(shaded_regions(μ_assumption, qs), [0, 0.1])), colour=:red, label = "μ violated", alpha=0.3)


text:
include("help.jl")
function perimetr!(r::Robot)
    num_vert = moves!(r, Sud)
    num_hor = moves!(r, West)

    for side in (Nord, Ost, Sud, West)
        putmarkers!(r, side)
    end

    moves!(r, Nord, num_vert)
    moves!(r, Ost, num_hor)
end

max_stars_repo_name: MCPavlo/mirea-progs

max_issues_repo_name: MCPavlo/mirea-progs

max_forks_repo_name: MCPavlo/mirea-progs


text:
# p27.jl - Solve KdV eq. u_t + uu_x + u_xxx = 0 on [-pi,pi] by
#         FFT with integrating factor v = exp(-ik^3t)*u-hat.

# Set up grid and two-soliton initial data:
N = 256; dt = .4/N^2; x = (2*pi/N)*(-N/2:N/2-1);
A = 25; B = 16; clf();
u = @. 3*A^2*sech(.5*(A*(x+2)))^2 + 3*B^2*sech(.5*(B*(x+1)))^2;
v = fft(u); k = [0:N/2-1;0;-N/2+1:-1]; ik3 = 1im*k.^3;

# Solve PDE and plot results:
tmax = 0.006; nplt = floor(Int,(tmax/25)/dt); nmax = round(Int,tmax/dt);
udata = u; tdata = [0.0];
for n = 1:nmax
    t = n*dt; g = -.5im*dt*k;
    E = exp.(dt*ik3/2); E2 = E.^2;
    a = g.*fft(real( ifft(     v    ) ).^2);
    b = g.*fft(real( ifft(E.*(v+a/2)) ).^2);     # 4th-order
    c = g.*fft(real( ifft(E.*v + b/2) ).^2);     # Runge-Kutta
    d = g.*fft(real( ifft(E2.*v+E.*c) ).^2);
    v = E2.*v + (E2.*a + 2*E.*(b+c) + d)/6;
    if mod(n,nplt) == 0
        u = real(ifft(v));
        udata = [udata u]; tdata = [tdata;t];
    end
end
mesh(x,tdata,udata',ccount=0,rcount=N); view(-20,25);
xlabel("x"); ylabel("y"); grid(true);
xlim(-pi,pi); ylim(0,tmax); zlim(0,12000);|
zticks([0,2000]);


text:
#
# errors.jl -
#
# Handling of errors for Julia interface to ActiveSilicon Phoenix (PHX)
# library.
#
#------------------------------------------------------------------------------
#
# This file is part of the `Phoenix.jl` package which is licensed under the MIT
# "Expat" License.
#
# Copyright (C) 2017-2019, Éric Thiébaut (https://github.com/emmt/Phoenix.jl).
# Copyright (C) 2016, Éric Thiébaut & Jonathan Léger.
#

Base.showerror(io::IO, e::PHXError) =
    print(io, "PHXError: ", geterrormessage(e.status), " [",
          geterrorsymbol(e.status), "]")

"""
   geterrormessage(code)

yields the error message corresponding to `code`.

See also: [`geterrorsymbol`](@ref), [`printerror`](@ref).

"""
function geterrormessage(code::Integer)
    buf = zeros(UInt8, 512)
    ccall(_PHX_ErrCodeDecode[], Nothing, (Ptr{UInt8}, Status), buf, code)
    return unsafe_string(pointer(buf))
end

"""
   geterrorsymbol(code)

yields the symbol corresponding to status value `code`.

See also: [`geterrormessage`](@ref), [`printerror`](@ref).

"""
geterrorsymbol(code::Integer) = get(_ERRSYM, code, :PHX_UNKNOWN_STATUS)

const _ERRSYM = Dict{Status, Symbol}()
for s in (:PHX_OK, :PHX_ERROR_BAD_HANDLE, :PHX_ERROR_BAD_PARAM,
          :PHX_ERROR_BAD_PARAM_VALUE, :PHX_ERROR_READ_ONLY_PARAM,
          :PHX_ERROR_OPEN_FAILED, :PHX_ERROR_INCOMPATIBLE,
          :PHX_ERROR_HANDSHAKE, :PHX_ERROR_INTERNAL_ERROR,
          :PHX_ERROR_OVERFLOW, :PHX_ERROR_NOT_IMPLEMENTED,
          :PHX_ERROR_HW_PROBLEM, :PHX_ERROR_NOT_SUPPORTED,
          :PHX_ERROR_OUT_OF_RANGE, :PHX_ERROR_MALLOC_FAILED,
          :PHX_ERROR_SYSTEM_CALL_FAILED, :PHX_ERROR_FILE_OPEN_FAILED,
          :PHX_ERROR_FILE_CLOSE_FAILED, :PHX_ERROR_FILE_INVALID,
          :PHX_ERROR_BAD_MEMBER, :PHX_ERROR_HW_NOT_CONFIGURED,
          :PHX_ERROR_INVALID_FLASH_PROPERTIES, :PHX_ERROR_ACQUISITION_STARTED,
          :PHX_ERROR_INVALID_POINTER, :PHX_ERROR_LIB_INCOMPATIBLE,
          :PHX_ERROR_SLAVE_MODE, :PHX_ERROR_DISPLAY_CREATE_FAILED,
          :PHX_ERROR_DISPLAY_DESTROY_FAILED, :PHX_ERROR_DDRAW_INIT_FAILED,
          :PHX_ERROR_DISPLAY_BUFF_CREATE_FAILED,
          :PHX_ERROR_DISPLAY_BUFF_DESTROY_FAILED,
          :PHX_ERROR_DDRAW_OPERATION_FAILED, :PHX_ERROR_WIN32_REGISTRY_ERROR,
          :PHX_ERROR_PROTOCOL_FAILURE,
          :PHX_WARNING_TIMEOUT, :PHX_WARNING_FLASH_RECONFIG,
          :PHX_WARNING_ZBT_RECONFIG, :PHX_WARNING_NOT_PHX_COM,
          :PHX_WARNING_NO_PHX_BOARD_REGISTERED, :PHX_WARNING_TIMEOUT_EXTENDED)
    _ERRSYM[@eval($s)] = s
end

"""

    printerror() -> curmode

yields the current error mode which indicates whether Phoenix error messages
are printed by the default error handler.  To choose the behavior, call:

    printerror(newmode) -> oldmode

which set the error mode to be `newmode` and yields the previous setting.

See also: [`geterrormessage`](@ref).

"""
printerror() = _PRNTERRMODE_REF[]
printerror(newmode::Bool) =
    (oldmode = _PRNTERRMODE_REF[]; _PRNTERRMODE_REF[] = newmode; oldmode)

# `_PRNTERRMODE_REF` is a global reference (hence non-volatile) to store whether or
# not print error messages.  In case of error, if `_PRNTERRMODE_REF[]` is true, the
# error handler `_errorhandler` will call Phoenix default error handler (to
# avoid using any Julia i/o routines); otherwise nothing is printed.  This
# mecanism is intended to be able to toggle printing of error messages while
# being thread safe because nothing from Julia is used.
#
const _PRNTERRMODE_REF = Ref{Bool}(true)
const _LASTERRFUNC_REF = Ref{Ptr{Cchar}}(C_NULL)
const _LASTERRCODE_REF = Ref{Status}(PHX_OK)
const _LASTERRMESG_REF = Ref{Ptr{Cchar}}(C_NULL)

"""
    _errorhandler(func, code, mesg)

error handler for the Phoenix library.  Temporarily stores last error and, if
`_PRNTERRMODE_REF[]` is true, immediately calls `_printlasterror()`; otherwise the
error message can be printed a bit later (however prior to calling any Phoenix
functions) by calling `_printlasterror()`.  This method is thread safe (it does
not use Julia engine) as you can check with:

    code_native(Phoenix._errorhandler, (Ptr{Cchar}, Cint, Ptr{Cchar}))

"""
function _errorhandler(func::Ptr{Cchar}, code::Status, mesg::Ptr{Cchar})
    _LASTERRFUNC_REF[] = func
    _LASTERRCODE_REF[] = code
    _LASTERRMESG_REF[] = mesg
    if _PRNTERRMODE_REF[]
        _printlasterror()
    end
end

"""
    _printlasterror()

prints last error (if and only if `_PRNTERRMODE_REF[]` is true) using the
default error handler of the Phoenix library and clear memorized error.  This
method is intended to be called right after an error occured.  This method is
thread safe (it does not use Julia engine) as you can check with:

    code_native(Phoenix._printlasterror, ())

"""
function _printlasterror()
    if _PRNTERRMODE_REF[] && _LASTERRFUNC_REF[] != C_NULL && _LASTERRMESG_REF[] != C_NULL
        ccall(_PHX_ErrHandlerDefault[], Cvoid, (Ptr{Cchar}, Status, Ptr{Cchar}),
              _LASTERRFUNC_REF[], _LASTERRCODE_REF[], _LASTERRMESG_REF[])
    end
    _LASTERRFUNC_REF[] = C_NULL
    _LASTERRCODE_REF[] = PHX_OK
    _LASTERRMESG_REF[] = C_NULL
    return nothing
end


text:
using Parameters

@with_kw struct Stokes_Para
    # mesh
    Nintervals::Int = 32        # number of intervals for your mesh
    xmin::Float64 = -2.0        # minimum x-value for mesh
    xmax::Float64 = 2.0         # minimum x-value for mesh
    ymin::Float64 = -2.0        # minimum x-value for mesh
    ymax::Float64 = 1.0         # minimum x-value for mesh

    # parameters
    μ::Float64 = 1.0            # viscosity

    # boundary conditions
    lBCval::Float64 = 1.0       # strength of inflow

    # visualization
    solFolder::String = "solutions"         # folder for solution to go in
    solFile::String = "stokes_example"     # name of vtk solutions output
end

# test comment

max_stars_repo_name: pseastham/eFEMpart

max_issues_repo_name: pseastham/eFEMpart


text:
export RedisString

struct RedisString <: AbstractRedisCollection
    conn::AbstractRedisConnection
    key::String
end

RedisString(key) = RedisString(default_connection, key)

function getindex(rs::RedisString, x::Integer, y::Integer)
    # it is always not null: non-exist keys are treated as empty string by redis
    exec(rs.conn, "getrange", rs.key, zero_index(x), zero_index(y)) |> String
end

function getindex(rs::RedisString, x::UnitRange{T}) where T<:Integer
    rs[x.start, x.stop]
end

function getindex(rs::RedisString, x::Integer) # do we really need this?
    rs[x, x][1]
end

function getindex(rs::RedisString, ::Colon=Colon())
    String(@some(exec(rs.conn, "get", rs.key)))
end

function setindex!(rs::RedisString, value, offset::Integer)
    exec(rs.conn, "setrange", rs.key, zero_index(offset), value)
    rs
end

function setindex!(rs::RedisString, value, x::UnitRange{T}) where T<:Integer
    value = value |> string |> String
    x.stop - x.start + 1 == sizeof(value) || throw(DimensionMismatch())
    rs[x.start] = value
    rs
end

function setindex!(rs::RedisString, value, ::Colon=Colon())
    exec(rs.conn, "set", rs.key, value)
    rs
end

"this methods should always be used in the form of `rs += xx`"
function (+)(rs::RedisString, value)
    exec(rs.conn, "append", rs.key, value)
    rs
end

"NOTE: Redis count length by byte, while julia count by char"
function length(rs::RedisString)
    exec(rs.conn, "strlen", rs.key)::Int64
end

function lastindex(rs::RedisString)
    length(rs)
end

function string(rs::RedisString)
    rs[:]
end


text:
"""
    simulate_continuous(i, t, η)

Simulate a simple (linear) birth and death process starting at time zero with
`i` units. The stochastic process is observed continuously for `t` time units.
Rates are given by parameter `η = [λ, μ]` where `λ` is the birth rate and `μ` is
the death rate.
"""
function simulate_continuous(
  i::I,
  t::F,
  η::Vector{F}
) where {
  I <: Integer,
  F <: AbstractFloat
}
  ω = η[1] + η[2]
  p = η[1] / ω

  if t < floatmin(F)
    return (zero(Int), zero(Int), zero(F), zero(Float64), zero(Int), [], Int(i),
            [])
  end

  tot_births = 0
  tot_deaths = 0
  integrated_jump = zero(F)
  sum_log_n = zero(Float64)

  # pre-allocate vectors. If needed, we will grow the vector by doubling its
  # size each time.
  len = 100
  waiting_time = zeros(F, len)
  increment = zeros(Int, len)

  popul_size = i
  h = 0
  total_time = zero(F)
  while true
    # sample next event time from a negative exponential distribution with rate
    # n * (λ + μ)
    τ = - log(rand(F)) / (popul_size * ω)
    total_time += τ

    if total_time > t
      integrated_jump += popul_size * (t - total_time + τ)
      break
    end

    if h == len - 1
      # not enough space to store this event
      len *= 2
      resize!(waiting_time, len)
      resize!(increment, len)
    end

    h += 1

    integrated_jump += popul_size * τ
    sum_log_n += log(popul_size)
    waiting_time[h] = τ

    # sample event type with probability (λ / (λ + μ), μ / (λ + μ))
    if rand(F) <= p
      tot_births += 1
      increment[h] = 1
      popul_size += 1
    else
      tot_deaths += 1
      increment[h] = -1
      popul_size -= 1

      if popul_size == 0
        break
      end
    end
  end

  resize!(waiting_time, h)
  resize!(increment, h)

  (tot_births, tot_deaths, integrated_jump, sum_log_n, h, waiting_time, Int(i),
   increment)
end

"""
    simulate_discrete(i, k, u, η)

Simulate a simple (linear) process starting at time zero with `i` units. The
stochastic process is observed every `u` time units until it reaches `u * k`
time. Rates are given by parameter `η = [λ, μ]` where `λ` is the birth rate and
`μ` is the death rate.
"""
function simulate_discrete(
  i::I,
  k::I,
  u::F,
  η::Vector{F}
)::Vector{I} where {
  I <: Integer,
  F <: AbstractFloat
}
  if k == 0
    return [i]
  end

  ω = η[1] + η[2]
  p = η[1] / ω

  event_time = zero(F)
  popul_size = i

  sim_size = zeros(I, k + 1)
  sim_size[1] = popul_size

  counter = 1
  t = u
  while true
    # sample event time from a negative exponential distribution with rate
    # n * (λ + μ)
    event_time -= log(rand(F)) / F(popul_size * ω)

    # new event might be after various null observations
    while (event_time > t) && (counter <= k)
      counter += 1
      t = counter * u
      sim_size[counter] = popul_size
    end

    # sample event type with probability (λ / (λ + μ), μ / (λ + μ))
    popul_size += (rand(F) <= p) ? I(1) : -I(1)

    if (counter > k) || (popul_size == 0)
      break
    end
  end

  sim_size
end

"""
    rand_continuous(n, i, t, η)

Simulate `n` simple (linear) birth and death processes with the same rate
parameter `η = [λ, μ]`, where `λ` is the birth rate and `μ` is the death rate.
Starting at time `0` with `i` units, observe continuously each process until
time `t` is reached.

For simulating only one process, use `rand_continuous(i, t, η)`.
"""
function rand_continuous(
  n::Integer,
  i::Integer,
  t::Real,
  η::Vector{T}
) where {
  T <: Real
}
  if n < 1
    error("Number of simulations 'n' is less than or equal to zero.")
  end

  if i < 1
    error("Initial population size 'i' is less than or equal to zero.")
  end

  if t < 0
    error("Time 't' is negative.")
  end

  if η[1] < 0
    error("Birth rate 'η[1]' is negative.")
  end

  if η[2] < 0
    error("Death rate 'η[2]' is negative.")
  end

  if (η[1] < floatmin(float(η[1]))) && (η[2] < floatmin(float(η[2])))
    error("Birth rate 'η[1]' and death rate 'η[2]' are jointly zero.")
  end

  F = typeof(float((η[1] - η[2]) * t))

  [ObservationContinuousTime{F}(
   simulate_continuous(i, F(t), [F(η[1]), F(η[2])])...) for s = 1:n]
end

"""
    rand_continuous(i, t, η)

Simulate one simple (linear) birth and death process with rate parameter
`η = [λ, μ]`, where `λ` is the birth rate and `μ` is the death rate. Starting at
time `0` with `i` units, observe continuously the process until time `t` is
reached.

For simulating `n` processes, use `rand_continuous(n, i, t, η)`.
"""
function rand_continuous(
  i::Integer,
  t::Real,
  η::Vector{T}
) where {
  T <: Real
}

  if i < 1
    error("Initial population size 'i' is less than or equal to zero.")
  end

  if t < 0
    error("Time 't' is negative.")
  end

  if η[1] < 0
    error("Birth rate 'η[1]' is negative.")
  end

  if η[2] < 0
    error("Death rate 'η[2]' is negative.")
  end

  if (η[1] < floatmin(float(η[1]))) && (η[2] < floatmin(float(η[2])))
    error("Birth rate 'η[1]' and death rate 'η[2]' are jointly zero.")
  end

  F = typeof(float((η[1] - η[2]) * t))

  ObservationContinuousTime{F}(
  simulate_continuous(i, F(t), [F(η[1]), F(η[2])])...)
end

"""
    rand_discrete(n, i, k, u, η)

Simulate `n` simple (linear) birth and death processes with the same rate
parameter `η = [λ, μ]`, where `λ` is the birth rate and `μ` is the death rate.
Starting at time `0` with `i` units, observe each process for a total of `k`
times, each separated by the same lag `u`.

For simulating only one process, use `rand_discrete(i, k, u, η)`.
"""
function rand_discrete(
  n::Integer,
  i::Integer,
  k::Integer,
  u::Real,
  η::Vector{T}
) where {
  T <: Real
}
  if n < 1
    error("Number of simulations 'n' is less than or equal to zero.")
  end

  if i < 1
    error("Initial population size 'i' is less than or equal to zero.")
  end

  if k < 0
    error("Number of observations 'k' is negative.")
  end

  if u < floatmin(float(u))
    error("Time lag 'u' is less than or equal to zero.")
  end

  if η[1] < 0
    error("Birth rate 'η[1]' is negative.")
  end

  if η[2] < 0
    error("Death rate 'η[2]' is negative.")
  end

  if (η[1] < floatmin(float(η[1]))) && (η[2] < floatmin(float(η[2])))
    error("Birth rate 'η[1]' and death rate 'η[2]' are jointly zero.")
  end

  I = typeof(i)
  F = typeof(float((η[1] - η[2]) * u))

  state = zeros(I, k + 1, n)

  for s = 1:n
    state[:, s] = simulate_discrete(i, I(k), F(u), [F(η[1]), F(η[2])])
  end

  ObservationDiscreteTimeEven{F, I}(Int(n), Int(k), F(u), state)
end

"""
    rand_discrete(i, k, u, η)

Simulate one simple (linear) birth and death process with rate parameter
`η = [λ, μ]`, where `λ` is the birth rate and `μ` is the death rate. Starting at
time `0` with `i` units, observe each process for a total of `k` times, each
separated by the same lag `u`.

For simulating `n` processes, use `rand_discrete(n, i, k, u, η)`.
"""
function rand_discrete(
  i::Integer,
  k::Integer,
  u::Real,
  η::Vector{T}
) where {
  T <: Real
}
  rand_discrete(1, i, k, u, η)
end


text:
using Query
using NDSparseData

source_ndsparsearray1 = NDSparse([fill("New York",3); fill("Boston",3)],
                            repmat(Date(2016,7,6):Date(2016,7,8), 2),
                            [91,89,91,95,83,76])

q = @from i in source_ndsparsearray1 begin
    @where i.index[1]=="Boston"
    @select i.value
    @collect
end

println(q)

source_ndsparsearray2 = NDSparse(Columns(city = [fill("New York",3); fill("Boston",3)],
                            date = repmat(Date(2016,7,6):Date(2016,7,8), 2)),
                            [91,89,91,95,83,76])

q = @from i in source_ndsparsearray2 begin
    @where i.index.city=="New York"
    @select i.value
    @collect
end

println(q)


text:
################################################################################
######################### LowLevelSolution Constructors ########################
################################################################################
export
    get_initial_solution,
    get_infeasible_solution,
    default_solution

"""
    `get_initial_solution`
"""
function get_initial_solution(mapf::MAPF{E,S,G}) where {S,A,G,T,C<:AbstractCostModel{T},E<:AbstractLowLevelEnv{S,A,C}}
    LowLevelSolution{S,A,T,C}(
        paths = map(i->Path{S,A,T}(
                    s0=get_start(mapf,i),
                    cost=get_initial_cost(mapf.env)),
                1:num_agents(mapf)),
        cost_model = get_cost_model(mapf.env),
        costs = Vector{T}(map(a->get_initial_cost(mapf.env),1:num_agents(mapf))),
        cost = get_initial_cost(mapf.env),
    )
end
"""
    `get_infeasible_solution`
"""
function get_infeasible_solution(mapf::MAPF{E,S,G}) where {S,A,G,T,C<:AbstractCostModel{T},E<:AbstractLowLevelEnv{S,A,C}}
    LowLevelSolution{S,A,T,C}(
        paths = Vector{Path{S,A,T}}(map(a->Path{S,A,T}(cost=get_initial_cost(mapf.env)),1:num_agents(mapf))),
        cost_model = get_cost_model(mapf.env),
        costs = Vector{T}(map(a->get_initial_cost(mapf.env),1:num_agents(mapf))),
        cost = get_infeasible_cost(mapf.env),
    )
end
"""
    `default_solution(solver, mapf::AbstractMAPF)`

    Defines what is returned by the solver in case of failure to find a feasible
    solution.
"""
function default_solution(mapf::M) where {M<:AbstractMAPF}
    return get_infeasible_solution(mapf), get_infeasible_cost(mapf.env)
end

export add_to_path!

"""
    add_to_path!(path,env,s,a,sp)

Adds the new (s,a,sp) tuple and its cost (under env) to path.
"""
function add_to_path!(path::Path,env,s,a,sp)
    push!(path,PathNode(s,a,sp))
    set_cost!(path,accumulate_cost(env,get_cost(path),
        get_transition_cost(env,s,a,sp)))
    path
end

################################################################################
################################# Extend Paths #################################
################################################################################
export
    extend_path!

function extend_path!(env::E,path::P,T::Int) where {E<:AbstractLowLevelEnv,P<:Path}
    # while get_index_from_time(path,get_end_index(path)) < T
    while get_end_index(path) < T
        s = get_final_state(path)
        a = wait(env,s)
        sp = get_next_state(env,s,a)
        push!(path,PathNode(s,a,sp))
        set_cost!(path, accumulate_cost(env, get_cost(path), get_transition_cost(env,s,a,sp)))
    end
    return path
end
function extend_path!(env::E,path::P,T::Float64) where {E<:AbstractLowLevelEnv,P<:Path}
    @assert abs(T - Int(round(T))) < 0.01
    extend_path!(env,path,Int(round(T)))
end
function extend_path(env::E,path::P,args...) where {E<:AbstractLowLevelEnv,P<:Path}
    new_path = copy(path)
    extend_path!(new_path,args...)
    return new_path
end

get_transition_cost(env,n::PathNode) = get_transition_cost(env,get_s(n),get_a(n),get_sp(n))

export recompute_cost

"""
    recompute_cost

Recompute the cost of `path` (according to `env`), beginning from initial cost
`c0`.
"""
function recompute_cost(env,path,c0=get_initial_cost(env))
    cost = c0
    for n in path.path_nodes
        cost = accumulate_cost(env,cost,get_transition_cost(env,n))
    end
    return cost
end

export replace_cost_model

export trim_path!

"""
    trim_path!(env,path,T)

Modify `path` to terminate at time step `T`. If `length(path) < T`, the path
will be extended to time step `T`.
"""
function trim_path!(env,path,T)
    while length(path) > T
        # deleteat!(path.path_nodes,T)
        pop!(path.path_nodes)
    end
    if length(path) < T
        extend_path!(env,path,T)
    end
    set_cost!(path,recompute_cost(env,path))
    @assert length(path) == T
    return path
end

export trim_solution!

"""
    trim_solution!

Modify `solution` so that all paths terminate at time step `T`.
"""
function trim_solution!(env,solution,T)
    for (agent_id,path) in enumerate(get_paths(solution))
        trim_path!(env,path,T)
        set_path_cost!(solution,get_cost(path),agent_id)
    end
    set_cost!(solution,
        aggregate_costs(get_cost_model(env),get_path_costs(solution)))
    return solution
end

export sorted_actions

"""
    sorted_actions(env,s)

Return a vector of actions sorted lowest cost to highest cost.
"""
function sorted_actions(env,s)
    f = (s,a,sp)->compute_heuristic_cost(env,get_transition_cost(env,s,a,sp),sp)
    sort(
        collect(get_possible_actions(env,s)),
        by=a->f(s,a,get_next_state(env,s,a))
    )
end

################################################################################
################################## Utilities ###################################
################################################################################
export
    is_consistent,
    assert_valid,
    is_valid

"""
    is_consistent(solution,mapf)

    Check if solution satisfies start and end constraints
"""
function is_consistent(path::P,start::S,goal::G) where {S,G,P<:Path}
    return (states_match(get_initial_state(path), start)
        && states_match(get_final_state(path), goal))

end
function is_consistent(paths::Vector{P},starts::Vector{S},goals::Vector{G}) where {S,G,P <: Path}
    for (i,path) in enumerate(paths)
        if !is_consistent(path,starts[i],goals[i])
            return false
        end
    end
    return true
end
function is_consistent(solution::L,starts::Vector{S},goals::Vector{G}) where {S,G,L<:LowLevelSolution}
    is_consistent(get_paths(solution),starts,goals)
end
function is_consistent(solution::L,mapf::M) where {L<:LowLevelSolution, M<:AbstractMAPF}
    is_consistent(solution,get_starts(mapf),get_goals(mapf))
end
assert_valid(args...) = @assert(is_consistent(args...))
function is_valid(r)
    try
        assert_valid(r)
    catch e
        if isa(e,AssertionError)
           showerror(stdout,e)
           return false
        else
           rethrow(e)
        end
    end
    return true
end

################################################################################
######################### Visualization and Debugging ##########################
################################################################################

export
    convert_to_vertex_lists

function convert_to_vertex_lists(path) end
function convert_to_vertex_lists(solution::S) where {S<:LowLevelSolution}
    return [convert_to_vertex_lists(path) for path in get_paths(solution)]
end

export
    sprint_route_plan

function sprint_route_plan(io::IO,solution::LowLevelSolution;kwargs...)
    sprint_indexed_list_array(io::IO,convert_to_vertex_lists(solution);kwargs...)
end
function sprint_route_plan(route_plan::LowLevelSolution;kwargs...)
    buffer = IOBuffer()
    sprint_route_plan(buffer,route_plan;kwargs...)
    String(take!(buffer))
end
function Base.show(io::IO,route_plan::LowLevelSolution)
    # buffer = IOBuffer()
    # print(buffer,"LowLevelSolution:\n")
    print(io,"LowLevelSolution:\n")
    # sprint_route_plan(buffer,route_plan;leftaligned=true)
    sprint_route_plan(io,route_plan;leftaligned=true)
    # print(io,String(take!(buffer)))
end

export
    summarize

summarize(env) = ""


text:
tests = ["mxarray",
         "matfile",
         "mstatements",
         "engine"]

println("Testing MATLAB.jl")
for t in tests
    fp = joinpath("test", "$(t).jl")
    println("* running $fp ...")
    include(fp)
end


text:
module LiquidsStructureGrids

### Imports
using Reexport
@reexport using ApproximationGrids
@reexport using LiquidsStructure

import ApproximationGrids: ApproximationGrid

### Exports
export StructureFactorGrid

### Implemetation
struct StructureFactorGrid{S <: StructureFactor, G <: ApproximationGrid}
    f::S
    grid::G
end

function Base.show(io::IO, ::MIME"text/plain", s::StructureFactorGrid)
    print(io, "$(s.f),$(s.grid)")
end

end # module


text:
using IntervalArithmetic, IntervalRootFinding, StaticArrays
using Base.Test

@testset "1D roots" begin
    rts = roots(sin, -5..5)
    @test length(rts) == 4
    @test length(find(x->x==:unique, [root.status for root in rts])) == 2

    rts = roots(sin, -5..6, Bisection)
    @test length(rts) == 3
    rts = roots(sin, rts, Newton)
    @test all(root.status == :unique for root in rts)
end

@testset "2D roots" begin
    f(x, y) = SVector(x^2 + y^2 - 1, y - 2x)
    f(X) = f(X...)

    rts = roots(f, (-6..6) × (-6..6), Bisection, 1e-3)
    @test length(rts) == 4

    rts = roots(f, rts, Newton)
    @test length(rts) == 2
end

@testset "Complex roots" begin
    x = -5..6
    Xc = Complex(x, x)
    f(z) = z^3 - 1

    rts = roots(f, Xc, Bisection, 1e-3)
    @test length(rts) == 7
    rts = roots(f, rts, Newton)
    @test length(rts) == 3
    rts = roots(f, Xc)
    @test length(rts) == 3
end

# From R docs:
# https://www.rdocumentation.org/packages/pracma/versions/1.9.9/topics/broyden

@testset "3D roots" begin
    function g(x)
        (x1, x2, x3) = x
        SVector(    x1^2 + x2^2 + x3^2 - 1,
                    x1^2 + x3^2 - 0.25,
                    x1^2 + x2^2 - 4x3
                )
    end

    X = (-5..5)
    rts = roots(g, IntervalBox(X, 3))
    @test length(rts) == 4
end

@testset "Stationary points" begin
    f(xx) = ( (x, y) = xx; sin(x) * sin(y) )
    rts = roots(∇(f), IntervalBox(-5..6, 2), Newton, 1e-5)
    @test length(rts) == 25
end


text:
module Helper

export initialize, importData, config!, unpackStackedArray, allocMat!,
  unpackStackedArray!, collectResult, assignFilename, updateIVmatrix!,
  convert2share, initδ!, importθ̂₂
using TypeDef
using DelimitedFiles, Statistics, Printf, Random, StatsBase

function importData(
  model  ::Int64,   # ID for demand model to estimate
  iv     ::String,  # IV approach (opt,blp,diff-quad,diff-local)
  nPerson::Int64,   # Number of simulation draws for random coefficients
  ver    ::String,  # Version of estimation data (824 for final version)
  impFile::String;  # Type of import files (corresponding to model & IV)
  output ::Int64=1) # Type of function output

  workPath = "$(homedir())/work/kantar/brand/workfiles/"

  if (impFile=="base") | (impFile=="optIVtest") | (impFile=="optIV2")
    fileD = string(workPath,"demand",ver,".csv")
    fileX = string(workPath,"Xinput",ver,".csv")
    if iv=="blp"
      fileZ = string(workPath,"ZinputBlp",ver,".csv")
    else
      fileZ = string(workPath,"ZinputBlp",ver,"core.csv")
    end
    fileDIV  = string(workPath,"DiffIVinput",ver,"reduced.csv")
  elseif impFile=="extra"
    fileD = string(workPath,"demand",ver,"extra.csv")
    fileX = string(workPath,"Xinput",ver,"extra.csv")
    fileZ = string(workPath,"ZinputBlp",ver,"coreextra.csv")
    fileDIV  = string(workPath,"DiffIVinput",ver,"reduced2.csv")
  elseif impFile=="ms15"
    fileD = string(workPath,"demand",ver,"ms15.csv")
    fileX = string(workPath,"Xinput",ver,".csv")
    fileZ = string(workPath,"ZinputBlp",ver,"core.csv")
    fileDIV  = string(workPath,"DiffIVinput",ver,"reduced2.csv")
  elseif impFile=="noAllow"
    fileD = string(workPath,"demand",ver,"NoAllow.csv")
    fileX = string(workPath,"Xinput",ver,"NoAllow.csv")
    fileZ = string(workPath,"ZinputBlp",ver,"core.csv")
    fileDIV  = string(workPath,"DiffIVinput",ver,"reduced2.csv")
  end

  impD   = readdlm(fileD,',';header=true)
  impX   = readdlm(fileX,',';header=true)
  impZ   = readdlm(fileZ,',';header=true)
  impDIV = readdlm(fileDIV,',';header=true)

  qtr     = convert(Vector{Int64},impD[1][:,1])
  region  = convert(Vector{Int64},impD[1][:,2])
  network = convert(Vector{Int64},impD[1][:,3])
  brand   = convert(Vector{Int64},impD[1][:,4])
  pop     = convert(Vector{Float64},impD[1][:,5])
  ms1d    = convert(Vector{Float64},impD[1][:,6])
  msOO1d  = convert(Vector{Float64},impD[1][:,7])
  price   = convert(Vector{Float64},impD[1][:,8])
  # p̂       = convert(Vector{Float64},impD[1][:,9]) # Generated by setIV()
  rawInc  = convert(Vector{Float64},impD[1][:,10])
  sdrInc  = convert(Vector{Float64},impD[1][:,11])

  indata  = convert(BitArray{1},impD[1][:,12])
  prepost = convert(Vector{Int64},impD[1][:,13])
  idObs   = convert(Vector{Int64},impD[1][:,14])
  idLag   = convert(Vector{Int64},impD[1][:,15])

  logInc  = convert(Vector{Float64},impD[1][:,16])
  sdlInc  = convert(Vector{Float64},impD[1][:,17])

  Xinput  = impX[1]
  Zinput  = impZ[1]
  Dinput  = impDIV[1]
  product = sortslices(unique([brand prepost]; dims=1); dims=1)

  # ms = Vector{Vector{Float64}}(undef,0)
  dat = datType(region,qtr,network,brand,pop,ms1d,msOO1d,price,prepost)
  set = setType(unique(region),unique(qtr),unique(network),unique(brand),
                unique(prepost),product)

  # Model setup
  Xexo,Xᵢ,nXᵢ,nθ₂,nDᵢ,nνᵢ,nRᵢ,labΣᵢᵢ,labL,diagL,indxL,indxeθ,ixP =
    setModel(model,Xinput,dat)

  # Dimensions
  nMarket = size(unique([dat.region dat.qtr];dims=1),1)
  nθ₁ = size(Xexo,2) # Excludes price in X
  nProdt = Array{Float64,1}(undef,0)
  dim = dimType(length(set.region),length(set.qtr),length(set.network),
                length(set.brand),length(set.prepost),size(set.product,1),
                nPerson,nMarket,length(indata),count(indata.==true),
                nProdt,nθ₂,nθ₁,nθ₂+nθ₁,nXᵢ,nDᵢ,nνᵢ,nRᵢ)

  # Markt index
  mkt = BitArray(undef,dim.obs,0)
  for g in eachindex(set.region)
    for t in eachindex(set.qtr)
      mask = (dat.region.==set.region[g]).&(dat.qtr.==set.qtr[t])
      mkt = [mkt mask]
      push!(dim.prodt, sum(mask))
    end
  end
  if size(mkt,2)!=dim.market
    error("Mask dimension not matching unique markets.")
  end
  indx = indxType(mkt,diagL,indxL,idObs,idLag,indxeθ)

  ξ    = Matrix{Float64}(undef,dim.obs,1)
  Lξ   = Matrix{Float64}(undef,dim.obs,1)
  Eξ   = Matrix{Float64}(undef,dim.obs,1)
  ν    = zeros(nνᵢ,nPerson)
  y    = zeros(dim.obs,dim.person)
  dμdL = [zeros(Float64,dim.prodt[t],dim.θ₂,dim.person) for t=1:dim.market]
  dμ̂dL = [zeros(Float64,dim.prodt[t],dim.θ₂,dim.person) for t=1:dim.market]
  L    = Array{Float64,2}(undef,nXᵢ,nXᵢ)
  Lν   = zeros(Float64,dim.Xᵢ,dim.person)
  ∂δ̂∂θ = zeros(Float64,dim.obs,dim.θ₂)
  maskBrand = BitArray{2}(undef,dim.obs,dim.brand)
  for i=1:dim.brand
    @. maskBrand[:,i] = (~indata) & (dat.brand==i)
  end
  inc  = [zeros(Float64,2) for t=1:dim.market]
  linc = [zeros(Float64,2) for t=1:dim.market]
  Dᵢ   = [zeros(Float64,dim.Dᵢ,dim.prodt[t],dim.person) for t=1:dim.market]

  Z, X̂ᵢ, p̂ = setIV(Xᵢ,Xexo,Zinput,Dinput,price,ixP,dim,indx,iv)
  X = Xexo; # excludes price in X
  XZ = X'*Z; invPhi = inv(Z'*Z)
  XZPhiZXZ = (XZ*invPhi*XZ')\(XZ*invPhi*Z')
  #X̂ = Vector{Matrix{Float64}}(undef,0);
  X̂  = [zeros(Float64,dim.prodt[t],dim.θ₁) for t=1:dim.market]
  ms = [zeros(Float64,dim.prodt[t],) for t=1:dim.market]
  mat = matType(X,X̂,Xᵢ,X̂ᵢ,Xexo,Z,invPhi,XZ,XZPhiZXZ,p̂,ξ,Lξ,Eξ,ν,y,
                dμdL,dμ̂dL,L,Lν,∂δ̂∂θ,indata,maskBrand,inc,linc,Dᵢ,ms)

  for t=1:dim.market
    id = findfirst(view(indx.mkt,:,t).==true)
    mat.income[t]    .= [rawInc[id], sdrInc[id]]
    mat.logIncome[t] .= [logInc[id], sdlInc[id]]
  end

  # Labels
  labθ₁ = []
  if dim.Dᵢ==0
    labθ₁  = ["Price/income"]
  end
  for k in eachindex(impX[2])
    labθ₁ = [labθ₁; String(impX[2][k])]
  end
  labθ₁ = [labθ₁; "Constant"]
  lab = labType(labθ₁,labΣᵢᵢ,labL)
  if output==1
    return dat, dim, lab, mat, set, indx
  elseif output==0
    return dim, lab, indx
  end
end

function config!(
  model     ::Int64,
  dat       ::datType,
  dim       ::dimType,
  mat       ::matType,
  set       ::setType,
  indx      ::indxType,
  simEmpiric::Bool=false,
  lowBoundOn::Bool=false,
  truncate  ::Bool=false)

  workPath = "$(homedir())/work/kantar/brand/workfiles/simdraws/"
  fileName = string("lhsnorm",dim.Rᵢ,"d",dim.person,".csv")
  csvname  = string(workPath*fileName)
  impRand  = readdlm(csvname,',';header=false)
  randDraw = impRand[1:dim.person,:]'
  lowBound = -Inf64
  if lowBoundOn==true
    lowBound = -1.5615
  end
  z = max.(randDraw[1:1,:],lowBound)
  ν = randDraw[2:end,:]
  if (dim.Rᵢ < (1 + dim.νᵢ))
    ν = vcat(ν, repeat(ν[end:end,:], 1+dim.νᵢ-dim.Rᵢ,1))
    display("text/plain",ν[:,1:5])
    # exit()
  end
  println("Sim. draws imported from: ",fileName)

  # Import random draws for truncated log normal
  if truncate==true
    rng = MersenneTwister(999)
  end

  if simEmpiric==true
    impEmp = readdlm(string(workPath*"incomeEMdraws.csv"),',';header=true)
    region = convert(Vector{Int64},impEmp[1][:,1])
    income = convert(Vector{Float64},impEmp[1][:,2])/100
    yInv   = convert(Vector{Float64},impEmp[1][:,3])
  end

  # Variables of heterogeneous taste
  if dim.Dᵢ>0
    Dᵢ = zero(mat.y)

    if simEmpiric==false
      for t=1:dim.market
        if (model>=100) # Log normal income
          ymean = mat.logIncome[t][1]
          ystd  = mat.logIncome[t][2]
          ysim  = ymean .+ ystd*z[1:1,:]
          yᵣ = max.(exp.(ysim),300)/100
          if truncate==true
            for j=1:size(yᵣ,2)
              while yᵣ[j]>80
                ysim = ymean + ystd*randn(rng)
                yᵣ[1,j] = max(exp(ysim),300)/100
              end
            end
          end
        else            # Normal income
          yᵣ = max.(mat.income[t][1] .+ mat.income[t][2]*z[1:1,:],300)/100
          ### Test of lower cutoff
          # for j=1:size(yᵣ,2)
          #   while yᵣ[j]==3
          #     yᵣ[1,j]=max.(mat.income[t][1].+mat.income[t][2]*randn(rng),300)/100
          #   end
          # end
          ### end of test
        end
        # println("Market: ", t, " Minimum income: ",minimum(yᵣ))
        mask = view(indx.mkt,:,t)
        mat.y[mask,:] = repeat(yᵣ,dim.prodt[t],1)
        if (model==101) | (model==116)  # Non-inverse
          Dᵢ[mask,:] = mat.y[mask,:]/100
        else                            # Inverse function
          Dᵢ[mask,:] = 1 ./ mat.y[mask,:]
        end
        for k=1:dim.Dᵢ
          @. mat.Dᵢ[t][k,:,:] = Dᵢ[mask,:]
        end
      end
    elseif simEmpiric==true
      for r in set.region
        maskObs = findall(dat.region.==r)
        maskInc = findall(region.==r)
        mat.y[maskObs,:] = repeat(income[maskInc]',length(maskObs),1)
        Dᵢ[maskObs,:] = 1 ./ mat.y[maskObs,:]
      end
      for t=1:dim.market
        mask = view(indx.mkt,:,t)
        for k=1:dim.Dᵢ
          @. mat.Dᵢ[t][k,:,:] = Dᵢ[mask,:]
        end
      end
    end
  end
  # println("Mean: ",round(mean(Dᵢ[:]),digits=6)," Median: ",median(Dᵢ[:]))
  # println("Min: ",minimum(Dᵢ[:])," Max: ",maximum(Dᵢ[:]))
  # println("Perc 95%: ",percentile(Dᵢ[:],95))
  # println("Perc 97.5%: ",percentile(Dᵢ[:],97.5))

  Xᵢ = mat.Xᵢ
  X̂ᵢ = mat.X̂ᵢ
  dμdL = Array{Float64,2}(undef,dim.obs,0)
  dμ̂dL = Array{Float64,2}(undef,dim.obs,0)

  for k=1:dim.Dᵢ
    dμdL = [dμdL Xᵢ[:,k:k].*Dᵢ]
    dμ̂dL = [dμ̂dL X̂ᵢ[:,k:k].*Dᵢ]
  end
  for k=1:dim.νᵢ
    l = dim.Dᵢ + k
    dμdL = [dμdL kron(Xᵢ[:,l:l],ν[k:k,:])]
    dμ̂dL = [dμ̂dL kron(X̂ᵢ[:,l:l],ν[k:k,:])]
  end

  # println([dim.Dᵢ,dim.θ₂,dim.νᵢ,dim.Rᵢ])
  dμdL = reshape(dμdL, (dim.obs,dim.person,dim.θ₂))
  dμ̂dL = reshape(dμ̂dL, (dim.obs,dim.person,dim.θ₂))

  @.mat.νᵢ = ν # Normal random draws
  dμdL3d   = permutedims(dμdL, [1,3,2])
  dμ̂dL3d   = permutedims(dμ̂dL, [1,3,2])

  # Predictions
  s    = [zeros(Float64,dim.prodt[t],) for t=1:dim.market]
  ms   = [zeros(Float64,dim.prodt[t],) for t=1:dim.market]
  mss  = [zeros(Float64,dim.prodt[t],dim.prodt[t]) for t=1:dim.market]
  num  = [zeros(Float64,dim.prodt[t],) for t=1:dim.market]
  den  = [zeros(Float64,dim.person,) for t=1:dim.market]
  adj  = [zeros(Float64,dim.person,) for t=1:dim.market]
  ∂S∂θ = [zeros(Float64,dim.prodt[t],dim.θ₂) for t=1:dim.market]
  pred = predType(s,ms,mss,num,den,adj,∂S∂θ)

  # Fixed point
  eδ1d = dat.ms1d./dat.msOO1d; eδ̂1d = copy(eδ1d)
  eδ = [zeros(Float64,dim.prodt[t],) for t=1:dim.market]

  for t=1:dim.market
    mask = view(indx.mkt,:,t)
    eδ[t] .= eδ1d[mask]
    mat.ms[t] .= dat.ms1d[mask]
    if dim.Dᵢ>0
      mat.X̂[t] .= mat.Xexo[mask,:]
    # else
      # mat.X̂[t] .= [-mat.p̂ʸ[mask] mat.Xexo[mask,:]]
    end
    mat.dμdL[t] .= dμdL3d[mask,:,:]
    mat.dμ̂dL[t] .= dμ̂dL3d[mask,:,:]
  end
  eδ̂ = deepcopy(eδ)
  eμ = [zeros(Float64,dim.prodt[t],dim.person) for t=1:dim.market]
  eμ̂ = [zeros(Float64,dim.prodt[t],dim.person) for t=1:dim.market]
  ∂δ∂θ = [zeros(Float64,dim.prodt[t],dim.θ₂) for t=1:dim.market]
  ∂δ̂∂θ = [zeros(Float64,dim.prodt[t],dim.θ₂) for t=1:dim.market]
  fxp = fxpType(eδ,eδ̂,eδ1d,eδ̂1d,eμ,eμ̂,∂δ∂θ,∂δ̂∂θ)

  # Diags type
  iter     ::Int64 = 0
  iterBLP  ::Int64 = 0
  tolBLP   ::Float64 = 1e-14
  tolOptIV ::Float64 = 1e-8
  gtolSimplex ::Float64 = 1e-6
  gtolBFGS    ::Float64 = 1e-8
  diags = diagsType(iter,iterBLP,tolBLP,tolOptIV,gtolSimplex,gtolBFGS)

  return fxp, pred, diags

end

function setModel(
  model ::Int64,
  Xinput::Matrix{Float64},
  dat   ::datType)

  free    = convert(Vector{Float64},dat.brand.==4)
  prepaid = convert(Vector{Float64},dat.prepost.==1)
  fightbr = convert(Vector{Float64},(dat.brand.>=5) .& (dat.brand.<=7))
  ant4gs  = sum(Xinput[:,5:6],dims=2)
  ant4g   = sum(Xinput[:,5:5],dims=2)
  lowcost = convert(Vector{Float64},(free.==1) .| (fightbr.==1))
  postpay = convert(Vector{Float64},(dat.prepost.==2))
  fbloque = convert(Vector{Float64},(dat.prepost.==3))
  price   = dat.price
  commitl = convert(Vector{Float64},(prepaid==1) .| (lowcost.==1))
  allcall = Xinput[:,9]
  alldata = Xinput[:,10]

  nDᵢ    = 0
  nνᵢ    = 0
  nRᵢ    = 0
  ixP    = [1]
  Xexo   = [Xinput ones(size(Xinput,1))]
  o1     = ones(length(price))
  indxeθ = Vector{Int64}(undef,0)

  if model==0
    Xᵢ = [-price fbloque*10 -ant4gs]
    labΣᵢᵢ = ["Price/y";"Fbloque/y";"Ant4G/y"]
    nDᵢ = 3; nRᵢ = 1
  elseif model==15
    Xᵢ = [-price -ant4gs o1 fbloque*10 prepaid*10]
    labΣᵢᵢ = ["Price/y";"Ant4G/y";"Intercept/y";"Fbloque/y";"Prepaid/y"]
    nDᵢ = 5; nRᵢ = 1
  elseif model==27
    Xᵢ = [-price -ant4gs fbloque    prepaid    o1]
    labΣᵢᵢ = ["Price/y";"Ant4G*ν";"Fbloque*ν";"Prepaid*ν";"Intercept*ν"]
    nDᵢ = 1; nνᵢ = 4; nRᵢ = 5
  else
    error("No model number matched.")
  end

  ### Setting dimensions
  nθ₁   = size(Xexo,2)
  nXᵢ   = size(Xᵢ,2)
  labL  = labΣᵢᵢ
  diagL = collect(1:nXᵢ)
  indxL = [diagL diagL]
  nθ₂   = length(labL)
  nObs  = size(Xᵢ,1)
  nθ    = nθ₁ + nθ₂

  return Xexo,Xᵢ,nXᵢ,nθ₂,nDᵢ,nνᵢ,nRᵢ,labΣᵢᵢ,labL,diagL,indxL,indxeθ,ixP
end

function initialize(
  model   ::Int64,
  iv      ::String,
  nPerson ::Int64,
  ver     ::String,
  impFile ::String,
  simEmp  ::Bool=false;
  lowBound::Bool=false,
  truncate::Bool=false)

  dat, dim, lab, mat, set, indx = importData(model,iv,nPerson,ver,impFile)
  fxp,pred,diags = config!(model,dat,dim,mat,set,indx,simEmp,lowBound,truncate)
  initδ!(fxp.eδ,dat,dim,mat,indx)
  for t=1:dim.market
    if any(isnan.(fxp.eδ[t])) | any(isinf.(fxp.eδ[t]))
      error(" Error: fxp.eδ NaN inside initialize()")
    end
  end
  # println("fxp.eδ initialized inside initialize()")
  return dat, dim, fxp, lab, mat, set, indx, pred, diags
end

function initδ!(
  eδ::Vector{Vector{Float64}},
  dat::datType,
  dim::dimType,
  mat::matType,
  indx::indxType)

  for t=1:dim.market
    mask = view(indx.mkt,:,t)
    @. eδ[t] = mat.ms[t]/dat.msOO1d[mask]
  end
end

function unpackStackedArray(stackedArray::Vector{Matrix{Float64}},
  nRow::Int64, nCol::Int64, nStack::Int64, mask::BitArray{2})

  arr2d = Array{Float64,2}(undef,nRow,nCol)
  for t=1:nStack
    arr2d[view(mask,:,t),:] = stackedArray[t]
  end
  return arr2d
end

function unpackStackedArray!(array2d::Matrix{Float64},
  stackedArray::Vector{Matrix{Float64}}, nStack::Int64, mask::BitArray{2})

  for t=1:nStack
    array2d[view(mask,:,t),:] = stackedArray[t]
  end
end

function unpackStackedArray!(array1d::Vector{Float64},
  stackedArray::Vector{Vector{Float64}}, nStack::Int64, mask::BitArray{2})

  for t=1:nStack
    array1d[view(mask,:,t)] = stackedArray[t]
  end
end

function allocMat!(mat::matType, dim::dimType)
  mat.Z        = Matrix{Float64}(undef,dim.obs,dim.θ)
  mat.XZ       = Matrix{Float64}(undef,dim.θ₁,dim.θ)
  mat.invPhi   = Matrix{Float64}(undef,dim.θ,dim.θ)
  mat.XZPhiZXZ = Matrix{Float64}(undef,dim.θ₁,dim.obs)
end

function updateIVmatrix!(mat::matType)
  mat.Z        .= [mat.∂δ̂∂θ mat.Xexo]
  mat.XZ       .= mat.X'*mat.Z
  mat.invPhi   .= inv(mat.Z'*mat.Z)
  mat.XZPhiZXZ .= (mat.XZ*mat.invPhi*mat.XZ')\(mat.XZ*mat.invPhi*mat.Z')
end

function collectResult(θ₂::Matrix{Float64}, θ₁::Matrix{Float64},
  seθ::Matrix{Float64}, fval::Vector{Float64}, exitflag::BitArray{1})

  # Jstat, idMin = findmin(fval .+ .!(exitflag)*Inf)
  # r = idMin[1]
  # param = paramType(θ₂[:,r],θ₁[:,r],seθ[:,r],Jstat)
  param = paramType(θ₂[:,end],θ₁[:,end],seθ[:,end],fval[end])
  return param
end

function assignFilename(
  fileName::fileNameType,
  model   ::Int64,
  stage   ::Int64,
  iv      ::String,
  nPerson ::Int64,
  impFile ::String)

  pathName = string("out/",impFile,"/m",model,"s",stage,iv,nPerson)
  if ispath(pathName)==true
    rm(pathName, recursive=true)
  end
  mkpath(pathName)

  if stage!=1
    fJLD = string(pathName,"/",fileName.JLD2)
    fCSV = string(pathName,"/",fileName.CSV2)
    fINT = string(pathName,"/",fileName.INT2)
  elseif stage==1
    fJLD = string(pathName,"/",fileName.JLD1)
    fCSV = string(pathName,"/",fileName.CSV1)
    fINT = string(pathName,"/",fileName.INT1)
  end
  fPAR = string(pathName,"/gmmParam.csv")

  return fJLD, fCSV, fINT, fPAR
end

function convert2share(A::Array{Vector{Float64},1}, N::Int64)
  B = Array{SharedArray,1}(N,)
  for n=1:N
    B[n] = convert(SharedArray,A[n])
  end
  return B
end
function convert2share(A::Array{Matrix{Float64},1}, N::Int64)
  B = Array{SharedArray,1}(N,)
  for n=1:N
    B[n] = convert(SharedArray,A[n])
  end
  return B
end

function setIV(
  Xᵢ    ::Matrix{Float64},
  Xexo  ::Matrix{Float64},
  Zinput::Matrix{Float64},
  Dinput::Matrix{Float64},
  price ::Vector{Float64},
  ixP   ::Vector{Int64},
  dim   ::dimType,
  indx  ::indxType,
  iv    ::String)

  z  = hcat(Xexo, Zinput)
  p̂  = z*inv(z'*z)*(z'*price)
  X̂ᵢ = copy(Xᵢ)
  for k in ixP
    X̂ᵢ[:,k] = -p̂ .* (Xᵢ[:,k]./(-price))
  end

  if iv=="opt"
    Z = hcat(randn(dim.obs,dim.θ₂), Xexo)
  elseif iv=="blp"
    Z = z
    @printf("Degrees of freedom: %d\n",size(Z,2)-dim.θ-1)
  elseif (iv=="diff-quad") | (iv=="diff-local")
    # Z = hcat(Xexo, p̂)
    V = hcat(Dinput, p̂)
    D = zero(V)
    for l=1:size(D,2)
      binary = (Set(V[:,l])==Set([0.0,1.0]))
      for t=1:dim.market
        M = view(indx.mkt,:,t)
        x = V[M,l]
        d = x .- x'
        if iv=="diff-quad"
          D[M,l] = sum(d.^2,dims=2)
        elseif iv=="diff-local"
          if binary==false
            dl = d .* (abs.(d) .< std(x))
          else
            dl = d .* (abs.(d) .< std(x))
          end
          D[M,l] = sum(dl,dims=2)
        end
      end
    end
    col0 = vec(sum(abs.(D),dims=1) .!= 0) # Drop zero columns from diff IV
    if iv=="diff-quad"
      Z = hcat(Xexo,D[:,col0])
    elseif iv=="diff-local"
      Z = hcat(Xexo,D[:,col0],p̂)
    end
    @printf("Degrees of freedom: %d\n",size(Z,2)-dim.θ)
    for i in findall(col0.==false)
      @printf("Column %d dropped from Dinput.\n",i)
    end
  end
  return Z, X̂ᵢ, p̂

end

function importθ̂₂(
  model  ::Int64,
  stage  ::Int64,
  mcid   ::Int64,
  iv     ::String,
  nPerson::Int64,
  impF   ::String,
  dim    ::dimType,
  update ::Bool=true)

  if (update==false) & (stage==2)
    pathName = string("out/",impF,"/m",model)
    if iv=="opt"
      if impF=="optIVtest"
        fnameStage1 = string(pathName,"s0opt",nPerson,"/gmmStage2.csv")
      elseif (impF=="optIV1") & (impF=="base")
        fnameStage1 = string(pathName,"s2opt",nPerson,"/gmmStage2.csv")
      elseif (impF=="optIV2") & (impF=="base")
        # fnameStage1 = string("./out/optIV1/m",model,"s2opt","/gmmStage2.csv")
        fnameStage1 = string(pathName,"s2opt",nPerson,"/gmmStage2.csv")
      else
        fnameStage1 = string(pathName,"s1blp",nPerson,"/gmmStage1.csv")
      end
    else
      fnameStage1 = string(pathName,"s1",iv,nPerson,"/gmmStage1.csv")
    end

    imp = readdlm(fnameStage1,',')
    θ̂  = convert(Vector{Float64},imp[mcid,2:(end-2)])
    θ̂₂ = θ̂[1:dim.θ₂]
    println("First stage θ̂₂: ",round.(θ̂₂,digits=4))
  else
    θ̂₂ = zeros(Float64,0)
  end
  return θ̂₂
end

end # end of module

max_stars_repo_name: yutec/FightBrand

max_issues_repo_name: yutec/FightBrand

max_forks_repo_name: yutec/FightBrand


text:
# Julia wrapper for the sgp4 Python library:
# https://pypi.python.org/pypi/sgp4/

module SGP4

using PyCall, Dates

import Base.getindex

const sgp4io = PyNULL()
const earth_gravity = PyNULL()
const sgp4_propagation = PyNULL()

function __init__()
    copy!(sgp4io, pyimport_conda("sgp4.io", "sgp4", "conda-forge"))
    copy!(earth_gravity, pyimport("sgp4.earth_gravity"))
    copy!(sgp4_propagation, pyimport("sgp4.propagation"))
end

export GravityModel,
       twoline2rv,
       propagate

struct GravityModel
    model::PyObject # can be any of {wgs72old, wgs72, wgs84}
end

mutable struct SGP4Sat
    s::PyObject
end
getindex(sat::SGP4Sat, sym::Symbol) = sat.s[sym]

GravityModel(ref::AbstractString) = GravityModel(earth_gravity[ref])

# sgp4.io convenience functions
function twoline2rv(line1::String, line2::String, grav::GravityModel)
    return SGP4Sat(sgp4io["twoline2rv"](line1,line2,grav.model))
end

function sgp4( sat::SGP4Sat,
               dtmin::Real )
    r,v = sgp4_propagation["sgp4"](sat.s, dtmin)
    return vcat(r...), vcat(v...)
end

"""
Propagate the satellite by `dtmin` minutes.

Returns (position, velocity)
"""
function propagate( sat::SGP4Sat,
                    dtmin::Real )
    sgp4(sat,dtmin)
end

function propagate(sats::Vector{SGP4Sat},
                   dtmin::Real)
    f = x->propagate(x, dtmin)
    f.(sats)
end

"""
Propagate the satellite from its epoch to the date/time specified

Returns (position, velocity) at the specified time
"""
function propagate( sat::SGP4Sat,
                    year::Real,
                    month::Real,
                    day::Real,
                    hour::Real,
                    min::Real,
                    sec::Real )
    (pos, vel) = sat.s[:propagate](year,month,day,hour,min,sec)

    # check for errors
    if sat.s[:error] != 0
        println(sat.s[:error_message])
    end

    return ([pos...],[vel...])
end

function propagate( sat::SGP4Sat,
                    t::DateTime )
    propagate(sat,
              Dates.year(t),
              Dates.month(t),
              Dates.day(t),
              Dates.hour(t),
              Dates.minute(t),
              Dates.second(t) + Dates.millisecond(t)/1000)
end

function propagate(sat::SGP4Sat,
                   t::AbstractVector{DateTime})
    pos = zeros(3, length(t))
    vel = zeros(3, length(t))

    for (idx, ti) in enumerate(t)
        pos[:,idx],vel[:,idx] = propagate(sat, ti)
    end
    return (pos,vel)
end

"Generate a satellite ephemeris"
function propagate( sat::SGP4Sat,
                    tstart::DateTime,
                    tstop::DateTime,
                    tstep::Dates.TimePeriod )
    propagate(sat, tstart:tstep:tstop)
end

"tstep specified in seconds"
function propagate( sat::SGP4Sat,
                    tstart::DateTime,
                    tstop::DateTime,
                    tstep::Real )
    propagate(sat,tstart,tstop,Dates.Second(tstep))
end

"Propagate many satellites to a common time"
function propagate( sats::Vector{SGP4Sat},
                    year::Real,
                    month::Real,
                    day::Real,
                    hour::Real,
                    min::Real,
                    sec::Real )
    f = x->propagate(x, year, month, day, hour, min, sec)
    f.(sats)
end

function propagate( sats::Vector{SGP4Sat},
                    t::DateTime )
    f = x->propagate(x, Dates.year(t), Dates.month(t), Dates.day(t), Dates.hour(t), Dates.minute(t), Dates.second(t) + Dates.millisecond(t)/1000)
    f.(sats)
end

end #module