text:
if !@isdefined(TabLevel)
    TabLevel = ""
end
println(TabLevel*"Open OpInnerProd.jl")
TabLevel=TabLevel*"    "

function OpInnerProd(State1,MiddleOp,State2)
    OpInnerProd(State1,MiddleOp*State2)
end

function OpInnerProd(State1::Array,State2::Array)
    ###Work out the desired scale
    MyType=typeof(one(OpScaleType(State1))+one(OpScaleType(State2)))
    IPMat=zeros(MyType,length(State1),length(State2))
    for I1 in 1:length(State1)
        for I2 in 1:length(State2)
            IPMat[I1,I2]=OpInnerProd(State1[I1],State2[I2])
        end
    end
    return IPMat
end

function OpInnerProd(x::Number,OP::MajoranaOp)
    if x == 0
        return 0
    else
        OpInnerProd(MajScale(x),OP)
    end
end
function OpInnerProd(OP::MajoranaOp,x::Number)
    if x == 0
        return 0
    else
        OpInnerProd(OP,MajScale(x))
    end
end

function OpInnerProd(State2::MajoranaSp,OP1::MajoranaOp)
    ##Complex conjugate
    return conj(OpInnerProd(OP1,State2))
end

function OpInnerProd(OP1::MajoranaOp,State2::MajoranaSp)
    ###Sum over all the inner product terms
    sum=0im
    for OP2 in State2.SpList
        sum+=OpInnerProd(OP1,OP2)
    end
    return sum
end

function OpInnerProd(State1::MajoranaSp,State2::MajoranaSp)
    ###Sum over all the inner product terms
    sum=0im
    for OP1 in State1.SpList
        sum+=OpInnerProd(OP1,State2)
    end
    return sum
end

function OpInnerProd(State1::MajoranaOp,State2::MajoranaOp)
    if State1.scale == 0 || State2.scale == 0
        return 0  ###If one is zero the inner product can safely be taken to be zero.
    elseif !State1.KetState ||  !State2.KetState
        ###The states need to be ket-states
        throw(BoundsError("Enclosing elements of the inner product should be kets!"))
    end
    if State1.scale==0 || State2.scale==0
        return 0
    end
    ### Peel of the last element of state1, reverse the order and conjugate the scale
    #print("State1: ")
    #OpPrintln(State1)
    ConjState1=MajoranaOp(conj(State1.scale),State1.OpList[end:-1:1])
    #print("ConjState1: ")
    #OpPrintln(ConjState1)
    ReducedState=ConjState1*State2
    #print("ReducedState: ")
    #OpPrintln(ReducedState)
    if length(ReducedState.OpList)==0
        ###If all the elements cancel then there is an inner product
        return ReducedState.scale
    else
        ### Otherwise not
        return 0
    end
end

TabLevel=TabLevel[1:end-4]
println(TabLevel*"Close OpInnerProd.jl")


text:
function plot_traj!(plt::Plots.Plot, traj::SystemTrajectory, g::AbstractGame,
                    player_colors::AbstractArray, args...; kwargs...)
    return  plot_traj!(plt, traj, xyindex(g), player_colors, args...; kwargs...)
end

function plot_traj(traj::SystemTrajectory)
    return plot_traj(traj, [:black], tuple(@S(1:length(eltype(traj.u)))))
end
function plot_traj(traj::SystemTrajectory, player_colors::Vector, uids)
    plots = []

    if !isnothing(uids)
        push!(plots, plot_inputs(traj, player_colors, uids))
    end

    push!(plots, plot_states(traj))

    return plot(plots...)
end

function plot_inputs(traj::SystemTrajectory)
    return plot_inputs(traj, [:black], tuple(@S(1:length(eltype(traj.u)))))
end
function plot_inputs(traj::SystemTrajectory, player_colors, uids; kwargs...)
    # names for each input
    nu = length(eltype(traj.u))
    input_labels = reshape([latexstring("u_$i") for i in 1:nu], 1, nu)
    # find the player color for each input
    input_colors = reshape([player_colors[findfirst(in.(i, uids))] for i in 1:nu], 1, nu)
    return plot(hcat(traj.u...)'; layout=(nu, 1),
                legend=:none,
                xlabel="time step",
                ylabel=input_labels,
                seriescolor=input_colors, kwargs...)
end

function plot_states(traj::SystemTrajectory)
    nx = length(eltype(traj.x))
    state_labels = reshape(["x$i" for i in 1:nx], 1, nx)
    return plot(hcat(traj.x...)'; layout=(nx, 1), label=state_labels,
                seriescolor=:black)
end

function plot_traj!(plt::Plots.Plot, traj::SystemTrajectory, xy_ids,
                    player_colors::AbstractArray,
                    uids=nothing,
                    plot_attributes=NamedTuple()
                    ; k::Int=1, kp=nothing)
    # buffer for all the plots
    default_plot_attributes = (legend=:none, seriesalpha=1.,
                               xlabel=L"p_x [m]", ylabel=L"p_y [m]",
                               xlims=(-3.5, 3.5), ylims=(-3.5, 3.5),
                               aspect_ratio=:equal)
    plot_attributes = merge(default_plot_attributes, plot_attributes)

    nu = length(eltype(traj.u))
    pu = plot(; layout=(nu, 1))

    if !isnothing(uids)
        pu = plot_inputs(traj, player_colors, uids)
    end

    if length(xy_ids) == 1
        player_colors = [:black for p in player_colors]
    end

    function mark_position!(p::Plots.Plot, px, py, marker)
        # marker at the current time step
        scatter!(plt, [px], [py]; marker=marker, legend=:none)
    end

    for (i, xy_i) in enumerate(xy_ids)
        # the trajectory
        pxs = collect(x[first(xy_i)] for x in traj.x)
        pys = collect(x[last(xy_i)] for x in traj.x)
        plot!(plt, pxs, pys; seriescolor=player_colors[i], plot_attributes...)
        if k > 0
            mark_position!(plt, pxs[k], pys[k], (:diamond, player_colors[i]))
        end
        if !isnothing(kp) && kp > 0
            mark_position!(plt, pxs[kp], pys[kp], (:circle, player_colors[i]))
        end
    end

    return isnothing(uids) ? plt : plot(pu, plt)
end

plot_traj(args...; kwargs...) = begin p = plot(); plot_traj!(p, args...; kwargs...) end

function scatter_positions(args...; kwargs...)
    plt = plot()
    return scatter_positions!(plt, args; kwargs...)
end
function scatter_positions!(plt, x, g, player_colors, marker=(:circle))
    for (i, xy_i) in enumerate(xyindex(g))
        scatter!(plt, [x[xy_i][1]], [x[xy_i][2]], seriescolor=player_colors[i],
                 marker=marker)
    end
    return plt
end

# cost plots
function plot_cost(g::AbstractGame, op::SystemTrajectory, dims, i::Int=1,
                   st::Symbol=:contour; k::Int=1)
    lqg = lq_approximation(g, op)
    nx = n_states(g)
    nu = n_controls(g)
    t = k * samplingtime(g)

    offset2vec(Δd1, Δd2) = begin
        Δx = zeros(nx)
        Δx[dims[1]] = Δd1
        if length(dims) == 2
            Δx[dims[2]] = Δd2
        end
        return SVector{nx}(Δx)
    end

    projected_cost(Δd1, Δd2=0) = begin
        Δx = offset2vec(Δd1, Δd2)
        return player_costs(g)[i](g, op.x[k]+Δx, op.u[k], t)
    end

    projected_cost_approx(Δd1, Δd2=0) = begin
        Δx = offset2vec(Δd1, Δd2)
        c0 = player_costs(g)[i](g, op.x[k], op.u[k], t)
        return player_costs(lqg)[k][i](Δx, (@SVector zeros(nu))) + c0
    end

    Δd1_range = Δd2_range = -10:0.1:10

    if length(dims) == 1
        p = plot(Δd1_range, projected_cost, label="g")
        plot!(p, Δd1_range, projected_cost_approx, label="lqg")
        return p
    elseif length(dims) == 2
        return plot(Δd1_range, Δd2_range, projected_cost, st=st)
    end

    @assert false "Can only visualize one or two dimensions"
end

function animate_plot(plot_frame::Function, plot_args...;
                      k_range::UnitRange, frame_sample::Int=1, fps::Int=10,
                      filename::String="$(@__DIR__)/../debug_out/test.gif")
    anim = @animate for k in k_range
        plot_frame(plot_args...; k=k)
    end every frame_sample

    return gif(anim, filename; fps=fps)
end

ext: jl


text:
using Distributions
using PyPlot
using MAT

srand(0) #seed

function generate_gaussian_data(N::Int64, D::Int64, K::Int64)
	x = randn(D,N)
	tpi = rand(Dirichlet(ones(K)))
	tzn = rand(Multinomial(N,tpi))  #number of points in each cluster, e.g. [2, 2, 6]
	tz = zeros(N)

	tmean = zeros(D,K)
	tcov = zeros(D,D,K)

	ind = 1
	println(tzn)
	for i=1:length(tzn)
		indices = ind:ind+tzn[i]-1
		tz[indices] = i
		tmean[:,i] = rand(MvNormal(zeros(D), 100*eye(D)))  #sample cluster mean
		tcov[:,:,i] = rand(InverseWishart(D+2, eye(D)))    #sample cluster covariance
		# T = chol(slice(tcov,:,:,i))
		# x[:,indices] = broadcast(+, T'*x[:,indices], tmean[:,i]);
		d = MvNormal(tmean[:,i], tcov[:,:,i])  #form cluster distribution
		for j=indices
			x[:,j] = rand(d)  #generate cluster point
		end
		ind += tzn[i]
	end
	x, tz, tmean, tcov
end

N = 1000
D = 2
K = 6
latent, latent_z, latent_mu, latent_cov = generate_gaussian_data(N, D, K)

figure()
scatter(latent[1,:], latent[2,:], c=latent_z)

file = matopen("clusters.mat", "w")
write(file, "data", latent)
write(file, "labels", latent_z)
close(file)


text:
#! /usr/bin/julia

# Rosetta Code, MD5

using Nettle

function md5sum(s::String)
    bytes2hex(md5_hash(s))
end

function Base.trunc(s::String, n::Integer)
    0 < n || throw(DomainError())
    len = length(s)
    len > n || return s
    3 < n || return s[1:n]
    return s[1:n-3]*"..."
end

tests = ["" => "d41d8cd98f00b204e9800998ecf8427e",
         "a" => "0cc175b9c0f1b6a831c399e269772661",
         "abc" => "900150983cd24fb0d6963f7d28e17f72",
         "message digest" => "f96b697d7cb7938d525a2f31aaf161d0",
         "abcdefghijklmnopqrstuvwxyz" => "c3fcd3d76192e4007dfb496cca67e13b",
         "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789" =>
         "d174ab98d277d9f5a5611c2c9f419d9f",
         "12345678901234567890123456789012345678901234567890123456789012345678901234567890" =>
         "57edf4a22be3c955ac49da2e2107b67a",
         "foobad" => "3858f62230ac3c915f300c664312c63f"]

println("Testing Julia's MD5 hash against RFC 1321.")
for k in sort(collect(keys(tests)), by=(x)->length(x))
    mysum = md5sum(k)
    print(@sprintf("    %15s => ", trunc(k, 15)), mysum)
    if mysum == tests[k]
        println(" MD5 OK")
    else
        println(" MD5 Bad")
        println("    The sum should be    ", tests[k])
    end
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
# Generate quadrature information for B-spline product.
# Returns weights and nodes for integration in [-1, 1].
#
# See https://en.wikipedia.org/wiki/Gaussian_quadrature.
#
# Some notes:
#
# - On each interval between two neighbouring knots, each B-spline is a
#   polynomial of degree (k - 1). Hence, the product of two B-splines has degree
#   (2k - 2).
#
# - The Gauss--Legendre quadrature rule is exact for polynomials of degree
#   <= (2n - 1), where n is the number of nodes and weights.
#
# - Conclusion: on each knot interval, `k` nodes should be enough to get an
#   exact integral. (I verified that results don't change when using more than
#   `k` nodes.)
#
# Here, p is the polynomial order (p = 2k - 2 for the product of two B-splines).
function _quadrature_prod(::Val{p}) where {p}
    n = cld(p + 1, 2)
    _gausslegendre(Val(n))
end

# Compile-time computation of Gauss--Legendre nodes and weights.
# If the @generated branch is taken, then the computation is done at compile
# time, making sure that no runtime allocations are performed.
function _gausslegendre(::Val{n}) where {n}
    if @generated
        vecs = _gausslegendre_impl(Val(n))
        :( $vecs )
    else
        _gausslegendre_impl(Val(n))
    end
end

function _gausslegendre_impl(::Val{n}) where {n}
    data = gausslegendre(n)
    map(SVector{n}, data)
end

# Metric for integration on interval [a, b].
# This is to transform from integration on the interval [-1, 1], in which
# quadratures are defined.
# See https://en.wikipedia.org/wiki/Gaussian_quadrature#Change_of_interval
struct QuadratureMetric{T}
    α :: T
    β :: T

    function QuadratureMetric(a::T, b::T) where {T}
        α = (b - a) / 2
        β = (a + b) / 2
        new{T}(α, β)
    end
end

# Apply metric to normalised coordinate (x ∈ [-1, 1]).
Base.:*(M::QuadratureMetric, x::Real) = M.α * x + M.β
Broadcast.broadcastable(M::QuadratureMetric) = Ref(M)

max_stars_repo_name: jipolanco/BasisSplines

max_issues_repo_name: jipolanco/BasisSplines

max_forks_repo_name: jipolanco/BasisSplines


text:
using DynamicRemoteCall
using Distributed

addprocs(3)

@everywhere function testfunc()
    r = rand()/10
    sleep(r)
    return r
end

results = zeros(100)

begin
    starttime = time()
    dynamicremotecall!(testfunc, results)
    endtime = time()
    (endtime-starttime)*nworkers()/sum(results)
end

max_stars_repo_name: wangl-cc/DynamicRemoteCall-jl

max_issues_repo_name: wangl-cc/DynamicRemoteCall-jl

max_forks_repo_name: wangl-cc/DynamicRemoteCall-jl


text:
using DarkSky
using Test
using Dates

if haskey(ENV, "DARKSKY_API_KEY")
    @info("Dark Sky API key was found in your system environment variables, running tests...")
    @test isa(forecast(42.3601, -71.0589), DarkSky.DarkSkyResponse)
    @test isa(forecast(42.3601, -71.0589, exclude=["currently"], extend="hourly"), DarkSky.DarkSkyResponse)
    response = forecast(42.3601, -71.0589, DateTime(2018, 3, 7, 0, 0, 0));
    @test isa(response, DarkSky.DarkSkyResponse)
    @test latitude(response) === 42.3601
    @test longitude(response) === -71.0589
    @test timezone(response) == "America/New_York"
    @test isa(currently(response), Dict{String, Any})
    @test isa(minutely(response), Nothing)
    @test isa(daily(response), Dict{String, Any})
    @test isa(alerts(response), Array{Any, 1}) | isa(alerts(response), Nothing)
    @test isa(flags(response), Dict{String, Any})
    @test isa(print(response), Nothing)
    dict_response = forecast(42.3601, -71.0589, DateTime(2018, 3, 7, 0, 0, 0), out_type="");
    @test Dict(response) == dict_response
    @test Dict(response) == convert(Dict, response)
    @test isa(forecast(42.3601, -71.0589, DateTime(2018, 1, 1, 0, 0, 0), lang="es", units="si", exclude=["minutely"]), DarkSky.DarkSkyResponse)
else
    @warn "Dark Sky API key was not found in your system environment variables, skipping tests..."
end


text:
include("../src/star.jl")
using DistributedArrays

function exclusive_cumsum(x)
  sum = 0
  for i = 1 : length(x)
    temp = x[i]
    x[i] = sum
    sum = sum + temp
  end
  x
end

function cumsum(initial_value, x)
  if length(x) > 0
    x[1] = initial_value + x[1]
  end
  cumsum!(x, x)
end

# computes cumsum on a distributed array x
cumsum_star(x) = Star.scan(cumsum, exclusive_cumsum(Star.reduce(sum, x)), x)

function cumsum_star_test(x)
  x_ = convert(Array{Float64, 1}, x)
  t = @elapsed cumsum_star(x)
  out_ = convert(Array{Float64, 1}, x)
  t1 = @elapsed cumsum!(x_, x_)
  println("parallel time ", t, " serial time ", t1)
  println("rel norm ", norm(out_-x_)/abs(x_[end]))
end

function cumsum_no_star(x)
  ref = Array(Any, nworkers())
  for i = 1 : nworkers()
    ref[i] = @spawnat workers()[i] sum(localpart(x))
  end
  y = [fetch(r) for r in ref]
  exclusive_cumsum(y)
  @sync begin
    for i = 1 : nworkers()
      @spawnat workers()[i] cumsum(y[i], localpart(x))
    end
  end
end


text:
module ACL

export
    floorsum

include("math/floorsum.jl")

end


text:
module TestConversions

using Test, DataFrames, InteractiveUtils
using DataStructures: OrderedDict, SortedDict
const ≅ = isequal

df = DataFrame()
df[:A] = 1:5
df[:B] = [:A, :B, :C, :D, :E]
@test isa(convert(Matrix, df), Matrix{Any})
@test isa(convert(Matrix{Any}, df), Matrix{Any})

df = DataFrame()
df[:A] = 1:5
df[:B] = 1.0:5.0
@test isa(convert(Matrix, df), Matrix{Float64})
@test isa(convert(Matrix{Any}, df), Matrix{Any})
@test isa(convert(Matrix{Float64}, df), Matrix{Float64})
@test isa(Matrix(df), Matrix{Float64})
@test isa(Matrix{Any}(df), Matrix{Any})
@test isa(Matrix{Float64}(df), Matrix{Float64})

df = DataFrame()
df[:A] = Vector{Union{Float64, Missing}}(1.0:5.0)
df[:B] = Vector{Union{Float64, Missing}}(1.0:5.0)
a = convert(Matrix, df)
aa = convert(Matrix{Any}, df)
ai = convert(Matrix{Int}, df)
@test isa(a, Matrix{Union{Float64, Missing}})
@test a == convert(Array, convert(Matrix{Union{Float64, Missing}}, df))
@test a == convert(Matrix, df)
@test a == Matrix(df)
@test isa(aa, Matrix{Any})
@test aa == convert(Matrix{Any}, df)
@test aa == Matrix{Any}(df)
@test isa(ai, Matrix{Int})
@test ai == convert(Matrix{Int}, df)
@test ai == Matrix{Int}(df)

df[1,1] = missing
@test_throws ErrorException convert(Matrix{Float64}, df)
na = convert(Matrix{Union{Float64, Missing}}, df)
naa = convert(Matrix{Any}, df)
nai = convert(Matrix{Union{Int, Missing}}, df)
@test isa(na, Matrix{Union{Float64, Missing}})
@test na ≅ convert(Matrix, df)
@test na ≅ Matrix(df)
@test isa(naa, Matrix{Union{Any, Missing}})
@test naa ≅ convert(Matrix{Any}, df)
@test naa ≅ Matrix{Any}(df)
@test isa(nai, Matrix{Union{Int, Missing}})
@test nai ≅ convert(Matrix{Union{Int, Missing}}, df)
@test nai ≅ Matrix{Union{Int, Missing}}(df)

a = Union{Float64, Missing}[1.0,2.0]
b = Union{Float64, Missing}[-0.1,3]
c = Union{Float64, Missing}[-3.1,7]
di = Dict("a"=>a, "b"=>b, "c"=>c)

df = convert(DataFrame, di)
@test isa(df, DataFrame)
@test names(df) == [Symbol(x) for x in sort(collect(keys(di)))]
@test df[:a] == a
@test df[:b] == b
@test df[:c] == c

od = OrderedDict("c"=>c, "a"=>a, "b"=>b)
df = convert(DataFrame,od)
@test isa(df, DataFrame)
@test names(df) == [Symbol(x) for x in keys(od)]
@test df[:a] == a
@test df[:b] == b
@test df[:c] == c

sd = SortedDict("c"=>c, "a"=>a, "b"=>b)
df = convert(DataFrame,sd)
@test isa(df, DataFrame)
@test names(df) == [Symbol(x) for x in keys(sd)]
@test df[:a] == a
@test df[:b] == b
@test df[:c] == c

a = [1.0]
di = Dict("a"=>a, "b"=>b, "c"=>c)
@test_throws DimensionMismatch convert(DataFrame,di)

end # module


text:
using ProfileView
using Statistics

using SubmodularMaximization

num_agents = 4
num_sensors = 6
nominal_area = 2.0

sensor_radius = sqrt(nominal_area / (num_agents * pi))
station_radius = 2 * sensor_radius

agent_specification = CircleAgentSpecification(sensor_radius, station_radius,
                                         num_sensors)

agents = generate_agents(agent_specification, num_agents)

f(x) = mean_area_coverage(x, 100)
problem = ExplicitPartitionProblem(f, agents)

solution = solve_optimal(problem)

times = Float64[]

function foo(n)
  global times = Float64[]

  for i in 1:n
    agents = generate_agents(agent_specification, num_agents)
    problem = ExplicitPartitionProblem(f, agents)

    time = @elapsed solution = solve_optimal(problem)

    @show solution.value

    push!(times, time)
  end
end

@profview foo(1)
ProfileView.closeall()
@profview foo(50)

println("Mean time: $(mean(times))")

PyPlot.plt.hist(times, 20)

nothing

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
using Epidemia
using Documenter

DocMeta.setdocmeta!(Epidemia, :DocTestSetup, :(using Epidemia); recursive=true)

makedocs(;
    modules=[Epidemia],
    authors="Andreas Koher",
    repo="https://github.com/andreaskoher/Epidemia.jl/blob/{commit}{path}#{line}",
    sitename="Epidemia.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://andreaskoher.github.io/Epidemia.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/andreaskoher/Epidemia.jl",
)


text:
include("./ZMQ_utils.jl")
include("./test_functions.jl")
# include("./user_functions.jl")

"""
    server_0mq4lv(fns=(;); initOK=false)

`server_0mq4lv` (meaning zero MQ for LabView) is the top-level function of the package.
User functions must be supplied as a NamedTuple or Dict. Start ZMQ socket, listen to
requests, parse them, execute the requested user function, send response. Repeat.
# Arguments
- `initOK`: set it to `true` to skip check of the user script
# Examples
```julia-repl

julia> server_0mq4lv((;foo=foo, bar, baz))
```
"""
function server_0mq4lv(fns=(;); initOK=false)

    context = Context()
    socket = Socket(context, REP)
    ZMQ.bind(socket, "tcp://*:5555")
    version = string(PkgVersion.Version(LabView0mqJl))
    if isempty(fns)
        fnlist = "built-in test functions only"
    else
        fnlist = join(["built-in test functions", keys(fns)...], ", ")
    end

    external_fns(;) = (;ext_fns=keys(fns))  # external_fns(;kwargs...)

    # println(external_fns())
    # println(external_fns)
    # println(eval(:external_fns))

    global scriptexists
    global scriptOK

    if ! initOK
        if isempty(fns)
            initOK = true # package-internal functions only
        else
            initOK = scriptexists && scriptOK
        end
    end

    if initOK
        println("starting ZMQ server, Julia for LabVIEW: LabView0mqJl v=$version")
        println("available functions:")
        println(fnlist)
    else
        println("server initialisation failed")
        @show scriptexists scriptOK
    end

    try
        while true
            # Wait for next request from client

            bytesreceived = ZMQ.recv(socket)

            cmnd = parse_cmnd(bytesreceived)
            response = UInt8[]
            if cmnd.command == :stop
                response = UInt8.([1, PROTOC_V])
            elseif !initOK
                if !isnothing(scriptexcep)
                    excep, stack_trace = scriptexcep
                else
                    stack_trace = []
                    excep = "Julia error on initialisation script"
                end
                err = nothing
                try
                    err = build_err_info(;
                        err = true,
                        errcode = 5235817,
                        source = @__FILE__,
                        stack_trace = err_stack_trace,
                        excep = excep,
                    )
                catch
                    err = build_err_info(;
                        err = true,
                        errcode = 5235817,
                        source = @__FILE__,
                        excep = "Julia error on initialisation script",
                    )
                end
                response = puttogether(; err = err, returncode = 3)
            elseif cmnd.command == :ping
                response = UInt8.([2, PROTOC_V])
            elseif cmnd.command == :callfun
                try
                    pr = parse_REQ(bytesreceived)
                    fn = pr.fun2call
                    if fn == :external_fns
                        f = external_fns
                    elseif haskey(fns, fn)
                        f = fns[fn]
                    else
                        f = eval(fn)
                    end
                    y = f(; pr.args...)
                    response = puttogether(; y = y)
                catch excep
                    err_stack_trace = stacktrace(catch_backtrace())
                    # # https://docs.julialang.org/en/v1/manual/stacktraces/#Error-handling-1
                    err = build_err_info(;
                        err = true,
                        errcode = 5235805,
                        source = @__FILE__,
                        stack_trace = err_stack_trace,
                        excep = excep,
                    )
                    response = puttogether(; err = err, returncode = 3)
                end
            else
                response = UInt8.([255, PROTOC_V])
            end

            # Send reply back to client
            ZMQ.send(socket, response)

            cmnd.command == :stop && break
        end

    finally
        # clean up
        ZMQ.close(socket)
        ZMQ.close(context)
    end # try

end

# server_0mq4lv()


text:
interval_grids = [EquispacedGrid,
    PeriodicEquispacedGrid, MidpointEquispacedGrid,
    UnitEquispacedGrid, UnitMidpointEquispacedGrid,
    UnitPeriodicEquispacedGrid,
    ChebyshevNodes, ChebyshevExtremae,
    ChebyshevUNodes, LegendreNodes]

function test_grids(T)
    ## Equispaced grids
    len = 21
    a = -T(1.2)
    b = T(3.5)
    g1 = EquispacedGrid(len, Interval(a, b))
    g2 = PeriodicEquispacedGrid(len, Interval(a-1, b+1))
    g3 = g1 × g2
    g4 = g1 × g3

    test_generic_grid(g1)
    test_generic_grid(g2)
    test_generic_grid(g3)
    test_generic_grid(g4)

    # Test a subgrid
    g5 = g1[10:20]
    @test g5[1] == g1[10]
    @test g5[11] == g1[20]
    @test length(g5) == 20-10+1
    test_generic_grid(g5)
    g6 = g1[10:2:20]
    @test g6[2] == g1[12]
    @test length(g6) == 6

    g = EquispacedGrid(len, a, b)
    idx = 5
    @test g[idx] ≈ a + (idx-1) * (b-a)/(len-1)
    @test g[len] ≈ b
    @test_throws BoundsError g[len+1] == b

    # Test iterations
    (l,s) = grid_iterator1(g)
    @assert s ≈ len * (a+b)/2

    ## Periodic equispaced grids
    len = 120
    a = -T(1.2)
    b = T(3.5)
    g = PeriodicEquispacedGrid(len, a, b)

    idx = 5
    @test g[idx] ≈ a + (idx-1) * (b-a)/len
    @test g[len] ≈ b - step(g)
    @test_throws BoundsError g[len+1] == b

    (l,s) = grid_iterator1(g)
    @test s ≈ (len-1)*(a+b)/2 + a

    (l,s) = grid_iterator2(g)
    @test s ≈ (len-1)*(a+b)/2 + a

    ## Tensor product grids
    len = 11
    g1 = PeriodicEquispacedGrid(len, -one(T), one(T))
    g = g1^2
    @test isperiodic(g1)
    @test isperiodic(g)
    @test size(g) == (length(g1),length(g1))
    @test productgrid(g1) ==g1

    g2 = EquispacedGrid(len, -one(T), one(T))
    @test !isperiodic(g2)
    g = g1 × g2
    @test !isperiodic(g)
    @test length(g) == length(g1) * length(g2)
    @test size(g) == (length(g1),length(g2))
    @test size(g,1) == length(g1)

    @test component(g, 1) == g1
    @test component(g, 2) == g2
    @test covering(g) ≈ covering(g1)×covering(g2)

    idx1 = 5
    idx2 = 9
    x1 = g1[idx1]
    x2 = g2[idx2]
    x = g[idx1,idx2]
    @test x[1] ≈ x1
    @test x[2] ≈ x2

    (l,s) = grid_iterator1(g)
    @test s ≈ -len

    (l,s) = grid_iterator2(g)
    @test s ≈ -len

    # Test a tensor of a tensor
    g3 = g × g2
    idx1 = 5
    idx2 = 7
    idx3 = 4
    x = g3[idx1,idx2,idx3]
    x1 = g1[idx1]
    x2 = g2[idx2]
    x3 = g2[idx3]
    @test x[1] ≈ x1
    @test x[2] ≈ x2
    @test x[3] ≈ x3

    # Test a mapped grid
    m = mapto(T(0)..T(1), T(2)..T(3))
    # Make a MappedGrid by hand because map_grid would simplify
    mg1 = MappedGrid(m, PeriodicEquispacedGrid(30, T(0), T(1)))
    test_generic_grid(mg1)
    # Does map_grid simplify?
    mg2 = map_grid(m, PeriodicEquispacedGrid(30, T(0), T(1)))
    @test typeof(mg2) <: PeriodicEquispacedGrid
    @test infimum(covering(mg2)) ≈ T(2)
    @test supremum(covering(mg2)) ≈ T(3)

    # Apply a second map and check whether everything simplified
    m2 = mapto(T(2)..T(3), T(4)..T(5))
    mg3 = map_grid(m2, mg1)
    @test infimum(covering(mg3)) ≈ T(4)
    @test supremum(covering(mg3)) ≈ T(5)
    @test mg3 isa PeriodicEquispacedGrid

    # Scattered grid
    pts = rand(T, 10)
    sg = ScatteredGrid(pts)
    test_generic_grid(sg)
end

function test_laguerre(T)
    grid = LaguerreNodes(10,rand(T))
    @test infimum(covering(LaguerreNodes(0.,rand(10)))) == 0
    @test supremum(covering(LaguerreNodes(0.,rand(10)))) == Inf
    test_generic_grid(grid)
end

function test_hermite(T)
    grid = HermiteNodes(10)
    @test DomainSets.FullSpace{T}()== covering(HermiteNodes(rand(T,10)))
    test_generic_grid(grid)
end

function test_jacobi(T)
    grid = JacobiNodes(10,rand(T),rand(T))
    @test covering(JacobiNodes(T(0),T(0),rand(T,10))) == ChebyshevInterval{T}()
    test_generic_grid(grid)
    @test JacobiNodes(10,zero(T),zero(T)) ≈ LegendreNodes(10)
    @test JacobiNodes(10,T(1//2),T(1//2)) ≈ ChebyshevUNodes(10)
    @test JacobiNodes(10,T(-1//2),T(-1//2)) ≈ ChebyshevTNodes(10)
end

for T in types
    delimit(string(T))
    for GRID in interval_grids
        @testset "$(string(GRID))" begin
            if GRID <: GridArrays.AbstractEquispacedRangeGrid
                g = GRID(10,UnitInterval{T}())
            else
                g = GRID(10)
            end
            test_interval_grid(g)
        end
    end

    for grid in (JacobiNodes(10,rand(T),rand(T)),)
        @testset "JacobiNodes" begin
            test_interval_grid(grid)
        end
    end

    @testset "HermiteNodes" begin
        test_hermite(T)
    end

    @testset "Laguerre" begin
        test_laguerre(T)
    end

    @testset "Jacobi" begin
        test_jacobi(T)
    end

    @testset "Specific grid tests" begin
        test_grids(T)
    end
end


text:
# Canonical test models

"""
    lotka!(du, u, p, t)

lotka-volterra, inplace
"""
function lotka!(du, u, p, t)
	x, y = u
	a, b, c, d = p
	du[1] = a*x - b*x*y
	du[2] = -c*y + d*x*y
	nothing
end

# Murray Biomolecular Feedback Systems http://www.cds.caltech.edu/~murray/BFSwiki/index.php/Main_Page
# Ch 3 examples

"Example 3.1 (Transcriptional component)"
function transcriptional_component!(du, u, p, t)
    ctrl = 5. # idk
    α, Κ, δ, κ, γ, n = p
    du[1] = α / (1 + (ctrl / Κ)^n) - δ * u[1]
    du[2] = κ * u[1] - γ * u[2]
    nothing
end

"Example 3.2 (Bistable gene circuit) - Murray"
function bistable_gene_circuit!(du, u, p, t)
    a, b, c, d, e, f = p
    du[1] = a / (1 + (u[2] / b)^c) - d * u[1]
    du[2] = e / (1 + (du[1] / f)^c) - d * u[2]
    nothing
end

"Example 3.4 (Negative autoregulation) - Murray"
function negative_autoregulation!(du, u, p, t)
    du[1] = 1 / (1 + u[2]) - u[1]
    du[2] = u[1] - u[2]
    nothing
end

""" Example 3.5 (Phosphorylation cycle) - Murray

Z + X -> Z + X*
Y + X* -> Y + X
X + X* = Xtot
Ytot = c

Notice that in this example there is a fixed amount of Xtot and Ytot.
"""
function phosphorylation_cycle!(out, du, u, p, t)
    Z, X, X2, Y = u
    Xtot, Ytot = p
    out[1] = - 0.04u[1]              + 1e4 * u[2] * u[3] - du[1]
    out[2] = + 0.04u[1] - 3e7 * u[2]^2 - 1e4 * u[2] * u[3] - du[2]
    out[3] = X + X2 - Xtot
end


text:
using DrWatson
@quickactivate "VelocityMTD"

using PyCall

pushfirst!(PyVector(pyimport("sys")."path"), srcdir())
pypsd = pyimport("PSD")
pyimport("importlib")."reload"(pypsd)

"""
Plot the PSD and its fit on `ax`. Input is velocity timeseries (or the notes
to calculate the timeseries from, by averaging velocity).

`kwargs` are propagated into `Python.psd`, while `minsub, tpq` are
"""
function psdfit!(ax, notes::Notes, middle = nothing;
    minsub = 12, tpq = 960, kwargs...)

    t, v = timeseries(notes, :velocity, mean, 0:(1//minsub):1)
    v = [ismissing(vi) ? 0 : vi for vi in v]
    psdfit!(ax, v, middle; minsub = minsub, tpq = tpq, kwargs...)
end

function psdfit!(ax, v::Vector, middle = nothing;
    minsub = 12, tpq = 960, kwargs...)

    mindur = tpq÷minsub

    f2b = (f) -> @. mindur/f/tpq
    b2f = (b) -> @. mindur/b/tpq

    f, psd = pypsd.psd(v; kwargs...)
    lo, up, N = pypsd.chisquare_ci(get(kwargs, :overlap_factor, 20) - 1)

    beats = f2b.(f[2:end])

    ax.loglog(beats, psd[2:end], c= "C0", lw = 1.5, basex=2)
    ax.fill_between(beats, psd[2:end] .* up ./ N, psd[2:end] .* lo ./ N, color = "C0", alpha = 0.25)

    if isnothing(middle)
        fs1, sl1, slope1 = pypsd.fit_powerlaw(f, psd, mindur, tpq; upper = maximum(beats), lower = minimum(beats))
        ax.plot(f2b.(fs1), sl1, label = "\$\\beta\$ = $(round(slope1; digits=2))", c = "C5")
        slope2 = 0.0
    else
        fs1, sl1, slope1 = pypsd.fit_powerlaw(f, psd, mindur, tpq; upper = middle, lower = minimum(beats))
        fs2, sl2, slope2 = pypsd.fit_powerlaw(f, psd, mindur, tpq; upper = maximum(beats), lower = middle)
        ax.plot(f2b.(fs1), sl1, label = "\$\\delta\$ = $(round(slope1; digits=2))", c = "C2")
        ax.plot(f2b.(fs2), sl2, label = "\$\\beta\$ = $(round(slope2; digits=2))", c = "C5")
    end

    ax.legend()
    ax.set_xlabel("period (beats)")
    ax.grid()
    return slope1, slope2
end

max_stars_repo_name: Datseris/VelocityMTD

max_issues_repo_name: Datseris/VelocityMTD

max_forks_repo_name: Datseris/VelocityMTD


text:
using Test
using Random
using Base.Order

using BinaryHeaps
using BinaryHeaps:
    walk,
    heapparent, heapleft, heapright,
    heapify!, bubble!,
    PermutationTracker,
    reduce_subtrees!, Reducer!

function isheap(v,order)
    n = length(v)

    # Nodes with two children
    for i = 1:heapparent(n-1)
        if lt(order, v[heapleft(i)], v[i]) ||
           lt(order, v[heapright(i)], v[i])
            return false
        end
    end

    # Node with one child
    for i = heapparent(n+1) : heapparent(n)
        if lt(order, v[heapleft(i)], v[i])
            return false
        end
    end

    return true
end

macro return_if_false(expr)
    return quote
        if !$(esc(expr))
            return false
        end
    end
end

function isreduction(v,r,reduce)
    n = length(v)
    if n == 0; return true; end

    # Check no-child nodes
    for i = n:-1:heapparent(n)+1
        @return_if_false r[i] == reduce(v[i])
    end

    # Check one-child node, if present
    for i = heapparent(n+1) : heapparent(n)
        @return_if_false r[i] == reduce(v[i],r[n])
    end

    # Reduce two-children nodes
    for i = heapparent(n)-1:-1:1
        @return_if_false r[i] == reduce(v[i], r[heapleft(i)], r[heapright(i)])
    end

    return true
end

include("walk.jl")
include("binary_heaps.jl")
include("permutations.jl")
include("reduction.jl")

include("MutableBinaryHeapsWithReduction.jl")


text:
p1 = codom(Open([1], PetriNet(1), [1]))
p2 = codom(Open([1,2], PetriNet(2), [1,2]))
pc1 = id(p1)
pc2 = id(p2)

@test dom(pc1 ⊗ pc2) == dom(pc1) ⊗ dom(pc2)
@test codom(pc1 ⊗ pc2) == codom(pc1) ⊗ codom(pc2)
@test dom(braid(p1,p2)) == p1 ⊗ p2
@test codom(braid(p1,p2)) == p2 ⊗ p1

@test dom(mcopy(p1)) == p1
@test codom(mcopy(p1)) == p1 ⊗ p1
@test dom(delete(p1)) == p1
@test codom(delete(p1)) == munit(OpenPetriNetOb)

max_stars_repo_name: AlgebraicJulia/AlgebraicPetri

max_issues_repo_name: AlgebraicJulia/AlgebraicPetri

max_forks_repo_name: AlgebraicJulia/AlgebraicPetri


text:
export RawAcquisitionData, EncodingCounters, AcquisitionHeader

"""
Encoding counters used in each Profile of a RawAcquisitionData object.
"""
Base.@kwdef struct EncodingCounters
  kspace_encode_step_1::Int16 = 0
  kspace_encode_step_2::Int16 = 0
  average::Int16 = 0
  slice::Int16 = 0
  contrast::Int16 = 0
  phase::Int16 = 0
  repetition::Int16 = 0
  set::Int16 = 0
  segment::Int16 = 0
  user::NTuple{8,Int16} = ntuple(i->Int16(0),8)
end

"""
Encoding counters used in each Profile of a RawAcquisitionData object.
"""
Base.@kwdef struct AcquisitionHeader
  version::Int16 = 0
  flags::Int64 = 0
  measurement_uid::Int32 = 0
  scan_counter::Int32 = 0
  acquisition_time_stamp::Int32 = 0
  physiology_time_stamp::NTuple{3,Int32} = ntuple(i->Int32(0),3)
  number_of_samples::Int16 = 0
  available_channels::Int16 = 0
  active_channels::Int16 = 0
  channel_mask::NTuple{16,Int64} = ntuple(i->Int64(0),16)
  discard_pre::Int16 = 0
  discard_post::Int16 = 0
  center_sample::Int16 = 0
  encoding_space_ref::Int16 = 0
  trajectory_dimensions::Int16 = 0
  sample_time_us::Float32 = 0.0
  position::NTuple{3,Float32} = ntuple(i->Float32(0),3)
  read_dir::NTuple{3,Float32} = ntuple(i->Float32(0),3)
  phase_dir::NTuple{3,Float32} = ntuple(i->Float32(0),3)
  slice_dir::NTuple{3,Float32} = ntuple(i->Float32(0),3)
  patient_table_position::NTuple{3,Float32} = ntuple(i->Float32(0),3)
  idx::EncodingCounters = EncodingCounters()
  user_int::NTuple{8,Int32} = ntuple(i->Int32(0),8)
  user_float::NTuple{8,Float32} = ntuple(i->Float32(0),8)
end

"""
Struct to describe data from one profile of a trajectory (from all coils).
"""
mutable struct Profile
  head::AcquisitionHeader
  traj::Array{Float32,2}
  data::Array{Complex{Float32},2}
end

"""
RawAcquisitionData object.

# Fields
* `params::Dict{String, Any}` - Dict containing the information of the XML header in ISMRMRD
* `profiles::Vector{Profile}` - Vector containing all the profiles of the acquisition

"""
mutable struct RawAcquisitionData
  params::Dict{String, Any}
  profiles::Vector{Profile}
end

"""
  trajectory(f::RawAcquisitionData; slice::Int=1, contrast::Int=1)

returns the `Trajectory` for given `slice` and `contrast` of a `RawAcquisitionData` object.
"""
function trajectory(f::RawAcquisitionData; slice::Int=1, contrast::Int=1)
  name = get(f.params, "trajectory","cartesian")
  if lowercase(name) == "cartesian"
    dim = ( maximum(encSteps2(f))==1 ? 2 : 3)
    if dim==3
      tr = trajectory("Cartesian3D", f.params["encodedSize"][2],
                                     f.params["encodedSize"][1],
                                     numSlices=f.params["encodedSize"][3])
    else
      tr = trajectory("Cartesian", f.params["encodedSize"][2],
                                   f.params["encodedSize"][1])
    end
  else
    name  = f.params["trajectory"]
    encSt1 = encSteps1(f)
    encSt2 = encSteps2(f)
    sl = slices(f)
    rep = repetitions(f)

    numSl = length(unique(sl))
    numRep = length(unique(rep))
    numEncSt1 = length(unique(encSt1))
    numEncSt2 = length(unique(encSt2))

    # assume constant number of samplings per profile
    numSampPerProfile = size(f.profiles[1].data,1)
    numChan = size(f.profiles[1].data,2)
    D = Int(f.profiles[1].head.trajectory_dimensions)

    # remove data that should be discarded
    i1 = f.profiles[1].head.discard_pre + 1
    i2 = numSampPerProfile - f.profiles[1].head.discard_post

    traj = zeros(Float32, D, length(i1:i2), numEncSt1, numEncSt2, numSl, numRep)
    times = zeros(Float32, length(i1:i2), numEncSt1, numEncSt2, numSl, numRep)

    for l=1:length(f.profiles)
      if f.profiles[l].head.idx.slice+1 != slice || f.profiles[l].head.idx.contrast+1 != contrast
        continue
      end
      traj[:, :, encSt1[l], encSt2[l], sl[l]-minimum(sl)+1, rep[l]-minimum(rep)+1] .= f.profiles[l].traj[:,i1:i2]
      dt = f.profiles[l].head.sample_time_us*1e-6
      times[:, encSt1[l], encSt2[l], sl[l]-minimum(sl)+1, rep[l]-minimum(rep)+1] .= 0:dt:(length(i1:i2)-1)*dt
    end

    traj_ = reshape(traj[:,:,:,:,1,1], D, :)
    # tr = Trajectory(traj_, size(traj_,3), size(traj_,2), circular=true)
    tr = Trajectory(traj_, size(traj_,3), size(traj_,2), circular=true, times=vec(times[:,:,:,1,1]))
  end
  return tr
end

function sequence(f::RawAcquisitionData)
  # TODO
end

"""
    numChannels(f::RawAcquisitionData)

returns the number of channels in a `RawAcquisitionData` object.
"""
function numChannels(f::RawAcquisitionData)
  return f.profiles[1].head.active_channels
end

encSteps1(f::RawAcquisitionData) =
   [f.profiles[l].head.idx.kspace_encode_step_1+1 for l=1:length(f.profiles)]
encSteps2(f::RawAcquisitionData) =
   [f.profiles[l].head.idx.kspace_encode_step_2+1 for l=1:length(f.profiles)]
slices(f::RawAcquisitionData) =
   [f.profiles[l].head.idx.slice+1 for l=1:length(f.profiles)]
repetitions(f::RawAcquisitionData) =
   [f.profiles[l].head.idx.repetition+1 for l=1:length(f.profiles)]
contrasts(f::RawAcquisitionData) =
  [f.profiles[l].head.idx.contrast+1 for l=1:length(f.profiles)]

  """
    subsampleIndices(f::RawAcquisitionData; slice::Int=1, contrast::Int=1)

  returns the sampled indices for a given `slice` and `contrast` in a `RawAcquisitionData` object.
  """
function subsampleIndices(f::RawAcquisitionData; slice::Int=1, contrast::Int=1, estimateProfileCenter::Bool=false)
  idx = Int64[]
  encSt1 = encSteps1(f)
  encSt2 = encSteps2(f)
  numEncSamp, numProf, numSl = f.params["encodedSize"]
  for i=1:length(f.profiles)
    # only consider data for the specified slice, contrast and repetition
    if f.profiles[i].head.idx.slice+1 != slice || f.profiles[i].head.idx.contrast+1 != contrast || f.profiles[i].head.idx.repetition != 0
      continue
    end
    numSamp = size(f.profiles[i].data,1)  # number of measured points
    if estimateProfileCenter
      center_sample = div(numSamp,2)+1
    else
      center_sample = f.profiles[i].head.center_sample+1 # idx of the center trajectory in the measurement
    end
    tr_center_idx = div(numEncSamp,2)+1 # center of the encoded trajectory

    # start and end indices of sampled profile
    i1 = tr_center_idx - center_sample + 1 + f.profiles[i].head.discard_pre
    i2 = tr_center_idx - center_sample + numSamp - f.profiles[i].head.discard_post
    # convert to linear index
    lineIdx = collect(i1:i2) .+ numSamp*((encSt2[i]-1)*numProf + (encSt1[i]-1))
    append!(idx, lineIdx)
  end

  return unique(idx)
end

"""
    rawdata(f::RawAcquisitionData)

returns the rawdata contained `RawAcquisitionData` object.
The output is an `Array{Matrix{ComplexF64},3}`, which can be stored in a `AcquisitionData` object.
"""
function rawdata(f::RawAcquisitionData; slice::Int=1, contrast::Int=1, repetition::Int=1)

  # find profiles corresponding to the given slice contrast and repitiion
  sl = slices(f)
  rep = repetitions(f)
  contr = contrasts(f)
  idx_sl = findall(x->x==slice,sl)
  idx_contr = findall(x->x==contrast,contr)
  idx_rep = findall(x->x==repetition,rep)
  idx = intersect(idx_sl,idx_contr,idx_rep)

  # number of unique combination of encoding statuses
  encSt1 = encSteps1(f)[idx]
  encSt2 = encSteps2(f)[idx]
  # find first occurence of each combination of encoding statuses (as `unique` returns)
  idx_unique = uniqueidx(hcat(encSt1,encSt2))
  numEncSt = length(idx_unique)

  # allocate space for k-space data
  # assume the same number of samples for all profiles
  numSampPerProfile, numChan = size(f.profiles[idx[1]].data)
  numSampPerProfile -= (f.profiles[idx[1]].head.discard_pre+f.profiles[idx[1]].head.discard_post)
  kdata = zeros(ComplexF64, numSampPerProfile, numEncSt, numChan)

  # store one profile (kspace data) for each unique encoding status
  cnt = 1
  for l in idx[idx_unique] # TODO: set of profiles (with unique encoding status)
    i1 = f.profiles[l].head.discard_pre + 1
    i2 = i1+numSampPerProfile-1
    kdata[:,cnt,:] .= ComplexF64.(f.profiles[l].data[i1:i2, :])
    cnt += 1
  end

  return reshape(kdata, :, numChan)
end

"""
    AcquisitionData(f::RawAcquisitionData; estimateProfileCenter::Bool=false)

converts `RawAcquisitionData` into the equivalent `AcquisitionData` object.
"""
function AcquisitionData(f::RawAcquisitionData; estimateProfileCenter::Bool=false)
  contrs = sort(unique(contrasts(f)))
  sls = sort(unique(slices(f)))
  reps = sort(unique(repetitions(f)))
  numContr = length(unique(contrasts(f)))
  numSl = length(unique(slices(f)))
  numRep = length(unique(repetitions(f)))
  # tr = [trajectory(f,contrast=contr) for contr=1:numContr]
  # subsampleIdx = [subsampleIndices(f,contrast=contr,estimateProfileCenter=estimateProfileCenter) for contr=1:numContr]
  # kdata = [rawdata(f; contrast=i, slice=j, repetition=k) for i=1:numContr, j=1:numSl, k=1:numRep]
  tr = [trajectory(f,contrast=contr) for contr=contrs]
  subsampleIdx = [subsampleIndices(f,contrast=contr,estimateProfileCenter=estimateProfileCenter) for contr=contrs]
  kdata = [rawdata(f; contrast=i, slice=j, repetition=k) for i=contrs, j=sls, k=reps]

  return AcquisitionData(tr, kdata,
                          idx=subsampleIdx,
                          encodingSize=collect(f.params["encodedSize"]),
                          fov = collect(f.params["encodedFOV"]) )
end

"""
    RawAcquisitionData(acqData::AcquisitionData)

converts `acqData` into the equivalent `RawAcquisitionData` object.
"""
function RawAcquisitionData(acqData::AcquisitionData)
  # XML header
  params = minimalHeader(Tuple(acqData.encodingSize),Tuple(acqData.fov),tr_name=string(trajectory(acqData,1)))
  # acquisition counter
  counter = 1
  # profiles
  profiles = Vector{Profile}()
  for rep = 1:numRepititions(acqData)
    for slice=1:numSlices(acqData)
      for contr = 1:numContrasts(acqData)
        tr = trajectory(acqData,contr)
        profIdx = profileIdx(acqData,contr)
        numSamp = numSamplingPerProfile(tr)
        numProf = div(length(acqData.subsampleIndices[contr]),numSamp)
          for prof_tr = 1:length(profIdx)
            for slice_tr = 1:numSlices(tr)
              head = AcquisitionHeader(acqData,rep,slice,slice_tr,contr,profIdx[prof_tr],counter)
              nodes = profileNodes(tr,prof_tr,slice_tr)
              if dims(tr)==2
                data = profileData(acqData,contr,slice,rep,prof_tr)
              else
                data = profileData(acqData,contr,slice_tr,rep,prof_tr)
              end
              profile = Profile( head, Float32.(nodes), Complex{Float32}.(data) )
              push!( profiles, profile )
              counter += 1
            end
          end
      end
    end
  end

  return RawAcquisitionData(params,profiles)
end

function AcquisitionHeader(acqData::AcquisitionData, rep::Int, slice::Int, slice_tr::Int
                            , contr::Int, prof_tr::Int, counter::Int
                            ; phase::Int=1, set::Int=1, segment::Int=1, repetition::Int=1
                            , user::NTuple{8,Int16} = ntuple(i->Int16(0),8), kargs...)
  numChan = numChannels(acqData)
  tr = trajectory(acqData,contr)
  numSamples = numSamplingPerProfile(tr)
  isCartesian(tr) ? center_sample=div(numSamples,2) : center_sample=0 # needs to be fixed for partial Fourier,etc
  tr_dims = dims(tr)
  sample_time_us = acqTimePerProfile(tr)/numSamples*1.e6

  idx = EncodingCounters( Int16(prof_tr-1)
                        , Int16(slice_tr-1)
                        , Int16(rep-1)
                        , Int16(slice-1)
                        , Int16(contr-1)
                        , Int16(phase-1)
                        , Int16(repetition-1)
                        , Int16(set-1)
                        , Int16(segment-1)
                        , user )

  return AcquisitionHeader(; scan_counter=Int32(counter)
                          , number_of_samples=Int16(numSamples)
                          , available_channels=Int16(numChan)
                          , active_channels=Int16(numChan)
                          , center_sample=Int16(center_sample)
                          , trajectory_dimensions=Int16(tr_dims)
                          , sample_time_us=Float32(sample_time_us)
                          , read_dir=Float32.((1,0,0))
                          , phase_dir=Float32.((0,1,0))
                          , slice_dir=Float32.((0,0,1))
                          , idx=idx )
end

# get indices of the first occurence of unique elements in the matrix x
function uniqueidx(x::Matrix{T}) where T
  uniqueset = Set{Vector{T}}()
  ex = eachindex(x[:,1])
  idxs = Vector{eltype(ex)}()
  for i in ex
      xi = x[i,:]
      if !(xi in uniqueset)
          push!(idxs, i)
          push!(uniqueset, xi)
      end
  end
  idxs
end


text:
module SKLearners

using PyCall

# standard included modules
using DataFrames
using Random
using ..AbsTypes
using ..Utils

import ..AbsTypes: fit, fit!, transform, transform!
export fit, fit!, transform, transform!
export SKLearner, sklearners

const learner_dict = Dict{String,PyObject}()
const ENS   = PyNULL()
const LM    = PyNULL()
const DA    = PyNULL()
const NN    = PyNULL()
const SVM   = PyNULL()
const TREE  = PyNULL()
const ANN   = PyNULL()
const GP    = PyNULL()
const KR    = PyNULL()
const NB    = PyNULL()
const ISO   = PyNULL()

function __init__()
   copy!(ENS , pyimport_conda("sklearn.ensemble","scikit-learn"))
   copy!(LM  , pyimport_conda("sklearn.linear_model","scikit-learn"))
   copy!(DA  , pyimport_conda("sklearn.discriminant_analysis","scikit-learn"))
   copy!(NN  , pyimport_conda("sklearn.neighbors","scikit-learn"))
   copy!(SVM , pyimport_conda("sklearn.svm","scikit-learn"))
   copy!(TREE, pyimport_conda("sklearn.tree","scikit-learn"))
   copy!(ANN , pyimport_conda("sklearn.neural_network","scikit-learn"))
   copy!(GP  , pyimport_conda("sklearn.gaussian_process","scikit-learn"))
   copy!(KR  , pyimport_conda("sklearn.kernel_ridge","scikit-learn"))
   copy!(NB  , pyimport_conda("sklearn.naive_bayes","scikit-learn"))
   copy!(ISO , pyimport_conda("sklearn.isotonic","scikit-learn"))

   # Available scikit-learn learners.
   learner_dict["AdaBoostClassifier"]             = ENS
   learner_dict["BaggingClassifier"]              = ENS
   learner_dict["ExtraTreesClassifier"]           = ENS
   learner_dict["VotingClassifier"]               = ENS
   learner_dict["GradientBoostingClassifier"]     = ENS
   learner_dict["RandomForestClassifier"]         = ENS
   learner_dict["QuadraticDiscriminantAnalysis"] = DA
   learner_dict["LinearDiscriminantAnalysis"]     = DA
   learner_dict["LogisticRegression"]             = LM
   learner_dict["PassiveAggressiveClassifier"]    = LM
   learner_dict["RidgeClassifier"]                = LM
   learner_dict["RidgeClassifierCV"]              = LM
   learner_dict["SGDClassifier"]                  = LM
   learner_dict["KNeighborsClassifier"]           = NN
   learner_dict["RadiusNeighborsClassifier"]      = NN
   learner_dict["NearestCentroid"]                = NN
   learner_dict["SVC"]                            = SVM
   learner_dict["LinearSVC"]                      = SVM
   learner_dict["NuSVC"]                          = SVM
   learner_dict["MLPClassifier"]                  = ANN
   learner_dict["GaussianProcessClassifier"]      = GP
   learner_dict["DecisionTreeClassifier"]         = TREE
   learner_dict["GaussianNB"]                     = NB
   learner_dict["MultinomialNB"]                  = NB
   learner_dict["ComplementNB"]                   = NB
   learner_dict["BernoulliNB"]                    = NB
   learner_dict["SVR"]                            = SVM
   learner_dict["Ridge"]                          = LM
   learner_dict["RidgeCV"]                        = LM
   learner_dict["Lasso"]                          = LM
   learner_dict["ElasticNet"]                     = LM
   learner_dict["Lars"]                           = LM
   learner_dict["LassoLars"]                      = LM
   learner_dict["OrthogonalMatchingPursuit"]      = LM
   learner_dict["BayesianRidge"]                  = LM
   learner_dict["ARDRegression"]                  = LM
   learner_dict["SGDRegressor"]                   = LM
   learner_dict["PassiveAggressiveRegressor"]     = LM
   learner_dict["KernelRidge"]                    = KR
   learner_dict["KNeighborsRegressor"]            = NN
   learner_dict["RadiusNeighborsRegressor"]       = NN
   learner_dict["GaussianProcessRegressor"]       = GP
   learner_dict["DecisionTreeRegressor"]          = TREE
   learner_dict["RandomForestRegressor"]          = ENS
   learner_dict["ExtraTreesRegressor"]            = ENS
   learner_dict["AdaBoostRegressor"]              = ENS
   learner_dict["GradientBoostingRegressor"]      = ENS
   learner_dict["IsotonicRegression"]             = ISO
   learner_dict["MLPRegressor"]                   = ANN
end

"""
    SKLearner(learner::String, args::Dict=Dict())

A Scikitlearn wrapper to load the different machine learning models.
Invoking `sklearners()` will list the available learners. Please
consult Scikitlearn documentation for arguments to pass.

Implements `fit!` and `transform!`.
"""
mutable struct SKLearner <: Learner
   name::String
   model::Dict{Symbol,Any}

   function SKLearner(args=Dict{Symbol,Any}())
      default_args=Dict{Symbol,Any}(
         :name => "sklearner",
         :output => :class,
         :learner => "LinearSVC",
         :impl_args => Dict{Symbol,Any}()
      )
      cargs = nested_dict_merge(default_args, args)
      cargs[:name] = cargs[:name]*"_"*randstring(3)
      skl = cargs[:learner]
      if !(skl in keys(learner_dict))
         println("$skl is not supported.")
         println()
         sklearners()
         error("Argument keyword error")
      end
      new(cargs[:name],cargs)
   end
end

function SKLearner(learner::String, args::Dict)
   SKLearner(Dict(:learner => learner,:name=>learner, args...))
end

function SKLearner(learner::String; args...)
   SKLearner(Dict(:learner => learner,:name=>learner,:impl_args=>Dict(pairs(args))))
end

function (skl::SKLearner)(;objargs...)
   skl.model[:impl_args] = Dict(pairs(objargs))
   skname = skl.model[:learner]
   skobj = getproperty(learner_dict[skname],skname)
   newskobj = skobj(;objargs...)
   skl.model[:sklearner] = newskobj
   return skl
end

"""
    function sklearners()

List the available scikitlearn machine learners.
"""
function sklearners()
  learners = keys(learner_dict) |> collect |> x-> sort(x,lt=(x,y)->lowercase(x)<lowercase(y))
  println("syntax: SKLearner(name::String, args::Dict=Dict())")
  println("where 'name' can be one of:")
  println()
  [print(learner," ") for learner in learners]
  println()
  println()
  println("and 'args' are the corresponding learner's initial parameters.")
  println("Note: Consult Scikitlearn's online help for more details about the learner's arguments.")
end

function fit!(skl::SKLearner, xx::DataFrame, y::Vector)::Nothing
  x = xx |> Array
  impl_args  = copy(skl.model[:impl_args])
  learner    = skl.model[:learner]
  py_learner = getproperty(learner_dict[learner],learner)

  # Assign CombineML-specific defaults if required
  if learner == "RadiusNeighborsClassifier"
    if get(impl_args, :outlier_label, nothing) == nothing
      impl_options[:outlier_label] = labels[rand(1:size(labels, 1))]
    end
  end

  # Train
  modelobj = py_learner(;impl_args...)
  modelobj.fit(x,y)
  skl.model[:sklearner] = modelobj
  skl.model[:impl_args] = impl_args
  return nothing
end

function fit(skl::SKLearner, xx::DataFrame, y::Vector)::SKLearner
   fit!(skl,xx,y)
   return deepcopy(skl)
end

function transform!(skl::SKLearner, xx::DataFrame)::Vector
	x = deepcopy(xx) |> Array
  #return collect(skl.model[:predict](x))
  sklearner = skl.model[:sklearner]
  return collect(sklearner.predict(x))
end

transform(skl::SKLearner, xx::DataFrame)::Vector = transform!(skl,xx)

end


text:
export opt_linear_fit!

"""
    opt_linear_fit!(
        graph,
        objfun,
        discr,
        linear_cref;
        input = :A,
        errtype = :abserr,
        linlsqr = :backslash,
        droptol = 0,
    )

Linear fitting of a `graph` of the form

    c_1 g_1(x) + c_2 g_2(x) + … + c_n g_n(x)

to the values of `objfun`, in the points `discr`. Reference to the coefficients
`c_1,…,c_n` should be given `linear_cref`.

The variable `graph` is modified during the iterations and the function has no
return value.

See [`opt_gauss_newton!`](@ref) for a description the kwarg `errtype` and
`input`, and [`solve_linlsqr!`](@ref) for the kwargs `linlsqr and `droptol`.
"""
function opt_linear_fit!(
    graph,
    objfun,
    discr,
    linear_cref;
    input = :A,
    errtype = :abserr,
    linlsqr = :backslash,
    droptol = 0,
)
    objfun_vals = objfun.(discr)

    vals = init_vals_eval_graph!(graph, discr, nothing, input)
    eval_graph(graph, discr, vals = vals, input = input)

    n = length(linear_cref)
    T = eltype(valtype(vals))
    A = ones(T, size(discr, 1), n)
    for k = 1:n
        parent = graph.parents[linear_cref[k][1]][linear_cref[k][2]]
        A[:, k] = vals[parent]
    end

    adjust_for_errtype!(A, objfun_vals, objfun_vals, errtype)

    c = solve_linlsqr!(A, objfun_vals, linlsqr, droptol)
    set_coeffs!(graph, c, linear_cref)

    return nothing
end


text:
#=
Masking operations
==================

Copyright (c) 2017 Ben J. Ward & Luis Yanes

This Source Code Form is subject to the terms of the Mozilla Public
License, v. 2.0. If a copy of the MPL was not distributed with this
file, You can obtain one at http://mozilla.org/MPL/2.0/.
=#

"""
    nibble_mask{T<:Unsigned}(value::T, x::T)

Create a mask for the nibbles (aligned 4 bit segments) in an unsigned integer
`x` that filter nibbles matching the corresponding nibble in `value`.
"""
@inline function nibble_mask(value::T, x::T) where {T<:Unsigned}
    # XOR with the desired values. So matching nibbles will be 0000.
    x = x ⊻ value
    # Horizontally OR the nibbles.
    x |= (x >>> 1)
    x |= (x >>> 2)
    # AND removes junk, we then widen x by multiplication and return
    # the inverse.
    x &= repeatpattern(T, 0x11)
    x *= 0x0F
    return ~x
end

"""
    mask{T<:Unsigned}(::Type{T}, n::Integer)

Creates a bit mask for given number of bits `n`.

The mask starts from the least significant bit, and end at bit `n`.

e.g:

```jldoctest
julia> Twiddle.mask(UInt64, 8)
0x00000000000000ff
```
"""
@inline mask(::Type{T}, n::Integer) where {T<:Unsigned} = (T(1) << n) - 0x1
@inline mask(n::Integer) = mask(UInt64, n)


text:
"""
    nlsv_obj!(obs::WSVarLmmObs, β, τ, Lγ, needgrad::Bool)
    nlsv_obj!(m::WSVarLmmModel; needgrad::Bool)

Evaluate the nonlinear least squares (NLS) criterion for variance estimation at
the given data and parameter values. Gradient is calculated if `needgrad=true`.
Expected Hessian is calculated if `needhess=true`.  If `updateres=true`, update
mean level residuals first. If `m.iswtnls[1]=true`, evaluate the weighted
`nlsv_obj!()` function. `update_wtmat!(m)` should be called to update the
weight matrix componentsprior to using the weighted version of `nlsv_obj!()`.
"""
function nlsv_obj!(
    obs       :: WSVarLmmObs{T},
    β         :: Vector{T},
    τ         :: Vector{T},
    Lγ        :: Matrix{T},  # must be lower triangular
              :: Val{false}, # un-weighted fitting
    needgrad  :: Bool = true,
    needhess  :: Bool = true,
    updateres :: Bool = false
    ) where T <: BlasReal
    (q, n) = size(obs.Zt)
    l = size(obs.Wt, 1)
    ###########
    # objective
    ###########
    # update the residual vector ri = y_i - Xi β
    updateres && update_res!(obs, β)
    # obs.storage_qq = (Z' * Z) * L
    copyto!(obs.storage_qq, obs.ztz)
    BLAS.trmm!('R', 'L', 'N', 'N', T(1), Lγ, obs.storage_qq)
    # ∇Lγ = (Z' * Z) * L * L' * (Z' * Z) for now, needed for gradient later
    needgrad && BLAS.syrk!('U', 'N', T(1), obs.storage_qq, T(0), obs.∇Lγ)
    # HLγLγ = C'(L'Z'Z ⊗ Z'ZL)KC for now, needed for hessian later
    needhess && Ct_At_kron_A_KC!(fill!(obs.HLγLγ, 0), obs.storage_qq)
    # obs.storage_qq = L' * (Z' * Z) * L
    BLAS.trmm!('L', 'L', 'T', 'N', T(1), Lγ, obs.storage_qq)
    # storage_qn = L' * Z'
    copyto!(obs.storage_qn, obs.Zt)
    BLAS.trmm!('L', 'L', 'T', 'N', T(1), Lγ, obs.storage_qn)
    # storage_q◺n = Cq' * (L'Z' ⊙ Z'), needed for hessian later
    needhess && Ct_A_kr_B!(fill!(obs.storage_q◺n, 0), obs.storage_qn, obs.Zt)
    # storage_q1 = L' * Z' * res
    copyto!(obs.storage_q1, obs.ztres)
    BLAS.trmv!('L', 'T', 'N', Lγ, obs.storage_q1)
    # update W * τ
    mul!(obs.expwτ, transpose(obs.Wt), τ)
    obs.obj[1]  = (1//2) * (abs2(norm(obs.storage_qq)) + abs2(obs.resnrm2[1]))
    obs.obj[1] -= abs2(norm(obs.storage_q1))
    map!(exp, obs.expwτ, obs.expwτ)
    @inbounds for j in 1:n
        obs.obj[1] += (1//2) * abs2(obs.expwτ[j])
        obs.zlltzt_dg[j] = 0
        for i in 1:q
            obs.zlltzt_dg[j] += abs2(obs.storage_qn[i, j])
        end
        obs.obj[1] += obs.expwτ[j] * (obs.zlltzt_dg[j] - obs.res2[j])
    end
    ###########
    # gradient
    ###########
    if needgrad
        # wrt τ
        @inbounds for j in 1:n
            Rjj = obs.res2[j] - obs.expwτ[j] - obs.zlltzt_dg[j]
            obs.storage_n1[j] = Rjj * obs.expwτ[j]
        end
        BLAS.gemv!('N', T(-1), obs.Wt, obs.storage_n1, T(0), obs.∇τ)
        # wrt Lγ
        # ∇Lγ = (Z' * R * Z) * Lγ
        # obs.storage_qn = obs.Zt * Diagonal(sqrt.(obs.expwτ))
        @inbounds for j in 1:n
            sqrtej = sqrt(obs.expwτ[j])
            for i in 1:q
                obs.storage_qn[i, j] = sqrtej * obs.Zt[i, j]
            end
        end
        # ∇Lγ = (Z' * Z) * L * L' * (Z' * Z) was computed earlier
        # ∇Lγ += storage_qn * storage_qn' - ztres * ztres'
        BLAS.syrk!('U', 'N',  T(1), obs.storage_qn, T(1), obs.∇Lγ)
        BLAS.syrk!('U', 'N', T(-1),      obs.ztres, T(1), obs.∇Lγ)
        copytri!(obs.∇Lγ, 'U')
        # so far ∇Lγ holds ∇Σγ, now ∇Lγ = ∇Σγ * Lγ
        # obs.∇Lγ = obs.∇Σγ, collect all then multiply by Lγ at model level for ∇Lγ.
        # BLAS.trmm!('R', 'L', 'N', 'N', T(2), Lγ, obs.∇Lγ)
    end
    ###########
    # hessian
    ###########
    if needhess
        # Hττ = W' * Diagonal(expwτ.^2) * W
        # storage_ln = W' * Diagonal(expwτ)
        @inbounds for j in 1:n
            ej = obs.expwτ[j]
            for i in 1:l
                obs.storage_ln[i, j] = ej * obs.Wt[i, j]
            end
        end
        BLAS.syrk!('U', 'N', T(1), obs.storage_ln, T(0), obs.Hττ)
        copytri!(obs.Hττ, 'U')
        # HτLγ = 2 W' * Diagonal(expwτ) * (L'Z' ⊙ Z')' * Cq
        # storage_ln = W' * Diagonal(expwτ) was computed above
        # storage_q◺n = Cq' * (L'Z' ⊙ Z') was computed earlier
        BLAS.gemm!('N', 'T', T(2), obs.storage_ln, obs.storage_q◺n, T(0), obs.HτLγ)
        # HLγLγ = 2 [ C'(L'Z'ZL ⊗ Z'Z)C + C'(L'Z'Z ⊗ Z'ZL)KC ]
        # HLγLγ = C'(L'Z'Z ⊗ Z'ZL)KC was calcualted earlier
        Ct_A_kron_B_C!(obs.HLγLγ, obs.storage_qq, obs.ztz)
        lmul!(2, obs.HLγLγ)
    end
    obs.obj[1]
end

function nlsv_obj!(
    obs       :: WSVarLmmObs{T},
    β         :: Vector{T},
    τ         :: Vector{T},
    Lγ        :: Matrix{T}, # must be lower triangular
              :: Val{true}, # weighted fitting
    needgrad  :: Bool = true,
    needhess  :: Bool = true,
    updateres :: Bool = false,
    ) where T <: BlasReal
    (q, n) = size(obs.Zt)
    l = size(obs.Wt, 1)
    # update the residual vector ri = y_i - Xi β
    updateres && update_res!(obs, β)
    # Precompute and fill necessary objects
    if needgrad
        fill!(obs.∇τ, 0)
    end
    if needhess
        fill!(obs.Hττ, 0)
        fill!(obs.HτLγ, 0)
        fill!(obs.storage_q◺n, 0)
    end
    # obs.expwτ = exp.(Wτ)
    mul!(obs.expwτ, transpose(obs.Wt), τ)
    # for ∇Lγ need storage_n1 = diag(exp.(0.5Wτ))
    obs.expwτ .= exp.(obs.expwτ)
    needgrad && (obs.storage_n1 .= sqrt.(obs.expwτ))
    # terms to compute and store
    mul!(obs.Lt_Zt_Dinv_r, transpose(Lγ), obs.Zt_Dinv_r)
    mul!(obs.Dinv_Z_L, transpose(obs.Zt_Dinv), Lγ)
    mul!(obs.UUt_Z_L, transpose(obs.Zt_UUt), Lγ)
    # Lt_Zt_Dinv_Z_L =  Zt_Dinv_Z * L
    mul!(obs.Lt_Zt_Dinv_Z_L, obs.Zt_Dinv_Z, Lγ)
    BLAS.trmm!('L', 'L', 'T', 'N', T(1), Lγ, obs.Lt_Zt_Dinv_Z_L)
    # Lt_Zt_UUt_Z_L =  Zt_UUt_Z_L
    mul!(obs.Lt_Zt_UUt_Z_L, obs.Zt_UUt_Z, Lγ)
    BLAS.trmm!('L', 'L', 'T', 'N', T(1), Lγ, obs.Lt_Zt_UUt_Z_L)
    # Evalute objective function and precompute diag(Vinv*R*Vinv) for ∇τ
    #############
    # objective #
    #############
    # objective function sums
    obs.obj[1]  = zero(T)
    obs.obj[1] += abs2(obs.rt_Dinv_r[1]) #1
    obs.obj[1] -= 2 * obs.rt_UUt_r[1] * obs.rt_Dinv_r[1] #2
    obs.obj[1] += abs2(obs.rt_UUt_r[1]) #7
    # calculate Ut * D * U
    fill!(obs.Ut_D_U, 0)
    @inbounds for k in 1:n, j in 1:q, i in j:q
        obs.Ut_D_U[i, j] += obs.Ut[j, k] * obs.Ut[i, k] * obs.expwτ[k]
    end
    copytri!(obs.Ut_D_U, 'L')
    if needgrad
        # storage_qn = Ut_D_U_Ut
        mul!(obs.storage_qn, obs.Ut_D_U, obs.Ut)
        fill!(obs.diagDVRV, 0)
        @inbounds for j in 1:n, i in 1:q  #avx speeds up
            obs.diagDVRV[j] -= obs.storage_qn[i, j] * obs.Ut[i, j] * obs.expwτ[j] #τ 8
        end
    end
    @inbounds @simd for j in 1:n
        obs.obj[1] += (-2 * abs2(obs.Dinv_r[j]) #3
            + 4 * obs.rt_UUt[j] * obs.Dinv_r[j]  #4
            - 2 * abs2(obs.rt_UUt[j]) # 8
            + abs2(obs.Dinv[j]) * obs.expwτ[j] #11
            - 2 * obs.diagUUt_Dinv[j] * obs.expwτ[j]) * obs.expwτ[j] #12
    end
    if needgrad
        @inbounds @simd for j in 1:n
            obs.diagDVRV[j] += (abs2(obs.Dinv_r[j])   #τ 1
                - 2 * obs.Dinv_r[j] * obs.rt_UUt[j] #τ 2
                - abs2(obs.Dinv[j]) * obs.expwτ[j]  #τ 3
                + 2 * obs.diagUUt_Dinv[j] * obs.expwτ[j]  #τ 4
                + abs2(obs.rt_UUt[j])) * obs.expwτ[j]  #τ 7
        end
    end
    @inbounds for j in 1:q  #j-i looping for memory access
        for i in 1:n
            obs.obj[1] += 2 * (abs2(obs.Dinv_Z_L[i, j]) # 13
                - obs.UUt_Z_L[i, j] * obs.Dinv_Z_L[i, j] # 14
                - obs.Dinv_Z_L[i, j] * obs.UUt_Z_L[i, j] # 16
                + abs2(obs.UUt_Z_L[i, j])) * obs.expwτ[i]# 17
        end
    end
    if needgrad
        @inbounds for j in 1:q
            for i in 1:n
                obs.diagDVRV[i] += (-abs2(obs.Dinv_Z_L[i, j]) #τ 5
                    + 2 * obs.Dinv_Z_L[i, j] * obs.UUt_Z_L[i, j] #τ 6
                    - abs2(obs.UUt_Z_L[i, j])) * obs.expwτ[i] #τ 9
            end
        end
    end
    #obs.storage_qq = L' * Z' * UU' * rr' * Dinv * Z * L
    copyto!(obs.storage_qq, obs.Zt_UUt_rrt_Dinv_Z)
    BLAS.trmm!('L', 'L', 'T', 'N', one(T), Lγ, obs.storage_qq)
    BLAS.trmm!('R', 'L', 'N', 'N', one(T), Lγ, obs.storage_qq) #for #9
    #use ∇Lγ as temporary storage
    #obs.∇Lγ = L' Zt_UUt_rrt_UUt_Z * L
    copyto!(obs.∇Lγ, obs.Zt_UUt_rrt_UUt_Z)
    BLAS.trmm!('L', 'L', 'T', 'N', T(1), Lγ, obs.∇Lγ)
    BLAS.trmm!('R', 'L', 'N', 'N', T(1), Lγ, obs.∇Lγ) #for #10
    @inbounds for j in 1:q
        obs.obj[1] -= 2 * abs2(obs.Lt_Zt_Dinv_r[j]) #5
        obs.obj[1] += 2 * obs.storage_qq[j, j] #9
        obs.obj[1] -= 2 * obs.∇Lγ[j, j] # 10
        for i in 1:q
            obs.obj[1] += obs.Ut_D_U[i, j]^2 #15
            obs.obj[1] += 2 * obs.rt_UUt_Z[i] * Lγ[i, j] * obs.Lt_Zt_Dinv_r[j] #6
            obs.obj[1] += abs2(obs.Lt_Zt_Dinv_Z_L[i, j]) #18
            obs.obj[1] -= 2 * obs.Lt_Zt_Dinv_Z_L[i, j] * obs.Lt_Zt_UUt_Z_L[i,j] #19
            obs.obj[1] += abs2(obs.Lt_Zt_UUt_Z_L[i, j]) #20
        end
    end
    obs.obj[1] *= (1//2)
    ############
    # Gradient #
    ############
    if needgrad
        #wrt τ
        #∇τ = -W' * diag(D * Vinv * R * Vinv)
        BLAS.gemv!('N', T(-1), obs.Wt, obs.diagDVRV, T(0), obs.∇τ)
        #wrt Lγ
        # ∇Lγ = -2(Z' * Vinv * R *  Vinv * Z * L
        # obs.storage_qq = (Z' * Vinv * Z) * L
        copyto!(obs.storage_qq, obs.Zt_Vinv_Z)
        BLAS.trmm!('R', 'L', 'N', 'N', T(1), Lγ, obs.storage_qq)
        #∇Lγ = Z' * Vinv * Z * L * L' * Z' * Vinv * Z
        BLAS.syrk!('U', 'N', T(1), obs.storage_qq, T(0), obs.∇Lγ)
        # obs.storage_qn = Zt_Vinv * sqrt.(D)
        # storage_n1 = diag(exp^{1/2 Wt})
        @inbounds for j in 1:n, i in 1:q
            obs.storage_qn[i, j] = obs.storage_n1[j] * obs.Zt_Vinv[i, j]
        end
        # obs.∇Lγ += Zt_Vinv_D_Vinv_Z
        BLAS.syrk!('U', 'N',  T(1), obs.storage_qn, T(1), obs.∇Lγ)
        # obs.∇Lγ += Zt_Vinv_rrt_Vinv_Z
        BLAS.syr!('U', T(-1), obs.Zt_Vinv_r, obs.∇Lγ)
        copytri!(obs.∇Lγ, 'U')
        # obs.∇Lγ = obs.∇Σγ, collect all then multiply by Lγ at model level for ∇Lγ.
    end
    ###########
    # Hessian #
    ###########
    if needhess
        #wrt τ
        # Hττ = W' * D * Vinv .* Vinv * D  * W
        # storage_ln = Wt * D
        @inbounds @simd for j in 1:n
            for i in 1:l
                obs.storage_ln[i, j] = obs.Wt[i, j] * obs.expwτ[j]
            end
        end
        @inbounds for j in 1:n, i in 1:l
            obs.Wt_D_Dinv[i, j] = obs.Dinv[j] * obs.storage_ln[i, j]
            obs.Wt_D_sqrtdiagDinv_UUt[i, j] = obs.storage_ln[i, j] * obs.sqrtDinv_UUt[j]
        end
        BLAS.syrk!('U', 'N', T(1), obs.Wt_D_Dinv, T(0), obs.Hττ) #first term
        BLAS.syrk!('U', 'N', T(-2), obs.Wt_D_sqrtdiagDinv_UUt, T(1), obs.Hττ) #second term
        mul!(obs.Wt_D_Ut_kr_Utt, obs.storage_ln, transpose(obs.Ut_kr_Ut))
        BLAS.syrk!('U', 'N', T(1), obs.Wt_D_Ut_kr_Utt, T(1), obs.Hττ) #third term
        copytri!(obs.Hττ, 'U')
        #wrt HτLγ
        # HτLγ = 2 W' * Diagonal(expwτ) * (L'Z'(V^-1) ⊙ Z'(V^-1))' * Cq
        # storage_ln = W' * Diagonal(expwτ) was computed above
        # storage_q◺n = Cq' * (L'Z'(V^-1) ⊙ Z'(V^-1))
        copyto!(obs.storage_qn, obs.Zt_Vinv)
        BLAS.trmm!('L', 'L', 'T', 'N', one(T), Lγ, obs.storage_qn)
        Ct_A_kr_B!(obs.storage_q◺n, obs.storage_qn, obs.Zt_Vinv)
        BLAS.gemm!('N', 'T', T(2), obs.storage_ln, obs.storage_q◺n, zero(T), obs.HτLγ)
        #wrt HLγLγ
        # HLγLγ = 2 [ C'(L'Z'(V^-1)ZL ⊗ Z'(V^-1)Z)C + C'(L'Z'(V^-1)Z ⊗ Z'(V^-1)ZL)KC ]
        # obs.storage_qq = (Z' (V^-1) Z) * L
        copyto!(obs.storage_qq, obs.Zt_Vinv_Z)
        BLAS.trmm!('R', 'L', 'N', 'N', T(1), Lγ, obs.storage_qq)
        # HLγLγ = C'(L'Z' (V^-1) Z ⊗ Z' (V^-1) ZL)KC first
        Ct_At_kron_A_KC!(fill!(obs.HLγLγ, 0), obs.storage_qq)
        # obs.storage_qq = L' * (Z' (V^-1) Z) * L
        BLAS.trmm!('L', 'L', 'T', 'N', T(1), Lγ, obs.storage_qq)
        Ct_A_kron_B_C!(obs.HLγLγ, obs.storage_qq, obs.Zt_Vinv_Z)
        lmul!(2, obs.HLγLγ)
    end
    obs.obj[1]
end

"""
    nlsv_obj!(m::WSVarLmmModel, needgrad::Bool, needhess:Bool, updateres::Bool)

Calculate the objective function of a `WSVarLmmModel` object and optionally the
gradient and hessian.
"""
function nlsv_obj!(
    m         :: WSVarLmmModel{T},
    needgrad  :: Bool = true,
    needhess  :: Bool = true,
    updateres :: Bool = false
    ) where T <: BlasReal
    # accumulate obj and gradient
    obj = zero(T)
    if needgrad
        fill!(   m.∇β, 0)
        fill!(   m.∇τ, 0)
        fill!(  m.∇Lγ, 0)
    end
    if needhess
        fill!(m.Hττ  , 0)
        fill!(m.HτLγ , 0)
        fill!(m.HLγLγ, 0)
    end
    if m.ismthrd[1]
        Threads.@threads for obs in m.data
            nlsv_obj!(obs, m.β, m.τ, m.Lγ, Val(m.iswtnls[1]),
                needgrad, needhess, updateres)
        end
        for obs in 1:length(m.data)
            wtobs = isempty(m.obswts) ? one(T) : T(m.obswts[obs])
            obj += wtobs * m.data[obs].obj[1]
            if needgrad
                BLAS.axpy!(wtobs, m.data[obs].∇τ   , m.∇τ )
                BLAS.axpy!(wtobs, m.data[obs].∇Lγ  , m.∇Lγ)
            end
            if needhess
                BLAS.axpy!(wtobs, m.data[obs].Hττ  , m.Hττ  )
                BLAS.axpy!(wtobs, m.data[obs].HτLγ , m.HτLγ )
                BLAS.axpy!(wtobs, m.data[obs].HLγLγ, m.HLγLγ)
            end
        end
    else
        for obs in 1:length(m.data)
            wtobs = isempty(m.obswts) ? one(T) : T(m.obswts[obs])
            obj += wtobs * nlsv_obj!(m.data[obs], m.β, m.τ, m.Lγ,
                Val(m.iswtnls[1]), needgrad, needhess, updateres)
            if needgrad
                BLAS.axpy!(wtobs, m.data[obs].∇τ   , m.∇τ )
                BLAS.axpy!(wtobs, m.data[obs].∇Lγ  , m.∇Lγ)
            end
            if needhess
                BLAS.axpy!(wtobs, m.data[obs].Hττ  , m.Hττ  )
                BLAS.axpy!(wtobs, m.data[obs].HτLγ , m.HτLγ )
                BLAS.axpy!(wtobs, m.data[obs].HLγLγ, m.HLγLγ)
            end
        end
    end
    # multiply m.∇Lγ by Lγ once instead of at observation level
    needgrad && BLAS.trmm!('R', 'L', 'N', 'N', T(2), m.Lγ, m.∇Lγ)
    obj
end

"""
    update_res!(obs, β)

Update the residual vector of `obs::WSVarLmmObs` according to `β`.
"""
function update_res!(
    obs :: WSVarLmmObs{T},
    β   :: Vector{T}
    ) where T <: BlasReal
    BLAS.gemv!('T', T(-1), obs.Xt, β, T(1), copyto!(obs.res, obs.y))
    obs.res2      .= abs2.(obs.res)
    obs.resnrm2[1] = sum(obs.res2)
    mul!(obs.ztres, obs.Zt, obs.res)
    obs.res
end

"""
    update_res!(m::WSVarLmmModel)

Update residual vector of each observation in `m` according to `m.β`.
"""
function update_res!(m::WSVarLmmModel{T}) where T <: BlasReal
    for obs in m.data
        update_res!(obs, m.β)
    end
    nothing
end

"""
    update_wtmat!(m::WSVarLmmModel)

Update the observation weight matrix according to the parameter values
`m.τ` and `m.Lγ`. Update `m.β` by WLS and update the residuals accordingly.
Precompute and store various objects needed to  evaluate the objective function,
gradient, and Hessian. At return,
- `m.data[i].∇β  = Xi' inv(Vi) ri`
- `m.data[i].Hββ = Xi' inv(Vi) Xi`
- `m.∇β  = sum_i Xi' inv(Vi) ri`
- `m.Hββ = sum_i Xi' inv(Vi) Xi`
"""
function update_wtmat!(m::WSVarLmmModel{T}) where T <: BlasReal
    p, q = m.p, m.q
    # update Dinv and U such that Vinv = Dinv - U * U'
    # accumulate quantities for updating β
    fill!(m.Hββ, 0)
    fill!(m.∇β , 0)
    for (idx, obs) in enumerate(m.data)
        # Form Dinv - UU^T (Woodbury structure of Vi inverse)
        n = length(obs.y)
        # Step 1: assemble Iq + Lt Zt diag(e^{-η}) Z L
        # storage_qn = Lt Zt Diagonal(e^{-0.5η})
        # storage_pn = Xt Diagonal(e^{-0.5η})
        # storage_n1 = Diagonal(e^{-0.5η}) * y
        copyto!(obs.storage_qn, obs.Zt)
        BLAS.trmm!('L', 'L', 'T', 'N', T(1), m.Lγ, obs.storage_qn)
        mul!(obs.expwτ, transpose(obs.Wt), m.τ)
        @inbounds for j in 1:n
            invsqrtj    = exp(-(1//2)obs.expwτ[j])
            obs.Dinv[j] = abs2(invsqrtj)
            for i in 1:q
                obs.storage_qn[i, j] *= invsqrtj
            end
            for i in 1:p
                obs.storage_pn[i, j] = invsqrtj * obs.Xt[i, j]
            end
            obs.storage_n1[j] = invsqrtj * obs.y[j]
        end
        # storage_qq = Iq + Lt Zt diag(e^{-η}) Z L
        BLAS.syrk!('U', 'N', T(1), obs.storage_qn, T(0), obs.storage_qq)
        @inbounds for i in 1:q
            obs.storage_qq[i, i] += 1
        end
        # Step 2: Cholesky: (Iq + Lt Zt diag(e^{-η}) Z L) = R'R
        LAPACK.potrf!('U', obs.storage_qq)
        # storage_qn = inv(R') Lt Zt diag(exp(-0.5η))
        BLAS.trsm!('L', 'U', 'T', 'N', T(1), obs.storage_qq, obs.storage_qn)
        # storage_qp = inv(R') Lt Zt diag(exp(-η)) X = U' X
        mul!(obs.storage_qp, obs.storage_qn, transpose(obs.storage_pn))
        # Step 3: accumulate X' Vinv X = X' Dinv X - X' U U' X
        BLAS.syrk!('U', 'N',  T(1), obs.storage_pn, T(0), obs.Hββ)
        BLAS.syrk!('U', 'T', T(-1), obs.storage_qp, T(1), obs.Hββ)
        copytri!(obs.Hββ, 'U')
        wtobs = isempty(m.obswts) ? one(T) : T(m.obswts[idx])
        BLAS.axpy!(wtobs, obs.Hββ, m.Hββ)
        # Step 4: accumulate X' Vinv y = X' Dinv y - X' U U' y
        mul!(obs.∇β, obs.storage_pn, obs.storage_n1)
        mul!(obs.storage_q1, obs.storage_qn, obs.storage_n1)
        BLAS.gemv!('T', T(-1), obs.storage_qp, obs.storage_q1, T(1), obs.∇β)
        BLAS.axpy!(wtobs, obs.∇β, m.∇β)
    end
    # update β by WLS
    copytri!(m.Hββ, 'U')
    copyto!(m.data[1].storage_pp, m.Hββ) # m.data[1].storage_pp as scratch space
    _, info = LAPACK.potrf!('U', m.data[1].storage_pp)
    info > 0 && throw("sum_i Xi' Vi^{-1} Xi is not positive definite")
    LAPACK.potrs!('U', m.data[1].storage_pp, copyto!(m.β, m.∇β))
    # update residuals according to new β
    update_res!(m)
    # precompute quantities for obj, grad and Hess evaluations
    for obs in m.data
        n = length(obs.y)
        # update obs.∇β = X' Vinv (y - Xβ) = X' Vinv y - X' Vinv X β
        BLAS.gemv!('N', T(-1), obs.Hββ, m.β, T(1), obs.∇β)
        # Update: Dinv_r, rt_Dinv_r, Ut = R'^{-1/2} Lt Zt diag(exp(-η))
        # storange_qn = R'^{-1/2} Lt Zt diag(exp(-0.5η)) from 1st loop
        copyto!(obs.Ut, obs.storage_qn)
        obs.rt_Dinv_r[1] = 0
        @inbounds for j in 1:n
            obs.Dinv_r[j]     = obs.Dinv[j] * obs.res[j]
            obs.rt_Dinv_r[1] += obs.Dinv[j] * obs.res2[j]
            invsqrtj          = sqrt(obs.Dinv[j])
            for i in 1:q
                obs.Ut[i, j] *= invsqrtj
            end
        end
        # Zt_Dinv_r
        mul!(obs.Zt_Dinv_r, obs.Zt, obs.Dinv_r)
        # rt_U and rt_UUt_r
        mul!(obs.rt_U, transpose(obs.res), transpose(obs.Ut))
        obs.rt_UUt_r[1] = abs2(norm(obs.rt_U))
        # rt_UUt
        mul!(obs.rt_UUt, obs.rt_U, obs.Ut)
        # rt_UUt_Z
        mul!(obs.rt_UUt_Z, obs.rt_UUt, transpose(obs.Zt))
        # Zt_Dinv_r
        mul!(obs.Zt_Dinv_r, obs.Zt, obs.Dinv_r)
        # storage_qq = Ut * Z
        mul!(obs.storage_qq, obs.Ut, transpose(obs.Zt))
        mul!(obs.Zt_UUt, transpose(obs.storage_qq), obs.Ut)
        # storage_qn = Zt * Dinv
        fill!(obs.diagUUt_Dinv, 0)
        @inbounds for j in 1:n
            invj = obs.Dinv[j]
            for i in 1:q
                obs.diagUUt_Dinv[j] += abs2(obs.Ut[i, j])
                obs.storage_qn[i, j] = obs.Zt[i, j] * invj
            end
            obs.diagUUt_Dinv[j] *= invj
        end
        copyto!(obs.Zt_Dinv, obs.storage_qn)
        # obs.Zt_Vinv = Z' * Dinv
        copyto!(obs.Zt_Vinv, obs.storage_qn)
        mul!(obs.Zt_Dinv_Z, obs.storage_qn, transpose(obs.Zt))
        # storage_qq = Zt * U
        mul!(obs.storage_qq, obs.Zt, transpose(obs.Ut))
        BLAS.syrk!('U', 'N', T(1), obs.storage_qq, T(0), obs.Zt_UUt_Z)
        copytri!(obs.Zt_UUt_Z, 'U')
        # for gradient wrt Lγ
        mul!(obs.Zt_UUt_rrt_Dinv_Z, transpose(obs.rt_UUt_Z), transpose(obs.Zt_Dinv_r))
        BLAS.syrk!('U', 'T', T(1), obs.rt_UUt_Z, T(0), obs.Zt_UUt_rrt_UUt_Z)
        copytri!(obs.Zt_UUt_rrt_UUt_Z, 'U')
        # for Hessian wrt τ
        map!(sqrt, obs.sqrtDinv_UUt, obs.diagUUt_Dinv)
        fill!(obs.Ut_kr_Ut, 0)
        kr_axpy!(obs.Ut, obs.Ut, obs.Ut_kr_Ut)
        # From earlier, obs.Zt_Vinv = Z' * Dinv, continue forming it
        BLAS.axpy!(T(-1), obs.Zt_UUt, obs.Zt_Vinv)
        copyto!(obs.Zt_Vinv_Z, obs.Zt_Dinv_Z)
        BLAS.axpy!(T(-1), obs.Zt_UUt_Z, obs.Zt_Vinv_Z)
        mul!(obs.Zt_Vinv_r, obs.Zt_Vinv, obs.res)
    end
    nothing
end


text:
type MultiGrid{T,A<:AbstractVector,B<:AbstractVector,C<:AbstractVector} <: AbstractGrid3d{T,1}
    c::A  # coords
    w::B  # weights
    p::C  # partition function
    atomptr::Vector{Int}  # atom i contained in atomptr[i]:atomptr[i+1]-1
    function MultiGrid(c::AbstractVector{Point3{T}},w::AbstractVector{T},
            p::AbstractVector{T},atomptr::Vector{Int})
        @assert size(c) == size(w)
        new(c,w,p,atomptr)
    end
end

MultiGrid{T}(c::AbstractVector{Point3{T}},w::AbstractVector{T},
        p::AbstractVector{T},atomptr::Vector{Int}) =
    MultiGrid{eltype(w),typeof(c),typeof(w),typeof(p)}(c,w,p,atomptr)

# empty grid
MultiGrid() = MultiGrid(Point3f[],Float64[],Float64[],Int[])

eachatom(g::MultiGrid) = [g.atomptr[i]:g.atomptr[i+1]-1 for i in 1:length(g.atomptr)-1]

function multigrid(atoms::Atoms)
    c = vcat([vec(at[:g]) for at in atoms]...)
    p = vcat([vec(at[:p](at[:g])) for at in atoms]...)
    w = p.*vcat([vec(at[:g].w) for at in atoms]...)
    atomptr = [1; cumsum([length(at[:g]) for at in atoms])+1...]
    MultiGrid(c,w,p,atomptr)
end

max_stars_repo_name: johncsnyder/Element


text:
mutable struct NDFTPlan{T,D} <: AbstractNFFTPlan{T,D,1}
  N::NTuple{D,Int64}
  M::Int64
  x::Matrix{T}
end

AbstractNFFTs.size_in(p::NDFTPlan) = p.N
AbstractNFFTs.size_out(p::NDFTPlan) = (p.M,)

mutable struct NDCTPlan{T,D} <: AbstractNFCTPlan{T,D,1}
  N::NTuple{D,Int64}
  M::Int64
  x::Matrix{T}
end

AbstractNFFTs.size_in(p::NDCTPlan) = p.N
AbstractNFFTs.size_out(p::NDCTPlan) = (p.M,)

mutable struct NDSTPlan{T,D} <: AbstractNFSTPlan{T,D,1}
  N::NTuple{D,Int64}
  M::Int64
  x::Matrix{T}
end

AbstractNFFTs.size_in(p::NDSTPlan) = p.N .- 1
AbstractNFFTs.size_out(p::NDSTPlan) = (p.M,)

mutable struct NNDFTPlan{T} <: AbstractNNFFTPlan{T,1,1}
  N::Int64
  M::Int64
  x::Matrix{T}
  y::Matrix{T}
end

AbstractNFFTs.size_in(p::NNDFTPlan) = (p.N,)
AbstractNFFTs.size_out(p::NNDFTPlan) = (p.M,)

### constructors ###

function NDFTPlan(x::Matrix{T}, N::NTuple{D,Int}; kwargs...) where {T,D}

  if D != size(x,1)
    throw(ArgumentError("Nodes x have dimension $(size(x,1)) != $D"))
  end

  M = size(x, 2)

  return NDFTPlan{T,D}(N, M, x)
end

function NDCTPlan(x::Matrix{T}, N::NTuple{D,Int}; kwargs...) where {T,D}

  if D != size(x,1)
    throw(ArgumentError("Nodes x have dimension $(size(x,1)) != $D"))
  end

  M = size(x, 2)

  return NDCTPlan{T,D}(N, M, x)
end

function NDSTPlan(x::Matrix{T}, N::NTuple{D,Int}; kwargs...) where {T,D}

  if D != size(x,1)
    throw(ArgumentError("Nodes x have dimension $(size(x,1)) != $D"))
  end

  M = size(x, 2)

  return NDSTPlan{T,D}(N, M, x)
end

function NNDFTPlan(x::Matrix{T}, y::Matrix{T}; kwargs...) where {T}

  M = size(x, 2)
  N = size(y, 2)

  return NNDFTPlan{T}(N, M, x, y)
end

### ndft functions ###

ndft(x::AbstractArray, f::AbstractArray; kargs...) =
    NDFTPlan(x, size(f); kargs...) * f

ndft_adjoint(x, N, fHat::AbstractVector; kargs...) =
    adjoint(NDFTPlan(x, N; kargs...)) * fHat

ndct(x::AbstractArray, f::AbstractArray; kargs...) =
    NDCTPlan(x, size(f); kargs...) * f

ndct_transposed(x, N, fHat::AbstractVector; kargs...) =
    transpose(NDCTPlan(x, N; kargs...)) * fHat

function LinearAlgebra.mul!(g::AbstractArray{Tg}, plan::NDFTPlan{Tp,D}, f::AbstractArray{T,D}) where {D,Tp,T,Tg}

    plan.N == size(f) ||
        throw(DimensionMismatch("Data f is not consistent with NDFTPlan"))
    plan.M == length(g) ||
        throw(DimensionMismatch("Output g is inconsistent with NDFTPlan"))

    g .= zero(Tg)

    for l=1:prod(plan.N)
        idx = CartesianIndices(plan.N)[l]

        for k=1:plan.M
            arg = zero(T)
            for d=1:D
                arg += plan.x[d,k] * ( idx[d] - 1 - plan.N[d] / 2 )
            end
            g[k] += f[l] * cis(-2*pi*arg)
        end
    end

    return g
end

function LinearAlgebra.mul!(g::AbstractArray{Tg,D}, pl::Adjoint{Complex{Tp},<:NDFTPlan{Tp,D}}, fHat::AbstractVector{T}) where {D,Tp,T,Tg}
  p = pl.parent

  p.M == length(fHat) ||
      throw(DimensionMismatch("Data f inconsistent with NDFTPlan"))
  p.N == size(g) ||
      throw(DimensionMismatch("Output g inconsistent with NDFTPlan"))

  g .= zero(Tg)

  for l=1:prod(p.N)
      idx = CartesianIndices(p.N)[l]

      for k=1:p.M
          arg = zero(T)
          for d=1:D
              arg += p.x[d,k] * ( idx[d] - 1 - p.N[d] / 2 )
          end
          g[l] += fHat[k] * cis(2*pi*arg)
      end
  end

  return g
end

function LinearAlgebra.mul!(g::AbstractArray{Tg}, p::NDCTPlan{Tp,D}, f::AbstractArray{T,D}) where {D,Tp,T,Tg}

    p.N == size(f) ||
        throw(DimensionMismatch("Data f is not consistent with NDCTPlan"))
    p.M == length(g) ||
        throw(DimensionMismatch("Output g is inconsistent with NDCTPlan"))

    g .= zero(Tg)

    for l=1:prod(p.N)
        idx = CartesianIndices(p.N)[l]

        for k=1:p.M
            arg = one(T)
            for d=1:D
                arg *= cos( 2 * pi * p.x[d,k] * ( idx[d] - 1 ) )
            end
            g[k] += f[l] * arg
        end
    end

    return g
end

function LinearAlgebra.mul!(g::AbstractArray{Tg,D}, pl::Transpose{Tp,<:NDCTPlan{Tp,D}}, fHat::AbstractVector{T}) where {D,Tp,T,Tg}
  p = pl.parent

  p.M == length(fHat) ||
      throw(DimensionMismatch("Data f inconsistent with NDCTPlan"))
  p.N == size(g) ||
      throw(DimensionMismatch("Output g inconsistent with NDCTPlan"))

  g .= zero(Tg)

  for l=1:prod(p.N)
      idx = CartesianIndices(p.N)[l]

      for k=1:p.M
          arg = one(T)
          for d=1:D
              arg *= cos( 2 * pi * p.x[d,k] * ( idx[d] - 1 ) )
          end
          g[l] += fHat[k] * arg
      end
  end

  return g
end

function LinearAlgebra.mul!(g::AbstractArray{Tg}, p::NDSTPlan{Tp,D}, f::AbstractArray{T,D}) where {D,Tp,T,Tg}

    p.N == size(f) .+ 1 ||
        throw(DimensionMismatch("Data f is not consistent with NDSTPlan"))
    p.M == length(g) ||
        throw(DimensionMismatch("Output g is inconsistent with NDSTPlan"))

    g .= zero(Tg)

    for l=1:prod(p.N .- 1)
        idx = CartesianIndices(p.N .- 1)[l]

        for k=1:p.M
            arg = one(T)
            for d=1:D
                arg *= sin( 2 * pi * p.x[d,k] * ( idx[d] ) )
            end
            g[k] += f[l] * arg
        end
    end

    return g
end

function LinearAlgebra.mul!(g::AbstractArray{Tg,D}, pl::Transpose{Tp,<:NDSTPlan{Tp,D}}, fHat::AbstractVector{T}) where {D,Tp,T,Tg}
  p = pl.parent

  p.M == length(fHat) ||
      throw(DimensionMismatch("Data f inconsistent with NDSTPlan"))
  p.N == size(g) .+ 1 ||
      throw(DimensionMismatch("Output g inconsistent with NDSTPlan"))

  g .= zero(Tg)

  for l=1:prod(p.N .- 1)
      idx = CartesianIndices(p.N .- 1)[l]

      for k=1:p.M
          arg = one(T)
          for d=1:D
              arg *= sin( 2 * pi * p.x[d,k] * ( idx[d] ) )
          end
          g[l] += fHat[k] * arg
      end
  end

  return g
end

function LinearAlgebra.mul!(g::AbstractArray{Tg}, p::NNDFTPlan{Tp}, f::AbstractArray{T}) where {Tp,T,Tg}

  g .= zero(Tg)

  for l=1:p.N
      for k=1:p.M
          arg = zero(T)
          for d=1:size(p.x,1)
              arg += p.x[d,k] * p.y[d,l]
          end
          g[k] += f[l] * cis(-2*pi*arg)
      end
  end

  return g
end

function LinearAlgebra.mul!(g::AbstractArray{Tg}, pl::Adjoint{Complex{Tp},<:NNDFTPlan{Tp}}, fHat::AbstractVector{T}) where {Tp,T,Tg}
  p = pl.parent
  g .= zero(Tg)

  for l=1:p.N
      for k=1:p.M
          arg = zero(T)
          for d=1:size(p.x,1)
            arg += p.x[d,k] * p.y[d,l]
          end
          g[l] += fHat[k] * cis(2*pi*arg)
      end
  end

  return g
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule marble_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("marble")
JLLWrappers.@generate_main_file("marble", UUID("678d7417-9a84-558b-a975-2deb8d71bebc"))
end  # module marble_jll


text:
using Test, ModelingToolkit

@parameters t

@connector function Foo(;name)
    @variables x(t)
    ODESystem(Equation[], t, [x], [], defaults=Dict(x=>1.0))
end

@connector function Goo(;name)
    @variables x(t)
    @parameters p
    ODESystem(Equation[], t, [x], [p], defaults=Dict(x=>1.0, p=>1.0))
end

function ModelingToolkit.connect(::Type{<:Foo}, ss...)
    n = length(ss)-1
    eqs = Vector{Equation}(undef, n)
    for i in 1:n
        eqs[i] = ss[i].x ~ ss[i+1].x
    end
    eqs
end

@named f1 = Foo()
@named f2 = Foo()
@named f3 = Foo()
@named f4 = Foo()
@named g = Goo()

@test isequal(connect(f1, f2), [f1.x ~ f2.x])
@test_throws ArgumentError connect(f1, g)

# Note that since there're overloadings, these tests are not re-runable.
ModelingToolkit.promote_connect_rule(::Type{<:Foo}, ::Type{<:Goo}) = Foo
@test isequal(connect(f1, g), [f1.x ~ g.x])
@test isequal(connect(f1, f2, g), [f1.x ~ f2.x; f2.x ~ g.x])
@test isequal(connect(f1, f2, g, f3), [f1.x ~ f2.x; f2.x ~ g.x; g.x ~ f3.x])
@test isequal(connect(f1, f2, g, f3, f4), [f1.x ~ f2.x; f2.x ~ g.x; g.x ~ f3.x; f3.x ~ f4.x])
ModelingToolkit.promote_connect_rule(::Type{<:Goo}, ::Type{<:Foo}) = Foo
@test isequal(connect(f1, g), [f1.x ~ g.x])
# test conflict
ModelingToolkit.promote_connect_rule(::Type{<:Goo}, ::Type{<:Foo}) = Goo
@test_throws ArgumentError connect(f1, g)


text:
using Pkg
Pkg.activate(".")

using Documenter
using Literate
using ProBeTune

# # generate examples
# examples = [
#     joinpath(@__DIR__, "..", "examples", "simple_graph_example.jl")
#     #joinpath(@__DIR__, "..", "examples", "swing_equation.jl")
# ]
# OUTPUT = joinpath(@__DIR__, "src/generated")
# isdir(OUTPUT) && rm(OUTPUT, recursive=true)
# mkpath(OUTPUT)

# for ex in examples
#     Literate.markdown(ex, OUTPUT)
#     Literate.script(ex, OUTPUT)
# end

# makedocs(;
#     modules=[BAC],
#     authors="Frank Hellmann, Ekaterina Zolotarevskaia",
#     repo="https://github.com/FHell/BAC.jl",
#     sitename="BAC.jl",
#     format=Documenter.HTML(;
#         prettyurls=get(ENV, "CI", "false") == "true",
#         canonical="https://fhell.github.io/BAC.jl",
#         assets=String[],
#     ),
#     pages=[
#         "Home" => "index.md",
#         "Examples" => ["Simple graph example" => "generated/simple_graph_example.md"]
#     ],
# )

# deploydocs(;
#     repo="github.com/FHell/BAC.jl.git",
#     devbranch="master",
#     # push_preview=true,
# )


text:
module Intensity

using Distributions, Random

export InhomogeneousExponential

EPSILON_DIGIT = 10

AVR = AbstractVector{<:Real}
AV = AbstractVector

abstract type AbstractTermStructure end

abstract type AbstractInterExtrapolation end
abstract type PWL <: AbstractInterExtrapolation end
abstract type Constant <: AbstractInterExtrapolation end

"""
    TermStructure <: AbstractTermStructure

basic struct to store scalar quantities possessing a nontrivial term structure.
    Terms are stored as floating point numbers.
    An interpolation scheme has to be specified

    Fields:
    - terms::AbstractVector{<:Real}
    - values::AbstractVector{<:Real}
    - interpolation::Type{<:AbstractInterExtrapolation}

    _N.B._ : terms have to be chronologically ordered
"""
struct TermStructure <: AbstractTermStructure
    terms::AVR
    values::AVR
    interpolation::Type{<:AbstractInterExtrapolation}

    function TermStructure(terms::AVR, values::AVR, interpolation::Type{<:AbstractInterExtrapolation})
        @assert sort(terms) == terms "Terms should be time-ordered"
        return new(terms, values, interpolation)
    end
end

function bracket_term(terms::AVR, target::Float64)
    lowidxs = findall(terms .<= target)
    return isempty(lowidxs) ? 0 : maximum(lowidxs)
end

"""
interextrapolate
    arguments
        - ::Type{<:AbstractInterExtrapolation},
        - terms::AV,
        - qtys::AV,
        - term::Float64,

    given the terms and the quantities corresponding to each of them, determine
    the interpolated/extrapolated value according to the desired scheme.

    _Notes_
    - Unless otherwise specified, the initial time is 0.0.
    - For PWL, if the first time peg (0.0) is omitted, the initial value is assumed to
    be zero
    - For Constant, the interpolation is continuous on the left (i.e. values
    are assumed to be on (a_i, a_{i+1}])
"""
function interextrapolate(
    ::Type{<:AbstractInterExtrapolation},
    terms::AV,
    qtys::AV,
    term::Float64)
    return error("interextrapolate not implemented.")
end

function interextrapolate(::Type{<:PWL}, terms::AV, qtys::AV, term::Float64)
    lowidx = bracket_term(terms, term)
    if lowidx == 0
        t_low = 0.0
        t_high = terms[1]
        qty_low = 0
        qty_high = qtys[1]
    else
        t_low = terms[lowidx]
        qty_low = qtys[lowidx]
        if lowidx < length(terms)
            t_high = terms[lowidx + 1]
            qty_high = qtys[lowidx + 1]
        else
            t_high = NaN
            qty_high = qtys[end]
        end
    end
    delta = qty_high - qty_low
    if delta == 0.0
        return qty_low
    else
        slope = delta / (t_high - t_low)
        dt = term - t_low
        return delta / (t_high - t_low) * (term - t_low) + qty_low
    end
end

function interextrapolate(::Type{Constant}, terms::AV, qtys::AV, term::Float64)
    highidxs = findall(terms .>= term)
    return isempty(highidxs) ? qtys[end] : qtys[minimum(highidxs)]
end

function interextrapolate(ts::AbstractTermStructure, term::Float64)
    @assert term >= 0
    return interextrapolate(ts.interpolation, ts.terms, ts.values, term)
end

function interextrapolate(ts::AbstractTermStructure, terms::AVR)
    return map(x -> interextrapolate(ts.interpolation, ts.terms, ts.values, x), terms)
end

"""
    SurvivalProbability

    A simple TermStructure enforcing the semantics of survival probabilities,
    without interpolation
"""
struct SurvivalProbability <: AbstractTermStructure
    terms::AVR
    probabilities::AVR

    function SurvivalProbability(terms::AVR, probs::AVR)
        sort(terms) != terms && throw(DomainError("Terms must be sorted"))
        length(terms) != length(probs) && throw(DomainError("Terms and probabilities have different length"))
        any(probs[2:end] .> probs[1:end-1]) && throw(DomainError("Survival probabilities must be decreasing"))
        any(probs .< 0.0) && throw(DomainError("Survival probabilities must be positive"))
        any(probs .> 1.0) && throw(DomainError("Survival probabilities must be not greater than one"))
        return new(terms, probs)
    end
end

"""
    DeterministicIntensity

    implement deterministic intensity semantics, with interpolation/extrapolation
    method specified
"""
mutable struct DeterministicIntensity <: AbstractTermStructure
    terms::AVR
    values::AVR
    interpolation::Type{<:AbstractInterExtrapolation}

    function DeterministicIntensity(
        terms::AVR,
        intensities::AVR,
        interpolation::Type{<:AbstractInterExtrapolation},
    )
        any(intensities .< 0) && throw(DomainError("All intensities have to be positive"))
        length(terms) != length(intensities) && throw(DomainError("Terms and intensities have different lengths"))
        return new(terms, intensities, interpolation)
    end
end

DeterministicIntensity(terms, intensities) = DeterministicIntensity(terms, intensities, PWL)

"""
    get_rolling_probabilities

    for each term of the DeterministicIntensity, compute the survival probability
"""
function get_rolling_probabilities(di::DeterministicIntensity)
    p = 1.0
    terms = vcat([0.0], di.terms)
    deltas = terms[2:end] - terms[1:end-1]
    intensities = interextrapolate(di, terms)
    if di.interpolation === Constant
        areas = deltas .* intensities[2:end]
    elseif di.interpolation === PWL
        areas = deltas .* (intensities[2:end] + intensities[1:end-1]) * 0.5
    end
    probs = zeros(size(di.terms))
    for (i, area) in enumerate(areas)
        p *= exp(-area)
        probs[i] = p
    end
    return probs
end

"""
    get_survival_probability

    args
    - di::DeterministicIntensity
    - term::Float64
"""
function get_survival_probability(di::DeterministicIntensity, term::Float64)
    lowidx = bracket_term(di.terms, term)
    probs = get_rolling_probabilities(di)
    if lowidx == 0
        base = 0.0
        base_int = 0.0
        p = 1.0
    else
        base = di.terms[lowidx]
        base_int = di.values[lowidx]
        p = probs[lowidx]
    end
    high_int = interextrapolate(di, term)
    delta = term - base
    height = di.interpolation == Constant ? high_int : (base_int + high_int) * 0.5
    area = delta * height
    return p * exp(-area)
end

function get_survival_probability(di::DeterministicIntensity, terms::AbstractArray{Float64})
    return map(x -> get_survival_probability(di, x), terms)
end

"""
    calibrate_intensity(sp::SurvivalProbability, interpolation::Type{<:AbstractInterExtrapolation})

    args
    - sp::SurvivalProbability,
    - interpolation
"""
function calibrate_intensity(sp::SurvivalProbability, interpolation::Type{<:AbstractInterExtrapolation})
    terms = sp.terms
    probs = sp.probabilities
    intensities = Float64[]
    new_terms = Float64[]
    if terms[1] == 0
        @assert probs[1] == 1 "survival probability at time zero must be one"
        terms = terms[2:end]
        probs = probs[2:end]
    end
    logp_prev = 0.0
    base = 0.0
    base_int = 0.0
    for i in 1:length(terms)
        term = terms[i]
        delta = term - base
        area = logp_prev - log(probs[i])
        height = area / delta
        if interpolation == PWL
            new_int = 2 * height - base_int
        else
            new_int = height
        end
        new_int = round(new_int; digits=EPSILON_DIGIT)
        @assert new_int >= 0.0 "Intensity has to be positive: $new_int, $height, $base_int, $base"
        base_int = new_int
        push!(new_terms, term)
        push!(intensities, new_int)
        base = term
        logp_prev = log(probs[i])
    end
    return DeterministicIntensity(new_terms, intensities, interpolation)
end

function make_grid(start::Float64, stop::Float64, npts::Integer; dates=nothing)
    grid = range(start, stop; length=npts)
    if dates != nothing
        grid = sort(union(grid, dates))
    end
    return grid
end

struct InhomogeneousExponential <: Sampleable{Univariate, Continuous}
    intensity::DeterministicIntensity
    grid::AVR
end

"""
    rand(rng::AbstractRNG, ie::InhomogeneousExponential)

Sample using Lewis'[1] thinning algorithm

References:
[1] Lewis, P. A. W. and G. S. Shedler 1979.
    `Simulation of nonhomogeneous Poisson processes by thinning.' Naval Res. Logistics Quart, 26:403– 413.
"""
function Base.rand(rng::AbstractRNG, ie::InhomogeneousExponential)
    t = []
    s = [0.0]
    λs = interextrapolate(ie.intensity, ie.grid)
    λmax = maximum(λs)
    while length(t) == 0
        ξ = - log(rand(rng)) / λmax
        s_new = s[end] + ξ
        push!(s, s_new)
        λnew = interextrapolate(ie.intensity, s_new)
        θ = λnew / λmax
        ζ = rand(rng)
        ζ <= θ && push!(t, s_new)
    end
    return first(t)
end

end

max_stars_repo_name: lsindoni/Intensity

max_issues_repo_name: lsindoni/Intensity

max_forks_repo_name: lsindoni/Intensity


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct RenewableDispatch <: RenewableGen
        name::String
        available::Bool
        bus::Bus
        activepower::Float64
        reactivepower::Float64
        tech::TechRenewable
        op_cost::TwoPartCost
        _forecasts::InfrastructureSystems.Forecasts
        internal::InfrastructureSystemsInternal
    end

# Arguments
- `name::String`
- `available::Bool`
- `bus::Bus`
- `activepower::Float64`
- `reactivepower::Float64`
- `tech::TechRenewable`
- `op_cost::TwoPartCost`
- `_forecasts::InfrastructureSystems.Forecasts`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct RenewableDispatch <: RenewableGen
    name::String
    available::Bool
    bus::Bus
    activepower::Float64
    reactivepower::Float64
    tech::TechRenewable
    op_cost::TwoPartCost
    _forecasts::InfrastructureSystems.Forecasts
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function RenewableDispatch(name, available, bus, activepower, reactivepower, tech, op_cost, _forecasts=InfrastructureSystems.Forecasts(), )
    RenewableDispatch(name, available, bus, activepower, reactivepower, tech, op_cost, _forecasts, InfrastructureSystemsInternal())
end

function RenewableDispatch(; name, available, bus, activepower, reactivepower, tech, op_cost, _forecasts=InfrastructureSystems.Forecasts(), )
    RenewableDispatch(name, available, bus, activepower, reactivepower, tech, op_cost, _forecasts, )
end

# Constructor for demo purposes; non-functional.

function RenewableDispatch(::Nothing)
    RenewableDispatch(;
        name="init",
        available=false,
        bus=Bus(nothing),
        activepower=0.0,
        reactivepower=0.0,
        tech=TechRenewable(nothing),
        op_cost=TwoPartCost(nothing),
        _forecasts=InfrastructureSystems.Forecasts(),
    )
end

"""Get RenewableDispatch name."""
get_name(value::RenewableDispatch) = value.name
"""Get RenewableDispatch available."""
get_available(value::RenewableDispatch) = value.available
"""Get RenewableDispatch bus."""
get_bus(value::RenewableDispatch) = value.bus
"""Get RenewableDispatch activepower."""
get_activepower(value::RenewableDispatch) = value.activepower
"""Get RenewableDispatch reactivepower."""
get_reactivepower(value::RenewableDispatch) = value.reactivepower
"""Get RenewableDispatch tech."""
get_tech(value::RenewableDispatch) = value.tech
"""Get RenewableDispatch op_cost."""
get_op_cost(value::RenewableDispatch) = value.op_cost
"""Get RenewableDispatch _forecasts."""
get__forecasts(value::RenewableDispatch) = value._forecasts
"""Get RenewableDispatch internal."""
get_internal(value::RenewableDispatch) = value.internal