text:
using Isoband
using Test

# let's use this quadruple bar thing as the replacement for setequal() from the R tests for convenience
≣(s1::Set, s2::Set) = s1 == s2
≣(s1, s2) = Set(s1) == Set(s2)

@testset "elementary polygons get merged" begin
    # two connected polygons get merged
    z = [0 0 1;
         1 1 1]

    out = isobands(1:3, 2:-1:1, z, 0.5, 1.5)

    @test 1000 * out.x + out.y ≣ 1000 * [3.0, 2.0, 1.0, 1.0, 2.0, 2.5, 3.0] + [1.0, 1.0, 1.0, 1.5, 1.5, 2.0, 2.0]
    @test out.id == ones(7)

    # two unconnected polygons don't get merged
    z = [1 2 1;
         1 2 2]
    out = isobands(1:3, 2:-1:1, z, 0.5, 1.5)
    @test 1000 * out.x + out.y ≣ 1000 * [3.0, 2.5, 3.0, 1.0, 1.5, 1.5, 1.0] + [1.5, 2.0, 2.0, 2.0, 2.0, 1.0, 1.0]
    @test out.id ≣ [1, 2]
    @test length(out.id) == 7

    # two separate bands get merged in second row
    z = [1 2 1;
         1 2 1;
         0 0 0]

    out = isobands(1:3, 3:-1:1, z, 0.5, 1.5)
    @test 1000 * out.x + out.y ≣ 1000 * [3.0, 2.0, 1.0, 1.0, 1.0, 1.5, 1.5, 2.0, 2.5, 2.5, 3.0, 3.0] + [1.50, 1.25, 1.50, 2.00, 3.00, 3.00, 2.00, 1.75, 2.00, 3.00, 3.00, 2.00]
    @test out.id == ones(12)

    # circle gets closed
    z = [1 1 1;
         1 2 1;
         1 1 1]

    out = isobands(1:3, 3:-1:1, z, 0.5, 1.5)
    @test 1000 * out.x + out.y ≣ 1000 * [3.0, 2.0, 1.0, 1.0, 1.0, 2.0, 3.0, 3.0, 2.0, 2.5, 2.0, 1.5] + [1.0, 1.0, 1.0, 2.0, 3.0, 3.0, 3.0, 2.0, 1.5, 2.0, 2.5, 2.0]
    @test out.id ≣ [1, 2]
    @test length(out.id) == 12
end

@testset "NaNs are handled correctly" begin
    z = [NaN 1 1;
         1 1 1;
         1 1 1]

    out = isobands(1:3, 3:-1:1, z, 0.5, 1.5)
    @test 1000 * out.x + out.y ≣ 1000 * [3.0, 2, 1, 1, 2, 2, 3, 3] + [1.0, 1, 1, 2, 2, 3, 3, 2]
    @test out.id == ones(8)

    z = [NaN 1 1;
         1 1 1;
         1 1 NaN]

    out = isobands(1:3, 3:-1:1, z, 0.5, 1.5)
    @test 1000 * out.x + out.y ≣ 1000 * [1.0, 1, 2, 2, 2, 2, 3, 3] + [1, 2, 2, 1, 2, 3, 3, 2]
    @test out.id ≣ [1, 2]
    @test length(out.id) == 8
end

@testset "All elementary shapes are calculated correctly" begin
  # a matrix that requires all elementary shapes for isobanding
  z =  [0 0 0 1 0 2 1 0 1 1 1 2 2 0 2 1 2 2;
        0 0 0 1 0 2 1 0 1 1 1 2 2 0 2 1 2 2;
        0 1 0 2 1 0 1 1 1 2 2 0 2 1 2 2 0 0;
        1 0 1 1 1 2 2 0 2 1 2 2 0 0 0 1 0 2;
        2 0 2 1 2 2 0 0 0 1 0 2 1 0 1 1 1 2;
        0 1 0 2 1 0 1 1 1 2 2 0 2 1 2 2 0 0;
        2 0 2 1 2 2 0 0 0 1 0 2 1 0 1 1 1 2;
        1 0 1 1 1 2 2 0 2 1 2 2 0 0 0 1 0 2;
        0 1 0 2 1 0 1 1 1 2 2 0 2 1 2 2 0 0;
        0 0 0 1 0 2 1 0 1 1 1 2 2 0 2 1 2 2;
        0 0 0 1 0 2 1 0 1 1 1 2 2 0 2 1 2 2]

    out = isobands(1:18, 11:-1:1, z, 0.5, 1.5)

    @test 1000 * out.x + out.y ≣ 1000 * [
        16.00, 15.50, 15.50, 16.00, 16.25, 16.00, 15.00, 14.50, 14.75, 14.75, 14.25, 14.25,
        14.00, 13.75, 13.75, 13.25, 13.25, 13.50, 13.00, 12.75, 12.00, 11.50, 11.50, 11.00,
        10.00,  9.00,  8.50,  8.50,  8.00,  7.50,  7.50,  7.00,  6.50,  6.50,  6.00,  5.75,
        5.75,  5.25,  5.25,  5.00,  4.50,  4.50,  4.00,  3.50,  3.50,  3.25,  3.00,  2.50,
        2.00,  1.50,  1.00,  1.00,  1.00,  1.25,  1.00,  1.00,  1.50,  1.00,  1.00,  1.25,
        1.00,  1.00,  1.00,  1.50,  2.00,  2.50,  3.00,  3.25,  3.50,  3.50,  4.00,  4.50,
        4.50,  5.00,  5.25,  5.25,  5.75,  5.75,  6.00,  6.50,  6.50,  7.00,  7.50,  7.50,
        8.00,  8.50,  8.50,  9.00, 10.00, 11.00, 11.50, 11.50, 12.00, 12.75, 13.00, 13.50,
        13.25, 13.25, 13.75, 13.75, 14.00, 14.25, 14.25, 14.75, 14.75, 14.50, 15.00, 16.00,
        16.25, 16.00, 15.50, 15.50, 16.00, 16.50, 16.50, 17.00, 18.00, 18.00, 17.00, 16.75,
        16.50, 17.00, 17.25, 18.00, 18.00, 17.75, 17.50, 18.00, 18.00, 17.00, 16.75, 17.00,
        18.00, 18.00, 17.50, 17.75, 18.00, 18.00, 17.25, 17.00, 16.50, 16.75, 17.00, 18.00,
        18.00, 17.00, 16.50, 16.50, 12.00, 12.25, 12.00, 11.75, 11.00, 11.25, 12.00, 12.25,
        12.50, 12.00, 11.75, 11.00, 10.50, 10.00,  9.50,  9.00,  8.75,  9.00,  9.50, 10.00,
        6.00,  6.50,  6.00,  5.50,  4.00,  4.50,  4.00,  3.75, 13.00, 14.00, 15.00, 15.50,
        15.00, 14.50, 14.00, 13.50, 13.00, 12.75,  8.00,  8.25,  9.00,  9.50,  9.00,  8.00,
        7.00,  6.75,  7.00,  7.75,  7.00,  7.25,  7.00,  6.25,  6.00,  5.00,  4.50,  5.00,
        6.00,  6.25,  7.00,  7.25,  7.00,  6.00,  5.50,  5.00,  4.50,  4.00,  3.75,  4.00,
        4.50,  5.00,  5.50,  6.00,  2.00,  2.50,  2.25,  2.00,  1.75,  1.50, 11.00, 11.25,
        11.00, 10.50, 13.00, 13.50, 13.00, 12.75, 12.00, 12.25, 12.00, 11.75, 11.00, 11.25,
        11.00, 10.00,  9.50, 10.00,  9.00,  8.25,  8.00,  7.75,  7.00,  6.75,  7.00,  8.00,
        9.00,  9.50, 16.25, 16.00, 15.00, 14.50, 15.00, 16.00, 15.00, 14.00, 13.00, 12.75,
        13.00, 13.50, 14.00, 14.50, 15.00, 15.50, 12.00, 12.50, 12.25, 12.00, 11.25, 11.00,
        10.00,  9.50,  9.00,  8.75,  9.00,  9.50, 10.00, 10.50, 11.00, 11.75, 11.00, 11.25,
        11.00, 10.50,  6.00,  5.50,  6.00,  6.50,  3.00,  2.50,  3.00,  3.25,  3.00,  3.50,
        3.00,  2.75,  2.00,  2.25,  2.50,  2.00,  1.50,  1.75,  3.00,  2.75,  3.00,  3.50,
        4.00,  3.75,  4.00,  4.50, 12.00, 12.25, 12.00, 11.75,  6.00,  6.50,  6.00,  5.50] +

        [
        1.00,  1.00,  2.00,  2.50,  3.00,  3.50,  3.25,  3.00,  2.00,  1.00,  1.00,  2.00,
        2.50,  2.00,  1.00,  1.00,  2.00,  3.00,  3.25,  3.00,  2.25,  2.00,  1.00,  1.00,
        1.00,  1.00,  1.00,  2.00,  2.50,  2.00,  1.00,  1.00,  1.00,  2.00,  2.25,  2.00,
        1.00,  1.00,  2.00,  2.50,  2.00,  1.00,  1.00,  1.00,  2.00,  3.00,  3.50,  3.00,
        2.50,  3.00,  3.50,  4.00,  4.50,  5.00,  5.25,  5.75,  6.00,  6.25,  6.75,  7.00,
        7.50,  8.00,  8.50,  9.00,  9.50,  9.00,  8.50,  9.00, 10.00, 11.00, 11.00, 11.00,
       10.00,  9.50, 10.00, 11.00, 11.00, 10.00,  9.75, 10.00, 11.00, 11.00, 11.00, 10.00,
        9.50, 10.00, 11.00, 11.00, 11.00, 11.00, 11.00, 10.00,  9.75,  9.00,  8.75,  9.00,
       10.00, 11.00, 11.00, 10.00,  9.50, 10.00, 11.00, 11.00, 10.00,  9.00,  8.75,  8.50,
        9.00,  9.50, 10.00, 11.00, 11.00, 11.00, 10.00,  9.75,  9.75,  9.25,  9.25,  9.00,
        8.00,  7.50,  8.00,  8.75,  8.25,  8.00,  7.00,  6.75,  6.25,  6.50,  6.00,  5.50,
        5.75,  5.25,  5.00,  4.00,  3.75,  3.25,  4.00,  4.50,  4.00,  3.00,  2.75,  2.75,
        2.25,  2.25,  2.00,  1.00,  2.75,  3.00,  3.25,  3.00,  2.50,  3.00,  3.75,  4.00,
        5.00,  5.25,  5.00,  4.25,  4.00,  3.50,  4.00,  4.25,  4.00,  3.50,  3.00,  2.50,
        2.75,  3.00,  3.25,  3.00,  2.50,  3.00,  3.50,  3.00,  3.75,  3.50,  3.75,  4.00,
        4.50,  5.00,  5.50,  5.00,  4.50,  4.00,  3.50,  4.00,  4.75,  5.00,  5.50,  5.50,
        5.50,  5.00,  4.75,  4.00,  3.50,  4.00,  4.25,  5.00,  5.25,  5.50,  6.00,  6.50,
        6.75,  7.00,  7.75,  8.00,  8.50,  8.25,  8.00,  7.50,  7.00,  6.50,  6.00,  5.50,
        5.00,  4.50,  4.00,  3.75,  3.50,  4.00,  5.00,  5.50,  5.00,  4.00,  4.75,  5.00,
        5.25,  5.00,  5.50,  6.00,  6.50,  6.00,  5.75,  6.00,  6.25,  6.00,  5.75,  6.00,
        6.25,  6.50,  6.00,  5.50,  7.25,  8.00,  8.50,  8.00,  7.25,  7.00,  6.50,  6.50,
        6.50,  7.00,  6.00,  6.50,  6.50,  6.00,  5.50,  5.50,  8.25,  8.50,  8.25,  8.00,
        7.50,  7.00,  6.50,  7.00,  7.50,  8.00,  6.75,  7.00,  8.00,  8.25,  9.00,  9.50,
        9.50,  9.00,  8.50,  8.00,  7.75,  8.00,  8.50,  8.00,  7.75,  7.00,  6.75,  7.00,
        7.25,  7.00,  6.25,  6.00,  5.75,  6.00,  6.25,  6.00,  5.75,  6.00,  6.75,  7.00,
        7.50,  7.00,  6.50,  7.00,  8.00,  8.50,  8.00,  7.00,  5.25,  5.00,  4.50,  5.00,
        9.50,  9.00,  8.50,  9.00,  8.75,  9.00,  9.25,  9.00,  8.75,  9.00,  9.25,  9.00
     ]

  @test out.id ≣ 1:26
  @test length(out.id) == 324
end

@testset "Six-sided saddles" begin
    # a matrix that contains all six-sided saddles
    z = [0 1 1 2;
        1 0 2 1;
        0 1 1 2]

    # midpoint outside the band
    out = isobands(1:4, 3:-1:1, z, 0.6, 1.4)
    @test 1000 * out.x + out.y ≣ 1000 * [3.0, 2.0, 1.6, 2.0, 2.3, 2.0,
        1.6, 2.0, 3.0, 3.4, 3.0, 2.7,
        3.0, 3.4, 1.0, 1.0, 1.0, 1.4,
        4.0, 4.0, 4.0, 3.6] +
     [1.0, 1.0, 1.0, 1.4, 2.0, 2.6,
        3.0, 3.0, 3.0, 3.0, 2.6, 2.0,
        1.4, 1.0, 1.6, 2.0, 2.4, 2.0,
        2.4, 2.0, 1.6, 2.0]
    @test out.id ≣ 1:3
    @test length(out.id) == 22

    # midpoint inside the band
    out = isobands(1:4, 3:-1:1, z, 0.4, 1.6)
    @test 1000 * out.x + out.y ≣ 1000 * [3.0, 2.0, 1.4, 1.0, 1.0, 1.0,
        1.4, 2.0, 3.0, 3.6, 4.0, 4.0,
        4.0, 3.6, 3.0, 3.4, 3.0, 2.8,
        2.0, 2.2, 2.0, 1.6] +
    [1.0, 1.0, 1.0, 1.4, 2.0, 2.6,
        3.0, 3.0, 3.0, 3.0, 2.6, 2.0,
        1.4, 1.0, 1.6, 2.0, 2.4, 2.0,
        1.6, 2.0, 2.4, 2.0]
    @test out.id ≣ 1:3
    @test length(out.id) == 22
end

@testset "Seven-sided saddles" begin
    # a matrix that contains all seven-sided saddles
    z = [0 1 0 1 2 1;
        2 0 2 2 0 2;
        0 1 0 1 2 1]

    # midpoint inside the band
    out = isobands(1:6, 3:-1:1, z, 0.5, 1.5)
    @test 1000 * out.x + out.y ≣ 1000 * [6.00, 6.00, 5.50, 5.00, 4.50, 4.00,
        3.50, 3.00, 2.50, 2.00, 1.50, 1.00,
        1.00, 1.25, 1.00, 1.00, 1.50, 2.00,
        2.50, 3.00, 3.50, 4.00, 4.50, 5.00,
        5.50, 6.00, 6.00, 5.75, 4.00, 4.25,
        4.00, 3.00, 2.75, 3.00, 2.00, 2.25,
        2.00, 1.75, 5.00, 5.25, 5.00, 4.75] +
    [1.50, 1.00, 1.00, 1.25, 1.00, 1.00,
        1.00, 1.25, 1.00, 1.00, 1.00, 1.25,
        1.75, 2.00, 2.25, 2.75, 3.00, 3.00,
        3.00, 2.75, 3.00, 3.00, 3.00, 2.75,
        3.00, 3.00, 2.50, 2.00, 1.50, 2.00,
        2.50, 2.25, 2.00, 1.75, 1.50, 2.00,
        2.50, 2.00, 1.75, 2.00, 2.25, 2.00]

    @test out.id ≣ 1:4
    @test length(out.id) == 42

    # midpoint outside the band
    out = isobands(1:6, 3:-1:1, z, 0.8, 1.2)
    @test 1000 * out.x + out.y ≣ 1000 * [6.0, 6.0, 5.8, 5.0, 4.4, 5.0,
        5.6, 4.2, 4.0, 3.8, 3.0, 2.4,
        3.0, 3.8, 4.0, 4.2, 4.0, 3.0,
        2.6, 3.0, 4.0, 2.2, 2.0, 1.8,
        2.0, 1.0, 1.4, 1.0, 1.0, 1.6,
        1.0, 5.0, 5.4, 5.0, 4.6, 6.0,
        5.8, 6.0, 2.2, 2.0, 1.8, 2.0] +
    [1.2, 1.0, 1.0, 1.4, 2.0, 2.6,
        2.0, 1.0, 1.0, 1.0, 1.4, 2.0,
        2.6, 3.0, 3.0, 3.0, 2.8, 2.4,
        2.0, 1.6, 1.2, 1.0, 1.0, 1.0,
        1.2, 1.6, 2.0, 2.4, 2.6, 2.0,
        1.4, 1.6, 2.0, 2.4, 2.0, 2.8,
        3.0, 3.0, 3.0, 2.8, 3.0, 3.0]
    @test out.id ≣ 1:8
    @test length(out.id) == 42
end

@testset "Eight-sided saddles" begin
    # a matrix that contains all eight-sided saddles
    z = [0 2 0;
         2 0 2]

    # midpoint above the band
    out = isobands(1:3, 2:-1:1, z, 0.5, 0.8)
    @test 1000 * out.x + out.y ≣ 1000 * [3.00, 3.00, 2.60, 2.75, 2.00,
        2.40, 2.25, 2.00, 1.75, 1.60,
        1.00, 1.00, 1.25, 1.40] +
    [1.75, 1.60, 2.00, 2.00, 1.40,
        1.00, 1.00, 1.25, 1.00, 1.00,
        1.60, 1.75, 2.00, 2.00]
    @test out.id ≣ 1:3
    @test length(out.id) == 14

    # midpoint inside the band
    out = isobands(1:3, 2:-1:1, z, 0.5, 1.5)
    @test 1000 * out.x + out.y ≣ 1000 * [2.75, 2.25, 2.00, 1.75, 1.25,
        1.00, 1.00, 1.25, 1.75, 2.00,
        2.25, 2.75, 3.00, 3.00] +
    [1.00, 1.00, 1.25, 1.00, 1.00,
        1.25, 1.75, 2.00, 2.00, 1.75,
        2.00, 2.00, 1.75, 1.25]
    @test out.id == ones(14)

    # midpoint below the band
    out = isobands(1:3, 2:-1:1, z, 1.2, 1.5)
    @test 1000 * out.x + out.y ≣ 1000 * [3.00, 3.00, 2.75, 2.60, 2.40,
        2.00, 1.60, 1.75, 2.00, 2.25,
        1.40, 1.25, 1.00, 1.00] +
    [1.40, 1.25, 1.00, 1.00, 2.00,
        1.60, 2.00, 2.00, 1.75, 2.00,
        1.00, 1.00, 1.25, 1.40]
    @test out.id ≣ 1:3
    @test length(out.id) == 14
end

@testset "Inconsistent numbers of isoband levels cause an error" begin
  m = [0 0 1 1;
        0 1 1 1;
        1 1 0 0;
        0 0 1 0]

  # error, multiple values are not recycled
  @test_throws ErrorException isobands(1:4, 1:4, m, [0.5, 1.5, 2.5], [0.5, 1.5])
  @test_throws ErrorException isobands(1:4, 1:4, m, [0.5, 1.5], [0.5, 1.5, 2.5])
end

#######################################################
# isolines

@testset "line segments get merged" begin
    # two connected line segments get merged
    z = [0 0 1;
        1 1 1]
    out = isolines(1:3, 2:-1:1, z, 0.5)
    @test 1000 * out.x + out.y ≣ 1000 * [1, 2, 2.5] + [1.5, 1.5, 2.0]
    @test out.id == ones(3)

    # two unconnected line segments don't get merged
    z = [0 1 0;
        0 1 1]
    out = isolines(1:3, 2:-1:1, z, 0.5)
    @test 1000 * out.x + out.y ≣ 1000 * [2.5, 3.0, 1.5, 1.5] + [2.0, 1.5, 2.0, 1.0]
    @test out.id ≣ 1:2
    @test length(out.id) == 4

    # two separate lines get merged in second row
    z = [0 1 0;
        0 1 0;
        0 0 0]
    out = isolines(1:3, 3:-1:1, z, 0.5)
    @test 1000 * out.x + out.y ≣ 1000 * [2.5, 2.5, 2.0, 1.5, 1.5] + [3.0, 2.0, 1.5, 2.0, 3.0]
    @test out.id == ones(5)

    # circle gets closed
    z = [0 0 0;
        0 1 0;
        0 0 0]
    out = isolines(1:3, 3:-1:1, z, 0.5)
    # circle is closed
    @test out.x[1] == out.x[5]
    @test out.y[1] == out.y[5]
    # coords are correct
    @test 1000 * out.x + out.y ≣ 1000 * [2.5, 2.0, 1.5, 2.0, 2.5] + [2.0, 2.5, 2.0, 1.5, 2.0]
    @test out.id == ones(5)
end

@testset "NAs are handled correctly" begin
    z = [NaN 0 0;
        0 1 1;
        0 1 1]
    out = isolines(1:3, 3:-1:1, z, 0.5)
    @test 1000 * out.x + out.y ≣ 1000 * [1.5, 1.5, 2.0, 3.0] + [2.0, 1.0, 2.5, 2.5]
    @test out.id ≣ 1:2
    @test length(out.id) == 4
end

@testset "All elementary segments are calculated correctly" begin
    # a matrix that requires all elementary segments for isolines
    z = [0 0 0 1 1 0 1 1;
        0 0 0 1 1 0 1 1;
        0 1 1 0 1 1 0 0;
        1 1 0 0 0 1 1 0;
        1 0 1 1 0 0 0 1]
    out = isolines(1:8, 5:-1:1, z, 0.5)

    @test 1000 * out.x + out.y ≣ 1000 * [7.5, 7.0, 6.0, 5.5, 5.0, 4.5, 4.0, 3.5,
            3.0, 2.5, 3.0, 4.0, 4.5, 2.5, 2.0, 1.5,
            8.0, 7.0, 6.5, 7.0, 7.5, 8.0, 6.5, 6.5,
            6.0, 5.5, 5.5, 3.5, 3.5, 3.0, 2.0, 1.5,
            1.0] +
            [1.0, 1.5, 1.5, 2.0, 2.5, 3.0, 3.5, 3.0,
            2.5, 2.0, 1.5, 1.5, 1.0, 1.0, 1.5, 1.0,
            3.5, 3.5, 3.0, 2.5, 2.0, 1.5, 5.0, 4.0,
            3.5, 4.0, 5.0, 5.0, 4.0, 3.5, 3.5, 3.0,
            2.5]

    @test out.id ≣ 1:5
    @test length(out.id) == 33
end

@testset "Saddles" begin
  # a matrix that contains all saddles (there are only two)
  z = [0 1 0;
       1 0 1]

  out = isolines(1:3, 2:-1:1, z, 0.5)

  @test 1000 * out.x + out.y ≣ 1000 * [2.5, 3.0, 2.5, 2.0, 1.5, 1.5, 1.0] + [2.0, 1.5, 1.0, 1.5, 1.0, 2.0, 1.5]
  @test out.id ≣ 1:3
  @test length(out.id) == 7

  out = isolines(1:3, 2:-1:1, z, 0.6)

  @test 1000 * out.x + out.y ≣ 1000 * [1.6, 2.0, 2.4, 3.0, 2.6, 1.0, 1.4] + [2.0, 1.6, 2.0, 1.4, 1.0, 1.4, 1.0]
  @test out.id ≣ 1:3
  @test length(out.id) == 7
end


text:
using Pkg
Pkg.activate(".")
Pkg.instantiate()

using IJulia
notebook(dir=pwd())

ext: jl


text:
using Flux, Flux.Tracker
using Flux.@epochs
using PyPlot
using CSV

pygui(false)

debugging = false
raw_data = CSV.read("logs/mturk_trusted_condition_1_mod.csv")

# set random seed for repeatability
#  srand(12345)
# 4 is good
# 21 is good
# 39 is beter, with 0.3 dropout and 10 hidden nodes
srand_num = 39
srand(srand_num)

rand_order = randperm(size(raw_data,1))
x1 = [float(q) for q in raw_data[:tprob][rand_order]]
x2 = [float(q) for q in raw_data[:N][rand_order]]
y = [float(q) for q in raw_data[:X3_1][rand_order]]

nx = length(x1)
train_pct = 0.90
train_limit = convert(Int,round(nx*train_pct))
train_set = 1:train_limit

valid_pct = 1.00
valid_limit = convert(Int,round(nx*valid_pct))
valid_set = train_limit+1:valid_limit

#  test_set = valid_limit+1:length(x1)

x1_train = x1[train_set]
x2_train = x2[train_set]
y_train = y[train_set]'

x_train = [x1_train x2_train]'
#  x_train = x1_train'

data = zip(x_train,y_train)

x1_valid = x1[valid_set]
x2_valid = x2[valid_set]
y_valid = y[valid_set]'

x_valid = [x1_valid x2_valid]'
#  x_valid = x1_valid'

#  x1_test = x1[test_set]
#  x2_test = x2[test_set]
#  y_test = y[test_set]'
#
#  x_test = [x1_test x2_test]'

num_hidden = 10
dropout_rate = 0.30

net = Chain(
            Dense(size(x_train,1),num_hidden,relu),
            Dropout(dropout_rate),
            Dense(num_hidden,num_hidden,relu),
            Dropout(dropout_rate),
            Dense(num_hidden,1)
           )

display(params(net))

loss(x,y) = Flux.mse(net(x),y)

opt = ADAM(params(net) , 0.001; β1 = 0.9, β2 = 0.999, ϵ = 1e-08, decay = 0)

evalcb() = @show(loss(tx,ty))

# array to store the validation error per batch
res = []
tr = []
tst = []

function check_err(x,y,ary)
    Flux.testmode!(net)
    append!(ary,loss(x,y).data)
    Flux.testmode!(net,false)
end

function valid_cb()
    # validate data
    check_err(x_valid,y_valid,res)
end
function train_cb()
    # validate data
    check_err(x_train,y_train,tr)
end
function test_cb()
    # validate data
    check_err(x_test,y_test,tst)
end

function cb_fun()
    valid_cb()
    train_cb()
    #  test_cb()
end

num_epocs = 50

@epochs num_epocs Flux.train!(loss,data,opt,cb = cb_fun)
#  Flux.train!(loss,data,opt,cb = Flux.throttle(evalcb(),3))
#  Flux.train!(loss,data,opt,cb = ()=>println("test"))

fig,ax = PyPlot.subplots(1,1)
fig[:set_size_inches](6,4)

epoc_mark = collect(0+1/length(y_train):1/length(y_train):num_epocs)
#  println(length(res))
#  println(length(epoc_mark))
#  display(epoc_mark)
#  readline()

ax[:plot](epoc_mark,res,color=:black,label="validation (N=$(length(y_valid)))")
ax[:plot](epoc_mark,tr,color=:blue,label="training (N=$(length(y_train)))")
#  ax[:plot](epoc_mark,tst,color=:orange,label="test (N=$(length(y_test)))")
if debugging
    ax[:set_title]("Num Hidden: $num_hidden, its: $num_epocs, srand: $srand_num, dropout: $dropout_rate")
else
    ax[:set_title]("Loss vs Epoch")
end

ax[:set_xlabel]("epoc")
ax[:set_ylabel]("loss (MSE)")
ax[:ticklabel_format](axis="y",style="sci",scilimits=(0,1),useOffset=false)
ax[:legend]()

#  Flux.testmode!(net)
#  test_err = loss(x_test,y_test)
#  Flux.testmode!(net,false)

#  ax[:axhline](y=test_err.data,color=:red)

output_fldr = "sqmodel_convergence"

savefig("$(output_fldr)/$(num_epocs)_$(num_hidden)_$(srand_num).pdf")
savefig("$(output_fldr)/current_fit.pdf")

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
@testset "StaticString" begin
    test_strings = ["foo, bar, buzz!", "T = Σ(δτ)", "Cars: 🚗🚕🚙"]
    test_UInt = [UInt8, UInt16, UInt32]

    @testset "Construction, Conversion and Comparison" begin
        for (i, s) in enumerate(test_strings)

            N = length(s)
            T = test_UInt[i]

            @test String(StaticString{N, T}(s))                 == s  # basic
            @test String(StaticString{N, T}(NTuple{N, T}(s)))   == s  # tuple
            @test String(StaticString{N, T}(s...))              == s  # Vararg Chars
            @test String(StaticString{N, T}(Vector{UInt8}(s)))  == s  # Vector{UInt8}
            @test String(StaticString{N, T}(collect(T, s)...))  == s  # Vararg UInt
            @test String(StaticString{N, T}(Base.CodeUnits(s))) == s  # CodeUnits
            @test String(@static_str(s))                        == s  # macro
            @test String(StaticString{0, T}(""))                == "" # empty 1
            @test String(StaticString{0, T}())                  == "" # empty 2
            @test String(StaticString{N, UInt64}(s))            == s  # larger UInt type

            if T <: Union{UInt16, UInt32}
                @test_throws InexactError StaticString{N, test_UInt[i-1]}(s)  # smaller UInt type
            end

            @test_throws DimensionMismatch StaticString{N-1, T}(s)  # too small
            @test_throws DimensionMismatch StaticString{N+1, T}(s)  # too large

            @test String(StaticString{N, T}(StaticString{N, T}(s))) == s

            # Conversion
            @test convert(StaticString{N, T}, s)::StaticString{N, T} == StaticString{N, T}(s)
            @test convert(String, StaticString{N, T}(s))::String == s

            # Equality
            @test s == StaticString{N, T}(s)
            @test StaticString{N, T}(s)      == s
            @test StaticString{N, T}(s)      == StaticString{N, T}(s)
            @test StaticString{N, UInt64}(s) == StaticString{N, T}(s)

            s_prime = let v = collect(Char, s)
                v2 = shuffle(v)
                while v2 == v; shuffle!(v2); end
                String(v2)
            end

            @test StaticString{N, T}(s_prime) != StaticString{N, T}(s)
            @test StaticString{N, T}(s) != StaticString{N, T}(s_prime)
            @test StaticString{N, T}(s_prime) != s
            @test s != StaticString{N, T}(s_prime)
        end

        @test static"bar" == StaticString{3, UInt8}("bar")

        @test isless("abc", StaticString{3, UInt8}("acc"))
        @test isless(StaticString{3, UInt16}("αβγ"), "αβδ")
    end

    @testset "AbstractString Interface" begin
        for (i, s) in enumerate(test_strings)

            N = length(s)
            T = test_UInt[i]

            stat = StaticString{N, T}(s)

            @test ncodeunits(stat) == N

            @test codeunit(stat) == T

            @test isvalid(stat, rand(1:N)) == true
            @test isvalid(stat, N+1) == false
            @test isvalid(stat, 0) == false

            @test sizeof(stat) == N * sizeof(T)

            s0 = StaticString{0, T}()
            @test isempty(s0) == true
            @test isempty(stat) == false
        end
    end

    @testset "Iteration and Indexing" begin
        for (i, s) in enumerate(test_strings)

            N = length(s)
            T = test_UInt[i]

            ss = StaticString{N, T}(s)
            v = collect(Char, s)

            @test firstindex(ss) == 1
            @test lastindex(ss) == N

            # Not sure how to test that inbounds actually works
            inbound_codeunit(j) = (local s = ss; @inbounds codeunit(s, j))
            @test_broken inbound_codeunit(0) isa UInt8

            inbound_getindex(j) = (local s = ss; @inbounds s[j])
            @test_broken inbound_getindex(0) isa Char

            for j in 1:N
                @test codeunit(ss, j) === T(v[j])
                @test inbound_codeunit(j) === T(v[j])

                @test ss[j] === v[j]
                @test inbound_getindex(j) === v[j]

                @test ss[1:j] == String(v[1:j])
                @test ss[collect(j:N)] == String(v[collect(j:N)])
            end

            let v = collect(Char, s)
                for c in ss
                    @test c isa Char
                    @test c == popfirst!(v)
                end
                @test length(v) == 0
            end

            @test Base.IteratorSize(typeof(ss)) == Base.HasLength()
            @test length(ss) == length(typeof(ss)) == N

            @test Base.IteratorEltype(typeof(ss)) == Base.HasEltype()
            @test eltype(ss) == eltype(typeof(ss)) == Char
        end
    end
end


text:
using Clang.Generators

import LAMMPS_jll

header_dir = joinpath(LAMMPS_jll.artifact_dir, "include", "lammps")
isdir(header_dir) || error("$header_dir does not exist")

const LAMMPS_INCLUDE = normpath(header_dir)

args = get_default_args()
push!(args, "-I$LAMMPS_INCLUDE")

options = load_options(joinpath(@__DIR__, "wrap.toml"))

ctx = create_context(joinpath(LAMMPS_INCLUDE, "library.h"), args, options)
build!(ctx)


text:
##### classic_degeneracy_hunter.jl
# Created by Alexander Dowling

##############
##### Degeneracy Hunter features

function degeneracyHunter(m::Model, mySolver;
	includeBounds::Bool = false,
	includeWeaklyActive::Bool = true,
	removeFixedVar::Bool = true,
	onlyCandidateSearch::Bool = true,
	epsiActive::Float64 = 1E-6,
	epsiLambda::Float64 = 1E-6,
	lambdaM::Float64 = 1E5,
	printThresholdSmallestSingularVector = 0.0,
	f=STDOUT)

	ds = DegenSettings(includeBounds, includeWeaklyActive, removeFixedVar,
		onlyCandidateSearch, epsiActive, epsiLambda, lambdaM, mySolver,
		printThresholdSmallestSingularVector)

	return degeneracyHunter(m, ds, f)

end

function degeneracyHunter(m::Model, ds::DegenSettings, f)

	println(f,"******************************************")
	println(f,"Welcome to Degeneracy Hunter!")

	println(f," ")
	println(f,"Options: ")
	println(f,"includeBounds = ",string(ds.includeBounds))
	println(f,"includeWeaklyActive = ",string(ds.includeWeaklyActive))
	println(f,"removeFixedVar = ",string(ds.removeFixedVar))
	println(f,"epsiActive = ",string(ds.epsiActive))
	println(f,"epsiLambda = ",string(ds.epsiLambda))
	println(f,"lambdaM = ",string(ds.lambdaM))
	println(f," ")

### Step 0: Setup

# Setup main data structure, evaluate full Jacobian

	dd = DegenData(m, f)

	checkJacobianScaling(m, dd, f)

### Step 1: Assemble "active" Jacobian

# Determine which constraints and bounds are active

	# Grab all equality constraints and inequality constraints that are active
	gLow = (abs.(dd.gLB - dd.g) .< ds.epsiActive)
	gUp = (abs.(dd.gUB - dd.g) .< ds.epsiActive)
	gEql = (abs.(dd.gUB - dd.gLB) .< ds.epsiActive)

	gActive = gLow .| gUp .| gEql

	# Tip: "$" was the bitwise xor operator; ⊻ is the new operator.
	actIneql = sum(xor.(gLow, gUp))
	eql = sum(gEql)
	totalIneql = length(dd.g) - eql
	inactiveIneql = totalIneql - actIneql

	println(f," ")
	println(f,"(Weakly) Active Inequality Constraints: ",string(actIneql))
	println(f,"Inactive Inequality Constraints: ",string(inactiveIneql))
	println(f,"Equality Constraints: ",string(eql))
	println(f," ")

	if(ds.includeWeaklyActive)
		# Do nothing
	else
		# Query multipliers

		# Need to implement

	end

	vFixed = abs.(m.colLower - m.colUpper) .< ds.epsiActive

	bActive = (abs.(m.colLower - dd.x) .< ds.epsiActive) .| (abs.(m.colUpper - dd.x) .< ds.epsiActive)

	println(f,"Variables: ",string(length(dd.x)))
	println(f,"(Weakly) Active Variable Bounds: ",string(sum(bActive)))
	println(f,"Fixed Variables: ",string(sum(vFixed)))
	println(f," ")

	if(ds.includeBounds)

		if(ds.includeWeaklyActive)
			# Do nothing
		else
			# Query multiplier

			# Need to implement
		end
	else
		bActive[:] = false
	end

	# Remove fixed variables from bActive
	if(ds.removeFixedVar)
		bActive = bActive .& .!vFixed
	end

# Add only active bounds to Jacobian
	print(f,"Adding Jacobian elements for bounds... ")
	tic()

	# Number of constraints
	n = length(dd.g)

	dd.bMap = find(bActive)
	for i = 1:length(dd.bMap)
		j = dd.bMap[i]
		push!(dd.iR, i + n)
		push!(dd.jC, j)
		push!(dd.J, 1.0)
	end

	tm = toq()
	println(f,string(tm," seconds"))

	print(f,"Assembling J_sparse... ")
	tic()
	J_sparse = sparse(dd.iR, dd.jC, dd.J, n + length(dd.bMap), dd.nVar)
	tm = toq()
	println(f,string(tm," seconds"))

# Remove inactive constraints and variables from the Jacobian

	dd.gMap = find(gActive)

	if(ds.removeFixedVar)
		dd.vMap = find(.!vFixed)
		cols = .!vFixed
	else
		dd.vMap = 1:nVar
		cols = trues(nVar)
	end

	rows = [gActive; ones(length(dd.bMap)) .> 0]

	print(f,"Assembling J_active... ")
	tic()
	dd.J_active = J_sparse[rows, cols]
	tm = toq()
	println(f, string(tm," seconds"))

	dd.nLambda = sum(rows)
	dd.nVarActive = sum(cols)

# Step 3: Check if any equations/bounds are degenerate themselves
# 			of if any Jacobian entries are NaN or Inf

	checkActiveJacobianRows(m, dd, f)
	checkActiveJacobianEntries(m, dd, f)

	flush(f)

# Step 3.5: Perform SVD.

	print("Computing SVD... ")
	tic()
	# This is ugly, as it converts dd.J_active from a sparse to dense matrix for SVD.
	# TODO: Use a sparse SVD implementation.

	# U*diagm(S)*V' restores matrix
	# Recall for J_active:
	#   row = constraints
	#   columns = variables
	F = LinAlg.svdfact(Matrix(dd.J_active))
	tm = toq()
	println(f, string(tm," seconds"))

	println(" ")
	println(f,"Largest Singular Value: ",F[:S][1])
	println(f,"Smallest Singular Value: ",F[:S][end])

	# Sort elements of U[:,end] from largest to smallest
	if ds.printSmallestSVD > 0

		println("Considered left singular vector for smallest singular value...")
		println("Sorted absolute value of elements from largest to smallest...")
		sigma_sorted = sortperm(abs.(F[:U][:,end]),rev=true)

		k = 1
		norm_total = 0
		norm_total_cutoff = (ds.printSmallestSVD)^2

		while k <= length(sigma_sorted) && norm_total <= norm_total_cutoff
			u_temp = F[:U][sigma_sorted[k],end]
			println("\nExamining U[",k,",end] = ",u_temp)
			printEquation(m, dd, sigma_sorted[k], f)
			printVariablesInEquation(m, dd, sigma_sorted[k], true, true, f)

			k = k + 1
			norm_total = norm_total + u_temp^2

			println("Fraction singular vector examined = ",round(100*sqrt(norm_total),2),"%")

		end

		println(" ")
	end

# Step 4: Identify candidate equations

	cand, ids = findCandidates(dd, ds, true)
	println(f,"Identified ",length(cand)," candidate constraints/bounds.")

	if(length(cand) > 0)

# Check for each candidate
		println(f,"Degenerate set from candidate search:")
		printIDS(ids, m, dd, ds, f)

		println(f,"")
		println(f,"nVar = ",string(dd.nVarActive))
		println(f,"nLambda = ",string(dd.nLambda))
		println(f,"\t Constraints = ",string(length(dd.gMap)))
		println(f,"\t Variable Bounds = ",string(length(dd.bMap)))
		println(f,"")

		# This is probably not the most efficient approach.
		# To Do: Extract row and column index vectors from J_active
		# and try working with those
		# Where is this used?
		# sparsityPattern = J_active .!= 0

		if(!ds.onlyCandidateSearch)

			sets = Array(IrreducibleDegenerateSet,length(cand))

			print(f,"Setting up MILP... ")
			tic()
			dh = setupMILP(dd, ds)
			tm = toq()
			println(f,string(tm," seconds"))

			for c = 1:length(cand)
				sets[c] = solveMILP(dh, dd, ds, c)
				printIDS(sets[c], m, dd, ds, f);
			end

		else
			sets = ids
		end

	else
		sets = Array{IrreducibleDegenerateSet}(0)
	end

	println(f," ")
	println(f,"Leaving Degeneracy Hunter.")
	println(f,"***************************")
	println(f," ")

	return sets, dd, ds

end

function setupMILP(dd::DegenData, ds::DegenSettings)

	M = ds.lambdaM

	#m2 = Model(solver=GurobiSolver(IntFeasTol=1E-6, NumericFocus=3))
	m2 = Model(solver=mySolver)
	L = 1:dd.nLambda
	@variable(m2,y[L],Bin)
	@variable(m2, -M <= lambda[L] <= M)

	@constraint(m2, degen[j=1:dd.nVarActive], sum(dd.J_active[i,j]*lambda[i] for i in L) == 0)

	@constraint(m2, lower[i=L], -M*y[i] <= lambda[i])
	@constraint(m2, upper[i=L],  lambda[i] <= M*y[i])

	@objective(m2, Min, sum(y[i] for i in L))

	return m2

end

function solveMILP(m2::Model, dd::DegenData, ds::DegenSettings, c::Int64)

	M = ds.lambdaM

	lambda = getvariable(m2, :lambda)
	L = 1:dd.nLambda

	for i = L
		setlowerbound(lambda[i],-M)
		setupperbound(lambda[i],M)
	end

	setupperbound(lambda[c],1.0)
	setlowerbound(lambda[c],1.0)

	status = solve(m2)

	if(status == :Optimal)

		lmbd_ = getvalue(lambda)
		lmbd = zeros(dd.nLambda)
		for i = L
			lmbd[i] = lmbd_[i]
		end

		y = getvariable(m2, :y)

		y_ = getvalue(y)
		ySln = zeros(dd.nLambda)
		for i = L
			ySln[i] = y_[i]
		end

		return IrreducibleDegenerateSet(c, sum(ySln), lmbd, ySln)

	else

		return IrreducibleDegenerateSet(c, 0, [], [])

	end

end

function findCandidates(dd::DegenData, ds::DegenSettings, explicit::Bool)

	if(explicit)
		mSmall = 1E-4
	else
		mSmall = 1.0
	end

	M = ds.lambdaM

	#m2 = Model(solver=Solver(IntFeasTol=1E-6, NumericFocus=3))
	m2 = Model(solver=ds.mySolver)
	L = 1:dd.nLambda
	@variable(m2, -M - mSmall <= lambda[L] <= M + mSmall)
	@variable(m2, yPos[i=L], Bin)
	@variable(m2, yNeg[i=L], Bin)
	@constraint(m2, logic[i=L], yPos[i] + yNeg[i] <= 1)

	if(explicit)

		target = 1

		@variable(m2, ySelect[i=L], Bin)
		@constraint(m2, selectLow[i=L], -(M+target)*(1-ySelect[i]) + target <= lambda[i])
		@constraint(m2, selectHigh[i=L], lambda[i] <= (M-target)*(1-ySelect[i]) + target)

		@constraint(m2, positiveUp[i=L], lambda[i] <= M*yPos[i])
		@constraint(m2, positiveLow[i=L], mSmall*yPos[i] - M*(1-yPos[i]) <= lambda[i])

		@constraint(m2, negativeLow[i=L], -M*yNeg[i] <= lambda[i])
		@constraint(m2, negativeUp[i=L], lambda[i] <= -mSmall*yNeg[i] + M*(1-yNeg[i]))

		@constraint(m2, sum(ySelect[i] for i in L) >= 1)

		@objective(m2, Min, sum(yPos[i] + yNeg[i] for i in L))

	else
		@variable(m2, 0 <= bound[i=L] <= M + mSmall)

		@constraint(m2, positiveUp[i=L], lambda[i] <= M*yPos[i] + mSmall)
		@constraint(m2, negativeLow[i=L], -M*yNeg[i] - mSmall <= lambda[i])

		# When yPos = 1, lambda >= mSmall
		# When yPos = 0, lambda >= -mSmall
		@constraint(m2, positiveLow[i=L], -mSmall + 2*mSmall*yPos[i] <= lambda[i])

		# When yNeg = 1, lambda <= -mSmall
		# When yNeg = 0, lambda <= mSmall
		@constraint(m2, negativeUp[i=L], lambda[i] <= mSmall - 2*mSmall*yNeg[i])

		@constraint(m2, lower[i=L], -bound[i] <= lambda[i])
		@constraint(m2, upper[i=L], lambda[i] <= bound[i])

		# At least one constraint must be in the degenerate set
		@constraint(m2, sum(yPos[i] + yNeg[i] for i in L) >= 1)

		@objective(m2, Min, sum(bound[i] for i in L))
	end

	@constraint(m2, degen[j=1:dd.nVarActive], sum(dd.J_active[i,j]*lambda[i] for i in L) == 0)

	status = solve(m2)

	c = Array{Int64}(0)

	if(status == :Optimal)

		lmbd_ = getvalue(lambda)
		lmbd = zeros(dd.nLambda)
		for i = L
			lmbd[i] = lmbd_[i]
		end

		abs_lmbd = abs.(lmbd)

		lmax = 0
		for i = L
			if(lmax < abs_lmbd[i])
				lmax = abs_lmbd[i]
			end
		end

		for i = L
			if(abs_lmbd[i] > 0.1*lmax)
				push!(c,i)
			end
		end

		y = 1.0*(abs_lmbd .> 1/M)

		ids = IrreducibleDegenerateSet(0, sum(y), lmbd, y)

	else
		ids = IrreducibleDegenerateSet(0, 0.0, zeros(0), zeros(0))
	end

	return c, ids

end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     cell_metadata_filter: -all
#     formats: jl:light,ipynb
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.3.2
#   kernelspec:
#     display_name: Julia 1.1.0
#     language: julia
#     name: julia-1.1
# ---

import DrWatson: quickactivate
quickactivate(@__DIR__, "Chemostat_Rath2017")
Base.active_project() |> println

# +
using Distributed
using Serialization
using Dates

nworkers = 1 # length(Sys.cpu_info()) - 1
length(workers()) < nworkers && addprocs(nworkers; exeflags = "--project")
println("Working in: ", workers())
# -
@everywhere begin
    using Serialization
    using Distributed

    import Chemostat
    Ch = Chemostat
end

# ### Initializing everywhere

@everywhere begin

    const obj_ider = "biom";
    const cost_ider ="tot_cost"

    const ep_alpha = 1e9
    const ep_epsconv = 1e-9 # This is an exageration!!
    const ep_ochlen = 100

    const hr_b0_nsamples = 100_000 #10_000_000
    const hr_nsamples = 1000 # 10_000
    const intake_info = Dict("gt" => Dict("ub" => 100.0, "c" => 10.0));

    const ξs = Ch.Utils.logspace(-1,1,20);
    const ξs_str = round.(ξs, digits = 2);

    const β0 = 0.0 # It is important to include beta = 0 if doing HR
    const βs = [β0; Ch.Utils.logspace(-1, 3.8, 20)]
    const βs_str = round.(βs, digits = 2);
end

println("ep_alpha: ", ep_alpha)
println("ep_epsconv: ", ep_epsconv)
println("βs: ", βs_str)
println("ξs: ", ξs_str)

# RES IDS
# Collect all the computed results ids for bundling
const chnl = RemoteChannel() do
    Channel{Any}(10)
end
const res_ids = []
const collector = @async while true
    id = take!(chnl)
    push!(res_ids, id)
end;

@everywhere begin
    const base_model = Ch.Test.toy_model();
    function prepare_model(ξ)
        model = base_model |> deepcopy
        model = Ch.SteadyState.apply_bound!(model, ξ, intake_info |> deepcopy)
        model = Ch.LP.fva_preprocess(model; verbose = false)
    end
end

# cache dir
cache_dir = "cache"
!isdir(cache_dir) && mkdir(cache_dir)
@assert isdir(cache_dir)
Ch.Utils.set_cache_dir(cache_dir)

# ### EP

# SIMULATION
# Any of the loops can be parallelized by just
# changing one of the 'map' functions
map(ξs) do ξ

    # HASH SIM
    sim_hash = (:TOY_SIMULATION, ξ)

    # FBA EP
    dat = Ch.SimulationUtils.cached_simulation(;
        epochlen = ep_ochlen,
        verbose = true,
        sim_id = sim_hash,
        get_model = () -> prepare_model(ξ),
        objider = obj_ider,
        beta_info = [(obj_ider, βs)],
        costider = cost_ider,
        clear_cache = false,
        use_seed = true,
        epmodel_kwargs = Dict(:alpha => ep_alpha),
        epconv_kwargs = Dict(:epsconv => ep_epsconv)
    )

    ## HR
    model = prepare_model(ξ)
    hrout_β0 = nothing # seed
    Ch.Utils.tagprintln_inmw("DOING HR ")
    for (βi, β) in enumerate(βs)

        Ch.Utils.tagprintln_inmw("DOING BETA ", βi)
        if β == β0 # hrout seed
                hrout_β0 = Ch.MaxEntHR.maxent_hr(model,
                            nsamples = hr_b0_nsamples,
                            drop_samples = false, verbose = false);
            hrout = Ch.Utils.HRout(hrout_β0, drop_samples = true) # a copy without samples
        else
            isnothing(hrout_β0) && (error("The minimum β must be cero!!!"); return)

            hrout = Ch.MaxEntHR.maxent_hr(model, hrout_β0, obj_ider,
                β * 10.0^(-sqrt(ℯ)); # I have no idea why these work! TODO find it!
                nsamples = hr_nsamples,
                maxiter = 1e16,
                verbose = false
            )
        end
        dat[(:hr, βi)] = hrout
    end

    ## SAVING DATA
    res_id = (:RESULT, sim_hash)
    Ch.Utils.save_cache(res_id, (ξ, βs, model, dat);
        headline = "CATCHING RESULTS\n")

    ## PASSING ID TO MASTER
    put!(chnl, res_id)

    GC.gc()
    return nothing
end; # map(ξs) do ξ

for ξ in ξs
    m = prepare_model(ξ)
    Ch.LP.fba(m, obj_ider, cost_ider)
end

# ### Parallel loop

# COLLECTING RESULTS
Ch.Utils.tagprintln_inmw("COLLECTING RESULTS ")
sleep(1) # wait for collector to get all ids
const bundle = Ch.Utils.ChstatBundle()
for id in res_ids

    ξ_, βs_, model_, dat_ = Ch.Utils.load_cache(id; verbose = false)

    bundle[ξ_, :net] = model_
    bundle[ξ_, :fba] = dat_[:fba]

    for (βi, β_) in βs_ |> enumerate
        bundle[ξ_, β_, :ep] = dat_[(:ep, βi)]
        bundle[ξ_, β_, :hr] = dat_[(:hr, βi)]
    end
end

# ### Saving results

### Catching
cache_file = "toy_model_cache.jls"
serialize(cache_file, bundle)
println(relpath(cache_file), " created!!!")

# ### Deleting temporal caches

Ch.Utils.tagprintln_inmw("DELETING TEMPORAL CACHE ")
Ch.Utils.delete_temp_caches()

max_stars_repo_name: josePereiro/Chemostat

max_issues_repo_name: josePereiro/Chemostat

max_forks_repo_name: josePereiro/Chemostat


text:
# Definition of an orthonormal basis
"""
    OrthonormalBasis{T} <: Basis{T}

A list of vector like objects of type `T` that are mutually orthogonal and normalized to
one, representing an orthonormal basis for some subspace (typically a Krylov subspace). See
also [`Basis`](@ref)

Orthonormality of the vectors contained in an instance `b` of `OrthonormalBasis`
(i.e. `all(dot(b[i],b[j]) == I[i,j] for i=1:length(b), j=1:length(b))`) is not checked when
elements are added; it is up to the algorithm that constructs `b` to guarantee
orthonormality.

One can easily orthogonalize or orthonormalize a given vector `v` with respect to a
`b::OrthonormalBasis` using the functions
[`w, = orthogonalize(v,b,...)`](@ref orthogonalize) or
[`w, = orthonormalize(v,b,...)`](@ref orthonormalize). The resulting vector `w` of the
latter can then be added to `b` using `push!(b, w)`. Note that in place versions
[`orthogonalize!(v, b, ...)`](@ref orthogonalize) or
[`orthonormalize!(v, b, ...)`](@ref orthonormalize) are also available.

Finally, a linear combination of the vectors in `b::OrthonormalBasis` can be obtained by
multiplying `b` with a `Vector{<:Number}` using `*` or `mul!` (if the output vector is
already allocated).
"""
struct OrthonormalBasis{T} <: Basis{T}
    basis::Vector{T}
end
OrthonormalBasis{T}() where {T} = OrthonormalBasis{T}(Vector{T}(undef, 0))

# Iterator methods for OrthonormalBasis
Base.IteratorSize(::Type{<:OrthonormalBasis}) = Base.HasLength()
Base.IteratorEltype(::Type{<:OrthonormalBasis}) = Base.HasEltype()

Base.length(b::OrthonormalBasis) = length(b.basis)
Base.eltype(b::OrthonormalBasis{T}) where {T} = T

Base.iterate(b::OrthonormalBasis) = Base.iterate(b.basis)
Base.iterate(b::OrthonormalBasis, state) = Base.iterate(b.basis, state)

Base.getindex(b::OrthonormalBasis, i) = getindex(b.basis, i)
Base.setindex!(b::OrthonormalBasis, i, q) = setindex!(b.basis, i, q)
Base.firstindex(b::OrthonormalBasis) = firstindex(b.basis)
Base.lastindex(b::OrthonormalBasis) = lastindex(b.basis)

Base.first(b::OrthonormalBasis) = first(b.basis)
Base.last(b::OrthonormalBasis) = last(b.basis)

Base.popfirst!(b::OrthonormalBasis) = popfirst!(b.basis)
Base.pop!(b::OrthonormalBasis) = pop!(b.basis)
Base.push!(b::OrthonormalBasis{T}, q::T) where {T} = (push!(b.basis, q); return b)
Base.empty!(b::OrthonormalBasis) = (empty!(b.basis); return b)
Base.sizehint!(b::OrthonormalBasis, k::Int) = (sizehint!(b.basis, k); return b)
Base.resize!(b::OrthonormalBasis, k::Int) = (resize!(b.basis, k); return b)

# Multiplication methods with OrthonormalBasis
function Base.:*(b::OrthonormalBasis, x::AbstractVector)
    y = zero(eltype(x)) * first(b)
    return mul!(y, b, x)
end
LinearAlgebra.mul!(y, b::OrthonormalBasis, x::AbstractVector) = unproject!(y, b, x, 1, 0)

const BLOCKSIZE = 4096

"""
    project!(y::AbstractVector, b::OrthonormalBasis, x,
        [α::Number = 1, β::Number = 0, r = Base.OneTo(length(b))])

For a given orthonormal basis `b`, compute the expansion coefficients `y` resulting from
projecting the vector `x` onto the subspace spanned by `b`; more specifically this computes

```
    y[j] = β*y[j] + α * dot(b[r[j]], x)
```

for all ``j ∈ r``.
"""
function project!(
    y::AbstractVector,
    b::OrthonormalBasis,
    x,
    α::Number = true,
    β::Number = false,
    r = Base.OneTo(length(b))
)
    # no specialized routine for IndexLinear x because reduction dimension is large dimension
    length(y) == length(r) || throw(DimensionMismatch())
    if get_num_threads() > 1
        @sync for J in splitrange(1:length(r), get_num_threads())
            Threads.@spawn for j in $J
                @inbounds begin
                    if β == 0
                        y[j] = α * dot(b[r[j]], x)
                    else
                        y[j] = β * y[j] + α * dot(b[r[j]], x)
                    end
                end
            end
        end
    else
        for j in 1:length(r)
            @inbounds begin
                if β == 0
                    y[j] = α * dot(b[r[j]], x)
                else
                    y[j] = β * y[j] + α * dot(b[r[j]], x)
                end
            end
        end
    end
    return y
end

"""
    unproject!(y, b::OrthonormalBasis, x::AbstractVector,
        [α::Number = 1, β::Number = 0, r = Base.OneTo(length(b))])

For a given orthonormal basis `b`, reconstruct the vector-like object `y` that is defined by
expansion coefficients with respect to the basis vectors in `b` in `x`; more specifically
this computes

```
    y = β*y + α * sum(b[r[i]]*x[i] for i = 1:length(r))
```
"""
function unproject!(
    y,
    b::OrthonormalBasis,
    x::AbstractVector,
    α::Number = true,
    β::Number = false,
    r = Base.OneTo(length(b))
)
    if y isa AbstractArray && IndexStyle(y) isa IndexLinear && get_num_threads() > 1
        return unproject_linear_multithreaded!(y, b, x, α, β, r)
    end
    # general case: using only vector operations, i.e. axpy! (similar to BLAS level 1)
    length(x) == length(r) || throw(DimensionMismatch())
    if β == 0
        rmul!(y, false) # should be hard zero
    elseif β != 1
        rmul!(y, β)
    end
    @inbounds for (i, ri) in enumerate(r)
        y = axpy!(α * x[i], b[ri], y)
    end
    return y
end
function unproject_linear_multithreaded!(
    y::AbstractArray,
    b::OrthonormalBasis{<:AbstractArray},
    x::AbstractVector,
    α::Number = true,
    β::Number = false,
    r = Base.OneTo(length(b))
)
    # multi-threaded implementation, similar to BLAS level 2 matrix vector multiplication
    m = length(y)
    n = length(r)
    length(x) == n || throw(DimensionMismatch())
    for rj in r
        length(b[rj]) == m || throw(DimensionMismatch())
    end
    if n == 0
        return β == 1 ? y : β == 0 ? fill!(y, 0) : rmul!(y, β)
    end
    let m = m, n = n, y = y, x = x, b = b, blocksize = prevpow(2, div(BLOCKSIZE, n))
        @sync for II in splitrange(1:blocksize:m, get_num_threads())
            Threads.@spawn for I in $II
                unproject_linear_kernel!(y, b, x, I:min(I + blocksize - 1, m), α, β, r)
            end
        end
    end
    return y
end
function unproject_linear_kernel!(
    y::AbstractArray,
    b::OrthonormalBasis{<:AbstractArray},
    x::AbstractVector,
    I,
    α::Number,
    β::Number,
    r
)
    @inbounds begin
        if β == 0
            @simd for i in I
                y[i] = zero(y[i])
            end
        elseif β != 1
            @simd for i in I
                y[i] *= β
            end
        end
        for (j, rj) in enumerate(r)
            xj = x[j] * α
            Vj = b[rj]
            @simd for i in I
                y[i] += Vj[i] * xj
            end
        end
    end
end

"""
    rank1update!(b::OrthonormalBasis, y, x::AbstractVector,
        [α::Number = 1, β::Number = 1, r = Base.OneTo(length(b))])

Perform a rank 1 update of a basis `b`, i.e. update the basis vectors as

```
    b[r[i]] = β*b[r[i]] + α * y * conj(x[i])
```

It is the user's responsibility to make sure that the result is still an orthonormal basis.
"""
@fastmath function rank1update!(
    b::OrthonormalBasis,
    y,
    x::AbstractVector,
    α::Number = true,
    β::Number = true,
    r = Base.OneTo(length(b))
)
    if y isa AbstractArray && IndexStyle(y) isa IndexLinear && Threads.nthreads() > 1
        return rank1update_linear_multithreaded!(b, y, x, α, β, r)
    end
    # general case: using only vector operations, i.e. axpy! (similar to BLAS level 1)
    length(x) == length(r) || throw(DimensionMismatch())
    @inbounds for (i, ri) in enumerate(r)
        if β == 1
            b[ri] = axpy!(α * conj(x[i]), y, b[ri])
        elseif β == 0
            b[ri] = mul!(b[ri], α * x[i], y)
        else
            b[ri] = axpby!(α * x[i], y, β, b[ri])
        end
    end
    return b
end
@fastmath function rank1update_linear_multithreaded!(
    b::OrthonormalBasis{<:AbstractArray},
    y::AbstractArray,
    x::AbstractVector,
    α::Number,
    β::Number,
    r
)
    # multi-threaded implementation, similar to BLAS level 2 matrix vector multiplication
    m = length(y)
    n = length(r)
    length(x) == n || throw(DimensionMismatch())
    for rj in r
        length(b[rj]) == m || throw(DimensionMismatch())
    end
    if n == 0
        return b
    end
    let m = m, n = n, y = y, x = x, b = b, blocksize = prevpow(2, div(BLOCKSIZE, n))
        @sync for II in splitrange(1:blocksize:m, get_num_threads())
            Threads.@spawn for I in $II
                @inbounds begin
                    for (j, rj) in enumerate(r)
                        xj = α * conj(x[j])
                        Vj = b[rj]
                        if β == 0
                            @simd for i in I:min(I + blocksize - 1, m)
                                Vj[i] = zero(Vj[i])
                            end
                        elseif β != 1
                            @simd for i in I:min(I + blocksize - 1, m)
                                Vj[i] *= β
                            end
                        end
                        if I + blocksize - 1 <= m
                            @simd for i in Base.OneTo(blocksize)
                                Vj[I-1+i] += y[I-1+i] * xj
                            end
                        else
                            @simd for i in I:m
                                Vj[i] += y[i] * xj
                            end
                        end
                    end
                end
            end
        end
    end
    return b
end

"""
    basistransform!(b::OrthonormalBasis, U::AbstractMatrix)

Transform the orthonormal basis `b` by the matrix `U`. For `b` an orthonormal basis,
the matrix `U` should be real orthogonal or complex unitary; it is up to the user to ensure
this condition is satisfied. The new basis vectors are given by

```
    b[j] ← b[i] * U[i,j]
```

and are stored in `b`, so the old basis vectors are thrown away. Note that, by definition,
the subspace spanned by these basis vectors is exactly the same.
"""
function basistransform!(b::OrthonormalBasis{T}, U::AbstractMatrix) where {T} # U should be unitary or isometric
    if T <: AbstractArray && IndexStyle(T) isa IndexLinear && get_num_threads() > 1
        return basistransform_linear_multithreaded!(b, U)
    end
    m, n = size(U)
    m == length(b) || throw(DimensionMismatch())

    let b2 = [similar(b[1]) for j in 1:n]
        if get_num_threads() > 1
            @sync for J in splitrange(1:n, get_num_threads())
                Threads.@spawn for j in $J
                    mul!(b2[j], b[1], U[1, j])
                    for i in 2:m
                        axpy!(U[i, j], b[i], b2[j])
                    end
                end
            end
        else
            for j in 1:n
                mul!(b2[j], b[1], U[1, j])
                for i in 2:m
                    axpy!(U[i, j], b[i], b2[j])
                end
            end
        end
        for j in 1:n
            b[j] = b2[j]
        end
    end
    return b
end

function basistransform_linear_multithreaded!(
    b::OrthonormalBasis{<:AbstractArray},
    U::AbstractMatrix
) # U should be unitary or isometric
    m, n = size(U)
    m == length(b) || throw(DimensionMismatch())
    K = length(b[1])

    blocksize = prevpow(2, div(BLOCKSIZE, m))
    let b2 = [similar(b[1]) for j in 1:n], K = K, m = m, n = n
        @sync for II in splitrange(1:blocksize:K, get_num_threads())
            Threads.@spawn for I in $II
                @inbounds for j in 1:n
                    b2j = b2[j]
                    @simd for i in I:min(I + blocksize - 1, K)
                        b2j[i] = zero(b2j[i])
                    end
                    for k in 1:m
                        bk = b[k]
                        Ukj = U[k, j]
                        @simd for i in I:min(I + blocksize - 1, K)
                            b2j[i] += bk[i] * Ukj
                        end
                    end
                end
            end
        end
        for j in 1:n
            b[j] = b2[j]
        end
    end
    return b
end

# function basistransform2!(b::OrthonormalBasis, U::AbstractMatrix) # U should be unitary or isometric
#     m, n = size(U)
#     m == length(b) || throw(DimensionMismatch())
#
#     # apply basis transform via householder reflections
#     for j = 1:size(U,2)
#         h, ν = householder(U, j:m, j)
#         lmul!(h, view(U, :, j+1:n))
#         rmul!(b, h')
#     end
#     return b
# end

# Orthogonalization of a vector against a given OrthonormalBasis
orthogonalize(v, args...) = orthogonalize!(true*v, args...)

function orthogonalize!(v::T, b::OrthonormalBasis{T}, alg::Orthogonalizer) where {T}
    S = promote_type(eltype(v), eltype(T))
    c = Vector{S}(undef, length(b))
    return orthogonalize!(v, b, c, alg)
end

function orthogonalize!(
    v::T,
    b::OrthonormalBasis{T},
    x::AbstractVector,
    ::ClassicalGramSchmidt
) where {T}
    x = project!(x, b, v)
    v = unproject!(v, b, x, -1, 1)
    return (v, x)
end
function reorthogonalize!(
    v::T,
    b::OrthonormalBasis{T},
    x::AbstractVector,
    ::ClassicalGramSchmidt
) where {T}
    s = similar(x) ## EXTRA ALLOCATION
    s = project!(s, b, v)
    v = unproject!(v, b, s, -1, 1)
    x .+= s
    return (v, x)
end
function orthogonalize!(
    v::T,
    b::OrthonormalBasis{T},
    x::AbstractVector,
    ::ClassicalGramSchmidt2
) where {T}
    (v, x) = orthogonalize!(v, b, x, ClassicalGramSchmidt())
    return reorthogonalize!(v, b, x, ClassicalGramSchmidt())
end
function orthogonalize!(
    v::T,
    b::OrthonormalBasis{T},
    x::AbstractVector,
    alg::ClassicalGramSchmidtIR
) where {T}
    nold = norm(v)
    orthogonalize!(v, b, x, ClassicalGramSchmidt())
    nnew = norm(v)
    while eps(one(nnew)) < nnew < alg.η * nold
        nold = nnew
        (v, x) = reorthogonalize!(v, b, x, ClassicalGramSchmidt())
        nnew = norm(v)
    end
    return (v, x)
end

function orthogonalize!(
    v::T,
    b::OrthonormalBasis{T},
    x::AbstractVector,
    ::ModifiedGramSchmidt
) where {T}
    for (i, q) in enumerate(b)
        s = dot(q, v)
        v = axpy!(-s, q, v)
        x[i] = s
    end
    return (v, x)
end
function reorthogonalize!(
    v::T,
    b::OrthonormalBasis{T},
    x::AbstractVector,
    ::ModifiedGramSchmidt
) where {T}
    for (i, q) in enumerate(b)
        s = dot(q, v)
        v = axpy!(-s, q, v)
        x[i] += s
    end
    return (v, x)
end
function orthogonalize!(
    v::T,
    b::OrthonormalBasis{T},
    x::AbstractVector,
    ::ModifiedGramSchmidt2
) where {T}
    (v, x) = orthogonalize!(v, b, x, ModifiedGramSchmidt())
    return reorthogonalize!(v, b, x, ModifiedGramSchmidt())
end
function orthogonalize!(
    v::T,
    b::OrthonormalBasis{T},
    x::AbstractVector,
    alg::ModifiedGramSchmidtIR
) where {T}
    nold = norm(v)
    (v, x) = orthogonalize!(v, b, x, ModifiedGramSchmidt())
    nnew = norm(v)
    while eps(one(nnew)) < nnew < alg.η * nold
        nold = nnew
        (v, x) = reorthogonalize!(v, b, x, ModifiedGramSchmidt())
        nnew = norm(v)
    end
    return (v, x)
end

# Orthogonalization of a vector against a given normalized vector
orthogonalize!(v::T, q::T, alg::Orthogonalizer) where {T} = _orthogonalize!(v, q, alg)
# avoid method ambiguity on Julia 1.0 according to Aqua.jl

function _orthogonalize!(
    v::T,
    q::T,
    alg::Union{ClassicalGramSchmidt,ModifiedGramSchmidt}
) where {T}
    s = dot(q, v)
    v = axpy!(-s, q, v)
    return (v, s)
end
function _orthogonalize!(
    v::T,
    q::T,
    alg::Union{ClassicalGramSchmidt2,ModifiedGramSchmidt2}
) where {T}
    s = dot(q, v)
    v = axpy!(-s, q, v)
    ds = dot(q, v)
    v = axpy!(-ds, q, v)
    return (v, s + ds)
end
function _orthogonalize!(
    v::T,
    q::T,
    alg::Union{ClassicalGramSchmidtIR,ModifiedGramSchmidtIR}
) where {T}
    nold = norm(v)
    s = dot(q, v)
    v = axpy!(-s, q, v)
    nnew = norm(v)
    while eps(one(nnew)) < nnew < alg.η * nold
        nold = nnew
        ds = dot(q, v)
        v = axpy!(-ds, q, v)
        s += ds
        nnew = norm(v)
    end
    return (v, s)
end

"""
    orthogonalize(v, b::OrthonormalBasis, [x::AbstractVector,] alg::Orthogonalizer]) -> w, x
    orthogonalize!(v, b::OrthonormalBasis, [x::AbstractVector,] alg::Orthogonalizer]) -> w, x

    orthogonalize(v, q, algorithm::Orthogonalizer]) -> w, s
    orthogonalize!(v, q, algorithm::Orthogonalizer]) -> w, s

Orthogonalize vector `v` against all the vectors in the orthonormal basis `b` using the
orthogonalization algorithm `alg` of type [`Orthogonalizer`](@ref), and return the resulting
vector `w` and the overlap coefficients `x` of `v` with the basis vectors in `b`.

In case of `orthogonalize!`, the vector `v` is mutated in place. In both functions, storage
for the overlap coefficients `x` can be provided as optional argument `x::AbstractVector`
with `length(x) >= length(b)`.

One can also orthogonalize `v` against a given vector `q` (assumed to be normalized), in
which case the orthogonal vector `w` and the inner product `s` between `v` and `q` are
returned.

Note that `w` is not normalized, see also [`orthonormalize`](@ref).

For more information on possible orthogonalization algorithms, see [`Orthogonalizer`](@ref)
and its concrete subtypes [`ClassicalGramSchmidt`](@ref), [`ModifiedGramSchmidt`](@ref),
[`ClassicalGramSchmidt2`](@ref), [`ModifiedGramSchmidt2`](@ref),
[`ClassicalGramSchmidtIR`](@ref) and [`ModifiedGramSchmidtIR`](@ref).
"""
orthogonalize, orthogonalize!

# Orthonormalization: orthogonalization and normalization
orthonormalize(v, args...) = orthonormalize!(true*v, args...)

function orthonormalize!(v, args...)
    out = orthogonalize!(v, args...) # out[1] === v
    β = norm(v)
    v = rmul!(v, inv(β))
    return tuple(v, β, Base.tail(out)...)
end

"""
    orthonormalize(v, b::OrthonormalBasis, [x::AbstractVector,] alg::Orthogonalizer]) -> w, β, x
    orthonormalize!(v, b::OrthonormalBasis, [x::AbstractVector,] alg::Orthogonalizer]) -> w, β, x

    orthonormalize(v, q, algorithm::Orthogonalizer]) -> w, β, s
    orthonormalize!(v, q, algorithm::Orthogonalizer]) -> w, β, s

Orthonormalize vector `v` against all the vectors in the orthonormal basis `b` using the
orthogonalization algorithm `alg` of type [`Orthogonalizer`](@ref), and return the resulting
vector `w` (of norm 1), its norm `β` after orthogonalizing and the overlap coefficients `x`
of `v` with the basis vectors in `b`, such that `v = β * w + b * x`.

In case of `orthogonalize!`, the vector `v` is mutated in place. In both functions, storage
for the overlap coefficients `x` can be provided as optional argument `x::AbstractVector`
with `length(x) >= length(b)`.

One can also orthonormalize `v` against a given vector `q` (assumed to be normalized), in
which case the orthonormal vector `w`, its norm `β` before normalizing and the inner product
`s` between `v` and `q` are returned.

See [`orthogonalize`](@ref) if `w` does not need to be normalized.

For more information on possible orthogonalization algorithms, see [`Orthogonalizer`](@ref)
and its concrete subtypes [`ClassicalGramSchmidt`](@ref), [`ModifiedGramSchmidt`](@ref),
[`ClassicalGramSchmidt2`](@ref), [`ModifiedGramSchmidt2`](@ref),
[`ClassicalGramSchmidtIR`](@ref) and [`ModifiedGramSchmidtIR`](@ref).
"""
orthonormalize, orthonormalize!


text:
## Parser combinator from strings to Julia expressions

reader_table = Dict{Symbol, Function}()

"Parser combinator for Lisp syntax."
lisp         = Delayed()
floaty_dot   = p"[-+]?[0-9]*\.[0-9]+([eE][-+]?[0-9]+)?[Ff]" > (x -> parse(Float32, x[1:end-1]))
floaty_nodot = p"[-+]?[0-9]*[0-9]+([eE][-+]?[0-9]+)?[Ff]" > (x -> parse(Float32, x[1:end-1]))
floaty       = floaty_dot | floaty_nodot
white_space  = p"(([\s\n\r]*(?<!\\);[^\n\r$]+[\n\r\s$]*)+|[\s\n\r]+)"
opt_ws       = white_space | e""

doubley      = p"[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?[dD]" > (x -> parse(Float64, x[1:end-1]))

inty         = p"[-+]?\d+" > (x -> parse(Int, x))

uchary       = p"\\(u[\da-fA-F]{4})" > (x -> first(unescape_string(x)))
achary       = p"\\[0-7]{3}" > (x -> unescape_string(x)[1])
chary        = p"\\." > (x -> x[2])

stringy      = p"(?<!\\)\".*?(?<!\\)\"" > (x -> x[2:end-1]) #_0[2:end-1] } #r"(?<!\\)\".*?(?<!\\)"
booly        = p"(true|false)" > (x -> x == "true" ? true : false)
symboly      = p"[^\d():\?{}#'`,@~;~\[\]^\s][^\s:\?()#'`,@~;^{}~\[\]]*" > Symbol
macrosymy    = p"@[^\d():\?{}#'`,@~;~\[\]^\s][^\s:\?()#'`,@~;^{}~\[\]]*" > Symbol

sexpr        = E"(" + ~opt_ws + Repeat(lisp + ~opt_ws) + E")" |> (x -> x)
hashy        = E"#{" + ~opt_ws + Repeat(lisp + ~opt_ws) + E"}" |> (x -> Set(x))
curly        = E"{" + ~opt_ws + Repeat(lisp + ~opt_ws) + E"}" |> (x -> Dict(x[i] => x[i+1] for i = 1:2:length(x)))
dispatchy    = E"#" + symboly + ~opt_ws + lisp |> (x -> reader_table[x[1]](x[2]))
bracket      = E"[" + ~opt_ws + Repeat(lisp + ~opt_ws) + E"]" |> (x -> x)

# Additional combinators to handle PDDL-specific syntax
vary         = p"\?[^\d():\?{}#'`,@~;~\[\]^\s][^\s():\?#'`,@~;^{}~\[\]]*" > (x -> Var(Symbol(uppercasefirst(x[2:end]))))
keywordy     = p":[^\d():\?{}#'`,@~;~\[\]^\s][^\s():\?#'`,@~;^{}~\[\]]*" > (x -> Keyword(Symbol(x[2:end])))

lisp.matcher = doubley | floaty | inty | uchary | achary | chary | stringy | booly |
               vary | keywordy | symboly | macrosymy |
               dispatchy | sexpr | hashy | curly | bracket

const top_level = Repeat(~opt_ws + lisp) + ~opt_ws + Eos()

parse_string(str::AbstractString) = parse_one(lowercase(str), top_level)[1]


text:
macro typenames()
  esc(quote
    uname = options[:longform] ? ((FSS < 7) ? "UnumSmall" : "UnumLarge") : "Unum"
    gname = options[:longform] ? ((FSS < 7) ? "GnumSmall" : "GnumLarge") : "Gnum"
    bname = options[:longform] ? ((FSS < 7) ? "UboundSmall" : "UboundLarge") : "Ubound"
  end)
end

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
function stmt_effect_free(@nospecialize(stmt), src, mod::Module)
    isa(stmt, Union{PiNode, PhiNode}) && return true
    isa(stmt, Union{ReturnNode, GotoNode, GotoIfNot}) && return false
    return effect_free(stmt, src, mod, true)
end

function abstract_eval_ssavalue(s::SSAValue, src::IRCode)
    return src.types[s.id]
end

function abstract_eval_ssavalue(s::SSAValue, src::IncrementalCompact)
    return types(src)[s]
end

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
module FIGlet

using Pkg

if isdefined(Pkg, :Artifacts)
    using Pkg.Artifacts
    @eval fontsroot = artifact"fonts"
else
    fontsroot = normpath(@__DIR__, "..", "deps")
end
const FONTSDIR = abspath(normpath(joinpath(fontsroot, "FIGletFonts-0.5.0", "fonts")))
const UNPARSEABLES = [
              "nvscript.flf",
             ]

abstract type FIGletError <: Exception end

"""
Width is not sufficient to print a character
"""
struct CharNotPrinted <: FIGletError end

"""
Font can't be located
"""
struct FontNotFoundError <: FIGletError
    msg::String
end

Base.showerror(io::IO, e::FontNotFoundError) = print(io, "FontNotFoundError: $(e.msg)")

"""
Problem parsing a font file
"""
struct FontError <: FIGletError
    msg::String
end

Base.showerror(io::IO, e::FontError) = print(io, "FontError: $(e.msg)")

"""
Color is invalid
"""
struct InvalidColorError <: FIGletError end

Base.@enum(Layout,
    FullWidth                   =       -1,
    HorizontalSmushingRule1     =        1,
    HorizontalSmushingRule2     =        2,
    HorizontalSmushingRule3     =        4,
    HorizontalSmushingRule4     =        8,
    HorizontalSmushingRule5     =       16,
    HorizontalSmushingRule6     =       32,
    HorizontalFitting           =       64,
    HorizontalSmushing          =      128,
    VerticalSmushingRule1       =      256,
    VerticalSmushingRule2       =      512,
    VerticalSmushingRule3       =     1024,
    VerticalSmushingRule4       =     2048,
    VerticalSmushingRule5       =     4096,
    VerticalFitting             =     8192,
    VerticalSmushing            =    16384,
)

struct FIGletHeader
    hardblank::Char
    height::Int
    baseline::Int
    max_length::Int
    old_layout::Int
    comment_lines::Int
    print_direction::Int
    full_layout::Int
    codetag_count::Int

    function FIGletHeader(
                          hardblank,
                          height,
                          baseline,
                          max_length,
                          old_layout,
                          comment_lines,
                          print_direction=0,
                          full_layout=Int(HorizontalSmushingRule2),
                          codetag_count=0,
                          args...,
                      )
        length(args) >0 && @warn "Received unknown header attributes: `$args`."
        height < 1 && ( height = 1 )
        max_length < 1 && ( max_length = 1 )
        print_direction < 0 && ( print_direction = 0 )
        # max_length += 100 # Give ourselves some extra room
        new(hardblank, height, baseline, max_length, old_layout, comment_lines, print_direction, full_layout, codetag_count)
    end
end

function FIGletHeader(
                      hardblank,
                      height::AbstractString,
                      baseline::AbstractString,
                      max_length::AbstractString,
                      old_layout::AbstractString,
                      comment_lines::AbstractString,
                      print_direction::AbstractString="0",
                      full_layout::AbstractString="2",
                      codetag_count::AbstractString="0",
                      args...,
                     )
    return FIGletHeader(
                        hardblank,
                        parse(Int, height),
                        parse(Int, baseline),
                        parse(Int, max_length),
                        parse(Int, old_layout),
                        parse(Int, comment_lines),
                        parse(Int, print_direction),
                        parse(Int, full_layout),
                        parse(Int, codetag_count),
                        args...,
                       )
end

struct FIGletChar
    ord::Char
    thechar::Matrix{Char}
end

struct FIGletFont
    header::FIGletHeader
    font_characters::Dict{Char,FIGletChar}
    version::VersionNumber
end

Base.show(io::IO, ff::FIGletFont) = print(io, "FIGletFont(n=$(length(ff.font_characters)))")

function readmagic(io)
    magic = read(io, 5)
    magic[1:4] != UInt8['f', 'l', 'f', '2'] && throw(FontError("File is not a valid FIGlet Lettering Font format. Magic header values must start with `flf2`."))
    magic[5] != UInt8('a') && @warn "File may be a FLF format but not flf2a."
    return magic # File has valid FIGlet Lettering Font format magic header.
end

function readfontchar(io, ord, height)

    s = readline(io)
    width = length(s)-1
    width == -1 && throw(FontError("Unable to find character `$ord` in FIGlet Font."))
    thechar = Matrix{Char}(undef, height, width)

    for (w, c) in enumerate(s)
        w > width && break
        thechar[1, w] = c
    end

    for h in 2:height
        s = readline(io)
        for (w, c) in enumerate(s)
            w > width && break
            thechar[h, w] = c
        end
    end

    return FIGletChar(ord, thechar)
end

Base.show(io::IO, fc::FIGletChar) = print(io, "FIGletChar(ord='$(fc.ord)')")

function readfont(s::AbstractString)
    name = s
    if !isfile(name)
        name = abspath(normpath(joinpath(FONTSDIR, name)))
        if !isfile(name)
            name = "$name.flf"
            !isfile(name) && throw(FontNotFoundError("Cannot find font `$s`."))
        end
    end

    font = open(name) do f
        readfont(f)
    end
    return font
end

function readfont(io)
    magic = readmagic(io)

    header = split(readline(io))
    fig_header = FIGletHeader(
                           header[1][1],
                           header[2:end]...,
                          )

    for i in 1:fig_header.comment_lines
        discard = readline(io)
    end

    fig_font = FIGletFont(
                          fig_header,
                          Dict{Char, FIGletChar}(),
                          v"2.0.0",
                         )

    for c in ' ':'~'
        fig_font.font_characters[c] = readfontchar(io, c, fig_header.height)
    end

    for c in ['Ä', 'Ö', 'Ü', 'ä', 'ö', 'ü', 'ß']
        if bytesavailable(io) > 1
            fig_font.font_characters[c] = readfontchar(io, c, fig_header.height)
        end
    end

    while bytesavailable(io) > 1
        s = readline(io)
        strip(s) == "" && continue
        s = split(s)[1]
        c = if '-' in s
            Char(-(parse(UInt16, strip(s, '-'))))
        else
            Char(parse(Int, s))
        end
        fig_font.font_characters[c] = readfontchar(io, c, fig_header.height)
    end

    return fig_font
end

function availablefonts(substring)
    fonts = String[]
    for (root, dirs, files) in walkdir(FONTSDIR)
        for file in files
            if !(file in UNPARSEABLES)
                if occursin(lowercase(substring), lowercase(file)) || substring == ""
                    push!(fonts, replace(file, ".flf"=>""))
                end
            end
        end
    end
    sort!(fonts)
    return fonts
end

"""
    availablefonts() -> Vector{String}
    availablefonts(substring::AbstractString) -> Vector{String}

Returns all available fonts.
If `substring` is passed, returns available fonts that contain the case insensitive `substring`.

Example:

    julia> availablefonts()
    680-element Array{String,1}:
     "1943____"
     "1row"
     ⋮
     "zig_zag_"
     "zone7___"

    julia> FIGlet.availablefonts("3d")
    5-element Array{String,1}:
     "3D Diagonal"
     "3D-ASCII"
     "3d"
     "Henry 3D"
     "Larry 3D"

    julia>
"""
availablefonts() = availablefonts("")

function render(io, text::AbstractString, font::FIGletFont)

    iob = IOBuffer()

    fig_chars = FIGletChar[]

    for c in text
        fc = font.font_characters[c]
        push!(fig_chars, fc)
    end

    for r in 1:font.header.height
        for fc in fig_chars
            print(iob, join(fc.thechar[r, :]))
        end
        print(iob, '\n')
    end

    print(io, replace(String(take!(iob)), font.header.hardblank=>' '))
end

render(io, text::AbstractString, font::AbstractString) = render(io, text, readfont(font))

"""
    render(text::AbstractString, font::Union{AbstractString, FIGletFont})

Renders `text` using `font` to `stdout`

Example:

    render("hello world", "standard")
    render("hello world", readfont("standard"))
"""
render(text::AbstractString, font) = render(stdout, text, font)

end # module

max_stars_repo_name: wookay/TestFIGlet

max_issues_repo_name: wookay/TestFIGlet

max_forks_repo_name: wookay/TestFIGlet


text:
# Example based on https://github.com/LLNL/sundials/blob/master/examples/arkode/C_serial/ark_twowaycouple_mri.c
#=
/* ----------------------------------------------------------------
 * Programmer(s): David J. Gardner @ LLNL
 * ----------------------------------------------------------------
 * Based a linear example program by Rujeko Chinomona @ SMU.
 * ----------------------------------------------------------------
 * SUNDIALS Copyright Start
 * Copyright (c) 2002-2020, Lawrence Livermore National Security
 * and Southern Methodist University.
 * All rights reserved.
 *
 * See the top-level LICENSE and NOTICE files for details.
 *
 * SPDX-License-Identifier: BSD-3-Clause
 * SUNDIALS Copyright End
 * ----------------------------------------------------------------
 * Example problem:
 *
 * This example simulates an ODE system with 3 components,
 * Y = [u,v,w], given by the equations,
 *
 *   du/dt =  100v+w
 *   dv/dt = -100u
 *   dw/dt = -w+u
 *
 * for t in the interval [0.0, 2.0] with intial conditions
 * u(0)=9001/10001, v(0)=-1e-5/10001, and w(0)=1000. In this problem
 * the slow (w) and fast (u and v) components depend on one another.
 *
 * This program solves the problem with the MRI stepper. Outputs are
 * printed at equal intervals of 0.1 and run statistics are printed
 * at the end.
 * ----------------------------------------------------------------*/
=#

using Sundials# Test

function ff(t, y_nv, ydot_nv, user_data)
    y = convert(Vector, y_nv)
    ydot = convert(Vector, ydot_nv)
    ydot[1] = 100.0 * y[2]
    ydot[2] = -100.0 * y[1]
    ydot[3] = y[1]
    return Sundials.ARK_SUCCESS
end

function fs(t, y_nv, ydot_nv, user_data)
    y = convert(Vector, y_nv)
    ydot = convert(Vector, ydot_nv)
    ydot[1] = y[3]
    ydot[2] = 0.0
    ydot[3] = -y[3]
    return Sundials.ARK_SUCCESS
end

T0 = 0.0
Tf = 2.0
dTout = 0.1
Neq = 3
Nt = ceil(Tf / dTout)
hs = 0.001
hf = 0.00002
y0 = [0.90001, -9.999, 1000.0]

# Fast Integration portion
_mem_ptr = Sundials.ARKStepCreate(ff, C_NULL, T0, y0);
inner_arkode_mem = Sundials.Handle(_mem_ptr)
Sundials.@checkflag Sundials.ARKStepSetTableNum(
    inner_arkode_mem,
    -1,
    Sundials.KNOTH_WOLKE_3_3,
)
Sundials.@checkflag Sundials.ARKStepSetFixedStep(inner_arkode_mem, hf)

# Slow integrator portion
_arkode_mem_ptr = Sundials.MRIStepCreate(fs, T0, y0, inner_arkode_mem)
arkode_mem = Sundials.Handle(_arkode_mem_ptr)
Sundials.@checkflag Sundials.MRIStepSetFixedStep(arkode_mem, hs)

t = [T0]
tout = T0 + dTout
res = Dict(0 => y0)
for i in 1:Nt
    y = similar(y0)
    global retval = Sundials.MRIStepEvolve(arkode_mem, tout, y, t, Sundials.ARK_NORMAL)
    global tout += dTout
    global tout = (tout > Tf) ? Tf : tout
    res[i] = y
end

for i in 1:3
    sol_1 = [-0.927671 -8.500060 904.786828]
    sol_end = [0.547358 -0.523577 135.169441]
    @test isapprox(res[1][i], sol_1[i], atol = 1e-3)
    @test isapprox(res[Nt][i], sol_end[i], atol = 1e-3)
end


text:
module AssigmentProblem
using JuMP
using GLPK

"""
Problema de asignación. (Vanderbei LP, pg. 227). Dado un conjunto S de m personas, un conjunto D de m tareas, y para cada iS,jD un costo cij de asignar la tarea j a la persona i, el problema de asignación consiste en asignar a cada persona una y solo una tarea, de manera tal que se minimize el costo total. Formule este problema como un problema de programación lineal con enteros. Escriba una función para resolverlo, que tome como entrada una matriz de costos cij arbitraria.
"""

#https://github.com/josePereiro/AssigmentProblem.jl/

include("BruteForce.jl")
include("LinearProgramming.jl")
include("Common.jl");

function hello()
	println("hello");
end

end # module


text:
import LinearAlgebra.BLAS: asum, dot, blascopy!, nrm2, scal, scal!, gemm, gemm!, gemv, gemv!,
    syrk, symm, symm!, symv, symv!, trmm, trsm, trmv, trsv, trsv!, ger!

const SA = StridedArray

# # `syrk` sensitivity implementations.
# @explicit_intercepts(
#     syrk,
#     Tuple{Char, Char, ∇Scalar, StridedVecOrMat{<:∇Scalar}},
#     [false, false, true, true],
# )
# function ∇(::typeof(syrk), ::Type{Arg{3}}, p, Y, Ȳ,
#     uplo::Char,
#     trans::Char,
#     α::∇Scalar,
#     A::StridedVecOrMat{<:∇Scalar},
# )
#     g! = uppercase(uplo) == 'L' ? tril! : triu!
#     return sum(g!(Ȳ .* Y)) / α
# end
# function ∇(::typeof(syrk), ::Type{Arg{4}}, p, Y, Ȳ,
#     uplo::Char,
#     trans::Char,
#     α::∇Scalar,
#     A::StridedVecOrMat{<:∇Scalar},
# )
#     triȲ = uppercase(uplo) == 'L' ? tril(Ȳ) : triu(Ȳ)
#     out = gemm('N', trans, α, triȲ .+ triȲ', A)
#     return uppercase(trans) == 'N' ? out : out'
# end
# function ∇(Ā::StridedVecOrMat{T}, ::typeof(syrk), ::Type{Arg{4}}, p, Y, Ȳ,
#     uplo::Char,
#     trans::Char,
#     α::∇Scalar,
#     A::StridedVecOrMat{T},
# ) where T<:∇Scalar
#     triȲ = uppercase(uplo) == 'L' ? tril(Ȳ) : triu(Ȳ)
#     out = gemm('N', trans, α, triȲ .+ triȲ', A)
#     return broadcast!((ā, δā)->ā+δā, Ā, Ā, uppercase(trans) == 'N' ? out : out')
# end

# # `syrk` sensitivity implementations for `α=1`.
# @explicit_intercepts(
#     syrk,
#     Tuple{Char, Char, StridedVecOrMat{<:∇Scalar}},
#     [false, false, true],
# )
# ∇(::typeof(syrk), ::Type{Arg{3}}, p, Y, Ȳ,
#     uplo::Char,
#     trans::Char,
#     A::StridedVecOrMat{<:∇Scalar},
# ) = ∇(syrk, Arg{4}, p, Y, Ȳ, uplo, trans, one(eltype(A)), A)
# ∇(Ā::StridedVecOrMat{T}, ::typeof(syrk), ::Type{Arg{4}}, p, Y, Ȳ,
#     uplo::Char,
#     trans::Char,
#     A::StridedVecOrMat{T},
# ) where T<:∇Scalar = ∇(Ā, syrk, Arg{4}, p, Y, Ȳ, uplo, char, one(eltype(A)), A)

# `symm` sensitivity implementations.
@explicit_intercepts(
    symm,
    Tuple{Char, Char, T, StridedMatrix{T}, StridedMatrix{T}} where T<:∇Scalar,
    [false, false, true, true, true],
)
∇(::typeof(symm), ::Type{Arg{3}}, p, Y, Ȳ,
    side::Char,
    ul::Char,
    α::T,
    A::StridedMatrix{T},
    B::StridedVecOrMat{T},
) where T<:∇Scalar = sum(Ȳ .* Y) / α
function ∇(::typeof(symm), ::Type{Arg{4}}, p, Y, Ȳ,
    side::Char,
    ul::Char,
    α::T,
    A::StridedMatrix{T},
    B::StridedVecOrMat{T},
) where T<:∇Scalar
    tmp = uppercase(side) == 'L' ? Ȳ * B' : B'Ȳ
    g! = uppercase(ul) == 'L' ? tril! : triu!
    return α * g!(tmp + tmp' - Diagonal(tmp))
end
function ∇(Ā::StridedMatrix{T}, ::typeof(symm), ::Type{Arg{4}}, p, Y, Ȳ,
    side::Char,
    ul::Char,
    α::T,
    A::StridedMatrix{T},
    B::StridedVecOrMat{T},
) where T<:∇Scalar
    tmp = uppercase(side) == 'L' ? Ȳ * B' : B'Ȳ
    g! = uppercase(ul) == 'L' ? tril! : triu!
    return broadcast!((ā, δā)->ā + δā, Ā, Ā, α * g!(tmp + tmp' - Diagonal(tmp)))
end
∇(::typeof(symm), ::Type{Arg{5}}, p, Y, Ȳ,
    side::Char,
    ul::Char,
    α::T,
    A::StridedMatrix{T},
    B::StridedVecOrMat{T},
) where T<:∇Scalar = symm(side, ul, α, A, Ȳ)
∇(B̄::StridedMatrix{T}, ::typeof(symm), ::Type{Arg{5}}, p, Y, Ȳ,
    side::Char,
    ul::Char,
    α::T,
    A::StridedMatrix{T},
    B::StridedVecOrMat{T},
) where T<:∇Scalar = symm!(side, ul, α, A, Ȳ, 1.0, B̄)

# `symm` sensitivity implementations for `α=1`.
@explicit_intercepts(
    symm,
    Tuple{Char, Char, StridedMatrix{T}, StridedMatrix{T}} where T<:∇Scalar,
    [false, false, true, true],
)
∇(::typeof(symm), ::Type{Arg{3}}, p, Y, Ȳ,
    side::Char,
    ul::Char,
    A::StridedMatrix{T},
    B::StridedVecOrMat{T},
) where T<:∇Scalar = ∇(symm, Arg{4}, p, Y, Ȳ, side, ul, one(T), A, B)
∇(Ā::StridedMatrix{T}, ::typeof(symm), ::Type{Arg{3}}, p, Y, Ȳ,
    side::Char,
    ul::Char,
    A::StridedMatrix{T},
    B::StridedVecOrMat{T},
) where T<:∇Scalar = ∇(Ā, symm, Arg{4}, p, Y, Ȳ, side, ul, one(T), A, B)
∇(::typeof(symm), ::Type{Arg{4}}, p, Y, Ȳ,
    side::Char,
    ul::Char,
    A::StridedMatrix{T},
    B::StridedVecOrMat{T},
) where T<:∇Scalar = ∇(symm, Arg{5}, p, Y, Ȳ, side, ul, one(T), A, B)
∇(B̄::StridedMatrix{T}, ::typeof(symm), ::Type{Arg{4}}, p, Y, Ȳ,
    side::Char,
    ul::Char,
    A::StridedMatrix{T},
    B::StridedVecOrMat{T},
) where T<:∇Scalar = ∇(B̄, symm, Arg{5}, p, Y, Ȳ, side, ul, one(T), A, B)

# `symv` sensitivity implementations.
@explicit_intercepts(
    symv,
    Tuple{Char, T, StridedMatrix{T}, StridedVector{T}} where T<:∇Scalar,
    [false, true, true, true],
)
∇(::typeof(symv), ::Type{Arg{2}}, p, y, ȳ,
    ul::Char,
    α::T,
    A::StridedMatrix{T},
    x::StridedVector{T},
) where T<:∇Scalar = dot(ȳ, y) / α
∇(::typeof(symv), ::Type{Arg{3}}, p, y, ȳ,
    ul::Char,
    α::T,
    A::StridedMatrix{T},
    x::StridedVector{T},
) where T<:∇Scalar = ∇(symm, Arg{4}, p, y, ȳ, 'L', ul, α, A, x)
∇(Ā::StridedMatrix{T}, ::typeof(symv), ::Type{Arg{3}}, p, y, ȳ,
    ul::Char,
    α::T,
    A::StridedMatrix{T},
    x::StridedVector{T},
) where T<:∇Scalar = ∇(Ā, symm, Arg{4}, p, y, ȳ, 'L', ul, α, A, x)
∇(::typeof(symv), ::Type{Arg{4}}, p, y, ȳ,
    ul::Char,
    α::T,
    A::StridedMatrix{T},
    x::StridedVector{T},
) where T<:∇Scalar = symv(ul, α, A, ȳ)
∇(x̄::StridedVector{T}, ::typeof(symv), ::Type{Arg{4}}, p, y, ȳ,
    ul::Char,
    α::T,
    A::StridedMatrix{T},
    x::StridedVector{T},
) where T<:∇Scalar = symv!(ul, α, A, ȳ, one(T), x̄)

# `symv` sensitivity implementations for `α=1`.
@explicit_intercepts(
    symv,
    Tuple{Char, StridedMatrix{T}, StridedVector{T}} where T<:∇Scalar,
    [false, true, true],
)
∇(::typeof(symv), ::Type{Arg{2}}, p, y, ȳ,
    ul::Char,
    A::StridedMatrix{T},
    x::StridedVector{T},
) where T<:∇Scalar = ∇(symv, Arg{3}, p, y, ȳ, ul, one(T), A, x)
∇(Ā::StridedMatrix{T}, ::typeof(symv), ::Type{Arg{2}}, p, y, ȳ,
    ul::Char,
    A::StridedMatrix{T},
    x::StridedVector{T},
) where T<:∇Scalar = ∇(Ā, symv, Arg{3}, p, y, ȳ, ul, one(T), A, x)
∇(::typeof(symv), ::Type{Arg{3}}, p, y, ȳ,
    ul::Char,
    A::StridedMatrix{T},
    x::StridedVector{T},
) where T<:∇Scalar = ∇(symv, Arg{4}, p, y, ȳ, ul, one(T), A, x)
∇(B̄::StridedVector{T}, ::typeof(symv), ::Type{Arg{3}}, p, y, ȳ,
    ul::Char,
    A::StridedMatrix{T},
    x::StridedVector{T},
) where T<:∇Scalar = ∇(B̄, symv, Arg{4}, p, y, ȳ, ul, one(T), A, x)

# `trmm` sensitivity implementations.
@explicit_intercepts(
    trmm,
    Tuple{Char, Char, Char, Char, T, StridedMatrix{T}, StridedMatrix{T}} where T<:∇Scalar,
    [false, false, false, false, true, true, true],
)
∇(::typeof(trmm), ::Type{Arg{5}}, p, Y, Ȳ,
    side::Char, ul::Char, ta::Char, dA::Char,
    α::T,
    A::StridedMatrix{T},
    B::StridedVecOrMat{T},
) where T<:∇Scalar = sum(Ȳ .* Y) / α
function ∇(::typeof(trmm), ::Type{Arg{6}}, p, Y, Ȳ,
    side::Char, ul::Char, ta::Char, dA::Char,
    α::T,
    A::StridedMatrix{T},
    B::StridedVecOrMat{T},
) where T<:∇Scalar
    Ā_full = uppercase(side) == 'L' ?
        uppercase(ta) == 'N' ?
            gemm('N', 'T', α, Ȳ, B) :
            gemm('N', 'T', α, B, Ȳ) :
        uppercase(ta) == 'N' ?
            gemm('T', 'N', α, B, Ȳ) :
            gemm('T', 'N', α, Ȳ, B)
    dA == 'U' && fill!(view(Ā_full, diagind(Ā_full)), zero(T))
    return (uppercase(ul) == 'L' ? tril! : triu!)(Ā_full)
end
∇(::typeof(trmm), ::Type{Arg{7}}, p, Y, Ȳ,
    side::Char, ul::Char, ta::Char, dA::Char,
    α::T,
    A::StridedMatrix{T},
    B::StridedVecOrMat{T},
) where T<:∇Scalar =
    uppercase(side) == 'L' ?
        uppercase(ta) == 'N' ?
            trmm('L', ul, 'T', dA, α, A, Ȳ) :
            trmm('L', ul, 'N', dA, α, A, Ȳ) :
        uppercase(ta) == 'N' ?
            trmm('R', ul, 'T', dA, α, A, Ȳ) :
            trmm('R', ul, 'N', dA, α, A, Ȳ)

# `trmv` sensitivity implementations.
@explicit_intercepts(
    trmv,
    Tuple{Char, Char, Char, StridedMatrix{T}, StridedVector{T}} where T<:∇Scalar,
    [false, false, false, true, true],
)
function ∇(::typeof(trmv), ::Type{Arg{4}}, p, y, ȳ,
    ul::Char, ta::Char, dA::Char,
    A::StridedMatrix{T},
    b::StridedVector{T},
) where T<:∇Scalar
    Ā = (uppercase(ul) == 'L' ? tril! : triu!)(uppercase(ta) == 'N' ? ȳ * b' : b * ȳ')
    dA == 'U' && fill!(view(Ā, diagind(Ā)), zero(T))
    return Ā
end
∇(::typeof(trmv), ::Type{Arg{5}}, p, y, ȳ,
    ul::Char, ta::Char, dA::Char,
    A::StridedMatrix{T},
    b::StridedVector{T},
) where T<:∇Scalar = trmv(ul, uppercase(ta) == 'N' ? 'T' : 'N', dA, A, ȳ)

# `trsm` sensitivity implementations.
@explicit_intercepts(
    trsm,
    Tuple{Char, Char, Char, Char, T, StridedMatrix{T}, StridedMatrix{T}} where T<:∇Scalar,
    [false, false, false, false, true, true, true],
)
∇(::typeof(trsm), ::Type{Arg{5}}, p, Y, Ȳ,
    side::Char, ul::Char, ta::Char, dA::Char,
    α::T,
    A::StridedMatrix{T},
    X::StridedMatrix{T},
) where T<:∇Scalar = sum(Ȳ .* Y) / α
function ∇(::typeof(trsm), ::Type{Arg{6}}, p, Y, Ȳ,
    side::Char, ul::Char, ta::Char, dA::Char,
    α::T,
    A::StridedMatrix{T},
    X::StridedVecOrMat{T},
) where T<:∇Scalar
    Ā_full = uppercase(side) == 'L' ?
        uppercase(ta) == 'N' ?
            trsm('L', ul, 'T', dA, -1.0, A, Ȳ * Y') :
            trsm('R', ul, 'T', dA, -1.0, A, Y * Ȳ') :
        uppercase(ta) == 'N' ?
            trsm('R', ul, 'T', dA, -1.0, A, Y'Ȳ) :
            trsm('L', ul, 'T', dA, -1.0, A, Ȳ'Y)
    dA == 'U' && fill!(view(Ā_full, diagind(Ā_full)), zero(T))
    return (uppercase(ul) == 'L' ? tril! : triu!)(Ā_full)
end
∇(::typeof(trsm), ::Type{Arg{7}}, p, Y, Ȳ,
    side::Char, ul::Char, ta::Char, dA::Char,
    α::T,
    A::StridedMatrix{T},
    X::StridedMatrix{T},
) where T<:∇Scalar =
    uppercase(side) == 'L' ?
        uppercase(ta) == 'N' ?
            trsm('L', ul, 'T', dA, α, A, Ȳ) :
            trsm('L', ul, 'N', dA, α, A, Ȳ) :
        uppercase(ta) == 'N' ?
            trsm('R', ul, 'T', dA, α, A, Ȳ) :
            trsm('R', ul, 'N', dA, α, A, Ȳ)

# `trsv` sensitivity implementations.
@explicit_intercepts(
    trsv,
    Tuple{Char, Char, Char, StridedMatrix{T}, StridedVector{T}} where T<:∇Scalar,
    [false, false, false, true, true],
)
function ∇(::typeof(trsv), ::Type{Arg{4}}, p, y, ȳ,
    ul::Char, ta::Char, dA::Char,
    A::StridedMatrix{T},
    x::StridedVector{T},
) where T<:∇Scalar
    Ā = ∇(trsm, Arg{6}, p, y, ȳ, 'L', ul, ta, dA, one(T), A, x)
    dA == 'U' && fill!(view(Ā, diagind(Ā)), zero(T))
    return Ā
end
∇(::typeof(trsv), ::Type{Arg{5}}, p, y, ȳ,
    ul::Char, ta::Char, dA::Char,
    A::StridedMatrix{T},
    x::StridedVector{T},
) where T<:∇Scalar = trsv(ul, uppercase(ta) == 'N' ? 'T' : 'N', dA, A, ȳ)

# # # TODO: Banded matrix operations.
# # # gbmv
# # # sbmv


text:
using CmdStan, DiffEqBayes, OrdinaryDiffEq, ParameterizedFunctions,
      ModelingToolkit, RecursiveArrayTools, Distributions, Random, Test

# Uncomment for local testing only, make sure MCMCChains and StatsPlots are available
using MCMCChains, StatsPlots

#Random.seed!(123)
cd(@__DIR__)
isdir("tmp") && rm("tmp", recursive=true)

include("../../stan_inference.jl")

println("\nFour parameter case\n")
f1 = @ode_def begin
  dx = a*x - b*x*y
  dy = -c*y + d*x*y
end a b c d
u0 = [1.0,1.0]
tspan = (0.0,10.0)
p = [1.5,1.0,3.0,1.0]
prob1 = ODEProblem(f1,u0,tspan,p)
sol = solve(prob1,Tsit5())
t = collect(range(1,stop=10,length=10))
randomized = VectorOfArray([(sol(t[i]) + 0.5randn(2)) for i in 1:length(t)])
data = convert(Array,randomized)
priors = [truncated(Normal(1.0,1),0.1,2),truncated(Normal(1.5,0.5),0.1,1.5),
          truncated(Normal(2.0,1),0.1,4),truncated(Normal(1.3,0.5),0.1,2)]

bayesian_result = stan_inference(prob1,t,data,priors;num_samples=2000, nchains=4,
  num_warmup=1000,vars =(DiffEqBayes.StanODEData(),InverseGamma(4,1)))
sdf  = CmdStan.read_summary(bayesian_result.model)
@test sdf[sdf.parameters .== :theta1, :mean][1] ≈ 1.5 atol=5e-1
@test sdf[sdf.parameters .== :theta2, :mean][1] ≈ 1.0 atol=5e-1
@test sdf[sdf.parameters .== :theta3, :mean][1] ≈ 3.0 atol=5e-1
@test sdf[sdf.parameters .== :theta4, :mean][1] ≈ 1.0 atol=5e-1

# Uncomment for local chain inspection
chn = CmdStan.convert_a3d(bayesian_result.chains, bayesian_result.cnames, Val(:mcmcchains))
plot(chn)
savefig("$(@__DIR__)/four_par_case.png")

v1 = [sdf[sdf.parameters .== Symbol("u_hat[$i,1]"), :mean][1] for i in 1:10]
v2 = [sdf[sdf.parameters .== Symbol("u_hat[$i,2]"), :mean][1] for i in 1:10]
qa_prey = zeros(10, 3);
qa_pred = zeros(10, 3);
achn = Array(chn);
for i in 1:10
  qa_prey[i, :] = quantile(achn[:, 10+i], [0.055, 0.5, 0.945])
  qa_pred[i, :] = quantile(achn[:, 20+i], [0.055, 0.5, 0.945])
end
p1 = plot(1:10, v1; ribbon=(qa_prey[:, 1], qa_prey[:, 3]), color=:lightgrey, leg=false)
title!("Prey u_hat 89% quantiles")
plot!(v1, lab="prey", xlab="time", ylab="prey", color=:darkred)

p2 = plot(1:10, v2; ribbon=(qa_pred[:, 1], qa_pred[:, 3]), color=:lightgrey, leg=false)
title!("Preditors u_hat 89% quantiles")
plot!(v2, lab="pred", xlab="time", ylab="preditors", color=:darkblue)

plot(p1, p2, layout=(2,1))
savefig("$(@__DIR__)/four_par_pred_prey.png")


text:
struct Parameters
	N                        ::Int
	d                        ::Int
	m                        ::Float64
	L                        ::Float64
	T                        ::Float64
	τ                        ::Float64
	ε_LJ                     ::Float64
	σ_LJ                     ::Float64
	k_B                      ::Float64
	seed                     ::Int
	n_iter                   ::Int
	ensemble_type            ::Symbol
	rescale_velocity_interval::Int
	rng                      ::AbstractRNG
	diag_indices             ::Vector{CartesianIndex{3}}

	function Parameters(params_dict)
		new(
			params_dict[:N],
			params_dict[:d],
			params_dict[:m],
			params_dict[:L],
			params_dict[:T],
			params_dict[:τ],
			params_dict[:ε_LJ],
			params_dict[:σ_LJ],
			params_dict[:k_B],
			params_dict[:seed],
			params_dict[:n_iter],
			params_dict[:ensemble_type],
			params_dict[:rescale_velocity_interval],
			MersenneTwister(params_dict[:seed]),
			[CartesianIndex((i, 1, i)) for i in 1:params_dict[:N]]
		)
	end
end

mutable struct State
	x                        ::Matrix{Float64}
	v                        ::Matrix{Float64}
	a                        ::Matrix{Float64}
	a_prev                   ::Matrix{Float64}
	r_nearest_image          ::Array{Float64, 3}
	r_nearest_image_magnitude::Array{Float64, 3}
	E_kin                    ::Float64
	E_pot                    ::Float64

	function State(P::Parameters, params_dict::Dict)
		x = get(params_dict, :x, initialize_random_positions(P))
		v = get(params_dict, :v, initialize_random_velocities(P))
		a                         = zeros((P.N, P.d))
		a_prev                    = zeros((P.N, P.d))
		r_nearest_image           = zeros((P.N, P.d, P.N))
		r_nearest_image_magnitude = zeros((P.N, 1, P.N))
		E_kin                     = 0
		E_pot                     = 0

		new(
			x,
			v,
			a,
			a_prev,
			r_nearest_image,
			r_nearest_image_magnitude,
			E_kin,
			E_pot,
		)
	end
end

mutable struct History
	x    ::Array{Float64, 3}
	v    ::Array{Float64, 3}
	E_kin::Vector{Float64}
	E_pot::Vector{Float64}

	function History(P::Parameters)
		new(
			zeros((P.N, P.d, P.n_iter + 1)),
			zeros((P.N, P.d, P.n_iter + 1)),
			zeros(P.n_iter + 1),
			zeros(P.n_iter + 1),
		)
	end
end

mutable struct MDSimulation
	params  ::Parameters
	state   ::State
	history ::History

	function MDSimulation(params_dict)
		params  = Parameters(params_dict)
		state   = State(params, params_dict)
		history = History(params)

		mdsim = new(params, state, history)
		compute_accelerations!(mdsim)
		compute_energies!(mdsim)
		update_history!(mdsim, 1)

		return mdsim
	end
end


text:
source = ARGS[1]
plugin = ARGS[2]

include(source)

mod = eval(Meta.parse(plugin))
println(mod.entryPoint())

max_stars_repo_name: ELIFE-ASU/dglmoore

max_issues_repo_name: ELIFE-ASU/dglmoore

max_forks_repo_name: ELIFE-ASU/dglmoore


text:
"""
This module reimplements models from the paper:

    Spoerer, C.J., McClure, P. and Kriegeskorte, N. (2017).
    Recurrent convolutional neural networks: a better model of biological object recognition.
    Frontiers in Psychology 8, 1551.
"""
module Spoerer2017
using Flux
using ...Splitters
using ...Mergers
using ...FeedbackChains
using ..LRNs
using ..Flatten

export spoerer_model_b, spoerer_model_bf, spoerer_model_bk,
       spoerer_model_bl, spoerer_model_bt, spoerer_model_blt

"""
    spoerer_model_fw(T; channels=1, inputsize=(28, 28), kernel=(3,3), features=32, classes=10)

Generate one of the forward models (B, B-K, B-F) from the paper:

    Spoerer, C.J., McClure, P. & Kriegeskorte, N. (2017).
    Recurrent convolutional neural networks: a better model of biological object recognition.
    Frontiers in Psychology 8, 1551.
"""
function spoerer_model_fw(T; channels=1, inputsize=(28, 28), kernel=(3,3), features=32, classes=10)
    return Chain(
        Conv(kernel, channels=>features, relu, pad=map(x -> x ÷ 2, kernel)),
        LRN(T(1.0), T(0.0001), T(0.5), 5),
        MaxPool((2,2), stride=(2,2)),
        Conv(kernel, features=>features, relu, pad=map(x -> x ÷ 2, kernel)),
        LRN(T(1.0), T(0.0001), T(0.5), 5),
        MaxPool(map(x -> x ÷ 2, inputsize)),
        flatten,
        Dense(features, classes, σ)
    )
end # function spoerer_model_fw

"""
    spoerer_model_b(T; channels=1, inputsize=(28, 28), classes=10)

Generate the bottom-up (B) convolutional neural network from:

    Spoerer, C.J., McClure, P. & Kriegeskorte, N. (2017).
    Recurrent convolutional neural networks: a better model of biological object recognition.
    Frontiers in Psychology 8, 1551.
"""
spoerer_model_b(T; channels=1, inputsize=(28,28), classes=10) =
    spoerer_model_fw(T, channels=channels, inputsize=inputsize, classes=classes)

"""
    spoerer_model_bk(T; channels=1, inputsize=(28, 28), classes=10)

Generate the convolutional neural network with increased kernel size (BK) from:

    Spoerer, C.J., McClure, P. & Kriegeskorte, N. (2017).
    Recurrent convolutional neural networks: a better model of biological object recognition.
    Frontiers in Psychology 8, 1551.
"""
spoerer_model_bk(T; channels=1, inputsize=(28,28), classes=10) =
    spoerer_model_fw(T, channels=channels, inputsize=inputsize, classes=classes, kernel=(5,5))

"""
    spoerer_model_bf(T; channels=1, inputsize=(28, 28), classes=10)

Generate the convolutional neural network with additional feature maps (BF) from:

    Spoerer, C.J., McClure, P. & Kriegeskorte, N. (2017).
    Recurrent convolutional neural networks: a better model of biological object recognition.
    Frontiers in Psychology 8, 1551.
"""
spoerer_model_bf(T; channels=1, inputsize=(28,28), classes=10) =
    spoerer_model_fw(T, channels=channels, inputsize=inputsize, classes=classes, features=64)

# TODO: currently, the forward, lateral and backward convolutions all have their
#       own biases. This should not make the model more powerful, as they are
#       combined additively before the non-linearity, but it wastes some resources.
"""
    spoerer_model_bl(T; channels=1, inputsize=(28, 28), kernel=(3,3), features=32, classes=10)

Generate the convolutional neural network with lateral recurrence (BL) from:

    Spoerer, C.J., McClure, P. & Kriegeskorte, N. (2017).
    Recurrent convolutional neural networks: a better model of biological object recognition.
    Frontiers in Psychology 8, 1551.
"""
function spoerer_model_bl(T; channels=1, inputsize=(28, 28), kernel=(3,3), features=32, classes=10)
    return FeedbackChain(
        Conv(kernel, channels=>features, pad=map(x -> x ÷ 2, kernel)),
        Merger("l1", ConvTranspose(kernel, features=>features, pad=map(x -> x ÷ 2, kernel)), +),
        x -> relu.(x),
        Splitter("l1"),
        LRN(T(1.0), T(0.0001), T(0.5), 5),
        MaxPool((2,2), stride=(2,2)),
        Conv(kernel, features=>features, pad=map(x -> x ÷ 2, kernel)),
        Merger("l2", ConvTranspose(kernel, features=>features, pad=map(x -> x ÷ 2, kernel)), +),
        x -> relu.(x),
        Splitter("l2"),
        LRN(T(1.0), T(0.0001), T(0.5), 5),
        MaxPool(map(x -> x ÷ 2, inputsize)),
        flatten,
        Dense(features, classes, σ)
    )
end # function spoerer_model_bl

"""
    spoerer_model_bt(T; channels=1, inputsize=(28, 28), kernel=(3,3), features=32, classes=10)

Generate the convolutional neural network with top-down recurrence (BT) from:

    Spoerer, C.J., McClure, P. & Kriegeskorte, N. (2017).
    Recurrent convolutional neural networks: a better model of biological object recognition.
    Frontiers in Psychology 8, 1551.
"""
function spoerer_model_bt(T; channels=1, inputsize=(28, 28), kernel=(3,3), features=32, classes=10)
    return FeedbackChain(
        Conv(kernel, channels=>features, pad=map(x -> x ÷ 2, kernel)),
        Merger("l2",  ConvTranspose((2,2), features=>features, stride=2), +),
        x -> relu.(x),
        LRN(T(1.0), T(0.0001), T(0.5), 5),
        MaxPool((2,2), stride=(2,2)),
        Conv(kernel, features=>features, relu, pad=map(x -> x ÷ 2, kernel)),
        Splitter("l2"),
        LRN(T(1.0), T(0.0001), T(0.5), 5),
        MaxPool(map(x -> x ÷ 2, inputsize)),
        flatten,
        Dense(features, classes, σ)
    )
end # function spoerer_model_bt

"""
    spoerer_model_blt(T; channels=1, inputsize=(28, 28), kernel=(3,3), features=32, classes=10)

Generate the convolutional neural network with lateral and top-down recurrence (BLT) from:

    Spoerer, C.J., McClure, P. & Kriegeskorte, N. (2017).
    Recurrent convolutional neural networks: a better model of biological object recognition.
    Frontiers in Psychology 8, 1551.
"""
function spoerer_model_blt(T; channels=1, inputsize=(28, 28), kernel=(3,3), features=32, classes=10)
    return FeedbackChain(
        Conv(kernel, channels=>features, pad=map(x -> x ÷ 2, kernel)),
        Merger("l1", ConvTranspose(kernel, features=>features, pad=map(x -> x ÷ 2, kernel)), +),
        Merger("l2", ConvTranspose((2,2), features=>features, stride=2), +),
        x -> relu.(x),
        Splitter("l1"),
        LRN(T(1.0), T(0.0001), T(0.5), 5),
        MaxPool((2,2), stride=(2,2)),
        Conv(kernel, features=>features, pad=map(x -> x ÷ 2, kernel)),
        Merger("l2", ConvTranspose(kernel, features=>features, pad=map(x -> x ÷ 2, kernel)), +),
        x -> relu.(x),
        Splitter("l2"),
        LRN(T(1.0), T(0.0001), T(0.5), 5),
        MaxPool(map(x -> x ÷ 2, inputsize)),
        flatten,
        Dense(features, classes, σ)
    )
end # function spoerer_model_blt

end # module Spoerer2017


text:
Base.@kwdef mutable struct RectangularLayout <: AbstractLayout
    pos::ImVec2
    width::Cfloat
    height::Cfloat
end

function get_pos(l::RectangularLayout)
    l.pos
end

function set_pos!(l::RectangularLayout, pos::ImVec2)
    l.pos = pos
end

function get_width(l::RectangularLayout)
    l.width
end

function set_width!(l::RectangularLayout, width::Cfloat)
    l.width = width
end

function get_height(l::RectangularLayout)
    l.height
end

function set_height!(l::RectangularLayout, height::Cfloat)
    l.height = height
end