text:
type GroupLassoData{T<:AbstractFloat, I}
  groups::Vector{Int64}
  numActive::Int64
  groupToIndex::Vector{I}
  xold::Vector{T}
  Asvd::Vector{Any}
  maxGroupLen::Int64
  btild::Vector{T}      # vector of max group len
  x1::Vector{T}         # vector of max group len
  Utb::Vector{T}         # vector of max group len
end

function GroupLassoData{T<:AbstractFloat, I}(
    x::StridedVector{T},
    groupToIndex::Vector{I};
    options::LassoOptions = LassoOptions()
    )
  zero_thr = options.zerothr
  numElem = length(groupToIndex)
  groups = [1:numElem;]
  numActive = 0
  maxGroupLen = 0
  for gi = 1:numElem
    group_len = length(groupToIndex[gi])
    if maxGroupLen < group_len
      maxGroupLen = group_len
    end
    if vecnorm(sub(x, groupToIndex[gi])) >= zero_thr
      numActive += 1
      groups[numActive], groups[gi] = groups[gi], groups[numActive]
    end
  end
  GroupLassoData{T, I}(groups, numActive, groupToIndex,
                       zeros(length(x)), Array(Any, numElem),
                       maxGroupLen,
                       zeros(maxGroupLen),zeros(maxGroupLen),zeros(maxGroupLen)
                       )
end

function _add_violator!{T<:AbstractFloat}(
    gld::GroupLassoData,
    x::StridedVector{T},
    A::StridedMatrix{T},
    b::StridedVector{T},
    λ::StridedVector{T};
    options::LassoOptions   =   LassoOptions()
    )
  changed = false

  zerothr = options.zerothr
  groups = gld.groups
  numActive = gld.numActive
  groupToIndex = gld.groupToIndex
  numElem = length(groups)
  # check for things to be removed from the active set
  i = 0
  @inbounds while i < numActive
    i = i + 1
    t = groups[i]
    if ~isdefined(gld.Asvd, t)
      gld.Asvd[t] = svdfact(sub(A, groupToIndex[t], groupToIndex[t]))
    end
    if _group_norm(x, gld, t) < zerothr
      _fill_zero!(x, gld, t)
      changed = true
      groups[numActive], groups[i] = groups[i], groups[numActive]
      numActive -= 1
      i = i - 1
    end
  end
  gld.numActive = numActive

  btild = gld.btild
  I = 0
  V = zero(T)
  @inbounds for i=numActive+1:numElem
    t = groups[i]
    group_norm = _Axkpb!(btild, A, b, x, t, gld)
    nV = group_norm - λ[t]
    if V < nV
      I = i
      V = nV
    end
  end
  if I > 0 && V > options.gradtol
    t = groups[I]
    if ~isdefined(gld.Asvd, t)
      gld.Asvd[t] = svdfact(sub(A, groupToIndex[t], groupToIndex[t]))
    end
    changed = true
    numActive += 1
    groups[numActive], groups[I] = groups[I], groups[numActive]
  end
  gld.numActive = numActive
  changed
end

function _group_norm{T<:AbstractFloat}(
    x::StridedVector{T},
    gld::GroupLassoData{T},
    j::Int64
    )
  r = zero(T)
  @inbounds for i in gld.groupToIndex[j]
    r += x[i]^2.
  end
  sqrt(r)
end

function _fill_zero!{T<:AbstractFloat, I}(
    x::StridedVector{T},
    gld::GroupLassoData{T, I},
    j::Int64
    )
  @inbounds for i in gld.groupToIndex[j]
    x[i] = zero(T)
  end
  x
end

# computes x[-gj]'A[-gj, gj] + b[gj]'
# returns norm of res
function _Axkpb!{T<:AbstractFloat}(
    res::StridedVector{T},
    A::StridedMatrix{T},
    b::StridedVector{T},
    x::StridedVector{T},
    j::Int64,
    gld::GroupLassoData
    )
  normRes = zero(T)

  groups = gld.groups
  numActive = gld.numActive
  groupToIndex = gld.groupToIndex

  group_j = groupToIndex[j]
  indC = 0
  @inbounds for column_ind in group_j
    indC += 1
    res[indC] = zero(T)
    for i=1:numActive
      gr = groups[i]
      if gr == j
        continue
      end
      for ind in groupToIndex[gr]
        res[indC] += x[ind] * A[ind, column_ind]
      end
    end
    res[indC] += b[column_ind]
    normRes += res[indC]^2.
  end
  sqrt(normRes)
end

function _tgemv!{T<:AbstractFloat}(
    β::T,
    A::StridedMatrix{T},
    b::StridedVector{T},
    α::T,
    out::StridedVector{T})

  nr, nc = size(A)
  @inbounds for c=1:nc
    v = zero(T)
    for r=1:nr
      v += b[r] * A[r, c]
    end
    out[c] = out[c] * α + β*v
  end
  out
end

function _min_one_group!{T<:AbstractFloat}(
    x::StridedVector{T},
    A::SVD{T,T},
    b::StridedVector{T},
    gld::GroupLassoData,
    j::Int64,
    λ::T;
    options::LassoOptions   =   LassoOptions()
    )

  gjToIndex = gld.groupToIndex[j]
  x1 = gld.x1
  Utb = gld.Utb
  U = A[:U]
  D = A[:S]
  Vt = A[:Vt]

  group_len = size(U, 1)
  gn = zero(T)
  @inbounds for i=1:group_len
    gn += b[i]^2.
  end
  gn = sqrt(gn) - λ
  lb = gn / D[1]
  ub = gn / D[end]
  #
  _tgemv!(-1., U, b, 0., Utb)
  @inbounds while true
    ##
    for i=1:group_len
      x1[i] = Utb[i] / (D[i] + λ / gn)
    end
    #     At_mul_B!(x, Vt, x1)
    en = zero(T)
    for i=1:group_len
      xind = gjToIndex[i]
      x[xind] = zero(T)
      for j=1:group_len
        x[xind] += Vt[j, i] * x1[j]
      end
      en += x[xind]^2.
    end
    en = sqrt(en)
    ##
    if abs(gn - en) < options.xtol || ub - lb < options.xtol
      break
    else
      if en > gn
        lb = gn
      else
        ub = gn
      end
      gn = (ub + lb) / 2.
    end
  end
  x
end

function _group_lasso!{T<:AbstractFloat}(
    x::StridedVector{T},
    A::StridedMatrix{T},
    b::StridedVector{T},
    λ::StridedVector{T},
    gld::GroupLassoData;
    options::LassoOptions   =   LassoOptions()
    )

  numActive = gld.numActive
  groups = gld.groups
  groupToIndex = gld.groupToIndex
  btild = gld.btild
  xold = gld.xold

  iter = 0
  @inbounds while true
    iter += 1
    maxUpdate = zero(T)
    for i=1:numActive
      j = groups[i]
      group_norm = _Axkpb!(btild, A, b, x, j, gld)
      if group_norm <= λ[j]
        _fill_zero!(x, gld, j)
      else
        _min_one_group!(x, gld.Asvd[j], btild, gld, j, λ[j]; options=options)
      end
      # compute norm of change
      norm_diff = zero(T)
      for k in groupToIndex[j]
        norm_diff += (x[k] - xold[k])^2.
      end
      norm_diff = sqrt(norm_diff)
      if norm_diff > maxUpdate
        maxUpdate = norm_diff
      end
    end
    if iter > options.max_inner_iter || maxUpdate < options.xtol
      break
    end
    copy!(xold, x)
  end
  x
end

group_lasso!{T<:AbstractFloat, I}(
    x::StridedVector{T},
    A::StridedMatrix{T},
    b::StridedVector{T},
    λ::T,
    groupToIndex::Vector{I};
    options::LassoOptions   =   LassoOptions(),
    gld                     =   nothing
    ) = group_lasso!(x, A, b, fill(λ, length(groupToIndex)), groupToIndex; options=options, gld=gld)

function group_lasso!{T<:AbstractFloat, I}(
    x::StridedVector{T},
    A::StridedMatrix{T},
    b::StridedVector{T},
    λ::StridedVector{T},
    groupToIndex::Vector{I};
    options::LassoOptions   =   LassoOptions(),
    gld                     =   nothing
    )

  maxoutiter = options.max_outer_iter
  if is(gld, nothing)
    gld = GroupLassoData(x, groupToIndex; options=options)
  end
  _add_violator!(gld, x, A, b, λ; options=options)

  for outiter=1:maxoutiter
    _group_lasso!(x, A, b, λ, gld; options=options)
    if ~_add_violator!(gld, x, A, b, λ; options=options)
      break
    end
  end
  x, gld
end

max_stars_repo_name: mlakolar/CDLasso

max_issues_repo_name: mlakolar/CDLasso

max_forks_repo_name: mlakolar/CDLasso


text:
function GPSBlockIIA(arg0::jdouble, arg1::jdouble, arg2::AbsoluteDate, arg3::AbsoluteDate, arg4::ExtendedPVCoordinatesProvider, arg5::Frame)
    return GPSBlockIIA((jdouble, jdouble, AbsoluteDate, AbsoluteDate, ExtendedPVCoordinatesProvider, Frame), arg0, arg1, arg2, arg3, arg4, arg5)
end

function get_default_yaw_rate(::Type{GPSBlockIIA}, arg0::jint)
    return jcall(GPSBlockIIA, "getDefaultYawRate", jdouble, (jint,), arg0)
end

ext: jl


text:
function part1(input)
    and_mask = 2^36 - 1
    or_mask = 0
    memory = Dict{Int, Int}()
    for line in readlines(input)
        if startswith(line, "mask")
            mask = last(split(line, " = "))
            and_mask = parse(Int, replace(mask, "X" => "1"), base = 2)
            or_mask = parse(Int, replace(mask, "X" => "0"), base = 2)
        else
            address_string, value_string = match(r"mem\[(\d+)\] = (\d+)",
                                                 line).captures
            address = parse(Int, address_string)
            value = parse(Int, value_string)
            memory[address] = (value & and_mask) | or_mask
        end
    end
    return sum(values(memory))
end

function part2(input)
    and_mask = 2^36 - 1
    or_mask = 0
    float_bits = Int[]
    memory = Dict{Int, Int}()
    for line in readlines(input)
        if startswith(line, "mask")
            mask = last(split(line, " = "))
            or_mask = parse(Int, replace(mask, "X" => "0"), base = 2)
            and_mask = parse(Int, replace(replace(mask, "0" => "1"),
                                          "X" => "0"), base = 2)
            float_bits = 36 .- findall(split(mask, "") .== "X")
        else
            address_string, value_string = match(r"mem\[(\d+)\] = (\d+)",
                                                 line).captures
            address = parse(Int, address_string)
            value = parse(Int, value_string)
            address |= or_mask
            for i = 0:((1 << length(float_bits)) - 1)
                address &= and_mask
                for j = 1:length(float_bits)
                    if (i & (1 << (j - 1))) != 0
                        address |= 1 << (float_bits[j])
                    end
                end
                memory[address] = value
            end
        end
    end
    return sum(values(memory))
end


text:
using DiagonalSBP
using PGFPlotsX: @pgf, Axis, Plot, pgfsave, Table, LogLogAxis
using ForwardDiff: derivative
using OrdinaryDiffEq: ODEProblem, solve, Tsit5
using SparseArrays: sparsevec
using Printf: @sprintf

# This uses the fourth-order, low-storage, Runge--Kutta scheme of
# Carpenter and Kennedy (1994) (in their notation (5,4) 2N-Storage RK
#                               scheme).
#
# @TECHREPORT{CarpenterKennedy1994,
#             author = {M.~H. Carpenter and C.~A. Kennedy},
#             title = {Fourth-order {2N-storage} {Runge-Kutta}
#                      schemes},
#             institution = {National Aeronautics and Space
#                            Administration},
#             year = {1994},
#             number = {NASA TM-109112},
#             address = {Langley Research
#                        Center, Hampton, VA},
#            }
function timestep!(q, f!, p, dt, (t0, t1), Δq = similar(q), Δq2 = similar(q))
    T = eltype(q)

    RKA = (
        T(0),
        T(-567301805773 // 1357537059087),
        T(-2404267990393 // 2016746695238),
        T(-3550918686646 // 2091501179385),
        T(-1275806237668 // 842570457699),
    )

    RKB = (
        T(1432997174477 // 9575080441755),
        T(5161836677717 // 13612068292357),
        T(1720146321549 // 2090206949498),
        T(3134564353537 // 4481467310338),
        T(2277821191437 // 14882151754819),
    )

    RKC = (
        T(0),
        T(1432997174477 // 9575080441755),
        T(2526269341429 // 6820363962896),
        T(2006345519317 // 3224310063776),
        T(2802321613138 // 2924317926251),
    )

    nstep = ceil(Int, (t1 - t0) / dt)
    dt = (t1 - t0) / nstep

    fill!(Δq, 0)
    fill!(Δq2, 0)
    for step in 1:nstep
        t = t0 + (step - 1) * dt
        for s in 1:length(RKA)
            f!(Δq2, q, p, t + RKC[s] * dt)
            Δq .+= Δq2
            q .+= RKB[s] * dt * Δq
            Δq .*= RKA[s % length(RKA) + 1]
        end
    end

    nothing
end

function bisection(F, a::T, b::T, tol = 100 * eps(T)) where T
    (a, b) = a < b ? (a, b) : (b, a)

    Fa = F(a)
    Fb = F(b)
    if abs(Fa) < abs(Fb)
        abs(Fa) < eps(T) && return a
    else
        abs(Fb) < eps(T) && return b
    end
    @assert Fa * Fb < 0

    V = (a + b) / 2
    while b - a > tol
        FV = F(V)
        if FV == 0
            break
        elseif FV * Fa < 0
            b = V
        else
            a = V
        end
        V = (a + b) / 2
    end

    return V
end

function semianalytic_solution(u0, v0, σ0, tend::T, F) where T

    # Characteristics propagating into interface
    w_m(x, t) = v0(x - t) - σ0(x - t)
    w_p(x, t) = v0(x + t) + σ0(x + t) # zero

    # Find V by root finding
    η(t) = w_p(0, t) - w_m(0, t)
    g(V, t) = 2F(V) + V - η(t)
    V(t) = t > 0 ? bisection(V -> g(V, t), zero(T), η(t)) : zero(t)

    # Set Characteristics propagating out of interface
    Q_m(t) = w_m(0, t) + 2F(V(t))
    Q_p(t) = w_p(0, t) - 2F(V(t))

    # Integrate variables in order to compute displacement
    ∫Q_m_ = solve(ODEProblem((_, _, t) -> Q_m(t), zero(T), (zero(T), tend)),
                  Tsit5(),
                  dtmax = 0.1,
                  abstol=1e-12,
                  reltol=1e-12)
    ∫Q_m(t) = t > 0 ? ∫Q_m_(t) : zero(T)
    ∫Q_p_ = solve(ODEProblem((_, _, t) -> Q_p(t), zero(T), (zero(T), tend)),
                  Tsit5(),
                  dtmax = 0.1,
                  abstol=1e-12,
                  reltol=1e-12)
    ∫Q_p(t) = t > 0 ? ∫Q_p_(t) : zero(T)

    ∫q_m(x, t) = ∫Q_m(t + x)
    ∫q_p(x, t) = ∫Q_p(t - x)

    q_m(x, t) = Q_m(t + x)
    q_p(x, t) = Q_p(t - x)

    u_m(x, t) = u0(x - t) + ∫q_m(x, t) / 2
    v_m(x, t) = (q_m(x, t) + w_m(x, t)) / 2
    σ_m(x, t) = (q_m(x, t) - w_m(x, t)) / 2

    u_p(x, t) = u0(x + t) + ∫q_p(x, t) / 2
    v_p(x, t) =  (q_p(x, t) + w_p(x, t)) / 2
    σ_p(x, t) = -(q_p(x, t) - w_p(x, t)) / 2

    return u_m, u_p, v_m, v_p, σ_m, σ_p
end

function noncharacteristic_operator!(dq, q, D, b0, bN, HI, F)
    Nq = size(D, 1)
    rngs = ntuple(j-> ((j-1) * Nq + 1) : (j * Nq), 4)
    v_m, v_p, u_m, u_p = view.(Ref(q), rngs)

    V_m = v_p[1] - v_m[Nq]
    us_m, us_p = u_m[Nq], u_p[1]
    τs_m = F(V_m)
    τs_p = -τs_m

    operator!(dq, q, D, b0, bN, HI, τs_m, τs_p, us_m, us_p)
    nothing
end

function characteristic_operator!(dq, q, D, b0, bN, HI, F, γ)
    Nq = size(D, 1)
    rngs = ntuple(j-> ((j-1) * Nq + 1) : (j * Nq), 4)
    v_m, v_p, u_m, u_p = view.(Ref(q), rngs)

    us_m, us_p = q[4Nq + 1], q[4Nq + 2]

    # Normal
    n0, nN = -1, 1

    τ_m = nN * bN' * u_m + γ * (us_m - u_m[Nq])
    τ_p = n0 * b0' * u_p + γ * (us_p - u_p[ 1])

    w_m = v_m[Nq] - τ_m
    w_p = v_p[ 1] - τ_p

    η_m = w_p - w_m
    g(V_m) = 2F(V_m) + V_m - η_m
    V_m = bisection(g, zero(η_m), η_m)
    τs_m = F(V_m)
    τs_p = -τs_m

    operator!(dq, q, D, b0, bN, HI, τs_m, τs_p, us_m, us_p)

    dq[4Nq + 1] = w_m + τs_m
    dq[4Nq + 2] = w_p + τs_p

    nothing
end

function operator!(dq, q, D, b0, bN, HI, τs_m, τs_p, us_m, us_p)
    Nq = size(D, 1)
    rngs = ntuple(j-> ((j-1) * Nq + 1) : (j * Nq), 4)
    dv_m, dv_p, du_m, du_p = view.(Ref(dq), rngs)
    v_m, v_p, u_m, u_p = view.(Ref(q), rngs)

    # Set the displacement rates
    du_m .= v_m
    du_p .= v_p

    # Set the velocity rates
    dv_m .= D * u_m
    dv_p .= D * u_p

    # Normal
    n0, nN = -1, 1

    # Set the Neumann bcs
    dv_m[ 1] -= n0 * HI[ 1, 1] * b0' * u_m
    dv_p[Nq] -= nN * HI[Nq,Nq] * bN' * u_p

    # Set the interface values
    dv_m[Nq] += HI[Nq, Nq] * (τs_m - nN * bN' * u_m)
    dv_p[ 1] += HI[ 1,  1] * (τs_p - n0 * b0' * u_p)
    dv_m .-= nN * HI * bN * (us_m - u_m[Nq])
    dv_p .-= n0 * HI * b0 * (us_p - u_p[ 1])
    nothing
end

function convergence_test(p, Ns, β, cfl_c, cfl_nc, output = false)
    T = Float64

    x_m = range(T(-1), stop = 0, length=100)
    x_p = range(T(0), stop = 1, length=100)

    # Pulse center
    μ = -1 / 2

    # Pulse width
    w = 1 / 15

    u0(x) = exp(-((x - μ) / w)^2)
    v0(x) = derivative(t->u0(x-t), 0)
    σ0(x) = derivative(x->u0(x), x)

    tend::T = 1

    F(V) = β * asinh(V)

    ue_m, ue_p, ve_m, ve_p, σe_m, σe_p = semianalytic_solution(u0, v0, σ0, tend, F)

    # storage for error
    err_nc  = zeros(length(Ns))
    err_c   = zeros(length(Ns))

    for (iter, N) = enumerate(Ns)
        if output
            @info @sprintf " level %d with N = %4d" iter N
        end
        Nq = N + 1

        # Setup the SBP operators
        (D, S0, SN, HI, H, x_p) = DiagonalSBP.D2(p, N; xc = (0, 1))
        b0, bN = S0[1, :], SN[end, :]
        θ_R = DiagonalSBP.D2_remainder_parameters(p).θ_R
        θ_H = DiagonalSBP.D2_remainder_parameters(p).θ_H
        γ = N * (1 / θ_R + 1 / θ_H)

        # Create the grid
        x_m = x_p .- 1
        x = [x_m; x_p]

        # Create the solution
        q0 = zeros(T, 4Nq)
        q0[1:2Nq] = v0.(x)
        q0[2Nq .+ (1:2Nq)] = u0.(x)
        qnc = copy(q0)

        fnc!(Δq, q, _, _) = noncharacteristic_operator!(Δq, q, D, b0, bN, HI, F)

        qc = [q0; u0(0); u0(0)]
        fc!(Δq, q, _, _) = characteristic_operator!(Δq, q, D, b0, bN, HI, F, γ)

        h = 1 / N
        dtnc = cfl_nc * h
        dtc = cfl_c * h

        timestep!(qnc, fnc!, nothing, dtnc, (T(0), tend))
        timestep!(qc, fc!, nothing, dtc, (T(0), tend))

        rngs = ntuple(j-> ((j-1) * Nq + 1) : (j * Nq), 4)
        _, _, unc_m, unc_p = view.(Ref(qnc), rngs)
        _, _, uc_m, uc_p = view.(Ref(qc), rngs)
        err_c[iter] = sqrt(sum(H * ((uc_m - ue_m.(x_m, tend)).^2 +
                                    (uc_p - ue_p.(x_p, tend)).^2)))
        err_nc[iter] = sqrt(sum(H * ((unc_m - ue_m.(x_m, tend)).^2 +
                                     (unc_p - ue_p.(x_p, tend)).^2)))
        if output
            if iter > 1
                rate = ((log(err_nc[iter-1]) - log(err_nc[iter])) /
                        (log(Ns[iter]) - log(Ns[iter-1])))
                @info @sprintf """ Non-Characteristic
                error = %.2e
                rate  = %.2e""" err_nc[iter] rate
            else
                @info @sprintf """ Non-Characteristic
                error = %.2e""" err_nc[iter]
            end
            if iter > 1
                rate = ((log(err_c[iter-1]) - log(err_c[iter])) /
                        (log(Ns[iter]) - log(Ns[iter-1])))
                @info @sprintf """ Characteristic
                error = %.2e
                rate  = %.2e""" err_c[iter] rate
            else
                @info @sprintf """ Characteristic
                error = %.2e""" err_c[iter]
            end
        end
    end
    return err_c, err_nc
end

let
    Ns = 17 * 2 .^ (0:5)

    hs = 1 ./ Ns
    tri_x = [hs[end-1], hs[end], hs[end-1], hs[end-1]]
    tri_y = [hs[end-1], hs[end], hs[end  ], hs[end-1]] ./ hs[end-1]

    cfl_nc = Dict()
    cfl_nc[2, 1] = 1 // 2
    cfl_nc[4, 1] = 1 // 2
    cfl_nc[6, 1] = 1 // 2

    cfl_nc[2, 2] = 1 // 2
    cfl_nc[4, 2] = 1 // 4
    cfl_nc[6, 2] = 1 // 4

    cfl_nc[2, 32] = 1 // 32
    cfl_nc[4, 32] = 1 // 64
    cfl_nc[6, 32] = 1 // 64

    cfl_nc[2, 64] = 1 // 64
    cfl_nc[4, 64] = 1 // 128
    cfl_nc[6, 64] = 1 // 128

    cfl_nc[2, 128] = 1 // 128
    cfl_nc[4, 128] = 1 // 256
    cfl_nc[6, 128] = 1 // 256
    for β in (32, 64, 128)
        @show β

        C2, NC2 = convergence_test(2, Ns, β, 1 // 2, cfl_nc[2, β])
        C4, NC4 = convergence_test(4, Ns, β, 1 // 2, cfl_nc[4, β])
        C6, NC6 = convergence_test(6, Ns, β, 1 // 4, cfl_nc[6, β])

        @pgf a = LogLogAxis(
                            {
                             ymax = 1,
                             ymin = 1e-11,
                             xmin = hs[end]/2,
                             xmax = 2hs[1],
                             xlabel = "\$h\$",
                             ylabel = "\$\\|\\Delta u\\|_{H}\$",
                            }
                           )
        tc2 = 0.6*NC2[end-1] * tri_y.^2
        @pgf push!(a, Plot(
                           {
                            color = "black",
                           },
                           Table(tri_x, tc2)
                          )
                  )
        tc4 = 0.6*NC4[end-1] * tri_y.^4
        @pgf push!(a, Plot(
                           {
                            color = "black",
                           },
                           Table(tri_x, 0.6*NC4[end-1] * tri_y.^4)
                          )
                  )
        tc6 = 0.6*NC6[end-1] * tri_y.^6
        @pgf push!(a, Plot(
                           {
                            color = "black",
                           },
                           Table(tri_x, 0.6*NC6[end-1] * tri_y.^5)
                          )
                  )
        #=
        @pgf push!(a,
        [@sprintf("""
        \\node[anchor=north] () at (axis cs:%e, %e){\$s = 2\$};
        \\node[anchor=north] () at (axis cs:%e, %e){\$s = 4\$};
        \\node[anchor=north] () at (axis cs:%e, %e){\$s = 5\$};
        """,
        (hs[end-1] + hs[end]) / 2, tc2[2],
        (hs[end-1] + hs[end]) / 2, tc4[2],
        (hs[end-1] + hs[end]) / 2, tc6[2],
        )])
        =#
        @pgf push!(a,
                   [@sprintf("""
                             \\node[anchor=west] () at (axis cs:%e, %e){\$2\$};
                             \\node[anchor=west] () at (axis cs:%e, %e){\$4\$};
                             \\node[anchor=west] () at (axis cs:%e, %e){\$5\$};
                             """,
                             hs[end-1], exp((log(tc2[1]) + log(tc2[2]))/2),
                             hs[end-1], exp((log(tc4[1]) + log(tc4[2]))/2),
                             hs[end-1], exp((log(tc6[1]) + log(tc6[2]))/2),
                            )])

        @pgf push!(a, Plot(
                           {
                            color = "red",
                            mark = "o",
                           },
                           Table(hs, C2)
                          )
                  )
        @pgf push!(a, Plot(
                           {
                            color = "red",
                            mark = "x",
                           },
                           Table(hs, NC2)
                          )
                  )
        @pgf push!(a, Plot(
                           {
                            color = "blue",
                            mark = "o",
                           },
                           Table(hs, C4)
                          )
                  )
        @pgf push!(a, Plot(
                           {
                            color = "blue",
                            mark = "x",
                           },
                           Table(hs, NC4)
                          )
                  )
        @pgf push!(a, Plot(
                           {
                            color = "green",
                            mark = "o",
                           },
                           Table(hs, C6)
                          )
                  )
        @pgf push!(a, Plot(
                           {
                            color = "green",
                            mark = "x",
                           },
                           Table(hs, NC6)
                          )
                  )
        pgfsave("convergence_interface_$β.tex", a)
        pgfsave("convergence_interface_$β.pdf", a)
    end
end


text:
using Distributed

#Need to remove all procs before running these tests
addprocs(4)

@everywhere using DistributedArrays, DistributedJets, JSON, Jets, LinearAlgebra, Statistics, Test

@everywhere JopFoo_df!(d,m;diagonal,kwargs...) = d .= diagonal .* m
@everywhere function JopFoo(diag)
    spc = JetSpace(Float64, size(diag))
    JopLn(;df! = JopFoo_df!, df′! = JopFoo_df!, dom = spc, rng = spc, s = (diagonal=diag,))
end

@everywhere JopClose_df!(d,m;diagonal,kwargs...) = d .= diagonal .* m
@everywhere function JopClose(diag)
    spc = JetSpace(Float64, size(diag))
    file = touch(tempname())
    JopLn(;df! = JopClose_df!, dom = spc, rng = spc, s = (diagonal=diag, file=file))
end
@everywhere Base.close(J::Jet{D,R,typeof(JopClose_df!)}) where {D,R} = rm(state(J).file)

@everywhere JopBar_f!(d,m) = d .= m.^2
@everywhere JopBar_df!(δd,δm;mₒ,kwargs...) = δd .= 2 .* mₒ .* δm
@everywhere function JopBar(n)
    spc = JetSpace(Float64, n)
    JopNl(f! = JopBar_f!, df! = JopBar_df!, df′! = JopBar_df!, dom = spc, rng = spc)
end
@everywhere Jets.perfstat(J::T) where {D,R,T<:Jets.Jet{D,R,typeof(JopBar_f!)}} = Float64(π)

@everywhere JopBaz_df!(d,m;A,kwargs...) = d .= A*m
@everywhere JopBaz_df′!(m,d;A,kwargs...) = m .= A'*d
@everywhere function JopBaz(A)
    dom = JetSpace(eltype(A), size(A,2))
    rng = JetSpace(eltype(A), size(A,1))
    JopLn(;df! = JopBaz_df!, df′! = JopBaz_df′!, dom = dom, rng = rng, s = (A=A,))
end

@testset "DArray irregular construction, T=$T" for T in (Float32,Float64,Complex{Float32},Complex{Float64})
    A = DArray(I->myid()*ones(T,length(I[1]),length(I[2])), workers()[1:2], [1:2,3:10], [1:2])
    @test size(A) == (10,2)
    @test A.indices[1] == (1:2, 1:2)
    @test A.indices[2] == (3:10, 1:2)
    @test indices(A,1) == [1:2,3:10]
    @test indices(A,2) == [1:2]
    @test fetch(@spawnat procs(A)[1] all(a->a≈T(myid()), localpart(A)))
    @test fetch(@spawnat procs(A)[2] all(a->a≈T(myid()), localpart(A)))
    A = DArray(I->myid()*ones(length(I[1])), workers(), [1:2,3:10])
    @test size(A) == (10,)
    @test A.indices[1] == (1:2,)
    @test A.indices[2] == (3:10,)
    @test indices(A,1) == [1:2,3:10]
    @test indices(A) == indices(A,1)
end

@testset "JetDSpace construction from JetBSpace array" begin
    _blkspaces = [JetBSpace([JetSpace(Float64,2),JetSpace(Float64,2)]), JetBSpace([JetSpace(Float64,2),JetSpace(Float64,2)])]
    blkspaces = distribute(_blkspaces; procs=workers()[1:2])
    R = JetDSpace(blkspaces)
    @test length(R) == 8
    @test nblocks(R) == 4
    @test blockmap(R) == [1:2,3:4]
    @test R.indices == [1:4,5:8]
end

@testset "JetDSpace construction, 1D arrays" begin
    A = @blockop DArray(I->[JopFoo(rand(2)) for i in I[1], j in I[2]], (2,1), workers()[1:2])
    R = range(A)
    @test size(R) == (4,)
    @test length(R) == 4
    @test eltype(R) == Float64
    @test eltype(typeof(R)) == Float64
    @test ndims(R) == 1
    @test indices(R) == [1:2,3:4]
    @test procs(R) == workers()[1:2]
    @test nprocs(R) == 2
    @test nprocs(A) == 2
    @test nblocks(R) == 2
end

@testset "JetDSpace construction, 2D arrays" begin
    A = @blockop DArray(I->[JopFoo(rand(2,3)) for i in I[1], j in I[2]], (2,1), workers()[1:2])
    R = range(A)
    @test size(R) == (12,)
    @test length(R) == 12
    @test eltype(R) == Float64
    @test eltype(typeof(R)) == Float64
    @test ndims(R) == 1
    @test indices(R) == [1:6,7:12]
    @test procs(R) == workers()[1:2]
    @test nprocs(R) == 2
    @test nprocs(A) == 2
    @test nblocks(R) == 2
end

@testset "JetDSpace construction from DBArray" begin
    A = @blockop DArray(I->[JopFoo(rand(2,3)) for i in I[1], j in I[2]], (2,1), workers()[1:2])
    R = range(A)
    @test R == space(rand(R))

    A = @blockop DArray(I->[JopFoo(rand(2)) for i in I[1], j in I[2]], (2,1), workers()[1:2])
    R = range(A)
    @test R == space(rand(R))
end

@testset "JetDSpace operations, 1D arrays" begin
    A = @blockop DArray(I->[JopFoo(rand(2)) for i in I[1], j in I[2]], (2,1), workers()[1:2])
    R = range(A)
    @test dzeros((4,), workers()[1:2]) ≈ zeros(R)
    @test dones((4,), workers()[1:2]) ≈ ones(R)
    d = rand(R)
    _d = drand((4,), workers()[1:2])
    @test size(d) == size(_d)
    @test d.darray.cuts == _d.cuts
    @test d.darray.indices == _d.indices
    d = Array(R)
    @test size(d) == size(_d)
    @test d.darray.cuts == _d.cuts
    @test d.darray.indices == _d.indices
    @test indices(d) == [1:2,3:4]

    y = similar(d)
    @test isa(y, DistributedJets.DBArray{Float64})
    z = similar(d, Float32)
    @test isa(z, DistributedJets.DBArray{Float32})
    _z = similar(d, Float32, length(y))
    @test isa(_z, DistributedJets.DBArray{Float32})
    __z = similar(d, Float32, (length(y),))
    @test isa(__z, DistributedJets.DBArray{Float32})

    @test remotecall_fetch(localindices, workers()[1], R) == 1:2
    @test remotecall_fetch(localindices, workers()[2], R) == 3:4
    @test remotecall_fetch(localblockindices, workers()[1], R) == 1:1
    @test remotecall_fetch(localblockindices, workers()[2], R) == 2:2

    x = getblock(d,1)
    x .= π
    setblock!(d,1,x)
    @test d[1:2] ≈ [π,π]
    x .= 0
    @test remotecall_fetch(getblock!, workers()[1], d, 1, x) ≈ [π,π]
end

@testset "JetDSpace operations, 2D arrays" begin
    A = @blockop DArray(I->[JopFoo(rand(2,3)) for i in I[1], j in I[2]], (4,1), workers()[1:2])
    R = range(A)
    @test dzeros((24,), workers()[1:2]) ≈ zeros(R)
    @test dones((24,), workers()[1:2]) ≈ ones(R)
    d = rand(R)
    _d = drand((24,), workers()[1:2])
    @test size(d) == size(_d)
    @test d.darray.cuts == _d.cuts
    @test d.darray.indices == _d.indices
    d = Array(R)
    @test size(d) == size(_d)
    @test d.darray.cuts == _d.cuts
    @test d.darray.indices == _d.indices
    @test blockmap(d) == [1:2,3:4]

    @test remotecall_fetch(localindices, workers()[1], R) == 1:12
    @test remotecall_fetch(localindices, workers()[2], R) == 13:24
    @test remotecall_fetch(localblockindices, workers()[1], R) == 1:2
    @test remotecall_fetch(localblockindices, workers()[2], R) == 3:4

    x = getblock(d,1)
    x .= π
    setblock!(d,1,x)
    @test d[1:6] ≈ [π, π, π, π, π, π]
    x .= 0
    @test remotecall_fetch(getblock!, workers()[1], d, 1, x) ≈ [π π π; π π π]
end

@testset "DBArray, construction" begin
    @everywhere foo(i) = rem(i,2) == 0 ? i*π*ones(2) : i*π*ones(3)
    nblks = 7
    A = DBArray(foo, (nblks,), workers(), [4])
    B = DBArray(foo, (nblks,), workers())
    C = DBArray(foo, (nblks,))
    D = DBArray(foo, blockmap(A))
    E = DBArray(foo, [1:1,2:2,3:3,4:7])

    blockmap(E)

    E.blkindices
    E.indices

    @test A ≈ B
    @test A ≈ C
    @test A ≈ D
    for i = 1:18
        @test A[i] ≈ E[i]
    end

    @test A.blkindices == B.blkindices
    @test A.blkindices == C.blkindices
    @test A.blkindices == D.blkindices
    @test A.blkindices != E.blkindices
    @test A.indices == B.indices
    @test A.indices == C.indices
    @test A.indices == D.indices
    @test A.indices != E.indices

    for iblock in 1:nblocks(A)
        a = getblock(A, iblock)
        b = getblock(B, iblock)
        c = getblock(C, iblock)
        d = getblock(D, iblock)
        e = getblock(E, iblock)
        n = length(a)
        m = rem(iblock,2) == 0 ? 2 : 3
        @test n == m
        @test length(b) == n
        @test length(c) == n
        @test length(d) == n
        @test length(e) == n
        @test a ≈ iblock*π*ones(m)
        @test b ≈ a
        @test c ≈ a
        @test d ≈ a
        @test e ≈ a
    end
end

@testset "DBArray, inner product" begin
    A = @blockop DArray(I->[JopFoo(rand(2,3)) for i in I[1], j in I[2]], (4,1), workers()[1:2])
    d₁ = rand(range(A))
    d₂ = rand(range(A))
    @test dot(d₁, d₂) ≈ dot(convert(Array, d₁), convert(Array, d₂))
end

@testset "DBArray, norm" begin
    A = @blockop DArray(I->[JopFoo(rand(2,3)) for i in I[1], j in I[2]], (4,1))
    d = rand(range(A))
    @test norm(d) ≈ norm(convert(Array, d))
    @test norm(d, 0) ≈ norm(convert(Array, d), 0)
    @test norm(d, Inf) ≈ norm(convert(Array, d), Inf)

    _d = similar(d, Float32)
    _d .= d
    @test typeof(norm(_d, 0.5)) == Float32
end

@testset "DBArray, mapreduce" begin
    A = @blockop DArray(I->[JopFoo(rand(2,3)) for i in I[1], j in I[2]], (4,1))
    d = rand(range(A))
    @test mapreduce(abs, +, d) ≈ mapreduce(abs, +, convert(Array, d))
end

@testset "DBArray, f-maximum" begin
    A = @blockop DArray(I->[JopFoo(rand(2,3)) for i in I[1], j in I[2]], (4,1))
    d = rand(range(A))
    @test maximum(abs, d) ≈ maximum(abs, convert(Array, d))
end

@testset "DBArray, maximum" begin
    A = @blockop DArray(I->[JopFoo(rand(2,3)) for i in I[1], j in I[2]], (4,1))
    d = rand(range(A))
    @test maximum(d) ≈ maximum(convert(Array, d))
end

@testset "DBArray, f-minimum" begin
    A = @blockop DArray(I->[JopFoo(rand(2,3)) for i in I[1], j in I[2]], (4,1))
    d = rand(range(A))
    @test maximum(abs, d) ≈ maximum(abs, convert(Array, d))
end

@testset "DBArray, minimum" begin
    A = @blockop DArray(I->[JopFoo(rand(2,3)) for i in I[1], j in I[2]], (4,1))
    d = rand(range(A))
    @test minimum(d) ≈ minimum(convert(Array, d))
end

@testset "DBArray, extrema" begin
    A = @blockop DArray(I->[JopFoo(rand(2,3)) for i in I[1], j in I[2]], (4,1), workers()[1:2])
    d = rand(range(A))
    mn,mx = extrema(d)
    _mn,_mx = extrema(convert(Array,d))
    @test mn ≈ _mn
    @test mx ≈ _mx
end

@testset "DBArray, mean" begin
    A = @blockop DArray(I->[JopFoo(rand(2,3)) for i in I[1], j in I[2]], (4,1), workers()[1:2])
    d = rand(range(A))
    μ = mean(d)
    _μ = mean(convert(Array, d))
    @test μ ≈ _μ
end

@testset "DBArray, variance" begin
    A = @blockop DArray(I->[JopFoo(rand(2,3)) for i in I[1], j in I[2]], (4,1), workers()[1:2])
    d = rand(range(A))

    σ² = var(d)
    _σ² = var(convert(Array, d))
    @test σ² ≈ _σ²

    σ² = var(d; mean=100.0)
    _σ² = var(convert(Array, d); mean=100.0)
    @test σ² ≈ _σ²

    σ² = var(d; corrected=false)
    _σ² = var(convert(Array, d); corrected=false)
    @test σ² ≈ _σ²
end

@testset "DBArray broadcasting, 1D arrays" begin
    A = @blockop DArray(I->[JopFoo(rand(2)) for i in I[1], j in I[2]], (7,1), workers()[1:2], [2,1])
    R = range(A)
    R.blkindices
    R.blkspaces[1]
    d₁ = ones(R)
    d₂ = ones(R)
    d₃ = ones(R)
    d = ones(R)

    α₁ = rand(Float64)
    α₂ = rand(Float64)
    α₃ = rand(Float64)

    d = α₃*d₃

    d .= α₁*d₁ .+ α₂*d₂ .+ α₃*d₃
    α = α₁ + α₂ + α₃

    for i = 1:length(d)
        @test d[i] ≈ α
    end

    d₁ .= 3.14

    for i = 1:length(d₁)
        @test d₁[i] ≈ 3.14
    end
end

@testset "DBArray broadcasting, 2D arrays" begin
    A = @blockop DArray(I->[JopFoo(rand(2,3)) for i in I[1], j in I[2]], (7,1), workers()[1:2], [2,1])
    R = range(A)
    R.blkindices
    R.blkspaces[1]
    d₁ = ones(R)
    d₂ = ones(R)
    d₃ = ones(R)
    d = ones(R)

    α₁ = rand(Float64)
    α₂ = rand(Float64)
    α₃ = rand(Float64)

    d = α₃*d₃

    d .= α₁*d₁ .+ α₂*d₂ .+ α₃*d₃
    α = α₁ + α₂ + α₃

    for i = 1:length(d)
        @test d[i] ≈ α
    end

    d₁ .= 3.14
    for i = 1:length(d₁)
        @test d₁[i] ≈ 3.14
    end
end

@everywhere function myblocks(i,j)
    if i ∈ (1,2,3) && j ∈ (2,3,1)
        return JopBar(10)
    elseif i == (2,3) && j == 4
        JopZeroBlock(JetSpace(Float64,10),JetSpace(Float64,10))
    elseif i == 2 && j == 4
        JopBaz(rand(10,10)) ∘ JopBar(10)
    else
        return JopBaz(rand(10,10))
    end
end

@testset "JopDBlock, heterogeneous, tall and skinny" begin
    _F = DArray(I->[myblocks(i,j) for i in I[1], j in I[2]], (3,4), workers()[1:2], [2,1])
    F = @blockop _F

    @test nblocks(F,1) == 3
    @test nblocks(F,2) == 4
    @test nblocks(F) == (3,4)
    @test blockmap(F)[1,1] == (1:2,1:4)
    @test blockmap(F)[2,1] == (3:3,1:4)

    @test isa(F, JopNl{<:Jet{<:JetBSpace,<:DistributedJets.JetDSpace,typeof(DistributedJets.JetDBlock_f!)}})

    F₁ = remotecall_fetch(localpart, workers()[1], F)
    F₂ = remotecall_fetch(localpart, workers()[2], F)

    @test remotecall_fetch(localblockindices, workers()[1], F) == (1:2,1:4)
    @test remotecall_fetch(localblockindices, workers()[2], F) == (3:3,1:4)
    @test remotecall_fetch(localblockindices, workers()[1], F, 1) == 1:2
    @test remotecall_fetch(localblockindices, workers()[1], F, 2) == 1:4

    @test isa(F₁, JopNl{<:Jet{<:JetBSpace,<:JetBSpace,typeof(Jets.JetBlock_f!)}})

    _G = [_F[i,j] for i in 1:3, j in 1:4]
    G = @blockop _G

    @test procs(F) == procs(state(F).ops)

    @test ones(range(F)) ≈ DArray(I->ones(length(I[1])), procs(_F), [1:20,21:30])
    @test ones(domain(F)) ≈ ones(40)
    @test zeros(range(F)) ≈ DArray(I->zeros(length(I[1])), procs(_F), [1:20,21:30])
    @test zeros(domain(F)) ≈ zeros(40)
    @test size(rand(domain(F))) == (40,)
    x = rand(range(F))
    @test size(x) == (30,)
    @test x.darray.cuts == [[1,21,31]]
    x = Array(range(F))
    @test size(x) == (30,)
    @test x.darray.cuts == [[1,21,31]]

    m = rand(domain(F))

    @test collect(F*m) ≈ G*m
    @test collect(F*m) ≈ [F₁*m ; F₂*m]

    J = jacobian!(F, m)
    _J = jacobian!(G, m)

    δm = rand(domain(J))
    @test collect(J*δm) ≈ _J*δm

    δd = rand(range(J))
    @test J'*δd ≈ _J'*collect(δd)

    F₃₂ = getblock(JopNl,F,3,2)
    G₃₂ = _G[3,2]
    @test F₃₂*m[21:30] ≈ G₃₂*m[21:30]
    A₁₄ = getblock(JopLn,F,1,4)
    B₁₄ = _G[1,4]
    @test A₁₄*m[31:40] ≈ B₁₄*m[31:40]
end

@testset "JopDBlock, homogeneous, tall and skinny" begin
    _F = DArray(I->[JopBar(10) for i in I[1], j in I[2]], (3,4), workers()[1:2], [2,1])
    F = @blockop _F

    _G = [_F[i,j] for i in 1:3, j in 1:4]
    G = @blockop _G

    @test isa(F, JopNl{<:Jet{<:JetBSpace,<:DistributedJets.JetDSpace,typeof(DistributedJets.JetDBlock_f!)}})

    @test ones(range(F)) ≈ DArray(I->ones(length(I[1])), procs(_F), [1:20,21:30])
    @test ones(domain(F)) ≈ ones(40)
    @test zeros(range(F)) ≈ DArray(I->zeros(length(I[1])), procs(_F), [1:20,21:30])
    @test zeros(domain(F)) ≈ zeros(40)
    @test size(rand(domain(F))) == (40,)
    x = rand(range(F))
    @test size(x) == (30,)
    @test x.darray.cuts == [[1,21,31]]
    x = Array(range(F))
    @test size(x) == (30,)
    @test x.darray.cuts == [[1,21,31]]

    m = rand(domain(F))
    @test collect(F*m) ≈ G*m

    J = jacobian!(F, m)
    _J = jacobian!(G, m)

    δm = rand(domain(J))
    @test collect(J*δm) ≈ _J*δm

    δd = rand(range(J))
    @test J'*δd ≈ _J'*collect(δd)
end

@testset "JopDBlock, distributed->distributed, block diagonal, square blocks" begin
    _F = DArray(I->[i==j ? JopBar(10) : JopZeroBlock(JetSpace(Float64,10),JetSpace(Float64,10)) for i in I[1], j in I[2]], (4,4), workers(), [4,1])
    F = @blockop _F isdiag=true

    _G = [_F[i,j] for i in 1:4, j in 1:4]
    G = @blockop _G

    @test isa(F, JopNl{<:Jet{<:DistributedJets.JetDSpace,<:DistributedJets.JetDSpace,typeof(DistributedJets.JetDBlock_f!)}})

    @test ones(range(F)) ≈ DArray(I->ones(length(I[1])), procs(_F), [1:10,11:20,21:30,31:40])
    @test ones(domain(F)) ≈ DArray(I->ones(length(I[1])), procs(_F), [1:10,11:20,21:30,31:40])

    m = rand(domain(F))
    @test collect(F*m) ≈ G*collect(m)

    J = jacobian!(F, m)
    _J = jacobian!(G, collect(m))

    δm = rand(domain(F))
    @test collect(J*δm) ≈ _J*collect(δm)
    δd = rand(range(J))
    @test collect(J'*δd) ≈ _J'*collect(δd)
end

@testset "JopDBlock, distributed->distributed, block diagonal, tall blocks" begin
    _A = DArray(I->[i==j ? JopBaz(rand(10,5)) : JopZeroBlock(JetSpace(Float64,5),JetSpace(Float64,10)) for i in I[1], j in I[2]], (4,4), workers(), [4,1])
    A = @blockop _A isdiag=true

    _B = [_A[i,j] for i in 1:4, j in 1:4]
    B = @blockop _B

    @test isa(A, JopLn{<:Jet{<:DistributedJets.JetDSpace,<:DistributedJets.JetDSpace,typeof(DistributedJets.JetDBlock_f!)}})

    @test ones(range(A)) ≈ DArray(I->ones(length(I[1])), procs(_A), [1:10,11:20,21:30,31:40])
    @test ones(domain(A)) ≈ DArray(I->ones(length(I[1])), procs(_A), [1:5,6:10,11:15,16:20])

    m = rand(domain(A))
    @test collect(A*m) ≈ B*collect(m)
    d = rand(range(A))
    @test collect(A'*d) ≈ B'*collect(d)
end

@testset "JopDBlock, distributed->distributed, block diagonal, number of blocks different from number of procs" begin
    _A = DArray(I->[i==j ? JopBaz(rand(10,5)) : JopZeroBlock(JetSpace(Float64,5),JetSpace(Float64,10)) for i in I[1], j in I[2]], (10,10), workers(), [4,1])
    A = @blockop _A isdiag=true

    _B = [_A[i,j] for i in 1:10, j in 1:10]
    B = @blockop _B

    @test isa(A, JopLn{<:Jet{<:DistributedJets.JetDSpace,<:DistributedJets.JetDSpace,typeof(DistributedJets.JetDBlock_f!)}})

    m = rand(domain(A))
    @test collect(A*m) ≈ B*collect(m)
    d = rand(range(A))
    @test collect(A'*d) ≈ B'*collect(d)
end

@testset "JopDBlock, heterogeneous, distributed->distributed" begin
    _F = DArray(I->Jop[myblocks(i,j) for i in I[1], j in I[2]], (3,5), workers(), [2,2])
    F = @blockop _F

    @test nblocks(F,1) == 3
    @test nblocks(F,2) == 5
    @test nblocks(F) == (3,5)
    @test blockmap(F)[1,1] == (1:2,1:3)
    @test blockmap(F)[2,1] == (3:3,1:3)
    @test blockmap(F)[1,2] == (1:2,4:5)
    @test blockmap(F)[2,2] == (3:3,4:5)

    @test isa(F, JopNl{<:Jet{<:DistributedJets.JetDSpace,<:DistributedJets.JetDSpace,typeof(DistributedJets.JetDBlock_f!)}})

    @test procs(F) == [workers()[1] workers()[3]; workers()[2] workers()[4]]

    @test remotecall_fetch(localblockindices, workers()[1], F) == (1:2,1:3)
    @test remotecall_fetch(localblockindices, workers()[2], F) == (3:3,1:3)
    @test remotecall_fetch(localblockindices, workers()[3], F) == (1:2,4:5)
    @test remotecall_fetch(localblockindices, workers()[4], F) == (3:3,4:5)

    @test remotecall_fetch(localblockindices, workers()[1], F, 1) == 1:2
    @test remotecall_fetch(localblockindices, workers()[1], F, 2) == 1:3

    F₁₁ = remotecall_fetch(localpart, workers()[1], F)
    @test isa(F₁₁, JopNl{<:Jet{<:JetBSpace,<:JetBSpace,typeof(Jets.JetBlock_f!)}})

    _G = [_F[i,j] for i in 1:3, j in 1:5]
    G = @blockop _G

    @test procs(F) == procs(state(F).ops)

    @test ones(range(F)) ≈ DArray(I->ones(length(I[1])), procs(_F), [1:20,21:30])
    @test ones(domain(F)) ≈ DArray(I->ones(length(I[1])), procs(_F), [1:30,31:50])
    @test zeros(range(F)) ≈ DArray(I->zeros(length(I[1])), procs(_F), [1:20,21:30])
    @test zeros(domain(F)) ≈ DArray(I->zeros(length(I[1])), procs(_F), [1:30,31:50])
    @test size(rand(domain(F))) == (50,)
    x = rand(range(F))
    @test size(x) == (30,)
    @test x.darray.cuts == [[1,21,31]]
    x = Array(range(F))
    @test size(x) == (30,)
    @test x.darray.cuts == [[1,21,31]]

    m = rand(domain(F))
    @test collect(F*m) ≈ G*collect(m)

    J = jacobian!(F, m)
    _J = jacobian!(G, m)

    δm = rand(domain(J))
    @test collect(J*δm) ≈ _J*collect(δm)

    δd = rand(range(J))
    @test J'*δd ≈ _J'*collect(δd)

    F₃₂ = getblock(JopNl,F,3,2)
    G₃₂ = _G[3,2]
    @test F₃₂*m[21:30] ≈ G₃₂*m[21:30]
    A₁₄ = getblock(JopLn,F,1,4)
    B₁₄ = _G[1,4]
    @test A₁₄*m[31:40] ≈ B₁₄*m[31:40]
end

@testset "JetDBlock, localpart" begin
    F = @blockop DArray(I->[JopBar(10) for i in I[1], j in I[2]], (3,4), workers()[1:2], [2,1])
    m = rand(domain(F))
    J = jacobian!(F, m)
    _J = remotecall_fetch(localpart, workers()[1], J)
    @test isa(_J, JopLn)
    for jblock = 1:nblocks(_J,2), iblock = 1:nblocks(_J,1)
        @test isa(getblock(_J,iblock,jblock), JopLn)
    end

    __J = @blockop [jacobian(JopBar(10), m[1:10]) jacobian(JopBar(10), m[11:20]) jacobian(JopBar(10), m[21:30]) jacobian(JopBar(10), m[31:40]);
           jacobian(JopBar(10), m[1:10]) jacobian(JopBar(10), m[11:20]) jacobian(JopBar(10), m[21:30]) jacobian(JopBar(10), m[31:40])]
    @test __J*m ≈ _J*m

    JT = J'
    _JT = remotecall_fetch(localpart, workers()[1], JT)

    typeof(_JT)
    getblock(JT, 1, 1)

    for jblock = 1:nblocks(JT,2), iblock = 1:nblocks(JT,1)
        @test isa(getblock(JT,iblock,jblock), JopAdjoint)
    end
end

@testset "JopDBlock, statistics reporting" begin
    rm("stats.json", force=true)

    _F = DArray(I->[JopBar(10) for i in I[1], j in I[2]], (3,4), workers(), [2,1])
    F = @blockop _F perfstatfile="stats.json"

    m = rand(domain(F))
    d = F*m

    s = JSON.parse(read("stats.json", String))
    @test s["step"][1]["pid"][1]["localblock"] ≈ π*ones(8)
    @test s["step"][1]["pid"][2]["localblock"] ≈ π*ones(4)
    @test s["step"][1]["pid"][1]["operation"] == "f"
    @test s["step"][1]["pid"][2]["operation"] == "f"

    J = jacobian(F, m)
    d = J*m

    s = JSON.parse(read("stats.json", String))
    @test s["step"][2]["pid"][1]["localblock"] ≈ π*ones(8)
    @test s["step"][2]["pid"][2]["localblock"] ≈ π*ones(4)
    @test s["step"][2]["pid"][1]["operation"] == "df"
    @test s["step"][2]["pid"][2]["operation"] == "df"

    m = J'*d

    s = JSON.parse(read("stats.json", String))
    @test s["step"][3]["pid"][1]["localblock"] ≈ π*ones(8)
    @test s["step"][3]["pid"][2]["localblock"] ≈ π*ones(4)
    @test s["step"][3]["pid"][1]["operation"] == "df′"
    @test s["step"][3]["pid"][2]["operation"] == "df′"

    rm("stats.json", force=true)
end

@testset "block operator composed with operator" begin
    A₁ = JopFoo(rand(2))
    A₂ = @blockop DArray(I->[JopBar(2) for i in I[1], j in I[2]], (2,1))
    A = A₂ ∘ A₁

    m = rand(domain(A))
    d = A*m
    A₁₁ = getblock(A, 1, 1)
    A₂₁ = getblock(A, 2, 1)
    @test getblock(d,1) ≈ A₂₁ * m
    @test getblock(d,2) ≈ A₂₁ * m
end

@testset "block operators with close" begin
    A = @blockop DArray(I->[JopClose(2) for i in I[1], j in I[2]], (2,1))
    A₁ = getblock(A,1,1)
    A₂ = getblock(A,2,1)
    @test isfile(state(A₁).file)
    @test isfile(state(A₂).file)
    close(A)
    @test !isfile(state(A₁).file)
    @test !isfile(state(A₂).file)
end

@testset "localblockindices for composite operators" begin
    A = @blockop DArray(I->[JopBar(2) for i in I[1], j in I[2]], (2,1))
    B = JopBar(2)
    C = A ∘ B
    @test_throws RemoteException remotecall_fetch(localblockindices, workers()[1], C)
end

@testset "vectorized operators" begin
    A = @blockop DArray(I->[JopFoo(rand(2,3)) for i in I[1], j in I[2]], (2,1))
    x = rand(domain(A))

    @test size(vec(domain(A))) == (2*3,)
    @test size(vec(range(A))) == (2*2*3,)
    @test size(range(A)) == size(vec(range(A)))

    B = vec(A)
    @test isa(jet(B), Jet{<:JetAbstractSpace,<:JetAbstractSpace,typeof(Jets.JetVec_f!)})

    d = A * x
    _d = vec(A) * vec(x)
    @test d ≈ _d

    a = A' * d
    _a = vec(A') * d
    __a = reshape(_a, domain(A))
    @test a ≈ __a
end

rmprocs(workers())


text:
using Base.Test
include("../SDE.jl")
using SDE
#is called with arg "travis" from travis
test_type = length(ARGS) == 1 ? ARGS[1] : ""

include("testsde.jl")
include("testlinproc.jl")
include("testtc0.jl")
include("testtc.jl")


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

using JuliaFEM
using JuliaFEM.Testing

ALL_ELEMENTS = [
    Seg2, Seg3,
    Tri3, Tri6, Tri7,
    Quad4, Quad8, Quad9,
    Tet4, Tet10,
    Wedge6,
    Hex8, Hex20, Hex27
]

@testset "Evaluating basis" begin
    for T in ALL_ELEMENTS
        el = Element(T)
        nnodes = length(el)
        for (i, X) in enumerate(get_reference_coordinates(T))
            Ni = vec(el(X))
            expected = zeros(nnodes)
            expected[i] = 1.0
            @test isapprox(Ni, expected)
        end
    end
end

function get_volume{T<:AbstractElement}(::Type{T})
    X = get_reference_coordinates(T)
    element = Element(T)
    update!(element, "geometry", X)
    V = 0.0
    for ip in get_integration_points(element)
        V += ip.weight*element(ip, 0.0, Val{:detJ})
    end
    return V
end

@testset "Calculate reference element length/area/volume" begin
    @test isapprox(get_volume(Seg2),  2.0)
    @test isapprox(get_volume(Seg3),  2.0)
    @test isapprox(get_volume(Tri3),  0.5)
    @test isapprox(get_volume(Tri6),  0.5)
    @test isapprox(get_volume(Tri7),  0.5)
    @test isapprox(get_volume(Quad4), 2.0^2)
    @test isapprox(get_volume(Quad8), 2.0^2)
    @test isapprox(get_volume(Quad9), 2.0^2)
    @test isapprox(get_volume(Tet4),  1/6)
    @test isapprox(get_volume(Tet10), 1/6)
    @test isapprox(get_volume(Wedge6), 1.0)
    @test isapprox(get_volume(Hex8),  2.0^3)
    @test isapprox(get_volume(Hex20), 2.0^3)
    @test isapprox(get_volume(Hex27), 2.0^3)
end


text:
abstract type AbstractCondition{T} end

struct ConditionPrice <: AbstractCondition{PRICE}
  conjunction::String
  is_more::Bool
  value::Float64
  conId::Int
  exchange::String
  triggerMethod::Int
end

struct ConditionTime <: AbstractCondition{TIME}
  conjunction::String
  is_more::Bool
  value::String
end

struct ConditionMargin <: AbstractCondition{MARGIN}
  conjunction::String
  is_more::Bool
  value::Int
end

struct ConditionExecution <: AbstractCondition{EXECUTION}
  conjunction::String
  secType::String
  exchange::String
  symbol::String
end

struct ConditionVolume <: AbstractCondition{VOLUME}
  conjunction::String
  is_more::Bool
  value::Int
  conId::Int
  exchange::String
end

struct ConditionPercentChange <: AbstractCondition{PERCENTCHANGE}
  conjunction::String
  is_more::Bool
  value::Float64
  conId::Int
  exchange::String
end

const condition_map = Dict(PRICE         => ConditionPrice,
                           TIME          => ConditionTime,
                           MARGIN        => ConditionMargin,
                           EXECUTION     => ConditionExecution,
                           VOLUME        => ConditionVolume,
                           PERCENTCHANGE => ConditionPercentChange)


text:
# Contains all helper functions for Tmatrix module

using LinearAlgebra

"""
    Calculate dot product for two vectors.

# I made this function because the "dot" function in LinearAlgebra package doesn't work as expected for complex vectors.

Parameters

A, B : two vectors

# Returns

scalar dot product
"""
function vector_dot_product(A, B)
    return sum(A .* B)
end

"""
    Create a meshgrid, to be used with numerical integrals.

The meshgrid created is following this rule:
x_grid[dim1,dim2]
dim1 corresponds to x index
dim2 corresponds to y index
"""
function meshgrid(x, y)
    x_grid = repeat(x, 1, length(y))
    y_grid = repeat(y', length(x), 1)
    return x_grid, y_grid
end

"""
    Create meshgrid of θ and ϕ
"""
function meshgrid_θ_ϕ(
    n_θ_points,
    n_ϕ_points;
    min_θ = 1e-16,
    min_ϕ = 1e-16,
    rotationally_symmetric = false,
)
    #TODO can avoid min values w full algebra
    θ_1D_array = LinRange(min_θ, π, n_θ_points)
    ϕ_1D_array = LinRange(min_ϕ, 2π, n_ϕ_points)
    if rotationally_symmetric
        θ_array = collect(θ_1D_array)
        ϕ_array = zeros(size(θ_array))
    else
        θ_array, ϕ_array = meshgrid(θ_1D_array, ϕ_1D_array)
    end
    return θ_array, ϕ_array
end

"""
Create a single index from m and n
We fill the index like this:
`idx = 0 for n = 1:n_max     for m = -n:n global idx idx += 1 end end`
"""
function single_index_from_m_n(m::Int, n::Int)
    return n * (n + 1) + m
end

"""
    Get the maximum single index, given the maximum n.
"""
function get_max_single_index_from_n_max(n_max::Int)
    return single_index_from_m_n(n_max, n_max)
end

"""
    Multiply the integrand by the `dS` element, which equals r²sin(θ)
"""
function surface_integrand(
    integrand::AbstractVecOrMat{C},
    r_array::AbstractVecOrMat{R},
    θ_array::AbstractVecOrMat{R},
) where {R <: Real, C <: Complex{R}}
    return integrand .* r_array .^ 2 .* sin.(θ_array)
end

"""
    Multiply the integrand by the `dS` element, which equals r²sin(θ)
"""
function surface_integrand(
    integrand::AbstractVecOrMat{R},
    r_array::AbstractVecOrMat{R},
    θ_array::AbstractVecOrMat{R},
) where {R <: Real}
    return integrand .* r_array .^ 2 .* sin.(θ_array)
end

"""
    Multiply the integrand by the `dS` element, which equals r²sin(θ)
"""
function surface_integrand(
    integrand::AbstractVecOrMat{Matrix{C}},
    r_array::AbstractVecOrMat{R},
    θ_array::AbstractVecOrMat{R},
) where {R <: Real, C <: Complex{R}}
    return integrand .* r_array .^ 2 .* sin.(θ_array)
end

"""
    Multiply the integrand by the `dS` element, which equals r²sin(θ)
"""
function surface_integrand(
    integrand::AbstractVecOrMat{Matrix{R}},
    r_array::AbstractVecOrMat{R},
    θ_array::AbstractVecOrMat{R},
) where {R <: Real, C <: Complex{R}}
    return integrand .* r_array .^ 2 .* sin.(θ_array)
end

"""
    Save T-matrix to HDF5 file, with fields "Tmatrix_real_CELES_convention" and "Tmatrix_imag_CELES_convention"
"""
function save_Tmatrix_to_HDF5_file(T, HDF5_filename)
    h5write(HDF5_filename, "Tmatrix_real_CELES_convention", Float64.(real(T))) # TODO: can we store BigFloat in HDF5? I converted to `Float64` because I couldn't read the file when I stored BigFloat
    h5write(HDF5_filename, "Tmatrix_imag_CELES_convention", Float64.(imag(T)))
end

# getting indices of T-matrix
"""
    Get 2D array, first and second columns are n and m values, respectively.

m = -n : 1 : +n
n will be the repeated value of input
"""
function get_n_m_array_given_n(n)
    return hcat(repeat([n], 2 * n + 1), (-n):n)
end

"""
    Get 2D array, first and second columns are n and m values, respectively.

This is looping `get_n_m_array_given_n` over all n = 1:1:n_max
"""
function get_n_m_array_given_n_max(n_max)
    return vcat(get_n_m_array_given_n.(1:n_max)...)
end

"""
    Get 2D array, first, second and third columns are n, m, and idx values, respectively.

The idx value is just the index of the row.
"""
function get_n_m_idx_array_given_n_max(n_max)
    nm_array = get_n_m_array_given_n_max(n_max)
    return hcat(nm_array, 1:size(nm_array, 1))
end

"""
    Get 4 matrices for m, n, m_, n_, corresponding to rank and order of incident and scattered VSWF, represented by elements of T-matrix

## Arguments

    n_max : Int, determine the size of the T-matrix

## Returns

    m_matrix, n_matrix, m__matrix, n__matrix : each is a square matrix
"""
function get_m_n_m__n__matrices_for_T_matrix(n_max)
    n_m_idx = get_n_m_idx_array_given_n_max(n_max)
    n_array = n_m_idx[:, 1]
    m_array = n_m_idx[:, 2]
    idx_array = n_m_idx[:, 3]

    idx_max = idx_array[end]
    idx_matrix = repeat(idx_array, 1, idx_max)
    idx__matrix = repeat(idx_array', idx_max, 1)

    n_matrix = n_array[idx_matrix]
    m_matrix = m_array[idx_matrix]

    n__matrix = n_array[idx__matrix]
    m__matrix = m_array[idx__matrix]

    return m_matrix, n_matrix, m__matrix, n__matrix
end

"""
    using for loops, Get 4 matrices for m, n, m_, n_, corresponding to rank and order of incident and scattered VSWF, represented by elements of T-matrix
"""
function get_m_n_m__n__matrices_for_T_matrix_with_forloop(n_max)
    idx_max = Tmatrix.get_max_single_index_from_n_max(n_max)
    n_matrix = zeros(Int, idx_max, idx_max)
    m_matrix = zeros(Int, idx_max, idx_max)
    n__matrix = zeros(Int, idx_max, idx_max)
    m__matrix = zeros(Int, idx_max, idx_max)
    idx = 0
    for n in 1:n_max
        for m in (-n):n
            idx += 1
            idx_ = 0
            for n_ in 1:n_max
                for m_ in (-n_):n_
                    idx_ += 1
                    m_matrix[idx, idx_] = m
                    n_matrix[idx, idx_] = n
                    m__matrix[idx, idx_] = m_
                    n__matrix[idx, idx_] = n_
                end
            end
        end
    end
    return m_matrix, n_matrix, m__matrix, n__matrix
end

"""
    Create m,n,m_,n_ matrices with for loop, and compare the result with `get_m_n_m__n__matrices_for_T_matrix`
"""
function validate_get_m_n_m__n__matrices_for_T_matrix(n_max)
    m_matrix, n_matrix, m__matrix, n__matrix =
        get_m_n_m__n__matrices_for_T_matrix_with_forloop(n_max)
    _m_matrix, _n_matrix, _m__matrix, _n__matrix =
        get_m_n_m__n__matrices_for_T_matrix(n_max)

    println(_n_matrix == n_matrix)
    println(_m_matrix == m_matrix)
    println(_n__matrix == n__matrix)
    println(_m__matrix == m__matrix)
end

#######################
# TOBE moved to a new package for trapezoidal integraion, that will be compatible with autodiff
"""
    1D numerical integral using trapezoidal rule
x and y are 1D arrays
"""
#TODO rename to not elzouka lol
function trapz_ELZOUKA(
    x::AbstractVector{R},
    y::AbstractVector{N},
) where {R <: Real, N <: Number}
    # TODO: small error if compared with Trapz.trapz
    base = x[2:end] - x[1:(end - 1)]
    av_height = (y[2:end] + y[1:(end - 1)]) / 2
    areas = base .* av_height
    total_area = sum(areas)
    return total_area
end

"""
    2D numerical integral using trapezoidal rule

x and y are 1D arrays, z is 2D array
"""
function trapz_ELZOUKA(
    x::AbstractVector{R},
    y::AbstractVector{R},
    z::AbstractMatrix{N},
) where {R <: Real, N <: Number}
    integrand_wrt_x = trapz_ELZOUKA.(eachcol(repeat(x, 1, size(z, 2))), eachcol(z))
    return trapz_ELZOUKA(y, integrand_wrt_x)
end

#######################
# Calculate orientation-averaged sscattering, extinction and absorption cross sections.
"""
    calculate orientation averaged scattering cross section given a T-matrix and a wavevector.

Input T-matrix can be a complex square matrix or a concatenation of separate real and imag parts of T-matrix
"""
function get_OrentationAv_scattering_CrossSection_from_Tmatrix(
    T::AbstractMatrix,
    k1::C,
) where {C <: Complex}
    if size(T)[1] == size(T)[2] # if T-matrix is a square matrix, then T-matrix is complex
        return real(2 * pi / k1^2 * sum(T .* conj(T)))
    elseif size(T)[1] == size(T)[2] / 2 # if T-matrix has number of columns double the number of rows, then T-matrix is hcat() of real and imag parts of Tmatrix
        return get_OrentationAv_scattering_CrossSection_from_Tmatrix_SeparateRealImag(
            T,
            real(k1),
            imag(k1),
        )
    end
end

"""
    calculate orientation averaged scattering cross section given a T-matrix, wavelength and optical properties of surrounding.

Input T-matrix can be a complex square matrix or a concatenation of separate real and imag parts of T-matrix
"""
function get_OrentationAv_scattering_CrossSection_from_Tmatrix(
    T::AbstractMatrix,
    wl_or_freq_input::R,
    input_unit::String,
    Eps_r_1::C,
    Mu_r_1::C,
) where {C <: Complex, R <: Real}
    return get_OrentationAv_scattering_CrossSection_from_Tmatrix(
        T,
        get_WaveVector(
            wl_or_freq_input;
            input_unit = input_unit,
            Eps_r = Eps_r_1,
            Mu_r = Mu_r_1,
        ),
    )
end

"""
    calculate orientation averaged scattering cross section given a T-matrix and a wavevector. This works for "_SeparateRealImag"
"""
function get_OrentationAv_scattering_CrossSection_from_Tmatrix_SeparateRealImag(
    T::AbstractMatrix{R},
    k1_r::R,
    k1_i::R,
) where {R <: Real}
    T_r, T_i = Tmatrix.separate_real_imag(T)
    T_by_conjT_sum = sum(complex_multiply.(T_r, T_i, T_r, -T_i))
    k1_squared = complex_multiply(k1_r, k1_i, k1_r, k1_i)
    two_pi_over_k1_squared = complex_divide(2 * pi, 0, k1_squared[1], k1_squared[2])
    return complex_multiply(
        two_pi_over_k1_squared[1],
        two_pi_over_k1_squared[2],
        T_by_conjT_sum[1],
        T_by_conjT_sum[2],
    )[1]
end

"""
    calculate orientation averaged extinction cross section given a T-matrix and a wavevector.
"""
function get_OrentationAv_extinction_CrossSection_from_Tmatrix(
    T::AbstractMatrix,
    k1::C,
) where {C <: Complex}
    if size(T)[1] == size(T)[2] # if T-matrix is a square matrix, then T-matrix is complex
        return real(-2 * pi / k1^2 * tr(real(T)))
    elseif size(T)[1] == size(T)[2] / 2 # if T-matrix has number of columns double the number of rows, then T-matrix is hcat() of real and imag parts of Tmatrix
        return get_OrentationAv_extinction_CrossSection_from_Tmatrix_SeparateRealImag(
            T,
            real(k1),
            imag(k1),
        )
    end
end

"""
    calculate orientation averaged extinction cross section given a T-matrix and a wavevector. This works for "_SeparateRealImag"
"""
function get_OrentationAv_extinction_CrossSection_from_Tmatrix_SeparateRealImag(
    T::AbstractMatrix{R},
    k1_r::R,
    k1_i::R,
) where {R <: Real}
    T_r, T_i = Tmatrix.separate_real_imag(T)
    k1_squared = complex_multiply(k1_r, k1_i, k1_r, k1_i)
    two_pi_over_k1_squared = complex_divide(2 * pi, 0, k1_squared[1], k1_squared[2])
    return -two_pi_over_k1_squared[1] * tr(T_r)
end

"""
    calculate orientation averaged absorption cross section given a T-matrix and a wavevector.
"""
function get_OrentationAv_absorption_CrossSection_from_Tmatrix(
    T::AbstractMatrix,
    k1::C,
) where {C <: Complex}
    return get_OrentationAv_extinction_CrossSection_from_Tmatrix(T, k1) -
           get_OrentationAv_scattering_CrossSection_from_Tmatrix(T, k1)
end

"""
    calculate orientation averaged absorption cross section given a T-matrix and a wavevector. This works for "_SeparateRealImag"
"""
function get_OrentationAv_absorption_CrossSection_from_Tmatrix_SeparateRealImag(
    T::AbstractMatrix{R},
    k1_r::R,
    k1_i::R,
) where {R <: Real}
    return get_OrentationAv_extinction_CrossSection_from_Tmatrix_SeparateRealImag(
        T,
        k1_r,
        k1_i,
    ) - get_OrentationAv_scattering_CrossSection_from_Tmatrix_SeparateRealImag(
        T,
        k1_r,
        k1_i,
    )
end

"""
    calculate orientation averaged emissivity cross section given a T-matrix and a wavevector.

`particle_surface_area` has to be in units consistent with wavevector units. If waveevctor unit is per m, then the particle surface area has to be in m^2
"""
function get_OrentationAv_emissivity_from_Tmatrix(
    T::AbstractMatrix,
    k1::C,
    particle_surface_area::R,
) where {R <: Real, C <: Complex}
    return get_OrentationAv_absorption_CrossSection_from_Tmatrix(T, k1) * 4 /
           particle_surface_area
end

"""
    calculate orientation averaged emissivity cross section given a T-matrix and a wavevector. This works for "_SeparateRealImag"

`particle_surface_area` has to be in units consistent with wavevector units. If waveevctor unit is per m, then the particle surface area has to be in m^2
"""
function get_OrentationAv_emissivity_from_Tmatrix_SeparateRealImag(
    T::AbstractMatrix{R},
    k1_r::R,
    k1_i::R,
    particle_surface_area::R,
) where {R <: Real}
    return get_OrentationAv_absorption_CrossSection_from_Tmatrix_SeparateRealImag(
        T,
        k1_r,
        k1_i,
    ) * 4 / particle_surface_area
end

"""
    calculate orientation averaged emissivity cross section given a T-matrix, wavelength and optical properties of surrounding. This works for "_SeparateRealImag"

`particle_surface_area` has to be in units consistent with wavevector units. If waveevctor unit is per m, then the particle surface area has to be in m^2
"""
function get_OrentationAv_emissivity_from_Tmatrix_SeparateRealImag(
    T::AbstractMatrix{R},
    wl_or_freq_input::R,
    input_unit::String,
    Eps_r_r_1::R,
    Eps_r_i_1::R,
    Mu_r_r_1::R,
    Mu_r_i_1::R,
    particle_surface_area::R,
) where {R <: Real}
    k1_complex = get_WaveVector(
        wl_or_freq_input;
        input_unit = input_unit,
        Eps_r = Complex(Eps_r_r_1, Eps_r_i_1),
        Mu_r = Complex(Mu_r_r_1, Mu_r_i_1),
    )
    return get_OrentationAv_absorption_CrossSection_from_Tmatrix_SeparateRealImag(
        T,
        k1_r,
        k1_i,
    ) * 4 / particle_surface_area
end


text:
module TestStatifier

using Test
using TSML

function test_statifier()

  Random.seed!(123)
  dt=[missing;rand(1:10,3);missing;missing;missing;rand(1:5,3)]
  dat = DataFrame(Date= DateTime(2017,12,31,1):Dates.Hour(1):DateTime(2017,12,31,10) |> collect,
		  Value = dt)
  statfier = Statifier(Dict(:processmissing=>false))
  fit!(statfier,dat)
  res=transform!(statfier,dat)
  @test ncol(res) == 20
  @test res[1,3:end]|>Vector|>sum |> x->round(x,digits=5) == 16.29114

  statfier = Statifier(Dict(:processmissing=>true))
  fit!(statfier,dat)
  res=transform!(statfier,dat)
  @test ncol(res) == 26
  @test res[1,3:end]|>Vector|>sum |> x->round(x,digits=5) == 28.29114

  fname = joinpath(dirname(pathof(TSML)),"../data/testdata.csv")
  csvfilter = CSVDateValReader(Dict(:filename=>fname,:dateformat=>"dd/mm/yyyy HH:MM"))
  valgator = DateValgator(Dict(:dateinterval=>Dates.Hour(1)))
  valnner = DateValNNer(Dict(:dateinterval=>Dates.Hour(1)))

  stfier = Statifier(Dict(:processmissing=>true))
  mpipeline0 = @pipeline csvfilter |> stfier
  respipe0 = fit_transform!(mpipeline0)
  @test round(respipe0[1,:sfreq],digits=2) ==  0.18
  vals = respipe0[1,3:end] |> Vector
  @test (vals[(!isnan).(vals)] |> sum |> x->round(x,sigdigits=2)) == -1.3e6

  stfier = Statifier(Dict(:processmissing=>true))
  mpipeline1 = @pipeline csvfilter |> valgator |> stfier
  fit!(mpipeline1)
  respipe1 = transform!(mpipeline1)
  @test respipe1[1,3:end] |> Vector |> sum |> x->round(x,digits=5) == -102779.88226

  mpipeline2 = @pipeline csvfilter |> valgator |> valnner |> stfier
  fit!(mpipeline2)
  respipe2 = transform!(mpipeline2)
  res2 = respipe2[1,3:end] |> Vector
  @test (res2[(!isnan).(res2)] |> sum |> x->round(x,sigdigits=2)) == -230000.00

  stfier = Statifier(Dict(:processmissing=>false))
  mpipeline1 = @pipeline csvfilter |> valgator |> stfier
  respipe1 = fit_transform!(mpipeline1)
  res1 = respipe1[1,3:end] |> Vector
  @test res1 |> sum |> x->round(x,digits=5) == -105185.44114

  stfier = Statifier(Dict(:processmissing=>true))
  mpipeline2 = @pipeline csvfilter |> valgator |> valnner |> stfier
  respipe2 = fit_transform!(mpipeline2)
  res2 = respipe2[1,3:end] |> Vector
  @test res2[(!isnan).(res2)] |> sum |> x->round(x,digits=5) == -227824.85354

end
@testset "Statifier: readcsv |> valgator |> valnner |> stfier" begin
  test_statifier()
end

end


text:
function myzip(xs, ys)
    n = min(length(xs),length(ys))
    return map(i->(xs[i],ys[i]), [i for i in 1:n])
end
function myunzip(xys)
    xs = map(pair -> pair[1], xys)
    ys = map(pair -> pair[2], xys)
    return xs,ys
end

max_stars_repo_path: II rok/I semestr/ANM/Pracownia 3/prog/funkcje.jl

max_stars_repo_name: antonitomaszewski/Studia

max_issues_repo_name: antonitomaszewski/Studia

max_forks_repo_name: antonitomaszewski/Studia


text:
########################################################################################################################
module GenerateTimeSeries
########################################################################################################################
# imports and exports
using Debugger, JLD2, Distributed
using Random, Distributions, LinearAlgebra
using ..DeSolvers, ..L96, ..IEEE39bus
export L96_time_series, IEEE39bus_time_series

########################################################################################################################
# generate timeseries based on the model, solver and parameters

function L96_time_series(args::Tuple{Int64,Int64,Float64,Int64,Int64,Float64,Float64})

    # time the experiment
    t1 = time()

    # unpack the experiment parameters determining the time series
    seed, state_dim, tanl, nanl, spin, diffusion, F = args

    # define the model
    dx_dt = L96.dx_dt
    dx_params = Dict{String, Array{Float64}}("F" => [8.0])

    # define the integration scheme
    if diffusion == 0.0
        # generate the observations with the Runge-Kutta scheme
        step_model! = DeSolvers.rk4_step!

        # parameters for the Runge-Kutta scheme
        h = 0.05

    else
        # generate the observations with the strong Taylor scheme
        step_model! = L96.l96s_tay2_step!

        # parameters for the order 2.0 strong Taylor scheme
        h = 0.005
        p = 1
        α = L96.α(p)
        ρ = L96.ρ(p)
    end

    # set the number of discrete integrations steps between each observation time
    f_steps = convert(Int64, tanl/h)

    # set storage for the ensemble timeseries
    obs = Array{Float64}(undef, state_dim, nanl)

    # define the integration parameters in the kwargs dict
    kwargs = Dict{String, Any}(
                               "h" => h,
                               "diffusion" => diffusion,
                               "dx_params" => dx_params,
                               "dx_dt" => dx_dt,
                              )
    if diffusion != 0.0
        kwargs["p"] = p
        kwargs["α"] = α
        kwargs["ρ"] = ρ
    end

    # seed the random generator
    Random.seed!(seed)
    x = rand(Normal(), state_dim)

    # spin the model onto the attractor
    for j in 1:spin
        for k in 1:f_steps
            step_model!(x, 0.0, kwargs)
        end
    end

    # save the model state at timesteps of tanl
    for j in 1:nanl
        for k in 1:f_steps
            step_model!(x, 0.0, kwargs)
        end
        obs[:, j] = x
    end

    data = Dict{String, Any}(
                             "h" => h,
                             "diffusion" => diffusion,
                             "dx_params" => dx_params,
                             "tanl" => tanl,
                             "nanl" => nanl,
                             "spin" => spin,
                             "state_dim" => state_dim,
                             "obs" => obs,
                             "model" => "L96"
                            )

    name = "L96_time_series_seed_" * lpad(seed, 4, "0") *
           "_dim_" * lpad(state_dim, 2, "0") *
           "_diff_" * rpad(diffusion, 5, "0") *
           "_F_" * lpad(F, 4, "0") *
           "_tanl_" * rpad(tanl, 4, "0") *
           "_nanl_" * lpad(nanl, 5, "0") *
           "_spin_" * lpad(spin, 4, "0") *
           "_h_" * rpad(h, 5, "0") *
           ".jld2"

    path = joinpath(@__DIR__, "../data/time_series/")
    save(path * name, data)
    print("Runtime " * string(round((time() - t1)  / 60.0, digits=4))  * " minutes\n")
end

########################################################################################################################
# generate timeseries based on the parameter tuple

function IEEE39bus_time_series(args::Tuple{Int64,Float64,Int64,Int64,Float64})

    # time the experiment
    t1 = time()

    # unpack the experiment parameters determining the time series
    seed, tanl, nanl, spin, diffusion = args
    Random.seed!(seed)

    # define the model
    dx_dt = IEEE39bus.dx_dt
    state_dim = 20

    # define the model parameters
    input_data = joinpath(@__DIR__, "../models/IEEE39bus_inputs/NE_EffectiveNetworkParams.jld2")
    tmp = load(input_data)
    dx_params = Dict{String, Array{Float64}}(
                                             "A" => tmp["A"],
                                             "D" => tmp["D"],
                                             "H" => tmp["H"],
                                             "K" => tmp["K"],
                                             "γ" => tmp["γ"],
                                             "ω" => tmp["ω"]
                                            )

    # define the integration scheme
    step_model! = DeSolvers.rk4_step!
    h = 0.01

    # define the diffusion coefficient structure matrix
    # notice that the perturbations are only applied to the frequencies
    # based on the change of variables derivation
    # likewise, the diffusion parameter is applied separately as an amplitude
    # in the Runge-Kutta scheme
    diff_mat = zeros(20,20)
    diff_mat[LinearAlgebra.diagind(diff_mat)[11:end]] = tmp["ω"][1] ./ (2.0 * tmp["H"])

    # set the number of discrete integrations steps between each observation time
    f_steps = convert(Int64, tanl/h)

    # set storage for the ensemble timeseries
    obs = Array{Float64}(undef, state_dim, nanl)

    # define the integration parameters in the kwargs dict
    kwargs = Dict{String, Any}(
              "h" => h,
              "diffusion" => diffusion,
              "dx_params" => dx_params,
              "dx_dt" => dx_dt,
              "diff_mat" => diff_mat
             )

    # load the steady state, generated by long simulation without noise
    x = tmp["synchronous_state"]

    # spin the model onto the attractor
    for j in 1:spin
        for k in 1:f_steps
            step_model!(x, 0.0, kwargs)
            # set phase angles mod 2pi
            x[1:10] .= rem2pi.(x[1:10], RoundNearest)
        end
    end

    # save the model state at timesteps of tanl
    for j in 1:nanl
        for k in 1:f_steps
            step_model!(x, 0.0, kwargs)
            # set phase angles mod 2pi
            x[1:10] .= rem2pi.(x[1:10], RoundNearest)
        end
        obs[:, j] = x
    end

    data = Dict{String, Any}(
                             "h" => h,
                             "diffusion" => diffusion,
                             "diff_mat" => diff_mat,
                             "dx_params" => dx_params,
                             "tanl" => tanl,
                             "nanl" => nanl,
                             "spin" => spin,
                             "obs" => obs,
                             "model" => "IEEE39bus"
                            )

    name = "IEEE39bus_time_series_seed_" * lpad(seed, 4, "0") *
           "_diff_" * rpad(diffusion, 5, "0") *
           "_tanl_" * rpad(tanl, 4, "0") *
           "_nanl_" * lpad(nanl, 5, "0") *
           "_spin_" * lpad(spin, 4, "0") *
           "_h_" * rpad(h, 5, "0") *
           ".jld2"

    path = joinpath(@__DIR__, "../data/time_series/")
    save(path * name, data)
    print("Runtime " * string(round((time() - t1)  / 60.0, digits=4))  * " minutes\n")
end

########################################################################################################################
# end module

end

max_stars_repo_name: sukhreens/DataAssimilationBenchmarks

max_issues_repo_name: sukhreens/DataAssimilationBenchmarks

max_forks_repo_name: sukhreens/DataAssimilationBenchmarks


text:
using VCFTools

chr = 20
truth  = "target.chr$chr.full.vcf.gz"
mendel = "mendel.imputed.chr$chr.vcf.gz"
beagle = "beagle.imputed.chr$chr.vcf.gz"
mmac4  = "minimac4.chr$chr.result.dose.vcf.gz"

isfile(truth) && isfile(mendel) && isfile(beagle) && isfile(mmac4) || error("file doens't exist!")

Xtrue = convert_gt(UInt8, truth, msg="importing Xtrue")
Xmendel = convert_gt(UInt8, mendel, msg="importing Xmendel")
Xbeagle = convert_gt(UInt8, beagle, msg="importing Xbeagle")
Xmmac4 = convert_gt(UInt8, mmac4, msg="importing Xmmac")
n, p = size(Xtrue)

# need to import each chunk of impute5 result separately
chunks = chr == 20 ? 3 : 7
Ximp5 = Vector{Matrix{UInt8}}(undef, chunks)
for i in 1:chunks
    Ximp5[i] = convert_gt(UInt8, "impute5.chr$chr.result.chunk$i.vcf.gz", msg="importing Ximp5 $i / $chunks")
end
Ximp5 = hcat(Ximp5...)

println("error mendel = $(sum(Xmendel .!= Xtrue) / n / p)")
println("error beagle = $(sum(Xbeagle .!= Xtrue) / n / p)")
println("error minimac4 = $(sum(Xmmac4 .!= Xtrue) / n / p)")
println("error impute5 = $(sum(Ximp5 .!= Xtrue) / n / p)")

max_stars_repo_name: Hua-Zhou/MendelImpute


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

args = ARGS

if length(args) < 1
    println(stdout, "Specify pkg name.")
    exit(1)
end

Pkg.add(args[1])

max_stars_repo_name: mkanenobu/myToolBox

max_issues_repo_name: mkanenobu/myToolBox

max_forks_repo_name: mkanenobu/toolbox


text:
src_path = joinpath(splitdir(dirname(@__FILE__))[1], "src")
push!(LOAD_PATH, src_path)

using Documenter, Slomo

makedocs(
    sitename = "Slomo.jl",
    modules = [Slomo],
    doctest = false,
    pages = [
	"Home" => "index.md",
	"Installation" => "install.md",
        "Getting started" => "start.md",
        "Examples" => "examples.md",
        "Public API" => "public_api.md",
        "Private API" => "private_api.md"
    ]
)

deploydocs(
    repo = "github.com/adwasser/Slomo.jl.git"
)


text:
using Distributions
using LinearAlgebra
using Random
using Test

using EnsembleKalmanProcesses
using EnsembleKalmanProcesses.ParameterDistributions
using EnsembleKalmanProcesses.Localizers
import EnsembleKalmanProcesses: construct_mean, construct_cov, construct_sigma_ensemble
const EKP = EnsembleKalmanProcesses

@testset "Localization" begin

    # Seed for pseudo-random number generator
    rng_seed = 42
    rng = Random.MersenneTwister(rng_seed)

    # Linear problem with d == p >> N_ens - Section 6.1 of Tong and Morzfeld (2022)
    G(u) = u
    N_ens = 10
    p = 50
    N_iter = 20
    # Generate random truth
    y = 10.0 * rand(p)
    Γ = 1.0 * I

    #### Define prior information on parameters
    priors = ParameterDistribution[] #empty PD-array
    for i in 1:p
        push!(priors, ParameterDistribution(Parameterized(Normal(0.0, 0.5)), no_constraint(), string("u", i)))
    end
    prior = combine_distributions(priors)
    #prior_distns = repeat([Parameterized(Normal(0.0, 0.5))], p)
    #constraints = repeat([[no_constraint()]], p)
    #prior_names = [string("u_", i) for i in 1:p]
    #prior = ParameterDistribution(prior_distns, constraints, prior_names)

    initial_ensemble = EKP.construct_initial_ensemble(rng, prior, N_ens)

    # Solve problem without localization
    ekiobj_vanilla = EKP.EnsembleKalmanProcess(initial_ensemble, y, Γ, Inversion(); rng = rng)
    for i in 1:N_iter
        g_ens_vanilla = G(get_u_final(ekiobj_vanilla))
        EKP.update_ensemble!(ekiobj_vanilla, g_ens_vanilla)
    end
    nonlocalized_error = get_error(ekiobj_vanilla)[end]

    # Test different localizers
    loc_methods = [Delta(), RBF(1.0), RBF(0.1), BernoulliDropout(0.1)]

    for loc_method in loc_methods
        ekiobj =
            EKP.EnsembleKalmanProcess(initial_ensemble, y, Γ, Inversion(); rng = rng, localization_method = loc_method)
        @test isa(ekiobj.localizer, Localizer)

        for i in 1:N_iter
            g_ens = G(get_u_final(ekiobj))
            EKP.update_ensemble!(ekiobj, g_ens, deterministic_forward_map = true)
        end

        # Test that localized version does better in the setting p >> N_ens
        @test get_error(ekiobj)[end] < nonlocalized_error

        # Test Schur product theorem
        u_final = get_u_final(ekiobj)
        g_final = get_g_final(ekiobj)
        cov_ug = cov(u_final, g_final, dims = 2, corrected = false)
        kernel = ekiobj.localizer.kernel
        @test rank(cov_ug) < rank(kernel .* cov_ug)
    end

end