text:
##
# Parameter layout: (κ₀, τ₀, ϕ₀, κ₁, τ₁, ϕ₁, ν)
# The optimization: note that this model has data-determined nonlinear contraints
#
#   0.0 > ϕ₀ - τ₀*iΨ(-1.0/κ₀) - minimum(data) **IF κ₀ < 0.0**
#   0.0 < ϕ₁ + τ₁*iΨ(-1.0/κ₁) - maximum(data) **IF κ₁ < 0.0**
#   κ₀/ν > -0.5  # two derivatives at endpoint
#   κ₁/ν > -0.5  # two derivatives at endpoint
#
# The way that Ipopt does handles these is by putting box constraints on the
# nonlinear argument functions. The functions constr1 and constr2 give those
# constraints in mathematical language. The code representation for Ipopt is
# much harder to interpret.
#
##

# This function is constrained to be NEGATIVE. It is not differentiable at κ=0,
# and so the jacobian should break there.
#
# In math:
# 0.0 > ϕ₀ - τ₀*iΨ(-1.0/κ₀) - minimum(data) **IF κ₀ < 0.0**
function constr1(p, mindat)
  (κ₀, τ₀, ϕ₀, κ₁, τ₁, ϕ₁, ν) = p
  κ₀ ≥ 0.0 && return -1.0 # If κ₀>0, return a very feasible number.
  ϕ₀ - τ₀*iΨ(-1.0/κ₀) - mindat
end

# This function is constrained to be POSITIVE. It is not differentiable at κ=0,
# and so the jacobian should break there.
#
# In math:
# 0.0 < ϕ₁ + τ₁*iΨ(-1.0/κ₁) - maximum(data) **IF κ₁ < 0.0**
function constr2(p, maxdat)
  (κ₀, τ₀, ϕ₀, κ₁, τ₁, ϕ₁, ν) = p
  κ₁ ≥ 0.0 && return 1.0 # If κ₁>0, return a very feasible number.
  ϕ₁ + τ₁*iΨ(-1.0/κ₁) - maxdat
end

# Tail indices:
@inline constr3(p) = p[1]/p[7]
@inline constr4(p) = p[4]/p[7]

function Distributions.fit_mle(::Type{BulkAndTailsDist}, data::AbstractArray;
                               init=[1.0, 1.0, 0.0, 1.0, 1.0, 0.0, 2.0],
                               tol=1.0e-5, maxit::Int=2000, print_level=0)

  # Box constraints (lowerb and upperb) and inequality constraints for the way
  # that I've implemented the nonlinear constraints.
  lowerb = [-1e22,  0.0,  -1e22, -1e22,  0.0,  -1e22, 0.0]
  upperb = [1e22,  1e22,   1e22,  1e22, 1e22,   1e22, 1e22]
  cons_low = [-1e22, 0.0, -0.5, -0.5]
  cons_upp = [0.0,  1e22, 1e22, 1e22]

  # objective, gradient, and Hessian, but with a cache to save some intermediate
  # results and cut down on repeated calculations.
  obj_cache  = Dict{UInt64, DiffResults.MutableDiffResult}()
  # objective function:
  obj = p->-sum(x->batslogpdf(x, p), data)
  # a method to add keys to the cache, which effectively computes the hessian
  # and gradient of obj and stores them for later access.
  addkey! = p -> begin
    haskey(obj_cache, hash(p)) && return
    res = DiffResults.HessianResult(p)
    ForwardDiff.hessian!(res, obj, p)
    obj_cache[hash(p)] = res
  end
  # a gradient for obj, which either gets it from the cache or adds it to the
  # cache and then gets it from the cache.
  grad! = (p,g) -> begin
    if haskey(obj_cache, hash(p))
      g .= DiffResults.gradient(obj_cache[hash(p)])
      return nothing
    end
    addkey!(p)
    grad!(p,g)
  end
  # analogous functionality for the Hessian.
  hess = p -> begin
    haskey(obj_cache, hash(p)) && return DiffResults.hessian(obj_cache[hash(p)])
    addkey!(p)
    hess(p)
  end

  # Writing the constraints in several forms, all of which are necessary for the
  # various ways that we use the constraints. Note that these constraints also
  # are involved in what ultimately gets given to Ipopt as the "hessian",
  # because the objective function gets a scaled constraint function added to it.
  #
  # (i):  one argument form
  # (ii): a tuple of one-argument functions
  # (iii): a one-argument vector-valued function
  _constr1(p) = constr1(p, minimum(data)) # (i)
  _constr2(p) = constr2(p, maximum(data)) # (i)
  constrf     = (_constr1, _constr2, constr3, constr4) # (ii)
  constr(p)   = [cons_f(p) for cons_f in constrf] # (iii)

  # Write the constraint in the final format that ipopt actually wants.
  ipopt_constr(p,g) = (g .= constr(p))

  # A properly reshaped jacobian (not very pretty sorry).
  g_jac(p) = vec(Matrix(ForwardDiff.jacobian(constr, p)'))

  # Prepare the final hessian and constraint jacobian functions.
  nconstr  = length(constrf)
  constrfh = [p->ForwardDiff.hessian(cj, p) for cj in constrf]
  ipopt_hess = (x,r,c,o,l,v) -> ipopt_hessian(x,r,c,o,l,v,hess,constrfh,nconstr)
  ipopt_jac_constr = (x,r,c,v) -> ipopt_constr_jac(x,r,c,v,g_jac,nconstr)

  # Set up the problem:
  prob = CreateIpoptProblem(
                       length(init), # number of parameters
                       lowerb, upperb, # box for parameter constraints
                       length(constrf), # number of nonlinear constraints
                       cons_low, cons_upp, # box for nonlinear constraint fxns
                       length(init)*length(constrf), # size of CONSTRAINT JACOBIAN
                       div(length(init)*(length(init)+1), 2), # size of ltri of HESSIAN
                       obj,
                       ipopt_constr,
                       grad!,
                       ipopt_jac_constr,
                       ipopt_hess)
  AddIpoptStrOption(prob, "sb", "yes")
  AddIpoptNumOption(prob, "tol", tol)
  AddIpoptIntOption(prob, "max_iter", maxit)
  AddIpoptIntOption(prob, "print_level", Int(print_level))
  prob.x = deepcopy(init)
  status = IpoptSolve(prob)

  # If a solution was found, return it. Otherwise throw an error.
  iszero(status) && return (mle = BulkAndTailsDist(prob.x), hess = Matrix(hess(prob.x)), nll = obj(prob.x))
  throw(error("Optimization failed with Ipopt return code $status."))
end

# for R users:
function fitbats(data; kwargs...)
  (mle, hes, nll) = fit_mle(BulkAndTailsDist, data; kwargs...)
  [getfield(mle, f) for f in fieldnames(BulkAndTailsDist)], hes, nll
end

## version for covariates

function bats_negloglikelihood_covariates(
  data::AbstractArray,
  τmatrix::AbstractArray,
  ϕmatrix::AbstractArray,
  parms
)

  m = size(data, 1)
  nϕ = size(ϕmatrix, 2)
  nτ = size(τmatrix, 2)

  @assert ((isequal(m,size(τmatrix,1))) && isequal(m,size(ϕmatrix,1))) "Dimension mismatch"

  #bit of messy indexing
  κ₀ = parms[1]
  τ₀vec = parms[2:(nτ+1)]
  ϕ₀vec = parms[(nτ+2):(nτ+1+nϕ)]
  κ₁ = parms[nτ+2+nϕ]
  τ₁vec = parms[(nτ+3+nϕ):(2*nτ+2+nϕ)]
  ϕ₁vec = parms[(2*nτ+3+nϕ):(2*nτ+2+2*nϕ)]
  ν = parms[2*nτ+3+2*nϕ]

  eltp = eltype(parms)
  output = zero(eltp)

  @inbounds for i in 1:m
    ϕ₀ = zero(eltp)
    ϕ₁ = zero(eltp)
    logτ₀ = zero(eltp)
    logτ₁ = zero(eltp)

    @inbounds for j in 1:nϕ
      ϕ₀ += ϕ₀vec[j] * ϕmatrix[i, j]
      ϕ₁ += ϕ₁vec[j] * ϕmatrix[i, j]
    end

    @inbounds for j in 1:nτ
      logτ₀ += τ₀vec[j] * τmatrix[i, j]
      logτ₁ += τ₁vec[j] * τmatrix[i, j]
    end

    output -= logpdf(BulkAndTailsDist(κ₀,exp(logτ₀),ϕ₀,κ₁,exp(logτ₁),ϕ₁,ν),data[i])
  end

  return output

end

# fit
function fit_bats_mle_covariates(data::AbstractArray, τmatrix::AbstractArray, ϕmatrix::AbstractArray, init;
  tol=1.0e-5, maxit::Int=1000, print_level=0)

  # Box constraints (lowerb and upperb) and inequality constraints for the way
  # that I've implemented the nonlinear constraints.
  lowerb = [repeat([-1e22], length(init)-1); 0.0]
  upperb = repeat([1e22], length(init))
  cons_low = [-0.5, -0.5]
  cons_upp = [1e22, 1e22]

  # objective function:
  obj = p -> bats_negloglikelihood_covariates(
    data,
    τmatrix,
    ϕmatrix,
    p
  )

  #cache functionality seems to be slower here, so just calling ForwardDiff directly
  grad!(p, g) = ForwardDiff.gradient!(g, obj, p)
  hess(p) = ForwardDiff.hessian(obj, p)

  # Writing the constraints in several forms, all of which are necessary for the
  # various ways that we use the constraints. Note that these constraints also
  # are involved in what ultimately gets given to Ipopt as the "hessian",
  # because the objective function gets a scaled constraint function added to it.
  #
  # (i):  one argument form
  # (ii): a tuple of one-argument functions
  # (iii): a one-argument vector-valued function
  nϕ = size(ϕmatrix,2)
  nτ = size(τmatrix,2)
  lenparms = length(init)
  function constr1_covariate(p, lenparms) #(i)
    p[1]/p[lenparms]
  end
  function constr2_covariate(p, nτ, nϕ, lenparms) #(i)
    p[2+nτ+nϕ]/p[lenparms]
  end
  _constr1(p) = constr1_covariate(p, lenparms) # (i)
  _constr2(p) = constr2_covariate(p, nτ, nϕ, lenparms) # (i)
  constrf     = (_constr1, _constr2) # (ii)
  constr(p)   = [cons_f(p) for cons_f in constrf] # (iii)

  # Write the constraint in the final format that ipopt actually wants.
  ipopt_constr(p,g) = (g .= constr(p))

  # A properly reshaped jacobian (not very pretty sorry).
  g_jac(p) = vec(Matrix(ForwardDiff.jacobian(constr, p)'))
  #g_jac(p) = vec(Matrix(wrapper(:jacobian, p)'))

  # Prepare the final hessian and constraint jacobian functions.
  nconstr  = length(constrf)
  constrfh = [p->ForwardDiff.hessian(cj, p) for cj in constrf]
  ipopt_hess = (x,r,c,o,l,v) -> ipopt_hessian(x,r,c,o,l,v,hess,constrfh,nconstr)
  ipopt_jac_constr = (x,r,c,v) -> ipopt_constr_jac(x,r,c,v,g_jac,nconstr)

  # Set up the problem:
  prob = CreateIpoptProblem(
  length(init), # number of parameters
  lowerb, upperb, # box for parameter constraints
  length(constrf), # number of nonlinear constraints
  cons_low, cons_upp, # box for nonlinear constraint fxns
  length(init)*length(constrf), # size of CONSTRAINT JACOBIAN
  div(length(init)*(length(init)+1), 2), # size of ltri of HESSIAN
  obj,
  ipopt_constr,
  grad!,
  ipopt_jac_constr,
  ipopt_hess)
  AddIpoptStrOption(prob, "sb", "yes")
  AddIpoptNumOption(prob, "tol", tol)
  AddIpoptIntOption(prob, "max_iter", maxit)
  AddIpoptIntOption(prob, "print_level", Int(print_level))
  prob.x = deepcopy(init)
  status = IpoptSolve(prob)

  return (status = status, nll = obj(prob.x), κ₀ = prob.x[1], τ₀ = prob.x[2:(nτ+1)], ϕ₀ = prob.x[(nτ+2):(nτ+1+nϕ)], κ₁ = prob.x[nτ+2+nϕ], τ₁ = prob.x[(nτ+3+nϕ):(2*nτ+2+nϕ)], ϕ₁ = prob.x[(2*nτ+3+nϕ):(2*nτ+2+2*nϕ)], ν = prob.x[2*nτ+3+2*nϕ])

end

#for R users, making the output a bit more accessible as a list without unicode characters
function fitbats_covariates(data, τmatrix, ϕmatrix, init; kwargs...)
  x = fit_bats_mle_covariates(data, τmatrix, ϕmatrix, init; kwargs...)
  return Dict("status" => x.status, "nll" => x.nll, "kappa0" => x.κ₀, "tau0" => x.τ₀, "phi0" => x.ϕ₀, "kappa1" => x.κ₁, "tau1" => x.τ₁, "phi1" => x.ϕ₁, "nu" => x.ν)
end


text:
mutable struct A
  a1::Int64
  a2::Int64
end

mutable struct B
  b1::Int64
  b2::A
end

mutable struct C
  c1::A
  c2::B
end

a1 = A(2, 3)
a2 = A(3, 4)
b = B(5, a1)
c = C(a2, b)

println(c)

ext: jl

max_stars_repo_name: judy-vscode/judy

max_issues_repo_name: judy-vscode/judy

max_forks_repo_name: judy-vscode/judy


text:
# Autogenerated wrapper script for GeographicLibWrapper_jll for x86_64-linux-musl-cxx03-julia_version+1.6.0
export libGeographicLibWrapper

using GeographicLib_jll
using libcxxwrap_julia_jll
JLLWrappers.@generate_wrapper_header("GeographicLibWrapper")
JLLWrappers.@declare_library_product(libGeographicLibWrapper, "libGeographicLibWrapper.so")
function __init__()
    JLLWrappers.@generate_init_header(GeographicLib_jll, libcxxwrap_julia_jll)
    JLLWrappers.@init_library_product(
        libGeographicLibWrapper,
        "lib/libGeographicLibWrapper.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
struct ElementContext
  default :: Default
  units :: Units
  autoname :: Int
end
ElementContext() = ElementContext(Default(), Units(), 0)

const ContextDict = Dict{Element,ElementContext}

struct Context
  dict :: ContextDict
  firstunits :: Units
  title :: String
end
function Context(element::Element)
  recursioncheck(element)
  Context(contextdict(element),firstunits(element),title(element))
end

const NamedElements = Union{Node,Segment,UniformPlane}

function firstunits(element::Element)
  result = Units()
  for e in traverseTree(element)
    if firstunits_(e) != nothing
      result =  firstunits_(e)
      break
    end
  end
  return result
end
firstunits_(::Element) = nothing
firstunits_(units::Units) = units

function contextdict(element::Element)
  cd = ContextDict()
  previouselementcontext = ElementContext()
  for e in traverseTree(element)
    previouselementcontext = appendelementcontext!(cd,previouselementcontext,e)
  end
  return cd
end

function appendelementcontext!(cd::ContextDict, pec::ElementContext, x::Default)
  newelementcontext = ElementContext(x, pec.units, pec.autoname)
  cd[x] = newelementcontext
end
function appendelementcontext!(cd::ContextDict, pec::ElementContext, x::Units)
  newelementcontext = ElementContext(pec.default, x, pec.autoname)
end
function appendelementcontext_!(cd::ContextDict, pec::ElementContext, x::NamedElements)
  if ~haskey(cd,x)
    if x.name == Symbol("")
      newelementcontext = ElementContext(pec.default, pec.units, pec.autoname+1)
    else
      newelementcontext = pec
    end
    cd[x] = newelementcontext
    return newelementcontext
  else
    return pec
  end
end
function appendelementcontext!(cd::ContextDict, pec::ElementContext, x::Node)
  appendelementcontext_!(cd,pec,x)
end
function appendelementcontext!(cd::ContextDict, pec::ElementContext, x::Segment)
  pec = appendelementcontext_!(cd,pec,x)
  pec = appendelementcontext_!(cd,pec,x.node1)
  pec = appendelementcontext_!(cd,pec,x.node2)# store default wxyz here?
  if x.wxwywz.isdefault
    node1_xyz1 = x.node1.xyz1
    node2_xyz1 = x.node2.xyz1 * segmentnodescaleratio(x,cd)
    v1 = node1_xyz1[1:3] - node2_xyz1[1:3]
    normalize!(v1)
    v2 = [0.0, 0.0, 1.0]
    w = cross(v1,v2)
    if norm(w) < 1e-10  # close to 0
      v2 = [0.0, 1.0, 0.0]
      w = cross(v1,v2)
    end
    x.wxwywz.xyz[1:3] = normalize(w)
  end
  return pec
end
function appendelementcontext!(cd::ContextDict, pec::ElementContext, x::UniformPlane)
  pec = appendelementcontext_!(cd,pec,x)
  for node in x.nodes
    pec = appendelementcontext_!(cd,pec,node)
  end
  pec
end
function appendelementcontext!(cd::ContextDict, pec::ElementContext, x::Equiv)
  for node in x.nodes
    pec = appendelementcontext_!(cd,pec,node)
  end
  pec
end
function appendelementcontext!(cd::ContextDict, pec::ElementContext, x::External)
  pec = appendelementcontext_!(cd,pec,x.node1)
  pec = appendelementcontext_!(cd,pec,x.node2)
end
appendelementcontext!(::ContextDict, pec::ElementContext, ::Element) = pec

# methods using Context return values in the first unit of their group

const tometers = Dict("km"  =>1e3,
                      "m"   =>1.0,
                      "cm"  =>1e-2,
                      "mm"  =>1e-3,
                      "um"  =>1e-6,
                      "in"  =>2.54e-2,
                      "mils"=>2.54e-5,
                      ""    =>1.0)  # meters is default

function scaletofirstunits(x::Element, context::Context)
  tometers[context.dict[x].units.unitname] /
    tometers[context.firstunits.unitname]
end

function xyz1(node::Node, scale::Float64)
  result = Array{Float64}(undef,4)
  result[1:3] = node.xyz1[1:3]*scale
  result[4] = 1.0
  return result
end

function xyz1(node::Node, context::Context)
  scale = scaletofirstunits(node,context)
  xyz1(node, scale)
end
function segmentnodescaleratio(segment::Segment, cd::ContextDict)
  if haskey(cd,segment.node1)
    scale1 = tometers[cd[segment.node1].units.unitname]
  else
    scale1 = tometers[cd[segment].units.unitname]
  end
  if haskey(cd,segment.node2)
    scale2 = tometers[cd[segment.node2].units.unitname]
  else
    scale2 = tometers[cd[segment].units.unitname]
  end
  return scale2/scale1
end
function nodes_xyz1(segment::Segment, context::Context)
  scale = scaletofirstunits(segment,context)
  if haskey(context.dict,segment.node1)
    node1_xyz1 = xyz1(segment.node1,context)
  else
    node1_xyz1 = xyz1(segment.node1,scale)
  end
  if haskey(context.dict,segment.node2)
    node2_xyz1 = xyz1(segment.node2,context)
  else
    node2_xyz1 = xyz1(segment.node2,scale)
  end
  return (node1_xyz1, node2_xyz1)
end

function width_height(segment::Segment, context::Context)
  scale = scaletofirstunits(segment,context)
  if isnan(segment.wh.w) || isnan(segment.wh.h)
    scaledefault = scaletofirstunits(context.dict[segment].default,context)
  end
  if isnan(segment.wh.w)
    width = context.dict[segment].default.wh.w * scaledefault
  else
    width = segment.wh.w * scale
  end
  if isnan(segment.wh.h)
    height = context.dict[segment].default.wh.h * scaledefault
  else
    height = segment.wh.h * scale
  end
  return (width, height)
end

#wxyz(segment::Segment, ::Context) = segment.wxwywz.xyz
# todo: wxyz should return a normalized vector perp to length vector in context
# just a warning?
# todoing...

function wxyz(segment::Segment, context::Context)
  wxyz_ = segment.wxwywz.xyz
  if ~segment.wxwywz.isdefault
    (n1,n2) = nodes_xyz1(segment,context)
    lxyz = n2[1:3] - n1[1:3]
    if dot(lxyz,segment.wxwywz.xyz)>1e-9
      wxyz_ = cross(cross(lxyz,segment.wxwywz.xyz),lxyz)
      normalize!(wxyz_)
      @warn "wx,wy,wz not perpendicular to segment length" is=segment.wxwywz.xyz should_be=wxyz_
    end
  end
  return wxyz_
end

function corners_xyz1_thick(uniformplane::UniformPlane, context::Context)
  scale = scaletofirstunits(uniformplane,context)
  c1 = Array{Float64}(undef,4)
  c2 = similar(c1)
  c3 = similar(c1)
  c1[1:3] = uniformplane.corner1[1:3] * scale
  c2[1:3] = uniformplane.corner2[1:3] * scale
  c3[1:3] = uniformplane.corner3[1:3] * scale
  c1[4] = 1.0
  c2[4] = 1.0
  c3[4] = 1.0
  thick = uniformplane.thick * scale
  (c1,c2,c3,thick)
end

function nodes_xyz1(uniformplane::UniformPlane, context::Context)
  if length(uniformplane.nodes)!=0
    scale = scaletofirstunits(uniformplane,context)
  end
  f(i) = xyz1(uniformplane.nodes[i],scale) #plane nodes should not be in context
  ntuple(f, length(uniformplane.nodes))
end

function title(e1::Element)
  for e2 in traverseTree(e1)
    return title_(e2)
  end
  return ""
end
title_(::Element) = ""
title_(x::Comment) = x.text

recursioncheck(element::Element) = recursioncheck!(element, Dict{Group,Group}())
recursioncheck!(::Element,::Dict{Group,Group}) = nothing
function recursioncheck!(group::Group, recursioncheckdict::Dict{Group,Group})
  if haskey(recursioncheckdict,group)
    throw(ErrorException("recursion not allowed"))
  end
  recursioncheckdict[group] = group
  for element in group.elements
    recursioncheck!(element, recursioncheckdict)
  end
  delete!(recursioncheckdict, group)
  return nothing
end


text:
module ESDLTools
using Distributed
import ..ESDL: ESDLdir
export mypermutedims!, freshworkermodule, passobj, @everywhereelsem,
toRange, getiperm, CItimes, CIdiv, @loadOrGenerate, PickAxisArray
struct PickAxisArray{T,N,AT<:AbstractArray,P,NCOL}
    parent::AT
end

function PickAxisArray(parent, indmask; ncol=0)
    @assert sum(indmask)+ncol==ndims(parent)
    f  = findall(indmask)
    PickAxisArray{eltype(parent),length(indmask),typeof(parent),(f...,),ncol}(parent)
end
indmask(p::PickAxisArray{<:Any,<:Any,<:Any,i,<:Any}) where i = i
getncol(p::PickAxisArray{<:Any,<:Any,<:Any,<:Any,NCOL}) where NCOL = NCOL
function Base.view(p::PickAxisArray, i...)
    inew = map(j->i[j],indmask(p))
    cols = ntuple(_ -> Colon(), getncol(p))
    view(p.parent,cols...,inew...)
end
function Base.getindex(p::PickAxisArray, i...)
    inew = map(j->i[j],indmask(p))
    cols = ntuple(_ -> Colon(), getncol(p))
    getindex(p.parent,cols...,inew...)
end
Base.getindex(p::PickAxisArray,i::CartesianIndex) = p[i.I...]
function newparent(a::PickAxisArray{T,N,AT,P,NCOL},parent) where {T,N,AT,P,NCOL}
  eltype(parent) == T || error("Types do not match")
  size(parent) == size(a.parent) || error("Sizes do not match")
  PickAxisArray{T,N,AT,P,NCOL}(parent)
end

function getiperm(perm)
    iperm = Array{Int}(undef,length(perm))
    for i = 1:length(perm)
        iperm[perm[i]] = i
    end
    return ntuple(i->iperm[i],length(iperm))
end

using Base.Cartesian
@generated function mypermutedims!(dest::AbstractArray{T,N},src::AbstractArray{S,N},perm::Type{Q}) where {Q,T,S,N}
    ind1=ntuple(i->Symbol("i_",i),N)
    ind2=ntuple(i->Symbol("i_",perm.parameters[1].parameters[1][i]),N)
    ex1=Expr(:ref,:src,ind1...)
    ex2=Expr(:ref,:dest,ind2...)
    quote
        @nloops $N i src begin
            $ex2=$ex1
        end
    end
end

@generated function CIdiv(index1::CartesianIndex{N}, index2::CartesianIndex{N}) where N
    I = index1
    args = [:(Base.div(index1[$d],index2[$d])) for d = 1:N]
    :($I($(args...)))
end
@generated function CItimes(index1::CartesianIndex{N}, index2::CartesianIndex{N}) where N
    I = index1
    args = [:(.*(index1[$d],index2[$d])) for d = 1:N]
    :($I($(args...)))
end

@generated function Base.getindex(t::NTuple{N},p::NTuple{N,Int}) where N
    :(@ntuple $N d->t[p[d]])
end

toRange(r::CartesianIndices)=map(colon,r.start.I,r.stop.I)
toRange(c1::CartesianIndex,c2::CartesianIndex)=map(colon,c1.I,c2.I)

function passobj(src::Int, target::Vector{Int}, nm::Symbol;
                 from_mod=Main, to_mod=Main)
    r = RemoteChannel(src)
    @spawnat(src, put!(r, getfield(from_mod, nm)))
    @sync for to in target
        @spawnat(to, Core.eval(to_mod, Expr(:(=), nm, fetch(r))))
    end
    nothing
end

function passobj(src::Int, target::Int, nm::Symbol; from_mod=Main, to_mod=Main)
    passobj(src, [target], nm; from_mod=from_mod, to_mod=to_mod)
end

function passobj(src::Int, target, nms::Vector{Symbol};
                 from_mod=Main, to_mod=Main)
    for nm in nms
        passobj(src, target, nm; from_mod=from_mod, to_mod=to_mod)
    end
end

function sendto(p::Int; args...)
    for (nm, val) in args
        @spawnat(p, Core.eval(Main, Expr(:(=), nm, val)))
    end
end

function sendto(ps::Vector{Int}; args...)
    for p in ps
        sendto(p; args...)
    end
end

getfrom(p::Int, nm::Symbol; mod=Main) = fetch(@spawnat(p, getfield(mod, nm)))

function freshworkermodule()
    in(:PMDATMODULE,names(Main)) || Core.eval(Main,:(module PMDATMODULE
        using ESDL
        using Distributed
    end))
    Core.eval(Main,quote
      using Distributed
      rs=Future[]
      for pid in workers()
        n=remotecall_fetch(()->in(:PMDATMODULE,names(Main)),pid)
        if !n
          r1=remotecall(()->(Core.eval(Main,:(using ESDL, Distributed));nothing),pid)
          r2=remotecall(()->(Core.eval(Main,:(module PMDATMODULE
          using ESDL
          using Distributed
        end));nothing),pid)
          push!(rs,r1)
          push!(rs,r2)
        end
      end
      [wait(r) for r in rs]
  end)

  nothing
end

import Distributed: extract_imports, remotecall_eval
macro everywhereelsem(ex)
  modu=:(Main.PMDATMODULE)
  procs = :(workers())
  imps = [Expr(:import, m) for m in extract_imports(ex)]
  return quote
    p = $(esc(procs))
    if p!=[1]
      $(isempty(imps) ? nothing : Expr(:toplevel, imps...)) # run imports locally first
      let ex = $(Expr(:quote, ex)), procs = p, mo=$(esc(modu))
        remotecall_eval(mo, procs, ex)
      end
    end
  end
end

"""
    macro loadOrGenerate(x...,expression)

Takes a list of `variablename=>"Storage Name"` pairs. Checks if all datasets can be found
on disk and loads them. If not, the datasets will be regenerated by evaluating the given expression.

To force recalculation, call `ESDL.recalculate(true)` before evaluating the macro.

### Example

The following lines will check if cubes with the names "Filled" and "Normalized"
exist on disk, load them and assign the variable names `cube_filled` and `cube_norm`.
If the datasets to not exist on disk, they are generated and saved under the given names.

````julia
@loadOrGenerate cube_filled=>"Filled" cube_norm=>"Normalized" begin
cube_filled = mapCube(gapFillMSC,d)
cube_norm   = mapCube(normalize_TS,d)
end

````
"""
macro loadOrGenerate(x...)
  code=x[end]
  x=x[1:end-1]
  x2=map(x) do i
    isa(i,Symbol) ? (i,string(i)) : (i.head==:call && i.args[1]==:(=>)) ? (i.args[2],i.args[3]) : error("Wrong Argument type")
  end
  xnames=map(i->i[2],x2)
  loadEx=map(x2) do i
    :($(i[1]) = loadCube($(i[2])))
  end
  loadEx=Expr(:block,loadEx...)
  saveEx=map(x2) do i
    :(saveCube($(i[1]),$(i[2])))
  end
  saveEx=Expr(:block,saveEx...)
  rmEx=map(x2) do i
    :(rmCube($(i[2])))
  end
  rmEx=Expr(:block,rmEx...)
  esc(quote
    if !ESDL.recalculate() && all(i->isdir(joinpath(ESDLdir(),i)),$xnames)
      $loadEx
    else
      $rmEx
      $code
      $saveEx
    end
  end)
end
end


text:
using CSFML, CSFML.LibCSFML

include("audio.jl")

const scaleFactor = 4.0

sfColor(val::UInt32) = CSFML.sfColor(((val >> 16) & 0xFF) % UInt8, ((val >> 8) & 0xFF) % UInt8, (val & 0xFF) % UInt8, 0xFF)

function window(name::String, engine)

  game = engine.loadgame(name)

  texture = sfTexture_create(256, 240)
  @assert texture != C_NULL

  sprite = sfSprite_create()
  sfSprite_setTexture(sprite, texture, sfTrue)
  sfSprite_setScale(sprite, sfVector2f(scaleFactor, scaleFactor))

  image = sfImage_create(256, 240)

  mode = sfVideoMode(Int(256 * scaleFactor), Int(240 * scaleFactor), 32)

  window = sfRenderWindow_create(mode, name, sfResize | sfClose, C_NULL)
  @assert window != C_NULL
  sfWindow_setFramerateLimit(window, 60)
  updatescreen!(window, texture, image, sprite, engine, game)

  previousSound = C_NULL
  previousBuffer = C_NULL

  event_ref = Ref(sfEvent(sfEvtClosed, ntuple(_ -> UInt32(0), 20)))

  starttime = time()
  totalticks = 0
  while Bool(sfRenderWindow_isOpen(window))
    # process events
    while Bool(sfRenderWindow_pollEvent(window, event_ref))
      # close window : exit
      if event_ref.x.type == sfEvtClosed
        sfRenderWindow_close(window)
      end
    end
    # Read input
    input = readinput()
    engine.setbuttons1!(game, input)

    # Step a single frame
    totalticks += engine.stepframe!(game)

    # Play audio.
    samples = engine.audiosamples(game)
    currentBuffer = sfSoundBuffer_createFromSamples(samples, length(samples), 1, 44100)
    currentSound = sfSound_create()
    sfSound_setBuffer(currentSound, currentBuffer)
    sfSound_play(currentSound)

    # Render the screen
    updatescreen!(window, texture, image, sprite, engine, game)

    # Swap out sound buffers
    sfSound_destroy(previousSound)
    sfSoundBuffer_destroy(previousBuffer)
    previousSound = currentSound
    previousBuffer = currentBuffer
  end

  sfSound_destroy(previousSound)
  sfSoundBuffer_destroy(previousBuffer)
  sfSprite_destroy(sprite)
  sfTexture_destroy(texture)
  sfRenderWindow_destroy(window)
end

function updatescreen!(window, texture, image, sprite, engine, game)
    s = engine.screen(game)
    for i = 1:256, j = 1:240
      color = sfColor(s[j, i].color)
      sfImage_setPixel(image, i-1, j-1, color)
    end
    sfTexture_updateFromImage(texture, image, 0, 0)
    # draw the sprite
    sfRenderWindow_drawSprite(window, sprite, C_NULL)
    # update the window
    sfRenderWindow_display(window)
end

function readinput()
  input = 0x00
  if Bool(sfKeyboard_isKeyPressed(sfKeyC))
    input |= 0x01
  end
  if Bool(sfKeyboard_isKeyPressed(sfKeyX))
    input |= 0x02
  end
  if Bool(sfKeyboard_isKeyPressed(sfKeyTab))
    input |= 0x04
  end
  if Bool(sfKeyboard_isKeyPressed(sfKeyEnter))
    input |= 0x08
  end
  if Bool(sfKeyboard_isKeyPressed(sfKeyUp))
    input |= 0x10
  end
  if Bool(sfKeyboard_isKeyPressed(sfKeyDown))
    input |= 0x20
  end
  if Bool(sfKeyboard_isKeyPressed(sfKeyLeft))
    input |= 0x40
  end
  if Bool(sfKeyboard_isKeyPressed(sfKeyRight))
    input |= 0x80
  end
  input
end


text:
macro custom(expr)
    (Meta.isexpr(expr, :call, 3) && expr.args[1] === :~) || error("incorrect macro usage")
    quote
        $(esc(expr.args[2])) = 0.0
    end
end

macro mymodel1(ex)
    # check if expression was modified by the DynamicPPL "compiler"
    if ex == :(y ~ Uniform())
        return esc(:(x ~ Normal()))
    else
        return esc(:(z ~ Exponential()))
    end
end

struct MyModelStruct{T}
    x::T
end
Base.:~(x, y::MyModelStruct) = y.x
macro mymodel2(ex)
    # check if expression was modified by the DynamicPPL "compiler"
    if ex == :(y ~ Uniform())
        # Just returns 42
        return :(4 ~ MyModelStruct(42))
    else
        return :(return -1)
    end
end

# Used to test sampling of immutable types.
struct MyCoolStruct{T}
    a::T
end

@testset "compiler.jl" begin
    @testset "model macro" begin
        @model function testmodel_comp(x, y)
            s ~ InverseGamma(2, 3)
            m ~ Normal(0, sqrt(s))

            x ~ Normal(m, sqrt(s))
            y ~ Normal(m, sqrt(s))

            return x, y
        end
        @test length(methods(testmodel_comp)) == 2
        testmodel_comp(1.0, 1.2)

        # check if drawing from the prior works
        @model function testmodel01(x=missing)
            x ~ Normal()
            return x
        end
        @test length(methods(testmodel01)) == 3
        f0_mm = testmodel01()
        @test mean(f0_mm() for _ in 1:1000) ≈ 0.0 atol = 0.1

        # Test #544
        @model function testmodel02(x=missing)
            if x === missing
                x = Vector{Float64}(undef, 2)
            end
            x[1] ~ Normal()
            x[2] ~ Normal()
            return x
        end
        @test length(methods(testmodel02)) == 3
        f0_mm = testmodel02()
        @test all(x -> isapprox(x, 0; atol=0.1), mean(f0_mm() for _ in 1:1000))

        @model function testmodel03(x=missing)
            x ~ Bernoulli(0.5)
            return x
        end
        f01_mm = testmodel03()
        @test length(methods(testmodel03)) == 3
        @test mean(f01_mm() for _ in 1:1000) ≈ 0.5 atol = 0.1

        # test if we get the correct return values
        @model function testmodel1(x1, x2)
            s ~ InverseGamma(2, 3)
            m ~ Normal(0, sqrt(s))

            x1 ~ Normal(m, sqrt(s))
            x2 ~ Normal(m, sqrt(s))

            return x1, x2
        end
        @test length(methods(testmodel1)) == 2
        f1_mm = testmodel1(1.0, 10.0)
        @test f1_mm() == (1, 10)

        # alternatives with keyword arguments
        testmodel1kw(; x1, x2) = testmodel1(x1, x2)
        f1_mm = testmodel1kw(; x1=1.0, x2=10.0)
        @test f1_mm() == (1, 10)

        @model function testmodel2(; x1, x2)
            s ~ InverseGamma(2, 3)
            m ~ Normal(0, sqrt(s))

            x1 ~ Normal(m, sqrt(s))
            x2 ~ Normal(m, sqrt(s))

            return x1, x2
        end
        @test length(methods(testmodel2)) == 2
        f1_mm = testmodel2(; x1=1.0, x2=10.0)
        @test f1_mm() == (1, 10)

        @info "Testing the compiler's ability to catch bad models..."

        # Test for assertions in observe statements.
        @model function brokentestmodel_observe1(x1, x2)
            s ~ InverseGamma(2, 3)
            m ~ Normal(0, sqrt(s))

            x1 ~ Normal(m, sqrt(s))
            x2 ~ x1 + 2

            return x1, x2
        end

        btest = brokentestmodel_observe1(1.0, 2.0)
        @test_throws ArgumentError btest()

        @model function brokentestmodel_observe2(x)
            s ~ InverseGamma(2, 3)
            m ~ Normal(0, sqrt(s))

            x = Vector{Float64}(undef, 2)
            x ~ [Normal(m, sqrt(s)), 2.0]

            return x
        end

        btest = brokentestmodel_observe2([1.0, 2.0])
        @test_throws ArgumentError btest()

        # Test for assertions in assume statements.
        @model function brokentestmodel_assume1()
            s ~ InverseGamma(2, 3)
            m ~ Normal(0, sqrt(s))

            x1 ~ Normal(m, sqrt(s))
            x2 ~ x1 + 2

            return x1, x2
        end

        btest = brokentestmodel_assume1()
        @test_throws ArgumentError btest()

        @model function brokentestmodel_assume2()
            s ~ InverseGamma(2, 3)
            m ~ Normal(0, sqrt(s))

            x = Vector{Float64}(undef, 2)
            x ~ [Normal(m, sqrt(s)), 2.0]

            return x
        end

        btest = brokentestmodel_assume2()
        @test_throws ArgumentError btest()

        # Test missing input arguments
        @model function testmodel_missing1(x)
            x ~ Bernoulli(0.5)
            return x
        end
        @test_throws MethodError testmodel_missing1()

        # Test missing initialization for vector observation turned parameter
        @model function testmodel_missing2(x)
            x[1] ~ Bernoulli(0.5)
            return x
        end
        @test_throws MethodError testmodel_missing2(missing)()

        # Test use of internal names
        @model function testmodel_missing3(x)
            x[1] ~ Bernoulli(0.5)
            global varinfo_ = __varinfo__
            global sampler_ = __context__.sampler
            global model_ = __model__
            global context_ = __context__
            global rng_ = __context__.rng
            global lp = getlogp(__varinfo__)
            return x
        end
        model = testmodel_missing3([1.0])
        varinfo = VarInfo(model)
        @test getlogp(varinfo) == lp
        @test varinfo_ isa AbstractVarInfo
        @test model_ === model
        @test context_ isa SamplingContext
        @test rng_ isa Random.AbstractRNG

        # disable warnings
        @model function testmodel_missing4(x)
            x[1] ~ Bernoulli(0.5)
            global varinfo_ = __varinfo__
            global sampler_ = __context__.sampler
            global model_ = __model__
            global context_ = __context__
            global rng_ = __context__.rng
            global lp = getlogp(__varinfo__)
            return x
        end false
        lpold = lp
        model = testmodel_missing4([1.0])
        varinfo = VarInfo(model)
        @test getlogp(varinfo) == lp == lpold

        # test DPPL#61
        @model function testmodel_missing5(z)
            m ~ Normal()
            z[1:end] ~ MvNormal(fill(m, length(z)), I)
            return m
        end
        model = testmodel_missing5(rand(10))
        @test all(z -> isapprox(z, 0; atol=0.2), mean(model() for _ in 1:1000))

        # test Turing#1464
        @model function gdemo(x)
            s ~ InverseGamma(2, 3)
            m ~ Normal(0, sqrt(s))
            for i in eachindex(x)
                x[i] ~ Normal(m, sqrt(s))
            end
        end
        x = [1.0, missing]
        VarInfo(gdemo(x))
        @test ismissing(x[2])

        # https://github.com/TuringLang/Turing.jl/issues/1464#issuecomment-731153615
        vi = VarInfo(gdemo(x))
        @test haskey(vi.metadata, :x)
        vi = VarInfo(gdemo(x))
        @test haskey(vi.metadata, :x)

        # Non-array variables
        @model function testmodel_nonarray(x, y)
            s ~ InverseGamma(2, 3)
            m ~ Normal(0, √s)
            for i in 1:(length(x.a) - 1)
                x.a[i] ~ Normal(m, √s)
            end

            # Dynamic indexing
            x.a[end] ~ Normal(100.0, 1.0)

            # Immutable set
            y.a ~ Normal()

            # Dotted
            z = Vector{Float64}(undef, 3)
            z[1:2] .~ Normal()
            z[end:end] .~ Normal()

            return (; s=s, m=m, x=x, y=y, z=z)
        end

        m_nonarray = testmodel_nonarray(
            MyCoolStruct([missing, missing]), MyCoolStruct(missing)
        )
        result = m_nonarray()
        @test !any(ismissing, result.x.a)
        @test result.y.a !== missing
        @test result.x.a[end] > 10

        # Ensure that we can work with `Vector{Real}(undef, N)` which is the
        # reason why we're using `BangBang.prefermutation` in `src/compiler.jl`
        # rather than the default from Setfield.jl.
        # Related: https://github.com/jw3126/Setfield.jl/issues/157
        @model function vdemo()
            x = Vector{Real}(undef, 10)
            for i in eachindex(x)
                x[i] ~ Normal(0, sqrt(4))
            end

            return x
        end
        x = vdemo()()
        @test all((isassigned(x, i) for i in eachindex(x)))
    end
    @testset "nested model" begin
        function makemodel(p)
            @model function testmodel(x)
                x[1] ~ Bernoulli(p)
                global lp = getlogp(__varinfo__)
                return x
            end
            return testmodel
        end
        model = makemodel(0.5)([1.0])
        varinfo = VarInfo(model)
        @test getlogp(varinfo) == lp
    end
    @testset "user-defined variable name" begin
        @model f1() = x ~ NamedDist(Normal(), :y)
        @model f2() = x ~ NamedDist(Normal(), @varname(y[2][:, 1]))
        @model f3() = x ~ NamedDist(Normal(), @varname(y[1]))
        vi1 = VarInfo(f1())
        vi2 = VarInfo(f2())
        vi3 = VarInfo(f3())
        @test haskey(vi1.metadata, :y)
        @test vi1.metadata.y.vns[1] == @varname(y)
        @test haskey(vi2.metadata, :y)
        @test vi2.metadata.y.vns[1] == @varname(y[2][:, 1])
        @test haskey(vi3.metadata, :y)
        @test vi3.metadata.y.vns[1] == @varname(y[1])
    end
    @testset "custom tilde" begin
        @model demo() = begin
            $(@custom m ~ Normal())
            return m
        end
        model = demo()
        @test all(iszero(model()) for _ in 1:1000)
    end
    @testset "docstring" begin
        "This is a test"
        @model function demo(x)
            m ~ Normal()
            return x ~ Normal(m, 1)
        end

        s = @doc(demo)
        @test string(s) == "This is a test\n"

        # Verify that adding docstring didn't completely break execution of model
        m = demo(0.0)
        @test m() isa Float64
    end
    @testset "type annotations" begin
        @model function demo_without(x)
            return x ~ Normal()
        end
        @test isempty(VarInfo(demo_without(0.0)))

        @model function demo_with(x::Real)
            return x ~ Normal()
        end
        @test isempty(VarInfo(demo_with(0.0)))
    end

    @testset "macros within model" begin
        # Macro expansion
        @model function demo1()
            @mymodel1(y ~ Uniform())
        end

        @test haskey(VarInfo(demo1()), @varname(x))

        # Interpolation
        # Will fail if:
        # 1. Compiler expands `y ~ Uniform()` before expanding the macros
        #    => returns -1.
        # 2. `@mymodel` is expanded before entire `@model` has been
        #    expanded => errors since `MyModelStruct` is not a distribution,
        #    and hence `tilde_observe` errors.
        @model function demo2()
            return $(@mymodel2(y ~ Uniform()))
        end
        @test demo2()() == 42
    end

    @testset "submodel" begin
        # No prefix, 1 level.
        @model function demo1(x)
            return x ~ Normal()
        end
        @model function demo2(x, y)
            @submodel demo1(x)
            return y ~ Uniform()
        end
        # No observation.
        m = demo2(missing, missing)
        vi = VarInfo(m)
        ks = keys(vi)
        @test @varname(x) ∈ ks
        @test @varname(y) ∈ ks

        # Observation in top-level.
        m = demo2(missing, 1.0)
        vi = VarInfo(m)
        ks = keys(vi)
        @test @varname(x) ∈ ks
        @test @varname(y) ∉ ks

        # Observation in nested model.
        m = demo2(1000.0, missing)
        vi = VarInfo(m)
        ks = keys(vi)
        @test @varname(x) ∉ ks
        @test @varname(y) ∈ ks

        # Observe all.
        m = demo2(1000.0, 0.5)
        vi = VarInfo(m)
        ks = keys(vi)
        @test isempty(ks)

        # Check values makes sense.
        @model function demo3(x, y)
            @submodel demo1(x)
            return y ~ Normal(x)
        end
        m = demo3(1000.0, missing)
        # Mean of `y` should be close to 1000.
        @test abs(mean([VarInfo(m)[@varname(y)] for i in 1:10]) - 1000) ≤ 10

        # Prefixed submodels and usage of submodel return values.
        @model function demo_return(x)
            x ~ Normal()
            return x
        end

        @model function demo_useval(x, y)
            @submodel prefix = "sub1" x1 = demo_return(x)
            @submodel prefix = "sub2" x2 = demo_return(y)

            return z ~ Normal(x1 + x2 + 100, 1.0)
        end
        m = demo_useval(missing, missing)
        vi = VarInfo(m)
        ks = keys(vi)
        @test VarName{Symbol("sub1.x")}() ∈ ks
        @test VarName{Symbol("sub2.x")}() ∈ ks
        @test @varname(z) ∈ ks
        @test abs(mean([VarInfo(m)[@varname(z)] for i in 1:10]) - 100) ≤ 10

        # AR1 model. Dynamic prefixing.
        @model function AR1(num_steps, α, μ, σ, ::Type{TV}=Vector{Float64}) where {TV}
            η ~ MvNormal(zeros(num_steps), I)
            δ = sqrt(1 - α^2)

            x = TV(undef, num_steps)
            x[1] = η[1]
            @inbounds for t in 2:num_steps
                x[t] = @. α * x[t - 1] + δ * η[t]
            end

            return @. μ + σ * x
        end

        @model function demo(y)
            α ~ Uniform()
            μ ~ Normal()
            σ ~ truncated(Normal(), 0, Inf)

            num_steps = length(y[1])
            num_obs = length(y)
            @inbounds for i in 1:num_obs
                @submodel prefix = "ar1_$i" x = AR1(num_steps, α, μ, σ)
                y[i] ~ MvNormal(x, 0.01 * I)
            end
        end

        ys = [randn(10), randn(10)]
        m = demo(ys)
        vi = VarInfo(m)

        for k in [:α, :μ, :σ, Symbol("ar1_1.η"), Symbol("ar1_2.η")]
            @test VarName{k}() ∈ keys(vi)
        end
    end

    @testset "check_tilde_rhs" begin
        @test_throws ArgumentError DynamicPPL.check_tilde_rhs(randn())

        x = Normal()
        @test DynamicPPL.check_tilde_rhs(x) === x

        x = [Laplace(), Normal(), MvNormal(zeros(3), I)]
        @test DynamicPPL.check_tilde_rhs(x) === x
    end
    @testset "isliteral" begin
        @test DynamicPPL.isliteral(:([1.0]))
        @test DynamicPPL.isliteral(:([[1.0], 1.0]))
        @test DynamicPPL.isliteral(:((1.0, 1.0)))

        @test !(DynamicPPL.isliteral(:([x])))
        @test !(DynamicPPL.isliteral(:([[x], 1.0])))
        @test !(DynamicPPL.isliteral(:((x, 1.0))))
    end

    @testset "array literals" begin
        # Verify that we indeed can parse this.
        @test @model(function array_literal_model()
            # `assume` and literal `observe`
            m ~ MvNormal(zeros(2), I)
            return [10.0, 10.0] ~ MvNormal(m, 0.25 * I)
        end) isa Function

        @model function array_literal_model2()
            # `assume` and literal `observe`
            m ~ MvNormal(zeros(2), I)
            return [10.0, 10.0] ~ MvNormal(m, 0.25 * I)
        end

        @test array_literal_model2()() == [10.0, 10.0]
    end

    # https://github.com/TuringLang/DynamicPPL.jl/issues/260
    @testset "anonymous function" begin
        error = ArgumentError("anonymous functions without name are not supported")
        @test_throws LoadError(@__FILE__, (@__LINE__) + 1, error) @macroexpand begin
            @model function (x)
                return x ~ Normal()
            end
        end
        @test_throws LoadError(@__FILE__, (@__LINE__) + 1, error) @macroexpand begin
            model = @model(x -> (x ~ Normal()))
        end
    end

    @testset "dispatching with model" begin
        f(x) = false

        @model demo() = x ~ Normal()
        @test !f(demo())
        f(::Model{typeof(demo)}) = true
        @test f(demo())

        # Leads to re-definition of `demo` and trait is not affected.
        @test length(methods(demo)) == 2
        @model demo() = x ~ Normal()
        @test length(methods(demo)) == 2
        @test f(demo())

        # Ensure we can specialize on arguments.
        @model demo(x) = x ~ Normal()
        @test length(methods(demo)) == 4
        @test f(demo(1.0))
        f(::Model{typeof(demo),(:x,)}) = false
        @test !f(demo(1.0))
        @test f(demo()) # should still be `true`

        # Set it to `false` again.
        f(::Model{typeof(demo),()}) = false
        @test !f(demo())
    end

    @testset "return value" begin
        # Make sure that a return-value of `x = 1` isn't combined into
        # an attempt at a `NamedTuple` of the form `(x = 1, __varinfo__)`.
        @model empty_model() = return x = 1
        empty_vi = VarInfo()
        retval_and_vi = DynamicPPL.evaluate!!(empty_model(), empty_vi, SamplingContext())
        @test retval_and_vi isa Tuple{Int,typeof(empty_vi)}

        # Even if the return-value is `AbstractVarInfo`, we should return
        # a `Tuple` with `AbstractVarInfo` in the second component too.
        @model demo() = return __varinfo__
        retval, svi = DynamicPPL.evaluate!!(demo(), SimpleVarInfo(), SamplingContext())
        @test svi == SimpleVarInfo()
        if Threads.nthreads() > 1
            @test retval isa DynamicPPL.ThreadSafeVarInfo{<:SimpleVarInfo}
            @test retval.varinfo == svi
        else
            @test retval == svi
        end

        # We should not be altering return-values other than at top-level.
        @model function demo()
            # If we also replaced this `return` inside of `f`, then the
            # final `return` would be include `__varinfo__`.
            f(x) = return x^2
            return f(1.0)
        end
        retval, svi = DynamicPPL.evaluate!!(demo(), SimpleVarInfo(), SamplingContext())
        @test retval isa Float64

        @model demo() = x ~ Normal()
        retval, svi = DynamicPPL.evaluate!!(demo(), SimpleVarInfo(), SamplingContext())

        # Return-value when using `@submodel`
        @model inner() = x ~ Normal()
        # Without assignment.
        @model outer() = @submodel inner()
        @test outer()() isa Real

        # With assignment.
        @model outer() = @submodel x = inner()
        @test outer()() isa Real
    end

    @testset "issue #368: hasmissing dispatch" begin
        @test !DynamicPPL.hasmissing(typeof(Union{}[]))

        # (nested) arrays with `Missing` eltypes
        @test DynamicPPL.hasmissing(Vector{Union{Missing,Float64}})
        @test DynamicPPL.hasmissing(Matrix{Union{Missing,Real}})
        @test DynamicPPL.hasmissing(Vector{Matrix{Union{Missing,Float32}}})

        # no `Missing`
        @test !DynamicPPL.hasmissing(Vector{Float64})
        @test !DynamicPPL.hasmissing(Matrix{Real})
        @test !DynamicPPL.hasmissing(Vector{Matrix{Float32}})
    end

    @testset "issue #393: anonymous argument with type parameter" begin
        @model f_393(::Val{ispredict}=Val(false)) where {ispredict} = ispredict ? 0 : 1
        @test f_393()() == 1
        @test f_393(Val(true))() == 0
    end
end

max_stars_repo_name: TuringLang/MicroPPL

max_issues_repo_name: TuringLang/MicroPPL

max_forks_repo_name: TuringLang/MicroPPL


text:
#####
##### Tests copied from IterativeSolvers.jl
##### https://github.com/JuliaLinearAlgebra/IterativeSolvers.jl/blob/v0.9.2/src/lsmr.jl
#####

# Type used in Dampenedtest
# solve (A'A + diag(v).^2 ) x = A'b
# using LSMR in the augmented space Ã = [A ; diag(v)] b̃ = [b; zeros(size(A, 2)]
struct DampenedMatrix{Tv,TA<:AbstractMatrix{Tv},TD<:AbstractVector{Tv}} <: AbstractMatrix{Tv}
    A::TA
    diagonal::TD
end

function Base.size(A::DampenedMatrix)
    m, n = size(A.A)
    l = length(A.diagonal)
    (m + l, n)
end

function Base.size(A::DampenedMatrix, dim::Integer)
    m, n = size(A.A)
    l = length(A.diagonal)
    dim == 1 ? (m + l) : (dim == 2 ? n : 1)
end

function LinearAlgebra.mul!(y::AbstractVector{Tv}, mw::DampenedMatrix, x::AbstractVector{Tv}) where {Tv}
    m₁ = size(mw.A, 1)
    m₂ = size(mw, 1)
    mul!(view(y, 1:m₁), mw.A, x)
    y[m₁+1:m₂] .= mw.diagonal .* x
    return y
end

function LinearAlgebra.mul!(y::AbstractVector, mw::Adjoint{Tv,<:DampenedMatrix}, x::AbstractVector) where {Tv}
    m₁ = size(mw.parent.A, 1)
    m₂ = size(mw.parent, 1)
    mul!(y, adjoint(mw.parent.A), view(x, 1:m₁))
    y .+= mw.parent.diagonal .* view(x, m₁+1:m₂)
    return y
end

"""
Produces the m × n submatrix from
A = [ 1
      1 2
        2 3
          3 4
            ...
              n ]
suitably padded by zeros.
"""
function sol_matrix(m, n)
    mn = min(m, n)
    I, J, V = SparseArrays.spdiagm_internal(-1 => 1.0 : mn - 1, 0 => 1.0 : mn)
    sparse(I, J, V, m, n)
end

@testset "Small dense matrix" for T = (Float32, Float64, ComplexF32, ComplexF64)
    Random.seed!(501)
    A = rand(T, 10, 5)
    b = rand(T, 10)
    x = QSM.lsmr(A, b)
    @test norm(x - A\b) ≤ √eps(real(T))
end

@testset "SOL test" for (m, n, damp) = ((10, 10, 0), (20, 10, 0), (20, 10, 0.1))
    # Test adapted from the BSD-licensed Matlab implementation at
    #    http://www.stanford.edu/group/SOL/software/lsqr.html
    #              Michael Saunders, Systems Optimization Laboratory,
    #              Dept of MS&E, Stanford University.
    #-----------------------------------------------------------------------
    # 11 Apr 1996: First version for distribution with lsqr.m.
    #              Michael Saunders, Dept of EESOR, Stanford University.

    A = sol_matrix(m, n)
    x = float(n : -1 : 1)
    b = A * x
    x_lsmr = QSM.lsmr(A, b, atol = 1e-7, btol = 1e-7, conlim = 1e10, maxit = 10n)
    @test norm(b - A * x) ≤ 1e-4
end

@testset "Dampened test" for (m, n) = ((10, 10), (20, 10))
    Random.seed!(501)
    # Test used to make sure A, b can be generic matrix / vector
    b = rand(m)
    A = rand(m, n)
    v = rand(n)
    A′ = DampenedMatrix(A, v)
    b′ = [b; zeros(n)]
    x = QSM.lsmr(A′, b′)
    @test norm((A'A + Matrix(Diagonal(v)) .^ 2)x - A'b) ≤ 1e-3
end


text:
using transD_GP

# Forward setttings
# dummy physics struct defined in src/FooPhysics.jl
F = transD_GP.FooPhysics.Foo(25, rand(10))

# some dummy forward
@info transD_GP.FooPhysics.returnphysics!(F, rand(3))

# struct for data and physics operator defined in src/FooPhysicsInversion.jl
FOp = transD_GP.FooPhysicsInversion.FooInversion(rand(10), F)

## Inverse settings
# prior bounds and setings
fbounds = [0 1.]
xall = collect(permutedims(0:0.1:1))
xbounds = permutedims([extrema(xall)...])
λ = vec(diff(xbounds, dims=2))
nmin, nmax = 2, 5
demean, sampledc = false, true
# McMC proposals
δ = 0.1*diff(fbounds, dims=2)[1]
sdev_pos = abs.(vec(diff(xbounds, dims=2)))
sdev_prop = abs.(vec(0.05*diff(fbounds, dims=2)))
K = transD_GP.GP.OrstUhn()

# make MCMC options using prior and proposals
opt = transD_GP.OptionsStat(nmin = nmin,
                        nmax = nmax,
                        xbounds = xbounds,
                        fbounds = fbounds,
                        xall = xall,
                        λ = λ,
                        δ = δ,
                        demean = demean,
                        sampledc = sampledc,
                        sdev_prop = sdev_prop,
                        sdev_pos = sdev_pos,
                        quasimultid = false,
                        K = K
                        )

## Put forward and inverse together
# initialise a random model
m = transD_GP.init(opt)

# compute misfit
transD_GP.get_misfit(m, opt, FOp)


text:
# ---
# title: 165. Compare Version Numbers
# id: problem165
# author: Indigo
# date: 2021-01-19
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/compare-version-numbers/description/>
# hidden: true
# ---
#
# Given two version numbers, `version1` and `version2`, compare them.
#
# Version numbers consist of **one or more revisions** joined by a dot `'.'`.
# Each revision consists of **digits**  and may contain leading **zeros**. Every
# revision contains **at least one character**. Revisions are **0-indexed from
# left to right** , with the leftmost revision being revision 0, the next
# revision being revision 1, and so on. For example `2.5.33` and `0.1` are valid
# version numbers.
#
# To compare version numbers, compare their revisions in **left-to-right
# order**. Revisions are compared using their  **integer value ignoring any
# leading zeros**. This means that revisions `1` and `001` are considered
# **equal**. If a version number does not specify a revision at an index, then
# **treat the revision as  `0`**. For example, version `1.0` is less than
# version `1.1` because their revision 0s are the same, but their revision 1s
# are `0` and `1` respectively, and `0 < 1`.
#
# _Return the following:_
#
#   * If `version1 < version2`, return `-1`.
#   * If `version1 > version2`, return `1`.
#   * Otherwise, return `0`.
#
#
#
# **Example 1:**
#
#
#
#     Input: version1 = "1.01", version2 = "1.001"
#     Output: 0
#     Explanation: Ignoring leading zeroes, both "01" and "001" represent the same integer "1".
#
#
# **Example 2:**
#
#
#
#     Input: version1 = "1.0", version2 = "1.0.0"
#     Output: 0
#     Explanation: version1 does not specify revision 2, which means it is treated as "0".
#
#
# **Example 3:**
#
#
#
#     Input: version1 = "0.1", version2 = "1.1"
#     Output: -1
#     Explanation:  version1's revision 0 is "0", while version2's revision 0 is "1". 0 < 1, so version1 < version2.
#
#
# **Example 4:**
#
#
#
#     Input: version1 = "1.0.1", version2 = "1"
#     Output: 1
#
#
# **Example 5:**
#
#
#
#     Input: version1 = "7.5.2.4", version2 = "7.5.3"
#     Output: -1
#
#
#
#
# **Constraints:**
#
#   * `1 <= version1.length, version2.length <= 500`
#   * `version1` and `version2` only contain digits and `'.'`.
#   * `version1` and `version2`  **are valid version numbers**.
#   * All the given revisions in `version1` and `version2` can be stored in a  **32-bit integer**.
#
#
## @lc code=start
using LeetCode

function compare_version(version1::String, version2::String)
    s_v1 = parse.(Int, split(version1, "."))
    s_v2 = parse.(Int, split(version2, "."))
    if length(s_v1) > length(s_v2)
        append!(s_v2, fill(0, length(s_v1) - length(s_v2)))
    else
        append!(s_v1, fill(0, length(s_v2) - length(s_v1)))
    end
    for i in 1:length(s_v1)
        if s_v1 > s_v2
            return 1
        elseif s_v1 < s_v2
            return -1
        end
    end
    return 0
end
## @lc code=end


text:
using FixedPointAcceleration
using Test

# Run tests
println("Test putting together iterates with jumps")
@time @test include("TestPuttingInputsAndOutputsTogether.jl")
println("Test simple Scalar function")
@time @test include("SimpleScalarFunction.jl")
println("Test simple vector function")
@time @test include("SimpleVectorFunction.jl")
println("Test Chaining together calls with cos functions.")
@time @test include("Chaining.jl")
println("Testing Bounds Error")
@time @test include("BoundsError.jl")
println("Testing ReplaceInvalidsTest")
@time @test include("ReplaceInvalidsTesting.jl")
println("CommonIncrementsTest")
@time @test include("CommonIncrementsTest.jl")
# Note consumption smoothing problem is not tested for time and dependency reasons.
# Note that the Autodiff extention is not included for the same reason.
# And test_side_products for the same reason.


text:
using Parameters

@with_kw mutable struct Constants
    CellTyp::String = "Doyle_94"
    T::Float64 = 298.15 # Cell Temperature
    T_ref::Float64 = 298.15 # Reference Temperature
    t_plus::Float64 = 0.2594 # Inital Transference Number
    tpf::Function = ce -> -0.1287*(ce/1000)^3+0.4106*(ce/1000)^2-0.4717*(ce/1000)+0.4492 # Transference Number Function - Requires ce in dm-3
    De::Float64 = 1.769e-10   # Inital Electrolyte Diffusivity
    Def::Function = ce -> 8.794e-11*(ce/1000)^2-3.972e-10*(ce/1000)+4.862e-10 # Electrolyte Diffusivity Function - Requires ce in dm-3
    SOC::Float64 = 0.747 # Initial State of Charge
    ce0::Float64 = 1000 # Initial Electrolyte Concentration
    #dln::Float64 = 3.
    Ea_κ = 0.
    Ea_De::Float64 = 0
    CC_A::Float64 =  0.1027  #Electrode Plate Area
    κ::Float64 = 0.9487
    κf::Function = ce -> 0.1297*(ce/1000)^3-2.51*(ce/1000)^1.5+3.329*(ce/1000) #Requires ce in dm-3
    Uocp::Function = (Electrode, θ) ->
        if Electrode == "Neg"
            Uocp = @. 1.97938*2.7182818284*exp(-39.3631*θ) + 0.2482 - 0.0909*tanh(29.8538*(θ - 0.1234)) - 0.04478*tanh(14.9159*(θ - 0.2769)) - 0.0205*tanh(30.4444*(θ - 0.6103))
        else
            Uocp = @. -0.8090*θ + 4.4875 - 0.0428*tanh(18.5138*(θ - 0.5542)) - 17.7326*tanh(15.7890*(θ - 0.3117)) + 17.5842*tanh(15.9308*(θ - 0.3120))
        end
    ∂Uocp::Function = (Electrode,θ) ->
        if Electrode == "Neg"
            #∂Uocp = @. 0.667934002(tanh(14.9159*θ - 4.13021271)^2) + 2.71371042(tanh(29.8538*θ - 3.68395892)^2) + 0.6241102(tanh(30.4444*θ - 18.58021732)^2) - 4.005754622 - (77.91453287630758(2.7182818284^(-39.3631*θ)))
            ∂Uocp = @. -0.62411*((sech(18.5802 - 30.4444*θ))^2 + 4.34813*(sech(3.68396 - 29.8538*θ))^2 + 1.07022*(sech(4.13021 - 14.9159*θ))^2) - 211.794*exp(-39.3631*θ)
        else
            ∂Uocp = @. 279.9800214*(tanh(15.789*θ - 4.9214313)^2) + 0.79239064*(tanh(18.5138*θ - 10.26034796)^2) - 1.4510386800000594 - (280.13037336*(tanh(15.9308*θ - 4.9704096)^2))
        end
    Ce_M::Int64 = 1     #(Rewritten)
    D1::Float64 = 1.0
    D2::Float64 = 1.0
    D3::Float64 = 1.0
    Ltot::Float64 = 0.
    Lnegsep::Float64 = 0.
end

@with_kw mutable struct Negative
    L::Float64 = 85.2e-6    # Negative Electrode Length
    Rs::Float64 = 5.86e-6   # Particle Radius [m]
    Ds::Float64 = 3.3e-14   # Solid Diffusivity [m^2/s]
    Ea_σ::Float64 = 0     # Activation Energy Solid Conductivity
    Ea_Ds::Float64 = 0    # Activation Energy Solid Diffusivity
    ϵ_s::Float64 = 0.75     # Active Material Volume Fraction
    ϵ_e::Float64 = 0.25     # Porosity of Negative Electrode
    De_brug::Float64 = 1.5  # Bruggeman Diffusion Exponent
    κ_brug::Float64 = 1.5   # Bruggeman Electrolyte Conductivity Exponent
    σ::Float64 = 215    # Solid Phase Conductivity
    σ_brug::Float64 = 1.5   # Bruggeman Solid Conductivity Exponent
    θ_100::Float64 = 0.910612#0.9014 # Theta @ 100% Lithium Concentration
    θ_0::Float64 = 0.0263473 #0.0279   # Theta @ 0% Lithium Concentration
    cs_max::Float64 = 33133 # Max Electrode Concentration
    α::Float64 = 0.5    # Alpha Factor
    k_norm::Float64 = 7.264265E-06 #6.48e-6 #7.226781E-06 #6.8973799e-13 #4.1580e-8 #2.12e-10 #Initial Reaction Rate
    Ea_κ::Float64 = 35000   # Activation Energy
    RFilm::Float64 = 0 # Film Resistance
    D1::Float64 = 1.   # Init Value
    D1f::Function = De -> De * ϵ_e^De_brug # Effective Diffusivity
    as::Float64 = 3.0*ϵ_s/Rs # Specific Interfacial Surf. Area
end

@with_kw mutable struct Positive
    L::Float64 = 75.6e-6    # Positive Electrode Length
    Rs::Float64 = 5.22e-6    # Particle radius [m]
    Ds::Float64 = 4e-15   # Solid diffusivity [m^2/s]
    Ea_σ::Float64 = 0
    Ea_Ds::Float64 = 0
    ϵ_s::Float64 = 0.665    # Active Material Volume Fraction
    ϵ_e::Float64 = 0.335    # Porosity of positive electrode
    De_brug::Float64 = 1.5  # Bruggeman Diffusivity Exponent
    κ_brug::Float64 = 1.5   # Bruggeman Electrolyte Conductivity Exponent
    σ::Float64 = 0.18   # Solid Phase Conductivity
    σ_brug::Float64 = 1.5   # Bruggeman Solid Conductivity Exponent
    θ_100::Float64 = 0.263849 #0.27 # Theta @ 100% Lithium Concentration
    θ_0::Float64 = 0.853974 #0.9084   # Theta @ 0% Lithium Concentration
    cs_max::Float64 = 63104 # Max Electrode Concentration
    α::Float64 = 0.5    # Alpha Factor
    k_norm::Float64 =  7.26426E-05 #9E-05 #3.640283886203905e-12 #3.5954e-7 #3.42e-6 #1.12e-9  #Initial Reaction Rate
    Ea_κ::Float64 = 17800   # Activation Energy
    RFilm::Float64 = 0 # Film Resistance
    D3::Float64 = 1   # Init Value
    D3f::Function = De -> De * ϵ_e^De_brug
    as::Float64 = 3.0*ϵ_s/Rs # Specific interfacial surf. area
end

@with_kw mutable struct Seperator
    L::Float64 = 12e-6  # Seperator Length
    ϵ_e::Float64 = 0.47    # Porosity of separator
    De_brug::Float64 = 1.5  # Bruggeman Diffusivity Factor
    κ_brug::Float64 = 1.5   # Bruggeman Electrolyte Conductivity Factor
    D2::Float64 = 1   # Init Value
    D2f::Function = De -> De * ϵ_e^De_brug
end

@with_kw mutable struct RealisationAlgorthim
    Fs::Float64 = 8    # Sampling Frequency of Transfer Functions [Hz]
    SamplingT::Float64 = 0.25     # Final Model Sampling Time [s]
    M::Int64 = 6    # Model Order
    N::Int64 = 1    # Number of Inputs
    Tlen::Int64 = 18000 #64800 #21600 #16200 #28800 #43200 #65536 #131072 #1048576 #2097152 #262144 #32768 #24    #Transfer Function Response Length [s] (Change to min)
    H1::Array{Int64,1} = 0:2000 #4000 #4612     # Hankel Dimensions 1
    H2::Array{Int64,1} = 0:2000 #4000 #4612     # Hankel Dimensions 2
    Outs::Int64 = 1    # Number of Outputs (Rewritten)
    Nfft::Float64 = 0
    f::Array{Float64,1} = [0]
    s::Array{ComplexF64} = [0 - 0.0im]
    Nfft!::Function = (Fs,Tlen) -> ceil(2^(log2(Fs*Tlen))) #2^(ceil(log2(Fs*Tlen))) Old way (constains to 2^ values)
    f!::Function = (Nfft) -> 0:Nfft-1
    s!::Function = (Fs,Nfft,f) -> transpose(((2im.*Fs)*tan.(pi.*f./Nfft)))
end

@with_kw mutable struct TransferFun
    #6 - 4
    #tfs =   [[C_e, Phi_e, C_se, Phi_s, Phi_se, Flux, C_se, Phi_s, Flux, Phi_se] ["Na", "Na", "Pos", "Pos", "Pos", "Pos", "Neg", "Neg", "Neg", "Neg"] [Float64[0.00,4.26e-5,8.52e-5,9.72e-5,1.35e-04,1.728e-4], Float64[4.26e-5,8.52e-5,9.72e-5,1.35e-04,1.728e-4], Float64[0,0.333,0.666,1],Float64[1],Float64[0,0.333,0.666,1],Float64[0,0.333,0.666,1],Float64[0,0.333,0.666,1],Float64[1],Float64[0,0.333,0.666,1],Float64[0,0.333,0.666,1]]]

    #4 - 4
    #tfs =   [[C_e, Phi_e, C_se, Phi_s, Phi_se, Flux, C_se, Phi_s, Flux, Phi_se] ["Na", "Na", "Pos", "Pos", "Pos", "Pos", "Neg", "Neg", "Neg", "Neg"] [Float64[0.0, 8.52e-5,9.72e-5,0.0001728], Float64[8.52e-5,9.72e-5,0.0001728], Float64[0,0.333,0.666,1],Float64[1],Float64[0,0.333,0.666,1],Float64[0,0.333,0.666,1],Float64[0,0.333,0.666,1],Float64[1],Float64[0,0.333,0.666,1],Float64[0,0.333,0.666,1]]]

    #8 - 4
    #tfs =   [[C_e, Phi_e, C_se, Phi_s, Phi_se, Flux, C_se, Phi_s, Flux, Phi_se] ["Na", "Na", "Pos", "Pos", "Pos", "Pos", "Neg", "Neg", "Neg", "Neg"] [Float64[0,2.84e-5,5.68e-5,8.52e-5,9.72e-5,0.0001224,1.224e-4,1.476e-4,1.728e-4], Float64[2.84e-5,5.68e-5,8.52e-5,9.72e-5,0.0001224,1.224e-4,1.476e-4,1.728e-4], Float64[0,0.333,0.666,1],Float64[1],Float64[0,0.333,0.666,1],Float64[0,0.333,0.666,1],Float64[0,0.333,0.666,1],Float64[1],Float64[0,0.333,0.666,1],Float64[0,0.333,0.666,1]]]

    # 9 - 6
    #tfs =   [[C_e, Phi_e, C_se, Phi_s, Phi_se, Flux, C_se, Phi_s, Flux, Phi_se] ["Na", "Na", "Pos", "Pos", "Pos", "Pos", "Neg", "Neg", "Neg", "Neg"] [Float64[0,2.84e-5,5.68e-5,8.52e-5,9.12e-5,9.72e-5,0.0001224,0.0001476,0.0001728], Float64[2.84e-5,5.68e-5,8.52e-5,9.12e-5,9.72e-5,0.0001224,0.0001476,0.0001728], Float64[0,0.2,0.4,0.6,0.8,1], Float64[1],Float64[0,0.2,0.4,0.6,0.8,1],Float64[0,0.2,0.4,0.6,0.8,1],Float64[0,0.2,0.4,0.6,0.8,1],Float64[1],Float64[0,0.2,0.4,0.6,0.8,1],Float64[0,0.2,0.4,0.6,0.8,1]]]

    # 4 - 2
    tfs =   [[C_e, Phi_e, C_se, Phi_s, Phi_se, Flux, C_se, Phi_s, Flux, Phi_se] ["Na", "Na", "Pos", "Pos", "Pos", "Pos", "Neg", "Neg", "Neg", "Neg"] [Float64[0, 8.52E-05, 9.72E-05, 1.728E-04], Float64[8.52E-05, 9.72E-05, 1.728E-04], Float64[0,1], Float64[1], Float64[0,1], Float64[0,1], Float64[0,1], Float64[1], Float64[0,1], Float64[0,1]]]

    # 4 - 6
    ##tfs =   [[C_e, Phi_e, C_se, Phi_s, Phi_se, Flux, C_se, Phi_s, Flux, Phi_se] ["Na", "Na", "Pos", "Pos", "Pos", "Pos", "Neg", "Neg", "Neg", "Neg"] [Float64[0.00,4.26e-5,8.52e-5,9.72e-5,1.35e-04,1.728e-4], Float64[4.26e-5,8.52e-5,9.72e-5,1.35e-04,1.728e-4], Float64[0,0.2,0.4,0.6,0.8,1], Float64[1],Float64[0,0.2,0.4,0.6,0.8,1],Float64[0,0.2,0.4,0.6,0.8,1],Float64[0,0.2,0.4,0.6,0.8,1],Float64[1],Float64[0,0.2,0.4,0.6,0.8,1],Float64[0,0.2,0.4,0.6,0.8,1]]]

    # 6 - 2
    #tfs =   [[C_e, Phi_e, C_se, Phi_s, Phi_se, Flux, C_se, Phi_s, Flux, Phi_se] ["Na", "Na", "Pos", "Pos", "Pos", "Pos", "Neg", "Neg", "Neg", "Neg"] [Float64[0.00,4.26e-5,8.52e-5,9.72e-5,1.35e-04,1.728e-4], Float64[4.26e-5,8.52e-5,9.72e-5,1.35e-04,1.728e-4], Float64[0,1], Float64[1],Float64[0,1],Float64[0,1],Float64[0,1],Float64[1],Float64[0,1],Float64[0,1]]]

    # 6 - 6
    #tfs =   [[C_e, Phi_e, C_se, Phi_s, Phi_se, Flux, C_se, Phi_s, Flux, Phi_se] ["Na", "Na", "Pos", "Pos", "Pos", "Pos", "Neg", "Neg", "Neg", "Neg"] [Float64[0.00,4.26e-5,8.52e-5,9.72e-5,1.35e-04,1.728e-4], Float64[4.26e-5,8.52e-5,9.72e-5,1.35e-04,1.728e-4], Float64[0,0.2,0.4,0.6,0.8,1], Float64[1],Float64[0,0.2,0.4,0.6,0.8,1],Float64[0,0.2,0.4,0.6,0.8,1],Float64[0,0.2,0.4,0.6,0.8,1],Float64[1],Float64[0,0.2,0.4,0.6,0.8,1],Float64[0,0.2,0.4,0.6,0.8,1]]]

    #Entries::Function = (x) -> for i in 1:size(tfs[:,3],1) return [x;tfs[i,3]] end
end

@with_kw mutable struct Params
    Const::Constants
    Neg::Negative
    Pos::Positive
    Sep::Seperator
    RA::RealisationAlgorthim
    Transfer::TransferFun
end

Cell = Params(Constants(),Negative(),Positive(),Seperator(),RealisationAlgorthim(),TransferFun())
Cell.Const.Lnegsep, Cell.Const.Ltot = Cell.Neg.L+Cell.Sep.L,Cell.Neg.L+Cell.Sep.L+Cell.Pos.L
Cell.Const.D1 = Cell.Const.De*Cell.Neg.ϵ_e^Cell.Neg.De_brug
Cell.Const.D2 = Cell.Const.De*Cell.Sep.ϵ_e^Cell.Sep.De_brug
Cell.Const.D3 = Cell.Const.De*Cell.Pos.ϵ_e^Cell.Pos.De_brug
Cell.Const.Ce_M = size(Cell.Transfer.tfs[1,3],1)
Cell.RA.Outs = sum([size(Cell.Transfer.tfs[i,3],1) for i in 1:size(Cell.Transfer.tfs[:,1],1)])
#Cell.RA.Outs = size(Cell.Transfer.Entries,1)
# Cell.RA.Nfft = 2^(ceil(log2(Cell.RA.Fs*Cell.RA.Tlen)))
# Cell.RA.f = 0:Nfft-1
# Cell.RA.s = transpose(((2im.*Cell.RA.Fs)*tan.(pi.*f./Nfft)))

max_stars_repo_name: BradyPlanden/LIBRA

max_issues_repo_name: BradyPlanden/LIBRA

max_forks_repo_name: BradyPlanden/LIBRA


text:
#=
author: mingchuno
simple permutation brute fource
=#
array = [0 ,2 ,3, 4, 5, 6, 7 ,8, 9]

for arr in permutations(array)
  a = arr[1]
  if a == 0 continue end # fail fast

  b = arr[2]

  c = arr[3]
  if c == 0 continue end # fail fast

  d = arr[4]

  e = arr[5]
  if e == 0 continue end # fail fast

  f = arr[6]

  g = arr[7]
  if g == 0 continue end # fail fast

  h = arr[8]

  ab = a*10 + b
  cd = c*10 + d
  ef = e*10 + f
  gh = g*10 + h

  if ab - cd == ef && ef + gh == 111
    println("($(ab) - $(cd) = $(ef)) + $(gh) = 111")
  end
end

max_stars_repo_name: jackhftang/ABCDEFGHPPP

max_issues_repo_name: jackhftang/ABCDEFGHPPP

max_forks_repo_name: jackhftang/ABCDEFGHPPP


text:
module lotka_volterra
using ReachabilityAnalysis
@taylorize function lotka_volterra!(du, u, p, t)
    local α, β, γ, δ = 1.5, 1.0, 3.0, 1.0
    du[1] = u[1] * (α - β*u[2])
    du[2] = -u[2] * (γ - δ*u[1])
    return du
end

function model(X0)
    S = @system(x' = lotka_volterra!(x), dim:2)
    return IVP(S, X0)
end
end # module

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
export has_symmetry

"""
    hassym = has_symmetry(relax_ctx, pb::Problem, symtype::T) where T

    Check the `pb` problem for general symmetry `sym`, by examining the objective polynomial and then all constraints bodys.
"""
function has_symmetry(relax_ctx, pb::Problem, sym::T) where T
    # Check objective function for symmetry sym
    if !has_symmetry(relax_ctx, pb.objective, sym)
        return false
    end

    # Check constraints for symmetry sym
    for (cstrname, cstr) in pb.constraints
        if !has_symmetry(relax_ctx, cstr.p, sym)
            return false
        end
    end
    return true
end

"""
    hassym = has_symmetry(relax_ctx, p::Polynomial, sym::T) where T

    Check the `p` polynomial for general symmetry `sym`, by checking each exponent for `sym` symmetry.
"""
function has_symmetry(relax_ctx, p::Polynomial, sym::T) where T
    for (expo, λ) in p
        if !has_symmetry(relax_ctx, expo, sym)
            return false
        end
    end
    return true
end

## Phase Invariance symmetry
"""
    has_symmetry(relax_ctx, expo::Exponent, sym::T) where T<:Type{PhaseInvariance}

    Check whether `expo` is phase invariant, i.e. the exponent value is the same when evaluated at a complex point or a phase shifted complex point (*-1 for real variables).
"""
function has_symmetry(relax_ctx, expo::Exponent, sym::T) where T<:Type{PhaseInvariance}
    explsum, conjsum = get_sumdegs(expo)
    if relax_ctx.hierarchykind == :Real
        (conjsum != 0) && error(LOGGER, "is_homogeneous(): Exponent $expo has conjugated variables for a real hierarchy.")
        return explsum % 2 == 0
    elseif relax_ctx.hierarchykind == :Complex
        return explsum == conjsum
    else
        error(LOGGER, "is_homogeneous(expo, kind): kind should be either :Real or :Complex ($kind here).")
    end
end


text:
using AtomicMiscellany
using Test

@testset "AtomicMiscellany.jl" begin
    include("nuclearmodels.jl")
    include("hydrogenic_energies.jl")
end


text:
# Autogenerated wrapper script for alsa_jll for x86_64-linux-gnu
export libatopology, libasound

## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"

# Relative path to `libatopology`
const libatopology_splitpath = ["lib", "libatopology.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libatopology_path = ""

# libatopology-specific global declaration
# This will be filled out by __init__()
libatopology_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libatopology = "libatopology.so.2"

# Relative path to `libasound`
const libasound_splitpath = ["lib", "libasound.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libasound_path = ""

# libasound-specific global declaration
# This will be filled out by __init__()
libasound_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libasound = "libasound.so.2"

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"alsa")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # We first need to add to LIBPATH_list the libraries provided by Julia
    LIBPATH_list = [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)]
    global libatopology_path = normpath(joinpath(artifact_dir, libatopology_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libatopology_handle = dlopen(libatopology_path)
    push!(LIBPATH_list, dirname(libatopology_path))

    global libasound_path = normpath(joinpath(artifact_dir, libasound_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libasound_handle = dlopen(libasound_path)
    push!(LIBPATH_list, dirname(libasound_path))

    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(LIBPATH_list, ':')

    # Add each element of LIBPATH to our DL_LOAD_PATH (necessary on platforms
    # that don't honor our "already opened" trick)
    #for lp in LIBPATH_list
    #    push!(DL_LOAD_PATH, lp)
    #end
end  # __init__()


text:
# Also see related functions in users.jl and follow.jl

"""
    top_tracks(;offset=0, limit=20, time_range="medium_term")

**Summary**: Get the current user's top tracks based on calculated affinity.

# Optional keywords
- `offset` _Optional_: The index of the first tracks to return. Default 0.\n
- `limit` _Optional_: The number of tracks to return. Default 20. \n
- `time_range` _Optional_: Over what time frame the affinities are computed.
                            Valid Options:
                            * `long_term` : Calculated from several years of data including all new data as it becomes available
                            * `medium_term` : Approx. last 6 months
                            * `short_term` : Approx. last 4 weeks
"""
function top_tracks(;offset=0, limit=20, time_range="medium_term")

    url1 = "me/top/tracks?offset=$offset&limit=$limit"
    url2 = "&time_range=$time_range"

    return spotify_request(url1 * url2; scope = "user-top-read")

end

"""
    top_artists(;offset=0, limit=20, time_range="medium_term")

**Summary**: Get the current user's top artists based on calculated affinity.

# Optional keywords
- `offset` _Optional_: The index of the first tracks to return. Default 0.\n
- `limit` _Optional_: The number of tracks to return. Default 20. \n
- `time_range` _Optional_: Over what time frame the affinities are computed.
                            Valid Options:
                            * `long_term` : Calculated from several years of data including all new data as it becomes available
                            * `medium_term` : Approx. last 6 months
                            * `short_term` : Approx. last 4 weeks
"""
function top_artists(;offset=0, limit=20, time_range="medium_term")

    url1 = "me/top/artists?offset=$offset&limit=$limit"
    url2 = "&time_range=$time_range"

    return spotify_request(url1 * url2; scope = "user-top-read")

end


text:
# # ZEUS interface
#
# Here is a short demonstration of how to access the ZEUS transfer matrix and bin interface.

using PartonDensity, CSV, DelimitedFiles

# ## Transfer matrix

eMPp = 1 # e+/e- switch 0/1

# Read in an example integrated cross section
numbers_from_file = readdlm("data/HERAPDF20_NNLO_EIG_ePp.txt")

# List of integrated cross section values in 429 bins
integ_xsec = numbers_from_file[:,3]

# Corresponding list of expected event numbers
prediction = get_pred_N(integ_xsec, eMPp);

integ_xsec[153]

prediction[151]

# ## Bins in detector space
#
# The transfer matrix projects into a list of 153 bins with irregular x and Q2 spacing.
# The bin edges are provided in `data/zeus_bin_edges.jl` but can also be quickly and
# easily accessed with the helper function `get_bin_info` as shown below.

get_bin_info(10)


text:
struct CompositeChebyshevGrid
    degree::Int # Chebyshev degree
    x::Vector{Float64} # Chebyshev nodes
    w::Vector{Float64} # Chebyshev node weights

    np::Int # number of panel
    panel::Vector{Float64}

    ngrid::Int # size of the grid = (np-1)*degree
    grid::Vector{Float64}  # fine grid
    function CompositeChebyshevGrid(degree, panel)
        # fill each panel with N order Chebyshev nodes
        np = length(panel) # number of panels break points
        xc, wc = Interp.barychebinit(degree)
        fineGrid = zeros(Float64, (np - 1) * degree) # np break points have np-1 panels
        for i = 1:np-1
            a, b = panel[i], panel[i+1]
            fineGrid[(i-1)*degree+1:i*degree] = a .+ (b - a) .* (xc .+ 1.0) ./ 2
        end
        return new(degree, xc, wc, np, panel, (np - 1) * degree, fineGrid)
    end
end

function ωChebyGrid(dlrGrid, degree, print = true)
    Λ, rtol = dlrGrid.Λ, dlrGrid.rtol

    npo = Int(ceil(log(Λ) / log(2.0))) # subintervals on [0,lambda] in omega space (subintervals on [-lambda,lambda] is 2*npo)

    if dlrGrid.symmetry == :ph || dlrGrid.symmetry == :pha
        # Panel break points for the real frequency ∈ [0, Λ]
        # get exponentially dense near 0⁺
        pbpo = zeros(Float64, npo + 1)
        pbpo[1] = 0.0
        for i = 1:npo
            pbpo[i+1] = Λ / 2^(npo - i)
        end
    else #τ in (0, 1)
        ############ ω discretization ##################
        # Panel break points for the real frequency ∈ [-Λ, Λ]
        # get exponentially dense near 0⁻ and 0⁺
        pbpo = zeros(Float64, 2npo + 1)
        pbpo[npo+1] = 0.0
        for i = 1:npo
            pbpo[npo+i+1] = Λ / 2^(npo - i)
        end
        pbpo[1:npo] = -pbpo[2npo+1:-1:npo+2]
    end
    return CompositeChebyshevGrid(degree, pbpo)
end

function τChebyGrid(dlrGrid, degree, print = true)
    Λ, rtol = dlrGrid.Λ, dlrGrid.rtol

    npt = Int(ceil(log(Λ) / log(2.0))) - 2 # subintervals on [0,1/2] in tau space (# subintervals on [0,1] is 2*npt)

    if dlrGrid.symmetry == :ph || dlrGrid.symmetry == :pha
        ############# Tau discretization ##############
        # Panel break points for the imaginary time ∈ (0, 1)
        # get exponentially dense near 0⁺
        pbpt = zeros(Float64, npt + 1)
        pbpt[1] = 0.0
        for i = 1:npt
            pbpt[i+1] = 0.5 / 2^(npt - i)
        end

    else #τ in (0, 1)
        ############# Tau discretization ##############
        # Panel break points for the imaginary time ∈ (0, 1)
        # get exponentially dense near 0⁺ and 1⁻
        pbpt = zeros(Float64, 2npt + 1)
        pbpt[1] = 0.0
        for i = 1:npt
            pbpt[i+1] = 1.0 / 2^(npt - i + 1)
        end
        pbpt[npt+2:2npt+1] = 1 .- pbpt[npt:-1:1]
    end

    return CompositeChebyshevGrid(degree, pbpt)
end

"""
function preciseKernelT(dlrGrid, τ, ω, print::Bool = true)

    Calculate the kernel matrix(τ, ω) for given τ, ω grids

# Arguments
- τ: a CompositeChebyshevGrid struct or a simple one-dimensional array
- ω: a CompositeChebyshevGrid struct or a simple one-dimensional array
- print: print information or not
"""
function preciseKernelT(dlrGrid, τ, ω, print::Bool = true)
    # Assume τ.grid is particle-hole symmetric!!!
    # if (τ isa CompositeChebyshevGrid)
    #     @assert (τ.np - 1) * τ.degree == τ.ngrid
    # end
    τGrid = (τ isa CompositeChebyshevGrid) ? τ.grid : τ
    ωGrid = (ω isa CompositeChebyshevGrid) ? ω.grid : ω
    kernel = zeros(Float64, (length(τGrid), length(ωGrid)))
    symmetry = dlrGrid.symmetry

    if symmetry == :none && (τ isa CompositeChebyshevGrid) && (ω isa CompositeChebyshevGrid)
        #symmetrize K(τ, ω)=K(β-τ, -ω) for τ>0
        @assert isodd(τ.np) #symmetrization is only possible for odd τ panels
        halfτ = ((τ.np - 1) ÷ 2) * τ.degree
        kernel[1:halfτ, :] = Spectral.kernelT(Val(true), Val(symmetry), τ.grid[1:halfτ], ω.grid, 1.0, true)
        kernel[end:-1:(halfτ+1), :] = Spectral.kernelT(Val(true), Val(symmetry), τ.grid[1:halfτ], ω.grid[end:-1:1], 1.0, true)
        # use the fermionic kernel for both the fermionic and bosonic propagators
    else
        kernel = Spectral.kernelT(Val(dlrGrid.isFermi), Val(symmetry), τGrid, ωGrid, 1.0, true)
    end

    # print && println("=====  Kernel Discretization =====")
    # print && println("fine grid points for τ     = ", τGrid)
    # print && println("fine grid points for ω     = ", ωGrid)
    return kernel
end

function testInterpolation(dlrGrid, τ, ω, kernel, print = true)
    ############# test interpolation accuracy in τ #######
    τGrid = (τ isa CompositeChebyshevGrid) ? τ.grid : τ
    ωGrid = (ω isa CompositeChebyshevGrid) ? ω.grid : ω
    if τ isa CompositeChebyshevGrid
        τ2 = CompositeChebyshevGrid(τ.degree * 2, τ.panel)
        kernel2 = preciseKernelT(dlrGrid, τ2, ω, print)
        err = 0.0
        for ωi = 1:length(ωGrid)
            tmp = 0.0
            for i = 1:τ2.np-1
                for k = 1:τ2.degree
                    τidx = (i - 1) * τ2.degree + k
                    kaccu = kernel2[τidx, ωi]
                    kinterp = Interp.barycheb(τ.degree, τ2.x[k], kernel[(i-1)*τ.degree+1:i*τ.degree, ωi], τ.w, τ.x)
                    tmp = max(tmp, abs(kaccu - kinterp))
                end
            end
            err = max(err, tmp / maximum(kernel[:, ωi]))
        end
        print && println("Max relative L∞ error of kernel discretization in τ = ", err)
        @assert err < dlrGrid.rtol "Discretization error is too big! $err >= $(dlrGrid.rtol). Increase polynominal degree."
    end

    if ω isa CompositeChebyshevGrid
        ω2 = CompositeChebyshevGrid(ω.degree * 2, ω.panel)
        kernel2 = preciseKernelT(dlrGrid, τ, ω2)
        err = 0.0
        for τi = 1:length(τGrid)
            tmp = 0.0
            for i = 1:ω2.np-1
                for k = 1:ω2.degree
                    idx = (i - 1) * ω2.degree + k
                    kaccu = kernel2[τi, idx]
                    kinterp = Interp.barycheb(ω.degree, ω2.x[k], kernel[τi, (i-1)*ω.degree+1:i*ω.degree], ω.w, ω.x)
                    tmp = max(tmp, abs(kaccu - kinterp))
                end
            end
            err = max(err, tmp / maximum(kernel[τi, :]))
        end
        print && println("Max relative L∞ error of kernel discretization in ω = ", err)
        @assert err < dlrGrid.rtol "Discretization error is too big! $err >= $(dlrGrid.rtol). Increase polynominal degree."
    end
end

function preciseKernelΩn(dlrGrid, ω, print::Bool = true)
    function Freq2Index(isFermi, ωnList)
        if isFermi
            # ωn=(2n+1)π
            return [Int(round((ωn / π - 1) / 2)) for ωn in ωnList]
        else
            # ωn=2nπ
            return [Int(round(ωn / π / 2)) for ωn in ωnList]
        end
    end

    function nGrid(isFermi, symmetry, Λ)
        # generate n grid from a logarithmic fine grid
        degree = 100
        np = Int(round(log(10 * Λ) / log(2)))
        xc = [(i - 1) / degree for i = 1:degree]
        panel = [2^(i - 1) - 1 for i = 1:(np+1)]
        nGrid = zeros(Int, np * degree)
        for i = 1:np
            a, b = panel[i], panel[i+1]
            nGrid[(i-1)*degree+1:i*degree] = Freq2Index(isFermi, a .+ (b - a) .* xc)
        end
        unique!(nGrid)
        return symmetry == :none ? vcat(-nGrid[end:-1:2], nGrid) : nGrid
    end

    ωGrid = (ω isa CompositeChebyshevGrid) ? ω.grid : ω
    symmetry = dlrGrid.symmetry
    n = nGrid(dlrGrid.isFermi, symmetry, dlrGrid.Λ)

    nkernelFermi = Spectral.kernelΩ(Val(true), Val(symmetry), n, Float64.(ωGrid), 1.0, true)
    nkernelBose = Spectral.kernelΩ(Val(false), Val(symmetry), n, Float64.(ωGrid), 1.0, true)

    return n, nkernelFermi, nkernelBose
end


text:
module Oulu

greet() = print("Hello World!")

function myfunc(x)
    return x+2
end

end # module


text:
### Trait definitions
#####################
using Traits

# simple
@traitdef Tr1{X} begin
    fun1(X) -> Number   # this means a method with signature fun1(::X)
                        # returning a Number
end
@traitdef Tr2{X,Y} begin
    fun2(X,Y) -> Number
end
# subtrait
@traitdef Tr3{X,Y} <: Tr1{X}, Tr2{X,Y} begin
    fun3(X,Y,Int)
end
# with additional constraint on the types
@traitdef Tr4{X,Y} begin
    fun4(X,Y)
    @constraints begin
        # both Types need to start with the same letter:
        string(X.name)[1]==string(Y.name)[1]
    end
end

### Trait implementation
########################

# manual, i.e. just define the functions
fun1(x::Int) = 5x
@assert istrait(Tr1{Int})

# using @traitimpl
@traitimpl Tr1{Float64} begin
    fun1(x::Float64) = 7x # the explicit "::Float64" is needed at the moment
end
@assert istrait(Tr1{Float64})

# wrong usage of @traitimpl
try
    @traitimpl Tr1{Float32} begin
        fun1(x::Float64) = 7x # if the explicit type is wrong, it may error
    end
catch e
    println(e)
end

# This gives an error because supertypes have not been defined yet:
try
    eval(:(
    @traitimpl Tr3{Int, Int} begin
        fun3(x::Int, y::Int, t::Int) = x+y+t
    end))
catch e
    println(e) # TraitException("Not all supertraits of Tr3{Int64,Int64} are implemented.\nImplement them first.")
end

# this works:
@traitimpl Tr2{Int, Int} begin
    fun2(x::Int, y::Int) = x+y
end
@traitimpl Tr3{Int, Int} begin
    fun3(x::Int, y::Int, t::Int) = x+y+t
end
@traitimpl Tr4{Int, Int} begin
    fun4(x::Int, y::Int) = x+y
end

# This gives an error because constraints are not satisfied
try
    eval(:(
    @traitimpl Tr4{Int, Float64} begin
        fun4(x::Int, y::Int) = x+y
    end))
catch e
    println(e)  # ErrorException("assertion failed: istrait(Tr4{Int,Float64})")
end

### Trait functions
###################

@traitfn tf1{X, Y; Tr1{X}, Tr1{Y}}(a::X, b::Y) = fun1(a) + fun1(b)             # I
@traitfn tf1{X, Y; Tr1{X}, Tr1{Y}}(a::X, b::Y, c::Int) = fun1(a) + fun1(b) + c # II
@traitfn tf1{X, Y; Tr2{X,Y}}(a::X, b::Y) = fun2(a,b)                           # III
# Note that all the type-parameters are in the {} and that all
# arguments need a type parameter (a limitation of the
# macro-parser). This doesn't work:
#
# julia> @traitfn ttt1{X, Y; Tr1{X}, Tr1{Y}}(a::X, b::Y, c) = fun1(a) + fun1(b) + c
# ERROR: type Symbol has no field args
#
# But this works:
#
# julia> @traitfn ttt1{X, Y, Z; Tr1{X}, Tr1{Y}}(a::X, b::Y, c::Z) = fun1(a) + fun1(b) + c
# ttt1 (generic function with 6 methods)

# tf1 now dispatches on traits
@assert tf1(5.,6.)==77. # -> 77 ; dispatches to I because istrait(Tr1{Float64})
                        #         but not istrait(Tr2{Float64,Float64})
@assert tf1(5.,6.,77)==154. # -> 154. ; dispatches to II because of the extra argument

# Errors because of dispatch ambiguity:
try
    tf1(5,6)  # istrait(Tr1{Int}) and istrait(Tr2{Int,Int}) are both true!
catch e
    println(e)
end

# Implementing Tr1 for a type will make it work with tf1:
type MyType
    a::Int
end
try
    tf1(MyType(8), 9) # not implemented yet
catch e
    println(e)
end
@traitimpl Tr1{MyType} begin
    fun1(x::MyType) = x.a+9
end
@assert tf1(MyType(8), 9)==62 # -> 62 ; dispatches to I

### Generated machine code
##########################
f(x,y) = 7x + 7y
@code_llvm f(5.,6.)
@code_llvm tf1(5.,6.)
# It's the same.


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*--------------------------------------------------------------------
#* Author: Paula B. Staudt and Rafael P. Soares
#* $Id$
#*------------------------------------------------------------------
type batch_basic
	batch_basic()=begin
		NComp=outers.NComp
		new(
			DanaInteger(),
			DanaInteger(),
			fill(DanaReal(Dict{Symbol,Any}(
				:Brief=>"Stoichiometric Matrix"
			)),(NComp, NReac)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet stream",
				:PosX=>0,
				:PosY=>0,
				:Symbol=>"_{in}"
			)),
			fill(conc_mol(Dict{Symbol,Any}(
				:Brief=>"Components concentration"
			)),(NComp)),
			fill(reaction_mol(Dict{Symbol,Any}(
				:Brief=>"Reaction rates"
			)),(NReac)),
			temperature(Dict{Symbol,Any}(
				:Brief=>"Reactor temperature"
			)),
			pressure(Dict{Symbol,Any}(
				:Brief=>"Reactor pressure"
			)),
			volume(Dict{Symbol,Any}(
				:Brief=>"Reacting Volume"
			)),
			[
				:(diff(C*Vr) = Inlet.F*Inlet.z + sumt(stoic*r)*Vr),
			],
			[
				"Component Molar Balance",
			],
			[:NComp,:NReac,:stoic,],
			[:Inlet,:C,:r,:T,:P,:Vr,]
		)
	end
	NComp::DanaInteger
	NReac::DanaInteger
	stoic::Array{DanaReal}
	Inlet::stream
	C::Array{conc_mol}
	r::Array{reaction_mol}
	T::temperature
	P::pressure
	Vr::volume
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export batch_basic
function setEquationFlow(in::batch_basic)
	addEquation(1)
end
function atributes(in::batch_basic,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/batch"
	fields[:Brief]="Model of a batch reactor"
	fields[:Info]="== Assumptions ==
* isothermic
"
	drive!(fields,_)
	return fields
end
batch_basic(_::Dict{Symbol,Any})=begin
	newModel=batch_basic()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using NGWP, MTSG, Plots, LightGraphs, MultivariateStats
using LaTeXStrings

Nx, Ny = 7, 3
G = LightGraphs.grid([Nx, Ny]); N = nv(G);
L = Matrix(laplacian_matrix(G))
Q = incidence_matrix(G; oriented = true)
𝛌, 𝚽 = eigen(L); 𝚽 = 𝚽.*sign.(𝚽[1,:])'; # sign of DCT

∇𝚽 = Q' * 𝚽
W = 1.0 * adjacency_matrix(G)

grid2eig_ind = [1,2,3,6,8,12,15,4,5,7,9,13,16,18,10,11,14,17,19,20,21];
eig2grid_ind = sortperm(grid2eig_ind);
eig2dct = Array{Int64,3}(undef, Nx, Ny, 2);
for i = 1:Nx; for j = 1:Ny; eig2dct[i,j,1] = i-1; eig2dct[i,j,2] = j-1; end; end
eig2dct = reshape(eig2dct, (N, 2)); eig2dct = eig2dct[eig2grid_ind, :];

function grid7x3_mds_heatmaps(E, 𝚽; Nx = 7, Ny = 3, backend = :pyplot,
                              annotate_ind = 1:N, plotOrder = 1:N)
    # set up all heatmap plots' positions
    max_x = maximum(E[1, :]); min_x = minimum(E[1, :])
    width_x = max_x - min_x
    max_y = maximum(E[2, :]); min_y = minimum(E[2, :])
    width_y = max_y - min_y
    dx = 0.005 * width_x; dy = dx;
    xej = zeros(Nx, N); yej=zeros(Ny, N);
    a = 5.0; b = 7.0;
    for k = 1:N
        xej[:,k] = LinRange(E[1,k] - Ny * a * dx, E[1, k] + Ny * a * dx, Nx)
        yej[:,k] = LinRange(E[2,k] - a * dy, E[2, k] + a * dy, Ny)
    end

    # generate Grid7x3 2D MDS heatmaps plot
    if backend == :gr
        gr(dpi = 200)
    elseif backend == :pyplot
        pyplot(dpi = 200)
    elseif backend == :plotlyjs
        plotlyjs(dpi = 200)
    else
        @error("backend does not support $(backend)!")
    end

    plot()
    for k in plotOrder
        if k in annotate_ind
            heatmap!(xej[:, k], yej[:, k], reshape(𝚽[:, k], (Nx, Ny))', c = :viridis,
                     colorbar = false, ratio = 1, annotations = (xej[4, k],
                     yej[3, k] + b*dy, text(latexstring("\\varphi_{",
                     string(eig2dct[k, 1]), ",", string(eig2dct[k, 2]), "}"), 10)))
        else
            heatmap!(xej[:, k], yej[:, k], reshape(𝚽[:, k], (Nx, Ny))', c = :viridis,
                     colorbar = false, ratio = 1)
        end
    end
    plt = plot!(xlim = [min_x - 0.12 * width_x, max_x + 0.12 * width_x],
                ylim = [min_y - 0.16 * width_y, max_y + 0.16 * width_y],
                grid = false, clims = (-0.4, 0.4),
                xlab = "X₁", ylab = "X₂")
    return plt
end

# display(vcat((1:N)', eig2dct'))


text:
module Examples

using ..Radish

include("aob.jl")
include("rcp.jl")

end


text:
∑ = sum
using StatsFuns:logistic
sigmoid = logistic
using CSV:File
using DataFrames
using Printf:@sprintf

# DataTypes
using LinearAlgebra:Transpose
Numeric = Union{Int64,Float64}
NumericV = Union{Array{Int64,1},Array{Float64,1}}
NumericM = Union{Array{Int64,2},Array{Float64,2},Transpose{Float64,Array{Float64,2}},Transpose{Int64,Array{Int64,2}}}
arrayNumericV = Union{Array{Array{Int64,1},1},Array{Array{Float64,1},1}}
TupleNumericM = Tuple{NumericM, NumericM};

function maskBiasUnits(array::NumericV)::NumericV
    mask = ones(size(array));
    mask[1,1] = 0;
    maskedArray = array.*mask;
end

function maskBiasUnits(array::NumericM)::NumericM
    mask = ones(size(array));
    mask[:,1] .= 0;
    maskedArray = array.*mask;
end

function powerElements(array::Union{NumericV,NumericM}, power::Int)::NumericM
    # function defined to be broadcasted to elements of a tuple
    array = array.^power
end

function getLabels(y::NumericV)::NumericV
    labels = collect(BitSet(y));
    sort!(labels)
end

function sortArrayByCol(array::Array{Any,2}, col::Int)
    array = DataFrame(array);
    sort!(array, col)
    array = Matrix(array)
end

function getLayersQ(all_θ)::Int
    θ_matricesQ = length(collect(all_θ));
    layersQ = θ_matricesQ+1;
end

function nnClasificationFeedForward(all_θ, X, y, labels, layersQ::Int=getLayersQ(all_θ))
    labelsQ = length(labels);
    hx_and_indices = reshape([],0,labelsQ+1);
    for lb in labels
        ly = 1;
        lb_idx = findall(==(lb), y)
        a_ly = X[lb_idx,:];
        ly += 1;
        while ly ≤ layersQ
            # add bias units to a_ly_lb
            a_ly = [ones(size(a_ly)[1],1) a_ly];
            a_previous_ly = a_ly;
            θ_previous_ly = all_θ[ly-1];
            z_ly = a_previous_ly*θ_previous_ly';
            a_ly = sigmoid.(z_ly);
            ly += 1;
        end
        ly = 1;
        hx_lb = a_ly;
        hx_and_indices_lb = hcat(hx_lb, lb_idx);
        hx_and_indices = vcat(hx_and_indices, hx_and_indices_lb);
    end
    hx_and_indices = sortArrayByCol(hx_and_indices, labelsQ+1);
    hx = hx_and_indices[:,1:labelsQ];
end

function nnCostClasification(all_θ::TupleNumericM, hx, y::NumericV, labels::NumericV, λ::Numeric)::Numeric
    J = 0;
    for lb in labels
        lb_idx = findall(==(lb), y)
        hx_lb = hx[lb_idx,:];
        y_lb = vec(zeros(size(labels))); y_lb[lb] = 1;
        J += ∑(-((log.(hx_lb))*y_lb)-(log.(1 .- hx_lb)*(1 .- y_lb)));
    end
    casesQ = size(hx)[1];
    J = (1/casesQ)*J;
    # add regularization cost
    J += (λ/(2*casesQ))* ∑(∑.(powerElements.(maskBiasUnits.(all_θ),2)));
end

function calcSigmoidDerivate(z::Numeric)::Numeric
    sigmoid(z)*(1-sigmoid(z));
end

function calcSigmoidDerivate(z::NumericM)::NumericM
    sigmoid.(z).*(1 .- sigmoid.(z));
end

calcSigmoidGradient = calcSigmoidDerivate;

function calcNiceInitϵ(layerInputsQ::Int, layerOutputsQ::Int)::Numeric
    initϵ = √6/√(layerInputsQ+layerOutputsQ)
end

function genRandomInitialθ(layerInputsQ::Int, layerOutputsQ::Int, initϵ::Numeric=calcNiceInitϵ(layerInputsQ,layerOutputsQ))::NumericM
    initϵInt = floor(initϵ*100);
    randInitialθ = rand(-initϵInt:initϵInt,layerOutputsQ,layerInputsQ+1);
    randInitialθ = randInitialθ./100;
end

function substractElements(sub_array::Array{Any,1},vector::NumericV)
    sub_array .- vectors
end

function gen_initial_z(all_θ, labelsQ::Int)
    initial_z_singleLy = [];
    for θ in all_θ
        initial_z_singleLy = [initial_z_singleLy, ones(size(θ))];
    end
    initial_z_singleLy[1]=initial_z_singleLy[1][2];

    initial_z = [];
    for lb in labelsQ
        initial_z = [initial_z, initial_z_singleLy];
    end
    initial_z = initial_z[2]
end

function nnClasificationBackPropagation(hx, all_θ, layersQ::Int, y, labels)
    # calculate the gradients -i.e. partial derivates with respect to θ- of the cost (J)
    nodesQ = labelsQ = length(labels);
    casesQ = size(hx)[1];
    ly = layersQ;
    idx_v = collect(1:labelsQ);
    outputLayer = map(idx -> (vec(hx[:,idx])), idx_v); #Matrix of subarrays where each of them correspond to each row of hx'
    δ_ly =  map(idx -> allCasesOutputLayer[idx] .- y, idx_v);
    ly -= 1;
    while ly ≥ 2
        δ_previous_ly = δ_ly;
        θ_ly = all_θ[ly];
        z_ly =

end

#parameters
# csv files -> DataFrames -> arrays
X =  Matrix{Float64}(File("./X.csv", header=false)
    |> DataFrame);
y = vec(Matrix(File("./y.csv", header=false) |> DataFrame));
θ1 = Matrix{Float64}(File("./theta1.csv", header=false)
    |> DataFrame);
θ2 = Matrix{Float64}(File("./theta2.csv", header=false) |> DataFrame);
λ = 1;

#implementation
all_θ = (θ1, θ2);
labels = getLabels(y);
layersQ = getLayersQ(all_θ);
hx = nnClasificationFeedForward(all_θ, X, y, labels, layersQ);
J = nnCostClasification(all_θ, hx, y, labels, λ);
allCasesδ_ly = nnClasificationBackPropagation(hx, layersQ, y, labels);

max_stars_repo_name: luis-perez-one/courseStandfordML

max_issues_repo_name: luis-perez-one/courseStandfordML

max_forks_repo_name: luis-perez-one/courseStandfordML


text:
using Surrogates
using Flux
using Flux: @epochs

#1D
a = 0.0
b = 10.0
obj_1D = x -> 2*x+3
x = sample(10,0.0,10.,SobolSample())
y = obj_1D.(x);
model = Chain(Dense(1,1), first)
loss(x, y) = Flux.mse(model(x), y)
opt = Descent(0.01)
n_echos = 1
my_neural = NeuralSurrogate(x,y,a,b,model,loss,opt,n_echos)
add_point!(my_neural,8.5,20.0)
add_point!(my_neural,[3.2,3.5],[7.4,8.0])
val = my_neural(5.0)

#ND

lb = [0.0,0.0]
ub = [5.0,5.0]
x = sample(5,lb,ub, SobolSample())
obj_ND_neural(x) = x[1]*x[2];
y = obj_ND_neural.(x)
model = Chain(Dense(2,1), first)
loss(x, y) = Flux.mse(model(x), y)
opt = Descent(0.01)
n_echos = 1
my_neural = NeuralSurrogate(x,y,lb,ub,model,loss,opt,n_echos)
my_neural((3.5, 1.49))
my_neural([3.4,1.4])
add_point!(my_neural,(3.5,1.4),4.9)
add_point!(my_neural,[(3.5,1.4),(1.5,1.4),(1.3,1.2)],[1.3,1.4,1.5])

# Multi-output #98
f  = x -> [x^2, x]
lb = 1.0
ub = 10.0
x  = sample(5, lb, ub, SobolSample())
push!(x, 2.0)
y  = f.(x)
model = Chain(Dense(1,2))
loss(x, y) = Flux.mse(model(x), y)
surrogate = NeuralSurrogate(x,y,lb,ub,model,loss,opt,n_echos)

f  = x -> [x[1], x[2]^2]
lb = [1.0, 2.0]
ub = [10.0, 8.5]
x  = sample(20, lb, ub, SobolSample())
push!(x, (1.0, 2.0))
y  = f.(x)
model = Chain(Dense(2,2))
loss(x, y) = Flux.mse(model(x), y)
surrogate = NeuralSurrogate(x,y,lb,ub,model,loss,opt,n_echos)
surrogate((1.0, 2.0))
x_new = (2.0, 2.0)
y_new = f(x_new)
add_point!(surrogate, x_new, y_new)


text:
using Test
using FlightMechanics.Models
using FlightMechanics

llh = [deg2rad(49.996908), 0.000000, 9907.31]
xyz_ecef = [4114496.258, 0.0, 4870157.031]
xyz_earth = [0., 0., -9907.31]
llh_ref = LLHPosition(llh[1], llh[2], 0.0)

pllh = LLHPosition(llh...)

@test isapprox(get_llh(pllh), llh)
@test isapprox(
    get_xyz_earth(pllh, ref_point=llh_ref),
    xyz_earth
)
@test isapprox(get_xyz_ecef(pllh), xyz_ecef)
@test isapprox(get_height(pllh), llh[3])

xyz_earth = [0., 0., -9907.31]
pearth = EarthPosition(xyz_earth..., llh_ref)

@test isapprox(get_llh(pearth), llh)
@test isapprox(get_xyz_earth(pearth, ref_point=llh_ref), xyz_earth)
@test isapprox(get_xyz_ecef(pearth), xyz_ecef)
@test isapprox(get_height(pearth), -xyz_earth[3])

pecef = ECEFPosition(xyz_ecef...)

@test isapprox(get_llh(pecef), llh, rtol=1e-5)
@test isapprox(get_xyz_earth(pecef, ref_point=llh_ref), xyz_earth, rtol=1e-5)
@test isapprox(get_xyz_ecef(pecef), xyz_ecef)
@test isapprox(get_height(pecef), llh[3], atol=0.17)

# Start Earth come back from ECEF and LLH
ref_llh = LLHPosition(π/4, π/3, 0.0)
earth = EarthPosition(1000.0, 300.0, 5000.0, ref_llh)

llh = convert(LLHPosition(earth.ref_point.ellipsoid), earth)
earth_from_llh = convert(earth, llh)
@test isapprox(earth, earth_from_llh)

ecef = convert(ECEFPosition, earth)
earth_from_ecef = convert(earth, ecef)
@test isapprox(earth, earth_from_ecef)

# Start LLH come back from Earth and ECEF
llh = LLHPosition(π/4, π/3, 1000.0)

earth = convert(EarthPosition(ref_llh), llh)
llh_from_earth = convert(LLHPosition(WGS84), earth)
@test isapprox(llh, llh_from_earth)

ecef = convert(ECEFPosition, llh)
llh_from_ecef = convert(llh, ecef)
@test isapprox(llh, llh_from_ecef)

# Start ECEF come back from Earth and LLH
ecef = ECEFPosition(2.25e6, 3.91e6, 4.48e6)

llh = convert(LLHPosition(WGS84), ecef)
ecef_from_llh = convert(ECEFPosition, llh)
@test isapprox(ecef, ecef_from_llh)

earth = convert(EarthPosition(ref_llh), ecef)
ecef_from_earth = convert(ECEFPosition, earth)
@test isapprox(ecef, ecef_from_earth)

# Test getters


text:
using Test
using Records
using SparseArrays

include("test_listrecords.jl")
include("test_queuerecords.jl")
include("test_frames.jl")
include("test_conversions.jl")
include("test_io.jl")

max_stars_repo_name: tawheeler/Records

max_issues_repo_name: tawheeler/Records

max_forks_repo_name: tawheeler/Records


text:
"""
    AnovaResult <: AnovaEffect

A set of values for an Anova factor which has been tested

`name` - the name of this factor

`ss` - sum of squares

`df` - degrees of freedom

`ms` - mean square

`f` - the F statistic

`p` - the probability of a Type I error (incorrect rejection of null hypothesis)

`ω²` - the effect size
"""
struct AnovaResult <: AnovaEffect
    name::String
    ss::Float64
    df::Float64
    ms::Float64
    f::Float64
    p::Float64
    ω²::Float64
end

AnovaResult(factor, f, p) = AnovaResult(factor.name, factor.ss, factor.df, factor.ms, f, p, 0)

AnovaResult(result, ω²) = AnovaResult(result.name, result.ss, result.df, result.ms, result.f, result.p, ω²)

import Base.isapprox
isapprox(x::AnovaResult, y::AnovaResult) =
    x.name == y.name &&
    x.ss ≈ y.ss &&
    x.df == y.df &&
    x.ms ≈ y.ms &&
    x.f ≈ y.f &&
    x.p ≈ y.p &&
    (isnan(x.ω²) && isnan(y.ω²) || x.ω² ≈ y.ω²)