text:
using MPFR_wrap
using Test
using Base.MPFR

@testset "assignment" begin
    prec, atol = 512, 1e-100
    c = BigFloat(precision=prec)

    @test set!(c, BigFloat("3", precision=prec)) == BigFloat("3", precision=prec)
    @test set!(c, UInt32(3)) == BigFloat("3", precision=prec)
    @test set!(c, Int32(-3)) == BigFloat("-3", precision=prec)
    @test set!(c, -3.0) == BigFloat("-3", precision=prec)
end

@testset "basic arithmetic functions" begin
    prec, atol = 512, 1e-100
    c = BigFloat("1.0", precision=prec)

    add!(c, BigFloat("1.2", precision=prec), BigFloat("1.3", precision=prec))
    @test isapprox(c, BigFloat("2.5", precision=prec), atol=atol)
    add!(c, BigFloat("1.2", precision=prec), UInt32(2))
    @test isapprox(c, BigFloat("3.2", precision=prec), atol=atol)
    add!(c, BigFloat("1.2", precision=prec), Int32(-2))
    @test isapprox(c, BigFloat("-0.8", precision=prec), atol=atol)
    add!(c, BigFloat("1.2", precision=prec), -2.0)
    @test isapprox(c, BigFloat("-0.8", precision=prec), atol=atol)
    add!(c, BigFloat("1.2", precision=prec), BigInt(-4))
    @test isapprox(c, BigFloat("-2.8", precision=prec), atol=atol)

    sub!(c, BigFloat("1.2", precision=prec), BigFloat("1.3", precision=prec))
    @test isapprox(c, BigFloat("-0.1", precision=prec), atol=atol)
    sub!(c, BigFloat("1.2", precision=prec), UInt32(2))
    @test isapprox(c, BigFloat("-0.8", precision=prec), atol=atol)
    sub!(c, BigFloat("1.2", precision=prec), Int32(-2))
    @test isapprox(c, BigFloat("3.2", precision=prec), atol=atol)
    sub!(c, BigFloat("1.2", precision=prec), -2.0)
    @test isapprox(c, BigFloat("3.2", precision=prec), atol=atol)
    sub!(c, BigFloat("1.2", precision=prec), BigInt(-4))
    @test isapprox(c, BigFloat("5.2", precision=prec), atol=atol)
    sub!(c, UInt32(2), BigFloat("1.2", precision=prec))
    @test isapprox(c, BigFloat("0.8", precision=prec), atol=atol)
    sub!(c, Int32(-2), BigFloat("1.2", precision=prec))
    @test isapprox(c, BigFloat("-3.2", precision=prec), atol=atol)
    sub!(c, -2.0, BigFloat("1.2", precision=prec))
    @test isapprox(c, BigFloat("-3.2", precision=prec), atol=atol)
    sub!(c, BigInt(-4), BigFloat("1.2", precision=prec))
    @test isapprox(c, BigFloat("-5.2", precision=prec), atol=atol)

    sqr!(c, BigFloat("2", precision=prec))
    @test isapprox(c, BigFloat("4", precision=prec), atol=atol)

    sqrt!(c, BigFloat("4", precision=prec))
    @test isapprox(c, BigFloat("2", precision=prec), atol=atol)
    @test_throws DomainError sqrt!(c, BigFloat("-4", precision=prec))

    rec_sqrt!(c, BigFloat("4", precision=prec))
    @test isapprox(c, BigFloat("0.5", precision=prec), atol=atol)

    cbrt!(c, BigFloat("-27", precision=prec))
    @test isapprox(c, BigFloat("-3", precision=prec), atol=atol)

    rootn!(c, BigFloat("81", precision=prec), UInt32(4))
    @test isapprox(c, BigFloat("3", precision=prec), atol=atol)
    @test_throws DomainError rootn!(c, BigFloat("-4"), UInt32(4))

    pow!(c, BigFloat("3", precision=prec), BigFloat("3", precision=prec))
    @test isapprox(c, BigFloat("27", precision=prec), atol=atol)
    pow!(c, BigFloat("3", precision=prec), UInt32(3))
    @test isapprox(c, BigFloat("27", precision=prec), atol=atol)
    pow!(c, BigFloat("3", precision=prec), Int32(3))
    @test isapprox(c, BigFloat("27", precision=prec), atol=atol)
    pow!(c, BigFloat("3", precision=prec), BigInt(3))
    @test isapprox(c, BigFloat("27", precision=prec), atol=atol)
    pow!(c, UInt32(3), UInt32(3))
    @test isapprox(c, BigFloat("27", precision=prec), atol=atol)
    pow!(c, UInt32(3), BigFloat("3", precision=prec))
    @test isapprox(c, BigFloat("27", precision=prec), atol=atol)

    neg!(c, BigFloat("-27", precision=prec))
    @test isapprox(c, BigFloat("27", precision=prec), atol=atol)

    abs!(c, BigFloat("-27", precision=prec))
    @test isapprox(c, BigFloat("27", precision=prec), atol=atol)

    dim!(c, BigFloat("27", precision=prec), BigFloat("20", precision=prec))
    @test isapprox(c, BigFloat("7", precision=prec), atol=atol)
    @test_throws DomainError dim!(c, BigFloat("27"), BigFloat(NaN))

    mul_2!(c, BigFloat("3", precision=prec), UInt32(2))
    @test isapprox(c, BigFloat("12", precision=prec), atol=atol)
    mul_2!(c, BigFloat("3", precision=prec), Int32(2))
    @test isapprox(c, BigFloat("12", precision=prec), atol=atol)

    div_2!(c, BigFloat("12", precision=prec), UInt32(2))
    @test isapprox(c, BigFloat("3", precision=prec), atol=atol)
    div_2!(c, BigFloat("12", precision=prec), Int32(2))
    @test isapprox(c, BigFloat("3", precision=prec), atol=atol)
end

@testset "comparison" begin
    @test cmp_2exp(BigFloat("8"), UInt32(2), 4) < 0
    @test cmp_2exp(BigFloat("8"), Int32(2), 2) == 0

    @test number_p(BigFloat(3))

    @test !regular_p(BigFloat(NaN))

    @test unordered_p(BigFloat(3.0), BigFloat(NaN))
end

@testset "special functions" begin
    prec, atol = 512, 1e-100
    c, d = BigFloat(precision=prec), BigFloat(precision=prec)

    log!(c, BigFloat(ℯ, precision=prec))
    @test isapprox(c, BigFloat("1.0", precision=prec), atol=atol)

    y0!(c, BigFloat("1.0", precision=64))
    @test isapprox(c, BigFloat("0.088256964215676957"), atol=1e-15)

    y1!(c, BigFloat("1.0", precision=64))
    @test isapprox(c, BigFloat("-0.781212821300288716"), atol=1e-15)

    @test isapprox(j0!(c, BigFloat("1.0", precision=prec)),
                   jn!(d, Int32(0), BigFloat("1.0", precision=prec)), atol=atol)

    @test isapprox(y0!(c, BigFloat("1.0", precision=prec)),
                   yn!(d, Int32(0), BigFloat("1.0", precision=prec)), atol=atol)

    fma!(c, BigFloat("3", precision=prec), BigFloat("2", precision=prec),
         BigFloat("1", precision=prec))
    @test isapprox(c, BigFloat("7", precision=prec), atol=atol)

    hypot!(c, BigFloat("3", precision=prec), BigFloat("4", precision=prec))
    @test isapprox(c, BigFloat("5", precision=prec), atol=atol)
end

@testset "integer and reminder functions" begin
    prec, atol = 512, 1e-100
    c, d = BigFloat(precision=prec), BigFloat(precision=prec)

    ceil!(c, BigFloat("3.4", precision=prec))
    @test isapprox(c, BigFloat("4.0", precision=prec), atol=atol)

    floor!(c, BigFloat("3.4", precision=prec))
    @test isapprox(c, BigFloat("3.0", precision=prec), atol=atol)

    round!(c, BigFloat("3.4", precision=prec))
    @test isapprox(c, BigFloat(3.0, precision=prec), atol=atol)

    roundeven!(c, BigFloat(3.4, precision=prec))
    @test isapprox(c, BigFloat(3.0, precision=prec), atol=atol)

    rint!(c, BigFloat("3.4", precision=prec), Base.MPFR.MPFRRoundUp)
    @test isapprox(c, BigFloat("4.0", precision=prec), atol=atol)

    frac!(c, BigFloat("3.4", precision=prec), Base.MPFR.MPFRRoundUp)
    @test isapprox(c, BigFloat("0.4", precision=prec), atol=atol)

    modf!(c, d, BigFloat("3.4", precision=prec), Base.MPFR.MPFRRoundUp)
    @test isapprox(c, BigFloat("3.0", precision=prec), atol=atol)
    @test isapprox(d, BigFloat("0.4", precision=prec), atol=atol)

    @test integer_p(BigFloat("4.0"))
end


text:
function calc_grad( pw::PWGrid, Potentials, Focc, psi::Array{ComplexF64,2} )

    Npoints = size(psi)[1]
    Nstates = size(psi)[2]
    Ω = pw.Ω
    Ns = pw.Ns
    #
    grad = zeros( ComplexF64, Npoints, Nstates )

    H_psi = op_H( pw, Potentials, psi )
    for i = 1:Nstates
        grad[:,i] = H_psi[:,i]
        for j = 1:Nstates
            grad[:,i] = grad[:,i] - dot( psi[:,j], H_psi[:,i] ) * psi[:,j]
        end
        grad[:,i] = Focc[i]*grad[:,i]
    end
    return grad

end

ext: jl


text:
### parsing.jl
#
# Copyright (C) 2018 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: uncertainty, error propagation, physics
#
# This file is a part of Measurements.jl.
#
# License is MIT "Expat".
#
### Commentary:
#
# This file defines the methods to represent `Measurement` objects in various places.
#
### Code:

function truncated_print(io::IO, m::Measurement, error_digits::Int;
                         atbeg = "", atend = "", pm = "±")
    print(io, atbeg,
          (iszero(m.err) || isnan(m.err)) ? m.val : round(m.val, digits =
                                                          -Base.hidigit(m.err, 10) + error_digits),
          get(io, :compact, false) ? pm : " $pm ",
          round(m.err, sigdigits=error_digits), atend)
end

full_print(io::IO, measure::Measurement) =
    print(io, measure.val, get(io, :compact, false) ? "±" : " ± ", measure.err)

function Base.show(io::IO, m::Measurement)
    error_digits = get(io, :error_digits, 2)
    if error_digits > 0
            truncated_print(io, m, error_digits)
    elseif error_digits == 0
        full_print(io, m)
    else
        error("`error_digits` must be non-negative")
    end # if
end

function Base.show(io::IO, ::MIME"text/latex", measure::Measurement)
    error_digits = get(io, :error_digits, 2)
    truncated_print(io, measure, error_digits, atbeg = "\$", atend = "\$", pm = "\\pm")
end

for mime in (MIME"text/x-tex", MIME"text/x-latex")
    function Base.show(io::IO, ::mime, measure::Measurement)
        error_digits = get(io, :error_digits, 2)
        truncated_print(io, measure, error_digits, pm = "\\pm")
    end
end

# Representation of complex measurements.  Print something that is easy to
# understand and that can be meaningfully copy-pasted into the REPL, at least
# for standard numeric types.
for mime in (MIME"text/plain", MIME"text/x-tex", MIME"text/x-latex")
    function Base.show(io::IO, mtype::mime, measure::Complex{<:Measurement})
        r, i = reim(measure)
        compact = get(io, :compact, false)
        print(io, "(")
        show(io, mtype, r)
        print(io, ")")
        if signbit(i) && !isnan(i)
            i = -i
            print(io, compact ? "-" : " - ")
        else
            print(io, compact ? "+" : " + ")
        end
        print(io, "(")
        show(io, mtype, i)
        print(io, ")im")
    end
end

# This is adapted from base/show.jl for Complex type.
function Base.alignment(io::IO, measure::Measurement)
    m = match(r"^(.*[\±])(.*)$", sprint(show, measure, context=io, sizehint=0))
    m === nothing ? (length(sprint(show, measure, context=io, sizehint=0)), 0) :
        (length(m.captures[1]), length(m.captures[2]))
end

### Juno pretty printing
@require Juno="e5e0dc1b-0480-54bc-9374-aad01c23163d" begin
    Juno.render(i::Juno.Inline, measure::Measurement) =
        Juno.render(i, Juno.Row(measure.val, Text(" ± "), measure.err))

    Juno.Row(measure.val, Text(" ± "), measure.err)

    function Juno.render(ji::Juno.Inline, cm::Complex{<:Measurement})
        r, i = reim(cm)
        if signbit(i) && !isnan(i)
            i = -i
            sss = " - "
        else
            sss = " + "
        end
        Juno.render(ji, Juno.Row("(", Juno.render(ji, r), ")", sss,
                                 "(", Juno.render(ji, i), ")im"))
    end
end


text:
function interpolate_rgb(c1::Colorant, c2::Colorant, w::Real)
  rgb1 = RGBA(c1)
  rgb2 = RGBA(c2)
  r = interpolate(rgb1.r, rgb2.r, w)
  g = interpolate(rgb1.g, rgb2.g, w)
  b = interpolate(rgb1.b, rgb2.b, w)
  a = interpolate(rgb1.alpha, rgb2.alpha, w)
  RGBA(r, g, b, a)
end

function interpolate(v1::Real, v2::Real, w::Real)
  (1-w) * v1 + w * v2
end

# --------------------------------------------------------------

# Methods to automatically generate gradients for color selection based on
# background color and a short list of seed colors

# here are some magic constants that could be changed if you really want
const _lightness_darkbg = [80.0]
const _lightness_lightbg = [60.0]
const _lch_c_const = [60]

function adjust_lch(color, l, c)
    lch = convert(LCHab, color)
    convert(RGBA{Float64}, LCHab(l, c, lch.h))
end

function lightness_from_background(bgcolor)
	bglight = convert(LCHab, bgcolor).l
	bglight < 50.0 ? _lightness_darkbg[1] : _lightness_lightbg[1]
end

function gradient_from_list(cs)
    zvalues = get_zvalues(length(cs))
    indices = sortperm(zvalues)
    sorted_colors = map(RGBA, cs[indices])
    sorted_zvalues = zvalues[indices]
    ColorGradient(sorted_colors, sorted_zvalues)
end

function generate_colorgradient(bgcolor = plot_color(:white);
                               color_bases = plot_color([colorant"steelblue",colorant"orangered"]),
                               lightness = lightness_from_background(bgcolor),
                               chroma = _lch_c_const[1],
                               n = 17)
    seed_colors = vcat(bgcolor, map(c -> adjust_lch(c, lightness, chroma), color_bases))
	seed_colors = convert(Vector{RGB{Float64}}, seed_colors)
	colors = distinguishable_colors(
		n,
		seed_colors,
		lchoices=Float64[lightness],
		cchoices=Float64[chroma],
		hchoices=range(0; stop=340, length=20)
	)[2:end]
	gradient_from_list(colors)
end

function get_color_palette(palette, bgcolor::Colorant, numcolors::Integer)
	grad = if palette == :auto
		generate_colorgradient(bgcolor)
	else
		cgrad(palette)
	end
	zrng = get_zvalues(numcolors)
	RGBA{Float64}[grad[z] for z in zrng]
end

get_color_palette(palette::Vector{C}, bgcolor::Colorant, numcolors::Integer) where {C<:Colorant} = palette

# ----------------------------------------------------------------------------------

function getpctrange(n::Int)
    n > 0 || error()
    n == 1 && return zeros(1)
    zs = [0.0, 1.0]
    for i in 3:n
        sorted = sort(zs)
        diffs = diff(sorted)
        widestj = 0
        widest = 0.0
        for (j,d) in enumerate(diffs)
            if d > widest
                widest = d
                widestj = j
            end
        end
        push!(zs, sorted[widestj] + 0.5 * diffs[widestj])
    end
    zs
end

function get_zvalues(n::Int)
    offsets = getpctrange(ceil(Int,n/4)+1)/4
    offsets = vcat(offsets[1], offsets[3:end])
    zvalues = Float64[]
    for offset in offsets
        append!(zvalues, offset .+ [0.0, 0.5, 0.25, 0.75])
    end
    vcat(zvalues[1], 1.0, zvalues[2:n-1])
end

# ----------------------------------------------------------------------------------

function darken(c, v=0.1)
    rgba = convert(RGBA, c)
    r = max(0, min(rgba.r - v, 1))
    g = max(0, min(rgba.g - v, 1))
    b = max(0, min(rgba.b - v, 1))
    RGBA(r,g,b,rgba.alpha)
end
function lighten(c, v=0.3)
    darken(c, -v)
end

# isbackgrounddark(bgcolor::Color) = Lab(bgcolor).l < 0.5

# move closer to lighter/darker depending on background value
function adjust_away(val, bgval, vmin=0., vmax=100.)
  if bgval < 0.5 * (vmax+vmin)
    tmp = max(val, bgval)
    return 0.5 * (tmp + max(tmp, vmax))
  else
    tmp = min(val, bgval)
    return 0.5 * (tmp + min(tmp, vmin))
  end
end

# borrowed from http://stackoverflow.com/a/1855903:
lightness_level(c::Colorant) = 0.299 * red(c) + 0.587 * green(c) + 0.114 * blue(c)

isdark(c::Colorant) = lightness_level(c) < 0.5
islight(c::Colorant) = !isdark(c)

function Base.convert(::Type{RGB}, h::Unsigned)
  mask = 0x0000FF
  RGB([(x & mask) / 0xFF for x in  (h >> 16, h >> 8, h)]...)
end

make255(x) = round(Int, 255 * x)

function rgb_string(c::Colorant)
  @sprintf("rgb(%d, %d, %d)", make255(red(c)), make255(green(c)), make255(blue(c)))
end

function rgba_string(c::Colorant)
	@sprintf("rgba(%d, %d, %d, %1.3f)", make255(red(c)), make255(green(c)), make255(blue(c)), alpha(c))
end


text:
using Test, Polyhedra

@testset "./src/LocalPolytope/facets.jl" begin

using BellScenario

@testset "facets()" begin
    @testset "3-simplex normalized" begin
        vertices = [[1,0],[0,1],[0,0]]
        facet_dict = LocalPolytope.facets(vertices)

        @test facet_dict["equalities"] == Vector{Vector{Int64}}[]
        @test facet_dict["facets"] == [[-1,0,0],[0,-1,0],[1,1,1]]
        @test length(keys(facet_dict)) == 2
    end

    @testset "3-simplex generalized" begin
        vertices = [[1,0,0],[0,1,0],[0,0,1]]
        facet_dict = LocalPolytope.facets(vertices)

        @test facet_dict["equalities"] == [[1,1,1,1]]
        @test facet_dict["facets"] == [[0,-1,0,0],[0,0,-1,0],[0,1,1,1]]
        @test length(keys(facet_dict)) == 2
    end

    @testset "CH inequalities" begin
        scenario = BipartiteNonSignaling(2,2,2,2)
        vertices = LocalPolytope.vertices(scenario)
        facet_dict = LocalPolytope.facets(vertices)

        @test facet_dict["equalities"] == Vector{Vector{Int64}}[]
        @test facet_dict["facets"] == [
            [0, 0, 0, 0, -1, 0, 0, 0, 0], [0, 0, 0, 0, 0, -1, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, -1, 0, 0], [0, 0, 0, 0, 0, 0, 0, -1, 0],
            [-1, 0, 0, 0, 0, 1, 0, 0, 0], [-1, 0, 0, 0, 1, 0, 0, 0, 0],
            [0, -1, 0, 0, 0, 0, 0, 1, 0], [0, -1, 0, 0, 0, 0, 1, 0, 0],
            [0, 0, -1, 0, 0, 0, 1, 0, 0], [0, 0, -1, 0, 1, 0, 0, 0, 0],
            [0, 0, 0, -1, 0, 0, 0, 1, 0], [0, 0, 0, -1, 0, 1, 0, 0, 0],
            [-1, 0, -1, 0, 1, 1, 1, -1, 0], [-1, 0, 0, -1, 1, 1, -1, 1, 0],
            [0, -1, -1, 0, 1, -1, 1, 1, 0], [0, -1, 0, -1, -1, 1, 1, 1, 0],
            [0, 1, 0, 1, 0, 0, 0, -1, 1], [0, 1, 1, 0, 0, 0, -1, 0, 1],
            [1, 0, 0, 1, 0, -1, 0, 0, 1], [1, 0, 1, 0, -1, 0, 0, 0, 1],
            [0, 1, 0, 1, 1, -1, -1, -1, 1], [0, 1, 1, 0, -1, 1, -1, -1, 1],
            [1, 0, 0, 1, -1, -1, 1, -1, 1], [1, 0, 1, 0, -1, -1, -1, 1, 1]
        ]
        @test length(keys(facet_dict)) == 2
    end
end

@testset "facets(::Polyhedron)" begin
    scenario = LocalSignaling(3,3,2)
    local_poly = LocalPolytope.vrep(scenario)

    @test !hrepiscomputed(local_poly)

    facets = LocalPolytope.facets(local_poly)
    @test length(facets) == 15

    @test hrepiscomputed(local_poly)

    match_vertices = LocalPolytope.vertices(scenario)
    match_facets = LocalPolytope.facets(match_vertices)["facets"]

    @test facets == match_facets
end

end


text:
module Beep
using WAV

export beep

include("notifyme.jl")

end


text:
module TestDSP
using Test
include("simulator.jl")

p = Parameters(0,2,5,15)

const delta_frac_bits = Fix_frac_bits - 15
const maxval = int2fix(((1<<15)-1)*(1 << delta_frac_bits))
const minval = Fix(-1)

triangle(x) = abs.(mod.(x .+ 1.5,2).-1)*2 .- 1

function test_waveform(f, inst, imm; tol = 0)
	t = int2fix.((-32768:32767)*(1 << delta_frac_bits))
	dp = DataPath(p)

	y = eval_inst!.((dp,),(inst,), t, imm)
	y_expected = f.(t)
	if tol == 0
		@test y == y_expected
	else
		@test maximum(abs.(y-y_expected)) <= tol
	end
end
test_waveform(f, inst) = test_waveform(f, inst, Fix(0))

function test_waveforms()
	triangle_i = Instruction(inst_triangle, IContext(0,-1))
	test_waveform(t->(min.(maxval, triangle(t))), triangle_i)

	test_waveform(t->(t >= 0 ? maxval : minval), inst_pulse, Fix(0))
	t0 = mem_trunc(p, Fix(0.73981))
	test_waveform(t->(t >= t0 ? maxval : minval), inst_pulse, t0)

	for k=0:12
		inst = Instruction(inst_pulse_exp2, IContext(0,-k))
		t0 = Fix(-1 + 1/2^k)
		test_waveform(t->(t >= t0 ? maxval : minval), inst)
	end

	test_waveform(t->sin(pi*t), inst_sin, Fix(0); tol = 6.3e-3)
end

function test_lfsr()
	initial = mem_trunc(p, int2fix(1<<delta_frac_bits))
	program = [
		(inst_imm,Fix(0.5)),
		(Instruction(inst_acc_to_imm, IContext(0,2)), Fix(0)), # force overflow every loop with imm += 2^2 * 0.5
		(Instruction(inst_lfsr, IContext(0)), initial),
		(inst_output_a,0)
	]
	dsp = create_dsp(p, program)
	y = run_loops!(dsp, 65535)
	@test dsp.data[3] == initial

	program = [
		(inst_imm,Fix(0.5)),
		(Instruction(inst_acc_to_imm, IContext(0)), Fix(0)), # only overflow every 4th loop
		(Instruction(inst_lfsr, IContext(0)), initial),
		(inst_output_a,0)
	]
	dsp = create_dsp(p, program)
	y2 = run_loops!(dsp, 65535*4)
	@test dsp.data[3] == initial

	x = 1
	y_expected = zeros(Fix, 65535)
	y2_expected = zeros(Fix, 65535*4)
	for k=1:65536
		x = ((x << 1) | lfsr16_bit(x)) & 65535
		xf = int2fix((((x .+ 32768) .& 65535) .- 32768) << delta_frac_bits)
		if (k < 65536)
			y_expected[k] = xf

			y2_expected[(k-1)*4+1] = xf
			y2_expected[(k-1)*4+2] = xf
		end
		if (k > 1)
			y2_expected[(k-1)*4-1] = xf
			y2_expected[(k-1)*4-0] = xf
		end
	end
	@test y == y_expected
	@test y2 == y2_expected
end

function test_bias()
	dp = DataPath(p)

	imm = Fix(0)
	biases = [0, 0x10000, 0x08000, 0x18000]
	for k=1:4
		imm = exec!(dp, inst_dither_update, imm)
		@test fix2int(imm) == 2*k
		@test dp.bias == biases[k]
		@test madd(p, Fix(2^-15), Fix((0x1fffc - dp.bias)*2^-17), Fix(0), dp.bias) == 0
		@test madd(p, Fix(2^-15), Fix((0x20000 - dp.bias)*2^-17), Fix(0), dp.bias) == Fix(2^-15)
	end
end

function test_dither_acc()
	# Test that dithering can give us freq_exp bits of extra frequency accuracy over a period of freq_exp cycles
	freq = reg_trunc(p, Fix(1/7))
	for freq_exp=1:12
		for k=1:2
			program = [
				(inst_imm, freq),
				(Instruction(inst_acc_to_imm, IContext(1, -freq_exp)), Fix(0)),
				(inst_output_a, Fix(0)),
			]
			if (k == 2)
				push!(program, (inst_dither_update, Fix(2^-15)))
			end
			dsp = create_dsp(p, program)
			y = run_loops!(dsp, 2^freq_exp)
			if (k == 1)
				@test y[end] != freq
			else
				@test y[end] == freq
			end
		end
	end
end

test_waveforms()
test_lfsr()

@test bit_reverse(1+8, 4) == 8+1
@test bit_reverse(1+8, 6) == 32+4
test_bias()
test_dither_acc()

end

max_stars_repo_name: toivoh/fmicrosynth

max_issues_repo_name: toivoh/fmicrosynth

max_forks_repo_name: toivoh/fmicrosynth


text:
findreplace(ex::Symbol, dict) = get(dict, ex, ex)
function findreplace(ex::Expr, dict)
    Expr(ex.head, map(x->findreplace(x, dict), ex.args)...)
end
findreplace(ex, dict) = ex

"""
```julia
ode_def_opts(name::Symbol,opts::Dict{Symbol,Bool},curmod,ex::Expr,params...;depvar=:t)
```

The core internal. Users should only interact with this through the `@ode_def_*` macros.

Options are self-explanatory by name mapping to `ODEFunction`:

- build_tgrad
- build_jac
- build_expjac
- build_invjac
- build_invW
- build_invW_t
- build_hes
- build_invhes
- build_dpfuncs

`depvar` sets the symbol for the dependent variable.

Example:

```julia
opts = Dict{Symbol,Bool}(
      :build_tgrad => true,
      :build_jac => true,
      :build_expjac => false,
      :build_invjac => true,
      :build_invW => true,
      :build_invW_t => true,
      :build_hes => false,
      :build_invhes => false,
      :build_dpfuncs => true)
```
"""
function ode_def_opts(name::Symbol,opts::Dict{Symbol,Bool},curmod,ex::Expr,params...;depvar=:t)
  # depvar is the dependent variable. Defaults to t
  # M is the mass matrix in RosW, must be a constant!

  origex = copy(ex) # Save the original expression

  if !(eltype(params) <: Symbol)
    error("The syntax for ParameterizedFunctions has changed. Simply list the parameters at the end, i.e. `a b c d`, instead of `a=5.0 b=>3.0 ...`. Parameters are defined in the problem type. See the documentation for more information.")
  end
  params = Symbol[params...]

  ## Build independent variable dictionary
  indvar_dict,syms = build_indvar_dict(ex,depvar)
  ####

  t = Symbolics.unwrap((@variables t)[1])
  vars = Symbolics.unwrap.([(@variables $x(t))[1] for x in syms])
  params = Symbolics.unwrap.([(@parameters $x)[1] for x in Symbol[params...]])

  vars_dict = Dict(x=>Symbol(v) for (x, v) in zip(syms, vars))

  # replace x with x(t) if it's a var
  ex = findreplace(ex, vars_dict)

  # Build the Expressions

  # Run find replace to make the function expression
  symex = copy(ex) # Different expression for symbolic computations
  #ode_findreplace(ex,symex,indvar_dict,params)
  funcs = build_component_funcs(symex)
  mtk_ops = modelingtoolkitize_expr.(funcs,([t;vars;params],),(curmod,))

  D = ModelingToolkit.Differential(t)

  mtk_diffeqs = [D(vars[i]) ~ mtk_ops[i] for i in 1:length(vars)]

  sys = ODESystem(mtk_diffeqs,t,vars,params,name=gensym(:Parameterized))

  f_ex_oop,f_ex_iip = ModelingToolkit.generate_function(sys,vars,params)

  if opts[:build_tgrad]
    try
      tgrad_ex_oop,tgrad_ex_iip = ModelingToolkit.generate_tgrad(sys,vars,params)
    catch
      @warn "tgrad construction failed"
      tgrad_ex_oop,tgrad_ex_iip = nothing,nothing
    end
  else
    tgrad_ex_oop,tgrad_ex_iip = nothing,nothing
  end

  if opts[:build_jac]
    try
      J_ex_oop,J_ex_iip = ModelingToolkit.generate_jacobian(sys,vars,params)
    catch
      @warn "Jacobian construction failed"
      J_ex_oop,J_ex_iip = nothing,nothing
    end
  else
    J_ex_oop,J_ex_iip = nothing,nothing
  end

  if opts[:build_invW] && length(mtk_diffeqs) < 4
    try
      W_exs = ModelingToolkit.generate_factorized_W(sys,vars,params,false)
      W_ex_oop,W_ex_iip = W_exs[1]
      W_t_ex_oop,W_t_ex_iip = W_exs[2]
    catch
      @warn "W-expression construction failed"
      W_ex_oop,W_ex_iip = (nothing,nothing)
      W_t_ex_oop,W_t_ex_iip = (nothing,nothing)
    end
  else
    W_ex_oop,W_ex_iip = (nothing,nothing)
    W_t_ex_oop,W_t_ex_iip = (nothing,nothing)
  end

  fname = gensym(:ParameterizedDiffEqFunction)
  tname = gensym(:ParameterizedTGradFunction)
  jname = gensym(:ParameterizedJacobianFunction)
  Wname = gensym(:ParameterizedWFactFunction)
  W_tname = gensym(:ParameterizedW_tFactFunction)
  funcname = gensym(:ParameterizedODEFunction)

  if tgrad_ex_oop !== nothing
    full_tex = quote
      $tname($(tgrad_ex_oop.args[1].args...)) = $(tgrad_ex_oop.args[2])
      $tname($(tgrad_ex_iip.args[1].args...)) = $(tgrad_ex_iip.args[2])
    end
  else
    full_tex = quote
      $tname = nothing
    end
  end

  if J_ex_oop !== nothing
    full_jex = quote
      $jname($(J_ex_oop.args[1].args...)) = $(J_ex_oop.args[2])
      $jname($(J_ex_iip.args[1].args...)) = $(J_ex_iip.args[2])
    end
  else
    full_jex = quote
      $jname = nothing
    end
  end

  if W_ex_oop !== nothing
    full_wex = quote
      $Wname($(W_ex_oop.args[1].args...)) = $(W_ex_oop.args[2])
      $Wname($(W_ex_iip.args[1].args...)) = $(W_ex_iip.args[2])
      $W_tname($(W_t_ex_oop.args[1].args...)) = $(W_t_ex_oop.args[2])
      $W_tname($(W_t_ex_iip.args[1].args...)) = $(W_t_ex_iip.args[2])
    end
  else
    full_wex = quote
      $Wname = nothing
      $W_tname = nothing
    end
  end

  quote
    struct $name{F,TG,TJ,TW,TWt,S} <: ParameterizedFunctions.DiffEqBase.AbstractParameterizedFunction{true}
      f::F
      mass_matrix::ParameterizedFunctions.LinearAlgebra.UniformScaling{Bool}
      analytic::Nothing
      tgrad::TG
      jac::TJ
      jvp::Nothing
      vjp::Nothing
      jac_prototype::Nothing
      sparsity::Nothing
      Wfact::TW
      Wfact_t::TWt
      paramjac::Nothing
      syms::Vector{Symbol}
      indepvar::Symbol
      colorvec::Nothing
      sys::S
    end

    (f::$name)(args...) = f.f(args...)

    $fname($(f_ex_oop.args[1].args...)) = $(f_ex_oop.args[2])
    $fname($(f_ex_iip.args[1].args...)) = $(f_ex_iip.args[2])
    $full_tex
    $full_jex
    $full_wex

    $name($fname,ParameterizedFunctions.LinearAlgebra.I,nothing,$tname,$jname,nothing,nothing,
          nothing,nothing,$Wname,$W_tname,nothing,$syms,$(Meta.quot(depvar)),nothing,$sys)
  end |> esc
end

max_stars_repo_name: JuliaDiffEq/ParameterizedFunctions

max_issues_repo_name: JuliaDiffEq/ParameterizedFunctions

max_forks_repo_name: JuliaDiffEq/ParameterizedFunctions


text:
using Lattices

include("chain.jl")
include("square.jl")


text:
function shifted_by(obj::FieldTimeShiftable, arg0::RealFieldElement)
    return jcall(obj, "shiftedBy", FieldTimeInterpolable, (RealFieldElement,), arg0)
end

function shifted_by(obj::FieldTimeShiftable, arg0::jdouble)
    return jcall(obj, "shiftedBy", FieldTimeInterpolable, (jdouble,), arg0)
end


text:
"""
    module PyExprIO

Functions to read and write Python expressions.
"""
module PyExprIO

export readpyexpr

using ..PyObjects

mutable struct State
    buf::String
    pos::Int
    posend::Int
end

State(buf::String) = State(buf, firstindex(buf), lastindex(buf))

"""
    readpyexpr(io_or_filename; simplify=false)

Reads a single expression from a .py file.
"""
function readpyexpr(io::IO; simplify=nothing)
    ans = parse_any(State(read(io, String)))
    if simplify !== nothing
        ans = PyObjects.simplify(ans)
        ans = PyObjects.walk(simplify, ans)
    end
    return ans
end
readpy(fn::AbstractString; kw...) = open(io->readpy(io; kw...), fn)

function inc(s::State)
    s.pos = nextind(s.buf, s.pos)
    return
end

function char(s::State)
    return s.buf[s.pos]
end

function atend(s::State)
    return s.pos > s.posend
end

function skip_space(s::State)
    while isspace(char(s))
        inc(s)
    end
end

function parse_any(s::State)
    i = skip_space(s)
    c = char(s)
    if ('a' ≤ c ≤ 'z') || ('A' ≤ c ≤ 'Z') || (c == '_')
        ident = parse_ident(s)
        if ident == "True"
            return true
        elseif ident == "False"
            return false
        elseif ident == "None"
            return nothing
        else
            error("not implemented: ident $(repr(ident))")
        end
    elseif c == ''' || c == '"'
        return parse_str(s)
    elseif c == '['
        return parse_list(s)
    elseif c == '('
        return parse_tuple(s)
    elseif c == '{'
        return parse_dict_or_set(s)
    elseif ('0' ≤ c ≤ '9') || (c == '-') || (c == '+') || (c == '.')
        return parse_number(s)
    else
        error("unexpected char: $(repr(c))")
    end
end

function parse_ident(s::State)
    i0 = i1 = s.pos
    inc(s)
    while !atend(s)
        c = char(s)
        if isletter(c) || isnumeric(c) || (c == '_')
            i1 = s.pos
            inc(s)
        else
            break
        end
    end
    return SubString(s.buf, i0:i1)
end

function parse_str(s::State)
    # TODO: triple-quotes
    c0 = char(s)
    inc(s)
    io = IOBuffer()
    while true
        c = char(s)
        inc(s)
        if c == c0
            break
        elseif c == '\\'
            error("escapes not implemented")
        else
            write(io, c)
        end
    end
    return String(take!(io))
end

function parse_list(s::State)
    xs = Any[]
    inc(s)
    while true
        skip_space(s)
        c = char(s)
        if c == ']'
            inc(s)
            break
        else
            x = parse_any(s)
            push!(xs, x)
            skip_space(s)
            c = char(s)
            if c == ','
                inc(s)
            elseif c == ']'
                inc(s)
                break
            else
                error("unexpected $(repr(c)), expecting ',' or ']'")
            end
        end
    end
    return PyList(xs)
end

function parse_tuple(s::State)
    xs = Any[]
    inc(s)
    while true
        skip_space(s)
        c = char(s)
        if c == ')'
            inc(s)
            break
        else
            x = parse_any(s)
            push!(xs, x)
            skip_space(s)
            c = char(s)
            if c == ','
                inc(s)
            elseif c == ')'
                inc(s)
                break
            else
                error("unexpected $(repr(c)), expecting ',' or ')'")
            end
        end
    end
    return PyTuple(xs)
end

function parse_dict_or_set(s::State)
    inc(s)
    skip_space(s)
    c = char(s)
    if c == '}'
        inc(s)
        return PyDict()
    else
        x0 = parse_any(s)
        skip_space(s)
        c = char(s)
        if c == '}'
            inc(s)
            return PySet(Any[x0])
        elseif c == ','
            inc(s)
            xs = Any[x0]
            while true
                skip_space(s)
                c = char(s)
                if c == '}'
                    inc(s)
                    break
                else
                    x = parse_any(s)
                    push!(xs, x)
                    skip_space(s)
                    c = char(s)
                    if c == ','
                        inc(s)
                    elseif c == '}'
                        inc(s)
                        break
                    else
                        error("unexpected $(repr(c)), expecting ',' or '}'")
                    end
                end
            end
            return PySet(xs)
        elseif c == ':'
            inc(s)
            k0 = x0
            skip_space(s)
            v0 = parse_any(s)
            xs = Pair{Any,Any}[Pair{Any,Any}(k0, v0)]
            skip_space(s)
            c = char(s)
            if c == ','
                inc(s)
                while true
                    skip_space(s)
                    c = char(s)
                    if c == '}'
                        inc(s)
                        break
                    else
                        k = parse_any(s)
                        skip_space(s)
                        c = char(s)
                        c == ':' || error("unexpected $(repr(c)), expecting ':'")
                        inc(s)
                        v = parse_any(s)
                        push!(xs, Pair{Any,Any}(k, v))
                        skip_space(s)
                        c = char(s)
                        if c == ','
                            inc(s)
                        elseif c == '}'
                            inc(s)
                            break
                        else
                            error("unexpected $(repr(c)), expecting ',' or '}'")
                        end
                    end
                end
            elseif c == '}'
                inc(s)
            else
                error("unexpected $(repr(c)), expecting ',' or '}'")
            end
            return PyDict(xs)
        else
            error("unexpected $(repr(c)), expecting ',', '}' or ':'")
        end
    end
end

function parse_number(s::State)
    i0 = i1 = s.pos
    inc(s)
    while !atend(s)
        c = char(s)
        if ('0' ≤ c ≤ '9') || (c == '-') || (c == '+') || (c == '.') || ('a' ≤ c ≤ 'f') || ('A' ≤ c ≤ 'F') || (c == 'x') || (c == 'X') || (c == 'b') || (c == 'B') || (c == 'o') || (c == 'O')
            i1 = s.pos
            inc(s)
        else
            break
        end
    end
    x = SubString(s.buf, i0:i1)
    if startswith(x, "0x") || startswith(x, "0o") || startswith(x, "0b")
        y = parse(BigInt, x)
    elseif ('.' in x) || ('e' in x) || ('f' in x)
        y = parse(Float64, x)
    else
        y = parse(BigInt, x)
    end
    if y isa BigInt
        y1 = mod(y, Int)
        if y1 == y
            y = y1
        end
    end
    return y
end

end


text:
"""
This benchmark is for demonstrating that `@simd ivdep` can be
beneficial for some cases.  It is constructed such that it is hard for
LLVM to insert run-time memory aliasing check and include a vectorized
version of the loop.
"""
module BenchForeachSeqSumMany

using BenchmarkTools
using ThreadsX

@inline cumsumargs(x, xs...) = (x, _cumsumargs(x, xs...)...)
@inline _cumsumargs(acc, x, xs...) = (acc + x, _cumsumargs(acc + x, xs...)...)
@inline _cumsumargs(acc) = ()

@noinline function manual_sum_many!(vecs::NTuple{N,AbstractArray}) where {N}
    @simd ivdep for i in eachindex(vecs...)
        sums = cumsumargs(map(a -> (@inbounds a[i]), vecs)...)
        ntuple(Val(N)) do j
            Base.@_inline_meta
            @inbounds vecs[j][i] = sums[j]
        end
    end
end

@noinline foreach_sum_many!(vecs::NTuple{N,AbstractArray}, simd) where {N} =
    let nvecs = Val(N)
        ThreadsX.Implementations.foreach_linear_seq(eachindex(vecs...), simd) do i
            Base.@_inline_meta
            sums = cumsumargs(map(a -> (@inbounds a[i]), vecs)...)
            ntuple(nvecs) do j
                Base.@_inline_meta
                @inbounds vecs[j][i] = sums[j]
            end
        end
    end

suite = BenchmarkGroup()

for nvecs in [8]
    m = 2^10
    vecs = Tuple([randn(m) for _ in 1:nvecs])

    s0 = suite[:nvecs=>nvecs] = BenchmarkGroup()
    s0["man"] = @benchmarkable(manual_sum_many!(vecs), setup = (vecs = map(copy, $vecs)))
    s1 = s0["tx"] = BenchmarkGroup()
    for simd in [false, true, :ivdep]
        s1[:simd=>simd] = @benchmarkable(
            foreach_sum_many!(vecs, $(Val(simd))),
            setup = (vecs = map(copy, $vecs)),
        )
    end
end

end  # module
BenchForeachSeqSumMany.suite


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "libiio"
version = v"0.23.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/analogdevicesinc/libiio.git", "92d6a35f3d8d721cda7d6fe664b435311dd368b4")
]

dependencies = [
    Dependency("libusb_jll", compat="~1.0.24"),
    Dependency("XML2_jll", compat="~2.9.12"),
    Dependency("libaio_jll", compat="~0.3.112"),
]

# Bash recipe for building across all platforms
script = raw"""
cd libiio

mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DHAVE_DNS_SD=OFF \
    -DENABLE_IPV6=OFF \
    ..
make -j${nproc}
if [[ "${target}" == *-apple-* ]]; then
    cp -r iio.framework ${libdir}
else
    make install
fi
"""

include("../../fancy_toys.jl")

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(p -> !Sys.isapple(p), supported_platforms(;experimental=true))
platforms_macos = filter!(p -> Sys.isapple(p), supported_platforms(;experimental=true))

# The products that we will ensure are always built
products = [
    LibraryProduct("libiio", :libiio)
]

products_apple = [
    FrameworkProduct("iio", :libiio)
]

# Build the tarballs, and possibly a `build.jl` as well.
if any(should_build_platform.(triplet.(platforms_macos)))
    build_tarballs(ARGS, name, version, sources, script, platforms_macos, products_apple, dependencies; julia_compat="1.6")
end
if any(should_build_platform.(triplet.(platforms)))
    build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")
end

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"Create the matrix S which will be inverted to find the scattering coefficients."
function scattering_matrix(medium::PhysicalMedium, particles::AbstractParticles, t_matrices::Vector, ω::T, basis_order::Integer)::Matrix{Complex{T}} where T
    # Generate response for one specific k
    # Number of particles
    P = length(particles)

    # Length of scattering basisfor each particle
    N = basisorder_to_basislength(typeof(medium),basis_order)

    # No particles means no scattering
    if P == 0
        return Matrix{Complex{T}}(undef,0,0)
    end

    # Faire: this could potentially return an MMatrix
    function S_block(j,l)
        if j == l
            return zeros(Complex{T}, N, N)
        else
            x_lj = origin(particles[j]) .- origin(particles[l])
            U = outgoing_translation_matrix(medium, basis_order, ω, x_lj)
            return - transpose(U) * t_matrices[l]
        end
    end

    # Evaluates all the blocks in the big matrix
    S_blocks = [S_block(j,l) for j in 1:P, l in 1:P]

    # Reshape S_blocks into big matrix
    S = zeros(Complex{T}, P*N, P*N)
    for i in 1:P
        for j in 1:P
            S[((i-1)*N+1):(i*N), ((j-1)*N+1):(j*N)] .= S_blocks[i,j]
        end
    end

    return S
end


text:
using StanSample, ParetoSmooth, NamedTupleTools
using StatisticalRethinking
import StatisticalRethinking: pk_plot
import ParetoSmooth: loo_compare

df = CSV.read(sr_datadir("WaffleDivorce.csv"), DataFrame);
scale!(df, [:Marriage, :MedianAgeMarriage, :Divorce])
data = (N=size(df, 1), D=df.Divorce_s, A=df.MedianAgeMarriage_s,
    M=df.Marriage_s)

stan5_1 = "
data {
    int < lower = 1 > N; // Sample size
    vector[N] D; // Outcome
    vector[N] A; // Predictor
}
parameters {
    real a; // Intercept
    real bA; // Slope (regression coefficients)
    real < lower = 0 > sigma;    // Error SD
}
transformed parameters {
    vector[N] mu;               // mu is a vector
    for (i in 1:N)
        mu[i] = a + bA * A[i];
}
model {
    a ~ normal(0, 0.2);         //Priors
    bA ~ normal(0, 0.5);
    sigma ~ exponential(1);
    D ~ normal(mu , sigma);     // Likelihood
}
generated quantities {
    vector[N] log_lik;
    for (i in 1:N)
        log_lik[i] = normal_lpdf(D[i] | mu[i], sigma);
}
";

stan5_2 = "
data {
    int N;
    vector[N] D;
    vector[N] M;
}
parameters {
    real a;
    real bM;
    real<lower=0> sigma;
}
transformed parameters {
    vector[N] mu;
    for (i in 1:N)
        mu[i]= a + bM * M[i];

}
model {
    a ~ normal( 0 , 0.2 );
    bM ~ normal( 0 , 0.5 );
    sigma ~ exponential( 1 );
    D ~ normal( mu , sigma );
}
generated quantities {
    vector[N] log_lik;
    for (i in 1:N)
        log_lik[i] = normal_lpdf(D[i] | mu[i], sigma);
}
";

stan5_3 = "
data {
  int N;
  vector[N] D;
  vector[N] M;
  vector[N] A;
}
parameters {
  real a;
  real bA;
  real bM;
  real<lower=0> sigma;
}
transformed parameters {
    vector[N] mu;
    for (i in 1:N)
        mu[i] = a + bA * A[i] + bM * M[i];
}
model {
  a ~ normal( 0 , 0.2 );
  bA ~ normal( 0 , 0.5 );
  bM ~ normal( 0 , 0.5 );
  sigma ~ exponential( 1 );
  D ~ normal( mu , sigma );
}
generated quantities{
    vector[N] log_lik;
    for (i in 1:N)
        log_lik[i] = normal_lpdf(D[i] | mu[i], sigma);
}
";

function loo_compare(models::Vector{SampleModel};
    loglikelihood_name="log_lik", model_names=nothing, sort_models=true)

    if isnothing(model_names)
        mnames = [models[i].name for i in 1:length(models)]
    end

    nmodels = length(models)

    ka = Vector{KeyedArray}(undef, nmodels)
    ll = Vector{Array{Float64, 3}}(undef, nmodels)

    for i in 1:length(models)
        ka[i] = read_samples(models[i], :keyedarray)
        ll[i] = permutedims(Array(matrix(ka[i], loglikelihood_name)), [1, 3, 2])
    end

    loo_compare(ll; model_names=mnames, sort_models)
end

function waic_compare(models::Vector{SampleModel};
    loglikelihood_name="log_lik", model_names=nothing, sort_models=true)

    if isnothing(model_names)
        mnames = [models[i].name for i in 1:length(models)]
    end

    nmodels = length(models)

    ka = Vector{KeyedArray}(undef, nmodels)
    ll = Vector{Array{Float64, 3}}(undef, nmodels)
    llp = Vector{Array{Float64, 3}}(undef, nmodels)
    llpr = Vector{Array{Float64, 3}}(undef, nmodels)

    for i in 1:length(models)
        ka[i] = read_samples(models[i], :keyedarray)
        ll[i] = permutedims(Array(matrix(ka[i], loglikelihood_name)), [1, 3, 2])
        llp[i] = permutedims(ll[1],[ 1, 3, 2])
        if i == 1
            ndraws, nchains, nobs = size(llp[1])
        end
        llpr[i] = reshape(llp, ndraws*nchains, nobs)
    end

    compare(llpr, Val(:waic); mnames=mnames)
end

m5_1s = SampleModel("m5.1s", stan5_1)
rc5_1s = stan_sample(m5_1s; data)

m5_2s = SampleModel("m5.2s", stan5_2)
rc5_2s = stan_sample(m5_2s; data)

m5_3s = SampleModel("m5.3s", stan5_3)
rc5_3s = stan_sample(m5_3s; data)

if success(rc5_1s) && success(rc5_2s) && success(rc5_3s)

    p5_1s = read_samples(m5_1s, :particles)
    NamedTupleTools.select(p5_1s, (:a, :bA, :sigma)) |> display
    p5_2s = read_samples(m5_2s, :particles)
    NamedTupleTools.select(p5_2s, (:a, :bM, :sigma)) |> display
    p5_3s = read_samples(m5_3s, :particles)
    NamedTupleTools.select(p5_3s, (:a, :bA, :bM, :sigma)) |> display

    models = [m5_1s, m5_2s, m5_3s]
    loglikelihood_name = :log_lik
    loo_comparison = loo_compare(models)
    println()
    for (i, psis) in enumerate(loo_comparison.psis)
        psis |> display
        pk_plot(psis.pointwise(:pareto_k))
        savefig(joinpath(@__DIR__, "m5.$(i)s.png"))
    end
    println()
    loo_comparison |> display
end
#=
With SR/ulam():
```
       PSIS    SE dPSIS  dSE pPSIS weight
m5.1u 126.0 12.83   0.0   NA   3.7   0.67
m5.3u 127.4 12.75   1.4 0.75   4.7   0.33
m5.2u 139.5  9.95  13.6 9.33   3.0   0.00
```
=#


text:
using RubyCall
using Base.Test

# write your own tests here
@test_approx_eq rbModule(:Math)[:sqrt](3) sqrt(3)
a=RubyCall.RbValue([13,7,3,9])
@test a[:sort]()[1] == 3


text:
export GuildCreate,
    GuildUpdate,
    GuildDelete,
    GuildBanAdd,
    GuildBanRemove,
    GuildEmojisUpdate,
    GuildIntegrationsUpdate,
    GuildMemberAdd,
    GuildMemberRemove,
    GuildMemberUpdate,
    GuildMembersChunk,
    GuildRoleCreate,
    GuildRoleUpdate,
    GuildRoleDelete

"""
Sent when a new [`Guild`](@ref) is created.
"""
struct GuildCreate <: AbstractEvent
    guild::Guild
end
@boilerplate GuildCreate :docs
GuildCreate(d::Dict{String, Any}) = GuildCreate(Guild(d))

"""
Sent when a [`Guild`](@ref) is updated.
"""
struct GuildUpdate <: AbstractEvent
    guild::Guild
end
@boilerplate GuildUpdate :docs
GuildUpdate(d::Dict{String, Any}) = GuildUpdate(Guild(d))

"""
Sent when a [`Guild`](@ref) is deleted.
"""
struct GuildDelete <: AbstractEvent
    guild::AbstractGuild  # Supposed to be an UnavailableGuild.
end
@boilerplate GuildDelete :docs
GuildDelete(d::Dict{String, Any}) = GuildDelete(AbstractGuild(d))

"""
Sent when a [`User`](@ref) is banned from a [`Guild`](@ref).
"""
struct GuildBanAdd <: AbstractEvent
    guild_id::Snowflake
    user::User
end
@boilerplate GuildBanAdd :dict :docs

"""
Sent when a [`User`](@ref) is unbanned from a [`Guild`](@ref).
"""
struct GuildBanRemove <: AbstractEvent
    guild_id::Snowflake
    user::User
end
@boilerplate GuildBanRemove :dict :docs

"""
Sent when a [`Guild`](@ref) has its [`Emoji`](@ref)s updated.
"""
struct GuildEmojisUpdate <: AbstractEvent
    guild_id::Snowflake
    emojis::Vector{Emoji}
end
@boilerplate GuildEmojisUpdate :dict :docs

"""
Sent when a [`Guild`](@ref) has its [`Integration`](@ref)s updated.
"""
struct GuildIntegrationsUpdate <: AbstractEvent
    guild_id::Snowflake
end
@boilerplate GuildIntegrationsUpdate :dict :docs

"""
Sent when a [`Member`](@ref) is added to a [`Guild`](@ref).
"""
struct GuildMemberAdd <: AbstractEvent
    guild_id::Snowflake
    member::Member
end
@boilerplate GuildMemberAdd :docs
GuildMemberAdd(d::Dict{String, Any}) = GuildMemberAdd(snowflake(d["guild_id"]), Member(d))

"""
Sent when a [`Member`](@ref) is removed from a [`Guild`](@ref).
"""
struct GuildMemberRemove <: AbstractEvent
    guild_id::Snowflake
    user::User
end
@boilerplate GuildMemberRemove :dict :docs

"""
Sent when a [`Member`](@ref) is updated in a [`Guild`](@ref).
"""
struct GuildMemberUpdate <: AbstractEvent
    guild_id::Snowflake
    roles::Vector{Snowflake}
    user::User
    nick::Union{String, Nothing}  # Not supposed to be nullable.
end
@boilerplate GuildMemberUpdate :dict :docs

"""
Sent when the [`Client`](@ref) requests guild members with [`request_guild_members`](@ref).
"""
struct GuildMembersChunk <: AbstractEvent
    guild_id::Snowflake
    members::Vector{Member}
end
@boilerplate GuildMembersChunk :dict :docs

"""
Sent when a new [`Role`](@ref) is created in a [`Guild`](@ref).
"""
struct GuildRoleCreate <: AbstractEvent
    guild_id::Snowflake
    role::Role
end
@boilerplate GuildRoleCreate :dict :docs

"""
Sent when a [`Role`](@ref) is updated in a [`Guild`](@ref).
"""
struct GuildRoleUpdate <: AbstractEvent
    guild_id::Snowflake
    role::Role
end
@boilerplate GuildRoleUpdate :dict :docs

"""
Sent when a [`Role`](@ref) is deleted from a [`Guild`](@ref).
"""
struct GuildRoleDelete <: AbstractEvent
    guild_id::Snowflake
    role_id::Snowflake
end
@boilerplate GuildRoleDelete :dict :docs


text:
_nextpow2exp(n) = (8 * sizeof(n) - leading_zeros(n) - (count_ones(n) == 1))

_get_bit_depth(img::BitArray) = 8  # TODO: write 1 bit-depth images
_get_bit_depth(img::AbstractArray{C}) where {C<:Colorant} = __get_bit_depth(eltype(C))
_get_bit_depth(img::AbstractArray{T}) where {T<:Normed} = __get_bit_depth(T)
# __get_bit_depth(::Type{Normed{T,1}}) where T = 1  # TODO: write 1 bit-depth images
# __get_bit_depth(::Type{Normed{T,2}}) where T = 2  # TODO: write 2 bit-depth images
# __get_bit_depth(::Type{Normed{T,4}}) where T = 4  # TODO: write 4 bit-depth images
__get_bit_depth(::Type{Bool}) = 8  # TODO: write 1 bit-depth images
__get_bit_depth(::Type{Normed{T,8}}) where T = 8
__get_bit_depth(::Type{Normed{T,16}}) where T = 16
__get_bit_depth(::Type{Normed{T,N}}) where {T,N} = ifelse(N <= 8, 8, 16)
__get_bit_depth(::Type{<:AbstractFloat}) = 8
_get_bit_depth(img::AbstractArray{T}) where {T<:AbstractFloat} = 8
_get_bit_depth(img::AbstractArray{<:Bool}) = 8  # TODO: write 1 bit-depth images
_get_bit_depth(img::AbstractArray{<:UInt8}) = 8
_get_bit_depth(img::AbstractArray{<:UInt16}) = 16

function _inspect_png_read(fpath, gamma::Union{Nothing,Float64}=nothing)
    fp = open_png(fpath)
    png_ptr = create_read_struct()
    info_ptr = create_info_struct(png_ptr)
    png_init_io(png_ptr, fp)
    png_set_sig_bytes(png_ptr, PNG_BYTES_TO_CHECK)
    png_read_info(png_ptr, info_ptr)

    width = png_get_image_width(png_ptr, info_ptr)
    height = png_get_image_height(png_ptr, info_ptr)
    color_type_orig = png_get_color_type(png_ptr, info_ptr)
    color_type = color_type_orig
    bit_depth_orig = png_get_bit_depth(png_ptr, info_ptr)
    bit_depth = bit_depth_orig
    num_channels = png_get_channels(png_ptr, info_ptr)
    interlace_type = png_get_interlace_type(png_ptr, info_ptr)

    backgroundp = png_color_16p()
    if png_get_bKGD(png_ptr, info_ptr, Ref(backgroundp)[]) != 0
        png_set_background(png_ptr, backgroundp, PNG_BACKGROUND_GAMMA_FILE, 1, 1.0)
    end

    if color_type == PNG_COLOR_TYPE_PALETTE
        png_set_palette_to_rgb(png_ptr)
        color_type = PNG_COLOR_TYPE_RGB
    end

    if color_type == PNG_COLOR_TYPE_GRAY && bit_depth < 8
        png_set_expand_gray_1_2_4_to_8(png_ptr)
        png_set_packing(png_ptr)
        bit_depth = UInt8(8)
    end

    if png_get_valid(png_ptr, info_ptr, PNG_INFO_tRNS) != 0
        png_set_tRNS_to_alpha(png_ptr)
        if color_type == PNG_COLOR_TYPE_GRAY || color_type == PNG_COLOR_TYPE_RGB
            color_type |= PNG_COLOR_MASK_ALPHA
        end
    end

    screen_gamma = PNG_DEFAULT_sRGB
    image_gamma = Ref{Cdouble}(0.0)
    intent = Ref{Cint}(-1)
    white_x = Ref{Cdouble}(0.0)
    white_y = Ref{Cdouble}(0.0)
    red_x = Ref{Cdouble}(0.0)
    red_y = Ref{Cdouble}(0.0)
    green_x = Ref{Cdouble}(0.0)
    green_y = Ref{Cdouble}(0.0)
    blue_x = Ref{Cdouble}(0.0)
    blue_y = Ref{Cdouble}(0.0)
    if isnothing(gamma)
        if png_get_valid(png_ptr, info_ptr, PNG_INFO_sRGB) != 0
            if png_get_sRGB(png_ptr, info_ptr, intent) != 0
                png_set_gamma(png_ptr, screen_gamma, PNG_DEFAULT_sRGB);
            else
                if png_get_gAMA(png_ptr, info_ptr, image_gamma) != 0
                    png_set_gamma(png_ptr, screen_gamma, image_gamma[])
                else
                    image_gamma[] = 0.45455
                    png_set_gamma(png_ptr, screen_gamma, image_gamma[])
                end
            end
        elseif png_get_valid(png_ptr, info_ptr, PNG_INFO_gAMA) != 0
            if png_get_gAMA(png_ptr, info_ptr, image_gamma) != 0
                png_set_gamma(png_ptr, screen_gamma, image_gamma[])
            else
                image_gamma[] = 0.45455
                png_set_gamma(png_ptr, screen_gamma, image_gamma[])
            end
        end
    elseif gamma != 1
        image_gamma[] = gamma
        png_set_gamma(png_ptr, screen_gamma, image_gamma[])
    end

    buffer_eltype = _buffer_color_type(color_type, bit_depth)
    bit_depth == 16 && png_set_swap(png_ptr)
    n_passes = png_set_interlace_handling(png_ptr)
    png_read_update_info(png_ptr, info_ptr)
    chunk_info(chunk) = png_get_valid(png_ptr, info_ptr, chunk)

    @show(fpath,
        chunk_info(PNG_INFO_sRGB),
        chunk_info(PNG_INFO_iCCP),
        chunk_info(PNG_INFO_gAMA),
        chunk_info(PNG_INFO_tRNS),
        chunk_info(PNG_INFO_cHRM),
        chunk_info(PNG_INFO_PLTE),
        chunk_info(PNG_INFO_sPLT),
        chunk_info(PNG_INFO_hIST),
        chunk_info(PNG_INFO_sBIT),
        chunk_info(PNG_INFO_bKGD),
        chunk_info(PNG_INFO_pCAL),
        (height, width, num_channels),
        (color_type_orig, color_type),
        (bit_depth_orig, bit_depth),
        (image_gamma[], screen_gamma),
        intent[],
        (white_x[], white_y[], red_x[], red_y[], green_x[], green_y[], blue_x[], blue_y[]),
        (interlace_type, n_passes),
        buffer_eltype,
    )

    png_destroy_read_struct(Ref{Ptr{Cvoid}}(png_ptr), Ref{Ptr{Cvoid}}(info_ptr), C_NULL)
    close_png(fp)
end

mutable struct _png_color_16_struct
    index::png_byte
    red::png_uint_16
    green::png_uint_16
    blue::png_uint_16
    gray::png_uint_16
end
_png_color_16_struct(c::AbstractGray{<:AbstractFloat}) = _png_color_16_struct(0x00, gray(Gray{N0f8}(c)).i, gray(Gray{N0f8}(c)).i, gray(Gray{N0f8}(c)).i, gray(Gray{N0f8}(c)).i)
_png_color_16_struct(c::AbstractGray{<:Normed}) = _png_color_16_struct(0x00, gray(c).i, gray(c).i, gray(c).i, gray(c).i)
_png_color_16_struct(c::UInt8) = _png_color_16_struct(c, 0x0000, 0x0000, 0x0000, 0x0000)
function _png_color_16_struct(c::AbstractRGB{<:AbstractFloat})
    n = RGB{N0f8}(c)
    _png_color_16_struct(0x00, red(n).i, green(n).i, blue(n).i, 0x0000)
end
_png_color_16_struct(c::AbstractRGB{<:Normed}) =  _png_color_16_struct(0x00, red(c).i, green(c).i, blue(c).i, 0x0000)

function _png_get_bKGD(png_ptr, info_ptr, background)
    ccall((:png_get_bKGD, libpng), png_uint_32, (png_const_structrp, png_inforp, Ptr{Ptr{_png_color_16_struct}}), png_ptr, info_ptr, background)
end

function _png_set_background(png_ptr, background_color, background_gamma_code, need_expand, background_gamma)
    ccall((:png_set_background, libpng), Cvoid, (png_structrp, Ptr{_png_color_16_struct}, Cint, Cint, Cdouble), png_ptr, background_color, background_gamma_code, need_expand, background_gamma)
end

function process_background(png_ptr, info_ptr, background::Bool)
    if background
        bg = _png_color_16_struct(0, 0, 0, 0, 0)
        GC.@preserve bg begin
            bgpp = Ref(Ptr{_png_color_16_struct}(pointer_from_objref(bg)))
            if _png_get_bKGD(png_ptr, info_ptr, bgpp) != 0
                _png_set_background(png_ptr, bgpp[], PNG_BACKGROUND_GAMMA_FILE, 1, 1.0)
            end
        end
        return bg
    end
    return
end

function process_background(png_ptr, info_ptr, background::Union{AbstractRGB,AbstractGray})
    bg = _png_color_16_struct(background)
    GC.@preserve bg begin
        bgp = Ptr{_png_color_16_struct}(pointer_from_objref(bg))
        _png_set_background(png_ptr, bgp, PNG_BACKGROUND_GAMMA_FILE, 1, 1.0)
        # _png_set_background(png_ptr, bgp, PNG_BACKGROUND_GAMMA_SCREEN, 0, 1.0)
    end
    return bg
end

function process_background(png_ptr, info_ptr, background::UInt8)
    bg = _png_color_16_struct(background)
    GC.@preserve bg begin
        bgp = Ptr{_png_color_16_struct}(pointer_from_objref(bg))
        _png_set_background(png_ptr, bgp, PNG_BACKGROUND_GAMMA_FILE, 1, 1.0)
        # _png_set_background(png_ptr, bgp, PNG_BACKGROUND_GAMMA_SCREEN, 0, 1.0)
    end
    return bg
end

_check_background_load(background::Bool, is_gray, has_pallette, as_paletted) = background
_check_background_load(background::AbstractGray, is_gray, has_pallette, as_paletted) = true
function _check_background_load(background::UInt8, is_gray, has_pallette, as_paletted)
    has_pallette || throw(ArgumentError("Only paletted images can use a paletted index (background::UInt8)."))
    !as_paletted || throw(ArgumentError("In order to use a (non-default) palette index to specify background, the `expand_paletted` flag must be set."))
    return true
end
function _check_background_load(background::AbstractRGB, is_gray, has_pallette, as_paletted)
    is_gray && throw(ArgumentError("Gray Images cannot use RGB background."))
    return true
end

function _check_background_save(background::AbstractGray, color_type)
    color_type == PNG_COLOR_TYPE_GRAY || color_type == PNG_COLOR_TYPE_GRAY_ALPHA ||
        throw(ArgumentError("Only gray scale images can use a Gray for a background."))
end
function _check_background_save(background::AbstractRGB, color_type)
   (color_type & PNG_COLOR_MASK_COLOR) != 0 ||
        throw(ArgumentError("Only true color images can use a RGB for background."))
end
function _check_background_save(background::UInt8, color_type)
    (color_type & PNG_COLOR_MASK_PALETTE) != 0 ||
        throw(ArgumentError("Only paletted images can use a paletted index for background."))
end

function _png_set_bKGD(png_ptr, info_ptr, background)
    bg = _png_color_16_struct(background)
    GC.@preserve bg begin
        bgp = Ptr{_png_color_16_struct}(pointer_from_objref(bg))
        png_set_bKGD(png_ptr, info_ptr, bgp)
    end
    return
end

_adjust_background_bitdepth(x, bit_depth) = x
function _adjust_background_bitdepth(x::AbstractRGB{T}, bit_depth) where {T}
    __get_bit_depth(T) == bit_depth && return x
    return bit_depth == 8 ? RGB{N0f8}(x) : RGB{N0f16}(x)
end
function _adjust_background_bitdepth(x::AbstractGray{T}, bit_depth) where {T}
    __get_bit_depth(T) == bit_depth && return x
    return bit_depth == 8 ? Gray{N0f8}(x) : Gray{N0f16}(x)
end


text:
using PowerSimulationsDynamics
using Sundials

"""
Case 25:
This case study a three-bus system, with two Marconato generators (in buses 1 and 2), and a constant impedance load in bus 3.
The perturbation increase the reference of mechanical power of generator-2 from 0.8 to 0.9 at t=1.0s.
"""

##################################################
############### LOAD DATA ########################
##################################################

include(joinpath(dirname(@__FILE__), "data_tests/test25.jl"))

##################################################
############### SOLVE PROBLEM ####################
##################################################

#time span
tspan = (0.0, 40.0);

#PSCAD benchmark data
csv_file = joinpath(dirname(@__FILE__), "benchmarks/pscad/Test25/Test25_v102.csv")
t_offset = 49.0

#Define Fault using Callbacks
gen2 = get_dynamic_injector(get_component(Generator, sys, "generator-102-1"));
Pref_change = ControlReferenceChange(1.0, gen2, PSID.P_ref_index, 0.9);

@testset "Test 25 Marconato with Dynamic Lines ImplicitModel" begin
    path = (joinpath(pwd(), "test-25"))
    !isdir(path) && mkdir(path)
    try
        #Define Simulation Problem
        sim = Simulation!(ImplicitModel, sys, path, tspan, Pref_change)

        small_sig = small_signal_analysis(sim)
        @test small_sig.stable

        #Solve problem in equilibrium
        execute!(sim, Sundials.IDA(), dtmax = 0.01, saveat = 0.01)

        #Obtain voltage magnitude data
        series = get_voltage_magnitude_series(sim, 102)
        t = series[1]
        v = series[2]

        #Obtain benchmark data from PSCAD
        M = get_csv_data(csv_file)
        t_pscad = M[:, 1] .- t_offset
        v_pscad = M[:, 2]

        diff = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in test25_x0_init
            diff[1] += LinearAlgebra.norm(res[k] - v)
        end
        @test (diff[1] < 1e-3)
        @test sim.solution.retcode == :Success
        #relaxed constraint to account for mismatch in damping
        @test LinearAlgebra.norm(v - v_pscad) <= 0.05
        @test LinearAlgebra.norm(t - round.(t_pscad, digits = 3)) == 0.0
    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end

@testset "Test 25 Marconato with Dynamic Lines MassMatrixModel" begin
    path = (joinpath(pwd(), "test-25"))
    !isdir(path) && mkdir(path)
    try
        #Define Simulation Problem
        sim = Simulation!(MassMatrixModel, sys, path, tspan, Pref_change)

        small_sig = small_signal_analysis(sim)
        @test small_sig.stable

        #Solve problem in equilibrium
        execute!(sim, Rodas5(autodiff = false), dtmax = 0.01, saveat = 0.01)

        #Obtain voltage magnitude data
        series = get_voltage_magnitude_series(sim, 102)
        t = series[1]
        v = series[2]

        #Obtain benchmark data from PSCAD
        M = get_csv_data(csv_file)
        t_pscad = M[:, 1] .- t_offset
        v_pscad = M[:, 2]

        diff = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in test25_x0_init
            diff[1] += LinearAlgebra.norm(res[k] - v)
        end
        @test (diff[1] < 1e-3)
        @test sim.solution.retcode == :Success
        #relaxed constraint to account for mismatch in damping
        @test LinearAlgebra.norm(v - v_pscad) <= 0.05
        @test LinearAlgebra.norm(t - round.(t_pscad, digits = 3)) == 0.0
    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end


text:
## ExoplanetsSysSim/src/ExoplanetsSysSim.jl
## (c) 2015 Eric B. Ford
module ExoplanetsSysSim

# Packages to be used                 # Note: Tried to keep to a minimum in core package to help with maintainability
import Compat: @compat, readstring, is_windows
if VERSION < v"0.7"
  import Compat: UTF8String, ASCIIString
end
if VERSION >= v"0.5-"
  using Combinatorics
end
using Distributions
# using StatsFuns
using DataFrames
if VERSION < VersionNumber(0,3,20)
  using Docile
end

using CORBITS

# Includes & associated  exports for types, then generic functions, then demo functions
include("constants.jl")
#require(joinpath(Pkg.dir("ExoplanetsSysSim"), "src", "constants.jl"))
export SimParam
export SimulationParameters
export add_param_fixed, add_param_active, set_active, set_inactive, is_active, update_param, get, get_real, get_int, get_function, get_any, haskey, make_vector_of_sim_param, update_sim_param_from_vector!, make_vector_of_active_param_keys, get_range_for_sim_param
export setup_sim_param_demo
include("simulation_parameters.jl")
using ExoplanetsSysSim.SimulationParameters
export Orbit
include("orbit.jl")
export Planet
include("planet.jl")
export StarAbstract, Star, SingleStar, BinaryStar
include("star.jl")
export flux, mass
export generate_stars
export PlanetarySystemAbstract, PlanetarySystemSingleStar, PlanetarySystem
include("limb_darkening.jl")
include("planetary_system.jl")
#include("corbits.jl")
#export prob_of_transits_approx
include("window_function.jl")
export WindowFunction
export setup_window_function, get_window_function_data, get_window_function_id, eval_window_function
include("stellar_table.jl")
export  StellarTable
export setup_star_table, star_table, num_usable_in_star_table, set_star_table, star_table_has_key
export KeplerTarget
export num_planets, generate_kepler_target_from_table, generate_kepler_target_simple
include("target.jl")
export KeplerTargetObs
include("transit_observations.jl")
export semimajor_axis
#export setup_koi_table, koi_table, num_koi
export KeplerPhysicalCatalog, KeplerObsCatalog
include("kepler_catalog.jl")
export generate_kepler_physical_catalog, observe_kepler_targets_single_obs, observe_kepler_targets_sky_avg, simulated_read_kepler_observations, setup_actual_planet_candidate_catalog, read_koi_catalog
export CatalogSummaryStatistics
include("summary_statistics.jl")
export calc_summary_stats_obs_demo, calc_summary_stats_sim_pass_one_demo, calc_summary_stats_sim_pass_two_demo
include("abc_distance.jl")
export dist_L1_fractional, dist_L1_abs, dist_L2_fractional, dist_L2_abs, calc_scalar_distance, combine_scalar_distances, distance_poisson_likelihood, distance_poisson_draw, distance_sum_of_bernoulli_draws
export calc_distance_vector_demo
export TestEvalModel
include("eval_model.jl")   # Also includes macros to help write eval model using different variables for closures
export test_eval_model
using ExoplanetsSysSim.TestEvalModel
export SysSimIO
include("io.jl")
using ExoplanetsSysSim.SysSimIO
export save_sim_param, save_sim_results, load_sim_param, load_distances, load_summary_stats

end # module


text:
function GLMakie.convert_arguments(::PointBased, contres::AbstractBranchResult, vars = nothing)
	ind1, ind2 = getPlotVars(contres, vars)
	return ([Point2f0(i, j) for (i, j) in zip(getproperty(contres.branch, ind1), getproperty(contres.branch, ind2))],)
end

function plotBranchCont(contres::ContResult,
		sol::BorderedArray,
		contparms,
		plotuserfunction;
		plotfold = false,
		plotspecialpoints = true,
		filterspecialpoints = false,
		plotcirclesbif = true)
	if length(contres) == 0; return ;end

	ind1, ind2 = getPlotVars(contres, nothing)
	xlab, ylab = getAxisLabels(ind1, ind2, contres)

	fig = Figure(resolution = (1200, 700))
	ax1 = fig[1, 1] = Axis(fig, xlabel = String(xlab), ylabel = String(ylab), tellheight = true)
	lines!(ax1, contres, color = :blue)

	ax2 = fig[1, 2] = Axis(fig, xlabel = "step", ylabel = String(xlab))
	lines!(ax2, contres.step, contres.param, color = :blue)

	if isnothing(plotuserfunction) == false
		ax3, line3 = plotuserfunction(fig[2:3, 1:2], sol.u, sol.p)
	end
	#
	# put arrow to indicate the order of computation
	if length(contres) > 1
		# plot!([contres.branch[end-1:end].param], [getproperty(contres.branch,1)[end-1:end]], label = "", arrow = true, subplot = 1)
		x = contres.branch[end].param
		y = getproperty(contres.branch,1)[end]
		# u = contres.branch[end].param - contres.branch[end-1].param
		# v = getproperty(contres.branch,1)[end]-getproperty(contres.branch,1)[end-1]
		scatter!(ax1, [x], [y], marker = :cross )
	end

	# display bifurcation points
	bifpt = filter(x -> (x.type != :none) && (plotfold || x.type != :fold), contres.specialpoint)
	if length(bifpt) >= 1 && plotspecialpoints #&& (ind1 == :param)
		if filterspecialpoints == true
			bifpt = filterBifurcations(bifpt)
		end
		scatter!(ax1, map(x -> getproperty(x, ind1), bifpt), map(x -> getproperty(x.printsol, ind2), bifpt), marker = map(x -> (x.status == :guess) && (plotcirclesbif==false) ? :rect : :circle, bifpt), markersize = 7, color = map(x -> colorbif[x.type], bifpt))
	end

	display(fig)
	fig
end

function plotBranch(contres::AbstractBranchResult; plotfold = false, plotspecialpoints = true, filterspecialpoints = false, plotcirclesbif = true, linewidthunstable = 1.0, linewidthstable = 2linewidthunstable)
	if length(contres) == 0; return ;end

	# Special case labels when vars = (:p,:y,:z) or (:x) or [:x,:y] ...
	ind1, ind2 = getPlotVars(contres, vars)
	xlab, ylab = getAxisLabels(ind1, ind2, contres)

	fig = Figure(resolution = (1200, 700))
	ax1 = fig[1, 1] = Axis(fig, xlabel = String(xlab), ylabel = String(ylab), tellheight = true)

	#TODO use band to have similar  result to Plots
	colors = hasstability(contres) ? map(x -> isodd(x) ? :green : :red, contres.stable) : :black
	lines!(ax1, contres, color = colors)

	#
	# put arrow to indicate the order of computation
	if length(contres) > 1
		x = xplot[end]
		y = yplot[end]
		scatter!(ax1, [x], [y], marker = :cross )
	end

	# display bifurcation points
	bifpt = filter(x -> (x.type != :none) && (plotfold || x.type != :fold) && (x.idx <= length(contres)-1), contres.specialpoint)

	if length(bifpt) >= 1 && plotspecialpoints #&& (ind1 == :param)
		if filterspecialpoints == true
			bifpt = filterBifurcations(bifpt)
		end
		scatter!(ax1, [applytoX(getproperty(contres[pt.idx], ind1)) for pt in bifpt],
					  [applytoY(getproperty(contres[pt.idx], ind2)) for pt in bifpt],
			 marker = map(x -> (x.status != :converged) && (plotcirclesbif==false) ? :rect : :circle, bifpt), markersize = 7, color = map(x -> colorbif[x.type], bifpt))
	end

	display(fig)
	fig
end

function plotBranch(contresV::AbstractBranchResult...; plotfold = false, plotspecialpoints = true, filterspecialpoints = false, plotcirclesbif = true, linewidthunstable = 1.0, linewidthstable = 2linewidthunstable)
	if length(contresV) == 0; return ;end
	fig = Figure(resolution = (1200, 700))
	ax1 = fig[1, 1] = Axis(fig)

	for contres in contresV

		ind1, ind2 = getPlotVars(contres, nothing)
		xlab, ylab = getAxisLabels(ind1, ind2, contres)

		colors = hasstability(contres) ? map(x -> isodd(x) ? :green : :red, contres.stable) : :black
		# lines!(ax1, contres, color = colors)

		h = lines!(ax1, contres, color = colors)
		ax1.xlabel = xlab; ax1.ylabel = ylab

		# display bifurcation points
		bifpt = filter(x -> (x.type != :none) && (plotfold || x.type != :fold), contres.specialpoint)
		if length(bifpt) >= 1 && plotspecialpoints #&& (ind1 == :param)
			if filterspecialpoints == true
				bifpt = filterBifurcations(bifpt)
			end
			scatter!(ax1, map(x -> getproperty(x, ind1), bifpt), map(x -> getproperty(x.printsol, ind2), bifpt), marker = map(x -> (x.status == :guess) && (plotcirclesbif==false) ? :rect : :circle, bifpt), markersize = 7, color = map(x -> colorbif[x.type], bifpt))
		end
	end

	display(fig)
	fig
end


text:
"""
    stream_velocity(r, horseshoes, Γs, V, Ω)

Evaluate the total induced velocity at a point `r` given Horseshoes, vortex strengths `Γ`s, rotation rates `Ω`, and freestream flow vector `freestream` in the global reference frame.
"""
stream_velocity(r, horseshoes, Γs, V, Ω) = sum(x -> velocity(r, x[1], x[2], V / norm(V)), zip(horseshoes, Γs)) + V + Ω × r

"""
    streamlines(point, freestream :: Freestream, horseshoes, Γs, length, num_steps)

Compute the streamlines from a given starting point, a Freestream, Horseshoes and their associated strengths Γs with a specified length of the streamline and number of evaluation points.
"""
function streamlines(V, Ω, horseshoes, Γs, point, length, num_steps :: Integer)
    streamlines = fill(point, num_steps)
    f(x) = stream_velocity(x, horseshoes, Γs, V, Ω)
    for i ∈ 2:num_steps
        update = f(streamlines[i-1])
        streamlines[i] = streamlines[i-1] + (update / norm(update) * length / num_steps)
    end
    streamlines
end

streamlines(fs :: Freestream, refs :: References, vortices, Γs, points, length, num_steps :: Integer) = mapreduce(pt -> streamlines(refs.speed * velocity(fs), fs.omega, vortices, Γs, pt, length, num_steps), hcat, points)

streamlines(system :: VortexLatticeSystem, points, length, num_steps :: Integer) = streamlines(system.freestream, system.reference, system.vortices, system.circulations, points, length, num_steps)


text:
using Test

@testset "MarketData.jl" begin
  include("utils.jl")
  include("const.jl")
  include("downloads.jl")
end


text:
"""
    fr = metagauss(f)

Generate a metagaussian field multiplying the spectralfield `f` with random phases and performing an inverse FFT transform to real space.
"""
function metagauss(f)

    (ns, ns, nt) = size(f)
    #	phases as fft of a gaussian noise random field
    #   ph=zeros(ns,ns,nt,'single')

    ph = zeros(ns, ns, nt)
    for i = 1:nt
        ph[:,:,i] = randn(ns,ns)
    end
    ph = fft(ph)
    ph = ph./abs.(ph); ph[1,1,1] = 0
    ph = f.*ph
    ph = real(ifft(ph))
end


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

function genpotential(foldername, filename, u, xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, pot)

n = length(u);
gen = 1;
ifile = 0;
stropu = "";
strcpu = "";
strgpu = "";
potnum = 0;

for i = 1:n
    ui = u[i];
    if ui != 0
        potnum = potnum + 1;
        fn = filename * string(potnum);
        if pot==1
            tmpopu, tmpcpu, tmpgpu = gensingle(foldername, fn, ui, xi, qi, ti, ai, mu, eta, kappa, gen, ifile);
        elseif pot==2
            tmpopu, tmpcpu, tmpgpu = genpair(foldername, fn, ui, xij, qi, qj, ti, tj, ai, aj, mu, eta, kappa, gen, ifile);
        elseif pot==3
            tmpopu, tmpcpu, tmpgpu = gentriplet(foldername, fn, ui, xij, xik, qi, qj, qk, ti, tj, tk, ai, aj, ak, mu, eta, kappa, gen, ifile);
        elseif pot==4
            tmpopu, tmpcpu, tmpgpu = genquadruplet(foldername, fn, ui, xij, xik, xil, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, gen, ifile);
        end
        tmpopu = modifystr1(tmpopu);
        tmpcpu = modifystr1(tmpcpu);
        tmpgpu = modifystr1(tmpgpu);
        stropu = stropu * tmpopu * "\n";
        strcpu = strcpu * tmpcpu * "\n";
        strgpu = strgpu * tmpgpu * "\n";
    end
end

if potnum==0
    nopotential(foldername, filename, u, xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, pot);
else
    sp0, sp1, sp2, sp3 = getpotstr(pot);
    sp0 = replace(sp0, "int ng)" => "int ng, int potnum)");
    sp1 = replace(sp1, "int);" => "int, int);");
    sp2 = replace(sp2, "int);" => "int, int);");
    tmpopu, tmpcpu, tmpgpu = genpotnum(filename, potnum, sp0, sp1, sp2, sp3);
    stropu = stropu * tmpopu;
    strcpu = strcpu * tmpcpu;
    strgpu = strgpu * tmpgpu;

    # here
    if pot==1
        tmpgpu = replace(tmpgpu, "(T *u," => "Gradient(T *u, T *du, T *u_xi,");
        tmpgpu = replace(tmpgpu, "(u," => "Gradient(u, du, u_xi,");
        tmpgpu = replace(tmpgpu, "(double *" => "Gradient(double *, double *, double*");
        tmpgpu = replace(tmpgpu, "(float *" => "Gradient(float *, float *, float*");
    elseif pot==2
        tmpgpu = replace(tmpgpu, "(T *u," => "Gradient(T *u, T *du, T *u_xij,");
        tmpgpu = replace(tmpgpu, "(u," => "Gradient(u, du, u_xij,");
        tmpgpu = replace(tmpgpu, "(double *" => "Gradient(double *, double *, double*");
        tmpgpu = replace(tmpgpu, "(float *" => "Gradient(float *, float *, float*");
    elseif pot==3
        tmpgpu = replace(tmpgpu, "(T *u," => "Gradient(T *u, T *du, T *u_xij, T *u_xik,");
        tmpgpu = replace(tmpgpu, "(u," => "Gradient(u, du, u_xij, u_xik,");
        tmpgpu = replace(tmpgpu, "(double *" => "Gradient(double *, double *, double*, double*");
        tmpgpu = replace(tmpgpu, "(float *" => "Gradient(float *, float *, float*, float*");
    elseif pot==4
        tmpgpu = replace(tmpgpu, "(T *u," => "Gradient(T *u, T *du, T *u_xij, T *u_xik, T *u_xil,");
        tmpgpu = replace(tmpgpu, "(u," => "Gradient(u, du, u_xij, u_xik, u_xil,");
        tmpgpu = replace(tmpgpu, "(double *" => "Gradient(double *, double *, double*, double*, double*");
        tmpgpu = replace(tmpgpu, "(float *" => "Gradient(float *, float *, float*, float*, float*");
    end
    #strgpu = strgpu * "\n" * tmpgpu;
    strgpu = strgpu * "\n" * "#ifdef _ENZYME\n" * tmpgpu * "#endif\n";

    cppfiles(foldername, filename, stropu, strcpu, strgpu, 0);
end

end

function modifystr1(oldstr)

newstr = oldstr;
s = "u[" * string(0) * " + i*" * string(1) * "]";
if contains(newstr, s)
    news = "u[i]";
    newstr = replace(newstr, s => news);
end

return newstr

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP