text:
function slope_intercept(wn_dict::Dict{Any,Any}, wn::PyCall.PyObject, link_results::Dict{Any,Any}, node_results::Dict{Any,Any}, num_timesteps::Int64)
    slopes = Dict{String,Any}()
    intercepts = Dict{String, Any}()
    pump_off = Dict{String, Any}()
    slope_sum = 0
    intercept_sum = 0
    count = 0
    for (key, pump) in wn_dict["pumps"]
        name = pump["name"]
        flow = link_results["flowrate"][name][:values][1:num_timesteps]
        headloss = link_results["headloss"][name][:values][1:num_timesteps]
        efficiency = nothing
        try
            efficiency = pump["efficiency"][:points]
        catch
            wn_dict["options"]["energy"]["global_efficiency"] != nothing ? efficiency = wn_dict["options"]["energy"]["global_efficiency"]: efficiency = 0.65
        end
        if length(efficiency) > 1
            energy = [((1000.0 * 9.81 * headloss[i] * flow[i])/efficiency[1][1]) for i = 1:num_timesteps]
        else
            energy = [((1000.0 * 9.81 * headloss[i] * flow[i])/efficiency) for i = 1:num_timesteps]
        end
        flows = Array{Float64}(0)
        energies = Array{Float64}(0)
        for i= 1:length(energy)
            if energy[i] != 0.0 && flow[i] != 0.0
                flows = vcat(flows,flow[i])
                energies = vcat(energies, energy[i])
            end
        end
        if length(energies) != 0.0
            count = count +1
            intercept, slope = LeastSquares(flows, energies)
            slopes[key] = slope
            intercepts[key] = intercept
            slope_sum = slope_sum + slope
            intercept_sum = intercept_sum + intercept
        else
            pump_off[key] = key
        end
    end

    avg_slope = slope_sum/count
    avg_intercept = intercept_sum/count
    for (key, name) in pump_off
        slopes[key] = avg_slope
        intercepts[key] = avg_intercept
    end
    return slopes, intercepts
end

function LeastSquares(flows::Array{Float64}, energies::Array{Float64})
    #if pump never turns on, delete control that inhibits pump, calulation energy, then replace control
        avg_energy = mean(energies)
        a, b = linear_fit(flows, energies)
        error = Array{Float64}(0)
        for i = 1:length(energies)
            push!(error, energies[i] - (avg_energy +b *flows[i]))
        end
        avg_error = mean(error)
        intercept = avg_energy + avg_error
        return intercept, b
    # total_error = energies - (intercept + b * flows)
end


text:
using Compat

abstract type WeaveChunk end
abstract type Inline end

mutable struct WeaveDoc
    source::AbstractString
    basename::AbstractString
    path::AbstractString
    chunks::Array{WeaveChunk}
    cwd::AbstractString
    format
    doctype::AbstractString
    header_script::String
    header::Dict
    template::AbstractString
    css::AbstractString
    highlight_theme
    fig_path::AbstractString
    function WeaveDoc(source, chunks, header)
        path, fname = splitdir(abspath(source))
        basename = splitext(fname)[1]
        new(source, basename, path, chunks, "", nothing, "", "", header,
          "", "", Highlights.Themes.DefaultTheme, "")
    end
end

struct ChunkOutput
    code::AbstractString
    stdout::AbstractString
    displayed::AbstractString
    rich_output::AbstractString
    figures::Array{AbstractString}
end

mutable struct CodeChunk <: WeaveChunk
    content::AbstractString
    number::Int
    result_no::Int
    start_line::Int
    optionstring::AbstractString
    options::Dict{Symbol, Any}
    output::AbstractString
    rich_output::AbstractString
    figures::Array{AbstractString}
    result::Array{ChunkOutput}
    function CodeChunk(content, number, start_line, optionstring, options)
        new(rstrip(content) * "\n", number, 0, start_line, optionstring, options, "","", AbstractString[], ChunkOutput[])
    end
end

mutable struct DocChunk <: WeaveChunk
    content::Array{Inline}
    number::Int
    start_line::Int
    function DocChunk(text::AbstractString, number::Int, start_line::Int, inline_regex = nothing)
        chunks = parse_inline(text, inline_regex)
        new(chunks, number, start_line)
    end
end

mutable struct InlineText <: Inline
    content::AbstractString
    si::Int
    ei::Int
    number::Int
end

mutable struct InlineCode <: Inline
    content::AbstractString
    si::Int
    ei::Int
    number::Int
    output::AbstractString
    rich_output::AbstractString
    figures::Array{AbstractString}
    function InlineCode(content, si, ei, number)
        new(content, si, ei, number, "", "", AbstractString[])
    end
end

mutable struct TermResult
end

mutable struct ScriptResult
end

mutable struct CollectResult
end

ext: jl


text:
# will need to test with parallel capabilities enabled
addprocs(3)

@everywhere using NumericExtensions: sumsq
@everywhere using PLINK
@everywhere using Base.Test

# use this epsilon for comparisons of approximate equality
const my_eps = 4.0*eps()

#############################
### test the constructors ###
#############################

# we rely on testing data bundled with standard PLINK 1.07 download
# it has been precompressed into binary format
# first initialize what the data *should* look like
const num_case = 6
const num_snps = 3
const block    = iceil(num_case/4)
const tblock   = iceil(num_snps/4)
const x        = convert(SharedArray{Int8,1}, [-36; 15; -25; 15; 107; 1])
const xt       = convert(SharedArray{Int8,1}, [60; 39; 41; 31; 31; 15])

# now use all four constructors
const x1 = BEDFile(x, xt, num_case, num_snps, block, tblock)
const x2 = BEDFile("../data/test_binary.bed", "../data/test_binary_t.bed")
const x3 = BEDFile("../data/test_binary.bed", "../data/test_binary_t.bed", num_case, num_snps)
const x4 = BEDFile("../data/test_binary.bed", "../data/test_binary_t.bed", num_case, num_snps, block, tblock)

# did inner constructor work?
@test isequal(x1.x, x)
@test isequal(x1.p, num_snps)
@test isequal(x1.n, num_case)
@test isequal(x1.blocksize, block)

# did outer constructors yield same result?
@test isequal(x1,x2)
@test isequal(x1,x3)
@test isequal(x1,x4)

##################################
### test the utility functions ###
##################################

# test size functions
@test isequal(size(x1), (num_case, num_snps))
@test isequal(size(x1,1), num_case)
@test isequal(size(x1,2), num_snps)

###################################
### test decompression routines ###
###################################

# first initialize the correct results
Y = convert(SharedArray{Float64,2},
	[0.0 2.0 2.0
     2.0 NaN 1.0
 	 NaN 1.0 1.0
	 2.0 2.0 NaN
	 2.0 2.0 NaN
	 2.0 2.0 0.0])

const means = convert(SharedArray{Float64,1}, [1.6, 1.8, 1.0])
const invstds = convert(SharedArray{Float64,1}, [1.1180339887498947,2.2360679774997894,1.224744871391589])
for j = 1:num_snps
	for i = 1:num_case
		Y[i,j] = (Y[i,j] - means[j]) * invstds[j]
	end
end
Y[isnan(Y)] = 0.0
const Yt = Y'

# this tests the full decompression routine
Y1 = SharedArray(Float64, num_case, num_snps)
decompress_genotypes!(Y1,x1, means=means, invstds=invstds)

Y2 = SharedArray(Float64,0,0)
colidx = trues(num_snps)
rowidx = trues(num_case)
# this tests the indexed decompression
for j = 0:num_snps
	fill!(colidx, true)
	colidx[1:j] = false
	Y2          = SharedArray(Float64, x1.n, sum(colidx))
	decompress_genotypes!(Y2,x1,colidx,means=means,invstds=invstds)
	@test isequal(Y2, Y[:, colidx])
	Y2 = false
	gc()
end

x5 = copy(x1)
# now test several permutations of subsetting and decompression
for i = 0:num_case, j = 0:num_snps
	println("\n\t***testing subsetting for $(num_case-i) rows, $(num_snps-j) columns***\n")
	fill!(rowidx, true)
	fill!(colidx, true)
	rowidx[1:i] = false
	colidx[1:j] = false
	x5          = x1[rowidx, colidx]
	Y2          = SharedArray(Float64, x5.n, x5.p)

#	decompress_genotypes!(Y2,x5,colidx, means=means[colidx], invstds=invstds[colidx])
	decompress_genotypes!(Y2,x5, means=means[colidx], invstds=invstds[colidx])
	@test_approx_eq_eps(Y2, Y[rowidx, colidx], my_eps)
	Y2 = false
	gc()
end

####################################
### test linear algebra routines ###
####################################

# create temporary arrays for linear algebra routines
y       = SharedArray(Float64, x1.n, init = S -> S[localindexes(S)] = 1.0)
b       = SharedArray(Float64, x1.p, init = S -> S[localindexes(S)] = 1.0 ./ localindexes(S))
r       = SharedArray(Float64, x1.n, init = S -> S[localindexes(S)] = 0.0)
indices = trues(num_snps)

const XtY_out   = Yt*y
const Xb_out    = Y*b
const sumsq_out = convert(SharedArray, vec(NumericExtensions.sumsq(sdata(Y),1)))
const res_out   = y - Xb_out

XtY = xty(x1,y)
Xb  = xb(x1, b, indices, x1.p)
ssq = PLINK.sumsq(x1)
PLINK.update_partial_residuals!(r, y, x1, indices, b, x1.p, Xb=Xb)

@test_approx_eq_eps(XtY_out, XtY, my_eps)
@test_approx_eq_eps(Xb_out, Xb, my_eps)
@test_approx_eq_eps(sumsq_out, ssq, my_eps)
@test_approx_eq_eps(res_out, r, my_eps)

max_stars_repo_name: klkeys/PLINK

max_issues_repo_name: klkeys/PLINK


text:
function test_lmp(scenario)
    @unpack tx, rx, bfield, species, ground = scenario
    waveguide = HomogeneousWaveguide(bfield, species, ground)

    E, amplitude, phase = propagate(waveguide, tx, rx)

    lmp_d = collect(distance(rx, tx)/1e3)
    lmp_a = amplitude
    lmp_p = rad2deg.(phase)

    return lmp_d, lmp_a, lmp_p
end

function test_lmp_segmented(scenario)
    @unpack tx, rx, bfield, species, ground, distances = scenario

    waveguide = SegmentedWaveguide([HomogeneousWaveguide(bfield[i], species[i], ground[i],
                                                         distances[i]) for i in 1:2])

    E, amplitude, phase = propagate(waveguide, tx, rx)

    lmp_d = collect(distance(rx, tx)/1e3)
    lmp_a = amplitude
    lmp_p = rad2deg.(phase)

    return lmp_d, lmp_a, lmp_p
end

function compare(lwpc_file, lmp_d, lmp_a, lmp_p)
    lwpc_d, lwpc_a, lwpc_p = readlog(lwpc_file)

    distmask = lwpc_d .> 300
    lm_p = mod.(lmp_p[distmask], 360)  # modulo because we're not interested in wrapping
    lw_p = mod.(lwpc_p[distmask], 360)

    @test lmp_d ≈ lwpc_d
    @test meanabsdiff(lmp_a[distmask], lwpc_a[distmask]) < 0.4
    @test meanabsdiff(lm_p, lw_p) < 4.0
end

@testset "LWPC comparisons" begin
    @info "Comparing to LWPC..."

    files = Dict(verticalB_scenario=>"verticalB.log", resonant_scenario=>"resonant.log",
        nonresonant_scenario=>"nonresonant.log",
        resonant_elevatedrx_scenario=>"resonant_elevatedrx.log",
        resonant_horizontal_scenario=>"resonant_horizontal.log",
        segmented_scenario=>"segmented.log")

    @info "  Running:"
    for scn in (verticalB_scenario, resonant_scenario, nonresonant_scenario,
        resonant_elevatedrx_scenario, resonant_horizontal_scenario,)

        @info "    "*files[scn]
        compare(joinpath(LWPC_PATH, files[scn]), test_lmp(scn)...)
    end
    for scn in (segmented_scenario,)
        @info "    "*files[scn]
        compare(joinpath(LWPC_PATH, files[scn]), test_lmp_segmented(scn)...)
    end
end


text:
# Problem 4 - Project Euler
# http://projecteuler.net/index.php?section=problems&id=4

is_palindromic(x) = string(x) == reverse(string(x))

function search()
    pmax = 0
    for i in 100:999
        for j in i:999
            p = i * j
            if is_palindromic(p) && p > pmax
                pmax = p
            end
        end
    end
    pmax
end

println(search())

max_stars_repo_name: emergent/ProjectEuler

max_issues_repo_name: emergent/ProjectEuler

max_forks_repo_name: emergent/ProjectEuler


text:
######## ISENTROPIC RELATIONS #########

function TemperatureRatio(pratio::PressureRatio,::Type{Isentropic};gas::PerfectGas=DefaultPerfectGas)
    γ = SpecificHeatRatio(gas)
    TemperatureRatio(pratio^((γ-1)/γ))
end

function TemperatureRatio(ρratio::DensityRatio,::Type{Isentropic};gas::PerfectGas=DefaultPerfectGas)
    γ = SpecificHeatRatio(gas)
    TemperatureRatio(ρratio^(γ-1))
end

function PressureRatio(Tratio::TemperatureRatio,::Type{Isentropic};gas::PerfectGas=DefaultPerfectGas)
    γ = SpecificHeatRatio(gas)
    PressureRatio(Tratio^(γ/(γ-1)))
end

function T0OverT(M::MachNumber,::Type{Isentropic};gas::PerfectGas=DefaultPerfectGas)
    γ = SpecificHeatRatio(gas)
    return TemperatureRatio(1+0.5*(γ-1)*M^2)
end

T0OverT(M::MachNumber;gas::PerfectGas=DefaultPerfectGas) = T0OverT(M,Isentropic,gas=gas)

function P0OverP(M::MachNumber,::Type{Isentropic};gas::PerfectGas=DefaultPerfectGas)
    γ = SpecificHeatRatio(gas)
    Tratio = T0OverT(M,Isentropic,gas=gas)
    PressureRatio(Tratio,Isentropic,gas=gas)
end

function DensityRatio(Tratio::TemperatureRatio,::Type{Isentropic};gas::PerfectGas=DefaultPerfectGas)
    γ = SpecificHeatRatio(gas)
    DensityRatio(Tratio^(1/(γ-1)))
end

function ρ0Overρ(M::MachNumber,::Type{Isentropic};gas::PerfectGas=DefaultPerfectGas)
    γ = SpecificHeatRatio(gas)
    Tratio = T0OverT(M,Isentropic,gas=gas)
    DensityRatio(Tratio,Isentropic,gas=gas)
end

# Temperature relations
function StagnationTemperature(T::Temperature,M::MachNumber,::Type{Isentropic};gas::PerfectGas=DefaultPerfectGas)
    StagnationTemperature(T*T0OverT(M,Isentropic,gas=gas))
end

StagnationTemperature(T::Temperature,M::MachNumber;gas::PerfectGas=DefaultPerfectGas) = StagnationTemperature(T,M,Isentropic,gas=gas)

Temperature(T0::StagnationTemperature,M::MachNumber,::Type{Isentropic};gas::PerfectGas=DefaultPerfectGas) = Temperature(T0/T0OverT(M,gas=gas))

Temperature(T0::StagnationTemperature,M::MachNumber;gas::PerfectGas=DefaultPerfectGas) = Temperature(T0,M,Isentropic,gas=gas)

function MachNumber(T_over_T0::TemperatureRatio,::Type{Isentropic};gas::PerfectGas=DefaultPerfectGas)
    value(T_over_T0) <= 1.0 || error("T/T0 must be 1 or smaller. Maybe you need to supply the inverse?")
    γ = SpecificHeatRatio(gas)
    M2 = ((1.0/T_over_T0)-1)*2/(γ-1)
    MachNumber(sqrt(M2))
end

MachNumber(T_over_T0::TemperatureRatio;gas::PerfectGas=DefaultPerfectGas) = MachNumber(T_over_T0,Isentropic,gas=gas)

MachNumber(T::Temperature,T0::StagnationTemperature,::Type{Isentropic};gas::PerfectGas=DefaultPerfectGas) = MachNumber(TemperatureRatio(T/T0),gas=gas)

MachNumber(T::Temperature,T0::StagnationTemperature;gas::PerfectGas=DefaultPerfectGas) = MachNumber(T,T0,Isentropic,gas=gas)

# Pressure relations
function StagnationPressure(p::Pressure,M::MachNumber,::Type{Isentropic};gas::PerfectGas=DefaultPerfectGas)
    StagnationPressure(p*PressureRatio(T0OverT(M,gas=gas),Isentropic,gas=gas))
end
function Pressure(p0::StagnationPressure,M::MachNumber,::Type{Isentropic};gas::PerfectGas=DefaultPerfectGas)
    Pressure(p0/PressureRatio(T0OverT(M,gas=gas),Isentropic,gas=gas))
end

function MachNumber(p_over_p0::PressureRatio,::Type{Isentropic};gas::PerfectGas=DefaultPerfectGas)
    value(p_over_p0) <= 1.0 || error("p/p0 must be 1 or smaller. Maybe you need to supply the inverse?")
    γ = SpecificHeatRatio(gas)
    MachNumber(TemperatureRatio(p_over_p0,Isentropic,gas=gas),gas=gas)
end

MachNumber(p::Pressure,p0::StagnationPressure,::Type{Isentropic};gas::PerfectGas=DefaultPerfectGas) = MachNumber(PressureRatio(p/p0),Isentropic,gas=gas)

# Density relations
function StagnationDensity(ρ::Density,M::MachNumber,::Type{Isentropic};gas::PerfectGas=DefaultPerfectGas)
    StagnationDensity(ρ*DensityRatio(T0OverT(M,gas=gas),Isentropic,gas=gas))
end
function Density(ρ0::StagnationDensity,M::MachNumber,::Type{Isentropic};gas::PerfectGas=DefaultPerfectGas)
    Density(ρ0/DensityRatio(T0OverT(M,gas=gas),Isentropic,gas=gas))
end

function MachNumber(ρ_over_ρ0::DensityRatio,::Type{Isentropic};gas::PerfectGas=DefaultPerfectGas)
    value(ρ_over_ρ0) <= 1.0 || error("ρ/ρ0 must be 1 or smaller. Maybe you need to supply the inverse?")
    γ = SpecificHeatRatio(gas)
    MachNumber(TemperatureRatio(ρ_over_ρ0,Isentropic,gas=gas))
end

MachNumber(ρ::Density,ρ0::StagnationDensity,::Type{Isentropic};gas::PerfectGas=DefaultPerfectGas) = MachNumber(DensityRatio(ρ/ρ0),gas=gas)

# Area-Mach number relation
function AOverAStar(M::MachNumber,::Type{Isentropic};gas::PerfectGas=DefaultPerfectGas)
    γ = SpecificHeatRatio(gas)
    return AreaRatio(1/M*(2/(γ+1)*(T0OverT(M,gas=gas)))^(0.5*(γ+1)/(γ-1)))
end
AOverAStar(M::MachNumber;gas::PerfectGas=DefaultPerfectGas) = AOverAStar(M,Isentropic,gas=gas)

function AStar(A::Area,M::MachNumber,::Type{Isentropic};gas::PerfectGas=DefaultPerfectGas)
    return Area(A/AOverAStar(M,gas=gas))
end
AStar(A::Area,M::MachNumber;gas::PerfectGas=DefaultPerfectGas) = AStar(A,M,Isentropic,gas=gas)

function SubsonicMachNumber(A_over_Astar::AreaRatio,::Type{Isentropic}; gas::PerfectGas=DefaultPerfectGas)
    value(A_over_Astar) >= 1.0 || error("A/A* must be 1 or larger. Maybe you need to supply the inverse?")
    Msub = find_zero(x -> AOverAStar(MachNumber(x),gas=gas)-A_over_Astar,(0,1),order=16)
    return MachNumber(Msub)
end
function SupersonicMachNumber(A_over_Astar::AreaRatio,::Type{Isentropic}; gas::PerfectGas=DefaultPerfectGas)
  value(A_over_Astar) >= 1.0 || error("A/A* must be 1 or larger. Maybe you need to supply the inverse?")
  Msup = find_zero(x -> AOverAStar(MachNumber(x),gas=gas)-A_over_Astar,(1,Inf),order=16)
  return MachNumber(Msup)
end

function SubsonicMachNumber(A::Area,Aloc::Area,ploc_over_p0::PressureRatio,::Type{Isentropic}; gas::PerfectGas=DefaultPerfectGas)
    Mloc = MachNumber(ploc_over_p0,Isentropic,gas=gas)
    Astar = AStar(Aloc,Mloc,gas=gas)
    SubsonicMachNumber(A,Astar,Isentropic,gas=gas)
end

function SupersonicMachNumber(A::Area,Aloc::Area,ploc_over_p0::PressureRatio,::Type{Isentropic}; gas::PerfectGas=DefaultPerfectGas)
    Mloc = MachNumber(ploc_over_p0,Isentropic,gas=gas)
    Astar = AStar(Aloc,Mloc,gas=gas)
    SupersonicMachNumber(A,Astar,Isentropic,gas=gas)
end

SubsonicMachNumber(A::Area,Astar::Area, ::Type{Isentropic}; gas::PerfectGas=DefaultPerfectGas) =
    SubsonicMachNumber(AreaRatio(A/Astar),Isentropic,gas=gas)

SupersonicMachNumber(A::Area,Astar::Area, ::Type{Isentropic}; gas::PerfectGas=DefaultPerfectGas) =
    SupersonicMachNumber(AreaRatio(A/Astar),Isentropic,gas=gas)

function MachNumber(A_over_Astar::AreaRatio,::Type{Isentropic};gas::PerfectGas=DefaultPerfectGas)
    return SubsonicMachNumber(A_over_Astar,Isentropic,gas=gas),
           SupersonicMachNumber(A_over_Astar,Isentropic,gas=gas)
end

MachNumber(A_over_Astar::AreaRatio;gas::PerfectGas=DefaultPerfectGas) = MachNumber(A_over_Astar,Isentropic,gas=gas)

function MachNumber(M1::MachNumber,A1::Area,A2::Area,::Type{Isentropic};gas::PerfectGas=DefaultPerfectGas)
    M2sub, M2sup = MachNumber(AreaRatio(A2/AStar(A1,M1,gas=gas)),gas=gas)
end
MachNumber(M1::MachNumber,A1::Area,A2::Area;gas::PerfectGas=DefaultPerfectGas) = MachNumber(M1,A2,A2,Isentropic,gas=gas)

function SubsonicPOverP0(A::Area,Aloc::Area,ploc_over_p0::PressureRatio,::Type{Isentropic};gas::PerfectGas=DefaultPerfectGas)
  Mloc = MachNumber(ploc_over_p0,Isentropic,gas=gas)
  Astar = AStar(Aloc,Mloc,gas=gas)
  SubsonicPOverP0(A,Astar,Isentropic,gas=gas)
end

function SupersonicPOverP0(A::Area,Aloc::Area,ploc_over_p0::PressureRatio,::Type{Isentropic};gas::PerfectGas=DefaultPerfectGas)
  Mloc = MachNumber(ploc_over_p0,Isentropic,gas=gas)
  Astar = AStar(Aloc,Mloc,gas=gas)
  SupersonicPOverP0(A,Astar,Isentropic,gas=gas)
end

function SubsonicPOverP0(A::Area,Astar::Area,::Type{Isentropic};gas::PerfectGas=DefaultPerfectGas)
  M = SubsonicMachNumber(A,Astar,Isentropic,gas=gas)
  PressureRatio(1/P0OverP(M,Isentropic,gas=gas))
end

function SupersonicPOverP0(A::Area,Astar::Area,::Type{Isentropic};gas::PerfectGas=DefaultPerfectGas)
  M = SupersonicMachNumber(A,Astar,Isentropic,gas=gas)
  PressureRatio(1/P0OverP(M,Isentropic,gas=gas))
end


text:
const INTERNALNAMES = (:__model__, :__context__, :__varinfo__)

"""
    isassumption(expr[, vn])

Return an expression that can be evaluated to check if `expr` is an assumption in the
model.

Let `expr` be `:(x[1])`. It is an assumption in the following cases:
    1. `x` is not among the input data to the model,
    2. `x` is among the input data to the model but with a value `missing`, or
    3. `x` is among the input data to the model with a value other than missing,
       but `x[1] === missing`.

When `expr` is not an expression or symbol (i.e., a literal), this expands to `false`.

If `vn` is specified, it will be assumed to refer to a expression which
evaluates to a `VarName`, and this will be used in the subsequent checks.
If `vn` is not specified, `AbstractPPL.drop_escape(varname(expr))` will be
used in its place.
"""
function isassumption(expr::Union{Expr,Symbol}, vn=AbstractPPL.drop_escape(varname(expr)))
    return quote
        if $(DynamicPPL.contextual_isassumption)(__context__, $vn)
            # Considered an assumption by `__context__` which means either:
            # 1. We hit the default implementation, e.g. using `DefaultContext`,
            #    which in turn means that we haven't considered if it's one of
            #    the model arguments, hence we need to check this.
            # 2. We are working with a `ConditionContext` _and_ it's NOT in the model arguments,
            #    i.e. we're trying to condition one of the latent variables.
            #    In this case, the below will return `true` since the first branch
            #    will be hit.
            # 3. We are working with a `ConditionContext` _and_ it's in the model arguments,
            #    i.e. we're trying to override the value. This is currently NOT supported.
            #    TODO: Support by adding context to model, and use `model.args`
            #    as the default conditioning. Then we no longer need to check `inargnames`
            #    since it will all be handled by `contextual_isassumption`.
            if !($(DynamicPPL.inargnames)($vn, __model__)) ||
                $(DynamicPPL.inmissings)($vn, __model__)
                true
            else
                $(maybe_view(expr)) === missing
            end
        else
            false
        end
    end
end

# failsafe: a literal is never an assumption
isassumption(expr, vn) = :(false)
isassumption(expr) = :(false)

"""
    contextual_isassumption(context, vn)

Return `true` if `vn` is considered an assumption by `context`.

The default implementation for `AbstractContext` always returns `true`.
"""
contextual_isassumption(::IsLeaf, context, vn) = true
function contextual_isassumption(::IsParent, context, vn)
    return contextual_isassumption(childcontext(context), vn)
end
function contextual_isassumption(context::AbstractContext, vn)
    return contextual_isassumption(NodeTrait(context), context, vn)
end
function contextual_isassumption(context::ConditionContext, vn)
    if hasvalue(context, vn)
        val = getvalue(context, vn)
        # TODO: Do we even need the `>: Missing`, i.e. does it even help the compiler?
        if eltype(val) >: Missing && val === missing
            return true
        else
            return false
        end
    end

    # We might have nested contexts, e.g. `ContextionContext{.., <:PrefixContext{..., <:ConditionContext}}`
    # so we defer to `childcontext` if we haven't concluded that anything yet.
    return contextual_isassumption(childcontext(context), vn)
end
function contextual_isassumption(context::PrefixContext, vn)
    return contextual_isassumption(childcontext(context), prefix(context, vn))
end

# If we're working with, say, a `Symbol`, then we're not going to `view`.
maybe_view(x) = x
maybe_view(x::Expr) = :(@views($x))

"""
    isliteral(expr)

Return `true` if `expr` is a literal, e.g. `1.0` or `[1.0, ]`, and `false` otherwise.
"""
isliteral(e) = false
isliteral(::Number) = true
isliteral(e::Expr) = !isempty(e.args) && all(isliteral, e.args)

"""
    check_tilde_rhs(x)

Check if the right-hand side `x` of a `~` is a `Distribution` or an array of
`Distributions`, then return `x`.
"""
function check_tilde_rhs(@nospecialize(x))
    return throw(
        ArgumentError(
            "the right-hand side of a `~` must be a `Distribution` or an array of `Distribution`s",
        ),
    )
end
check_tilde_rhs(x::Distribution) = x
check_tilde_rhs(x::AbstractArray{<:Distribution}) = x

"""
    unwrap_right_vn(right, vn)

Return the unwrapped distribution on the right-hand side and variable name on the left-hand
side of a `~` expression such as `x ~ Normal()`.

This is used mainly to unwrap `NamedDist` distributions.
"""
unwrap_right_vn(right, vn) = right, vn
unwrap_right_vn(right::NamedDist, vn) = unwrap_right_vn(right.dist, right.name)

"""
    unwrap_right_left_vns(right, left, vns)

Return the unwrapped distributions on the right-hand side and values and variable names on the
left-hand side of a `.~` expression such as `x .~ Normal()`.

This is used mainly to unwrap `NamedDist` distributions and adjust the indices of the
variables.

# Example
```jldoctest; setup=:(using Distributions, LinearAlgebra)
julia> _, _, vns = DynamicPPL.unwrap_right_left_vns(MvNormal(ones(2), I), randn(2, 2), @varname(x)); vns[end]
x[:,2]

julia> _, _, vns = DynamicPPL.unwrap_right_left_vns(Normal(), randn(1, 2), @varname(x)); vns[end]
x[1,2]

julia> _, _, vns = DynamicPPL.unwrap_right_left_vns(Normal(), randn(1, 2), @varname(x[:])); vns[end]
x[:][1,2]

julia> _, _, vns = DynamicPPL.unwrap_right_left_vns(Normal(), randn(3), @varname(x[1])); vns[end]
x[1][3]
```
"""
unwrap_right_left_vns(right, left, vns) = right, left, vns
function unwrap_right_left_vns(right::NamedDist, left, vns)
    return unwrap_right_left_vns(right.dist, left, right.name)
end
function unwrap_right_left_vns(
    right::MultivariateDistribution, left::AbstractMatrix, vn::VarName
)
    # This an expression such as `x .~ MvNormal()` which we interpret as
    #     x[:, i] ~ MvNormal()
    # for `i = size(left, 2)`. Hence the symbol should be `x[:, i]`,
    # and we therefore add the `Colon()` below.
    vns = map(axes(left, 2)) do i
        return vn ∘ Setfield.IndexLens((Colon(), i))
    end
    return unwrap_right_left_vns(right, left, vns)
end
function unwrap_right_left_vns(
    right::Union{Distribution,AbstractArray{<:Distribution}},
    left::AbstractArray,
    vn::VarName,
)
    vns = map(CartesianIndices(left)) do i
        return vn ∘ Setfield.IndexLens(Tuple(i))
    end
    return unwrap_right_left_vns(right, left, vns)
end

#################
# Main Compiler #
#################

"""
    @model(expr[, warn = false])

Macro to specify a probabilistic model.

If `warn` is `true`, a warning is displayed if internal variable names are used in the model
definition.

# Examples

Model definition:

```julia
@model function model(x, y = 42)
    ...
end
```

To generate a `Model`, call `model(xvalue)` or `model(xvalue, yvalue)`.
"""
macro model(expr, warn=false)
    # include `LineNumberNode` with information about the call site in the
    # generated function for easier debugging and interpretation of error messages
    return esc(model(__module__, __source__, expr, warn))
end

function model(mod, linenumbernode, expr, warn)
    modelinfo = build_model_info(expr)

    # Generate main body
    modelinfo[:body] = generate_mainbody(mod, modelinfo[:modeldef][:body], warn)

    return build_output(modelinfo, linenumbernode)
end

"""
    build_model_info(input_expr)

Builds the `model_info` dictionary from the model's expression.
"""
function build_model_info(input_expr)
    # Break up the model definition and extract its name, arguments, and function body
    modeldef = MacroTools.splitdef(input_expr)

    # Check that the function has a name
    # https://github.com/TuringLang/DynamicPPL.jl/issues/260
    haskey(modeldef, :name) ||
        throw(ArgumentError("anonymous functions without name are not supported"))

    # Print a warning if function body of the model is empty
    warn_empty(modeldef[:body])

    ## Construct model_info dictionary

    # Shortcut if the model does not have any arguments
    if !haskey(modeldef, :args) && !haskey(modeldef, :kwargs)
        modelinfo = Dict(
            :allargs_exprs => [],
            :allargs_syms => [],
            :allargs_namedtuple => NamedTuple(),
            :defaults_namedtuple => NamedTuple(),
            :modeldef => modeldef,
        )
        return modelinfo
    end

    # Ensure that all arguments have a name, i.e., are of the form `name` or `name::T`
    addargnames!(modeldef[:args])

    # Extract the positional and keyword arguments from the model definition.
    allargs = vcat(modeldef[:args], modeldef[:kwargs])

    # Split the argument expressions and the default values.
    allargs_exprs_defaults = map(allargs) do arg
        MacroTools.@match arg begin
            (x_ = val_) => (x, val)
            x_ => (x, NO_DEFAULT)
        end
    end

    # Extract the expressions of the arguments, without default values.
    allargs_exprs = first.(allargs_exprs_defaults)

    # Extract the names of the arguments.
    allargs_syms = map(allargs_exprs) do arg
        MacroTools.@match arg begin
            (name_::_) => name
            x_ => x
        end
    end

    # Build named tuple expression of the argument symbols and variables of the same name.
    allargs_namedtuple = to_namedtuple_expr(allargs_syms)

    # Extract default values of the positional and keyword arguments.
    default_syms = []
    default_vals = []
    for (sym, (expr, val)) in zip(allargs_syms, allargs_exprs_defaults)
        if val !== NO_DEFAULT
            push!(default_syms, sym)
            push!(default_vals, val)
        end
    end

    # Build named tuple expression of the argument symbols with default values.
    defaults_namedtuple = to_namedtuple_expr(default_syms, default_vals)

    modelinfo = Dict(
        :allargs_exprs => allargs_exprs,
        :allargs_syms => allargs_syms,
        :allargs_namedtuple => allargs_namedtuple,
        :defaults_namedtuple => defaults_namedtuple,
        :modeldef => modeldef,
    )

    return modelinfo
end

"""
    generate_mainbody(mod, expr, warn)

Generate the body of the main evaluation function from expression `expr` and arguments
`args`.

If `warn` is true, a warning is displayed if internal variables are used in the model
definition.
"""
generate_mainbody(mod, expr, warn) = generate_mainbody!(mod, Symbol[], expr, warn)

generate_mainbody!(mod, found, x, warn) = x
function generate_mainbody!(mod, found, sym::Symbol, warn)
    if warn && sym in INTERNALNAMES && sym ∉ found
        @warn "you are using the internal variable `$sym`"
        push!(found, sym)
    end

    return sym
end
function generate_mainbody!(mod, found, expr::Expr, warn)
    # Do not touch interpolated expressions
    expr.head === :$ && return expr.args[1]

    # Do we don't want escaped expressions because we unfortunately
    # escape the entire body afterwards.
    Meta.isexpr(expr, :escape) && return generate_mainbody(mod, found, expr.args[1], warn)

    # If it's a macro, we expand it
    if Meta.isexpr(expr, :macrocall)
        return generate_mainbody!(mod, found, macroexpand(mod, expr; recursive=true), warn)
    end

    # Modify dotted tilde operators.
    args_dottilde = getargs_dottilde(expr)
    if args_dottilde !== nothing
        L, R = args_dottilde
        return Base.remove_linenums!(
            generate_dot_tilde(
                generate_mainbody!(mod, found, L, warn),
                generate_mainbody!(mod, found, R, warn),
            ),
        )
    end

    # Modify tilde operators.
    args_tilde = getargs_tilde(expr)
    if args_tilde !== nothing
        L, R = args_tilde
        return Base.remove_linenums!(
            generate_tilde(
                generate_mainbody!(mod, found, L, warn),
                generate_mainbody!(mod, found, R, warn),
            ),
        )
    end

    return Expr(expr.head, map(x -> generate_mainbody!(mod, found, x, warn), expr.args)...)
end

function generate_tilde_literal(left, right)
    # If the LHS is a literal, it is always an observation
    @gensym value
    return quote
        $value, __varinfo__ = $(DynamicPPL.tilde_observe!!)(
            __context__, $(DynamicPPL.check_tilde_rhs)($right), $left, __varinfo__
        )
        $value
    end
end

"""
    generate_tilde(left, right)

Generate an `observe` expression for data variables and `assume` expression for parameter
variables.
"""
function generate_tilde(left, right)
    isliteral(left) && return generate_tilde_literal(left, right)

    # Otherwise it is determined by the model or its value,
    # if the LHS represents an observation
    @gensym vn isassumption value

    # HACK: Usage of `drop_escape` is unfortunate. It's a consequence of the fact
    # that in DynamicPPL we the entire function body. Instead we should be
    # more selective with our escape. Until that's the case, we remove them all.
    return quote
        $vn = $(AbstractPPL.drop_escape(varname(left)))
        $isassumption = $(DynamicPPL.isassumption(left, vn))
        if $isassumption
            $(generate_tilde_assume(left, right, vn))
        else
            # If `vn` is not in `argnames`, we need to make sure that the variable is defined.
            if !$(DynamicPPL.inargnames)($vn, __model__)
                $left = $(DynamicPPL.getvalue_nested)(__context__, $vn)
            end

            $value, __varinfo__ = $(DynamicPPL.tilde_observe!!)(
                __context__,
                $(DynamicPPL.check_tilde_rhs)($right),
                $(maybe_view(left)),
                $vn,
                __varinfo__,
            )
            $value
        end
    end
end

function generate_tilde_assume(left, right, vn)
    # HACK: Because the Setfield.jl macro does not support assignment
    # with multiple arguments on the LHS, we need to capture the return-values
    # and then update the LHS variables one by one.
    @gensym value
    expr = :($left = $value)
    if left isa Expr
        expr = AbstractPPL.drop_escape(
            Setfield.setmacro(BangBang.prefermutation, expr; overwrite=true)
        )
    end

    return quote
        $value, __varinfo__ = $(DynamicPPL.tilde_assume!!)(
            __context__,
            $(DynamicPPL.unwrap_right_vn)($(DynamicPPL.check_tilde_rhs)($right), $vn)...,
            __varinfo__,
        )
        $expr
        $value
    end
end

"""
    generate_dot_tilde(left, right)

Generate the expression that replaces `left .~ right` in the model body.
"""
function generate_dot_tilde(left, right)
    isliteral(left) && return generate_tilde_literal(left, right)

    # Otherwise it is determined by the model or its value,
    # if the LHS represents an observation
    @gensym vn isassumption value
    return quote
        $vn = $(AbstractPPL.drop_escape(varname(left)))
        $isassumption = $(DynamicPPL.isassumption(left, vn))
        if $isassumption
            $(generate_dot_tilde_assume(left, right, vn))
        else
            # If `vn` is not in `argnames`, we need to make sure that the variable is defined.
            if !$(DynamicPPL.inargnames)($vn, __model__)
                $left .= $(DynamicPPL.getvalue_nested)(__context__, $vn)
            end

            $value, __varinfo__ = $(DynamicPPL.dot_tilde_observe!!)(
                __context__,
                $(DynamicPPL.check_tilde_rhs)($right),
                $(maybe_view(left)),
                $vn,
                __varinfo__,
            )
            $value
        end
    end
end

function generate_dot_tilde_assume(left, right, vn)
    # We don't need to use `Setfield.@set` here since
    # `.=` is always going to be inplace + needs `left` to
    # be something that supports `.=`.
    @gensym value
    return quote
        $value, __varinfo__ = $(DynamicPPL.dot_tilde_assume!!)(
            __context__,
            $(DynamicPPL.unwrap_right_left_vns)(
                $(DynamicPPL.check_tilde_rhs)($right), $(maybe_view(left)), $vn
            )...,
            __varinfo__,
        )
        $left .= $value
        $value
    end
end

# Note that we cannot use `MacroTools.isdef` because
# of https://github.com/FluxML/MacroTools.jl/issues/154.
"""
    isfuncdef(expr)

Return `true` if `expr` is any form of function definition, and `false` otherwise.
"""
function isfuncdef(e::Expr)
    return if Meta.isexpr(e, :function)
        # Classic `function f(...)`
        true
    elseif Meta.isexpr(e, :->)
        # Anonymous functions/lambdas, e.g. `do` blocks or `->` defs.
        true
    elseif Meta.isexpr(e, :(=)) && Meta.isexpr(e.args[1], :call)
        # Short function defs, e.g. `f(args...) = ...`.
        true
    else
        false
    end
end

"""
    replace_returns(expr)

Return `Expr` with all `return ...` statements replaced with
`return ..., DynamicPPL.return_values(__varinfo__)`.

Note that this method will _not_ replace `return` statements within function
definitions. This is checked using [`isfuncdef`](@ref).
"""
replace_returns(e) = e
function replace_returns(e::Expr)
    if isfuncdef(e)
        return e
    end

    if Meta.isexpr(e, :return)
        # We capture the original return-value in `retval` and return
        # a `Tuple{typeof(retval),typeof(__varinfo__)}`.
        # If we don't capture the return-value separately, cases such as
        # `return x = 1` will result in `(x = 1, __varinfo__)` which will
        # mistakenly attempt to construct a `NamedTuple` (which fails on Julia 1.3
        # and is not our intent).
        @gensym retval
        return quote
            $retval = $(e.args...)
            return $retval, __varinfo__
        end
    end

    return Expr(e.head, map(replace_returns, e.args)...)
end

# If it's just a symbol, e.g. `f(x) = 1`, then we make it `f(x) = return 1`.
make_returns_explicit!(body) = Expr(:return, body)
function make_returns_explicit!(body::Expr)
    # If the last statement is a return-statement, we don't do anything.
    # Otherwise we replace the last statement with a `return` statement.
    if !Meta.isexpr(body.args[end], :return)
        body.args[end] = Expr(:return, body.args[end])
    end
    return body
end

const FloatOrArrayType = Type{<:Union{AbstractFloat,AbstractArray}}
hasmissing(::Type) = false
hasmissing(::Type{>:Missing}) = true
hasmissing(::Type{<:AbstractArray{TA}}) where {TA} = hasmissing(TA)
hasmissing(::Type{Union{}}) = false # issue #368

"""
    build_output(modelinfo, linenumbernode)

Builds the output expression.
"""
function build_output(modelinfo, linenumbernode)
    ## Build the anonymous evaluator from the user-provided model definition.
    evaluatordef = deepcopy(modelinfo[:modeldef])

    # Add the internal arguments to the user-specified arguments (positional + keywords).
    evaluatordef[:args] = vcat(
        [
            :(__model__::$(DynamicPPL.Model)),
            :(__varinfo__::$(DynamicPPL.AbstractVarInfo)),
            :(__context__::$(DynamicPPL.AbstractContext)),
        ],
        modelinfo[:allargs_exprs],
    )

    # Delete the keyword arguments.
    evaluatordef[:kwargs] = []

    # Replace the user-provided function body with the version created by DynamicPPL.
    # We use `MacroTools.@q begin ... end` instead of regular `quote ... end` to ensure
    # that no new `LineNumberNode`s are added apart from the reference `linenumbernode`
    # to the call site.
    # NOTE: We need to replace statements of the form `return ...` with
    # `return (..., __varinfo__)` to ensure that the second
    # element in the returned value is always the most up-to-date `__varinfo__`.
    # See the docstrings of `replace_returns` for more info.
    evaluatordef[:body] = MacroTools.@q begin
        $(linenumbernode)
        $(replace_returns(make_returns_explicit!(modelinfo[:body])))
    end

    ## Build the model function.

    # Extract the named tuple expression of all arguments and the default values.
    allargs_namedtuple = modelinfo[:allargs_namedtuple]
    defaults_namedtuple = modelinfo[:defaults_namedtuple]

    # Obtain or generate the name of the model to support functors:
    # https://github.com/TuringLang/DynamicPPL.jl/issues/367
    modeldef = modelinfo[:modeldef]
    if MacroTools.@capture(modeldef[:name], ::T_)
        name = gensym(:f)
        modeldef[:name] = Expr(:(::), name, T)
    elseif MacroTools.@capture(modeldef[:name], (name_::_ | name_))
    else
        throw(ArgumentError("unsupported format of model function"))
    end

    # Update the function body of the user-specified model.
    # We use `MacroTools.@q begin ... end` instead of regular `quote ... end` to ensure
    # that no new `LineNumberNode`s are added apart from the reference `linenumbernode`
    # to the call site
    modeldef[:body] = MacroTools.@q begin
        $(linenumbernode)
        return $(DynamicPPL.Model)($name, $allargs_namedtuple, $defaults_namedtuple)
    end

    return MacroTools.@q begin
        $(MacroTools.combinedef(evaluatordef))
        $(Base).@__doc__ $(MacroTools.combinedef(modeldef))
    end
end

function warn_empty(body)
    if all(l -> isa(l, LineNumberNode), body.args)
        @warn("Model definition seems empty, still continue.")
    end
    return nothing
end

"""
    matchingvalue(sampler, vi, value)
    matchingvalue(context::AbstractContext, vi, value)

Convert the `value` to the correct type for the `sampler` or `context` and the `vi` object.

For a `context` that is _not_ a `SamplingContext`, we fall back to
`matchingvalue(SampleFromPrior(), vi, value)`.
"""
function matchingvalue(sampler, vi, value)
    T = typeof(value)
    if hasmissing(T)
        _value = convert(get_matching_type(sampler, vi, T), value)
        if _value === value
            return deepcopy(_value)
        else
            return _value
        end
    else
        return value
    end
end
function matchingvalue(sampler::AbstractSampler, vi, value::FloatOrArrayType)
    return get_matching_type(sampler, vi, value)
end

function matchingvalue(context::AbstractContext, vi, value)
    return matchingvalue(NodeTrait(matchingvalue, context), context, vi, value)
end
function matchingvalue(::IsLeaf, context::AbstractContext, vi, value)
    return matchingvalue(SampleFromPrior(), vi, value)
end
function matchingvalue(::IsParent, context::AbstractContext, vi, value)
    return matchingvalue(childcontext(context), vi, value)
end
function matchingvalue(context::SamplingContext, vi, value)
    return matchingvalue(context.sampler, vi, value)
end

"""
    get_matching_type(spl::AbstractSampler, vi, ::Type{T}) where {T}

Get the specialized version of type `T` for sampler `spl`.

For example, if `T === Float64` and `spl::Hamiltonian`, the matching type is
`eltype(vi[spl])`.
"""
get_matching_type(spl::AbstractSampler, vi, ::Type{T}) where {T} = T
function get_matching_type(spl::AbstractSampler, vi, ::Type{<:Union{Missing,AbstractFloat}})
    return Union{Missing,floatof(eltype(vi, spl))}
end
function get_matching_type(spl::AbstractSampler, vi, ::Type{<:AbstractFloat})
    return floatof(eltype(vi, spl))
end
function get_matching_type(spl::AbstractSampler, vi, ::Type{<:Array{T,N}}) where {T,N}
    return Array{get_matching_type(spl, vi, T),N}
end
function get_matching_type(spl::AbstractSampler, vi, ::Type{<:Array{T}}) where {T}
    return Array{get_matching_type(spl, vi, T)}
end

floatof(::Type{T}) where {T<:Real} = typeof(one(T) / one(T))
floatof(::Type) = Real # fallback if type inference failed

max_stars_repo_name: TuringLang/MicroPPL

max_issues_repo_name: TuringLang/MicroPPL

max_forks_repo_name: TuringLang/MicroPPL


text:
function iterate!(job::BasicMCJob, ::Type{MALA}, ::Type{Univariate})
  if (isa(job.tuner, VanillaMCTuner) && job.tuner.verbose) || isa(job.tuner, AcceptanceRateMCTuner)
    job.sstate.tune.proposed += 1
  end

  job.sstate.μ = job.pstate.value+0.5*job.sstate.tune.step*job.pstate.gradlogtarget
  job.sstate.pstate.value = job.sstate.μ+sqrt(job.sstate.tune.step)*randn()

  job.parameter.uptogradlogtarget!(job.sstate.pstate)

  job.sstate.ratio = job.sstate.pstate.logtarget-job.pstate.logtarget

  job.sstate.ratio += 0.5*(abs2(job.sstate.μ-job.sstate.pstate.value)/job.sstate.tune.step)
  job.sstate.μ = job.sstate.pstate.value+0.5*job.sstate.tune.step*job.sstate.pstate.gradlogtarget
  job.sstate.ratio -= 0.5*(abs2(job.sstate.μ-job.pstate.value)/job.sstate.tune.step)

  if (job.sstate.ratio > 0 || (job.sstate.ratio > log(rand())))
    job.pstate.value = job.sstate.pstate.value

    job.pstate.gradlogtarget = job.sstate.pstate.gradlogtarget
    if in(:gradloglikelihood, job.outopts[:monitor]) && job.parameter.gradloglikelihood! != nothing
      job.pstate.gradloglikelihood = job.sstate.pstate.gradloglikelihood
    end
    if in(:gradlogprior, job.outopts[:monitor]) && job.parameter.gradlogprior! != nothing
      job.pstate.gradlogprior = job.sstate.pstate.gradlogprior
    end

    job.pstate.logtarget = job.sstate.pstate.logtarget
    if in(:loglikelihood, job.outopts[:monitor]) && job.parameter.loglikelihood! != nothing
      job.pstate.loglikelihood = job.sstate.pstate.loglikelihood
    end
    if in(:logprior, job.outopts[:monitor]) && job.parameter.logprior! != nothing
      job.pstate.logprior = job.sstate.pstate.logprior
    end

    if !isempty(job.sstate.diagnosticindices)
      if haskey(job.sstate.diagnosticindices, :accept)
        job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:accept]] = true
      end
    end

    if (isa(job.tuner, VanillaMCTuner) && job.tuner.verbose) || isa(job.tuner, AcceptanceRateMCTuner)
      job.sstate.tune.accepted += 1
    end
  else
    if !isempty(job.sstate.diagnosticindices)
      if haskey(job.sstate.diagnosticindices, :accept)
        job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:accept]] = false
      end
    end
  end

  if (isa(job.tuner, VanillaMCTuner) && job.tuner.verbose) || isa(job.tuner, AcceptanceRateMCTuner)
    if (job.sstate.tune.totproposed <= job.range.burnin && mod(job.sstate.tune.proposed, job.tuner.period) == 0)
      rate!(job.sstate.tune)

      if isa(job.tuner, AcceptanceRateMCTuner)
        tune!(job.sstate.tune, job.tuner)
      end

      if job.tuner.verbose
        println(
          "Burnin iteration ",
          job.fmt_iter(job.sstate.tune.totproposed),
          " of ",
          job.range.burnin,
          ": ",
          job.fmt_perc(100*job.sstate.tune.rate),
          " % acceptance rate"
        )
      end

      reset_burnin!(job.sstate.tune)
    end
  end
end

function iterate!(job::BasicMCJob, ::Type{MALA}, ::Type{Multivariate})
  if (isa(job.tuner, VanillaMCTuner) && job.tuner.verbose) || isa(job.tuner, AcceptanceRateMCTuner)
    job.sstate.tune.proposed += 1
  end

  job.sstate.μ[:] = job.pstate.value+0.5*job.sstate.tune.step*job.pstate.gradlogtarget
  job.sstate.pstate.value[:] = job.sstate.μ+sqrt(job.sstate.tune.step)*randn(job.pstate.size)

  job.parameter.uptogradlogtarget!(job.sstate.pstate)

  job.sstate.ratio = job.sstate.pstate.logtarget-job.pstate.logtarget

  job.sstate.ratio += sum(0.5*(abs2.(job.sstate.μ-job.sstate.pstate.value)/job.sstate.tune.step))
  job.sstate.μ[:] = job.sstate.pstate.value+0.5*job.sstate.tune.step*job.sstate.pstate.gradlogtarget
  job.sstate.ratio -= sum(0.5*(abs2.(job.sstate.μ-job.pstate.value)/job.sstate.tune.step))

  if (job.sstate.ratio > 0 || (job.sstate.ratio > log(rand())))
    job.pstate.value = copy(job.sstate.pstate.value)

    job.pstate.gradlogtarget = copy(job.sstate.pstate.gradlogtarget)
    if in(:gradloglikelihood, job.outopts[:monitor]) && job.parameter.gradloglikelihood! != nothing
      job.pstate.gradloglikelihood = copy(job.sstate.pstate.gradloglikelihood)
    end
    if in(:gradlogprior, job.outopts[:monitor]) && job.parameter.gradlogprior! != nothing
      job.pstate.gradlogprior = copy(job.sstate.pstate.gradlogprior)
    end

    job.pstate.logtarget = job.sstate.pstate.logtarget
    if in(:loglikelihood, job.outopts[:monitor]) && job.parameter.loglikelihood! != nothing
      job.pstate.loglikelihood = job.sstate.pstate.loglikelihood
    end
    if in(:logprior, job.outopts[:monitor]) && job.parameter.logprior! != nothing
      job.pstate.logprior = job.sstate.pstate.logprior
    end

    if !isempty(job.sstate.diagnosticindices)
      if haskey(job.sstate.diagnosticindices, :accept)
        job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:accept]] = true
      end
    end

    if (isa(job.tuner, VanillaMCTuner) && job.tuner.verbose) || isa(job.tuner, AcceptanceRateMCTuner)
      job.sstate.tune.accepted += 1
    end
  else
    if !isempty(job.sstate.diagnosticindices)
      if haskey(job.sstate.diagnosticindices, :accept)
        job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:accept]] = false
      end
    end
  end

  if (isa(job.tuner, VanillaMCTuner) && job.tuner.verbose) || isa(job.tuner, AcceptanceRateMCTuner)
    if (job.sstate.tune.totproposed <= job.range.burnin && mod(job.sstate.tune.proposed, job.tuner.period) == 0)
      rate!(job.sstate.tune)

      if isa(job.tuner, AcceptanceRateMCTuner)
        tune!(job.sstate.tune, job.tuner)
      end

      if job.tuner.verbose
        println(
          "Burnin iteration ",
          job.fmt_iter(job.sstate.tune.totproposed),
          " of ",
          job.range.burnin,
          ": ",
          job.fmt_perc(100*job.sstate.tune.rate),
          " % acceptance rate"
        )
      end

      reset_burnin!(job.sstate.tune)
    end
  end
end


text:
mutable struct Translational{T,N} <: Joint{T,N}
    V3::Adjoint{T,SVector{3,T}} # in body1's frame
    V12::SMatrix{2,3,T,6} # in body1's frame
    vertices::NTuple{2,SVector{3,T}} # in body1's & body2's frames

    spring::T
    damper::T

    Fτ::SVector{3,T}

    childid::Int64

    function Translational{T,N}(body1::AbstractBody, body2::AbstractBody;
            p1::AbstractVector = szeros(T,3), p2::AbstractVector = szeros(T,3), axis::AbstractVector = szeros(T,3), spring = zero(T), damper = zero(T)
        ) where {T,N}

        vertices = (p1, p2)
        V1, V2, V3 = orthogonalrows(axis)
        V12 = [V1;V2]

        Fτ = zeros(T,3)

        childid = body2.id

        new{T,N}(V3, V12, vertices, spring, damper, Fτ, childid), body1.id, body2.id
    end
end

Translational0 = Translational{T,0} where T
Translational1 = Translational{T,1} where T
Translational2 = Translational{T,2} where T
Translational3 = Translational{T,3} where T

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, constraint::Translational{T,N}) where {T,N}
    summary(io, constraint)
    println(io,"")
    println(io, " V3:       "*string(constraint.V3))
    println(io, " V12:      "*string(constraint.V12))
    println(io, " vertices: "*string(constraint.vertices))
end

### Constraints and derivatives
## Position level constraints (for dynamics)
@inline function g(joint::Translational, xa::AbstractVector, qa::UnitQuaternion, xb::AbstractVector, qb::UnitQuaternion)
    vertices = joint.vertices
    return vrotate(xb + vrotate(vertices[2], qb) - (xa + vrotate(vertices[1], qa)), inv(qa))
end
@inline function g(joint::Translational, xb::AbstractVector, qb::UnitQuaternion)
    vertices = joint.vertices
    return xb + vrotate(vertices[2], qb) - vertices[1]
end

## Derivatives NOT accounting for quaternion specialness
@inline function ∂g∂posa(joint::Translational, xa::AbstractVector, qa::UnitQuaternion, xb::AbstractVector, qb::UnitQuaternion)
    point2 = xb + vrotate(joint.vertices[2], qb)

    X = -VLᵀmat(qa) * RVᵀmat(qa)
    Q = 2 * VLᵀmat(qa) * (Lmat(UnitQuaternion(point2)) - Lmat(UnitQuaternion(xa)))

    return X, Q
end
@inline function ∂g∂posb(joint::Translational, xa::AbstractVector, qa::UnitQuaternion, xb::AbstractVector, qb::UnitQuaternion)
    X = VLᵀmat(qa) * RVᵀmat(qa)
    Q = 2 * VLᵀmat(qa) * Rmat(qa) * Rᵀmat(qb) * Rmat(UnitQuaternion(joint.vertices[2]))

    return X, Q
end
@inline function ∂g∂posb(joint::Translational, xb::AbstractVector, qb::UnitQuaternion)
    X = I
    Q = 2 * VRᵀmat(qb) * Rmat(UnitQuaternion(joint.vertices[2]))

    return X, Q
end

## vec(G) Jacobian (also NOT accounting for quaternion specialness in the second derivative: ∂(∂ʳg∂posx)∂y)
@inline function ∂2g∂posaa(joint::Translational{T}, xa::AbstractVector, qa::UnitQuaternion, xb::AbstractVector, qb::UnitQuaternion) where T
    Lpos = Lmat(UnitQuaternion(xb + vrotate(joint.vertices[2], qb) - xa))
    Ltpos = Lᵀmat(UnitQuaternion(xb + vrotate(joint.vertices[2], qb) - xa))

    XX = szeros(T, 9, 3)
    XQ = -kron(Vmat(T),VLᵀmat(qa))*∂R∂qsplit(T) - kron(VRᵀmat(qa),Vmat(T))*∂Lᵀ∂qsplit(T)
    QX = -kron(VLᵀmat(qa),2*VLᵀmat(qa))*∂L∂qsplit(T)[:,SA[2; 3; 4]]
    QQ = kron(Vmat(T),2*VLᵀmat(qa)*Lpos)*∂L∂qsplit(T) + kron(VLᵀmat(qa)*Ltpos,2*Vmat(T))*∂Lᵀ∂qsplit(T)

    return XX, XQ, QX, QQ
end
@inline function ∂2g∂posab(joint::Translational{T}, xa::AbstractVector, qa::UnitQuaternion, xb::AbstractVector, qb::UnitQuaternion) where T
    XX = szeros(T, 9, 3)
    XQ = szeros(T, 9, 4)
    QX = kron(VLᵀmat(qa),2*VLᵀmat(qa))*∂L∂qsplit(T)[:,SA[2; 3; 4]]
    QQ = kron(VLᵀmat(qa),2*VLᵀmat(qa)*Lmat(qb)*Lmat(UnitQuaternion(joint.vertices[2])))*∂Lᵀ∂qsplit(T) + kron(VLᵀmat(qa)*Lmat(qb)*Lᵀmat(UnitQuaternion(joint.vertices[2])),2*VLᵀmat(qa))*∂L∂qsplit(T)

    return XX, XQ, QX, QQ
end
@inline function ∂2g∂posba(joint::Translational{T}, xa::AbstractVector, qa::UnitQuaternion, xb::AbstractVector, qb::UnitQuaternion) where T
    XX = szeros(T, 9, 3)
    XQ = kron(Vmat(T),VLᵀmat(qa))*∂R∂qsplit(T) + kron(VRᵀmat(qa),Vmat(T))*∂Lᵀ∂qsplit(T)
    QX = szeros(T, 9, 3)
    QQ = kron(VLᵀmat(qb)*Rᵀmat(UnitQuaternion(joint.vertices[2]))*Rmat(qb),2*VLᵀmat(qa))*∂R∂qsplit(T) + kron(VLᵀmat(qb)*Rᵀmat(UnitQuaternion(joint.vertices[2]))*Rmat(qb)*Rᵀmat(qa),2*Vmat(T))*∂Lᵀ∂qsplit(T)

    return XX, XQ, QX, QQ
end
@inline function ∂2g∂posbb(joint::Translational{T}, xa::AbstractVector, qa::UnitQuaternion, xb::AbstractVector, qb::UnitQuaternion) where T
    XX = szeros(T, 9, 3)
    XQ = szeros(T, 9, 4)
    QX = szeros(T, 9, 3)
    QQ = kron(Vmat(T),2*VLᵀmat(qa)*Rmat(qa)*Rᵀmat(qb)*Rmat(UnitQuaternion(joint.vertices[2])))*∂L∂qsplit(T) + kron(VLᵀmat(qb)*Rᵀmat(UnitQuaternion(joint.vertices[2])),2*VLᵀmat(qa)*Rmat(qa))*∂Rᵀ∂qsplit(T)

    return XX, XQ, QX, QQ
end
@inline function ∂2g∂posbb(joint::Translational{T}, xb::AbstractVector, qb::UnitQuaternion) where T
    XX = szeros(T, 9, 3)
    XQ = szeros(T, 9, 4)
    QX = szeros(T, 9, 3)
    QQ = kron(Vmat(T),2*VRᵀmat(qb)*Rmat(UnitQuaternion(joint.vertices[2])))*∂L∂qsplit(T) + kron(VLᵀmat(qb)*Rᵀmat(UnitQuaternion(joint.vertices[2])),2*Vmat(T))*∂Rᵀ∂qsplit(T)

    return XX, XQ, QX, QQ
end

### Spring and damper
## Forces for dynamics
@inline function springforcea(joint::Translational, xa::AbstractVector, qa::UnitQuaternion, xb::AbstractVector, qb::UnitQuaternion)
    A = nullspacemat(joint)
    Aᵀ = zerodimstaticadjoint(A)
    distance = A * g(joint, xa, qa, xb, qb)
    force = Aᵀ * A * joint.spring * Aᵀ * distance  # Currently assumes same spring constant in all directions
    return [force;szeros(3)]
end
@inline function springforceb(joint::Translational, xa::AbstractVector, qa::UnitQuaternion, xb::AbstractVector, qb::UnitQuaternion)
    A = nullspacemat(joint)
    Aᵀ = zerodimstaticadjoint(A)
    distance = A * g(joint, xa, qa, xb, qb)
    force = - Aᵀ * A * joint.spring * Aᵀ * distance  # Currently assumes same spring constant in all directions
    return [force;szeros(3)]
end
@inline function springforceb(joint::Translational, xb::AbstractVector, qb::UnitQuaternion)
    A = nullspacemat(joint)
    Aᵀ = zerodimstaticadjoint(A)
    distance = A * g(joint, xb, qb)
    force = - Aᵀ * A * joint.spring * Aᵀ * distance  # Currently assumes same spring constant in all directions
    return [force;szeros(3)]
end

@inline function damperforcea(joint::Translational, xa::AbstractVector, va::AbstractVector, qa::UnitQuaternion, ωa::AbstractVector,
        xb::AbstractVector, vb::AbstractVector, qb::UnitQuaternion, ωb::AbstractVector)
    A = nullspacemat(joint)
    Aᵀ = zerodimstaticadjoint(A)
    velocity = A * (vb - va)
    force = Aᵀ * A * joint.damper * Aᵀ * velocity  # Currently assumes same damper constant in all directions
    return [force;szeros(3)]
end
@inline function damperforceb(joint::Translational, xa::AbstractVector, va::AbstractVector, qa::UnitQuaternion, ωa::AbstractVector,
    xb::AbstractVector, vb::AbstractVector, qb::UnitQuaternion, ωb::AbstractVector)
    A = nullspacemat(joint)
    Aᵀ = zerodimstaticadjoint(A)
    velocity = A * (vb - va)
    force = - Aᵀ * A * joint.damper * Aᵀ * velocity  # Currently assumes same damper constant in all directions
    return [force;szeros(3)]
end
@inline function damperforceb(joint::Translational, xb::AbstractVector, vb::AbstractVector, qb::UnitQuaternion, ωb::AbstractVector)
    A = nullspacemat(joint)
    Aᵀ = zerodimstaticadjoint(A)
    velocity = A * vb
    force = - Aᵀ * A * joint.damper * Aᵀ * velocity  # Currently assumes same damper constant in all directions
    return [force;szeros(3)]
end

## Damper velocity derivatives
@inline function diagonal∂damper∂ʳvel(joint::Translational{T}) where T
    A = nullspacemat(joint)
    AᵀA = zerodimstaticadjoint(A) * A
    Z = szeros(T, 3, 3)
    return [[-AᵀA * joint.damper * AᵀA; Z] [Z; Z]]
end
@inline function offdiagonal∂damper∂ʳvel(joint::Translational{T}, xa::AbstractVector, qa::UnitQuaternion, xb::AbstractVector, qb::UnitQuaternion) where T
    A = nullspacemat(joint)
    AᵀA = zerodimstaticadjoint(A) * A
    Z = szeros(T, 3, 3)
    return [[AᵀA * joint.damper * AᵀA; Z] [Z; Z]]
end
@inline function offdiagonal∂damper∂ʳvel(joint::Translational{T}, xb::AbstractVector, qb::UnitQuaternion) where T
    A = nullspacemat(joint)
    AᵀA = zerodimstaticadjoint(A) * A
    Z = szeros(T, 3, 3)
    return [[AᵀA * joint.damper * AᵀA; Z] [Z; Z]]
end

### Forcing
## Application of joint forces (for dynamics)
@inline function applyFτ!(joint::Translational{T}, statea::State, stateb::State, clear::Bool) where T
    F = joint.Fτ
    vertices = joint.vertices
    _, qa = posargsk(statea)
    _, qb = posargsk(stateb)

    Fa = vrotate(-F, qa)
    Fb = -Fa

    τa = vrotate(torqueFromForce(Fa, vrotate(vertices[1], qa)),inv(qa)) # in local coordinates
    τb = vrotate(torqueFromForce(Fb, vrotate(vertices[2], qb)),inv(qb)) # in local coordinates

    statea.Fk[end] += Fa
    statea.τk[end] += τa
    stateb.Fk[end] += Fb
    stateb.τk[end] += τb
    clear && (joint.Fτ = szeros(T,3))
    return
end
@inline function applyFτ!(joint::Translational{T}, stateb::State, clear::Bool) where T
    F = joint.Fτ
    vertices = joint.vertices
    _, qb = posargsk(stateb)

    Fb = F
    τb = vrotate(torqueFromForce(Fb, vrotate(vertices[2], qb)),inv(qb)) # in local coordinates

    stateb.Fk[end] += Fb
    stateb.τk[end] += τb
    clear && (joint.Fτ = szeros(T,3))
    return
end

## Forcing derivatives (for linearization)
# Control derivatives
@inline function ∂Fτ∂ua(joint::Translational, statea::State, stateb::State)
    vertices = joint.vertices
    _, qa = posargsk(statea)

    BFa = -VLmat(qa) * RᵀVᵀmat(qa)
    Bτa = -skew(vertices[1])

    return [BFa; Bτa]
end
@inline function ∂Fτ∂ub(joint::Translational, statea::State, stateb::State)
    vertices = joint.vertices
    xa, qa = posargsk(statea)
    xb, qb = posargsk(stateb)
    qbinvqa = qb\qa

    BFb = VLmat(qa) * RᵀVᵀmat(qa)
    Bτb = skew(vertices[2]) * VLmat(qbinvqa) * RᵀVᵀmat(qbinvqa)

    return [BFb; Bτb]
end
@inline function ∂Fτ∂ub(joint::Translational, stateb::State)
    vertices = joint.vertices
    _, qb = posargsk(stateb)

    BFb = I
    Bτb = skew(vertices[2]) * VLᵀmat(qb) * RVᵀmat(qb)

    return [BFb; Bτb]
end

# Position derivatives
@inline function ∂Fτ∂posa(joint::Translational{T}, statea::State, stateb::State) where T
    _, qa = posargsk(statea)
    _, qb = posargsk(stateb)
    F = joint.Fτ
    vertices = joint.vertices

    FaXa = szeros(T,3,3)
    FaQa = -2*VRᵀmat(qa)*Rmat(UnitQuaternion(F))*LVᵀmat(qa)
    τaXa = szeros(T,3,3)
    τaQa = szeros(T,3,3)
    FbXa = szeros(T,3,3)
    FbQa = 2*VRᵀmat(qa)*Rmat(UnitQuaternion(F))*LVᵀmat(qa)
    τbXa = szeros(T,3,3)
    τbQa = 2*skew(vertices[2])*VLᵀmat(qb)*Rmat(qb)*Rᵀmat(qa)*Rmat(UnitQuaternion(F))*LVᵀmat(qa)

    return FaXa, FaQa, τaXa, τaQa, FbXa, FbQa, τbXa, τbQa
end
@inline function ∂Fτ∂posb(joint::Translational{T}, statea::State, stateb::State) where T
    _, qa = posargsk(statea)
    _, qb = posargsk(stateb)
    F = joint.Fτ
    vertices = joint.vertices

    FaXb = szeros(T,3,3)
    FaQb = szeros(T,3,3)
    τaXb = szeros(T,3,3)
    τaQb = szeros(T,3,3)
    FbXb = szeros(T,3,3)
    FbQb = szeros(T,3,3)
    τbXb = szeros(T,3,3)
    τbQb = 2*skew(vertices[2])*VLᵀmat(qb)*Lmat(qa)*Lmat(UnitQuaternion(F))*Lᵀmat(qa)*LVᵀmat(qb)

    return FaXb, FaQb, τaXb, τaQb, FbXb, FbQb, τbXb, τbQb
end
@inline function ∂Fτ∂posb(joint::Translational{T}, stateb::State) where T
    xb, qb = posargsk(stateb)
    F = joint.Fτ
    vertices = joint.vertices

    FaXb = szeros(T,3,3)
    FaQb = szeros(T,3,3)
    τaXb = szeros(T,3,3)
    τaQb = szeros(T,3,3)
    FbXb = szeros(T,3,3)
    FbQb = szeros(T,3,3)
    τbXb = szeros(T,3,3)
    τbQb = 2*skew(vertices[2])*VLᵀmat(qb)*Lmat(UnitQuaternion(F))*LVᵀmat(qb)

    return FaXb, FaQb, τaXb, τaQb, FbXb, FbQb, τbXb, τbQb
end

### Minimal coordinates
## Position and velocity offsets
@inline function getPositionDelta(joint::Translational, body1::AbstractBody, body2::Body, x::SVector)
    Δx = zerodimstaticadjoint(nullspacemat(joint)) * x # in body1 frame
    return Δx
end
@inline function getVelocityDelta(joint::Translational, body1::AbstractBody, body2::Body, v::SVector)
    Δv = zerodimstaticadjoint(nullspacemat(joint)) * v # in body1 frame
    return Δv
end

## Minimal coordinate calculation
@inline function minimalCoordinates(joint::Translational, body1::Body, body2::Body)
    statea = body1.state
    stateb = body2.state
    return nullspacemat(joint) * g(joint, statea.xc, statea.qc, stateb.xc, stateb.qc)
end
@inline function minimalCoordinates(joint::Translational, body1::Origin, body2::Body)
    stateb = body2.state
    return nullspacemat(joint) * g(joint, stateb.xc, stateb.qc)
end
@inline function minimalVelocities(joint::Translational, body1::Body, body2::Body)
    statea = body1.state
    stateb = body2.state
    return nullspacemat(joint) * (stateb.vc - statea.vc)
end
@inline function minimalVelocities(joint::Translational, body1::Origin, body2::Body)
    stateb = body2.state
    return nullspacemat(joint) * stateb.vc
end


text:
"""
```
weight_kernel!(coeff_terms, log_e_1_terms, log_e_2_terms, φ_old, Ψ, y_t,
    s_t_nontemp, det_HH, inv_HH; initialize = false, parallel = false,
    poolmodel = false)
```

The outputs of the weight_kernel function are meant to speed up the adaptive φ
finding, so that we don't do the same matrix multiplication step in every
iteration of the root-solving algorithm.

The exponential terms are logged first and then exponentiated in the
incremental weight calculation so the problem is well-conditioned (i.e. not
exponentiating very large negative numbers).

This function modifies `coeff_terms`, `log_e_1_terms`, and `log_e_2_terms`.
"""
function weight_kernel!(coeff_terms::V, log_e_1_terms::V, log_e_2_terms::V,
                        φ_old::Float64, Ψ::Function, y_t::Vector{Float64},
                        s_t_nontemp::AbstractMatrix{Float64},
                        det_HH::Float64, inv_HH::Matrix{Float64};
                        initialize::Bool = false,
                        parallel::Bool = false,
                        poolmodel::Bool = false) where V<:AbstractVector{Float64}
    # Sizes
    n_particles = length(coeff_terms)
    n_obs = length(y_t)

    if poolmodel
        if parallel
            @sync @distributed for i in 1:n_particles
                if initialize
                    # Initialization step (using 2π instead of φ_old)
                    coeff_terms[i] = (2*pi)^(-n_obs/2) # this may need to be adjusted
                    log_e_1_terms[i] = 0.
                    log_e_2_terms[i] = log(Ψ(s_t_nontemp[:,i])) # pred dens passed -> log it
                else
                    coeff_terms[i] = (φ_old)^(-n_obs/2)
                    log_e_1_terms[i] = -φ_old * log(Ψ(s_t_nontemp[:,i]))
                    log_e_2_terms[i] = log(Ψ(s_t_nontemp[:,i]))
                end
            end
        else
            for i in 1:n_particles
                if initialize
                    # Initialization step (using 2π instead of φ_old)
                    coeff_terms[i] = (2*pi)^(-n_obs/2) # this may need to be adjusted
                    log_e_1_terms[i] = 0.
                    log_e_2_terms[i] = log(Ψ(s_t_nontemp[:,i]))
                else
                    coeff_terms[i] = (φ_old)^(-n_obs/2)
                    log_e_1_terms[i] = -φ_old * log(Ψ(s_t_nontemp[:,i]))
                    log_e_2_terms[i] = log(Ψ(s_t_nontemp[:,i]))
                end
            end
        end
    else
        if parallel
            @sync @distributed for i in 1:n_particles
                error    = y_t - Ψ(s_t_nontemp[:, i])
                sq_error = dot(error, inv_HH * error)

                if initialize
                    # Initialization step (using 2π instead of φ_old)
                    coeff_terms[i]   = (2*pi)^(-n_obs/2) * det_HH^(-1/2)
                    log_e_1_terms[i] = 0.
                    log_e_2_terms[i] = -1/2 * sq_error
                else
                    # Non-initialization step (tempering and final iteration)
                    coeff_terms[i]   = (φ_old)^(-n_obs/2)
                    log_e_1_terms[i] = -1/2 * (-φ_old) * sq_error
                    log_e_2_terms[i] = -1/2 * sq_error
                end
            end
        else
           for i in 1:n_particles
                error    = y_t - Ψ(s_t_nontemp[:, i])
                sq_error = dot(error, inv_HH * error)

                if initialize
                    # Initialization step (using 2π instead of φ_old)
                    coeff_terms[i]   = (2*pi)^(-n_obs/2) * det_HH^(-1/2)
                    log_e_1_terms[i] = 0.
                    log_e_2_terms[i] = -1/2 * sq_error
                else
                    # Non-initialization step (tempering and final iteration)
                    coeff_terms[i]   = (φ_old)^(-n_obs/2)
                    log_e_1_terms[i] = -1/2 * (-φ_old) * sq_error
                    log_e_2_terms[i] = -1/2 * sq_error
                end
            end
        end
    end
    return nothing
end

"""
```
next_φ(φ_old, coeff_terms, log_e_1_terms, log_e_2_terms, n_obs, r_star, stage;
    fixed_sched = [], findroot = bisection, xtol = 1e-3)
```

Return the next tempering factor `φ_new`. If `isempty(fixed_sched)`, adaptively
choose `φ_new` by setting InEff(φ_new) = r_star. Otherwise, return
`fixed_sched[stage]`.
"""
function next_φ(φ_old::Float64, coeff_terms::V, log_e_1_terms::V, log_e_2_terms::V,
                n_obs::Int, r_star::Float64, stage::Int;
                fixed_sched::Vector{Float64} = Float64[],
                findroot::Function = bisection,
                xtol::Float64 = 1e-3) where V<:AbstractVector{Float64}

    if isempty(fixed_sched)
        n_particles  = length(coeff_terms)
        inc_weights  = Vector{Float64}(undef, n_particles)
        norm_weights = Vector{Float64}(undef, n_particles)
        ineff0(φ) =
            ineff!(inc_weights, norm_weights, φ, coeff_terms, log_e_1_terms, log_e_2_terms, n_obs) - r_star

        if stage == 1 || (sign(ineff0(φ_old)) != sign(ineff0(1.0)))
            # Solve for optimal φ if either
            # 1. First stage
            # 2. Sign change from ineff0(φ_old) to ineff0(1.0)
            return findroot(ineff0, φ_old, 1.0, xtol = xtol)
        else
            # Otherwise, set φ = 1
            return 1.0
        end
    else
        return fixed_sched[stage]
    end
end

"""
```
correction!(inc_weights, norm_weights, φ_new, coeff_terms, log_e_1_terms,
     log_e_2_terms, n_obs)
```

Compute (and modify in-place) incremental weights w̃ₜʲ and normalized weights W̃ₜʲ:

        w̃ₜʲ(φₙ) = pₙ(yₜ|sₜʲ'ⁿ⁻¹) / pₙ₋₁(yₜ|sₜʲ'ⁿ⁻¹)
                = (φₙ/φₙ₋₁)^(d/2) exp{-1/2 (φₙ-φₙ₋₁) [yₜ-Ψ(sₜʲ'ⁿ⁻¹)]' Σᵤ⁻¹ [yₜ-Ψ(sₜʲ'ⁿ⁻¹)]}

        W̃ₜʲ(φₙ) = w̃ₜʲ(φₙ) / (1/M) ∑ w̃ₜʲ(φₙ)
"""
function correction!(inc_weights::Vector{Float64}, norm_weights::Vector{Float64},
                     φ_new::Float64, coeff_terms::V, log_e_1_terms::V, log_e_2_terms::V,
                     n_obs::Int) where V<:AbstractVector{Float64}
    # Compute incremental weights
    n_particles = length(inc_weights)
    for i = 1:n_particles
        inc_weights[i] =
            φ_new^(n_obs/2) * coeff_terms[i] * exp(log_e_1_terms[i]) * exp(φ_new*log_e_2_terms[i])
    end

    # Normalize weights
    norm_weights .= inc_weights ./ mean(inc_weights)

    return nothing
end

"""
```
ineff!(inc_weights, norm_weights, φ_new, coeff_terms, exp_1_terms, exp_2_terms, n_obs)
```

Compute and return InEff(φₙ), where:

        InEff(φₙ) = (1/M) ∑ (W̃ₜʲ(φₙ))²
"""
function ineff!(inc_weights::Vector{Float64}, norm_weights::Vector{Float64},
                φ_new::Float64, coeff_terms::V, exp_1_terms::V, exp_2_terms::V,
                n_obs::Int) where V<:AbstractVector{Float64}

    correction!(inc_weights, norm_weights, φ_new, coeff_terms, exp_1_terms,
                exp_2_terms, n_obs)
    return sum(norm_weights.^2) / length(norm_weights)
end


text:
"Returns type parameter for a `Val` type or instance."
val_param(::Val{T}) where {T} = T
val_param(::Type{Val{T}}) where {T} = T

"""
    args_tupletype_expr(arg_exprs, esc=identity)

Given `arg_exprs`, a list of positional argument expressions, such as
`[:(x::Int), :(y::Float64), :(z::T...)]`, return a `Tuple` type expression
whose parameters are the argument types, e.g. `Tuple{Int,Float64,Vararg{T}}`.
"""
function args_tupletype_expr(arg_exprs, esc=identity)
    ret = Expr(:curly)
    ret.args = map(arg_exprs) do arg
        # Detect splatting
        splatted = Meta.isexpr(arg, :(...), 1)
        if splatted
            arg = arg.args[1]
        end
        # Extract type expressions
        if Meta.isexpr(arg, :(::), 2)
            _, ty = arg.args
            ty = esc(ty)
        elseif arg isa Symbol
            ty = :Any
        else
            error("Unexpected form of argument: $arg")
        end
        # Wrap within Vararg if splatted
        if splatted
            ty = Expr(:curly, :Vararg, ty)
        end
        return ty
    end
    pushfirst!(ret.args, :Tuple)
    return ret
end
args_tupletype_expr(signature_def::Dict{Symbol}, esc=identity) =
    args_tupletype_expr(signature_def[:args], esc)

"Check whether an expression is a Vararg (including escaped Varargs)."
is_vararg_expr(e) = Meta.isexpr(e, [:where, :curly, :escape]) ?
    is_vararg_expr(e.args[1]) : (e == :Vararg)

"Fill-in unnamed arguments."
fill_unnamed(e) =
    Meta.isexpr(e, [:(::)], 1) ? Expr(:(::), gensym(), e.args[2]) : e

"Remove type annotation from expression, if present."
rm_type_annotation(e) =
    Meta.isexpr(e, [:(::)], 2) ? e.args[1] : e

"""
    generate_switch_stmt(cond_exprs, branch_exprs, default_expr=:nothing)

Generates a switch expression using `if` and `elseif`, where `cond_exprs`
are the expressions for each branch condition, `branch_exprs` are the
expressions for statements to be executed within each branch, and `default_expr`
is the default statement if all conditions fail.
"""
function generate_switch_stmt(cond_exprs, branch_exprs, default_expr=:nothing)
    @assert length(cond_exprs) == length(branch_exprs)
    expr = default_expr
    for i in length(cond_exprs):-1:1
        head = i == 1 ? :if : :elseif
        expr = Expr(head, cond_exprs[i], branch_exprs[i], expr)
    end
    return expr
end


text:
#=

Base.convert() methods to convert shapefile geometry to Luxor geometry.

Use `include()` to include this file if you've installed Shapefile.jl.

Luxor y coordinates increase downwards, but maps and shapefile coordinates increase upwards.
In the Northern hemisphere, 0 longitude is North Pole, so we negate y coordinates.

=#

# Convert a Shapefile rectangle to an array of points suitable for a Luxor rectangle
function Base.convert(::Type{Array{Luxor.Point, 1}}, R::Shapefile.Rect)
    return [R.left, R.top, R.right-R.left, R.bottom-R.top]
end

# Convert a Shapefile Point to a Luxor Point (flipping y-coordinates)
function  Base.convert(::Type{Luxor.Point}, pt::Shapefile.Point)
    return Luxor.Point(pt.x, -pt.y)
end

# Convert Shapefile PolyLine to a Luxor poly
function Base.convert(::Type{Array{Luxor.Point, 1}}, pl::Shapefile.Polyline)
    resultbbox = convert(Array{Luxor.Point, 1}, pl.MBR)
    resultpoly = Luxor.Point[]
    for pt in pl.points
        push!(resultpoly, convert(Luxor.Point, pt))
    end
    return resultbbox, resultpoly
end

# Convert Shapefile Polygon to a Luxor poly.
# some Shapefile polygons have multiple parts, but these aren't subpaths, just
# extra polygons. For example, France (in Europe) and French Guiana in (South America)
# belong to the same Shapefile.Polygon and the return value will include these and other
# polygons in the array.

function Base.convert(::Type{Array{Luxor.Point, 1}}, pl::Shapefile.Polygon)
    bbox = [Luxor.Point(pl.MBR.left, -pl.MBR.top), pl.MBR.right - pl.MBR.left, pl.MBR.top - pl.MBR.bottom]
    polygons = Array{Point, 1}[] # to hold all the polygons for a shape
    if length(pl.parts) >= 2
        points = Luxor.Point[]
        for i in 1:length(pl.parts) - 1
            subpolygonpoints = pl.points[pl.parts[i]+1: pl.parts[i + 1]]
            for pt in subpolygonpoints
                push!(points, convert(Luxor.Point, pt))
            end
            push!(polygons, points)
            points = Luxor.Point[]
        end
        # don't forget the rest
        subpolygonpoints = pl.points[pl.parts[end] + 1: end]
        for pt in subpolygonpoints
            push!(points, convert(Luxor.Point, pt))
        end
        push!(polygons, points)
    else
        points = Luxor.Point[]
        for pt in pl.points
            push!(points, convert(Luxor.Point, pt))
        end
        push!(polygons, points)
    end
    return polygons, bbox
end


text:
module ReversedStrings
using AutoHashEquals

export ReversedString, reversed
"""
    ReversedString{V}<:AbstractString

Lazy reversed `V` struct holding representation and caching `lastindex`.
"""
@auto_hash_equals struct ReversedString{V}<:AbstractString
    representation::V
    lastindex::Int
    ReversedString(x) =
        # if lastindex(x) == 1
        #     x
        # else
        new{typeof(x)}(x, lastindex(x))
        ##end
end

"""
    reversed(x::AbstractString)

A lazy implementation of `Base.reverse` for `String`s.
"""
reversed(x::ReversedString) = x.representation
reversed(x::AbstractString) = ReversedString(x)

"""
    reverse_index(x::ReversedString,i)

Return corresponding index in unreversed String `x.lastindex-i+1`.
"""
@inline function reverse_index(x::ReversedString,i)
    ri = x.lastindex-i+1
end

"""
    reverse_index(x::AbstractString,i)

Return original index `i`.
"""
@inline reverse_index(x::AbstractString,i) =
    i

"""
    Base.SubString(x::ReversedString,start::Int,stop::Int)

Create a reversed `SubString` using [`reverse_index`](@ref).
"""
@inline Base.SubString(x::ReversedString,start::Int,stop::Int) =
    reversed(SubString(x.representation, reverse_index(x, stop), reverse_index(x, start)))

"""
    Base.ncodeunits(x::ReversedString)

`ncodeunits(x.representation)`
"""
@inline Base.ncodeunits(x::ReversedString) = ncodeunits(x.representation)

"""
    Base.ncodeunits(x::ReversedString)

`codeunit(x.representation,`[`reverse_index`](@ref)`(x,i))`.
"""
@inline Base.codeunit(s::ReversedString, i::Integer) = codeunit(x.representation, reverse_index(x,i))

"""
    Base.isvalid(x::ReversedString,i::Int)

`isvalid(x.representation,`[`reverse_index`](@ref)`(x,i))`.
"""
@inline Base.isvalid(x::ReversedString,i::Int) =
    isvalid(x.representation,reverse_index(x,i))

"""
    Base.firstindex(x::ReversedString)

1
"""
@inline Base.firstindex(x::ReversedString) = 1

"""
    Base.lastindex(x::ReversedString)

cached `x.lastindex`.
"""
@inline Base.lastindex(x::ReversedString) = x.lastindex

"""
    Base.length(x::ReversedString)

`length(x.representation)`
"""
@inline Base.length(s::ReversedString) =
    length(s.representation)

"""
    Base.getindex(x::ReversedString,is::UnitRange{<:Integer})

`ReversedString` of `getindex` on representation.
!!! note
    `getindex` is creating a new `String` representation.
    See [`SubString`](@ref)
"""
@inline Base.getindex(x::ReversedString,is::UnitRange{<:Integer}) =
    reversed(getindex(x.representation,reverse_index(x,is.stop):reverse_index(x,is.start)))

"""
    Base.getindex(x::ReversedString,i::Integer)

`getindex(x.representation,`[`reverse_index`](@ref)`(x,i))`.
"""
@inline Base.getindex(x::ReversedString,i::Integer) =
    getindex(x.representation,reverse_index(x,i))

"""
    Base.thisind(x::ReversedString,i::Integer)

[`reverse_index`](@ref)`(thisind(x.representation,`[`reverse_index`](@ref)`(x,i)))`.
"""
@inline function Base.thisind(x::ReversedString,i::Int)
    ri = reverse_index(x, i)
    reverse_index(x, thisind(x.representation, ri))
end

"""
    Base.nextind(x::ReversedString,i[,n=1])

[`reverse_index`](@ref)`(nextind(x.representation,`[`reverse_index`](@ref)`(x,i))[,n])`.
"""
@inline function Base.nextind(x::ReversedString,i::Int)
    ri = reverse_index(x, i)
    reverse_index(x, prevind(x.representation, ri))
end

@inline function Base.nextind(x::ReversedString,i::Int,n::Int)
    ri = reverse_index(x, i)
    reverse_index(x, prevind(x.representation, ri, n))
end

"""
    Base.prevind(x::ReversedString,i[,n=1])

[`reverse_index`](@ref)`(prevind(x.representation,`[`reverse_index`](@ref)`(x,i))[,n])`.
"""
@inline function Base.prevind(x::ReversedString,i::Int)
    ri = reverse_index(x, i)
    reverse_index(x, nextind(x.representation, ri))
end

@inline function Base.prevind(x::ReversedString,i::Int,n::Int)
    ri = reverse_index(x, i)
    reverse_index(x,nextind(x.representation, ri, n))
end

"""
    Base.iterate(x::ReversedString[, state])

Iterate lazy reversed string.
"""
@inline Base.iterate(x::ReversedString) =
    iterate(x,1)

@inline Base.iterate(x::ReversedString,i::Int) =
    if reverse_index(x,i) >= 1 && reverse_index(x,i) <= lastindex(x) # ? <=
        x[i], nextind(x,i)
    else
        nothing
    end

end # module


text:
module IntroCrusher

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Intro Crusher" => Ahorn.EntityPlacement(
        Maple.IntroCrusher,
        "rectangle",
        Dict{String, Any}(),
        function(entity)
            entity.data["nodes"] = [(Int(entity.data["x"]) + Int(entity.data["width"]) + 8, Int(entity.data["y"]))]
            Ahorn.tileEntityFinalizer(entity)
        end,
    ),
)

Ahorn.editingOptions(entity::Maple.IntroCrusher) = Dict{String, Any}(
    "tiletype" => Ahorn.tiletypeEditingOptions()
)

Ahorn.nodeLimits(entity::Maple.IntroCrusher) = 1, 1

# The Intro Crusher functionality only works if the entity is >= 22 wide
Ahorn.minimumSize(entity::Maple.IntroCrusher) = 24, 8
Ahorn.resizable(entity::Maple.IntroCrusher) = true, true

function Ahorn.selection(entity::Maple.IntroCrusher)
    x, y = Ahorn.position(entity)
    nx, ny = Int.(entity.data["nodes"][1])

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    return [Ahorn.Rectangle(x, y, width, height), Ahorn.Rectangle(nx, ny, width, height)]
end

Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.IntroCrusher, room::Maple.Room) = Ahorn.drawTileEntity(ctx, room, entity, blendIn=false)

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.IntroCrusher, room::Maple.Room)
    x, y = Ahorn.position(entity)
    nodes = get(entity.data, "nodes", ())

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    if !isempty(nodes)
        nx, ny = Int.(nodes[1])
        cox, coy = floor(Int, width / 2), floor(Int, height / 2)

        material = get(entity.data, "tiletype", "3")[1]

        fakeTiles = Ahorn.createFakeTiles(room, nx, ny, width, height, material, blendIn=false)
        Ahorn.drawFakeTiles(ctx, room, fakeTiles, room.objTiles, true, nx, ny, clipEdges=true)
        Ahorn.drawArrow(ctx, x + cox, y + coy, nx + cox, ny + coy, Ahorn.colors.selection_selected_fc, headLength=6)
    end
end

end


text:
import JSON

using MetabolicFluxAnalysis

# load file
json_1_f = open("datasets/y76_no_labeling/y76Model.json")
json_1_s = readstring(json_1_f)

json_data = JSON.parse(json_1_s)

imp_s_result = MetabolicFluxAnalysis.import_json_species( json_data["metabolites"] )
imp_f_result = MetabolicFluxAnalysis.import_json_fluxes( json_data["reactions"] )

import_a = MetabolicFluxAnalysis.import_json( json_data )

net_a = import_a[1]
constraints_a = import_a[2]

tic()
fp  = MetabolicFluxAnalysis.get_flux_polyhedron(net_a , constraints_a ; max_val=Inf)
toc()
sx_fluxes = MetabolicFluxAnalysis.get_sx_fluxes(net_a)

fva_result = MetabolicFluxAnalysis.fva(net_a,constraints_a,collect(1:100))
fva_lb     = fva_result[1]
fva_ub     = fva_result[2]
#rfp = MetabolicFluxAnalysis.get_rounded_fulldimensional_parametrization( fp ; remove_red_constr=1)


text:
"""
`Bins{TS, TW, TBI, TT}`: A bin structure designed for WE

### Fields

* `Ω` - structure containing information for uniquely identifying each bin
* `n` - number of particles in each bin
* `target` - target number of particles in each bin
* `ν` - weight of each bin
* `d` - auxiliary bin data
"""
struct Bins{TS, TW<:AbstractFloat, TBI<:Integer, TT<:Real, TD} <: AbstractBins
   # structure which identifies the bins
   Ω::Vector{TS}
   # number of walkers in each bin
   n::Vector{TBI}
   # target number of walkers in each bin - specify if real/int with TT
   target::Vector{TT}
   # weight associated with each bin
   ν::Vector{TW}
   # auxiliary data structure for each bin
   d::Vector{TD}
end

function Base.push!(B::TB, Ω, n, target, ν, d) where {TB<:Bins}
   push!(B.Ω, Ω);
   push!(B.n, n);
   push!(B.target, target);
   push!(B.ν, ν);
   push!(B.d, d)
end

function Base.pop!(B::TB) where {TB<:Bins}
   Ω = pop!(B.Ω);
   n = pop!(B.n);
   target = pop!(B.target);
   ν = pop!(B.ν);
   d = pop!(B.d)
   #particles = pop!(B.particles)
   return Ω, n, target, ν, d
end

function Base.popfirst!(B::TB) where {TB<:Bins}
   Ω = popfirst!(B.Ω);
   n = popfirst!(B.n);
   target = popfirst!(B.target);
   ν = popfirst!(B.ν);
   d = popfirst!(B.d);
   #particles = pop!(B.particles)
   return Ω, n, target, ν, d
end

function Base.length(B::TB) where {TB<:Bins}
   return length(B.Ω)
end

function Base.eltype(B::TB) where {TB<:Bins}
   return typeof((B.Ω, B.n, B.target, B.ν, B.d))
end

function Base.iterate(B::TB, state = 1) where {TB<:Bins}

   if state > length(B)
      return nothing
   end
   return (B.Ω[state], B.n[state], B.target[state], B.ν[state], B.d[state]), state+1
end


text:
include("orm.jl")
using JLD
using Images
using FileIO
import FileIO: load,
               save
using StaticArrays
mutable struct Algorithm <: DBTable
    id::Union{Integer, Nothing}
    name::String
    params::Dict
end
function Algorithm(;id::Union{Integer, Nothing}=nothing,
                   name::String="",
                   params::Dict=Dict())
    return Algorithm(id, name, params)
end

mutable struct SegmentedImage <: DBTable
    id::Union{Integer, Nothing}
    dataset::String
    image::Integer
    file_path::String
    algorithm::Union{Algorithm, Integer}
    metrics::Dict
end
function SegmentedImage(;
                        id::Union{Integer, Nothing}=nothing,
                        dataset::String="",
                        image::Integer=0,
                        file_path::String="",
                        algorithm::Union{Algorithm, Integer}=-1,
                        metrics::Dict=Dict())
    return SegmentedImage(
        id,
        dataset,
        image,
        file_path,
        algorithm,
        metrics)
end
function save(segmented_image::SegmentedImage, img, segments, conn)
    insert(segmented_image, conn)
    img = map(i->segment_mean(segments, i), labels_map(segments))
    @info("Saving jld file")
    save(string(segmented_image.file_path, ".jld"),
         "segments", segments)

    @info("Saving png file")
    save(string(segmented_image.file_path, ".png"),
         colorview(RGB, img))
    @info("Checking that the file can be read")
    s = load(string(segmented_image.file_path, ".jld"))
end
function save(segmented_image::SegmentedImage, conn)
    update(segmented_image, conn)
end
function load(s::SegmentedImage)
    file_path = string(s.file_path, ".jld")
    @info("Opening file $file_path")
    loaded_file = load(file_path)
    @info("Loaded")
    return loaded_file["segments"]
end

mutable struct Result <: DBTable
    id::Union{Integer, Nothing}
    s::Union{SegmentedImage, Integer}
    results::Dict
end
function Result(;id::Union{Integer, Nothing}=nothing,
                s::Union{SegmentedImage, Integer}=-1,
                results::Dict=Dict())
    return Result(id, s, results)
end

function create_tables(con)
    execute(
        conn,
        create_table_sql(Algorithm))
    execute(
        conn,
        create_table_sql(SegmentedImage))
    execute(
        conn,
        create_table_sql(Result))
end

function sql_to_jsonb(d::Dict)
    return "'$(json(d))'::jsonb"
end

function filename(seg::SegmentedImage)
    base_filename = filename(seg.dataset, seg.image)
    return joinpath(seg.algorithm.name, base_filename)
end

function filename(dataset::RGBDDataset, image_number::Integer)
    return filename(name(dataset),image_number)
end
function filename(ds_name::String, image_number::Integer)
    return "$(ds_name)_$(image_number)_seg"
end

function exists(seg::SegmentedImage, conn; exclude::Array{Symbol}=[])
    do_exists = invoke(exists, Tuple{DBTable, Any}, seg, conn;exclude=exclude)
    if ! do_exists
        return do_exists
    else
        return isfile(string(seg.file_path, ".jld"))
    end
end
connection_params = Dict(:host=>"localhost",
                         :dbname=>"rgbd_segmentation_results",
                         :user=>"rgbd_segmentation",
                         :password=>"rgbd_seg")
connection = join([join(a,'=') for a in connection_params], ' ')
conn = LibPQ.Connection(connection)
create_tables(conn)

function tests()
    a = SegmentedImage(5,
                       "asa",
                       5,
                       "asdasd/aads",
                       Algorithm(15, "aaa", Dict("a"=>"b")))

    println(exists(a, conn))
    insert(a, conn)
    println(exists(a, conn))
end


text:
const Float = Float64

@enum EField Ex = 1 Ey Ez
@enum HField Hx = 1 Hy Hz


text:
@doc doc"""
Minimise le problème : ``min_{||s||< \delta_{k}} q_k(s) = s^{t}g + (1/2)s^{t}Hs``
                        pour la ``k^{ème}`` itération de l'algorithme des régions de confiance

# Syntaxe
```julia
sk = Gradient_Conjugue_Tronque(fk,gradfk,hessfk,option)
```

# Entrées :
   * **gradfk**           : (Array{Float,1}) le gradient de la fonction f appliqué au point xk
   * **hessfk**           : (Array{Float,2}) la Hessienne de la fonction f appliqué au point xk
   * **options**          : (Array{Float,1})
      - **delta**    : le rayon de la région de confiance
      - **max_iter** : le nombre maximal d'iterations
      - **tol**      : la tolérance pour la condition d'arrêt sur le gradient

# Sorties:
   * **s** : (Array{Float,1}) le pas s qui approche la solution du problème : ``min_{||s||< \delta_{k}} q(s)``

# Exemple d'appel:
```julia
gradf(x)=[-400*x[1]*(x[2]-x[1]^2)-2*(1-x[1]) ; 200*(x[2]-x[1]^2)]
hessf(x)=[-400*(x[2]-3*x[1]^2)+2  -400*x[1];-400*x[1]  200]
xk = [1; 0]
options = []
s = Gradient_Conjugue_Tronque(gradf(xk),hessf(xk),options)
```
"""
function Gradient_Conjugue_Tronque(gradfk,hessfk,options)

    "# Si option est vide on initialise les 3 paramètres par défaut"
    if options == []
        deltak = 2
        max_iter = 100
        tol = 1e-6
    else
        deltak = options[1]
        max_iter = options[2]
        tol = options[3]
    end

   n = length(gradfk)
   s = zeros(n)
   q(s) = gradfk' * s + (1 / 2) * transpose(s) * hessfk * s

   sj = zeros(n)
   gj = gradfk
   pj = -gradfk

   for i in 1 : max_iter
       kj = pj' * hessfk * pj

       if kj < 0
           a = norm(pj) ^ 2
           b = 2 * sj' * pj
           c = norm(sj) ^ 2 - deltak ^ 2
           d = sqrt(b ^ 2 - 4 * a * c)
           rho1, rho2 = (-b - d) / (2 * a), (-b + d) / (2 * a)

           if q(sj + rho1 * pj) > q(sj + rho2 * pj)
               rho = rho2
           else
               rho = rho1
           end
           s = sj + rho * pj
           break
       end

       alphaj = (gj' * gj) / kj

       if norm(sj + alphaj * pj) > deltak
           a = norm(pj) ^ 2
           b = 2 * sj' * pj
           c = norm(sj) ^ 2 - deltak ^ 2

           d = sqrt(b ^ 2 - 4 * a * c)
           alpha1, alpha2 = (-b - d) / (2 * a), (-b + d) / (2 * a)

           if alpha1 > 0
               alpha = alpha1
           elseif alpha2 > 0
               alpha = alpha2
               s = sj + alpha * pj
               break
           end
       end

       sj = sj + alphaj * pj
       gj_1 = gj + alphaj * hessfk * pj
       betaj = (gj_1' * gj_1) / (gj' * gj)
       pj = -gj_1 + betaj * pj

       if norm(gj_1) <= tol * norm(gradfk)
           s = sj
           break
       end
       gj = gj_1
   end

   return s
end

max_stars_repo_name: MOUDDENEHamza/ENSEEIHT

max_issues_repo_name: MOUDDENEHamza/ENSEEIHT

max_forks_repo_name: MOUDDENEHamza/ENSEEIHT


text:
using Interpolations
using Interpolations: tweight, tcoef, prefilter, copy_with_padding, prefilter!, BSplineInterpolation

A_x1 = 1 : 0.1 : 10
A_x2 = 1 : 0.5 : 20

f(x1, x2) = log(x1+x2)
A = [f(x1,x2) for x1 in A_x1, x2 in A_x2]

splinetype(r::AbstractRange) = BSpline(Quadratic(InPlace()))
splinetype(r::UnitRange) = BSpline(Constant())

b2 = BSpline(Constant())
b1 = BSpline(Quadratic(Free(OnGrid())))

it = (b1, b2,)
it isa NTuple{N,BSpline} where {N}

interpolate(A, it) # start here

# infer tw and tc, then
tw = tweight(A)
tc = tcoef(A)

# interpolate
@which interpolate(tw, tc, A, it)
# D:\software_libraries\julia\packages\Interpolations\0zZ6m\src\b-splines\b-splines.jl:132
# function interpolate(::Type{TWeights}, ::Type{TC}, A, it::IT) where {TWeights,TC,IT<:DimSpec{BSpline}}
#     Apad = prefilter(TWeights, TC, A, it)
#     BSplineInterpolation(TWeights, Apad, it, axes(A))
# end

# D:\software_libraries\julia\packages\Interpolations\0zZ6m\src\b-splines\prefiltering.jl:39
@which prefilter(tw, tc, A, it)
# function prefilter(::Type{TWeights}, ::Type{TC}, A::AbstractArray, it) where {TWeights,TC}
#     ret = copy_with_padding(TC, A, it)
#     prefilter!(TWeights, ret, it)
# end

# D:\software_libraries\julia\packages\Interpolations\0zZ6m\src\b-splines\prefiltering.jl:21
@which copy_with_padding(tc, A, it)
ret = copy_with_padding(tc, A, it)
size(ret)
# copy_with_padding(A, it) = copy_with_padding(eltype(A), A, it)
# function copy_with_padding(::Type{TC}, A, it::DimSpec{InterpolationType}) where {TC}
#     indsA = axes(A)
#     indspad = padded_axes(indsA, it)
#     coefs = padded_similar(TC, indspad)
#     if indspad == indsA
#         coefs = copyto!(coefs, A)
#     else
#         fill!(coefs, zero(TC))
#         ct!(coefs, indsA, A, indsA)
#     end
#     coefs
# end

@which BSplineInterpolation(tw, Apad, it, axes(A))
bsp = BSplineInterpolation(tw, Apad, it, axes(A))

Interpolations.padded_axes(axes(A), it)

function padded_coefs(::Type{TC}, A, it) where {TC}
    indsA = axes(A)
    indspad = Interpolations.padded_axes(indsA, it)
    return Interpolations.padded_similar(TC, indspad)
end
padded_coefs(A,it) = padded_coefs(tcoef(A), A, it)

sitp = scale(itp, A_x1, A_x2)

sitp(5., 10.) # exactly log(5 + 10)
sitp(5.6, 7.1) # approximately log(5.6 + 7.1)

struct EmaxCoefs{TC,IT<:DimSpec{BSpline},A1<:AbstractArray{T},A2<:AbstractArray{T}}
    EV::A1
    dEV::A2
    it::IT
    function EmaxCoefs(ev::A1, dev::A2, it::IT) where {TC,IT,A1<:AbstractArray{TC},A2<:AbstractArray{TC}}
        szdev = size(dev)
        (szdev[1:end-2]..., szdev[end]) == size(ev) || throw(DimensionMismatch())
        length(it)+1 == length(szdev) || throw(DimensionMismatch())
        new{T,IT,A1,A2}(ev, dev, it)
    end
end

EV(x::EmaxCoefs) = x.EV
dEV(x::EmaxCoefs) = x.dEV
ITEV(x::EmaxCoefs) = (x.it..., BSpline(Constant()))
ITdEV(x::EmaxCoefs) = (x.it..., BSpline(Constant()), BSpline(Constant()))

function EmaxCoefs(evs::DCDPEmax, it::Tuple)
    bspc = BSpline(Constant())
    ev_it = (it..., bspc, )
    dev_it = (it..., bspc, bspc,)

    ev_coef = padded_coefs(EV(evs), ev_it)
    dev_coef = padded_coefs(EV(evs), dev_it)

    return EmaxCoefs(ev_coef, dev_coef, it)
end

function copy_with_padding!(::Type{TC}, coefs, A, it::DimSpec{InterpolationType}) where {TC}
    indsA = axes(A)
    indspad = Interpolations.padded_axes(indsA, it)
    lengths.(indspad) == size(coefs)
    if indspad == indsA
        copyto!(coefs, A)
    else
        fill!(coefs, zero(TC))
        Interpolations.ct!(coefs, indsA, A, indsA)
    end
    return nothing
end

function update!(coefs::EmaxCoefs{TC}, evs::DCDPEmax, dograd) where {TC}
    TC == tcoef(EV(evs)) || throw(error())

    twev = tweights(EV(evs))
    copy_with_padding!(TC, EV(coefs), EV(evs), ITEV(coefs))
    prefilter!(twev, EV(coefs), it)

    if dograd
        twdev = tweights(dEV(evs))
        copy_with_padding!(TC, dEV(coefs), dEV(evs), ITdEV(coefs))
        prefilter!(twdev, dEV(coefs), dEV(evs), ITdEV(coefs))
    end
    return nothing
end

function interpolate(coefs::EmaxCoefs, evs::DCDPEmax)
    evtw = tweights(EV(evs))
    evax = axes(EV(evs))

    devtw = tweights(dEV(evs))
    devax = axes(dEV(evs))

    ev_itp = BSplineInterpolation(tw, EV(coefs), ITEV(coefs), ax)
    dev_itp = BSplineInterpolation(tw, dEV(coefs), ITdEV(coefs), ax)

    return ev_itp, dev_itp
end

function scale(coefs::EmaxCoefs, evs::DCDPEmax, ddm::DyanmicDrillintModel)
    ev, dev = interpolate(coefs, evs)
    ns = OneTo(length(statespace(ddm)))
    nk = OneTo(_nparm(reward(ddm)))
    sev = scale(ev, zspace(ddm)..., ns)
    sdev = scale(dev, zspace(ddm)..., nk, ns)
    return sev, sdev
end

# scaling 114
# indexint 32