text:
@timed_testset "data_constructors" begin
    @testset "data_response" begin
        @testset "NamedTuples" begin
            expected = [2, 3, 4, 5]
            f = @formula y_int ~ 0 + x_float + x_cat
            y = T.data_response(f, nt_str)
            @test y == expected
            y = T.data_response(f, nt_cat)
            @test y == expected
            f = @formula y_int ~ 1 + x_float + x_cat
            y = T.data_response(f, nt_str)
            @test y == expected
            y = T.data_response(f, nt_cat)
            @test y == expected

            expected = [2.3, 3.4, 4.5, 5.4]
            f = @formula y_float ~ 0 + x_float + x_cat
            y = T.data_response(f, nt_str)
            @test y == expected
            y = T.data_response(f, nt_cat)
            @test y == expected
            f = @formula y_float ~ 1 + x_float + x_cat
            @test y == expected
            y = T.data_response(f, nt_str)
            @test y == expected
            y = T.data_response(f, nt_cat)
            @test y == expected
        end
        @testset "DataFrames" begin
            expected = [2, 3, 4, 5]
            f = @formula y_int ~ 0 + x_float + x_cat
            y = T.data_response(f, df_str)
            @test y == expected
            y = T.data_response(f, df_cat)
            @test y == expected
            f = @formula y_int ~ 1 + x_float + x_cat
            y = T.data_response(f, df_str)
            @test y == expected
            y = T.data_response(f, df_cat)
            @test y == expected

            expected = [2.3, 3.4, 4.5, 5.4]
            f = @formula y_float ~ 0 + x_float + x_cat
            y = T.data_response(f, df_str)
            @test y == expected
            y = T.data_response(f, df_cat)
            @test y == expected
            f = @formula y_float ~ 1 + x_float + x_cat
            y = T.data_response(f, df_str)
            @test y == expected
            y = T.data_response(f, df_cat)
            @test y == expected
        end
    end
    @testset "data_fixed_effects" begin
        @testset "NamedTuples" begin
            expected = [
                1.1 0.0 0.0 0.0
                2.3 1.0 0.0 0.0
                3.14 0.0 1.0 0.0
                3.65 0.0 0.0 1.0
            ]

            f = @formula y_int ~ 0 + x_float + x_cat
            X = T.data_fixed_effects(f, nt_str)
            @test X == expected

            f = @formula y_int ~ 1 + x_float + x_cat
            X = T.data_fixed_effects(f, nt_str)
            @test X == expected

            f = @formula y_float ~ 0 + x_float + x_cat
            X = T.data_fixed_effects(f, nt_str)
            @test X == expected

            f = @formula y_float ~ 1 + x_float + x_cat
            X = T.data_fixed_effects(f, nt_str)
            @test X == expected

            f = @formula y_int ~ 0 + x_float + x_cat_ordered
            X = T.data_fixed_effects(f, nt_cat)
            @test X == expected

            f = @formula y_int ~ 1 + x_float + x_cat_ordered
            X = T.data_fixed_effects(f, nt_cat)
            @test X == expected

            f = @formula y_float ~ 0 + x_float + x_cat_ordered
            X = T.data_fixed_effects(f, nt_cat)
            @test X == expected

            f = @formula y_float ~ 1 + x_float + x_cat_ordered
            X = T.data_fixed_effects(f, nt_cat)
            @test X == expected

            # Interactions
            expected = [
                1.0 0.0 0.0 0.0 0.0 0.0 0.0
                2.0 1.0 0.0 0.0 2.0 0.0 0.0
                3.0 0.0 1.0 0.0 0.0 3.0 0.0
                4.0 0.0 0.0 1.0 0.0 0.0 4.0
            ]

            f = @formula y_float ~ 0 + x_int * x_cat
            X = T.data_fixed_effects(f, nt_str)
            @test X == expected

            f = @formula y_float ~ 1 + x_int * x_cat
            X = T.data_fixed_effects(f, nt_str)
            @test X == expected

            f = @formula y_float ~ 0 + x_int * x_cat
            X = T.data_fixed_effects(f, nt_cat)
            @test X == expected

            f = @formula y_float ~ 1 + x_int * x_cat
            X = T.data_fixed_effects(f, nt_cat)
            @test X == expected

            f = @formula y_float ~ 0 + x_int * x_cat_ordered
            X = T.data_fixed_effects(f, nt_cat)
            @test X == expected

            f = @formula y_float ~ 1 + x_int * x_cat_ordered
            X = T.data_fixed_effects(f, nt_cat)
            @test X == expected

            # Interactions coming first
            expected = [
                1.0 0.0 0.0 0.0 1.1 0.0 0.0 0.0
                2.0 1.0 0.0 0.0 2.3 2.0 0.0 0.0
                3.0 0.0 1.0 0.0 3.14 0.0 3.0 0.0
                4.0 0.0 0.0 1.0 3.65 0.0 0.0 4.0
            ]

            f = @formula y_float ~ 1 + x_int * x_cat + x_float
            X = T.data_fixed_effects(f, nt_str)
            @test X == expected
        end
        @testset "DataFrames" begin
            expected = [
                1.1 0.0 0.0 0.0
                2.3 1.0 0.0 0.0
                3.14 0.0 1.0 0.0
                3.65 0.0 0.0 1.0
            ]

            f = @formula(y_int ~ 0 + x_float + x_cat)
            X = T.data_fixed_effects(f, df_str)
            @test X == expected

            f = @formula(y_int ~ 1 + x_float + x_cat)
            X = T.data_fixed_effects(f, df_str)
            @test X == expected

            f = @formula(y_float ~ 0 + x_float + x_cat)
            X = T.data_fixed_effects(f, df_str)
            @test X == expected

            f = @formula(y_float ~ 1 + x_float + x_cat)
            X = T.data_fixed_effects(f, df_str)
            @test X == expected

            f = @formula(y_int ~ 0 + x_float + x_cat_ordered)
            X = T.data_fixed_effects(f, df_cat)
            @test X == expected

            f = @formula(y_int ~ 1 + x_float + x_cat_ordered)
            X = T.data_fixed_effects(f, df_cat)
            @test X == expected

            f = @formula(y_float ~ 0 + x_float + x_cat_ordered)
            X = T.data_fixed_effects(f, df_cat)
            @test X == expected

            f = @formula(y_float ~ 1 + x_float + x_cat_ordered)
            X = T.data_fixed_effects(f, df_cat)
            @test X == expected

            # Interactions
            expected = [
                1.0 0.0 0.0 0.0 0.0 0.0 0.0
                2.0 1.0 0.0 0.0 2.0 0.0 0.0
                3.0 0.0 1.0 0.0 0.0 3.0 0.0
                4.0 0.0 0.0 1.0 0.0 0.0 4.0
            ]

            f = @formula y_float ~ 0 + x_int * x_cat
            X = T.data_fixed_effects(f, df_str)
            @test X == expected

            f = @formula y_float ~ 1 + x_int * x_cat
            X = T.data_fixed_effects(f, df_str)
            @test X == expected

            f = @formula y_float ~ 0 + x_int * x_cat
            X = T.data_fixed_effects(f, df_cat)
            @test X == expected

            f = @formula y_float ~ 1 + x_int * x_cat
            X = T.data_fixed_effects(f, df_cat)
            @test X == expected

            f = @formula y_float ~ 0 + x_int * x_cat_ordered
            X = T.data_fixed_effects(f, df_cat)
            @test X == expected

            f = @formula y_float ~ 1 + x_int * x_cat_ordered
            X = T.data_fixed_effects(f, df_cat)
            @test X == expected

            # Interactions coming first
            expected = [
                1.0 0.0 0.0 0.0 1.1 0.0 0.0 0.0
                2.0 1.0 0.0 0.0 2.3 2.0 0.0 0.0
                3.0 0.0 1.0 0.0 3.14 0.0 3.0 0.0
                4.0 0.0 0.0 1.0 3.65 0.0 0.0 4.0
            ]

            f = @formula y_float ~ 1 + x_int * x_cat + x_float
            X = T.data_fixed_effects(f, df_str)
            @test X == expected
        end
    end

    @testset "data_random_effects" begin
        expected = nothing
        f = @formula y_float ~ 1 + x_int * x_cat + x_float
        Z = T.data_random_effects(f, nt_str)
        @test Z == expected
        Z = T.data_random_effects(f, nt_cat)
        @test Z == expected
        Z = T.data_random_effects(f, df_str)
        @test Z == expected
        Z = T.data_random_effects(f, df_cat)
        @test Z == expected

        f = @formula y_float ~ 1 + x_float + (1 | x_cat) + (1 | x_cat)
        Z = T.data_random_effects(f, nt_str)
        @test Z == expected
        Z = T.data_random_effects(f, nt_cat)
        @test Z == expected
        Z = T.data_random_effects(f, df_str)
        @test Z == expected
        Z = T.data_random_effects(f, df_cat)
        @test Z == expected

        expected = Dict("slope_x_int" => [1.0, 2.0, 3.0, 4.0])
        f = @formula y_float ~ 1 + (1 + x_int | x_cat)
        Z = T.data_random_effects(f, nt_str)
        @test Z == expected
        Z = T.data_random_effects(f, nt_cat)
        @test Z == expected
        Z = T.data_random_effects(f, df_str)
        @test Z == expected
        Z = T.data_random_effects(f, df_cat)
        @test Z == expected

        expected = Dict(
            "slope_x_float" => [1.1, 2.3, 3.14, 3.65], "slope_x_int" => [1.0, 2.0, 3.0, 4.0]
        )
        f = @formula y_float ~ 1 + (1 + x_int + x_float | x_cat)
        Z = T.data_random_effects(f, nt_str)
        @test Z == expected
        Z = T.data_random_effects(f, nt_cat)
        @test Z == expected
        Z = T.data_random_effects(f, df_str)
        @test Z == expected
        Z = T.data_random_effects(f, df_cat)
        @test Z == expected

        expected = Dict(
            "slope_x_float" => [1.1, 2.3, 3.14, 3.65], "slope_x_int" => [1.0, 2.0, 3.0, 4.0]
        )
        f = @formula y_float ~
            1 + (1 + x_int + x_float | x_cat) + (1 + x_int + x_float | group)
        Z = T.data_random_effects(f, nt_str)
        @test Z == expected
        Z = T.data_random_effects(f, nt_cat)
        @test Z == expected
        Z = T.data_random_effects(f, df_str)
        @test Z == expected
        Z = T.data_random_effects(f, df_cat)
        @test Z == expected
    end

    @testset "has_ranef" begin
        f = @formula y_float ~ 1 + x_int + x_cat
        @test T.has_ranef(f) == false

        f = @formula y_float ~ 1 + x_int + (1 | x_cat)
        @test T.has_ranef(f) == true

        f = @formula y_float ~ 0 + x_int + x_cat
        @test T.has_ranef(f) == false

        f = @formula y_float ~ 0 + x_int + (1 | x_cat)
        @test T.has_ranef(f) == true

        f = @formula y_float ~ x_int + x_cat
        @test T.has_ranef(f) == false

        f = @formula y_float ~ x_int + (1 | x_cat)
        @test T.has_ranef(f) == true
    end

    @testset "ranef" begin
        f = @formula y_float ~ 1 + x_int + x_cat
        @test T.ranef(f) === nothing

        f = @formula y_float ~ x_int + (1 | x_cat)
        @test T.ranef(f) isa Tuple{T.RandomEffectsTerm}

        f = @formula y_float ~ x_int + (1 + x_float | x_cat)
        @test T.ranef(f) isa Tuple{T.RandomEffectsTerm}
    end

    @testset "n_ranef" begin
        f = @formula y_float ~ 1 + x_int + x_cat
        @test T.n_ranef(f) == 0

        f = @formula y_float ~ x_int + (1 | x_cat)
        @test T.n_ranef(f) == 1

        f = @formula y_float ~ 1 + x_float + (1 | x_cat) + (1 | x_cat)
        @test T.n_ranef(f) == 1

        f = @formula y_float ~ 1 + x_float + (1 | x_cat) + (1 | group)
        @test T.n_ranef(f) == 2

        f = @formula y_float ~ x_int + (1 + x_float | x_cat)
        @test T.n_ranef(f) == 2

        f = @formula y_float ~ 1 + (1 + x_int + x_float | x_cat)
        @test T.n_ranef(f) == 3

        f = @formula y_float ~
            1 + (1 + x_int + x_float | x_cat) + (1 + x_int + x_float | group)
        @test T.n_ranef(f) == 6
    end

    @testset "intercept_per_ranef" begin
        f = @formula y_float ~ 1 + x_int + xcat + (1 | x_cat)
        @test T.intercept_per_ranef(T.ranef(f)) == ["x_cat"]

        f = @formula y_float ~ 1 + (1 + x_int + x_float | x_cat)
        @test T.intercept_per_ranef(T.ranef(f)) == ["x_cat"]

        f = @formula y_float ~ 1 + x_float + (1 | x_cat) + (1 | x_cat)
        @test T.intercept_per_ranef(T.ranef(f)) == ["x_cat"]

        f = @formula y_float ~ 1 + x_float + (1 | x_cat) + (1 | group)
        @test T.intercept_per_ranef(T.ranef(f)) == ["x_cat", "group"]

        f = @formula y_float ~ 1 + (1 + x_int + x_float | x_cat)
        @test T.intercept_per_ranef(T.ranef(f)) == ["x_cat"]

        f = @formula y_float ~
            1 + (1 + x_int + x_float | x_cat) + (1 + x_int + x_float | group)
        @test T.intercept_per_ranef(T.ranef(f)) == ["x_cat", "group"]
    end

    @testset "slope_per_ranef" begin
        f = @formula y_float ~ 1 + x_int + xcat + (1 | x_cat)
        @test T.slope_per_ranef(T.ranef(f)) == T.SlopePerRanEf()

        f = @formula y_float ~ 2 + (1 + x_int + x_float | x_cat)
        @test T.slope_per_ranef(T.ranef(f)) ==
            T.SlopePerRanEf(Dict("x_cat" => ["x_int", "x_float"]))

        f = @formula y_float ~ 1 + x_float + (1 | x_cat) + (1 | x_cat)
        @test T.slope_per_ranef(T.ranef(f)) == T.SlopePerRanEf()

        f = @formula y_float ~ 1 + x_float + (1 | x_cat) + (1 | group)
        @test T.slope_per_ranef(T.ranef(f)) == T.SlopePerRanEf()

        f = @formula y_float ~
            1 + (1 + x_int + x_float | x_cat) + (1 + x_int + x_float | group)
        @test T.slope_per_ranef(T.ranef(f)) == T.SlopePerRanEf(
            Dict("x_cat" => ["x_int", "x_float"], "group" => ["x_int", "x_float"])
        )
    end

    @testset "get_idx" begin
        expected = ([1, 2, 3, 4], Dict(1.1 => 1, 3.65 => 4, 2.3 => 2, 3.14 => 3))
        @test T.get_idx(T.term("x_float"), nt_str) == expected
        @test T.get_idx(T.term("x_float"), df_str) == expected
        @test T.get_idx(T.term("x_float"), nt_cat) == expected
        @test T.get_idx(T.term("x_float"), df_cat) == expected

        expected = ([1, 2, 3, 4], Dict("4" => 4, "1" => 1, "2" => 2, "3" => 3))
        @test T.get_idx(T.term("x_cat"), nt_str) == expected
        @test T.get_idx(T.term("x_cat"), df_str) == expected

        cv = categorical([1, 2, 3, 4])
        expected = (
            [1, 2, 3, 4],
            Dict(
                CategoricalValue(4, cv) => 4,
                CategoricalValue(2, cv) => 2,
                CategoricalValue(3, cv) => 3,
                CategoricalValue(1, cv) => 1,
            ),
        )
        @test T.get_idx(T.term("x_cat"), nt_cat) == expected
        @test T.get_idx(T.term("x_cat"), df_cat) == expected
    end

    @testset "get_var" begin
        expected = [1.1, 2.3, 3.14, 3.65]
        @test T.get_var(T.term("x_float"), nt_str) == expected
        @test T.get_var(T.term("x_float"), df_str) == expected

        expected = [1.1, 2.3, 3.14, 3.65]
        @test T.get_var(T.term("x_float"), nt_cat) == expected
        @test T.get_var(T.term("x_float"), df_cat) == expected

        expected = ["1", "2", "3", "4"]
        @test T.get_var(T.term("x_cat"), nt_str) == expected
        @test T.get_var(T.term("x_cat"), df_str) == expected

        expected = categorical([1, 2, 3, 4])
        @test T.get_var(T.term("x_cat"), nt_cat) == expected
        @test T.get_var(T.term("x_cat"), df_cat) == expected
    end
end


text:
using Pkg;
pkg"add BenchmarkTools, LaTeXStrings, LinearAlgebra, Plots, Printf, PyPlot, Random, Suppressor"

ext: jl


text:
#=
  From https://dtai.cs.kuleuven.be/problog/tutorial/mpe/01_bn.html
  (from ~/problog/bayesian_networks.pl )

  See ~/blog/bayesian_network.blog
      ~/psi/bayesian_network.psi
      ~/webppl/bayesian_network.wppl
=#

using Turing, StatsPlots, DataFrames
using ReverseDiff, Zygote, Tracker
# Turing.setadbackend(:reversediff)
# Turing.setadbackend(:zygote)
# Turing.setadbackend(:tracker)

include("jl_utils.jl")

@model function bayesian_network()

    burglary ~ flip(0.7)
    earthquake ~ flip(0.2)

    alarm ~ flip(0.5)

    if burglary && earthquake
        alarm ~ flip(0.9)
    elseif burglary && !earthquake
        alarm ~ flip(0.8)
    elseif !burglary && earthquake
        alarm ~ flip(0.1)
    else
        alarm ~ flip(0.0); # false
    end

    # x calls depending on the alarm
    # But using function don't seems to work!
    # function calls(x)
    #    return alarm ? flip(0.8) : flip(0.1)
    # end
    # calls(john) == false ||  begin Turing.@addlogprob! -Inf; return; end
    # calls(mary) == trye ||  begin Turing.@addlogprob! -Inf; return; end

    calls_john ~ alarm ? flip(0.8) : flip(0.1)
    calls_mary ~ alarm ? flip(0.8) : flip(0.1)

    true ~ Dirac(calls_john == false)
    true ~ Dirac(calls_mary == true)

    return
end

model = bayesian_network()
num_chains = 4
# chains = sample(model, Prior(), MCMCThreads(), 1000, num_chains)
chains = sample(model, MH(), MCMCThreads(), 40_000, num_chains)
# chains = sample(model, PG(20), MCMCThreads(), 10_000, num_chains)
# chains = sample(model, IS(), MCMCThreads(), 10_000, num_chains)

# chains = sample(model, NUTS(1000,0.65), MCMCThreads(), 40_000, num_chains)
# chains = sample(model, Gibbs(MH(:zlabels),NUTS(1000,0.65,:m,:b,:sigma)), MCMCThreads(), 40_000, num_chains)

display(chains)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
# This file is part of TaylorSeries.jl, MIT licensed
#

using TaylorSeries

using Test

@testset "Broadcasting with Taylor1 expansions" begin
    t = Taylor1(Int, 5)

    # @test t .= t
    @test t .== t
    @test t .≈ t
    @test t .!= (1 + t)

    @test 1.0 .+ t == 1.0 + t
    @test typeof(1.0 .+ t) == Taylor1{Float64}
    @test 1.0 .+ [t] == [1.0 + t]
    @test typeof(1.0 .+ [t]) == Vector{Taylor1{Float64}}
    @test 1.0 .+ [t, 2t] == [1.0 + t, 1.0 + 2t]
    @test [1.0,2.0] .+ [t 2t] == [1.0+t 1.0+2t; 2.0+t 2.0+2t]
    @test [1.0] .+ t == t .+ [1.0] == [1.0 + t]
    @test 1.0 .* t == t
    @test typeof(1.0 .* t) == Taylor1{Float64}

    st = sin(t)
    @test st .== st
    @test st == sin.(t)
    @test st.(pi/3) == evaluate(st, pi/3)
    @test st(pi/3) == evaluate.(st, pi/3)
    @test st.([0.0, pi/3]) == evaluate(st, [0.0, pi/3])

    @test typeof(Float32.(t)) == Taylor1{Float32}
    @test (Float32.(t))[1] == Float32(1.0)
    @test_throws MethodError Float32(t)

    # Nested Taylor1 tests
    t = Taylor1(Int, 3)
    ts = zero(t)
    ts .= t
    @test ts == t
    @. ts = 3 * t^2 - 1
    @test ts == 3 * t^2 - 1

    #  `tt` has to be `Taylor1{Taylor1{Float64}}` (instead of `Taylor1{Taylor1{Int}}`)
    # since the method a^n (n integer) is equivalent to `a^float(n).`
    tt = Taylor1([zero(1.0*t), one(t)], 2)
    tts = zero(tt)
    @test tt .== tt
    @. tts = 3 * tt^2 - 1
    @test tts == 3 * tt^2 - 1

    ttt = Taylor1([zero(tt), one(tt)])
    ttts = zero(ttt)
    @test ttt .≈ ttt
    @. ttts = 3 * ttt^1 - 1
    @test ttts == 3 * ttt^1 - 1
    @. ttts = 3 * ttt^3 - 1
    @test ttts == - 1.0
end

@testset "Broadcasting with HomogeneousPolynomial and TaylorN" begin
    x, y = set_variables("x y", order=3)
    xH = x[1]
    yH = y[1]

    @test xH .== xH
    @test yH .≈ yH
    @test xH .== xH
    @test x[2] .== y[2]

    xHs = zero(xH)
    xHs .= xH
    @test xHs == xH
    @. xHs = 2 * xH + yH
    @test xHs == 2 * xH + yH

    @test 1 .* xH == xH
    @test 1 .* [xH] == [xH]
    @test [1] .* xH == xH .* [1] == [xH]

    @test x .== x
    @test y .≈ y
    @test x .!= (1 + x)

    @test typeof(Float32.(x)) == TaylorN{Float32}
    @test (Float32.(x))[1] == HomogeneousPolynomial(Float32[1.0, 0.0])
    @test_throws MethodError Float32(x)

    p = zero(x)
    p .= x
    @test p == x
    @. p = 1 + 2*x + 3x^2 - x * y
    @test p == 1 + 2*x + 3*x^2 - x * y

    @test 1.0 .+ x == 1.0 + x
    @test y .+ x == y + x
    @test typeof(big"1.0" .+ x) == TaylorN{BigFloat}
    @test 1.0 .+ [x] == [1.0 + x]
    @test typeof(1.0 .+ [y]) == Vector{TaylorN{Float64}}
    @test 1.0 .+ [x, 2y] == [1.0 + x, 1.0 + 2y]
    @test [1.0,2.0] .+ [x 2y] == [1.0+x 1.0+2y; 2.0+x 2.0+2y]
    @test [1.0] .+ x == x .+ [1.0] == [1.0 + x]
    @test 1.0 .* y == y
    @test typeof(1.0 .* x .* y) == TaylorN{Float64}
end

@testset "Broadcasting with mixtures Taylor1{TalorN{T}}" begin
    x, y = set_variables("x", numvars=2, order=6)
    tN = Taylor1(TaylorN{Float64}, 3)

    @test tN .== tN
    @test tN .≈ tN
    @test tN .!= (1 + tN)

    @test 1.0 .+ tN == 1.0 + tN
    @test typeof(1.0 .+ tN) == Taylor1{TaylorN{Float64}}
    @test 1.0 .+ [tN] == [1.0 + tN]
    @test typeof(1.0 .+ [tN]) == Vector{Taylor1{TaylorN{Float64}}}
    @test 1.0 .+ [tN, 2tN] == [1.0 + tN, 1.0 + 2tN]
    @test [1.0,2.0] .+ [tN 2tN] == [1.0+tN 1.0+2tN; 2.0+tN 2.0+2tN]
    @test [1.0] .+ tN == tN .+ [1.0] == [1.0 + tN]
    @test 1.0 .* tN == 1.0 * tN
    @test typeof(1.0 .* tN) == Taylor1{TaylorN{Float64}}

    tNs = zero(tN)
    tNs .= tN
    @test tNs == tN
    @. tNs = y[1] * tN^2 - 1
    @test tNs == y[1] * tN^2 - 1
    @. tNs = y * tN^2 - 1
    @test tNs == y * tN^2 - 1
end

@testset "Broadcasting with mixtures TaylorN{Talor1{T}}" begin
    set_variables("x", numvars=2, order=6)
    t = Taylor1(3)
    xHt = HomogeneousPolynomial([one(t), zero(t)])
    yHt = HomogeneousPolynomial([zero(t), t])
    tN1 = TaylorN([HomogeneousPolynomial([t]),xHt,yHt^2])

    @test tN1 .== tN1
    @test tN1 .≈ tN1
    @test tN1 .!= (1 + tN1)

    @test 1.0 .+ tN1 == 1.0 + tN1
    @test typeof(1.0 .+ tN1) == TaylorN{Taylor1{Float64}}
    @test 1.0 .+ [tN1] == [1.0 + tN1]
    @test typeof(1.0 .+ [tN1]) == Vector{TaylorN{Taylor1{Float64}}}
    @test 1.0 .+ [tN1, 2tN1] == [1.0 + tN1, 1.0 + 2tN1]
    @test [1.0, 2.0] .+ [tN1 2tN1] == [1.0+tN1 1.0+2tN1; 2.0+tN1 2.0+2tN1]
    @test [1.0] .+ tN1 == tN1 .+ [1.0] == [1.0 + tN1]
    @test 1.0 .* tN1 == tN1
    @test typeof(1.0 .* tN1) == TaylorN{Taylor1{Float64}}

    tN1s = zero(tN1)
    tN1s .= tN1
    @test tN1s == tN1
    @. tN1s = t * tN1^2 - 1
    @test tN1s == t * tN1^2 - 1
end


text:
using Distributed
addprocs(2)

using Test
using Dagger

include("domain.jl")
include("array.jl")
include("scheduler.jl")
include("fault-tolerance.jl")
println(stderr, "tests done. cleaning up...")
Dagger.cleanup()
#include("cache.jl")
println(stderr, "all done.")


text:
#=
Similar to Batch Normalization, except online and without the rescaling/skew
    y = (a .- μ) ./ σ

    TODO: This is currently broken because OnlineStats.Variances no longer
    exists.
=#
type InputNorm{T,W<:Weight} <: Transformation
    n::Int  # maps n --> n
    vars::Variances{W}
    input::SumNode{T,1}    # a
    output::OutputNode{T,1}  # y
end

function InputNorm{T}(::Type{T}, n::Int,
                    lookback::Int = 100,
                    α::Float64 = NaN,
                    wgt::Weight =  BoundedEqualWeight(isnan(α) ? lookback : α)
                    )
    InputNorm(n,
        Variances(n, wgt),
        InputNode(T, n),
        OutputNode(T, n)
    )
end
InputNorm(n::Int, args...; kw...) = InputNorm(Float64, n, args...; kw...)

function Base.show(io::IO, t::InputNorm)
    print(io, "InputNorm{n=$(t.n)}")
end

function transform!{T}(layer::InputNorm{T})
    a = layer.input.val
    y = layer.output.val

    OnlineStats.fit!(layer.vars, a)
    μ = mean(layer.vars)
    σ = std(layer.vars)
    for i=1:layer.n
        y[i] = (a[i] - μ[i]) / max(σ[i], T(1e-10))
    end
    y
end

function grad!{T}(layer::InputNorm{T})
    ∇a = layer.input.∇
    ∇y = layer.output.∇

    σ = std(layer.vars)
    for i=1:layer.n
        ∇a[i] = ∇y[i] / max(σ[i], T(1e-10))
    end
    ∇a
end


text:
# # SIGNATURES

const DOC_SIGNATURES =
"""
A learning network *signature* is an intermediate object defined when
a user constructs a learning network machine, `mach`. They are named
tuples whose values are the nodes consitituting interface points
between the network and the machine.  Examples are

    (predict=yhat, )
    (transform=Xsmall,)
    (predict=yhat, transform=W, report=(loss=loss_node,))

where `yhat`, `Xsmall`, `W` and `loss_node` are nodes in the network.

If a key `k` is the name of an operation (such as `:predict`,
`:predict_mode`, `:transform`, `inverse_transform`) then `k(mach, X)`
returns `n(X)` where `n` is the corresponding node value.  Each such
node must have a unique origin (`length(origins(n)) === 1`).

The only other allowed key is `:report`, whose corresponding value
must be a named tuple

    (k1=n1, k2=n2, ...)

whose keys are arbitrary, and whose values are nodes of the
network. For each such key-value pair `k=n`, the value returned by
`n()` is included in the named tuple `report(mach)`, with
corresponding key `k`. So, in the third example above,
`report(mach).loss` will return the value of `loss_node()`.

"""

function _operation_part(signature)
    ops = filter(in(OPERATIONS), keys(signature))
    return NamedTuple{ops}(map(op->getproperty(signature, op), ops))
end
function _report_part(signature)
    :report in keys(signature) || return NamedTuple()
    return signature.report
end

_operations(signature) = keys(_operation_part(signature))

function _nodes(signature)
    return (values(_operation_part(signature))...,
            values(_report_part(signature))...)
end

function _call(nt::NamedTuple)
    _call(n) = deepcopy(n())
    _keys = keys(nt)
    _values = values(nt)
    return NamedTuple{_keys}(_call.(_values))
end

"""
    model_supertype(signature)

Return, if this can be inferred, which of `Deterministic`,
`Probabilistic` and `Unsupervised` is the appropriate supertype for a
composite model obtained by exporting a learning network with the
specified `signature`.

$DOC_SIGNATURES

If a supertype cannot be inferred, `nothing` is returned.

If the network with given `signature` is not exportable, this method
will not error but it will not a give meaningful return value either.

**Private method.**

"""
function model_supertype(signature)

    operations = _operations(signature)

    length(intersect(operations, (:predict_mean, :predict_median))) == 1 &&
        return Deterministic

    if :predict in operations
        node = signature.predict
        if node isa Source
            return Deterministic
        end
        if node.machine !== nothing
            model = node.machine.model
            model isa Deterministic && return Deterministic
            model isa Probabilistic && return Probabilistic
        end
    end

    return nothing

end

# # FITRESULTS FOR COMPOSITE MODELS

mutable struct CompositeFitresult
    signature
    glb
    report_additions
    function CompositeFitresult(signature)
        glb = MLJBase.glb(_nodes(signature)...)
        new(signature, glb)
    end
end
signature(c::CompositeFitresult) = getfield(c, :signature)
glb(c::CompositeFitresult) = getfield(c, :glb)
report_additions(c::CompositeFitresult) = getfield(c, :report_additions)

update!(c::CompositeFitresult) =
    setfield!(c, :report_additions, _call(_report_part(signature(c))))

# To accommodate pre-existing design (operations.jl) arrange
# that `fitresult.predict` returns the predict node, etc:
Base.propertynames(c::CompositeFitresult) = keys(signature(c))
Base.getproperty(c::CompositeFitresult, name::Symbol) =
    getproperty(signature(c), name)

# # LEARNING NETWORK MACHINES

surrogate(::Type{<:Deterministic})  = Deterministic()
surrogate(::Type{<:Probabilistic})  = Probabilistic()
surrogate(::Type{<:Unsupervised}) = Unsupervised()
surrogate(::Type{<:Static}) = Static()

caches_data_by_default(::Type{<:Surrogate}) = false

const ERR_MUST_PREDICT = ArgumentError(
    "You must specify at least `predict=<some node>`. ")
const ERR_MUST_TRANSFORM = ArgumentError(
    "You must specify at least `transform=<some node>`. ")
const ERR_MUST_OPERATE = ArgumentError(
    "You must specify at least one operation, as in `predict=<some node>`. ")
const ERR_MUST_SPECIFY_SOURCES = ArgumentError(
    "You must specify at least one source `Xs`, as in "*
    "`machine(surrogate_model, Xs, ...; kwargs...)`. ")
const ERR_BAD_SIGNATURE = ArgumentError(
    "Only the following keyword arguments are supported in learning network "*
    "machine constructors: `report` or one of: `$OPERATIONS`. ")
const ERR_EXPECTED_NODE_IN_SIGNATURE = ArgumentError(
    "Learning network machine constructor syntax error. "*
    "Did not enounter `Node` in place one was expected. ")

function check_surrogate_machine(::Surrogate, signature, _sources)
    isempty(_operations(signature)) && throw(ERR_MUST_OPERATE)
    isempty(_sources) && throw(ERR_MUST_SPECIFY_SOURCES)
    return nothing
end

function check_surrogate_machine(::Union{Supervised,SupervisedAnnotator},
                                 signature,
                                 _sources)
    isempty(_operations(signature)) && throw(ERR_MUST_PREDICT)
    length(_sources) > 1 || throw(err_supervised_nargs())
    return nothing
end

function check_surrogate_machine(::Union{Unsupervised},
                                 signature,
                                 _sources)
    isempty(_operations(signature)) && throw(ERR_MUST_TRANSFORM)
    length(_sources) < 2 || throw(err_unsupervised_nargs())
    return nothing
end

function machine(model::Surrogate, _sources::Source...; pair_itr...)

    # named tuple, such as `(predict=yhat, transform=W)`:
    signature = (; pair_itr...)

    # signature checks:
    isempty(_operations(signature)) && throw(ERR_MUST_OPERATE)
    for k in keys(signature)
        if k in OPERATIONS
            getproperty(signature, k) isa AbstractNode ||
                throw(ERR_EXPECTED_NODE_IN_SIGNATURE)
        elseif k === :report
            all(v->v isa AbstractNode, values(signature.report)) ||
                throw(ERR_EXPECTED_NODE_IN_SIGNATURE)
        else
            throw(ERR_BAD_SIGNATURE)
        end
    end

    check_surrogate_machine(model, signature, _sources)

    mach = Machine(model, _sources...)

    mach.fitresult = CompositeFitresult(signature)

    return mach

end

function machine(_sources::Source...; pair_itr...)

    signature = (; pair_itr...)

    T = model_supertype(signature)
    if T == nothing
        @warn "Unable to infer surrogate model type. \n"*
            "Using Deterministic(). To override specify "*
            "surrogate model, as in "*
        "`machine(Probabilistic(), ...)` or `machine(Interval(), ...)`"
        model = Deterministic()
    else
        model = surrogate(T)
    end

    return machine(model, _sources...; pair_itr...)

end

"""
    N = glb(mach::Machine{<:Union{Composite,Surrogate}})

A greatest lower bound for the nodes appearing in the signature of
`mach`.

$DOC_SIGNATURES

**Private method.**

"""
glb(mach::Machine{<:Union{Composite,Surrogate}}) = glb(mach.fitresult)

"""
    fit!(mach::Machine{<:Surrogate};
         rows=nothing,
         acceleration=CPU1(),
         verbosity=1,
         force=false))

Train the complete learning network wrapped by the machine `mach`.

More precisely, if `s` is the learning network signature used to
construct `mach`, then call `fit!(N)`, where `N` is a greatest lower
bound of the nodes appearing in the signature (values in the signature
that are not `AbstractNode` are ignored). For example, if `s =
(predict=yhat, transform=W)`, then call `fit!(glb(yhat, W))`.

See also [`machine`](@ref)

"""
function fit!(mach::Machine{<:Surrogate}; kwargs...)
    glb_node = glb(mach)
    fit!(glb_node; kwargs...)
    update!(mach.fitresult) # updates `report_additions`
    mach.state += 1
    mach.report = merge(report(glb_node), report_additions(mach.fitresult))
    return mach
end

MLJModelInterface.fitted_params(mach::Machine{<:Surrogate}) =
    fitted_params(glb(mach))

# # CONSTRUCTING THE RETURN VALUE FOR A COMPOSITE FIT METHOD

# Identify which properties of `model` have, as values, a model in the
# learning network wrapped by `mach`, and check that no two such
# properties have have identical values (#377). Return the property name
# associated with each model in the network (in the order appearing in
# `models(glb(mach))`) using `nothing` when the model is not
# associated with any property.
function network_model_names(model::M,
                             mach::Machine{<:Surrogate}) where M<:Model

    network_model_ids = objectid.(MLJBase.models(glb(mach)))

    names = propertynames(model)

    # intialize dict to detect duplicity a la #377:
    name_given_id = Dict{UInt64,Vector{Symbol}}()

    # identify location of properties whose values are models in the
    # learning network, and build name_given_id:
    for name in names
        id = objectid(getproperty(model, name))
        if id in network_model_ids
            if haskey(name_given_id, id)
                push!(name_given_id[id], name)
            else
                name_given_id[id] = [name,]
            end
        end
    end

    # perform #377 check:
    no_duplicates = all(values(name_given_id)) do name
        length(name) == 1
    end
    if !no_duplicates
        for (id, name) in name_given_id
            if length(name) > 1
                @error "The hyperparameters $name of "*
                    "$model have identical model "*
                "instances as values. "
            end
        end
        throw(ArgumentError(
        "Two distinct hyper-parameters of a "*
            "composite model that are both "*
            "associated with models in the underlying learning "*
            "network (eg, any two components of a `@pipeline` model) "*
            "cannot have identical values, although they can be `==` "*
            "(corresponding nested properties are `==`). "*
            "Consider constructing instances "*
            "separately or use `deepcopy`. "))
    end

    return map(network_model_ids) do id
        if id in keys(name_given_id)
            return name_given_id[id] |> first
        else
            return nothing
        end
    end

end

"""

    return!(mach::Machine{<:Surrogate}, model, verbosity)

The last call in custom code defining the `MLJBase.fit` method for a
new composite model type. Here `model` is the instance of the new type
appearing in the `MLJBase.fit` signature, while `mach` is a learning
network machine constructed using `model`. Not relevant when defining
composite models using `@pipeline` or `@from_network`.

For usage, see the example given below. Specifically, the call does
the following:

- Determines which hyper-parameters of `model` point to model
  instances in the learning network wrapped by `mach`, for recording
  in an object called `cache`, for passing onto the MLJ logic that
  handles smart updating (namely, an `MLJBase.update` fallback for
  composite models).

- Calls `fit!(mach, verbosity=verbosity)`.

- Moves any data in source nodes of the learning network into `cache`
  (for data-anonymization purposes).

- Records a copy of `model` in `cache`.

- Returns `cache` and outcomes of training in an appropriate form
  (specifically, `(mach.fitresult, cache, mach.report)`; see [Adding
  Models for General
  Use](https://alan-turing-institute.github.io/MLJ.jl/dev/adding_models_for_general_use/)
  for technical details.)

### Example

The following code defines, "by hand", a new model type `MyComposite`
for composing standardization (whitening) with a deterministic
regressor:

```
mutable struct MyComposite <: DeterministicComposite
    regressor
end

function MLJBase.fit(model::MyComposite, verbosity, X, y)
    Xs = source(X)
    ys = source(y)

    mach1 = machine(Standardizer(), Xs)
    Xwhite = transform(mach1, Xs)

    mach2 = machine(model.regressor, Xwhite, ys)
    yhat = predict(mach2, Xwhite)

    mach = machine(Deterministic(), Xs, ys; predict=yhat)
    return!(mach, model, verbosity)
end
```

"""
function return!(mach::Machine{<:Surrogate},
                 model::Union{Model,Nothing},
                 verbosity)

    _network_model_names = network_model_names(model, mach)

    verbosity isa Nothing || fit!(mach, verbosity=verbosity)

    # anonymize the data:
    sources = mach.args
    data = Tuple(s.data for s in sources)
    [MLJBase.rebind!(s, nothing) for s in sources]

    # record the current hyper-parameter values:
    old_model = deepcopy(model)

    cache = (sources = sources,
             data=data,
             network_model_names=_network_model_names,
             old_model=old_model)

    return mach.fitresult, cache, mach.report

end

network_model_names(model::Nothing, mach::Machine{<:Surrogate}) =
    nothing

## DUPLICATING AND REPLACING PARTS OF A LEARNING NETWORK MACHINE

"""
    replace(mach, a1=>b1, a2=>b2, ...; empty_unspecified_sources=false)

Create a deep copy of a learning network machine `mach` but replacing
any specified sources and models `a1, a2, ...` of the original
underlying network with `b1, b2, ...`.

If `empty_unspecified_sources=true` then any source nodes not
specified are replaced with empty source nodes.

"""
function Base.replace(mach::Machine{<:Surrogate},
                      pairs::Pair...; empty_unspecified_sources=false)

    signature = MLJBase.signature(mach.fitresult)
    operation_nodes = values(_operation_part(signature))
    report_nodes = values(_report_part(signature))

    W = glb(operation_nodes..., report_nodes...)

    # Note: We construct nodes of the new network as values of a
    # dictionary keyed on the nodes of the old network. Additionally,
    # there are dictionaries of models keyed on old models and
    # machines keyed on old machines. The node and machine
    # dictionaries must be built simultaneously.

    # build model dict:
    model_pairs = filter(collect(pairs)) do pair
        first(pair) isa Model
    end
    models_ = models(W)
    models_to_copy = setdiff(models_, first.(model_pairs))
    model_copy_pairs = [model=>deepcopy(model) for model in models_to_copy]
    newmodel_given_old = IdDict(vcat(model_pairs, model_copy_pairs))

    # build complete source replacement pairs:
    sources_ = mach.args # sources(W)
    specified_source_pairs = filter(collect(pairs)) do pair
        first(pair) isa Source
    end
    unspecified_sources = setdiff(sources_, first.(specified_source_pairs))
    unspecified_sources_wrapping_something =
        filter(s -> !isempty(s), unspecified_sources)
    if !isempty(unspecified_sources_wrapping_something) &&
        !empty_unspecified_sources
        @warn "No replacement specified for one or more non-empty source "*
        "nodes. Contents will be duplicated. "
    end
    if empty_unspecified_sources
        unspecified_source_pairs = [s => source() for
                                    s in unspecified_sources]
    else
        unspecified_source_pairs = [s => deepcopy(s) for
                                    s in unspecified_sources]
    end

    all_source_pairs = vcat(specified_source_pairs, unspecified_source_pairs)

    # drop source nodes from all nodes of network terminating at W:
    nodes_ = filter(nodes(W)) do N
        !(N isa Source)
    end
    isempty(nodes_) && error("All nodes in network are source nodes. ")
    # instantiate node and machine dictionaries:
    newnode_given_old =
        IdDict{AbstractNode,AbstractNode}(all_source_pairs)
    newsources = [newnode_given_old[s] for s in sources_]
    newoperation_node_given_old =
        IdDict{AbstractNode,AbstractNode}()
    newreport_node_given_old =
        IdDict{AbstractNode,AbstractNode}()
    newmach_given_old = IdDict{Machine,Machine}()

    # build the new network:
    for N in nodes_
       args = [newnode_given_old[arg] for arg in N.args]
         if N.machine === nothing
             newnode_given_old[N] = node(N.operation, args...)
         else
             if N.machine in keys(newmach_given_old)
                 m = newmach_given_old[N.machine]
             else
                 train_args = [newnode_given_old[arg] for arg in N.machine.args]
                 m = Machine(newmodel_given_old[N.machine.model],
                                train_args...)
                 newmach_given_old[N.machine] = m
             end
             newnode_given_old[N] = N.operation(m, args...)
         end
        if N in operation_nodes
            newoperation_node_given_old[N] = newnode_given_old[N]
        elseif N in report_nodes
            newreport_node_given_old[N] = newnode_given_old[N]
        end
    end

    newoperation_nodes = Tuple(newoperation_node_given_old[N] for N in
                          operation_nodes)
    newreport_nodes = Tuple(newreport_node_given_old[N] for N in
                            report_nodes)
    report_tuple =
        NamedTuple{keys(_report_part(signature))}(newreport_nodes)
    operation_tuple =
        NamedTuple{keys(_operation_part(signature))}(newoperation_nodes)

    newsignature = if isempty(report_tuple)
        operation_tuple
    else
        merge(operation_tuple, (report=report_tuple,))
    end

    return machine(mach.model, newsources...; newsignature...)

end