text:
using Base.Test
using MIPVerify
using MIPVerify: BlurringPerturbationFamily
isdefined(:TestHelpers) || include("../../../TestHelpers.jl")
using TestHelpers: batch_test_adversarial_example

@testset "FC + Softmax" begin

### Parameters for neural net
batch = 1
in1_height = 8
in1_width = 8
in1_channels = 2

A_height = 5
A_width = batch*in1_height*in1_width*in1_channels

B_height = 3
B_width = A_height

srand(5)
input_size = (batch, in1_height, in1_width, in1_channels)
x0 = rand(input_size)

kernelf1 = rand(-10:10, A_width, A_height)
biasf1 = rand(-10:10, A_height)

kernelf2 = rand(-10:10, B_width, B_height)
biasf2 = rand(-10:10, B_height)

nn = Sequential([Flatten(4),
        Linear(kernelf1, biasf1), ReLU(),
        Linear(kernelf2, biasf2)
    ],
    "tests.integration.generated_weights.fc+softmax"
)

@testset "Blurring perturbations accept multiple input channels." begin
    pp_blur = BlurringPerturbationFamily((5, 5))

    expected_objective_values = Dict(
        (1, pp_blur, Inf, 0) => 0,
        (2, pp_blur, Inf, 0) => 0.1936444023676,
        (3, pp_blur, Inf, 0) => 0.4151904104813
    )

    batch_test_adversarial_example(nn, x0, expected_objective_values)
end

end


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

"""
    g1 ∩ g2

Return the intersection of two geometries `g1` and `g2` as a new geometry.
"""
Base.intersect(g1::Geometry, g2::Geometry) = get(intersection(g1, g2))

"""
    intersection([f], g1, g2)

Compute the intersection of two geometries `g1` and `g2`
and apply function `f` to it. Default function is [`identity`](@ref).

## Examples

```julia
intersecttype(g1, g2) do I
  if I isa CrossingLines
    # do something
  else
    # do nothing
  end
end
```

### Notes

When a custom function `f` is used that reduces the number of
return types, Julia is able to optimize the branches of the code
and generate specialized code. This is not the case when
`f === identity`.
"""
intersection(f, g1, g2) = intersection(f, g2, g1)
intersection(g1, g2)    = intersection(identity, g1, g2)

"""
    IntersectionType

The different types of intersection that may occur between geometries.
Type `IntersectionType` in a Julia session to see the full list.
"""
@enum IntersectionType begin
  # no intersection
  NoIntersection

  # segment-segment intersection
  CrossingSegments
  MidTouchingSegments
  CornerTouchingSegments
  OverlappingSegments

  # box-box intersection
  OverlappingBoxes
  FaceTouchingBoxes
  CornerTouchingBoxes

  # line-line intersection
  CrossingLines
  OverlappingLines

  # ray-box intersection
  CrossingRayBox
  TouchingRayBox

  # segment-triangle intersection
  IntersectingSegmentTriangle

  # ray-triangle intersection
  IntersectingRayTriangle

  # segment-plane intersection
  CrossingSegmentPlane
  TouchingSegmentPlane
  OverlappingSegmentPlane
end

"""
    Intersection{G}

An intersection between geometries holding a geometry of type `G`.
"""
struct Intersection{GeometryType}
  type::IntersectionType
  geom::GeometryType
end

Intersection(type, geom) = Intersection{typeof(geom)}(type, geom)

"""
    type(intersection)

Return the type of intersection computed between geometries.
"""
type(I::Intersection) = I.type

"""
    get(intersection)

Return the underlying geometry stored in the intersection object.
"""
Base.get(I::Intersection) = I.geom

# helper macro for developers in case we decide to
# change the internal representation of Intersection
macro IT(type, geom, func)
  type = esc(type)
  geom = esc(geom)
  func = esc(func)
  :(Intersection($type, $geom) |> $func)
end

# ----------------
# IMPLEMENTATIONS
# ----------------

include("intersections/lines.jl")
include("intersections/segments.jl")
include("intersections/boxes.jl")
include("intersections/raybox.jl")
include("intersections/segmentplane.jl")
include("intersections/segmenttriangle.jl")
include("intersections/raytriangle.jl")
include("intersections/geompolygon.jl")

# ------------------------
# MISSING IMPLEMENTATIONS
# ------------------------

intersection(f, ::Ball, ::Polygon) = throw(ErrorException("not implemented"))

# -----------------
# TRUE/FALSE CHECK
# -----------------

"""
    hasintersect(g1, g2)

Return `true` if geometries `g1` and `g2` intersect and `false` otherwise.

The algorithm works with any geometry that has a well-defined [`supportfun`](@ref).

## References

* Gilbert, E., Johnson, D., Keerthi, S. 1988. [A fast
  Procedure for Computing the Distance Between Complex
  Objects in Three-Dimensional Space]
  (https://ieeexplore.ieee.org/document/2083)
"""
function hasintersect(g1::Geometry{Dim,T}, g2::Geometry{Dim,T}) where {Dim,T}
  # handle non-convex geometries
  if !isconvex(g1)
    d1 = discretize(g1, FIST())
    return hasintersect(d1, g2)
  elseif !isconvex(g2)
    d2 = discretize(g2, FIST())
    return hasintersect(g1, d2)
  end

  # initial direction
  c1, c2 = centroid(g1), centroid(g2)
  d = c1 == c2 ? rand(Vec{Dim,T}) : c2 - c1

  # first point in Minkowski difference
  P = minkowskipoint(g1, g2, d)

  # origin of coordinate system
  O = minkowskiorigin(Dim, T)

  # initialize simplex vertices
  points = [P]

  # move towards the origin
  d = O - P
  while true
    P = minkowskipoint(g1, g2, d)
    if (P - O) ⋅ d < zero(T)
      return false
    end
    push!(points, P)

    # line segment case
    if length(points) == 2
      B, A = points
      AB = B - A
      AO = O - A
      d  = AB × AO × AB
    else # simplex case
      C, B, A = points
      AB  = B - A
      AC  = C - A
      AO  = O - A
      ABᵀ = AC × AB × AB
      ACᵀ = AB × AC × AC
      if ABᵀ ⋅ AO > zero(T)
        popat!(points, 1) # pop C
        d = ABᵀ
      elseif ACᵀ ⋅ AO > zero(T)
        popat!(points, 2) # pop B
        d = ACᵀ
      else
        return true
      end
    end
  end
end

hasintersect(d1::Domain, g2::Geometry) =
  any(g1 -> hasintersect(g1, g2), d1)

hasintersect(g1::Geometry, d2::Domain) =
  hasintersect(d2, g1)

function hasintersect(d1::Domain, d2::Domain)
  for g1 in d1, g2 in d2
    hasintersect(g1, g2) && return true
  end
  return false
end

# support point in Minkowski difference
function minkowskipoint(g1::Geometry{Dim,T}, g2::Geometry{Dim,T}, d) where {Dim,T}
  n = Vec{Dim}(d[1:Dim])
  v = supportfun(g1, n) - supportfun(g2, -n)
  Point(ntuple(i-> i ≤ Dim ? v[i] : zero(T), max(Dim, 3)))
end

# origin of coordinate system
minkowskiorigin(Dim, T) = Point(ntuple(i->zero(T), max(Dim, 3)))

ext: jl


text:
###################################################################################################
#                                        Load Packages
###################################################################################################
cd(@__DIR__)
using Pkg
Pkg.activate("../..")
using Revise, ACTRSimulators
import_gui()
import ACTRSimulators: start!, press_key!, repaint!
include("../PVT_Example/PVT.jl")
include("../PVT_Example/PVT_Model.jl")
###################################################################################################
#                                        Run Model
###################################################################################################
scheduler = ACTRScheduler(;model_trace=true)
task = PVT(;scheduler, n_trials=4, visible=true, realtime=true)
n_models = 2
models = map(id -> init_model(scheduler, task, id), 1:n_models)
@time run!(models, task)


text:
"""
```julia
function Not(; name)
```

NOT gate in 9-level logic.

# Connectors
- `x`
  Input [`DigitalPin`](@ref)
- `y`
  Output [`DigitalPin`](@ref)
"""
function Not(; name)
    @named x = DigitalPin()
    @named y = DigitalPin()

    eqs = [
        x.i ~ y.i
        y.val ~ _not(x.val)
    ]
    ODESystem(eqs, t, [], [], systems=[x, y], name=name)
end

"""
```julia
function And(; name, N=2)
```

AND gate in 9-level logic, with `N` inputs

# Connectors
- `x1`, `x2`, ...
  `N` input [`DigitalPin`](@ref)s
- `y`
  Output [`DigitalPin`](@ref)
"""
function And(; name, N=2)
    x = map(1:N) do i
        DigitalPin(name=Symbol(:x, i))
    end
    @named y = DigitalPin()

    vals = [k.val for k in x]
    eqs = [
        y.val ~ _and(vals...)
        y.i ~ sum(k -> k.i, x)
    ]
    ODESystem(eqs, t, [], [], systems=[x..., y], name=name)
end

"""
```julia
function Nand(; name, N=2)
```

NAND gate in 9-level logic, with `N` inputs

# Connectors
- `x1`, `x2`, ...
  `N` input [`DigitalPin`](@ref)s
- `y`
  Output [`DigitalPin`](@ref)
"""
function Nand(; name, N=2)
    x = map(1:N) do i
        DigitalPin(name=Symbol(:x, i))
    end
    @named y = DigitalPin()

    vlist = [k.val for k in x]
    eqs = [
        y.val ~ _not(_and(vlist...))
        y.i ~ sum(k -> k.i, x)
    ]
    ODESystem(eqs, t, [], [], systems=[x..., y], name=name)
end

"""
```julia
function Or(; name, N=2)
```

OR gate in 9-level logic, with `N` inputs

# Connectors
- `x1`, `x2`, ...
  `N` input [`DigitalPin`](@ref)s
- `y`
  Output [`DigitalPin`](@ref)
"""
function Or(; name, N=2)
    x = map(1:N) do i
        DigitalPin(name=Symbol(:x, i))
    end
    @named y = DigitalPin()

    vals = [k.val for k in x]
    eqs = [
        y.val ~ _or(vals...)
        y.i ~ sum(k -> k.i, x)
    ]
    ODESystem(eqs, t, [], [], systems=[x..., y], name=name)
end

"""
```julia
function Nor(; name, N=2)
```

NOR gate in 9-level logic, with `N` inputs

# Connectors
- `x1`, `x2`, ...
  `N` input [`DigitalPin`](@ref)s
- `y`
  Output [`DigitalPin`](@ref)
"""
function Nor(; name, N=2)
    x = map(1:N) do i
        DigitalPin(name=Symbol(:x, i))
    end
    @named y = DigitalPin()

    vlist = [k.val for k in x]
    eqs = [
        y.val ~ _not(_or(vlist...))
        y.i ~ sum(k -> k.i, x)
    ]
    ODESystem(eqs, t, [], [], systems=[x..., y], name=name)
end

"""
```julia
function Xor(; name, N=2)
```

XOR gate in 9-level logic, with `N` inputs

# Connectors
- `x1`, `x2`, ...
  `N` input [`DigitalPin`](@ref)s
- `y`
  Output [`DigitalPin`](@ref)
"""
function Xor(; name, N=2)
    x = map(1:N) do i
        DigitalPin(name=Symbol(:x, i))
    end
    @named y = DigitalPin()

    vals = [k.val for k in x]
    eqs = [
        y.val ~ _xor(vals...)
        y.i ~ sum(k -> k.i, x)
    ]
    ODESystem(eqs, t, [], [], systems=[x..., y], name=name)
end

"""
```julia
function Xnor(; name, N=2)
```

XNOR gate in 9-level logic, with `N` inputs

# Connectors
- `x1`, `x2`, ...
  `N` input [`DigitalPin`](@ref)s
- `y`
  Output [`DigitalPin`](@ref)
"""
function Xnor(; name, N=2)
    x = map(1:N) do i
        DigitalPin(name=Symbol(:x, i))
    end
    @named y = DigitalPin()

    vlist = [k.val for k in x]
    eqs = [
        y.val ~ _not(_xor(vlist...))
        y.i ~ sum(k -> k.i, x)
    ]
    ODESystem(eqs, t, [], [], systems=[x..., y], name=name)
end


text:
################################
######    Polymorphism    ######
################################
# Expected number of polymorphism above frequency x from the standard diffusion theory
# f(x) = ∫(s) θs* (1/x(1-x)) * ( (ℯ^(Ns)*(1-ℯ^(-4Ns(1-x))) ) / (ℯ^(4Ns-1)))
# Convolote with the binomial to obtain the downsampled sfs
# E[P(x)] = ∑_x=x*→x=1 fB(x*)
############Neutral#############

"""

	DiscSFSNeutDown()

Expected rate of neutral allele frequency reduce by backgrou	nd selection. The spectrum depends on the number of individual []

```math
\\mathbb{E}[Ps_{(x)}] = \\sum{x^{*}=x}{x^{*}=1}f_{B}(x)
```

# Return:
 - `Array{Float64}`: expected rate of neutral alleles frequencies.
"""
function DiscSFSNeutDown(param::parameters,binom::SparseMatrixCSC{Float64,Int64})

	NN2 = convert(Int64,ceil(param.NN*param.B))

	# Allocating variables
	neutral_sfs(i::Int64) = 1.0/(i)

	x = collect(0:NN2)
	solved_neutral_sfs = neutral_sfs.(x)
	replace!(solved_neutral_sfs, Inf => 0.0)

	# subsetDict = get(param.bn,param.B,1)
	# subsetDict = binom
	out::Array{Float64,1} = param.B * (param.θ_coding) * 0.25 * (binom*solved_neutral_sfs)

	return out
end

############Positive############
# Variable gamma in function changed to s to avoid problem with exported SpecialFunctions.gamma
"""

	DiscSFSSelPosDown(s,p)

Expected rate of positive selected allele frequency reduce by background selection. The spectrum depends on the number of individuals.

# Arguments
 - `s::Int64`: selection strength.
 - `p::Float64`: positive selected alleles probabilty.
# Return:
 - `Array{Float64}`: expected positive selected alleles frequencies.
"""
function DiscSFSSelPosDown(param::parameters,s::Int64,p::Float64,binom::SparseMatrixCSC{Float64,Int64})

	if p == 0.0
		out = zeros(Float64,param.nn + 1)
		out = out[2:end-1]
	else

		# Solving sfs
		red_plus = Φ(param,s)
		NN2      = convert(Int64,ceil(param.NN*param.B))
		xa1      = collect(0:NN2)
		xa2      = xa1/(NN2)

		# Solving float precision performance using exponential rule. Only one BigFloat estimation.
		s_corrected = s*param.B

		s_exp1 = exp(s_corrected*2)
		s_exp2 = exp(s_corrected*-2)

		positiveSfs(i::Float64,g1::Float64=s_exp1,g2::Float64=s_exp2,p::Float64=p) = Float64(p*0.5*(g1*(1- g2^(1.0-i))/((g1-1.0)*i*(1.0-i))))

		# Original
		# p*0.5*(ℯ^(2*s_corrected)*(1-ℯ^(-2.0*s_corrected*(1.0-i)))/((ℯ^(2*s_corrected)-1.0)*i*(1.0-i)))

		# Allocating outputs
		solved_positive_sfs::Array{Float64,1} = (1.0/(NN2)) * (positiveSfs.(xa2))
		replace!(solved_positive_sfs, NaN => 0.0)

		# subsetDict = get(param.bn,param.B,1)
		# out               = (param.θ_coding)*red_plus*0.75*(subsetDict*solved_positive_sfs)
		out::Array{Float64,1} = param.θ_coding * red_plus * 0.75 * (binom*solved_positive_sfs)
		#=out::Array{Float64,1} = param.θᵣ .* red_plus .* 0.75 .* (binom*solved_positive_sfs)=#

	end

	return out
end

function DiscSFSSelPosDownArb(param::parameters,s::Int64,p::Float64,binom::SparseMatrixCSC{Float64,Int64})

	if p == 0.0
		out = zeros(Float64,param.nn + 1)
		out = out[2:end-1]
	else

		red_plus = Φ(param,s)

		# Solving sfs
		NN2  = convert(Int64,ceil(param.NN*param.B))
		xa1  = collect(0:NN2)
		xa2  = xa1/(NN2)

		# Solving float precision performance using exponential rule. Only one BigFloat estimation.
		s_corrected = s*param.B
		s_exp1::Quadmath.Float128 = exp(Quadmath.Float128(s_corrected*2))
		s_exp2::Quadmath.Float128 = exp(Quadmath.Float128(s_corrected*-2))

		positiveSfs(i::Float64,g1::Quadmath.Float128=s_exp1,g2::Quadmath.Float128=s_exp2,p::Float64=p) = Float64(p*0.5*(g1*(1- g2^(1.0-i))/((g1-1.0)*i*(1.0-i))))
		# Allocating outputs
		solved_positive_sfs::Array{Float64,1} = (1.0/(NN2)) * (positiveSfs.(xa2))
		replace!(solved_positive_sfs, NaN => 0.0)
		out::Array{Float64,1} = param.θ_coding * red_plus * 0.75 * (binom*solved_positive_sfs)
		#=out::Array{Float64,1} = param.θᵣ .* red_plus .* 0.75 .* (binom*solved_positive_sfs)=#

	end

	return out
end

######Slightly deleterious######
"""

	DiscSFSSelNegDown(param,p)

Expected rate of positive selected allele frequency reduce by background selection. Spectrum drawn on a gamma DFE. It depends on the number of individuals.

# Arguments
 - `p::Float64`: positive selected alleles probabilty.
# Return:
 - `Array{Float64}`: expected negative selected alleles frequencies.
"""
function DiscSFSSelNegDown(param::parameters,p::Float64,binom::SparseMatrixCSC{Float64,Int64})

	# subsetDict = get(param.bn,param.B,1)
	solved_negative = DiscSFSSelNeg(param,p)
	out = param.B .* (param.θ_coding) .* 0.75 .* (binom*solved_negative)

	return out
end

function DiscSFSSelNeg(param::parameters,p::Float64)

	beta     = param.scale/(1.0*param.B)
	NN2      = convert(Int64, ceil(param.NN*param.B))
	xa       = collect(0:NN2)/NN2

	solve_z   = similar(xa)

	z(x::Float64,p::Float64=p) = (1.0-p)*(2.0^-param.shape)*(beta^param.shape)*(-zeta(param.shape,x+beta/2.0) + zeta(param.shape,(2+beta)/2.0))/((-1.0+x)*x)

	solve_z   = z.(xa)

	if (solve_z[1] == Inf || isnan(solve_z[1]))
		solve_z[1] = 0.0
	end
	if (solve_z[lastindex(solve_z)] == Inf || isnan(solve_z[lastindex(solve_z)]))
		solve_z[lastindex(solve_z)] = 0.0
	end

	return 1.0/(NN2+0.0).*solve_z
end

"""
	cumulative_sfs(sfs_tmp)

Changing SFS considering all values above a frequency *x*. The original asymptotic-MK approach takes Pn(x) and Ps(x) as the number of polymorphic sites at frequency *x* rather than above *x*, but this approach scales poorly as sample size increases. We define the polymorphic spectrum as stated above since these quantities trivially have the same asymptote but are less affected by changing sample size.
"""
function cumulative_sfs(sfs_tmp::Array,freqs::Bool=true)

	out      = Array{Float64}(undef, size(sfs_tmp,1),size(sfs_tmp,2))

	if freqs
		idx = 2
	else
		idx = 1
	end

	out[1,idx:end] = sum(sfs_tmp[:,idx:end],dims=1)

	@simd for i in 2:(size(sfs_tmp)[1])

		#=app = view(out,i-1,:) .- view(sfs_tmp,i-1,:)=#
		app = out[i-1,idx:end] .- sfs_tmp[i-1,idx:end]

		if sum(app) > 0.0
			out[i,idx:end] = app
		else
			out[i,idx:end] = zeros(length(app))
		end
	end

	if freqs
		out[:,1] = sfs_tmp[:,1]
	end

	return out
end

"""
	reduce_sfs(sfs_tmp,bins)

Function to bin the SFS into a sample of N individuals.
"""
function reduce_sfs(sfs_tmp::Array,bins::Int64)

	n   = (bins*2) - 1
	f   = sfs_tmp[:,1]
	sfs = sfs_tmp[:,2:end]

	b    = collect(1/n:1/n:1)
	inds = searchsortedfirst.(Ref(b), f)
	out  = zeros((n,size(sfs_tmp,2)))
	out[:,1] = unique(inds)
	sfs_grouped = hcat(inds,sfs)

	for i in unique(inds)
		out[out[:,1] .== i,2:end] = sum(sfs_grouped[sfs_grouped[:,1] .== i,2:end],dims=1)
	end
	out[:,1] = b

	return(out)
end


text:
# NOTE in mcmc.jl

using ForwardDiff
using ForwardDiff: value, Dual

CTAG = ForwardDiff.Tag{Val{:custom_tag}, Float64}

#TODO need to pause langevin updates for a moment
"""
    findProposalLaw(xx, tt, P˟, P̃, Ls, Σs; dt=1/5000, timeChange=true,
                    solver::ST=Ralston3())

Initialise the object with proposal law and all the necessary containers needed
for the simulation of the guided proposals
"""
function findProposalLaw(::LangevinUpdt, ::Type{K}, xx, tt, P˟, P̃, Ls, Σs, τ; dt=1/5000,
                         solver::ST=Ralston3(),
                         changePt::ODEChangePt=NoChangePt()) where {K,ST}
    m = length(xx) - 1
    P = Array{ContinuousTimeProcess,1}(undef,m)
    params(Pˣ)
#    P˟_D = clone(Pˣ, )#TODO here is stop point
    for i in m:-1:1
        numPts = Int64(ceil((tt[i+1]-tt[i])/dt))+1
        t = τ(tt[i], tt[i+1]).( range(tt[i], stop=tt[i+1], length=numPts) )
        xx_D = Dual{CT}.(xx[i+1])
        L_D = Dual{CT}.(Ls[i])
        Σ_D = Dual{CT}.(Σs[i])

        P[i] = ( (i==m) ? GuidPropBridge(K, t, P˟, P̃[i], L_D, xx_D, Σ_D;
                                         changePt=changePt, solver=ST()) :
                          GuidPropBridge(K, t, P˟, P̃[i], L_D, xx_D, Σ_D,
                                         P[i+1].H[1], P[i+1].Hν[1], P[i+1].c[1];
                                         changePt=changePt, solver=ST()) )
    end
    P
end

function prepareLangevin(𝓦𝓢::Workspace{ObsScheme}, θ, ::UpdtIdx, y, m, yPr,
                         priors, ::ST, uidx) where {ObsScheme,UpdtIdx,ST}
    idxToUpdt = idx(UpdtIdx())
    function _ll(ϑ)
        XX, WW, P, fpt = 𝓦𝓢.XX, 𝓦𝓢.WW, 𝓦𝓢.P, 𝓦𝓢.fpt
        updateLaws!(P, ϑ)
        solveBackRec!(NoBlocking(), P, ST()) # changes nothing, but needed for ∇
        findPathFromWiener!(XX, y, WW, P, 1:m)

        ll = logpdf(yPr, y)
        ll += pathLogLikhd(ObsScheme(), XX, P, 1:m, fpt)
        ll += lobslikelihood(P[1], y)
        for prior in priors
            ll += logpdf(prior, θ)
        end
        ll
        ϑ[1]
    end
    ϑ = [θ[i] for i in idxToUpdt]
    chunkSize = 1
    result = 𝓦𝓢.result[uidx]
    cfg = ForwardDiff.GradientConfig(_ll, ϑ, ForwardDiff.Chunk{chunkSize}())
    ForwardDiff.gradient!(result, _ll, ϑ, cfg)
    DiffResults.value(result), DiffResults.gradient(result)
end

function postProcessLangevin(𝓦𝓢::Workspace{ObsScheme}, θᵒ, ::UpdtIdx, y, m, yPr,
                         priors, ::ST, uidx) where {ObsScheme,UpdtIdx,ST}
    idxToUpdt = idx(UpdtIdx())
    function _ll(ϑ)
        XXᵒ, WW, Pᵒ, fpt = 𝓦𝓢.XXᵒ, 𝓦𝓢.WW, 𝓦𝓢.Pᵒ, 𝓦𝓢.fpt
        for (i, ui) in enumerate(idxToUpdt)
            θᵒ[ui] = ϑ[i]
        end
        updateLaws!(Pᵒ, θᵒ)
        solveBackRec!(NoBlocking(), Pᵒ, ST()) # changes nothing, but needed for ∇
        findPathFromWiener!(XX, y, WW, P, 1:m)

        yPrᵒ = invStartPt(y, yPr, Pᵒ[1])

        ll = logpdf(yPrᵒ, y)
        ll += pathLogLikhd(ObsScheme(), XXᵒ, Pᵒ, 1:m, fpt)
        ll += lobslikelihood(Pᵒ[1], y)
        for prior in priors
            ll += logpdf(prior, θᵒ)
        end
        ll
        ϑ[2]
    end
    ϑ = [θᵒ[i] for i in idxToUpdt]
    chunkSize = 1
    result = 𝓦𝓢.resultᵒ[uidx]
    cfg = ForwardDiff.GradientConfig(_ll, ϑ, ForwardDiff.Chunk{chunkSize}())
    ForwardDiff.gradient!(result, _ll, ϑ, cfg)

    yPrᵒ = invStartPt(y, yPr, 𝓦𝓢.Pᵒ[1])
    DiffResults.value(result), DiffResults.gradient(result), yPrᵒ
end

function updateParam!(::LangevinUpdt, 𝔅::NoBlocking, tKern, θ,
                      ::UpdtIdx, yPr, 𝓦𝓢::Workspace{ObsScheme}, ll, priors,
                      recomputeODEs; solver::ST=Ralston3(), verbose=false,
                      it=NaN, uidx=NaN) where {ObsScheme <: AbstractObsScheme, ST, UpdtIdx}
    WW, Pᵒ, P, XXᵒ, XX, fpt = 𝓦𝓢.WW, 𝓦𝓢.Pᵒ, 𝓦𝓢.P, 𝓦𝓢.XXᵒ, 𝓦𝓢.XX, 𝓦𝓢.fpt
    m = length(WW)
    y = XX[1].yy[1]
    ll, ∇ll = prepareLangevin(𝓦𝓢, θ, UpdtIdx(), y, m, yPr, priors, ST(), uidx) # TODO pre-allocate ∇ll
    θᵒ = rand(tKern, θ, ∇ll, UpdtIdx())               # sample new parameter
    llᵒ, ∇llᵒ, yPrᵒ = postProcessLangevin(𝓦𝓢, θᵒ, UpdtIdx(), y, m, yPr, priors, ST(), uidx)

    printInfo(verbose, it, ll, llᵒ)

    llr = ( llᵒ - ll + logpdf(tKern, θᵒ, θ, ∇llᵒ, UpdtIdx()) - logpdf(tKern, θ, θᵒ, ∇ll, UpdtIdx()))

    # Accept / reject
    if acceptSample(llr, verbose)
        swap!(XX, XXᵒ, P, Pᵒ, 1:m)
        swap!(𝓦𝓢.resultᵒ, 𝓦𝓢.result, 1:m)
        return llᵒ, true, θᵒ, yPrᵒ
    else
        return ll, false, θ, yPr
    end
end


text:
isassign(frame) = isassign(frame, frame.pc)
isassign(frame, pc) = (pc in frame.framecode.used)

lookup_var(frame, val::SSAValue) = frame.framedata.ssavalues[val.id]
lookup_var(frame, ref::GlobalRef) = getfield(ref.mod, ref.name)
function lookup_var(frame, slot::SlotNumber)
    val = frame.framedata.locals[slot.id]
    val !== nothing && return val.value
    throw(UndefVarError(frame.framecode.src.slotnames[slot.id]))
end

function lookup_expr(frame, e::Expr)
    head = e.head
    head === :the_exception && return frame.framedata.last_exception[]
    if head === :static_parameter
        arg = e.args[1]::Int
        if isassigned(frame.framedata.sparams, arg)
            return frame.framedata.sparams[arg]
        else
            syms = sparam_syms(frame.framecode.scope)
            throw(UndefVarError(syms[arg]))
        end
    end
    head === :boundscheck && length(e.args) == 0 && return true
    error("invalid lookup expr ", e)
end

"""
    rhs = @lookup(frame, node)
    rhs = @lookup(mod, frame, node)

This macro looks up previously-computed values referenced as SSAValues, SlotNumbers,
GlobalRefs, QuoteNode, sparam or exception reference expression.
It will also lookup symbols in `moduleof(frame)`; this can be supplied ahead-of-time via
the 3-argument version.
If none of the above apply, the value of `node` will be returned.
"""
macro lookup(args...)
    length(args) == 2 || length(args) == 3 || error("invalid number of arguments ", length(args))
    havemod = length(args) == 3
    local mod
    if havemod
        mod, frame, node = args
    else
        frame, node = args
    end
    nodetmp = gensym(:node)  # used to hoist, e.g., args[4]
    if havemod
        fallback = quote
            isa($nodetmp, Symbol) ? getfield($(esc(mod)), $nodetmp) :
            $nodetmp
        end
    else
        fallback = quote
            $nodetmp
        end
    end
    quote
        $nodetmp = $(esc(node))
        isa($nodetmp, SSAValue) ? lookup_var($(esc(frame)), $nodetmp) :
        isa($nodetmp, GlobalRef) ? lookup_var($(esc(frame)), $nodetmp) :
        isa($nodetmp, SlotNumber) ? lookup_var($(esc(frame)), $nodetmp) :
        isa($nodetmp, QuoteNode) ? $nodetmp.value :
        isa($nodetmp, Symbol) ? getfield(moduleof($(esc(frame))), $nodetmp) :
        isa($nodetmp, Expr) ? lookup_expr($(esc(frame)), $nodetmp) :
        $fallback
    end
end

# This is used only for new struct/abstract/primitive nodes.
# The most important issue is that in these expressions, :call Exprs can be nested,
# and hence our re-use of the `callargs` field of Frame would introduce
# bugs. Since these nodes use a very limited repertoire of calls, we can special-case
# this quite easily.
function lookup_or_eval(@nospecialize(recurse), frame, @nospecialize(node))
    if isa(node, SSAValue)
        return lookup_var(frame, node)
    elseif isa(node, SlotNumber)
        return lookup_var(frame, node)
    elseif isa(node, Symbol)
        return getfield(moduleof(frame), node)
    elseif isa(node, QuoteNode)
        return node.value
    elseif isa(node, Expr)
        ex = Expr(node.head)
        for arg in node.args
            push!(ex.args, lookup_or_eval(recurse, frame, arg))
        end
        if ex.head === :call
            f = ex.args[1]
            if f === Core.svec
                return Core.svec(ex.args[2:end]...)
            elseif f === Core.apply_type
                return Core.apply_type(ex.args[2:end]...)
            elseif f === Core.typeof
                return Core.typeof(ex.args[2])
            elseif f === Base.getproperty
                return Base.getproperty(ex.args[2], ex.args[3])
            else
                error("unknown call f ", f)
            end
        else
            dump(ex)
            error("unknown expr ", ex)
        end
    elseif isa(node, Int) || isa(node, Number)   # Number is slow, requires subtyping
        return node
    elseif isa(node, Type)
        return node
    end
    return eval_rhs(recurse, frame, node)
end

function resolvefc(frame, @nospecialize(expr))
    if isa(expr, SlotNumber)
        expr = lookup_var(frame, expr)
    elseif isa(expr, SSAValue)
        expr = lookup_var(frame, expr)
        isa(expr, Symbol) && return QuoteNode(expr)
    end
    (isa(expr, Symbol) || isa(expr, String) || isa(expr, Ptr) || isa(expr, QuoteNode)) && return expr
    isa(expr, Tuple{Symbol,Symbol}) && return expr
    isa(expr, Tuple{String,String}) && return expr
    isa(expr, Tuple{Symbol,String}) && return expr
    isa(expr, Tuple{String,Symbol}) && return expr
    if isexpr(expr, :call)
        a = (expr::Expr).args[1]
        (isa(a, QuoteNode) && a.value === Core.tuple) || error("unexpected ccall to ", expr)
        return Expr(:call, GlobalRef(Core, :tuple), (expr::Expr).args[2:end]...)
    end
    error("unexpected ccall to ", expr)
end

function collect_args(frame::Frame, call_expr::Expr; isfc::Bool=false)
    args = frame.framedata.callargs
    resize!(args, length(call_expr.args))
    mod = moduleof(frame)
    args[1] = isfc ? resolvefc(frame, call_expr.args[1]) : @lookup(mod, frame, call_expr.args[1])
    for i = 2:length(args)
        args[i] = @lookup(mod, frame, call_expr.args[i])
    end
    return args
end

"""
    ret = evaluate_foreigncall(frame::Frame, call_expr)

Evaluate a `:foreigncall` (from a `ccall`) statement `callexpr` in the context of `frame`.
"""
function evaluate_foreigncall(frame::Frame, call_expr::Expr)
    head = call_expr.head
    args = collect_args(frame, call_expr; isfc = head==:foreigncall)
    for i = 2:length(args)
        arg = args[i]
        args[i] = isa(arg, Symbol) ? QuoteNode(arg) : arg
    end
    head === :cfunction && (args[2] = QuoteNode(args[2]))
    scope = frame.framecode.scope
    data = frame.framedata
    if !isempty(data.sparams) && scope isa Method
        sig = scope.sig
        args[2] = instantiate_type_in_env(args[2], sig, data.sparams)
        args[3] = Core.svec(map(args[3]) do arg
            instantiate_type_in_env(arg, sig, data.sparams)
        end...)
    end
    return Core.eval(moduleof(frame), Expr(head, args...))
end

# We have to intercept ccalls / llvmcalls before we try it as a builtin
function bypass_builtins(frame, call_expr, pc)
    if isassigned(frame.framecode.methodtables, pc)
        tme = frame.framecode.methodtables[pc]
        if isa(tme, Compiled)
            fargs = collect_args(frame, call_expr)
            f = to_function(fargs[1])
            fmod = parentmodule(f)
            if fmod === JuliaInterpreter.CompiledCalls || fmod === Core.Compiler
                return Some{Any}(Base.invokelatest(f, fargs[2:end]...))
            else
                return Some{Any}(f(fargs[2:end]...))
            end
        end
    end
    return nothing
end

function evaluate_call_compiled!(::Compiled, frame::Frame, call_expr::Expr; enter_generated::Bool=false)
    # @assert !enter_generated
    pc = frame.pc
    ret = bypass_builtins(frame, call_expr, pc)
    isa(ret, Some{Any}) && return ret.value
    ret = maybe_evaluate_builtin(frame, call_expr, false)
    isa(ret, Some{Any}) && return ret.value
    fargs = collect_args(frame, call_expr)
    f = fargs[1]
    popfirst!(fargs)  # now it's really just `args`
    return f(fargs...)
end

function evaluate_call_recurse!(@nospecialize(recurse), frame::Frame, call_expr::Expr; enter_generated::Bool=false)
    pc = frame.pc
    ret = bypass_builtins(frame, call_expr, pc)
    isa(ret, Some{Any}) && return ret.value
    ret = maybe_evaluate_builtin(frame, call_expr, true)
    isa(ret, Some{Any}) && return ret.value
    call_expr = ret
    fargs = collect_args(frame, call_expr)
    if fargs[1] === Core.eval
        return Core.eval(fargs[2], fargs[3])  # not a builtin, but worth treating specially
    elseif fargs[1] === Base.rethrow
        err = length(fargs) > 1 ? fargs[2] : frame.framedata.last_exception[]
        throw(err)
    end
    if fargs[1] === Core.invoke # invoke needs special handling
        f_invoked = which(fargs[2], fargs[3])
        fargs_pruned = [fargs[2]; fargs[4:end]]
        sig = Tuple{_Typeof.(fargs_pruned)...}
        ret = prepare_framecode(f_invoked, sig; enter_generated=enter_generated)
        isa(ret, Compiled) && invoke(fargs[2:end]...)
        framecode, lenv = ret
        lenv === nothing && return framecode  # this was a Builtin
        fargs = fargs_pruned
    else
        framecode, lenv = get_call_framecode(fargs, frame.framecode, frame.pc; enter_generated=enter_generated)
        if lenv === nothing
            if isa(framecode, Compiled)
                f = popfirst!(fargs)  # now it's really just `args`
                return Base.invokelatest(f, fargs...)
            end
            return framecode  # this was a Builtin
        end
    end
    newframe = prepare_frame_caller(frame, framecode, fargs, lenv)
    npc = newframe.pc
    shouldbreak(newframe, npc) && return BreakpointRef(newframe.framecode, npc)
    # if the following errors, handle_err will pop the stack and recycle newframe
    if recurse === finish_and_return!
        # Optimize this case to avoid dynamic dispatch
        ret = finish_and_return!(finish_and_return!, newframe, false)
    else
        ret = recurse(recurse, newframe, false)
    end
    isa(ret, BreakpointRef) && return ret
    frame.callee = nothing
    return_from(newframe)
    return ret
end

"""
    ret = evaluate_call!(Compiled(), frame::Frame, call_expr)
    ret = evaluate_call!(recurse,    frame::Frame, call_expr)

Evaluate a `:call` expression `call_expr` in the context of `frame`.
The first causes it to be executed using Julia's normal dispatch (compiled code),
whereas the second recurses in via the interpreter.
`recurse` has a default value of [`JuliaInterpreter.finish_and_return!`](@ref).
"""
evaluate_call!(::Compiled, frame::Frame, call_expr::Expr; kwargs...) = evaluate_call_compiled!(Compiled(), frame, call_expr; kwargs...)
evaluate_call!(@nospecialize(recurse), frame::Frame, call_expr::Expr; kwargs...) = evaluate_call_recurse!(recurse, frame, call_expr; kwargs...)
evaluate_call!(frame::Frame, call_expr::Expr; kwargs...) = evaluate_call!(finish_and_return!, frame, call_expr; kwargs...)

# The following come up only when evaluating toplevel code
function evaluate_methoddef(frame, node)
    f = node.args[1]
    if isa(f, Symbol)
        mod = moduleof(frame)
        f = isdefined(mod, f) ? getfield(mod, f) : Core.eval(moduleof(frame), Expr(:function, f))  # create a new function
    end
    length(node.args) == 1 && return f
    sig = @lookup(frame, node.args[2])::SimpleVector
    body = @lookup(frame, node.args[3])
    ccall(:jl_method_def, Cvoid, (Any, Any, Any), sig, body, moduleof(frame))
    return nothing
end

function structname(frame, node)
    name = node.args[1]
    if isa(name, GlobalRef)
        mod = name.module
        name = name.name
    else
        mod = moduleof(frame)
        name = name::Symbol
    end
    return name, mod
end

function set_structtype_const(mod::Module, name::Symbol)
    dt = Base.unwrap_unionall(getfield(mod, name))
    ccall(:jl_set_const, Cvoid, (Any, Any, Any), mod, dt.name.name, dt.name.wrapper)
end

function inplace_lookup!(ex, i, frame)
    a = ex.args[i]
    if isa(a, SSAValue) || isa(a, SlotNumber)
        ex.args[i] = lookup_var(frame, a)
    elseif isexpr(a, :call)
        for j = 1:length(a.args)
            inplace_lookup!(a, j, frame)
        end
    end
    return ex
end

function evaluate_structtype(@nospecialize(recurse), frame, node)
    grsvec!(ex::Expr) = (ex.args[1] = GlobalRef(Core, :svec); return ex)

    name, mod = structname(frame, node)
    supertype = lookup_or_eval(recurse, frame, node.args[4])::Type
    ismutable = node.args[6]
    ninit = node.args[7]
    newstructexpr = Expr(:struct_type, name, nothing, nothing, supertype, nothing, ismutable, ninit)
    for idx in (2, 3, 5)
        ex = newstructexpr.args[idx] = grsvec!(copy(node.args[idx]))
        for i = 2:length(ex.args)
            inplace_lookup!(ex, i, frame)
        end
    end
    Core.eval(mod, newstructexpr)
    VERSION < v"1.2.0-DEV.239" && set_structtype_const(mod, name)
end

function evaluate_abstracttype(@nospecialize(recurse), frame, node)
    name, mod = structname(frame, node)
    params = lookup_or_eval(recurse, frame, node.args[2])::SimpleVector
    supertype = lookup_or_eval(recurse, frame, node.args[3])::Type
    Core.eval(mod, Expr(:abstract_type, name, params, supertype))
    VERSION < v"1.2.0-DEV.239" && set_structtype_const(mod, name)
end

function evaluate_primitivetype(@nospecialize(recurse), frame, node)
    name, mod = structname(frame, node)
    params = lookup_or_eval(recurse, frame, node.args[2])::SimpleVector
    nbits = node.args[3]::Int
    supertype = lookup_or_eval(recurse, frame, node.args[4])::Type
    Core.eval(mod, Expr(:primitive_type, name, params, nbits, supertype))
    VERSION < v"1.2.0-DEV.239" && set_structtype_const(mod, name)
end

function do_assignment!(frame, @nospecialize(lhs), @nospecialize(rhs))
    code, data = frame.framecode, frame.framedata
    if isa(lhs, SSAValue)
        data.ssavalues[lhs.id] = rhs
    elseif isa(lhs, SlotNumber)
        counter = (frame.assignment_counter += 1)
        data.locals[lhs.id] = Some{Any}(rhs)
        data.last_reference[lhs.id] = counter
    elseif isa(lhs, GlobalRef)
        Core.eval(lhs.mod, :($(lhs.name) = $(QuoteNode(rhs))))
    elseif isa(lhs, Symbol)
        Core.eval(moduleof(code), :($lhs = $(QuoteNode(rhs))))
    end
end

function maybe_assign!(frame, @nospecialize(stmt), @nospecialize(val))
    pc = frame.pc
    if isexpr(stmt, :(=))
        lhs = stmt.args[1]
        do_assignment!(frame, lhs, val)
    elseif isassign(frame, pc)
        lhs = SSAValue(pc)
        do_assignment!(frame, lhs, val)
    end
    return nothing
end
maybe_assign!(frame, @nospecialize(val)) = maybe_assign!(frame, pc_expr(frame), val)

function eval_rhs(@nospecialize(recurse), frame, node::Expr)
    head = node.head
    if head === :new
        mod = moduleof(frame)
        args = [@lookup(mod, frame, arg) for arg in node.args]
        T = popfirst!(args)
        rhs = ccall(:jl_new_structv, Any, (Any, Ptr{Any}, UInt32), T, args, length(args))
        return rhs
    elseif head === :splatnew  # Julia 1.2+
        mod = moduleof(frame)
        rhs = ccall(:jl_new_structt, Any, (Any, Any), @lookup(mod, frame, node.args[1]), @lookup(mod, frame, node.args[2]))
        return rhs
    elseif head === :isdefined
        return check_isdefined(frame, node.args[1])
    elseif head === :call
        # here it's crucial to avoid dynamic dispatch
        isa(recurse, Compiled) && return evaluate_call_compiled!(recurse, frame, node)
        return evaluate_call_recurse!(recurse, frame, node)
    elseif head === :foreigncall || head === :cfunction
        return evaluate_foreigncall(frame, node)
    elseif head === :copyast
        val = (node.args[1]::QuoteNode).value
        return isa(val, Expr) ? copy(val) : val
    elseif head === :enter
        return length(frame.framedata.exception_frames)
    elseif head === :boundscheck
        return true
    elseif head === :meta || head === :inbounds || head == (@static VERSION >= v"1.2.0-DEV.462" ? :loopinfo : :simdloop) ||
           head === :gc_preserve_begin || head === :gc_preserve_end
        return nothing
    elseif head === :method && length(node.args) == 1
        return evaluate_methoddef(frame, node)
    end
    return lookup_expr(frame, node)
end

function check_isdefined(frame, @nospecialize(node))
    data = frame.framedata
    if isa(node, SlotNumber)
        return data.locals[node.id] !== nothing
    elseif isexpr(node, :static_parameter)
        return isassigned(data.sparams, node.args[1]::Int)
    elseif isa(node, GlobalRef)
        return isdefined(node.mod, node.name)
    elseif isa(node, Symbol)
        return isdefined(moduleof(frame), node)
    end
    error("unrecognized isdefined node ", node)
end

# For "profiling" where JuliaInterpreter spends its time. See the commented-out block
# in `step_expr!`
const _location = Dict{Tuple{Method,Int},Int}()

function step_expr!(@nospecialize(recurse), frame, @nospecialize(node), istoplevel::Bool)
    pc, code, data = frame.pc, frame.framecode, frame.framedata
    if !is_leaf(frame)
        show_stackloc(frame)
        @show node
    end
    @assert is_leaf(frame)
    local rhs
    # For debugging:
    # show_stackloc(frame)
    # @show node
    # For profiling:
    # location_key = (scopeof(frame), pc)
    # _location[location_key] = get(_location, location_key, 0) + 1
    try
        if isa(node, Expr)
            if node.head === :(=)
                lhs, rhs = node.args
                if isa(rhs, Expr)
                    rhs = eval_rhs(recurse, frame, rhs)
                else
                    rhs = istoplevel ? @lookup(moduleof(frame), frame, rhs) : @lookup(frame, rhs)
                end
                isa(rhs, BreakpointRef) && return rhs
                do_assignment!(frame, lhs, rhs)
            elseif node.head === :gotoifnot
                arg = @lookup(frame, node.args[1])
                if !isa(arg, Bool)
                    throw(TypeError(nameof(frame), "if", Bool, arg))
                end
                if !arg
                    return (frame.pc = node.args[2]::Int)
                end
            elseif node.head === :enter
                rhs = node.args[1]
                push!(data.exception_frames, rhs)
            elseif node.head === :leave
                for _ = 1:node.args[1]
                    pop!(data.exception_frames)
                end
            elseif node.head === :pop_exception
                n = lookup_var(frame, node.args[1])
                deleteat!(data.exception_frames, n+1:length(data.exception_frames))
            elseif node.head === :return
                return nothing
            elseif istoplevel
                if node.head === :method && length(node.args) > 1
                    evaluate_methoddef(frame, node)
                elseif node.head === :struct_type
                    evaluate_structtype(recurse, frame, node)
                elseif node.head === :abstract_type
                    evaluate_abstracttype(recurse, frame, node)
                elseif node.head === :primitive_type
                    evaluate_primitivetype(recurse, frame, node)
                elseif node.head === :module
                    error("this should have been handled by split_expressions")
                elseif node.head === :using || node.head === :import || node.head === :export
                    Core.eval(moduleof(frame), node)
                elseif node.head === :const
                    g = node.args[1]
                    if isa(g, GlobalRef)
                        mod, name = g.module, g.name
                    else
                        mod, name = moduleof(frame), g::Symbol
                    end
                    if VERSION >= v"1.2.0-DEV.239"  # depends on https://github.com/JuliaLang/julia/pull/30893
                        Core.eval(mod, Expr(:const, name))
                    end
                elseif node.head === :thunk
                    newframe = prepare_thunk(moduleof(frame), node)
                    if isa(recurse, Compiled)
                        finish!(recurse, newframe, true)
                    else
                        newframe.caller = frame
                        frame.callee = newframe
                        finish!(recurse, newframe, true)
                        frame.callee = nothing
                    end
                    return_from(newframe)
                elseif node.head === :global
                    Core.eval(moduleof(frame), node)
                elseif node.head === :toplevel
                    mod = moduleof(frame)
                    modexs, _ = split_expressions(mod, node)
                    rhs = Core.eval(mod, Expr(:toplevel,
                        :(for modex in $modexs
                              newframe = ($prepare_thunk)(modex)
                              newframe === nothing && continue
                              while true
                                  ($through_methoddef_or_done!)($recurse, newframe) === nothing && break
                              end
                              $return_from(newframe)
                          end)))
                elseif node.head === :error
                    error("unexpected error statement ", node)
                elseif node.head === :incomplete
                    error("incomplete statement ", node)
                else
                    rhs = eval_rhs(recurse, frame, node)
                end
            elseif node.head === :thunk || node.head === :toplevel
                error("this frame needs to be run at top level")
            else
                rhs = eval_rhs(recurse, frame, node)
            end
        elseif isa(node, GotoNode)
            return (frame.pc = node.label)
        elseif is_GotoIfNot(node)
            node = node::Core.GotoIfNot
            arg = @lookup(frame, node.cond)
            if !isa(arg, Bool)
                throw(TypeError(nameof(frame), "if", Bool, arg))
            end
            if !arg
                return (frame.pc = node.dest)
            end
        elseif is_ReturnNode(node)
            return nothing
        elseif isa(node, NewvarNode)
            # FIXME: undefine the slot?
        elseif istoplevel && isa(node, LineNumberNode)
        elseif istoplevel && isa(node, Symbol)
            rhs = getfield(moduleof(frame), node)
        else
            rhs = @lookup(frame, node)
        end
    catch err
        return handle_err(recurse, frame, err)
    end
    @isdefined(rhs) && isa(rhs, BreakpointRef) && return rhs
    if isassign(frame, pc)
        if !@isdefined(rhs)
            @show frame node
        end
        lhs = SSAValue(pc)
        do_assignment!(frame, lhs, rhs)
    end
    return (frame.pc = pc + 1)
end

"""
    pc = step_expr!(recurse, frame, istoplevel=false)
    pc = step_expr!(frame, istoplevel=false)

Execute the next statement in `frame`. `pc` is the new program counter, or `nothing`
if execution terminates, or a [`BreakpointRef`](@ref) if execution hits a breakpoint.

`recurse` controls call evaluation; `recurse = Compiled()` evaluates :call expressions
by normal dispatch. The default value `recurse = finish_and_return!` will use recursive
interpretation.

If you are evaluating `frame` at module scope you should pass `istoplevel=true`.
"""
step_expr!(@nospecialize(recurse), frame::Frame, istoplevel::Bool=false) =
    step_expr!(recurse, frame, pc_expr(frame), istoplevel)
step_expr!(frame::Frame, istoplevel::Bool=false) =
    step_expr!(finish_and_return!, frame, istoplevel)

"""
    loc = handle_err(recurse, frame, err)

Deal with an error `err` that arose while evaluating `frame`. There are one of three
behaviors:

- if `frame` catches the error, `loc` is the program counter at which to resume
  evaluation of `frame`;
- if `frame` doesn't catch the error, but `break_on_error[]` is `true`,
  `loc` is a `BreakpointRef`;
- otherwise, `err` gets rethrown.
"""
function handle_err(@nospecialize(recurse), frame, err)
    data = frame.framedata
    err_will_be_thrown_to_top_level = isempty(data.exception_frames) && !data.caller_will_catch_err
    if break_on_throw[] || (break_on_error[] && err_will_be_thrown_to_top_level)
        return BreakpointRef(frame.framecode, frame.pc, err)
    end
    if isempty(data.exception_frames)
        if !err_will_be_thrown_to_top_level
            return_from(frame)
        end
        # Check for world age errors, which generally indicate a failure to go back to toplevel
        if isa(err, MethodError)
            is_arg_types = isa(err.args, DataType)
            arg_types = is_arg_types ? err.args : Base.typesof(err.args...)
            if (err.world != typemax(UInt) &&
                hasmethod(err.f, arg_types) &&
                !hasmethod(err.f, arg_types, world = err.world))
                @warn "likely failure to return to toplevel, try JuliaInterpreter.split_expressions"
            end
        end
        rethrow(err)
    end
    data.last_exception[] = err
    return (frame.pc = data.exception_frames[end])
end

if isdefined(Core, :ReturnNode)
    lookup_return(frame, node::Core.ReturnNode) = @lookup(frame, node.val)
else
    lookup_return(frame, node::Expr) = @lookup(frame, node.args[1])
end

"""
    ret = get_return(frame)

Get the return value of `frame`. Throws an error if `frame.pc` does not point to a `return` expression.
`frame` must have already been executed so that the return value has been computed (see,
e.g., [`JuliaInterpreter.finish!`](@ref)).
"""
function get_return(frame)
    node = pc_expr(frame)
    is_return(node) || error("expected return statement, got ", node)
    return lookup_return(frame, node)
end
get_return(t::Tuple{Module,Expr,Frame}) = get_return(t[end])


text:
"""
	ShotProfileWEM(m,d,adj=true;<keyword arguments>)

Shot Profile Wave Equation Migration and Demigration of 3D isotropic data.*

# Arguments
- `m::AbstractString` : filename of image
- `d::AbstractString` : filename of data
- `adj=true` : flag for adjoint (migration), or forward (demigration) (default=true)
- `pspi=true` : flag for Phase Shift Plus Interpolation (default=true)
- `nref=5` : number of reference velocities to use if pspi is selected (default=5)
- `damping=1000.` : damping for deconvolution imaging condition
- `vel="vel"` : seis file containing the velocity (should have same x and z dimensions as the desired image)
- `angx="angx"` : seis file containing incidence angles in the x direction for each shot
- `angy="angy"` : seis file containing incidence angles in the y direction for each shot
- `wav="wav"` : seis file containing the source wavelet (in time domain)
- `sz=0.` : source depth (Dev: read this from source wavelet file for variable source depth)
- `gz=0.` : receiver depth (Dev: read this from data file for variable source depth (but then what to do in fwd op?))
- `nangx=1` : number of angle bins in x direction
- `oangx=0.` : min angle in x direction (angle between source incidence angle and reflector normal in Degrees)
- `dangx=1.` : angle increment in x direction
- `nangy=1` : number of angle bins in y direction
- `oangy=0.` : min angle in y direction (angle between source incidence angle and reflector normal in Degrees)
- `dangy=1.` : angle increment in y direction
- `fmin=0.` : min frequency to process (Hz)
- `fmax=80.` : max frequency to process (Hz)
- `padt=2` : pad factor for the time axis
- `padx=2` : pad factor for the spatial axes
- `verbose=false` : flag for error / debugging messages
- `sx=[0.]` : array of source X positions (meters)
- `sy=[0.]` : array of source Y positions (meters)

*Credits: AS, 2015*

"""
function ShotProfileWEM(m::AbstractString,d::AbstractString,adj=true;pspi=true,nref=5,damping=1000.,vel="vel",angx="angx",angy="angy",wav="wav",sz=0.,gz=0.,nangx=1,oangx=0,dangx=1,nangy=1,oangy=0,dangy=1,fmin=0,fmax=80,padt=2,padx=2,verbose=false,sx=[0],sy=[0])

	nshot = length(sx)
	v,h,e = SeisRead(vel)
	min_imx = h[1].imx
	max_imx = h[end].imx
	nx = max_imx - min_imx + 1
	ox = h[1].mx
	dx = h[2].mx - ox
	min_imy = h[1].imy
	max_imy = h[end].imy
	ny = max_imy - min_imy + 1
	oy = h[1].my
	dy = ny > 1 ? h[nx+1].my - oy : dx
	nz = h[1].n1
	dz = h[1].d1
	oz = h[1].o1
	w,h,e = SeisRead(wav)
	nt = h[1].n1
	dt = h[1].d1
	ot = h[1].o1
	nsx = length(sx)
	osx = sx[1]
	dsx = nsx > 1 ? sx[2] - sx[1] : 1.
	nsy = length(sy)
	osy = sy[1]
	dsy = nsy > 1 ? sy[2] - sy[1] : 1.
	dsx = dsx == 0. ? 1. : dsx
	dsy = dsy == 0. ? 1. : dsy

	shot_list = Array(Shot,nshot)
	for ishot = 1 : nshot
		shot_list[ishot] = Shot(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
		shot_list[ishot].m = join([m "_shot_" Int(floor(sx[ishot])) "_" Int(floor(sy[ishot]))])
		shot_list[ishot].d = join([d "_shot_" Int(floor(sx[ishot])) "_" Int(floor(sy[ishot]))])
		shot_list[ishot].angx = join([angx "_shot_" Int(floor(sx[ishot])) "_" Int(floor(sy[ishot]))])
		shot_list[ishot].angy = join([angy "_shot_" Int(floor(sx[ishot])) "_" Int(floor(sy[ishot]))])
		shot_list[ishot].vel = vel
		shot_list[ishot].wav = wav
		shot_list[ishot].damping = damping
		shot_list[ishot].sx = sx[ishot]
		shot_list[ishot].sy = sy[ishot]
		shot_list[ishot].sz = sz
		shot_list[ishot].gz = gz
		shot_list[ishot].fmin = fmin
		shot_list[ishot].fmax = fmax
		shot_list[ishot].padt = padt
		shot_list[ishot].padx = padx
		shot_list[ishot].adj = (adj == true) ? "y" : "n"
		shot_list[ishot].pspi = (pspi == true) ? "y" : "n"
		shot_list[ishot].nref = nref
		shot_list[ishot].verbose = (verbose == true) ? "y" : "n"
	end

	if (nangx != 1 || nangy != 1)
		for ishot = 1 : nshot
			SeisWindow(angx,shot_list[ishot].angx,key=["sx","sy"],minval=[sx[ishot],sy[ishot]],maxval=[sx[ishot],sy[ishot]])
			SeisWindow(angy,shot_list[ishot].angy,key=["sx","sy"],minval=[sx[ishot],sy[ishot]],maxval=[sx[ishot],sy[ishot]])
		end
	end

	if (adj == true)
		for ishot = 1 : nshot
			SeisWindow(d,shot_list[ishot].d,key=["sx","sy"],minval=[sx[ishot],sy[ishot]],maxval=[sx[ishot],sy[ishot]])
		end
		@sync @parallel for ishot = 1 : nshot
			a = shotwem(shot_list[ishot])
		end
		j = 1
		gather = zeros(Float32,nz,nangx*nangy)
		extent = Seismic.Extent(nz,nx,ny,nangx,nangy,
				oz,ox,oy,oangx,oangy,
				dz,dx,dy,dangx,dangy,
				"Depth","mx","my","angx","angy",
				"m","m","m","degrees","degrees",
				"")
		for imx = 1 : nx
			for imy = 1 : ny
				h = Header[]
				for iangx = 1 : nangx
					for iangy = 1 : nangy
						push!(h,Seismic.InitSeisHeader())
						h[(iangx-1)*nangy + iangy].tracenum = convert(typeof(h[1].tracenum),j + (iangx-1)*nangy + iangy)
						h[(iangx-1)*nangy + iangy].o1 = convert(typeof(h[1].o1),0)
						h[(iangx-1)*nangy + iangy].n1 = convert(typeof(h[1].n1),nz)
						h[(iangx-1)*nangy + iangy].d1 = convert(typeof(h[1].d1),dz)
						h[(iangx-1)*nangy + iangy].imx = convert(typeof(h[1].imx),imx-1 + min_imx)
						h[(iangx-1)*nangy + iangy].imy = convert(typeof(h[1].imy),imy-1 + min_imy)
						h[(iangx-1)*nangy + iangy].mx = convert(typeof(h[1].mx),dx*(imx-1 + min_imx))
						h[(iangx-1)*nangy + iangy].my = convert(typeof(h[1].my),dy*(imy-1 + min_imy))
						h[(iangx-1)*nangy + iangy].iang = convert(typeof(h[1].iang),iangx-1)
						h[(iangx-1)*nangy + iangy].ang = convert(typeof(h[1].ang),(iangx-1)*dangx + oangx)
						h[(iangx-1)*nangy + iangy].iaz = convert(typeof(h[1].iaz),iangy-1)
						h[(iangx-1)*nangy + iangy].az = convert(typeof(h[1].az),(iangy-1)*dangy + oangy)
					end
				end
				SeisWrite(m,gather,h,extent,itrace=j)
				j += nangx*nangy
			end
		end

		m_m = ParseDataName(m)
		m_h = ParseHeaderName(m)
		stream_m = open(m_m,"a+")
		stream_h = open(m_h,"a+")
		for ishot = 1 : nshot

			# println("reading ",shot_list[ishot].m)
			m_shot,h_shot  = SeisRead(shot_list[ishot].m)
			if (nangx != 1 || nangy != 1)
				angx_shot,h_ang,extent = SeisRead(shot_list[ishot].angx)
				angy_shot,h_ang,extent = SeisRead(shot_list[ishot].angy)
			else
				angx_shot = 0 .* m_shot
				angy_shot = 0 .* m_shot
			end
			nx_shot = size(m_shot,2)
			for ix = 1 : nx_shot
				itrace = (h_shot[ix].imx)*ny*nangx*nangy + (h_shot[ix].imy)*nangx*nangy
				position_m = 4*nz*itrace
				seek(stream_m,position_m)
				m = read(stream_m,Float32,nz*nangx*nangy)
				m = reshape(m,convert(Int64,nz),convert(Int64,nangx*nangy))
				for iz = 1 : nz
					# bilinear interpolation of angx and angy onto grid points
					angx = angx_shot[iz,ix]
					iangx = Int(floor((angx - oangx)/dangx)) + 1
					angx1 = (iangx-1)*dangx + oangx
					angx2 = angx1 + dangx
					angy = angy_shot[iz,ix]
					iangy = Int(floor((angy - oangy)/dangy)) + 1
					angy1 = (iangy-1)*dangy + oangy
					angy2 = angy1 + dangy
					w1 = (angx2-angx)*(angy2-angy)/((angx2-angx1)*(angy2-angy1))
					w2 = (angx-angx1)*(angy2-angy)/((angx2-angx1)*(angy2-angy1))
					w3 = (angx2-angx)*(angy-angy1)/((angx2-angx1)*(angy2-angy1))
					w4 = (angx-angx1)*(angy-angy1)/((angx2-angx1)*(angy2-angy1))
					if (iangx >= 1 && iangx <= nangx && iangy >= 1 && iangy <= nangy)
						m[iz,(iangx-1)*nangy + iangy]   += w1*m_shot[iz,ix]
						if (iangx < nangx)
							m[iz,(iangx)*nangy   + iangy]   += w2*m_shot[iz,ix]
						end
						if (nangy > 1)
							m[iz,(iangx-1)*nangy + iangy+1] += w3*m_shot[iz,ix]
							if (iang < nang)
								m[iz,(iangx)*nangy   + iangy+1] += w4*m_shot[iz,ix]
							end
						end
					end
				end
				seek(stream_m,position_m)
				write(stream_m,convert(Array{Float32,1},m[:]))
			end
			SeisRemove(shot_list[ishot].m)
			SeisRemove(shot_list[ishot].d)
			if (nangx != 1 || nangy != 1)
				SeisRemove(shot_list[ishot].angx)
				SeisRemove(shot_list[ishot].angy)
			end

		end
		close(stream_m)
		close(stream_h)

	else

		m_m = ParseDataName(m)
		m_h = ParseHeaderName(m)
		stream_m = open(m_m,"r")
		stream_h = open(m_h,"r")

		for ishot = 1 : nshot
			sx = shot_list[ishot].sx
			sy = shot_list[ishot].sy
			j = 1
			h_shot = Array(Header,nx*ny)
			for imx = 1 : nx
				for imy = 1 : ny
					h_shot[j] = Seismic.InitSeisHeader()
					h_shot[j].tracenum = convert(typeof(h_shot[1].tracenum),j)
					h_shot[j].o1 = convert(typeof(h_shot[1].o1),0)
					h_shot[j].n1 = convert(typeof(h_shot[1].n1),nz)
					h_shot[j].d1 = convert(typeof(h_shot[1].d1),dz)
					h_shot[j].imx = convert(typeof(h_shot[1].imx),imx-1 + min_imx)
					h_shot[j].imy = convert(typeof(h_shot[1].imy),imy-1 + min_imy)
					h_shot[j].mx = convert(typeof(h_shot[1].mx),dx*h_shot[j].imx)
					h_shot[j].my = convert(typeof(h_shot[1].my),dy*h_shot[j].imy)
					j += 1
				end
			end
			m_shot = zeros(nz,nx*ny)
			if (nangx != 1 || nangy != 1)
				angx_shot,h_ang,extent = SeisRead(shot_list[ishot].angx)
				angy_shot,h_ang,extent = SeisRead(shot_list[ishot].angy)
			else
				angx_shot = 0 .* m_shot
				angy_shot  = 0 .* m_shot
			end
			for imx = 1 : nx
				for imy = 1 : ny
					ix = (imx - 1)*ny + imy
					h_shot[ix].sx = convert(typeof(h[1].sx),shot_list[ishot].sx)
					h_shot[ix].sy = convert(typeof(h[1].sy),shot_list[ishot].sy)
					itrace = (imx-1 + min_imx)*ny*nangx*nangy + (imy-1 + min_imy)*nangx*nangy
					position_m = 4*nz*itrace
					seek(stream_m,position_m)
					m = read(stream_m,Float32,nz*nangx*nangy)
					m = reshape(m,convert(Int64,nz),convert(Int64,nangx*nangy))
					for iz = 1 : nz
						# bilinear interpolation of angx and angy onto grid points
						angx = angx_shot[iz,ix]
						iangx = Int(floor((angx - oangx)/dangx)) + 1
						angx1 = (iangx-1)*dangx + oangx
						angx2 = angx1 + dangx
						angy = angy_shot[iz,ix]
						iangy = Int(floor((angy - oangy)/dangy)) + 1
						angy1 = (iangy-1)*dangy + oangy
						angy2 = angy1 + dangy
						w1 = (angx2-angx)*(angy2-angy)/((angx2-angx1)*(angy2-angy1))
						w2 = (angx-angx1)*(angy2-angy)/((angx2-angx1)*(angy2-angy1))
						w3 = (angx2-angx)*(angy-angy1)/((angx2-angx1)*(angy2-angy1))
						w4 = (angx-angx1)*(angy-angy1)/((angx2-angx1)*(angy2-angy1))

						if (iangx >= 1 && iangx <= nangx && iangy >= 1 && iangy <= nangy)
							m_shot[iz,ix] += w1*m[iz,(iangx-1)*nangy + iangy]
							if (iangx < nangx)
								m_shot[iz,ix] += w2*m[iz,(iangx)*nangy   + iangy]
							end
							if (nangy > 1)
								m_shot[iz,ix] += w3*m[iz,(iangx-1)*nangy + iangy+1]
								if (iangx < nangx)
									m_shot[iz,ix] += w4*m[iz,(iangx)*nangy   + iangy+1]
								end
							end
						end
					end
				end
			end

			extent = Seismic.Extent(nz,nx,ny,1,1,
				oz,ox,oy,0.,0.,
				dz,dx,dy,1.,1.,
				"Depth","mx","my","","",
				"m","m","m","","",
				"")
			SeisWrite(shot_list[ishot].m,m_shot,h_shot,extent)
		end
		close(stream_m)
		close(stream_h)
		@sync @parallel for ishot = 1 : nshot
			a = shotwem(shot_list[ishot])
		end
		extent = Seismic.Extent(nt,nx,ny,nsx,nsy,
				ot,ox,oy,osx,osy,
				dt,dx,dy,dsx,dsy,
				"Time","gx","gy","sx","sy",
				"s","m","m","m","m",
				"")
		j = 1
		for ishot = 1 : nshot
			d_shot,h_shot,e = SeisRead(shot_list[ishot].d)
			SeisWrite(d,d_shot,h_shot,extent,itrace=j)
			j += size(d_shot,2)
			SeisRemove(shot_list[ishot].d)
			SeisRemove(shot_list[ishot].m)
			if (nangx != 1 || nangy != 1)
				SeisRemove(shot_list[ishot].angx)
				SeisRemove(shot_list[ishot].angy)
			end
		end
	end
end

struct Shot
	m
	d
	vel
	angx
	angy
	wav
	damping
	sx
	sy
	sz
	gz
	fmin
	fmax
	padt
	padx
	adj
	pspi
	nref
	verbose
end

function shotwem(shot)
	@compat if (Libdl.find_library(["shotwem"]) == "")
		error("Couldn't find shared library shotwem.so, make sure it is installed and check LD_LIBRARY_PATH environmental variable.")
	end
	argv = ["0",
	join(["adj=",shot.adj]),
	join(["pspi=",shot.pspi]),
	join(["nref=",shot.nref]),
	join(["d=",shot.d]), join(["m=",shot.m]), join(["vel=",shot.vel]),  join(["wav=",shot.wav]),
	join(["damping=",shot.damping]), join(["sx=",shot.sx]), join(["sy=",shot.sy]),  join(["sz=",shot.sz]),  join(["gz=",shot.gz]),
	join(["fmin=",shot.fmin]),  join(["fmax=",shot.fmax]),
	join(["padt=",shot.padt]),  join(["padx=",shot.padx]),
	join(["verbose=",shot.verbose]) ]
	@compat a = ccall((:main, "shotwem"), Int32, (Int32, Ptr{Ptr{UInt8}}), length(argv), argv)
	return(a)
end


text:
#Write a function that tests whether a string is a palindrome.

str = "toot"
str2 = "Alex"

function pal(x)
    x == reverse(x)
end

pal(str)

max_stars_repo_name: atott/Simple-Programming-Problems-Julia

max_issues_repo_name: atott/Simple-Programming-Problems-Julia

max_forks_repo_name: atott/Simple-Programming-Problems-Julia


text:
using ImmutableArrays
using Base.Test

typealias Vec2d Vector2{Float64}
typealias Vec3d Vector3{Float64}
typealias Vec4d Vector4{Float64}
typealias Vec3f Vector3{Float32}

v1 = Vec3d(1.0,2.0,3.0)
v2 = Vec3d(6.0,5.0,4.0)

# indexing
@test v1[1] == 1.0
@test v1[2] == 2.0
@test v1[3] == 3.0
@test_throws BoundsError v1[-1]
@test_throws BoundsError v1[0]
@test_throws BoundsError v1[4]

# negation
@test -v1 == Vec3d(-1.0,-2.0,-3.0)
@test isa(-v1,Vec3d)

# addition
@test v1+v2 == Vec3d(7.0,7.0,7.0)

# subtraction
@test v2-v1 == Vec3d(5.0,3.0,1.0)

# multiplication
@test v1.*v2 == Vec3d(6.0,10.0,12.0)

# division
@test v1 ./ v1 == Vec3d(1.0,1.0,1.0)

# scalar operations
@test 1.0 + v1 == Vec3d(2.0,3.0,4.0)
@test 1.0 .+ v1 == Vec3d(2.0,3.0,4.0)
@test v1 + 1.0 == Vec3d(2.0,3.0,4.0)
@test v1 .+ 1.0 == Vec3d(2.0,3.0,4.0)
@test 1 + v1 == Vec3d(2.0,3.0,4.0)
@test 1 .+ v1 == Vec3d(2.0,3.0,4.0)
@test v1 + 1 == Vec3d(2.0,3.0,4.0)
@test v1 .+ 1 == Vec3d(2.0,3.0,4.0)

@test v1 - 1.0 == Vec3d(0.0,1.0,2.0)
@test v1 .- 1.0 == Vec3d(0.0,1.0,2.0)
@test 1.0 - v1 == Vec3d(0.0,-1.0,-2.0)
@test 1.0 .- v1 == Vec3d(0.0,-1.0,-2.0)
@test v1 - 1 == Vec3d(0.0,1.0,2.0)
@test v1 .- 1 == Vec3d(0.0,1.0,2.0)
@test 1 - v1 == Vec3d(0.0,-1.0,-2.0)
@test 1 .- v1 == Vec3d(0.0,-1.0,-2.0)

@test 2.0 * v1 == Vec3d(2.0,4.0,6.0)
@test 2.0 .* v1 == Vec3d(2.0,4.0,6.0)
@test v1 * 2.0 == Vec3d(2.0,4.0,6.0)
@test v1 .* 2.0 == Vec3d(2.0,4.0,6.0)
@test 2 * v1 == Vec3d(2.0,4.0,6.0)
@test 2 .* v1 == Vec3d(2.0,4.0,6.0)
@test v1 * 2 == Vec3d(2.0,4.0,6.0)
@test v1 .* 2 == Vec3d(2.0,4.0,6.0)

@test v1 / 2.0 == Vec3d(0.5,1.0,1.5)
@test v1 ./ 2.0 == Vec3d(0.5,1.0,1.5)
@test v1 / 2 == Vec3d(0.5,1.0,1.5)
@test v1 ./ 2 == Vec3d(0.5,1.0,1.5)

@test 12.0 ./ v1 == Vec3d(12.0,6.0,4.0)
@test 12 ./ v1 == Vec3d(12.0,6.0,4.0)

@test v1.^2 == Vec3d(1.0,4.0,9.0)
@test v1.^2.0 == Vec3d(1.0,4.0,9.0)
@test 2.0.^v1 == Vec3d(2.0,4.0,8.0)
@test 2.^v1 == Vec3d(2.0,4.0,8.0)

# vector norm
@test norm(Vec3d(1.0,2.0,2.0)) == 3.0

# cross product
@test cross(v1,v2) == Vec3d(-7.0,14.0,-7.0)
@test isa(cross(v1,v2),Vec3d)
x = unit(Vec2d,1)
y = unit(Vec2d,2)
@test cross(x,y) == 1.0
@test cross(y,x) == -1.0

# basis vectors
e1 = unit(Vec4d,1)
e2 = unit(Vec4d,2)
e3 = unit(Vec4d,3)
e4 = unit(Vec4d,4)
@test e1 == Vec4d(1.0,0.0,0.0,0.0)
@test e2 == Vec4d(0.0,1.0,0.0,0.0)
@test e3 == Vec4d(0.0,0.0,1.0,0.0)
@test e4 == Vec4d(0.0,0.0,0.0,1.0)

# type conversion
@test isa(convert(Vec3f,v1),Vec3f)
@test Vector3([1.0,2.0,3.0]) == v1
@test convert(Vec3d,[1.0,2.0,3.0]) == v1
@test isa(convert(Vector{Float64},v1),Vector{Float64})
@test convert(Vector{Float64},v1) == [1.0,2.0,3.0]

# matrix operations

typealias Mat1d Matrix1x1{Float64}
typealias Mat2d Matrix2x2{Float64}
typealias Mat3d Matrix3x3{Float64}
typealias Mat4d Matrix4x4{Float64}

@test zero(Mat2d) == Mat2d(Vec2d(0.0,0.0),Vec2d(0.0,0.0))

v = Vec4d(1.0,2.0,3.0,4.0)
r = row(v)
c = column(v)

@test c*r == Mat4d(Vec4d(1.0,2.0,3.0,4.0),
                     Vec4d(2.0,4.0,6.0,8.0),
                     Vec4d(3.0,6.0,9.0,12.0),
                     Vec4d(4.0,8.0,12.0,16.0))
@test r*c == Matrix1x1(30.0)
@test r' == c
@test c' == r
@test row(r,1) == v
@test column(c,1) == v
@test row(r+c',1) == 2*v
@test sum(r) == sum(v)
@test prod(c) == prod(v)
@test eye(Mat3d) == Mat3d(Vec3d(1.0,0.0,0.0),Vec3d(0.0,1.0,0.0),Vec3d(0.0,0.0,1.0))
@test v*eye(Mat4d)*v == 30.0
@test -r == -1.0*r
@test diag(diagm(v)) == v

# type conversion
@test isa(convert(Matrix1x4{Float32},r),Matrix1x4{Float32})
jm = rand(4,4)
im = Matrix4x4(jm)
@test isa(im,Mat4d)
@test jm == im
im = convert(Mat4d,jm)
@test isa(im,Mat4d)
@test jm == im
jm2 = convert(Array{Float64,2},im)
@test isa(jm2,Array{Float64,2})
@test jm == jm2
@test isa(convert(Vector2{Float32},[4.;5.]).e1, Float32)

# Issue #44
let
    A = Matrix2x2(eye(2))
    @test isa(A^0, Matrix2x2)
    @test isa(A^1, Matrix2x2)
    @test isa(A^2, Matrix2x2)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiPolicyV1beta1PodSecurityPolicy <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    spec::Any # spec type: Union{ Nothing, IoK8sApiPolicyV1beta1PodSecurityPolicySpec } # spec name: spec

    function IoK8sApiPolicyV1beta1PodSecurityPolicy(;apiVersion=nothing, kind=nothing, metadata=nothing, spec=nothing)
        o = new()
        validate_property(IoK8sApiPolicyV1beta1PodSecurityPolicy, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiPolicyV1beta1PodSecurityPolicy, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiPolicyV1beta1PodSecurityPolicy, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiPolicyV1beta1PodSecurityPolicy, Symbol("spec"), spec)
        setfield!(o, Symbol("spec"), spec)
        o
    end
end # type IoK8sApiPolicyV1beta1PodSecurityPolicy

const _property_map_IoK8sApiPolicyV1beta1PodSecurityPolicy = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("spec")=>Symbol("spec"))
const _property_types_IoK8sApiPolicyV1beta1PodSecurityPolicy = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("spec")=>"IoK8sApiPolicyV1beta1PodSecurityPolicySpec")
Base.propertynames(::Type{ IoK8sApiPolicyV1beta1PodSecurityPolicy }) = collect(keys(_property_map_IoK8sApiPolicyV1beta1PodSecurityPolicy))
Swagger.property_type(::Type{ IoK8sApiPolicyV1beta1PodSecurityPolicy }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiPolicyV1beta1PodSecurityPolicy[name]))}
Swagger.field_name(::Type{ IoK8sApiPolicyV1beta1PodSecurityPolicy }, property_name::Symbol) =  _property_map_IoK8sApiPolicyV1beta1PodSecurityPolicy[property_name]

function check_required(o::IoK8sApiPolicyV1beta1PodSecurityPolicy)
    true
end

function validate_property(::Type{ IoK8sApiPolicyV1beta1PodSecurityPolicy }, name::Symbol, val)
end


text:
import Gaugefields:Traceless_antihermitian_add!

struct DomainwallFermiAction{Dim,Dirac,fermion,gauge} <: FermiAction{Dim,Dirac,fermion,gauge}
    hascovnet::Bool
    covneuralnet::Union{Nothing,CovNeuralnet{Dim}}
    diracoperator::Dirac
    _temporary_fermionfields::Vector{fermion}
    _temporary_gaugefields::Vector{gauge}

    function DomainwallFermiAction(D::Dirac_operator{Dim},hascovnet,covneuralnet) where Dim
        num = 8
        temps = get_temporaryvectors(D)
        x = temps[1]
        xtype = typeof(x)
        _temporary_fermionfields = Array{xtype,1}(undef,num)
        for i=1:num
            _temporary_fermionfields[i] = similar(x)
        end

        Utemp = D.U[1]
        Utype = typeof(Utemp)
        numU = 2
        _temporary_gaugefields = Array{Utype,1}(undef,numU)
        for i=1:numU
            _temporary_gaugefields[i] = similar(Utemp)
        end

        return new{Dim,typeof(D),xtype,Utype }(hascovnet,covneuralnet,D,
                        _temporary_fermionfields,
                        _temporary_gaugefields)

    end
end

function evaluate_FermiAction(fermi_action::DomainwallFermiAction{Dim,Dirac,fermion,gauge },U,ϕ::AbstractFermionfields) where {Dim,Dirac,fermion,gauge }
    W = fermi_action.diracoperator(U)
    η = fermi_action._temporary_fermionfields[1]
    solve_DinvX!(η,W',ϕ)
    Sf = dot(η,η)
    return real(Sf)
end

function calc_UdSfdU!(UdSfdU::Vector{<: AbstractGaugefields},fermi_action::DomainwallFermiAction{Dim,Dirac,fermion,gauge } ,U::Vector{<: AbstractGaugefields},ϕ::AbstractFermionfields) where  {Dim,Dirac,fermion,gauge }
    #println("------dd")
    QD5DW = fermi_action.diracoperator.D5DW(U)
    Q = D5DWdagD5DW_Wilson_operator(QD5DW)
    D5_PV = fermi_action.diracoperator.D5DW_PV(U)

    temps_dw =fermi_action._temporary_fermionfields[1]

    X = fermi_action._temporary_fermionfields[6]
    Y = fermi_action._temporary_fermionfields[5]

    mul!(temps_dw,D5_PV',ϕ)

    solve_DinvX!(X,Q,temps_dw)
    #set_wing_fermion!(X)
    set_wing_fermion!(ϕ)

    clear_U!(UdSfdU)

    calc_UdSfdU_fromX!(UdSfdU,Y,ϕ,fermi_action,U,X)

    #println("----aa--")
    set_wing_U!(UdSfdU)
end

function calc_UdSfdU_fromX!(UdSfdU::Vector{<: AbstractGaugefields},Y,ϕ,
        fermi_action::DomainwallFermiAction{Dim,Dirac,fermion,gauge } ,U,X;coeff = 1) where {Dim,Dirac,fermion,gauge }
    W = fermi_action.diracoperator.D5DW(U)
    temps_dw = fermi_action._temporary_fermionfields[2]
    mul!(temps_dw,W,X)
    clear_fermion!(Y)
    add_fermion!(Y,-1,ϕ,1,temps_dw)
    set_wing_fermion!(Y)

    temp0_g = fermi_action._temporary_gaugefields[1]

    κ = 1/2
    Dwilson = W.wilsonoperator

    L5 = fermi_action.diracoperator.D5DW.L5

    if L5 != X.L5
        @assert L5 % 2 == 0
        irange = Int64[]
        irange_out = Int64[]
        #irange = 1:L5
        #irange_out = (L5+1):X.L5

        for i5=1:X.L5
            if i5 <= div(L5,2) || i5 >= X.L5-div(L5,2)+1
                push!(irange,i5)
            else
                push!(irange_out,i5)
            end

        end

    else
        irange = 1:L5
    end

#    for i5=1:X.L5
    for i5 in irange

        temp0_f = fermi_action._temporary_fermionfields[1].w[i5] #F_field
        temp1_f = fermi_action._temporary_fermionfields[2].w[i5] #F_field

        for μ=1:Dim
            #!  Construct U(x,mu)*P1
            Xs = X.w[i5]
            Ys = Y.w[i5]

            # U_{k,μ} X_{k+μ}
            Xsplus = shift_fermion(Xs,μ)

            #@time mul!(temp0_f,U[μ],X)
            mul!(temp0_f,U[μ],Xsplus)

            # (r-γ_μ) U_{k,μ} X_{k+μ}
            mul!(temp1_f,view(Dwilson.rminusγ,:,:,μ),temp0_f)

            # κ (r-γ_μ) U_{k,μ} X_{k+μ}
            mul!(temp0_f,κ,temp1_f)

            # κ ((r-γ_μ) U_{k,μ} X_{k+μ}) ⊗ Y_k
            mul!(temp0_g,temp0_f,Ys')

            add_U!(UdSfdU[μ],-coeff,temp0_g)

            #!  Construct P2*U_adj(x,mu)
            # Y_{k+μ}^dag U_{k,μ}^dag
            Ysplus = shift_fermion(Ys,μ)
            mul!(temp0_f,Ysplus',U[μ]')

            # Y_{k+μ}^dag U_{k,μ}^dag*(r+γ_μ)
            mul!(temp1_f,temp0_f,view(Dwilson.rplusγ,:,:,μ))

            # κ Y_{k+μ}^dag U_{k,μ}^dag*(r+γ_μ)
            mul!(temp0_f,κ,temp1_f)

            # X_k ⊗ κ Y_{k+μ}^dag U_{k,μ}^dag*(r+γ_μ)
            #println(getvalue(temp0_g,1,1,1,1,1,1))
            #println(temp0_g[1,1,1,1,1,1])
            mul!(temp0_g,Xs,temp0_f)
            #println("temp0_g , ",getvalue(temp0_g,1,1,1,1,1,1))
            #println("temp0_g , ",temp0_g[1,1,1,1,1,1])
            #println("coeff ", coeff)
            #println("before ",getvalue(UdSfdU[μ],1,1,1,1,1,1))
            #println("before ",UdSfdU[μ][1,1,1,1,1,1])

            add_U!(UdSfdU[μ],coeff,temp0_g)
            #error("h")
            #println("after ",getvalue(UdSfdU[μ],1,1,1,1,1,1))
            #println("after ",UdSfdU[μ][1,1,1,1,1,1])

        end

    end

end

function calc_p_UdSfdU!(p,
    fermi_action::DomainwallFermiAction{Dim,Dirac,fermion,gauge } ,U::Vector{<: AbstractGaugefields},ϕ::AbstractFermionfields,coeff = 1) where {Dim,Dirac,fermion,gauge }
    #println("------dd")
    QD5DW = fermi_action.diracoperator.D5DW(U)
    Q = D5DWdagD5DW_Wilson_operator(QD5DW)
    D5_PV = fermi_action.diracoperator.D5DW_PV(U)

    temps_dw =fermi_action._temporary_fermionfields[1]

    X = fermi_action._temporary_fermionfields[end]
    Y = fermi_action._temporary_fermionfields[end-1]

    mul!(temps_dw,D5_PV',ϕ)

    solve_DinvX!(X,Q,temps_dw)

    #set_wing_fermion!(ϕ)

    calc_p_UdSfdU_fromX!(p,Y,ϕ,fermi_action,U,X,coeff= coeff)
    #println("----aa--")
    #set_wing_U!(UdSfdU)
end

function calc_p_UdSfdU_fromX!(p,Y,ϕ,
    fermi_action::DomainwallFermiAction{Dim,Dirac,fermion,gauge }  ,U,X;coeff = 1) where {Dim,Dirac,fermion,gauge }
    W = fermi_action.diracoperator.D5DW(U)
    temps_dw = fermi_action._temporary_fermionfields[1]
    mul!(temps_dw,W,X)
    clear_fermion!(Y)
    add_fermion!(Y,-1,ϕ,1,temps_dw)
    set_wing_fermion!(Y)

    temp0_g = fermi_action._temporary_gaugefields[1]

    κ = 1/2
    Dwilson = W.wilsonoperator
    for i5=1:X.L5

        temp0_f = fermi_action._temporary_fermionfields[1].w[i5] #F_field
        temp1_f = fermi_action._temporary_fermionfields[2].w[i5] #F_field

        for μ=1:Dim
            #!  Construct U(x,mu)*P1
            Xs = X.w[i5]
            Ys = Y.w[i5]

            # U_{k,μ} X_{k+μ}
            Xsplus = shift_fermion(Xs,μ)

            #@time mul!(temp0_f,U[μ],X)
            mul!(temp0_f,U[μ],Xsplus)

            # (r-γ_μ) U_{k,μ} X_{k+μ}
            mul!(temp1_f,view(Dwilson.rminusγ,:,:,μ),temp0_f)

            # κ (r-γ_μ) U_{k,μ} X_{k+μ}

            mul!(temp0_f,κ,temp1_f)

            # κ ((r-γ_μ) U_{k,μ} X_{k+μ}) ⊗ Y_k
            mul!(temp0_g,temp0_f,Ys')

            Traceless_antihermitian_add!(p[μ],-coeff,temp0_g)
            #println(p[μ][1,1,1,1,1])

            #!  Construct P2*U_adj(x,mu)
            # Y_{k+μ}^dag U_{k,μ}^dag
            Ysplus = shift_fermion(Ys,μ)
            mul!(temp0_f,Ysplus',U[μ]')

            # Y_{k+μ}^dag U_{k,μ}^dag*(r+γ_μ)
            mul!(temp1_f,temp0_f,view(Dwilson.rplusγ,:,:,μ))

            # κ Y_{k+μ}^dag U_{k,μ}^dag*(r+γ_μ)
            mul!(temp0_f,κ,temp1_f)

            # X_k ⊗ κ Y_{k+μ}^dag U_{k,μ}^dag*(r+γ_μ)
            mul!(temp0_g,Xs,temp0_f)

            Traceless_antihermitian_add!(p[μ],coeff,temp0_g)
        end
    end

end

function gauss_sampling_in_action!(η::AbstractFermionfields,U,fermi_action::DomainwallFermiAction{Dim,Dirac,fermion,gauge } ) where {Dim,Dirac,fermion,gauge }
    #gauss_distribution_fermion!(η)
    gauss_distribution_fermion!(η,rand)
end

using InteractiveUtils

function sample_pseudofermions!(ϕ::AbstractFermionfields,U,fermi_action::DomainwallFermiAction{Dim,Dirac,fermion,gauge } ,ξ::AbstractFermionfields) where {Dim,Dirac,fermion,gauge }
    W = fermi_action.diracoperator(U)
    mul!(ϕ,W',ξ)
    set_wing_fermion!(ϕ)
end


text:
# abstract type for matrix-valued kernels
abstract type MatrixKernel <: AbstractKernel end

# uniform scaling of scalar-valued kernel
struct UniformScalingKernel{T<:Real,K} <: MatrixKernel
    λ::T
    kernel::K

    function UniformScalingKernel{T,K}(λ::T, kernel::K) where {T,K}
        λ ≥ zero(λ) || throw(ArgumentError("`λ` must be non-negative"))

        new{T,K}(λ, kernel)
    end
end

UniformScalingKernel(λ::Real, kernel) =
    UniformScalingKernel{typeof(λ),typeof(kernel)}(λ, kernel)
UniformScalingKernel(kernel) = UniformScalingKernel(1, kernel)

(kernel::UniformScalingKernel)(x, y) = UniformScaling(kernel.λ * kernel.kernel(x, y))

kernel_result_type(kernel::UniformScalingKernel, x, y) =
    typeof(zero(kernel.λ) * kernel.kernel(zero(eltype(x)), zero(eltype(y))))

# scaling of scalar-valued kernel by diagonal matrix
struct DiagonalKernel{V<:AbstractVector{<:Real},K} <: MatrixKernel
    diag::V
    kernel::K

    function DiagonalKernel{V,K}(diag::V, kernel::K) where {V,K}
        all(x -> x ≥ zero(x), diag) ||
            throw(ArgumentError("all entries of `diag` must be non-negative"))

        new{V,K}(diag, kernel)
    end
end

DiagonalKernel(diag::AbstractVector{<:Real}, kernel) =
    DiagonalKernel{typeof(diag),typeof(kernel)}(diag, kernel)

(kernel::DiagonalKernel)(x, y) = Diagonal(kernel.diag .* kernel.kernel(x, y))

kernel_result_type(kernel::DiagonalKernel, x, y) =
    typeof(zero(eltype(kernel.diag)) * kernel.kernel(zero(eltype(x)), zero(eltype(y))))


text:
exp_name = "AL"

data_dirs = [
             # "ALOI", # too many observations
             "Annthyroid",
             # "Arrhythmia", # too many attributes
             "Cardiotocography",
             "Glass",
             "HeartDisease",
             "Hepatitis",
             "Ionosphere",
             # "InternetAds", # too many attributes
             # "KDDCup99", # too many observations
             "Lymphography",
             "PageBlocks",
             "Parkinson",
             "PenDigits",
             "Pima",
             "Shuttle",
             "SpamBase",
             "Stamps",
             "Waveform",
             "WBC",
             "WDBC",
             "WPBC"]

solver = with_optimizer(Gurobi.Optimizer; OutputFlag=0, Threads=1)
oracle_param = Dict{Symbol, Any}( :type => PoolOracle, :param => Dict{Symbol, Any}())

quality_metrics = Dict(:kappa => cohens_kappa, :f1 => f1_score)
selection_metric = :kappa

n_iter = 50
n_c_steps = 20

learning_methods = vcat(
    # Kernel Active Learning with custom selection strategy
    [ALKernel(nn_stability = nns,
             stability_window = 10,
             stability_threshold = st,
             C_steps = n_c_steps,
             max_iteration = n_iter,
             initial_n = 4,
             max_gamma = 100.0,
             adjusted = true) for nns in [1, 5, 10] for st in [0.0]],
    # Kernel Active Learning with random selection strategy
    [ALRandom(nn_stability = nns,
                       n_samples = n_iter,
                       max_gamma = 100.0,
                       adjusted = true,
                       C_steps = n_c_steps) for nns in [1, 5, 10]])

seeds = [42 * i for i in 1:5]

max_stars_repo_name: holtri/lama

max_issues_repo_name: holtri/lama

max_forks_repo_name: holtri/lama


text:
include("../models/bernoulliModel.jl")

setRNGs(12345)

θ = 0.4

model = makeBernoulliModel(θ)
ys = RecyclePermanent.simulateModel(model, 100)
model = makeBernoulliModel(θ, ys)

N = 1000*model.n

trueValue = sum(ys)*log(θ)+(model.n-sum(ys))*log(1-θ)

kmio = KMIO{model.Tx}(N, model.n)
@test kmEstimate(model, kmio) ≈ trueValue atol=0.1

rpio = RPIO{model.Tx}(N, model.n, 1)
@test simpleEstimate(model, rpio) ≈ trueValue atol=0.1
@test biasedEstimate(model, rpio) ≈ trueValue atol=0.1
@test biasedEstimateDiscrete(model, rpio) ≈ trueValue atol=0.1
@test rpEstimate(model, rpio) ≈ trueValue atol=0.1
@test rpEstimateDiscrete(model, rpio) ≈ trueValue atol=0.1

rpio = RPIO{model.Tx}(N, model.n, Threads.nthreads())
@test simpleEstimate(model, rpio) ≈ trueValue atol=0.1
@test biasedEstimate(model, rpio) ≈ trueValue atol=0.1
@test biasedEstimateDiscrete(model, rpio) ≈ trueValue atol=0.1
@test rpEstimate(model, rpio) ≈ trueValue atol=0.1
@test rpEstimateDiscrete(model, rpio) ≈ trueValue atol=0.1

N = 3*model.n
trials = 10000
vals = Vector{Float64}(undef, trials)

t1 = 1+(θ/θ^2-1)/N
t2 = 1+((1-θ)/(1-θ)^2-1)/N
trueVar = t1^sum(ys)*t2^(model.n - sum(ys)) - 1

rpio = RPIO{model.Tx}(N, model.n, Threads.nthreads())

for i in 1:trials
  vals[i] = exp(simpleEstimate(model, rpio) - trueValue)
end
@test mean(vals) ≈ 1.0 atol=0.1
@test var(vals) ≈ trueVar atol=0.1

for i in 1:trials
  vals[i] = exp(rpEstimate(model, rpio) - trueValue)
end
@test mean(vals) ≈ 1.0 atol=0.1

for i in 1:trials
  vals[i] = exp(rpEstimateDiscrete(model, rpio) - trueValue)
end
@test mean(vals) ≈ 1.0 atol=0.1

kmio = KMIO{model.Tx}(N, model.n)
for i in 1:trials
  vals[i] = exp(kmEstimate(model, kmio) - trueValue)
end
@test mean(vals) ≈ 1.0 atol=0.1


text:
# Copyright 2020 Google LLC

# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at

#     https://www.apache.org/licenses/LICENSE-2.0

# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
raw_input = readlines("day03input.txt")

@enum Square tree snow

function parse(lines)
    h = length(lines)
    w = length(first(lines))
    area = fill(snow, (w, h))
    for y in 1:h
        for x in 1:w
            if lines[y][x] == '.'
                area[x, y] = snow
            elseif lines[y][x] == '#'
                area[x, y] = tree
            end
        end
    end
    return area
end

input = parse(raw_input)

function count_trees(input, step::Tuple{Int,Int})
    w, h = size(input)
    dx, dy = step
    pos = (1, 1)
    n = 0
    while pos[2] <= h
        if input[(pos[1] - 1) % w + 1, pos[2]] == tree
            n += 1
        end
        pos = (pos[1] + dx, pos[2] + dy)
    end
    return n
end

part1(input) = count_trees(input, (3, 1))

function part2(input)
    slopes = [(1, 1), (3, 1), (5, 1), (7, 1), (1, 2)]
    reduce(
        *,
        map(slope -> count_trees(input, slope), slopes))
end

println(part1(input))

println(part2(input))


text:
using Many
using Base.Test

@test collect( Many.hours(DateTime(2016,04,04,15), 3)) |> length == 3
@test collect( Many.hours(DateTime(2016,04,04,15), Dates.Hour(3))) |> length == 3
@test first( collect( Many.hours(DateTime(2016,04,04,15), Dates.Hour(3))) ) == DateTime(2016,4,4,15)
@test last( collect(Many.hours(DateTime(2016,04,04,15), Dates.Day(1))))==DateTime(2016,4,5,14)
@test last(collect(Many.hours(Date(2016,04,04), Dates.Day(1)))) == DateTime(2016,4,4,23)

max_stars_repo_name: closescreen/Many

max_issues_repo_name: closescreen/Many

max_forks_repo_name: closescreen/Many


text:
using Test, DiffEqBase
using DiffEqBase: __has_jac, __has_tgrad, __has_Wfact, __has_Wfact_t,
  __has_paramjac, __has_syms, __has_analytic, __has_colorvec,  has_jac, has_tgrad,
  has_Wfact, has_Wfact_t, has_paramjac, has_syms, has_analytic, has_colorvec,
  AbstractDiffEqFunction

struct Foo <: AbstractDiffEqFunction{false}
  jac
  tgrad
  Wfact
  Wfact_t
  paramjac
  syms
  analytic
  colorvec
end

f = Foo(1,1,1,1,1,1,1,1)

@test __has_jac(f)
@test __has_tgrad(f)
@test __has_Wfact(f)
@test __has_Wfact_t(f)
@test __has_paramjac(f)
@test __has_syms(f)
@test __has_analytic(f)
@test __has_colorvec(f)

@test has_jac(f)
@test has_tgrad(f)
@test has_Wfact(f)
@test has_Wfact_t(f)
@test has_paramjac(f)
@test has_syms(f)
@test has_analytic(f)
@test has_colorvec(f)

struct Foo2 <: AbstractDiffEqFunction{false}
  jac
  tgrad
  Wfact
  Wfact_t
end

f2 = Foo2(1,1,nothing,nothing)

@test __has_jac(f2)
@test __has_tgrad(f2)
@test __has_Wfact(f2)
@test __has_Wfact_t(f2)
@test !__has_paramjac(f2)
@test !__has_syms(f2)
@test !__has_analytic(f2)
@test !__has_colorvec(f2)

@test has_jac(f2)
@test has_tgrad(f2)
@test !has_Wfact(f2)
@test !has_Wfact_t(f2)
@test !has_paramjac(f2)
@test !has_syms(f2)
@test !has_analytic(f2)
@test !has_colorvec(f2)


text:
using FlightSims
const FS = FlightSims
using UnPack
using Transducers
using Plots

function test()
    power_loop = PowerLoop()
    cg = Command(power_loop)
    @unpack t0, t_go_straight, t_loop = power_loop
    Δt = 0.01
    tf = t0 + 2*t_go_straight + t_loop + 5
    ts = t0:Δt:tf
    ps = ts |> Map(cg) |> collect
    ps_cat = hcat(ps...)'
    p_traj = plot3d(1;
                    xlim=(-20, 20),
                    ylim=(-20, 20),
                    zlim=(-20, 20),
                    aspect_ratio=:equal,
                    label="trajectory",
                    title="Power Loop")
    anim = Animation()
    for p in ps[1:10:end]
        push!(p_traj, p...)
        frame(anim)
    end
    gif(anim, "power_loop.gif")
end


text:
module GlobalSensitivity

using Statistics, RecursiveArrayTools, LinearAlgebra, Random
using QuasiMonteCarlo, ForwardDiff, KernelDensity, Trapz
using Parameters: @unpack
using FFTW, Distributions, StatsBase

abstract type GSAMethod end

include("morris_sensitivity.jl")
include("sobol_sensitivity.jl")
include("regression_sensitivity.jl")
include("DGSM_sensitivity.jl")
include("eFAST_sensitivity.jl")
include("delta_sensitivity.jl")
include("easi_sensitivity.jl")
include("rbd-fast_sensitivity.jl")
include("fractional_factorial_sensitivity.jl")

export gsa

export Sobol, Morris, RegressionGSA, DGSM, eFAST, DeltaMoment, EASI, FractionalFactorial

end # module


text:
using ImageCore, Colors, FixedPointNumbers, OffsetArrays
using Test, Random

@info "Deprecations are expected"
@testset "permuteddimsview" begin
    a = [1 3; 2 4]
    v = permuteddimsview(a, (1,2))
    @test v == a
    v = permuteddimsview(a, (2,1))
    @test v == a'
    a = rand(3,7,5)
    v = permuteddimsview(a, (2,3,1))
    @test v == permutedims(a, (2,3,1))
end


text:
"""
    RCLogit <: GMMModel

A random coefficients logit model with endogeneity.
An `RCLogit` model consists of outcomes, `y` ∈ (0,1),  regressors
`x`, instruments `z`, and random draws `ν ∼ N(0,I)`.  The moment condition is

``E[ξz] = 0``

where

`` y = ∫ exp(x(β + ν) + ξ)/(1 + exp(x(β + ν) + ξ)) dΦ(ν;Σ) ``

where Φ(ν;Σ) is the normal distribution with variance Σ.

The dimensions of `x`, `y`, `z`, and `ν` must be such that
`length(y) == size(x,1) == size(z,1) == size(ν,2)`
and
`size(ν,3) == size(x,2) ≤ size(z,2)`.
"""
struct RCLogit <: GMMModel
  x::Matrix{Float64}
  y::Vector{Float64}
  z::Matrix{Float64}
  ν::Array{Float64,3}
end

"""
    RCLogit(n::Integer, β::AbstractVector,
            π::AbstractMatrix, Σ::AbstractMatrix,
            ρ, nsim=100)

Simulates a RCLogit model.

# Arguments

- `n` number of observations
- `β` mean coefficients on `x`
- `π` first stage coefficients `x = z*π + e`
- `Σ` variance of random coefficients
- `ρ` correlation between x[:,1] and structural error.
- `nsim` number of draws of `ν` for monte-carlo integration
"""
function RCLogit(n::Integer, β::AbstractVector,
                 π::AbstractMatrix, Σ::AbstractMatrix,
                 ρ, nsim=100)
  z = randn(n, size(π)[1])
  endo = randn(n, length(β))
  x = z*π .+ endo
  ξ = rand(Normal(0,sqrt((1.0-ρ^2))),n).+endo[:,1]*ρ
  ν = randn(nsim,n,length(β))
  U = cholesky(Σ).U
  y = zeros(n)
  for s in 1:nsim
    y .= y .+ cdf.(Logistic(), x*β .+ sum(x.*(ν[s,:,:]*U),dims=2) .+ ξ)[:]
  end
  y  .= y./nsim
  #ν = randn(nsim÷10,n,length(β))
  RCLogit(x,y,z,ν)
end

number_parameters(model::RCLogit) = size(model.x,2) +
  size(model.x,2)*(size(model.x,2)-1) ÷ 2 + number_observations(model)
number_observations(model::RCLogit) = length(model.y)
number_moments(model::RCLogit) = size(model.z,2)

function get_gi(model::RCLogit)
  function(θ)
    n = number_observations(model)
    ξ = θ[(length(θ)-n+1):end]
    ξ.*model.z
  end
end

"""
    gmm_constraints(model::RCLogit)

Returns

`` c(θ) = ∫ exp(x(β + ν) + ξ)/(1 + exp(x(β + ν) + ξ)) dΦ(ν;Σ) - y``

where `θ = [β, uvec, ξ]` with `uvec = vector(cholesky(Σ).U)`.

The integral is computed by monte-carlo integration.
"""
function gmm_constraints(model::RCLogit)
  function(θ)
    (nsim,n,K) = size(model.ν)
    β = θ[1:K]
    uvec = θ[(K+1):(K+(K*(K+1))÷2)]
    U =  zeros(eltype(uvec),K,K)
    l = 1
    for i in 1:K
      for j in i:K
        U[i,j] = uvec[l]
        l = l+1
      end
    end
    ξ = θ[((K+(K*(K-1))÷2)+1):end]
    share = zeros(eltype(θ), n)
    for s in 1:nsim
      share += cdf.(Logistic(), model.x*β + sum(model.x.*(model.ν[s,:,:]*U),dims=2) .+ ξ)[:]
    end
    share /= nsim
    return(share - model.y)
  end
end


text:
using Revise
using EvoTrees
using DataFrames
using CSV
using MLJ

num_cols = [
    "ClientPeriod",
    "MonthlySpending",
    "TotalSpent"
];

cat_cols = [
    "Sex",
    "IsSeniorCitizen",
    "HasPartner",
    "HasChild",
    "HasPhoneService",
    "HasMultiplePhoneNumbers",
    "HasInternetService",
    "HasOnlineSecurityService",
    "HasOnlineBackup",
    "HasDeviceProtection",
    "HasTechSupportAccess",
    "HasOnlineTV",
    "HasMovieSubscription",
    "HasContractPhone",
    "IsBillingPaperless",
    "PaymentMethod"
];

all_feature_cols = [num_cols; cat_cols];
target_col = "Churn";

# ,types=Dict("Sex"=>CategoricalValue{String, UInt32})
# + tags=[]
df = DataFrame!(CSV.File("data/train.csv", pool=0.1, missingstrings=[" "]))
categorical!(df,[cat_cols;target_col]);
describe(df,:eltype,:nunique, :nmissing)
dropmissing!(df);
describe(df,:eltype,:nunique, :nmissing)

X_df = df[!, all_feature_cols];
y = df[!,target_col];

mach_x = machine(ContinuousEncoder(), X_df)
fit!(mach_x)
X = MLJ.transform(mach_x, X_df)

tree_model = EvoTreeClassifier(T=Float32, max_depth=9, λ=1, nrounds=2000, colsample=0.3, metric=:mlogloss)
@time fit_evotree(tree_model, Array(X), y, print_every_n=100);
tree_model = EvoTreeClassifier(T=Float64, max_depth=9, λ=1, nrounds=2000, colsample=0.3, metric=:mlogloss)
@time fit_evotree(tree_model, Array(X), y, print_every_n=100);

y_num = Float64.(y.refs) .- 1
tree_model = EvoTreeRegressor(T=Float32, loss=:logistic, max_depth=9, nrounds=2000, colsample=0.3, metric=:logloss)
@time fit_evotree(tree_model, Array(X), y_num, print_every_n=100);
tree_model = EvoTreeRegressor(T=Float64, loss=:logistic, max_depth=9, nrounds=2000, colsample=0.3, metric=:logloss)
@time fit_evotree(tree_model, Array(X), y_num, print_every_n=100);

mach = machine(tree_model, X, y)

using XGBoost

params_xgb = ["max_depth" => 8,
         "eta" => 0.1,
         "objective" => "multi:softmax",
         "subsample" => 1.0,
         "colsample_bytree" => 0.3,
         "tree_method" => "hist",
         "max_bin" => 64,
         "lambda" => 1.0,
         "num_class" => 2,
         "print_every_n" => 50]
metrics = ["mlogloss"]

Int.(y_num)
@time m_xgb = xgboost(Array(X), 2000, label=Int.(y_num), param=params_xgb, metrics=metrics, nthread=8, print_every_n=100, silent=0);


text:
# Test the branching of the SoilWaterParameterizations functions,
# test the instantiation of them works as expected, and that they
# return what we expect.
using MPI
using OrderedCollections
using StaticArrays
using Test

using CLIMAParameters
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

using ClimateMachine
using ClimateMachine.Land
using ClimateMachine.Land.SoilWaterParameterizations

@testset "Land water parameterizations" begin
    FT = Float64
    test_array = [0.5, 1.0]
    vg_model = vanGenuchten(FT;)
    mm = MoistureDependent{FT}()
    bc_model = BrooksCorey(FT;)
    hk_model = Haverkamp(FT;)

    #Use an array to confirm that extra arguments are unused.
    @test viscosity_factor.(Ref(ConstantViscosity{FT}()), test_array) ≈
          [1.0, 1.0]
    @test impedance_factor.(Ref(NoImpedance{FT}()), test_array) ≈ [1.0, 1.0]
    @test moisture_factor.(
        Ref(MoistureIndependent{FT}()),
        Ref(vg_model),
        test_array,
    ) ≈ [1.0, 1.0]

    viscosity_model = TemperatureDependentViscosity{FT}(; T_ref = FT(1.0))
    @test viscosity_factor(viscosity_model, FT(1.0)) == 1
    impedance_model = IceImpedance{FT}(; Ω = 2.0)
    @test impedance_factor(impedance_model, 0.5) == FT(0.1)

    imp_f = impedance_factor(impedance_model, 0.5)
    vis_f = viscosity_factor(viscosity_model, 1.0)
    m_f = moisture_factor(MoistureDependent{FT}(), vg_model, 1.0)
    Ksat = 1.0
    @test hydraulic_conductivity(Ksat, imp_f, vis_f, m_f) == FT(0.1)

    @test_throws Exception effective_saturation(0.5, -1.0, 0.0)
    @test effective_saturation(0.5, 0.25, 0.05) - 4.0 / 9.0 < eps(FT)

    n = vg_model.n
    m = 1.0 - 1.0 / n
    α = vg_model.α
    S_s = 0.001
    θ_r = 0.0
    ν = 1.0
    @test moisture_factor.(Ref(mm), Ref(vg_model), test_array) ==
          sqrt.(test_array) .*
          (FT(1) .- (FT(1) .- test_array .^ (FT(1) / m)) .^ m) .^ FT(2)
    @test moisture_factor.(Ref(mm), Ref(bc_model), test_array) ==
          test_array .^ (FT(2) * bc_model.m + FT(3))
    ψ =
        -(
            (test_array .^ (-FT(1) / hk_model.m) .- FT(1)) .*
            hk_model.α^(-hk_model.n)
        ) .^ (FT(1) / hk_model.n)
    @test moisture_factor.(Ref(mm), Ref(hk_model), test_array) ==
          hk_model.A ./ (hk_model.A .+ abs.(ψ) .^ hk_model.k)

    @test pressure_head.(
        Ref(vg_model),
        Ref(ν),
        Ref(S_s),
        Ref(θ_r),
        test_array,
        Ref(0.0),
    ) ≈ .-((-1 .+ test_array .^ (-1 / m)) .* α^(-n)) .^ (1 / n)
    #test branching in pressure head
    @test pressure_head(vg_model, ν, S_s, θ_r, 1.5, 0.0) == 500
    @test pressure_head.(
        Ref(hk_model),
        Ref(ν),
        Ref(S_s),
        Ref(θ_r),
        test_array,
        Ref(0.0),
    ) ≈ .-((-1 .+ test_array .^ (-1 / m)) .* α^(-n)) .^ (1 / n)
    m = FT(0.5)
    ψb = FT(0.1656)
    @test pressure_head(bc_model, ν, S_s, θ_r, 0.5, 0.0) ≈ -ψb * 0.5^(-m)
    @test volumetric_liquid_fraction.([0.5, 1.5], Ref(0.75)) ≈ [0.5, 0.75]

    @test inverse_matric_potential(
        bc_model,
        matric_potential(bc_model, FT(0.5)),
    ) - FT(0.5) < eps(FT)
    @test inverse_matric_potential(
        vg_model,
        matric_potential(vg_model, FT(0.5)),
    ) == FT(0.5)
    @test inverse_matric_potential(
        hk_model,
        matric_potential(hk_model, FT(0.5)),
    ) == FT(0.5)

    @test_throws Exception inverse_matric_potential(bc_model, 1)
    @test_throws Exception inverse_matric_potential(hk_model, 1)
    @test_throws Exception inverse_matric_potential(vg_model, 1)

    # Test that spatial dependence works properly
    vg_spatial = vanGenuchten(FT; n = (x) -> 2 * x, α = (x) -> x^2)
    @test supertype(typeof(vg_spatial.n)) == Function
    vg_point = vg_spatial(FT(2.0))
    @test typeof(vg_point.n) == FT
    @test vg_point.n == 4.0
    @test vg_point.α == 4.0
    @test vg_point.m == 0.75
end