text:
dottosym(x)::Symbol = x

function dottosym(x::Expr)::Symbol
  s1 = dottosym(x.args[1])
  xa2 = x.args[2]
  xa2 isa QuoteNode ? Symbol(s1, "###extractarray###", xa2.value) : Symbol(s1, "###extractarray###", xa2)
end
function extract_array_symbol_from_ref!(ls::LoopSet, ex::Expr, offset1::Int)::Symbol
  ar = ex.args[1 + offset1]
  if isa(ar, Symbol)
    return ar
  elseif isa(ar, Expr) && (ar.head === :(.) || ar.head === :ref)
    s = dottosym(ar)
    # pushpreamble!(ls, Expr(:(=), s, ar))
    pushprepreamble!(ls, Expr(:(=), s, ar))
    return s
  else
    throw("Indexing into the following expression was not recognized: $ar")
  end
end

function ref_from_expr!(ls, ex, offset1::Int, offset2::Int)
    ar = extract_array_symbol_from_ref!(ls, ex, offset1)
    ar, @view(ex.args[2 + offset2:end])
end
ref_from_ref!(ls::LoopSet, ex::Expr) = ref_from_expr!(ls, ex, 0, 0)
ref_from_getindex!(ls::LoopSet, ex::Expr) = ref_from_expr!(ls, ex, 1, 1)
ref_from_setindex!(ls::LoopSet, ex::Expr) = ref_from_expr!(ls, ex, 1, 2)
function ref_from_expr!(ls::LoopSet, ex::Expr)
  if ex.head === :ref
    ref_from_ref!(ls, ex)
  else#if ex.head === :call
    f = first(ex.args)::Symbol
    f === :getindex ? ref_from_getindex!(ls, ex) : ref_from_setindex!(ls, ex)
  end
end

add_vptr!(ls::LoopSet, op::Operation) = add_vptr!(ls, op.ref)
add_vptr!(ls::LoopSet, mref::ArrayReferenceMeta) = add_vptr!(ls, mref.ref.array, vptr(mref))
# using VectorizationBase: noaliasstridedpointer
presbufsym(array) = Symbol('#', array, "#preserve#buffer#")
function add_vptr!(ls::LoopSet, array::Symbol, vptrarray::Symbol, actualarray::Bool = true)
  if !includesarray(ls, array)
    push!(ls.includedarrays, array)
    actualarray && push!(ls.includedactualarrays, array)
    pushpreamble!(ls, Expr(:(=), Expr(:tuple, vptrarray, presbufsym(array)), Expr(:call, lv(:stridedpointer_preserve), array)))
  end
  nothing
end

@inline staticdims(::Any) = One()
@inline staticdims(::CartesianIndices{N}) where {N} = StaticInt{N}()

function append_loop_staticdims!(valcall::Expr, loop::Loop, constoffset::Int)
  isstaticloop(loop) && return constoffset + 1
  push!(valcall.args, Expr(:call, lv(:staticdims), loop.rangesym))
  constoffset
end
function subset_vptr!(ls::LoopSet, vptr::Symbol, indnum::Int, ind, previndices, loopindex, D::Int)
  subset = D == 0
  str_typ = subset ? "subset" : "index"
  subsetvptr = Symbol(vptr, "_##$(str_typ)##_$(indnum)_##with##_$(ind)_##")
  valcall = Expr(:call, +)
  constoffset = 1
  if indnum > 1
    offset = first(previndices) === DISCONTIGUOUS
    # valcall =  Expr(:call, :(+), valcall)
    for i ∈ 1:indnum-1
      if loopindex[i]
        loopdep = previndices[i+offset]
        if loopdep === CONSTANTZEROINDEX
          constoffset += 1
          continue
        end
      else
        # assumes all staticdims will be of equal length once expanded...
        # A[I + J, constindex], I and J may be CartesianIndices. This requires they all be of same number of dims
        loopdep = first(loopdependencies(ls.opdict[previndices[i+offset]]))
      end
      constoffset = append_loop_staticdims!(valcall, getloop(ls, loopdep), constoffset)
    end
  end
  # indm1 = ind isa Integer ? ind - 1 : Expr(:call, :-, ind, 1)
  f = lv(Core.ifelse(subset, :subsetview, :_gesp))
  constoffsetexpr = staticexpr(constoffset)
  if length(valcall.args) ≠ 1
    push!(valcall.args, constoffsetexpr)
    constoffsetexpr = valcall
  end
  fexpr = Expr(:call, f, vptr, constoffsetexpr, ind)
  subset || push!(fexpr.args, staticexpr(D))
  pushpreamble!(ls, Expr(:(=), subsetvptr, fexpr))
  subsetvptr
end

function gesp_const_offset!(ls::LoopSet, vptrarray::Symbol, ninds::Int, indices::Vector{Symbol}, loopedindex::Vector{Bool}, mlt::Integer, sym, D::Int)
  if isone(mlt)
    subset_vptr!(ls, vptrarray, ninds, sym, indices, loopedindex, D)
  else
    # want same name for arrays to be given the same name if possible
    mltsym = Symbol(sym, "##multiplied##by##", mlt)
    pushpreamble!(ls, Expr(:(=), mltsym, Expr(:call, :(*), mlt, sym)))
    # pushprepreamble!(ls, Expr(:(=), mltsym, Expr(:call, :(*), mlt, sym)))
    subset_vptr!(ls, vptrarray, ninds, mltsym, indices, loopedindex, D)
  end
end
function gesp_const_offsets!(ls::LoopSet, vptrarray::Symbol, ninds::Int, indices::Vector{Symbol}, loopedindex::Vector{Bool}, mltsyms::Vector{Tuple{Int,Symbol}}, D::Int)
  length(mltsyms) > 1 && sort!(mltsyms, by = last) # if multiple have same combination of syms, make sure they match even if order is different
  for (mlt,sym) ∈ mltsyms
    vptrarray = gesp_const_offset!(ls, vptrarray, ninds, indices, loopedindex, mlt, sym, D)
  end
  vptrarray
end

byterepresentable(x)::Bool = false
byterepresentable(x::Integer)::Bool = typemin(Int8) ≤ x ≤ typemax(Int8)
function _addoffset!(indices, offsets, strides, loopedindex, loopdependencies, ind, offset, stride)
    push!(indices, ind)
    push!(offsets, offset % Int8)
    # push!(offsets, (offset+stride-1) % Int8)
    push!(strides, stride % Int8)
    push!(loopedindex, true)
    push!(loopdependencies, ind)
    true
end

function addconstindex!(indices, offsets, strides, loopedindex, offset)
    push!(indices, CONSTANTZEROINDEX)
    # push!(offsets, (offset-1) % Int8)
    push!(offsets, (offset) % Int8)
    push!(strides, zero(Int8))
    push!(loopedindex, true)
    true
end
function addopindex!(
    parents::Vector{Operation}, loopdependencies::Vector{Symbol}, reduceddeps::Vector{Symbol},
    indices::Vector{Symbol}, offsets::Vector{Int8}, strides::Vector{Int8}, loopedindex::Vector{Bool}, indop::Operation,
    stride = one(Int8), offset = zero(Int8)
)
    pushparent!(parents, loopdependencies, reduceddeps, indop)
    push!(indices, name(indop));
    push!(offsets, offset)
    push!(strides, stride)
    push!(loopedindex, false)
    nothing
end

####################################################################################################
###################################### Index Parsing Helpers #######################################
####################################################################################################

function add_affine_index_expr!(ls::LoopSet, mult_syms::Vector{Tuple{Int,Symbol}}, constant::Base.RefValue{Int}, stride::Int, expr::Symbol)
    push!(mult_syms, (stride, expr))
    return nothing
end
function add_affine_index_expr!(ls::LoopSet, mult_syms::Vector{Tuple{Int,Symbol}}, constant::Base.RefValue{Int}, stride::Int, expr::Integer)
    constant[] += stride*expr
    return nothing
end
function add_affine_op!(ls::LoopSet, mult_syms::Vector{Tuple{Int,Symbol}}, constant::Base.RefValue{Int}, stride::Int, expr::Expr)
    parent = add_operation!(ls, gensym!(ls, "indexpr"), expr, sizeof(Int), length(ls.loopsymbols))
    add_affine_index_expr!(ls, mult_syms, constant, stride, name(parent))
    return nothing
end
function add_mul!(ls::LoopSet, mult_syms::Vector{Tuple{Int,Symbol}}, constant::Base.RefValue{Int}, stride::Int, arg1, arg2, expr)
    if arg1 isa Integer
        add_affine_index_expr!(ls, mult_syms, constant, stride * arg1, arg2)
    elseif arg2 isa Integer
        add_affine_index_expr!(ls, mult_syms, constant, stride * arg2, arg1)
    else
        add_affine_op!(ls, mult_syms, constant, stride, expr)
    end
    return nothing
end
function add_affine_index_expr!(ls::LoopSet, mult_syms::Vector{Tuple{Int,Symbol}}, constant::Base.RefValue{Int}, stride::Int, expr::Expr)
    expr.head === :call || return add_affine_op!(ls, mult_syms, constant, stride, expr)
    f = expr.args[1]
    if f === :(*)
        @assert length(expr.args) == 3
        add_mul!(ls, mult_syms, constant, stride, expr.args[2], expr.args[3], expr)
    elseif f === :(-)
        if length(expr.args) == 3
            add_affine_index_expr!(ls, mult_syms, constant, stride, expr.args[2])
        elseif length(expr.args) ≠ 2
            throw("Subtraction expression: $expr has an invalid number of arguments.")
        end
        add_affine_index_expr!(ls, mult_syms, constant, -stride, last(expr.args))
    elseif f === :(+)
        for i ∈ 2:length(expr.args)
            add_affine_index_expr!(ls, mult_syms, constant, stride, expr.args[i])
        end
    else
        add_affine_op!(ls, mult_syms, constant, stride, expr)
    end
    return nothing
end
function affine_index_expression(ls::LoopSet, expr)::Tuple{Int,Vector{Tuple{Int,Symbol}}}
    mult_syms = Tuple{Int,Symbol}[]
    constant = Ref(0)
    add_affine_index_expr!(ls, mult_syms, constant, 1, expr)
    return constant[], mult_syms
end

function muladd_index!(
    ls::LoopSet, parents, loopdependencies, reduceddeps, indices, offsets, strides, loopedindex,
    mlt::Int, sym::Symbol, offset::Int
)
    muladd_index!(ls, parents, loopdependencies, reduceddeps, indices, offsets, strides, loopedindex, mlt, getop(ls, sym, sizeof(Int)), offset)
end
function muladd_op!(ls::LoopSet, mlt::Int, sym::Symbol, offset::Int)
    muladd_op!(ls, mlt, getop(ls, sym, sizeof(Int)), offset)
end
function muladd_op!(ls::LoopSet, mlt::Int, symop::Operation, offset::Int)
    indsym = gensym!(ls, "affindexop")
    vparents = [symop]
    # vparents = [, mltop, offop]
    if mlt == -1
        f = :vsub_nsw
        if offset ≠ 0
            pushfirst!(vparents, add_constant!(ls, offset, sizeof(Int)))
        end
    elseif mlt == 1
        offset == 0 && return only(vparents)
        push!(vparents, add_constant!(ls, offset, sizeof(Int)))
        f = :vadd_nsw
    else
        push!(vparents, add_constant!(ls, mlt, sizeof(Int)))
        f = if offset == 0
            :vmul_nsw
        else
            push!(vparents, add_constant!(ls, offset, sizeof(Int)))
            :muladd
        end
    end
    add_compute!(ls, indsym, instruction(f), vparents, sizeof(Int))
end
function muladd_index!(
    ls::LoopSet, parents, loopdependencies, reduceddeps, indices, offsets, strides, loopedindex,
    mlt::Int, symop::Operation, offset::Int
)
    if byterepresentable(offset) & byterepresentable(mlt)
        addopindex!(
            parents, loopdependencies, reduceddeps, indices,
            offsets, strides, loopedindex, symop, mlt, offset
        )
    else
        indop = muladd_op!(ls, mlt, symop, offset)
        addopindex!(parents, loopdependencies, reduceddeps, indices, offsets, strides, loopedindex, indop)
    end
end

function getintconstant(ls::LoopSet, op::Operation)
  for (i, (v,bits,bool)) ∈ ls.preamble_symint
    i == identifier(op) || continue
    return byterepresentable(v), v%Int8
  end
  return false, zero(Int8)
end

function add_additive_index!(ls::LoopSet, opparents, vptrarray, sym, mop, ninds, indices, offsets, strides, loopedindex, loopdependencies, reduceddeps, offset, mlt, D)
  factor = Core.ifelse((instruction(mop).instr === :sub_fast), -1, 1)
  if length(parents(mop)) == 2
    sub1 = parents(mop)[1]
    sub2 = parents(mop)[2]
    if isloopvalue(sub1) & isconstant(sub2)
      isliteralint, literalval = getintconstant(ls, sub2)
      if isliteralint
        _addoffset!(indices, offsets, strides, loopedindex, loopdependencies, name(sub1), offset+literalval, mlt)
      else
        vptrarray = gesp_const_offset!(ls, vptrarray, ninds, indices, loopedindex, factor*mlt, name(sub2), D)
        _addoffset!(indices, offsets, strides, loopedindex, loopdependencies, name(sub1), offset, mlt)
      end
    elseif isloopvalue(sub2) & isconstant(sub1)
      isliteralint, literalval = getintconstant(ls, sub1)
      if isliteralint
        _addoffset!(indices, offsets, strides, loopedindex, loopdependencies, name(sub2), offset+literalval, factor*mlt)
      else
        vptrarray = gesp_const_offset!(ls, vptrarray, ninds, indices, loopedindex, mlt, name(sub1), D)
        _addoffset!(indices, offsets, strides, loopedindex, loopdependencies, name(sub2), offset, factor*mlt)
      end
    else
      muladd_index!(ls, opparents, loopdependencies, reduceddeps, indices, offsets, strides, loopedindex, mlt, sym, offset)
    end
  else
    muladd_index!(ls, opparents, loopdependencies, reduceddeps, indices, offsets, strides, loopedindex, mlt, sym, offset)
  end
  vptrarray
end

function checkforoffset!(
    ls::LoopSet, vptrarray::Symbol, ninds::Int, opparents::Vector{Operation}, indices::Vector{Symbol}, offsets::Vector{Int8}, strides::Vector{Int8},
    loopedindex::Vector{Bool}, loopdependencies::Vector{Symbol}, reduceddeps::Vector{Symbol}, ind::Expr, D::Int
)::Symbol

  offset, mult_syms = affine_index_expression(ls, ind)
  let deleted = 0, N = length(mult_syms)
    for n ∈ 1:N
      ntemp = n - deleted
      mlt, sym = mult_syms[ntemp]
      opm = get(ls.opdict, sym, nothing)
      opm === nothing && continue
      isconstant(opm) || continue
      found = false
      for (opid,(intval,intsz,signed)) ∈ ls.preamble_symint
        if opid == identifier(opm)
          offset += intval * mlt
          deleted += 1
          deleteat!(mult_syms, ntemp)
          found = true
          break
        end
      end
      found && continue
      for (opid,nt) ∈ ls.preamble_zeros
        if opid == identifier(opm)
          deleted += 1
          deleteat!(mult_syms, ntemp)
          break
        end
      end
    end
  end
  if !byterepresentable(offset)
    if length(mult_syms) == 1
      mlt,sym = only(mult_syms)
      if !byterepresentable(mlt)
        # this is so we don't unnecessarilly add a separate offset
        muladd_index!(ls, opparents, loopdependencies, reduceddeps, indices, offsets, strides, loopedindex, mlt, sym, offset)
        return vptrarray
      end
    end
    r = copysign(abs(offset) & 127, offset)
    vptrarray = gesp_const_offset!(ls, vptrarray, ninds, indices, loopedindex, 1, offset - r, D)
    offset = r
  end
  # (success && byterepresentable(offset)) || return false, vptrarray
  if length(mult_syms) == 0
    addconstindex!(indices, offsets, strides, loopedindex, offset)
    return vptrarray
  elseif length(mult_syms) == 1
    mlt, sym = mult_syms[1]
    if sym ∈ ls.loopsymbols
      if byterepresentable(mlt)
        _addoffset!(indices, offsets, strides, loopedindex, loopdependencies, sym, offset, mlt)
      else
        muladd_index!(ls, opparents, loopdependencies, reduceddeps, indices, offsets, strides, loopedindex, mlt, sym, offset)
      end
    elseif !byterepresentable(mlt)
      muladd_index!(ls, opparents, loopdependencies, reduceddeps, indices, offsets, strides, loopedindex, mlt, sym, offset)
    else
      mop = get(ls.opdict, sym, nothing)
      if mop === nothing
        vptrarray = gesp_const_offset!(ls, vptrarray, ninds, indices, loopedindex, mlt, sym, D)
        addconstindex!(indices, offsets, strides, loopedindex, offset)
      elseif isconstant(mop)
        vptrarray = gesp_const_offset!(ls, vptrarray, ninds, indices, loopedindex, mlt, name(mop), D)
        addconstindex!(indices, offsets, strides, loopedindex, offset)
      elseif (instruction(mop).instr === :add_fast) || (instruction(mop).instr === :sub_fast)
        vptrarray = add_additive_index!(ls, opparents, vptrarray, sym, mop, ninds, indices, offsets, strides, loopedindex, loopdependencies, reduceddeps, offset, mlt, D)
      else
        muladd_index!(ls, opparents, loopdependencies, reduceddeps, indices, offsets, strides, loopedindex, mlt, sym, offset)
      end
    end
    return vptrarray
  end
  loopsym_ind = 0
  operations = Operation[]
  operation_mults = Int[]
  deleteat_inds = Int[]
  for (i,(m,s)) ∈ enumerate(mult_syms)
    if s ∈ ls.loopsymbols
      if loopsym_ind == 0
        loopsym_ind = i
        push!(deleteat_inds, loopsym_ind)
      else
        push!(operations, add_loopvalue!(ls, s, sizeof(Int)))
        push!(deleteat_inds, i)
        push!(operation_mults, m)
      end
    else
      sop = get(ls.opdict, s, nothing)
      if sop === nothing
        push!(deleteat_inds, i)
        vptrarray = gesp_const_offset!(ls, vptrarray, ninds, indices, loopedindex, m, s, D)
      elseif isconstant(sop)
        push!(deleteat_inds, i)
        vptrarray = gesp_const_offset!(ls, vptrarray, ninds, indices, loopedindex, m, name(sop), D)
      else
        # @show sop
        # if instruction(sop).instr === :sub_fast
        #   if length(parents(sop)) == 1
        #     sub1 = parents(sop)[1]
        #   elseif length(parents(op)) == 2
        #     sub1 = parents(sop)[1]
        #     sub2 = parents(sop)[2]
        #     if isloopvalue(sub1)
        #     elseif isloopvalue(sub2)
        #       @show sub2
        #     end
        #   end
        # elseif instruction(op).instr === :add_fast
        # end
        push!(operations, getop(ls, s, 8))
        push!(operation_mults, m)
        push!(deleteat_inds, i)
      end
    end
  end
  if (length(operations) > 0) & (loopsym_ind > 0) # turn into an operation
    _m,_s = mult_syms[loopsym_ind]
    push!(operations, add_loopvalue!(ls, _s, sizeof(Int)))
    push!(operation_mults, _m)
  end
  if length(operations) == 0
    if loopsym_ind == 0
      addconstindex!(indices, offsets, strides, loopedindex, offset)
    else
      mlt, sym = mult_syms[loopsym_ind]
      _addoffset!(indices, offsets, strides, loopedindex, loopdependencies, sym, offset, mlt)
    end
    deleteat!(mult_syms, deleteat_inds)
    return gesp_const_offsets!(ls, vptrarray, ninds, indices, loopedindex, mult_syms, D)
  end
  deleteat!(mult_syms, deleteat_inds)
  vptrarray = gesp_const_offsets!(ls, vptrarray, ninds, indices, loopedindex, mult_syms, D)
  if length(operations) == 1
    _mlt = only(operation_mults)
    indop = muladd_op!(ls, Core.ifelse(byterepresentable(_mlt), 1, _mlt), only(operations), 0)
    addopindex!(opparents, loopdependencies, reduceddeps, indices, offsets, strides, loopedindex, indop, Core.ifelse(byterepresentable(_mlt), _mlt%Int8, one(Int8)), offset%Int8)
  else
    mlt1ind = findfirst(isone, operation_mults)
    opbase = if mlt1ind === nothing # if none of them have a multiplier of `1`, pick the base arbtirarility
      muladd_op!(ls, pop!(operation_mults), pop!(operations), 0)
    else # otherwise, we start accumulating with a 1-multiplier.
      deleteat!(operation_mults, mlt1ind)
      popat!(operations, mlt1ind)
    end
    for i ∈ eachindex(operations)
      _op = operations[i]
      _mlt = operation_mults[i]
      opbase = if _mlt == -1
        add_compute!(ls, gensym!(ls, "indexaccum"), instruction(:(-)), [opbase, _op], sizeof(Int))
      elseif _mlt == 1
        add_compute!(ls, gensym!(ls, "indexaccum"), instruction(:(+)), [opbase, _op], sizeof(Int))
      else
        add_compute!(ls, gensym!(ls, "indexaccum"), instruction(:muladd), [add_constant!(ls, _mlt, sizeof(Int)), _op, opbase], sizeof(Int))
      end
    end
    addopindex!(opparents, loopdependencies, reduceddeps, indices, offsets, strides, loopedindex, opbase, one(Int8), offset%Int8)
  end
  return vptrarray
end

function move_to_last!(x, i)
    i == length(x) && return
    xᵢ = x[i]
    deleteat!(x, i)
    push!(x, xᵢ)
    nothing
end
# TODO: Make this work with Cartesian Indices
function repeated_index!(ls::LoopSet, indices::Vector{Symbol}, vptr::Symbol, indnum::Int, firstind::Int)
    # Move ind to last position
    vptrrepremoved = Symbol(vptr, "##ind##", firstind, "##repeated##", indnum, "##")
    f = Expr(:(.), Expr(:(.), :LoopVectorization, QuoteNode(:VectorizationBase)), QuoteNode(:double_index))
    fiv = Expr(:call, Expr(:curly, :Val, firstind - 1))
    siv = Expr(:call, Expr(:curly, :Val, indnum - 1))
    pushpreamble!(ls, Expr(:(=), vptrrepremoved, Expr(:call, f, vptr, fiv, siv)))
    vptrrepremoved
end

function array_reference_meta!(ls::LoopSet, array::Symbol, rawindices, elementbytes::Int, var::Union{Nothing,Symbol} = nothing)
  vptrarray = vptr(array)
  add_vptr!(ls, array, vptrarray) # now, subset
  indices = Symbol[]
  offsets = Int8[]
  strides = Int8[]
  loopedindex = Bool[]
  parents = Operation[]
  loopdependencies = Symbol[]
  reduceddeps = Symbol[]
  loopset = ls.loopsymbols
  ninds = 1
  D = length(rawindices)
  for ind ∈ rawindices
    if ind isa Integer # subset
      if byterepresentable(ind)
        addconstindex!(indices, offsets, strides, loopedindex, ind)
        ninds += 1
      else
        vptrarray = subset_vptr!(ls, vptrarray, ninds, ind, indices, loopedindex, 0)
        length(indices) == 0 && push!(indices, DISCONTIGUOUS)
      end
    elseif ind isa Expr
      #FIXME: position (in loopnest) wont be length(ls.loopsymbols) in general
      vptrarray = checkforoffset!(
        ls, vptrarray, ninds, parents, indices, offsets, strides, loopedindex, loopdependencies, reduceddeps, ind, D
      )
      ninds += 1
    elseif ind isa Symbol
      if ind ∈ loopset
        ind_prev_index = findfirst(Base.Fix2(===,ind), indices)
        if ind_prev_index === nothing
          push!(indices, ind); ninds += 1
          push!(offsets, zero(Int8))
          push!(strides, one(Int8))
          push!(loopedindex, true)
          push!(loopdependencies, ind)
        else
          move_to_last!(indices, ind_prev_index)
          move_to_last!(offsets, ind_prev_index)
          move_to_last!(strides, ind_prev_index)
          move_to_last!(loopedindex, ind_prev_index)
          move_to_last!(loopdependencies, ind_prev_index)
          vptrarray = repeated_index!(ls, indices, vptrarray, ninds, ind_prev_index + (first(indices) === DISCONTIGUOUS))
          makediscontiguous!(indices)
        end
      else
        indop = get(ls.opdict, ind, nothing)
        if indop !== nothing  && !isconstant(indop)
          if iscompute(indop) && Base.sym_in(instruction(indop).instr, (:add_fast, :sub_fast))
            vptrarray = add_additive_index!(ls, parents, vptrarray, ind, indop, ninds, indices, offsets, strides, loopedindex, loopdependencies, reduceddeps, 0, 1, D)
            ninds += 1
          else
            pushparent!(parents, loopdependencies, reduceddeps, indop)
            push!(indices, name(indop)); ninds += 1
            push!(offsets, zero(Int8))
            push!(strides, one(Int8))
            push!(loopedindex, false)
          end
        else
          vptrarray = subset_vptr!(ls, vptrarray, ninds, ind, indices, loopedindex, 0)
          length(indices) == 0 && push!(indices, DISCONTIGUOUS)
        end
      end
    else
      throw("Unrecognized loop index: $ind.")
    end
  end
  mref = ArrayReferenceMeta(ArrayReference( array, indices, offsets, strides ), loopedindex, vptrarray)
  ArrayReferenceMetaPosition(mref, parents, loopdependencies, reduceddeps, var === nothing ? Symbol("") : var )
end
function tryrefconvert(ls::LoopSet, ex::Expr, elementbytes::Int, var::Union{Nothing,Symbol} = nothing)::Tuple{Bool,ArrayReferenceMetaPosition}
    ya, yinds = if ex.head === :ref
        ref_from_ref!(ls, ex)
    elseif ex.head === :call
        f = first(ex.args)
        if f === :getindex
            ref_from_getindex!(ls, ex)
        elseif f === :setindex!
            ref_from_setindex!(ls, ex)
        else
            return false, NOTAREFERENCEMP
        end
    else
        return false, NOTAREFERENCEMP
    end
    true, array_reference_meta!(ls, ya, yinds, elementbytes, var)
end


text:
# Scheduling of processes in GeMM

"""
    simulate!(world, settings, timesteps)

This is the central function of the model with the main event loop. It defines
the scheduling for all submodels and output functions.
"""
function simulate!(world::Array{Patch,1}, settings::Dict{String, Any}, timesteps::Int=1000, timeoffset::Int = 0)
    simlog("Starting simulation.")
    checkviability!(world, settings)
    for t in (timeoffset + 1):(timeoffset + timesteps)
        simlog("UPDATE $t")
        # ecological processes are outsourced to specialised methods below
        if settings["mode"] == "default"
            defaultexperiment(world, settings)
        elseif settings["mode"] == "invasion"
            invasionexperiment(world, settings, t)
        elseif settings["mode"] == "zosterops"
            zosteropsexperiment(world, settings)
        else
            simlog("Mode setting not recognised: $(settings["mode"])", 'e')
        end
        if settings["lineages"]
            recordstatistics(world, settings)
            recordlineages(world, settings, t)
        end
        if mod(t, settings["outfreq"]) == 0 && any([settings["fasta"] != "off", settings["raw"], settings["stats"]])
            writedata(world, settings, t)
        end
    end
end

"""
    defaultexperiment(world, settings)

The standard annual update procedure, designed primarily for plant communities.
"""
function defaultexperiment(world::Array{Patch,1}, settings::Dict{String, Any})
    establish!(world, settings["nniches"], settings["static"])
    survive!(world, settings)
    grow!(world, settings)
    compete!(world, settings["static"])
    reproduce!(world, settings)
    if settings["mutate"]
        mutate!(world, settings)
    end
    disperse!(world, settings["static"])
    checkviability!(world, settings)
    changehabitat!(world, settings) # model output
end

"""
    invasionexperiment(world, settings)

The annual update procedure for the invasion experiments.
"""
function invasionexperiment(world::Array{Patch,1}, settings::Dict{String, Any}, t::Int)
    establish!(world, settings["nniches"], settings["static"])
    survive!(world, settings)
    grow!(world, settings)
    compete!(world, settings["static"])
    reproduce!(world, settings)
    if 0 < settings["burn-in"] < t
        disturb!(world, settings)
        invade!(world, settings)
    end
    disperse!(world, settings["static"])
    checkviability!(world, settings)
end

"""
    zosteropsexperiment(world, settings)

The annual update procedure for the Zosterops experiments, this time for bird populations.
"""
function zosteropsexperiment(world::Array{Patch,1}, settings::Dict{String, Any})
    #TODO
    establish!(world, settings["nniches"], settings["static"])
    survive!(world, settings)
    grow!(world, settings)
    reproduce!(world, settings)
    if settings["mutate"]
        mutate!(world, settings)
    end
    zdisperse!(world, settings)
    checkviability!(world, settings)
    #changehabitat!(world, settings) # model output
end

# TODO add more modes for gradient/habitat change experiments?

ext: jl

max_stars_repo_name: lleiding/gemm

max_issues_repo_name: lleiding/gemm

max_forks_repo_name: lleiding/gemm


text:
#=
  Port of WebPPL code from https://github.com/probmods/ppaml2016/blob/gh-pages/chapters/5-election.md
  """
  Basic model
  Learning a state-wide preference from poll data
  """

  Clinton: 1 Trump: 2
  Distributions of variable (num:0)
  2.00000 =>    8750  (0.875000)
  1.00000 =>    1250  (0.125000)

  According to the WebPPL model, the probability that Clinton wins (won) is 0.12207
=#

using Turing
include("jl_utils.jl")

@model function sampleElection()
    function trueStatePref()
        pref ~ Beta(1,1)
        counts = Dict(:trump=>304, :clinton=>276)
        total = sum(values(counts))
        c ~ Binomial(total, pref)
        true ~ Dirac(c == counts[:clinton])
        return pref
    end
    clinton = 0
    trump = 1
    pref = trueStatePref()
    turnout = 2400000
    clintonVotes ~ Binomial(turnout,pref)
    trumpVotes ~ Dirac(turnout - clintonVotes)
    winner ~ Dirac(clintonVotes > trumpVotes ? clinton : trump)

    clintonWins ~ Dirac(winner == clinton)
    trumpWins ~ Dirac(winner == trump)
    return winner
end

model = sampleElection()

# chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(15), 10_000)
# chns = sample(model, IS(), 10_000)
chns = sample(model, SMC(), 10_000)
# chns = sample(model, SMC(), MCMCThreads(), 10_000, 4)

display(chns)

chns_params = Turing.MCMCChains.get_sections(chns, :parameters)
genq = generated_quantities(model, chns_params)
println("Clinton: 0 Trump: 1")
show_var_dist_pct(genq)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
type LmResp{T<:FP} <: ModResp  # response in a linear model
    mu::Vector{T}     # mean response
    offset::Vector{T} # offset added to linear predictor (may have length 0)
    wts::Vector{T}    # prior weights (may have length 0)
    y::Vector{T}      # response
    function LmResp(mu::Vector{T}, off::Vector{T}, wts::Vector{T}, y::Vector{T})
        n = length(y); length(mu) == n || error("mismatched lengths of mu and y")
        ll = length(off); ll == 0 || ll == n || error("length of offset is $ll, must be $n or 0")
        ll = length(wts); ll == 0 || ll == n || error("length of wts is $ll, must be $n or 0")
        new(mu,off,wts,y)
    end
end
LmResp{T<:FP}(y::Vector{T}) = LmResp{T}(zeros(T,length(y)), T[], T[], y)

function updatemu!{T<:FP}(r::LmResp{T}, linPr::Vector{T})
    n = length(linPr); length(r.y) == n || error("length(linPr) is $n, should be $(length(r.y))")
    length(r.offset) == 0 ? copy!(r.mu, linPr) : map!(Add(), r.mu, linPr, r.offset)
    deviance(r)
end
updatemu!{T<:FP}(r::LmResp{T}, linPr) = updatemu!(r, convert(Vector{T},vec(linPr)))

type WtResid <: TernaryFunctor end
evaluate{T<:FP}(::WtResid,wt::T,y::T,mu::T) = (y - mu)*sqrt(wt)
result_type{T<:FP}(::WtResid,wt::T,y::T,mu::T) = T

deviance(r::LmResp) = length(r.wts) == 0 ? sumsqdiff(r.y, r.mu) : wsumsqdiff(r.wts,r.y,r.mu)
residuals(r::LmResp)= length(r.wts) == 0 ? r.y - r.mu : map(WtResid(),r.wts,r.y,r.mu)

type LmMod <: LinPredModel
    fr::ModelFrame
    rr::LmResp
    pp::LinPred
    ff::Formula
end

function lm(f::Formula, df::AbstractDataFrame)
    mf = ModelFrame(f, df); mm = ModelMatrix(mf)
    rr = LmResp(model_response(mf)); pp = DensePredQR(mm)
    installbeta!(delbeta!(pp, rr.y)); updatemu!(rr, linpred(pp,0.))
    LmMod(mf, rr, pp, f)
end
lm(f::Expr, df::AbstractDataFrame) = lm(Formula(f), df)
lm(f::String, df::AbstractDataFrame) = lm(Formula(parse(f)[1]), df)

function lmc(f::Formula, df::AbstractDataFrame)
    mf = ModelFrame(f, df); mm = ModelMatrix(mf)
    rr = LmResp(model_response(mf)); pp = DensePredChol(mm)
    installbeta!(delbeta!(pp, rr.y)); updatemu!(rr, linpred(pp,0.))
    LmMod(mf, rr, pp, f)
end
lmc(f::Expr, df::AbstractDataFrame) = lmc(Formula(f), df)
lmc(f::String, df::AbstractDataFrame) = lmc(Formula(parse(f)[1]), df)

## scale(m) -> estimate, s, of the scale parameter
## scale(m,true) -> estimate, s^2, of the squared scale parameter
function scale(x::LmMod, sqr::Bool=false)
    ssqr = deviance(x.rr)/df_residual(x)
    sqr ? ssqr : sqrt(ssqr)
end

function coeftable(mm::LmMod)
    cc = coef(mm)
    se = stderr(mm)
    tt = cc ./ se
    DataFrame({cc, se, tt, ccdf(FDist(1, df_residual(mm)), tt .* tt)},
              ["Estimate","Std.Error","t value", "Pr(>|t|)"])
end

function predict(mm::LmMod, newx::Matrix)
    newx * coef(mm)
end

function confint(obj::LmMod, level::Real)
    cft = coeftable(obj)
    hcat(coef(obj),coef(obj)) + cft["Std.Error"] *
    quantile(TDist(df_residual(obj)), (1. - level)/2.) * [1. -1.]
end
confint(obj::LmMod) = confint(obj, 0.95)


text:
using ClimateMachine
const clima_dir = dirname(dirname(pathof(ClimateMachine)));

if parse(Bool, get(ENV, "CLIMATEMACHINE_PLOT_EDMF_COMPARISON", "false"))
    plot_dir = joinpath(clima_dir, "output", "sbl_edmf", "pycles_comparison")
else
    plot_dir = nothing
end

include(joinpath(@__DIR__, "compute_mse.jl"))

data_file = Dataset(joinpath(PyCLES_output_dataset_path, "Gabls.nc"), "r")

#! format: off
best_mse = OrderedDict()
best_mse["prog_ρ"] = 6.8041561724036673e-03
best_mse["prog_ρu_1"] = 6.2586216904022822e+03
best_mse["prog_ρu_2"] = 1.2965826326450741e-04
best_mse["prog_turbconv_environment_ρatke"] = 4.9035225536000081e+02
best_mse["prog_turbconv_environment_ρaθ_liq_cv"] = 8.7727377301948991e+01
best_mse["prog_turbconv_updraft_1_ρa"] = 1.8213581913998944e+01
best_mse["prog_turbconv_updraft_1_ρaw"] = 1.7800899665237452e-01
best_mse["prog_turbconv_updraft_1_ρaθ_liq"] = 1.3358308964295857e+01
#! format: on

computed_mse = compute_mse(
    solver_config.dg.grid,
    solver_config.dg.balance_law,
    time_data,
    dons_arr,
    data_file,
    "Gabls",
    best_mse,
    60,
    plot_dir,
)

@testset "SBL EDMF Solution Quality Assurance (QA) tests" begin
    #! format: off
    test_mse(computed_mse, best_mse, "prog_ρ")
    test_mse(computed_mse, best_mse, "prog_ρu_1")
    test_mse(computed_mse, best_mse, "prog_ρu_2")
    test_mse(computed_mse, best_mse, "prog_turbconv_environment_ρatke")
    test_mse(computed_mse, best_mse, "prog_turbconv_environment_ρaθ_liq_cv")
    test_mse(computed_mse, best_mse, "prog_turbconv_updraft_1_ρa")
    test_mse(computed_mse, best_mse, "prog_turbconv_updraft_1_ρaw")
    test_mse(computed_mse, best_mse, "prog_turbconv_updraft_1_ρaθ_liq")
    #! format: on
end


text: calc_T_k(reactor, R_0, B_0) = K_T(reactor) * B_0 * R_0


text:
"""
# Module API

Application Programming Interface (API) of JustSayIt.

#### Macros
- [`@voiceargs`](@ref)

#### FUNCTIONS
- [`is_next`](@ref)
- [`are_next`](@ref)
- [`pause_recording`](@ref)
- [`restart_recording`](@ref)

To see a description of a function, macro or module type `?<functionname>`, `?<macroname>` (including the `@`) or `?<modulename>`, respectively.
"""
module API
    import ..JustSayIt: @voiceargs, is_next, are_next, pause_recording, restart_recording
    export @voiceargs, is_next, are_next, pause_recording, restart_recording
end


text:
using FactCheck
using ValidatedNumerics

facts("DecoratedInterval tests") do
    a = DecoratedInterval(@interval(1, 2), com)
    @fact decoration(a) --> com

    b = sqrt(a)
    @fact interval_part(b) --> sqrt(interval_part(a))
    @fact decoration(b) --> com

    a = DecoratedInterval(@interval(-1, 1), com)
    b = sqrt(a)
    @fact interval_part(b) --> sqrt(Interval(0, 1))
    @fact decoration(b) --> trv

    d = DecoratedInterval(a, dac)
    @fact decoration(d) --> dac

    @fact decoration(DecoratedInterval(1.1)) --> com
    @fact decoration(DecoratedInterval(1.1, dac)) --> dac

    @fact decoration(DecoratedInterval(2, 0.1, com)) --> ill
    @fact decoration(DecoratedInterval(2, 0.1)) --> ill
    @fact isnai(interval_part(DecoratedInterval(2, 0.1))) --> true
    @fact decoration(@decorated(2, 0.1)) --> ill
    @fact decoration(DecoratedInterval(big(2), big(1))) --> ill
    @fact isnai(interval_part((DecoratedInterval(big(2), big(1))))) --> true
    @fact isnai(interval_part(@decorated(big(2), big(1)))) --> true

    # Disabling the following tests, because Julia 0.5 has some strange behaviour here
    # @fact_throws ArgumentError DecoratedInterval(BigInt(1), 1//10)
    # @fact_throws ArgumentError @decorated(BigInt(1), 1//10)
end


text:
# contains main ASWING calling routines

"""
    setgeom(asw::Configuration)

    Inputs geometry into ASWING
"""
function setgeom(asw::Configuration)

    # clear variable array sizes
    ASWING.NB[:, :] .= 0
    ASWING.NANGJ[:] .= 0
    npylo = 0 #pylon counter

    # NAME
    for i = 1:length(asw.namein)
        ASWING.NAME[i] = asw.namein[i]
    end

    # UNITS
    if asw.si_units
        UNITS.UNITL[1] = 1.0
        UNITS.UNCHL[:] = Vector{Char}(rpad("m",32))
        UNITS.UNITT[1] = 1.0
        UNITS.UNCHT[:] = Vector{Char}(rpad("s",32))
        UNITS.UNITM[1] = 1.0
        UNITS.UNCHM[:] = Vector{Char}(rpad("kg",32))
        UNITS.UNITF[1] = 0.0
        UNITS.UNCHF[:] = Vector{Char}(rpad(" ",32))
    else
        UNITS.UNITL[1] = 1.0
        UNITS.UNCHL[:] = Vector{Char}(rpad("ft",32))
        UNITS.UNITT[1] = 1.0
        UNITS.UNCHT[:] = Vector{Char}(rpad("s ",32))
        UNITS.UNITM[1] = 0.0
        UNITS.UNCHM[:] = Vector{Char}(rpad("  ",32))
        UNITS.UNITF[1] = 1.0
        UNITS.UNCHF[:] = Vector{Char}(rpad("lb",32))
    end

    # CONSTANTS BLOCK
    ASWING.GEEW[1] = asw.constants.gravity
    ASWING.RHOSL[1] = asw.constants.sea_level_air_density
    ASWING.VSOSL[1] = asw.constants.sea_level_speed_of_sound

    # REFERENCE BLOCK
    ASWING.SREF[1] = asw.refvals.area
    ASWING.CREF[1] = asw.refvals.chord
    ASWING.BREF[1] = asw.refvals.span
    ASWING.XYZREF[:,1] = asw.refvals.moments
    ASWING.XYZREF[:,2] = asw.refvals.accelerations
    ASWING.XYZREF[:,3] = asw.refvals.velocities

    # SENSOR BLOCK
    nsen = length(asw.sensors)
    for i = 1:nsen
        ASWING.KPTYPE[npylo+i] = 100 + asw.sensors[i].Ksens
        ASWING.KBPYLO[npylo+i] = asw.sensors[i].Nbeam
        ASWING.QPYLO[0,npylo+i] = asw.sensors[i].t
        ASWING.QPYLO[1:3,npylo+i] = asw.sensors[i].location
        ASWING.QPYLO[4:6,npylo+i] = asw.sensors[i].V
        ASWING.QPYLO[7:9,npylo+i] = asw.sensors[i].A
    end
    npylo += nsen

    # ENGINE BLOCK
    neng = length(asw.engines)
    for i = 1:neng
        keng = asw.engines[i].Keng
        ietyp = asw.engines[i].IEtyp
        ASWING.KBPYLO[npylo+i] = asw.engines[i].Nbeam
        ASWING.QPYLO[0,npylo+i] = asw.engines[i].t
        ASWING.QPYLO[1:3,npylo+i] = asw.engines[i].location
        ASWING.QPYLO[4:6,npylo+i] = asw.engines[i].thrust_direction
        ASWING.QPYLO[7,npylo+i] = asw.engines[i].dFdPe
        ASWING.QPYLO[8,npylo+i] = asw.engines[i].dMdPe
        ASWING.QPYLO[9,npylo+i] = asw.engines[i].Rdisk
        ASWING.QPYLO[10,npylo+i] = asw.engines[i].Omega
        ASWING.QPYLO[11,npylo+i] = asw.engines[i].cdA
        ASWING.QPYLO[12,npylo+i] = asw.engines[i].cl
        ASWING.QPYLO[13,npylo+i] = asw.engines[i].CLa
        ASWING.QPYLO[14,npylo+i] = asw.engines[i].S0
        ASWING.QPYLO[15,npylo+i] = asw.engines[i].C0
        ASWING.QPYLO[16,npylo+i] = asw.engines[i].S1
        ASWING.QPYLO[17,npylo+i] = asw.engines[i].C1
        ASWING.QPYLO[18,npylo+i] = asw.engines[i].S2
        ASWING.QPYLO[19,npylo+i] = asw.engines[i].C2
        ASWING.QPYLO[20,npylo+i] = asw.engines[i].S3
        ASWING.QPYLO[21,npylo+i] = asw.engines[i].C3

        ASWING.KPTYPE[npylo+i] = 10 + keng
        ke = max(1, min(NENGX, keng))
        ASWING.IENGTYP[ke] = ietyp
    end
    npylo += neng

    # WEIGHT BLOCK
    nwei = length(asw.weights)
    for i = 1:nwei
        ASWING.KPTYPE[npylo+i] = 1
        ASWING.KBPYLO[npylo+i] = asw.weights[i].Nbeam
        ASWING.QPYLO[0,npylo+i] = asw.weights[i].t
        ASWING.QPYLO[1:3,npylo+i] = asw.weights[i].location
        ASWING.QPYLO[4,npylo+i] = asw.weights[i].weight
        ASWING.QPYLO[5,npylo+i] = asw.weights[i].CDA
        ASWING.QPYLO[6,npylo+i] = asw.weights[i].vol
        ASWING.QPYLO[7:9,npylo+i] = asw.weights[i].Ho
        ASWING.QPYLO[10,npylo+i] = asw.weights[i].Ixx
        ASWING.QPYLO[11,npylo+i] = asw.weights[i].Iyy
        ASWING.QPYLO[12,npylo+i] = asw.weights[i].Izz
        ASWING.QPYLO[13,npylo+i] = asw.weights[i].Ixy
        ASWING.QPYLO[14,npylo+i] = asw.weights[i].Ixz
        ASWING.QPYLO[15,npylo+i] = asw.weights[i].Iyz
    end
    npylo += nwei

    # STRUT BLOCK
    nstr = length(asw.struts)
    for i = 1:nstr
        ASWING.KPTYPE[npylo+i] = 2
        ASWING.KBPYLO[npylo+i] = asw.struts[i].Nbeam
        ASWING.QPYLO[0,npylo+i] = asw.struts[i].t
        ASWING.QPYLO[1:3,npylo+i] = asw.struts[i].location
        ASWING.QPYLO[4:6,npylo+i] = asw.struts[i].anchor_location
        ASWING.QPYLO[7,npylo+i] = asw.struts[i].dLo
        ASWING.QPYLO[8,npylo+i] = asw.struts[i].EAw
    end
    npylo += nstr
    ASWING.NPYLO[1] = npylo

    # JOINT BLOCK
    njoin = length(asw.joints)
    ASWING.NJOIN[1] = njoin
    for i = 1:njoin
        ASWING.KBJOIN[1,i] = asw.joints[i].Nbeam1
        ASWING.KBJOIN[2,i] = asw.joints[i].Nbeam2
        ASWING.TJOIN[1,i] = asw.joints[i].t1
        ASWING.TJOIN[2,i] = asw.joints[i].t2
        ASWING.KJTYPE[i] = asw.joints[i].KJtype
    end

    # JANGLE BLOCK
    for i = 1:length(asw.jangles)
        kj = asw.jangles[i].Njoint
        nangj = length(asw.jangles[i].hinge_moment)
        ASWING.NANGJ[kj] = nangj
        ASWING.HJAX[1:3,kj] = asw.jangles[i].hinge_axis
        ASWING.MOMJ[1:nangj,kj] = asw.jangles[i].hinge_moment
        ASWING.ANGJ[1:nangj,kj] = asw.jangles[i].hinge_angle
    end

    # GROUNDS BLOCK
    ngrou = length(asw.grounds)
    ASWING.NGROU[1] = ngrou
    for i = 1:ngrou
        ASWING.KBGROU[i] = asw.grounds[i].Nbeam
        ASWING.TGROU[i] = asw.grounds[i].t
        ASWING.KGTYPE[i] = asw.grounds[i].KGtype
    end

    # BEAM DEFINITION BLOCKS
    nbeam = length(asw.beams)
    ASWING.NBEAM[1] = nbeam
    for ibeam = 1:nbeam
        ASWING.BNAME[ibeam][1:length(asw.beams[ibeam].name)] = Vector{Char}(asw.beams[ibeam].name)
        ASWING.KBNUM[ibeam] = asw.beams[ibeam].number
        ASWING.IBEAM[ibeam] = asw.beams[ibeam].physical_index
        for var in keys(beamvarval)
            j = beamvarval[var]
            varval = getfield(asw.beams[ibeam].spanwise_variables,Symbol(var))
            nb = size(varval,1)
            ASWING.NB[j, ibeam] = nb
            ASWING.TB[1:nb, j, ibeam] = varval[:,1]
            ASWING.QB[1:nb, j, ibeam] = varval[:,2]
        end
    end

    ccall((:input_, libaswing), Nothing, (Ref{Cchar},), ' ')
    ccall((:qmset_, libaswing), Nothing, ())
    ccall((:dscheck_, libaswing), Nothing, ())
    return nothing
end

"""
    setpnt(point)
    Inputs operating point(s) into ASWING
"""
function setpnt(point::AbstractArray{<:OperatingPoint,1})

    # set number of operating points
    npoint = length(point)
    ASWING.NPOINT[1] = npoint
    for ipnt = 1:npoint
        # set constraints
        setcons(point[ipnt].constraints, ipnt)
        # set parameter values
        setpars(point[ipnt].parameters, ipnt)
        # set other options and weights
        ASWING.LMACH[ipnt] = point[ipnt].mach_from_airspeed
        ASWING.IGRIM[ipnt] = point[ipnt].ground_image
        ASWING.WFLAP[:,ipnt] = point[ipnt].wflap
        ASWING.WPENG[:,ipnt] = point[ipnt].wpeng
        # solution is unconverged
        ASWING.LCONV[ipnt] = false
    end

    # set Mach number for Prandtl-Glauert correction
    for ipnt = 1:npoint
        for ip = 1:KPFREE
            ccall((:psp2par_, libaswing), Nothing, (Ref{Int32}, Ref{Int32}), ip, ipnt)
        end
        ASWING.PARAM[KPMACH, ipnt] =  point[ipnt].machpg
    end

    return nothing
end

"""
    setcons(pnt::dotpnt,ipnt=1)
    Inputs operating points contraints into ASWING
"""
function setcons(constraints, ipnt::Integer=1)
    ASWING.IPPAR[KPVACX:KPVACZ,ipnt] = constraints.linear_acceleration
    ASWING.IPPAR[KPRACX:KPRACZ,ipnt] = constraints.angular_acceleration
    ASWING.IPPAR[KPVELX:KPVELZ,ipnt] = constraints.velocity
    ASWING.IPPAR[KPROTX:KPROTZ,ipnt] = constraints.rotation_rate
    ASWING.IPPAR[KPPOSX:KPPOSZ,ipnt] = constraints.position
    ASWING.IPPAR[KPBANK,ipnt] = constraints.phi
    ASWING.IPPAR[KPELEV,ipnt] = constraints.theta
    ASWING.IPPAR[KPHEAD,ipnt] = constraints.psi
    ASWING.IPPAR[KPFLP1:KPFLP1+length(constraints.flap_defl_ctrl_var)-1,ipnt] = constraints.flap_defl_ctrl_var
    for ip = KPFLP1+length(constraints.flap_defl_ctrl_var):KPFLP1+NFLPX
        ASWING.IPPAR[ip,ipnt] = ip
    end
    ASWING.IPPAR[KPENG1:KPENG1+length(constraints.eng_pwr_ctrl_var)-1,ipnt] = constraints.eng_pwr_ctrl_var
    for ip = KPENG1+length(constraints.eng_pwr_ctrl_var):KPENG1+NENGX
        ASWING.IPPAR[ip,ipnt] = ip
    end
    for kf=1 : NFLPX
        kpflp1 = KPFLP1[1]
        ipflp1 = IPFLP1[1]
        if ASWING.LSFLAP[kf] == 1
            if ASWING.IPPAR[kpflp1-1+kf,ipnt] == 0
                ASWING.IPPAR[kpflp1-1+kf,ipnt] = ipflp1-1+kf
            end
        else
            ASWING.IPPAR[kpflp1-1+kf,ipnt] = 0
        end
    end
    for ke=1 : NENGX
        kpeng1 = KPENG1[1]
        ipeng1 = IPENG1[1]
        if ASWING.LSPENG[ke] == 1
            if ASWING.IPPAR[kpeng1-1+ke, ipnt] == 0
                ASWING.IPPAR[kpeng1-1+ke, ipnt] = ipeng1-1 + ke
            end
        else
            ASWING.IPPAR[kpeng1[1]-1+ke,ipnt] = 0
        end
    end
    ASWING.IPPAR[KPVIDT,ipnt] = constraints.err_int_Vinf
    ASWING.IPPAR[KPBEDT,ipnt] = constraints.err_int_beta
    ASWING.IPPAR[KPALDT,ipnt] = constraints.err_int_alpha
    ASWING.IPPAR[KPBADT,ipnt] = constraints.err_int_phi
    ASWING.IPPAR[KPELDT,ipnt] = constraints.err_int_theta
    ASWING.IPPAR[KPHEDT,ipnt] = constraints.err_int_psi
    ASWING.IPPAR[KPWDTX:KPWDTZ,ipnt] = constraints.err_int_ROTxyz
    ASWING.IPPAR[KPADTX:KPADTZ,ipnt] = constraints.err_int_VACxyz
    ASWING.LCONV[ipnt] = false
    return nothing
end

"""
    setpars(parameters, ipnt::Integer=1)
    Sets operating point parameters into ASWING
"""
function setpars(parameters, ipnt::Integer=1)

    ASWING.PSPEC[IPVACX:IPVACZ,ipnt] = parameters.linear_acceleration
    ASWING.PSPEC[IPRACX:IPRACZ,ipnt] = parameters.angular_acceleration
    ASWING.PSPEC[IPVIAS,ipnt] = parameters.velocity
    ASWING.PSPEC[IPBETA,ipnt] = parameters.beta
    ASWING.PSPEC[IPALFA,ipnt] = parameters.alpha
    ASWING.PSPEC[IPROTX:IPROTZ,ipnt] = parameters.rotation_rate
    ASWING.PSPEC[IPPOSX:IPPOSZ,ipnt] = parameters.position
    ASWING.PSPEC[IPBANK,ipnt] = parameters.phi
    ASWING.PSPEC[IPELEV,ipnt] = parameters.theta
    ASWING.PSPEC[IPHEAD,ipnt] = parameters.psi
    ASWING.PSPEC[IPFLP1:IPFLP1+length(parameters.flap_defl_ctrl_var)-1,ipnt] = parameters.flap_defl_ctrl_var
    for ip = IPFLP1+length(parameters.flap_defl_ctrl_var):IPFLP1+NFLPX
      ASWING.PSPEC[ip,ipnt] = 0.0
    end
    ASWING.PSPEC[IPENG1:IPENG1+length(parameters.eng_pwr_ctrl_var)-1,ipnt] = parameters.eng_pwr_ctrl_var
    for ip = IPENG1+length(parameters.eng_pwr_ctrl_var):IPENG1+NENGX
      ASWING.PSPEC[ip,ipnt] = 0.0
    end
    ASWING.PSPEC[IPFORX:IPFORZ,ipnt] = parameters.sum_force
    ASWING.PSPEC[IPMOMX:IPMOMZ,ipnt] = parameters.sum_mom
    ASWING.PSPEC[IPLIFT,ipnt] = parameters.lift
    ASWING.PSPEC[IPGAMM,ipnt] = parameters.climb_angle
    ASWING.PSPEC[IPRACC,ipnt] = parameters.radial_acceleration
    ASWING.PSPEC[IPUSR1,ipnt] = parameters.usr1
    ASWING.PSPEC[IPUSR2,ipnt] = parameters.usr2
    ASWING.PSPEC[IPLSQD,ipnt] = parameters.least_squared

    return nothing
end

"""
    solvesteady(;repeat::Integer=0)
    solvesteady(ipnt1::Integer; repeat::Integer=0)
    solvesteady(ipnt1::Integer, ipnt2::Integer; repeat::Integer=0)

    Quasi-steady ASWING solve. If no arguments are specified, all points are
    solved.
"""
function solvesteady(ipnt1::Integer, ipnt2::Integer; repeat::Integer=0)

    ASWING.STEADY[1] = true
    for ipnt in ipnt1:ipnt2
        # set operating point
        ASWING.ARGP1[1:4] .= Vector{Char}(rpad("$ipnt",4)) # hack to avoid string passing
        ccall((:operjl_, libaswing), Nothing, (Ref{Float64}, Ref{Int32}), 0.0, 0)
        for i = 1:repeat+1
            # converge operating point
            ASWING.ARGP1[1:4] .= Vector{Char}(rpad("X",4)) # hack to avoid string passing
            ccall((:operjl_, libaswing), Nothing, (Ref{Float64}, Ref{Int32}), 0.0, 0)
            if Bool(ASWING.LCONV[ipnt])
                break
            end
        end
    end

    # return formatted solution
    return getsolution(ipnt1,ipnt2)
end
solvesteady(; repeat::Integer=0) = solvesteady(1, ASWING.NPOINT[1], repeat=repeat)
solvesteady(ipnt1::Integer; repeat::Integer=0) = solvesteady(ipnt1, ipnt1, repeat=repeat)[ipnt1]

"""
    solveunsteady(deltat::Real, ntimes::Integer)
    Unsteady ASWING time-domain solve
"""
function solveunsteady(deltat::Real, ntimes::Integer, ipnt::Integer=1)

    # set operating point
    ASWING.ARGP1[1:4] .= Vector{Char}(rpad("$ipnt",4)) # hack to avoid string passing
    ccall((:operjl_, libaswing), Nothing, (Ref{Float64}, Ref{Int32}), 0.0, 0)

    # change to unsteady
    ASWING.STEADY[1] = false

    # solve unsteady
    ASWING.ARGP1[1:4] .= Vector{Char}(rpad("X",4)) # hack to avoid string passing
    ccall((:operjl_, libaswing), Nothing, (Ref{Float64}, Ref{Int32}), deltat, ntimes)

    # return formatted solution
    return getsolution(1, ntimes+1)
end

"""
    retainparam()
    retainparam(ipnt1::Integer)
    retainparam(ipnt1::Integer, ipnt2::Integer)
    Sets operating point parameters equivalent to those of the current state
"""
function retainparam(ipnt1::Integer, ipnt2::Integer=ipnt1)

    for ipnt in ipnt1:ipnt2
        # set operating point
        ASWING.ARGP1[1:4] .= Vector{Char}(rpad("$ipnt",4)) #hack to avoid string passing
        ccall((:operjl_, libaswing), Nothing, (Ref{Float64}, Ref{Int32}), 0.0, 0)
        # set parameters on operating point
        ASWING.ARGP1[1:4] .= Vector{Char}(rpad("R",4)) #hack to avoid string passing
        ccall((:operjl_, libaswing), Nothing, (Ref{Float64}, Ref{Int32}), 0.0, 0)
    end

    return nothing
end
retainparam() = retainparam(1, ASWING.NPOINT[1])

"""
    freeze(ipoint = 1)
    Freezes geometry at current state
"""
function freeze(ipoint = 1)

  ccall((:freeze_,libaswing), Nothing, (Ref{Int32},), ipoint)

  for i = 1 : ASWING.NPOINT[1]
    ASWING.LCONV[i] = false
  end

end

"""
    setgust(igust::Integer, params::AbstractArray{<:Real}=[-4.0,0.0,1.0,8.0,2.0,0.0])
    Sets gust type and corresponding parameters
"""
function setgust(igust::Integer, params::AbstractArray{<:Real}=[-4.0,0.0,1.0,8.0,2.0,0.0])

    gustparam = [-4.0,0.0,1.0,8.0,2.0,0.0]
    if length(params) > 6
        error("No more than six gust parameters may be specified")
    end
    gustparam[1:length(params)] = params
    ASWING.IGUST[1] = igust
    ASWING.VGCON[1:6] = gustparam

    return nothing
end

"""
    geteigs(neigen, zshift, ipnt; evmin = 1e-4, evtol = 1e-3)
    Returns eigenvalues of the system
"""
function geteigs(neig, zshift, ipnt; evmin = 1e-4, evtol = 1e-3)

    ASWING.STEADY[1] = false
    ASWING.LFGUST[1] = true

    # number of new modes to compute
    ASWING.NBEIGEN[1] = neig

    # clear current eigenmode counter, and start from new
    EIGEN.NEIGEN[1] = 0
    EIGEN.NROM[1] = 0

    # fix shift if too close to origin
    if abs(zshift) < 0.01
        zshift = -0.01 + 0.0im
    end

    # use shift as new default seed
    EIGEN.ZESEED[1] = zshift

    # compute only eigenvalues for specified operating points
    lrvec = true
    lgerr = false
    ccall((:evexec_, libaswing), Nothing, (Ref{Int32}, Ref{ComplexF64}, Ref{Int32},
        Ref{Int32}, Ref{Int32}, Ref{Float64}, Ref{Float64}, Ref{Int32}),
        ipnt, zshift, ASWING.NBEIGEN[1], lrvec, lgerr, evmin, evtol, 0)

    ccall((:wbcalc_, libaswing), Nothing, (Ref{Int32},), ipnt)

    zlamda = copy(EIGEN.ZLAMDA[1:EIGEN.NEIGEN[1]])

    return zlamda
end

"""
    stabilityderivatives(ipnt = 1, momref=Float64[])
    Calculates stability derivatives
"""
function stabilityderivatives(ipnt::Integer = 1; momref::AbstractArray{<:Real,1}=Float64[])

    if !isempty(momref)
        # assign new moment reference point
        ASWING.XYZREF[:,1] = momref
    end

    # get stability derivatives
    f_vba = zeros(Float64, 3, 3)
    m_vba = zeros(Float64, 3, 3)
    f_rot = zeros(Float64, 3, 3)
    m_rot = zeros(Float64, 3, 3)
    f_flap = zeros(Float64, 3, NFLPX)
    m_flap = zeros(Float64, 3, NFLPX)
    f_peng = zeros(Float64, 3, NENGX)
    m_peng = zeros(Float64, 3, NENGX)
    ccall((:dumpdm_, libaswing), Nothing,
    (Ref{Int32}, Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Float64},
    Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Int32}),
    6, f_vba, m_vba, f_rot, m_rot,
    f_flap, m_flap, f_peng, m_peng, ipnt)

    # limit output to defined flaps
    idx = findall(x->x!=0, ASWING.LSFLAP)
    f_flap = f_flap[:,idx]
    m_flap = m_flap[:,idx]
    # limit output to defined engines
    idx = findall(x->x!=0, ASWING.LSPENG)
    f_peng = f_peng[:,idx]
    m_peng = m_peng[:,idx]

  return StabilityDerivatives(f_vba, m_vba, f_rot, m_rot, f_flap, m_flap, f_peng, m_peng)
end

"""
    conclr()
    Clears all control laws
"""
function conclr()
    ASWING.CONLAW[1] = false
    ASWING.CONSET[1] = false
    ASWING.LDCON[:] .= false
    ASWING.LUCON[:] .= false
    ASWING.LQCON[:] .= false
    CONLAW.CU[:,:,:,:] .= 0.0
    CONLAW.CQ[:,:,:,:,:] .= 0.0
    CONLAW.DFREQ[:] .= 0.0
    CONLAW.DELTA[:] .= 0.0
    return nothing
end

"""
    conadd(i,j,ls,gfac,c)
    Sets up a proportional control law
"""
function conadd(i::Integer, j::Integer, ls::Integer, gfac::Real, c::Real)

    ASWING.CONLAW[1] = true
    if ls <= 0
        CONLAW.NUPAR1[i,j] = 2
        CONLAW.NUPAR2[i,j] = 2
        CONLAW.UPAR1[i,j,1:2] .= 0.0
        CONLAW.UPAR2[i,j,1:2] .= 0.0
        CONLAW.CU[i,j,1:2,1:2] .= c*gfac
        ASWING.LDCON[i] = true
        ASWING.LUCON[i] = true
    else
        CONLAW.NQPAR1[i,j,ls] = 2
        CONLAW.NQPAR2[i,j,ls] = 2
        CONLAW.QPAR1[i,j,ls,1:2] .= 0.0
        CONLAW.QPAR2[i,j,ls,1:2] .= 0.0
        CONLAW.CQ[i,j,ls,1:2,1:2] .= c*gfac
        ASWING.LDCON[i] = true
        ASWING.LQCON[j,ls] = true
    end
    ASWING.CONSET[1] = true

    return nothing
end

"""
    getintervals(ibeam::Integer)
    Gets number of panels between each break point
"""
function getintervals(ibeam::Integer)
    # get first beam index
    if ASWING.SBRK[1,ibeam] == ASWING.STIPL[ibeam]
        kbfrst = 1
    else
        kbfrst = 0
    end
    # get last beam index
    nbreak = ASWING.NBRK[ibeam]
    if ASWING.SBRK[nbreak,ibeam] == ASWING.STIPR[ibeam]
        kblast = nbreak
    else
        kblast = nbreak+1
    end
    # get spacing
    return ASWING.IBRK[kbfrst+1:kblast,ibeam].-ASWING.IBRK[kbfrst:kblast-1,ibeam].-1
end

"""
    nodeset(intervals::AbstractArray{<:AbstractArray{<:Integer, 1},1})

    Sets new number of nodes for ibeam.  A custom number of panels can be used
    between break points.
"""
function nodeset(intervals::AbstractArray{<:AbstractArray{<:Integer, 1},1})

    if length(intervals) != ASWING.NBEAM[1]
        error("length of intervals must equal ASWING.NBEAM")
    end

    for ibeam = 1:ASWING.NBEAM[1]
        # get first beam index
        if ASWING.SBRK[1,ibeam] == ASWING.STIPL[ibeam]
            kbfrst = 1
        else
            kbfrst = 0
        end

        # get last beam index
        nbreak = ASWING.NBRK[ibeam]
        if ASWING.SBRK[nbreak,ibeam] == ASWING.STIPR[ibeam]
            kblast = nbreak
        else
            kblast = nbreak+1
        end

        # assign new number of nodes to beam
        if !isempty(intervals[ibeam])
            nnodes = sum(intervals[ibeam])+length(intervals[ibeam])
            if nnodes < 2
                error("sum(intervals[$ibeam]) too small. nnodes must be >= 2")
            elseif nnodes > IIX
                error("sum(intervals[$ibeam]) too large. nnodes must be < IIX")
            end
            ASWING.II[ibeam] = nnodes
        end
        ccall((:iiset_, libaswing), Nothing, ())

        # assign new node distribution
        if !isempty(intervals[ibeam])
            THCALC.CUSTOMINTERVALS[1] = true
            THCALC.INTERVALS[kbfrst+1:kblast] = intervals[ibeam]
        else
            THCALC.CUSTOMINTERVALS[1] = false
        end
        ccall((:winit_, libaswing), Nothing, (Ref{Int32}, Ref{Int32}), ibeam, 0)
        ccall((:qinit_, libaswing), Nothing, (Ref{Int32}, Ref{Int32}, Ref{Float64}, Ref{Float64}),
            ASWING.IFRST[ibeam], ASWING.ILAST[ibeam], view(ASWING.Q,1:18,1:ASWING.IITOT[1],ASWING.IPOINT[1]), ASWING.Q0)
    end
    ccall((:spltset_, libaswing), Nothing, ())
    ccall((:micalc_, libaswing), Nothing, (Ref{Int32},), 1)
    ccall((:nnset_, libaswing), Nothing, ())

    # set to unconverged
    ASWING.LQINI[1:NPNTX] .= false
    ASWING.LCONV[1:NPNTX] .= false
    ASWING.LAINI[1:NPNTX] .= false

    # restore default node distribution settings
    THCALC.CUSTOMINTERVALS[1] = false
    return nothing
end

function getstabilityderivatives(ipnt::Integer=1)

    # no convergence failures yet
    fail = false

    # Get quasi-steady solution
    solution = solvesteady(ipnt)

    if solution.converged == false
        fail = true
    end

    # Retain parameters
    retainparam(ipnt)

    # Switch to anchored constraints
    setcons(Constraints("anchored"))

    # Get quasi-steady solution
    solution_qs = solvesteady(ipnt)

    # check if convergence failed
    if solution_qs.converged == false
        fail = true
    end

    # Get stability matrices using center of gravity as moment reference
    stabderivs_qs = stabilityderivatives(ipnt, momref=solution_qs.pnt_param.position_aircraft_cg)

    # rigidize aircraft at current deformed shape
    freeze(ipnt)

    # Get quasi-steady solution
    solution_rigid = solvesteady(ipnt)

    # check if convergence failed
    if solution_rigid.converged == false
        fail = true
    end

    # Get stability matrices using center of gravity as moment reference
    stabderivs_rigid = stabilityderivatives(ipnt, momref=solution_rigid.pnt_param.position_aircraft_cg)

    return stabderivs_qs, stabderivs_rigid, fail
end #getstabilityderivatives

"""
    getstaticmargin(ipnt::Integer=1)

    Returns quasi-steady static margin, rigidized static margin, and a solution
    failure flag.  Assumes constraints are set up for trimmed flight and follows
    the instructions in derivative_calc.txt. This function performs irreversible
    changes so you will need to re-input the geometry afterwards.
"""
function getstaticmargin(ipnt::Integer=1)

    stabderivs_qs, stabderivs_rigid, fail = getstabilityderivatives(ipnt)

    # calculate static margin
    sm_qs = -stabderivs_qs.Cma/stabderivs_qs.CLa*100

    # calculate static margin
    sm_rigid = -stabderivs_rigid.Cma/stabderivs_rigid.CLa*100

    return sm_qs, sm_rigid, fail
end #getstaticmargin


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

function cornerplot(sdata, vars=nothing;
                    quantiles=[0.25,0.50,0.75],
                    bandwidthx=100, bandwidthy=100,
                    kwargs...)
  # valid variables
  validvars = name.(variables(sdata))

  # plot all variables by default
  isnothing(vars) && (vars = validvars)
  @assert vars ⊆ validvars "invalid variable name"

  plts = []
  n = length(vars)
  for i in 1:n, j in 1:n
    xticks = i == n
    xguide = i == n ? vars[j] : ""
    yticks = i > 1 && j == 1
    yguide = i > 1 && j == 1 ? vars[i] : ""
    if i == j
      p = distplot1d(sdata, vars[i], quantiles=quantiles,
                     xticks=xticks, yticks=yticks,
                     xguide=xguide, yguide=yguide)
    elseif i > j
      p = distplot2d(sdata, vars[j], vars[i], quantiles=quantiles,
                     bandwidthx=bandwidthx, bandwidthy=bandwidthy,
                     xticks=xticks, yticks=yticks,
                     xguide=xguide, yguide=yguide)
    else
      p = RecipesBase.plot(framestyle=:none)
    end
    push!(plts, p)
  end

  RecipesBase.plot(plts...; layout=(n,n), kwargs...)
end


text:
"""
Lesson 4 - This is pretty Eventful

Close the window in response to an event.

Shows a smiley face with a transparent background over a tiled background of red and yellow text.

"""
module lesson4

using LibSDL2
const SDL = LibSDL2

sdlerror() = unsafe_string(SDL.GetError())

function checkptr(ptr)
    if ptr == C_NULL
        # The tutorial has us freeing and quitting and all sorts.
        # My philosophy: exit(1) is the simplest way of freeing
        # resources.
        error(sdlerror())
    else
        ptr
    end
end

function load_texture!(renderer, filename)
    surface = checkptr(SDL.LoadBMP(filename))
    texture = checkptr(SDL.CreateTextureFromSurface(renderer, surface))
    SDL.FreeSurface(surface)
    return texture
end

"""
    draw!(ren, tex::Ptr{SDL.Texture}, x, y, w, h)

Draw texture to renderer at position x, y, scaled to width & height w, h.

"""
function draw!(ren, tex::Ptr{SDL.Texture}, x, y, w, h)
    dst = SDL.Rect(x, y, w, h)
    # Don't know if I should be using Ref or pointer.
    # dst is a regular Julia struct, so it's lifetime is controlled by Julia,
    # so I'm going with Ref.
    SDL.RenderCopy(ren, tex, C_NULL, Ref(dst))
    return nothing
end

"""
    draw!(ren, tex::Ptr{SDL.Texture}, x, y)

Draw texture to renderer at position x, y, preserving texture's width and height

"""
function draw!(ren, tex::Ptr{SDL.Texture}, x, y)
    w, h = texture_wh(tex)
    draw!(ren, tex, x, y, w, h)
    return nothing
end

function query_texture(tex)
    # Is there a better way to do this?
    # Ref{Cint}(0) also works.
    format = UInt32[0]
    access = Cint[0]
    w = Cint[0]
    h = Cint[0]
    checkptr(SDL.QueryTexture(tex, format, access, w, h))
    return format[], access[], w[], h[]
end

"Return the width and height of `tex`"
function texture_wh(tex)
    return query_texture(tex)[3:4]
end

### Main

SDL.Event(::UndefInitializer) = SDL.Event(SDL.EventType(0), UInt32(0), ntuple(_ -> 0x00, 48))

function main()
    screen_width = 640
    screen_height = 480
    tile_size = 40

    checkptr(SDL.Init(SDL.INIT_VIDEO))

    window = checkptr(SDL.CreateWindow(
        "Lesson 4", 100, 100,
        screen_width, screen_height,
        SDL.WINDOW_SHOWN))

    # -1 is for the "index of the rendering driver to initialise, or -1
    # for the first one supporting the requested flags"
    renderer = checkptr(SDL.CreateRenderer(window, -1, SDL.RENDERER_ACCELERATED))

    background = load_texture!(renderer, "res/Lesson2/background.bmp")
    image = load_texture!(renderer, "res/Lesson2/image.bmp")

    quit = false
    while !quit
        event = SDL.Event(undef)
        while SDL.PollEvent(Ref(event)) != 0
            #= if event.type in (SDL.QUIT, SDL.KEYDOWN, SDL.MOUSEBUTTONDOWN) =#
                quit = true
            #= end =#
        end

        SDL.RenderClear(renderer)

        # Tiled background
        xtiles = screen_width ÷ tile_size
        ytiles = screen_height ÷ tile_size
        for i in 0:(xtiles * ytiles - 1)
            x = (i % xtiles) * tile_size
            y = (i ÷ xtiles) * tile_size
            draw!(renderer, background, x, y, tile_size, tile_size)
        end

        # Foreground
        iW, iH = texture_wh(image)
        draw!(renderer, image,
            screen_width ÷ 2 - iW ÷ 2,
            screen_height ÷ 2 - iH ÷ 2)

        SDL.RenderPresent(renderer)
        # Sleep 10 ms to avoid blowing the screen right off
        sleep(0.010)
    end

    SDL.DestroyTexture(background)
    SDL.DestroyTexture(image)
    SDL.DestroyRenderer(renderer)
    SDL.DestroyWindow(window)
    SDL.Quit()

    return true
end

end


text:
# number of visible keypoints = 13
# mask = [1 1 1 1;
#         1 0 1 1;
#         1 0 1 1;
#         1 1 1 0]
viewpoint_indices = [1, 2, 3, 4, 1, 3, 4, 1, 3, 4, 1, 2, 3]
point_indices     = [1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4]
indices = SBA.Indices(viewpoint_indices, point_indices)
N = length(viewpoint_indices)
A = randn(Float64, 2, 4, N)  # n_pose_params = 4
B = randn(Float64, 2, 3, N)  # n_point_params = 3
x_true = x_pred = randn(Float64, 2, N)

# n_rows(J) = 2 * n_visible_keypoints
#           = 2 * 16 = 26
# n_cols(J) = n _point_params * n_points + n_pose_params * n_viewpoints
#           = 3 * 4 + 4 * 4 = 12 + 16 = 28
# n_rows(J) < n_cols(J)
# J' * J can not be invertible. ArgumentError should be thrown
@test_throws ArgumentError SBA.sba(indices, x_true, x_pred, A, B)

# number of visible keypoints = 14
# mask = [1 1 1 1;
#         1 1 1 1;
#         1 1 1 0;
#         1 1 1 0]
viewpoint_indices = [1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 1, 2, 3]
point_indices     = [1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4]
indices = SBA.Indices(viewpoint_indices, point_indices)
N = length(viewpoint_indices)
A = randn(Float64, 2, 4, N)  # n_pose_params = 4
B = randn(Float64, 2, 3, N)  # n_point_params = 3
x_true = x_pred = randn(Float64, 2, N)

# n_rows(J) = 2 * n_visible_keypoints
#           = 2 * 14 = 28
# n_cols(J) = n _point_params * n_points + n_pose_params * n_viewpoints
#           = 3 * 4 + 4 * 4 = 12 + 16 = 28
# n_rows(J) = n_cols(J)
# J' * J can be invertible. Nothing should be thrown
SBA.sba(indices, x_true, x_pred, A, B)


text:
# Turing.jl version of model at https://github.com/stan-dev/blob/master/basic_estimators/normal_mixture.stan

using StatsFuns
using ForwardDiff: Dual

@model nmmodel(y) = begin
  N = length(y)

  theta ~ Uniform(0, 1)

  mu = tzeros(2) # mu is sampled by HMC
  for i = 1:2
    mu[i] ~ Normal(0, 10)
  end

  k = tzeros(Int64, N)
  for i = 1:N
    k[i] ~ Categorical([theta, 1.0 - theta])
  end

  for i = 1:N
    y[i] ~ Normal(mu[k[i]], 1.0)
  end

  # logtheta_p = map(yᵢ -> [log(theta) + logpdf(Normal(mu[1], 1.0), yᵢ), log(1 - theta) + logpdf(Normal(mu[2], 1.0), yᵢ)], y)
  # map!(logtheta_pᵢ -> logtheta_pᵢ - logsumexp(logtheta_pᵢ), logtheta_p)   # normalization
  # for i = 1:N
  #   k[i] ~ Categorical(exp.(logtheta_p[i]))
  # end
end


text:
const hp = 1  # huber parameter
const php = 1 # pseudo huber parameter
const sp = 3  # saturation parameter

function l1norm(resid)
    return @. sign(resid) * sqrt(abs(resid))
end

function l2norm(resid)
    return resid
end

function huber(resid)
    return [abs(r) < hp ? r/√(2) : sign(r) * sqrt(hp * (abs(r) - hp/2)) for r in resid]
end

function pseudo_huber(resid)
    return @. sign(resid) * php * sqrt(sqrt(1 + (resid/php)^2) - 1)
end

function saturate(resid)
    return [abs(r) < sp ? r : sp * sign(r) for r in resid]
end

function makesensible!(param::Vector{Float64})
    param[1] = (param[1] + 90) % 180
    if param[1] < 0
        param[1] += 180
    end
    param[1] = abs(param[1] - 90) # set inclination
    param[2] %= 180
    if param[2] < 0
        param[2] += 180
    end # set position angle
    for i in 1:length(param)÷3
        param[3*i] = abs(param[3*i])
    end
end

function fit(
    vis2::Vector{Vis2},
    model::Function;
    statistic::Function = l1norm,
    set_flux = 0,
    guess = nothing,
    parinfo = nothing,
    config = nothing,
    ritr = 0,
)
    lengths = length.(vis2)
    len = length(lengths)
    splits = zeros(Int64, len+1)
    for i in eachindex(lengths)
        splits[i+1] = splits[i] + lengths[i]
    end
    _uv = uv.(vis2)
    guessParam = if isnothing(guess)
        [45., 131.3, fill(0.5, 3*len)...]
    else
        s = 2 + 3 * len
        length(guess) ≠ s && error("`guess` should have $s parameters")
        guess
    end
    parinfo = CMPFit.Parinfo(length(guessParam))
    for i in 0:len-1
        guessParam[3+(3*i)] = 3.0
        parinfo[4+(3*i)].limited = (1, 1)
        parinfo[4+(3*i)].limits = (0, 1)
        parinfo[5+(3*i)].limited = (1, 1)
        parinfo[5+(3*i)].limits = (0, 1)
    end
    function cmpfit_callback(param::Vector{Float64})
        resid = vcat([(
            Measurements.value.(vis2[i].data) -
            model(_uv[i], param[1], param[2], param[3*i], param[3*i+1], param[3*i+2]).^2
        ) ./ Measurements.uncertainty.(vis2[i].data) for i in 1:len]...)
        return statistic(resid)
    end
    b = if ritr > 0
        function makefit()
            guessParam = [90, 180, vcat(fill([6, 1, 1], len)...)...] .* rand(3*len+2)
            parinfo = CMPFit.Parinfo(length(guessParam))
            for i in 0:len-1
                guessParam[3+(3*i)] *= 6.0
                parinfo[4+(3*i)].limited = (1, 1)
                parinfo[4+(3*i)].limits = (0, 1)
                parinfo[5+(3*i)].limited = (1, 1)
                parinfo[5+(3*i)].limits = (0, 1)
            end
            return cmpfit(cmpfit_callback, guessParam, parinfo=parinfo, config=config)
        end
        a = [makefit() for i in 1:ritr]
        a[argmin([a[i].bestnorm for i in 1:ritr])]
    else
        cmpfit(cmpfit_callback, guessParam, parinfo=parinfo, config=config)
    end
    makesensible!(b.param)
    return b
end

function fit(
    vis2::Vis2,
    model::Function;
    statistic::Function = l1norm,
    set_flux = 0,
    guess = nothing,
    parinfo = nothing,
    config = nothing,
    ritr = 0,
)
    fit(
        [vis2],
        model;
        statistic = statistic,
        set_flux = set_flux,
        guess = guess,
        parinfo = parinfo,
        config = config,
        ritr = ritr,
    )
end

function bootfit(
    vis2::Union{Vis2, Vector{Vis2}},
    model::Function;
    statistic::Function = l1norm,
    set_flux = 0,
    parinfo = nothing,
    config = nothing,
    ritr = 10,
    straps = 10,
)
    t0 = time()
    _vis2 = vis2 isa Vis2 ? [vis2] : vis2
    len = length(_vis2)
    sp = [Random.Sampler(Random.MersenneTwister, 1:length(v)) for v in _vis2]
    params = fill(NaN, 3*len+2, straps)
    for strap in 1:straps
        params[:, strap] .= fit(
            [_vis2[i][rand(sp[i], length(_vis2[i]))] for i in eachindex(_vis2)],
            model;
            statistic=statistic,
            set_flux=set_flux,
            parinfo=parinfo,
            config=config,
            ritr=ritr,
        ).param
        makesensible!(params[:, strap])
    end
    ret = fit(
        _vis2,
        model,
        statistic = statistic,
        set_flux = set_flux,
        parinfo = parinfo,
        config = config,
        ritr = ritr*10
    )
    ret.perror = [std(params[i,:]) for i in 1:size(params)[1]]
    ret.covar = collect(transpose(params))
    ret.elapsed = time() - t0
    return ret
end


text:
module test_h3_lib_regions

using Test
using H3.Lib

Lib.polyfill
Lib.maxPolyfillSize
Lib.h3SetToLinkedGeo
Lib.destroyLinkedPolygon

using .Lib: H3Index, GeoCoord, Geofence, GeoPolygon

# https://github.com/uber/h3/blob/master/src/apps/testapps/testPolyfill.c#L25

sfVerts = [GeoCoord(0.659966917655, -2.1364398519396),
           GeoCoord(0.6595011102219, -2.1359434279405),
           GeoCoord(0.6583348114025, -2.1354884206045),
           GeoCoord(0.6581220034068, -2.1382437718946),
           GeoCoord(0.6594479998527, -2.1384597563896),
           GeoCoord(0.6599990002976, -2.1376771158464)]
sfGeofence = Geofence(6, pointer(sfVerts))
sfGeoPolygon = GeoPolygon(sfGeofence, 0, C_NULL)
@test Lib.maxPolyfillSize(Ref(sfGeoPolygon), 9) == 3169

holeVerts = [GeoCoord(0.6595072188743, -2.1371053983433),
             GeoCoord(0.6591482046471, -2.1373141048153),
             GeoCoord(0.6592295020837, -2.1365222838402)]
holeGeofence = Geofence(3, pointer(holeVerts))
holeGeoPolygon = GeoPolygon(sfGeofence, 1, pointer_from_objref(Ref(holeGeofence)))
@test Lib.maxPolyfillSize(Ref(holeGeoPolygon), 9) == 3169

emptyVerts = [GeoCoord(0.659966917655, -2.1364398519394),
              GeoCoord(0.659966917655, -2.1364398519395),
              GeoCoord(0.659966917655, -2.1364398519396)]
emptyGeofence = Geofence(3, pointer(emptyVerts))
emptyGeoPolygon = GeoPolygon(emptyGeofence, 0, C_NULL)
@test Lib.maxPolyfillSize(Ref(emptyGeoPolygon), 9) == 1

verts = [GeoCoord(0.01, 0.01),
         GeoCoord(0.01, -0.01),
         GeoCoord(-0.01, -0.01),
         GeoCoord(-0.01, 0.01)]
fence = Geofence(4, pointer(verts))
polygon = GeoPolygon(fence, 0, C_NULL)

countActualHexagons(hexagons::Vector{H3Index})::Int = length(filter(!iszero, hexagons))

hexagons = Lib.polyfill(Ref(sfGeoPolygon), 9)
@test countActualHexagons(hexagons) == 1253

refpolygon = Lib.h3SetToLinkedGeo(H3Index[])
Lib.destroyLinkedPolygon(refpolygon)

end # module test_h3_lib_regions


text:
using ITensors
using ITensorMPOCompression
N=10; #10 sites
NNN=7; #Include up to 7th nearest neighbour interactions
sites = siteinds("S=1/2",N);
H=make_transIsing_MPO(sites,NNN);
pprint(H)
orthogonalize!(H,orth=right)
pprint(H)
truncate!(H,orth=left)
pprint(H)


text:
using JSON
Nproc=16
time_LL_grad=zeros(Nproc)
time_LL=zeros(Nproc)
time_PR=zeros(Nproc)

for i in 1:Nproc
    rmprocs(workers())
    addprocs(i)
    println("workers: ",length(workers()))
    @everywhere import Base.convert
    @everywhere using MAT
    @everywhere using ForwardDiff
    @everywhere using Optim
    @everywhere using Polynomials
    @everywhere using JSON
    @everywhere model_num="10"

    @everywhere Nstim = 80000

    @everywhere path="/home/hcli64/hcli64751/rats/fit_behaviour_julia/"
    #
    @everywhere include(path*"functions/functions_julia.jl")
    @everywhere include(path*"functions/parameter_functions.jl")
    @everywhere include(path*"functions/utils.jl")
    @everywhere rat_num = "1"
    @everywhere path_load=path*"scripts_rats/data/"
    @everywhere filename="processed_data_rat"*rat_num*".json"

    @everywhere data = Dict()
    open(path_load*filename, "r") do f
        global data
        dicttxt = readstring(f)  # file information to string
        data=JSON.parse(dicttxt)  # parse and transform data
    end

    stim_10 = data["stim_10"][1:Nstim,:]
    choices_10 = data["choices_10"][1:Nstim,:]

    @everywhere history_bias=zeros(Nstim)
    ###### check the gradient ##########

    @everywhere Nframes=10
    @everywhere Tframe=.05
    model_num="10"

    @everywhere args=[("c2",[1]),("c4",0),("st",[1,2]),("x0",0),("bias_d",0),("hbias_d",0),("sigma_a",[1])]
    @everywhere consts=[("c2",[2]),("c4",[1]),("st",[3,4]),("x0",[1,2]),("bias_d",[1]),("hbias_d",[1]),("sigma_a",0)]

    @everywhere z = [0,1,0,0,0,0,0,0]

    @everywhere lb=[-10.,0.0,0.0,0.1]
    @everywhere ub=[10.,10,10,10.]

    @everywhere model_f=Dict(:c2_f=>c2_urgency,:c4_f=>c4_const,:c6_f=>c6_zeros,:bias_d_f=>bias_d_const,
    :hbias_d_f=>history_bias_d_const,:st_f=>st_logistic_bias2,:x0_f=>initial_condition_bias_hbias,
    :ro_f=>read_out_perfect,:lambda=>0.0)

    x_ini=[ 0.6,4.72,6.38,1.52]
    LL = ComputeLL2(stim_10,choices_10,history_bias,Tframe, args, x_ini,consts,z,model_f)
    ll_t,grad_t=LL_function_grads(stim_10,choices_10,history_bias,Tframe, args, x_ini,consts,z,model_f)
    time_LL[i]=@elapsed ComputeLL2(stim_10,choices_10,history_bias,Tframe, args, x_ini,consts,z,model_f)
    time_LL_grad[i]=@elapsed ll_t,grad_t=LL_function_grads(stim_10,choices_10,history_bias,Tframe, args, x_ini,consts,z,model_f)
    param=make_dict(args,x_ini,consts,z)
    time_PR[i] = @elapsed PR=ProbRight2(stim_10,history_bias,Tframe,param,model_f)

end
@everywhere path_save="/home/hcli64/hcli64751/rats/fit_behaviour_julia/results/"
mkpath(path_save)
results = Dict("time_LL"=>time_LL,"time_LL_grad"=>time_LL_grad,"time_PR"=>time_PR,"workers"=>1:Nproc,"Nstim"=>Nstim)
# pass data as a json string (how it shall be displayed in a file)
stringdata = JSON.json(results)
# write the file with the stringdata variable information
open(path_save*"number_workers.json", "w") do f
        write(f, stringdata)
     end
println("done")

max_stars_repo_name: emmafdz/FinalDegreeProject

max_issues_repo_name: emmafdz/FinalDegreeProject

max_forks_repo_name: emmafdz/FinalDegreeProject


text:
macro workunit(label, ex)
    quote
        local start = time()
        local val = $(esc(ex))
        local stop = time()
        local label = $(esc(label))
        local tid = Threads.threadid()
        duration = stop-start
        lock(stdout)
        out = "[$label] thread = $tid, start = $start, stop = $stop, duration = $duration"
        println(stdout, out)
        unlock(stdout)
        val
    end
end

function plotgantt(out)
    lines = split(out, "\n")
    patt = r"\[(.*)\] thread = ([0-9]+), start = ([0-9\.e]+), stop = ([0-9\.e]+), duration = ([0-9\.e\-]+)"
    data::Vector{Dict{Symbol,Any}} = []
    for line in lines
        m = match(patt, line)
        isnothing(m) && continue
        label, thread, start, stop, duration = m.captures
        thread = parse(Int, thread)
        start, stop, duration = parse.(Float64, (start, stop, duration))
        push!(data, Dict(
                         :label => label,
                         :thread => thread,
                         :start => start,
                         :stop => stop,
                         :duration => duration,
                        ))
    end

    groups = sort(unique([d[:label] for d in data]))
    threads = sort(unique([d[:thread] for d in data]))
    gt0 = minimum([d[:start] for d in data])
    gt1 = maximum([d[:stop] for d in data])
    shapes::Vector{PlotlyJS.Shape} = []
    traces::Vector{GenericTrace} = []
    colors = PlotlyJS.colors.tab10
    usems, multiplier = false, 1
    if gt1-gt0 < 2
        usems, multiplier = true, 1e3
    end
    for (color,g) in zip(colors,groups)
        xs = Vector{Union{Float64,Missing}}()
        ys = Vector{Union{Float64,Missing}}()
        for d in data
            d[:label] != g && continue
            x0 = multiplier*(d[:start]-gt0)
            x1 = multiplier*(d[:stop]-gt0)
            y0 = d[:thread]-0.3
            y1 = d[:thread]+0.3
            # https://github.com/plotly/plotly.py/blob/0a833291f2b01b77cecf2a168936aa90e17acaed/packages/python/plotly/plotly/figure_factory/_gantt.py#L172
            append!(xs, [x0,x1,x1,x0,missing])
            append!(ys, [y0,y0,y1,y1,missing])
        end
        push!(traces, scatter(x=xs, y=ys, fill="toself", hoverinfo="name", mode="none", fillcolor=color, name=g))
    end

    totalwalltime = gt1-gt0
    markedwalltime = sum(d[:duration] for d in data)
    eff = (markedwalltime/length(threads)) / totalwalltime

    xaxis_attr = attr(title="walltime", ticksuffix=(usems ? "ms" : "s"), zeroline=false)
    yaxis_attr = attr(range=[0.5, length(threads)+0.5], autorange=false, zeroline=false, tickvals=threads, tickprefix="thread ")
    legend_attr = attr(orientation="h", yanchor="bottom", y=1.02, xanchor="right", x=1)
    title = "duration = $(round(totalwalltime; sigdigits=4))s, efficiency = $(round(100*eff; sigdigits=3))%"
    layout = Layout(;shapes=shapes,
                    yaxis=yaxis_attr,
                    xaxis=xaxis_attr,
                    legend=legend_attr,
                    title=title,
                    showlegend=true,
                   )
    p = plot(traces, layout)
    return p
end


text:
function eom_henon_triple(u, p, n)
    O = zeros(Float64, n + 3, 3)
    x₁, x₂, x₃ = (u...,)
    a, b, c = (p...,)

    # Propagate initial condition to the three first time steps.
    for i = 1:3
        O[i, 1] = x₁
        O[i, 2] = x₂
        O[i, 3] = x₃
    end
    for i = 4:n+3
        x₁1 = O[i-1, 1]
        x₁2 = O[i-2, 1]
        x₂1 = O[i-1, 2]
        x₂2 = O[i-2, 2]
        x₃1 = O[i-1, 3]
        x₃2 = O[i-2, 3]

        x₁new = a - x₁1^2 + b*x₁2
        x₂new = a - c*x₁1*x₂1 - (1 - c)*x₂1^2 + b*x₂2
        x₃new = a - c*x₂1*x₃1 - (1 - c)*x₃1^2 + b*x₃2

        O[i, 1] = x₁new
        O[i, 2] = x₂new
        O[i, 3] = x₃new
    end

    return O[4:end, :]
end

function henon_triple(uᵢ, a, b, c, n::Int, n_transient::Int)
    p = [a, b, c]
    o = eom_henon_triple(uᵢ, p, n + n_transient)
    x, y, z = o[n_transient+1:end, 1], o[n_transient+1:end, 2], o[n_transient+1:end, 3]
    return Dataset(x, y, z)
end

"""
    henon_triple(x, p, n) -> Function

Iterate a 3D discrete system consisting of coupled Henon maps where the coupling
is x1 -> x2 -> x3 [1]. This version allows for tweaking the parameters of the
equations.

The difference equations are:

```math
\\begin{aligned}
x_1(t+1) &= a - x_1(t)^2 + b x_1(t-2) \\
x_2(t+1) &= a - c x_1(t) x_2(t)- (1 - c) x_2(t)^2 + b x_2(t-1) \\
x_3(t+1) &= c x_2(t) x_3(t) - (1 - c) x_3(t)^2 + b x_3(t-1)
\\end{aligned}
```

Here ``c`` is the coupling constant. The system becomes completely synchronized
for ``c >= 0.7``.

# References
1. Papana, A., Kyrtsou, C., Kugiumtzis, D., & Diks, C. (2013). Simulation study of
direct causality measures in multivariate time series. Entropy, 15(7), 2635–2661.
"""
function henon_triple(;uᵢ = rand(3), a = 1.4, b = 0.3, c = 0.0,  n::Int = 100, n_transient::Int = 100)
    henon_triple(uᵢ, a, b, c, n, n_transient)
end


text:
using SparseArrays

function poisson_linear_system(ndofs)
        operator = spdiagm(
                -1 => 1.0 * ones(ndofs - 1),
                0 => -2.0 * ones(ndofs),
                1 => 1.0 * ones(ndofs - 1),
        )
        operator[1,1] = 1.0
        operator[1,2] = 0.0
        operator[ndofs,ndofs-1] = 0.0
        operator[ndofs,ndofs] = 1.0
        return operator
end

function solve_direct(operator, rhs, bcleft, bcright, stepsize)
    numpts = length(rhs)
    copyrhs = copy(rhs)
    copyrhs[1] = bcleft / stepsize^2
    copyrhs[numpts] = bcright / stepsize^2
    sol = stepsize^2 * (operator \ copyrhs)
    return sol
end

function solve_direct_all_rhs(operator, bcleft, bcright, stepsize)
    ndofs, ndofs = size(operator)
    rhs = diagm(ones(ndofs))[:, 2:ndofs-1]
    rhs[1, :] .= bcleft / stepsize^2
    rhs[ndofs, :] .= bcright / stepsize^2
    sol = stepsize^2 * (operator \ rhs)
    return sol
end

max_stars_repo_name: ArjunNarayanan/ExploreJacobi

max_issues_repo_name: ArjunNarayanan/ExploreJacobi

max_forks_repo_name: ArjunNarayanan/ExploreJacobi


text:
# # Packages
#
# Julia has over 2000 registered packages, making packages a huge part of the Julia ecosystem.
#
# Even so, the package ecosystem still has some growing to do. Notably, we have first class function calls  to other languages, providing excellent foreign function interfaces. We can easily call into python or R, for example, with `PyCall` or `Rcall`.
#
# This means that you don't have to wait until the Julia ecosystem is fully mature, and that moving to Julia doesn't mean you have to give up your favorite package/library from another language!
#
# To see all available packages, check out
#
# https://pkg.julialang.org/
# or
# https://juliaobserver.com/
#
# For now, let's learn how to use a package.

#-

# The first time you use a package on a given Julia installation, you need to use the package manager to explicitly add it:

using Pkg
Pkg.add("Example")

# Every time you use Julia (start a new session at the REPL, or open a notebook for the first time, for example), you load the package with the `using` keyword

using Example

# In the source code of `Example.jl` at
# https://github.com/JuliaLang/Example.jl/blob/master/src/Example.jl
# we see the following function declared
#
# ```
# hello(who::String) = "Hello, $who"
# ```
#
# Having loaded `Example`, we should now be able to call `hello`

hello("it's me. I was wondering if after all these years you'd like to meet.")

# Now let's play with the Colors package

Pkg.add("Colors")

#-

using Colors

# Let's create a palette of 100 different colors

palette = distinguishable_colors(100)

# and then we can create a randomly checkered matrix using the `rand` command

rand(palette, 3, 3)

# In the next notebook, we'll use a new package to plot datasets.

#-

# # What happened?
#
# Ask Julia for its status

Pkg.status()

# Examine the Project.toml
#
# The Project.toml holds all high-level dependencies and a unique identifier.

# Examine the Manifest.toml
#
# The Manifest holds _all_ dependencies and their exact versions.

# # Environments
#
# This is all happening in the global environment right now — but we can choose
# to use the environment in this folder. This guarantees reproducibility!

# # Creating your own package
#
# Best way to do this is through the PkgTemplates meta-package

Pkg.add("PkgTemplates")

# Now use the generate_interactive() function at the REPL to create your new package

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
# Test all IC specification methods
using DelimitedFiles

N = 3
t0 = 7257.93115525 - 7300.0
H = [N,ones(Int64, N - 1)...]

function get_elements_array_ic()
    fname = "elements.txt"
    elements = readdlm(fname, ',', comments=true)[1:N,:]
    ic = ElementsIC(t0, H, elements)
    return ic
end

function get_cartesian_array_ic()
    fname = "coordinates.txt"
    coords = readdlm(fname, ',', comments=true)[1:N,:]
    ic = CartesianIC(t0, H[1], coords)
    return ic
end

function get_elements_file_ic()
    fname = "elements.txt"
    ic = ElementsIC(t0, H, fname)
    return ic
end

function get_cartesian_file_ic()
    fname = "coordinates.txt"
    ic = CartesianIC(t0, H[1], fname)
    return ic
end

function get_elements_ic()
    fname = "elements.txt"
    elements = readdlm(fname, ',', comments=true)[1:N,:]
    elems = Elements{Float64}[]
    for i in 1:N
        push!(elems, Elements(elements[i,:]...))
    end
    ic = ElementsIC(t0, H, elems)
end

@testset "Initial Conditions" begin
    @testset "Elements" begin
        ic_file = get_elements_file_ic()
        ic_array = get_elements_array_ic()
        ic_elems = get_elements_ic()

        @test ic_file.elements == ic_array.elements
        @test ic_elems.elements == ic_array.elements

        @test ic_file.amat == ic_array.amat
        @test ic_elems.amat == ic_array.amat
    end

    @testset "Cartesian" begin
        ic_file = get_cartesian_file_ic()
        ic_array = get_cartesian_array_ic()

        @test ic_file.x == ic_array.x
        @test ic_file.v == ic_array.v
        @test ic_file.m == ic_array.m
    end
end


text:
using BNP_WMReg_Joint

using Statistics
using Plotly

using PDMats
using StatsBase

using BSON
using BSON: @load

println(pwd())

## load
using LinearAlgebra # along with all others above

whichdat = "pinkSalmon_first30_log_n25_K5_SigXdiag_VSlocal_gaminitall_snr5_1121191"

whichdat = "Ricker_single_..."
whichdat = "Ricker_single_lognormal_..."
whichdat = "Ricker_twolag_lognormal_..."

BSON.@load "postsimB/postsim_$(whichdat).bson" sims model init mesg

### posterior checks
println(model.state.llik)

println(model.state.iter)
println(model.state.accpt / float(model.state.iter))
println(counts(model.state.S, 1:model.H))

println(exp.(model.state.lω))
println(model.state.α)

iter_use = 1:length(sims) # all simulations
iter_use = 1:1000 # first thousand
iter_use = range(1, stop=length(sims), step=5) # by 5
length(collect(iter_use))

plot( [ sims[i][:llik] for i in iter_use ] )

plot( [ sims[i][:n_occup] for i in iter_use ] )

println( round.( mean( hcat([ sims[i][:Scounts] for i in iter_use ]...), dims=2 ), digits=1 ) )
println( round.( mean( hcat([ sims[i][:Scounts] .> 0 for i in iter_use ]...), dims=2 ), digits=1 ) )
println( sum( mean( hcat([ sims[i][:Scounts] .> 0 for i in iter_use ]...), dims=2 ) .> 0.0 ) )

h = 1
plot( [ sims[i][:Scounts][h] for i in iter_use ] )

### lag selection

## global
k = 1
plot( [ 1*sims[i][:γ][k] for i in iter_use ] )
println( mean( [ 1*sims[i][:γ][k] for i in iter_use, k = 1:model.K ], dims=1 ) )

## local
γoccup = permutedims( hcat( [ vec(sims[ii][:γ_occupied]) for ii in iter_use ]... ) )
pm_γoccup = mean( γoccup, dims=1 )[1,:]
println(pm_γoccup)

b = 0.05
γ_alt = [ ( sims[iter_use[ii]][:γ][h,k] && abs(sims[iter_use[ii]][:β_y][h,k]) > b )  for ii = 1:length(iter_use), h = 1:model.H, k = 1:model.K ]
γ_alt_occup = [ sims[iter_use[ii]][:Scounts]'γ_alt[ii,:,k] / model.n  for ii = 1:length(iter_use), k = 1:model.K ]
pm_γ_alt_occup = mean( γ_alt_occup, dims=1 )[1,:]
println(pm_γ_alt_occup)

sims_gam_glob = [ any( sims[i][:γ][ 1:maximum( findall( sims[i][:Scounts] .> 0 ) ) , k] ) for i in iter_use, k = 1:model.K ]
pm_gam_glob = mean(sims_gam_glob, dims=1)[1,:]

k = 1
plot( sims_gam_glob[:,k] )

mean( [ sims[i][:π_γ] for i in iter_use ] )

k = 3
plot( [ sims[i][:π_γ][k] for i in iter_use ] )

mean( [ 1.0*sims[i][:ξ] for i in iter_use ] )

k = 2
plot( [ sims[i][:ξ][k] for i in iter_use ] )

sims_gam_w = [ sum( exp.(sims[i][:lω][sims[i][:γ][:,k]]) ) for i=1:length(sims), k = 1:model.K ]

k = 3
plot( sims_gam_w[iter_use, k] )
plot( γoccup[iter_use, k] )

h = 1
k = 2
println( round.( mean( hcat([ sims[i][:Scounts] for i in iter_use ]...), dims=2 ), digits=1 ) )

plot( [ 1*sims[i][:γ][h,k] for i in iter_use ] )
for h = 1:model.H
    println(mean( hcat([ 1*sims[i][:γ][h,:] for i in iter_use ]...), dims=2 ))
end

maxScount_indx = [ findmax(sims[i][:Scounts])[2] for i in 1:length(sims) ]
plot( [ 1*sims[i][:γ][maxScount_indx[i],k] for i in iter_use ] )

### component kernel parameters

model.state.δ_x
h = 1
k = 1
plot( [sims[i][:δ_x][h,k] for i in iter_use] )
Statistics.median([sims[i][:δ_x][h,k] for i in iter_use])

model.state.μ_x
h = 2
k = 1
plot([sims[i][:μ_x][h,k] for i in iter_use])
mean([sims[i][:μ_x][h,k] for i in iter_use])
std(sims.μ_x, dims=1)[1,:,:]

model.state.β_x
h = 1
kk = 2
k = 1
plot( [sims[i][:β_x][kk][h,k] for i in iter_use] )

model.state.μ_y
h = 2
plot([sims[i][:μ_y][h] for i in iter_use])

model.state.β_y
h = 1
k = 1
plot( [sims[i][:β_y][h,k] for i in iter_use] )
Statistics.median([sims[i][:β_y][h,k] for i in iter_use])
Statistics.mean([sims[i][:β_y][h,k] for i in iter_use])
println(Statistics.std([sims[i][:β_y][h,k] for i in iter_use]))

Statistics.mean([sims[i][:β_y][h,k] for i in iter_use])
println(Statistics.std([sims[i][:β_y][h,k] for i in iter_use]))

## prior covariance for beta_y
inv(model.state.Λ0star_ηy.mat) * model.state.δ_y[h]
model.prior.s0_δy_s0 * model.prior.Λ0star_ηy_S0.mat

model.state.δ_y
h = 2
plot( [sims[i][:δ_y][h] for i in iter_use] )
median([sims[i][:δ_y][h] for i in iter_use])

# intercepts
h = 3

# local lag selection
plot( [sims[i][:μ_y][h] + sum( sims[i][:γ][h,:] .* sims[i][:β_y][h,:] .* sims[i][:μ_x][h,:]) for i in iter_use ] )

# global lag selection
plot( [sims[i][:μ_y][h] + sum( sims[i][:γ] .* sims[i][:β_y][h,:] .* sims[i][:μ_x][h,:]) for i in iter_use ] )
Statistics.mean([sims[i][:μ_y][h] + sum( sims[i][:γ] .* sims[i][:β_y][h,:] .* sims[i][:μ_x][h,:]) for i in iter_use ])
println(Statistics.std([sims[i][:μ_y][h] + sum( sims[i][:γ] .* sims[i][:β_y][h,:] .* sims[i][:μ_x][h,:]) for i in iter_use ]))

h = 4
plot(([sims[i][:lω][h] for i in iter_use]))
plot(exp.([sims[i][:lω][h] for i in iter_use]))

vec(mean( exp.( permutedims(hcat([sims[i][:lω] for i in iter_use]...)) ) , dims=1 ))
vec(std( exp.( permutedims(hcat([sims[i][:lω] for i in iter_use]...)) ) , dims=1 ))
counts(model.state.S, 1:model.H)

plot(exp.([sims[i][:lω][h] for i in iter_use]))
plot(exp.([sims[i][:lω][model.H] for i in iter_use]))

plot([sims[i][:α] for i in iter_use])

### Functionals

include("Rfuns.jl")

npred = 25

minx = minimum(model.X)
maxx = maximum(model.X)
ranx = maxx - minx

X_pred_vals = collect(range(minx - 0.01*ranx, length=npred,
                        stop=(maxx + 0.01*ranx)))

## γ_use is a boolean vector specifying which lags we want to examine
γ_use = deepcopy(model.state.γ) # global only
γ_use = [false, true, false, false, false]
γ_use = pm_γoccup .> 0.11 # local
γ_use = pm_gam_glob .> 0.4 # local
println(1*γ_use)

iter_use = sort(findall([ sims[ii][:γ] == γ_use for ii in 1:length(sims)])) # global only
iter_use = range(1, stop=length(sims), step=10)
iter_use = 1:(length(sims))

length(iter_use)

γindx_use = findall(γ_use)
Kuse = sum(γ_use)
X_pred00 = rcopy(xpndgrid([X_pred_vals for k = 1:Kuse]...))

X_pred0 = hcat([ fill(mean(model.X[:,k]), size(X_pred00,1)) for k = 1:model.K ]...)
X_pred0 = hcat([ rand(size(X_pred00,1)) .* (maximum(model.X[:,k]) - minimum(model.X[:,k])) .+ minimum(model.X[:,k])  for k = 1:model.K ]...) # draw uniformly over the range

X_pred0[:,findall(γ_use)] = deepcopy(X_pred00)
X_pred = deepcopy(X_pred0)

# X_pred = deepcopy(model.X) # for residual analysis

## use posterior samples to recreate mixture weights
ldw = ldensweight_mat(X_pred, sims[iter_use], model.state.γδc)
dw = exp.(ldw)
sum(dw[1,1,:])

i = 1
ℓ = 2
h = 2
plot([scatter(x=X_pred[:,ℓ], y=dw[i,:,h], mode="scatter")])

## examine weight function
mean_dw = mean(dw, dims=1)[1,:,:]

println( round.( mean( hcat([ sims[i][:Scounts] for i in iter_use ]...), dims=2 ), digits=1 ) )

h = 2
plot( [ sims[i][:Scounts][h] for i in iter_use ] )

mean_dw_mat = Matrix(reshape(mean_dw[:,h], fill(npred,Kuse)...)) # weight "surface"
trace1 = surface(Dict(
  :x => X_pred_vals,
  :y => X_pred_vals,
  :z => mean_dw_mat,
  # :z => q025_Q_mat_plot,
  :colorscale => "Viridis",
  :opacity => 0.7,
  :showscale => false,
  :type => "surface"
))
data = [trace1]
plot(data)

ℓ = 1
h = 3
plot([scatter(x=X_pred[:,ℓ], y=mean_dw[:,h], mode="scatter")]) # 1D

plot([scatter(x=X_pred[:,ℓ], y=mean_dw[1,:,h], mode="scatter")])
plot(scatter(x=X[:,1], y=y, mode="markers"))

## use weight function to create transition mean function (surface)

Ey = getEy(X_pred, dw, sims[iter_use])

# ii = 500
# plot([scatter(x=X_pred[:,1], y=Ey[ii,:], mode="scatter")])

mean_Ey = mean(Ey, dims=1)[1,:]

q025_Ey = [quantile(Ey[:,j], 0.025) for j = 1:size(Ey,2)]
q50_Ey = [quantile(Ey[:,j], 0.5) for j = 1:size(Ey,2)]
q975_Ey = [quantile(Ey[:,j], 0.975) for j = 1:size(Ey,2)]

resid = model.y - mean_Ey # ONLY use this if X_pred is the X matrix from the model
plotR(mean_Ey, resid, ylab="residual", xlab="posterior mean predicted value")

## transition quantile function (surface)

quant = 0.2
Q = getQuant(quant, X_pred, dw, sims[iter_use])

mean_Q = mean(Q, dims=1)[1,:]

q025_Q = [quantile(Q[:,j], 0.025) for j = 1:size(Q,2)]
q975_Q = [quantile(Q[:,j], 0.975) for j = 1:size(Q,2)]

### plot transition surfaces

### two dimensions

mean_Ey_mat = Matrix(reshape(mean_Ey, fill(npred,Kuse)...))
maximum(mean_Ey_mat)
minimum(mean_Ey_mat)
mean_Ey_mat_plot = [ abs(mean_Ey_mat[i,j]) < maximum(abs.([minx, maxx])) + 0.6*ranx ? mean_Ey_mat[i,j] : NaN for i = 1:npred, j=1:npred ]
mean_Ey_mat_plot = [ abs(mean_Ey_mat[i,j]) < maxx + 0.3*ranx ? mean_Ey_mat[i,j] : NaN for i = 1:npred, j=1:npred ]
# plot([scatter(x=X_pred[:,1], y=mean_Ey, mode="scatter")])

q025_Ey_mat = Matrix(reshape(q025_Ey, fill(npred,Kuse)...))
q025_Ey_mat_plot = [ abs(q025_Ey_mat[i,j]) < maxx + 0.1*ranx ? q025_Ey_mat[i,j] : NaN for i = 1:npred, j=1:npred ]

q975_Ey_mat = Matrix(reshape(q975_Ey, fill(npred,Kuse)...))
q975_Ey_mat_plot = [ abs(q975_Ey_mat[i,j]) < maxx + 0.1*ranx ? q975_Ey_mat[i,j] : NaN for i = 1:npred, j=1:npred ]

## quantile
mean_Q_mat = Matrix(reshape(mean_Q, fill(npred,Kuse)...))
maximum(mean_Q_mat)
mean_Q_mat_plot = [ abs(mean_Q_mat[i,j]) < maxx + 0.3*ranx ? mean_Q_mat[i,j] : NaN for i = 1:npred, j=1:npred ]

q025_Q_mat = Matrix(reshape(q025_Q, fill(npred,Kuse)...))
q025_Q_mat_plot = [ abs(q025_Q_mat[i,j]) < maxx + 0.1*ranx ? q025_Q_mat[i,j] : NaN for i = 1:npred, j=1:npred ]
q975_Q_mat = Matrix(reshape(q975_Q, fill(npred,Kuse)...))
q975_Q_mat_plot = [ abs(q975_Q_mat[i,j]) < maxx + 0.1*ranx ? q975_Q_mat[i,j] : NaN for i = 1:npred, j=1:npred ]

trace1 = surface(Dict(
  :x => X_pred_vals,
  :y => X_pred_vals,
  :z => q025_Ey_mat_plot,
  # :z => q025_Q_mat_plot, # uncomment to view quantile surfaces
  :colorscale => "Viridis",
  :opacity => 0.7,
  :showscale => false,
  :type => "surface"
))

trace2 = surface(Dict(
  :x => X_pred_vals,
  :y => X_pred_vals,
  :z => mean_Ey_mat_plot,
  # :z => mean_Q_mat_plot, # uncomment to view quantile surfaces
  :colorscale => "Viridis",
  :opacity => 0.9,
  :showscale => false,
  :type => "surface"
))

trace3 = surface(Dict(
  :x => X_pred_vals,
  :y => X_pred_vals,
  :z => q975_Ey_mat_plot,
  # :z => q975_Q_mat_plot, # uncomment to view quantile surfaces
  :colorscale => "Viridis",
  :opacity => 0.7,
  :showscale => false,
  :type => "surface"
))

trace0 = scatter3d(Dict(
    :x => model.X[:,γindx_use[1]],
    :y => model.X[:,γindx_use[2]],
    :z => model.y,
    :marker => Dict(
        :color => "rgba(0,0,0,0.70)",
        :size => 2
    ),
    :mode => "markers"
))

data = [trace2, trace0] # mean surface and data only
data = [trace1, trace2, trace3, trace0] # include intervals

layout = Layout(Dict(
  :hovermode => "closest",
  :margin => Dict(
    :r => 10,
    :t => 25,
    :b => 40,
    :l => 60
  ),
  :scene => Dict( :aspectmode => "manual", # doesn't seem to be working
    :xaxis => Dict(:title => "y[t-$(γindx_use[1])]",
                   # :range => [minx - 0.1*ranx, maxx + 0.1*ranx]),
                   :autorange => true),
    :yaxis => Dict(:title => "y[t-$(γindx_use[2])]",
                   # :range => [minx - 0.1*ranx, maxx + 0.1*ranx]),
                   :autorange => true),
    :zaxis => Dict(:title => "y[t]",
                   # :range => [minx - 0.1*ranx, maxx + 0.1*ranx])
                   :autorange => true),
  ),
  :showlegend => false
))

p = plot(data, layout)

### one dimension
size(Ey)
plotR(X_pred_vals, mean_Ey, ylab="", xlab="",
    ylim=[minimum(q025_Ey), maximum(q975_Ey)], lwd=2, typ="l")
plotR(X_pred_vals, mean_Ey, ylab="", xlab="",
    ylim=[minimum(y), maximum(y)], lwd=2, typ="l")

linesR(X_pred_vals, q025_Ey, lty=2)
linesR(X_pred_vals, q975_Ey, lty=2)
pointsR(model.X[:,γindx_use[1]], model.y)

R"title(main='Pink Salmon')"
R"title(ylab='log y[t]')"
R"title(xlab='log y[t-2] ')"
# linesR(X_pred_vals, mean_Q, lwd=2, typ="l", col="blue")
# linesR(X_pred_vals, p05_Q, lty=2, col="blue")
# linesR(X_pred_vals, p95_Q, lty=2, col="blue")

## individual iteration mean functions
for i in 1:5
    linesR(X_pred_vals, Ey[i,:], col=i, lwd=1)
end

R"par(mar=c(4,4,1,1)+0.1)"
plotR(X_pred_vals, mean_Ey, ylab="", xlab="",
    # ylim=[minimum(q025_Ey), maximum(q975_Ey)],
    ylim = [-1.0, 13.0], # for pink salmon
    lwd=2, typ="l")
polygonR(vcat(X_pred_vals, reverse(X_pred_vals)), vcat(q025_Ey, reverse(q975_Ey)), col="gray90", border="NA")
pointsR(model.X[:,γindx_use], model.y)
linesR(X_pred_vals, mean_Ey, lwd=2)

R"title(main='Pink Salmon');
    title(ylab='log y[t]');
    title(xlab='log y[t-2]')"
R"abline(0,1, lty=3)"

R"dev.off()"

plotR(model.y, typ="o", ylab="") # time series
plotR(model.X[:,2], model.y, ylab="", xlab="") # lag scatter plot

### Density estimation
## This section requires several variables assigned at the beginning of the 'Functionals' section.

using Distributions

n_y = 25
y_grid = collect(range(minx - 0.3*ranx, length=n_y,
                        stop=(maxx + 0.3*ranx)))
# y_grid = collect(range(minimum(model.y)-2.0, length=n_y, stop=maximum(model.y)+2.0))

n_star_Py = 4 # number of x locations at which to estimate the transition density

## single lag
n_star_Py = 3
X_star_Py = deepcopy(X_pred[1:n_star_Py,:])
X_star_Py[1, [γindx_use[1]] ] = [ 1.0 ]
X_star_Py[2, [γindx_use[1]] ] = [ 3.0 ]
X_star_Py[3, [γindx_use[1]] ] = [ 5.0 ]

## two lag
X_star_Py = deepcopy(X_pred[1:n_star_Py,:])
X_star_Py[1, [γindx_use[1], γindx_use[2]] ] = [ 1.0, 2.0 ]
X_star_Py[2, [γindx_use[1], γindx_use[2]] ] = [ 2.0, 1.0 ]
X_star_Py[3, [γindx_use[1], γindx_use[2]] ] = [ 2.5, 2.5 ]
X_star_Py[4, [γindx_use[1], γindx_use[2]] ] = [ 3.0, 2.5 ]

# iter_use = 1:length(sims)
iter_use = range(1, stop=length(sims), step=5)

X_star_Py

## construct weight functions and transition densities from posterior samples
ldwPy = ldensweight_mat(X_star_Py, sims[iter_use])
lPy, Ey_0 = getlogdens_EY(X_star_Py, y_grid, ldwPy, sims[iter_use])

Py = exp.(lPy)
size(Py)

mean_Py = permutedims(mean(Py, dims=1)[1,:,:])
q025_Py = [quantile(Py[:,j,i], 0.025) for i = 1:n_y, j = 1:n_star_Py]
# p50_Py = [quantile(Py[:,j,i], 0.5) for i = 1:n_y, j = 1:n_star_Py]
q975_Py = [quantile(Py[:,j,i], 0.975) for i = 1:n_y, j = 1:n_star_Py]

@rput mean_Py
R"library(lattice); levelplot(mean_Py)" # quick, discretized estimate of transtition density

xind = 1 # select row of n_star_Py for which to evaluate density
println(X_star_Py[xind,:])

## Transition densities for individual MCMC iterations
ii = 1 # posterior sample
j = 1 # row of n_star_Py
plotR(y_grid, Py[ii,j,:], typ="l", ylab="", xlab="")
for ii in range(1, stop=1000, step=3)
    linesR(y_grid, Py[ii,j,:], typ="l", ylab="", xlab="", col=ii)
end
R"dev.off()"

R"par(mfrow=c(2,2))"
for xind in 1:n_star_Py
    plotR(y_grid, q025_Py[:,xind],
        main="y[t-1] = $(round(X_star_Py[xind,γindx_use[1]], digits=2))",
        ylab="density", xlab="y[t] (minutes)", typ="l", ylim=[0.0, maximum(q975_Py[:,xind])], lty=2)
    linesR(y_grid, q975_Py[:,xind], lty=2)
    linesR(y_grid, mean_Py[:,xind], lty=1, lwd=2)
end
R"dev.off()"


text:
using Test

include("flatland.jl")
include("typesystem.jl")