text:
module Quadrature

using Requires, Reexport,  MonteCarloIntegration, QuadGK, HCubature
@reexport using DiffEqBase
using Zygote, ReverseDiff, ForwardDiff , LinearAlgebra

import ChainRulesCore
import ChainRulesCore: NoTangent
import ZygoteRules

struct QuadGKJL <: DiffEqBase.AbstractQuadratureAlgorithm end
struct HCubatureJL <: DiffEqBase.AbstractQuadratureAlgorithm end
struct VEGAS <: DiffEqBase.AbstractQuadratureAlgorithm
    nbins::Int
    ncalls::Int
end

VEGAS(;nbins = 100,ncalls = 1000) = VEGAS(nbins,ncalls)

abstract type AbstractCubaAlgorithm <: DiffEqBase.AbstractQuadratureAlgorithm end
struct CubaVegas <:AbstractCubaAlgorithm end
struct CubaSUAVE <: AbstractCubaAlgorithm end
struct CubaDivonne <: AbstractCubaAlgorithm end
struct CubaCuhre <: AbstractCubaAlgorithm end

abstract type AbstractCubatureJLAlgorithm <: DiffEqBase.AbstractQuadratureAlgorithm end
struct CubatureJLh <: AbstractCubatureJLAlgorithm end
struct CubatureJLp <: AbstractCubatureJLAlgorithm end

abstract type QuadSensitivityAlg end
struct ReCallVJP{V}
    vjp::V
end

abstract type QuadratureVJP end
struct ZygoteVJP end
struct ReverseDiffVJP
    compile::Bool
end

function scale_x!(_x,ub,lb,x)
    _x .= (ub .- lb) .* x .+ lb
    _x
end

function scale_x(ub,lb,x)
    (ub .- lb) .* x .+ lb
end

function v_inf(t)
	return t ./ (1 .- t.^2)
end

function v_semiinf(t , a , upto_inf)
	if upto_inf == true
		return a .+ (t ./ (1 .- t))
	else
		return a .+ (t ./ (1 .+ t))
	end
end
function transfor_inf_number(t , p , f , lb , ub)
	if lb == -Inf && ub == Inf
		j = (1 .+ t.^2 )/(1 .- t.^2).^2
		return f(v_inf(t) , p)*(j)
	elseif lb != -Inf && ub == Inf
		a = lb
		j = 1 ./ ((1 .- t).^2)
		return f(v_semiinf(t , a , 1) , p)*(j)
	elseif lb == -Inf && ub != Inf
		a = ub
		j = 1 ./ ((1 .+ t ).^2)
		return f(v_semiinf(t , a , 0) , p)*(j)
	end
end

function transform_inf(t , p , f , lb , ub)
	if lb isa Number && ub isa Number
		return transfor_inf_number(t , p , f , lb , ub)
	end

	lbb = lb .== -Inf
	ubb = ub .== Inf
	_none = .!lbb .& .!ubb
	_inf = lbb .& ubb
	semiup = .!lbb .& ubb
	semilw = lbb  .& .!ubb

	function v(t)
		return t.*_none + v_inf(t).*_inf + v_semiinf(t , lb , 1).*semiup + v_semiinf(t , ub , 0).*semilw
	end
    jac = Zygote.@ignore ForwardDiff.jacobian(x -> v(x), t)
	j = det(jac)
	f(v(t) , p)*(j)
end

function transformation_if_inf(prob)
	if (prob.lb isa Number && prob.ub isa Number && (prob.ub == Inf || prob.lb == -Inf))
		if prob.lb == -Inf && prob.ub == Inf
			g = prob.f
			h(t , p) = transform_inf(t , p , g , prob.lb , prob.ub)
			lb = -1.00
			ub = 1.00
	        _prob = QuadratureProblem(h , lb ,ub, prob.p , nout = prob.nout , batch = prob.batch , prob.kwargs...)
			return _prob
		elseif prob.lb != -Inf && prob.ub == Inf
			g = prob.f
			h_(t , p) = transform_inf(t , p , g , prob.lb , prob.ub)
			lb = 0.00
			ub = 1.00
	        _prob = QuadratureProblem(h_ , lb ,ub, prob.p , nout = prob.nout , batch = prob.batch , prob.kwargs...)
			return _prob
		elseif prob.lb == -Inf && prob.ub != Inf
			g = prob.f
			_h(t , p) = transform_inf(t , p , g , prob.lb , prob.ub)
			lb = -1.00
			ub = 0.00
	        _prob = QuadratureProblem(_h , lb ,ub, prob.p , nout = prob.nout , batch = prob.batch , prob.kwargs...)
			return _prob
		end
	end
	if -Inf in prob.lb || Inf in prob.ub
		lbb = prob.lb .== -Inf
		ubb = prob.ub .== Inf
		_none = .!lbb .& .!ubb
		_inf = lbb .& ubb
		_semiup = .!lbb .& ubb
		_semilw = lbb  .& .!ubb

		_lb = 0.00.*_semiup + -1.00.*_inf + -1.00.*_semilw +  _none.*prob.lb
		_ub = 1.00.*_semiup + 1.00.*_inf  + 0.00.*_semilw  + _none.*prob.ub
		g = prob.f
		hs(t , p) = transform_inf(t , p , g , prob.lb , prob.ub)
		_prob = QuadratureProblem(hs, _lb ,_ub, prob.p , nout = prob.nout , batch = prob.batch , prob.kwargs...)
		return _prob
	end
	return prob
end
function DiffEqBase.solve(prob::QuadratureProblem,::Nothing,sensealg,lb,ub,p,args...;
                          reltol = 1e-8, abstol = 1e-8, kwargs...)
    if lb isa Number
        __solve(prob,QuadGKJL();reltol=reltol,abstol=abstol,kwargs...)
    elseif length(lb) > 8 && reltol < 1e-4 || abstol < 1e-4
        __solve(prob,VEGAS();reltol=reltol,abstol=abstol,kwargs...)
    else
        __solve(prob,HCubatureJL();reltol=reltol,abstol=abstol,kwargs...)
    end
end

function DiffEqBase.solve(prob::QuadratureProblem,
                            alg::DiffEqBase.AbstractQuadratureAlgorithm,
                            args...; sensealg = ReCallVJP(ZygoteVJP()), kwargs...)
   prob = transformation_if_inf(prob)
  __solvebp(prob,alg,sensealg,prob.lb,prob.ub,prob.p,args...;kwargs...)
end

# Give a layer to intercept with AD
__solvebp(args...;kwargs...) = __solvebp_call(args...;kwargs...)

function __solvebp_call(prob::QuadratureProblem,::QuadGKJL,sensealg,lb,ub,p,args...;
                          reltol = 1e-8, abstol = 1e-8,
                          maxiters = typemax(Int),
                          kwargs...)
    if isinplace(prob) || lb isa AbstractArray || ub isa AbstractArray
        error("QuadGKJL only accepts one-dimensional quadrature problems.")
    end
    @assert prob.batch == 0
    @assert prob.nout == 1
    p = p
    f = x -> prob.f(x,p)
    val,err = quadgk(f, lb, ub,
                     rtol=reltol, atol=abstol,
                     kwargs...)
    DiffEqBase.build_solution(prob,QuadGKJL(),val,err,retcode = :Success)
end

function __solvebp_call(prob::QuadratureProblem,::HCubatureJL,sensealg,lb,ub,p,args...;
                          reltol = 1e-8, abstol = 1e-8,
                          maxiters = typemax(Int),
                          kwargs...)
    p = p

    if isinplace(prob)
        dx = zeros(prob.nout)
        f = (x) -> (prob.f(dx,x,p); dx)
    else
        f = (x) -> prob.f(x,p)
    end
    @assert prob.batch == 0

    if lb isa Number
        val,err = hquadrature(f, lb, ub;
                              rtol=reltol, atol=abstol,
                              maxevals=maxiters, kwargs...)
    else
        val,err = hcubature(f, lb, ub;
                            rtol=reltol, atol=abstol,
                            maxevals=maxiters, kwargs...)
    end
    DiffEqBase.build_solution(prob,HCubatureJL(),val,err,retcode = :Success)
end

function __solvebp_call(prob::QuadratureProblem,alg::VEGAS,sensealg,lb,ub,p,args...;
                          reltol = 1e-8, abstol = 1e-8,
                          maxiters = typemax(Int),
                          kwargs...)
    p = p
    @assert prob.nout == 1
    if prob.batch == 0
        if isinplace(prob)
          dx = zeros(prob.nout)
          f = (x) -> (prob.f(dx,x,p); dx)
        else
          f = (x) -> prob.f(x,p)
        end
    else
        if isinplace(prob)
          dx = zeros(prob.batch)
          f = (x) -> (prob.f(dx,x',p); dx)
        else
          f = (x) -> prob.f(x',p)
        end
    end
    val,err,chi = vegas(f, lb, ub, rtol=reltol, atol=abstol,
                        maxiter = maxiters, nbins = alg.nbins,
                        ncalls = alg.ncalls, batch=prob.batch != 0, kwargs...)
    DiffEqBase.build_solution(prob,alg,val,err,chi=chi,retcode = :Success)
end

function __init__()
    @require Cubature="667455a9-e2ce-5579-9412-b964f529a492" begin
        function __solvebp_call(prob::QuadratureProblem,
                                  alg::AbstractCubatureJLAlgorithm,
                                  sensealg, lb, ub, p, args...;
                                  reltol = 1e-8, abstol = 1e-8,
                                  maxiters = typemax(Int),
                                  kwargs...)
            prob = transformation_if_inf(prob) #intercept for infinite transformation
            nout = prob.nout
            if nout == 1
                if prob.batch == 0
                    if isinplace(prob)
                        dx = zeros(prob.nout)
                        f = (x) -> (prob.f(dx,x,p); dx[1])
                    else
                        f = (x) -> prob.f(x,p)[1]
                    end
                    if lb isa Number
                        if alg isa CubatureJLh
                            _val,err = Cubature.hquadrature(f, lb, ub;
                                                           reltol=reltol, abstol=abstol,
                                                           maxevals=maxiters)
                        else
                            _val,err = Cubature.pquadrature(f, lb, ub;
                                                           reltol=reltol, abstol=abstol,
                                                           maxevals=maxiters)
                        end
                        val = prob.f(lb,p) isa Number ? _val : [_val]
                    else
                        if alg isa CubatureJLh
                            _val,err = Cubature.hcubature(f, lb, ub;
                                                     reltol=reltol, abstol=abstol,
                                                     maxevals=maxiters)
                        else
                            _val,err = Cubature.pcubature(f, lb, ub;
                                                         reltol=reltol, abstol=abstol,
                                                         maxevals=maxiters)
                        end

                        if isinplace(prob) || !isa(prob.f(lb,p), Number)
                            val = [_val]
                        else
                            val = _val
                        end
                     end
                else
                    if isinplace(prob)
                        f = (x,dx) -> prob.f(dx',x,p)
                    elseif lb isa Number
                        if prob.f([lb ub], p) isa Vector
                            f = (x,dx) -> (dx .= prob.f(x',p))
                        else
                            f = function (x,dx)
                                dx[:] = prob.f(x',p)
                            end
                        end
                    else
                        if prob.f([lb ub], p) isa Vector
                            f = (x,dx) -> (dx .= prob.f(x,p))
                        else
                            f = function (x,dx)
                                dx .= prob.f(x,p)[:]
                            end
                        end
                    end
                    if lb isa Number
                        if alg isa CubatureJLh
                            _val,err = Cubature.hquadrature_v(f, lb, ub;
                                                             reltol=reltol, abstol=abstol,
                                                             maxevals=maxiters)
                        else
                            _val,err = Cubature.pquadrature_v(f, lb, ub;
                                                             reltol=reltol, abstol=abstol,
                                                             maxevals=maxiters)
                        end
                    else
                        if alg isa CubatureJLh
                            _val,err = Cubature.hcubature_v(f, lb, ub;
                                                           reltol=reltol, abstol=abstol,
                                                           maxevals=maxiters)
                        else
                            _val,err = Cubature.pcubature_v(f, lb, ub;
                                                           reltol=reltol, abstol=abstol,
                                                           maxevals=maxiters)
                        end
                     end
                     val = _val isa Number ? [_val] : _val
                 end
             else
                 if prob.batch == 0
                     if isinplace(prob)
                         f = (x,dx) -> (prob.f(dx,x,p); dx)
                     else
                         f = (x,dx) -> (dx .= prob.f(x,p))
                     end
                     if lb isa Number
                         if alg isa CubatureJLh
                             val,err = Cubature.hquadrature(nout, f, lb, ub;
                                                            reltol=reltol, abstol=abstol,
                                                            maxevals=maxiters)
                         else
                             val,err = Cubature.pquadrature(nout, f, lb, ub;
                                                            reltol=reltol, abstol=abstol,
                                                            maxevals=maxiters)
                         end
                     else
                         if alg isa CubatureJLh
                             val,err = Cubature.hcubature(nout, f, lb, ub;
                                                          reltol=reltol, abstol=abstol,
                                                          maxevals=maxiters)
                         else
                             val,err = Cubature.pcubature(nout, f, lb, ub;
                                                          reltol=reltol, abstol=abstol,
                                                          maxevals=maxiters)
                         end
                      end
                 else
                     if isinplace(prob)
                         f = (x,dx) -> prob.f(dx,x,p)
                     else
                         if lb isa Number
                             f = (x,dx) -> (dx .= prob.f(x',p))
                         else
                             f = (x,dx) -> (dx .= prob.f(x,p))
                         end
                     end

                     if lb isa Number
                         if alg isa CubatureJLh
                             val,err = Cubature.hquadrature_v(nout, f, lb, ub;
                                                              reltol=reltol, abstol=abstol,
                                                              maxevals=maxiters)
                         else
                             val,err = Cubature.pquadrature_v(nout, f, lb, ub;
                                                              reltol=reltol, abstol=abstol,
                                                              maxevals=maxiters)
                         end
                     else
                         if alg isa CubatureJLh
                             val,err = Cubature.hcubature_v(nout, f, lb, ub;
                                                            reltol=reltol, abstol=abstol,
                                                            maxevals=maxiters)
                         else
                             val,err = Cubature.pcubature_v(nout, f, lb, ub;
                                                            reltol=reltol, abstol=abstol,
                                                            maxevals=maxiters)
                         end
                      end
                  end
              end
              DiffEqBase.build_solution(prob,alg,val,err,retcode = :Success)
        end
    end

    @require Cuba="8a292aeb-7a57-582c-b821-06e4c11590b1" begin
        function __solvebp_call(prob::QuadratureProblem,alg::AbstractCubaAlgorithm,sensealg,
                                  lb,ub,p,args...;
                                  reltol = 1e-8, abstol = 1e-8,
                                  maxiters = alg isa CubaSUAVE ? 1000000 : typemax(Int),
                                  kwargs...)
          prob = transformation_if_inf(prob) #intercept for infinite transformation
          p = p
          if lb isa Number && prob.batch == 0
              _x = Float64[lb]
          elseif lb isa Number
              _x = zeros(length(lb),prob.batch)
          elseif prob.batch == 0
              _x = zeros(length(lb))
          else
              _x = zeros(length(lb),prob.batch)
          end
          ub = ub
          lb = lb

          if prob.batch == 0
              if isinplace(prob)
                  f = function (x,dx)
                      prob.f(dx,scale_x!(_x,ub,lb,x),p)
                      dx .*= prod((y)->y[1]-y[2],zip(ub,lb))
                  end
              else
                  f = function (x,dx)
                      dx .= prob.f(scale_x!(_x,ub,lb,x),p) .* prod((y)->y[1]-y[2],zip(ub,lb))
                  end
              end
          else
              if lb isa Number
                  if isinplace(prob)
                      f = function (x,dx)
                          #todo check scale_x!
                          prob.f(dx',scale_x!(view(_x,1:length(x)),ub,lb,x),p)
                          dx .*= prod((y)->y[1]-y[2],zip(ub,lb))
                      end
                  else
                      if prob.f([lb ub], p) isa Vector
                          f = function (x,dx)
                              dx .= prob.f(scale_x(ub,lb,x),p)' .* prod((y)->y[1]-y[2],zip(ub,lb))
                          end
                      else
                          f = function (x,dx)
                              dx .= prob.f(scale_x(ub,lb,x),p) .* prod((y)->y[1]-y[2],zip(ub,lb))
                          end
                      end
                  end
              else
                  if isinplace(prob)
                      f = function (x,dx)
                          prob.f(dx,scale_x(ub,lb,x),p)
                          dx .*= prod((y)->y[1]-y[2],zip(ub,lb))
                      end
                  else
                      if prob.f([lb ub], p) isa Vector
                          f = function (x,dx)
                              dx .= prob.f(scale_x(ub,lb,x),p)' .* prod((y)->y[1]-y[2],zip(ub,lb))
                          end
                      else
                          f = function (x,dx)
                              dx .= prob.f(scale_x(ub,lb,x),p) .* prod((y)->y[1]-y[2],zip(ub,lb))
                          end
                      end
                  end
              end
          end

          ndim = length(lb)

          nvec = prob.batch == 0 ? 1 : prob.batch

          if alg isa CubaVegas
              out = Cuba.vegas(f, ndim, prob.nout; rtol = reltol,
                               atol = abstol, nvec = nvec,
                               maxevals = maxiters, kwargs...)
          elseif alg isa CubaSUAVE
              out = Cuba.suave(f, ndim, prob.nout; rtol = reltol,
                               atol = abstol, nvec = nvec,
                               maxevals = maxiters, kwargs...)
          elseif alg isa CubaDivonne
              out = Cuba.divonne(f, ndim, prob.nout; rtol = reltol,
                                 atol = abstol, nvec = nvec,
                                 maxevals = maxiters, kwargs...)
          elseif alg isa CubaCuhre
              out = Cuba.cuhre(f, ndim, prob.nout; rtol = reltol,
                               atol = abstol, nvec = nvec,
                               maxevals = maxiters, kwargs...)
          end

          if isinplace(prob) || prob.batch != 0
              val = out.integral
          else
              if prob.nout == 1 && prob.f(lb, p) isa Number
                  val = out.integral[1]
              else
                  val = out.integral
              end
          end

          DiffEqBase.build_solution(prob,alg,val,out.error,
                         chi=out.probability,retcode = :Success)
        end
    end
end

function ChainRulesCore.rrule(::typeof(__solvebp),prob,alg,sensealg,lb,ub,p,args...;kwargs...)
    out = __solvebp_call(prob,alg,sensealg,lb,ub,p,args...;kwargs...)
    function quadrature_adjoint(Δ)
        y = typeof(Δ) <: Array{<:Number,0} ? Δ[1] : Δ
        if isinplace(prob)
            dx = zeros(prob.nout)
            _f = (x) -> prob.f(dx,x,p)
            if sensealg.vjp isa ZygoteVJP
                dfdp = function (dx,x,p)
                    _,back = Zygote.pullback(p) do p
                        _dx = Zygote.Buffer(x, prob.nout, size(x,2))
                        prob.f(_dx,x,p)
                        copy(_dx)
                    end

                    z = zeros(size(x,2))
                    for idx in 1:size(x,2)
                        z[1] = 1
                        dx[:,idx] = back(z)[1]
                        z[idx]=0
                    end
                end
            elseif sensealg.vjp isa ReverseDiffVJP
                error("TODO")
            end
        else
            _f = (x) -> prob.f(x,p)
            if sensealg.vjp isa ZygoteVJP
                if prob.batch > 0
                    dfdp = function (x,p)
                        _,back = Zygote.pullback(p->prob.f(x,p),p)

                        out = zeros(length(p),size(x,2))
                        z = zeros(size(x,2))
                        for idx in 1:size(x,2)
                            z[idx] = 1
                            out[:,idx] = back(z)[1]
                            z[idx]=0
                        end
                        out
                    end
                else
                    dfdp = function (x,p)
                        _,back = Zygote.pullback(p->prob.f(x,p),p)
                        back(y)[1]
                    end
                end

            elseif sensealg.vjp isa ReverseDiffVJP
                error("TODO")
            end
        end

        dp_prob = QuadratureProblem(dfdp,lb,ub,p;nout=length(p),batch = prob.batch,kwargs...)

        if p isa Number
            dp = __solvebp_call(dp_prob,alg,sensealg,lb,ub,p,args...;kwargs...)[1]
        else
            dp = __solvebp_call(dp_prob,alg,sensealg,lb,ub,p,args...;kwargs...).u
        end

        if lb isa Number
            dlb = -_f(lb)
            dub = _f(ub)
            return (NoTangent(),NoTangent(),NoTangent(),NoTangent(),dlb,dub,dp,ntuple(x->NoTangent(),length(args))...)
        else
            return (NoTangent(),NoTangent(),NoTangent(),NoTangent(),NoTangent(),NoTangent(),dp,ntuple(x->NoTangent(),length(args))...)
        end
    end
    out,quadrature_adjoint
end

ZygoteRules.@adjoint function ZygoteRules.literal_getproperty(
                        sol::DiffEqBase.QuadratureSolution,::Val{:u})
    sol.u,Δ->(DiffEqBase.build_solution(sol.prob,sol.alg,Δ,sol.resid),)
end

### Forward-Mode AD Intercepts

# Direct AD on solvers with QuadGK and HCubature
function __solvebp(prob,alg::QuadGKJL,sensealg,lb,ub,p::AbstractArray{<:ForwardDiff.Dual{T,V,P},N},args...;kwargs...) where {T,V,P,N}
    __solvebp_call(prob,alg,sensealg,lb,ub,p,args...;kwargs...)
end

function __solvebp(prob,alg::HCubatureJL,sensealg,lb,ub,p::AbstractArray{<:ForwardDiff.Dual{T,V,P},N},args...;kwargs...) where {T,V,P,N}
    __solvebp_call(prob,alg,sensealg,lb,ub,p,args...;kwargs...)
end

# Manually split for the pushforward
function __solvebp(prob,alg,sensealg,lb,ub,p::AbstractArray{<:ForwardDiff.Dual{T,V,P},N},args...;kwargs...) where {T,V,P,N}
    primal = __solvebp_call(prob,alg,sensealg,lb,ub,ForwardDiff.value.(p),args...;kwargs...)

    nout = prob.nout*P

    if isinplace(prob)
        dfdp = function (out,x,p)
            dualp = reinterpret(ForwardDiff.Dual{T,V,P}, p)
            if prob.batch > 0
                dx = similar(dualp, prob.nout, size(x,2))
            else
                dx = similar(dualp, prob.nout)
            end
            prob.f(dx,x,dualp)

            ys = reinterpret(ForwardDiff.Dual{T,V,P}, dx)
            idx=0
            for y in ys
                for p in ForwardDiff.partials(y)
                    out[idx+=1] = p
                end
            end
            return out
        end
    else
        dfdp = function (x,p)
            dualp = reinterpret(ForwardDiff.Dual{T,V,P}, p)
            ys = prob.f(x,dualp)
            if prob.batch > 0
                out = similar(p, V, nout, size(x,2))
            else
                out = similar(p, V, nout)
            end

            idx=0
            for y in ys
                for p in ForwardDiff.partials(y)
                    out[idx+=1] = p
                end
            end

            return out
        end
    end
    rawp = copy(reinterpret(V, p))

    dp_prob = QuadratureProblem(dfdp,lb,ub,rawp;nout=nout,batch = prob.batch,kwargs...)
    dual = __solvebp_call(dp_prob,alg,sensealg,lb,ub,rawp,args...;kwargs...)
    res = similar(p, prob.nout)
    partials = reinterpret(typeof(first(res).partials), dual.u)
    for idx in eachindex(res)
        res[idx] = ForwardDiff.Dual{T,V,P}(primal.u[idx], partials[idx])
    end
    if primal.u isa Number
        res = first(res)
    end
    DiffEqBase.build_solution(prob,alg,res,primal.resid)
end

export QuadGKJL, HCubatureJL, VEGAS
export CubatureJLh, CubatureJLp
export CubaVegas, CubaSUAVE, CubaDivonne, CubaCuhre
end # module


text:
using JLD2,MCMCChains
import KenyaCoVSD
using StatsPlots,Distributions,Optim

## Group fitted models by semiurban/rural and rural

Semi_urban_filenames = ["modelfits/Nakuru_model.jld2",
                        "modelfits/Uasin Gishu_model.jld2",
                        "modelfits/Embu_model.jld2",
                        "modelfits/Kisumu_model.jld2",
                        "modelfits/Nyeri_model.jld2",
                        "modelfits/Siaya_model.jld2",
                         "modelfits/Kisii_model.jld2"]

rural_filenames = ["modelfits/Kilifi_model.jld2",
                        "modelfits/Kwale_model.jld2"]

## Combine posterior distribution draws across groups
fit_dict = load(Semi_urban_filenames[1])
model = fit_dict[first(keys(fit_dict))]
semi_urban_chain = model.MCMC_results.chain

for filename in Semi_urban_filenames[2:end]
    fit_dict = load(filename)
    model = fit_dict[first(keys(fit_dict))]
    semi_urban_chain = vcat(semi_urban_chain,model.MCMC_results.chain)
end

fit_dict = load(rural_filenames[1])
model = fit_dict[first(keys(fit_dict))]
rural_chain = model.MCMC_results.chain
for filename in rural_filenames[2:end]
    fit_dict = load(filename)
    model = fit_dict[first(keys(fit_dict))]
    rural_chain = vcat(rural_chain,model.MCMC_results.chain)
end

## Define a MLE method for creating a univariate fit for detection rate parameters across groups
"""
    fit_univariate_distributions_for_detection(post)

Fit a collection of univariate Gamma distributions to posterior samples of detection rate parameters.
"""
function fit_univariate_distributions_for_detection(post)
    v = get(post,[:p_test₁,:p_test₂,:χ₁,:χ₂])
    d_ptest1 = fit_mle(Gamma,v.p_test₁)
    d_ptest2 = fit_mle(Gamma,v.p_test₂)
    d_χ1 = fit_mle(Gamma,v.χ₁)
    d_χ2 = fit_mle(Gamma,v.χ₂)

    return (d_ptest1 = d_ptest1,
            d_ptest2= d_ptest2,
            d_χ1 = d_χ1,
            d_χ2 = d_χ2)
end

semi_urban_rural_fits = fit_univariate_distributions_for_detection(semi_urban_chain)
rural_fits = fit_univariate_distributions_for_detection(rural_chain)

@save("data/semi_urban_rural_detection_fits.jld2",semi_urban_rural_fits)
@save("data/rural_fits.jld2",rural_fits)

ext: jl


text:
export GraphNode, SplitNode, LeafNode

abstract type GraphNode end

"""
	mutable struct SplitNode <: GraphNode

Represent a split in the tree of a given segment. It is an internal node in the graph of splits.
- `color` is the color of the segment that `SplitNode` refers to
- `conf` defines the split. It is an array of `Bool` of length `N` (number of leaves). All leaves `i` such that `conf[i]` are on one side of the split, the others are on the other side.
- `anc::Union{SplitNode,Nothing}`
- `child::Array{GraphNode,1}`
- `isroot::Bool`
"""
mutable struct SplitNode <: GraphNode
	anc::Union{SplitNode,Nothing}
	child::Array{GraphNode,1}
	color::Int64
	conf::Array{Int,1}
	isroot::Bool
	SplitNode(anc, child, color, conf, isroot) = new(anc, child, color, sort(conf), isroot)
end
function SplitNode(;
	anc=nothing,
	child=Array{GraphNode,1}(undef, 0),
	color=0,
	conf=Array{Int,1}(undef,0),
	isroot=true
)
	return  SplitNode(anc, child, color, conf, isroot)
end

"""
	mutable struct LeafNode <: GraphNode

Represent a leaf in the graph of splits.
"""
mutable struct LeafNode <: GraphNode
	anc::Array{SplitNode,1}
	conf::Array{Int,1}
	index::Int
	LeafNode(anc, conf, index::Int) = begin
		if length(conf) == 1 && conf[1] == index
			return new(anc, conf, index)
		else
			error("Can't create LeafNode from conf=$conf and index=$index")
		end
	end
end
LeafNode(anc, index::Int) = LeafNode(anc, [index], index)
function LeafNode(;anc = Array{SplitNode}(undef, 0), conf = Int[0], index=0)
	return LeafNode(anc, conf, index)
end

"""
	mutable struct Graph

Graph of genealogic splits for a set of sequences consisting of many segments. Above leaf level, this graph consists of `K` trees of splits, where `K` is the number of segments.
"""
mutable struct Graph
	nodes::Array{GraphNode,1}
	leaves::Array{LeafNode,1}
	lleaves::Dict{String, LeafNode}
	internals::Array{SplitNode,1}
	labels::Array{String,1}
	labels_to_int::Dict{Any,Int64}
	K::Int64
end
function Graph(; nodes=Array{GraphNode,1}(undef,0), leaves=Array{LeafNode,1}(undef,0), lleaves=Dict{String, LeafNode}(), internals=Array{SplitNode,1}(undef,0), labels=Array{String,1}(undef,0), labels_to_int=Dict{Any,Int64}(), K=0)
	return Graph(nodes, leaves, lleaves, internals, labels, labels_to_int, K)
end

max_stars_repo_name: PierreBarrat/RecombTools

max_issues_repo_name: PierreBarrat/TreeKnit

max_forks_repo_name: PierreBarrat/TreeKnit


text:
using DomainSets, BasisFunctions, FrameFun.FrameFunDomains, Test

function test_extended_domains_3d()
    domains = (atomium(),)
end

function test_extended_domains_2d()
    D1 = characteristic(x->(x[1]-x[2])<0,disk(1.0))
    @test in([0.5;1.0],D1)
    @test ~in([1.0,0.5],D1)
    # basis = ChebyshevT(21)⊗ChebyshevT(21)
    # domain = disk(1.0)
    # f = (x,y)->exp(x+y)
    # F = Fun(f, basis, domain)
    # F2 = F*2
    # D2 = F<2
    # D3 = F<F2
    # @test ~in([0.5;0.5],D2)
    # @test in([0.25;0.25],D2)
    # @test in([0.5;0.5],D3)
    D4 = mandelbrot()
    @test in([0.0;0.0],D4)
    @test ~in([-0.5;0.5],D4)
    D5 = juliaset()
    @test in([0.6;0.0],D5)
    @test ~in([-0.1;-0.3],D5)
    # D6 = polardomain(x->1+0.2*cos(5x),disk(2.0))
    # @test ~in([1.25;0],D6)
    # @test in([0.79;0],D6)
end

domn(point,domain::Domain,normalvector; tol=1e-14) = norm(normal(point,domain)-normalvector)<tol
domb(point,domain::Domain; tol=1e-14) = abs(distance(point,domain))<tol

function test_distances_and_normals()
    D1 = simplex(Val{2})
    @test domb([0.5,0.5],D1)
    @test domn([0.5,0.5],D1,[1/sqrt(2),1/sqrt(2)])
    @test domb([0.5,0],D1)
    @test domn([0.5,0.0],D1,[0,-1])

    # D6 = polardomain(x->1+0.2*cos(5x),disk(2.0))
    # @test domb([-0.8,0],D6)
    # @test domn([-0.8,0],D6,[-1,0],tol=1e-12)

    D2 = disk(sqrt(0.5))
    @test domb([0.5,0.5],D2)
    @test domn([0.5,0.5],D2,[1/sqrt(2),1/sqrt(2)])

    D2b = circle(sqrt(0.5))
    @test domb([0.5,0.5],D2b)
    @test domn([0.5,0.5],D2b,[1/sqrt(2),1/sqrt(2)])

    D3 = -1..1.0
    @test isapprox(distance(0.5,D3),0.5)
    @test isapprox(normal(1.0,D3),1)

    D5 = cube(Val{2})
    @test domb([0.4,0.0],D5)
    @test domn([0.4,0.0],D5,[0,-1])

    # D4 = UnionDomain(D6, disk())\D2
    # @test domb([sqrt(0.5),0.0],D4)
    # @test domb([-1.0,0.0],D4)
    # @test ~domb([-0.8,0.0],D4)
    # @test domn([-1.0,0],D4,[-1,0])
    # @test domn([0.5,0.0],D4,[-1.0,0])

end

@testset "FrameFun Domains" begin

    test_extended_domains_3d()

    test_extended_domains_2d()

    test_distances_and_normals()
end


text:
export
    ChannelsList,
    ChannelsListResponse,
    ChannelsInfo,
    ChannelsInfoResponse,
    ChannelsArchive,
    ChannelsArchiveResponse,
    ChannelsJoin,
    ChannelsJoinResponse

@slackmethod(ChannelsList, "channels.list",
    begin
        exclude_archived::Nullable{Int64}
    end,

    begin
        channels::Array{SlackChannel,1}
    end)

@slackmethod(ChannelsInfo, "channels.info",
    begin
        channel::ChannelId
    end,

    begin
        channel::SlackChannel
    end)

@slackmethod(ChannelsArchive, "channels.archive",
    begin
        channel::ChannelId
    end,

    # Note: ChannelsArchiveResponse is intentionally empty, because we get no more information back than
    # the status.
    begin
    end)

@slackmethod(ChannelsJoin, "channels.join",
    begin
        name::UTF8String
    end,

    begin
        already_in_channel::Nullable{Bool}
        # TODO: The "channel" object here might be a limited channel type.
    end)


text:
export learn, learn!
export arc_requirements_do_not_conflict
export DSL_BayesianSearch,          LearnParams_BayesianSearch
export DSL_GreedyThickThinning,     LearnParams_GreedyThickThinning
export DSL_NaiveBayes,              LearnParams_NaiveBayes
export DSL_PC,                      LearnParams_PC
export DSL_TreeAugmentedNaiveBayes, LearnParams_TreeAugmentedNaiveBayes

# -------

const DSL_K2   = unsafe_load(cglobal((:EDSL_K2,   LIB_SMILE), Cint))
const DSL_BDeu = unsafe_load(cglobal((:EDSL_BDeu, LIB_SMILE), Cint))

# -------

abstract DSL_LearningAlgorithm
type DSL_BayesianSearch          <: DSL_LearningAlgorithm end
type DSL_GreedyThickThinning     <: DSL_LearningAlgorithm end
type DSL_NaiveBayes              <: DSL_LearningAlgorithm end
type DSL_PC                      <: DSL_LearningAlgorithm end
type DSL_TreeAugmentedNaiveBayes <: DSL_LearningAlgorithm end

# -------

abstract DSL_LearnParams
type LearnParams_BayesianSearch <: DSL_LearnParams
		maxparents           :: Int            # limits the maximum number of parents a node can have
		maxsearchtime        :: Int            # maximum runtime of the algorithm [seconds?] (0 = infinite)
		niterations          :: Int            # number of searches (and indirectly number of random restarts)
		linkprobability      :: Float64        # defines the probabilty of having an arc between two nodes
		priorlinkprobability :: Float64        # defines a prior existence of an arc between two nodes
		priorsamplesize      :: Int            # influences the "strength" of prior link probability.
		seed                 :: Int            # random seed (0=none)
		forced_arcs          :: Set{Tuple{Int,Int}} # a list of (i->j) arcs which are forced to be in the network
		forbidden_arcs       :: Set{Tuple{Int,Int}} # a list of (i->j) arcs which are forbidden in the network
		tiers                :: Set{Tuple{Int,Int}} # a list of (i->tier) associating nodes with a particular tier (tier starts at 1, variable starts at 0)

		LearnParams_BayesianSearch() = new(5, 0, 20, 0.1, 0.001, 50, 0, Set{Tuple{Int,Int}}(), Set{Tuple{Int,Int}}(), Set{Tuple{Int,Int}}())
end
type LearnParams_GreedyThickThinning <: DSL_LearnParams
		maxparents           :: Int # limits the maximum number of parents a node can have
		priors               :: Int # either K2 or BDeu
		netWeight            :: Float64
		forced_arcs          :: Set{Tuple{Int,Int}} # a list of (i->j) arcs which are forced to be in the network
		forbidden_arcs       :: Set{Tuple{Int,Int}} # a list of (i->j) arcs which are forbidden in the network
		tiers                :: Set{Tuple{Int,Int}} # a list of (i->tier) associating nodes with a particular tier
		LearnParams_GreedyThickThinning() = new(5, DSL_K2, 1.0, Set{Tuple{Int,Int}}(), Set{Tuple{Int,Int}}(), Set{Tuple{Int,Int}}())
end
type LearnParams_NaiveBayes <: DSL_LearnParams
	classVariableId :: AbstractString # the variable ID of the column that is our class

	LearnParams_NaiveBayes() = new("class")
	LearnParams_NaiveBayes(var::AbstractString) = new(var)
end
type LearnParams_PC <: DSL_LearnParams
	maxcache      :: UInt64
	maxAdjacency  :: Cint
	maxSearchTime :: Cint
	significance  :: Float64
	forced_arcs     :: Set{Tuple{Int,Int}} # a list of (i->j) arcs which are forced to be in the network
		forbidden_arcs  :: Set{Tuple{Int,Int}} # a list of (i->j) arcs which are forbidden in the network
		tiers           :: Set{Tuple{Int,Int}} # a list of (i->tier) associating nodes with a particular tier

	LearnParams_PC() = new(2048, 8, 0, 0.05, Set{Tuple{Int,Int}}(), Set{Tuple{Int,Int}}(), Set{Tuple{Int,Int}}())
end
type LearnParams_TreeAugmentedNaiveBayes <: DSL_LearnParams
	classvar      :: AbstractString # Used to pick the class variable for the network. Is case-sensitive.
	maxSearchTime :: Cint   # maximum runtime for the algorithm
	seed          :: UInt32 # random seed (0 for none)
	maxcache      :: UInt64 # maximum cache size (units?)

	LearnParams_TreeAugmentedNaiveBayes() = new("class", 0, 0, 2048)
	LearnParams_TreeAugmentedNaiveBayes(class) = new(class, 0, 0, 2048)
end

# ------

function learn!( net::Network, dset::Dataset, ::Type{DSL_BayesianSearch},
	params::LearnParams_BayesianSearch = LearnParams_BayesianSearch()
	)

	nforcedarcs,    forcedarcs_arr    = _get_arcs_arr(params.forced_arcs)
	nforbiddenarcs, forbiddenarcs_arr = _get_arcs_arr(params.forbidden_arcs)
	lentiers,       tiers_arr         = _get_arcs_arr(params.tiers)

	retval = ccall( (:dataset_learnBayesianSearch, LIB_SMILE), Bool,
		(Ptr{Void},Ptr{Void},Int32,Int32,Int32,Float64,Float64,Int32,Int32,
			Ptr{Int32},Int32,Ptr{Int32},Int32,Ptr{Int32},Int32),
		dset.ptr, net.ptr, params.maxparents, params.maxsearchtime, params.niterations,
		params.linkprobability, params.priorlinkprobability, params.priorsamplesize,
		params.seed, forcedarcs_arr, nforcedarcs,
		forbiddenarcs_arr, nforbiddenarcs, tiers_arr, lentiers )
	if ( !retval )
		warn("Dataset::learn_bayesian_search: did not work")
	end
	return retval # true on success
end
function learn!( net::Network, dset::Dataset, ::Type{DSL_GreedyThickThinning},
	params::LearnParams_GreedyThickThinning = LearnParams_GreedyThickThinning()
	)

	nforcedarcs,    forcedarcs_arr    = _get_arcs_arr(params.forced_arcs)
	nforbiddenarcs, forbiddenarcs_arr = _get_arcs_arr(params.forbidden_arcs)
	lentiers,       tiers_arr         = _get_arcs_arr(params.tiers)

	retval = ccall( (:dataset_learnGreedyThickThinning, LIB_SMILE), Bool,
		(Ptr{Void},Ptr{Void},Int32,Int32,Float64,
			Ptr{Int32},Int32,Ptr{Int32},Int32,Ptr{Int32},Int32),
		dset.ptr, net.ptr, params.priors, params.maxparents,
		params.netWeight, forcedarcs_arr, nforcedarcs,
		forbiddenarcs_arr, nforbiddenarcs, tiers_arr, lentiers )
	if ( !retval )
		warn("Dataset::learn_greedy_thick_thinning: did not work")
	end
	return retval # true on success
end
function learn!( net::Network, dset::Dataset, ::Type{DSL_NaiveBayes},
	params::LearnParams_NaiveBayes = LearnParams_NaiveBayes()
	)

	if find_variable(dset, params.classVariableId) == -1
		warn("class variable $(params.classVariableId) could not be found")
	end

	retval = ccall( (:dataset_learnNaiveBayes, LIB_SMILE), Bool,
		(Ptr{Void},Ptr{Void},Ptr{UInt8}),
		dset.ptr, net.ptr, bytestring(params.classVariableId) )
	if ( !retval )
		warn("Dataset::learn_naive_bayes: did not work")
	end
	return retval # true on success
end
function learn!( net::Network, dset::Dataset, ::Type{DSL_TreeAugmentedNaiveBayes},
	params::LearnParams_TreeAugmentedNaiveBayes = LearnParams_TreeAugmentedNaiveBayes()
	)

	if find_variable(dset, params.classvar) == -1
		warn("class variable $(params.classvar) could not be found")
	end

	retval = ccall( (:dataset_learnTAN, LIB_SMILE), Bool,
		(Ptr{Void},Ptr{Void},Ptr{UInt8}, Cint, UInt32, UInt64),
		dset.ptr, net.ptr, bytestring(params.classvar), params.maxSearchTime, params.seed, params.maxcache )
	if ( !retval )
		warn("Dataset::learn_TAN: did not work")
	end
	return retval # true on success
end
function learn!( net::Network, dset::Dataset, ::Type{DSL_PC},
	params::LearnParams_PC = LearnParams_PC()
	)

	pat = Pattern()
	worked = learn!( pat, dset, DSL_PC, params )
	to_network(pat, dset, net)
	return worked
end
function learn!( pat::Pattern, dset::Dataset, ::Type{DSL_PC},
	params::LearnParams_PC = LearnParams_PC()
	)

	nforcedarcs,    forcedarcs_arr    = _get_arcs_arr(params.forced_arcs)
	nforbiddenarcs, forbiddenarcs_arr = _get_arcs_arr(params.forbidden_arcs)
	lentiers,       tiers_arr         = _get_arcs_arr(params.tiers)

	retval = ccall( (:dataset_learnPC, LIB_SMILE), Bool,
		(Ptr{Void},Ptr{Void},UInt64,Cint,Cint,Float64,
			Ptr{Int32},Int32,Ptr{Int32},Int32,Ptr{Int32},Int32),
		dset.ptr, pat.ptr, params.maxcache, params.maxAdjacency,
		params.maxSearchTime, params.significance,
		forcedarcs_arr, nforcedarcs, forbiddenarcs_arr,
		nforbiddenarcs, tiers_arr, lentiers)
	if ( !retval )
		warn("Dataset::dataset_learnPC: did not work")
	end
	return retval # true on success
end

learn!( net::Network, dset::Dataset, params::LearnParams_BayesianSearch ) =
	learn!(net, dset, DSL_BayesianSearch, params)
learn!( net::Network, dset::Dataset, params::LearnParams_GreedyThickThinning ) =
	learn!(net, dset, DSL_GreedyThickThinning, params)
learn!( net::Network, dset::Dataset, params::LearnParams_NaiveBayes ) =
	learn!(net, dset, DSL_NaiveBayes, params)
learn!( net::Network, dset::Dataset, params::LearnParams_TreeAugmentedNaiveBayes ) =
	learn!(net, dset, DSL_TreeAugmentedNaiveBayes, params)
learn!( net::Network, dset::Dataset, params::LearnParams_PC ) =
	learn!(net, dset, DSL_PC, params)
learn!( pat::Pattern, dset::Dataset, params::LearnParams_PC ) =
	learn!(pat, dset, DSL_PC, params)

function learn{S<:DSL_LearningAlgorithm}( dset::Dataset, alg::Type{S} )
	net = Network()
	worked = learn!(net, dset, alg)
	(net, worked)
end
function learn{S<:DSL_LearnParams}( dset::Dataset, params::S )
	net = Network()
	worked = learn!(net, dset, params)
	(net, worked)
end

# -------

function arc_requirements_do_not_conflict(
	forced_arcs    :: Set{Tuple{Int,Int}},
	forbidden_arcs :: Set{Tuple{Int,Int}}
	)
	# returns true if a forced arc is also a forbidden arc
	for tup in params.forbidden_arcs
		if in(tup, params.forced_arcs)
			return false
		end
	end
	true
end
function arc_requirements_do_not_conflict(
	forced_arcs    :: Set{Tuple{Int,Int}},
	forbidden_arcs :: Set{Tuple{Int,Int}},
	tiers          :: Set{Tuple{Int,Int}},
	)
	# converts tiers to arc requirements before checking conflict
	forced, forbidden = _tiers_to_arc_requirements(tiers)
	union!(forced, forced_arcs)
	union!(forbidden, forbidden_arcs)
	arc_requirements_do_not_conflict(forced, forbidden)
end
arc_requirements_do_not_conflict( params::DSL_BayesianSearch ) =
	arc_requirements_do_not_conflict(params.forced_arcs, params.forbidden_arcs, params.tiers)
arc_requirements_do_not_conflict( params::DSL_GreedyThickThinning ) =
	arc_requirements_do_not_conflict(params.forced_arcs, params.forbidden_arcs, params.tiers)
arc_requirements_do_not_conflict( params::DSL_BayesianSearch ) =
	arc_requirements_do_not_conflict(params.forced_arcs, params.forbidden_arcs, params.tiers)

# -------

function _get_arcs_arr(arcs::Set{Tuple{Int,Int}})
	narcs = length(arcs)
	arr = zeros(Int32, 2*narcs)
	for (i,tup) in enumerate(arcs)
		arr[2*i-1] = tup[1]
		arr[2*i] = tup[2]
	end
	(narcs, arr)
end
function _tiers_to_arc_requirements(
	tiers :: Set{Tuple{Int,Int}} # a list of (i->tier) associating nodes with a particular tier
	)

	# the node starts counting at 0
	# the tier starts counting at 1

	forced_arcs    :: Set{Tuple{Int,Int}}()
	forbidden_arcs :: Set{Tuple{Int,Int}}()

	for (i,tier_i) in tiers
		for (j,tier_j) in tiers
			if tier_i > tier_j
				push!(forbidden_arcs, (i,j))
			elseif tier_j > tier_i
				push!(forbidden_arcs, (j,i))
			end
		end
	end

	(forced_arcs, forbidden_arcs)
end


text:
"""

## invar

This subroutine forms the stress invariants in 2- or 3-d. See equations
6.3 and 6.4

### Function
```julia
invar(stress, sigm, dsbar, theta)
```

### Arguments
```julia
* stress::Vector{Float64}    : Stress vector
* sigm::Float64              : Invariant, eq 6.4 (Updated)
* dsbar::Float64             : Invariant, eq 6.4 (Updated)
* theta::Float64             : Invariant, eq 6.3 (Updated)
```

### REturn values
```julia
* stress::Vector{Float64}    : Stress vector
* sigm::Float64              : Invariant, eq 6.4 (Updated)
* dsbar::Float64             : Invariant, eq 6.4 (Updated)
```
"""
function invar(stress::Vector{Float64}, sigm::Float64, dsbar::Float64,
  theta::Float64)
  #
  # This subroutine forms the stress invariants in 2- or 3-d.
  #
  # Updates: sigm, dsbar, theta
  #
 nst = size(stress, 1)

 if nst == 4
   sx = stress[1]
   sy = stress[2]
   txy =stress[3]
   sz = stress[4]
   sigm =(sx+sy+sz)/3.0
   dsbar = sqrt((sx-sy)^2+(sy-sz)^2+(sz-sx)^2+6.0*txy^2)/sqrt(2.0)
   if dsbar < 1.0e-10
     theta  =0.0
   else
     dx = (2.0*sx-sy-sz)/3.0
     dy = (2.0*sy-sz-sx)/3.0
     dz = (2.0*sz-sx-sy)/3.0
     xj3 = dx*dy*dz-dz*txy^2
     sine = -13.5*xj3/dsbar^3
     if sine >= 1.0
       sine = 1.0
     end
     if sine < -1.0
       sine = -1.0
     end
     theta = asin(sine)/3.0
   end
 elseif nst == 6
   sq3 = sqrt(3.0)
   s1 = stress[1]
   s2 = stress[2]
   s3 = stress[3]
   s4 = stress[4]
   s5 = stress[5]
   s6 = stress[6]
   sigm = (s1+s2+s3)/3.0
   d2 = ((s1-s2)^2+(s2-s3)^2+(s3-s1)^2)/6.0+s4*s4+s5*s5+s6*s6
   ds1 = s1-sigm
   ds2 = s2-sigm
   ds3 = s3-sigm
   d3 = ds1*ds2*ds3-ds1*s5*s5-ds2*s6*s6-ds3*s4*s4+2.0*s4*s5*s6
   dsbar = sq3*sqrt(d2)
   if dsbar < 1.0e-10
     theta = 0.0
   else
     sine = -3.0*sq3*3.0/(2.0*sqrt(2.0)^3)
     if sine >= 1.0
       sine = 1.0
     end
     if sine < -1.0
       sine = -1.0
     end
     theta = asin(sine)/3.0
   end
 else
   println("Wrong size for nst in invar.")
 end
 (sigm, dsbar, theta)
end


text:
module CellMapsTests

using Test
using CellwiseValues
using TensorValues

using ..CellValuesMocks
using ..MapsMocks

l = 10

a = VectorValue(10,10)
b = VectorValue(15,20)
p1 = VectorValue(1,1)
p2 = VectorValue(2,2)
p3 = VectorValue(3,3)
p = [p1,p2,p3]
m = MockMap(a)
r = evaluate(m,p)

cm = TestIterCellValue(m,l)
cp = TestIterCellValue(p,l)
rm = [ CachedArray(r) for i in 1:l]
test_iter_cell_map(cm,cp,rm)

cm = TestIndexCellValue(m,l)
cp = TestIndexCellValue(p,l)
rm = [ CachedArray(r) for i in 1:l]
test_index_cell_map_with_index_arg(cm,cp,rm)

end # module


text:
module MTHE493

using Images
using Statistics
using LinearAlgebra

include("eigenfaces.jl")
include("classification.jl")
include("utils.jl")

export  load_image,
        load_images,
        images_to_image_matrix,
        image_matrix_to_images,
        tile_images,
        rescale,
        get_projection,
        get_eigenfaces,
        image_to_eigenface_rep,
        eigenface_rep_to_image,
        euclidean_distance,
        k_nearest_neighbour,
        test_model

end


text:
module QuantumOpsTests

using QuantumOps
using ReTest
using StaticArrays
import SparseArraysN
import LinearAlgebra
const SparseArrays = SparseArraysN
using QuantumOps: Paulis, FermiOps

using Aqua
Aqua.test_all(QuantumOps; ambiguities=false)
Aqua.test_ambiguities(QuantumOps)

include("./test_pauli.jl")
include("./test_jordan_wigner.jl")
include("./test_op_term_sum.jl")
include("./test_fermi.jl")
include("./test_sparse.jl")
include("./test_z4group.jl")

# DISABLE this because it confuses PackageCompiler.
# @testset "quantum_ops_intro" begin
#     include("../examples/quantum_ops_intro.jl")
# end

# @testset "jordan_wigner_example" begin
#     include("../examples/jordan_wigner_example.jl")
# end

end # module QuantumOpsTests


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

run polynomial utilities tests
=#

using Test
using Printf
include(joinpath(@__DIR__, "polyutils.jl"))

real_types = [
    Float64,
    Float32,
    BigFloat,
    ]

@testset "polynomial utilities tests" begin
@testset "$T" for T in real_types
    println("$T ...")
    test_time = @elapsed begin
        test_interp_domain(T)
        test_complex_interp(T)
        test_quadrature(T)
    end
    @printf("%8.2e seconds\n", test_time)
end
end
;


text:
using SparseArrays
using LinearAlgebra
include("LeituraDados.jl")

"""
Corrige os números dos graus de liberdade da notação do Python para a notação
de Julia. Soma 1 aos graus de liberdade não impedidos.
"""
function corrigir_graus_lib(x)
    if x != -1
        x + 1
    else
        x
    end
end

"""
Faz a montagem da matriz de rigidez da estrutura. Se retornar vetores for true, retorna
os três vetores que compõem a matriz esparsa. Caso contrário, retorna a matriz esparsa.
No caso da saída com vetores, o retorno é feito com uma lista com linhas, colunas e termos,
respectivamente.
"""
function matriz_rigidez_estrutura(kelems, dados::Dict, retornar_vetores::Bool=false)
    gls_elems = dados[ctes[2]]
    gls_est = dados[ctes[3]]
    apoios = dados[ctes[4]]

    # Tamanho dos vetores
    num_it = sum([size(ke)[1]^2 for ke in kelems])

    linhas = ones(Int64, num_it)
    colunas = ones(Int64, num_it)
    termos = zeros(num_it)

    # Número de graus de liberdade livres
    num_gls_livres = length(gls_est) - length(apoios)

    t = 1

    for i in range(1, stop=length(gls_elems))
        gle = gls_elems[i]
        ke = kelems[i]

        # Graus de liberdade do elemento com apoios
        gle_com_apoios = gls_est[gle]

        # Posições livres do elemento
        pos_livres = [i for (i, gl) in enumerate(gle_com_apoios) if gl != -1]

        # Graus de liberdade livres na forma B
        gle_livres = gle_com_apoios[pos_livres]

        # Matriz de rigidez com os graus de liberdade impedidos
        ke_livre = ke[pos_livres, pos_livres]

        num_gle_livres = length(gle_livres)

        if length(gle_livres) > 0
            for ii in range(1, stop=num_gle_livres)
                for jj in range(1, stop=num_gle_livres)
                    linhas[t] = gle_livres[ii]
                    colunas[t] = gle_livres[jj]
                    termos[t] = ke_livre[ii, jj]
                    t += 1
                end
            end
        end
    end
    k = sparse(linhas, colunas, termos)
    dropzeros!(k)
    if retornar_vetores
        return linhas, colunas, termos
    else
        return k
    end
end

"""
Resolve o sistema de equações e retorna os deslocamentos nodais.
"""
function deslocamentos(kelems, dados)
    k = matriz_rigidez_estrutura(kelems, dados)

    forcas = dados["forcas"]
    apoios = dados["apoios"]

    gls_livres = [i for (i, f) in enumerate(forcas) if i ∉ apoios]
    forcas_com_apoios = forcas[gls_livres]

    # Aplicação do RCM
    rcm = dados["RCM"]
    k_rcm = k[rcm, rcm]
    forcas_rcm = forcas_com_apoios[rcm]

    # Solução do sistema
    u0 = k_rcm \ forcas_rcm

    # Reversão do RCM
    rcm_rev = zeros(Int64, length(rcm))

    for (i, v) in enumerate(rcm)
        rcm_rev[v] = i
    end

    u1 = u0[rcm_rev]

    u = zeros(length(forcas))
    u[gls_livres] = u1

    return u
end

max_stars_repo_name: lorranfoliveira/otm

max_issues_repo_name: lorranfoliveira/otm

max_forks_repo_name: lorranfoliveira/otm


text:
module Testbyrow

using Test
using DataFrames
using DataFramesMeta
using Statistics

const ≅ = isequal

@testset "@transform with @byrow" begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :transform, missing])

    @test @transform(df, @byrow :n = :i + :g) ≅ @transform(df, :n = :i + :g)
    @test @transform(df, @byrow :n = :t * string(:y)) ≅ @transform(df, :n = :t .* string.(:y))
    @test @transform(df, @byrow :n = :g == 1 ? 100 : 500) ≅ @transform(df, :n = ifelse.(:g .== 1, 100, 500))
    @test @transform(df, @byrow :n = :g == 1 && :t == "a") ≅ @transform(df, :n = map((g, t) -> g == 1 && t == "a", :g, :t))
    @test @transform(df, @byrow :n = first(:g)) ≅ @transform(df, :n = first.(:g))

    d = @transform df @byrow begin
        :n1 = :i
        :n2 = :i * :g
    end
    @test d ≅ @transform(df, :n1 = :i, :n2 = :i .* :g)
    @test d ≅ @transform(df, @byrow(:n1 = :i), @byrow(:n2 = :i * :g))

    d = @transform df @byrow begin
        $:n1 = :i
        :n2 = $"i" * :g
    end
    @test d ≅ @transform(df, :n1 = :i, :n2 = :i .* :g)
    d = @transform df @byrow begin
        :n1 = $"i"
        $:n2 = :i * :g
    end
    @test d ≅ @transform(df, :n1 = :i, :n2 = :i .* :g)

    d = @transform df @byrow begin
        :n1 = begin
            :i
        end
        :n2 = :i * :g
    end
    @test d ≅ @transform(df, :n1 = :i, :n2 = :i .* :g)

    d = @transform df @byrow begin
        :n1 = :i * :g
        :n2 = :i * :g
    end
    @test d ≅ @transform(df, :n1 = :i .* :g, :n2 = :i .* :g)
end

@testset "@rtransform" begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :transform, missing])

    @test @rtransform(df, :n = :i + :g) ≅ @transform(df, :n = :i + :g)
    @test @rtransform(df, :n = :t * string(:y)) ≅ @transform(df, :n = :t .* string.(:y))
    @test @rtransform(df, :n = :g == 1 ? 100 : 500) ≅ @transform(df, :n = ifelse.(:g .== 1, 100, 500))
    @test @rtransform(df, :n = :g == 1 && :t == "a") ≅ @transform(df, :n = map((g, t) -> g == 1 && t == "a", :g, :t))
    @test @rtransform(df, :n = first(:g)) ≅ @transform(df, :n = first.(:g))

    d = @rtransform df begin
        :n1 = :i
        :n2 = :i * :g
    end
    @test d ≅ @transform(df, :n1 = :i, :n2 = :i .* :g)
    @test d ≅ @transform(df, @byrow(:n1 = :i), @byrow(:n2 = :i * :g))

    d = @rtransform df begin
        cols(:n1) = :i
        :n2 = cols(:i) * :g
    end
    @test d ≅ @transform(df, :n1 = :i, :n2 = :i .* :g)

    d = @rtransform df begin
        :n1 = cols(:i)
        cols(:n2) = :i * :g
    end
    @test d ≅ @transform(df, :n1 = :i, :n2 = :i .* :g)

    d = @rtransform df begin
        :n1 = begin
            :i
        end
        :n2 = :i * :g
    end
    @test d ≅ @transform(df, :n1 = :i, :n2 = :i .* :g)

    d = @rtransform df begin
        :n1 = :i * :g
        :n2 = :i * :g
    end
    @test d ≅ @transform(df, :n1 = :i .* :g, :n2 = :i .* :g)
end

@testset "@transform! with @byrow" begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :transform!, missing])

    @test @transform!(copy(df), @byrow :n = :i + :g) ≅ @transform!(copy(df), :n = :i + :g)
    @test @transform!(copy(df), @byrow :n = :t * string(:y)) ≅ @transform!(copy(df), :n = :t .* string.(:y))
    @test @transform!(copy(df), @byrow :n = :g == 1 ? 100 : 500) ≅ @transform!(copy(df), :n = ifelse.(:g .== 1, 100, 500))
    @test @transform!(copy(df), @byrow :n = :g == 1 && :t == "a") ≅ @transform!(copy(df), :n = map((g, t) -> g == 1 && t == "a", :g, :t))
    @test @transform!(copy(df), @byrow :n = first(:g)) ≅ @transform!(copy(df), :n = first.(:g))

    d = @transform! df @byrow begin
        :n1 = :i
        :n2 = :i * :g
    end
    @test d ≅ @transform!(copy(df), :n1 = :i, :n2 = :i .* :g)
    @test d ≅ @transform!(copy(df), @byrow(:n1 = :i), @byrow(:n2 = :i * :g))

    d = @transform! df @byrow begin
        $:n1 = :i
        :n2 = $"i" * :g
    end
    @test d ≅ @transform!(copy(df), :n1 = :i, :n2 = :i .* :g)
    d = @transform! df @byrow begin
        :n1 = $"i"
        $:n2 = :i * :g
    end
    @test d ≅ @transform!(copy(df), :n1 = :i, :n2 = :i .* :g)

    d = @transform! df @byrow begin
        :n1 = begin
            :i
        end
        :n2 = :i * :g
    end
    @test d ≅ @transform!(copy(df), :n1 = :i, :n2 = :i .* :g)

    d = @transform! df @byrow begin
        :n1 = :i * :g
        :n2 = :i * :g
    end
    @test d ≅ @transform!(copy(df), :n1 = :i .* :g, :n2 = :i .* :g)
end

@testset "@rtransform!" begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :transform!, missing])

    @test @rtransform!(copy(df), :n = :i + :g) ≅ @transform!(copy(df), :n = :i + :g)
    @test @rtransform!(copy(df), :n = :t * string(:y)) ≅ @transform!(copy(df), :n = :t .* string.(:y))
    @test @rtransform!(copy(df), :n = :g == 1 ? 100 : 500) ≅ @transform!(copy(df), :n = ifelse.(:g .== 1, 100, 500))
    @test @rtransform!(copy(df), :n = :g == 1 && :t == "a") ≅ @transform!(copy(df), :n = map((g, t) -> g == 1 && t == "a", :g, :t))
    @test @rtransform!(copy(df), :n = first(:g)) ≅ @transform!(copy(df), :n = first.(:g))

    df2 = copy(df)
    d = @rtransform! df2 begin
        :n1 = :i
        :n2 = :i * :g
    end
    @test d === df2
    @test d ≅ @transform!(copy(df), :n1 = :i, :n2 = :i .* :g)
    @test d ≅ @transform!(copy(df), @byrow(:n1 = :i), @byrow(:n2 = :i * :g))

    df2 = copy(df)
    d = @rtransform! df2 begin
        cols(:n1) = :i
        :n2 = cols(:i) * :g
    end
    @test d === df2
    @test d ≅ @transform!(copy(df), :n1 = :i, :n2 = :i .* :g)
    d = @rtransform! df begin
        :n1 = cols(:i)
        cols(:n2) = :i * :g
    end
    @test d ≅ @transform!(copy(df), :n1 = :i, :n2 = :i .* :g)

    d = @rtransform! df begin
        :n1 = begin
            :i
        end
        :n2 = :i * :g
    end
    @test d ≅ @transform!(copy(df), :n1 = :i, :n2 = :i .* :g)

    d = @rtransform! df begin
        :n1 = :i * :g
        :n2 = :i * :g
    end
    @test d ≅ @transform!(copy(df), :n1 = :i .* :g, :n2 = :i .* :g)
end

@testset "@select with @byrow" begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :select, missing])

    @test @select(df, @byrow :n = :i + :g) ≅ @select(df, :n = :i + :g)
    @test @select(df, @byrow :n = :t * string(:y)) ≅ @select(df, :n = :t .* string.(:y))
    @test @select(df, @byrow :n = :g == 1 ? 100 : 500) ≅ @select(df, :n = ifelse.(:g .== 1, 100, 500))
    @test @select(df, @byrow :n = :g == 1 && :t == "a") ≅ @select(df, :n = map((g, t) -> g == 1 && t == "a", :g, :t))
    @test @select(df, @byrow :n = first(:g)) ≅ @select(df, :n = first.(:g))

    d = @select df @byrow begin
        :n1 = :i
        :n2 = :i * :g
    end
    @test d ≅ @select(df, :n1 = :i, :n2 = :i .* :g)
    @test d ≅ @select(df, @byrow(:n1 = :i), @byrow(:n2 = :i * :g))

    d = @select df @byrow begin
        $:n1 = :i
        :n2 = $"i" * :g
    end
    @test d ≅ @select(df, :n1 = :i, :n2 = :i .* :g)
    d = @select df @byrow begin
        :n1 = $"i"
        $:n2 = :i * :g
    end
    @test d ≅ @select(df, :n1 = :i, :n2 = :i .* :g)

    d = @select df @byrow begin
        :n1 = begin
            :i
        end
        :n2 = :i * :g
    end
    @test d ≅ @select(df, :n1 = :i, :n2 = :i .* :g)

    d = @select df @byrow begin
        :n1 = :i * :g
        :n2 = :i * :g
    end
    @test d ≅ @select(df, :n1 = :i .* :g, :n2 = :i .* :g)
end

@testset "@rselect" begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :select, missing])

    @test @rselect(df, :n = :i + :g) ≅ @select(df, :n = :i + :g)
    @test @rselect(df, :n = :t * string(:y)) ≅ @select(df, :n = :t .* string.(:y))
    @test @rselect(df, :n = :g == 1 ? 100 : 500) ≅ @select(df, :n = ifelse.(:g .== 1, 100, 500))
    @test @rselect(df, :n = :g == 1 && :t == "a") ≅ @select(df, :n = map((g, t) -> g == 1 && t == "a", :g, :t))
    @test @rselect(df, :n = first(:g)) ≅ @select(df, :n = first.(:g))

    d = @rselect df begin
        :n1 = :i
        :n2 = :i * :g
    end
    @test d ≅ @select(df, :n1 = :i, :n2 = :i .* :g)
    @test d ≅ @select(df, @byrow(:n1 = :i), @byrow(:n2 = :i * :g))

    d = @rselect df begin
        cols(:n1) = :i
        :n2 = cols(:i) * :g
    end
    @test d ≅ @select(df, :n1 = :i, :n2 = :i .* :g)
    d = @rselect df begin
        :n1 = cols(:i)
        cols(:n2) = :i * :g
    end
    @test d ≅ @select(df, :n1 = :i, :n2 = :i .* :g)

    d = @rselect df begin
        :n1 = begin
            :i
        end
        :n2 = :i * :g
    end
    @test d ≅ @select(df, :n1 = :i, :n2 = :i .* :g)

    d = @rselect df begin
        :n1 = :i * :g
        :n2 = :i * :g
    end
    @test d ≅ @select(df, :n1 = :i .* :g, :n2 = :i .* :g)
end

@testset "@select! with @byrow" begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :select!, missing])

    @test @select!(copy(df), @byrow :n = :i + :g) ≅ @select!(copy(df), :n = :i + :g)
    @test @select!(copy(df), @byrow :n = :t * string(:y)) ≅ @select!(copy(df), :n = :t .* string.(:y))
    @test @select!(copy(df), @byrow :n = :g == 1 ? 100 : 500) ≅ @select!(copy(df), :n = ifelse.(:g .== 1, 100, 500))
    @test @select!(copy(df), @byrow :n = :g == 1 && :t == "a") ≅ @select!(copy(df), :n = map((g, t) -> g == 1 && t == "a", :g, :t))
    @test @select!(copy(df), @byrow :n = first(:g)) ≅ @select!(copy(df), :n = first.(:g))

    d = @select! copy(df) @byrow begin
        :n1 = :i
        :n2 = :i * :g
    end
    @test d ≅ @select!(copy(df), :n1 = :i, :n2 = :i .* :g)
    @test d ≅ @select!(copy(df), @byrow(:n1 = :i), @byrow(:n2 = :i * :g))

    d = @select! copy(df) @byrow begin
        $:n1 = :i
        :n2 = $"i" * :g
    end
    @test d ≅ @select!(copy(df), :n1 = :i, :n2 = :i .* :g)
    d = @select! copy(df) @byrow begin
        :n1 = $"i"
        $:n2 = :i * :g
    end
    @test d ≅ @select!(copy(df), :n1 = :i, :n2 = :i .* :g)

    d = @select! copy(df) @byrow begin
        :n1 = begin
            :i
        end
        :n2 = :i * :g
    end
    @test d ≅ @select!(copy(df), :n1 = :i, :n2 = :i .* :g)

    d = @select! copy(df) @byrow begin
        :n1 = :i * :g
        :n2 = :i * :g
    end
    @test d ≅ @select!(copy(df), :n1 = :i .* :g, :n2 = :i .* :g)
end

@testset "@rselect!" begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :select!, missing])

    @test @rselect!(copy(df), :n = :i + :g) ≅ @select!(copy(df), :n = :i + :g)
    @test @rselect!(copy(df), :n = :t * string(:y)) ≅ @select!(copy(df), :n = :t .* string.(:y))
    @test @rselect!(copy(df), :n = :g == 1 ? 100 : 500) ≅ @select!(copy(df), :n = ifelse.(:g .== 1, 100, 500))
    @test @rselect!(copy(df), :n = :g == 1 && :t == "a") ≅ @select!(copy(df), :n = map((g, t) -> g == 1 && t == "a", :g, :t))
    @test @rselect!(copy(df), :n = first(:g)) ≅ @select!(copy(df), :n = first.(:g))

    df2 = copy(df)
    d = @rselect! df2 begin
        :n1 = :i
        :n2 = :i * :g
    end
    @test d === df2
    @test d ≅ @select!(copy(df), :n1 = :i, :n2 = :i .* :g)
    @test d ≅ @select!(copy(df), @byrow(:n1 = :i), @byrow(:n2 = :i * :g))

    df2 = copy(df)
    d = @rselect! df2 begin
        cols(:n1) = :i
        :n2 = cols(:i) * :g
    end
    @test d === df2
    @test d ≅ @select!(copy(df), :n1 = :i, :n2 = :i .* :g)
    d = @rselect! copy(df) begin
        :n1 = cols(:i)
        cols(:n2) = :i * :g
    end
    @test d ≅ @select!(copy(df), :n1 = :i, :n2 = :i .* :g)

    d = @rselect! copy(df) begin
        :n1 = begin
            :i
        end
        :n2 = :i * :g
    end
    @test d ≅ @select!(copy(df), :n1 = :i, :n2 = :i .* :g)

    d = @rselect! copy(df) begin
        :n1 = :i * :g
        :n2 = :i * :g
    end
    @test d ≅ @select!(copy(df), :n1 = :i .* :g, :n2 = :i .* :g)
end

@testset "@with with @byrow" begin
    df = DataFrame(A = 1:3, B = [2, 1, 2])

    @test @with(df, @byrow :A * 1)   ==  df.A .* 1
    @test @with(df, @byrow :A * :B)  ==  df.A .* df.B

    t = @with df @byrow begin
        :A * 1
    end
    @test t == df.A .* 1

    t = @with df @byrow begin
        :A * :B
    end
    @test t == df.A .* df.B
end

@testset "@subset with @byrow" begin
    df = DataFrame(A = [1, 2, 3, missing], B = [2, 1, 2, 1])

    d = @subset df begin
        @byrow :A > 1
        @byrow :B > 1
    end
    @test d ≅ @subset(df, :A .> 1, :B .> 1)

    d = @subset df @byrow begin
        :A > 1
        :B > 1
    end
    @test d ≅ @subset(df, :A .> 1, :B .> 1)
end

@testset "@rsubset" begin
    df = DataFrame(A = [1, 2, 3, missing], B = [2, 1, 2, 1])

    d = @rsubset df begin
        :A > 1
        :B > 1
    end
    @test d ≅ @subset(df, :A .> 1, :B .> 1)

    d = @rsubset df begin
        :A > 1
        :B > 1
    end
    @test d ≅ @subset(df, :A .> 1, :B .> 1)
end

@testset "@subset! with @byrow" begin
    df = DataFrame(A = [1, 2, 3, missing], B = [2, 1, 2, 1])

    df2 = copy(df)
    d = @subset! df2 begin
        @byrow :A > 1
        @byrow :B > 1
    end
    @test d === df2
    @test d ≅ @subset(df, :A .> 1, :B .> 1)

    d = @subset! copy(df) @byrow begin
        :A > 1
        :B > 1
    end
    @test d ≅ @subset(df, :A .> 1, :B .> 1)
end

@testset "@rsubset!" begin
    df = DataFrame(A = [1, 2, 3, missing], B = [2, 1, 2, 1])

    df2 = copy(df)
    d = @rsubset! df2 begin
        :A > 1
        :B > 1
    end
    @test d === df2
    @test d ≅ @subset(df, :A .> 1, :B .> 1)

    d = @rsubset! copy(df) begin
        :A > 1
        :B > 1
    end
    @test d ≅ @subset(df, :A .> 1, :B .> 1)
end

@testset "@orderby with @byrow" begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :transform, missing]
        )

    d = @orderby df begin
        @byrow :c
        @byrow :g *  2
    end
    @test d ≅ @orderby(df, :c, :g .* 2)

    d = @orderby df @byrow begin
        :c
        :g *  2
    end
    @test d ≅ @orderby(df, :c, :g .* 2)
end

@testset "@rorderby " begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :transform, missing]
        )

    d = @rorderby df begin
        :c
        :g *  2
    end
    @test d ≅ @orderby(df, :c, :g .* 2)

    d = @rorderby df begin
        :c
        :g *  2
    end
    @test d ≅ @orderby(df, :c, :g .* 2)
end

@testset "@combine with @byrow" begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :combine, missing])

    gd = groupby(df, :g)

    @test @combine(gd, @byrow :n = :i + :g) ≅ @combine(gd, :n = :i + :g)
    @test @combine(gd, @byrow :n = :t * string(:y)) ≅ @combine(gd, :n = :t .* string.(:y))
    @test @combine(gd, @byrow :n = :g == 1 ? 100 : 500) ≅ @combine(gd, :n = ifelse.(:g .== 1, 100, 500))
    @test @combine(gd, @byrow :n = :g == 1 && :t == "a") ≅ @combine(gd, :n = map((g, t) -> g == 1 && t == "a", :g, :t))
    @test @combine(gd, @byrow :n = first(:g)) ≅ @combine(gd, :n = first.(:g))

    d = @combine gd @byrow begin
        :n1 = :i
        :n2 = :i * :g
    end
    @test d ≅ @combine(gd, :n1 = :i, :n2 = :i .* :g)
    @test d ≅ @combine(gd, @byrow(:n1 = :i), @byrow(:n2 = :i * :g))

    d = @combine gd @byrow begin
        $:n1 = :i
        :n2 = $"i" * :g
    end
    @test d ≅ @combine(gd, :n1 = :i, :n2 = :i .* :g)
    d = @combine gd @byrow begin
        :n1 = $"i"
        $:n2 = :i * :g
    end
    @test d ≅ @combine(gd, :n1 = :i, :n2 = :i .* :g)

    d = @combine gd @byrow begin
        :n1 = begin
            :i
        end
        :n2 = :i * :g
    end
    @test d ≅ @combine(gd, :n1 = :i, :n2 = :i .* :g)

    d = @combine gd @byrow begin
        :n1 = :i * :g
        :n2 = :i * :g
    end
    @test d ≅ @combine(gd, :n1 = :i .* :g, :n2 = :i .* :g)
end

@testset "@by with @byrow" begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :combine, missing])

    @test @by(df, :g, @byrow :n = :i + :g) ≅ @by(df, :g, :n = :i + :g)
    @test @by(df, :g, @byrow :n = :t * string(:y)) ≅ @by(df, :g, :n = :t .* string.(:y))
    @test @by(df, :g, @byrow :n = :g == 1 ? 100 : 500) ≅ @by(df, :g, :n = ifelse.(:g .== 1, 100, 500))
    @test @by(df, :g, @byrow :n = :g == 1 && :t == "a") ≅ @by(df, :g, :n = map((g, t) -> g == 1 && t == "a", :g, :t))
    @test @by(df, :g, @byrow :n = first(:g)) ≅ @by(df, :g, :n = first.(:g))

    d = @by df :g @byrow begin
        :n1 = :i
        :n2 = :i * :g
    end
    @test d ≅ @by(df, :g, :n1 = :i, :n2 = :i .* :g)
    @test d ≅ @by(df, :g, @byrow(:n1 = :i), @byrow(:n2 = :i * :g))

    d = @by df :g @byrow begin
        $:n1 = :i
        :n2 = $"i" * :g
    end
    @test d ≅ @by(df, :g, :n1 = :i, :n2 = :i .* :g)

    d = @by df :g @byrow begin
        :n1 = $"i"
        $:n2 = :i * :g
    end
    @test d ≅ @by(df, :g, :n1 = :i, :n2 = :i .* :g)

    d = @by df :g @byrow begin
        :n1 = begin
            :i
        end
        :n2 = :i * :g
    end
    @test d ≅ @by(df, :g, :n1 = :i, :n2 = :i .* :g)

    d = @by df :g @byrow begin
        :n1 = :i * :g
        :n2 = :i * :g
    end
    @test d ≅ @by(df, :g, :n1 = :i .* :g, :n2 = :i .* :g)
end

@testset "multi-argument r*, #280" begin
    df = DataFrame(a = [1, 2, 3], b = [4, 5, 6])
    d = @rtransform(df, :x = :a, :y = :a + :b)
    @test d.x == df.a

    d = @rtransform!(copy(df), :x = :a, :y = :a + :b)
    @test d.x == df.a

    d = @rselect(df, :x = :a, :y = :a + :b)
    @test d.x == df.a

    d = @rselect!(copy(df), :x = :a, :y = :a + :b)
    @test d.x == df.a

    d = @rsubset(df, :a > 1, :b > 3)
    @test d == DataFrame(a = [2, 3], b = [5, 6])

    d = @rsubset!(copy(df), :a > 1, :b > 3)
    @test d == DataFrame(a = [2, 3], b = [5, 6])

    d = @rorderby(df, -:a, -:b)
    @test d == DataFrame(a = [3, 2, 1], b = [6, 5, 4])

    # no `@orderby!` for some reason.
end

@testset "row operations modify" begin
    df = DataFrame(a = 1, b = 2)
    cols = ["a", "b"]
    for c in cols
        @rtransform! df $c = $c + 100
    end

    @test df == DataFrame(a = 101, b = 102)

    df = DataFrame(a = 1, b = 2)
    @rselect! df :x = :a
    @test names(df) == ["x"]
end

end


text:
module OptimWrapper

using JavaCall

const AbstractConvergenceChecker = @jimport org.hipparchus.optim.AbstractConvergenceChecker
const AbstractOptimizationProblem = @jimport org.hipparchus.optim.AbstractOptimizationProblem
const BaseMultiStartMultivariateOptimizer = @jimport org.hipparchus.optim.BaseMultiStartMultivariateOptimizer
const BaseMultivariateOptimizer = @jimport org.hipparchus.optim.BaseMultivariateOptimizer
const BaseOptimizer = @jimport org.hipparchus.optim.BaseOptimizer
const Class = @jimport java.lang.Class
const ConvergenceChecker = @jimport org.hipparchus.optim.ConvergenceChecker
const Enum = @jimport java.lang.Enum
const Incrementor = @jimport org.hipparchus.util.Incrementor
const InitialGuess = @jimport org.hipparchus.optim.InitialGuess
const JString = @jimport java.lang.String
const Locale = @jimport java.util.Locale
const LocalizedOptimFormats = @jimport org.hipparchus.optim.LocalizedOptimFormats
const MaxEval = @jimport org.hipparchus.optim.MaxEval
const MaxIter = @jimport org.hipparchus.optim.MaxIter
const Object = @jimport java.lang.Object
const OptimizationData = @jimport org.hipparchus.optim.OptimizationData
const OptimizationProblem = @jimport org.hipparchus.optim.OptimizationProblem
const Optional = @jimport java.util.Optional
const Pair = @jimport org.hipparchus.util.Pair
const PointValuePair = @jimport org.hipparchus.optim.PointValuePair
const PointVectorValuePair = @jimport org.hipparchus.optim.PointVectorValuePair
const RandomVectorGenerator = @jimport org.hipparchus.random.RandomVectorGenerator
const SimpleBounds = @jimport org.hipparchus.optim.SimpleBounds
const SimplePointChecker = @jimport org.hipparchus.optim.SimplePointChecker
const SimpleValueChecker = @jimport org.hipparchus.optim.SimpleValueChecker
const SimpleVectorValueChecker = @jimport org.hipparchus.optim.SimpleVectorValueChecker

export AbstractConvergenceChecker
export AbstractOptimizationProblem
export BaseMultiStartMultivariateOptimizer
export BaseMultivariateOptimizer
export BaseOptimizer
export ConvergenceChecker
export InitialGuess
export LocalizedOptimFormats
export MaxEval
export MaxIter
export OptimizationData
export OptimizationProblem
export PointValuePair
export PointVectorValuePair
export SimpleBounds
export SimplePointChecker
export SimpleValueChecker
export SimpleVectorValueChecker
export compare_to
export converged
export create
export describe_constable
export equals
export get_absolute_threshold
export get_class
export get_convergence_checker
export get_declaring_class
export get_evaluation_counter
export get_evaluations
export get_first
export get_initial_guess
export get_iteration_counter
export get_iterations
export get_key
export get_localized_string
export get_lower
export get_lower_bound
export get_max_eval
export get_max_evaluations
export get_max_iter
export get_max_iterations
export get_optima
export get_point
export get_point_ref
export get_relative_threshold
export get_second
export get_source_string
export get_start_point
export get_upper
export get_upper_bound
export get_value
export get_value_ref
export hash_code
export name
export notify
export notify_all
export optimize
export ordinal
export to_string
export unbounded
export unlimited
export value_of
export values
export wait

include("LinearWrapper/LinearWrapper.jl")
include("NonlinearWrapper/NonlinearWrapper.jl")
include("UnivariateWrapper/UnivariateWrapper.jl")
include("abstract_convergence_checker.jl")
include("abstract_optimization_problem.jl")
include("base_multi_start_multivariate_optimizer.jl")
include("base_multivariate_optimizer.jl")
include("base_optimizer.jl")
include("convergence_checker.jl")
include("initial_guess.jl")
include("localized_optim_formats.jl")
include("max_eval.jl")
include("max_iter.jl")
include("optimization_data.jl")
include("optimization_problem.jl")
include("point_value_pair.jl")
include("point_vector_value_pair.jl")
include("simple_bounds.jl")
include("simple_point_checker.jl")
include("simple_value_checker.jl")
include("simple_vector_value_checker.jl")

end


text:
module Alexya

using Reexport

@reexport using Colors, Luxor

import Base: convert, length, getindex, get
import Gtk, Luxor.Cairo

import Gtk: GtkScale,
    GtkButton,
    GtkLabel,
    GtkEntry,
    GtkSpinButton,
    GtkExpander,
    GtkSpinner,
    GtkTextView,
    GtkStatusbar,
    GtkSwitch,
    GtkLinkButton,
    GtkFontButton,
    GtkAppChooser,
    GtkColorButton,
    GtkProgressBar,
    GtkFileChooser,
    GtkCheckButton,
    GtkToggleButton,
    GtkComboBoxText,
    GtkVolumeButton,
    GtkCanvas,
    GtkBox,
    GtkGrid,
    GtkFrame,
    GtkPaned,
    GtkWindow,
    GtkMenuBar,
    GtkToolbar,
    GtkTreeView,
    GtkNotebook,
    GtkScrolledWindow,
    CssProvider,
    GtkWidget,
    GdkRGBA,
    GConstants,
    GtkAlign,
    set_gtk_property!,
    get_gtk_property,
    signal_connect,
    signal_handler_disconnect,
    bytestring,
    show, showall,
    destroy,
    GAccessor,
    width, height

include("exports.jl")
include("utils.jl")
include("vector.jl")
include("math.jl")
include("gtk.jl")
include("widgets.jl")
include("grid.jl")
include("events.jl")
include("app.jl")
include("layout.jl")
include("svg.jl")
include("sprites.jl")
include("controls.jl")

function runexamples()
    include(joinpath(@__DIR__, "..", "examples", "runexamples.jl"))
end

export runexamples, GtkAlign

if VERSION >= v"1.1"   # work around https://github.com/JuliaLang/julia/issues/34121
    include("precompile.jl")
    _precompile_()
end

end


text:
export  normview

"""
```julia
# Example

# The resolvent of a random Hermitian matrix is approximately diagonal,
# see Section 3 in https://arxiv.org/pdf/1903.10060.pdf
 normview(resolvent(randHermitian(1000,norm=true))(0.5+0.1im))
```
"""
function normview(M::AbstractMatrix)
    M = norm.(M)
    heatmap(M,
            yflip=true,
                c=cgrad(:binary, 50, categorical = true),
                ratio = 1,
                xlims=(0.5,size(M)[1]+0.5),ylims=(0.5,size(M)[2]+0.5)
                )
end


text:
macro inmodule(mod, expr)
  expr = Expr(:quote, expr)
  esc(:(eval($mod, $expr)))
end

Requires.@init @inmodule Base begin

  if VERSION < v"0.4-"
    # Display primitives

    export HTML, Text

    type HTML{T}
      content::T
    end

    function HTML(xs...)
      HTML() do io
        for x in xs
          writemime(io, MIME"text/html"(), x)
        end
      end
    end

    writemime(io::IO, ::MIME"text/html", h::HTML) = print(io, h.content)
    writemime(io::IO, ::MIME"text/html", h::HTML{Function}) = h.content(io)

    type Text{T}
      content::T
    end

    print(io::IO, t::Text) = print(io, t.content)
    print(io::IO, t::Text{Function}) = t.content(io)
    writemime(io::IO, ::MIME"text/plain", t::Text) = print(io, t)
  end

  # Add Julia to the path

  @unix_only begin
    export addtopath!, rmfrompath!

    function addtopath!(target = "/usr/local/bin/julia")
      source = joinpath(JULIA_HOME, "julia")
      isfile(target) && error("There is already a file at $target. Please call rmfrompath!() first.")
      run(`sudo ln -s $source $target`)
    end

    function rmfrompath!(target = "/usr/local/bin/julia")
      isfile(target) && run(`sudo rm $target`)
    end
  end

  # Profiler

  export profile
  profile() = Main.Jewel.ProfileView.fetch()

end


text:
#!/usr/bin/env julia

using DelimitedFiles: writedlm
using LinearAlgebra: Diagonal

candidates(mat, idx) = @views setdiff(BitSet(1:size(mat[], 1)),
                                      BitSet(mat[][idx[1], :]), BitSet(mat[][:, idx[2]]))

function recursion2(mat, K)
    i = findfirst(==(0), mat[])
    i === nothing && return true

    for c in candidates(mat, i)
        mat[][i] = c
        recursion2(mat, K) && return true
        mat[][i] = 0
    end

    false
end

function recursion(fakediag, K)
    trace = fakediag[] |> sum
    trace > K && return false, nothing

    if trace == K
        mat = fakediag[] |> Diagonal |> Matrix |> Ref
        recursion2(mat, K) && return true, mat
    end

    i = findfirst(==(0), fakediag[])
    i === nothing && return false, nothing

    for ii in 1:min(length(fakediag[]), (K - trace - length(fakediag[]) + i))
        fakediag[][i] = ii
        res, mat = recursion(fakediag, K)
        res && return true, mat
    end

    fakediag[][i] = 0
    return false, nothing
end

function process(case, N, K)
    fakediag = Ref(zeros(Int, N))

    res, mat = recursion(fakediag, K)
    if res
        println("Case #$(case): POSSIBLE")
        writedlm(stdout, mat[], ' ')
    else
        println("Case #$(case): IMPOSSIBLE")
    end
end

function main()
    T = parse(Int, readline(stdin))
    for case in 1:T
        N, K = parse.(Int, stdin |> readline |> split)
        process(case, N, K)
    end
end

isinteractive() || main()

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
typealias Key Vector{UInt8}
typealias Value Vector{UInt8}


text:
#--------------------------------------------------------------------
# Spacetime Discretization methods in Julia
# Soham 01-2018
#--------------------------------------------------------------------

function sconv(bnd1::Function, bnd2::Function, frhs::Function, asol::Function,
                                                                   Nx::Int, Ny::Int, M::Int)::Float64
    dbase = fdistribute(bnd1, bnd2, frhs, Nx, Ny, M)
    sumL2 = 0.0
    for m in 1:M, n in 1:M
        sPatch = interpolatePatch(dbase[[m,n]], 4Nx, 4Ny).value
        fPatch = Float64[asol(x,y) for x in chebgrid(4Nx, M, m), y in chebgrid(4Ny, M, n)]
        sumL2 += (L2norm(fPatch, sPatch, chebweights(4Nx)/M, chebweights(4Ny)/M))^2
    end
    return sqrt(sumL2)
end

function showconv(bnd1::Function, bnd2::Function, frhs::Function, analyticsol::Function,
                                                               maxmodes::Int, maxlevels::Int, h::Int)
    if maxlevels > 1
        L2 = 1
        p  = maxmodes
        println("----------------------------------------------------------------------------------")
        println("==> h-convergence")
        println("----------------------------------------------------------------------------------")
        for level in 0:maxlevels
            L2error = sconv(bnd1, bnd2, frhs, analyticsol, p, p, h^level)
            level == 0 ? R = 0 : R =  L2/L2error
            @printf("p = %s | np = %s | h^(p+1) = %s | R = %e | L2 error = %e \n",
                    lpad(p,4," "), lpad(h^level,4," "),  lpad(h^(p+1),4," "), R, L2error)
            L2 = L2error
        end
    else
        println("----------------------------------------------------------------------------------")
        println("==> p-convergence")
        println("----------------------------------------------------------------------------------")
        L2 = 1
        h  = 1
        level = 0
        for p in 1:2:maxmodes
            L2error = sconv(bnd1, bnd2, frhs, analyticsol, p, p, h^level)
            p == 1 ? R = 0 : R =  L2/L2error
            @printf("p = %s | np = %s | h^(p+1) = %s | R = %e | L2 error = %e \n",
                    lpad(p,4," "), lpad(h^level,4," "),  lpad(h^(p+1),4," "), R, L2error)
            L2 = L2error
        end
    end
end


text:
#####
##### Reductions of AbstractField
#####

struct Reduction{R, O, D}
    reduce! :: R
    operand :: O
    dims :: D
end

"""
    Reduction(reduce!, operand; dims)

Return a `Reduction` of `operand` with `reduce!`, along `dims`. Note that `Reduction`
expects `reduce!` to operate in-place.

Example
=======

```jldoctest
using Oceananigans

Nx, Ny, Nz = 3, 3, 3

grid = RectilinearGrid(size=(Nx, Ny, Nz), x=(0, 1), y=(0, 1), z=(0, 1),
                       topology=(Periodic, Periodic, Periodic))

c = CenterField(grid)

set!(c, (x, y, z) -> x + y + z)

max_c² = Field(Reduction(maximum!, c^2, dims=3))

compute!(max_c²)

max_c²[1:Nx, 1:Ny]

# output
3×3 Matrix{Float64}:
 1.36111  2.25     3.36111
 2.25     3.36111  4.69444
 3.36111  4.69444  6.25
```
"""
Reduction(reduce!, operand; dims) = Reduction(reduce!, operand, dims)

location(r::Reduction) = reduced_location(location(r.operand); dims=r.dims)

function Field(reduction::Reduction;
               data = nothing,
               indices = indices(reduction.operand),
               recompute_safely = false)

    operand = reduction.operand
    grid = operand.grid
    LX, LY, LZ = loc = location(reduction)
    indices = reduced_indices(indices; dims=reduction.dims)

    if isnothing(data)
        data = new_data(grid, loc, indices)
        recompute_safely = false
    end

    boundary_conditions = FieldBoundaryConditions(grid, loc, indices)
    status = recompute_safely ? nothing : FieldStatus()

    return Field(loc, grid, data, boundary_conditions, indices, reduction, status)
end

const ReducedComputedField = Field{<:Any, <:Any, <:Any, <:Reduction}

function compute!(field::ReducedComputedField, time=nothing)
    reduction = field.operand
    compute_at!(reduction.operand, time)
    reduction.reduce!(field, reduction.operand)
    return field
end

#####
##### show
#####

Base.show(io::IO, field::ReducedComputedField) =
    print(io, "$(summary(field))\n",
          "├── data: $(typeof(field.data)), size: $(size(field))\n",
          "├── grid: $(summary(field.grid))\n",
          "├── operand: $(summary(field.operand))\n",
          "└── status: $(summary(field.status))")

Base.summary(r::Reduction) = string(r.reduce!,
                                    " over dims ", r.dims,
                                    " of ", summary(r.operand))

Base.show(io::IO, r::Reduction) =
    print(io, "$(summary(r))\n",
          "└── operand: $(summary(r.operand))\n",
          "    └── grid: $(summary(r.operand.grid))")


text:
# copied and pasted directly from
# https://github.com/ChrisRackauckas/VectorizedRoutines.jl/blob/master/src/matlab.jl
#
# using VectorizedRoutines.jl directly causes Pkg versioning issues with SpecialFunctions.jl
# Not sure why, so I just put the meshgrid code directly in here.

"""
meshgrid(vx)
Computes an (x,y)-grid from the vectors (vx,vx).
For more information, see the MATLAB documentation.

Copied and pasted directly from [VectorizedRoutines.jl](https://github.com/ChrisRackauckas/VectorizedRoutines.jl/blob/master/src/matlab.jl).
Using VectorizedRoutines.jl directly causes Pkg versioning issues with SpecialFunctions.jl
"""
meshgrid(v::AbstractVector) = meshgrid(v, v)

"""
meshgrid(vx,vy)
Computes an (x,y)-grid from the vectors (vx,vy).
For more information, see the MATLAB documentation.

Copied and pasted directly from [VectorizedRoutines.jl](https://github.com/ChrisRackauckas/VectorizedRoutines.jl/blob/master/src/matlab.jl).
Using VectorizedRoutines.jl directly causes Pkg versioning issues with SpecialFunctions.jl
"""
function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where {T}
    m, n = length(vy), length(vx)
    return repeat(reshape(vx, 1, n),m,1), repeat(reshape(vy, m, 1),1,n)
end

"""
meshgrid(vx,vy,vz)
Computes an (x,y,z)-grid from the vectors (vx,vy,vz).
For more information, see the MATLAB documentation.

Copied and pasted directly from [VectorizedRoutines.jl](https://github.com/ChrisRackauckas/VectorizedRoutines.jl/blob/master/src/matlab.jl).
Using VectorizedRoutines.jl directly causes Pkg versioning issues with SpecialFunctions.jl
"""
function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T},
                     vz::AbstractVector{T}) where {T}
    m, n, o = length(vy), length(vx), length(vz)
    vx = reshape(vx, 1, n, 1)
    vy = reshape(vy, m, 1, 1)
    vz = reshape(vz, 1, 1, o)
    om = ones(Int, m)
    on = ones(Int, n)
    oo = ones(Int, o)
    (vx[om, :, oo], vy[:, on, oo], vz[om, on, :])
end


text:
###############################################################################
##
## Asses convergence
##
##############################################################################

function assess_convergence(δx,
                            x,
                            maxabs_gr,
                            ssr,
                            trial_ssr,
                            xtol::Real,
                            ftol::Real,
                            grtol::Real)

    x_converged, f_converged, gr_converged = false, false, false
    maxabs_x = maxabs(x)
    if abs(trial_ssr - ssr) <= ftol * (abs(ssr) + ftol)
        f_converged = true
    elseif maxabs(δx) <= xtol
            x_converged = true
    elseif maxabs_gr <= grtol
            gr_converged = true
    end
    converged = x_converged || f_converged || gr_converged
    return x_converged, f_converged, gr_converged, converged
end

###############################################################################
##
## Finite Exception
##
##############################################################################

type IsFiniteException <: Exception
  indices::Vector{Int}
end
Base.show(io::IO, e::IsFiniteException) = print(io,
  "During the resolution of the non-linear system, the evaluation" *
  " of the following equation(s) resulted in a non-finite number: $(e.indices)")

function check_isfinite(x::Vector)
    i = find(!isfinite(x))
    if !isempty(i)
        throw(IsFiniteException(i))
    end
end

function check_isfinite(x)
end

###############################################################################
##
## Trace
##
##############################################################################

immutable OptimizationState
    iteration::Int
    value::Float64
    g_norm::Float64
end

immutable OptimizationTrace
    states::Vector{OptimizationState}
end

OptimizationTrace() = OptimizationTrace(Vector{OptimizationState}[])
function update!(tr::OptimizationTrace,
                 iteration::Integer,
                 f_x::Real,
                 grnorm::Real,
                 store_trace::Bool,
                 show_trace::Bool,
                 show_every::Int = 1)
    os = OptimizationState(iteration, f_x, grnorm)
    if store_trace
        push!(tr, os)
    end
    if show_trace
        if iteration % show_every == 0
            show(os)
        end
    end
    return
end

function Base.show(io::IO, t::OptimizationTrace)
    @printf io "Iter     Function value   Gradient norm \n"
    @printf io "------   --------------   --------------\n"
    for state in t.states
        show(io, state)
    end
    return
end

function Base.show(io::IO, t::OptimizationState)
    @printf io "%6d   %14e   %14e\n" t.iteration t.value t.g_norm
    return
end

Base.push!(t::OptimizationTrace, s::OptimizationState) = push!(t.states, s)
Base.getindex(t::OptimizationTrace, i::Integer) = getindex(t.states, i)
##############################################################################
##
## Methods for sparse and dense matrices
##
##############################################################################

for (name, symbol) in ((:Ac_mul_B!, 'T'),
                       (:A_mul_B!, 'N'))
    @eval begin
        $name(α::Number, A::StridedVecOrMat, x::AbstractVector, β, y::AbstractVector) = BLAS.gemm!($symbol, 'N', convert(eltype(y), α), A, x, convert(eltype(y), β), y)
    end
end

function colsumabs2!(v::AbstractVector, A::StridedVecOrMat)
    length(v) == size(A, 2) || error("v should have length size(A, 2)")
    @inbounds for j in 1:length(v)
        v[j] = sumabs2(view(A, :, j))
    end
end

function colsumabs2!(v::AbstractVector, A::SparseMatrixCSC)
    length(v) == size(A, 2) || error("v should have length size(A, 2)")
    @inbounds for j in 1:length(v)
        v[j] = sumabs2(view(nonzeros(A), nzrange(A, j)))
    end
end

_zeros(x) = fill!(similar(x), 0)

function wdot(x::AbstractVector, y::AbstractVector, w::AbstractVector)
    (length(x) == length(y) && length(y) == length(w)) || error("vectors have not the same length")
    out = zero(one(eltype(x)) * one(eltype(y)) * one(eltype(w)))
    @inbounds for i in 1:length(x)
        out += w[i] * x[i] * y[i]
    end
    return out
end

# can be user written
wdot(x, y, w) = dot(x, y, w)
wsumabs2(x, w) = wdot(x, x, w)
wnorm(x, w) = sqrt(wsumabs2(x, w))


text:
"""
`InitialGuessODE`: Initial guess for ordinary differential equations

### Fields

* `int`: interpolation structure
* `v`:   vector  field
* `Δt`:  time step
* `s`:   number of extrapolation stages (for initialisation)
"""
mutable struct InitialGuessODE{DT, TT, VT, IT <: Interpolator}
    int::IT
    v::VT
    Δt::TT
    s::Int

    function InitialGuessODE{DT,TT,VT,IT}(interp, v, Δt) where {DT,TT,VT,IT}
        new(interp, v, Δt, get_config(:ig_extrapolation_stages))
    end
end

function InitialGuessODE(interp, equation::ODE{DT,TT,VT}, Δt::TT) where {DT,TT,VT}
    int = interp(zero(DT), one(DT), Δt, ndims(equation))
    InitialGuessODE{DT, TT, VT, interp}(int, equation.v, Δt)
end

function InitialGuessODE(interp, equation::IODE{DT,TT,ΘT,FT,GT,HT,VT}, Δt::TT) where {DT,TT,ΘT,FT,GT,HT,VT}
    int = interp(zero(DT), one(DT), Δt, ndims(equation))
    InitialGuessODE{DT, TT, VT, interp}(int, equation.v, Δt)
end

function InitialGuessODE(interp, equation::VODE{DT,TT,AT,FT,GT,VT}, Δt::TT) where {DT,TT,AT,FT,GT,VT}
    int = interp(zero(DT), one(DT), Δt, ndims(equation))
    InitialGuessODE{DT, TT, VT, interp}(int, equation.v, Δt)
end

"Initialise initial guess, i.e., given t₀, t₁, q₁ compute q₀, v₀, v₁."
function initialize!(ig::InitialGuessODE{DT,TT,VT,IT},
                t₁::TT,
                q₁::SolutionVector{DT},
                v₁::SolutionVector{DT},
                t₀::TT,
                q₀::SolutionVector{DT},
                v₀::SolutionVector{DT}) where {DT,TT,VT,IT}
    midpoint_extrapolation(ig.v, t₁, t₀, q₁, q₀, ig.s)
    ig.v(t₀, q₀, v₀)
    ig.v(t₁, q₁, v₁)
end

"compute vector field of new solution"
function update!(ig::InitialGuessODE{DT,TT}, t₁::TT, q₁::SolutionVector{DT}, v₁::Vector{DT}) where {DT,TT}
    ig.v(t₁, q₁, v₁)
end

function CommonFunctions.evaluate!(ig::InitialGuessODE{DT,TT},
                q₀::SolutionVector{DT},
                v₀::SolutionVector{DT},
                q₁::SolutionVector{DT},
                v₁::SolutionVector{DT},
                guess::SolutionVector{DT},
                c::TT=one(TT)) where {DT,TT}

    if q₀ == q₁
        @warn "q₀ and q₁ in initial guess are identical! Setting q=q₁."
        guess .= q₁
    else
        evaluate!(ig.int, q₀, q₁, v₀, v₁, one(TT)+c, guess)
    end
end

function CommonFunctions.evaluate!(ig::InitialGuessODE{DT,TT},
                q₀::SolutionVector{DT},
                v₀::SolutionVector{DT},
                q₁::SolutionVector{DT},
                v₁::SolutionVector{DT},
                guess_q::SolutionVector{DT},
                guess_v::SolutionVector{DT},
                c::TT=one(TT)) where {DT,TT}
                @assert length(guess_q) == length(guess_v)

    if q₀ == q₁
        @warn "q₀ and q₁ in initial guess are identical! Setting q=q₁ and v=0."
        guess_q .= q₁
        guess_v .= 0
    else
        evaluate!(ig.int, q₀, q₁, v₀, v₁, one(TT)+c, guess_q, guess_v)
    end
end


text:
using JSON: JSON

# add_format(format"JSON", "JSON", [".json"])

function load_truss_json(filepath::AbstractString)
    open(filepath) do io
        load_truss_json(io)
    end
end

function load_truss_json(io::IO) #io::Stream{format"JSON"})
    data = JSON.parse(io)
    ndim = data["dimension"]
    n = data["node_num"]
    m = data["element_num"]
    iT = Int
    T = Float64

    node_points = Dict{iT,SVector{ndim,T}}()
    for (i, ndata) in enumerate(data["nodes"])
        node_points[i] = convert(SVector{ndim,T}, ndata["point"])
        if "node_ind" in keys(ndata)
            @assert 1 + ndata["node_ind"] == i
        end
    end
    @assert length(node_points) == n

    elements = Dict{iT,Tuple{iT,iT}}()
    element_inds_from_tag = Dict()
    for (i, edata) in enumerate(data["elements"])
        elements[i] = (edata["end_node_inds"]...,) .+ 1
        if "elem_ind" in keys(edata)
            @assert 1 + edata["elem_ind"] == i
        end
        elem_tag = edata["elem_tag"]
        if elem_tag ∉ keys(element_inds_from_tag)
            element_inds_from_tag[elem_tag] = []
        end
        push!(element_inds_from_tag[elem_tag], i)
    end
    @assert length(elements) == m

    E_from_tag = Dict()
    ν_from_tag = Dict()
    for mat in data["materials"]
        # if elem_tag list has length 0, use tag `nothing`
        mat["elem_tags"] = length(mat["elem_tags"]) == 0 ? [nothing] : mat["elem_tags"]
        for e_tag in mat["elem_tags"]
            if e_tag in keys(E_from_tag)
                @warn "Multiple materials assigned to the same element tag |$(e_tag)|!"
            end
            E_from_tag[e_tag] = T(mat["E"])
            if "mu" ∈ keys(mat)
                ν_from_tag[e_tag] = T(mat["mu"])
            else
                ν_from_tag[e_tag] = 0.0
            end
        end
    end
    A_from_tag = Dict()
    for cs in data["cross_secs"]
        # if elem_tag list has length 0, use tag `nothing`
        cs["elem_tags"] = length(cs["elem_tags"]) == 0 ? [nothing] : cs["elem_tags"]
        for e_tag in cs["elem_tags"]
            if e_tag in keys(A_from_tag)
                @warn "Multiple cross secs assigned to the same element tag |$(e_tag)|!"
            end
            A_from_tag[e_tag] = T(cs["A"])
        end
    end
    mats = Array{TrussFEAMaterial}(undef, m)
    crosssecs = Array{TrussFEACrossSec}(undef, m)
    for (tag, e_ids) in element_inds_from_tag
        # @show tag, e_ids
        for ei in e_ids
            if !(tag ∈ keys(A_from_tag))
                # use default material (key `nothing`)
                A = A_from_tag[nothing]
            else
                A = A_from_tag[tag]
            end
            crosssecs[ei] = TrussFEACrossSec(A)

            if !(tag ∈ keys(E_from_tag))
                # use default material (key `nothing`)
                E = E_from_tag[nothing]
            else
                E = E_from_tag[tag]
            end
            if !(tag ∈ keys(ν_from_tag))
                # use default material (key `nothing`)
                ν = ν_from_tag[nothing]
            else
                ν = ν_from_tag[tag]
            end
            mats[ei] = TrussFEAMaterial(E, ν)
        end
    end

    # TODO only translation dof for now
    @assert(length(data["supports"]) > 0)
    fixities = Dict{iT,SVector{ndim,Bool}}()
    for sdata in data["supports"]
        supp_v = iT(sdata["node_ind"]) + 1
        fixities[supp_v] = sdata["condition"][1:ndim]
    end

    load_cases = Dict()
    for (lc_ind, lc_data) in data["loadcases"]
        nploads = length(lc_data["ploads"])
        @assert nploads > 0
        ploads = Dict{iT,SVector{ndim,T}}()
        for pl in lc_data["ploads"]
            load_v = pl["node_ind"] + 1
            ploads[load_v] = convert(SVector{ndim,T}, pl["force"])
        end
        load_cases[lc_ind] = ploads
    end

    return node_points, elements, mats, crosssecs, fixities, load_cases
end