text:
struct GTRrel <: GTR
  α::Float64
  β::Float64
  γ::Float64
  δ::Float64
  ϵ::Float64
  πA::Float64
  πC::Float64
  πG::Float64
  πT::Float64
  function GTRrel(α::Float64, β::Float64, γ::Float64,
                  δ::Float64, ϵ::Float64,
                  πA::Float64, πC::Float64, πG::Float64, πT::Float64)
    if α <= 0.
      error("GTR parameter α must be positive")
    elseif β <= 0.
      error("GTR parameter β must be positive")
    elseif γ <= 0.
      error("GTR parameter γ must be positive")
    elseif δ <= 0.
      error("GTR parameter δ must be positive")
    elseif ϵ <= 0.
      error("GTR parameter ϵ must be positive")
    elseif sum([πA,πC,πG,πT]) != 1.0
      error("GTR frequencies must sum to 1.0")
    elseif any([πA,πC,πG,πT] .<= 0.0)
      error("GTR frequencies must be positive")
    end
    new(α, β, γ, δ, ϵ, πA, πC, πG, πT)
  end
end

function show(io::IO, object::GTRrel)
  print(io, "\r\e[0m\e[1mG\e[0meneralised \e[1mT\e[0mime \e[1mR\e[0meversible model (relative rate form)
α = $(object.α), β = $(object.β), γ = $(object.γ), δ = $(object.δ), ϵ = $(object.ϵ), π = [$(object.πA), $(object.πC), $(object.πG), $(object.πT)]")
end

GTR(α, β, γ, δ, ϵ, πA, πC, πG, πT) = GTRrel(α, β, γ, δ, ϵ, πA, πC, πG, πT)

const _α(mod::GTRrel) = mod.α
const _β(mod::GTRrel) = mod.β
const _γ(mod::GTRrel) = mod.γ
const _δ(mod::GTRrel) = mod.δ
const _ϵ(mod::GTRrel) = mod.ϵ
const _η(mod::GTRrel) = 1.0


text:
"""
    kneaddata(inputfile, outputfile; kwargs...)

Run `kneaddata` command line tool on `inputfile`,
creating `outputfile`.
Requires `kneaddata` to be installed and accessible in the `PATH`
(see [Getting Started](@ref)).

`kneaddata` options can be passed via keyword arguments.
For example, if on the command line you would run:

```sh
\$ kneaddata -i some.fastq.gz -o test --n 8 --bypass-trim
```

using this function, you would write:

```julia
kneaddata("some.fastq.gz", "test"; n = 8, bypass_trim=true)
```

To pass multiple databases,
pass an array of paths to the `reference_db` argument

Conda installations of `trimmomatic` (a dependency of `kneaddata`)
don't work properly out of the box.
If you have installed `kneaddata` using commandline conda
(instead of `Conda.jl`),
use `trimmomatic = /path/to/trimmomatic`,
where `/path/to/trimmomatic` is something like
`/home/username/miniconda3/envs/biobakery3/share/trimmomatic`.
If you used `BiobakeryUtils.install_deps()`,
you don't need to worry about this.
"""
function kneaddata(inputs, output; kwargs...)
    check_for_install("kneaddata")
    c = ["kneaddata"]
    if inputs isa AbstractString
        append!(c, ["-i", inputs])
    else
        length(inputs) > 2 && error("for now, only 2 input files possible")
        append!(c, ["-i", inputs[1], "-i", inputs[2]])
    end

    append!(c, ["-o", output])

    if !haskey(kwargs, :trimmomatic)
        env = get(kwargs, :condaenv, :BiobakeryUtils)
        trimpath = normpath(joinpath(Conda.bin_dir(env), "..", "share"), "trimmomatic")

        append!(c, ["--trimmomatic", trimpath])
    end

    if haskey(kwargs, :reference_db)
        if kwargs[:reference_db] isa String
            append!(c, ["-db", kwargs[:reference_db]])
        else
            for db in kwargs[:reference_db]
                append!(c, ["-db", db])
            end
        end
    end

    add_cli_kwargs!(c, kwargs; optunderscores=false, skip=[:reference_db])

    @info "Running command: $(Cmd(c))"
    return run(Cmd(c))
end

"""
    kneaddata_database(db, kind, path)

See `kneaddata_database --help`

```julia
kneaddata_database("human_genome", "bowtie2", "/some/database/dir/")
```
"""
function kneaddata_database(db, kind, path)
    check_for_install("kneaddata_database")
    c = ["kneaddata_database", "--download", db, kind, path]

    @info "Running command: $(Cmd(c))"
    return run(Cmd(c))
end

"""
    kneaddata_read_count_table(input, output)

See `kneaddata_read_count_table --help`

```julia
kneaddata_read_count_table("human_genome", "bowtie2", "/some/database/dir/")
```
"""
function kneaddata_read_count_table(input, output)
    check_for_install("kneaddata_read_count_table")
    c = ["kneaddata_read_count_table", "--input", input, "--output", output]

    @info "Running command: $(Cmd(c))"
    return run(Cmd(c))
end

ext: jl

max_stars_repo_name: BioJulia/BiobakeryUtils

max_issues_repo_name: BioJulia/BiobakeryUtils

max_forks_repo_name: BioJulia/BiobakeryUtils


text:
using Random
using Printf

include("Algencan/AlgencanWrapper.jl")
include("sh.jl")

import Voronoi
import Covering

drawing = (haskey(ENV, "DRAW") && ENV["DRAW"] != "" && ENV["DRAW"] != "0")
drawingAll = (haskey(ENV, "DRAWALL") && ENV["DRAWALL"] != "" && ENV["DRAWALL"] != "0")

WIDTH = 1.0
HEIGHT = 1.0

A = Array{Covering.Section, 1}([
        Covering.Section(Array{Tuple{Real, Real}, 1}([
            (0.2, 0.5),
            (0.8, 0.5),
            (0.8, 0.75),
            (0.2, 0.75),
        ])),
        Covering.Section(Array{Tuple{Real, Real}, 1}([
            (0.4, 0.5),
            (0.4, 0.25),
            (0.6, 0.25),
            (0.6, 0.5),
        ]))
    ])

#A = Array{Covering.Section, 1}([
#        Covering.Section(Array{Tuple{Real, Real}, 1}([
#            (0, 0.125),
#            (0.75, 0.5),
#            (1 ,0.875),
#            (0.25, 0.5),
#        ])),
#        Covering.Section(Array{Tuple{Real, Real}, 1}([
#            (0, 1),
#            (0, 0.5),
#            (0.5, 1),
#        ]))
#    ])

function f()
    total_area = 0
    for Aⱼ in A
        t, _, _ = Covering.areaAndGradient(Aⱼ, 0)
        total_area += t
    end
    return total_area
end
total_area = f()

#function avoidRepeats(points)
#    println("LÁ VAMOS NÓS")
#    println("points antes: ", points)
#    repeat = []
#    for i in 1:length(points)
#        for j in 1:length(points)
#            if j == i
#                continue
#            end
#
#            if points[i][1] == points[j][1] && points[i][2] == points[j][2]
#                push!(repeat, points[j])
#                deleteat!(points, j)
#            end
#        end
#    end
#
#    println("points depois: ", points)
#    println("repeats depois: ", repeat)
#    return points, repeat
#end

function coverWithCircles(n, WIDTH, HEIGHT, r, points)
   function pack(r, points)
       x = zeros(2n + 1)
       x[1] = r
       for (i, p) in enumerate(points)
           x[2i] = p[1]
           x[2i + 1] = p[2]
       end

       return x
   end

   function unpack(x)
       r = x[1]

       points = Array{Tuple{Real, Real}}([])

       for i in 1:n
           append!(points, [(x[2i], x[2i + 1])])
       end

       return r, points
   end

   function draw(r, points)
       if !drawing
           return
       end

       V = Voronoi.Fortune.compute(points, WIDTH, HEIGHT)
       Voronoi.Intersect.intersect(V, Voronoi.Intersect.Rectangle(WIDTH, HEIGHT))

       Draw.init(WIDTH, HEIGHT)
       Draw.circles(r, points)

       for Aⱼ in A
           W = SH.intersect(V, Aⱼ, r)
           Draw.coveringPartition(W, "xkcd:pink")
           Draw.filledPolygon(Aⱼ, "xkcd:black")
       end

       Draw.commit()
       print("")
   end

   function f(x)
       r, points = unpack(x)

       if drawing && drawingAll
           draw(r, points)
       end

       #println("R: ", r)
       return r
   end

   function ∇f(x)
       grad = zeros(2n + 1)
       grad[1] = 1
       #println("GRAD: ", grad)
       return grad
   end

   function ∇∇f(x)
       return [], [], []
   end

   function c(ind, x)
       r, points = unpack(x)

       #points, repeats = avoidRepeats(points)

       V = Voronoi.Fortune.compute(points, WIDTH, HEIGHT)
       Voronoi.Intersect.intersect(V, Voronoi.Intersect.Rectangle(WIDTH, HEIGHT))

       uncovered_area = total_area
       for Aⱼ in A
           W = SH.intersect(V, Aⱼ, r)
           area, _, _ = Covering.areaAndGradient(W)
           uncovered_area -= area
       end

       #println("C: ", uncovered_area)
       return uncovered_area
   end

   function ∇c(ind, x)
       r, points = unpack(x)

       #points, repeats = avoidRepeats(points)

       V = Voronoi.Fortune.compute(points, WIDTH, HEIGHT)
       Voronoi.Intersect.intersect(V, Voronoi.Intersect.Rectangle(WIDTH, HEIGHT))

       gᵣ = 0
       gₛ = fill((0.0, 0.0), n)
       for Aⱼ in A
           W = SH.intersect(V, Aⱼ, r)
           _, gᵣⱼ, gₛⱼ = Covering.areaAndGradient(W)
           gᵣ += gᵣⱼ
           for i in 1:n
               gₛ[i] = (gₛ[i][1] + gₛⱼ[i][1], gₛ[i][2] + gₛⱼ[i][2])
           end
       end

       #for i in repeats
       #    push!(gₛ, (0, 0))
       #end
       #println("GRAD C.: ", pack(gᵣ, gₛ))
       #println("Aperte Enter pra continuar")
       #readline(stdin)

       return collect(1:2n+1), pack(gᵣ, gₛ)
   end

   l = zeros(2n + 1)
   l[1] = 1e-14

   u = zeros(2n + 1)
   u[1] = 10e20
   for i in 2:2n+1
       u[i] = 1.0
   end

   draw(r, points)
   #println("Aperte Enter pra continuar")
   #readline(stdin)

   try
       x, fx = AlgencanWrapper.optimize(n = 2n + 1, m = 1,

                                        f = f, g = ∇f,# h = ∇∇f,
                                        equatn = [1],
                                        c = c, jac = ∇c,

                                        x = pack(r, points),
                                        l = l,
                                        u = u,

                                        nvparam = 2,
                                        vparam = [
                                                  "ITERATIONS-OUTPUT-DETAIL 10",
                                                  "PENALTY-PARAMETER-INITIAL-VALUE 1000"
                                                 ],
                                        #checkder = 1,
                                        epsopt = 1.0e-6,
                                        epsfeas = 1.0e-6
                                        )

       r, points = unpack(x)
       draw(r, points)

       println("\n\n=================\nRESULTADO\n=================\n")
       println("r = ", r)
       println("s = ", 1/r)
       println("pontos: ")
       for p in points
           println("(", p[1], ", ", p[2], ")")
       end
       if drawing
           println("Digite ENTER para sair.")
           readline(stdin)
       end
   catch exc
       println("Exception: ", exc)
       return 1e10, []
   end
end

if length(ARGS) == 0
    println("USO: cover.sh N_DE_PONTOS [DESENHAR] [DESENHAR_TODA_AVALIAÇÃO]")
else
    if drawing
        import Draw
    end

    n = tryparse(Int64, ARGS[1])

    function randf(start, finish, n)
        v = rand(n)
        return map(x -> start + x*(finish-start), v)
    end

    Random.seed!(1)
    points = convert(Array{Tuple{Real, Real}}, collect(zip(randf(1, WIDTH-1, n), randf(1, HEIGHT-1, n))))

    r = (0.5 + (rand(1)[1]))/n

    function insideHalfSpace(point, p, q)
        pq = Voronoi.Geometry.subVector(q, p)
        v = Voronoi.Geometry.rotateVectorCCW(pq)

        function dot(a, b)
            return a[1]*b[1] + a[2]*b[2]
        end

        u = Voronoi.Geometry.subVector(point, p)

        return (dot(u, v) >= 0)
    end

    function pointIsInsideConvexPolygon(point, polygon)
        p = polygon.borderHead

        first = true
        while p != polygon.borderHead || first
            first = false
            q = p.next
            if !insideHalfSpace(point, p.origin, q.origin)
                return false
            end

            p = p.next
        end

        return true
    end

    function pointIsInsideSomeConvexPolygon(point, polygons)
        for polygon in polygons
            if pointIsInsideConvexPolygon(point, polygon)
                return true
            end
        end

        return false
    end

    for (i, p) in enumerate(points)
        while !pointIsInsideSomeConvexPolygon(points[i], A)
            #println("Ponto $i ainda está fora. Re-sorteando.")
            points[i] = (randf(1, WIDTH-1, 1)[1], randf(1, HEIGHT-1, 1)[1])
        end
    end

    try
        coverWithCircles(n, WIDTH, HEIGHT, r, points)
    catch e
        @error "Deu pau:" exception=(e, catch_backtrace())
    end
end

max_stars_repo_name: arthursantana/coveringWithCircles

max_issues_repo_name: arthursantana/coveringWithCircles

max_forks_repo_name: arthursantana/coveringWithCircles


text:
# he linearized forward modeling operator `F`, such that
# `Fx` can be computed without explicitly storing the operator matrix (
#  see `LinearMaps.jl`);
# the imaging/migration operator `F*`;

# These maps are the building blocks of iterative optimization schemes.

for scenario in [:downhole, :pizza]
	println("@@@@@@@@@@@@TESTING ", scenario)
	for rfields in [[:P], [:Vx], [:Vz]]
		pa=JG.SeisInvExpt(scenario, born_flag=true, rfields=rfields)

		F=JF.operator_Born(pa);

		x1=randn(size(F,2))
		x2=randn(size(F,2))
		x12=x1.+x2

		d12=F*x12
		d1=F*x1
		δmodtt1=copy(pa.paf.c.δmodtt)
		d2=F*x2
		δmodtt2=copy(pa.paf.c.δmodtt)

		d12new=d1.+d2

		f=Misfits.error_squared_euclidean!(nothing, d12, d12new, nothing, norm_flag=true)

		@test f<1e-25

		function adjtest()
			x=randn(size(F,2))
			y=randn(size(F,1))
			a=LinearAlgebra.dot(y,F*x)
			b=LinearAlgebra.dot(x,adjoint(F)*y)
			c=LinearAlgebra.dot(x, transpose(F)*F*x)
			println("adjoint test: ", a, "\t", b)
			@test isapprox(a,b,rtol=1e-5)
			println("must be positive: ", c)
			@test c>0.0
		end

		for i in 1:3
			adjtest()
		end
	end
end


text:
#=
# Description of Market data model:

1. The highest level are the market data containers, e.g. `FXForwardMarketDataContainer <: AbstractMarketDataContainer`. They store all market data for a whole product class (e.g. FX Forwards) and for the whole time series.
2. The 2nd level are the product-specific market data, e.g. `SpecificFXForwardMarketData <: AbstractProductSpecificMarketData`. This contains the market data for a specific product (here currency pairs) for the whole time series.
The function `get_specific_market_data(product, container)` takes a container (level 1) and returns a product-specific market data item. It is rather slow and should not be called too often.

3. The lowest level are the specific curves. They are stored as OrderedDicts, with the valuation_dates as keys, in the product-specific market data structures. Example: `market_data.fx_forward_curves[valuation_date]`.

The idea of this data model is that it should be fast to price a whole time series of a product, either for studying price evolution / backtesting P&Ls (using real market data time series) or for market risk scenario P&Ls (using and artificial time series).

=#

## Market Data time series

abstract type AbstractTimeSeries end
Base.Broadcast.broadcastable(q:: AbstractTimeSeries) = Ref(q) # treat it as a scalar in broadcasting

Base.getindex(ts:: AbstractTimeSeries, valuation_date:: Date) = ts.curves[valuation_date]
Base.length(ts:: AbstractTimeSeries) = length(ts.curves)
Base.keys(ts:: AbstractTimeSeries) = keys(ts.curves)

struct FXForwardTimeSeries{BaseCurrency <: Currency, QuoteCurrency <: Currency, C <: Curve} <: AbstractTimeSeries
    curves:: SortedDict{Date, C}
end

FXForwardTimeSeries(::BaseCurrency, ::QuoteCurrency, dates:: AbstractVector{Date}, curves:: AbstractVector{C}) where {BaseCurrency, QuoteCurrency, C} =
    FXForwardTimeSeries{BaseCurrency, QuoteCurrency, C}(SortedDict(dates .=> curves))

struct DiscountTimeSeries{BaseCurrency <: Currency, C <: Curve} <: AbstractTimeSeries
    curves:: SortedDict{Date, C}
end

DiscountTimeSeries(::BaseCurrency, dates:: AbstractVector{Date}, curves:: AbstractVector{C}) where {BaseCurrency, C} =
    DiscountTimeSeries{BaseCurrency, C}(SortedDict(dates .=> curves))

## Market data definitions per product

# The product specific market data contains only market data required for valuation of the products.
# No dynamic market data selection (e.g. for currencies, rate indices, etc.) required anymore.

abstract type AbstractProductSpecificMarketData end
Base.Broadcast.broadcastable(q:: AbstractProductSpecificMarketData) = Ref(q) # treat it as a scalar in broadcasting

"""
Container for the market data required to price FX Forwards for a specific currency pair.
The currency pair information is encoded in the type signature so that it is available at compile time.
"""
struct SpecificFXForwardMarketData{BaseCurrency <: Currency, QuoteCurrency <: Currency, C1 <: FXForwardTimeSeries, C2 <: DiscountTimeSeries} <: AbstractProductSpecificMarketData
    fx_forward_curves:: C1
    discount_curves:: C2
    flip_fx_direction:: Bool
end

SpecificFXForwardMarketData(::BaseCurrency, ::QuoteCurrency, fx_forward_curve:: C1, base_ccy_discount_curve:: C2; flip_fx_direction=false) where {BaseCurrency, QuoteCurrency, C1, C2} =
    SpecificFXForwardMarketData{BaseCurrency, QuoteCurrency, C1, C2}(fx_forward_curve, base_ccy_discount_curve, flip_fx_direction)

## Market data container for a whole product class (e.g. FX Forwards)

abstract type AbstractMarketDataContainer end
Base.Broadcast.broadcastable(q:: AbstractMarketDataContainer) = Ref(q) # treat it as a scalar in broadcasting

"""
Check if the time series are consistent for all market data
"""
function _check_scenario_consistency(x:: X) where {X <: AbstractMarketDataContainer}
    ref_length = -1
    local ref_dates
    for field in fieldnames(X)
        md_dict = getproperty(x, field)
        for md_series in values(md_dict)

            if ref_length == -1
                ref_length = length(md_series)
                ref_length > 0 || error("time series is empty")
                ref_dates = collect(keys(md_series))
            else
                length(md_series) == ref_length || error("inconsistent time series lengths")
                ref_dates == collect(keys(md_series)) || error("inconsistent scenario dates")
            end
        end
    end
end

struct FXForwardMarketDataContainer <: AbstractMarketDataContainer
    fx_forward_curves:: Dict{Tuple{Currency, Currency}, FXForwardTimeSeries}
    discount_curves:: Dict{Currency, DiscountTimeSeries}

    function FXForwardMarketDataContainer(fx_forward_curves, discount_curves)
        container = new(fx_forward_curves, discount_curves)
        _check_scenario_consistency(container)
        return container
    end
end

function FXForwardMarketDataContainer(fx_forward_curves:: AbstractVector{<: FXForwardTimeSeries}, discount_curves:: AbstractVector{<: DiscountTimeSeries})
    fx_forward_curve_dict = Dict{Tuple{Currency, Currency}, FXForwardTimeSeries}()
    for curve in fx_forward_curves
        base_currency = typeof(curve).parameters[1]()
        quote_currency = typeof(curve).parameters[2]()
        fx_forward_curve_dict[(base_currency, quote_currency)] = curve
    end

    discount_curve_dict = Dict{Currency, DiscountTimeSeries}()
    for curve in discount_curves
        base_currency = typeof(curve).parameters[1]()
        discount_curve_dict[base_currency] = curve
    end
    return FXForwardMarketDataContainer(fx_forward_curve_dict, discount_curve_dict)
end

get_spot_rates(market_data_container:: FXForwardMarketDataContainer) = Dict(k => interpolate(t"2D", last(v.curves).second) for  (k, v) in pairs(market_data_container.fx_forward_curves))
get_spot_rates(market_data_container:: FXForwardMarketDataContainer, valuation_date:: Date) = Dict(k => interpolate(t"2D", v.curves[valuation_date]) for  (k, v) in pairs(market_data_container.fx_forward_curves))

get_dates(market_data_container:: FXForwardMarketDataContainer) = collect(keys(first(market_data_container.fx_forward_curves).second.curves))

"""
Returns the market data items required for valuation of a specific product
"""
function get_specific_market_data(product:: FXForward, container:: FXForwardMarketDataContainer)
    quote_currency = product.quote_currency
    base_currency = typeof(product.notional_base_currency).parameters[2]()

    discount_curve_time_series = container.discount_curves[quote_currency]
    if (base_currency, quote_currency) ∈ keys(container.fx_forward_curves)
        fx_time_series = container.fx_forward_curves[(base_currency, quote_currency)]
        flip_fx_direction = false
    elseif (quote_currency, base_currency) ∈ keys(container.fx_forward_curves)
        fx_time_series = container.fx_forward_curves[(quote_currency, base_currency)]
        flip_fx_direction = true
    else
        error("currency pair $((quote_currency, base_currency)) not found")
    end

    return SpecificFXForwardMarketData(base_currency, quote_currency, fx_time_series, discount_curve_time_series; flip_fx_direction=flip_fx_direction)
end

max_issues_repo_name: lungben/MarketRiskPrototype

max_forks_repo_name: lungben/MarketRiskPrototype


text:
"""
Helper infrastructure to compile and sample models using `cmdstan`.

[`StanModel`](@ref) wraps a model definition (source code), while [`stan_sample`](@ref) can
be used to sample from it.

[`stan_compile`](@ref) can be used to pre-compile a model without sampling. A
[`StanModelError`](@ref) is thrown if this fails, which contains the error messages from
`stanc`.
"""
module StanSample

using Reexport

@reexport using StanBase

using DocStringExtensions: FIELDS, SIGNATURES, TYPEDEF

import StanRun: stan_sample, stan_cmd_and_paths, default_output_base
import StanSamples: read_samples
import StanBase: cmdline

include("stanmodel/sample_types.jl")
include("stanmodel/SampleModel.jl")
include("stanrun/cmdline.jl")
include("stanrun/stan_generate_quantities.jl")
include("stansamples/read_samples.jl")

export
  SampleModel,
  stan_generate_quantities

end # module


text:
struct BlankGeometry <: Gadfly.GeometryElement
end

"""
    Geom.blank

A blank geometry is drawn, and guides maybe drawn if aesthetics are provided.
"""
blank = BlankGeometry

element_aesthetics(::BlankGeometry) = [:x, :y, :size, :color, :shape]

function render(geom::BlankGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    ctx = context()
    return compose!(ctx, svgclass("geometry"))
end


text:
# NOT CURRENTLY WORKING !!! NEED TO ADD TIME STEPPING ROUTINES TO eFEM

using eFEMpart, JLD

using LinearAlgebra # for condition number calculation

function main()
  function computeNorms(dt,i)
    start = time()
    mesh = squareMesh([-2,2.0,-1.,1.],16,2);
    Tfinal = 0.5
    varname      = "temperature"
    OperatorType = :Poisson2D     # -Delta u = f

    # Define auxiliary information
    dNodes = Dirichlet(:left,:right)
    nNodes = Neumann(:bottom,:top)
    dBCf = Dirichlet((x,y) -> 0.0)
    nBCf = Neumann((x,y) -> 0.0)

    # required technicality for grouping auxiliary information
    Nodes = [dNodes,nNodes]; bcfun = [dBCf,nBCf]

    # define problem variable (helps solver)
    prob = Problem(mesh,Nodes,bcfun,varname,OperatorType)

    param = PoissonParam(1.0)
    xic = [mesh.xy[i].x for i=1:length(mesh.xy)]
    U0  = sin.(pi*xic/2)

    U = FEMForwardEuler(prob,mesh,param,U0,Tfinal,dt)
    sol = ScalarSolution(U)

    #fn = Path("HeatValidation$(i)")
    # print problem
    #vtksave(mesh,prob,param,sol,fn)

    elapsed = time()-start

    # generate exact solution
    Uexact(x,y) = exp(-(pi/2)^2*Tfinal)*sin(pi*x/2)
    UexactArr = [Uexact.(mesh.xy[i].x,mesh.xy[i].y) for i=1:length(mesh.xy)]

    # compute difference
    err = sol.u - UexactArr

    # compute Li norm
    L1   = DomainNorm(mesh.xy,mesh.cm,err;normID="1")
    L2   = DomainNorm(mesh.xy,mesh.cm,err;normID="2")
    Linf = DomainNorm(mesh.xy,mesh.cm,err;normID="Inf")

    return L1,L2,Linf,elapsed
  end

  start = 1e-4

  dtarr    = [start,start,start/2,start/4,start/8,start/16]
  N       = length(dtarr)
  L1arr   = zeros(N)
  L2arr   = zeros(N)
  LInfarr = zeros(N)
  timearr = zeros(N)

  for i=1:N
      L1arr[i],L2arr[i],LInfarr[i],timearr[i] = computeNorms(dtarr[i],i)
      println("completed i=",i,"/",N)
  end

  L1arr = L1arr[2:end]
  L2arr = L2arr[2:end]
  LInfarr = LInfarr[2:end]
  timearr = timearr[2:end]
  dtarr = dtarr[2:end]

  run(`mkdir -p data`)

  save("data/TEMPHeat_validation.jld", "L1arr", L1arr,
                                        "L2arr", L2arr,
                                        "LInfarr",LInfarr,
                                        "timearr",timearr,
                                        "dtarr",dtarr)
  nothing
end

main()

max_stars_repo_name: pseastham/eFEMpart

max_issues_repo_name: pseastham/eFEMpart


text:
#workspace()
#=
using EAGOParametricInterval
using IntervalArithmetic

opt1 = Param_Bisect_Opts() # sets bisection options to default
opt1.kmax_main = 1      # number of iterations for generalized bisection routine
opt1.kmax_cntr = 100       # number of iterations for contractor
opt1.style = "NewtonGS"    # specifies NewtonGS
opt1.DAGflag = false       # disable interval constraint prop

opt2 = Param_Bisect_Opts() # sets bisection options to default
opt2.kmax_main = 1      # number of iterations for generalized bisection routine
opt2.kmax_cntr = 100        # number of iterations for contractor
opt2.style = "KrawczykCW"    # specifies NewtonGS
opt2.DAGflag = false       # disable interval constraint prop

h1(z,p) = [z[1]^2+z[2]^2+p[1]*z[1]+4.0;
           z[1]+p[2]*z[2]]
h1j(z,p) = [2.0*z[1]+p[1]   2.0*z[2];
                      1.0   p[2]]
X1 = [Interval(-10.0,10.0),Interval(-2.0,2.0)]
P1 = [Interval(5.0,7.0),Interval(5.0,7.0)]

println("-------------------------------------------------")
println("----------- EXAMPLE 1 SUMMARY (NEWTON) ----------")
println("-------------------------------------------------")
println(" ")
Nsol1,Nin1,Niter1,Nhist1 = Generalized_Param_Bisection(X1,P1,h1,h1j,opt1)

println("-------------------------------------------------")
println("----------- EXAMPLE 1 SUMMARY (KRAWCZ) ----------")
println("-------------------------------------------------")
println(" ")
Ksol1,Kin1,Kiter1,Khist1 = Generalized_Param_Bisection(X1,P1,h1,h1j,opt2)

println("----------- EXAMPLE 2 ----------")
h2(z,p) = [(3.25-z[1])/p[1]-z[3];
           z[1]/p[2]-z[3];
           z[2]-(z[1]^2)/(1+z[1]^2)]
h2j(z,p) = [-1/p[1] 0 (-1);
            1/p[2] 0 (-1);
            -2*z[1]/(1+z[1]^2)^2 1 0]
X2 = [Interval(-30..30),Interval(-30..30),Interval(-30..30)]
P2 = [Interval(1800..2200),Interval(900..1100)]
# #=
println("-------------------------------------------------")
println("----------- EXAMPLE 2 SUMMARY (NEWTON) ----------")
println("-------------------------------------------------")
println(" ")
#Nsol2,Nin2,Niter2,Nhist2 = Generalized_Param_Bisection(X2,P2,h2,h2j,opt1)

println("-------------------------------------------------")
println("----------- EXAMPLE 2 SUMMARY (KRAWCZ) ----------")
println("-------------------------------------------------")
println(" ")
#Ksol2,Kin2,Kiter2,Khist2 = Generalized_Param_Bisection(X2,P2,h2,h2j,opt2)

h3(z,p) = [(exp(38.0*z[1]-1.0+log(10.0^(-9.0)))) + p[1]*z[1]-1.6722*z[2] + 0.6689*z[3]-8.0267
           (exp(38.0*z[2]-1.0)+log(1.98*10^(-9.0))) + 0.6622*z[1] + p[2]*z[2] + 0.6622*z[3] + 4.0535
           (exp(38.0*z[3]- 1.0+(10.0^(-9.0)))) + z[1] - z[2] + p[3]*z[3] - 6.0]
h3j(z,p) = [p[1]+exp(38.0*z[1]-1.0+log(38.0*(10.0^(-9.0)))) -1.6722  0.6689;
            0.6622 exp(38.0*z[2]-1.0+log(38*(1.98*10.0^(-9.0))))+p[2] 0.6622;
            1.0 -2.0 p[3]+exp(38.0*z[3]-1.0+log(38.0*(10.0^(-9.0))));]
X3 = [Interval(-30.0,30.0),Interval(-30.0,30.0),Interval(-30.0,30.0)]
P3 = [Interval(0.6020,0.7358),Interval(1.2110,1.4801),Interval(3.6,4.4)]
h3temp = h3(X3,P3)
hj3temp = h3j(X3,P3)
Nsol3,Nin3,Niter3,Nhist3 = Generalized_Param_Bisection(X3,P3,h3,h3j,opt1)
Ksol3,Kin3,Kiter3,Khist3 = Generalized_Param_Bisection(X3,P3,h3,h3j,opt2)

h4(z,p) = [(-z[1]^2+p[1])*z[1]]
h4j(z,p) = [-3*z[1]^2+p[1]]
X4 = [Interval(-10.0,10.0)]
P4 = [Interval(0.25..20)]
#Nsol4,Nin4,Niter4,Nhist4 = Generalized_Param_Bisection(X4,P4,h4,h4j,opt1)
#Ksol4,Kin4,Kiter4,Khist4 = Generalized_Param_Bisection(X4,P4,h4,h4j,opt2)

h5(z,p) = [p[1]*z[1]^5-25.2*z[1]^3+6*p[1]*z[1]-p[2]*z[2];
           2*p[2]*z[2]-p[1]*z[1]]
h5j(z,p) = [5*p[1]*z[1]^4-3*25.2*z[1]^2+6*p[1]  (-p[2]);
            -p[1]  2*p[2]]
X5 = [Interval(3.0,5.0),Interval(4.25,7.75)]
P5 = [Interval(-10.0,10.0),Interval(-10.0,10.0)]
#Nsol5,Nin5,Niter5,Nhist5 = Generalized_Param_Bisection(X5,P5,h5,h5j,opt1)
#Ksol5,Kin5,Kiter5,Khist5 = Generalized_Param_Bisection(X5,P5,h5,h5j,opt2)

println("---------- BEGIN SUMMARY OF TEST PROBLEMS --------------")

println("Test Problem #1: ")
println("Newton   Iterations Taken: $(Niter1) (Target = 169)")
println("Krawczyk Iterations Taken: $(Kiter1) (Target = 131)")
println("Newton   Boxes Taken: $(length(Nsol1)) (Target = 5)")
println("Krawczyk Boxes Taken: $(length(Ksol1)) (Target = 3)")

#=
println("Test Problem #2: ")
println("Newton   Iterations Taken: $(Niter2) (Target = 1)")
println("Krawczyk Iterations Taken: $(Kiter2) (Target = 1)")
println("Newton   Boxes Taken: $(length(Nsol2)) (Target = 1)")
println("Krawczyk Boxes Taken: $(length(Ksol2)) (Target = 1)")
=#
#=
println("Test Problem #3: ")
println("Newton   Iterations Taken: $(Niter3) (Target = 55)")
println("Krawczyk Iterations Taken: $(Kiter3) (Target = 329)")
println("Newton   Boxes Taken: $(length(Nsol3)) (Target = 1)")
println("Krawczyk Boxes Taken: $(length(Ksol3)) (Target = 1)")

println("Test Problem #4: ")
println("Newton   Iterations Taken: $(Niter4) (Target = 270)")
println("Krawczyk Iterations Taken: $(Kiter4) (Target = 295)")
println("Newton   Boxes Taken: $(length(Nsol4)) (Target = 15)")
println("Krawczyk Boxes Taken: $(length(Ksol4)) (Target = 15)")

println("Test Problem #5: ")
println("Newton   Iterations Taken: $(Niter5) (Target = 2392)")
println("Krawczyk Iterations Taken: $(Kiter5) (Target = 2599)")
println("Newton   Boxes Taken: $(length(Nsol5)) (Target = 80)")
println("Krawczyk Boxes Taken: $(length(Ksol5)) (Target = 80)")

println("---------- END SUMMARY OF TEST PROBLEMS --------------")

# generates directed graph contractor parameters
if (DAGflag)
  DAGr::Int64 = opt.DAGpass
  DAGpack::Vector{Interval{Float64}} = vcat(Xin,Pin)
  exprs::Vector{Expr} = vcat(opt.DAGh,opt.DAGg)
  hbnds::Vector{Float64} = zeros(Float64,length(Xin))
  if (length(opt.DAGgL)>0)
   gL::Vector{Float64} = vcat(hbnds,opt.DAGgL)
   gU::Vector{Float64} = vcat(hbnds,opt.DAGgU)
  else
   gL = vcat(hbnds)
   gU = vcat(hbnds)
  end
  npx::Int64 = length(Xin) + length(Pin)
  DAGparam = Generate_TapeList(exprs,npx,gL,gU)
end
=#

#=
if (~Eflag)
 if (DAGflag)
    #### Contracts box using interval contractors on DAG ####
    DAGpack = vcat(Xw,Pw)
    DAGContractor!(DAGpack,DAGparam,DAGr)
    for i=1:length(DAGpack)
      #### Fathoms by interval contractors ####
      if (isempty(DAGpack[i]))
       Eflag = true
       break
      end
    end
    for i=1:length(Xw)
      Xw[i] = DAGpack[i]
    end
    count = 1
    for i=(length(Xw)+1):(length(Xw)+length(Pw))
      Pw[count] = DAGpack[i]
      count += 1
    end
 end
end
=#


text:
# Part 1
grid = open(f -> readlines(f), "day03_input.txt")
function treecount(grid, right, down)
    pos = [1, 1]
    hits = 0
    while pos[2] <= length(grid)
        row = grid[pos[2]]
        hits += row[(pos[1] - 1) % length(row) + 1] == '#' ? 1 : 0
        pos += [right, down]
    end
    hits
end
println("There are $(treecount(grid, 3, 1)) trees encountered")

# Part 2
slopes = [[1, 1], [3, 1], [5, 1], [7, 1], [1, 2]]
println("Tree product: $(prod([treecount(grid, slope...) for slope in slopes]))")

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
using Test, Omniscape

d = run_omniscape("input/config2.ini")
a, b, c = run_omniscape("input/config.ini")

e = run_omniscape("input/config3.ini")

@test typeof(a) == Array{Float64,2}
@test typeof(b) == Array{Float64,2}
@test typeof(c) == Array{Float64,2}
@test typeof(d) == Array{Float64,2}
@test typeof(e) == Array{Float64,2}
@test b ≈ d

rm("test1_output", recursive = true)
rm("test2_output", recursive = true)
rm("test3_output", recursive = true)


text:
include("Batcher.jl")
include("ForwardBatcher.jl")

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
"""
    restrictedSequencingDD(model::SequencingModel, maxWidth::Int;rootNode::Union{RestrictedSequencingNode,RelaxedSequencingNode,RelaxedStateSequencingNode,Nothing}=nothing,bestKnownSolution::Union{Vector{Int},Nothing}=nothing,bestKnownValue::Union{T,Nothing}=nothing, rrbSetting::Union{RRBSetting,Nothing}=nothing, rrbPacket=nothing, rootDomain::BitVector=BitVector(),debugOn::Bool=false)where{T<:Real}

A quick restricted decision diagram for sequencing problems of fixed length
Returns: bestKnownSolution (Vector{Int}), bestKnownValue (Real), exact {Bool}
'exact' is true if the dd is exact

# Arguments
- `model::SequencingModel`: The sequencing problem
- `maxWidth::Int`: The max width allowed in a decision diagram layer
- `rootNode::Union{RestrictedSequencingNode,RelaxedSequencingNode,RelaxedStateSequencingNode,Nothing}`: The starting node for the DD
- `bestKnownSolution::Union{Vector{Int},Nothing}`: An optional parameter that gives the DD a starting solution
- `bestKnownValue::Union{T,Nothing}`: An optional parameter that gives the DD a starting bound
- `rrbSetting::Union{RRBSetting,Nothing}`: An optional parameter that chooses an rough relaxed bounding method
- `rrbPacket`: Data to be passed to the rrbFunction
- `rootDomain::BitVector`: The domain of the root node if it is a RelaxedStateSequencingNode
- `debugOn::Bool`: An optional parameter that can be used to turn on debug statements in other functions
"""
function restrictedSequencingDD(model::SequencingModel, maxWidth::Int;rootNode::Union{RestrictedSequencingNode,RelaxedSequencingNode,RelaxedStateSequencingNode,Nothing}=nothing,bestKnownSolution::Union{Vector{Int},Nothing}=nothing,bestKnownValue::Union{T,Nothing}=nothing, rrbSetting::Union{RRBSetting,Nothing}=nothing, rrbPacket=nothing, rootDomain::BitVector=BitVector(),debugOn::Bool=false)where{T<:Real}
    doSecondLayerDomainCheck = false
    if isnothing(rootNode)
        rootNode = createRestrictedRootNode(model)
        queue = Vector{RestrictedSequencingNode}([rootNode])
        pathLength = length(rootNode.path)
    elseif isa(rootNode,RelaxedSequencingNode)
        pathLength = length(rootNode)
        rootDomain = falses(length(model))
        for arc in rootNode.arcsOut
            rootDomain[arc.assignment] = true
        end
        queue = Vector{RestrictedSequencingNode}([RestrictedSequencingNode(copy(rootNode.bestSeq),rootNode.allDown, rootNode.lengthToRoot,rootDomain)])
        if model.precedenceConstraints
            doSecondLayerDomainCheck = true
        end
    elseif isa(rootNode,RelaxedStateSequencingNode)
        pathLength = length(rootNode)
        queue = Vector{RestrictedSequencingNode}([RestrictedSequencingNode(copy(rootNode.bestSeq),rootNode.allDown, rootNode.lengthToRoot,rootDomain)])
        if model.precedenceConstraints
            doSecondLayerDomainCheck = true
        end
    else
        queue = Vector{RestrictedSequencingNode}([rootNode])
        pathLength = length(rootNode.path)
    end

    #track exactness
    exact = true

    #iterate down through the layers
    for layerIndex in 1:length(model)-pathLength
        if layerIndex == 2 && doSecondLayerDomainCheck
            secondLayerDomainCheck(model, queue)
        end
        #initialize array for new layer
        newLayer = Vector{RestrictedSequencingNode}()
        #for each node in the queue
        for node in queue
            if length(node.path)>0
                state = getState(node)
            else
                state = nothing
            end
            #for each unvisited location of the node
            for i in 1:length(model)
                #if i is in the domain add the new node
                if inDomain(node, i)
                    newNode = makeDecision(model, node, state, i)

                    if !isnothing(rrbSetting) && layerIndex+pathLength<length(model)-2
                        if  !isBetterSolutionValue(model, sequencingRRB(rrbSetting, newNode, rrbPacket), bestKnownValue)
                            continue
                        end
                    end

                    #if worth adding, add it
                    if isnothing(bestKnownValue) || model.objective.type==maximization
                        addToSortedLayer!(model, newLayer, newNode)
                    elseif isBetterSolutionValue(model, newNode.value, bestKnownValue)
                        addToSortedLayer!(model, newLayer, newNode)
                    end
                end
            end#end domain loop
            if length(newLayer)>maxWidth
                if exact
                    exact = false
                end
                newLayer = newLayer[1:maxWidth]
            end
        end#node in queue loop
        #trim queue
        if length(newLayer)>maxWidth
            if exact
                exact = false
            end
            queue = newLayer[1:maxWidth]
        elseif isempty(newLayer)
            return bestKnownSolution, bestKnownValue, exact
        else
            queue = newLayer
        end
    end#iterate over layers loop

    #update result
    if isnothing(bestKnownValue) || isBetterSolutionValue(model, queue[1].value, bestKnownValue)
        bestKnownSolution = queue[1].path
        bestKnownValue = queue[1].value
    end
    return bestKnownSolution, bestKnownValue, exact
end

"""
    makeDecision(model::SequencingModel, rootNode::RestrictedSequencingNode, state::Union{Nothing, Int}, decision::Int)

A method that generates the child node resulting from adding an arc to an existing node

# Arguments
- `model::SequencingModel`: The Sequencing problem
- `rootNode::RestrictedSequencingNode`: The starting node for the decision
- `state::Union{Nothing, Int}`: The last node visited
- `decision::Int`: The node being visited
"""
function makeDecision(model::SequencingModel, rootNode::RestrictedSequencingNode, state::Union{Nothing, Int}, decision::Int)
    #update path
    newPath = copy(rootNode.path)
    push!(newPath, decision)

    #update visited
    newVisited = copy(rootNode.visited)
    newVisited[decision] = true

    #update domain
    newDomain = copy(rootNode.domain)
    newDomain[decision] = false
    #handle precedence constraints
    if model.precedenceConstraints
        if haskey(model.precedenceP2F,decision)
            for i in model.precedenceP2F[decision]
                #check each of the relevant precedence constraints
                #if satisfied, put in the domain
                if haskey(model.precedenceF2P,i)
                    if checkPrecedence(newVisited,model.precedenceF2P[i])
                        newDomain[i] = true
                    end
                end
            end
        end
    end

    #update value
    if isnothing(state)
        return RestrictedSequencingNode(newPath, newVisited, rootNode.value, newDomain)
    else
        return RestrictedSequencingNode(newPath, newVisited, rootNode.value + evaluateDecision(model.objective, state, decision), newDomain)
    end
end

"""
    addToSortedLayer!(model::SequencingModel, layer::Vector{RestrictedSequencingNode}, newNode::RestrictedSequencingNode

A method that generates the child node resulting from adding an arc to an existing node

# Arguments
- `model::SequencingModel`: The Sequencing problem
- `layer::Vector{RestrictedSequencingNode}`: The list to add to
- `newNode::RestrictedSequencingNode`: the node to insert
"""
function addToSortedLayer!(model::SequencingModel, layer::Vector{RestrictedSequencingNode}, newNode::RestrictedSequencingNode)
    if isempty(layer)
        push!(layer, newNode)
    else
        insert!(
            layer,
            searchsorted(layer, newNode, by = x -> valueSequencingNode(model, x)).start,
            newNode
        )
    end
end

"""
    valueSequencingNode(model::CVRPModel, node::RestrictedCVRPNode)

Return a quick heuristic value of a node, lower is better

# Arguments
- `model::SequencingModel`: The sequencing problem
- `newNode::RestrictedSequencingNode`: the node to value
"""
function valueSequencingNode(model::SequencingModel, node::RestrictedSequencingNode)
    if model.objective.type == minimization
        return node.value
    else
        return (-1 * node.value)
    end
end

"""
    secondLayerDomainCheck(model::SequencingModel, queue::Vector{RestrictedSequencingNode})

A method that fixes the domain of the first nodes generated from a relaxed node

# Arguments
- `model::SequencingModel`: The Sequencing problem
- `queue::Vector{RestrictedSequencingNode}`: The list of nodes to check
"""
function secondLayerDomainCheck(model::SequencingModel, queue::Vector{RestrictedSequencingNode})
    for node in queue
        for i in 1:length(node.domain)
            if !node.visited[i] && !node.domain[i] && hasKey(model.objective, getState(node), i)
                #handle precedence constraints
                if haskey(model.precedenceF2P,i)
                    if checkPrecedence(node.visited,model.precedenceF2P[i])
                        node.domain[i] = true
                    end
                else
                    node.domain[i] = true
                end
            end
        end
    end
end


text:
include("euler/euler.jl")
using .BigIntegers: BigInteger
using .Numbers: digits_sum
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(n::Int)::Int
    current, previous = BigInteger(1), BigInteger(1)
    for i ∈ 2:n + 1
        current, previous = (i % 3 == 1 ? 2i ÷ 3 : 1) * current + previous, current
    end
    return digits_sum(current)
end

compute(10)

compute(100)

@benchmark compute(100)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
using Reactive

## Basics

facts("Basic checks") do
    x = Signal(Float32)
    @fact isa(x, Signal{Type{Float32}}) --> true

    a = Signal(number(); name="a")
    b = map(x -> x*x, a; name="b")

    context("map") do

        # Lift type
        #@fact typeof(b) --> Reactive.Lift{Int}

        # type conversion
        push!(a, 1.0)
        step()
        @fact value(b) --> 1
        # InexactError to be precise
        push!(a, 2.1, (n,x,error_node,err) -> @fact n --> a)
        step()

        @fact value(b) --> 1

        push!(a, number())
        step()
        @fact value(b) --> value(a)^2

        push!(a, -number())
        step()
        @fact value(b) --> value(a)^2

        ## Multiple inputs to Lift
        d = Signal(number())
        c = map(+, a, b, d, typ=Int)
        @fact value(c) --> value(a) + value(b) + value(d)

        push!(a, number())
        step()
        @fact value(c) --> value(a) + value(b) + value(d)

        push!(d, number())
        step()
        @fact value(c) --> value(a) + value(b) + value(d)
    end

    context("merge") do
        ## Merge
        d = Signal(number(); name="d")
        e = merge(b, d, a; name="e")

        # precedence to d
        @fact value(e) --> value(d)

        push!(a, number())
        step()
        # precedence to b over a -- a is older.

        @fact value(e) --> value(b)

        c = map(identity, a) # Make a younger than b
        f = merge(d, c, b)
        push!(a, number())
        step()
        @fact value(f) --> value(c)
    end

    context("foldp") do

        ## foldl over time

        gc()
        push!(a, 0)
        step()
        f = foldp(+, 0, a)
        nums = round.([Int], rand(100)*1000)
        map(x -> begin push!(a, x); step() end, nums)

        @fact sum(nums) --> value(f)
    end

    context("filter") do
        # filter
        g = Signal(0)
        pred = x -> x % 2 != 0
        h = filter(pred, 1, g)
        j = filter(x -> x % 2 == 0, 1, g)

        @fact value(h) --> 1
        @fact value(j) --> 0

        push!(g, 2)
        step()
        @fact value(h) --> 1

        push!(g, 3)
        step()
        @fact value(h) --> 3

        g = Signal(0)
        pred = x -> x % 2 != 0
        h = filter(pred, g)

        push!(g, 2)
        step()
        @fact value(h) --> 0

        push!(g, 3)
        step()
        @fact value(h) --> 3
    end

    context("filter counts") do
        a = Signal(1; name="a")
        b = Signal(2; name="b")
        c = filter(value(a), a; name="c") do aval; aval > 1 end
        d = map(*,b,c)
        count = foldp((x, y) -> x+1, 0, d)
        @fact value(count) --> 0
        push!(a, 0)
        step()
        @fact value(count) --> 0
    end

    context("sampleon") do
        # sampleon
        g = Signal(0)
        nv = number()
        push!(g, nv)
        println("step 1")
        step()
        i = Signal(true)
        j = sampleon(i, g)
        # default value
        @fact value(j) --> value(g) # j == g == nv
        push!(g, value(g)-1)
        println("step 2")
        step()
        @fact value(j) --> value(g)+1 # g is nv - 1, j is unchanged on nv
        push!(i, true)
        println("step 3")
        step() # resample
        @fact value(j) --> value(g)
    end

    context("droprepeats") do
        # droprepeats
        count = s -> foldp((x, y) -> x+1, 0, s)

        k = Signal(1)
        l = droprepeats(k)

        @fact value(l) --> value(k)
        push!(k, 1)
        step()
        @fact value(l) --> value(k)
        push!(k, 0)
        step()
        #println(l.value, " ", value(k))
        @fact value(l) --> value(k)

        m = count(k)
        n = count(l)

        seq = [1, 1, 1, 0, 1, 0, 1, 0, 0]
        map(x -> begin push!(k, x); step() end, seq)

        @fact value(m) --> length(seq)
        @fact value(n) --> 6
    end

    context("filterwhen") do
        # filterwhen
        b = Signal(false)
        n = Signal(1)
        dw = filterwhen(b, 0, n)
        @fact value(dw) --> 0
        push!(n, 2)
        step()
        @fact value(dw) --> 0
        push!(b, true)
        step()
        @fact value(dw) --> 0
        push!(n, 1)
        step()
        @fact value(dw) --> 1
        push!(n, 2)
        step()
        @fact value(dw) --> 2
        dw = filterwhen(b, 0, n)
        @fact value(dw) --> 2
    end

    context("push! inside push!") do
        a = Signal(0)
        b = Signal(1)
        Reactive.preserve(map(x -> push!(a, x), b))

        @fact value(a) --> 0

        step()
        @fact value(a) --> 1

        push!(a, 2)
        step()
        @fact value(a) --> 2
        @fact value(b) --> 1

        push!(b, 3)
        step()
        @fact value(b) --> 3
        @fact value(a) --> 2

        step()
        @fact value(a) --> 3
    end

    context("previous") do
        x = Signal(0)
        y = previous(x)
        @fact value(y) --> 0

        push!(x, 1)
        step()

        @fact value(y) --> 0

        push!(x, 2)
        step()

        @fact value(y) --> 1

        push!(x, 3)
        step()

        @fact value(y) --> 2
        @fact queue_size() --> 0
    end

    context("delay") do
        x = Signal(0)
        y = delay(x)
        @fact value(y) --> 0

        push!(x, 1)
        step()

        @fact queue_size() --> 1
        @fact value(y) --> 0

        step()
        @fact value(y) --> 1
        @fact queue_size() --> 0
    end

    context("bind") do
        x = Signal(0; name="x")
        y = Signal(0; name="y")
        zpre_count = 0
        zpost_count = 0
        zpre = map(yval->(zpre_count+=1; 2yval), y; name="zpre")
        # map(...) runs the function once to get the init value on creation
        @fact zpre_count --> 1
        bind!(y, x)
        @fact zpre_count --> 2 # initialising the bind should cause zpre to run too
        zpost = map(yval->(zpost_count+=1; 2yval), y; name="zpost")

        @fact zpre_count --> 2
        @fact zpost_count --> 1

        @show queue_size()
        push!(x,1000)
        step()

        @fact value(y) --> 1000
        @fact value(zpre) --> 2000
        @fact value(zpost) --> 2000
        @fact zpre_count --> 3
        @fact zpost_count --> 2
        @fact bound_srcs(y) --> [x]
        @fact bound_dests(x) --> [y]

        unbind!(y,x)
        push!(x,0)
        step()

        @fact value(y) --> 1000
        @fact value(zpre) --> 2000
        @fact value(zpost) --> 2000

        # bind where dest is before src in node list
        a = Signal(1; name="a")
        b = map(x->2x, a; name="b")
        c = Signal(1; name="c")
        d = map(x->4x, c; name="d")
        bind!(a, d)
        @fact value(d) --> value(a)

        @fact queue_size() --> 0

        push!(c, 3)
        @fact queue_size() --> 1
        step()
        @fact value(c) --> 3
        @fact value(d) --> 4*3
        @fact value(a) --> 4*3
        @fact value(b) --> 2*4*3

    end

    context("bindmap") do

        src2dst(x) = x + 2
        dst2src(x) = x - 2

        # oneway with initiation
        src = Signal(3)
        dst = Signal(0)
        bindmap!(dst, src2dst, src)
        @fact value(dst) --> 5
        push!(src, 2)
        step()
        @fact value(dst) --> 4
        push!(dst, 5)
        step()
        @fact value(src) --> 2
        unbind!(dst,src)
        push!(src,1)
        step()
        @fact value(dst) --> 5

        # twoway with initiation
        src = Signal(3)
        dst = Signal(0)
        bindmap!(dst, src2dst, src, dst2src)
        @fact value(dst) --> 5
        push!(src, 2)
        step()
        @fact value(dst) --> 4
        push!(dst, 5)
        step()
        @fact value(src) --> 3
        unbind!(dst,src,false) # test oneway first
        push!(src,1)
        step()
        @fact value(dst) --> 5
        push!(dst,1)
        step()
        @fact value(src) --> -1

        # This should work but I think line 444 in operators.jl shouldn't exit if the twoway unbind! was applied after a oneway unbind!.
        # unbind!(dst,src) # test other way as well
        # push!(dst,2)
        # step()
        # @fact value(src) --> -1

        # oneway without initiation
        src = Signal(3)
        dst = Signal(0)
        bindmap!(dst, src2dst, src, initial=false)
        @fact value(dst) --> 0
        push!(src, 2)
        step()
        @fact value(dst) --> 4
        push!(dst, 5)
        step()
        @fact value(src) --> 2
        unbind!(dst,src)
        push!(src,1)
        step()
        @fact value(dst) --> 5

        # twoway without initiation
        src = Signal(3)
        dst = Signal(0)
        bindmap!(dst, src2dst, src, dst2src, initial=false)
        @fact value(dst) --> 0
        push!(src, 2)
        step()
        @fact value(dst) --> 4
        push!(dst, 5)
        step()
        @fact value(src) --> 3
        unbind!(dst,src)
        push!(src,1)
        step()
        @fact value(dst) --> 5
        push!(dst,2)
        step()
        @fact value(src) --> 1

    end

end


text:
using TropicalTensors
using ForwardDiff
using Test

@testset "spinglass" begin
    lt = Cylinder(10, 8)
    sg = Spinglass(lt, ones(Float32, 152), zeros(Float32, 80))
    res = solve(sg; usecuda=false)
    @test res.n == 152
end

@testset "counting tropical" begin
    lt = Cylinder(10, 8)
    sg = Spinglass(lt, ones(Float32, 152), zeros(Float32, 80))
    res = solve(CountingTropicalF32, sg; usecuda=false)
    @test res.n == 152
    @test res.c == 2
end

using BitBasis
function count_degeneracy_exact(sg::Spinglass{LT,T}) where {LT, T}
    elist = T[]
    bonds = sgbonds(sg.lattice)
    for b in basis(BitStr64{length(sg.lattice)})
        eng = zero(T)
        for ((i, j), J) in zip(bonds, sg.Js)
            eng += b[i] == b[j] ? J : -J
        end
        push!(elist, eng)
    end
    findall(==(minimum(elist)), elist) |> length
end

@testset "counting tropical 2" begin
    lt = Cylinder(4, 4)
    sg = Spinglass(lt, rand([-1, 1], 28), zeros(Int, 16))
    res = solve(CountingTropical{Int,Int}, sg; usecuda=false)
    @test res.c == count_degeneracy_exact(sg)
end

@testset "forwarddiff" begin
    L = 3
    Js = ones(Float32, 2*L^2-L)
    gs = ForwardDiff.gradient(x->solve(Cylinder(L, L), x, zeros(eltype(x), L^2); usecuda=false).n, Js)
    @test gs ≈ ones(Float32, 2*L^2-L)
end


text:
struct Camera3D <: AbstractCamera
    eyeposition::Observable{Vec3f}
    lookat::Observable{Vec3f}
    upvector::Observable{Vec3f}

    zoom_mult::Observable{Float32}
    fov::Observable{Float32} # WGLMakie compat
    near::Observable{Float32}
    far::Observable{Float32}
    pulser::Observable{Float64}

    attributes::Attributes
end

"""
    Camera3D(scene[; attributes...])

Creates a 3d camera with a lot of controls.

The 3D camera is (or can be) unrestricted in terms of rotations and translations. Both `cam3d!(scene)` and `cam3d_cad!(scene)` create this camera type. Unlike the 2D camera, settings and controls are stored in the `cam.attributes` field rather than in the struct directly, but can still be passed as keyword arguments. The general camera settings include

- `fov = 45f0` sets the "neutral" field of view, i.e. the fov corresponding to no zoom. This is irrelevant if the camera uses an orthographic projection.
- `near = automatic` sets the value of the near clip. By default this will be chosen based on the scenes bounding box. The final value is in `cam.near`.
- `far = automatic` sets the value of the far clip. By default this will be chosen based on the scenes bounding box. The final value is in `cam.far`.
- `rotation_center = :lookat` sets the default center for camera rotations. Currently allows `:lookat` or `:eyeposition`.
- `projectiontype = Perspective` sets the type of the projection. Can be `Orthographic` or `Perspective`.
- `fixed_axis = false`: If true panning uses the (world/plot) z-axis instead of the camera up direction.
- `zoom_shift_lookat = true`: If true attempts to keep data under the cursor in view when zooming.
- `cad = false`: If true rotates the view around `lookat` when zooming off-center.

The camera view follows from the position of the camera `eyeposition`, the point which the camera focuses `lookat` and the up direction of the camera `upvector`. These can be accessed as `cam.eyeposition` etc and adjusted via `update_cam!(scene, cameracontrols(scene), eyeposition, lookat[, upvector = Vec3f(0, 0, 1)])`. They can also be passed as keyword arguments when the camera is constructed.

The camera can be controlled by keyboard and mouse. The keyboard has the following available attributes

- `up_key   = Keyboard.r` sets the key for translations towards the top of the screen.
- `down_key = Keyboard.f` sets the key for translations towards the bottom of the screen.
- `left_key  = Keyboard.a` sets the key for translations towards the left of the screen.
- `right_key = Keyboard.d` sets the key for translations towards the right of the screen.
- `forward_key  = Keyboard.w` sets the key for translations into the screen.
- `backward_key = Keyboard.s` sets the key for translations out of the screen.

- `zoom_in_key   = Keyboard.u` sets the key for zooming into the scene (enlarge, via fov).
- `zoom_out_key  = Keyboard.o` sets the key for zooming out of the scene (shrink, via fov).
- `stretch_view_key  = Keyboard.page_up` sets the key for moving `eyepostion` away from `lookat`.
- `contract_view_key = Keyboard.page_down` sets the key for moving `eyeposition` towards `lookat`.

- `pan_left_key  = Keyboard.j` sets the key for rotations around the screens vertical axis.
- `pan_right_key = Keyboard.l` sets the key for rotations around the screens vertical axis.
- `tilt_up_key   = Keyboard.i` sets the key for rotations around the screens horizontal axis.
- `tilt_down_key = Keyboard.k` sets the key for rotations around the screens horizontal axis.
- `roll_clockwise_key        = Keyboard.e` sets the key for rotations of the screen.
- `roll_counterclockwise_key = Keyboard.q` sets the key for rotations of the screen.

- `keyboard_rotationspeed = 1f0` sets the speed of keyboard based rotations.
- `keyboard_translationspeed = 0.5f0` sets the speed of keyboard based translations.
- `keyboard_zoomspeed = 1f0` sets the speed of keyboard based zooms.
- `update_rate = 1/30` sets the rate at which keyboard based camera updates are evaluated.

and mouse interactions are controlled by

- `translation_button = Mouse.right` sets the mouse button for drag-translations. (up/down/left/right)
- `scroll_mod = true` sets an additional modifier button for scroll-based zoom. (true being neutral)
- `rotation_button = Mouse.left` sets the mouse button for drag-rotations. (pan, tilt)

- `mouse_rotationspeed = 1f0` sets the speed of mouse rotations.
- `mouse_translationspeed = 0.5f0` sets the speed of mouse translations.
- `mouse_zoomspeed = 1f0` sets the speed of mouse zooming (mousewheel).
- `circular_rotation = (true, true, true)` enables circular rotations for (fixed x, fixed y, fixed z) rotation axis. (This means drawing a circle with your mouse around the center of the scene will result in a continuous rotation.)

There are also a few generally applicable controls:

- `fix_x_key = Keyboard.x` sets the key for fixing translations and rotations to the (world/plot) x-axis.
- `fix_y_key = Keyboard.y` sets the key for fixing translations and rotations to the (world/plot) y-axis.
- `fix_z_key = Keyboard.z` sets the key for fixing translations and rotations to the (world/plot) z-axis.
- `reset = Keyboard.home` sets the key for fully resetting the camera. This equivalent to setting `lookat = Vec3f(0)`, `upvector = Vec3f(0, 0, 1)`, `eyeposition = Vec3f(3)` and then calling `center!(scene)`.

You can also make adjustments to the camera position, rotation and zoom by calling relevant functions:

- `translate_cam!(scene, v)` will translate the camera by the given world/plot space vector `v`.
- `rotate_cam!(scene, angles)` will rotate the camera around its axes with the corresponding angles. The first angle will rotate around the cameras "right" that is the screens horizontal axis, the second around the up vector/vertical axis or `Vec3f(0, 0, +-1)` if `fixed_axis = true`, and the third will rotate around the view direction i.e. the axis out of the screen. The rotation respects the the current `rotation_center` of the camera.
- `zoom!(scene, zoom_step)` will change the zoom level of the scene without translating or rotating the scene. `zoom_step` applies multiplicatively to `cam.zoom_mult` which is used as a multiplier to the fov (perspective projection) or width and height (orthographic projection).
"""
function Camera3D(scene; kwargs...)
    attr = merged_get!(:cam3d, scene, Attributes(kwargs)) do
        Attributes(
            # Keyboard controls
            # Translations
            up_key        = Keyboard.r,
            down_key      = Keyboard.f,
            left_key      = Keyboard.a,
            right_key     = Keyboard.d,
            forward_key   = Keyboard.w,
            backward_key  = Keyboard.s,
            # Zooms
            zoom_in_key   = Keyboard.u,
            zoom_out_key  = Keyboard.o,
            stretch_view_key  = Keyboard.page_up,
            contract_view_key = Keyboard.page_down,
            # Rotations
            pan_left_key  = Keyboard.j,
            pan_right_key = Keyboard.l,
            tilt_up_key   = Keyboard.i,
            tilt_down_key = Keyboard.k,
            roll_clockwise_key        = Keyboard.e,
            roll_counterclockwise_key = Keyboard.q,
            # Mouse controls
            translation_button = Mouse.right,
            scroll_mod         = true,
            rotation_button    = Mouse.left,
            # Shared controls
            fix_x_key = Keyboard.x,
            fix_y_key = Keyboard.y,
            fix_z_key = Keyboard.z,
            reset = Keyboard.home,
            # Settings
            keyboard_rotationspeed = 1f0,
            keyboard_translationspeed = 0.5f0,
            keyboard_zoomspeed = 1f0,
            mouse_rotationspeed = 1f0,
            mouse_translationspeed = 1f0,
            mouse_zoomspeed = 1f0,
            circular_rotation = (true, true, true),
            fov = 45f0, # base fov
            near = automatic,
            far = automatic,
            rotation_center = :lookat,
            update_rate = 1/30,
            projectiontype = Perspective,
            fixed_axis = true,
            zoom_shift_lookat = false, # doesn't really work with fov
            cad = false,
            # internal
            selected = true
        )
    end

    cam = Camera3D(
        pop!(attr, :eyeposition, Vec3f(3)),
        pop!(attr, :lookat,      Vec3f(0)),
        pop!(attr, :upvector,    Vec3f(0, 0, 1)),

        Observable(1f0),
        Observable(attr[:fov][]),
        Observable(attr[:near][] === automatic ? 0.1f0 : attr[:near][]),
        Observable(attr[:far][]  === automatic ? 100f0 : attr[:far][]),
        Observable(-1.0),

        attr
    )

    disconnect!(camera(scene))

    # Keyboard controls
    # ticks every so often to get consistent position updates.
    on(cam.pulser) do prev_time
        current_time = time()
        active = on_pulse(scene, cam, Float32(current_time - prev_time))
        @async if active && attr.selected[]
            sleep(attr.update_rate[])
            cam.pulser[] = current_time
        else
            cam.pulser.val = -1.0
        end
    end

    keynames = (
        :up_key, :down_key, :left_key, :right_key, :forward_key, :backward_key,
        :zoom_in_key, :zoom_out_key, :stretch_view_key, :contract_view_key,
        :pan_left_key, :pan_right_key, :tilt_up_key, :tilt_down_key,
        :roll_clockwise_key, :roll_counterclockwise_key
    )

    # Start ticking if relevant keys are pressed
    on(camera(scene), events(scene).keyboardbutton) do event
        if event.action in (Keyboard.press, Keyboard.repeat) && cam.pulser[] == -1.0 &&
            attr.selected[] && any(key -> ispressed(scene, attr[key][]), keynames)

            cam.pulser[] = time()
            return Consume(true)
        end
        return Consume(false)
    end

    # de/select plot on click outside/inside
    # also deselect other cameras
    deselect_all_cameras!(root(scene))
    on(camera(scene), events(scene).mousebutton, priority = 100) do event
        if event.action == Mouse.press
            attr.selected[] = is_mouseinside(scene)
        end
        return Consume(false)
    end

    # Mouse controls
    add_translation!(scene, cam)
    add_rotation!(scene, cam)

    # add camera controls to scene
    cameracontrols!(scene, cam)

    # Trigger updates on scene resize and settings change
    on(camera(scene), scene.px_area, attr[:fov], attr[:projectiontype]) do _, _, _
        update_cam!(scene, cam)
    end
    on(camera(scene), attr[:near], attr[:far]) do near, far
        near === automatic || (cam.near[] = near)
        far  === automatic || (cam.far[] = far)
        update_cam!(scene, cam)
    end

    # reset
    on(camera(scene), events(scene).keyboardbutton) do event
        if attr.selected[] && event.key == attr[:reset][] && event.action == Keyboard.release
            # center keeps the rotation of the camera so we reset that here
            # might make sense to keep user set lookat, upvector, eyeposition
            # around somewhere for this?
            cam.lookat[] = Vec3f(0)
            cam.upvector[] = Vec3f(0,0,1)
            cam.eyeposition[] = Vec3f(3)
            center!(scene)
            return Consume(true)
        end
        return Consume(false)
    end

    cam
end

# These imitate the old camera
function cam3d!(scene; zoom_shift_lookat = true, fixed_axis = true, kwargs...)
    Camera3D(scene, zoom_shift_lookat = zoom_shift_lookat, fixed_axis = fixed_axis; kwargs...)
end
function cam3d_cad!(scene; cad = true, zoom_shift_lookat = false, fixed_axis = false, kwargs...)
    Camera3D(scene, cad = cad, zoom_shift_lookat = zoom_shift_lookat, fixed_axis = fixed_axis; kwargs...)
end

function deselect_all_cameras!(scene)
    cam = cameracontrols(scene)
    cam isa Camera3D && (cam.attributes.selected[] = false)
    for child in scene.children
        deselect_all_cameras!(child)
    end
    nothing
end

function add_translation!(scene, cam::Camera3D)
    translationspeed = cam.attributes[:mouse_translationspeed]
    zoomspeed = cam.attributes[:mouse_zoomspeed]
    shift_lookat = cam.attributes[:zoom_shift_lookat]
    cad = cam.attributes[:cad]
    button = cam.attributes[:translation_button]
    scroll_mod = cam.attributes[:scroll_mod]

    last_mousepos = RefValue(Vec2f(0, 0))
    dragging = RefValue(false)

    function compute_diff(delta)
        if cam.attributes[:projectiontype][] == Orthographic
            aspect = Float32((/)(widths(scene.px_area[])...))
            aspect_scale = Vec2f(1f0 + aspect, 1f0 + 1f0 / aspect)
            return cam.zoom_mult[] * delta .* aspect_scale ./ widths(scene.px_area[])
        else
            viewdir = cam.lookat[] - cam.eyeposition[]
            return 0.002f0 * cam.zoom_mult[] * norm(viewdir) * delta
        end
    end

    # drag start/stop
    on(camera(scene), scene.events.mousebutton) do event
        if ispressed(scene, button[])
            if event.action == Mouse.press && is_mouseinside(scene) && !dragging[]
                last_mousepos[] = mouseposition_px(scene)
                dragging[] = true
                return Consume(true)
            end
        elseif event.action == Mouse.release && dragging[]
            mousepos = mouseposition_px(scene)
            diff = compute_diff(last_mousepos[] - mousepos)
            last_mousepos[] = mousepos
            dragging[] = false
            translate_cam!(scene, cam, translationspeed[] * Vec3f(diff[1], diff[2], 0f0))
            return Consume(true)
        end
        return Consume(false)
    end

    # in drag
    on(camera(scene), scene.events.mouseposition) do mp
        if dragging[] && ispressed(scene, button[])
            mousepos = screen_relative(scene, mp)
            diff = compute_diff(last_mousepos[] - mousepos)
            last_mousepos[] = mousepos
            translate_cam!(scene, cam, translationspeed[] * Vec3f(diff[1], diff[2], 0f0))
            return Consume(true)
        end
        return Consume(false)
    end

    on(camera(scene), scene.events.scroll) do scroll
        if is_mouseinside(scene) && ispressed(scene, scroll_mod[])
            zoom_step = (1f0 + 0.1f0 * zoomspeed[]) ^ -scroll[2]
            zoom!(scene, cam, zoom_step, shift_lookat[], cad[])
            return Consume(true)
        end
        return Consume(false)
    end
end

function add_rotation!(scene, cam::Camera3D)
    rotationspeed = cam.attributes[:mouse_rotationspeed]
    button = cam.attributes[:rotation_button]
    last_mousepos = RefValue(Vec2f(0, 0))
    dragging = RefValue(false)
    e = events(scene)

    # drag start/stop
    on(camera(scene), e.mousebutton) do event
        if ispressed(scene, button[])
            if event.action == Mouse.press && is_mouseinside(scene) && !dragging[]
                last_mousepos[] = mouseposition_px(scene)
                dragging[] = true
                return Consume(true)
            end
        elseif event.action == Mouse.release && dragging[]
            mousepos = mouseposition_px(scene)
            dragging[] = false
            rot_scaling = rotationspeed[] * (e.window_dpi[] * 0.005)
            mp = (last_mousepos[] - mousepos) * 0.01f0 * rot_scaling
            last_mousepos[] = mousepos
            rotate_cam!(scene, cam, Vec3f(-mp[2], mp[1], 0f0), true)
            return Consume(true)
        end
        return Consume(false)
    end

    # in drag
    on(camera(scene), e.mouseposition) do mp
        if dragging[] && ispressed(scene, button[])
            mousepos = screen_relative(scene, mp)
            rot_scaling = rotationspeed[] * (e.window_dpi[] * 0.005)
            mp = (last_mousepos[] .- mousepos) * 0.01f0 * rot_scaling
            last_mousepos[] = mousepos
            rotate_cam!(scene, cam, Vec3f(-mp[2], mp[1], 0f0), true)
            return Consume(true)
        end
        return Consume(false)
    end
end

function on_pulse(scene, cam, timestep)
    attr = cam.attributes

    # translation
    right = ispressed(scene, attr[:right_key][])
    left = ispressed(scene, attr[:left_key][])
    up = ispressed(scene, attr[:up_key][])
    down = ispressed(scene, attr[:down_key][])
    backward = ispressed(scene, attr[:backward_key][])
    forward = ispressed(scene, attr[:forward_key][])
    translating = right || left || up || down || backward || forward

    if translating
        # translation in camera space x/y/z direction
        translation = attr[:keyboard_translationspeed][] * timestep *
            Vec3f(right - left, up - down, backward - forward)
        viewdir = cam.lookat[] - cam.eyeposition[]
        _translate_cam!(scene, cam, cam.zoom_mult[] * norm(viewdir) * translation)
    end

    # rotation
    up = ispressed(scene, attr[:tilt_up_key][])
    down = ispressed(scene, attr[:tilt_down_key][])
    left = ispressed(scene, attr[:pan_left_key][])
    right = ispressed(scene, attr[:pan_right_key][])
    counterclockwise = ispressed(scene, attr[:roll_counterclockwise_key][])
    clockwise = ispressed(scene, attr[:roll_clockwise_key][])
    rotating = up || down || left || right || counterclockwise || clockwise

    if rotating
        # rotations around camera space x/y/z axes
        angles = attr[:keyboard_rotationspeed][] * timestep *
            Vec3f(up - down, left - right, counterclockwise - clockwise)

        _rotate_cam!(scene, cam, angles)
    end

    # zoom
    zoom_out = ispressed(scene, attr[:zoom_out_key][])
    zoom_in = ispressed(scene, attr[:zoom_in_key][])
    zooming = zoom_out || zoom_in

    if zooming
        zoom_step = (1f0 + attr[:keyboard_zoomspeed][] * timestep) ^ (zoom_out - zoom_in)
        _zoom!(scene, cam, zoom_step, false)
    end

    stretch = ispressed(scene, attr[:stretch_view_key][])
    contract = ispressed(scene, attr[:contract_view_key][])
    if stretch || contract
        zoom_step = (1f0 + attr[:keyboard_zoomspeed][] * timestep) ^ (stretch - contract)
        cam.eyeposition[] = cam.lookat[] + zoom_step * (cam.eyeposition[] - cam.lookat[])
    end
    zooming = zooming || stretch || contract

    # if any are active, update matrices, else stop clock
    if translating || rotating || zooming
        update_cam!(scene, cam)
        return true
    else
        return false
    end
end

function translate_cam!(scene::Scene, cam::Camera3D, t::VecTypes)
    _translate_cam!(scene, cam, t)
    update_cam!(scene, cam)
    nothing
end
function _translate_cam!(scene, cam, t)
    # This uses a camera based coordinate system where
    # x expands right, y expands up and z expands towards the screen
    lookat = cam.lookat[]
    eyepos = cam.eyeposition[]
    up = normalize(cam.upvector[])
    u_z = normalize(eyepos - lookat)
    u_x = normalize(cross(up, u_z))
    u_y = normalize(cross(u_z, u_x))

    trans = u_x * t[1] + u_y * t[2] + u_z * t[3]

    # apply world space restrictions
    fix_x = ispressed(scene, cam.attributes[:fix_x_key][])
    fix_y = ispressed(scene, cam.attributes[:fix_y_key][])
    fix_z = ispressed(scene, cam.attributes[:fix_z_key][])
    if fix_x || fix_y || fix_z
        trans = Vec3f(fix_x, fix_y, fix_z) .* trans
    end

    cam.eyeposition[] = eyepos + trans
    cam.lookat[] = lookat + trans
    return
end

function rotate_cam!(scene, cam::Camera3D, angles::VecTypes, from_mouse=false)
    _rotate_cam!(scene, cam, angles, from_mouse)
    update_cam!(scene, cam)
    nothing
end
function _rotate_cam!(scene, cam::Camera3D, angles::VecTypes, from_mouse=false)
    # This applies rotations around the x/y/z axis of the camera coordinate system
    # x expands right, y expands up and z expands towards the screen
    lookat = cam.lookat[]
    eyepos = cam.eyeposition[]
    up = cam.upvector[]         # +y
    viewdir = lookat - eyepos   # -z
    right = cross(viewdir, up)  # +x

    x_axis = right
    y_axis = cam.attributes[:fixed_axis][] ? Vec3f(0, 0, sign(up[3])) : up
    z_axis = -viewdir

    fix_x = ispressed(scene, cam.attributes[:fix_x_key][])
    fix_y = ispressed(scene, cam.attributes[:fix_y_key][])
    fix_z = ispressed(scene, cam.attributes[:fix_z_key][])
    cx, cy, cz = cam.attributes[:circular_rotation][]
    rotation = Quaternionf(0, 0, 0, 1)
    if !xor(fix_x, fix_y, fix_z)
        # if there are more or less than one restriction apply all rotations
        rotation *= qrotation(y_axis, angles[2])
        rotation *= qrotation(x_axis, angles[1])
        rotation *= qrotation(z_axis, angles[3])
    else
        # apply world space restrictions
        if from_mouse && ((fix_x && (fix_x == cx)) || (fix_y && (fix_y == cy)) || (fix_z && (fix_z == cz)))
            # recontextualize the (dy, dx, 0) from mouse rotations so that
            # drawing circles creates continuous rotations around the fixed axis
            mp = mouseposition_px(scene)
            past_half = 0.5f0 .* widths(scene.px_area[]) .> mp
            flip = 2f0 * past_half .- 1f0
            angle = flip[1] * angles[1] + flip[2] * angles[2]
            angles = Vec3f(-angle, angle, -angle)
            # only one fix is true so this only rotates around one axis
            rotation *= qrotation(
                Vec3f(fix_x, fix_z, fix_y) .* Vec3f(sign(right[1]), viewdir[2], sign(up[3])),
                dot(Vec3f(fix_x, fix_y, fix_z), angles)
            )
        else
            # restrict total quaternion rotation to one axis
            rotation *= qrotation(y_axis, angles[2])
            rotation *= qrotation(x_axis, angles[1])
            rotation *= qrotation(z_axis, angles[3])
            # the first three components are related to rotations around the x/y/z-axis
            rotation = Quaternionf(rotation.data .* (fix_x, fix_y, fix_z, 1))
        end
    end

    cam.upvector[] = rotation * up
    viewdir = rotation * viewdir

    # TODO maybe generalize this to arbitrary center?
    # calculate positions from rotated vectors
    if cam.attributes[:rotation_center][] == :lookat
        cam.eyeposition[] = lookat - viewdir
    else
        cam.lookat[] = eyepos + viewdir
    end
    return
end

"""
    zoom!(scene, zoom_step)

Zooms the camera in or out based on the multiplier `zoom_step`. A `zoom_step`
of 1.0 is neutral, larger zooms out and lower zooms in.

Note that this method only applies to Camera3D.
"""
zoom!(scene::Scene, zoom_step) = zoom!(scene, cameracontrols(scene), zoom_step, false, false)
function zoom!(scene::Scene, cam::Camera3D, zoom_step, shift_lookat = false, cad = false)
    _zoom!(scene, cam, zoom_step, shift_lookat, cad)
    update_cam!(scene, cam)
    nothing
end
function _zoom!(scene::Scene, cam::Camera3D, zoom_step, shift_lookat = false, cad = false)
    if cad
        # move exeposition if mouse is not over the center
        lookat = cam.lookat[]
        eyepos = cam.eyeposition[]
        up = cam.upvector[]         # +y
        viewdir = lookat - eyepos   # -z
        right = cross(viewdir, up)  # +x

        rel_pos = 2f0 * mouseposition_px(scene) ./ widths(scene.px_area[]) .- 1f0
        shift = rel_pos[1] * normalize(right) + rel_pos[2] * normalize(up)
        shifted = eyepos + 0.1f0 * sign(1f0 - zoom_step) * norm(viewdir) * shift
        cam.eyeposition[] = lookat + norm(viewdir) * normalize(shifted - lookat)
    elseif shift_lookat
        lookat = cam.lookat[]
        eyepos = cam.eyeposition[]
        up = normalize(cam.upvector[])
        viewdir = lookat - eyepos
        u_z = normalize(-viewdir)
        u_x = normalize(cross(up, u_z))
        u_y = normalize(cross(u_z, u_x))

        if cam.attributes[:projectiontype][] == Perspective
            # translate both eyeposition and lookat to more or less keep data
            # under the mouse in view
            fov = cam.attributes[:fov][]
            before = tan(clamp(cam.zoom_mult[] * fov, 0.01f0, 175f0) / 360f0 * Float32(pi))
            after  = tan(clamp(cam.zoom_mult[] * zoom_step * fov, 0.01f0, 175f0) / 360f0 * Float32(pi))

            aspect = Float32((/)(widths(scene.px_area[])...))
            rel_pos = 2f0 * mouseposition_px(scene) ./ widths(scene.px_area[]) .- 1f0
            shift = rel_pos[1] * u_x + rel_pos[2] * u_y
            shift = -(after - before) * norm(viewdir) * normalize(aspect .* shift)
        else
            mx, my = 2f0 * mouseposition_px(scene) ./ widths(scene.px_area[]) .- 1f0
            aspect = Float32((/)(widths(scene.px_area[])...))
            w = 0.5f0 * (1f0 + aspect) * cam.zoom_mult[]
            h = 0.5f0 * (1f0 + 1f0 / aspect) * cam.zoom_mult[]
            shift = (1f0 - zoom_step) * (mx * w * u_x + my * h * u_y)
        end

        cam.lookat[]      = lookat + shift
        cam.eyeposition[] = eyepos + shift
    end

    # apply zoom
    cam.zoom_mult[] = cam.zoom_mult[] * zoom_step

    return
end

function update_cam!(scene::Scene, cam::Camera3D)
    @extractvalue cam (lookat, eyeposition, upvector)

    near = cam.near[]; far = cam.far[]
    aspect = Float32((/)(widths(scene.px_area[])...))

    if cam.attributes[:projectiontype][] == Perspective
        fov = clamp(cam.zoom_mult[] * cam.attributes[:fov][], 0.01f0, 175f0)
        cam.fov[] = fov
        proj = perspectiveprojection(fov, aspect, near, far)
    else
        w = 0.5f0 * (1f0 + aspect) * cam.zoom_mult[]
        h = 0.5f0 * (1f0 + 1f0 / aspect) * cam.zoom_mult[]
        proj = orthographicprojection(-w, w, -h, h, near, far)
    end

    view = Makie.lookat(eyeposition, lookat, upvector)

    scene.camera.projection[] = proj
    scene.camera.view[] = view
    scene.camera.projectionview[] = proj * view
    scene.camera.eyeposition[] = cam.eyeposition[]
end

function update_cam!(scene::Scene, camera::Camera3D, area3d::Rect)
    @extractvalue camera (lookat, eyeposition, upvector)
    bb = Rect3f(area3d)
    width = widths(bb)
    half_width = width/2f0
    middle = maximum(bb) - half_width
    old_dir = normalize(eyeposition .- lookat)
    camera.lookat[] = middle
    neweyepos = middle .+ (1.2*norm(width) .* old_dir)
    camera.eyeposition[] = neweyepos
    camera.upvector[] = Vec3f(0,0,1)
    if camera.attributes[:near][] === automatic
        camera.near[] = 0.1f0 * norm(widths(bb))
    end
    if camera.attributes[:far][] === automatic
        camera.far[]  = 3f0 * norm(widths(bb))
    end
    if camera.attributes[:projectiontype][] == Orthographic
        camera.zoom_mult[] = 0.6 * norm(width)
    else
        camera.zoom_mult[] = 1f0
    end
    update_cam!(scene, camera)
    return
end

function update_cam!(scene::Scene, camera::Camera3D, eyeposition, lookat, up = Vec3f(0, 0, 1))
    camera.lookat[] = Vec3f(lookat)
    camera.eyeposition[] = Vec3f(eyeposition)
    camera.upvector[] = Vec3f(up)
    update_cam!(scene, camera)
    return
end


text:
#### functions on `ThirdOrderExpFamily`
function Base.length(d::ThirdOrderExpFamily)::Int
    return d.nnode
end

function Distributions._logpdf(d::ThirdOrderExpFamily, x::AbstractVector)
    # unnormalized log probability (negative free energy)
    # and L2 norm square is added to ensure the functions is normalizable
    return -potential(d, x) - d.alpha * sum(x.^2)^2
end

function Distributions.insupport(d::ThirdOrderExpFamily, x::AbstractVector)
    return true
end

function potential(d::ThirdOrderExpFamily, x::AbstractVector)
    # free energy
    res = 0.0
    # singleton
    for i in 1:d.nnode
        res += ef1(x[i], d.node_attr[i])
    end
    # higher order
    for i in 1:d.nfactor
        if d.edge_type[i] == 2
            res += ef2(x[d.edge_index[i]], d.edge_attr[i])
        elseif d.edge_type[i] == 3
            res += ef3(x[d.edge_index[i]], d.edge_attr[i])
        else
            error("edge type note supported")
        end
    end

    return res
end

function ef1(x::Real, attr::AbstractVector)
    ans = 0.0
    for (i, a) in enumerate(attr)
        ans += a * x^i
    end
    return ans
end

function ef2(x::AbstractVector, attr::AbstractVector)
    @assert length(x) == 2 "len error, should be 2"
    @assert length(attr) == 1 "len error, should be 1"
    return attr[1] * prod(x)
end

function ef3(x::AbstractVector, attr::AbstractVector)
    @assert length(x) == 3 "len error, should be 3"
    @assert length(attr) == 1 "len error, should be 1"
    return attr[1] * prod(x)
end

function dist2dict(d::ThirdOrderExpFamily)
    # convert a distribution to a dictionary for serialization
    dic = Dict()
    dic["node_attr"] = d.node_attr
    dic["edge_attr"] = d.edge_attr
    dic["edge_type"] = d.edge_type
    dic["edge_index"] = d.edge_index
    dic["nnode"] = d.nnode
    dic["nfactor"] = d.nfactor
    dic["alpha"] = d.alpha
    return dic
end

function dist2dict(d::DiscreteSpinGlass)
    dic = Dict()
    dic["node_attr"] = d.node_attr
    dic["edge_attr"] = d.edge_attr
    dic["edge_type"] = d.edge_type
    dic["edge_index"] = d.edge_index
    dic["nnode"] = d.nnode
    dic["nfactor"] = d.nfactor
    dic["beta"] = d.beta
    return dic
end

## -----------------------------   SingleThirdOrderExpFamily
function Base.length(d::SingleThirdOrderExpFamily)::Int
    return 1
end

function Distributions._logpdf(d::SingleThirdOrderExpFamily, x::AbstractVector)
    @assert length(x) == 1 "size incorrect"
    return -potential(d, x) - d.alpha * x[1]^4
end

function Distributions.insupport(d::SingleThirdOrderExpFamily, x::AbstractVector)
    @assert length(x) == 1 "size incorrect"
    return true
end

function potential(d::SingleThirdOrderExpFamily, x::AbstractVector)
    @assert length(x) == 1 "size incorrect"
    ans = 0.0
    for i in 1:length(d.node_attr)
        ans += d.node_attr[i] * x[1]^i
    end
    return ans
end

## -----------------------------   DiscreteSpinGlass
function Base.length(d::DiscreteSpinGlass)::Int
    return d.nnode
end

function Distributions._logpdf(d::DiscreteSpinGlass, x::AbstractVector)
    return -potential(d, x) * d.beta
end

function potential(d::DiscreteSpinGlass, x::AbstractVector)
    ans = 0.0
    if !has_singleton_factor(d)      # if no singleton factor
        ans += sum([ef1(x[i],attr) for (i,attr) in enumerate(d.node_attr)])
        # ans += sum(d.node_attr .* x)  # bias term calculated from node_attr
    end
    for i in 1:d.nfactor
        ans += d.edge_attr[i][1] * prod(x[d.edge_index[i]])
    end
    return ans
end

function Distributions.insupport(d::DiscreteSpinGlass, x::AbstractVector)
    return true
end

## ------------------------------ MCMC sampling function -----
function sampleSingleChain(nstep::Integer, d::ContinuousDistribution, burnin::Integer = 3000,
        target::Real = 0.8)
    # Sample a single chain in one process
    # Return: Array of shape(nsample,nvariable)
    nvar = length(d)
    logp = x->logpdf(d, x)
    gradlogp = x->ForwardDiff.gradient(logp, x)
    logfgrad = x->(logp(x), gradlogp(x))
    theta = NUTSVariate(rand(Uniform(-100.0, 100.0), nvar), logfgrad; target = target)
    ans = zeros(nstep - burnin, nvar)

    for i in 1:nstep
        sample!(theta, adapt = (i <= burnin))
        if i > burnin
            ans[i - burnin,:] = theta
        end
    end
    return ans
end

function sampleMultiChain(nchain::Integer, nstep::Integer, d::ContinuousDistribution, burnin::Integer = 3000,
        target::Real = 0.8; distributed = true)
    sim = Chains(nstep, length(d), start = (burnin + 1), chains = nchain)
    args = []
    results = pmap(sampleSingleChain, repeat([nstep], nchain), repeat([d], nchain), repeat([burnin], nchain), repeat([target], nchain);distributed = distributed)
    for i in 1:length(results)
        sim[:,:,i] = results[i]
    end
    return sim
end

function sampleMultiChain2(model::Model, nvar::Integer, nstep::Integer, burnin::Integer = 5000, nchain::Integer = 6)
    inits = [
        Dict(:x => rand(Uniform(-10, 10), nvar)) for _ in 1:nchain
    ]
    scheme = [NUTS([:x])]
    setsamplers!(model, scheme)
    sim = mcmc(model, Dict{Symbol,Any}(), inits, nstep, burnin = burnin, chains = nchain;verbose = false)
    return sim
end


text:
using Unitful

export BrukerRobot, BrukerRobotParams

export moveCenter, movePark

Base.@kwdef struct BrukerRobotParams <: DeviceParams
  connectionName::String = "RobotServer"
  axisRange::Vector{Vector{typeof(1.0u"mm")}} = [[-85.0,225.0], [-Inf, Inf], [-Inf, Inf]]u"mm"
  namedPositions::Dict{String, Vector{typeof(1.0u"mm")}} = Dict("origin" => [0,0,0]u"mm")
end

BrukerRobotParams(dict::Dict) = params_from_dict(BrukerRobotParams, dict)

Base.@kwdef mutable struct BrukerRobot <: Robot
  "Unique device ID for this device as defined in the configuration."
  deviceID::String
  "Parameter struct for this devices read from the configuration."
  params::BrukerRobotParams
  "Flag if the device is optional."
  optional::Bool = false
  "Flag if the device is present."
  present::Bool = false
  "Vector of dependencies for this device."
  dependencies::Dict{String, Union{Device, Missing}}
  state::RobotState = INIT
end

""" `BrukerCommand(command::String)` """
struct BrukerCommand
  command::String
end

# Useful links
# http://unix.stackexchange.com/questions/87831/how-to-send-keystrokes-f5-from-terminal-to-a-process
#
# https://github.com/JuliaLang/julia/pull/6948
#
# Kind of deprecated:
# http://blog.leahhanson.us/post/julia/julia-commands.html

const center = "center\n"
const park = "park\n"
const pos = "pos\n"
const quit = "quit\n"
const exit = "exit\n"
const err = "err?\n"

dof(rob::BrukerRobot) = 3
axisRange(rob::BrukerRobot) = rob.params.axisRange
defaultVelocity(rob::BrukerRobot) = nothing
_getPosition(rob::BrukerRobot) = sendCommand(rob, BrukerCommand(pos))
_isReferenced(rob::BrukerRobot) = true
_enable(rob::BrukerRobot) = nothing
_disable(rob::BrukerRobot) = nothing
_reset(rob::BrukerRobot) = nothing
_setup(rob::BrukerRobot) = nothing
_doReferenceDrive(rob::BrukerRobot) = nothing

""" Move Bruker Robot to center"""
function moveCenter(sd::BrukerRobot)
    sendCommand(sd, BrukerCommand(center))
end

""" Move Bruker Robot to park"""
function movePark(sd::BrukerRobot)
    sendCommand(sd, BrukerCommand(park))
end

function _moveAbs(rob::BrukerRobot, pos::Vector{<:Unitful.Length}, speed::Union{Vector{<:Unitful.Velocity},Nothing})
    if speed !== nothing
        @warn "BrukerRobot does not support setting velocities!"
    end
    cmd = BrukerCommand("goto $(ustrip(Float64, u"mm", pos[1])),$(ustrip(Float64, u"mm", pos[2])),$(ustrip(Float64, u"mm", pos[3]))\n")
    res = sendCommand(rob, cmd)
end

""" Not Implemented """
function _moveRel(rob::BrukerRobot, dist::Vector{<:Unitful.Length}, speed::Union{Vector{<:Unitful.Velocity},Nothing})
    error("BrukerRobot does not support moveRel")
end

""" Send Command `sendCommand(sd::BrukerRobot, brukercmd::BrukerCommand)`"""
function sendCommand(sd::BrukerRobot, brukercmd::BrukerCommand)

    (result, startmovetime, endmovetime) = _sendCommand(sd, brukercmd)

    if result == "0\n"
        return true
    elseif length(split(result, ",")) == 3
        @info "$(brukercmd.command) returned position: $(result)"
        return parse.(Float64, split(result, ",")) * u"mm"
    elseif result == "?\n"
        @warn "$(brukercmd.command) is unknown! Try again..."
        return false
    elseif result == "!\n"
        throw(ErrorException("Error during command $(brukercmd.command) execution."))
    else
        throw(ErrorException("$(brukercmd.command) has unexpected result $(result)"))
    end
end

function _sendCommand(sd::BrukerRobot, brukercmd::BrukerCommand)
    p = open(`$(sd.params.connectionName)`, "r+");
  # p = open(`cat`,"r+");
    startmovetime = now(Dates.UTC)
    writetask = write(p.in, brukercmd.command)
    writetaskexit = write(p.in, exit)
    readtask = @async readavailable(p.out)
    wait(readtask)
    endmovetime = now(Dates.UTC)
    if readtask.state == :done
        return (ascii(String(readtask.result)), startmovetime, endmovetime);
    else
        @error "end" readTask.state readTask.exception
        return readtask.exception
    end
end


text:
# test mesh adaptation

function test_adapt_2d(;parallel=false)

  @testset "Testing 2D mesh adaptation" begin

    if parallel
      smb_name = "parallel2.smb"
    else
      smb_name = "tri8l.smb"
    end
    order = 1
    opts = Dict{Any, Any}(
      "dmg_name" => ".null",
      "smb_name" => smb_name,
      "order" => order,
      "coloring_distance" => 2,
      "numBC" => 1,
      "BC1" =>  [0],
      "run_type" => 4)

    PdePumiInterface.set_defaults(opts)

    sbp = getTriSBPOmega(degree=order)
    vtx = sbp.vtx
    sbpface = TriFace{Float64}(order, sbp.cub, vtx)

    mesh = PumiMeshDG2(Float64, sbp, opts, sbpface, dofpernode=1)

    # write linear field to mesh
    u = zeros(mesh.numDof)
    for i=1:mesh.numEl
      for j=1:mesh.numNodesPerElement
        x = mesh.coords[1, j, i]; y = mesh.coords[2, j, i]
        u[mesh.dofs[1, j, i]] = x + y + 1
      end
    end

    el_sizes = getElementSizes(mesh)
    scale!(el_sizes, 0.5)

    numEl_initial = mesh.numEl

    newmesh, unew = adaptMesh(mesh, sbp, opts, el_sizes, u)

    @test newmesh.numEl == 4*numEl_initial

    for i=1:newmesh.numEl
      for j=1:newmesh.numNodesPerElement
        x = newmesh.coords[1, j, i]; y = newmesh.coords[2, j, i]
        @test isapprox(unew[newmesh.dofs[1, j, i]], x + y + 1) atol=1e-13
      end
    end

  end  # end testset

  return nothing
end

function test_adapt_3d(; parallel=false)

  @testset "Testing 3D mesh adaptation" begin
    degree = 1
    dmg_name = ".null"
    if parallel
      smb_name = "pcube10.smb"
    else
      smb_name = "unitcube.smb"
    end

    opts = PdePumiInterface.get_defaults()
    opts["smb_name"] = smb_name
    opts["dmg_name"] = dmg_name
    opts["order"] = degree
    opts["coloring_distance"] = 2
    opts["use_linear_metrics"] = true
    opts["numBC"] = 1
    opts["BC1"] = [0,1,2,3,4,5]

    sbp = getTetSBPOmega(degree=degree, Tsbp=Float64)
    ref_verts = sbp.vtx
    face_verts = SummationByParts.SymCubatures.getfacevertexindices(sbp.cub)
    topo = ElementTopology{3}(face_verts)
    sbpface = TetFace{Float64}(degree, sbp.cub, ref_verts)

    mesh = PumiMeshDG3(Float64, sbp, opts, sbpface, topo)

    # write linear field to mesh
    println("writing field")
    u = zeros(mesh.numDof)
    for i=1:mesh.numEl
      for j=1:mesh.numNodesPerElement
        x = mesh.coords[1, j, i]; y = mesh.coords[2, j, i]; z = mesh.coords[3, j, i]
        u[mesh.dofs[1, j, i]] = x + 2*y + 3*z + 1
      end
    end

    el_sizes = getElementSizes(mesh)
    scale!(el_sizes, 0.5)

    numEl_initial = mesh.numEl

    println("adapting mesh")
    newmesh, unew = adaptMesh(mesh, sbp, opts, el_sizes, u)

    @test newmesh.numEl == 8*numEl_initial

    # check that the linear field was interpolated exactly
    for i=1:newmesh.numEl
      for j=1:newmesh.numNodesPerElement
        x = newmesh.coords[1, j, i]; y = newmesh.coords[2, j, i]; z = newmesh.coords[3, j, i]
        @test isapprox(unew[newmesh.dofs[1, j, i]], x + 2*y + 3*z + 1) atol=1e-13
      end
    end

  end  # end testset

  return nothing
end


text:
module HigherOrderClustering

export hoccf_data, higher_order_ccfs, kcliques
export read_undir_graph_txt, load_example_data

include("io.jl")
include("cliques.jl")

"""
hoccf_data
----------

This is the data type returned by the function that computes higher-order clustering coefficients.

The field values are

order::Int64
    The order of the clustering coefficient (order = 2 is the classical
    clustering coefficient definition).

global_hoccf::Float64
    The global higher-order clustering coefficient.

avg_hoccf::Float64
    The average higher-order clustering coefficient (the mean is taken over
    nodes at the center of at least one wedge).

avg_hoccf2::Float64
     The average higher-order clustering coefficient, where the local clustering
    of a node not in any wedge is considered to be 0.

local_hoccfs::Vector{Float64}
    The vector of local higher-order clustering coefficients. If a node is not
    the center of at least one wedge, then the value is 0.

ho_wedge_counts::Vector{Int64}
    Higher-order wedge counts of nodes: ho_wedge_counts[v] is the number of higher-order
    wedges with node v at its center.

clique_counts::Vector{Int64}
    Clique counts of nodes: clique_counts[v] is the number of k-cliques containing
    node v where k = order + 1.
"""
immutable hoccf_data
    order::Int64
    global_hoccf::Float64
    avg_hoccf::Float64
    avg_hoccf2::Float64
    local_hoccfs::Vector{Float64}
    ho_wedge_counts::Vector{Int64}
    clique_counts::Vector{Int64}
end

"""
higher_order_ccfs
-----------------

This is the function that computes higher-order clustering coefficients.
It takes two inputs:

A::SparseMatrixCSC{Int64,Int64}
    The adjacency matrix of an undirected graph.

l::Int64
    The order of the clustering coefficient (l = 2 is the classical definition).
"""
function higher_order_ccfs(A::SparseMatrixCSC{Int64,Int64}, l::Int64)
    A = min.(A, 1)
    A -= spdiagm(diag(A))
    n = size(A, 1)
    # Get clique counts
    clique_counts1 = kcliques(A, l)
    clique_counts2 = kcliques(A, l + 1)
    degs = vec(sum(A, 2))
    # normalize
    wedge_counts = (degs - l + 1) .* clique_counts1
    nz_inds = find(wedge_counts .> 0)
    local_hoccfs = zeros(Float64, n)
    local_hoccfs[nz_inds] = l * clique_counts2[nz_inds] ./ wedge_counts[nz_inds]
    return hoccf_data(l, l * sum(clique_counts2) / sum(wedge_counts),
                      mean(local_hoccfs[nz_inds]), mean(local_hoccfs),
                      local_hoccfs, wedge_counts, clique_counts2)
end

end # module


text:
using SegmentIntersections, Test

@testset "create point" begin
    p = Point(1,2)
    @test p.x == 1.0
    @test typeof(p.x) == Float64
    @test p.y == 2.0
    @test typeof(p.y) == Float64
    p = Point(1.0,2.0)
    @test p.x == 1.0
    @test typeof(p.x) == Float64
    @test p.y == 2.0
    @test typeof(p.y) == Float64
end

@testset "Test point order" begin
    @testset "Test different y" begin
        p = Point(1,2)
        q = Point(3,4)
        @test q < p
        @test q <= p
        @test !(q > p)
        @test !(q >= p)
        @test !(q == p)
    end
    @testset "Test same y" begin
        p = Point(1,2)
        q = Point(3,2)
        @test !(q < p)
        @test !(q <= p)
        @test (q > p)
        @test (q >= p)
        @test !(q == p)
    end
end


text:
# sample sequences of ergms with different parameters' values from R package ergm
# and test the PseudoLikelihoodScoreDrivenERGM filter
using Utilities,AReg,StaticNets,JLD,MLBase,StatsBase,CSV, RCall
using PyCall; pygui(); using PyPlot

using JLD2,Utilities, GLM
## load R MCMC simulation and estimates and estimate sdergmTest

Nsample = 50
 dgpType = "steps"
 T = 50
 N = 50
 Nterms = 2
 Nsteps1 ,Nsteps2 = 2,1
 load_fold = "./data/estimatesTest/sdergmTest/R_MCMC_estimates/"
 @load(load_fold*"test_Nodes_$(N)_T_$(T)_Sample_$(Nsample)_Ns_" * dgpType * "_$(Nsteps1)_$(Nsteps2)_MPLE.jld",
             stats_T, changeStats_T,estParSS_T,sampledMat_T ,parDgpT,Nsample)

onlyTest = false
 targetAllTv = true
 tvParFromTest =true; pValTh = 0.05
 gasParEst = fill(fill(Float64[],2), Nsample)
 startPointEst = fill(fill(0.0,2), Nsample)
 staticEst = fill(fill(0.0,2), Nsample)
 filtPar_T_Nsample = zeros(Nterms,T,Nsample)
 pVals_Nsample = zeros(Nterms,Nsample,2)
 scoreAutoc_Nsample = zeros(Nterms,Nsample)
 convFlag = falses(Nsample)
 tmp = Array{Array{Float64,2},2}(T,Nsample); for t=1:T,s=1:Nsample tmp[t,s] =  changeStats_T[t][s];end;#changeStats_T = tmp
 #for
     n=1#:Nsample
     @show(n)
      pVals_Nsample[:,n,1],tmpInfo,staticEst[n] =
         pValStatic_SDERGM( SdErgmDirBinGlobalPseudo(tmp[:,n],fooGasPar,falses(Nterms),""))
      @show(tmpInfo)
     if tvParFromTest
         indTvPar =  pVals_Nsample[:,n] .< pValTh
     else
         indTvPar = BitArray([true,true])
     end
     if targetAllTv
         indTargPar = indTvPar
     else
          indTargPar =BitArray([false,false])
     end
    model = SdErgmDirBinGlobalPseudo(tmp[:,n],fooGasPar,indTvPar,"")

    if .!onlyTest
    estPar,convFlag[n],UM,startPoint = estimate(model;UM =meanSq(estParSS_T[n,:,:],1),indTvPar = indTvPar, indTargPar = indTargPar)

    gasParEst[n] = estPar # store gas parameters
    startPointEst[n] = startPoint
    gasParVec = zeros(sum(indTvPar)*3); for i=0:(sum(indTvPar)-1) gasParVec[1+3i : 3(i+1)] = estPar[indTvPar][i+1]; end
    constParVec = zeros(sum(!indTvPar)); for i=1:sum(.!indTvPar) constParVec[i] = estPar[.!indTvPar][i][1]; end
    gasFiltPar , pseudolike = score_driven_filter_or_dgp(model,gasParVec,indTvPar;
                              vConstPar = constParVec,ftot_0= startPoint )

    filtPar_T_Nsample[:,:,n] = gasFiltPar
    pVals_Nsample[:,n,2],tmpInfo,~ =
       pValStatic_SDERGM( SdErgmDirBinGlobalPseudo(tmp[:,n],fooGasPar,falses(Nterms),"");estPar = gasFiltPar)
    #end
 end
 close("all")

tmp = gasScoreSeries(model,gasFiltPar;matScal=false)
plot(tmp[1]')
autocor(Float64.(tmp[1])',[1])
pValStatic_SDERGM( SdErgmDirBinGlobalPseudo(tmp[:,n],fooGasPar,falses(Nterms),""))

pValStatic_SDERGM( SdErgmDirBinGlobalPseudo(tmp[:,n],fooGasPar,falses(Nterms),"");estPar = gasFiltPar)
  # Estimation
#Plotta info relative ai tests
for testInt =1:2
 figure()
 namePar1 = "Number of Links"
 namePar2 = "GWESP"
 subplot(1,2,1);    title(namePar1)
 subplot(1,2,2);  title(namePar2 );# legend(legTex)

 pVals_Nsample[isnan(pVals_Nsample)] = pValTh
  pValTh = 0.05
      subplot(2,2,1); parInd = 1;   plt[:hist]((pVals_Nsample[parInd,:,testInt]), bins=logspace(minimum(log10(pVals_Nsample[1,:,testInt])),maximum(log10(pVals_Nsample[parInd,:,testInt])), 20))
      if maximum((pVals_Nsample[parInd,:,testInt])) >pValTh
           axvspan((pValTh),maximum((pVals_Nsample[parInd,:,testInt])),color = "r",alpha = 0.1);
      end
           xscale("log");
      subplot(2,2,2);  parInd = 2;   plt[:hist]((pVals_Nsample[parInd,:,testInt]), bins=logspace(minimum(log10(pVals_Nsample[2,:,testInt])),maximum(log10(pVals_Nsample[parInd,:,testInt])), 20))
      if maximum((pVals_Nsample[parInd,:])) >pValTh
           axvspan((pValTh),maximum((pVals_Nsample[parInd,:,testInt])),color = "r",alpha = 0.1);
      end
           xscale("log");
      subplot(2,2,3); plot((pVals_Nsample[1,:,testInt]),(pVals_Nsample[2,:,testInt]),".")
      axvline((pValTh));axhline((pValTh))
      xscale("log");yscale("log")
      xlabel("P-Value " * namePar1)
      ylabel("P-Value " * namePar2)
end


text:
using ApproxFun, QuantumTimeSteppers, Plots, MTFun

F = Fourier(-4..4)
u0F = Fun(x -> exp(10im*x)*exp(-30*(x+2.0)^2),F,1024)
Δt = 0.01
ε = 0.01
V = x -> .5*x^2
T = 20
tF, uF = Strang_evolve(u0F,V,ε,Δt,T,0);
@gif for k=1:5:length(tF)
    plot(uF[k],ylims=(-2,2),title=string("t = ",tF[k]))
end

MT = MalmquistTakenaka()
u0MT = Fun(x -> exp(10im*x)*exp(-30*(x+2.0)^2),MT,1024)
Δt = 0.01
ε = 0.01
V = x -> .5*x^2
T = 20
tMT, uMT = Strang_evolve(u0MT,V,ε,Δt,T,6);
@gif for k=1:5:length(tMT)
    plot(uMT[k],xlims=(-4,4),ylims=(-2,2),title=string("t = ",tMT[k]))
end

H = HermiteFSE()
u0H = Fun(x -> exp(10im*x)*exp(-.5*(x+2.0)^2),H,512)
Δt = 0.01
ε = 0.01
V = x -> .5*x^2
T = 10
tH, uH = Strang_evolve(u0H,V,ε,Δt,T,6);
@gif for k=1:5:length(tH)
    plot(uH[k],xlims=(-4,4),ylims=(-2,2),title=string("t = ",tH[k]))
end


text:
using Test
using ShowSet

A = BitSet([1,3,9])
@test string(A) == "{1,3,9}"

B = Set{String}()
push!(B,"Hello")
push!(B,"Bye")
@test "{Bye,Hello}" == string(B)