text:
@muladd function perform_step!(integrator,cache::DRI1ConstantCache,f=integrator.f)
  @unpack a021,a031,a032,a121,a131,b021,b031,b121,b131,b221,b222,b223,b231,b232,b233,α1,α2,α3,c02,c03,c12,c13,beta11,beta12,beta13,beta22,beta23,beta31,beta32,beta33,beta42,beta43,NORMAL_ONESIX_QUANTILE = cache
  @unpack t,dt,uprev,u,W,p = integrator

  # define three-point distributed random variables
  dW_scaled = W.dW / sqrt(dt)
  _dW = map(x -> calc_threepoint_random(integrator, NORMAL_ONESIX_QUANTILE, x), dW_scaled)
  chi1 = map(x -> (x^2-dt)/2, _dW) # diagonal of Ihat2
  if !(typeof(W.dW) <: Number)
    m = length(W.dW)
    # define two-point distributed random variables
    _dZ = map(x -> calc_twopoint_random(integrator, x),  W.dZ)
    Ihat2 = zeros(eltype(W.dZ), m, m) # I^_(k,l)
    for k = 1:m
      for l = 1:m
        if k<l
          Ihat2[k, l] = (_dW[k]*_dW[l]-integrator.sqdt*_dZ[k])/2
        elseif l<k
          Ihat2[k, l] = (_dW[k]*_dW[l]+integrator.sqdt*_dZ[l])/2
        else k==l
          Ihat2[k, k] = chi1[k]
        end
      end
    end
  end

  # compute stage values
  k1 = integrator.f(uprev,p,t)
  g1 = integrator.g(uprev,p,t)

  # H_i^(0), stage 1
  # H01 = uprev
  # H_i^(0), stage 2
  if !is_diagonal_noise(integrator.sol.prob) || typeof(W.dW) <: Number
    H02 = uprev + a021*k1*dt + b021*g1*_dW
  else
    H02 = uprev + a021*k1*dt + b021*g1.*_dW
  end

  # H_i^(0), stage 3 (requires H_i^(k) stage 2 in general)
  k2 = integrator.f(H02,p,t+c02*dt)
  H03 = uprev + a032*k2*dt + a031*k1*dt
  if !is_diagonal_noise(integrator.sol.prob) || typeof(W.dW) <: Number
    H03 += b031*g1*_dW
  else
    H03 += b031*g1.*_dW
  end

  k3 = integrator.f(H03,p,t+c03*dt)

  # H_i^(k), stage 1
  # H11 = uprev
  # H_i^(k), stage 2
  if typeof(W.dW) <: Number
    H12 = uprev + a121*k1*dt + b121*g1*integrator.sqdt
  elseif is_diagonal_noise(integrator.sol.prob)
    H12 = Vector{typeof(uprev)}[uprev .+ a121*k1*dt .+ b121*g1[k]*integrator.sqdt for k=1:m]
  else
    H12 = Vector{typeof(uprev)}[uprev .+ a121*k1*dt .+ b121*g1[:,k]*integrator.sqdt for k=1:m]
  end

  # # H_i^(k), stage 3
  if typeof(W.dW) <: Number
    H13 = uprev + a131*k1*dt + b131*g1*integrator.sqdt
  elseif is_diagonal_noise(integrator.sol.prob)
    H13 = Vector{typeof(uprev)}[uprev .+ a131*k1*dt .+ b131*g1[k]*integrator.sqdt for k=1:m]
  else
    H13 = Vector{typeof(uprev)}[uprev .+ a131*k1*dt .+ b131*g1[:,k]*integrator.sqdt for k=1:m]
  end

  # H^_i^(k), stage 1
  # H21 = uprev
  # H^_i^(k), stage 2

  if typeof(W.dW) <: Number
    g2 = integrator.g(H12,p,t+c12*dt)
    g3 = integrator.g(H13,p,t+c13*dt)
    # for m=1:  H22 = uprev
  else
    g2 = [integrator.g(H12[k],p,t+c12*dt) for k=1:m]
    g3 = [integrator.g(H13[k],p,t+c13*dt) for k=1:m]
    H22 = [uprev for k=1:m]
    # add inbounds for speed if working properly
    for k=1:m
      for l=1:m
        if k == l
          continue
        end
        if is_diagonal_noise(integrator.sol.prob)
          @.. H22[k] += (b221*g1[l]+b222*g2[l][l]+b223*g3[l][l])*Ihat2[k,l]/integrator.sqdt
        else
          H22[k] += (b221*g1[:,l]+b222*g2[l][:,l]+b223*g3[l][:,l])*Ihat2[k,l]/integrator.sqdt
        end
      end
    end
  end

  # H^_i^(k), stage 3
  # for m=1: H23 = uprev
  if !(typeof(W.dW) <: Number)
    H23 = [uprev for k=1:m]
    # add inbounds for speed if working properly
    for k=1:m
      for l=1:m
        if k == l
          continue
        end
        if is_diagonal_noise(integrator.sol.prob)
          @.. H23[k] += (b231*g1[l]+b232*g2[l][l]+b233*g3[l][l])*Ihat2[k,l]/integrator.sqdt
        else
          H23[k] += (b231*g1[:,l]+b232*g2[l][:,l]+b233*g3[l][:,l])*Ihat2[k,l]/integrator.sqdt
        end
      end
    end
  end

  # add stages together Eq. (3)
  u = uprev + α1*k1*dt + α2*k2*dt + α3*k3*dt

  # add noise
  if typeof(W.dW) <: Number
    # lines 2 and 3
    u += g1*(_dW*beta11) + g2*(_dW*beta12+chi1*beta22/integrator.sqdt) + g3*(_dW*beta13+chi1*beta23/integrator.sqdt)
    # lines 4 and 5 are zero by construction
    # u += g1*(_dW*(beta31+beta32+beta33)+chi1*integrator.sqdt*(beta42+beta43))
  else
      if is_diagonal_noise(integrator.sol.prob)
        u += g1.*_dW*beta11
        for k=1:m
          u[k] += g2[k][k]*_dW[k]*beta12+g3[k][k]*_dW[k]*beta13
          u[k] += g2[k][k]*chi1[k]*beta22/integrator.sqdt + g3[k][k]*chi1[k]*beta23/integrator.sqdt
          tmpg = integrator.g(H22[k],p,t)
          u[k] = u[k] + g1[k]*_dW[k]*beta31 + tmpg[k]*(_dW[k]*beta32 + integrator.sqdt*beta42)
          tmpg = integrator.g(H23[k],p,t)
          u[k] = u[k] + tmpg[k]*(_dW[k]*beta33 + integrator.sqdt*beta43)
        end
      else
        # non-diag noise
        for k=1:m
          g1k = @view g1[:,k]
          g2k = @view g2[k][:,k]
          g3k = @view g3[k][:,k]
          @.. u = u + g1k*_dW[k]*(beta11+beta31)
          @.. u = u + g2k*_dW[k]*beta12 + g3k*_dW[k]*beta13
          @.. u = u + g2k*chi1[k]*beta22/integrator.sqdt + g3k*chi1[k]*beta23/integrator.sqdt
          tmpg = integrator.g(H22[k],p,t)
          @.. u = u + tmpgk*_dW[k]*beta32 + tmpgk*integrator.sqdt*beta42
          tmpg = integrator.g(H23[k],p,t)
          @.. u = u + tmpgk*_dW[k]*beta33 + tmpgk*integrator.sqdt*beta43
        end
      end
  end

  # Test EM scheme:
  # integrator.u = uprev + 1*k1*dt + 1*g1*_dW

  integrator.u = u

end

@muladd function perform_step!(integrator,cache::DRI1Cache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack _dW,_dZ,chi1,Ihat2,tab,g1,g2,g3,k1,k2,k3,H02,H03,H12,H13,H22,H23,tmp1,tmpg = cache
  @unpack a021,a031,a032,a121,a131,b021,b031,b121,b131,b221,b222,b223,b231,b232,b233,α1,α2,α3,c02,c03,c12,c13,beta11,beta12,beta13,beta22,beta23,beta31,beta32,beta33,beta42,beta43,NORMAL_ONESIX_QUANTILE = cache.tab

  m = length(W.dW)

  if typeof(W.dW) <: Union{SArray,Number}
    # tbd
    _dW = map(x -> calc_threepoint_random(integrator, NORMAL_ONESIX_QUANTILE, x), W.dW / sqrt(dt))
  else
    # define three-point distributed random variables
    @.. chi1 = W.dW / sqrt(dt)
    calc_threepoint_random!(_dW, integrator, NORMAL_ONESIX_QUANTILE, chi1)
    map!(x -> (x^2-dt)/2, chi1, _dW)

    # define two-point distributed random variables
    calc_twopoint_random!(_dZ,integrator,W.dZ)

    for k = 1:m
      for l = 1:m
        if k<l
          Ihat2[k, l] = (_dW[k]*_dW[l]-integrator.sqdt*_dZ[k])/2
        elseif l<k
          Ihat2[k, l] = (_dW[k]*_dW[l]+integrator.sqdt*_dZ[l])/2
        else k==l
          Ihat2[k, k] = chi1[k]
        end
      end
    end
  end

  # compute stage values
  integrator.f(k1,uprev,p,t)
  integrator.g(g1,uprev,p,t)

  # H_i^(0), stage 1
  # H01 = uprev
  # H_i^(0), stage 2
  if is_diagonal_noise(integrator.sol.prob)
    @.. H02 = uprev + a021*k1*dt + b021*g1*_dW
  else
    mul!(tmp1,g1,_dW)
    @.. H02 = uprev + dt*a021*k1 + b021*tmp1
  end

  integrator.f(k2,H02,p,t+c02*dt)
  if is_diagonal_noise(integrator.sol.prob)
    @.. H03 = uprev + a032*k2*dt + a031*k1*dt + b031*g1*_dW
  else
    @.. H03 = uprev + a032*k2*dt + a031*k1*dt + b031*tmp1
  end

  integrator.f(k3,H03,p,t+c03*dt)

  # H_i^(k), stages
  # H11 = uprev
  for k=1:m
    if is_diagonal_noise(integrator.sol.prob)
      @.. H12[k] = uprev + a121*k1*dt + b121*g1[k]*integrator.sqdt
      @.. H13[k] = uprev + a131*k1*dt + b131*g1[k]*integrator.sqdt
    else
      g1k = @view g1[:,k]
      @.. H12[k] = uprev + a121*k1*dt + b121*g1k*integrator.sqdt
      @.. H13[k] = uprev + a131*k1*dt + b131*g1k*integrator.sqdt
    end
  end
  # H^_i^(k), stages

  for k=1:m
    integrator.g(g2[k],H12[k],p,t+c12*dt)
    integrator.g(g3[k],H13[k],p,t+c13*dt)

    H22[k] = uprev
    H23[k] = uprev
  end

  for k=1:m
    for l=1:m
      if k == l
        continue
      end
      if is_diagonal_noise(integrator.sol.prob)
        @.. H22[k] = H22[k]+(b221*g1[l]+b222*g2[l][l]+b223*g3[l][l])*Ihat2[k,l]/integrator.sqdt
        @.. H23[k] = H23[k]+(b231*g1[l]+b232*g2[l][l]+b233*g3[l][l])*Ihat2[k,l]/integrator.sqdt
      else
        g1l = @view g1[:,l]
        g2l = @view g2[l][:,l]
        g3l = @view g3[l][:,l]
        @.. H22[k] = H22[k]+(b221*g1l+b222*g2l+b223*g3l)*Ihat2[k,l]/integrator.sqdt
        @.. H23[k] = H23[k]+(b231*g1l+b232*g2l+b233*g3l)*Ihat2[k,l]/integrator.sqdt
      end
    end
  end

  # add stages together Eq. (3)
  @.. u = uprev + α1*k1*dt + α2*k2*dt + α3*k3*dt

  # add noise
  if is_diagonal_noise(integrator.sol.prob)
      @.. u = u + g1*_dW*beta11
      for k=1:m
        u[k] = u[k] + g2[k][k]*_dW[k]*beta12 + g3[k][k]*_dW[k]*beta13
        u[k] = u[k]  + g2[k][k]*chi1[k]*beta22/integrator.sqdt + g3[k][k]*chi1[k]*beta23/integrator.sqdt
        if m > 1
          integrator.g(tmpg,H22[k],p,t)
          u[k] = u[k] + g1[k]*_dW[k]*beta31 + tmpg[k]*(_dW[k]*beta32 + integrator.sqdt*beta42)
          integrator.g(tmpg,H23[k],p,t)
          u[k] = u[k] + tmpg[k]*(_dW[k]*beta33 + integrator.sqdt*beta43)
        end
      end
  else
    for k=1:m
      g1k = @view g1[:,k]
      g2k = @view g2[k][:,k]
      g3k = @view g3[k][:,k]
      tmpgk = @view tmpg[:,k]
      @.. u = u + g1k*_dW[k]*(beta11+beta31)
      @.. u = u + g2k*_dW[k]*beta12 + g3k*_dW[k]*beta13
      @.. u = u + g2k*chi1[k]*beta22/integrator.sqdt + g3k*chi1[k]*beta23/integrator.sqdt
      integrator.g(tmpg,H22[k],p,t)
      @.. u = u + tmpgk*_dW[k]*beta32 + tmpgk*integrator.sqdt*beta42
      integrator.g(tmpg,H23[k],p,t)
      @.. u = u + tmpgk*_dW[k]*beta33 + tmpgk*integrator.sqdt*beta43
    end
  end

end

# Roessler SRK for first order weak approx
@muladd function perform_step!(integrator,cache::RDI1WMConstantCache,f=integrator.f)
  @unpack a021,b021,α1,α2,c02,beta11,NORMAL_ONESIX_QUANTILE = cache
  @unpack t,dt,uprev,u,W,p = integrator

  # define three-point distributed random variables
  dW_scaled = W.dW / sqrt(dt)
  _dW = map(x -> calc_threepoint_random(integrator, NORMAL_ONESIX_QUANTILE, x), dW_scaled)
  chi1 = map(x -> (x^2-dt)/2, _dW) # diagonal of Ihat2
  if !(typeof(W.dW) <: Number)
    m = length(W.dW)
    # define two-point distributed random variables
    _dZ = map(x -> calc_twopoint_random(integrator, x),  W.dZ)
    Ihat2 = zeros(eltype(W.dZ), m, m) # I^_(k,l)
    for k = 1:m
      for l = 1:m
        if k<l
          Ihat2[k, l] = (_dW[k]*_dW[l]-integrator.sqdt*_dZ[k])/2
        elseif l<k
          Ihat2[k, l] = (_dW[k]*_dW[l]+integrator.sqdt*_dZ[l])/2
        else k==l
          Ihat2[k, k] = chi1[k]
        end
      end
    end
  end

  # compute stage values
  k1 = integrator.f(uprev,p,t)
  g1 = integrator.g(uprev,p,t)

  # H_i^(0), stage 1
  # H01 = uprev
  # H_i^(0), stage 2
  if !is_diagonal_noise(integrator.sol.prob) || typeof(W.dW) <: Number
    H02 = uprev + a021*k1*dt + b021*g1*_dW
  else
    H02 = uprev + a021*k1*dt + b021*g1.*_dW
  end
  k2 = integrator.f(H02,p,t+c02*dt)

  # add stages together Eq. (3)
  u = uprev + α1*k1*dt + α2*k2*dt

  # add noise
  if typeof(W.dW) <: Number
    # lines 2 and 3
    u += g1*(_dW*beta11)
  else
      if is_diagonal_noise(integrator.sol.prob)
        u += g1.*_dW*beta11
      else
        # non-diag noise
        for k=1:m
          g1k = @view g1[:,k]
          @.. u = u + g1k*_dW[k]*(beta11)
        end
      end
  end

  integrator.u = u

end

@muladd function perform_step!(integrator,cache::RDI1WMCache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack _dW,_dZ,chi1,Ihat2,tab,g1,k1,k2,H02,tmp1 = cache
  @unpack a021,b021,α1,α2,c02,beta11,NORMAL_ONESIX_QUANTILE = cache.tab

  m = length(W.dW)

  if typeof(W.dW) <: Union{SArray,Number}
    # tbd
    _dW = map(x -> calc_threepoint_random(integrator, NORMAL_ONESIX_QUANTILE, x), W.dW / sqrt(dt))
  else
    # define three-point distributed random variables
    @.. chi1 = W.dW / sqrt(dt)
    calc_threepoint_random!(_dW, integrator, NORMAL_ONESIX_QUANTILE, chi1)
    map!(x -> (x^2-dt)/2, chi1, _dW)

    # define two-point distributed random variables
    calc_twopoint_random!(_dZ,integrator,W.dZ)

    for k = 1:m
      for l = 1:m
        if k<l
          Ihat2[k, l] = (_dW[k]*_dW[l]-integrator.sqdt*_dZ[k])/2
        elseif l<k
          Ihat2[k, l] = (_dW[k]*_dW[l]+integrator.sqdt*_dZ[l])/2
        else k==l
          Ihat2[k, k] = chi1[k]
        end
      end
    end
  end

  # compute stage values
  integrator.f(k1,uprev,p,t)
  integrator.g(g1,uprev,p,t)

  # H_i^(0), stage 1
  # H01 = uprev
  # H_i^(0), stage 2
  if is_diagonal_noise(integrator.sol.prob)
    @.. H02 = uprev + a021*k1*dt + b021*g1*_dW
  else
    mul!(tmp1,g1,_dW)
    @.. H02 = uprev + dt*a021*k1 + b021*tmp1
  end

  integrator.f(k2,H02,p,t+c02*dt)

  # add stages together Eq. (3)
  @.. u = uprev + α1*k1*dt + α2*k2*dt

  # add noise
  if is_diagonal_noise(integrator.sol.prob)
      @.. u = u + g1*_dW*beta11
  else
    for k=1:m
      g1k = @view g1[:,k]
      @.. u = u + g1k*_dW[k]*beta11
    end
  end
end

# Stratonovich sense

@muladd function perform_step!(integrator,cache::RSConstantCache,f=integrator.f)
  @unpack a021,a031,a032,a131,a141,b031,b032,b121,b131,b132,b141,b142,b143,b221,b231,b331,b332,b341,b342,α1,α2,α3,α4,c02,c03,c13,c14,beta11,beta12,beta13,beta14,beta22,beta23,NORMAL_ONESIX_QUANTILE = cache
  @unpack t,dt,uprev,u,W,p = integrator

  # define three-point distributed random variables
  dW_scaled = W.dW / sqrt(dt)
  _dW = map(x -> calc_threepoint_random(integrator, NORMAL_ONESIX_QUANTILE, x), dW_scaled)
  if !(typeof(W.dW) <: Number)
    m = length(W.dW)
    # define two-point distributed random variables
    _dZ = map(x -> calc_twopoint_random(integrator, x),  W.dZ)
    Ihat2 = zeros(eltype(W.dZ), m, m) # I^_(k,l)
    for k = 1:m
      for l = 1:k-1
        Ihat2[k, l] = _dW[k]*_dZ[l]
        Ihat2[l, k] = -_dW[k]*_dZ[l]
      end
    end
  end
  # compute stage values
  k1 = integrator.f(uprev,p,t)
  g1 = integrator.g(uprev,p,t)

  # H_1^(0)
  # H01 = uprev
  # H_2^(0)
  H02 = uprev + a021*k1*dt

  # H_1^(k)
  # H11 = uprev
  # H_2^(k), stage 2
  if typeof(W.dW) <: Number
    H12 = uprev + b121*g1*_dW
  elseif is_diagonal_noise(integrator.sol.prob)
    H12 = Vector{typeof(uprev)}[uprev .+ b121*g1[k]*_dW[k] for k=1:m]
  else
    H12 = Vector{typeof(uprev)}[uprev .+ b121*g1[:,k]*_dW[k] for k=1:m]
  end

  if typeof(W.dW) <: Number
    g2 = integrator.g(H12,p,t)
  else
    g2 = [integrator.g(H12[k],p,t) for k=1:m]
  end

  # H_3^(k)
  if typeof(W.dW) <: Number
    H13 = uprev + a131*k1*dt + b131*g1*_dW + b132*g2*_dW
  else
    H13 = [zero(typeof(uprev)) for k=1:m]
    for k=1:m
      H13[k] += uprev .+ a131*k1*dt
      if is_diagonal_noise(integrator.sol.prob)
        H13[k] += b131*g1[k]*_dW[k] .+ b132*g2[k]*_dW[k]
        for l=1:m
          if l!=k
            H13[k] += b331*g1[l]*_dW[l] .+ b332*g2[l][l]*_dW[l]
          end
        end
      else
        H13 += b131*g1[:,k]*_dW[k] .+ b132*g2[k][:,k]*_dW[k]
        for l=1:m
          if l!=k
            H13[k] += b331*g1[:,l]*_dW[l] .+ b332*g2[l][:,l]*_dW[l]
          end
        end
      end
    end
  end

  if typeof(W.dW) <: Number
    g3 = integrator.g(H13,p,t+c13*dt)
  else
    g3 = [integrator.g(H13[k],p,t+c13*dt) for k=1:m]
  end
  # H_3^(k)
  if typeof(W.dW) <: Number
    H14 = uprev + a141*k1*dt + b141*g1*_dW + b142*g2*_dW + b143*g3*_dW
  else
    H14 = [zero(typeof(uprev)) for k=1:m]
    for k=1:m
      H14[k] += uprev .+ a141*k1*dt
      if is_diagonal_noise(integrator.sol.prob)
        H14[k] += b141*g1[k]*_dW[k] .+ b142*g2[k]*_dW[k] .+ b143*g3[k]*_dW[k]
        for l=1:m
          if l!=k
            H14[k] += b341*g1[l]*_dW[l] .+ b342*g2[l][l]*_dW[l]
          end
        end
      else
        H14 += b141*g1[:,k]*_dW[k] .+ b142*g2[k][:,k]*_dW[k] .+ b143*g3[k][:,k]*_dW[k]
        for l=1:m
          if l!=k
            H14[k] += b341*g1[:,l]*_dW[l] .+ b342*g2[l][:,l]*_dW[l]
          end
        end
      end
    end
  end

  if typeof(W.dW) <: Number
    g4 = integrator.g(H14,p,t+c14*dt)
  else
    g4 = [integrator.g(H14[k],p,t+c14*dt) for k=1:m]
  end

  # H_3^(0) (requires H_2^(k))
  k2 = integrator.f(H02,p,t+c02*dt)
  H03 = uprev + a032*k2*dt + a031*k1*dt

  if !is_diagonal_noise(integrator.sol.prob) || typeof(W.dW) <: Number
    H03 += b031*g1*_dW
    if typeof(W.dW) <: Number
      H03 += b032*g2*_dW
    else
      for k=1:m
        H03 += b032*g2[k][:,k]*_dW[k]
      end
    end
  else
    H03 += b031*g1.*_dW
    for k=1:m
      H03 += b032*g2[k][k]*_dW[k]
    end
  end

  k3 = integrator.f(H03,p,t+c03*dt)

  # H_4^(0)
  # H04 = uprev; k4 = integrator.f(uprev,p,t+0*dt) = k1

  # H^_1^(k)
  # H21 = uprev
  # H^_2^(k) # H^_3^(k)

  if !(typeof(W.dW) <: Number)
    H22 = [uprev for k=1:m]
    H23 = [uprev for k=1:m]
    # add inbounds for speed if working properly
    for k=1:m
      for l=1:m
        if k == l
          continue
        end
        if is_diagonal_noise(integrator.sol.prob)
          @.. H22[k] += b221*g1[l]*Ihat2[k,l]/integrator.sqdt
          @.. H23[k] += b231*g1[l]*Ihat2[k,l]/integrator.sqdt
        else
          H22[k] += b221*g1[:,l]*Ihat2[k,l]/integrator.sqdt
          H23[k] += b231*g1[:,l]*Ihat2[k,l]/integrator.sqdt
        end
      end
    end
  end
  # H^_4^(k)
  # H24 = uprev

  # add stages together Eq. (5.1)
  u = uprev + (α1*k1 + α2*k2 + α3*k3 + α4*k1)*dt

  # add noise
  if typeof(W.dW) <: Number
    u += (g1*beta11+g2*beta12+g3*beta13+g4*beta14)*_dW # beta2 terms are zero by construction
  else
      if is_diagonal_noise(integrator.sol.prob)
        u += g1.*_dW*beta11
        for k=1:m
          u[k] += (g2[k][k]*beta12+g3[k][k]*beta13+g4[k][k]*beta14)*_dW[k]
          tmpg = integrator.g(H22[k],p,t)
          u[k] = u[k] + tmpg[k]*beta22*integrator.sqdt
          tmpg = integrator.g(H23[k],p,t)
          u[k] = u[k] + tmpg[k]*beta23*integrator.sqdt
        end
      else
        # non-diag noise
        for k=1:m
          g1k = @view g1[:,k]
          g2k = @view g2[k][:,k]
          g3k = @view g3[k][:,k]
          g4k = @view g4[k][:,k]
          @.. u = u + (g1k*beta11 + g2k*beta12 + g3k*beta13 + g4k*beta14)*_dW[k]
          tmpg = integrator.g(H22[k],p,t)
          @.. u = u + tmpg*beta22*integrator.sqdt
          tmpg = integrator.g(H23[k],p,t)
          @.. u = u + tmpg*beta23*integrator.sqdt
        end
      end
  end

  integrator.u = u
end

@muladd function perform_step!(integrator,cache::RSCache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack _dW,_dZ,chi1,Ihat2,tab,g1,g2,g3,g4,k1,k2,k3,H02,H03,H12,H13,H14,H22,H23,tmp1,tmpg = cache
  @unpack a021,a031,a032,a131,a141,b031,b032,b121,b131,b132,b141,b142,b143,b221,b231,b331,b332,b341,b342,α1,α2,α3,α4,c02,c03,c13,c14,beta11,beta12,beta13,beta14,beta22,beta23,NORMAL_ONESIX_QUANTILE = cache.tab

  m = length(W.dW)

  if typeof(W.dW) <: Union{SArray,Number}
    # tbd
    _dW = map(x -> calc_threepoint_random(integrator, NORMAL_ONESIX_QUANTILE, x), W.dW / sqrt(dt))
  else
    # define three-point distributed random variables
    @.. chi1 = W.dW / sqrt(dt)
    calc_threepoint_random!(_dW, integrator, NORMAL_ONESIX_QUANTILE, chi1)
    map!(x -> (x^2-dt)/2, chi1, _dW)

    # define two-point distributed random variables
    calc_twopoint_random!(_dZ,integrator,W.dZ)

    for k = 1:m
      for l = 1:k-1
        Ihat2[k, l] = _dW[k]*_dZ[l]
        Ihat2[l, k] = -_dW[k]*_dZ[l]
      end
    end
  end

  # compute stage values
  integrator.f(k1,uprev,p,t)
  integrator.g(g1,uprev,p,t)

  # H_1^(0)
  # H01 = uprev
  # H_2^(0)
  @.. H02 = uprev + a021*k1*dt

  # H_1^(k)
  # H11 = uprev
  # H_2^(k), stage 2

  if is_diagonal_noise(integrator.sol.prob)
    for k=1:m
      @.. H12[k] = uprev + b121*g1[k]*_dW[k]
    end
  else
    for k=1:m
      g1k = @view g1[:,k]
      @.. H12[k] = uprev .+ b121*g1k*_dW[k]
    end
  end

  for k=1:m
    integrator.g(g2[k],H12[k],p,t)
  end

  # H_3^(k)
  for k=1:m
    @.. H13[k] = uprev+a131*k1*dt
    if is_diagonal_noise(integrator.sol.prob)
      @.. H13[k] = H13[k] + b131*g1[k]*_dW[k] + b132*g2[k]*_dW[k]
      for l=1:m
        if l!=k
          @.. H13[k] = H13[k] + b331*g1[l]*_dW[l] + b332*g2[l]*_dW[l]
        end
      end
    else
      g1k = @view g1[:,k]
      g2k = @view g2[k][:,k]
      @.. H13[k] = H13[k] + b131*g1k*_dW[k] + b132*g2k*_dW[k]
      for l=1:m
        if l!=k
          g1l = @view g1[:,l]
          g2l = @view g2[l][:,l]
          @.. H13[k] = H13[k] + b331*g1l*_dW[l] + b332*g2l*_dW[l]
        end
      end
    end
  end

  for k=1:m
    integrator.g(g3[k],H13[k],p,t+c13*dt)
  end

  # H_3^(k)
  for k=1:m
    @.. H14[k] = uprev + a141*k1*dt
    if is_diagonal_noise(integrator.sol.prob)
      @.. H14[k] = H14[k] + b141*g1[k]*_dW[k] + b142*g2[k]*_dW[k] + b143*g3[k]*_dW[k]
      for l=1:m
        if l!=k
          @.. H14[k] = H14[k] + b341*g1[l]*_dW[l] .+ b342*g2[l]*_dW[l]
        end
      end
    else
      g1k = @view g1[:,k]
      g2k = @view g2[k][:,k]
      g3k = @view g3[k][:,k]

      @.. H14[k] = H14[k] + b141*g1k*_dW[k] + b142*g2k*_dW[k] + b143*g3k*_dW[k]
      for l=1:m
        if l!=k
          g1l = @view g1[:,l]
          g2l = @view g2[l][:,l]
          @.. H14[k] = H14[k] + b341*g1l*_dW[l] + b342*g2l*_dW[l]
        end
      end
    end
  end

  for k=1:m
    integrator.g(g4[k],H14[k],p,t+c14*dt)
  end

  # H_3^(0) (requires H_2^(k))
  integrator.f(k2,H02,p,t+c02*dt)
  @.. H03 = uprev + a032*k2*dt + a031*k1*dt

  if !is_diagonal_noise(integrator.sol.prob) || typeof(W.dW) <: Number
    mul!(tmp1,g1,_dW)
    @.. H03 = H03 + b031*tmp1
    for k=1:m
      g2k = @view g2[k][:,k]
      @.. H03 = H03 +  b032*g2k*_dW[k]
    end
  else
    @.. H03 = H03 + b031*g1*_dW
    for k=1:m
      @.. H03 = H03 + b032*g2[k][k]*_dW[k]
    end
  end

  integrator.f(k3,H03,p,t+c03*dt)

  # H_4^(0)
  # H04 = uprev; k4 = integrator.f(uprev,p,t+0*dt) = k1

  # H^_1^(k)
  # H21 = uprev
  # H^_2^(k) # H^_3^(k)

  if (!(typeof(W.dW) <: Number) || (m!=1))
    for k=1:m
     H22[k] = uprev
     H23[k] = uprev
    end
    for k=1:m
      for l=1:m
        if k == l
          continue
        end
        if is_diagonal_noise(integrator.sol.prob)
          @.. H22[k] = H22[k] + b221*g1[l]*Ihat2[k,l]/integrator.sqdt
          @.. H23[k] = H23[k] + b231*g1[l]*Ihat2[k,l]/integrator.sqdt
        else
          g1l = @view g1[:,l]
          @.. H22[k] = H22[k] + b221*g1l*Ihat2[k,l]/integrator.sqdt
          @.. H23[k] = H23[k] + b231*g1l*Ihat2[k,l]/integrator.sqdt
        end
      end
    end
  end
  # H^_4^(k)
  # H24 = uprev

  # add stages together Eq. (5.1)
  @.. u = uprev + (α1*k1 + α2*k2 + α3*k3 + α4*k1)*dt

  # add noise

  if is_diagonal_noise(integrator.sol.prob)
    @.. u = u + g1*_dW*beta11
    for k=1:m
      u[k] = u[k] + (g2[k][k]*beta12+g3[k][k]*beta13+g4[k][k]*beta14)*_dW[k]
      if m != 1
        integrator.g(tmpg,H22[k],p,t)
        u[k] = u[k] + tmpg[k]*beta22*integrator.sqdt
        integrator.g(tmpg,H23[k],p,t)
        u[k] = u[k] + tmpg[k]*beta23*integrator.sqdt
      end
    end
  else
    # non-diag noise
    for k=1:m
      g1k = @view g1[:,k]
      g2k = @view g2[k][:,k]
      g3k = @view g3[k][:,k]
      g4k = @view g4[k][:,k]
      tmpgk = @view tmpg[:,k]
      @.. u = u + (g1k*beta11 + g2k*beta12 + g3k*beta13 + g4k*beta14)*_dW[k]
      integrator.g(tmpg,H22[k],p,t)
      @.. u = u + tmpgk*beta22*integrator.sqdt
      integrator.g(tmpg,H23[k],p,t)
      @.. u = u + tmpgk*beta23*integrator.sqdt
    end
  end
end

ext: jl


text:
# using MonteCarloMeasurements

"""
    importanceSample(p(p_args), q(q_args), observed_data)

Sample from `q`, and weight the result to behave as if the sample were taken from `p`. For example,

```
julia> p = @model begin
    x ~ Normal()
    y ~ Normal(x,1) |> iid(5)
end;

julia> q = @model μ,σ begin
    x ~ Normal(μ,σ)
end;

julia> y = rand(p()).y;

julia> importanceSample(p(),q(μ=0.0, σ=0.5), (y=y,))
Weighted(-7.13971.4
,(x = -0.12280566635062592,)
````
"""
function importanceSample end

export importanceSample
@inline function importanceSample(p::ConditionalModel, q::ConditionalModel, _data)
    return _importanceSample(getmoduletypencoding(p.model), p.model, p.args, q.model, q.args, _data)
end

sourceImportanceSample(p::Model,q::Model,_data) = sourceImportanceSample(_data)(p::Model,q::Model)

export sourceImportanceSample
function sourceImportanceSample(_data)
    function(p::Model,q::Model)
        m = merge(p,q)

        _datakeys = getntkeys(_data)

        function proc(m, st::Sample)
            st.x ∈ _datakeys && return :(_ℓ += logdensity_def($(st.rhs), $(st.x)))

            if hasproperty(p.dists, st.x)
                pdist = getproperty(p.dists, st.x)
                qdist = getproperty(q.dists, st.x)
                @gensym ℓx
                result = @q begin
                    $ℓx = importanceSample($pdist, $qdist, _data)
                    _ℓ += $ℓx.ℓ
                    $(st.x) = $ℓx.val
                end
                return flatten(result)
            else return :($(st.x) = rand($(st.rhs)))
            end
            return flatten(result)
        end
        proc(m, st::Assign)     = :($(st.x) = $(st.rhs))
        proc(m, st::Return)  = :(return $(st.rhs))
        proc(m, st::LineNumber) = nothing

        body = buildSource(m, proc) |> MacroTools.flatten

        kwargs = arguments(p) ∪ arguments(q)
        kwargsExpr = Expr(:tuple,kwargs...)

        stochExpr = begin
            vals = map(sampled(q)) do x Expr(:(=), x,x) end
            Expr(:tuple, vals...)
        end

        wrap(kernel) = @q begin
            _ℓ = 0.0
            $body
            return Weighted(_ℓ, $stochExpr)
        end

        buildSource(m, proc, wrap) |> MacroTools.flatten
    end
end

@inline function importanceSample(p, q, _data)
    x = rand(q)
    ℓ = logdensity_def(p,x) - logdensity_def(q,x)
    Weighted(ℓ,x)
end

# export sourceParticleImportance
# function sourceParticleImportance(p,q;ℓ=:ℓ)
#     p = canonical(p)
#     q = canonical(q)
#     @gensym N
#     # This determines how to initialize a Particle for a given expression
#     vars(expr) = (bound(p) ∪ bound(q) ∪ stochastic(p) ∪ stochastic(q)) ∩ variables(expr)

#     procp(p, st::Follows) = :($ℓ += logdensity_def($(st.rhs), $(st.x)))
#     procp(p, st::Let)     = convert(Expr, st)
#     procp(p, st::Return)  = nothing
#     procp(p, st::LineNumber) = convert(Expr, st)

#     function procq(q, st::Follows)
#         if isempty(vars(st.rhs))
#             @q begin
#                 $(st.x) = Particles($N, $(st.rhs))
#                 $ℓ -= logdensity_def($(st.rhs), $(st.x))
#             end
#         else
#             @q begin
#                 $(st.x) = rand($(st.rhs))
#                 $ℓ -= logdensity_def($(st.rhs), $(st.x))
#             end
#         end
#     end
#     procq(q, st::Let)     = convert(Expr, st)
#     procq(q, st::Return)  = nothing
#     procq(q, st::LineNumber) = convert(Expr, st)

#     pbody = buildSource(p, procp) |> striplines
#     qbody = buildSource(q, procq) |> striplines

#     kwargs = freeVariables(q) ∪ arguments(p)
#     kwargsExpr = Expr(:tuple,kwargs...)

#     stochExpr = begin
#         vals = map(stochastic(q)) do x Expr(:(=), x,x) end
#         Expr(:tuple, vals...)
#     end

#     @gensym particleImportance
#     result = @q function $particleImportance($N, pars)
#         @unpack $kwargsExpr = pars
#         $ℓ = 0.0 * Particles($N, Uniform())
#         $qbody
#         $pbody
#         return ($ℓ, $stochExpr)
#     end

#     flatten(result)
# end

#     @gensym rand

#     flatten(@q (
#         function $rand(args...;kwargs...)
#             @unpack $argsExpr = kwargs
#             # kwargs = Dict(kwargs)
#             $body
#             $stochExpr
#         end
#     ))

# end

function merge_pqargs(src)
    @q begin
        _args = merge(_pargs, _qargs)
        $src
    end |> MacroTools.flatten
end

@gg function _importanceSample(M::Type{<:TypeLevel}, p::Model, _pargs, q::Model, _qargs, _data)
    p = type2model(p)
    q = type2model(q)

    body = sourceImportanceSample(_data)(p,q) |>
        loadvals(_qargs, _data) |>
        loadvals(_pargs, NamedTuple()) |>
        merge_pqargs

    @under_global from_type(_unwrap_type(M)) @q let M
        $body
    end
end


text:
##############################################################################
##
## Type with stored cholesky
##
##############################################################################

struct DenseCholeskyAllocatedSolver{Tc <: StridedMatrix} <: AbstractAllocatedSolver
    cholm::Tc
    function DenseCholeskyAllocatedSolver{Tc}(cholm) where {Tc <: StridedMatrix}
        size(cholm, 1) == size(cholm, 2) || throw(DimensionMismatch("chol must be square"))
        new(cholm)
    end
end
function DenseCholeskyAllocatedSolver(cholm::Tc) where {Tc <: StridedMatrix}
    DenseCholeskyAllocatedSolver{Tc}(cholm)
end

function AbstractAllocatedSolver(nls::LeastSquaresProblem{Tx, Ty, Tf, TJ, Tg}, optimizer::AbstractOptimizer{Cholesky}) where {Tx, Ty, Tf, TJ <: StridedVecOrMat, Tg}
    return DenseCholeskyAllocatedSolver(Array{eltype(nls.J)}(undef, length(nls.x), length(nls.x)))
end

##############################################################################
##
## solve J'J \ J'y by Cholesky
##
##############################################################################

function ldiv!(x::AbstractVector, J::StridedMatrix, y::AbstractVector, A::DenseCholeskyAllocatedSolver)
    cholm = A.cholm
    mul!(cholm, J',  J)
    mul!(x, J',  y)
    ldiv!(cholesky!(Symmetric(cholm), Val(true)), x)
    return x, 1
end

##############################################################################
##
## solve (J'J + diagm(damp)) \ J'y by Cholesky
##
##############################################################################

function ldiv!(x::AbstractVector, J::StridedMatrix, y::AbstractVector,
            damp::AbstractVector, A::DenseCholeskyAllocatedSolver)
    cholm = A.cholm

    # update cholm as J'J + λdtd
    mul!(cholm, J', J)
    # transform dammp
    size(cholm, 1) == length(damp) || throw(DimensionMismatch("size(chol, 1) should equal length(damp)"))
   for i in 1:size(cholm, 1)
        cholm[i, i] += damp[i]
    end

    # solve
    mul!(x, J', y)
    ldiv!(cholesky!(Symmetric(cholm)), x)
    return x, 1
end


text:
function getsaturate(g, v)::Vector{Vector{Int}}
    return collect(combinations(neighbors(g, v), 2))
end

function argmindeg(g)
    return argmin(degree(g))
end

function add_edges!(g, edges::Vector{Vector{T}}) where T
    added = Vector{Bool}(undef, 0)
    for e in edges
        push!(added, add_edge!(g, e[1], e[2]))
    end
    return added
end

function vertices_meta(g::MetaGraph)
    return Iterators.Stateful([get_prop(g, v, :v) for v in vertices(g)])
end

function edges_meta(g::MetaGraph)
    return Iterators.Stateful([edges_m2s(e) for e in edges(g)])
end

function vertex_m2s(g::MetaGraph, v)
    return get_prop(g, v, :v)
end

function vertex_s2m(g::MetaGraph, v)
    # Here the index are integers so g[v, :v] wouldn't return the node index
    # There seems to be no function to retrieve the index of an integer indexed prop
    # so I access the field `metaindex` directly
    return g.metaindex[:v][v]
end

function edge_m2s(g::MetaGraph, src, dst)
    return [get_prop(g, src, :v), get_prop(g, dst, :v)]
end

function edge_s2m(g::MetaGraph, src, dst)
    return [(g[src, :v], g[dst, :v])]
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
module GraphSearchZero
using DataStructures

export SearchProblem
# Describes how problems need to be specified in order to work with this module.
include("problem_interface.jl")
# A simple but fairly fast implementation of a search node.
export SearchNode
include("search_node.jl")
# Implementations of some graphs search algorithms.
export
    InfeasibleSearchProblemError,
    generic_graph_search,
    astar_search,
    weighted_astar_search
include("search_algorithms.jl")

end # module


text:
@testset "plot recipes" begin
    gr()

    Z = randn(5, 100)
    labels = ["a", "asdf", "b", "c", "d"]
    latentboxplot(Z,labels)

    annotatedheatmap(rand(3,3),)

    h = MVHistory()
    for _ in 1:100
        push!(h, :a, rand(3))
        push!(h, :b, rand(3))
        push!(h, :c, rand(3))
        push!(h, :c, rand(3))
        push!(h, :xrec, ones(3))
    end

    plothistory(h)

    plotreconstruction(rand(100,2), rand(100,2))
    @test true
end


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
using DataFrames

@testset "two increments, linear solver" begin
    X = Dict{Int, Vector{Float64}}(
        1 => [0.0,0.0],
        2 => [1.0,0.0],
        3 => [1.0,1.0],
        4 => [0.0,1.0])
    element = Element(Quad4, [1, 2, 3, 4])
    update!(element, "geometry", X)
    update!(element, "temperature thermal conductivity", 6.0)
    update!(element, "temperature load", 0.0 => 12.0)
    update!(element, "temperature load", 1.0 => 24.0)
    problem = Problem(Heat, "one element heat problem", 1)
    problem.properties.formulation = "2D"
    push!(problem, element)
    boundary_element = Element(Seg2, [1, 2])
    update!(boundary_element, "geometry", X)
    update!(boundary_element, "temperature 1", 0.0)
    bc = Problem(Dirichlet, "fixed", 1, "temperature")
    push!(bc, boundary_element)
    solver = Solver(Linear, problem, bc)

    solver.time = 0.0

    empty!(problem.assembly)
    solver()
    @test isapprox(solver("temperature", 0.0)[3], 1.0)

    empty!(problem.assembly)
    solver()
    @test isapprox(solver("temperature", 0.0)[3], 1.0)

    solver.time = 1.0

    empty!(problem.assembly)
    solver()
    @test isapprox(solver("temperature", 1.0)[3], 2.0)

    empty!(problem.assembly)
    solver()
    @test isapprox(solver("temperature", 1.0)[3], 2.0)

end

@testset "two increments, nonlinear solver" begin
    X = Dict{Int, Vector{Float64}}(
        1 => [0.0,0.0],
        2 => [1.0,0.0],
        3 => [1.0,1.0],
        4 => [0.0,1.0])
    element = Element(Quad4, [1, 2, 3, 4])
    update!(element, "geometry", X)
    update!(element, "temperature thermal conductivity", 6.0)
    update!(element, "temperature load", 0.0 => 12.0)
    update!(element, "temperature load", 1.0 => 24.0)
    problem = Problem(Heat, "one element heat problem", 1)
    problem.properties.formulation = "2D"
    push!(problem, element)
    boundary_element = Element(Seg2, [1, 2])
    update!(boundary_element, "geometry", X)
    update!(boundary_element, "temperature 1", 0.0)
    bc = Problem(Dirichlet, "fixed", 1, "temperature")
    push!(bc, boundary_element)
    solver = Solver(Nonlinear, problem, bc)

    solver.time = 0.0

    empty!(problem.assembly)
    solver()
    @test isapprox(solver("temperature", 0.0)[3], 1.0)

    empty!(problem.assembly)
    solver()
    @test isapprox(solver("temperature", 0.0)[3], 1.0)

    solver.time = 1.0

    empty!(problem.assembly)
    solver()
    @test isapprox(solver("temperature", 1.0)[3], 2.0)

    empty!(problem.assembly)
    solver()
    @test isapprox(solver("temperature", 1.0)[3], 2.0)

end


text:
mutable struct Scene
    nodes::Union{Nothing,Vector{Int}}
    name::Union{Nothing,String}
    extensions::Union{Nothing,Dict}
    extras::Union{Nothing,Dict}
    function Scene(; nodes=Set(), name=nothing, extensions=Dict(), extras=nothing)
        obj = new()
        obj.nodes = nothing
        obj.name = nothing
        obj.extensions = nothing
        obj.extras = nothing
        obj
    end
end

function Base.setproperty!(obj::Scene, sym::Symbol, x)
    if sym === :nodes && x !== nothing
        all(n->n ≥ 0, x) || throw(ArgumentError("the indices of each root node should be ≥ 0"))
        allunique(x) || throw(ArgumentError("nodes must be unique."))
    end
    setfield!(obj, sym, x)
end

JSON3.StructType(::Type{Scene}) = JSON3.Mutable()
StructTypes.omitempties(::Type{Scene}) = (:nodes, :name, :extensions, :extras)


text:
module Metrics

using Photon: Tensor
using Statistics
import Knet

export Metric, BinaryAccuracy, OneHotBinaryAccuracy

abstract type Metric end

"""
Binary accuracy calculation
"""
struct BinaryAccuracy
    threshold

    BinaryAccuracy(;threshold=0.5) = new(threshold)
end

function (a::BinaryAccuracy)(y_pred::Tensor, y_true::Tensor)
    y_pred = y_pred .> a.threshold
    return mean(y_true .== y_pred)
end

"""
Binary accuracy for a onehot classification.
"""
struct OneHotBinaryAccuracy end

function (a::OneHotBinaryAccuracy)(y_pred::Tensor, y_true::Tensor)
    y_pred, y_true = Knet.mat(y_pred), Knet.mat(y_true)

    y_pred = mapslices(argmax,Knet.value(Array(y_pred)),dims=1)
    y_true = Array(y_true)

    mean(y_pred .== y_true)
end

end


text:
using Wflow
using Documenter

pages = [
    "Introduction" => [
        "index.md",
        "intro/use_cases.md",
        "intro/publications.md",
        "changelog.md"
    ],
    "User guide" => [
        "user_guide/intro.md",
        "user_guide/install.md",
        "How to use wflow" => [
            "user_guide/step1_requirements.md",
            "user_guide/step2_settings_file.md",
            "user_guide/step3_input_data.md",
            "user_guide/step4_running.md",
            "user_guide/step5_output.md",
            "user_guide/additional_options.md",
            "user_guide/sample_data.md",
            ],
        "user_guide/model-setup.md",

    ],
    "Model documentation" => [
        "model_docs/intro.md",
        "model_docs/model_configurations.md",
        "Vertical concepts" => [
            "model_docs/vertical/sbm.md",
            "model_docs/vertical/hbv.md",
            "model_docs/vertical/flextopo.md",
            "model_docs/vertical/sediment.md",
            "model_docs/shared_concepts.md",
        ],
        "Lateral concepts" => [
            "model_docs/lateral/gwf.md",
            "model_docs/lateral/kinwave.md",
            "model_docs/lateral/local-inertial.md",
            "model_docs/lateral/sediment_flux.md",
        ],
        "Model parameters" => [
            "model_docs/params_vertical.md",
            "model_docs/params_lateral.md",
            "model_docs/structures.md",
        ]
    ],
]

makedocs(;
    modules = [Wflow],
    authors = "Deltares and contributors",
    repo = "https://github.com/Deltares/Wflow.jl/blob/{commit}{path}#L{line}",
    sitename = "Wflow.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://deltares.github.io/Wflow.jl",
        assets = String[],
        collapselevel = 2,
    ),
    pages = pages,
)

deploydocs(; repo = "github.com/Deltares/Wflow.jl", push_preview = true)


text:
"""
Initialize thermal measurements on the given CPU.
"""
function init_thermal(cpuid::Integer)
    LibLikwid.thermal_init(cpuid)
    return true
end

"""
Read the current temperature of the given CPU in degrees Celsius.
"""
function get_temperature(cpuid::Integer)
    init_thermal(cpuid)
    data = Ref(zero(UInt32))
    LibLikwid.thermal_read(cpuid, data)
    return Int(data[]) * u"°C"
end


text:
using BranchAndBound
using PowerModels
using JuMP
using Ipopt

const BB = BranchAndBound
const PM = PowerModels

mutable struct NodeWRMModel <: PM.AbstractWRMModel PM.@pm_fields end

function PM.constraint_model_voltage(pm::NodeWRMModel, n::Int)

    WR = var(pm, n)[:WR]
    WI = var(pm, n)[:WI]

    nw = pm.cnw
    bus_ids = ids(pm, nw, :bus)
    lookup_w_index = Dict((bi,i) for (i,bi) in enumerate(bus_ids))

    # create voltage variables
    PM.variable_bus_voltage_real(pm, bounded=true)
    PM.variable_bus_voltage_imaginary(pm, bounded=true)

    vr = var(pm, n)[:vr]
    vi = var(pm, n)[:vi]

    # constraints on diagonal
    for (i, bus) in ref(pm, nw, :bus)
        w_idx = lookup_w_index[i]

        # vr
        JuMP.@constraint(pm.model, -bus["vmax"] * vr[i] - WR[i,i] <= 0)
        JuMP.@constraint(pm.model, bus["vmax"] * vr[i] - WR[i,i] >= 0)

        # vi
        JuMP.@constraint(pm.model, -bus["vmax"] * vi[i] - WI[i,i] <= 0)
        JuMP.@constraint(pm.model, bus["vmax"] * vi[i] - WI[i,i] >= 0)
    end

    # constraints on off-diagonal
    for (i,j) in ids(pm, nw, :buspairs)
        wi_idx = lookup_w_index[i]
        wj_idx = lookup_w_index[j]
        bus = ref(pm, nw, :bus)[i]

        # impose skew symmetric
        JuMP.@constraint(pm.model, WI[i,j] + WI[j,i] == 0)

        # vr
        JuMP.@constraint(pm.model, -bus["vmax"] * vr[j] - WR[i,j] <= 0)
        JuMP.@constraint(pm.model, bus["vmax"] * vr[j] - WR[i,j] >= 0)

        # vi
        JuMP.@constraint(pm.model, -bus["vmax"] * vi[j] - WI[i,j] <= 0)
        JuMP.@constraint(pm.model, bus["vmax"] * vi[j] - WI[i,j] >= 0)
    end
end

pm_data = PM.parse_file("../data/case9.m")
pm = PM.instantiate_model(pm_data, NodeWRMModel, PM.build_opf)

vr = var(pm, pm.cnw)[:vr]
vi = var(pm, pm.cnw)[:vi]
WR = var(pm, pm.cnw)[:WR]
WI = var(pm, pm.cnw)[:WI]

pm_data = PM.parse_file("../data/case9.m")
pm = PM.instantiate_model(pm_data, NodeWRMModel, PM.build_opf)

function BB.branch(node::BB.Node)
    vr_val = BB.node_solution.(node, vr)
    vi_val = BB.node_solution.(node, vi)
    WR_val = BB.node_solution.(node, WR)
    WI_val = BB.node_solution.(node, WI)

    for (i, bus) in ref(pm, pm.cnw, :bus)
        @show abs(vr_val[i]^2 - WR_val[i,i])
        @show abs(vi_val[i]^2 - WI_val[i,i])
    end
end


text:
#
# Example script to wrap PETSc
#

using Clang.cindex
using Clang.wrap_c

PETSC_INCLUDE = "/usr/include/petsc"
MPI_INCLUDE = "/usr/include/openmpi"

petsc_header = [joinpath(PETSC_INCLUDE, "petsc.h")]

# Set up include paths
clang_includes = ASCIIString[]
push!(clang_includes, PETSC_INCLUDE)
push!(clang_includes, MPI_INCLUDE)

# Clang arguments
clang_extraargs = ["-v"]
# clang_extraargs = ["-D", "__STDC_LIMIT_MACROS", "-D", "__STDC_CONSTANT_MACROS"]

# Callback to test if a header should actually be wrapped (for exclusion)
function wrap_header(top_hdr::ASCIIString, cursor_header::ASCIIString)
    return startswith(dirname(cursor_header), PETSC_INCLUDE)
end

lib_file(hdr::ASCIIString) = "petsc"
output_file(hdr::ASCIIString) = "PETSc.jl"

function wrap_cursor(name::ASCIIString, cursor)
    exc = false
    exc |= contains(name, "MPI")
    return !exc
end

const wc = wrap_c.init(;
                        headers = petsc_header,
                        output_file = "libPETSc_h.jl",
                        common_file = "libPETSc_common.jl",
                        clang_includes      = clang_includes,
                        clang_args          = clang_extraargs,
                        header_wrapped      = wrap_header,
                        header_library      = lib_file,
                        header_outputfile   = output_file,
                        cursor_wrapped      = wrap_cursor,
                        clang_diagnostics = true)

run(wc)


text:
#
# Functor pattern
#

abstract Functor

function get{T <: Functor}(f::T) end

#
# Functor utilities
#

_map(x::Type{Val{:done}}, f::Function) = x

_map(x::Symbol, f::Function) = _map(Val{x}, f)

_map{T <: Functor}(x::T, f::Function) = T(f(get(x)))

_filter(x::Type{Val{:done}}, f::Function) = x

_filter(x::Symbol, f::Function) = _filter(Val{x}, f)

_filter{T <: Functor}(x::T, f::Function) = (f(get(x)))? x : T(nothing)

_reduce{T <: Functor}(x::T, y::Type{Val{:done}}, f::Function) = x

_reduce{T <: Functor}(x::T, y::Symbol, f::Function) = _reduce(x, Val{y}, f)

_reduce{T <: Functor}(x::T, y::T, f::Function) = T(f(get(x), get(y)))

#
# Subject
#

immutable Subject <: Functor
    value
end

function get(subject::Subject)
    return subject.value
end

function init_task(f::Function, args::Vararg)
    x = @task f(args...)
    consume(x)
    return x
end

function (subject::Subject)(f::Function, g::Function)
    return subject(init_task(f, g))
end

function (subject::Subject)(t::Task)
    return consume(t, subject)
end


text:
"""
FeatureRegistries.jl is a Julia package to create registries of a
package’s features, letting package users discover and use them more easily.
FeatureRegistries.jl lets you create tables that list package features so that

- users can explore and discover features,
- users can search features through a consistent interface; and
- developers of third-party packages can transparently add functionality

## Usage

### Package users

Package users of packages can use registries through:

- [`find`](#)/`filter`/`findfirst` to subset the list of features based on filters,
- [`info`](#)`(registry)` to get more information on a specific registry,
- `registry[id]` to inspect an entry in a registry,
- [`load`](#)`(registry)` to access a feature

See also [how to use registries](/documents/docs/using.md).

### Package developers

Developers can create registries using [`Registry`](#) and [`Field`](#), `push!(registry)`.
"""
module FeatureRegistries

using Base64
using InlineTest
using Markdown
using PrettyTables
using StructArrays

import ImageShow

# see https://github.com/JuliaImages/ImageShow.jl/pull/51
if isdefined(ImageShow, :enable_html_render)
    ImageShow.enable_html_render()
end
import Setfield

include("registry.jl")
export load

include("find.jl")
export find

include("info.jl")
export info

include("printing.jl")
include("formats.jl")

include("test.jl")

end  # module


text:
# Copyright 2021 Massachusetts Institute of Technology
# @file R2RANGE example
# @author: Qiangqiang Huang

# constructor for variables and factors
struct my_var
  name
  type
  manifold
  gt_arr
end

struct f_unary
  var
  f_type
  obs
  cov_arr
end

struct f_bin
  var1
  var2
  f_type
  obs
  cov_arr
end

struct f_kary
  var1
  var2s
  binary
  weights
  obs
  cov_arr
  f_type
end

struct incVarFactor
  vars
  factors
end

function get_kw_idx(str_arr, idx, keyword)
  kw_idx = 0
  for i = idx:length(str_arr)
    if cmp(str_arr[i], keyword) == 0
      kw_idx = i
      break
    end
  end
  return kw_idx
end

# strings for parser
str_var = "Variable"
str_factor = "Factor"
str_rbt = "Pose"
str_lmk = "Landmark"
str_f_p2p = "SE2RelativeGaussianLikelihoodFactor"
str_f_range = "SE2R2RangeGaussianLikelihoodFactor"
str_f_ada = "AmbiguousDataAssociationFactor"
str_f_prior = "UnarySE2ApproximateGaussianPriorFactor"

using Caesar
using Gadfly

# using RoME, Distributions
# using GraphPlot
using RoMEPlotting
Gadfly.set_default_plot_size(35cm,20cm)

# using Rebugger
using NPZ
using MAT

# # add multiproc support
using Distributed
addprocs(8)
@everywhere using Caesar, RoME

# # 90% precompile
for i in 1:5
  warmUpSolverJIT()
end

kde_sample = 100
post_sample = 1000
inc_step = 1
case_f = "/home/chad/codebase/NF-iSAM/example/slam/problems_for_paper/small_range_gaussian_problem/data_association/test_caesar"

fg_path = string(case_f, "/factor_graph.fg")
println(fg_path)
cnt=1
output_dir = string(case_f,"/caesar",cnt)
while true
    if !isdir(output_dir)
        mkpath(output_dir)
        break
    else
        println(output_dir)
        cnt = cnt + 1
        output_dir = string(case_f,"/caesar",cnt)
    end
end

name2var = Dict()
lmknames = []
rbtnames = []

# read factor graph
lines = readlines(fg_path)

println(lines)

factor_arr = []
p2p_idx = []
range_idx = []
ada_idx = []
nh_idx = []
prior_idx = []

factor_idx = 0
for l in lines
    str_arr = split(l,' ')
    if cmp(str_arr[1], str_var) == 0
    # parsing vriables
    gt_arr = [parse(Float64, str_arr[i]) for i in 5:length(str_arr)]
    tmp_var = my_var(str_arr[4],str_arr[2],str_arr[3],gt_arr)
    merge!(name2var, Dict(tmp_var.name=>tmp_var))
    if cmp(tmp_var.type, str_rbt) == 0
        push!(rbtnames, tmp_var.name)
    elseif cmp(tmp_var.type, str_lmk) == 0
        push!(lmknames, tmp_var.name)
    end
    elseif cmp(str_arr[1], str_factor) == 0
    factor_idx += 1
    n_f = 0
    if cmp(str_arr[2], str_f_p2p) == 0
        first_obs_idx = 5
        cov_idx = get_kw_idx(str_arr, 4, "covariance")
        obs_arr = [parse(Float64, str_arr[i]) for i in first_obs_idx:cov_idx-1]
        cov_arr = [parse(Float64, str_arr[i]) for i in cov_idx+1:length(str_arr)]
        n_f = f_bin(str_arr[3],str_arr[4], str_arr[2], obs_arr, cov_arr)
        push!(p2p_idx, factor_idx)
    elseif cmp(str_arr[2], str_f_prior) == 0
        cov_idx = get_kw_idx(str_arr, 4, "covariance")
        obs_arr = [parse(Float64, str_arr[i]) for i in 4:cov_idx-1]
        cov_arr = [parse(Float64, str_arr[i]) for i in cov_idx+1:length(str_arr)]
        n_f = f_unary(str_arr[3], str_f_prior, obs_arr, cov_arr)
        push!(prior_idx, factor_idx)
    elseif cmp(str_arr[2], str_f_range) == 0
        first_obs_idx = 5
        obs_arr = [parse(Float64, str_arr[first_obs_idx])]
        cov_arr = [parse(Float64, str_arr[first_obs_idx+1])^2]
        n_f = f_bin(str_arr[3],str_arr[4], str_arr[2], obs_arr, cov_arr)
        push!(range_idx, factor_idx)
    elseif cmp(str_arr[2], str_f_ada) == 0
        va1_idx = get_kw_idx(str_arr, 1, "Observer")
        va2_idx = get_kw_idx(str_arr, 1, "Observed")
        w_idx = get_kw_idx(str_arr, 1, "Weights")
        f_idx = get_kw_idx(str_arr, 1, "Binary")
        obs_idx = get_kw_idx(str_arr, 1, "Observation")
        sig_idx = get_kw_idx(str_arr, 1, "Sigma")
        var2_arr = [str_arr[i] for i=va2_idx+1:w_idx-1]
        w_arr = [parse(Float64, str_arr[i]) for i=w_idx+1:f_idx-1]
        obs_arr = [parse(Float64, str_arr[i]) for i=obs_idx+1:sig_idx-1]
        if cmp(str_arr[f_idx+1], str_f_range) == 0
        cov_arr = [parse(Float64, str_arr[i])^2 for i=sig_idx+1:length(str_arr)]
        else
        cov_arr = [parse(Float64, str_arr[i]) for i=sig_idx+1:length(str_arr)]
        end
        n_f = f_kary(str_arr[va1_idx+1],var2_arr, str_arr[f_idx+1],w_arr, obs_arr, cov_arr, str_f_ada)
        push!(ada_idx, factor_idx)
    else
        throw(ErrorException("Unknown factor type"))
    end
    push!(factor_arr, n_f)
    end
end

# group factors and variables incrementally
# //return values and indices of factors
incVarFactorPairs = []
# vector<pair<Values, vector<int>>> res;

if (inc_step == 0 || inc_step > length(rbtnames))
    inc_step = rbtnames.size();
    println("Reset inc_step to ", inc_step)
end

newVars = []
newFactors = []
addedRbts = []
addedLmks = []
# incremental updates with robot vars
for i in 1:length(rbtnames)
    push!(newVars, name2var[rbtnames[i]])
    push!(addedRbts, rbtnames[i])

    tmp_factor_idx = []
    for j in prior_idx
        if cmp(factor_arr[j].var, rbtnames[i]) == 0
            push!(tmp_factor_idx, j)
            println("Push a prior factor at ", rbtnames[i])
        end
    end
    prior_idx = setdiff(prior_idx, tmp_factor_idx)
    union!(newFactors, tmp_factor_idx)

    tmp_factor_idx = []
    for j in p2p_idx
        tmp_vars = [factor_arr[j].var1, factor_arr[j].var2]
        if issubset(tmp_vars, addedRbts)
            push!(tmp_factor_idx, j)
            println("Push a p2p factor between ", factor_arr[j].var1, factor_arr[j].var2)
        end
    end
    if length(tmp_factor_idx) == 0 && length(addedRbts) > 1
        throw(ErrorException("No pose2pose factors for the newly added robot variable."))
    else
        p2p_idx = setdiff(p2p_idx, tmp_factor_idx)
        union!(newFactors, tmp_factor_idx)
    end

    tmp_factor_idx = []
    for j in range_idx
        lmk_var = setdiff([factor_arr[j].var1, factor_arr[j].var2], [rbtnames[i]])
        len_lmk_var = length(lmk_var)
        if len_lmk_var == 1
            lmk_var = lmk_var[1]
            if ~(lmk_var in addedLmks)
                push!(newVars, name2var[lmk_var])
                push!(addedLmks, lmk_var)
            end
            push!(tmp_factor_idx, j)
            println("Push a range factor between ", factor_arr[j].var1, factor_arr[j].var2)
        elseif len_lmk_var == 2
            ;
        else
            throw(ErrorException("Not a binary factor."))
        end
    end
    range_idx = setdiff(range_idx, tmp_factor_idx)
    union!(newFactors, tmp_factor_idx)

    tmp_factor_idx = []
    for j in ada_idx
        if cmp(factor_arr[j].var1, rbtnames[i]) == 0
            if ~(issubset(factor_arr[j].var2s, addedRbts) || issubset(factor_arr[j].var2s, addedLmks))
                throw(ErrorException("Unobserved var is associated."))
            end
            push!(tmp_factor_idx, j)
            println("Push an ADA factor between ", factor_arr[j].var1, factor_arr[j].var2s)
        end
    end
    ada_idx = setdiff(ada_idx, tmp_factor_idx)
    union!(newFactors, tmp_factor_idx)

    # TODO: add null hypothesis

    if( i%inc_step == 0)
        inc_data = incVarFactor(deepcopy(newVars), deepcopy(newFactors))
        push!(incVarFactorPairs, inc_data)
        empty!(newVars)
        empty!(newFactors)
        println("New batch loaded.")
    end
end
println("There are ", length(incVarFactorPairs), " pairs of vars and factors.")

N = post_sample
kde_N = kde_sample
# start with an empty factor graph object
fg = initfg()
init_tree = false
tree = 1
# tree = wipeBuildNewTree!(fg)
# fieldnames(typeof(getSolverParams(fg)))
# getSolverParams(fg).graphinit = true  # init on graph first then solve on tree (default)
# getSolverParams(fg).graphinit = false # init and solve on tree
# getSolverParams(fg).useMsgLikelihoods = true
getSolverParams(fg).N = kde_N

timing = zeros(0)
added_syms = Vector{Symbol}(undef, 0)
for i in 1:length(incVarFactorPairs)
    cur_vars = incVarFactorPairs[i].vars
    cur_factors = incVarFactorPairs[i].factors
    # Add at a fixed location Prior to pin :x0 to a starting location (0,0,0)
    s_time = time_ns()
    for var in cur_vars
        if ~(var in added_syms)
        push!(added_syms, Symbol(var.name))
        if cmp(var.manifold, "SE2") == 0
            addVariable!(fg, Symbol(var.name), Pose2)
        elseif cmp(var.manifold, "R2") == 0
            addVariable!(fg, Symbol(var.name), Point2)
        else
            throw(ErrorException("Unknown manifold"))
        end
        end
    end

    for j in cur_factors
        f = factor_arr[j]
        if cmp(f.f_type, str_f_prior) == 0
            addFactor!(fg, [Symbol(f.var)], PriorPose2( MvNormal(f.obs, reshape(f.cov_arr,3,3)) ))
        elseif cmp(f.f_type, str_f_range) == 0
            addFactor!(fg, [Symbol(f.var1), Symbol(f.var2)], Pose2Point2Range(MvNormal(f.obs, diagm(f.cov_arr))))
        elseif cmp(f.f_type, str_f_p2p) == 0
            addFactor!(fg, [Symbol(f.var1), Symbol(f.var2)], Pose2Pose2(MvNormal(f.obs, reshape(f.cov_arr,3,3))))
        elseif cmp(f.f_type, str_f_ada) == 0
            symbols = [Symbol(f.var1)]
            hypos = [1.0]
            for (v_i, v) in enumerate(f.var2s)
                push!(symbols, Symbol(v))
                push!(hypos, f.weights[v_i])
            end
            # push!(symbols, Symbol("o$(j)"))
            # push!(hypos, 1.0)
            if cmp(f.binary, str_f_p2p) == 0
                # addVariable!(fg, Symbol("o$(j)"), Pose2)
                addFactor!(fg, symbols, Pose2Pose2(MvNormal(f.obs, reshape(f.cov_arr,3,3))), multihypo=hypos)
            elseif cmp(f.binary, str_f_range) == 0
                # addVariable!(fg, Symbol("o$(j)"), ContinuousScalar)
                addFactor!(fg, symbols, Pose2Point2Range(MvNormal(f.obs, diagm(f.cov_arr))), multihypo=hypos)
            else
                throw(ErrorException("Unknown binary type."))
            end
        else
            throw(ErrorException("Unknown factor types."))
        end
    end

    ## Perform inference
    if !init_tree
        tree, smt, hist = solveTree!(fg)
        init_tree = true
    else
        tree, smt, hist = solveTree!(fg, tree)
    end
    e_time = time_ns()
    elapsed = (e_time - s_time)/10^9
    append!(timing, elapsed)
    println("elapsed time: ", elapsed)

    # TODO: save variable order
    pl = plotKDE(fg, added_syms,dims=[1;2], levels=4)
    Gadfly.draw(PDF(joinpath(output_dir,"step$(i-1).pdf"), 20cm, 10cm),pl)
    res = rand(getBelief(fg, added_syms[1]) ,N)
    var_names = [string(added_syms[1])]
    if length(added_syms)>1
        for s in added_syms[2:end]
        res = vcat(res, rand(getBelief(fg, s) ,N))
        push!(var_names, string(s))
        end
    end
    npzwrite(joinpath(output_dir,"step$(i-1).npz"),res)
    outfile = joinpath(output_dir,"step$(i-1)_ordering")

    open(outfile, "w") do io
        for i in var_names # or for note in notes
            println(io, i)
        end
    end
end
npzwrite(joinpath(output_dir,"timing.npz"),timing)

max_stars_repo_name: MarineRoboticsGroup/NF-iSAM

max_issues_repo_name: MarineRoboticsGroup/NF-iSAM

max_forks_repo_name: MarineRoboticsGroup/NF-iSAM


text:
using MakieTeX
using MakieTeX.Makie, MakieTeX.CairoMakie

using Downloads

using Test

example_path = joinpath(@__DIR__, "test_images")
mkpath(example_path)

function render_texample(url)

    fig = Figure()

    lt = LTeX(fig[1, 1], CachedTeX(TeXDocument(read(Downloads.download(url), String))))

    @test true

    resize_to_layout!(fig)

    filename = splitdir(splitext(url)[1])[2]

    save(joinpath(example_path, "$filename.png"), fig; px_per_unit=3)
    save(joinpath(example_path, "$filename.pdf"), fig; px_per_unit=1)
    save(joinpath(example_path, "$filename.svg"), fig; px_per_unit=0.75)

    @test true

end
@testset "MakieTeX.jl" begin

    @testset "texample.net" begin
        names = [
            "rotated-triangle",
            "city",
            "planets",
            "model-physics",
            "smart-description",
            "or-gate",
            "polar-plot",
            "dominoes",
            "cielab"
        ]

        for name in names

            @testset "$name" begin
                render_texample("https://texample.net/media/tikz/examples/TEX/$name.tex")
            end

        end

        @testset "mandala" begin

            fig = Figure()

            @test_warn r"The PDF has more than 1 page!  Choosing the first page." LTeX(fig[1, 1], CachedTeX(TeXDocument(read(Downloads.download("https://texample.net/media/tikz/examples/TEX/mandala.tex"), String))))

            resize_to_layout!(fig)

            filename = "mandala"

            save(joinpath(example_path, "$filename.png"), fig; px_per_unit=3)
            save(joinpath(example_path, "$filename.pdf"), fig; px_per_unit=1)

        end

    end

    @testset "Rendering math and text" begin

        fig = Figure()

        @test_nowarn LTeX(fig[1, 1], "This is Lorem Ipsum")

        @test_nowarn LTeX(fig[1, 2], L"\iiint_a^{\mathbb{R}} \mathfrak D ~dt = \textbf{Poincar\'e quotient}")

        save(joinpath(example_path, "plaintex.png"), fig; px_per_unit=3)
        save(joinpath(example_path, "plaintex.pdf"), fig; px_per_unit=1)
        save(joinpath(example_path, "plaintex.svg"), fig; px_per_unit=0.75)
    end

    @testset "Layouting" begin

        @testset "Logo" begin
            fig = Figure(figure_padding = 1, resolution = (1, 1))
            @test_nowarn LTeX(fig[1, 1], "Makie\\TeX.jl")
            @test_nowarn resize_to_layout!(fig)

            save(joinpath(example_path, "logo.png"), fig; px_per_unit=3)
            save(joinpath(example_path, "logo.pdf"), fig; px_per_unit=1)
            save(joinpath(example_path, "logo.svg"), fig; px_per_unit=0.75)

            @test true
        end
    end
end


text:
const GAMMA = 6.796381511118120480641144794583e-5;

function Ham!(ham,L,jt,jcouplings)
    @views begin
        for i in 1:Int(L*(L-1)/2)
            ham .= ham .+ jcouplings[i].*jt[:,:,i]
        end
        nothing
    end
end

function incorporateNoise!(j0s, γs, τs, sigmas, disGam, jSWAP, j0)
    if sigmas[2] > 0.0
        γdist = truncated(Normal(disGam,sigmas[2]*disGam),0.0,Inf) # Prob. distribution of γ. σγ ∝ J^1
        for i in eachindex(γs)
            γs[i] = rand(γdist)
        end
    end
    for i in eachindex(j0s)
        j0s[i] = j0 * (1.0 + randn()*sigmas[1]) # This isn't truncated because σJ is so small ~0.01
    end
    for i in eachindex(τs)
        τs[i] = pi/(4.0 * jSWAP) + randn().*sigmas[3]
    end
end

function diagexp!(matrix)
    for i in 1:length(@view matrix[:,1])
        matrix[i,i] = exp(matrix[i,i])
    end
    nothing
end

function zeros!(A)
    for i in eachindex(A)
        A[i] = 0.0
    end
    nothing
end


text:
"""
```julia
Q_R2V(tau::T, delta::S) where {T,S <: Real}
```
Computes the transformation Quaternion from the R-frame to the V-frame.
"""
function Q_R2V(tau::T, delta::S) where {T,S <: Real}
  return inverse(Q_V2R(tau, delta))
end

"""
```julia
Q_V2R(tau::T, delta::S) where {T,S <: Real}
```
Computes the transformation Quaternion from the V-frame to the R-frame.
"""
function Q_V2R(tau::T, delta::S) where {T,S <: Real}
  return angleaxis2quat(pi/2 + delta,[0,1,0])*angleaxis2quat(-tau,[0,0,1])
end

#=
warn("The Q_V2B transformation has to be redone with the correct simulation state vector")
"""
```julia
Q_V2B(VehState::Vector{Float64}, GMST::Float64)
```
Computes the transformation matrix from the V-frame to the B-frame, using the
vehicle state and the GMST as inputs.
"""
function Q_V2B(VehState::Vector{T}, GMST::S) where {T,S <: Real}
  QI2R = Q_I2R(GMST)
  PosSph = PositionCart2Sph(QI2R*VehState[4:6])
  return inverse(Quaternion(VehState[10],VehState[11:13]))*(QI2R*Q_R2V(PosSph[2], PosSph[3])))
end
=#

"""
```julia
[V_air, γ, χ] = VelocityCart2Sph(VelCart::Vector{Float64})
```
Computes the spherical velocity components from its cartesian counterparts (
  in V-frame ).
"""
function VelocityCart2Sph(VelCart::Vector{T}) where {T <: Real}
  Va = norm(VelCart)
  if Va < 1e-13
    return 0, 0, 0
  else
    return [Va -asin(VelCart[3]/Va) atan2(VelCart[2], VelCart[1])]
  end
end

"""
```julia
[R, τ, δ] = PositionCart2Sph(PosCart::Vector{Float64})
```
Computes the spherical position from its cartesian components.
"""
function PositionCart2Sph(PosCart::Vector{T}) where {T <: Real}
  return [norm(PosCart), atan2(PosCart[2], PosCart[1]),
    atan2(PosCart[3], norm(PosCart))]
end

"""
```julia
PosCart = PositionSph2Cart(PosSph::Vector{Float64})
```
Computes the cartesian position from its spherical components.
"""
function PositionSph2Cart(PosSph::Vector{T}) where {T <: Real}
  return PosSph[1].*[cos(PosSph[3])*cos(PosSph[2]) cos(PosSph[3])*sin(PosSph[2])
    sin(PosSph[2])]
end

"""
```julia
[VA, γ_a, χ_a], Q_V2TA =
  AirVelocity(VG_NED::Vector{Float64}, VW_NED::Vector{Float64} = zeros(3))
```
Computes the airspeed, flight path angle, the heading, and the vertical to
airspeed based trajectory reference frame transformation matrix.
"""
function AirVelocity(VG_NED::Vector{T}, VW_NED::Vector{S} = [0,0,0]) where {T,S <: Real}
  VA_NED = VG_NED - VW_NED
  VA = norm(VA_NED)
  if VA < eps()
    γ_a = χ_a = 0
  else
    γ_a = -asin(VA_NED[3]/VA)
    χ_a = atan2(VA_NED[2], VA_NED[1])
  end
  return [VA, γ_a, χ_a], Q_V2TA(γ_a, χ_a)
end


text:
using DataFrames

df = DataFrame("event_id" => [1,2,3,4,5],
               "index" => [10,20,30,40,50],
               "data" => [100.,200.,300.,400.,500.])

@testset "util" begin
    x = collect(1:100)
    y = collect(6:9)
    ## Tests for in_range function
    xr = ATools.in_range(x, 5, 10) # interval is ( )

    @test all(y .== xr)
    @test all(ATools.in_range(x, 5, 10, ATools.ClosedBound) .== collect(5:10))
    @test all(ATools.in_range(x, 5, 10, ATools.LeftClosed ) .== collect(5:9 ))
    @test all(ATools.in_range(x, 5, 10, ATools.RightClosed) .== collect(6:10))

    @test nrow(ATools.select_values(df, x -> x < 5, "event_id")) == 4
    @test ATools.select_by_index(df,
    "index", 1) == ATools.select_by_column_value(df, "index", 1)

    df3 = ATools.select_event(df, 3)
    @test df3.index[1] == 30
    @test df3.data[1] == 300.0
    @test ATools.select_by_column_value(df, "data", 100.0).index[1]==10
    @test ATools.select_by_column_value_lt(df, "data", 200.0).index[1]==10
    @test ATools.select_by_column_value_gt(df, "data", 400.0).index[1] == 50
    @test ATools.select_by_column_value_interval(df, "data", 200.0, 400.0).index[1] ==30

    @test ATools.find_max_xy(df, "event_id", "data") == (5, 5, 500.0)
end

max_stars_repo_name: andLaing/ATools

max_issues_repo_name: andLaing/ATools

max_forks_repo_name: andLaing/ATools


text:
using VarianceComponentsHDFE
using Random, BenchmarkTools, Test
using JLD, SparseArrays, LinearAlgebra
using AlgebraicMultigrid, IterativeSolvers
# Make sure LDLFactorizations is version 0.5.0 and that the `multiple-rhs` branch is checked out
using DataDeps

register(DataDep("VarianceComponentsHDFE",
       "Downloading medium_nocontrols_pruned.csv, medium_controls_pruned.csv, large_nocontrols_pruned.csv, and large_controls_pruned.csv, huge_pruned_main.csv, large_controls_main.csv",
       ["https://vchdfe.s3-us-west-2.amazonaws.com/medium_pruned_main.tar.gz",
       "https://vchdfe.s3-us-west-2.amazonaws.com/medium_controls_pruned_main.tar.gz",
       "https://vchdfe.s3-us-west-2.amazonaws.com/large_pruned_main.tar.gz",
       "https://vchdfe.s3-us-west-2.amazonaws.com/large_controls_pruned_main.tar.gz",
       "https://vchdfe.s3-us-west-2.amazonaws.com/huge_pruned_main.tar.gz",
       "https://vchdfe.s3-us-west-2.amazonaws.com/large_controls_main.tar.gz"],
       ["e01474ebd8d7df6d645b0c865079a108c9fd33e052b7ab163bb2f5749bc58387",
       "a896b114eb952e60212196ecd1dffe2ea8bded263590a5b93e3b30b069feb0e1",
       "583ec55959b16ef93db210736c93daa8c8eab9097d4492347fe442237f6b3c2c",
       "e6ca435c7ac3bf7a8e425ab58581f36877845c099ad4da086d1990a3c270ac88",
       "5743ed199dd6e6180f6fd721d3232bec0573a05d33f358c9a228a930e517ec12",
       "82a95b368b2f6be381f10228026ebf2ffffc7e0d80bf5580ee905adf29b26128"];
       post_fetch_method= [unpack, unpack, unpack, unpack, unpack, unpack]
   ))

pkg_dir = pkgdir(VarianceComponentsHDFE)

# Obtain environment variables
use_matlab_CMG = get(ENV, "USE_MATLAB_CMG", "false") == "true" ? true : false
run_large_benchmark = get(ENV, "VCHDFE_LARGE_BENCHMARK", "false")  == "true" ? true : false
num_rhs_large = try
        parse(Int, get(ENV, "BENCHMARK_NUM_RHS_LARGE", "2"))
    catch
        2
    end

if use_matlab_CMG
    using Laplacians, Pkg
    include(string(Pkg.dir("Laplacians") , "/src/matlabSolvers.jl"))
end

# Entry point for the PkgBenchmarks call.  Can split into different files later.
include("prepare_benchmark_data.jl")
# NOTE: Suite below can assume that the `benchmark/data/...` has been filled

Random.seed!(1234)

μ = sqrt(eps())
medium_data = load(pkg_dir*"/benchmark/data/medium_main.jld")
X = medium_data["X_GroundedLaplacian"]
m,k = size(X)
X̃ = [sparse(1.0I, m,m ) X; X' spzeros(k, k)]

S_xx = medium_data["S_xx"]
X̃_regularized = [sparse(1.0I, m,m) X; X' sparse(-μ*I,k, k)]

max_rhs = 200
R_p = convert(Array{Float64,2}, bitrand(max_rhs,m))
rademacher!(R_p)

const SUITE = BenchmarkGroup()

# Prepare for direct method benchmarks
luX̃ = lu(X̃)
qrX̃ = qr(X̃)
Rz = zeros(m+k)
# RHS_aug is the first column of the mxm identity augmented by k zeros
RHS_aug = zeros(m+k)
RHS_aug[1] = 1.0

# Direct methods on the augmented system X̃
SUITE["X_tilde inplace direct solve: LU"] = @benchmarkable ldiv!($Rz, $luX̃, $RHS_aug)

# Non-inplace direct methods on the augmented system X̃
SUITE["X_tilde direct solve"] = @benchmarkable \($X̃, $RHS_aug)
SUITE["X_tilde direct solve: LU"] = @benchmarkable \($luX̃, $RHS_aug)
SUITE["X_tilde direct solve: QR"] = @benchmarkable \($qrX̃, $RHS_aug)

# Non-inplace factorizations of the augmented system X̃
SUITE["X_tilde factorization: LU"] = @benchmarkable lu($X̃)
SUITE["X_tilde factorization: QR"] = @benchmarkable qr($X̃)

# Prepare for direct method benchmarks
P = aspreconditioner(ruge_stuben(S_xx))
RHS = SparseMatrixCSC{Float64,Int64}(X[1,:])
z = 0.1.*ones(length(RHS))

# Iterative methods on the original system S_xx
SUITE["S_xx iterative solve: AMG"] = @benchmarkable cg!($z, $S_xx, $RHS, Pl = $P , log=true, maxiter=300)
if use_matlab_CMG
    sparseS_xx = sparse(S_xx)
    SUITE["S_xx iterative solve: CMG"] = @benchmarkable matlabCmgSolver($sparseS_xx, $RHS; tol=1e-6, maxits=300)
end

# Computation of the preconditioner
SUITE["S_xx precondition: AMG ruge_stuben"] = @benchmarkable aspreconditioner(ruge_stuben($S_xx))

# Prepapre for direct method (augmented system) benchmarks
ldltX̃_reg = ldlt(X̃_regularized)
luX̃_reg = lu(X̃_regularized)
qrX̃_reg = qr(X̃_regularized)
ldlX̃_reg = ldl(X̃_regularized)

# Direct methods on the regularized augmented X̃_regularized
SUITE["X_tilde_reg direct solve: LDLT"] = @benchmarkable \($ldltX̃_reg, $RHS_aug)
SUITE["X_tilde_reg direct solve: LU"] = @benchmarkable \($luX̃_reg, $RHS_aug)
SUITE["X_tilde_reg direct solve: QR"] = @benchmarkable \($qrX̃_reg, $RHS_aug)
SUITE["X_tilde_reg direct solve: LDL"] = @benchmarkable \($ldlX̃_reg, $RHS_aug)

SUITE["X_tilde_reg inplace direct solve: LDL"] = @benchmarkable ldiv!($Rz, $ldlX̃_reg, $RHS_aug)
SUITE["X_tilde_reg inplace direct solve: LU"] = @benchmarkable ldiv!($Rz, $luX̃_reg, $RHS_aug)

# Non-inplace factorizations of the regularized augmented system X̃
SUITE["X_tilde_reg factorization: LDLT"] = @benchmarkable ldlt($X̃_regularized)
SUITE["X_tilde_reg factorization: LU"] = @benchmarkable lu($X̃_regularized)
SUITE["X_tilde_reg factorization: QR"] = @benchmarkable qr($X̃_regularized)
SUITE["X_tilde_reg factorization: LDL"] = @benchmarkable ldl($X̃_regularized)

# TODO: Move around factorizations/preconditioners since they are independent of the RHS

## JLA problem with a single RHS
# JLA_RHS_aug is the first column of the mxm identity augmented by k zeros
JLA_RHS_aug = zeros(m+k)
JLA_RHS_aug[1:m] = R_p[1,:]

# Direct methods on the augmented system X̃
SUITE["JLA: X_tilde inplace direct solve: LU"] = @benchmarkable ldiv!($Rz, $luX̃, $JLA_RHS_aug)

# Non-inplace direct methods on the augmented system X̃
SUITE["JLA: X_tilde direct solve"] = @benchmarkable \($X̃, $JLA_RHS_aug)
SUITE["JLA: X_tilde direct solve: LU"] = @benchmarkable \($luX̃, $JLA_RHS_aug)
SUITE["JLA: X_tilde direct solve: QR"] = @benchmarkable \($qrX̃, $JLA_RHS_aug)

# Prepare for direct method benchmarks
# JLA_RHS is not too sparse, so we benchmark both the sparse and non-sparse representation
JLA_RHS = (R_p*X)[1,:]
JLA_RHS_sparse =  SparseMatrixCSC{Float64,Int64}(sparse(JLA_RHS))
z = 0.1.*ones(length(JLA_RHS))

# Iterative methods on the original system S_xx
SUITE["JLA: S_xx iterative solve: AMG"] = @benchmarkable cg!($z, $S_xx, $JLA_RHS, Pl = $P , log=true, maxiter=300)
SUITE["JLA: S_xx iterative solve: AMG, sparse RHS"] = @benchmarkable cg!($z, $S_xx, $JLA_RHS_sparse, Pl = $P , log=true, maxiter=300)
if use_matlab_CMG
    sparseS_xx = sparse(S_xx)
    SUITE["JLA: S_xx iterative solve: CMG"] = @benchmarkable matlabCmgSolver($sparseS_xx, $JLA_RHS; tol=1e-6, maxits=300)
    SUITE["JLA: S_xx iterative solve: CMG, sparse RHS"] = @benchmarkable matlabCmgSolver($sparseS_xx, $JLA_RHS_sparse; tol=1e-6, maxits=300)
end

# Direct methods on the regularized augmented X̃_regularized
SUITE["JLA: X_tilde_reg direct solve: LDLT"] = @benchmarkable \($ldltX̃_reg, $JLA_RHS_aug)
SUITE["JLA: X_tilde_reg direct solve: LU"] = @benchmarkable \($luX̃_reg, $JLA_RHS_aug)
SUITE["JLA: X_tilde_reg direct solve: QR"] = @benchmarkable \($qrX̃_reg, $JLA_RHS_aug)
SUITE["JLA: X_tilde_reg direct solve: LDL"] = @benchmarkable \($ldlX̃_reg, $JLA_RHS_aug)

SUITE["JLA: X_tilde_reg inplace direct solve: LDL"] = @benchmarkable ldiv!($Rz, $ldlX̃_reg, $JLA_RHS_aug)
SUITE["JLA: X_tilde_reg inplace direct solve: LU"] = @benchmarkable ldiv!($Rz, $luX̃_reg, $JLA_RHS_aug)

## LDL Factorization (regularized augmented system) with multiple right hand sides
JLA_RHS_aug_m = zeros(m+k,2)
JLA_RHS_aug_m[1:m,1:2] = R_p[1:2,:]'
Rz = zeros(size(JLA_RHS_aug_m))
SUITE["JLA: X_tilde_reg inplace direct solve: LDL, multiple RHS: 2"] = @benchmarkable ldiv!($Rz, $ldlX̃_reg, $JLA_RHS_aug_m)
JLA_RHS_aug_m = zeros(m+k,4)
JLA_RHS_aug_m[1:m,1:4] = R_p[1:4,:]'
Rz = zeros(size(JLA_RHS_aug_m))
SUITE["JLA: X_tilde_reg inplace direct solve: LDL, multiple RHS: 4"] = @benchmarkable ldiv!($Rz, $ldlX̃_reg, $JLA_RHS_aug_m)
JLA_RHS_aug_m = zeros(m+k,8)
JLA_RHS_aug_m[1:m,1:8] = R_p[1:8,:]'
Rz = zeros(size(JLA_RHS_aug_m))
SUITE["JLA: X_tilde_reg inplace direct solve: LDL, multiple RHS: 8"] = @benchmarkable ldiv!($Rz, $ldlX̃_reg, $JLA_RHS_aug_m)
JLA_RHS_aug_m = zeros(m+k,16)
JLA_RHS_aug_m[1:m,1:16] = R_p[1:16,:]'
Rz = zeros(size(JLA_RHS_aug_m))
SUITE["JLA: X_tilde_reg inplace direct solve: LDL, multiple RHS: 16"] = @benchmarkable ldiv!($Rz, $ldlX̃_reg, $JLA_RHS_aug_m)
JLA_RHS_aug_m = zeros(m+k,32)
JLA_RHS_aug_m[1:m,1:32] = R_p[1:32,:]'
Rz = zeros(size(JLA_RHS_aug_m))
SUITE["JLA: X_tilde_reg inplace direct solve: LDL, multiple RHS: 32"] = @benchmarkable ldiv!($Rz, $ldlX̃_reg, $JLA_RHS_aug_m)
JLA_RHS_aug_m = zeros(m+k,64)
JLA_RHS_aug_m[1:m,1:64] = R_p[1:64,:]'
Rz = zeros(size(JLA_RHS_aug_m))
SUITE["JLA: X_tilde_reg inplace direct solve: LDL, multiple RHS: 64"] = @benchmarkable ldiv!($Rz, $ldlX̃_reg, $JLA_RHS_aug_m)
JLA_RHS_aug_m = zeros(m+k,200)
JLA_RHS_aug_m[1:m,1:200] = R_p[1:200,:]'
Rz = zeros(size(JLA_RHS_aug_m))
SUITE["JLA: X_tilde_reg inplace direct solve: LDL, multiple RHS: 200"] = @benchmarkable ldiv!($Rz, $ldlX̃_reg, $JLA_RHS_aug_m)

# TODO JLA problem loop using iterative solver to compare against multiple RHS at the same time

##Medium Data with controls

medium_controls_data = load(pkg_dir*"/benchmark/data/medium_controls_main.jld")
Xcontrols = medium_controls_data["Xcontrols"]
m,k = size(Xcontrols)
X̃_controls = [sparse(1.0I, m,m ) Xcontrols; Xcontrols' spzeros(k, k)]
S_xx_controls = medium_controls_data["S_xx"]
X̃_regularized_controls = [sparse(1.0I, m,m) Xcontrols; Xcontrols' sparse(-μ*I,k, k)]

# Prepare for direct method benchmarks
luX̃ = lu(X̃_controls)
# TODO qr is taking too long locally, take a look at it more closely
# qrX̃ = qr(X̃_controls)
Rz = zeros(m+k)
# RHS_aug is the first column of the kxk identity augmented by m zeros
RHS_aug = zeros(m+k)
RHS_aug[1] = 1.0

# Direct methods on the augmented system X̃
SUITE["X_tilde_controls inplace direct solve: LU"] = @benchmarkable ldiv!($Rz, $luX̃, $RHS_aug)

# Non-inplace direct methods on the augmented system X̃
SUITE["X_tilde_controls direct solve"] = @benchmarkable \($X̃_controls, $RHS_aug)
SUITE["X_tilde_controls direct solve: LU"] = @benchmarkable \($luX̃, $RHS_aug)
# SUITE["X_tilde_controls direct solve: QR"] = @benchmarkable \($qrX̃, $RHS_aug)

# Non-inplace factorizations of the augmented system X̃
SUITE["X_tilde_controls factorization: LU"] = @benchmarkable lu($X̃_controls)
# SUITE["X_tilde_controls factorization: QR"] = @benchmarkable qr($X̃_controls)

# Prepare for direct method benchmarks
P = aspreconditioner(ruge_stuben(S_xx_controls))
RHS = SparseMatrixCSC{Float64,Int64}(Xcontrols[1,:])
z = 0.1.*ones(length(RHS))

# Iterative methods on the original system S_xx
SUITE["S_xx_controls iterative solve: AMG"] = @benchmarkable cg!($z, $S_xx_controls, $RHS, Pl = $P , log=true, maxiter=300)

# Computation of the preconditioner
SUITE["S_xx_controls precondition: AMG ruge_stuben"] = @benchmarkable aspreconditioner(ruge_stuben($S_xx_controls))

# Prepapre for direct method (augmented system) benchmarks
ldltX̃_reg = ldlt(X̃_regularized_controls)
luX̃_reg = lu(X̃_regularized_controls)
# qrX̃_reg = qr(X̃_regularized_controls)
ldlX̃_reg = ldl(X̃_regularized_controls)

# Direct methods on the regularized augmented X̃_regularized
SUITE["X_tilde_reg_controls direct solve: LDLT"] = @benchmarkable \($ldltX̃_reg, $RHS_aug)
SUITE["X_tilde_reg_controls direct solve: LU"] = @benchmarkable \($luX̃_reg, $RHS_aug)
# SUITE["X_tilde_reg_controls direct solve: QR"] = @benchmarkable \($qrX̃_reg, $RHS_aug)
SUITE["X_tilde_reg_controls direct solve: LDL"] = @benchmarkable \($ldlX̃_reg, $RHS_aug)

SUITE["X_tilde_reg_controls inplace direct solve: LDL"] = @benchmarkable ldiv!($Rz, $ldlX̃_reg, $RHS_aug)
SUITE["X_tilde_reg_controls inplace direct solve: LU"] = @benchmarkable ldiv!($Rz, $luX̃_reg, $RHS_aug)

# Non-inplace factorizations of the regularized augmented system X̃
SUITE["X_tilde_reg_controls factorization: LDLT"] = @benchmarkable ldlt($X̃_regularized_controls)
SUITE["X_tilde_reg_controls factorization: LU"] = @benchmarkable lu($X̃_regularized_controls)
# SUITE["X_tilde_reg_controls factorization: QR"] = @benchmarkable qr($X̃_regularized_controls)
SUITE["X_tilde_reg_controls factorization: LDL"] = @benchmarkable ldl($X̃_regularized_controls)

if run_large_benchmark

    large_data = load(pkg_dir*"/benchmark/data/large_main.jld")
    X = large_data["X_GroundedLaplacian"]
    m,k = size(X)
    X̃ = Symmetric([sparse(1.0I, m,m ) X; spzeros(k, k+m)])

    S_xx = large_data["S_xx"]
    X̃_regularized = Symmetric([sparse(1.0I, m,m) X; spzeros(k,m) sparse(-μ*I,k, k)])

    R_p = convert(Array{Float64,2}, bitrand(num_rhs_large,m))
    rademacher!(R_p)

    # Prepare for direct method benchmarks
    Rz = zeros(m+k)
    # RHS_aug is the first column of the mxm identity augmented by k zeros
    RHS_aug = zeros(m+k)
    RHS_aug[1] = 1.0

    # Prepare for iterative method benchmarks
    P = aspreconditioner(ruge_stuben(S_xx))
    RHS = SparseMatrixCSC{Float64,Int64}(X[1,:])
    z = 0.1.*ones(length(RHS))

    # Iterative methods on the original system S_xx
    SUITE["Large: S_xx iterative solve: AMG"] = @benchmarkable cg!($z, $S_xx, $RHS, Pl = $P , log=true, maxiter=300)
    if use_matlab_CMG
        sparseS_xx = sparse(S_xx)
        SUITE["Large: S_xx iterative solve: CMG"] = @benchmarkable matlabCmgSolver($sparseS_xx, $RHS; tol=1e-6, maxits=300)
    end

    # Computation of the preconditioner
    SUITE["Large: S_xx precondition: AMG ruge_stuben"] = @benchmarkable aspreconditioner(ruge_stuben($S_xx))

    # Prepapre for direct method (augmented system) benchmarks
    ldltX̃_reg = ldlt(X̃_regularized)
    ldlX̃_reg = ldl(X̃_regularized)

    # Direct methods on the regularized augmented X̃_regularized
    SUITE["Large: X_tilde_reg direct solve: LDLT"] = @benchmarkable \($ldltX̃_reg, $RHS_aug)
    # SUITE["Large: X_tilde_reg direct solve: LDL"] = @benchmarkable \($ldlX̃_reg, $RHS_aug)

    SUITE["Large: X_tilde_reg inplace direct solve: LDL"] = @benchmarkable ldiv!($Rz, $ldlX̃_reg, $RHS_aug)

    # Non-inplace factorizations of the regularized augmented system X̃
    SUITE["Large: X_tilde_reg factorization: LDLT"] = @benchmarkable ldlt($X̃_regularized)
    SUITE["Large: X_tilde_reg factorization: LDL"] = @benchmarkable ldl($X̃_regularized)

    ## JLA problem with a single RHS
    # JLA_RHS_aug is the first column of the mxm identity augmented by k zeros
    JLA_RHS_aug = zeros(m+k)
    JLA_RHS_aug[1:m] = R_p[1,:]

    # JLA_RHS is not too sparse, so we benchmark both the sparse and non-sparse representation
    JLA_RHS = (R_p*X)[1,:]
    JLA_RHS_sparse =  SparseMatrixCSC{Float64,Int64}(sparse(JLA_RHS))
    z = 0.1.*ones(length(JLA_RHS))

    # Iterative methods on the original system S_xx
    SUITE["Large: JLA: S_xx iterative solve: AMG"] = @benchmarkable cg!($z, $S_xx, $JLA_RHS, Pl = $P , log=true, maxiter=300)
    z = 0.1.*ones(length(JLA_RHS))
    SUITE["Large: JLA: S_xx iterative solve: AMG, sparse RHS"] = @benchmarkable cg!($z, $S_xx, $JLA_RHS_sparse, Pl = $P , log=true, maxiter=300)
    if use_matlab_CMG
        sparseS_xx = sparse(S_xx)
        SUITE["Large: JLA: S_xx iterative solve: CMG"] = @benchmarkable matlabCmgSolver($sparseS_xx, $JLA_RHS; tol=1e-6, maxits=300)
        SUITE["Large: JLA: S_xx iterative solve: CMG, sparse RHS"] = @benchmarkable matlabCmgSolver($sparseS_xx, $JLA_RHS_sparse; tol=1e-6, maxits=300)
    end

    # Direct methods on the regularized augmented X̃_regularized
    SUITE["Large: JLA: X_tilde_reg direct solve: LDLT"] = @benchmarkable \($ldltX̃_reg, $JLA_RHS_aug)
    # SUITE["Large: JLA: X_tilde_reg direct solve: LDL"] = @benchmarkable \($ldlX̃_reg, $JLA_RHS_aug)

    SUITE["Large: JLA: X_tilde_reg inplace direct solve: LDL"] = @benchmarkable ldiv!($Rz, $ldlX̃_reg, $JLA_RHS_aug)

    ## LDL Factorization (regularized augmented system) with multiple right hand sides
    JLA_RHS_aug_m = zeros(m+k,num_rhs_large)
    JLA_RHS_aug_m[1:m,1:num_rhs_large] = R_p[1:num_rhs_large,:]'
    Rz = zeros(size(JLA_RHS_aug_m))
    SUITE["JLA: X_tilde_reg inplace direct solve: LDL, multiple RHS"] = @benchmarkable ldiv!($Rz, $ldlX̃_reg, $JLA_RHS_aug_m)

    # TODO Add back the controls later

end


text:
"""
    $(TYPEDEF)

In-memory representation of Tiff Tags, which are essentially key value pairs.
The `data` field can either be a `String`, a `Number`, an Array of bitstypes, or
a [`RemoteData`](@ref) type.

$FIELDS
"""
struct Tag{T}
    tag::UInt16
    data::T
end

Tag(tag::TiffTag, data) = Tag(UInt16(tag), data)
Tag(tag::UInt16, data::String) = Tag{String}(tag, data)
Tag(tag::UInt16, data::T) where {T <: Enum} = Tag{UInt16}(tag, UInt16(data))
function Tag(tag::UInt16, data::AbstractVector{T}) where {T}
    if length(data) == 1
        Tag{T}(tag, first(data))
    else
        Tag{Vector{T}}(tag, data)
    end
end

function Base.getproperty(t::Tag{<: AbstractString}, f::Symbol)
    if f == :data
        # strip NUL terminators from the end of strings
        strip(getfield(t, :data), '\0')
    else
        getfield(t, f)
    end
end

Base.length(t::Tag{<: AbstractVector}) = length(t.data)
Base.length(t::Tag{<: AbstractString}) = (endswith(t.data, '\0') ? length(t.data) : length(t.data) + 1)
Base.length(t::Tag) = 1

Base.eltype(::Tag{T}) where {T} = T
Base.eltype(::Tag{<: AbstractVector{T}}) where {T} = T
Base.eltype(t::Tag{RemoteData}) = t.data.datatype

load(tf::TiffFile, t::Tag) = t

function load(tf::TiffFile{O}, t::Tag{RemoteData{O}}) where {O <: Unsigned}
    T = t.data.datatype
    N = t.data.count
    rawdata = Vector{UInt8}(undef, bytes(T)*N)

    pos = position(tf.io)
    seek(tf, t.data.position)
    read!(tf, rawdata)

    # if this datatype is comprised of multiple bytes and this file needs to be
    # bitswapped then we'll need to reverse the byte order inside each datatype
    # unit
    if tf.need_bswap && bytes(T) >= 2
        reverse!(rawdata)
        data = Array{T}(reverse(reinterpret(T, rawdata)))
    elseif T == String
        data = String(rawdata)
    elseif T == Any
        data = Array{Any}(rawdata)
    else
        data = Array{T}(reinterpret(T, rawdata))
    end

    if N == 1
        data = first(data)
    end
    seek(tf, pos)

    Tag(t.tag, data)
end

bytes(x::Type) = sizeof(x)
bytes(::Type{Any}) = 1
bytes(::Type{String}) = 1
bytes(::Type{RemoteData{O}}) where {O} = one(O)
bytes(::Type{<: AbstractVector{T}}) where {T} = bytes(T)

function Base.read(tf::TiffFile{O}, ::Type{Tag}) where O <: Unsigned
    tag = read(tf, UInt16)
    datatype = read(tf, UInt16)
    count = read(tf, O)
    data = Vector{UInt8}(undef, sizeof(O))
    read!(tf, data)

    T = Any
    if datatype in keys(tiff_to_julian)
        T = tiff_to_julian[datatype]
    end

    nbytes = bytes(T) * count
    if nbytes <= sizeof(O)
        if tf.need_bswap
            reverse!(view(data, 1:nbytes))
        end
        if T == String
            return Tag(tag, String(data))
        elseif T == Any
            return Tag(tag, Array{Any}(data))
        elseif count == 1
            return Tag(tag, first(reinterpret(T, data)))
        else
            return Tag(tag, Array(reinterpret(T, data)[1:Int(count)]))
        end
    else
        (tf.need_bswap) && reverse!(data)
        return Tag(tag, RemoteData(first(reinterpret(O, data)), T, count))
    end
end

_showdata(io::IO, t::Tag{RemoteData{O}}) where {O} = print(io, "REMOTE@$(t.data.position) $(t.data.datatype)[] len=$(t.data.count)")
_showdata(io::IO, t::Tag{<: AbstractString}) = print(io, "\"", first(t.data, 20), (length(t.data) > 20) ? "..." : "", "\"")
_showdata(io::IO, t::Tag{<: AbstractVector}) = print(io, "$(eltype(t.data))[", join(t.data[1:min(5, end)], ", "), (length(t.data) > 5) ? ", ..." : "", "]")
_showdata(io::IO, t::Tag) = print(io, t.data)

function Base.show(io::IO, t::Tag)
    print(io, "Tag(")
    try
        print(io, TiffTag(t.tag), ", ")
    catch
        print(io, "UNKNOWN($(Int(t.tag))), ")
    end
    if t.tag == Int(COMPRESSION)
        print(io, CompressionType(t.data))
    else
        _showdata(io, t)
    end
    print(io, ")")
end

"""
    write(tf, t)

Write tag `t` to the tiff file `tf`. Returns `true` if the tag data fit
entirely in the IFD space and was written to disk. Otherwise it returns false.
"""
function Base.write(tf::TiffFile{O}, t::Tag{T}) where {O <: Unsigned, T}
    # if the data are too large to fit then we'll need to skip writing this tag
    # for now until we know the length of the entire IFD
    if length(t)*bytes(T) > sizeof(O)
        _writeblank(tf)
        return false
    end

    # add NUL terminator to the end of Strings that don't have it already
    data = (T == String && !endswith(t.data, '\0')) ? t.data * "\0" : t.data

    write(tf, t.tag)
    write(tf, julian_to_tiff[eltype(t)])
    write(tf, O(length(t)))
    nbytes = write(tf.io, data)

    # write padding
    if nbytes < sizeof(O)
        write(tf, fill(0x00, sizeof(O) - nbytes))
    end
    true
end

function Base.write(tf::TiffFile{O}, t::Tag{RemoteData{O}}) where {O <: Unsigned}
    write(tf, t.tag)
    write(tf, julian_to_tiff[t.data.datatype])
    write(tf, t.data.count)
    write(tf, t.data.position)
    true
end

function Base.write(tf::TiffFile{O}, t::Tag, offset) where {O <: Unsigned}
    write(tf, Tag(t.tag, RemoteData(O(offset), eltype(t), O(length(t)))))
end

# Base.write(tf::TiffFile, t::Tag) = error("Tag offsets must agree with file offsets")

function _writeblank(tf::TiffFile{O}) where O
    write(tf, UInt32(0))
    write(tf, zero(O))
    write(tf, zero(O))
end

const tagfields = fieldnames(Tag)
function Base.:(==)(t1::Tag{O1}, t2::Tag{O2}) where {O1, O2}
    return O1 == O2 && getproperty.(Ref(t1), tagfields) == getproperty.(Ref(t2), tagfields)
end


text:
export get_rt_d, plot_kinx, plot_kinxt, plot_kiny, plot_kinyt, plot_kin_binx, plot_kin_biny

#=
Reaction Time Calculations
=#

function get_rt_d(kin::Array{Kinematic,1},thres::Float64)

    rt=zeros(Int64,length(kin))

    for i=1:length(kin)
        c1=(mean(kin[i].px[1:10]),mean(kin[i].py[1:10]))

        t_rt=0
        for j=10:length(kin[i].px)
            if (kin[i].px[j]-c1[1])^2 + (kin[i].py[j]-c1[2])^2 > thres*thres
                t_rt=j
                break
            end
        end

        rt[i]=t_rt
    end
    rt
end

#=
Velocity Calculations
=#

function kin_ci(y::Array{Array{Float64,1},1},l::Int64)

    myci=zeros(Float64,l,2)
    mymean=zeros(Float64,l)
    for i=1:l
        test=zeros(Float64,0)
        for j=1:length(y)
            if length(y[j])>i
                push!(test,y[j][i])
            end
        end
        #temp=ci(OneSampleTTest(test),tail = :both)
        temp=quantile(test,[.25,.75])
        myci[i,1]=temp[1]
        myci[i,2]=temp[2]
        mymean[i]=mean(test)
    end
    (myci,mymean)
end

function get_kinx(kin::Array{Kinematic,1},tp::Array{Int64,1},myind::Array{Int64,1})
    vx=[kin[i].vx[tp[i]-9:end] for i=1:length(kin[myind])]
end

function get_kiny(kin::Array{Kinematic,1},tp::Array{Int64,1},myind::Array{Int64,1})
    vy=[kin[i].vy[tp[i]-9:end] for i=1:length(kin[myind])]
end

function kin_hist(y::Array{Array{Float64,1},1},l::Int64,bins::Array{Float64,1})

    test=zeros(Float64,0)
    for j=1:length(y)
        if length(y[j])>l
            push!(test,y[j][l])
        end
    end
    hist(test,bins)[2]
end

#=
General Plotting
=#

function shaded_line(y::Array{Float64,1},myci::Array{Float64,2},ax)
    ax[:plot](collect(1:length(y)),y,color="blue")
    ax[:fill_between](collect(1:length(y)),myci[:,1],myci[:,2],color="blue",alpha=0.5)
    nothing
end

#=
Velocity Plotting
=#

function plot_kinx(kin::Array{Kinematic,1},tp::Array{Int64,1},myind::Array{Int64,1},l::Int64,ax)
    y=get_kinx(kin,tp,myind)
    (myci,mymean)=kin_ci(y,l)
    shaded_line(mymean,myci,ax)
end

function plot_kinxt(kin::Array{Kinematic,1},tp::Array{Int64,1},myind::Array{Int64,1},l::Int64,ax)

    x=[collect(1:length(tp[i]-9:length(kin[i].vx))) for i=1:length(kin[myind])]
    y=get_kinx(kin,tp,myind)
    xy=[hcat(x[i],y[i]) for i=1:length(x)]

    c1=C.LineCollection(xy,color="black",linewidth=1.0,alpha=.2)
    ax[:add_collection](c1)
    ax[:set_xlim]([0.0, l])
    ax[:set_ylim]([-1.0, 1.0])

    nothing
end

function plot_kiny(kin::Array{Kinematic,1},tp::Array{Int64,1},myind::Array{Int64,1},l::Int64,ax)
    y=get_kiny(kin,tp,myind)
    (myci,mymean)=kin_ci(y,l)
    shaded_line(mymean,myci,ax)
end

function plot_kinyt(kin::Array{Kinematic,1},tp::Array{Int64,1},myind::Array{Int64,1},l::Int64,ax)

    x=[collect(1:length(tp[i]-9:length(kin[i].vy))) for i=1:length(kin[myind])]
    y=get_kiny(kin,tp,myind)
    xy=[hcat(x[i],y[i]) for i=1:length(x)]

    c1=C.LineCollection(xy,color="black",linewidth=1.0,alpha=.2)
    ax[:add_collection](c1)
    ax[:set_xlim]([0.0, l])
    ax[:set_ylim]([-1.0, 1.0])

    nothing
end

function plot_kin_binx(kin::Array{Kinematic,1},tp::Array{Int64,1},myind::Array{Int64,1},l::Int64,bins::Array{Float64,1},ax)
    y=get_kinx(kin,tp,myind)
    myhist=kin_hist(y,l,bins)
    ax[:bar](bins[2:end],myhist,bins[2]-bins[1], color="black")
    ax[:set_xlim]([bins[1], bins[end]])
end

function plot_kin_biny(kin::Array{Kinematic,1},tp::Array{Int64,1},myind::Array{Int64,1},l::Int64,bins::Array{Float64,1},ax)
    y=get_kiny(kin,tp,myind)
    myhist=kin_hist(y,l,bins)
    ax[:bar](bins[2:end],myhist,bins[2]-bins[1], color="black")
    ax[:set_xlim]([bins[1], bins[end]])
end


text:
################################################################################
#                                     DCMs
################################################################################

#                              Introduced in v0.4
# ==============================================================================
@deprecate   dcm2angle   dcm_to_angle
@deprecate   dcm2quat    dcm_to_quat

################################################################################
#                             Euler Angle and Axis
################################################################################

#                              Introduced in v0.4
# ==============================================================================
@deprecate   angleaxis2quat     angleaxis_to_quat

################################################################################
#                                 Euler Angles
################################################################################

#                              Introduced in v0.4
# ==============================================================================
@deprecate   angle2dcm         angle_to_dcm
@deprecate   smallangle2dcm    smallangle_to_dcm
@deprecate   angle2quat        angle_to_quat
@deprecate   smallangle2quat   smallangle_to_quat
@deprecate   angle2rot         angle_to_rot
@deprecate   smallangle2rot    smallangle_to_rot

################################################################################
#                                 Quaternions
################################################################################

#                              Introduced in v0.4
# ==============================================================================
@deprecate   quat2angle       quat_to_angle
@deprecate   quat2angleaxis   quat_to_angleaxis
@deprecate   quat2dcm         quat_to_dcm


text:
@everywhere function spmd_test1()
    barrier(;tag=:b1)

    if myid() == 1
        @assert SPMD.recvfrom(2) == "Hello from 2"
        println("SPMD: Passed send/recv")
    elseif myid() == 2
        data = "Hello from 2"
        sendto(1, data)
    end

    stime = rand(1:5)
#    println("Sleeping for $stime seconds")
    sleep(stime)
    barrier(;tag=:b2)

    bcast_val = nothing
    if myid() == 1
        bcast_val = rand(2)
    end

    bcast_val = bcast(bcast_val, 1)

    if myid() == 1
        @assert bcast_val == SPMD.recvfrom(2)
        println("SPMD: Passed broadcast")
    elseif myid() == 2
        sendto(1, bcast_val)
    end

    barrier()

    scatter_data = nothing
    if myid() == 1
        scatter_data = rand(Int8, nprocs())
    end
    lp = scatter(scatter_data, 1, tag=1)

    if myid() == 1
        @assert scatter_data[2:2] == SPMD.recvfrom(2)
        println("SPMD: Passed scatter 1")
    elseif myid() == 2
        sendto(1, lp)
    end

    scatter_data = nothing
    if myid() == 1
        scatter_data = rand(Int8, nprocs()*2)
    end
    lp = scatter(scatter_data, 1, tag=2)

    if myid() == 1
        @assert scatter_data[3:4] == SPMD.recvfrom(2)
        println("SPMD: Passed scatter 2")
    elseif myid() == 2
        sendto(1, lp)
    end

    gathered_data = gather(myid(), 1, tag=3)
    if myid() == 1
        @assert gathered_data == procs()
        println("SPMD: Passed gather 1")
    end

    gathered_data = gather([myid(), myid()], 1, tag=4)
    if myid() == 1
        @assert gathered_data == [[p,p] for p in procs()]
        println("SPMD: Passed gather 2")
    end
end

spmd(spmd_test1)

# Test running only on the workers using the spmd function.

# define the function everywhere
@everywhere function foo_spmd(d_in, d_out, n)
    pids=sort(vec(procs(d_in)))
    pididx = findfirst(isequal(myid()), pids)
    mylp = localpart(d_in)
    localsum = 0

    # Have each node exchange data with its neighbors
    n_pididx = pididx+1 > length(pids) ? 1 : pididx+1
    p_pididx = pididx-1 < 1 ? length(pids) : pididx-1

#    println(p_pididx, " p", pids[p_pididx], " ", n_pididx, " p", pids[n_pididx])
#    println(mylp)

    for i in 1:n
        sendto(pids[n_pididx], mylp[2])
        sendto(pids[p_pididx], mylp[1])

        mylp[2] = SPMD.recvfrom(pids[p_pididx])
        mylp[1] = SPMD.recvfrom(pids[n_pididx])

#        println(mylp)

        barrier(;pids=pids)
        localsum = localsum + mylp[1] + mylp[2]
    end

    # finally store the sum in d_out
    d_out[:L] = localsum
end

# run foo_spmd on all workers, many of them, all concurrently using implictly different contexts.
in_arrays = map(x->DArray(I->fill(myid(), (map(length,I)...,)), (nworkers(), 2), workers(), [nworkers(),1]), 1:8)
out_arrays = map(x->ddata(), 1:8)

@sync for i in 1:8
    @async spmd(foo_spmd, in_arrays[i], out_arrays[i], nworkers(); pids=workers())
end
for i in 1:8
    @test Any[sum(workers())*2 for i in 1:nworkers()] == gather(out_arrays[i])
end

println("SPMD: Passed testing of spmd function run concurrently")

# run concurrently with explictly different contexts

# define the function everywhere
@everywhere function foo_spmd2(d_in, d_out, n)
    pids=sort(vec(procs(d_in)))
    pididx = findfirst(isequal(myid()), pids)
    mylp = localpart(d_in)

    # see if we have a value in the local store.
    store = context_local_storage()

    localsum = get!(store, :LOCALSUM, 0)

    # Have each node exchange data with its neighbors
    n_pididx = pididx+1 > length(pids) ? 1 : pididx+1
    p_pididx = pididx-1 < 1 ? length(pids) : pididx-1

    for i in 1:n
        sendto(pids[n_pididx], mylp[2])
        sendto(pids[p_pididx], mylp[1])

        mylp[2] = SPMD.recvfrom(pids[p_pididx])
        mylp[1] = SPMD.recvfrom(pids[n_pididx])

        barrier(;pids=pids)
        localsum = localsum + mylp[1] + mylp[2]
    end

    # finally store the sum in d_out
    d_out[:L] = localsum
    store[:LOCALSUM] = localsum
end

in_arrays = map(x->DArray(I->fill(myid(), (map(length,I)...,)), (nworkers(), 2), workers(), [nworkers(),1]), 1:8)
out_arrays = map(x->ddata(), 1:8)
contexts = map(x->context(workers()), 1:8)

@sync for i in 1:8
    @async spmd(foo_spmd2, in_arrays[i], out_arrays[i], nworkers(); pids=workers(), context=contexts[i])
end
# Second run will add the value stored in the previous run.
@sync for i in 1:8
    @async spmd(foo_spmd2, in_arrays[i], out_arrays[i], nworkers(); pids=workers(), context=contexts[i])
end

for i in 1:8
    @test Any[2*sum(workers())*2 for i in 1:nworkers()] == gather(out_arrays[i])
end

# verify localstores with appropriate context store values exist.
@everywhere begin
    if myid() != 1
        local n = 0
        for (k,v) in DistributedArrays.SPMD.map_ctxts
            store = v.store
            localsum = store[:LOCALSUM]
            if localsum != 2*sum(workers())*2
                println("localsum ", localsum, " != $(2*sum(workers())*2)")
                error("localsum mismatch")
            end
            n += 1
        end
        @assert n == 8
    end
end

# close the contexts
foreach(x->close(x), contexts)

# verify that the localstores have been deleted.
@everywhere begin
    @assert isempty(DistributedArrays.SPMD.map_ctxts)
end

println("SPMD: Passed spmd function with explicit context run concurrently")


text:
module BenchmarkTools

using Compat

# `show` compatibility for pre-JuliaLang/julia#16354 builds
if VERSION < v"0.5.0-dev+4305"
    Base.get(io::IO, setting::Symbol, default::Bool) = default
end

##############
# Parameters #
##############

include("parameters.jl")

export loadparams!

##############
# Trial Data #
##############

include("trials.jl")

export gctime,
       memory,
       allocs,
       params,
       ratio,
       judge,
       isinvariant,
       isregression,
       isimprovement,
       rmskew!,
       rmskew,
       trim

##################
# Benchmark Data #
##################

include("groups.jl")

export BenchmarkGroup,
       invariants,
       regressions,
       improvements,
       @tagged,
       addgroup!,
       leaves

######################
# Execution Strategy #
######################

include("execution.jl")

export tune!,
       warmup,
       @benchmark,
       @benchmarkable

##########################################
# Plotting Facilities (loaded on demand) #
##########################################

loadplotting() = include(joinpath(dirname(@__FILE__), "plotting.jl"))

end # module BenchmarkTools


text:
module Day8

run_part2(s) = "not done"

SEGMENT_COUNTS = Dict(2 => [1],
                      3 => [7],
                      4 => [4],
                      5 => [2, 3, 5],
                      6 => [0, 6, 9],
                      7 => [8])

possible_digits_by_segment_count(s) = SEGMENT_COUNTS[length(s)]

output_values(s) = split(s, "|")[2] |>
    strip |>
    y -> split(y, " ")

run_part1(data_dir) = joinpath(data_dir, "day8.txt") |>
    readlines |>
    ls -> output_values.(ls) |>
    Base.Iterators.flatten |>
    collect |>
    ovs -> possible_digits_by_segment_count.(ovs) |>
    ds -> length.(ds) |>
    ds -> filter(y -> y == 1, ds) |>
    sum

DISPLAY = Dict(Set("abcefg") => 0,
               Set("cf") => 1,
               Set("acdeg") => 2,
               Set("acdfg") => 3,
               Set("bcdf") => 4,
               Set("abdfg") => 5,
               Set("abdefg") => 6,
               Set("acf") => 7,
               Set("abcdefg") => 8,
               Set("abcdfg") => 9)

struct Permuation
    a
    b
    c
    d
    e
    f
    g
end

end