text:
"Parameters for right-hand side function of fully implicit Runge-Kutta methods."
mutable struct ParametersMidpointImplicit{DT, TT, D, S, ET <: NamedTuple} <: Parameters{DT,TT}
    equs::ET
    tab::Tableau{TT}
    Δt::TT

    t::TT
    q::Vector{DT}
    θ::Vector{DT}

    function ParametersMidpointImplicit{DT,D}(equs::ET, tab::Tableau{TT}, Δt::TT) where {DT, TT, D, ET <: NamedTuple}
        new{DT, TT, D, tab.s, ET}(equs, tab, Δt, zero(TT), zeros(DT,D), zeros(DT,D))
    end
end

@doc raw"""
Fully implicit Runge-Kutta integrator cache.

### Fields

* `q̃`: initial guess of solution
* `ṽ`: initial guess of vector field
* `s̃`: holds shift due to periodicity of solution
* `Q`: internal stages of solution
* `V`: internal stages of vector field
* `Θ`: implicit function of internal stages
* `F`: vector field of implicit function
"""
struct IntegratorCacheMidpointImplicit{DT,D,S} <: ODEIntegratorCache{DT,D}
    q::Vector{DT}
    v::Vector{DT}
    θ::Vector{DT}

    Q::Vector{DT}
    V::Vector{DT}
    F::Vector{DT}

    function IntegratorCacheMidpointImplicit{DT,D,S}() where {DT,D,S}
        new(zeros(DT,D), zeros(DT,D), zeros(DT,D),
            zeros(DT,D), zeros(DT,D), zeros(DT,D))
    end
end

function IntegratorCache{ST}(params::ParametersMidpointImplicit{DT,TT,D,S}; kwargs...) where {ST,DT,TT,D,S}
    IntegratorCacheMidpointImplicit{ST,D,S}(; kwargs...)
end

@inline CacheType(ST, params::ParametersMidpointImplicit{DT,TT,D,S}) where {DT,TT,D,S} = IntegratorCacheMidpointImplicit{ST,D,S}

"""
Implicit midpoint Runge-Kutta integrator.
"""
struct IntegratorMidpointImplicit{DT, TT, D, S, PT <: ParametersMidpointImplicit{DT,TT},
                                    ST <: NonlinearSolver{DT},
                                    IT <: InitialGuessODE{DT,TT}} <: AbstractIntegratorRK{DT,TT}
    params::PT
    solver::ST
    iguess::IT
    caches::CacheDict{PT}

    function IntegratorMidpointImplicit(params::ParametersMidpointImplicit{DT,TT,D,S}, solver::ST, iguess::IT, caches) where {DT,TT,D,S,ST,IT}
        new{DT, TT, D, S, typeof(params), ST, IT}(params, solver, iguess, caches)
    end

    function IntegratorMidpointImplicit{DT,D}(equations::NamedTuple, tableau::Tableau{TT}, Δt::TT) where {DT,TT,D}
        # get number of stages
        S = tableau.s

        @assert S == 1

        # create params
        params = ParametersMidpointImplicit{DT,D}(equations, tableau, Δt)

        # create cache dict
        caches = CacheDict(params)

        # create solver
        solver = create_nonlinear_solver(DT, D*S, params, caches)

        # create initial guess
        iguess = InitialGuessODE{DT,D}(get_config(:ig_extrapolation), equations[:v̄], Δt)

        # create integrator
        IntegratorMidpointImplicit(params, solver, iguess, caches)
    end

    # function IntegratorMidpointImplicit{DT,D}(v::Function, tableau::Tableau{TT}, Δt::TT; kwargs...) where {DT,TT,D}
    #     IntegratorMidpointImplicit{DT,D}(NamedTuple{(:v,)}((v,)), tableau, Δt; kwargs...)
    # end

    # function IntegratorMidpointImplicit{DT,D}(v::Function, h::Function, tableau::Tableau{TT}, Δt::TT; kwargs...) where {DT,TT,D}
    #     IntegratorMidpointImplicit{DT,D}(NamedTuple{(:v,:h)}((v,h)), tableau, Δt; kwargs...)
    # end

    function IntegratorMidpointImplicit(equation::Union{IODE{DT}, LODE{DT}}, tableau::Tableau{TT}, Δt::TT; kwargs...) where {DT,TT}
        IntegratorMidpointImplicit{DT, ndims(equation)}(get_functions(equation), tableau, Δt; kwargs...)
    end
end

@inline Base.ndims(::IntegratorMidpointImplicit{DT,TT,D,S}) where {DT,TT,D,S} = D

function initialize!(int::IntegratorMidpointImplicit, sol::AtomicSolutionODE)
    sol.t̄ = sol.t - timestep(int)

    equations(int)[:v](sol.t, sol.q, sol.v)

    initialize!(int.iguess, sol.t, sol.q, sol.v,
                            sol.t̄, sol.q̄, sol.v̄)
end

function update_params!(int::IntegratorMidpointImplicit, sol::AtomicSolutionPODE)
    # set time for nonlinear solver and copy previous solution
    int.params.t  = sol.t
    int.params.q .= sol.q
    equations(int)[:ϑ](sol.t, sol.q, sol.v, int.params.θ)
end

function initial_guess!(int::IntegratorMidpointImplicit{DT,TT}, sol::AtomicSolutionPODE{DT,TT},
                        cache::IntegratorCacheMidpointImplicit{DT}=int.caches[DT]) where {DT,TT}

    # compute initial guess for solution
    evaluate!(int.iguess, sol.q̄, sol.v̄, sol.q, sol.v, cache.q, cache.v, one(TT))
    int.solver.x .= cache.q
end

function compute_stages!(x::Vector{ST}, Q::Vector{ST}, V::Vector{ST}, F::Vector{ST},
                                        q::Vector{ST}, v::Vector{ST}, θ::Vector{ST},
                                        params::ParametersMidpointImplicit{DT,TT,D,S}) where {ST,DT,TT,D,S}

    # copy x to V
    q .= x

    # compute Q and V
    Q .= (q .+ params.q) ./ 2
    V .= (q .- params.q) ./ params.Δt

    # compute F = f(Q,V)
    params.equs[:f](params.t + params.Δt / 2, Q, V, F)

    # compute Θ = ϑ(q̄)
    params.equs[:ϑ](params.t + params.Δt, q, v, θ)
end

function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersMidpointImplicit{DT,TT,D,S},
                          caches::CacheDict) where {ST,DT,TT,D,S}
    # get cache for internal stages
    cache = caches[ST]

    # compute stages from nonlinear solver solution x
    compute_stages!(x, cache.Q, cache.V, cache.F, cache.q, cache.v, cache.θ, params)

    # compute b
    b .= cache.θ .- params.θ .- params.Δt .* cache.F
end

function integrate_step!(int::IntegratorMidpointImplicit{DT,TT}, sol::AtomicSolutionPODE{DT,TT},
                         cache::IntegratorCacheMidpointImplicit{DT}=int.caches[DT]) where {DT,TT}

    # update nonlinear solver parameters from atomic solution
    update_params!(int, sol)

    # compute initial guess
    initial_guess!(int, sol, cache)

    # reset atomic solution
    reset!(sol, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute vector field at internal stages
    compute_stages!(int.solver.x, cache.Q, cache.V, cache.F, cache.q, cache.v, cache.θ, int.params)

    # compute final update
    sol.q .= cache.q
    sol.p .= cache.θ

    # compute vector field for initial guess
    equations(int)[:v](sol.t, sol.q, sol.v)
    # update_vector_fields!(int.iguess, sol.t, sol.q, sol.v)
end


text:
export fit, fit!, Params

### PARAMETERS TYPE
abstract AbstractParams
Params(args...; kwargs...) = ProxGradParams(args...; kwargs...)

# default in-place fitting uses proximal gradient method
@compat function fit!(glrm::AbstractGLRM; kwargs...)
    kwdict = Dict(kwargs)
    if :params in keys(kwdict)
        return fit!(glrm, kwdict[:params]; kwargs...)
    else
        if isa(glrm.A,SparseMatrixCSC)
            # Default to sparse algorithm for a sparse dataset
            return fit!(glrm, SparseProxGradParams(); kwargs...)
        else
            # Classic proximal gradient method for non-sparse data
            return fit!(glrm, ProxGradParams(); kwargs...)
        end
    end
end

# fit without modifying the glrm object
function fit(glrm::AbstractGLRM, args...; kwargs...)
    X0 = Array(Float64, size(glrm.X))
    Y0 = Array(Float64, size(glrm.Y))
    copy!(X0, glrm.X); copy!(Y0, glrm.Y)
    X,Y,ch = fit!(glrm, args...; kwargs...)
    copy!(glrm.X, X0); copy!(glrm.Y, Y0)
    return X',Y,ch
end

ext: jl


text:
struct NullData end
const DEFAULT_DATA = Iterators.cycle((NullData(),))
Base.iterate(::NullData, i=1) = nothing
Base.length(::NullData) = 0

get_maxiters(data) = Iterators.IteratorSize(typeof(DEFAULT_DATA)) isa Iterators.IsInfinite ||
                     Iterators.IteratorSize(typeof(DEFAULT_DATA)) isa Iterators.SizeUnknown ?
                     typemax(Int) : length(data)

maybe_with_logger(f, logger) = logger === nothing ? f() : Logging.with_logger(f, logger)

function default_logger(logger)
    Logging.min_enabled_level(logger) ≤ ProgressLogging.ProgressLevel && return nothing
    if Sys.iswindows() || (isdefined(Main, :IJulia) && Main.IJulia.inited)
          progresslogger = ConsoleProgressMonitor.ProgressLogger()
    else
          progresslogger = TerminalLoggers.TerminalLogger()
    end
    logger1 = LoggingExtras.EarlyFilteredLogger(progresslogger) do log
        log.level == ProgressLogging.ProgressLevel
    end
    logger2 = LoggingExtras.EarlyFilteredLogger(logger) do log
        log.level != ProgressLogging.ProgressLevel
    end
    LoggingExtras.TeeLogger(logger1, logger2)
end

macro withprogress(progress, exprs...)
    quote
      if $progress
          $maybe_with_logger($default_logger($Logging.current_logger())) do
              $ProgressLogging.@withprogress $(exprs...)
            end
      else
        $(exprs[end])
      end
    end |> esc
end

decompose_trace(trace) = trace

function _map_optimizer_args(prob::OptimizationProblem, opt; kwargs...)
    __map_optimizer_args(prob, opt; kwargs...)
end

function _check_and_convert_maxiters(maxiters)
    if !(isnothing(maxiters)) && maxiters <= 0.0
        error("The number of maxiters has to be a non-negative and non-zero number.")
    elseif !(isnothing(maxiters))
        return convert(Int, maxiters)
    end
end

function _check_and_convert_maxtime(maxtime)
  if !(isnothing(maxtime)) && maxtime <= 0.0
      error("The maximum time has to be a non-negative and non-zero number.")
  elseif !(isnothing(maxtime))
      return convert(Float32, maxtime)
  end
end

function check_pkg_version(pkg::String,ver::String; branch::Union{String, Nothing}=nothing)
    deps = Pkg.dependencies()
    pkg_info = Dict{String, Pkg.Types.PackageInfo}()
    for (uuid, dep) in deps
        dep.is_direct_dep || continue
        dep.version === nothing && continue
        pkg_info[dep.name] = dep
    end

    return (isnothing(branch) | (pkg_info[pkg].git_revision == branch)) ? pkg_info[pkg].version >= VersionNumber(ver) : pkg_info[pkg].version > VersionNumber(ver)
end


text:
immutable FEMHeatIntegrator{T1,T2,T3,F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType}
  N::Int
  NT::Int
  dt::tType
  t::tType
  Minv::MinvType
  D::DiffType
  A::Array{uElType,2}
  freenode::Vector{Int}
  f::F1
  gD::F2
  gN::F3
  u::uType
  node::nodeType
  elem::Array{Int,2}
  area::AType
  bdnode::Vector{Int}
  mid::Array{uElType,3}
  dirichlet::Array{Int,2}
  neumann::Array{Int,2}
  islinear::Bool
  numvars::Int
  sqrtdt::Float64
  σ::F4
  noisetype::Symbol
  numiters::Int
  save_everystep::Bool
  timeseries::Vector{uType}
  ts::Vector{tType}
  solver::Symbol
  autodiff::Bool
  method::Symbol
  show_trace::Bool
  iterations::Int
  timeseries_steps::Int
  progressbar::Bool
  progress_steps::Int
  progressbar_name::String
end

@def femheat_footer begin
  u[bdnode] = gD(i*dt,node)[bdnode]
  if save_everystep && i%timeseries_steps==0
    push!(timeseries,copy(u))
    push!(ts,t)
  end
  if progressbar && i%progress_steps==0
    Juno.msg(prog,"dt="*string(dt))
    Juno.progress(prog,i/numiters)
  end
end

@def femheat_deterministicimplicitlinearsolve begin
  if solver==:Direct || solver==:Cholesky || solver==:QR || solver==:LU || solver==:SVD
    u[freenode,:] = lhs\(Dinv.*rhs(i,u))
  elseif solver==:CG
    for j=1:size(u,2)
      u[freenode,j],ch = cg!(u[freenode,j],lhs,Dinv.*rhs(i,u)[:,j]) # Requires Vector, need to change rhs
    end
  elseif solver==:GMRES
    for j=1:size(u,2)
      u[freenode,j],ch = gmres!(u[freenode,j],lhs,Dinv.*rhs(i,u)[:,j]) # Requires Vector, need to change rhs
    end
  end
end

@def femheat_stochasticimplicitlinearsolve begin
  dW = next(rands)
  if solver==:Direct || solver==:Cholesky || solver==:QR || solver==:LU || solver==:SVD
    u[freenode,:] = lhs\(Dinv.*rhs(i,u,dW))
  elseif solver==:CG
    for j=1:size(u,2)
      u[freenode,j],ch = cg!(u[freenode,j],lhs,Dinv.*rhs(i,u,dW)[:,j]) # Requires Vector, need to change rhs
    end
  elseif solver==:GMRES
    for j=1:size(u,2)
      u[freenode,j],ch = gmres!(u[freenode,j],lhs,Dinv.*rhs(i,u,dW)[:,j]) # Requires Vector, need to change rhs
    end
  end
end

@def femheat_implicitpreamble begin
  Dinv = D.^(-1)
  if solver==:Cholesky
    lhs = cholfact(lhs) # Requires positive definite, may be violated
  elseif solver==:LU
    lhs = lufact(lhs)
  elseif solver==:QR
    lhs = qrfact(lhs) #More stable, slower than LU
  elseif solver==:SVD
    lhs = svdfact(lhs)
  end
end

@def femheat_nonlinearsolvestochasticloop begin
  u = vec(u)
  uOld = copy(u)
  dW = next(rands)
  nlres = NLsolve.nlsolve((resid,u)->rhs!(i,u,resid,dW,uOld),uOld,autodiff=autodiff,method=method,show_trace=show_trace,iterations=iterations)
  u = nlres.zero
  if numvars > 1
    u = reshape(u,N,numvars)
  end
end

@def femheat_nonlinearsolvedeterministicloop begin
  u = vec(u)
  uOld = copy(u)
  nlres = NLsolve.nlsolve((resid,u)->rhs!(i,u,resid,uOld),uOld,autodiff=autodiff,method=method,show_trace=show_trace,iterations=iterations)
  u = nlres.zero
  if numvars > 1
    u = reshape(u,N,numvars)
  end
end

@def femheat_nonlinearsolvepreamble begin
  uOld = similar(vec(u))
end

@def femheat_deterministicpreamble begin
  @unpack N,NT,dt,t,Minv,D,A,freenode,f,gD,gN,u,node,elem,area,bdnode,mid,dirichlet,neumann,islinear,numvars,numiters,save_everystep,timeseries,ts,solver,autodiff,method,show_trace,iterations,timeseries_steps,progressbar,progress_steps, progressbar_name = integrator
  progressbar && (prog = Juno.ProgressBar(name=progressbar_name))
end

@def femheat_stochasticpreamble begin
  @unpack sqrtdt,σ,noisetype = integrator
  rands = getNoise(u,node,elem,noisetype=noisetype)
end

@def femheat_postamble begin
  progressbar && Juno.done(prog)
  if ts[end] != t
    push!(timeseries,copy(u))
    push!(ts,t)
  end
  u,timeseries,ts
end

function femheat_solve{F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType}(integrator::FEMHeatIntegrator{true,FEMDiffEqHeatEuler,false,F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType})
  @femheat_deterministicpreamble
  K = eye(N) - dt*Minv*D*A #D okay since numVar = 1 for linear
  @inbounds for i=1:numiters
    u[freenode,:] = K[freenode,freenode]*u[freenode,:] + (Minv*dt*quadfbasis((x)->f(t,x),(x)->gD(t,x),(x)->gN(t,x),
                A,u,node,elem,area,bdnode,mid,N,NT,dirichlet,neumann,islinear,numvars))[freenode,:]
    t += dt
    @femheat_footer
  end
  @femheat_postamble
end

function femheat_solve{F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType}(integrator::FEMHeatIntegrator{true,FEMDiffEqHeatEuler,true,F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType})
  @femheat_deterministicpreamble
  @femheat_stochasticpreamble
  K = eye(N) - dt*Minv*D*A #D okay since numVar = 1 for linear
  @inbounds for i=1:numiters
    dW = next(rands)
    u[freenode,:] = K[freenode,freenode]*u[freenode,:] + (Minv*dt*quadfbasis((x)->f(t,x),(x)->gD(t,x),(x)->gN(t,x),
                A,u,node,elem,area,bdnode,mid,N,NT,dirichlet,neumann,islinear,numvars))[freenode,:] +
                (sqrtdt.*dW.*Minv*quadfbasis((x)->σ(t,x),(x)->gD(t,x),(x)->gN(t,x),
                            A,u,node,elem,area,bdnode,mid,N,NT,dirichlet,neumann,islinear,numvars))[freenode,:]
    t += dt
    @femheat_footer
  end
  @femheat_postamble
end

function femheat_solve{F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType}(integrator::FEMHeatIntegrator{false,FEMDiffEqHeatEuler,false,F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType})
  @femheat_deterministicpreamble
  @inbounds for i=1:numiters
    u[freenode,:] = u[freenode,:] - D.*(dt*Minv[freenode,freenode]*A[freenode,freenode]*u[freenode,:]) + (Minv*dt*quadfbasis((x,u)->f(t,x,u),(x)->gD(t,x),(x)->gN(t,x),
            A,u,node,elem,area,bdnode,mid,N,NT,dirichlet,neumann,islinear,numvars))[freenode,:]
    t += dt
    @femheat_footer
  end
  @femheat_postamble
end

function femheat_solve{F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType}(integrator::FEMHeatIntegrator{false,FEMDiffEqHeatEuler,true,F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType})
  @femheat_deterministicpreamble
  @femheat_stochasticpreamble
  @inbounds for i=1:numiters
    dW = next(rands)
    u[freenode,:] = u[freenode,:] - D.*(dt*Minv[freenode,freenode]*A[freenode,freenode]*u[freenode,:]) + (Minv*dt*quadfbasis((x,u)->f(t,x,u),(x)->gD(t,x),(x)->gN(t,x),
                A,u,node,elem,area,bdnode,mid,N,NT,dirichlet,neumann,islinear,numvars))[freenode,:] +
                (sqrtdt.*dW.*(Minv*quadfbasis((x,u)->σ(t,x,u),(x)->gD(t,x),(x)->gN(t,x),
                            A,u,node,elem,area,bdnode,mid,N,NT,dirichlet,neumann,islinear,numvars)))[freenode,:]
    t += dt
    @femheat_footer
  end
  @femheat_postamble
end

function femheat_solve{F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType}(integrator::FEMHeatIntegrator{true,FEMDiffEqHeatImplicitEuler,true,F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType})
  @femheat_deterministicpreamble
  @femheat_stochasticpreamble
  K = eye(N) + dt*Minv*D*A #D okay since numVar = 1 for linear
  lhs = K[freenode,freenode]
  rhs(i,u,dW) = u[freenode,:] + (Minv*dt*quadfbasis((x)->f((i)*dt,x),(x)->gD((i)*dt,x),(x)->gN((i)*dt,x),A,u,node,elem,area,bdnode,mid,N,NT,dirichlet,neumann,islinear,numvars))[freenode,:] +
              (sqrtdt.*dW.*(Minv*quadfbasis((x)->σ((i-1)*dt,x),(x)->gD((i-1)*dt,x),(x)->gN((i-1)*dt,x),
                          A,u,node,elem,area,bdnode,mid,N,NT,dirichlet,neumann,islinear,numvars)))[freenode,:]
  @femheat_implicitpreamble
  @inbounds for i=1:numiters
    dW = next(rands)
    t += dt
    @femheat_stochasticimplicitlinearsolve
    @femheat_footer
  end
  @femheat_postamble
end

function femheat_solve{F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType}(integrator::FEMHeatIntegrator{true,FEMDiffEqHeatImplicitEuler,false,F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType})
  @femheat_deterministicpreamble
  K = eye(N) + dt*Minv*D*A #D okay since numVar = 1 for linear
  lhs = K[freenode,freenode]
  rhs(i,u) = u[freenode,:] + (Minv*dt*quadfbasis((x)->f((i)*dt,x),(x)->gD((i)*dt,x),(x)->gN((i)*dt,x),A,u,node,elem,area,bdnode,mid,N,NT,dirichlet,neumann,islinear,numvars))[freenode,:]
  @femheat_implicitpreamble
  @inbounds for i=1:numiters
    t += dt
    @femheat_deterministicimplicitlinearsolve
    @femheat_footer
  end
  @femheat_postamble
end

function femheat_solve{F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType}(integrator::FEMHeatIntegrator{true,FEMDiffEqHeatCrankNicholson,true,F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType})
  @femheat_deterministicpreamble
  @femheat_stochasticpreamble
  Km = eye(N) - dt*Minv*D*A/2 #D okay since numVar = 1 for linear
  Kp = eye(N) + dt*Minv*D*A/2 #D okay since numVar = 1 for linear
  lhs = Kp[freenode,freenode]
  rhs(i,u,dW) = Km[freenode,freenode]*u[freenode,:] + (Minv*dt*quadfbasis((x)->f((i-.5)*dt,x),(x)->gD((i-.5)*dt,x),(x)->gN((i-.5)*dt,x),A,u,node,elem,area,bdnode,mid,N,NT,dirichlet,neumann,islinear,numvars))[freenode,:] +
              (sqrtdt.*dW.*(Minv*quadfbasis((x)->σ((i-1)*dt,x),(x)->gD((i-1)*dt,x),(x)->gN((i-1)*dt,x),
                          A,u,node,elem,area,bdnode,mid,N,NT,dirichlet,neumann,islinear,numvars)))[freenode,:]
  @femheat_implicitpreamble
  @inbounds for i=1:numiters
    dW = next(rands)
    t += dt
    @femheat_stochasticimplicitlinearsolve
    @femheat_footer
  end
  @femheat_postamble
end

function femheat_solve{F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType}(integrator::FEMHeatIntegrator{true,FEMDiffEqHeatCrankNicholson,false,F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType})
  @femheat_deterministicpreamble
  Km = eye(N) - dt*Minv*D*A/2 #D okay since numVar = 1 for linear
  Kp = eye(N) + dt*Minv*D*A/2 #D okay since numVar = 1 for linear
  lhs = Kp[freenode,freenode]
  rhs(i,u) = Km[freenode,freenode]*u[freenode,:] + (Minv*dt*quadfbasis((x)->f((i-.5)*dt,x),(x)->gD((i-.5)*dt,x),(x)->gN((i-.5)*dt,x),A,u,node,elem,area,bdnode,mid,N,NT,dirichlet,neumann,islinear,numvars))[freenode,:]
  @femheat_implicitpreamble
  @inbounds for i=1:numiters
    t += dt
    @femheat_deterministicimplicitlinearsolve
    @femheat_footer
  end
  @femheat_postamble
end

function femheat_solve{F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType}(integrator::FEMHeatIntegrator{false,FEMDiffEqHeatSemiImplicitEuler,false,F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType}) #Incorrect for system with different diffusions
  @femheat_deterministicpreamble
  Dinv = D.^(-1)
  K = eye(N) + dt*Minv*A
  lhs = K[freenode,freenode]
  rhs(i,u) = u[freenode,:] + (Minv*dt*quadfbasis((x,u)->f((i)*dt,x,u),(x)->gD((i)*dt,x),(x)->gN((i)*dt,x),
              A,u,node,elem,area,bdnode,mid,N,NT,dirichlet,neumann,islinear,numvars))[freenode,:]
  @femheat_implicitpreamble
  @inbounds for i=1:numiters
    t += dt
    @femheat_deterministicimplicitlinearsolve
    @femheat_footer
  end
  @femheat_postamble
end

function femheat_solve{F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType}(integrator::FEMHeatIntegrator{false,FEMDiffEqHeatSemiImplicitEuler,true,F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType}) #Incorrect for system with different diffusions
  @femheat_deterministicpreamble
  @femheat_stochasticpreamble
  Dinv = D.^(-1)
  K = eye(N) + dt*Minv*A
  lhs = K[freenode,freenode]
  rhs(i,u,dW) = u[freenode,:] + (Minv*dt*quadfbasis((x,u)->f((i)*dt,x,u),(x)->gD((i)*dt,x),(x)->gN((i)*dt,x),
              A,u,node,elem,area,bdnode,mid,N,NT,dirichlet,neumann,islinear,numvars))[freenode,:] +
              (sqrtdt.*dW.*(Minv*quadfbasis((x,u)->σ((i-1)*dt,x,u),(x)->gD((i-1)*dt,x),(x)->gN((i-1)*dt,x),
                          A,u,node,elem,area,bdnode,mid,N,NT,dirichlet,neumann,islinear,numvars)))[freenode,:]
  @femheat_implicitpreamble
  @inbounds for i=1:numiters
    dW = next(rands)
    t += dt
    @femheat_stochasticimplicitlinearsolve
    @femheat_footer
  end
  @femheat_postamble
end

function femheat_solve{F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType}(integrator::FEMHeatIntegrator{false,FEMDiffEqHeatSemiImplicitCrankNicholson,false,F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType}) #Incorrect for system with different diffusions
  @femheat_deterministicpreamble
  Dinv = D.^(-1)
  Km = eye(N) - dt*Minv*A/2
  Kp = eye(N) + dt*Minv*A/2
  lhs = Kp[freenode,freenode]
  rhs(i,u) = Km[freenode,freenode]*u[freenode,:] + (Minv*dt*quadfbasis((x,u)->f((i-.5)*dt,x,u),(x)->gD((i-.5)*dt,x),(x)->gN((i-.5)*dt,x),
              A,u,node,elem,area,bdnode,mid,N,NT,dirichlet,neumann,islinear,numvars))[freenode,:]
  @femheat_implicitpreamble
  @inbounds for i=1:numiters
    t += dt
    @femheat_deterministicimplicitlinearsolve
    @femheat_footer
  end
  @femheat_postamble
end

function femheat_solve{F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType}(integrator::FEMHeatIntegrator{false,FEMDiffEqHeatSemiImplicitCrankNicholson,true,F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType}) #Incorrect for system with different diffusions
  @femheat_deterministicpreamble
  @femheat_stochasticpreamble
  Dinv = D.^(-1)
  Km = eye(N) - dt*Minv*A/2
  Kp = eye(N) + dt*Minv*A/2
  lhs = Kp[freenode,freenode]
  rhs(i,u,dW) = Km[freenode,freenode]*u[freenode,:] + (Minv*dt*quadfbasis((x,u)->f((i-.5)*dt,x,u),(x)->gD((i-.5)*dt,x),(x)->gN((i-.5)*dt,x),
              A,u,node,elem,area,bdnode,mid,N,NT,dirichlet,neumann,islinear,numvars))[freenode,:] +
              (sqrtdt.*dW.*(Minv*quadfbasis((x,u)->σ((i-1)*dt,x,u),(x)->gD((i-1)*dt,x),(x)->gN((i-1)*dt,x),
                          A,u,node,elem,area,bdnode,mid,N,NT,dirichlet,neumann,islinear,numvars)))[freenode,:]
  @femheat_implicitpreamble
  @inbounds for i=1:numiters
    dW = next(rands)
    t += dt
    @femheat_stochasticimplicitlinearsolve
    @femheat_footer
  end
  @femheat_postamble
end

function femheat_solve{F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType}(integrator::FEMHeatIntegrator{false,FEMDiffEqHeatImplicitEuler,false,F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType})
  @femheat_deterministicpreamble
  function rhs!(i,u,resid,uOld)
    u = reshape(u,N,numvars)
    uOld = reshape(uOld,N,numvars)
    resid = reshape(resid,N,numvars)
    resid[freenode,:] = u[freenode,:] - uOld[freenode,:] + D.*(dt*Minv[freenode,freenode]*A[freenode,freenode]*u[freenode,:]) -
    (Minv*dt*quadfbasis((x,u)->f((i)*dt,x,u),(x)->gD((i)*dt,x),(x)->gN((i)*dt,x),A,u,node,elem,area,bdnode,mid,N,NT,dirichlet,neumann,islinear,numvars))[freenode,:]
    u = vec(u)
    resid = vec(resid)
  end
  @femheat_nonlinearsolvepreamble
  @inbounds for i=1:numiters
    t += dt
    @femheat_nonlinearsolvedeterministicloop
    @femheat_footer
  end
  @femheat_postamble
end

function femheat_solve{F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType}(integrator::FEMHeatIntegrator{false,FEMDiffEqHeatImplicitEuler,true,F1,F2,F3,F4,uElType,uType,nodeType,AType,tType,DiffType,MinvType,StiffType})
  @femheat_deterministicpreamble
  @femheat_stochasticpreamble
  function rhs!(i,u,resid,dW,uOld)
    u = reshape(u,N,numvars)
    resid = reshape(resid,N,numvars)
    resid[freenode,:] = u[freenode,:] - uOld[freenode,:] + D.*(dt*Minv[freenode,freenode]*A[freenode,freenode]*u[freenode,:]) - (Minv*dt*quadfbasis((x,u)->f((i)*dt,x,u),(x)->gD((i)*dt,x),(x)->gN((i)*dt,x),A,u,node,elem,area,bdnode,mid,N,NT,dirichlet,neumann,islinear,numvars))[freenode,:] -(sqrtdt.*dW.*(Minv*quadfbasis((x,u)->σ((i)*dt,x,u),(x)->gD((i)*dt,x),(x)->gN((i)*dt,x),
                A,u,node,elem,area,bdnode,mid,N,NT,dirichlet,neumann,islinear,numvars)))[freenode,:]
    u = vec(u)
    resid = vec(resid)
  end
  @femheat_nonlinearsolvepreamble
  @inbounds for i=1:numiters
    t += dt
    @femheat_nonlinearsolvestochasticloop
    @femheat_footer
  end
  @femheat_postamble
end


text:
struct Simplex{T<:Number, U}
  vertices::Vector{Vertex{T,U}}
  function Simplex(vertices::Vector{Vertex{T,U}}) where {T<:Number, U}
    sort!(vertices, by=value)
    return new{T,U}(vertices)
  end
end
function Simplex(f::T, ic::AbstractVector{U}, initial_step::AbstractVector{V}
    ) where {T<:Function, U<:Number, V<:Number}
  if length(ic) != length(initial_step)
    throw(ArgumentError("ic, $ic must be same length as initial_step
                        $initial_step"))
  end
  dim = length(ic)
  positions = Vector{Vector{promote_type(U,V)}}()
  for i ∈ 1:dim+1
    x = [ic[j] + ((j == i) ? initial_step[j] : zero(V)) for j ∈ 1:dim]
    push!(positions, x)
  end
  return Simplex(f, positions)
end
function Simplex(f::T, positions::U
    ) where {T<:Function, W<:Number, V<:AbstractVector{W}, U<:AbstractVector{V}}
  if length(unique(length.(positions))) != 1
    throw(ArgumentError("All entries in positions $positions must be the same
                        length"))
  end
  dim = length(positions) - 1
  vertex = Vertex(positions[1], f(positions[1]))
  vertices = Vector{typeof(vertex)}()
  push!(vertices, vertex)
  map(i->push!(vertices, Vertex(positions[i], f(positions[i]))), 2:dim+1)
  return Simplex(vertices)
end

import Base: length, iterate, push!, iterate, getindex
import Base: eachindex, sort!, hash
Base.length(s::Simplex) = length(s.vertices)
Base.push!(s::Simplex, v::Vertex) = push!(s.vertices, v)
Base.iterate(s::Simplex) = iterate(s.vertices)
Base.iterate(s::Simplex, counter) = iterate(s.vertices, counter)
Base.getindex(s::Simplex, index) = s.vertices[index]
Base.sort!(s::Simplex; kwargs...) = sort!(s.vertices; kwargs...)
Base.hash(s::Simplex) = hash(s, hash(:Simplex))
Base.hash(s::Simplex, h::UInt64) = hash(hash.(s), h)

dimensionality(s::Simplex) = length(s) - 1

remove!(s::Simplex, v::Vertex) = filter!(x -> !≡(x, v), s.vertices)
remove!(s::Simplex, x::Vector) = deleteat!(s.vertices, x)

function getvertex(s::Simplex, i::Int)
  @assert 1 <= i <= length(s)
  @assert issorted(s)
  return s[i]
end

bestvertex(s::Simplex) = getvertex(s, 1)
worstvertex(s::Simplex) = getvertex(s, length(s))
secondworstvertex(s::Simplex) = getvertex(s, length(s)-1)

function centroidposition(s::Simplex)
  g(v) = ≡(v, worstvertex(s)) ? zero(position(v)) : position(v)
  return mapreduce(g, +, s) / (length(s) - 1)
end
function swapworst!(s::Simplex, forthis::Vertex)
  lengthbefore = length(s)
  remove!(s, worstvertex(s))
  push!(s, forthis)
  sort!(s, by=v->value(v))
  @assert length(s) == lengthbefore
  return nothing
end

function assessconvergence(simplex, xtol_abs, xtol_rel, ftol_abs, ftol_rel, stopval)
  value(bestvertex(simplex)) <= stopval && return :STOPVAL_REACHED

  allxtol = true
  for (vi, v) ∈ enumerate(simplex)
    for qi ∈ vi+1:length(simplex)
      q = getvertex(simplex, qi)
      thisxtol = true
      for (i, (pv, pq)) ∈ enumerate(zip(position(v), position(q)))
        thisxtol &= isapprox(pv, pq, rtol=xtol_rel[i], atol=xtol_abs[i])
      end
      allxtol &= thisxtol
    end
  end
  allxtol && return :FTOL_REACHED

  allftol = true
  for (vi, v) ∈ enumerate(simplex)
    for qi ∈ vi+1:length(simplex)
      q = getvertex(simplex, qi)
      allftol &= all(isapprox(value(v), value(q), rtol=ftol_rel, atol=ftol_abs))
    end
  end
  allftol && return :FTOL_REACHED

  return :CONTINUE
end


text:
# tests for ApproxManifoldProducts.jl

using ApproxManifoldProducts
using Test

include("basics.jl")
include("ex_1D.jl")
include("ex_2D_rot.jl")
include("testManifoldPartial.jl")
include("testManiProductBigSmall.jl")
include("testBasicManiProduct.jl")
include("testMarginalProducts.jl")
include("testMMD.jl")
include("testPartialProductSE2.jl")
include("basic_se3.jl")
include("testSymmetry.jl")

#


text:
using PyPlot, Images
using PyCall
@pyimport numpy

function DIC_sine_transform(img_name, lambda)

    # load in RGB data, convert to greyscale array
    raw_img = PyPlot.imread(img_name);
    img = convert(Image{Images.Gray},raw_img[:,:,1:3]);
    M = convert(Array,img);
    M = float(M);

    # take gradient of image, then process
    gX, gY = numpy.gradient(M);
    d1 = FFTW.r2r(gY,FFTW.RODFT00)
    nx, ny = size(d1)
    d3=zeros(nx,ny)
    d3=[(nx/pi)^2 * d1[j,i]/(lambda*j^2 + i^2) for i in 1:nx, j in 1:ny];
    dst = convert(Array{Float64,2},d3)
    e1 = FFTW.r2r(dst,FFTW.RODFT00)
    e1 = e1 / (2*nx)^2

    # clamp image range to 0 1, output as .png
    maxi = maximum(e1)
    mini = minimum(e1)
    clamped = 1 +  (-e1 + maxi) * 1 / ( mini-maxi);
    out_name = split(image_name,".")[1] * "_sin_l=$(@sprintf("%.2f", lambda))" * ".png"
    imwrite(grayim(clamped),out_name)

    return clamped

end

max_stars_repo_name: bcollyer/DICkit

max_issues_repo_name: bcollyer/DICkit

max_forks_repo_name: bcollyer/DICkit


text:
module SnFFT

#load the files included
include("Convolution.jl")
include("Correlation.jl")
include("Counter.jl")
include("Element.jl")
include("Examples.jl")
include("FFT.jl")
include("FFT_BL.jl")
include("FFT_SP.jl")
include("Function.jl")
include("Hook.jl")
include("IFFT.jl")
include("IFFT_BL.jl")
include("IFFT_P.jl")
include("Partitions.jl")
include("Partitions_BL.jl")
include("PartitionTree.jl")
include("PartitionTree_BL.jl")
include("YamanouchiSymbols.jl")
include("YamanouchiSymbols_BL.jl")
include("YoungsOrthogonalRepresentations.jl")
include("YoungsOrthogonalRepresentations_BL.jl")

export
	#core functions
	yor,
	yor_bl,
	sn_fft,
	sn_fft_bl,
	sn_fft_sp,
	sn_ifft,
	sn_ifft_bl,
	sn_ifft_p,
	sn_convolution,
	sn_correlation,

	#examples
	example1,
	example2,
	example3,
	example4,
	example5,
	example6,
	example7,
	example8,
	example_clustering,

	#functions from SnElement.jl
	sn_multiply,
	sn_inverse,
	sn_p,
	sn_cc,
	sn_at,
	sn_t,
	permutation_ccf,
	ccf_index,
	permutation_index,
	index_ccf,
	ccf_permutation,
	index_permutation,
	permutation_atf,
	yor_permutation,

	#functions fromSnFunction.jl
	snf,
	snf_bl,
	snf_sp,

	#other fucntions Examples.jl
	mallowsdistribution,
	preferencematrix,
	kendalldistance,
	permutation_string,
	partition_string,

	#other
	partitions

end


text:
struct EventRates{T <: EpidemicModel}
  exposure::Union{Nothing, Vector{Float64}}
  infection::Union{Nothing, Vector{Float64}}
  removal::Union{Nothing, Vector{Float64}}
  individuals::Int64

  function EventRates{T}(n::Int64) where T <: SEIR
    return new{T}(fill(0.0, n), fill(0.0, n), fill(0.0, n), n)
  end

  function EventRates{T}(n::Int64) where T <: SEI
    return new{T}(fill(0.0, n), fill(0.0, n), nothing, n)
  end

  function EventRates{T}(n::Int64) where T <: SIR
    return new{T}(nothing, fill(0.0, n), fill(0.0, n), n)
  end

  function EventRates{T}(n::Int64) where T <: SI
    return new{T}(nothing, fill(0.0, n), nothing, n)
  end
end

function Base.getindex(x::EventRates{T}, new_state::DiseaseState) where T <: EpidemicModel
  if new_state == State_E
    return x.exposure
  elseif new_state == State_I
    return x.infection
  elseif new_state == State_R
    return x.removal
  else
    @error "Unrecognized indexing disease state"
  end
end

function Base.getindex(x::EventRates{T}, new_states::Vector{DiseaseState}) where T <: EpidemicModel
  y = x[new_states[1]]
  for i=2:length(new_states)
    y = hcat(y, x[new_states[i]])
  end
  return y
end

function Base.sum(x::EventRates{T}) where T <: EpidemicModel
  return sum([sum(x[i]) for i in _state_progressions[T][2:end]])
end


text:
abstract type Animal end

function getName(animal::Animal)
    animal.name
end

struct Dog <: Animal
    name::String
end

struct Cat <: Animal
    name::String
end

println(Dog("Pochi") |> getName)
println(Cat("Tama") |> getName)
# Dog is Animal's subtype?
println(Dog <: Animal)

max_stars_repo_name: mkanenobu/trashbox

max_issues_repo_name: mkanenobu/trashbox

max_forks_repo_name: mkanenobu/trashbox


text:
"""
    NOCB(; context=Context())

Next observation carried backward (NOCB) iterates backwards through the `data` and fills
missing data with the next existing observation.

See also:
- [LOCF](@ref): Last Observation Carried Forward

WARNING: missing elements at the tail of the array may not be imputed if there is no
existing observation to carry backward. As a result, this method does not guarantee
that all missing values will be imputed.

# Keyword Arguments
* `context::AbstractContext`: A context which keeps track of missing data
  summary information

# Example
"""
struct NOCB <: Imputor
    context::AbstractContext
end

# TODO: Switch to using Base.@kwdef on 1.1
NOCB(; context=Context()) = NOCB(context)

function impute!(data::AbstractVector, imp::NOCB)
    imp.context() do c
        end_idx = findlast(c, data) - 1
        for i in end_idx:-1:firstindex(data)
            if ismissing!(c, data[i])
                data[i] = data[i+1]
            end
        end

        return data
    end
end


text:
using Distributed
using Images
using FileIO
using Mmap

# calculate a distributed fractal
width = 40000
height = 20000

rmin = -2.5
rmax = 1.5

imin = -1.25
imax = 1.25

iter = 500
epsilon = 0.25

mapfile = "zset-$(width)-$(height)"
s = open(mapfile)
n = read(s, Int)
println("size of zset is $n")
if n != (width * height)
	@warn("zset doesn't match dims")
	exit()
end

zset = Mmap.mmap(s, Array{ComplexF64,1}, n)
close(s)

num_procs =  0
hosts = ["robert@mbp2018.local"]

for i in 1:num_procs
	println("started proc $(i)")
	sshflags=`-i /Users/robert/src/julia/juliaSnippets/distributed/ms-keys`

	addprocs(hosts, tunnel=false, sshflags=sshflags, max_parallel=10, dir="/Users/robert/tmp", topology=:all_to_all)

end

# define mandel everywhere
include("mandel.jl")

println("starting pmap on $(workers())")
#
#generate image from results
# time measurements
print("starting...\n")
tStart=time()
#zset = Array{ComplexF64,1}(undef,0)
m_set = mandelbrot_set(zset, rmin, rmax, imin, imax, width, height, iter, epsilon)
tStop = time()

# write the image-file
img = colorview(RGB, m_set)
save("mandel.jpg", img)

print("done. took ", tStop-tStart, " seconds\n");

max_stars_repo_name: rabarar/juliaSnippets

max_issues_repo_name: rabarar/juliaSnippets

max_forks_repo_name: rabarar/juliaSnippets


text:
import gmsh

gmsh.initialize(ARGS)
gmsh.model.add("sphere_cut")
R = 1
R1 = 0.95
sph = gmsh.model.occ.addSphere(0,0,0, R, -1, 0, pi/2, pi/2)
b1 = gmsh.model.occ.addBox(R1,0,0, R,R,R)
b2 = gmsh.model.occ.addBox(0,R1,0, R,R,R)
b3 = gmsh.model.occ.addBox(0,0,R1, R,R,R)
gmsh.model.occ.cut([(3,sph)], [(3,b1), (3,b2), (3,b3)])
gmsh.model.occ.synchronize()

gmsh.model.removeEntities([(3,sph)])
gmsh.model.removeEntities([(2,2), (2,4), (2,6)], true)
gmsh.fltk.run()
gmsh.finalize()

max_stars_repo_name: Poofee/fastFEM

max_issues_repo_name: Christophe-Foyer/pyFEA

max_forks_repo_name: Christophe-Foyer/pyFEA


text:
using Documenter, SpeedyWeather

makedocs(
    format = Documenter.HTML(
    prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename="SpeedyWeather.jl",
    authors="M Klöwer",
    modules=[SpeedyWeather],
    pages=["Home"=>"index.md",
            "How to run SpeedyWeather.jl"=>"how_to_run_speedy.md",
            "Dynamical core"=>"dynamical_core.md",
            "Parameterizations"=>"parameterizations.md",
            "Boundary conditions"=>"boundary_conditions.md",
            "New model setups"=>"new_model_setups.md",
            "Function and type index"=>"functions.md"]
)

deploydocs(
    repo = "github.com/milankl/SpeedyWeather.jl.git",
    devbranch = "main"
)


text:
# copied and modified from https://github.com/tkf/IPython.jl/blob/master/test/install_dependencies.jl
# and https://github.com/JuliaPy/SymPy.jl/blob/master/test/install_dependencies.jl

# Adding Pkg in test/REQUIRE would be an error in 0.6.  Using
# Project.toml still has some gotchas.  So:

# Let PyCall.jl use Python interpreter from Conda.jl
# See: https://github.com/JuliaPy/PyCall.jl
ENV["PYTHON"] = "/opt/hostedtoolcache/Python/3.10.1/x64/bin/python"
Pkg.build("PyCall")


text:
module Param

export parse_parameters, SliceIterator

import YAML
using  DataStructures

"""
    parse_parameters(filename::String;
                     max_parameters::Union{Int, Nothing}=nothing)

Parse the parameters yml file to read information about partition parameters and output
sampling method.

Example Parameter file
======================
output:
  method: List
  params:
    bitstrings:
      - "01000"
      - "01110"
      - "10101"
"""
function parse_parameters(filename::String)
    param_dict = YAML.load_file(filename, dicttype=OrderedDict{String, Any})

    # parse the output method section of the parameter file
    method_params = OrderedDict{Symbol, Any}(Symbol(x[1]) => x[2] for x in param_dict["output"])
    method_params[:params] = OrderedDict{Symbol, Any}(Symbol(x[1]) => x[2] for x in method_params[:params])

    method_params
end

end


text:
#=
   Fair coin from a biased coin.

   http://cplint.eu/e/von_neumann_trick.swinb
   """
   If you have a biased coin, how to obtain a fair coin, i.e., from a coin that lands heads with
   probability p with p≠0.5, how to generate samples from {heads,tails} with P(heads)=P(tails)=0.5?

   John von Neuamnn gave the following solution

    1) Toss the coin twice.
    2) If you get HT return H, if you get TH return T
    3) Otherwise, return to 1 and start over

   In fact, if p is the probability of the biased coin landing heads, then the outcomes HT and TH are
   equiprobable with probability p(1−p). However, with probability pp+(1−p)(1−p)
   these results are not obtained. In this case, we simply repeat the process. The probability that we
   never get HT or TH is 0 so this procedure ends with certainty.

   See
   - https://en.wikipedia.org/wiki/Fair_coin#Fair_results_from_a_biased_coin
   - von Neumann, John (1951). "Various techniques used in connection with random digits".
     National Bureau of Standards Applied Math Series. 12: 36.
   """

   Distributions of variable biased_coin1 (num:0)
   1.00000 =>    7992  (0.799200)
   2.00000 =>    2008  (0.200800)

   Distributions of variable coin1 (num:0)
   1.00000 =>    5004  (0.500400)
   2.00000 =>    4996  (0.499600)

   Cf ~/blog/fair_coin_from_a_biased_coin.blog
      ~/webppl/fair_coin_from_a_biased_coin.wppl
=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function fair_coin_from_a_biased_coin()
    head = 1
    tail = 2

    n = 3
    function biased_coin(i)
        Categorical([0.8,0.2]) # [head,tail]
    end

    function coin(i)
        # StackOverflowError unless we limit the recursion.
        # But that don't help since we got a lot of 0's as the result....
        # This is fixed below...
        if i > 100
           return Dirac(0)
        end
        bi ~ biased_coin(i)
        bi1 ~ biased_coin(i+1)
        if bi != bi1
            Dirac(bi)
        else
            ci2 ~ coin(i+2)
            Dirac(ci2)
        end
    end

    biased_coin1 ~ biased_coin(1)
    coin1 ~ coin(1)
    # This is the fix. Not pretty though!
    true ~ Dirac(coin1 != 0)

end

model = fair_coin_from_a_biased_coin()

num_chains = 4

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

show_var_dist_pct(chns, :biased_coin1)
show_var_dist_pct(chns, :coin1)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using Documenter
using SimpleANOVA

DocMeta.setdocmeta!(
    SimpleANOVA,
    :DocTestSetup,
    :(using SimpleANOVA);
    recursive=true
)

makedocs(
    sitename = "SimpleANOVA.jl",
    pages = [
        "SimpleANOVA" => "index.md",
        "API" => "api.md"
    ],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    modules = [SimpleANOVA],
    strict = true,
    checkdocs = :exports
)

deploydocs(repo="github.com/BioTurboNick/SimpleANOVA.jl")


text:
using Trixi, OrdinaryDiffEq

dg = DGMulti(polydeg = 3, element_type = Quad(), approximation_type = GaussSBP(),
             surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs),
             volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha))

equations = CompressibleEulerEquations2D(1.4)

# Same as the default Trixi example, but v2 is a non-symmetric perturbation
function initial_condition_khi_nonsymmetric_perturbation(x, t, equations::CompressibleEulerEquations2D)
    # change discontinuity to tanh
    # typical resolution 128^2, 256^2
    # domain size is [-1,+1]^2
    slope = 15
     # discontinuous function with value 2 for -.5 <= x <= .5 and 0 elsewhere
    B = tanh(slope * x[2] + 7.5) - tanh(slope * x[2] - 7.5)
    rho = 0.5 + 0.75 * B
    v1 = 0.5 * (B - 1)
    # scale by (1 + .01 * sin(pi*x) * sin(pi*y)) to break symmetry
    v2 = 0.1 * sin(2 * pi * x[1]) * (1.0 + 1e-2 * sin(pi * x[1]) * sin(pi * x[2]))
    p = 1.0
    return prim2cons(SVector(rho, v1, v2, p), equations)
end
initial_condition = initial_condition_khi_nonsymmetric_perturbation

mesh = DGMultiMesh(dg, cells_per_dimension=(32, 32), periodicity=true)
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg)

tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()
alive_callback = AliveCallback(alive_interval=10)
analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg))
callbacks = CallbackSet(summary_callback,
                        analysis_callback,
                        alive_callback)

###############################################################################
# run the simulation

tsave = LinRange(tspan..., 100)
sol = solve(ode, RDPK3SpFSAL49(), abstol = 1.0e-7, reltol = 1.0e-7,
            dt = 1e-4,  save_everystep = false, saveat = tsave, callback = callbacks)

summary_callback() # print the timer summary


text:
# Estimation of Gaussian mixture models using EM

println("Importing modules ...")

using NumericExtensions
using Distributions
using MixtureModels
using Winston

### config

theta = [0,3,5,7] * (pi/4)
K = length(theta)

println("Generating data ...")

function gen_cluster(t, n)
	ct = cos(t)
	st = sin(t)

	x0 = randn(n) + 5.0
	y0 = randn(n) * 0.3

	x = x0 * ct - y0 * st
	y = x0 * st + y0 * ct

	[x y]'
end

obs = hcat([gen_cluster(t, 500) for t in theta]...)

### fit model

println("Run 1 ......")
r = fit_fmm(MvNormal{PDMat}, obs, K, fmm_em(;display=:iter))

for t in 2 : 5
	println("Run $t ......")
	rt = fit_fmm(MvNormal{PDMat}, obs, K, fmm_em(;display=:iter))
	if rt.objective > r.objective
		r = rt
	end
	println()
end

comps = r.mixture.components

### visualization

function gauss_circle(g::MvNormal, r::Float64)
	μ = g.μ

	t = linspace(0, 2pi, 300)
	x = [cos(t)'; sin(t)'] * r
	unwhiten!(g.Σ, x)

	Curve(x[1,:] + μ[1], x[2,:] + μ[2], "color", "red")
end

p = FramedPlot()
add(p, Points(obs[1,:], obs[2,:], "type", "dot", "color", "blue"))
for k in 1 : K
	add(p, gauss_circle(comps[k], 2.0))
end
Winston.display(p)

# block & wait for the signal to exit

println("Press enter to exit.")
readline(STDIN)


text:
# This file replicates the testsuite provided by MINPACK for the hybrj
# function. It consists in 14 functions with their analytical Jacobian.
# Julia replication taken from NLSOLVE.jl package
using LeastSquaresOptim, Printf,  LinearAlgebra, SparseArrays, Test

function rosenbrock()
    name = "rosenbrock"
    function f!(fvec, x)
        fvec[1] = 1 - x[1]
        fvec[2] = 10(x[2]-x[1]^2)
    end
    function g!(fjac, x)
        fjac[1,1] = -1
        fjac[1,2] = 0
        fjac[2,1] = -20x[1]
        fjac[2,2] = 10
    end
    x = [-1.2; 1.]
    return name, f!, g!, x
end

function powell_singular()
    name = "powell_singular"
    function f!(fvec, x)
        fvec[1] = x[1] + 10x[2]
        fvec[2] = sqrt(5)*(x[3] - x[4])
        fvec[3] = (x[2] - 2x[3])^2
        fvec[4] = sqrt(10)*(x[1] - x[4])^2
    end
    function g!(fjac, x)
        fill!(fjac, 0)
        fjac[1,1] = 1
        fjac[1,2] = 10
        fjac[2,3] = sqrt(5)
        fjac[2,4] = -fjac[2,3]
        fjac[3,2] = 2(x[2] - 2x[3])
        fjac[3,3] = -2fjac[3,2]
        fjac[4,1] = 2sqrt(10)*(x[1] - x[4])
        fjac[4,4] = -fjac[4,1]
    end
    x = [3.; -1.; 0.; 1.]
    return name, f!, g!, x
end

function powell_badly_scaled()
    name = "powell_badly_scaled"
    c1 = 1e4
    c2 = 1.0001
    function f!(fvec, x)
        fvec[1] = c1*x[1]*x[2] - 1
        fvec[2] = exp(-x[1]) + exp(-x[2]) - c2
    end
    function g!(fjac, x)
        fjac[1,1] = c1*x[2]
        fjac[1,2] = c1*x[1]
        fjac[2,1] = -exp(-x[1])
        fjac[2,2] = -exp(-x[2])
    end
    x = [0.; 1.]
    return name, f!, g!, x
end

function wood()
    name = "wood"
    c3 = 2e2
    c4 = 2.02e1
    c5 = 1.98e1
    c6 = 1.8e2

    function f!(fvec, x)
        temp1 = x[2] - x[1]^2
        temp2 = x[4] - x[3]^2
        fvec[1] = -c3*x[1]*temp1 - (1 - x[1])
        fvec[2] = c3*temp1 + c4*(x[2] - 1) + c5*(x[4] - 1)
        fvec[3] = -c6*x[3]*temp2 - (1 - x[3])
        fvec[4] = c6*temp2 + c4*(x[4] - 1) + c5*(x[2] - 1)
    end

    function g!(fjac, x)
        fill!(fjac, 0)
        temp1 = x[2] - 3x[1]^2
        temp2 = x[4] - 3x[3]^2
        fjac[1,1] = -c3*temp1 + 1
        fjac[1,2] = -c3*x[1]
        fjac[2,1] = -2*c3*x[1]
        fjac[2,2] = c3 + c4
        fjac[2,4] = c5
        fjac[3,3] = -c6*temp2 + 1
        fjac[3,4] = -c6*x[3]
        fjac[4,2] = c5
        fjac[4,3] = -2*c6*x[3]
        fjac[4,4] = c6 + c4
    end
    x = [-3.; -1; -3; -1]
    return name, f!, g!, x
end

function helical_valley()
    name = "helical_valley"
    tpi = 8*atan(1)
    c7 = 2.5e-1
    c8 = 5e-1

    function f!(fvec, x)
        if x[1] > 0
            temp1 = atan(x[2]/x[1])/tpi
        elseif x[1] < 0
            temp1 = atan(x[2]/x[1])/tpi + c8
        else
            temp1 = c7*sign(x[2])
        end
        temp2 = sqrt(x[1]^2+x[2]^2)
        fvec[1] = 10(x[3] - 10*temp1)
        fvec[2] = 10(temp2 - 1)
        fvec[3] = x[3]
    end

    function g!(fjac, x)
        temp = x[1]^2 + x[2]^2
        temp1 = tpi*temp
        temp2 = sqrt(temp)
        fjac[1,1] = 100*x[2]/temp1
        fjac[1,2] = -100*x[1]/temp1
        fjac[1,3] = 10
        fjac[2,1] = 10*x[1]/temp2
        fjac[2,2] = 10*x[2]/temp2
        fjac[2,3] = 0
        fjac[3,1] = 0
        fjac[3,2] = 0
        fjac[3,3] = 1
    end
    x = [-1.; 0; 0]
    return name, f!, g!, x
end

function watson(n::Integer)
    name = "watson"
    c9 = 2.9e1

    function f!(fvec, x)
        fill!(fvec, 0)
        for i = 1:29
            ti = i/c9
            sum1 = 0.0
            temp = 1.0
            for j = 2:n
                sum1 += (j-1)*temp*x[j]
                temp *= ti
            end
            sum2 = 0.0
            temp = 1.0
            for j = 1:n
                sum2 += temp*x[j]
                temp *= ti
            end
            temp1 = sum1 - sum2^2 - 1
            temp2 = 2*ti*sum2
            temp = 1/ti
            for k = 1:n
                fvec[k] += temp*(k-1-temp2)*temp1
                temp *= ti
            end
        end
        temp = x[2] - x[1]^2 - 1
        fvec[1] += x[1]*(1-2temp)
        fvec[2] += temp
    end

    function g!(fjac, x)
        fill!(fjac, 0)
        for i = 1:29
            ti = i/c9
            sum1 = 0.0
            temp = 1.0
            for j = 2:n
                sum1 += (j-1)*temp*x[j]
                temp *= ti
            end
            sum2 = 0.0
            temp = 1.0

            for j = 1:n
                sum2 += temp*x[j]
                temp *= ti
            end
            temp1 = 2(sum1 - sum2^2 - 1)
            temp2 = 2*sum2
            temp = ti^2
            tk = 1.0

            for k = 1:n
                tj = tk
                for j = k:n
                    fjac[k,j] += tj*((Float64(k-1)/ti - temp2)*(Float64(j-1)/ti - temp2) - temp1)
                    tj *= ti
                end
                tk *= temp
            end
        end
        fjac[1,1] += 6x[1]^2 - 2x[2] + 3
        fjac[1,2] -= 2x[1]
        fjac[2,2] += 1
        for k = 1:n
            for j = k:n
                fjac[j,k] = fjac[k,j]
            end
        end
    end
    x = zeros(n)
    return name, f!, g!, x
end

function chebyquad(n::Integer)
    name = "chebyquad"
    tk = 1/n

    function f!(fvec, x)
        fill!(fvec, 0)
        for j = 1:n
            temp1 = 1.0
            temp2 = 2x[j]-1
            temp = 2temp2
            for i = 1:n
                fvec[i] += temp2
                ti = temp*temp2 - temp1
                temp1 = temp2
                temp2 = ti
            end
        end
        iev = -1.0
        for k = 1:n
            fvec[k] *= tk
            if iev > 0
                fvec[k] += 1/(k^2-1)
            end
            iev = -iev
        end
    end

    function g!(fjac, x)
        for j = 1:n
            temp1 = 1.
            temp2 = 2x[j] - 1
            temp = 2*temp2
            temp3 = 0.0
            temp4 = 2.0
            for k = 1:n
                fjac[k,j] = tk*temp4
                ti = 4*temp2 + temp*temp4 - temp3
                temp3 = temp4
                temp4 = ti
                ti = temp*temp2 - temp1
                temp1 = temp2
                temp2 = ti
            end
        end
    end
    x = collect(1:n)/(n+1)
    return name, f!, g!, x
end

function brown_almost_linear(n::Integer)
    name = "brown_almost_linear"
    function f!(fvec, x)
        sum1 = sum(x) - (n+1)
        for k = 1:(n-1)
            fvec[k] = x[k] + sum1
        end
        fvec[n] = prod(x) - 1
    end

    function g!(fjac, x)
        fill!(fjac, 1)
        fjac[diagind(fjac)] .= 2
        prd = prod(x)
        for j = 1:n
            if x[j] == 0.0
                fjac[n,j] = 1.
                for k = 1:n
                    if k != j
                        fjac[n,j] *= x[k]
                    end
                end
            else
                fjac[n,j] = prd/x[j]
            end
        end
    end
    x = 0.5*ones(n)
    return name, f!, g!, x
end

function discrete_boundary_value(n::Integer)
    name = "discrete_boundary_value"
    h = 1/(n+1)

    function f!(fvec, x)
        for k = 1:n
            temp = (x[k] + k*h + 1)^3
            if k != 1
                temp1 = x[k-1]
            else
                temp1 = 0.0
            end
            if k != n
                temp2 = x[k+1]
            else
                temp2 = 0.0
            end
            fvec[k] = 2x[k] - temp1 - temp2 + temp*h^2/2
        end
    end

    function g!(fjac, x)
        for k = 1:n
            temp = 3*(x[k]+k*h+1)^2
            for j = 1:n
                fjac[k,j] = 0
            end
            fjac[k,k] = 2 + temp*h^2/2
            if k != 1
                fjac[k,k-1] = -1
            end
            if k != n
                fjac[k,k+1] = -1
            end
        end
    end
    x = collect(1:n)*h
    x = x .* (x .- 1)
    return name, f!, g!, x
end

function discrete_integral_equation(n::Integer)
    name = "discrete_integral_equation"
    h = 1/(n+1)

    function f!(fvec, x)
        for k = 1:n
            tk = k*h
            sum1 = 0.0
            for j = 1:k
                tj = j*h
                sum1 += tj*(x[j] + tj + 1)^3
            end
            sum2 = 0.0
            kp1 = k+1
            if n >= kp1
                for j = kp1:n
                    tj = j*h
                    sum2 += (1-tj)*(x[j] + tj + 1)^3
                end
            end
            fvec[k] = x[k] + h*((1-tk)*sum1 + tk*sum2)/2
        end
    end

    function g!(fjac, x)
        for k = 1:n
            tk = k*h
            for j = 1:n
                tj = j*h
                fjac[k,j] = h*min(tj*(1-tk), tk*(1-tj))*3(x[j] + tj + 1)^2/2
            end
            fjac[k,k] += 1
        end
    end

    x = collect(1:n)*h
    x = x .* (x .- 1)

    return name, f!, g!, x
end

function trigonometric(n::Integer)
    name = "trigonometric"
    function f!(fvec, x)
        for j = 1:n
            fvec[j] = cos(x[j])
        end
        sum1 = sum(fvec)
        for k = 1:n
            fvec[k] = n+k-sin(x[k]) - sum1 - k*fvec[k]
        end
    end

    function g!(fjac, x)
        for j = 1:n
            temp = sin(x[j])
            for k = 1:n
                fjac[k,j] = temp
            end
            fjac[j,j] = (j+1)*temp - cos(x[j])
        end
    end
    x = ones(n)/n
    return name, f!, g!, x
end

function variably_dimensioned(n::Integer)
    name = "variably_dimensioned"
    function f!(fvec, x)
        sum1 = 0.0
        for j = 1:n
            sum1 += j*(x[j]-1)
        end
        temp = sum1*(1+2sum1^2)
        for k = 1:n
            fvec[k] = x[k] - 1 + k*temp
        end
    end

    function g!(fjac, x)
        sum1 = 0.0
        for j = 1:n
            sum1 += j*(x[j]-1)
        end
        temp = 1 + 6sum1^2
        for k = 1:n
            for j = k:n
                fjac[k,j] = k*j*temp
                fjac[j,k] = fjac[k,j]
            end
            fjac[k,k] += 1
        end
    end
    x = collect(1:n)/n
    return name, f!, g!, x
end

function broyden_tridiagonal(n::Integer)
    name = "broyden_tridiagonal"
    function f!(fvec, x)
        for k = 1:n
            temp = (3-2x[k])*x[k]
            if k != 1
                temp1 = x[k-1]
            else
                temp1 = 0.0
            end
            if k != n
                temp2 = x[k+1]
            else
                temp2 = 0.0
            end
            fvec[k] = temp - temp1 - 2temp2 + 1
        end
    end

    function g!(fjac, x)
        fill!(fjac, 0)
        for k = 1:n
            fjac[k,k] = 3-4x[k]
            if k != 1
                fjac[k,k-1] = -1
            end
            if k != n
                fjac[k,k+1] = -2
            end
        end
    end
    x = -ones(n)
    return name, f!, g!, x
end

function broyden_banded(n::Integer)
    name = "broyden_banded"
    ml = 5
    mu = 1

    function f!(fvec, x)
        for k = 1:n
            k1 = max(1, k-ml)
            k2 = min(k+mu, n)
            temp = 0.0
            for j = k1:k2
                if j != k
                    temp += x[j]*(1+x[j])
                end
            end
            fvec[k] = x[k]*(2+5x[k]^2) + 1 - temp
        end
    end

    function g!(fjac, x)
        fill!(fjac, 0)
        for k = 1:n
            k1 = max(1, k-ml)
            k2 = min(k+mu, n)
            for j = k1:k2
                if j != k
                    fjac[k,j] = -(1+2x[j])
                end
            end
            fjac[k,k] = 2+15x[k]^2
        end
    end
    x = -ones(n)
    return name, f!, g!, x
end

for matrix in (:dense, :sparse)
    for (solver, solver_abbr) in ((LeastSquaresOptim.QR(), :qr), (LeastSquaresOptim.LSMR(), :lsmr))
        if (matrix == :sparse && typeof(solver) <: LeastSquaresOptim.QR)
            # tests below
            continue
        else
            for (optimizer, optimizer_abbr) in ((LeastSquaresOptim.Dogleg, :dl), (LeastSquaresOptim.LevenbergMarquardt, :lm))
                alltests = [rosenbrock();
                powell_singular(); powell_badly_scaled();
                wood();
                helical_valley();
                watson(6); watson(9);
                chebyquad(5); chebyquad(6); chebyquad(7); chebyquad(9);
                brown_almost_linear(10); brown_almost_linear(30); brown_almost_linear(40);
                discrete_boundary_value(10);
                discrete_integral_equation(1); discrete_integral_equation(10);
                trigonometric(10); variably_dimensioned(10);
                broyden_tridiagonal(10); broyden_banded(10)]
                for (name, f!, g!, x) in alltests
                    fcur = similar(x)
                    J = Array{Float64}(undef, length(x), length(x))
                    if matrix == :sparse
                        J = sparse(J)
                    end
                    nls = LeastSquaresProblem(x = x, y = fcur, f! = f!, J = J, g! = g!)
                    r = optimize!(nls, optimizer(solver))
                    @printf("%-6s %4s %2s %30s %5d %5d   %5d   %10e\n", matrix, solver_abbr, optimizer_abbr, name, r.iterations, r.f_calls, r.g_calls, r.ssr)
                    @test r.ssr <= 1e-3
                end
            end
        end
    end
end

# test sparse factorization
alltests =  [rosenbrock();
powell_singular(); powell_badly_scaled();
wood();
helical_valley();
watson(6); watson(9);
chebyquad(5); chebyquad(6); chebyquad(7); chebyquad(9);
brown_almost_linear(10); brown_almost_linear(30); brown_almost_linear(40);
discrete_boundary_value(10);
discrete_integral_equation(1); discrete_integral_equation(10);
trigonometric(10); variably_dimensioned(10);
broyden_tridiagonal(10); broyden_banded(10)]

for (name, f!, g!, x) in alltests
    fcur = similar(x)
    J = ones(length(x), length(x))
    J = sparse(J)
    fill!(J.nzval, 0)
    g!(J, x)
    # try because g! may change symbolic factorizations (see Julia issue #9906)
    try
        nls = LeastSquaresProblem(x = x, y = fcur, f! = f!, J = J, g! = g!)
        r = optimize!(nls, optimizer(solver))
        @test r.converged
        @test r.ssr <= 1e-3
        @printf("%-6s %4s %2s %30s %5d %5d   %5d   %10e\n", :sparse, :fact, :dl, name, r.iterations, r.f_calls, r.g_calls, r.ssr)
    catch y
        if !isa(y, Exception)
            throw(y)
        end
    end
end

alltests = [rosenbrock();
    powell_singular(); powell_badly_scaled();
    wood();
    helical_valley();
    watson(6);
    chebyquad(5); chebyquad(6); chebyquad(7); chebyquad(9);
    brown_almost_linear(10);# brown_almost_linear(30);
    discrete_boundary_value(10);
    discrete_integral_equation(1); discrete_integral_equation(10);
    trigonometric(10); variably_dimensioned(10);
    broyden_tridiagonal(10); broyden_banded(10)]
# test dense cholesky
for (optimizer, optimizer_abbr) in ((LeastSquaresOptim.Dogleg, :dl), (LeastSquaresOptim.LevenbergMarquardt, :lm))
    for (name, f!, g!, x) in alltests
        fcur = similar(x)
        J = Array{Float64}(undef, length(x), length(x))
        nls = LeastSquaresProblem(x = x, y = fcur, f! = f!, J = J, g! = g!)
        r = optimize!(nls, optimizer(LeastSquaresOptim.Cholesky()))
        @printf("%-6s %4s %2s %30s %5d %5d   %5d   %10e\n", :dense, :chol, optimizer_abbr, name, r.iterations, r.f_calls, r.g_calls, r.ssr)
        @test r.converged
        @test r.ssr <= 1e-3
    end
end

alltests = [rosenbrock();
  powell_singular(); powell_badly_scaled();
  wood();
  helical_valley();
  watson(6); watson(9);
  chebyquad(5); chebyquad(6); chebyquad(7); chebyquad(9);
  brown_almost_linear(10); brown_almost_linear(30); brown_almost_linear(40);
  discrete_boundary_value(10);
  discrete_integral_equation(1); discrete_integral_equation(10);
  trigonometric(10); variably_dimensioned(10);
  broyden_tridiagonal(10); broyden_banded(10)]
#test forwarddiff
for (optimizer, optimizer_abbr) in ((LeastSquaresOptim.Dogleg, :dl), (LeastSquaresOptim.LevenbergMarquardt, :lm))
    for (name, f!, g!, x) in alltests
        nls = LeastSquaresProblem(x = x, f! = f!, output_length = length(x))
        r = optimize!(nls, optimizer())
        @printf("%-10s %4s %2s %30s %5d %5d   %5d   %10e\n", :forwarddiff, :fact, optimizer_abbr, name, r.iterations, r.f_calls, r.g_calls, r.ssr)
        @test r.converged
        @test r.ssr <= 1e-3
    end
end

# check default are correctly
name, f!, g!, x = wood()
fcur = similar(x)
J = ones(length(x), length(x))
nls = LeastSquaresProblem(x = x, y = fcur, f! = f!, J = J, g! = g!)
result = optimize!(nls, show_trace = true)
@test result.optimizer == "Dogleg"
nls = LeastSquaresProblem(x = x, y = fcur, f! = f!, J = sparse(J), g! = g!)
result = optimize!(nls, show_trace = true)
@test result.optimizer == "LevenbergMarquardt"

# bigfloat
function newf(x)
    out = similar(x)
    out[1] = 1 - x[1]
    out[2] = 10(x[2]-x[1]^2)
    return out
end
nls = optimize(newf, [BigFloat(0.0), BigFloat(0.0)], Dogleg())
nls = optimize(newf, [BigFloat(0.0), BigFloat(0.0)], LevenbergMarquardt())

max_stars_repo_name: matthieugomez/LeastSquares

max_issues_repo_name: matthieugomez/LeastSquares

max_forks_repo_name: matthieugomez/LeastSquares


text:
# A server and a client help each other counting to 10.
using WebSockets
import WebSockets.with_logger
import WebSockets.string_with_env_ws
const COUNTTO = 10
const PORT = 8090
addsubproto("count")
# This overloading only affects the
# first argument to logging macros, and only when
# logging with WebSocketLogger
function string_with_env_ws(env, req::WebSockets.Request)
    iob = IOBuffer()
    ioc = IOContext(iob, env)
    printstyled(ioc, "Request: ", color= :yellow)
    printstyled(ioc, req.method, color = :bold)
    print(ioc, " ")
    printstyled(ioc, req.target, color = :cyan)
    subprot = subprotocol(req)
    if subprot != ""
        print(ioc, "\n")
        orig = WebSockets.origin(req)
        print(ioc, "\nOrigin: " * orig * " Subprotocol: " * subprot)
    end
    if length(req.headers) > 0
        print(ioc, " ")
        printstyled(ioc, "\tHeaders: ", color = :cyan)
        for (name, value) in req.headers
            print(ioc, "$name => $value \t")
        end
    end
    String(take!(iob))
end

httphandler(req::WebSockets.Request) = WebSockets.Response(200, "OK")

function coroutine_count(ws)
    @wslog ("Enter coroutine, ", ws)
    success = true
    protocolfollowed = true
    counter = 0
    # Server sends 1
    if ws.server
        counter += 1
        if !writeguarded(ws, Array{UInt8}([counter]))
            @warn ws, " could not write first message"
            protocolfollowed = false
        end
    end

    while isopen(ws) && protocolfollowed && counter < COUNTTO && success
        protocolfollowed = false
        data, success = readguarded(ws)
        if success
            OUTDATA = data
            if length(data) == 1
                if data[1] == counter + 1
                    protocolfollowed = true
                    counter += 2
                    @wslog ws, " Counter: ", counter
                    writeguarded(ws, Array{UInt8}([counter]))
                else
                    @error ws, " unexpected: ", counter
                end
            else
                protocolfollowed = false
                @warn (ws, " wrong message length: "), length(data)
                @wslog "Data is type ", typeof(data)
                @wslog data
            end
        end
    end
    if protocolfollowed
        @wslog (ws, " finished counting to $COUNTTO, exiting")
    else
        @wslog "Exiting ", ws, " Counter: ", counter
    end
end

function gatekeeper(req, ws)
    @wslog req
    if subprotocol(req) != ""
        coroutine_count(ws)
    else
        @wslog "No subprotocol"
    end
end

function serve_task(logger= WebSocketLogger(stderr, WebSockets.Logging.Debug))
    server = WebSockets.ServerWS(httphandler, gatekeeper)
    task = @async   with_logger(logger) do
        WebSockets.serve(server, port = PORT)
    end
    @info "http://localhost:$PORT"
    return server, task
end

server, sertask = serve_task()

clitask = @async with_logger(WebSocketLogger(stderr, WebSockets.Logging.Debug)) do
        WebSockets.open(coroutine_count, "ws://localhost:$PORT", subprotocol = "count")
    end

@async begin
    sleep(5)
    println("Time out 5 s, closing server")
    put!(server.in, "Close")
    nothing
end


text:
module TestDoctest

import FGenerators
using Documenter: doctest
using Test

@testset "doctest" begin
    doctest(FGenerators; manual = false)
end

end  # module


text:
# ordered dict

import Base: haskey, get, get!, getkey, delete!, push!, pop!, empty!,
             setindex!, getindex, length, isempty, start,
             next, done, keys, values, setdiff, setdiff!,
             union, union!, intersect, isequal, filter, filter!,
             hash, eltype

# This is just a simple wrapper around a HashDict, which is a modified
# implementation of the Dict implementation in Base allowing ordering
# information to be maintained.

# In particular, the HashDict stored in an OrderedDict is a
# HashDict{K,V,Ordered}

# A HashDict{K,V,Unordered} would be equivalent to Base.Dict

immutable OrderedDict{K,V} <: Associative{K,V}
    d::HashDict{K,V,Ordered}

    OrderedDict() = new(HashDict{K,V,Ordered}())
    OrderedDict(kv) = new(HashDict{K,V,Ordered}(kv))
    if VERSION >= v"0.4.0-dev+980"
        OrderedDict(ps::Pair{K,V}...) = new(HashDict{K,V,Ordered}(ps...))
    end
    #OrderedDict(ks,vs) = new(HashDict{K,V,Ordered}(ks,vs))
end

OrderedDict() = OrderedDict{Any,Any}()
OrderedDict(kv::@compat Tuple{}) = OrderedDict()
OrderedDict(kv) = ordered_dict_with_eltype(kv, eltype(kv))

ordered_dict_with_eltype{K,V}(kv, ::Type{@compat Tuple{K,V}}) = OrderedDict{K,V}(kv)
ordered_dict_with_eltype(kv, t) = OrderedDict{Any,Any}(kv)

if VERSION >= v"0.4.0-dev+980"
    OrderedDict{K,V}(ps::Pair{K,V}...)              = OrderedDict{K,V}(ps...)
    OrderedDict{K,V}(kv::@compat Tuple{Vararg{Pair{K,V}}})         = OrderedDict{K,V}(kv)
    OrderedDict{K}(kv::@compat Tuple{Vararg{Pair{K}}})             = OrderedDict{K,Any}(kv)
    OrderedDict{V}(kv::@compat Tuple{Vararg{Pair{TypeVar(:K),V}}}) = OrderedDict{Any,V}(kv)
    OrderedDict(kv::@compat Tuple{Vararg{Pair}})                   = OrderedDict{Any,Any}(kv)

    OrderedDict{K,V}(kv::AbstractArray{Pair{K,V}}) = OrderedDict{K,V}(kv)

    ordered_dict_with_eltype{K,V}(kv, ::Type{Pair{K,V}}) = OrderedDict{K,V}(kv)
end

copy(d::OrderedDict) = OrderedDict(d)

## Functions

## Most functions are simply delegated to the wrapped HashDict

@delegate OrderedDict.d [ haskey, get, get!, getkey, delete!, pop!,
                         empty!, setindex!, getindex,
                         length, isempty, start, next, done, keys,
                         values ]

sizehint(d::OrderedDict, sz::Integer) = (sizehint(d.d, sz); d)

if VERSION >= v"0.4.0-dev+980"
    push!(d::OrderedDict, kv::Pair) = (push!(d.d, kv); d)
    push!(d::OrderedDict, kv::Pair, kv2::Pair) = (push!(d.d, kv, kv2); d)
    push!(d::OrderedDict, kv::Pair, kv2::Pair, kv3::Pair...) = (push!(d.d, kv, kv2, kv3...); d)
end

push!(d::OrderedDict, kv) = (push!(d.d, kv); d)
push!(d::OrderedDict, kv, kv2...) = (push!(d.d, kv, kv2...); d)

similar{K,V}(d::OrderedDict{K,V}) = OrderedDict{K,V}()
in{T<:OrderedDict}(key, v::Base.KeyIterator{T}) = key in keys(v.dict.d.d)

max_stars_repo_name: ambujone/DIDebug

max_issues_repo_name: ambujone/DIDebug

max_forks_repo_name: ambujone/DIDebug


text:
using Roadnet_MDP
using POMDPToolbox
using MCTS

function visualize(mdp::roadnet_with_pursuer,s::roadnet_pursuer_state,policy::MCTS.MCTSPlanner;ms::Int64=10,fldr::String="default_movie_folder",fname::String="frame",max_steps::Int64=1000)
    hist = HistoryRecorder(max_steps=max_steps)
    hist = simulate(hist, mdp, policy, s)

    ## check for directory
    if !(fldr in readdir())
        mkdir(fldr)
    end
    if !("images" in readdir(fldr))
        mkdir(joinpath(fldr,"images"))
    end

    img_path = string(joinpath(fldr,"images"))

    #clear out old files
    if !(isempty(readdir(img_path)))
        display("removing pre-existing files")
        run(`rm $img_path/$(readdir(img_path))`)
    end

    it = 1
    println("processing $(length(hist)) steps")
    for (s,a,sp) in eachstep(hist, "s,a,sp")
        #  a = action(policy,s)
        a_int = POMDPs.action_index(mdp,a)
        if POMDPs.action_index(mdp,a) <= length(neighbors(mdp.road_net,s.node))
            a_loc = neighbors(mdp.road_net,s.node)[a_int]
        else
            a_loc = empty!([1])
        end
        println("taking action: $a, from node $(s.node) to get to node: $(sp.node) pursuer going to: $(sp.pnode)")
        #  display_network(mdp.road_net,evader_locs=[s.node],pursuer_locs=[s.pnode],action_locs=[a_loc],fname=string(joinpath(img_path,fname),it))
        #  display_network(mdp.road_net,evader_locs=[s.node],pursuer_locs=[s.pnode],fname=string(joinpath(img_path,fname),prepend_zeros("$it",tot_chars=length(hist))))
        display_network(mdp.road_net,evader_locs=[s.node],pursuer_locs=[s.pnode],fname=string(joinpath(img_path,fname),prepend_zeros("$it",tot_chars=1)))

        it += 1
        #  return a
    end;
    # original command line:
    # for file in *.svg; do inkscape $file -z -d=400 -e ${file%svg}png; done
    # -d is dpi, -z means no gui, -e is export as png
    temp_log = mktemp()
    for file in searchdir(img_path,".svg")
        f_split = splitext(file)
        # pipeline function allows to pipe output of command to a temp file
        run(pipeline(`inkscape $(joinpath(img_path,file)) -d=400 -z -e $(joinpath(img_path,f_split[1])).png`, stdout= "$(temp_log[1])"));
    end
    #  MAKE SURE TO SORT THESE BEFORE MAKING VIDEO
    run(pipeline(`ffmpeg -r 0.5 -y -i $img_path/frame%d.png -crf 25 $fldr/movie.avi`, stdout="$(temp_log[1])",stderr="$(temp_log[1])"))
    display("wrote movie to $(fldr)/movie.avi")
end

searchdir(path,key) = filter(x->contains(x,key),readdir(path)) # found this function on a forum

function prepend_zeros(s::String;tot_chars::Int64=5)
    # Help keep the files in alphanumerical order
    current_length = length(s)
    num_zeros = tot_chars-current_length
    if num_zeros > 0
        zeros = repeat("0",tot_chars-current_length)
    else
        zeros = ""
    end
    return zeros*s
end

function print_hist(mdp::roadnet_with_pursuer,hist::MDPHistory)
    hist_log = ""
    it = 1
    for (s,a,r,sp) in eachstep(hist, "s,a,r,sp")
        #  a = action(policy,s)
        a_int = POMDPs.action_index(mdp,a)
        if POMDPs.action_index(mdp,a) <= length(neighbors(mdp.road_net,s.node))
            a_loc = neighbors(mdp.road_net,s.node)[a_int]
        else
            a_loc = empty!([1])
        end
        h_str = "$it --- taking action: $a, from node $(s.node) to get to node: $(sp.node) pursuer going to: $(sp.pnode) rwd: $r"
        hist_log = "$hist_log \n $h_str"
        it += 1
    end
    return "$hist_log \nreward: $(undiscounted_reward(hist)), discounted: $(discounted_reward(hist))"
end

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
using RegexTools
using Test

@testset "RegexTools.jl" begin

    @testset "hex_escape" begin
        for char in Char.(0:126)
            rawstr = string(char)
            regstr = RegexTools.hex_escape(rawstr)

            reg = Regex(regstr)
            m = match(reg, rawstr)
            @test !isnothing(m)
        end
    end

end


text:
export Mesh

"""
    `Mesh(args)`

Constructs a mesh of collocation points and associated Fourier modes.

# Arguments
Can be either
- `xmin , xmax , N `; or
- `L , N ` (same as above with `xmin=-L` and `xmax=L`); or
- `param :: NamedTuple`, param contains `N` and `L`, then same as above; or
- `x` a vector of regularly spaced collocation points`.

The mesh as `N` collocation points regularly spaced between `xmin` (included) and `xmax` (excluded)

# Return values
`m=Mesh(args)` is of parametric type and offers
with
- `m.N `: number of collocation points and Fourier modes;
- `m.xmin`: minimum of the mesh;
- `m.xmax`: maximum of the mesh;
- `m.dx`: distance between two collocation points;
- `m.x`: the vector of collocation points;
- `m.kmin`: minimum of Fourier modes;
- `m.kmax`: maximum of Fourier modes;
- `m.dk`: distance between two Fourier modes;
- `m.k`: the vector of Fourier modes.

"""
struct Mesh

    N    :: Int64
    xmin :: Float64
    xmax :: Float64
    dx   :: Float64
    x    :: Vector{Float64}
    kmin :: Float64
    kmax :: Float64
    dk   :: Float64
    k    :: Vector{Float64}

    function Mesh( xmin , xmax , N :: Int64)

        dx   = (xmax-xmin)/N
        x    = zeros(Float64, N)
        x   .= range(xmin, stop=xmax, length=N+1)[1:end-1]
        dk   = 2π/(N*dx)
        kmin = -N/2*dk
        kmax = (N/2-1)*dk
        k    = zeros(Float64, N)
        k   .= dk .* vcat(0:(N-1)÷2, -N÷2:-1)

        new( N, xmin, xmax, dx, x, kmin, kmax, dk, k)

    end

    function Mesh(L , N :: Int64)

        xmin = - L
        xmax =   L
        N    =   N

        Mesh( xmin, xmax, N)

    end

    function Mesh( x )

        N    =   length(x)
        dx   =   (x[end]-x[1])/(N-1)
        xmin = x[1]
        xmax = x[end]+dx

        Mesh( xmin, xmax, N)

    end

    function Mesh(param :: NamedTuple)

        if :xmin in keys(param) && :xmax in keys(param)
            xmin = param.xmin
            xmax = param.xmax
        else
            xmin = - param.L
            xmax =   param.L
        end
        N    =   param.N

        Mesh( xmin, xmax, N)

    end
end


text:
"""
Perdew, Burke, Ernzerhof (DOI: 10.1103/PhysRevLett.77.3865)
"""
function energy_per_particle(::Val{:gga_x_pbe}, ρ, σ)
    β = 0.06672455060314922
    # μ eqn: below (12) in PBE paper
    energy_per_particle_x_pbe(ρ, σ; κ=0.8040, μ=β/3 * π^2)
end

"""
Zhang, Yang 1998 (DOI 10.1103/physrevlett.80.890)
"""
function energy_per_particle(::Val{:gga_x_pbe_r}, ρ, σ)
    β = 0.06672455060314922  # like in PBE
    energy_per_particle_x_pbe(ρ, σ; κ=1.245, μ=β/3 * π^2)
end

"""
Xu, Goddard 2004 (DOI 10.1063/1.1771632)
"""
function energy_per_particle(::Val{:gga_x_xpbe}, ρ, σ)
    energy_per_particle_x_pbe(ρ, σ; κ=0.91954, μ=0.23214)  # Table 1
end

"""
Perdew, Ruzsinszky, Csonka and others 2008 (DOI 10.1103/physrevlett.100.136406)
"""
function energy_per_particle(::Val{:gga_x_pbe_sol}, ρ, σ)
    energy_per_particle_x_pbe(ρ, σ; κ=0.8040, μ=10/81)  # μ given below equation (2)
end

"""
Constantin, Fabiano, Laricchia 2011 (DOI 10.1103/physrevlett.106.186406)
"""
function energy_per_particle(::Val{:gga_x_apbe}, ρ, σ)
    energy_per_particle_x_pbe(ρ, σ; κ=0.8040, μ=0.260)  # p. 1, right column, bottom
end

"""
del Campo, Gazqez, Trickey and others 2012 (DOI 10.1063/1.3691197)
"""
function energy_per_particle(::Val{:gga_x_pbe_mol}, ρ, σ)
    energy_per_particle_x_pbe(ρ, σ; κ=0.8040, μ=0.27583)  # p. 4, left column, bottom
end

"""
Sarmiento-Perez, Silvana, Marques 2015 (DOI 10.1021/acs.jctc.5b00529)
"""
function energy_per_particle(::Val{:gga_x_pbefe}, ρ, σ)
    energy_per_particle_x_pbe(ρ, σ; κ=0.437, μ=0.346)  # Table 1
end

function energy_per_particle_x_pbe(ρ, σ; κ, μ) # eqns refer to PBE paper
    pbe_x_f(s²) = 1 + κ - κ^2 / (κ + μ * s²)   # (14)
    # rₛ = cbrt(3 / (4π  * ρ))                 # page 2, left column, top
    # kF = cbrt(3π^2 * ρ)                      # page 2, left column, top
    # s  = sqrt(σ) / (2kF * ρ)                 # below (9)
    s² = σ / ( ρ^(4/3) * 2cbrt(3π^2) )^2

    energy_per_particle(Val(:lda_x), ρ) * pbe_x_f(s²)  # (10)
end


text:
doc"""
    sum(x::Var, dim::Int)

Returns the sum over the given dimension.
"""
function Base.sum(x::Var, dim::Int)
    throw("Not implemented yet.")
end


text:
module Services
include("Response.jl")

import Requests
export api_request, api_key, api_secret, api_resource_path, api_savefile_path, api_baseurl

outputtype = Dict(
:CSV  => "csv",
:TSV  => "tsv",
:JSON => "json",
:HTML => "html")

mutable struct ApiConfig
  default_outputtype::String
  resource_path::String
  savefile_path::String
  baseurl::String
  user_agent::String
  key::String
  secret::String
  ApiConfig() = new(EMPTY,EMPTY,EMPTY,EMPTY,EMPTY,EMPTY,EMPTY)
  ApiConfig(key,secret) = new(outputtype[:CSV],defaultresourcepath(),defaultsavefilepath(),defaultbaseurl(),defaultuseragent(),key,secret)
  ApiConfig(resource_path,key,secret) = new(outputtype[:CSV],resource_path,defaultsavepath(),defaultbaseurl(),defaultuseragent(),key,secret)
  ApiConfig(resource_path,savefile_path,key,secret) = new(outputtype[:CSV],resource_path,savefile_path,defaultbaseurl(),defaultuseragent(),key,secret)
end

function defaultresourcepath()
  if haskey(ENV,"CHARM_RESOURCE_PATH")
    return ENV["CHARM_RESOURCE_PATH"]
  else
    root = Sys.is_windows() ? "c:/" : homedir()
    return joinpath(root,"clarusft","data","test")
  end
end

function defaultsavefilepath()
root = Sys.is_windows() ? "c:/" : homedir()
return joinpath(root,"clarusft","data","test")
end

function defaultbaseurl()
  if haskey(ENV,"CHARM_API_URL")
    return ENV["CHARM_API_URL"]
  else
    return "@eval.clarusft.com/api/rest/v1/"
  end
end

function defaultuseragent()
  return "Clarus.jl"*string(Pkg.installed("Clarus"))*" "*"julia-requests"*string(Pkg.installed("Requests"))
end

function keypath()
  root = Sys.is_windows() ? "c:/" : homedir()
  return joinpath(root,"clarusft","keys")
end

const EMPTY           = ""
credentials           = ApiConfig()
credentials_init      = false
const KEYFILE         = "API-Key.txt"
const SECRETFILE      = "API-Secret.txt"
const UTIL_SERVICE    = "Util"

function __init_credentials!__(credentials_init)
  if credentials_init == true
    return
  else
    global credentials = ApiConfig(EMPTY,EMPTY)
    credentials_init=true
    _api_key!(credentials)
    _api_secret!(credentials)
  end
end
function _api_key!(c::ApiConfig)
  if length(c.key) == 0
    if haskey(ENV,"CHARM_API_KEY")
      c.key = ENV["CHARM_API_KEY"]
    else
      c.key = readCredentialFile(KEYFILE)
    end
  end
  return c.key
end

function _api_secret!(c::ApiConfig)
  if length(c.secret) == 0
    if haskey(ENV,"CHARM_API_SECRET")
      c.secret = ENV["CHARM_API_SECRET"]
    else
      c.secret = readCredentialFile(SECRETFILE)
    end
  end
  return c.secret
end

function readCredentialFile(filename::String)
  if length(filename) == 0
    return EMPTY
  end
  filepath = joinpath(keypath(),filename)
  if isfile(filepath)
    f = open(filepath)
    try
      x = strip(readstring(f))
      return x
    finally
      close(f)
    end
  end
  return EMPTY
end

function api_key(x)
  credentials.key = x
end

function api_secret(x)
  credentials.secret = x
end

function api_resource_path(x)
  credentials.resource_path = x
end

function api_savefile_path(x)
  credentials.savefile_path = x
end

function api_baseurl(x)
  credentials.baseurl = x
end

function requesterrormessage(r)
  errormessage   = get(r.headers,MESSAGES,"")  #If Kwargs are blank,
  errormessage   = string(errormessage,"\n",String(r.data)) #If function name is wrong.
  return errormessage
end

function url(category,functionName,output=credentials.default_outputtype)
  # urlBase = "https://" * _api_key!(credentials) * ":" * _api_secret!(credentials) * credentials.baseurl
  urlBase = "https://" * credentials.baseurl
  restUrl  =  urlBase * category * "/" * functionName * "."*output
  return restUrl
end

function api_request(category, functionName; params...)
  __init_credentials!__(credentials_init)
  restUrl  = url(category,functionName)
  headers = Dict("Authorization"=>"Basic $(base64encode(_api_key!(credentials)*":"*_api_secret!(credentials)))")
  r = Requests.post(restUrl,headers=headers, json=Dict(params))
  if Requests.statuscode(r)!=200
    error("Request to " * category * "/" * functionName * " failed with status code: " * string(Requests.statuscode(r))* "\n"*requesterrormessage(r))
  end
  return Response(r)
end

end #MODULE-END


text:
# This file is a part of EncodedArrays.jl, licensed under the MIT License (MIT).

__precompile__(true)

module EncodedArrays

using ArraysOfArrays
using BitOperations
using FillArrays
using StructArrays

include("encoded_array.jl")
include("varlen_io.jl")
include("varlen_diff_codec.jl")

end # module


text:
## from Juno

# Super-simple pandoc interface.
function pandoc(infn, infmt::AbstractString, outfmt::AbstractString, args::AbstractString...)
    cmd = String["pandoc",
                     "--from=$(infmt)",
                     "--to=$(outfmt)"]
    for arg in args
        push!(cmd, arg)
    end

    read(pipeline(infn, Cmd(cmd)), String)
end

"""
Take a string in markdown and covert to LaTeX via pandoc
"""
function markdown_to_latex(txt)
    infn = tempname() * ".md"
    io = open(infn, "w"); print(io, txt); close(io)
    out = pandoc(infn, "markdown", "latex")
    rm(infn)
    out
end


text:
using ClangCompiler
using Test

@testset "Types" begin
    src = joinpath(@__DIR__, "code", "main.cpp")
    args = get_compiler_args()
    haskey(ENV, "CI") && push!(args, "-v")
    irgen = IncrementalIRGenerator(src, args)
    ctx = ClangCompiler.getASTContext(get_instance(irgen))

    @test clty_to_jlty(jlty_to_clty(Cvoid, ctx)) == Cvoid
    @test clty_to_jlty(jlty_to_clty(Bool, ctx)) == Bool
    @test clty_to_jlty(jlty_to_clty(Int8, ctx)) == Int8
    @test clty_to_jlty(jlty_to_clty(Int16, ctx)) == Int16
    @test clty_to_jlty(jlty_to_clty(Int32, ctx)) == Int32
    @test clty_to_jlty(jlty_to_clty(Int64, ctx)) == Int64
    @test clty_to_jlty(jlty_to_clty(Int128, ctx)) == Int128
    @test clty_to_jlty(jlty_to_clty(UInt8, ctx)) == UInt8
    @test clty_to_jlty(jlty_to_clty(UInt16, ctx)) == UInt16
    @test clty_to_jlty(jlty_to_clty(UInt32, ctx)) == UInt32
    @test clty_to_jlty(jlty_to_clty(UInt64, ctx)) == UInt64
    @test clty_to_jlty(jlty_to_clty(UInt128, ctx)) == UInt128
    @test clty_to_jlty(jlty_to_clty(Float16, ctx)) == Float16
    @test clty_to_jlty(jlty_to_clty(Float32, ctx)) == Float32
    @test clty_to_jlty(jlty_to_clty(Float64, ctx)) == Float64
    @test clty_to_jlty(jlty_to_clty(Ptr{Cvoid}, ctx)) == Ptr{Cvoid}

    dispose(irgen)
end


text:
using BinaryBuilder

# Collection of sources required to build MbedTLS
sources = [
    "https://github.com/ARMmbed/mbedtls.git" =>
    "d1236a790ea0c685e136e9a72b4666882119d3f6",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd mbedtls/
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain -DUSE_SHARED_MBEDTLS_LIBRARY=On
make -j${nproc} && make install
if [ $target == "x86_64-w64-mingw32" ]; then
    cp $prefix/lib/*.dll $prefix/bin/.
elif [ $target == "i686-w64-mingw32" ]; then
    cp $prefix/lib/*.dll $prefix/bin/.
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, :glibc),
    Linux(:x86_64, :glibc),
    Linux(:aarch64, :glibc),
    Linux(:armv7l, :glibc),
    Linux(:powerpc64le, :glibc),
    Windows(:x86_64),
    Windows(:i686),
    MacOS()
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libmbedcrypto", :libmbedcrypto),
    LibraryProduct(prefix, "libmbedtls", :libmbedtls),
    LibraryProduct(prefix, "libmbedx509", :libmbedx509)
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "MbedTLS", sources, script, platforms, products, dependencies)

max_stars_repo_name: staticfloat/MbedTLSBuilder

max_issues_repo_name: staticfloat/MbedTLSBuilder

max_forks_repo_name: staticfloat/MbedTLSBuilder


text:
include("../src/ekp.jl") # for ekpobj
include("../src/transforms.jl") # for transform_prior_to_real, transform_real_to_prior
include("../src/units.jl") # for apply_units_transform

using Plots
using Plots.PlotMeasures
using LinearAlgebra
using JLD2
using Statistics

using .EKP

"""
Plotting script for the ensemble on a 2D parameter space, It plots:
- iterations of EKP in parameter space
- the location of the true parameters in parameter space
- first iteration in dark grey, intermediate iterations in light grey, final iteration in pink

3 panels        'transformed' (transformed parameters coordinates)
              'untransformed' (real parameter coordinates)
         'zoom_untransformed' (zoom-in around true parameters of 'untransformed')

requires:
--------
- EKP object
- min and max iteration of EKP
- true parameter coordinates (atm just given as numbers should be stored somewhere)
- transformations and inverse transformations real->prior space (if required)
- units transform to add units (if required)
- limits for x/y coordinates etc. and other plot preferences

"""

######################################################################
function main()

    #object
    @load "ekp.jld2" ekpobj

    #number of EKP iterations
    min_ekp_it = 1
    max_ekp_it = 10

    #true parameters
    true_params_raw = [0,0]

    #obtain parameters
    u = ekpobj.u[min_ekp_it:max_ekp_it]#it x [enssize x param]
    u_flat = cat(u...,dims=1)#[(itxens) x param]
    ens_size = ekpobj.J

    # Transformations
    u_raw = transform_prior_to_real(u_flat)
    true_params = transform_real_to_prior(true_params)

    # Specific units
    u_raw = apply_units_transform(u_raw)
    true_params_raw = apply_units_transform(true_params_raw)

    #plot - transformed
    gr(size=(500,500))

    circ=Shape(Plots.partialcircle(0, 2π))

    plot(u_flat[(min_ekp_it - 1) * ens_size + 1 : min_ekp_it * ens_size, 1],
         u_flat[(min_ekp_it - 1) * ens_size + 1 : min_ekp_it * ens_size, 2],
         seriestype=:scatter,
         markershape=circ,
         markercolor=:grey,
         markersize=5,
         msw=0,
         legend=false,
         grid=false,
         xlims=[-3.0,3.0],
         ylims=[7.5, 13],
         annotations=(-2.725,12.725,"(b)"),
         dpi=300,
         framestyle=:box,
         left_margin=50px,
         top_margin=25px,
         bottom_margin=50px)

    for i = min_ekp_it+1 : max_ekp_it
        if i<max_ekp_it
            plot!(u_flat[(i - 1) * ens_size + 1 : i * ens_size , 1],
                  u_flat[(i - 1) * ens_size + 1 : i * ens_size , 2],
                  seriestype=:scatter,
                  markershape=circ,
                  markercolor=:grey,
                  markeralpha=0.3*(max_ekp_it-i)/max_ekp_it,
                  msw=0,
                  markersize=5 )
       else
            plot!(u_flat[(i - 1) * ens_size + 1 : i * ens_size , 1],
                  u_flat[(i - 1) * ens_size + 1 : i * ens_size , 2],
                  seriestype=:scatter,
                  markershape=circ,
                  markercolor=:lightpink,
                  msw=0,
                  markersize=5)
       end

    end

    vline!([true_params[1]],color=:blue,linealpha=0.5,linestyle=:dash)
    hline!([true_params[2]],color=:blue,linealpha=0.5,linestyle=:dash)

    xlabel!("logit-relative humidity parameter")
    ylabel!("log-timescale parameter (seconds)")
    savefig(outdir*"ensemble_transformed.pdf")

    #plot - untransformed
    plot(u_raw[(min_ekp_it - 1) * ens_size + 1 : min_ekp_it * ens_size, 1],
         u_raw[(min_ekp_it - 1) * ens_size + 1 : min_ekp_it * ens_size, 2],
         seriestype=:scatter,
         markershape=circ,
         markercolor=:grey,
         msw=0,
         markersize=5,
         legend=false,
         grid=false,
         annotations=(0.05,62,"(a)"),
         xlims=[0,1],
         ylims=[0,65],
         dpi=300,
         framestyle=:box,
         left_margin=50px,
         top_margin=25px,
         bottom_margin=50px)

    for i = min_ekp_it+1 : max_ekp_it
        if i<max_ekp_it
            plot!(u_raw[(i - 1) * ens_size + 1 : i * ens_size , 1],
                  u_raw[(i - 1) * ens_size + 1 : i * ens_size , 2],
                  seriestype=:scatter,
                  markershape=circ,
                  markercolor=:grey,
                  markeralpha=0.3*(max_ekp_it-i)/max_ekp_it,
                  msw=0,
                  markersize=5 )
       else
            plot!(u_raw[(i - 1) * ens_size + 1 : i * ens_size , 1],
                  u_raw[(i - 1) * ens_size + 1 : i * ens_size , 2],
                  seriestype=:scatter,
                  markershape=circ,
                  markercolor=:lightpink,
                  msw=0,
                  markersize=5)
       end

    end
    vline!([true_params_raw[1]],color=:blue,linealpha=0.5,linestyle=:dash)
    hline!([true_params_raw[2]],color=:blue,linealpha=0.5,linestyle=:dash)

    xlabel!("relative humidity parameter")
    ylabel!("timescale parameter (hours)")
    savefig(outdir*"ensemble_untransformed.pdf")

    #zoomed in - removes points outside of zoomed in domain area
    zoom_rhum_bd = [0.6,0.8]
    zoom_tau_bd = [1.0,4.0]
    zoom_rhum_in_bd = (u_raw[:,1] .> zoom_rhum_bd[1]) .& (u_raw[:,1] .< zoom_rhum_bd[2])
    zoom_tau_in_bd = (u_raw[:,2] .> zoom_tau_bd[1]) .& (u_raw[:,2] .< zoom_tau_bd[2])
    onetoN = 1:ens_size*(max_ekp_it - min_ekp_it+1)
    idx_in_bd = onetoN[zoom_rhum_in_bd .& zoom_tau_in_bd]

    u_raw_in_bd = u_raw[idx_in_bd,:]

    it_in_bd = (idx_in_bd .>= (min_ekp_it - 1) * ens_size + 1 ) .& (idx_in_bd .<= min_ekp_it * ens_size)
    plot(u_raw_in_bd[it_in_bd, 1],
         u_raw_in_bd[it_in_bd, 2],
         seriestype=:scatter,
         markershape=circ,
         markercolor=:grey,
         msw=0,
         markersize=5,
         legend=false,
         grid=false,
         annotations=(0.05,62,"(a)"),
         xlims=zoom_rhum_bd,
         ylims=zoom_tau_bd,
         dpi=300,
         framestyle=:box,
         left_margin=50px,
         top_margin=25px,
         bottom_margin=50px)

    for i = min_ekp_it+1 : max_ekp_it
        it_in_bd = (idx_in_bd .>= (i - 1) * ens_size + 1 ) .& (idx_in_bd .<= i * ens_size)
        if i<max_ekp_it
            plot!(u_raw_in_bd[it_in_bd, 1],
                  u_raw_in_bd[it_in_bd, 2],
                  seriestype=:scatter,
                  markershape=circ,
                  markercolor=:grey,
                  markeralpha=0.3*(max_ekp_it-i)/max_ekp_it,
                  msw=0,
                  markersize=5 )
       else
            plot!(u_raw_in_bd[it_in_bd, 1],
                  u_raw_in_bd[it_in_bd, 2],
                  seriestype=:scatter,
                  markershape=circ,
                  markercolor=:lightpink,
                  msw=0,
                  markersize=5)
       end

    end
    vline!([true_params_raw[1]],color=:blue,linealpha=0.5,linestyle=:dash)
    hline!([true_params_raw[2]],color=:blue,linealpha=0.5,linestyle=:dash)

    xlabel!("relative humidity parameter")
    ylabel!("timescale parameter (hours)")
    savefig(outdir*"ensemble_zoom_untransformed.pdf")
end

main()


text:
using NiLang
using Test

@testset "identity" begin
    x, y = 0.2, 0.5
    @instr x += identity(y)
    @test x == 0.7 && y==0.5
end

@testset "*, /" begin
    x, y, out = 2.0, 2.0, 1.0
    @instr out += x * y
    @test x == 2.0 && y == 2.0 && out == 5.0
    x, y, out = 2.0, 2.0, 1.0
    @instr out += x / y
    @test x == 2.0 && y == 2.0 && out == 2.0

    out = Fixed43(0.0)
    x = 1
    @instr out += x/2
    @test out === Fixed43(0.5)
    @instr out -= x/2
    @test out === Fixed43(0.0)
end

@testset "SWAP" begin
    x, y = 1, 2
    @instr SWAP(x, y)
    @test x == 2 && y == 1
end

@testset "NEG" begin
    x = 0.3
    @instr NEG(x)
    @test x == -0.3
    @test check_inv(NEG, (x,))
end

@testset "INV" begin
    x = 0.2
    @instr INV(x)
    @test x == 5.0
    @test check_inv(INV, (x,))
end

@testset "AddConst" begin
    x = 0.3
    @instr AddConst(4.0)(x)
    @test x == 4.3
    @test check_inv(AddConst(4.0), (x,))

    x = 0.3
    @instr SubConst(4.0)(x)
    @test x == -3.7
    @test check_inv(SubConst(4.0), (x,))
end

@testset "FLIP" begin
    x = false
    @instr FLIP(x)
    @test x == true
    @test check_inv(FLIP, (x,))
end

@testset "ROT" begin
    x, y, θ = 0.0, 1.0, π
    @test check_inv(ROT, (x, y, θ); verbose=true)
    @test check_inv(IROT, (x, y, θ); verbose=true)
end

@testset "INC, DEC" begin
    x = Int32(2)
    @instr INC(x)
    @test x === Int32(3)
    @instr DEC(x)
    @test x === Int32(2)
end

@testset "HADAMARD" begin
    x = 0.5
    y = 0.8
    @test check_inv(HADAMARD, (x, y))
end

@testset "dataviews" begin
    @i function f(z, y, x)
        y += cos(x |> INV)
        z += tan(y |> AddConst(4.0))
        z += y * (x |> NEG |> SubConst(0.5) |> INV)
        z += sin(x |> INV)
    end
    @test check_inv(f, (0.2, 0.5, 0.8))
end

@testset "fixed point arithmetics" begin
    for op in [exp, log, sin, sinh, asin, cos, cosh, acos, tan, tanh, atan]
        x, y = Fixed43(2.0), Fixed43(0.5)
        @instr x += op(y)
        @test x ≈ 2.0 + op(0.5)
    end
    for op in [SWAP, HADAMARD]
        x, y = Fixed43(2.0), Fixed43(0.5)
        @instr op(x, y)
        @test x ≈ op(2.0, 0.5)[1]
        @test y ≈ op(2.0, 0.5)[2]
    end
    for op in [NEG, INC, DEC]
        x = Fixed43(2.0)
        @instr op(x)
        @test x ≈ op(2.0)
    end
    for op in [^, /]
        x, y, z = Fixed43(2.0), Fixed43(0.5), Fixed43(0.8)
        @instr PlusEq(op)(x, y, z)
        @test x ≈ 2.0 + op(0.5, 0.8)
    end
end

@testset "additive identity" begin
    struct TestAdd{T}
        x::T
        y::Vector{T}
    end
    @test getfield.(PlusEq(identity)(TestAdd(1, [2]), TestAdd(10, [2])), :x) == (TestAdd(11, [4]).x, TestAdd(10, [2]).x)
    @test getfield.(PlusEq(identity)(TestAdd(1, [2]), TestAdd(10, [2])), :y) == (TestAdd(11, [4]).y, TestAdd(10, [2]).y)
end