text:
# functions for doing the parallel communication required for reverse mode
# calculations
# this is different than the communication in parallel_rev2.jl because it
# sends both q and res_bar at the beginning, and therefore does not require
# a parallel communication at the end to return q_bar to the process that
# ownes it.  Unfortunately, this makes the code the physics modules must
# implement a little bit different in reverse mode.  The idea is:
#
# Primal: q_local + q_remote -> res_local
# Reverse: q_local + q_remote + res_bar_local + res_bar_remote -> q_bar_local
#
# Notice that only q_bar_local needs to be computed because q_bar_remote will
# be computed by the process that ownes it.  It is able to do this because
# it has access to res_bar_local and res_bar_remote

"""
  Similar to [`startSolutionExchange`](@ref), this function starts
  communciation of the solution `eqn.res_bar` and possibly `eqn.q`.
  See [`finishExchangeData_rev`](@ref) for a way to finish communication
  for reverse-mode calculations. `eqn.shared_data` and `eqn.shared_data_bar`
  are used for the communications.

  **Inputs**

   * mesh: an AbstractMesh
   * sbp: an SBP operator
   * eqn: an AbstractSolutionData
   * opts: options dictionary

  **Keyword arguments**

   * send_q: if true, start communication for the solution `eqn.q`.
             Default true.  Callers should carefully consider if this can
             be set to false, which may improve performance.
   * wait: wait for sends and receives to finish before exiting, default false

  **Implementation Notes**

  This function is slightly sub-optimial in that it posts all the sends and
  receives for the solution and then posts the sends and receives for
  `res_bar`.  It would be slightly better to post all receives first and then
  all sends.
"""
function startSolutionExchange_rev2(mesh::AbstractMesh, sbp::AbstractOperator,
                                  eqn::AbstractSolutionData, opts;
                                  send_q=true, wait=false)

  if mesh.npeers == 0
    return nothing
  end

  if send_q
    startSolutionExchange(mesh, sbp, eqn, opts, wait=wait)
  end

  pdata = eqn.shared_data_bar[1].pdata  # assume all are same
  if pdata == PARALLEL_DATA_FACE
    populate_buffer = getSendDataFace_resbar
  elseif pdata == PARALLEL_DATA_ELEMENT
    populate_buffer = getSendDataElement_resbar
  else
    throw(ErrorException("unsupported parallel_type = $(getParallelDataString(pdata))"))
  end

  exchangeData(mesh, sbp, eqn, opts, eqn.shared_data_bar, populate_buffer,
               wait=wait)

  return nothing
end

"""
  Similar to [`finishExchangeData`](@ref), but accepts two vectors of
  `SharedFaceData` objects, one for the solution data and the other for some
  `_bar` quantity.  The `calc_func` must have the signature:

  ```
    calc_func(mesh, sbp, eqn, opts, data::SharedFaceData,
              data_bar::SharedFaceData)
  ```

  **Inputs**

   * mesh: an AbstractMesh
   * sbp: an SBPOperator
   * eqn: an AbstractSolutionData
   * opts: the options dictonary
   * calc_func: function that does calculations for a set of shared faces
               described by a single SharedFaceData.  It must have the signature
               calc_func(mesh, sbp, eqn, opts, data::SharedFaceData)

  **Inputs/Outputs**

   * shared_data: vector of `SharedFaceData` for solution variables
   * shared_data_bar: vector of `SharedFaceData` for `_bar` variables

  In reality, whether or not `shared_data` contains the solution and
  `shared_data_bar` contains `_bar` variables of some kind does not matter.
  This function is simply the extensions of the `finishExchangeData` where
  two vector of `SharedFaceData` are provided.

  **Implementation Notes**

  The implementation uses MPI's Waitany function on `shared_data_bar`.  This
  is because the more likely case for reverse-mode calculations is that the
  solution was sent previously and only `shared_face_bar` needs to be
  waited on.

"""
function finishExchangeData_rev2(mesh, sbp, eqn, opts,
                                shared_data::Vector{SharedFaceData{T}},
                                shared_data_bar::Vector{SharedFaceData{T2}},
                                calc_func::Function) where {T, T2}

  npeers = length(shared_data)
  @assert npeers == length(shared_data_bar)
  val = assertReceivesConsistent(shared_data)
  val_bar = assertReceivesConsistent(shared_data_bar)

  for i=1:npeers

    #TODO: there is a slightly better algorithm for waiting on two sets of
    #      Requestions: Waitany on the first, check if the corresponding
    #      entry in the second set has finished.  If not, Waitany on the
    #      second set.  Continue ping-ponging back and forth unitil a matched
    #      set is ready
    if val_bar == 0  # request have not been waited on previously
      eqn.params.time.t_wait += @elapsed idx = waitAnyReceive(shared_data_bar)
    else
      idx = i
    end

    data_bar_idx = shared_data_bar[idx]
    if data_bar_idx.pdata == PARALLEL_DATA_FACE && val_bar == 0
      # permute the received nodes to be in the elementR orientation
      permuteinterface!(mesh.sbpface, data_bar_idx.interfaces, data_bar_idx.q_recv)
    end

    data_idx = shared_data[idx]
    if val == 0
      waitReceive(data_idx)
    end

    if data_idx.pdata == PARALLEL_DATA_FACE && val == 0
      # permute the received nodes to be in the elementR orientation
      permuteinterface!(mesh.sbpface, data_idx.interfaces, data_idx.q_recv)
    end

    calc_func(mesh, sbp, eqn, opts, data_idx, data_bar_idx)

  end

  # wait on all sends to complete.  This is not strictly necessary, but
  # avoids a problem with exchangeData_rev where it expects the sends
  # to already have been waited on
  waitAllSends(shared_data_bar)

  return nothing
end

"""
  Like [`getSendDataFace`](@ref), but takes the data from `eqn.res_bar`
  and puts it in the send buffer.
"""
function getSendDataFace_resbar(mesh::AbstractMesh, sbp::AbstractOperator,
                         eqn::AbstractSolutionData, opts, data::SharedFaceData)

  idx = data.peeridx
  bndryfaces = data.bndries_local
  boundaryinterpolate!(mesh.sbpface, bndryfaces, eqn.res_bar, data.q_send)

  return nothing
end

"""
  This function populates the send buffer from eqn.q for
  `PARALLEL_DATA_ELEMENT`

  Inputs:

    mesh: a mesh
    sbp: an SBP operator
    eqn: an AbstractSolutionData
    opts: options dictonary

  Inputs/Outputs:

    data: a SharedFaceData.  data.q_send will be overwritten
"""
function getSendDataElement_resbar(mesh::AbstractMesh, sbp::AbstractOperator,
                         eqn::AbstractSolutionData, opts, data::SharedFaceData)

  # copy data into send buffer
  idx = data.peeridx
  local_els = mesh.local_element_lists[idx]
  send_buff = data.q_send
  for j=1:length(local_els)
    el_j = local_els[j]
    for k=1:size(eqn.q, 2)
      for p=1:size(eqn.q, 1)
        send_buff[p, k, j] = eqn.res_bar[p, k, el_j]
      end
    end
  end

  return nothing
end


text:
# Autogenerated wrapper script for OpenSSL_jll for i686-w64-mingw32
export libcrypto, libssl, openssl

JLLWrappers.@generate_wrapper_header("OpenSSL")
JLLWrappers.@declare_library_product(libcrypto, "libcrypto-1_1.dll")
JLLWrappers.@declare_library_product(libssl, "libssl-1_1.dll")
JLLWrappers.@declare_executable_product(openssl)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libcrypto,
        "bin\\libcrypto-1_1.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libssl,
        "bin\\libssl-1_1.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        openssl,
        "bin\\openssl.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
export LBFGSOperator, InverseLBFGSOperator, diag, diag!

"A data type to hold information relative to LBFGS operators."
mutable struct LBFGSData{T,I<:Integer}
  mem::I
  scaling::Bool
  scaling_factor::T
  damped::Bool
  σ₂::T
  σ₃::T
  s::Vector{Vector{T}}
  y::Vector{Vector{T}}
  ys::Vector{T}
  α::Vector{T}
  a::Vector{Vector{T}}
  b::Vector{Vector{T}}
  insert::I
  Ax::Vector{T}
end

function LBFGSData(
  T::DataType,
  n::I;
  mem::I = 5,
  scaling::Bool = true,
  damped::Bool = false,
  inverse::Bool = true,
  σ₂::Float64 = 0.99,
  σ₃::Float64 = 10.0,
) where {I<:Integer}
  LBFGSData{T,I}(
    max(mem, 1),
    scaling,
    convert(T, 1),
    damped,
    convert(T, σ₂),
    convert(T, σ₃),
    [zeros(T, n) for _ = 1:mem],
    [zeros(T, n) for _ = 1:mem],
    zeros(T, mem),
    inverse ? zeros(T, mem) : zeros(T, 0),
    inverse ? Vector{T}(undef, 0) : [zeros(T, n) for _ = 1:mem],
    inverse ? Vector{T}(undef, 0) : [zeros(T, n) for _ = 1:mem],
    1,
    Vector{T}(undef, n),
  )
end

LBFGSData(n::I; kwargs...) where {I<:Integer} = LBFGSData(Float64, n; kwargs...)

"A type for limited-memory BFGS approximations."
mutable struct LBFGSOperator{T,S,I<:Integer,F,Ft,Fct} <: AbstractLinearOperator{T}
  nrow::I
  ncol::I
  symmetric::Bool
  hermitian::Bool
  prod!::F    # apply the operator to a vector
  tprod!::Ft    # apply the transpose operator to a vector
  ctprod!::Fct   # apply the transpose conjugate operator to a vector
  Mv::S # storage vector for prod!
  Mtu::S # storage vector for tprod!
  Maw::S # storage vector for ctprod!
  inverse::Bool
  data::LBFGSData{T,I}
  nprod::I
  ntprod::I
  nctprod::I
end

LBFGSOperator{T}(
  nrow::I,
  ncol::I,
  symmetric::Bool,
  hermitian::Bool,
  prod!::F,
  tprod!::Ft,
  ctprod!::Fct,
  Mv::S,
  Mtu::S,
  Maw::S,
  inverse::Bool,
  data::LBFGSData{T,I},
) where {T,S,I<:Integer,F,Ft,Fct} =
  LBFGSOperator{T,S,I,F,Ft,Fct}(nrow, ncol, symmetric, hermitian, prod!, tprod!, ctprod!, Mv, Mtu, Maw, inverse, data, 0, 0, 0)

"""
    InverseLBFGSOperator(T, n, [mem=5; scaling=true])
    InverseLBFGSOperator(n, [mem=5; scaling=true])
Construct a limited-memory BFGS approximation in inverse form. If the type `T`
is omitted, then `Float64` is used.
"""
function InverseLBFGSOperator(T::DataType, n::I; kwargs...) where {I<:Integer}
  kwargs = Dict(kwargs)
  delete!(kwargs, :inverse)
  lbfgs_data = LBFGSData(T, n; inverse = true, kwargs...)

  function lbfgs_multiply(res::AbstractVector, data::LBFGSData, x::AbstractArray, αm, βm)
    # Multiply operator with a vector.
    # See, e.g., Nocedal & Wright, 2nd ed., Procedure 7.4, p. 178.

    q = data.Ax # tmp vector
    q .= x
    # q .= αm .* x

    for i = 1:(data.mem)
      k = mod(data.insert - i - 1, data.mem) + 1
      if data.ys[k] != 0
        αk = dot(data.s[k], q) / data.ys[k]
        data.α[k] = αk
        for j ∈ eachindex(q)
          q[j] -= αk * data.y[k][j]
        end
      end
    end

    data.scaling && (q .*= data.scaling_factor)

    for i = 1:(data.mem)
      k = mod(data.insert + i - 2, data.mem) + 1
      if data.ys[k] != 0
        αk = data.α[k]
        β = αk - dot(data.y[k], q) / data.ys[k]
        for j ∈ eachindex(q)
          q[j] += β * data.s[k][j]
        end
      end
    end
    if βm != zero(T)
      res .= αm .* q .+ βm .* res
    else
      res .= αm .* q
    end
  end

  prod! = @closure (res, x, α, β) -> lbfgs_multiply(res, lbfgs_data, x, α, β)
  Mv = similar(lbfgs_data.Ax)
  return LBFGSOperator{T}(n, n, true, true, prod!, prod!, prod!, Mv, Mv, Mv, true, lbfgs_data)
end

InverseLBFGSOperator(n::Int; kwargs...) = InverseLBFGSOperator(Float64, n; kwargs...)

"""
    LBFGSOperator(T, n; [mem=5, scaling=true])
    LBFGSOperator(n; [mem=5, scaling=true])
Construct a limited-memory BFGS approximation in forward form. If the type `T`
is omitted, then `Float64` is used.
"""
function LBFGSOperator(T::DataType, n::I; kwargs...) where {I<:Integer}
  kwargs = Dict(kwargs)
  delete!(kwargs, :inverse)
  lbfgs_data = LBFGSData(T, n; inverse = false, kwargs...)

  function lbfgs_multiply(res::AbstractVector, data::LBFGSData, x::AbstractArray, α, β)
    # Multiply operator with a vector.
    # See, e.g., Nocedal & Wright, 2nd ed., Procedure 7.6, p. 184.

    q = data.Ax
    q .= x

    data.scaling && (q ./= data.scaling_factor)

    # B = B₀ + Σᵢ (bᵢbᵢ' - aᵢaᵢ').
    for i = 1:(data.mem)
      k = mod(data.insert + i - 2, data.mem) + 1
      if data.ys[k] != 0
        ax = dot(data.a[k], x)
        bx = dot(data.b[k], x)
        for j ∈ eachindex(q)
          q[j] += bx * data.b[k][j] - ax * data.a[k][j]
        end
      end
    end
    if β != zero(T)
      res .= α .* q .+ β .* res
    else
      res .= α .* q
    end
  end

  prod! = @closure (res, x, α, β) -> lbfgs_multiply(res, lbfgs_data, x, α, β)
  Mv = similar(lbfgs_data.Ax)
  return LBFGSOperator{T}(n, n, true, true, prod!, prod!, prod!, Mv, Mv, Mv, false, lbfgs_data)
end

LBFGSOperator(n::I; kwargs...) where {I<:Integer} = LBFGSOperator(Float64, n; kwargs...)

"""
    push!(op, s, y)
    push!(op, s, y, α, g)
Push a new {s,y} pair into a L-BFGS operator.
The second calling sequence is used in inverse LBFGS updating in conjunction with damping,
where α is the most recent steplength and g the gradient used when solving `d=-Hg`.
In forward updating with damping, it is not necessary to supply α and g.
"""
function push!(op::LBFGSOperator, s::Vector, y::Vector, α::Real = 1.0, g::Vector = eltype(y)[])
  ys = dot(y, s)
  σ₂ = op.data.σ₂
  σ₃ = op.data.σ₃

  if op.data.damped
    # Powell's damped update strategy
    Bs = op.inverse ? (-α * g) : (op * s)
    sBs = dot(s, Bs)
    damp = false
    if ys < (1 - σ₂) * sBs
      θ = σ₂ * sBs / (sBs - ys)
      damp = true
    elseif ys > (1 + σ₃) * sBs
      θ = σ₃ * sBs / (ys - sBs)
      damp = true
    end
    if damp
      y = θ * y + (1 - θ) * Bs  # damped y
      ys = θ * ys + (1 - θ) * sBs
    end
  else
    if ys <= eps(eltype(op))
      # op.counters.rejects +=1
      return op
    end
  end

  # op.counters.updates += 1
  data = op.data
  insert = data.insert

  data.s[insert] .= s
  data.y[insert] .= y
  data.ys[insert] = ys
  op.data.scaling && (op.data.scaling_factor = ys / dot(y, y))

  # Update arrays a and b used in forward products.
  if !op.inverse
    data.b[insert] .= y / sqrt(ys)

    for i = 1:(data.mem)
      k = mod(insert + i - 1, data.mem) + 1
      if data.ys[k] != 0
        data.a[k] .= data.s[k] / data.scaling_factor  # B₀ = I / γ.

        for j = 1:(i - 1)
          l = mod(insert + j - 1, data.mem) + 1
          if data.ys[l] != 0
            data.a[k] .+= dot(data.b[l], data.s[k]) * data.b[l]
            data.a[k] .-= dot(data.a[l], data.s[k]) * data.a[l]
          end
        end
        data.a[k] /= sqrt(dot(data.s[k], data.a[k]))
      end
    end
  end

  op.data.insert = mod(insert, data.mem) + 1
  return op
end

"""
    diag(op)
    diag!(op, d)
Extract the diagonal of a L-BFGS operator in forward mode.
"""
function diag(op::LBFGSOperator{T}) where {T}
  d = Vector{T}(undef, op.nrow)
  diag!(op, d)
end

function diag!(op::LBFGSOperator{T}, d) where {T}
  op.inverse && throw(
    LinearOperatorException("only the diagonal of a forward L-BFGS approximation is available"),
  )
  data = op.data

  fill!(d, 1)
  data.scaling && (d ./= data.scaling_factor)

  for i = 1:(data.mem)
    k = mod(data.insert + i - 2, data.mem) + 1
    if data.ys[k] != 0
      for j = 1:(op.nrow)
        d[j] = d[j] + data.b[k][j]^2 - data.a[k][j]^2
      end
    end
  end
  return d
end

"""
    reset!(data)
Resets the given LBFGS data.
"""
function reset!(data::LBFGSData{T,I}, inverse::Bool) where {T,I<:Integer}
  for i = 1:(data.mem)
    fill!(data.s[i], 0)
    fill!(data.y[i], 0)
    if !inverse
      fill!(data.a[i], 0)
      fill!(data.b[i], 0)
    end
  end
  fill!(data.ys, 0)
  fill!(data.α, 0)
  data.scaling_factor = T(1)
  data.insert = 1
  return data
end

"""
    reset!(op)
Resets the LBFGS data of the given operator.
"""
function reset!(op::LBFGSOperator)
  reset!(op.data, op.inverse)
  op.nprod = 0
  op.ntprod = 0
  op.nctprod = 0
  return op
end


text:
function AUSM_plusup(Qbase, Qcon, cellxmax, cellymax, vecAx, vecAy, specific_heat_ratio, Minf, volume, nval)
    E_adv_hat = zeros(cellxmax+1,   cellymax, nval)
    F_adv_hat = zeros(  cellxmax, cellymax+1, nval)
    g         = specific_heat_ratio

    for i in 2:cellxmax+1 -1
        for j in 2:cellymax -1
            # i-1セル
            cell_vecAx = zeros(2)
            cell_vecAx[1] = 0.5*(vecAx[i-1,j,1]+vecAx[i,j,1])
            cell_vecAx[2] = 0.5*(vecAx[i-1,j,2]+vecAx[i,j,2])

            rhoL = Qbase[i-1,j,1]
            UL = Qcon[i-1,j,2] / Qcon[i-1,j,1]*cell_vecAx[1] + Qcon[i-1,j,3] / Qcon[i-1,j,1]*cell_vecAx[2]
            pL = Qbase[i-1,j,4]

            # iセル
            cell_vecAx = zeros(2)
            cell_vecAx[1] = 0.5*(vecAx[i,j,1]+vecAx[i+1,j,1])
            cell_vecAx[2] = 0.5*(vecAx[i,j,2]+vecAx[i+1,j,2])

            rhoR = Qbase[i,j,1]
            #UR = Qbase[i,j,2]*cell_vecAx[1] + Qbase[i,j,3]*cell_vecAx[2]
            UR = Qcon[i,j,2] / Qcon[i,j,1]*cell_vecAx[1] + Qcon[i,j,3] / Qcon[i,j,1]*cell_vecAx[2]
            pR = Qbase[i,j,4]

            # AUSM+up
            mdot, ph = AUSM_plusup_half(rhoL, rhoR, UL, UR, pL, pR, Minf, g)
#=
            if rhoL == 0
                println(" L  ")
                println(i)
                println(j)
            end
            if rhoR == 0
                println("   ")
                println(i)
                println(j)
            end
=#
            # flux half
            temp_vecX = zeros(nval)
            temp_vecX[2] = vecAx[i,j,1]
            temp_vecX[3] = vecAx[i,j,2]

            Lpsi = zeros(nval)
            Rpsi = zeros(nval)

            for l in 1:nval
                Lpsi[l] = Qcon[i-1,j,l] / Qcon[i-1,j,1]
                Rpsi[l] = Qcon[i,j,l]   / Qcon[i,j,1]
            end
            Lpsi[4] = (Qcon[i-1,j,4] + Qbase[i-1,j,4]) / Qcon[i-1,j,1]
            Rpsi[4] = (Qcon[i,j,4]   + Qbase[i,j,4])   / Qcon[i,j,1]

            if mdot > 0
                for l in 1:nval
                    E_adv_hat[i,j,l] = mdot * Lpsi[l] + ph * temp_vecX[l]
                end
            else
                for l in 1:nval
                    E_adv_hat[i,j,l] = mdot * Rpsi[l] + ph * temp_vecX[l]
                end
            end
        end
    end

    for i in 2:cellxmax -1
        for j in 2:cellymax+1 -1

            # j-1セル
            cell_vecAy = zeros(2)
            cell_vecAy[1] = 0.5*(vecAy[i,j-1,1]+vecAy[i,j,1])
            cell_vecAy[2] = 0.5*(vecAy[i,j-1,2]+vecAy[i,j,2])

            rhoL = Qbase[i,j-1,1]
            VL = Qbase[i,j-1,2]*cell_vecAy[1] + Qbase[i,j-1,3]*cell_vecAy[2]
            pL = Qbase[i,j-1,4]

            # jセル
            cell_vecAy = zeros(2)
            cell_vecAy[1] = 0.5*(vecAy[i,j,1]+vecAy[i,j+1,1])
            cell_vecAy[2] = 0.5*(vecAy[i,j,2]+vecAy[i,j+1,2])

            rhoR = Qbase[i,j,1]
            VR = Qbase[i,j,2]*cell_vecAy[1] + Qbase[i,j,3]*cell_vecAy[2]
            pR = Qbase[i,j,4]

            # AUSM+up
            mdot, ph = AUSM_plusup_half(rhoL, rhoR, VL, VR, pL, pR, Minf, g)
            #=
            if i == 100 && j == 50
                println(" y ")
                println(mdot)
                println(ph)
            end
            =#

            # flux half
            temp_vecY = zeros(nval)
            temp_vecY[2] = vecAy[i,j,1]
            temp_vecY[3] = vecAy[i,j,2]

            Lpsi = zeros(nval)
            Rpsi = zeros(nval)

            for l in 1:nval
                Lpsi[l] = Qcon[i,j-1,l] / Qcon[i,j-1,1]
                Rpsi[l] = Qcon[i,j,l]   / Qcon[i,j,1]
            end
            Lpsi[4] = (Qcon[i,j-1,4] + Qbase[i,j-1,4]) / Qcon[i,j-1,1]
            Rpsi[4] = (Qcon[i,j,4]   + Qbase[i,j,4])   / Qcon[i,j,1]

            if mdot > 0
                for l in 1:nval
                    F_adv_hat[i,j,l] = mdot * Lpsi[l] + ph * temp_vecY[l]
                end
            else
                for l in 1:nval
                    F_adv_hat[i,j,l] = mdot * Rpsi[l] + ph * temp_vecY[l]
                end
            end

            #=
            if i == 35 && j ==2
                println(" y ")
                println(mdot)
                println(temp_vecY)
                println(Lpsi)
                println(Rpsi)
            end
            =#

        end
    end

    return E_adv_hat, F_adv_hat
end

function AUSM_plusup_half(rhoL, rhoR, UL, UR, pL, pR, Minf, g)
    # param
    beta  = 1/8
    kp    = 0.25
    sigma = 1.0
    ku    = 0.75

    # L, R
    aL = (g * pL / rhoL)^0.5
    aR = (g * pR / rhoR)^0.5

    # half
    ah   = 0.5*( aL + aR )
    rhoh = 0.5*( rhoL + rhoR )

    Mbar = (( UL^2 + UR^2 ) / ( 2 * ah^2 ))^0.5
    Mo   = (min(1,max( Mbar^2, Minf^2 )))^0.5
    fa   = Mo * (2-Mo)

    alpha = 3/16 * ( -4 + 5*fa )

    ML = UL/ah
    MR = UR/ah

    M_p4 = 0
    M_m4 = 0
    p_p5 = 0
    p_m5 = 0
    if abs(ML) >= 1
        M_p4 = 0.5*(ML + abs(ML))
        p_p5 = 0.5*(ML + abs(ML)) / ML
    else
        Mtp  = 0.25*(ML + 1)^2
        Mtm  = -0.25*(ML - 1)^2
        #M2m1 = (ML^2 - 1)^2
        M_p4 = Mtp * (1 - 16*beta*Mtm)
        #M_p4 = (1-ML)* (Mtp + beta*M2m1) + ML*0.5*(ML + abs(ML))
        p_p5 = Mtp * ((2-ML) - 16*alpha*ML*Mtm)
    end

    if abs(MR) >= 1
        M_m4 = 0.5*(MR - abs(MR))
        p_m5 = 0.5*(MR - abs(MR)) / MR
    else
        Mtp  = 0.25*(MR + 1)^2
        Mtm  = -0.25*(MR - 1)^2
        M_m4 = Mtm * (1 + 16*beta*Mtp)
        p_m5 = Mtm * ((-2-MR) + 16*alpha*MR*Mtp)
    end

    #=
    Mtp  = 0.25*(ML + 1)^2
    Mtm  = -0.25*(ML - 1)^2
    M2m1 = (ML^2 - 1)^2
    M_p4 = (1-ML)* (Mtp + beta*M2m1) + ML*0.5*(ML + abs(ML))

    Mtp  = 0.25*(MR + 1)^2
    Mtm  = -0.25*(MR - 1)^2
    M2m1 = (MR^2 - 1)^2
    M_m4 = (1-MR)* (Mtm - beta*M2m1) + MR*0.5*(MR - abs(MR))
    =#

    # M half
    Mp = -kp * max( 1 - sigma*Mbar^2, 0) * (pR - pL)/(rhoh*ah^2)

    # Mh = M_p4 + M_m4 + Mp/fa
    Mh = M_p4 + M_m4

    # mdot half
    mdot = ah * Mh
    if Mh > 0
        mdot = mdot * rhoL
    else
        mdot = mdot * rhoR
    end

    # p half
    pu = -ku * p_p5 * p_m5 * (rhoL + rhoR) * ah *(UR - UL)

    ph = p_p5*pL + p_m5*pR + fa * pu
    return mdot, ph, ah, Mh
end

function AUSM_plus(Qbase, Qcon, cellxmax, cellymax, vecAx, vecAy, specific_heat_ratio, volume, nval)
    E_adv_hat = zeros(cellxmax+1,   cellymax, nval)
    F_adv_hat = zeros(  cellxmax, cellymax+1, nval)
    g         = specific_heat_ratio

    for i in 2:cellxmax+1 -1
        for j in 2:cellymax -1
            # i-1セル
            cell_vecAx = zeros(2)
            cell_vecAx[1] = 0.5*(vecAx[i-1,j,1]+vecAx[i,j,1])
            cell_vecAx[2] = 0.5*(vecAx[i-1,j,2]+vecAx[i,j,2])

            rhoL = Qbase[i-1,j,1]
            UL = Qcon[i-1,j,2] / Qcon[i-1,j,1]*cell_vecAx[1] + Qcon[i-1,j,3] / Qcon[i-1,j,1]*cell_vecAx[2]
            pL = Qbase[i-1,j,4]

            # iセル
            cell_vecAx = zeros(2)
            cell_vecAx[1] = 0.5*(vecAx[i,j,1]+vecAx[i+1,j,1])
            cell_vecAx[2] = 0.5*(vecAx[i,j,2]+vecAx[i+1,j,2])

            rhoR = Qbase[i,j,1]
            #UR = Qbase[i,j,2]*cell_vecAx[1] + Qbase[i,j,3]*cell_vecAx[2]
            UR = Qcon[i,j,2] / Qcon[i,j,1]*cell_vecAx[1] + Qcon[i,j,3] / Qcon[i,j,1]*cell_vecAx[2]
            pR = Qbase[i,j,4]

            # AUSM+
            mdot, ph = AUSM_plus_half(rhoL, rhoR, UL, UR, pL, pR, g)

            # flux half
            temp_vecX = zeros(nval)
            temp_vecX[2] = vecAx[i,j,1]
            temp_vecX[3] = vecAx[i,j,2]

            Lpsi = zeros(nval)
            Rpsi = zeros(nval)

            for l in 1:nval
                Lpsi[l] = Qcon[i-1,j,l] / Qcon[i-1,j,1]
                Rpsi[l] = Qcon[i,j,l]   / Qcon[i,j,1]
            end
            Lpsi[4] = (Qcon[i-1,j,4] + Qbase[i-1,j,4]) / Qcon[i-1,j,1]
            Rpsi[4] = (Qcon[i,j,4]   + Qbase[i,j,4])   / Qcon[i,j,1]

            if mdot > 0
                for l in 1:nval
                    E_adv_hat[i,j,l] = mdot * Lpsi[l] + ph * temp_vecX[l]
                end
            else
                for l in 1:nval
                    E_adv_hat[i,j,l] = mdot * Rpsi[l] + ph * temp_vecX[l]
                end
            end
        end
    end

    for i in 2:cellxmax -1
        for j in 2:cellymax+1 -1

            # j-1セル
            cell_vecAy = zeros(2)
            cell_vecAy[1] = 0.5*(vecAy[i,j-1,1]+vecAy[i,j,1])
            cell_vecAy[2] = 0.5*(vecAy[i,j-1,2]+vecAy[i,j,2])

            rhoL = Qbase[i,j-1,1]
            VL = Qbase[i,j-1,2]*cell_vecAy[1] + Qbase[i,j-1,3]*cell_vecAy[2]
            pL = Qbase[i,j-1,4]

            # jセル
            cell_vecAy = zeros(2)
            cell_vecAy[1] = 0.5*(vecAy[i,j,1]+vecAy[i,j+1,1])
            cell_vecAy[2] = 0.5*(vecAy[i,j,2]+vecAy[i,j+1,2])

            rhoR = Qbase[i,j,1]
            VR = Qbase[i,j,2]*cell_vecAy[1] + Qbase[i,j,3]*cell_vecAy[2]
            pR = Qbase[i,j,4]

            # AUSM+up
            mdot, ph = AUSM_plus_half(rhoL, rhoR, VL, VR, pL, pR, g)

            # flux half
            temp_vecY = zeros(nval)
            temp_vecY[2] = vecAy[i,j,1]
            temp_vecY[3] = vecAy[i,j,2]

            Lpsi = zeros(nval)
            Rpsi = zeros(nval)

            for l in 1:nval
                Lpsi[l] = Qcon[i,j-1,l] / Qcon[i,j-1,1]
                Rpsi[l] = Qcon[i,j,l]   / Qcon[i,j,1]
            end
            Lpsi[4] = (Qcon[i,j-1,4] + Qbase[i,j-1,4]) / Qcon[i,j-1,1]
            Rpsi[4] = (Qcon[i,j,4]   + Qbase[i,j,4])   / Qcon[i,j,1]

            if mdot > 0
                for l in 1:nval
                    F_adv_hat[i,j,l] = mdot * Lpsi[l] + ph * temp_vecY[l]
                end
            else
                for l in 1:nval
                    F_adv_hat[i,j,l] = mdot * Rpsi[l] + ph * temp_vecY[l]
                end
            end
        end
    end

    return E_adv_hat, F_adv_hat
end

function AUSM_plus_half(rhoL, rhoR, UL, UR, pL, pR, g)
    # param
    beta  = 1/8
    alpha = 3/16

    # L, R
    aL = (g * pL / rhoL)^0.5
    aR = (g * pR / rhoR)^0.5

    # half
    ah   = 0.5*( aL + aR )
    rhoh = 0.5*( rhoL + rhoR )

    ML = UL/ah
    MR = UR/ah

    M_p4 = 0
    M_m4 = 0
    p_p5 = 0
    p_m5 = 0
    if abs(ML) >= 1
        M_p4 = 0.5*(ML + abs(ML))
        p_p5 = 0.5*(ML + abs(ML)) / ML
    else
        M_p4 = 0.5*(ML+1)^2 + beta*(ML^2-1)^2
        p_p5 = 0.25*(ML+1)^2 * (2-ML) + alpha*ML*(ML^2-1)^2
    end

    if abs(MR) >= 1
        M_m4 = 0.5*(MR - abs(MR))
        p_m5 = 0.5*(MR - abs(MR)) / MR
    else
        M_m4 = -0.5*(MR-1)^2 - beta*(MR^2-1)^2
        p_m5 = 0.25*(MR-1)^2 * (2+MR) - alpha*MR*(MR^2-1)^2
    end

    # Mh = M_p4 + M_m4 + Mp/fa
    Mh = M_p4 + M_m4

    # mdot half
    mdot = ah * Mh
    if Mh > 0
        mdot = mdot * rhoL
    else
        mdot = mdot * rhoR
    end

    # p half
    ph = p_p5*pL + p_m5*pR

    return mdot, ph
end

max_stars_repo_name: hide-dog/test

max_issues_repo_name: hide-dog/test

max_forks_repo_name: hide-dog/test


text:
__precompile__()
module AppleAccelerate

const libacc = "/System/Library/Frameworks/Accelerate.framework/Accelerate"

if !isfile(libacc)
    error("Accelerate framework not found at $(libacc)")
end

include("Array.jl")
include("DSP.jl")
include("Util.jl")

end # module


text:
using Revise
using BridgeLandmarks
const BL = BridgeLandmarks
using Distributions
using LinearAlgebra
using Random
# using FileIO
# using JLD2
using DelimitedFiles

Random.seed!(9)

################ set directories ##########################################
workdir = @__DIR__
cd(workdir)
include("outdirpath.jl")

dir1 = joinpath(outdirpath,"millerjoshi/ms")
outdir_ms = mkpath(dir1)

dir2 = joinpath(outdirpath,"millerjoshi/ahs")
outdir_ahs = mkpath(dir2)

################ read data ##########################################
workdir = @__DIR__
cd(workdir)
d1 = readdlm("revised_experiments/millerjoshi/dataset1.csv",',')
d2 = readdlm("revised_experiments/millerjoshi/dataset2.csv",',')
n = 12
# center and scale
center = PointF(40.0, 20.0)
xobs0 = map(i->PointF(d1[i,:]) - center, 1:n)/10.0
xobsT = map(i->PointF(d2[i,:]) - center, 1:n)/10.0

################ settings and mcmc #################################
ups = [:innov, :mala_mom, :parameter]
skip_saveITER = 1000
printskip = 1000
adaptskip = ITER = 20_000

p_ms = Pars_ms(γinit=1.0/√n, aprior=Pareto(1.0, 0.1), dt = 0.01,
                adaptskip=adaptskip, skip_saveITER=skip_saveITER,
                ρlowerbound=0.9)

@time landmarksmatching(xobs0,xobsT; ITER=ITER, pars=p_ms, updatescheme=ups, printskip=printskip, outdir=outdir_ms)

p_ahs = Pars_ahs(db=[3.0, 2.0], stdev=.5, γinit=.1, aprior=Pareto(1.0, 0.1), dt = 0.01,
                                adaptskip=adaptskip, skip_saveITER=skip_saveITER,
                                ρlowerbound=0.9)
@time landmarksmatching(xobs0,xobsT; ITER=ITER, pars=p_ahs, updatescheme=ups, printskip=printskip, outdir=outdir_ahs)

# using ProfileView
# using Profile
#Profile.init() @profview


text:
using GBIF
using Test

global anyerrors = false

tests = [
  "species retrieval" => "taxon.jl",
  "single occurrence functions" => "occurrence.jl",
  "multiple occurrences" => "occurrences.jl",
  "paging" => "paging.jl",
  "iteration" => "iteration.jl",
  "methods" => "methods.jl",
  "data frame support" => "dataframes.jl",
  "query support" => "query.jl"
]

for test in tests
  try
    include(test.second)
    println("\033[1m\033[32m✓\033[0m\t$(test.first)")
  catch e
    global anyerrors = true
    println("\033[1m\033[31m×\033[0m\t$(test.first)")
    println("\033[1m\033[38m→\033[0m\ttest/$(test.second)")
    showerror(stdout, e, backtrace())
    println()
    break
  end
end

if anyerrors
  throw("Tests failed")
end


text:
"""
     enkf(state_ens, d_matrix, hx_matrix)

Return updated states
"""
function enkf(state_ens, d_matrix, hx_matrix)

  # Variables (Mandel)

  X  = state_ens;
  D  = d_matrix;
  HX = hx_matrix;

  # Subtract ensemble mean (Mandel)

  (n, N) = size(X);
  (m, N) = size(D);

  A    = X - 1 / N * (X * ones(Float64, N, 1)) * ones(Float64, 1, N);
  HA   = HX - 1 / N * (HX * ones(Float64, N, 1)) * ones(Float64, 1, N);
  Dtmp = D - 1 / N * (D * ones(Float64, N, 1)) * ones(Float64, 1, N);

  # Observation error variance (Mandel-theoretic, Evensen-sample)

  R_sample = Dtmp * Dtmp' / (N-1);

  # Variance of predicted observations (DeChant and Mandel)

  C_YY = 1 / (N-1) * HA * HA';

  # Covariance between states ensemble and predicted observations (DeChant and Mandel)

  C_XY = 1 / (N-1) * A * HA';

  # Compute kalman gain (DeChant)

  K = C_XY / (C_YY + R_sample);

  if any(isnan(K))
    println("R_sample = $R_sample")
    println("C_YY = $C_YY")
    println("C_XY = $C_XY")
    error("Nans in Kalman gain")
  end

  # Update states (DeChant and Mandel)

  Xhat = X + K*(D-HX);

  return(Xhat);

end

"""
     enkf!(state_ens, d_matrix, hx_matrix)

Update states in place
"""
function enkf!(state_ens, d_matrix, hx_matrix)

  # Variables (Mandel)

  X  = state_ens;
  D  = d_matrix;
  HX = hx_matrix;

  # Subtract ensemble mean (Mandel)

  (n, N) = size(X);
  (m, N) = size(D);

  A    = X - 1 / N * (X * ones(Float64, N, 1)) * ones(Float64, 1, N);
  HA   = HX - 1 / N * (HX * ones(Float64, N, 1)) * ones(Float64, 1, N);
  Dtmp = D - 1 / N * (D * ones(Float64, N, 1)) * ones(Float64, 1, N);

  # Observation error variance (Mandel-theoretic, Evensen-sample)

  R_sample = Dtmp * Dtmp' / (N-1);

  # Variance of predicted observations (DeChant and Mandel)

  C_YY = 1 / (N-1) * HA * HA';

  # Covariance between states ensemble and predicted observations (DeChant and Mandel)

  C_XY = 1 / (N-1) * A * HA';

  # Compute kalman gain (DeChant)

  K = C_XY / (C_YY + R_sample);

  if any(isnan(K))
    println("R_sample = $R_sample")
    println("C_YY = $C_YY")
    println("C_XY = $C_XY")
    error("Nans in Kalman gain")
  end

  # Update states (DeChant and Mandel)

  X = X + K*(D-HX);

  nothing

end


text:
# Autogenerated wrapper script for boost_jll for powerpc64le-linux-gnu-cxx03
export libboost_atomic, libboost_chrono, libboost_container, libboost_context, libboost_contract, libboost_coroutine, libboost_date_time, libboost_filesystem, libboost_graph, libboost_iostreams, libboost_log, libboost_log_setup, libboost_prg_exec_monitor, libboost_program_options, libboost_random, libboost_regex, libboost_serialization, libboost_stacktrace_basic, libboost_stacktrace_noop, libboost_system, libboost_thread, libboost_timer, libboost_type_erasure, libboost_unit_test_framework, libboost_wave, libboost_wserialization

JLLWrappers.@generate_wrapper_header("boost")
JLLWrappers.@declare_library_product(libboost_atomic, "libboost_atomic.so.1.76.0")
JLLWrappers.@declare_library_product(libboost_chrono, "libboost_chrono.so.1.76.0")
JLLWrappers.@declare_library_product(libboost_container, "libboost_container.so.1.76.0")
JLLWrappers.@declare_library_product(libboost_context, "libboost_context.so.1.76.0")
JLLWrappers.@declare_library_product(libboost_contract, "libboost_contract.so.1.76.0")
JLLWrappers.@declare_library_product(libboost_coroutine, "libboost_coroutine.so.1.76.0")
JLLWrappers.@declare_library_product(libboost_date_time, "libboost_date_time.so.1.76.0")
JLLWrappers.@declare_library_product(libboost_filesystem, "libboost_filesystem.so.1.76.0")
JLLWrappers.@declare_library_product(libboost_graph, "libboost_graph.so.1.76.0")
JLLWrappers.@declare_library_product(libboost_iostreams, "libboost_iostreams.so.1.76.0")
JLLWrappers.@declare_library_product(libboost_log, "libboost_log.so.1.76.0")
JLLWrappers.@declare_library_product(libboost_log_setup, "libboost_log_setup.so.1.76.0")
JLLWrappers.@declare_library_product(libboost_prg_exec_monitor, "libboost_prg_exec_monitor.so.1.76.0")
JLLWrappers.@declare_library_product(libboost_program_options, "libboost_program_options.so.1.76.0")
JLLWrappers.@declare_library_product(libboost_random, "libboost_random.so.1.76.0")
JLLWrappers.@declare_library_product(libboost_regex, "libboost_regex.so.1.76.0")
JLLWrappers.@declare_library_product(libboost_serialization, "libboost_serialization.so.1.76.0")
JLLWrappers.@declare_library_product(libboost_stacktrace_basic, "libboost_stacktrace_basic.so.1.76.0")
JLLWrappers.@declare_library_product(libboost_stacktrace_noop, "libboost_stacktrace_noop.so.1.76.0")
JLLWrappers.@declare_library_product(libboost_system, "libboost_system.so.1.76.0")
JLLWrappers.@declare_library_product(libboost_thread, "libboost_thread.so.1.76.0")
JLLWrappers.@declare_library_product(libboost_timer, "libboost_timer.so.1.76.0")
JLLWrappers.@declare_library_product(libboost_type_erasure, "libboost_type_erasure.so.1.76.0")
JLLWrappers.@declare_library_product(libboost_unit_test_framework, "libboost_unit_test_framework.so.1.76.0")
JLLWrappers.@declare_library_product(libboost_wave, "libboost_wave.so.1.76.0")
JLLWrappers.@declare_library_product(libboost_wserialization, "libboost_wserialization.so.1.76.0")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libboost_atomic,
        "lib/libboost_atomic.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libboost_chrono,
        "lib/libboost_chrono.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libboost_container,
        "lib/libboost_container.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libboost_context,
        "lib/libboost_context.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libboost_contract,
        "lib/libboost_contract.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libboost_coroutine,
        "lib/libboost_coroutine.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libboost_date_time,
        "lib/libboost_date_time.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libboost_filesystem,
        "lib/libboost_filesystem.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libboost_graph,
        "lib/libboost_graph.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libboost_iostreams,
        "lib/libboost_iostreams.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libboost_log,
        "lib/libboost_log.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libboost_log_setup,
        "lib/libboost_log_setup.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libboost_prg_exec_monitor,
        "lib/libboost_prg_exec_monitor.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libboost_program_options,
        "lib/libboost_program_options.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libboost_random,
        "lib/libboost_random.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libboost_regex,
        "lib/libboost_regex.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libboost_serialization,
        "lib/libboost_serialization.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libboost_stacktrace_basic,
        "lib/libboost_stacktrace_basic.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libboost_stacktrace_noop,
        "lib/libboost_stacktrace_noop.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libboost_system,
        "lib/libboost_system.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libboost_thread,
        "lib/libboost_thread.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libboost_timer,
        "lib/libboost_timer.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libboost_type_erasure,
        "lib/libboost_type_erasure.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libboost_unit_test_framework,
        "lib/libboost_unit_test_framework.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libboost_wave,
        "lib/libboost_wave.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libboost_wserialization,
        "lib/libboost_wserialization.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function tcm_createRandomDistinctAssortments(M, N)
	temp = randperm(2^(N-1))[1:M] - 1

	assortments = zeros(M, N)
	for m in 1:M
		# @show bin(temp[m],N-1)
		assortments[m,1:(N-1)] = map(x -> parse(string(x)), [bin(temp[m],N-1)...] )
	end

	assortments[:,N] = 1

	return assortments
end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
#= src/SparseFFT.jl

References:
  F. Woolfe, E. Liberty, V. Rokhlin, M. Tygert. A fast randomized algorithm for
    the approximation of matrices. Appl. Comput. Harmon. Anal. 25: 335-366,
    2008.
=#

module SparseFFT

using FFTW:
  fftwReal, fftwComplex, fftwNumber, FORWARD, BACKWARD,
  FFTWPlan, plan_fft!, plan_bfft!, plan_rfft, plan_brfft

using LinearAlgebra: mul!, transpose!

export
  cSpFFTPlan1, rSpFFTPlan1, brSpFFTPlan1,
  cSpFFTPlan2, rSpFFTPlan2, brSpFFTPlan2,
  plan_spfft, spfft_f2s!, spfft_s2f!,
  plan_spbfft, spbfft_f2s!, spbfft_s2f!,
  plan_sprfft, sprfft_f2s!,
  plan_spbrfft, spbrfft_s2f!

# common

const SpFFTReal    = fftwReal
const SpFFTComplex = fftwComplex
const SpFFTNumber  = fftwNumber
const FFTPlan      = FFTWPlan

# for complex-to-real transpose (copy from master 2421ff2)
const transposebaselength=64
function transpose_f!(f,B::AbstractMatrix,A::AbstractMatrix)
    inds = axes(A)
    axes(B,1) == inds[2] && axes(B,2) == inds[1] || throw(DimensionMismatch(string(f)))

    m, n = length(inds[1]), length(inds[2])
    if m*n<=4*transposebaselength
        @inbounds begin
            for j = inds[2]
                for i = inds[1]
                    B[j,i] = f(A[i,j])
                end
            end
        end
    else
        transposeblock!(f,B,A,m,n,first(inds[1])-1,first(inds[2])-1)
    end
    return B
end
function transposeblock!(f,B::AbstractMatrix,A::AbstractMatrix,m::Int,n::Int,offseti::Int,offsetj::Int)
    if m*n<=transposebaselength
        @inbounds begin
            for j = offsetj .+ (1:n)
                for i = offseti .+ (1:m)
                    B[j,i] = f(A[i,j])
                end
            end
        end
    elseif m>n
        newm=m>>1
        transposeblock!(f,B,A,newm,n,offseti,offsetj)
        transposeblock!(f,B,A,m-newm,n,offseti+newm,offsetj)
    else
        newn=n>>1
        transposeblock!(f,B,A,m,newn,offseti,offsetj)
        transposeblock!(f,B,A,m,n-newn,offseti,offsetj+newn)
    end
    return B
end

# source files

include("spfft1.jl")
include("spfft2.jl")

end  # module


text:
InverseFunctions.inverse(::typeof(normcdf)) = norminvcdf
InverseFunctions.inverse(::typeof(norminvcdf)) = normcdf

InverseFunctions.inverse(::typeof(normccdf)) = norminvccdf
InverseFunctions.inverse(::typeof(norminvccdf)) = normccdf

InverseFunctions.inverse(::typeof(normlogcdf)) = norminvlogcdf
InverseFunctions.inverse(::typeof(norminvlogcdf)) = normlogcdf

InverseFunctions.inverse(::typeof(normlogccdf)) = norminvlogccdf
InverseFunctions.inverse(::typeof(norminvlogccdf)) = normlogccdf


text:
function preinstantiate(path=".")
  benchmark_dirs = ["Cartesian_vs_curved", "Gauss", "inlining", "numerical_fluxes",
                    "overintegration", "precompute_logs", "primitive_variables"]
  current_dir = pwd()

  for dir in benchmark_dirs
    cd(joinpath(path, dir))
    println("#"^100)
    println("### Instantiating $dir... ###")
    run(`julia -e 'import Pkg; Pkg.activate("."); Pkg.instantiate()'`)
    cd(current_dir)
  end
  println("#"^100)

  cd(current_dir)
end


text:
apply_covariance = (x,M) -> sum(dot.(x,M*x))

optional_posZ(x::NTuple{1,T}, offset::NTuple{1,T}) where {T,N} = 1
optional_posZ(x::NTuple{2,T}, offset::NTuple{2,T}) where {T,N} = 1
optional_posZ(x::NTuple{N,T}, offset::NTuple{N,T}) where {T,N} = x[3]-offset[3]

# List of functions and names we want to predefine
function generate_functions_expr()
    # offset and scale is already wrapped in the generator function
    x_expr = :(scale .* (x .- offset))
    x_expr1 = :(scale[1] .* (x[1] .- offset[1]))
    x_expr2 = :(scale[2] .* (x[2] .- offset[2]))
    x_expr3 = :(scale[3] .* (x[3] .- offset[3]))
    x_expr4 = :(scale[4] .* (x[4] .- offset[4]))
    x_expr5 = :(scale[5] .* (x[5] .- offset[5]))
    x_expr14 = :(maximum(abs.(scale .* (x .- offset)))) # useful for rectangular things
    x_expr15 = :(minimum(abs.(scale .* (x .- offset))))
    x_expr6 = :(prod(x .== offset))
    x_expr7 = :(cis(dot((x .- offset), scale)))
    x_expr8 = :(exp(.- sum(abs2.(x .- offset).*scale))) # scale is 1/(2 sigma^2)
    x_expr9 = :(exp(.- sum(abs2.(x .- offset).*scale)) ./ prod(sqrt.(pi ./ scale)))
    x_expr10 = :(exp(.- apply_covariance((x .- offset), scale) ))
    x_expr11 = :(exp(.- apply_covariance((x .- offset), scale)) ./ prod(sqrt.(pi ./ abs2(scale))))
    x_expr12 = :(abs2(scale[1] .* (x[1] .- offset[1]))+abs2(scale[2] .* (x[2] .- offset[2])))
    x_expr13 = :((x[1] .- offset[1]).* scale[1] .+ (x[2] .- offset[2]) .* scale[2])

    functions = [
        (:(rr2), :(x -> T(sum(abs2.($x_expr))))),
        (:(rr),  :(x -> T(sqrt.(sum(abs2.($x_expr)))))),
        (:(xx),  :(x -> T($x_expr1))),
        (:(yy),  :(x -> T($x_expr2))),
        (:(zz),  :(x -> T($x_expr3))),
        (:(ee),  :(x -> T($x_expr4))),
        (:(tt),  :(x -> T($x_expr5))),
        (:(idx_min),  :(x -> T($x_expr15))),
        (:(idx_max),  :(x -> T($x_expr14))),
        (:(delta),  :(x -> T($x_expr6))),
        (:(phiphi), :(x -> T(atan.($x_expr2, $x_expr1)))),  # this is the arcus tangens of y/x yielding a spiral phase ramp
        (:(phase_kz), :(x -> T(optional_posZ(x,offset).*sqrt.(max(1 .- $x_expr12,0))))),  # can be used for constucting a free-space propagator in optics
        (:(phase_kxy), :(x -> T($x_expr13))),  # useful for xy shifting
        (:(exp_is),  :(x -> T($x_expr7))),  # exp(2pi i s (x-o)) # by modifying s, this becomes exp(i kx)
        (:(exp_sqr),  :(x -> T($x_expr8))),  # maximum-normalized Gaussian
        (:(exp_sqr_norm),  :(x -> T($x_expr9))),  # integral-normalized Gaussian (over infinite ROI)
        (:(exp_sqr_cov),  :(x -> T($x_expr10))),  # Gaussian based on the covariance matrix
        (:(exp_sqr_cov_norm),  :(x -> T($x_expr11))),  # normalized Gaussian based on the covariance matrix
    ]
    return functions
end

function wrap_fkt(f)  # encapsulates the expression in a proper function with index, offset and scale argument
    return (idx, offset, scale, weight) -> f(idx)  # Not that f contains offset and scale
end

function wrap_zipped(g) # encapsulates a wrapped function as called with a sinle tuple of offset and scale
    return (idx, mytuple) -> g(idx,mytuple...)
end

curry(f, x) = (xs...) -> f(x, xs...)   # just a shorthand to remove x

## These functions ensure that also numbers can be iterated and zipped
function cast_iter(vals::Matrix)
    Tuple(Tuple(vals[:,n]) for n in 1:size(vals,2))
    # Tuple(vals[:,axes(vals, 2)])
end

function cast_iter(vals::Vector{<:Number})
    Tuple(vals)   # makes the matrix iterable, interpreting it as a series of vectors
end

function cast_iter(vals::IterType)
    vals
end

function cast_iter(vals)
    repeated(vals) # during the zip operation this type always yields the same value
end

function cast_number_iter(vals::Vector{<:Number})
    Tuple(vals)   # makes the matrix iterable, interpreting it as a series of vectors
end

function cast_number_iter(vals::NTuple{N,<:Number} where N)
    vals
end

function cast_number_iter(vals::Number)
    repeated(vals) # during the zip operation this type always yields the same value
end

function optional_mat_to_iter(vals)  # only for matrices
    vals
end

function optional_mat_to_iter(vals::Matrix)
    cast_iter(vals)
end

mat_to_tvec(v) = [Tuple(v[:,n]) for n in 1:size(v,2)] # converts a 2d matrix to a Vector of Tuples

function apply_dims(scale2, dims, N)  # replaces scale entries not in dims with zeros
    ntuple(i -> i ∈ dims ? scale2[i] : zero(eltype(scale2)), N)  # zero(scale2[1])
end

function apply_dims(scales::IterType, dims, N)  # replaces scale entries not in dims with zeros
    Tuple(ntuple(i -> i ∈ dims ? scale3[i] : zero(eltype(scale3)), N)  for scale3 in  scales)
end

# we automatically generate the functions for rr2, rr, ...
# We set the types for the arguments correctly in the default cases
for F in generate_functions_expr()

    # default functions with offset and scaling behavior. This version allows no list of numbers or tuples
    @eval function $(Symbol(:_, F[1]))(::Type{T}, size::NTuple{N, Int},
                           offset::Union{Type{<:Ctr}, Number, NTuple{N,Number}, Vector{<:Number}},
                           scale::Union{Type{<:Sca}, Number, NTuple{N,Number}, Vector{<:Number}},
                           dims,
                           accumulator,
                           weight::Number) where{N, T}
        offset = get_offset(size, offset)
        scale_init = get_scale(size, scale)
        scale = apply_dims(scale_init, dims, N) # replaces scale entries not in dims with zeros

        g(x) = weight .* ($(F[2])(x)) #
        IndexFunArray(T,g, size)
    end

    @eval function $(Symbol(:_, F[1]))(::Type{T}, size::NTuple{N, Int},
        offset,   # a version that supports an iterable collection of tuples for offset and scale
        scale,
        dims=ntuple(+, N),
        accumulator = sum,
        weight=1) where{N, T}

        # its important to use different variable names due to a julia bug
        offset_ = get_offset(size, offset)
        scale_init = get_scale(size, scale)
        scale_ = apply_dims(scale_init, dims, N)

        offset2 = cast_iter(offset_)
        # @show offset2
        scale2 = cast_iter(scale_)
        weight2 = cast_number_iter(weight)
        w = (x, offset, scale, weight3) -> weight3 .* $(F[2])(x) # adds offset and scale as requires parameters to the expression
        # the line below makes a function only depending on the index position but iterating over both, offsets and scales
        fkt = (idx) -> accumulator(curry(wrap_zipped(w),idx),zip(offset2,scale2,weight2))

        IndexFunArray(T, fkt, size)
    end
    # default functions with offset and scaling behavior
    @eval function $(F[1])(::Type{T}, size::NTuple{N, Int};
                           offset=CtrFT,
                           scale=ScaUnit,
                           dims=ntuple(+, N), accumulator=sum, weight=1) where{N, T}
        # $(F[1])(T, size, offset, scale, dims)
        # print("Dispatch T3\n")
        $(Symbol(:_, F[1]))(T, size, offset, scale, dims, accumulator, weight)
    end

    # change order of offset and scale
    @eval function $(F[1])(size::NTuple{N, Int};
                           offset=CtrFT, scale=ScaUnit, dims=ntuple(+, N), accumulator=sum, weight=1) where {N}
        T = $DEFAULT_T
        $(F[1])(T, size, scale=scale, offset=offset, dims=dims, accumulator=accumulator, weight=weight)
    end

    # convenient wrapper to provide an array as input
    @eval function $(F[1])(arr::AbstractArray{T, N};
                           offset=CtrFT, scale=ScaUnit, dims=ntuple(+, N), accumulator=sum, weight=1) where {T, N}
        T2 = default_type(T, $DEFAULT_T)
        $(F[1])(T2, size(arr), scale=scale, offset=offset, dims=dims, accumulator=accumulator, weight=weight)
    end

    @eval export $(F[1])
end

export axes1d

function axes1d(::Type{T}, size::NTuple{N, Int}; offset=CtrFT, scale=ScaUnit, dims=ntuple(+, N), keepdims=false) where {T,N}
    offset = get_offset(size, offset)
    scale_init = get_scale(size, scale)
    scale = apply_dims(scale_init, dims, N) # replaces scale entries not in dims with zeros
    if keepdims
        @show offset
        (ramp(T, d, size[d]; offset=offset[d],scale=scale[d]) for d in dims)
    else
        (xx(T, (size[d],); offset=offset[d], scale=scale[d]) for d in dims)
    end
end

function axes1d(size::NTuple{N, Int}; offset=CtrFT, scale=ScaUnit, dims=ntuple(+, N), keepdims=false) where {N}
    T = DEFAULT_T
    axes1d(T, size; offset=offset, scale=scale, dims=dims, keepdims=keepdims)
end

function axes1d(arr::AbstractArray{T, N}; offset=CtrFT, scale=ScaUnit, dims=ntuple(+, N), keepdims=false) where {T, N}
    axes1d(size(arr); offset=offset, scale=scale, dims=dims, keepdims=keepdims)
end