text:
"""
    TRMM3Hourly{ST<:AbstractString, DT<:TimeType} <: TRMMDataset

Object containing information on 3-Hourly TRMM datasets to be downloaded
"""
struct TRMM3Hourly{ST<:AbstractString, DT<:TimeType} <: TRMMDataset
	npdID :: ST
	lname :: ST
	doi   :: ST
    dtbeg :: DT
    dtend :: DT
    sroot :: ST
    hroot :: ST
    fpref :: ST
    fsuff :: ST
end

"""
    TRMM3Hourly(
        ST = String,
        DT = Date;
        dtbeg :: TimeType,
        dtend :: TimeType,
        sroot :: AbstractString = homedir(),
    ) -> npd :: TRMM3Hourly{ST,DT}

Creates a `TRMM3Hourly` dataset `npd` to retrieve datasets from the Final post-processing runs for 3-Hourly output

Keyword Arguments
=================
- `dtbeg` : Date at which download / analysis of the dataset begins
- `dtend` : Date at which download / analysis of the dataset ends
- `sroot` : The directory in which the folder `trmm3hourly` will be created for data downloads, storage and analysis, default is the home directoy called by `homedir()`

The following fields in `npd` will be fixed as below:
- `npdID` : trmm3hourly
- `lname` : Final TRMM 3-Hourly
- `doi`   : 10.5067/TRMM/TMPA/3H/7
- `hroot` : https://disc2.gesdisc.eosdis.nasa.gov/opendap/TRMM_L3/TRMM_3B42.7
- `fpref` : 3B42
- `fsuff` : 7.HDF
"""
function TRMM3Hourly(
    ST = String,
    DT = Date;
    dtbeg :: TimeType,
    dtend :: TimeType,
    sroot :: AbstractString = homedir(),
)

	@info "$(modulelog()) - Setting up data structure containing information on Final TRMM 3-Hourly data to be downloaded"

    fol = joinpath(sroot,"trmm3hourly"); if !isdir(fol); mkpath(fol) end

	dtbeg = Date(year(dtbeg),month(dtbeg),1)
	dtend = Date(year(dtend),month(dtend),daysinmonth(dtend))

    return TRMM3Hourly{ST,DT}(
		"trmm3hourly", "Final TRMM 3-Hourly", "10.5067/TRMM/TMPA/3H/7",
        dtbeg, dtend,
		joinpath(sroot,"trmm3hourly"),
        "https://disc2.gesdisc.eosdis.nasa.gov/opendap/TRMM_L3/TRMM_3B42.7",
        "3B42", "7.HDF",
    )

end

"""
    TRMM3HourlyNRT(
        ST = String,
        DT = Date;
        dtbeg :: TimeType,
        dtend :: TimeType,
        sroot :: AbstractString = homedir(),
    ) -> npd :: TRMM3Hourly{ST,DT}

Creates a `TRMM3Hourly` dataset `npd` to retrieve datasets from the Near Real-Time runs for 3-Hourly output

Keyword Arguments
=================
- `dtbeg` : Date at which download / analysis of the dataset begins
- `dtend` : Date at which download / analysis of the dataset ends
- `sroot` : The directory in which the folder `trmm3hourlynrt` will be created for data downloads, storage and analysis, default is the home directoy called by `homedir()`

The following fields in `npd` will be fixed as below:
- `npdID` : trmm3hourlynrt
- `lname` : Near Real-Time TRMM 3-Hourly
- `doi`   : 10.5067/TRMM/TMPA/3H-E/7
- `hroot` : https://disc2.gesdisc.eosdis.nasa.gov/opendap/TRMM_RT/TRMM_3B42RT.7
- `fpref` : 3B42RT
- `fsuff` : 7R2.nc4
"""
function TRMM3HourlyNRT(
    ST = String,
    DT = Date;
    dtbeg :: TimeType,
    dtend :: TimeType,
    sroot :: AbstractString = homedir(),
)

	@info "$(modulelog()) - Setting up data structure containing information on Near Real-Time TRMM 3-Hourly data to be downloaded"

    fol = joinpath(sroot,"trmm3hourlynrt"); if !isdir(fol); mkpath(fol) end

	dtbeg = Date(year(dtbeg),month(dtbeg),1)
	dtend = Date(year(dtend),month(dtend),daysinmonth(dtend))

    return TRMM3Hourly{ST,DT}(
		"trmm3hourlynrt", "Near Real-Time TRMM 3-Hourly", "10.5067/TRMM/TMPA/3H-E/7",
        dtbeg, dtend,
		joinpath(sroot,"trmm3hourlynrt"),
        "https://disc2.gesdisc.eosdis.nasa.gov/opendap/TRMM_RT/TRMM_3B42RT.7",
        "3B42RT", "7R2.nc4",
    )

end

function show(io::IO, npd::TRMM3Hourly{ST,DT}) where {ST<:AbstractString, DT<:TimeType}
    print(
		io,
		"The NASA Precipitation Dataset {$ST,$DT} has the following properties:\n",
		"    Dataset ID      (npdID) : ", npd.npdID, '\n',
		"    Logging Name    (lname) : ", npd.lname, '\n',
		"    DOI URL          (doi)  : ", npd.doi,   '\n',
		"    Data Directory  (sroot) : ", npd.sroot, '\n',
		"    Date Begin      (dtbeg) : ", npd.dtbeg, '\n',
		"    Date End        (dtend) : ", npd.dtend, '\n',
		"    Timestep                : 3 Hourly\n",
        "    Data Resolution         : 0.25º\n",
        "    Data Server     (hroot) : ", npd.hroot, '\n',
	)
end

max_stars_repo_name: natgeo-wong/NASAPrecipitation

max_issues_repo_name: natgeo-wong/NASAPrecipitation

max_forks_repo_name: natgeo-wong/NASAPrecipitation


text:
using Test
using aoc2019.Fuel
using aoc2019.computer
using aoc2019.wiring
using aoc2019.password
using aoc2019.star_chart
using aoc2019.thrusters
using aoc2019.sif
using aoc2019.asteroids
using aoc2019.Nanofactory

@testset "the tyranny of the rocket equation -- part one                          " begin
    @test Fuel.total([12]) == 2
    @test Fuel.total([14]) == 2
    @test Fuel.total([1969]) == 654
    @test Fuel.total([100756]) == 33583
end

@testset "the tyranny of the rocket equation -- part two                          " begin
    @test Fuel.correct_total([14]) == 2
    @test Fuel.correct_total([1969]) == 966
    @test Fuel.correct_total([100756]) == 50346
end

@testset "1202 program alarm -- part one                                          " begin
    @test computer.run([1,0,0,0,99]) == [2,0,0,0,99]
    @test computer.run([2,3,0,3,99]) == [2,3,0,6,99]
    @test computer.run([2,4,4,5,99,0]) == [2,4,4,5,99,9801]
    @test computer.run([1,1,1,4,99,5,6,0,99]) == [30,1,1,4,2,5,6,0,99]
end

@testset "crossed wires -- part one                                               " begin
    @test wiring.manhattan_distance(
        ["R75","D30","R83","U83","L12","D49","R71","U7","L72"],
        ["U62","R66","U55","R34","D71","R55","D58","R83"]
    ) == 159

    @test wiring.manhattan_distance(
        ["R98","U47","R26","D63","R33","U87","L62","D20","R33","U53","R51"],
        ["U98","R91","D20","R16","D67","R40","U7","R15","U6","R7"]
    ) == 135
end

@testset "crossed wires -- part two                                               " begin
    @test wiring.signal_delay(
        ["R75","D30","R83","U83","L12","D49","R71","U7","L72"],
        ["U62","R66","U55","R34","D71","R55","D58","R83"]
    ) == 610

    @test wiring.signal_delay(
        ["R98","U47","R26","D63","R33","U87","L62","D20","R33","U53","R51"],
        ["U98","R91","D20","R16","D67","R40","U7","R15","U6","R7"]
    ) == 410
end

@testset "secure container -- part one                                            " begin
    @test password.is_candidate(111111) == true
    @test password.is_candidate(223450) == false
    @test password.is_candidate(123789) == false
end

@testset "secure container -- part two                                            " begin
    @test password.is_candidate(112233) && password.is_candidate_strict(112233) == true
    @test password.is_candidate(123444) && password.is_candidate_strict(123444)== false
    @test password.is_candidate(111122) && password.is_candidate_strict(1111122)== true
end

@testset "sunny with a chance of asteroids -- part one                            " begin
    @test computer.run([1002,4,3,4,33]) == [1002,4,3,4,99]
end

@testset "sunny with a chance of asteroids -- part two                            " begin
    # equals tests
    stdin, stdout = io()
    put!(stdin, 8)
    @async computer.run([3,9,8,9,10,9,4,9,99,-1,8], stdin, stdout)
    @test take!(stdout) == 1

    stdin, stdout = io()
    put!(stdin, 9)
    @async computer.run([3,9,8,9,10,9,4,9,99,-1,8], stdin, stdout)
    @test take!(stdout) == 0

    stdin, stdout = io()
    put!(stdin, 8)
    @async computer.run([3,3,1108,-1,8,3,4,3,99], stdin, stdout)
    @test take!(stdout) == 1

    stdin, stdout = io()
    put!(stdin, 9)
    @async computer.run([3,3,1108,-1,8,3,4,3,99], stdin, stdout)
    @test take!(stdout) == 0

    # less than tests
    stdin, stdout = io()
    put!(stdin, 7)
    @async computer.run([3,9,7,9,10,9,4,9,99,-1,8], stdin, stdout)
    @test take!(stdout) == 1

    stdin, stdout = io()
    put!(stdin, 8)
    @async computer.run([3,9,7,9,10,9,4,9,99,-1,8], stdin, stdout)
    @test take!(stdout) == 0

    stdin, stdout = io()
    put!(stdin, 9)
    @async computer.run([3,9,7,9,10,9,4,9,99,-1,8], stdin, stdout)
    @test take!(stdout) == 0

    stdin, stdout = io()
    put!(stdin, 7)
    @async computer.run([3,3,1107,-1,8,3,4,3,99], stdin, stdout)
    @test take!(stdout) == 1

    stdin, stdout = io()
    put!(stdin, 8)
    @async computer.run([3,3,1107,-1,8,3,4,3,99], stdin, stdout)
    @test take!(stdout) == 0

    stdin, stdout = io()
    put!(stdin, 9)
    @async computer.run([3,3,1107,-1,8,3,4,3,99], stdin, stdout)
    @test take!(stdout) == 0

    # jump tests
    stdin, stdout = io()
    put!(stdin, 0)
    @async computer.run([3,12,6,12,15,1,13,14,13,4,13,99,-1,0,1,9], stdin, stdout)
    @test take!(stdout) == 0

    stdin, stdout = io()
    put!(stdin, 99)
    @async computer.run([3,12,6,12,15,1,13,14,13,4,13,99,-1,0,1,9], stdin, stdout)
    @test take!(stdout) == 1

    stdin, stdout = io()
    put!(stdin, 0)
    @async computer.run([3,3,1105,-1,9,1101,0,0,12,4,12,99,1], stdin, stdout)
    @test take!(stdout) == 0

    # other
    stdin, stdout = io()
    put!(stdin, 7)
    @async computer.run(
        [3,21,1008,21,8,20,1005,20,22,107,8,21,20,1006,20,31,
            1106,0,36,98,0,0,1002,21,125,20,4,20,1105,1,46,104,
            999,1105,1,46,1101,1000,1,20,4,20,1105,1,46,98,99
        ],
        stdin,
        stdout
    )
    @test take!(stdout) == 999

    stdin, stdout = io()
    put!(stdin, 8)
    @async computer.run(
        [3,21,1008,21,8,20,1005,20,22,107,8,21,20,1006,20,31,
            1106,0,36,98,0,0,1002,21,125,20,4,20,1105,1,46,104,
            999,1105,1,46,1101,1000,1,20,4,20,1105,1,46,98,99
        ],
        stdin,
        stdout
    )
    @test take!(stdout) == 1000

    stdin, stdout = io()
    put!(stdin, 9)
    @async computer.run(
        [3,21,1008,21,8,20,1005,20,22,107,8,21,20,1006,20,31,
            1106,0,36,98,0,0,1002,21,125,20,4,20,1105,1,46,104,
            999,1105,1,46,1101,1000,1,20,4,20,1105,1,46,98,99
        ],
        stdin,
        stdout
    )
    @test take!(stdout) == 1001
end

@testset "universal orbit map -- part one                                         " begin
    @test star_chart.total_orbits([
        "COM)B", "B)C", "C)D", "D)E", "E)F", "B)G", "G)H", "D)I", "E)J", "J)K", "K)L"
    ]) == 42
end

@testset "universal orbit map -- part two                                         " begin
    @test star_chart.total_transfers([
        "COM)B", "B)C", "C)D", "D)E", "E)F", "B)G", "G)H", "D)I",
        "E)J", "J)K", "K)L", "K)YOU", "I)SAN"
    ]) == 4
end

@testset "amplification circuit -- part one                                       " begin
    @test thrusters.output(
        [3,15,3,16,1002,16,10,16,1,16,15,15,4,15,99,0,0],
        [4,3,2,1,0],
        0
    ) == 43210

    @test thrusters.output(
        [3,23,3,24,1002,24,10,24,1002,23,-1,23,101,5,23,23,1,24,23,23,4,23,99,0,0],
        [0,1,2,3,4],
        0
    ) == 54321

    @test thrusters.output(
        [3,31,3,32,1002,32,10,32,1001,31,-2,31,1007,31,0,33,
            1002,33,7,33,1,33,31,31,1,32,31,31,4,31,99,0,0,0
        ],
        [1,0,4,3,2],
        0
    ) == 65210
end

@testset "amplification circuit -- part two                                       " begin
    @test thrusters.feedback_loop(
        [3,26,1001,26,-4,26,3,27,1002,27,2,27,1,27,26,
            27,4,27,1001,28,-1,28,1005,28,6,99,0,0,5
        ],
        [9,8,7,6,5],
        0
    ) == 139629729

    @test thrusters.feedback_loop(
        [3,52,1001,52,-5,52,3,53,1,52,56,54,1007,54,5,55,1005,55,26,1001,
            54,-5,54,1105,1,12,1,53,54,53,1008,54,0,55,1001,55,1,55,2,53,
            55,53,4,53,1001,56,-1,56,1005,56,6,99,0,0,0,0,10
        ],
        [9,7,8,5,6],
        0
    ) == 18216
end

@testset "space image format -- part one                                          " begin
    @test sif.checksum(sif.image([1,2,3,4,5,6,7,8,9,0,1,2], 3, 2)) == 1
    @test sif.checksum(sif.image([0,0,0,0,0,1,1,1,1,1,1,2], 2, 2)) == 3
end

@testset "space image format -- part two                                          " begin
    @test sif.resolve(sif.image([0,2,2,2,1,1,2,2,2,2,1,2,0,0,0,0], 2, 2)) == [0 1; 1 0]
end

@testset "sensor boost -- part one                                                " begin
    program = [109,1,204,-1,1001,100,1,100,1008,100,16,101,1006,101,0,99]
    stdin, stdout = computer.io()
    task = @async computer.run(copy(program), stdin, stdout)
    while !istaskdone(task); yield(); end
    close(stdout)
    @test [o for o in stdout] == program

    program = [1102,34915192,34915192,7,4,7,99,0]
    stdin, stdout = computer.io()
    task = @async computer.run(program, stdin, stdout)
    while !istaskdone(task); yield(); end
    @test length(digits(take!(stdout))) == 16

    program = [104,1125899906842624,99]
    stdin, stdout = computer.io()
    task = @async computer.run(program, stdin, stdout)
    while !istaskdone(task); yield(); end
    @test take!(stdout) == 1125899906842624
end

@testset "monitoring station -- part one                                          " begin
    map = asteroids.chart([".#..#", ".....", "#####", "....#", "...##"])
    (coords, visible) = asteroids.maxima(map)
    @test coords == asteroids.Point(4, 5)
    @test visible == 8

    map = asteroids.chart([
        "......#.#.",
        "#..#.#....",
        "..#######.",
        ".#.#.###..",
        ".#..#.....",
        "..#....#.#",
        "#..#....#.",
        ".##.#..###",
        "##...#..#.",
        ".#....####"
    ])
    (coords, visible) = asteroids.maxima(map)
    @test coords == asteroids.Point(6, 9)
    @test visible == 33

    map = asteroids.chart([
        "#.#...#.#.",
        ".###....#.",
        ".#....#...",
        "##.#.#.#.#",
        "....#.#.#.",
        ".##..###.#",
        "..#...##..",
        "..##....##",
        "......#...",
        ".####.###."
    ])
    (coords, visible) = asteroids.maxima(map)
    @test coords == asteroids.Point(2, 3)
    @test visible == 35

    map = asteroids.chart([
        ".#..##.###...#######",
        "##.############..##.",
        ".#.######.########.#",
        ".###.#######.####.#.",
        "#####.##.#.##.###.##",
        "..#####..#.#########",
        "####################",
        "#.####....###.#.#.##",
        "##.#################",
        "#####.##.###..####..",
        "..######..##.#######",
        "####.##.####...##..#",
        ".#####..#.######.###",
        "##...#.##########...",
        "#.##########.#######",
        ".####.#.###.###.#.##",
        "....##.##.###..#####",
        ".#.#.###########.###",
        "#.#.#.#####.####.###",
        "###.##.####.##.#..##"
    ])
    (coords, visible) = asteroids.maxima(map)
    @test coords == asteroids.Point(12, 14)
    @test visible == 210
end

@testset "monitoring station -- part two                                          " begin
    tau = 2 * atan(0, -1)
    map = asteroids.chart([
        ".#....#####...#..",
        "##...##.#####..##",
        "##...#...#.#####.",
        "..#.....X...###..",
        "..#.#.....#....##"
    ])
    asts = asteroids.sweep(map, Point(9, 4), atan(-1, 0) + tau, atan(-1, 0))
    @test asts[1] == Point(9, 2)
    @test asts[9] == Point(16, 2)
end

@testset "space stoichiometry -- part one                                         " begin
    formulae = Nanofactory.parse_equations([
        "10 ORE => 10 A"
        "1 ORE => 1 B"
        "7 A, 1 B => 1 C"
        "7 A, 1 C => 1 D"
        "7 A, 1 D => 1 E"
        "7 A, 1 E => 1 FUEL"
    ])
    inventory = Dict("FUEL" => -1)
    Nanofactory.produce(inventory, formulae)
    ore = abs(inventory["ORE"])
    @test ore == 31

    formulae = Nanofactory.parse_equations([
        "9 ORE => 2 A"
        "8 ORE => 3 B"
        "7 ORE => 5 C"
        "3 A, 4 B => 1 AB"
        "5 B, 7 C => 1 BC"
        "4 C, 1 A => 1 CA"
        "2 AB, 3 BC, 4 CA => 1 FUEL"
    ])
    inventory = Dict("FUEL" => -1)
    Nanofactory.produce(inventory, formulae)
    ore = abs(inventory["ORE"])
    @test ore == 165

    formulae = Nanofactory.parse_equations([
        "157 ORE => 5 NZVS"
        "165 ORE => 6 DCFZ"
        "44 XJWVT, 5 KHKGT, 1 QDVJ, 29 NZVS, 9 GPVTF, 48 HKGWZ => 1 FUEL"
        "12 HKGWZ, 1 GPVTF, 8 PSHF => 9 QDVJ"
        "179 ORE => 7 PSHF"
        "177 ORE => 5 HKGWZ"
        "7 DCFZ, 7 PSHF => 2 XJWVT"
        "165 ORE => 2 GPVTF"
        "3 DCFZ, 7 NZVS, 5 HKGWZ, 10 PSHF => 8 KHKGT"
    ])
    inventory = Dict("FUEL" => -1)
    Nanofactory.produce(inventory, formulae)
    ore = abs(inventory["ORE"])
    @test ore == 13312

    formulae = Nanofactory.parse_equations([
        "2 VPVL, 7 FWMGM, 2 CXFTF, 11 MNCFX => 1 STKFG"
        "17 NVRVD, 3 JNWZP => 8 VPVL"
        "53 STKFG, 6 MNCFX, 46 VJHF, 81 HVMC, 68 CXFTF, 25 GNMV => 1 FUEL"
        "22 VJHF, 37 MNCFX => 5 FWMGM"
        "139 ORE => 4 NVRVD"
        "144 ORE => 7 JNWZP"
        "5 MNCFX, 7 RFSQX, 2 FWMGM, 2 VPVL, 19 CXFTF => 3 HVMC"
        "5 VJHF, 7 MNCFX, 9 VPVL, 37 CXFTF => 6 GNMV"
        "145 ORE => 6 MNCFX"
        "1 NVRVD => 8 CXFTF"
        "1 VJHF, 6 MNCFX => 4 RFSQX"
        "176 ORE => 6 VJHF"
    ])
    inventory = Dict("FUEL" => -1)
    Nanofactory.produce(inventory, formulae)
    ore = abs(inventory["ORE"])
    @test ore == 180697

    formulae = Nanofactory.parse_equations([
        "171 ORE => 8 CNZTR"
        "7 ZLQW, 3 BMBT, 9 XCVML, 26 XMNCP, 1 WPTQ, 2 MZWV, 1 RJRHP => 4 PLWSL"
        "114 ORE => 4 BHXH"
        "14 VRPVC => 6 BMBT"
        "6 BHXH, 18 KTJDG, 12 WPTQ, 7 PLWSL, 31 FHTLT, 37 ZDVW => 1 FUEL"
        "6 WPTQ, 2 BMBT, 8 ZLQW, 18 KTJDG, 1 XMNCP, 6 MZWV, 1 RJRHP => 6 FHTLT"
        "15 XDBXC, 2 LTCX, 1 VRPVC => 6 ZLQW"
        "13 WPTQ, 10 LTCX, 3 RJRHP, 14 XMNCP, 2 MZWV, 1 ZLQW => 1 ZDVW"
        "5 BMBT => 4 WPTQ"
        "189 ORE => 9 KTJDG"
        "1 MZWV, 17 XDBXC, 3 XCVML => 2 XMNCP"
        "12 VRPVC, 27 CNZTR => 2 XDBXC"
        "15 KTJDG, 12 BHXH => 5 XCVML"
        "3 BHXH, 2 VRPVC => 7 MZWV"
        "121 ORE => 7 VRPVC"
        "7 XCVML => 6 RJRHP"
        "5 BHXH, 4 VRPVC => 5 LTCX"
    ])
    inventory = Dict("FUEL" => -1)
    Nanofactory.produce(inventory, formulae)
    ore = abs(inventory["ORE"])
    @test ore == 2210736
end

@testset "space stoichiometry -- part two                                         " begin
    formulae = Nanofactory.parse_equations([
        "157 ORE => 5 NZVS"
        "165 ORE => 6 DCFZ"
        "44 XJWVT, 5 KHKGT, 1 QDVJ, 29 NZVS, 9 GPVTF, 48 HKGWZ => 1 FUEL"
        "12 HKGWZ, 1 GPVTF, 8 PSHF => 9 QDVJ"
        "179 ORE => 7 PSHF"
        "177 ORE => 5 HKGWZ"
        "7 DCFZ, 7 PSHF => 2 XJWVT"
        "165 ORE => 2 GPVTF"
        "3 DCFZ, 7 NZVS, 5 HKGWZ, 10 PSHF => 8 KHKGT"
    ])
    low = 1
    high = 1000000000000
    while low < high
        guess = div(low + high + 1, 2)
        inventory = Dict("FUEL" => -1 * guess)
        Nanofactory.produce(inventory, formulae)
        ore = abs(inventory["ORE"])
        if ore < 1000000000000
            low = guess
        else
            high = guess - 1
        end
    end
    @test low == 82892753

    formulae = Nanofactory.parse_equations([
        "2 VPVL, 7 FWMGM, 2 CXFTF, 11 MNCFX => 1 STKFG"
        "17 NVRVD, 3 JNWZP => 8 VPVL"
        "53 STKFG, 6 MNCFX, 46 VJHF, 81 HVMC, 68 CXFTF, 25 GNMV => 1 FUEL"
        "22 VJHF, 37 MNCFX => 5 FWMGM"
        "139 ORE => 4 NVRVD"
        "144 ORE => 7 JNWZP"
        "5 MNCFX, 7 RFSQX, 2 FWMGM, 2 VPVL, 19 CXFTF => 3 HVMC"
        "5 VJHF, 7 MNCFX, 9 VPVL, 37 CXFTF => 6 GNMV"
        "145 ORE => 6 MNCFX"
        "1 NVRVD => 8 CXFTF"
        "1 VJHF, 6 MNCFX => 4 RFSQX"
        "176 ORE => 6 VJHF"
    ])
    low = 1
    high = 1000000000000
    while low < high
        guess = div(low + high + 1, 2)
        inventory = Dict("FUEL" => -1 * guess)
        Nanofactory.produce(inventory, formulae)
        ore = abs(inventory["ORE"])
        if ore < 1000000000000
            low = guess
        else
            high = guess - 1
        end
    end
    @test low == 5586022

    formulae = Nanofactory.parse_equations([
        "171 ORE => 8 CNZTR"
        "7 ZLQW, 3 BMBT, 9 XCVML, 26 XMNCP, 1 WPTQ, 2 MZWV, 1 RJRHP => 4 PLWSL"
        "114 ORE => 4 BHXH"
        "14 VRPVC => 6 BMBT"
        "6 BHXH, 18 KTJDG, 12 WPTQ, 7 PLWSL, 31 FHTLT, 37 ZDVW => 1 FUEL"
        "6 WPTQ, 2 BMBT, 8 ZLQW, 18 KTJDG, 1 XMNCP, 6 MZWV, 1 RJRHP => 6 FHTLT"
        "15 XDBXC, 2 LTCX, 1 VRPVC => 6 ZLQW"
        "13 WPTQ, 10 LTCX, 3 RJRHP, 14 XMNCP, 2 MZWV, 1 ZLQW => 1 ZDVW"
        "5 BMBT => 4 WPTQ"
        "189 ORE => 9 KTJDG"
        "1 MZWV, 17 XDBXC, 3 XCVML => 2 XMNCP"
        "12 VRPVC, 27 CNZTR => 2 XDBXC"
        "15 KTJDG, 12 BHXH => 5 XCVML"
        "3 BHXH, 2 VRPVC => 7 MZWV"
        "121 ORE => 7 VRPVC"
        "7 XCVML => 6 RJRHP"
        "5 BHXH, 4 VRPVC => 5 LTCX"
    ])
    low = 1
    high = 1000000000000
    while low < high
        guess = div(low + high + 1, 2)
        inventory = Dict("FUEL" => -1 * guess)
        Nanofactory.produce(inventory, formulae)
        ore = abs(inventory["ORE"])
        if ore < 1000000000000
            low = guess
        else
            high = guess - 1
        end
    end
    @test low == 460664
end

ext: jl


text:
"""
    array(M::MPS; reverse::Bool = true)

Generate the full dense vector from an MPS.
"""
function array(M::MPS; reverse::Bool=true)
  # check if it is a vectorized MPO
  if length(siteinds(M, 1)) == 2
    s = []
    for j in 1:length(M)
      push!(s, firstind(M[j]; tags="Input"))
      push!(s, firstind(M[j]; tags="Output"))
    end
  else
    s = siteinds(M)
  end
  if reverse
    s = Base.reverse(s)
  end
  C = combiner(s...)
  Mvec = prod(M) * dag(C)
  return ITensors.array(Mvec)
end

"""
    array(M::MPO; reverse::Bool = true)
    array(L::LPDO; reverse::Bool = true)

Generate the full dense matrix from an MPO or LPDO.
"""
function array(M::MPO; reverse::Bool=true)
  if length(siteinds(M, 1)) == 4
    s = []
    for j in 1:length(M)
      push!(s, firstind(M[j]; tags="Input", plev=0))
      push!(s, firstind(M[j]; tags="Output", plev=0))
    end
  else
    s = firstsiteinds(M; plev=0)
  end
  if reverse
    s = Base.reverse(s)
  end
  C = combiner(s...)
  Mmat = prod(M) * dag(C) * C'
  c = combinedind(C)
  return ITensors.array(permute(Mmat, c', c))

end

function array(L::LPDO{MPO}; reverse::Bool=true)
  !ischoi(L) && return PastaQ.array(MPO(L); reverse=reverse)
  M = MPO(L)
  s = []
  for j in 1:length(M)
    push!(s, firstind(M[j]; tags="Input", plev=0))
    push!(s, firstind(M[j]; tags="Output", plev=0))
  end
  if reverse
    s = Base.reverse(s)
  end
  C = combiner(s...)
  Mmat = prod(M) * dag(C) * C'
  c = combinedind(C)
  return ITensors.array(permute(Mmat, c', c))
end

is_operator(T::ITensor) =
  !isempty(inds(T, plev = 1))

"""
    PastaQ.array(T::ITensor; kwargs...)

Transform a prod(MPS/MPO) into a dense vector/matrix
"""
PastaQ.array(T::ITensor; kwargs...) =
  (is_operator(T) ? tomatrix(T; kwargs...) : tovector(T; kwargs...))

function tovector(M::ITensor; reverse::Bool = true)
  if length(inds(M,tags="n=1")) > 1
    error("Cannot transform a density matrix into a vector")
  end
  length(inds(M)) == 1 && return ITensors.array(M)
  s = []
  for j in 1:length(inds(M,plev=0))
    push!(s,firstind(M,tags="n=$(j)",plev=0))
  end
  if reverse
    s = Base.reverse(s)
  end
  C = combiner(s...)
  return ITensors.array(M * C)
end

function tomatrix(M::ITensor; reverse::Bool = true)
  if length(inds(M,tags="n=1")) == 1
    error("Cannot transform a wavefunctionm into a matrix")
  end
  length(inds(M)) == 2 && return ITensors.array(M)
  s = []
  if (hastags(M, "Input")) && (hastags(M, "Output"))
    for j in 1:length(inds(M,plev=0))÷2
      push!(s, firstind(M,tags="Input,n=$(j)",plev=0))
      push!(s, firstind(M,tags="Output,n=$(j)",plev=0))
    end
  else
    for j in 1:length(inds(M,plev=0))
      push!(s,firstind(M,tags="n=$(j)",plev=0))
    end
  end
  if reverse
    s = Base.reverse(s)
  end
  C = combiner(s...)
  Mmat = M * dag(C) * C'
  c = combinedind(C)
  return ITensors.array(permute(Mmat, c', c))
end

function itensor(M::AbstractMatrix, sites::Vector{<:Index}; reverse::Bool = true)
  sites = reverse ? Base.reverse(sites) : sites
  return ITensors.itensor(M, sites', ITensors.dag(sites))
end

function itensor(v::AbstractVector, sites::Vector{<:Index}; reverse::Bool = true)
  sites = reverse ? Base.reverse(sites) : sites
  return ITensors.itensor(v, sites)
end


text:
@testset "6.7 Miscellaneous" begin
    include("6.7.1 Hyperbolic functions.jl")
end

max_stars_repo_path: test/6 Hyperbolic functions/6.7 Miscellaneous/.jl


text:
# fDMRG_tfIsing.jl
# Using finite DMRG to calculate the ground state energy of TransverseFieldIsing model.

using TensorLabXD
using MPSLabXD
using KrylovKit

# Step 1: Initialize the MPS
V_virtual = ℤ₂Space(0=>5,1=>5)
V_physical = ℤ₂Space(0=>1,1=>1)
V_boundary = ℤ₂Space(0=>1)
T_left_boundary = TensorMap(randn, V_boundary*V_physical, V_virtual)
T_bulk = TensorMap(randn,V_virtual*V_physical,V_virtual)
T_right_boundary = TensorMap(randn,V_virtual*V_physical,V_boundary)
L = 10
s = append!([isomorphism(V_boundary,V_boundary)],
            [isomorphism(V_virtual,V_virtual) for i in 1:L-1],
            [isomorphism(V_boundary,V_boundary)])
mps = FiniteMPS(append!([T_left_boundary],[T_bulk for i in 1:L-2],
                            [T_right_boundary]),s,[:non_canonical])
make_right_canonical(mps)

# step 2: Define the model Hamiltonian as mpo
J = 1.0
g = 0.0
mpo, final_position = TransverseFieldIsing(L,J,g)

# step 2: Do finite DMRG
DMRG_parameters = Dict("max_E_error"=>10^-8,
                       "max_S_error"=>10^-3,
                       "max_sweep"=>100,
                       "min_sweep"=>5)

truncation_parameters = Dict("chi_max"=>100,
                             "svd_min"=>10^-15,
                             "truncation_error"=>10^-15)

finiteDMRG(mps, mpo, final_position, DMRG_parameters, truncation_parameters)


text:
"""model definition of arm22"""

struct StateArm22 <: AbsState
    theta1::Float64
    theta2::Float64
end

function get_mid_status(p::StateArm22, q::StateArm22)::StateArm22
    return StateArm22(
        diff_angles(q.theta1, p.theta1) / 2 + p.theta1,
        diff_angles(q.theta2, p.theta2) / 2 + p.theta2,
    )
end

function dist(q1::StateArm22, q2::StateArm22)
    return norm([
        diff_angles(q1.theta1, q2.theta1) / π,
        diff_angles(q1.theta2, q2.theta2) / π,
    ])
end

function get_arm22_positions(q::StateArm22, pos::Vector{Float64}, rad::Float64)
    a = pos
    b = a + [cos(q.theta1), sin(q.theta1)] * rad
    c = b + [cos(q.theta2), sin(q.theta2)] * rad
    return [a, b, c]
end

function gen_connect(
    q::StateArm22,
    obstacles::Vector{CircleObstacle2D},
    positions::Vector{Vector{Float64}},
    rads::Vector{Float64};
    step_dist::Float64 = STEP_DIST,
    max_dist::Union{Nothing,Float64} = nothing,
    safety_dist::Float64 = SAFETY_DIST_LINE,
)::Function

    # check: q \in C_free
    f(q::StateArm22, i::Int64)::Bool = begin
        a, b, c = get_arm22_positions(q, positions[i], rads[i])

        any(x -> (x < 0 || 1 < x), vcat(a, b, c)) && return false

        any(
            o -> (dist(a, b, [o.x, o.y]) < o.r || dist(b, c, [o.x, o.y]) < o.r),
            obstacles,
        ) && return false

        # self collision
        dist(a, b, c) < safety_dist && return false

        return true
    end

    f(q_from::StateArm22, q_to::StateArm22, i::Int64)::Bool = begin
        D = dist(q_from, q_to)
        !isnothing(max_dist) && D > max_dist && return false

        dt1 = diff_angles(q_to.theta1, q_from.theta1)
        dt2 = diff_angles(q_to.theta2, q_from.theta2)

        a = positions[i]

        for e in vcat(collect(0:step_dist:D) / D, 1.0)
            t1 = q_from.theta1 + e * dt1
            t2 = q_from.theta2 + e * dt2

            b = a + rads[i] * [cos(t1), sin(t1)]
            c = b + rads[i] * [cos(t2), sin(t2)]

            # outside
            any(x -> (x < 0 || 1 < x), vcat(a, b, c)) && return false

            # obstacles
            any(
                o -> (dist(a, b, [o.x, o.y]) < o.r || dist(b, c, [o.x, o.y]) < o.r),
                obstacles,
            ) && return false

            # self collision
            dist(a, b, c) < safety_dist && return false
        end

        return true
    end
    return f
end

function gen_collide(
    q::StateArm22,
    positions::Vector{Vector{Float64}},
    rads::Vector{Float64};
    step_dist::Float64 = STEP_DIST,
    safety_dist::Float64 = SAFETY_DIST_LINE,
)::Function

    N = length(rads)

    f(
        q_i_from::StateArm22,
        q_i_to::StateArm22,
        q_j_from::StateArm22,
        q_j_to::StateArm22,
        i::Int64,
        j::Int64,
    ) = begin
        # check each pair of step
        D_i = dist(q_i_from, q_i_to)
        D_j = dist(q_j_from, q_j_to)

        dt1_i = diff_angles(q_i_to.theta1, q_i_from.theta1)
        dt2_i = diff_angles(q_i_to.theta2, q_i_from.theta2)
        dt1_j = diff_angles(q_j_to.theta1, q_j_from.theta1)
        dt2_j = diff_angles(q_j_to.theta2, q_j_from.theta2)

        for e_i in vcat(collect(0:step_dist:D_i) / D_i, 1.0)
            # intermediate angles & positions for agent-i
            t1_i = q_i_from.theta1 + e_i * dt1_i
            t2_i = q_i_from.theta2 + e_i * dt2_i
            a_i = positions[i]
            b_i = rads[i] * [cos(t1_i), sin(t1_i)] + a_i
            c_i = rads[i] * [cos(t2_i), sin(t2_i)] + b_i

            for e_j in vcat(collect(0:step_dist:D_j) / D_j, 1.0)
                # intermediate angles & positions for agent-j
                t1_j = q_j_from.theta1 + e_j * dt1_j
                t2_j = q_j_from.theta2 + e_j * dt2_j
                a_j = positions[j]
                b_j = rads[j] * [cos(t1_j), sin(t1_j)] + a_j
                c_j = rads[j] * [cos(t2_j), sin(t2_j)] + b_j

                # check_collision
                any(
                    e -> dist(e...) < safety_dist,
                    [
                        [a_i, b_i, a_j, b_j],
                        [a_i, b_i, b_j, c_j],
                        [b_i, c_i, a_j, b_j],
                        [b_i, c_i, b_j, c_j],
                    ],
                ) && return true
            end
        end
        return false
    end

    f(q_i::StateArm22, q_j::StateArm22, i::Int64, j::Int64) = begin
        a_i, b_i, c_i = get_arm22_positions(q_i, positions[i], rads[i])
        a_j, b_j, c_j = get_arm22_positions(q_j, positions[j], rads[j])

        # check_collision
        return any(
            e -> dist(e...) < safety_dist,
            [
                [a_i, b_i, a_j, b_j],
                [a_i, b_i, b_j, c_j],
                [b_i, c_i, a_j, b_j],
                [b_i, c_i, b_j, c_j],
            ],
        )
    end

    f(Q_from::Vector{Node{StateArm22}}, Q_to::Vector{Node{StateArm22}}) = begin
        for i = 1:N, j = i+1:N
            if f(Q_from[i].q, Q_to[i].q, Q_from[j].q, Q_to[j].q, i, j)
                return true
            end
        end
        return false
    end

    f(C::Vector{StateArm22}, q_i_to::StateArm22, i::Int64) = begin
        q_i_from = C[i]
        D_i = dist(q_i_from, q_i_to)

        dt1_i = diff_angles(q_i_to.theta1, q_i_from.theta1)
        dt2_i = diff_angles(q_i_to.theta2, q_i_from.theta2)

        for e_i in vcat(collect(0:step_dist:D_i) / D_i, 1.0)
            # intermediate angles & positions for agent-i
            t1_i = q_i_from.theta1 + e_i * dt1_i
            t2_i = q_i_from.theta2 + e_i * dt2_i
            a_i = positions[i]
            b_i = rads[i] * [cos(t1_i), sin(t1_i)] + a_i
            c_i = rads[i] * [cos(t2_i), sin(t2_i)] + b_i

            for j in filter(k -> k != i, 1:N)
                a_j, b_j, c_j = get_arm22_positions(C[j], positions[j], rads[j])

                # check collision
                any(
                    e -> dist(e...) < safety_dist,
                    [
                        [a_i, b_i, a_j, b_j],
                        [a_i, b_i, b_j, c_j],
                        [b_i, c_i, a_j, b_j],
                        [b_i, c_i, b_j, c_j],
                    ],
                ) && return true
            end
        end

        return false
    end

    f(Q::Vector{Node{StateArm22}}, q_i_to::StateArm22, i::Int64) = begin
        return f(map(v -> v.q, Q), q_i_to, i)
    end

    return f
end

function gen_uniform_sampling(q::StateArm22)::Function
    () -> begin
        return StateArm22((rand(2) * 2π)...)
    end
end

function plot_motion!(
    q_from::StateArm22,
    q_to::StateArm22,
    pos::Vector{Float64},
    rad::Float64,
    params,
)
    p1 = get_arm22_positions(q_from, pos, rad)[end]
    p2 = get_arm22_positions(q_to, pos, rad)[end]
    plot!([p1[1], p2[1]], [p1[2], p2[2]]; params...)
end

function plot_start_goal!(
    q_init::StateArm22,
    q_goal::StateArm22,
    pos::Vector{Float64},
    rad::Float64,
    params,
)
    pos_init = get_arm22_positions(q_init, pos, rad)[end]
    pos_goal = get_arm22_positions(q_goal, pos, rad)[end]
    plot!([pos_init[1]], [pos_init[2]]; markershape = :hex, params...)
    plot!([pos_goal[1]], [pos_goal[2]]; markershape = :star, params...)
end

function plot_agent!(q::StateArm22, pos::Vector{Float64}, rad::Float64, color::String)
    _, pos1, pos2 = get_arm22_positions(q, pos, rad)
    plot!(
        [pos[1], pos1[1], pos2[1]],
        [pos[2], pos1[2], pos2[2]],
        color = color,
        lw = 5,
        markershape = :circle,
        markersize = 3,
        label = nothing,
    )
    scatter!(
        [pos[1]],
        [pos[2]],
        markersize = 8,
        markershape = :rect,
        color = color,
        label = nothing,
    )
end

function gen_random_instance_StateArm22(;
    N_min::Int64 = 2,
    N_max::Int64 = 8,
    N::Int64 = rand(N_min:N_max),
    num_obs_min::Int64 = 0,
    num_obs_max::Int64 = 10,
    num_obs::Int64 = rand(num_obs_min:num_obs_max),
    rad::Float64 = 0.025,
    rad_obs::Float64 = 0.05,
    rad_min::Float64 = rad,
    rad_max::Float64 = rad,
    rad_obs_min::Float64 = rad_obs,
    rad_obs_max::Float64 = rad_obs,
    TIME_LIMIT::Float64 = 0.5,
)::Tuple{
    Vector{StateArm22},  # initial configuration
    Vector{StateArm22},  # goal configuration
    Vector{CircleObstacle2D},
    Vector{Vector{Float64}},  # positions
    Vector{Float64},  # radius
}

    while true
        t_s = now()
        timeover() = elapsed_sec(t_s) > TIME_LIMIT

        # generate obstacles
        obstacles = gen_obstacles(2, num_obs, rad_obs_min, rad_obs_max)

        # determine positions
        positions = map(e -> rand(2), 1:N)

        # determine rads
        rads = map(e -> rand() * (rad_max - rad_min) + rad_min, 1:N)

        # generate starts & goals
        _q = StateArm22(0, 0)
        connect = gen_connect(_q, obstacles, positions, rads)
        collide = gen_collide(_q, positions, rads)
        config_init, config_goal = gen_config_init_goal(_q, N, connect, collide, timeover)

        timeover() && continue
        return (config_init, config_goal, obstacles, positions, rads)
    end
end


text:
# abstract type
abstract type  AdaptationAlgorithm end

# typs
"""
    FixedKernel

Fixed proposal distribution.

Parameters:

* `Cov::Array{Real}` covaraince matrix
"""
type FixedKernel <: AdaptationAlgorithm
  Cov::Array{Real}
end

"""
    AMUpdate

Adaptive tuning of the proposal distribution. Sources: *A tutorial on adaptive MCMC* [https://link.springer.com/article/10.1007/s11222-008-9110-y](https://link.springer.com/article/10.1007/s11222-008-9110-y),
and *Exploring the common concepts of adaptive MCMC and Covariance Matrix Adaptation schemes* [http://drops.dagstuhl.de/opus/volltexte/2010/2813/pdf/10361.MuellerChristian.Paper.2813.pdf](http://drops.dagstuhl.de/opus/volltexte/2010/2813/pdf/10361.MuellerChristian.Paper.2813.pdf)

Parameters:

* `C_0::Array{Real}`
* `r_cov_m::Real`
* `gamma_0::Real`
* `k::Real`
* `t_0::Integer`
"""
type AMUpdate <: AdaptationAlgorithm
  C_0::Array{Real}
  r_cov_m::Real
  gamma_0::Real
  k::Real
  t_0::Integer
end

"""
    AMUpdate

Adaptive tuning of the proposal distribution. Sources: *A tutorial on adaptive MCMC* [https://link.springer.com/article/10.1007/s11222-008-9110-y](https://link.springer.com/article/10.1007/s11222-008-9110-y),
and *Exploring the common concepts of adaptive MCMC and Covariance Matrix Adaptation schemes* [http://drops.dagstuhl.de/opus/volltexte/2010/2813/pdf/10361.MuellerChristian.Paper.2813.pdf](http://drops.dagstuhl.de/opus/volltexte/2010/2813/pdf/10361.MuellerChristian.Paper.2813.pdf)

Parameters:

* `C_0::Array{Float64}`
* `r_cov_m::Float64`
* `P_star::Float64`
* `gamma_0::Float64`
* `k::Float64`
* `t_0::Int64`
"""
type AMUpdate_gen <: AdaptationAlgorithm
  C_0::Array{Float64}
  r_cov_m::Float64
  P_star::Float64
  gamma_0::Float64
  k::Float64
  t_0::Int64
end

# functions

# set up functions
function set_adaptive_alg_params(algorithm::FixedKernel, nbr_of_unknown_parameters::Integer, theta::Vector,R::Integer)

  return (algorithm.Cov, NaN)

end

function set_adaptive_alg_params(algorithm::AMUpdate, nbr_of_unknown_parameters::Integer, theta::Vector,R::Integer)

  # define m_g m_g_1
  m_g = zeros(nbr_of_unknown_parameters,1)
  m_g_1 = zeros(nbr_of_unknown_parameters,1)

  return (algorithm.C_0, algorithm.gamma_0, algorithm.k, algorithm.t_0, [algorithm.r_cov_m], m_g, m_g_1)

end

function set_adaptive_alg_params(algorithm::AMUpdate_gen, nbr_of_unknown_parameters::Int64, theta::Vector,R::Int64)

  Cov_m = algorithm.C_0
  log_r_cov_m = log(algorithm.r_cov_m)
  log_P_star = log(algorithm.P_star)
  gamma_0 = algorithm.gamma_0
  k = algorithm.k
  t_0 = algorithm.t_0
  vec_log_r_cov_m = zeros(1,R)
  vec_log_r_cov_m[1:t_0] = log_r_cov_m
  # define m_g m_g_1
  m_g = zeros(nbr_of_unknown_parameters,1)
  m_g_1 = zeros(nbr_of_unknown_parameters,1)
  #diff_a_log_P_start = 1.
  #factor = 1.
  #gain = 1.

  return (Cov_m, log_P_star, gamma_0, k, t_0, vec_log_r_cov_m, [log_r_cov_m], m_g, m_g_1)

end

# return covariance functions
function return_covariance_matrix(algorithm::FixedKernel, adaptive_update_params::Tuple,r::Integer)

  return adaptive_update_params[1]

end

function return_covariance_matrix(algorithm::AMUpdate, adaptive_update_params::Tuple,r::Integer)
    r_cov_m = adaptive_update_params[end-2][1]
    Cov_m = adaptive_update_params[1]
    return r_cov_m^2*Cov_m
end

function return_covariance_matrix(algorithm::AMUpdate_gen, adaptive_update_params::Tuple,r::Int64)
    log_r_cov_m = adaptive_update_params[end-2][1]
    Cov_m = adaptive_update_params[1]
    return exp(log_r_cov_m)^2*Cov_m
end

# print_covariance functions
function print_covariance(algorithm::FixedKernel, adaptive_update_params::Tuple,r::Integer)

  println(adaptive_update_params[1])

end

function print_covariance(algorithm::AMUpdate, adaptive_update_params::Tuple,r::Integer)
    r_cov_m = adaptive_update_params[end-2][1]
    Cov_m = adaptive_update_params[1]
    println(r_cov_m^2*Cov_m)
end

function print_covariance(algorithm::AMUpdate_gen, adaptive_update_params::Tuple,r::Int64)
    log_r_cov_m = adaptive_update_params[end-2][1]
    Cov_m = adaptive_update_params[1]
    println(exp(log_r_cov_m)^2*Cov_m)
end

# get_covariance functions
function get_covariance(algorithm::FixedKernel, adaptive_update_params::Tuple,r::Integer)

  return adaptive_update_params[1]

end

function get_covariance(algorithm::AMUpdate, adaptive_update_params::Tuple,r::Integer)
    r_cov_m = adaptive_update_params[end-2][1]
    Cov_m = adaptive_update_params[1]
    return r_cov_m^2*Cov_m
end

function get_covariance(algorithm::AMUpdate_gen, adaptive_update_params::Tuple,r::Int64)
    log_r_cov_m = adaptive_update_params[end-2][1]
    Cov_m = adaptive_update_params[1]
    return exp(log_r_cov_m)^2*Cov_m
end

# gaussian random walk functions
function gaussian_random_walk(algorithm::FixedKernel, adaptive_update_params::Tuple, Theta::Vector, r::Integer)

  return rand(MvNormal(Theta, 1.0*adaptive_update_params[1])), zeros(size(Theta))

end

function gaussian_random_walk(algorithm::AMUpdate, adaptive_update_params::Tuple, Theta::Vector, r::Integer)
  r_cov_m = adaptive_update_params[end-2][1]
  Cov_m = adaptive_update_params[1]
  return rand(MvNormal(Theta, r_cov_m^2*Cov_m)), zeros(size(Theta))
end

function gaussian_random_walk(algorithm::AMUpdate_gen, adaptive_update_params::Tuple, Theta::Vector, r::Int64)
  log_r_cov_m = adaptive_update_params[end-2][1]
  Cov_m = adaptive_update_params[1]
  return rand(MvNormal(Theta, exp(log_r_cov_m)^2*Cov_m)), zeros(size(Theta))
end

# functions for adaptation of parameters
function adaptation(algorithm::FixedKernel, adaptive_update_params::Tuple, Theta::Array, r::Integer,a_log::Real)

  # do nothing

end

function adaptation(algorithm::AMUpdate, adaptive_update_params::Tuple, Theta::Array, r::Integer,a_log::Real)

  Cov_m = adaptive_update_params[1]
  m_g = adaptive_update_params[end-1]
  m_g_1 = adaptive_update_params[end]
  k = adaptive_update_params[3]
  gamma_0 = adaptive_update_params[2]
  t_0 = adaptive_update_params[4]

  g = r-1
  if r-1 == t_0
      m_g = mean(Theta[:,1:r-1],2)
      m_g_1 = m_g + gamma_0/( (g+1)^k ) * (Theta[:,g+1] - m_g)
      adaptive_update_params[1][:] = Cov_m + ( gamma_0/( (g+1)^k ) ) * ( (Theta[:, g+1] - m_g)*(Theta[:,g+1] - m_g)'     -  Cov_m)
      adaptive_update_params[end-1][:] = m_g_1
  elseif r-1 > t_0
      m_g_1 = m_g + gamma_0/( (g+1)^k ) * (Theta[:,g+1] - m_g)
      adaptive_update_params[1][:] = Cov_m + ( gamma_0/( (g+1)^k ) ) * ( (Theta[:, g+1] - m_g)*(Theta[:,g+1] - m_g)'     -  Cov_m)
      adaptive_update_params[end-1][:] = m_g_1
  end

end

function adaptation(algorithm::AMUpdate_gen, adaptive_update_params::Tuple, Theta::Array, r::Int64, a_log::Float64)

  Cov_m = adaptive_update_params[1]
  m_g = adaptive_update_params[end-1]
  m_g_1 = adaptive_update_params[end]
  log_P_star = adaptive_update_params[2]
  log_r_cov_m = adaptive_update_params[end-2][1]
  k = adaptive_update_params[4]
  gamma_0 = adaptive_update_params[3]
  t_0 = adaptive_update_params[5]

  g = r-1;
  if r-1 >= t_0
    diff_a_log_P_start = abs(min(1, exp(a_log)) - exp(log_P_star)) #abs( min(log(1), a_log) - log_P_star)
    factor = 1 #diff_a_log_P_start/abs(log_r_cov_m)
    gain = diff_a_log_P_start/factor
    if min(1, exp(a_log)) < exp(log_P_star) #min(log(1), a_log) < log_P_star
      sign_val = -1
    else
      sign_val = 1
    end
    adaptive_update_params[end-2][1] = log_r_cov_m + sign_val * ( gamma_0/( (g+1)^k ) ) * gain
    adaptive_update_params[6][g+1] = log_r_cov_m
    if r-1 == t_0
      m_g = mean(Theta[:,1:r-1],2)
      m_g_1 = m_g + ( gamma_0/( (g+1)^k ) ) * (Theta[:,g+1] - m_g)
      adaptive_update_params[1][:] = Cov_m + ( gamma_0/( (g+1)^k ) ) * ( (Theta[:,g+1] - m_g)*(Theta[:,g+1] - m_g)'     -  Cov_m)
      adaptive_update_params[end-1][:] = m_g_1
    elseif r-1 > t_0
      m_g_1 = m_g + gamma_0/( (g+1)^k ) * (Theta[:,g+1] - m_g)
      adaptive_update_params[1][:] = Cov_m + ( gamma_0/( (g+1)^k ) ) * ( (Theta[:,g+1] - m_g)*(Theta[:,g+1] - m_g)'     -  Cov_m)
      adaptive_update_params[end-1][:] = m_g_1
    end
  end

end


text:
###############################################################################
#
# Segmented Hydraulic system
#
###############################################################################
#= AbstractHydraulicOrgan type tree
---> LeafHydraulics
---> RootHydraulics
---> StemHydraulics
=#
"""
    abstract type AbstractHydraulicOrgan{FT}

Hierarchy of AbstractHydraulicOrgan
- [`LeafHydraulics`](@ref)
- [`RootHydraulics`](@ref)
- [`StemHydraulics`](@ref)
"""
abstract type AbstractHydraulicOrgan{FT<:AbstractFloat} end

"""
    mutable struct LeafHydraulics{FT}

A struct that contains leaf hydraulics information.

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct LeafHydraulics{FT} <: AbstractHydraulicOrgan{FT}
    # number of xylem slices
    N::Int = 5

    # leaf hydraulic parameters
    "Leaf area `[m²]`"
    area ::FT = FT(1500)
    "Maximal extra-xylary hydraulic conductance `[mol s⁻¹ MPa⁻¹ m⁻²]`"
    k_ox ::FT = FT(100)
    "Maximal leaf hydraulic conductance per leaf area `[mol s⁻¹ MPa⁻¹ m⁻²]`"
    k_sla::FT = FT(0.04)
    "Vulnerability curve"
    vc::AbstractXylemVC{FT} = WeibullSingle{FT}()
    "Critical xylem pressure `[MPa]`"
    p_crt::FT = -vc.b * log(FT(1000)) ^ (1/vc.c)

    # flows and pressures (need to be updated with time)
    "Flow rate in the xylem `[mol s⁻¹]`"
    flow  ::FT = FT(0)
    "Leaf xylem water pressure at the downstream end of leaf xylem `[MPa]`"
    p_dos ::FT = FT(0.0)
    "Leaf end water pressure at air-water interface `[MPa]`"
    p_leaf::FT = FT(0.0)
    "Leaf xylem water pressure at the leaf base (upstream) `[MPa]`"
    p_ups ::FT = FT(0.0)

    # pressure, k, and p_history profile
    "List of leaf k_max per element `[mol s⁻¹ MPa⁻¹ m⁻²]`"
    k_element::Array{FT,1} =  ones(FT,N) .* N .* k_sla
    "List of leaf kr history per element"
    k_history::Array{FT,1} =  ones(FT,N)
    "List of xylem water pressure `[MPa]`"
    p_element::Array{FT,1} = zeros(FT,N)
    "List of xylem water pressure history (normalized to 298.15 K) `[MPa]`"
    p_history::Array{FT,1} = zeros(FT,N)

    # temperature related (uniform leaf temperature), update with time
    "Relative surface tension"
    f_st ::FT = FT(1)
    "Relative viscosity"
    f_vis::FT = FT(1)
    "Temperature memory `[K]`"
    T_old::FT = K_25(FT)
    "Upstream sap temperature `[K]`"
    T_sap::FT = K_25(FT)

    # capacitance
    "Pressure volume curve for storage"
    pv::AbstractCapacity{FT} = PVCurveSegmented{FT}()
    "Pressure of storage"
    p_storage::FT = 0
    "Total capaciatance at Ψ = 0 `[mol m⁻²]`"
    v_maximum::FT = 20
    "Current capaciatance at Ψ_leaf `[mol m⁻²]`"
    v_storage::FT = 20
    "Flow rate into the tissue (used for non-steady state) `[mol m⁻² s⁻¹]`"
    q_in     ::FT = 0
    "Flow rate out of the tissue (used for non-steady state) `[mol m⁻² s⁻¹]`"
    q_out    ::FT = 0
end

"""
    mutable struct RootHydraulics{FT}

A struct that contains root hydraulics information.

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct RootHydraulics{FT} <: AbstractHydraulicOrgan{FT}
    # number of xylem slices
    N::Int = 5

    # root hydraulic parameters
    "Root cross-section area `[m²]`"
    area ::FT = FT(1)
    "Maximal hydraulic conductance `[mol s⁻¹ MPa⁻¹]`"
    k_max::FT = FT(25)
    "Maximal xylem hydraulic conductivity `[mol s⁻¹ MPa⁻¹ m⁻²]`"
    k_s  ::FT = FT(250)
    "Vulnerability curve"
    vc::AbstractXylemVC{FT} = WeibullSingle{FT}()
    "Root z difference `[m]`"
    Δh   ::FT = FT(1.0)

    # soil parameters
    "Rhizosphere  conductance `[mol s⁻¹ MPa⁻¹]`"
    k_rhiz   ::FT     = FT(5e14)
    "Soil hydraulics"
    sh::AbstractSoilVC{FT} = BrooksCorey{FT}()

    # flows and pressures (need to be updated with time)
    "Flow rate in the xylem `[mol s⁻¹]`"
    flow  ::FT = FT(0.0)
    "Xylem water pressure at the downstream end of xylem `[MPa]`"
    p_dos ::FT = FT(0.0)
    "Xylem-rhizosphere interface water pressure `[MPa]`"
    p_rhiz::FT = FT(0.0)
    "Soil matrix potential `[MPa]`"
    p_ups ::FT = FT(0.0)
    "Soil osmotic potential at 298.15 K `[MPa]"
    p_osm ::FT = FT(0.0)

    # pressure, k, and p_history profile
    "List of k_max per element `[mol s⁻¹ MPa⁻¹ m⁻²]`"
    k_element::Array{FT,1} =  ones(FT,N) .* N .* k_max
    "List of kr history per element"
    k_history::Array{FT,1} =  ones(FT,N)
    "List of xylem water pressure `[MPa]`"
    p_element::Array{FT,1} = zeros(FT,N)
    "List of pressure drop caused by gravity `[MPa]`"
    p_gravity::Array{FT,1} = zeros(FT,N) .+ ρg_MPa(FT) .* Δh ./ N;
    "List of xylem water pressure history (normalized to 298.15 K) `[MPa]`"
    p_history::Array{FT,1} = zeros(FT,N)

    # temperature related (uniform temperature), update with time
    "Relative surface tension"
    f_st ::FT = FT(1)
    "Relative viscosity"
    f_vis::FT = FT(1)
    "Temperature memory `[K]`"
    T_old::FT = K_25(FT)
    "Upstream sap temperature `[K]`"
    T_sap::FT = K_25(FT)

    # capacitance
    "Pressure volume curve for storage"
    pv::AbstractCapacity{FT} = PVCurveLinear{FT}()
    "Pressure of storage per element"
    p_storage::Array{FT,1} = zeros(FT,N)
    "Maximal storage per element `[mol]`"
    v_maximum::Array{FT,1} = area * Δh / N * 6000 * ones(FT,N)
    "Storage per element `[mol]`"
    v_storage::Array{FT,1} = area * Δh / N * 6000 * ones(FT,N)
    "List of xylem water flow `[mol m⁻²]`"
    q_element::Array{FT,1} = zeros(FT,N)
    "List of buffer water flow `[mol m⁻²]`"
    q_buffer ::Array{FT,1} = zeros(FT,N)
    "List of diiferntial water flow `[mol m⁻²]`"
    q_diff   ::Array{FT,1} = zeros(FT,N)
    "Flow rate into the tissue (used for non-steady state) `[mol s⁻¹]`"
    q_in ::FT = 0
    "Flow rate out of the tissue (used for non-steady state) `[mol s⁻¹]`"
    q_out::FT = 0
end

"""
    mutable struct StemHydraulics{FT}

A struct that contains stem hydraulics information.

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct StemHydraulics{FT} <: AbstractHydraulicOrgan{FT}
    # number of xylem slices
    N::Int = 5

    # stem hydraulic parameters
    "Stem cross-section area `[m²]`"
    area ::FT = FT(1)
    "Maximal hydraulic conductance `[mol s⁻¹ MPa⁻¹]`"
    k_max::FT = FT(50)
    "Maximal xylem hydraulic conductivity `[mol s⁻¹ MPa⁻¹ m⁻²]`"
    k_s  ::FT = FT(250)
    "Vulnerability curve"
    vc::AbstractXylemVC{FT} = WeibullSingle{FT}()
    "Stem height difference `[m]`"
    Δh   ::FT = FT(5.0)

    # flows and pressures (need to be updated with time)
    "Flow rate in the xylem `[mol s⁻¹]`"
    flow  ::FT = FT(0.0)
    "Xylem water pressure at the downstream end of xylem `[MPa]`"
    p_dos::FT = FT(0.0)
    "Xylem water pressure at the base (upstream) `[MPa]`"
    p_ups::FT = FT(0.0)

    # pressure, k, and p_history profile
    "List of k_max per element `[mol s⁻¹ MPa⁻¹ m⁻²]`"
    k_element::Array{FT,1} =  ones(FT,N) .* N .* k_max
    "List of kr history per element"
    k_history::Array{FT,1} =  ones(FT,N)
    "List of xylem water pressure `[MPa]`"
    p_element::Array{FT,1} = zeros(FT,N)
    "List of pressure drop caused by gravity `[MPa]`"
    p_gravity::Array{FT,1} = zeros(FT,N) .+ ρg_MPa(FT) .* Δh ./ N;
    "List of xylem water pressure history (normalized to 298.15 K) `[MPa]`"
    p_history::Array{FT,1} = zeros(FT,N)

    # temperature related (uniform temperature), update with time
    "Relative surface tension"
    f_st ::FT = FT(1)
    "Relative viscosity"
    f_vis::FT = FT(1)
    "Temperature memory `[K]`"
    T_old::FT = K_25(FT)
    "Upstream sap temperature `[K]`"
    T_sap::FT = K_25(FT)

    # capacitance
    "Pressure volume curve for storage"
    pv::AbstractCapacity{FT} = PVCurveLinear{FT}()
    "Pressure of storage per element"
    p_storage::Array{FT,1} = zeros(FT,N)
    "Maximal storage per element `[mol]`"
    v_maximum::Array{FT,1} = area * Δh / N * 6000 * ones(FT,N)
    "Storage per element `[mol]`"
    v_storage::Array{FT,1} = area * Δh / N * 6000 * ones(FT,N)
    "List of xylem water flow `[mol m⁻²]`"
    q_element::Array{FT,1} = zeros(FT,N)
    "Flow rate into the tissue (used for non-steady state) `[mol s⁻¹]`"
    q_in ::FT = 0
    "Flow rate out of the tissue (used for non-steady state) `[mol s⁻¹]`"
    q_out::FT = 0
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
### A Pluto.jl notebook ###
# v0.10.2

using Markdown
using InteractiveUtils

# ╔═╡ a0f0ae42-9ad8-11ea-2475-a735df64aa28
using Plots

# ╔═╡ 7b93882c-9ad8-11ea-0288-0941e163f9d5
md"""
# Plotting in Pluto

Pluto is an excellent enviroment to visualise your data. This notebook shows a few examples of using plots in Pluto using the `Plots` package.

I wrote this notebook so you can understand it even if you have never used `Plots`. However, it is not intend as a complete tutorial for the package. If you want to start making your own plots, I recommend looking at the package [documentation](https://docs.juliaplots.org/latest/) for a full tutorial as well.

Let's start by importing the `Plots` package. You can just put a cell like this anywhere in your notebook.
"""

# ╔═╡ 5ae65950-9ad9-11ea-2e14-35119d369acd
md"""
## The basics
Now we can get started. We can begin by making some X and Y values. I typed out some Y values so you can alter them and see the effect.
"""

# ╔═╡ aaa805d4-9ad8-11ea-21c2-3b20580fea0e
years = 2001:2010

# ╔═╡ c02bc44a-9ad8-11ea-117d-3997e1f3cab0
apples = [15, 25, 80, 75, 50, 30, 35, 15, 25, 35]

# ╔═╡ a405ae4c-9ad9-11ea-0008-f763e098846d
md"""
Great, let's plot them! The basic syntax of of `Plots` is very simple.
"""

# ╔═╡ 12a8c222-9ad9-11ea-2544-355bd080367f
plot(years, apples)

# ╔═╡ 19660da0-9ada-11ea-1942-a9d4ddd7753b
md"""
You can change the properties of the plot by giving different arguments to the function call.
"""

# ╔═╡ 5410faba-9ada-11ea-2fc5-35f0bdd40d96
plot(years, apples, legend = false, title = "Number of apples per year")

# ╔═╡ 90c47cae-9ada-11ea-10da-91f4dcf0f2d6
md"""
Alternatively, you can use `plot!()`. This function takes the same arguments, but alters our existing plot.
"""

# ╔═╡ b2727d2e-9ada-11ea-23d3-bddab9ec8de1
begin
	plot(years, apples)
	plot!(legend = false)
	plot!(title =  "Number of apples per year")
end

# ╔═╡ eb218bea-9adc-11ea-2f3a-298c41caf8a1
md"""
## More options

We can easily plot multiple series.
"""

# ╔═╡ bd30bdf4-9add-11ea-26d4-75ac3f0b8610
oranges = rand(1:100, 10)

# ╔═╡ dc50d480-9add-11ea-2821-018a98fab262
begin
	plot(years, apples, label = "apples")
	plot!(years, oranges, label = "oranges")
end

# ╔═╡ c9a64e98-9b67-11ea-1d6a-3b1f7548cb8d
md"""
Different plot types have their own functions.
"""

# ╔═╡ 13428812-9ade-11ea-3865-1d293993203d
scatter(apples, oranges, legend = false)

# ╔═╡ d76fd0f6-9b67-11ea-0a73-e9ddb60ca291
md"""
You can also plot a function instead of an output array!
"""

# ╔═╡ 69fa7744-9b67-11ea-26df-3b4446dfb3ea
begin
	function myfunction(x)
		return sin(x/2)
	end

	plot(years, sin, label = "sine")
	plot!(years, myfunction, label = "my function")
end

# ╔═╡ af31a428-9adf-11ea-280a-45e2b891c4f9
md"""
## Naming plots

Plots are also just a type of variable, so you can give them a name.
"""

# ╔═╡ 5ed7c444-9adf-11ea-2925-8573f018e820
apples_plot = plot(years, apples,
	legend = false, xlabel = "year", ylabel = "apples")

# ╔═╡ 932b6b3c-a64d-11ea-2aba-2bf577ce3cf8
apples_vs_oranges = scatter(apples, oranges,
	legend = false, xlabel = "apples", ylabel = "oranges")

# ╔═╡ 030c1d10-9b69-11ea-11b7-ff04c60892b9
md"""
This way, you can show them again
"""

# ╔═╡ d7bdab72-9b65-11ea-18c0-5353a214749d
apples_plot

# ╔═╡ 22992cec-9adb-11ea-078e-c7c21489c053
md"""
## Combining subplots

We can combine multiple plots in a single image using the `layout` argument. One option is to provide a matrix of values and call a single `plot`.
"""

# ╔═╡ 5d723220-9ade-11ea-37db-c110dda60a15
pears = rand(1:100, 10)

# ╔═╡ 632cb26c-9ade-11ea-1bc7-39b8207190e2
bananas = rand(1:100, 10)

# ╔═╡ a01bf430-9ade-11ea-2ae4-23287ba11406
fruits = hcat(apples, oranges, pears, bananas)

# ╔═╡ 69e045b0-9ade-11ea-0c55-05fa0da5893c
plot(years, fruits, layout = (4,1), legend = false)

# ╔═╡ 054e0c4e-9adf-11ea-3deb-4daf2b6a2548
md"""
This works when we have a nice output matrix. For more complicated options, named plots are usually a good solution. When we have created a few plots, we can put them together in a single plot.
"""

# ╔═╡ 91b4eedc-9adf-11ea-37aa-250ceec26640
plot(apples_plot, apples_vs_oranges, layout = (1, 2))

# ╔═╡ a7abf64e-9adb-11ea-2ae4-47886125fe60
md"""
## The Pluto way of plotting

For the most part, plotting in Pluto is not different from anywhere else. However, there are a few things to keep in mind.

The `plot!()` function alters an existing plot. In reactive programming, you are not supposed to alter the value of a variable you defined in a different cell. I strongly recommend that you only use `plot!()` to alter plots you initialised in the same cell.

For the sake of demonstration, here is what happens if you use `plot!()` in its own cell. Let's start by making a new plot.
"""

# ╔═╡ 2617a7b6-9ae2-11ea-1461-d914c7cfdb09
pears_plot = plot(years, pears, label = "pears")

# ╔═╡ 46388876-9ae2-11ea-2c5f-073698e45d44
md"""
The cell below adds a new series to the plot. It shows the altered plot, while our original plot still looks intact.

Maybe that's not so bad. It seems like we could use this to build a plot one step at a time. But try changing the plot below so that it plots oranges instead of bananas.
"""

# ╔═╡ 2598a29a-9ae2-11ea-02cc-21f3fa5a9472
plot!(pears_plot, years, bananas, label = "bananas")

# ╔═╡ a9cdba00-9ae2-11ea-3bf9-41d07e7a785f
md"""
Did you get three series? It's because the `plot` still has the `bananas` data added to it. So the bananas are still plotted, even though you now don't have any code doing that. This is not great for coding: it means that if you close your notebook, your plots might not look the same when you open it again.

Try running the cell a few times. You keep adding series! Also not great when working on a plot.
"""

# ╔═╡ f22b6fda-9ae3-11ea-3ff2-ebcf1374f595
md"""
So what if you do want to add things to a plot one at a time, instead of using a single cell? My recommendation is to define your plot in a function.
"""

# ╔═╡ 8819edf0-9ae4-11ea-15f4-f17e9e9db8ea
function plot_pears()
	myplot = plot(years, pears, label = "pears")
	return myplot
end

# ╔═╡ 9e4d1462-9ae4-11ea-3c0f-774d68a671c3
plot_pears()

# ╔═╡ a7906786-9ae4-11ea-0175-4f860c05a8a2
begin
	plot_pears()
	plot!(years, bananas, label = "bananas")
end

# ╔═╡ d3cbc548-9ae4-11ea-261f-7fb956839e53
md"""
Try changing the function for the base plot, or the addition in the bottom cell. No weird effects!
"""

# ╔═╡ Cell order:
# ╟─7b93882c-9ad8-11ea-0288-0941e163f9d5
# ╠═a0f0ae42-9ad8-11ea-2475-a735df64aa28
# ╟─5ae65950-9ad9-11ea-2e14-35119d369acd
# ╠═aaa805d4-9ad8-11ea-21c2-3b20580fea0e
# ╠═c02bc44a-9ad8-11ea-117d-3997e1f3cab0
# ╟─a405ae4c-9ad9-11ea-0008-f763e098846d
# ╠═12a8c222-9ad9-11ea-2544-355bd080367f
# ╟─19660da0-9ada-11ea-1942-a9d4ddd7753b
# ╠═5410faba-9ada-11ea-2fc5-35f0bdd40d96
# ╟─90c47cae-9ada-11ea-10da-91f4dcf0f2d6
# ╠═b2727d2e-9ada-11ea-23d3-bddab9ec8de1
# ╟─eb218bea-9adc-11ea-2f3a-298c41caf8a1
# ╠═bd30bdf4-9add-11ea-26d4-75ac3f0b8610
# ╠═dc50d480-9add-11ea-2821-018a98fab262
# ╟─c9a64e98-9b67-11ea-1d6a-3b1f7548cb8d
# ╠═13428812-9ade-11ea-3865-1d293993203d
# ╟─d76fd0f6-9b67-11ea-0a73-e9ddb60ca291
# ╠═69fa7744-9b67-11ea-26df-3b4446dfb3ea
# ╟─af31a428-9adf-11ea-280a-45e2b891c4f9
# ╠═5ed7c444-9adf-11ea-2925-8573f018e820
# ╠═932b6b3c-a64d-11ea-2aba-2bf577ce3cf8
# ╟─030c1d10-9b69-11ea-11b7-ff04c60892b9
# ╠═d7bdab72-9b65-11ea-18c0-5353a214749d
# ╟─22992cec-9adb-11ea-078e-c7c21489c053
# ╠═5d723220-9ade-11ea-37db-c110dda60a15
# ╠═632cb26c-9ade-11ea-1bc7-39b8207190e2
# ╠═a01bf430-9ade-11ea-2ae4-23287ba11406
# ╠═69e045b0-9ade-11ea-0c55-05fa0da5893c
# ╟─054e0c4e-9adf-11ea-3deb-4daf2b6a2548
# ╠═91b4eedc-9adf-11ea-37aa-250ceec26640
# ╟─a7abf64e-9adb-11ea-2ae4-47886125fe60
# ╠═2617a7b6-9ae2-11ea-1461-d914c7cfdb09
# ╟─46388876-9ae2-11ea-2c5f-073698e45d44
# ╠═2598a29a-9ae2-11ea-02cc-21f3fa5a9472
# ╟─a9cdba00-9ae2-11ea-3bf9-41d07e7a785f
# ╟─f22b6fda-9ae3-11ea-3ff2-ebcf1374f595
# ╠═8819edf0-9ae4-11ea-15f4-f17e9e9db8ea
# ╠═9e4d1462-9ae4-11ea-3c0f-774d68a671c3
# ╠═a7906786-9ae4-11ea-0175-4f860c05a8a2
# ╟─d3cbc548-9ae4-11ea-261f-7fb956839e53


text:
using EconPDEs

Base.@kwdef  struct CampbellCochraneModel
    # consumption process parameters
    μ::Float64 = 0.0189
    σ::Float64 = 0.015
    # utility
    γ::Float64 = 2.0
    ρ::Float64 = 0.116
    # habit
    κs::Float64 = 0.138
    b::Float64 = 0.0
end
# I choose persistence so that monthly simulation of the model matches processes in CC (1999)
# ρ = 12 * (1 - 0.89^(1/12))
# κs = 12 * (1 - 0.87^(1/12))

function initialize_stategrid(m::CampbellCochraneModel; sn = 1000)
    μ = m.μ ; σ = m.σ ; γ = m.γ ; ρ = m.ρ ; κs = m.κs ; b = m.b
    Sbar = σ * sqrt(γ / (κs - b / γ))
    sbar = log.(Sbar)
    smax =  sbar + 0.5 * (1 - Sbar^2)
    # corresponds to Grid 3 in Wachter (2005)
    shigh = log.(range(0.0, exp(smax), length = div(sn, 10)))
    slow = range(-300.0, shigh[2], length = sn - div(sn, 10))
    OrderedDict(:s => vcat(slow[1:(end-1)], shigh[2:end]))
end

function (m::CampbellCochraneModel)(state::NamedTuple, y::NamedTuple)
    (; μ, σ, γ, ρ, κs, b) = m
    (; s) = state
    (; p, ps_up, ps_down, pss) = y

    # drift and volatility of  s and p
    Sbar = σ * sqrt(γ / (κs - b / γ))
    sbar = log(Sbar)
    λ = 1 / Sbar * sqrt(1 - 2 * (s - sbar)) - 1
    μs = - κs * (s - sbar)
    ps = (μs >= 0) ? ps_up : ps_down
    σs = λ * σ
    σp = ps / p * σs
    μp = ps / p * μs + 0.5 * pss / p * σs^2

    # market price of risk κ
    κ = γ * (σ + σs)

    # risk free rate  r
    r = ρ + γ * μ - (γ * κs - b) / 2 + b * (sbar - s)

    # PDE
    pt = - p * (1 / p + μ + μp + σp * σ - r - κ * (σ + σp))
    return (; pt)
end

# Campbell Cochrane (1999)
m = CampbellCochraneModel()
stategrid = initialize_stategrid(m)
yend = OrderedDict(:p => ones(length(stategrid[:s])))
result = pdesolve(m, stategrid, yend)
@assert result.residual_norm <= 1e-5

# Wachter (2005) calibration
# m = CampbellCochraneModel(μ = 0.022, σ = 0.0086, γ = 2.0, ρ = 0.073, κs = 0.116, b = 0.011)
# stategrid = initialize_stategrid(m)
# yend = OrderedDict(:p => ones(length(stategrid[:s])))
# y, result, distance = pdesolve(m, stategrid, yend)


text:
#=
pfile.jl
read GE MRI k-space data from .p file
based on https://gitlab.com/fMRI/toppe/+toppe/+utils/loadpfile.m
=#

export loadpfile

# using MIRT: read_rdb_hdr

"""
`(dat, rdb_hdr) = loadpfile(pfile ; ...)`

Load data for one or more echoes from GE MRI scan Pfile.

in
- `pfile::String` filename

option
- `coils::AbstractVector{<:Integer}`
	only get data for these coils; default: all coils
- `echoes::AbstractVector{<:Integer}`
	only get data for these echoes; default: all echoes
- `slices::AbstractVector{<:Integer}`
	only get data for these slices; default: `2:nslices` (NB!)
	because first slice (dabslice=0 slot) may contain corrupt data.
- `views::AbstractVector{<:Integer}`
	only get data for these views; default: all views
- `quiet::Bool`	non-verbosity, default `false`

out
- `dat::Array{Complex{Int16}}` `[ndat, ncoil, nslice, necho, nview]`
- `rdb_hdr::NamedTuple`	header information

Note that to save memory the output type is complex-valued Int16.

Comments from Matlab version at
https://gitlab.com/fMRI/toppe/+toppe/+utils/loadpfile.m

This file was derived from part of the TOPPE development environment
for platform-independent MR pulse programming.

TOPPE is free software: you can redistribute it and/or modify
it under the terms of the GNU Library General Public License as published by
the Free Software Foundation version 2.0 of the License.

TOPPE is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See the GNU Library General Public License for more details.

You should have received a copy of the GNU Library General Public License
along with TOPPE.
If not, see <http://www.gnu.org/licenses/old-licenses/lgpl-2.0.html>.

(c) 2016 The Regents of the University of Michigan
Jon-Fredrik Nielsen, jfnielse@umich.edu

2019-05-22 Julia version by Jeff Fessler
"""
function loadpfile(pfile::String ;
		coils::AbstractVector{<:Integer} = empty([], Integer),
		echoes::AbstractVector{<:Integer} = empty([], Integer),
		slices::AbstractVector{<:Integer} = empty([], Integer),
		views::AbstractVector{<:Integer} = empty([], Integer),
		quiet::Bool = false,
	)

	fid = open(pfile, "r") # open pfile
	rdb_hdr = read_rdb_hdr(fid) # read header

	# Header parameters
	ndat = rdb_hdr.frame_size
	nslices = rdb_hdr.nslices
	ptsize = rdb_hdr.point_size # 2: data stored in short int format, int16)
	if ptsize != 2 # 4 (extended precision) unsupported here
		throw("ptsize = $ptsize unsupported")
	end
	nechoes = rdb_hdr.nechoes
	nviews = rdb_hdr.nframes
	ncoils = rdb_hdr.dab[2] - rdb_hdr.dab[1] + 1

	# Calculate size of data chunks
	# See pfilestruct.jpg, and rhrawsize calculation in .e file.
	# number of data points per 'echo' loaddab slot. Includes baseline (0) view:
	echores  = ndat * (nviews+1)
	sliceres = nechoes * echores # number of data points per 'slice'
	coilres  = nslices * sliceres # number of data points per receive coil

	# This size should match the Pfile size exactly:
	pfilesize = rdb_hdr.off_data +
		2 * ptsize * ncoils * nslices * nechoes * (nviews+1) * ndat

	# File size check
	if pfilesize != filesize(pfile)
		throw("Expected $(pfilesize/1e6) MB but see $(filesize(pfile)/1e6) MB file")
	#	fprintf('Press enter to continue anyway...');
	#	input('');
	end

	# Check if second view is empty
	# This happens when there's only one view, but the scanner sets nviews = 2
	if nviews == 2
		# Seek to view slice 2, echo 1, view 2, coil 1:
		seek(fid, rdb_hdr.off_data + (2 * ptsize * (sliceres + 2*ndat)))
		view2tmp = Array{Int16}(undef, 2*ndat)
		view2tmp = read!(fid, view2tmp)
		if all(view2tmp .== Int16(0))
			nviews = 1 # set nviews to 1 so we don't read in all the empty data
			!quiet && @warn("View 2 appears to be empty, only loading view 1.")
		end
	end

	# Determine which portions to load
	# NB: default is to skip first slice (sometimes contains corrupted data)
	isempty(coils) && (coils = 1:ncoils)
	isempty(echoes) && (echoes = 1:nechoes)
	isempty(slices) && (slices = 2:nslices)
	isempty(views) && (views = 1:nviews)

	# output dimensions:
	dims = (ndat, length(coils), length(slices), length(echoes), length(views))

	# Warn if output size will exceed system memory.
	# Exceeding free RAM may cause a system freeze due to attempted mem swap
	memneeded = 4 * prod(dims) # 4 bytes per complex value because Int16
	memfree = Sys.total_memory() # system total memory
	mempercentuse = 100 * memneeded / memfree
	if mempercentuse > 100
		@warn("Loading data ($(memneeded/1e9) GB) may exceed available RAM")
		@warn("This could freeze computer!")
	#	fprintf('Press enter to continue anyway...');
	#	input('');
	elseif mempercentuse > 90 # Warn if we will we use 90% of memory
		@warn("Loading data ($(memneeded/1e9) GB) will use $mempercentuse % of your available RAM. Proceed with caution!")
	end

	if !quiet
		@info("ndat = $ndat, memory = $(memneeded/1e9) GB")
		@info("slices: $(slices[1])-$(slices[end]) / $nslices")
		@info("echoes: $(echoes[1])-$(echoes[end]) / $nechoes")
		@info("views: $(views[1])-$(views[end]) / $nviews")
		@info("coils: $(coils[1])-$(coils[end]) / $ncoils")
	end

	# Read data from file
	# Julia stores complex data as real1,imag1, real2,imag2, ... unlike matlab!
	# This simplifies the IO here!
	data = Array{Complex{Int16}}(undef, dims)
	dtmp = Array{Complex{Int16}}(undef, ndat) # one readout
#	!quiet && textprogressbar('Loading data: ')
	for icoil = coils
		@show icoil
	#	!quiet && textprogressbar(icoil/ncoils*100)
		for islice = slices
			for iecho = echoes
				for iview = views
					offsetres = (icoil-1)*coilres + (islice-1)*sliceres + (iecho-1)*echores + iview*ndat
					offsetbytes = 2 * ptsize * offsetres
					seek(fid, rdb_hdr.off_data + offsetbytes)
					read!(fid, dtmp)
					data[:, icoil-coils[1]+1, islice-slices[1]+1,
						iecho-echoes[1]+1, iview-views[1]+1] .= dtmp
				end
			end
		end
	end

	close(fid)
#	!quiet && textprogressbar(' done.')

	return (data, rdb_hdr)
end

"""
`loadpfile(file, echo::Integer ; ...)`

load a single echo
"""
function loadpfile(pfile::String, echo::Integer ; kwarg...)
	return loadpfile(pfile, echoes=[echo], kwarg...)
end

"""
`loadpfile(:test)`

does nothing now because tests require huge data files
"""
function loadpfile(test::Symbol)
	true
end


text:
@testset "Model run" begin
  param_file = "params.yml"
  @test_nowarn model = Science_ABM.initialize_model(param_file)
  model = Science_ABM.initialize_model(param_file)
  if model.properties["replicates"] <= 1
    @test_nowarn adata, mdata = run!(model, dummystep, Science_ABM.model_step!, model.properties["years"], adata=Science_ABM.adata, when =  model.properties["save_gens"])
  else
    @test_nowarn adata, mdata, models = ensemblerun!([Science_ABM.initialize_model(param_file) for i in model.properties["replicates"]], dummystep, Science_ABM.model_step!, model.properties["years"], adata=Science_ABM.adata, when =  model.properties["save_gens"])
  end
end


text:
using Oceananigans.Advection: boundary_buffer
using Oceananigans.Operators: ℑxᶠᵃᵃ, ℑxᶜᵃᵃ, ℑyᵃᶠᵃ, ℑyᵃᶜᵃ, ℑzᵃᵃᶠ, ℑzᵃᵃᶜ
using Oceananigans.TurbulenceClosures: AbstractTurbulenceClosure, AbstractTimeDiscretization

const ATC = AbstractTurbulenceClosure
const ATD = AbstractTimeDiscretization

struct RasterDepthMask end

struct GridFittedBoundary{M, S} <: AbstractImmersedBoundary
    mask :: S
    mask_type :: M
end

GridFittedBoundary(mask; mask_type=nothing) = GridFittedBoundary(mask, mask_type)

@inline is_immersed(i, j, k, underlying_grid, ib::GridFittedBoundary) = ib.mask(node(c, c, c, i, j, k, underlying_grid)...)
@inline is_immersed(i, j, k, underlying_grid, ib::GridFittedBoundary{<:RasterDepthMask}) = ib.mask(i, j)

const IBG = ImmersedBoundaryGrid

const c = Center()
const f = Face()

#####
##### ImmersedBoundaryGrid
#####

@inline solid_cell(i, j, k, ibg) = is_immersed(i, j, k, ibg.grid, ibg.immersed_boundary)

@inline solid_node(LX, LY, LZ, i, j, k, ibg) = solid_cell(i, j, k, ibg) # fallback (for Center or Nothing LX, LY, LZ)

@inline solid_node(::Face, LX, LY, i, j, k, ibg) = solid_cell(i, j, k, ibg) || solid_cell(i-1, j, k, ibg)
@inline solid_node(LX, ::Face, LZ, i, j, k, ibg) = solid_cell(i, j, k, ibg) || solid_cell(i, j-1, k, ibg)
@inline solid_node(LX, LY, ::Face, i, j, k, ibg) = solid_cell(i, j, k, ibg) || solid_cell(i, j, k-1, ibg)

@inline solid_node(::Face, ::Face, LZ, i, j, k, ibg) = solid_node(c, f, c, i, j, k, ibg) || solid_node(c, f, c, i-1, j, k, ibg)
@inline solid_node(::Face, LY, ::Face, i, j, k, ibg) = solid_node(c, c, f, i, j, k, ibg) || solid_node(c, c, f, i-1, j, k, ibg)
@inline solid_node(LX, ::Face, ::Face, i, j, k, ibg) = solid_node(c, f, c, i, j, k, ibg) || solid_node(c, f, c, i, j, k-1, ibg)

@inline solid_node(::Face, ::Face, ::Face, i, j, k, ibg) = solid_node(c, f, f, i, j, k, ibg) || solid_node(c, f, f, i-1, j, k, ibg)

@inline conditional_flux_ccc(i, j, k, ibg::IBG{FT}, grid, flux, args...) where FT = ifelse(solid_node(c, c, c, i, j, k, ibg), zero(FT), flux(i, j, k, grid, args...))
@inline conditional_flux_ffc(i, j, k, ibg::IBG{FT}, grid, flux, args...) where FT = ifelse(solid_node(f, f, c, i, j, k, ibg), zero(FT), flux(i, j, k, grid, args...))
@inline conditional_flux_fcf(i, j, k, ibg::IBG{FT}, grid, flux, args...) where FT = ifelse(solid_node(f, c, f, i, j, k, ibg), zero(FT), flux(i, j, k, grid, args...))
@inline conditional_flux_cff(i, j, k, ibg::IBG{FT}, grid, flux, args...) where FT = ifelse(solid_node(c, f, f, i, j, k, ibg), zero(FT), flux(i, j, k, grid, args...))

@inline conditional_flux_fcc(i, j, k, ibg::IBG{FT}, grid, flux, args...) where FT = ifelse(solid_node(f, c, c, i, j, k, ibg), zero(FT), flux(i, j, k, grid, args...))
@inline conditional_flux_cfc(i, j, k, ibg::IBG{FT}, grid, flux, args...) where FT = ifelse(solid_node(c, f, c, i, j, k, ibg), zero(FT), flux(i, j, k, grid, args...))
@inline conditional_flux_ccf(i, j, k, ibg::IBG{FT}, grid, flux, args...) where FT = ifelse(solid_node(c, c, f, i, j, k, ibg), zero(FT), flux(i, j, k, grid, args...))

#####
##### Advective fluxes
#####

# ccc, ffc, fcf
 @inline viscous_flux_ux(i, j, k, ibg::IBG, disc::ATD, clo::ATC, args...) = conditional_flux_ccc(i, j, k, ibg, ibg.grid, viscous_flux_ux, disc, clo, args...)
 @inline viscous_flux_uy(i, j, k, ibg::IBG, disc::ATD, clo::ATC, args...) = conditional_flux_ffc(i, j, k, ibg, ibg.grid, viscous_flux_uy, disc, clo, args...)
 @inline viscous_flux_uz(i, j, k, ibg::IBG, disc::ATD, clo::ATC, args...) = conditional_flux_fcf(i, j, k, ibg, ibg.grid, viscous_flux_uz, disc, clo, args...)

 # ffc, ccc, cff
 @inline viscous_flux_vx(i, j, k, ibg::IBG, disc::ATD, clo::ATC, args...) = conditional_flux_ffc(i, j, k, ibg, ibg.grid, viscous_flux_vx, disc, clo, args...)
 @inline viscous_flux_vy(i, j, k, ibg::IBG, disc::ATD, clo::ATC, args...) = conditional_flux_ccc(i, j, k, ibg, ibg.grid, viscous_flux_vy, disc, clo, args...)
 @inline viscous_flux_vz(i, j, k, ibg::IBG, disc::ATD, clo::ATC, args...) = conditional_flux_cff(i, j, k, ibg, ibg.grid, viscous_flux_vz, disc, clo, args...)

 # fcf, cff, ccc
 @inline viscous_flux_wx(i, j, k, ibg::IBG, disc::ATD, clo::ATC, args...) = conditional_flux_fcf(i, j, k, ibg, ibg.grid, viscous_flux_wx, disc, clo, args...)
 @inline viscous_flux_wy(i, j, k, ibg::IBG, disc::ATD, clo::ATC, args...) = conditional_flux_cff(i, j, k, ibg, ibg.grid, viscous_flux_wy, disc, clo, args...)
 @inline viscous_flux_wz(i, j, k, ibg::IBG, disc::ATD, clo::ATC, args...) = conditional_flux_ccc(i, j, k, ibg, ibg.grid, viscous_flux_wz, disc, clo, args...)

# fcc, cfc, ccf
@inline diffusive_flux_x(i, j, k, ibg::IBG, disc::ATD, clo::ATC, args...) = conditional_flux_fcc(i, j, k, ibg, ibg.grid, diffusive_flux_x, disc, clo, args...)
@inline diffusive_flux_y(i, j, k, ibg::IBG, disc::ATD, clo::ATC, args...) = conditional_flux_cfc(i, j, k, ibg, ibg.grid, diffusive_flux_y, disc, clo, args...)
@inline diffusive_flux_z(i, j, k, ibg::IBG, disc::ATD, clo::ATC, args...) = conditional_flux_ccf(i, j, k, ibg, ibg.grid, diffusive_flux_z, disc, clo, args...)

#####
##### Advective fluxes
#####

@inline _advective_momentum_flux_Uu(i, j, k, ibg::IBG, args...) = conditional_flux_ccc(i, j, k, ibg, ibg, advective_momentum_flux_Uu, args...)
@inline _advective_momentum_flux_Vu(i, j, k, ibg::IBG, args...) = conditional_flux_ffc(i, j, k, ibg, ibg, advective_momentum_flux_Vu, args...)
@inline _advective_momentum_flux_Wu(i, j, k, ibg::IBG, args...) = conditional_flux_fcf(i, j, k, ibg, ibg, advective_momentum_flux_Wu, args...)

@inline _advective_momentum_flux_Uv(i, j, k, ibg::IBG, args...) = conditional_flux_ffc(i, j, k, ibg, ibg, advective_momentum_flux_Uv, args...)
@inline _advective_momentum_flux_Vv(i, j, k, ibg::IBG, args...) = conditional_flux_ccc(i, j, k, ibg, ibg, advective_momentum_flux_Vv, args...)
@inline _advective_momentum_flux_Wv(i, j, k, ibg::IBG, args...) = conditional_flux_cff(i, j, k, ibg, ibg, advective_momentum_flux_Wv, args...)

@inline _advective_momentum_flux_Uw(i, j, k, ibg::IBG, args...) = conditional_flux_fcf(i, j, k, ibg, ibg, advective_momentum_flux_Uw, args...)
@inline _advective_momentum_flux_Vw(i, j, k, ibg::IBG, args...) = conditional_flux_cff(i, j, k, ibg, ibg, advective_momentum_flux_Vw, args...)
@inline _advective_momentum_flux_Ww(i, j, k, ibg::IBG, args...) = conditional_flux_ccc(i, j, k, ibg, ibg, advective_momentum_flux_Ww, args...)

@inline _advective_tracer_flux_x(i, j, k, ibg::IBG, args...) = conditional_flux_fcc(i, j, k, ibg, ibg, advective_tracer_flux_x, args...)
@inline _advective_tracer_flux_y(i, j, k, ibg::IBG, args...) = conditional_flux_cfc(i, j, k, ibg, ibg, advective_tracer_flux_y, args...)
@inline _advective_tracer_flux_z(i, j, k, ibg::IBG, args...) = conditional_flux_ccf(i, j, k, ibg, ibg, advective_tracer_flux_z, args...)

#####
##### "Boundary-aware" interpolation
#####
##### Don't interpolate dead cells.
#####

@inline near_x_boundary(i, j, k, ibg, buffer) = any(ntuple(δ -> solid_cell(i - buffer - 1 + δ, j, k, ibg), Val(2buffer + 1)))
@inline near_y_boundary(i, j, k, ibg, buffer) = any(ntuple(δ -> solid_cell(i, j - buffer - 1 + δ, k, ibg), Val(2buffer + 1)))
@inline near_z_boundary(i, j, k, ibg, buffer) = any(ntuple(δ -> solid_cell(i, j, k - buffer - 1 + δ, ibg), Val(2buffer + 1)))

for bias in (:symmetric, :left_biased, :right_biased)

    for (d, ξ) in enumerate((:x, :y, :z))

        code = [:ᵃ, :ᵃ, :ᵃ]

        for loc in (:ᶜ, :ᶠ)
            code[d] = loc
            second_order_interp = Symbol(:ℑ, ξ, code...)
            interp = Symbol(bias, :_interpolate_, ξ, code...)
            alt_interp = Symbol(:_, interp)

            near_boundary = Symbol(:near_, ξ, :_boundary)

            # Conditional high-order interpolation in Bounded directions
            @eval begin
                import Oceananigans.Advection: $alt_interp
                using Oceananigans.Advection: $interp

                @inline $alt_interp(i, j, k, ibg::ImmersedBoundaryGrid, scheme, ψ) =
                    ifelse($near_boundary(i, j, k, ibg, boundary_buffer(scheme)),
                           $second_order_interp(i, j, k, ibg.grid, ψ),
                           $interp(i, j, k, ibg.grid, scheme, ψ))
            end
        end
    end
end


text:
function DSSTTesseral(arg0::Frame, arg1::jdouble, arg2::UnnormalizedSphericalHarmonicsProvider)
    return DSSTTesseral((Frame, jdouble, UnnormalizedSphericalHarmonicsProvider), arg0, arg1, arg2)
end

function DSSTTesseral(arg0::Frame, arg1::jdouble, arg2::UnnormalizedSphericalHarmonicsProvider, arg3::jint, arg4::jint, arg5::jint, arg6::jint, arg7::jint, arg8::jint, arg9::jint)
    return DSSTTesseral((Frame, jdouble, UnnormalizedSphericalHarmonicsProvider, jint, jint, jint, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function get_events_detectors(obj::DSSTTesseral)
    return jcall(obj, "getEventsDetectors", Vector{EventDetector}, ())
end

function get_field_events_detectors(obj::DSSTTesseral, arg0::Field)
    return jcall(obj, "getFieldEventsDetectors", Vector{FieldEventDetector}, (Field,), arg0)
end

function get_mean_element_rate(obj::DSSTTesseral, arg0::FieldSpacecraftState, arg1::FieldAuxiliaryElements, arg2::Vector{RealFieldElement})
    return jcall(obj, "getMeanElementRate", Vector{RealFieldElement}, (FieldSpacecraftState, FieldAuxiliaryElements, Vector{RealFieldElement}), arg0, arg1, arg2)
end

function get_mean_element_rate(obj::DSSTTesseral, arg0::SpacecraftState, arg1::AuxiliaryElements, arg2::Vector{jdouble})
    return jcall(obj, "getMeanElementRate", Vector{jdouble}, (SpacecraftState, AuxiliaryElements, Vector{jdouble}), arg0, arg1, arg2)
end

function get_parameters_drivers(obj::DSSTTesseral)
    return jcall(obj, "getParametersDrivers", Vector{ParameterDriver}, ())
end

function initialize(obj::DSSTTesseral, arg0::AuxiliaryElements, arg1::PropagationType, arg2::Vector{jdouble})
    return jcall(obj, "initialize", List, (AuxiliaryElements, PropagationType, Vector{jdouble}), arg0, arg1, arg2)
end

function initialize(obj::DSSTTesseral, arg0::FieldAuxiliaryElements, arg1::PropagationType, arg2::Vector{RealFieldElement})
    return jcall(obj, "initialize", List, (FieldAuxiliaryElements, PropagationType, Vector{RealFieldElement}), arg0, arg1, arg2)
end

function register_attitude_provider(obj::DSSTTesseral, arg0::AttitudeProvider)
    return jcall(obj, "registerAttitudeProvider", void, (AttitudeProvider,), arg0)
end

function update_short_period_terms(obj::DSSTTesseral, arg0::Vector{RealFieldElement}, arg1::Vector{FieldSpacecraftState})
    return jcall(obj, "updateShortPeriodTerms", void, (Vector{RealFieldElement}, Vector{FieldSpacecraftState}), arg0, arg1)
end

function update_short_period_terms(obj::DSSTTesseral, arg0::Vector{jdouble}, arg1::Vector{SpacecraftState})
    return jcall(obj, "updateShortPeriodTerms", void, (Vector{jdouble}, Vector{SpacecraftState}), arg0, arg1)
end


text:
"""
TLeakyGlaescher: delta-rule updates of the transition matrix
as in Glaescher et al., 2010
"""
struct TLeakyGlaescher <: TPs1a0s0
    nr_states::Int
    nr_actions::Int
    η::Float64
    SPE::Array{Float64, 1}
    Ps1a0s0::Array{Dict{Tuple{Int, Int}, Float64}, 1}
    terminalstates::Array{Int,1}
end
function TLeakyGlaescher(; nr_states = 10, nr_actions = 4, η = .9)
    SPE = [0.]
    Ps1a0s0 = [Dict{Tuple{Int, Int}, Float64}() for _ in 1:nr_states]
    [Ps1a0s0[sprime][(a, s)] = 1. /nr_states for sprime in 1:nr_states for a in 1:nr_actions for s in 1:nr_states]
    TLeakyGlaescher(nr_states, nr_actions, η, SPE, Ps1a0s0, Int[])
end
export TLeakyGlaescher
""" X[t] = etaleak * X[t-1] + etaleak * I[.] """
function updatet!(learnerT::TLeakyGlaescher, s0, a0, s1, done)
    calcSPE!(learnerT, s0, a0, s1)
    leaka0s0!(learnerT, s0, a0, s1, done)
    computeterminalPs1a0s0!(learnerT, s1, done)
end
export updatet!
function leaka0s0!(learnerT::TLeakyGlaescher, s0, a0, s1, done)
    tsasprime = learnerT.Ps1a0s0[s1][(a0, s0)] + learnerT.η * learnerT.SPE[1]
    for sprime in 1:learnerT.nr_states
        learnerT.Ps1a0s0[sprime][(a0, s0)] *= (1. - learnerT.η)
    end
    learnerT.Ps1a0s0[s1][(a0, s0)] = copy(tsasprime)
end
function calcSPE!(learnerT::TLeakyGlaescher, s0, a0, s1)
    learnerT.SPE[1] = 1. - learnerT.Ps1a0s0[s1][(a0, s0)]
end


text:
abstract type Game end

const RoomGames = Dict{String, Game}()

# Returns list of users in game g
# To be implemented by all Game subtypes
function inGameUsers(g::Game)
    @warn("Unimplemented inGameUsers method of $(typeof(g)) : $g")
    return []
end

# Clean-up Game References
function removeGameRef(roomname)
    if haskey(RoomGames, roomname)
        delete!(RoomGames, roomname)
    end

    if haskey(Disconnected, roomname)
        delete!(Disconnected, roomname)
    end

    msg = Dict{String, Any}()
    msg["responsetype"] = "game-stopped"
    return broadcastmsg(roomname, JSON.json(msg))
end


text:
const printbranch = false

Cassette.@context HasBranchingCtx

function Cassette.overdub(ctx::HasBranchingCtx, f, args...)
    if Cassette.canrecurse(ctx, f, args...)
        return Cassette.recurse(ctx, f, args...)
    else
        return Cassette.fallback(ctx, f, args...)
    end
end

for (mod, f, n) in DiffRules.diffrules()
    isdefined(@__MODULE__, mod) || continue
    @eval Cassette.overdub(::HasBranchingCtx, f::Core.Typeof($mod.$f), x::Vararg{Any, $n}) = f(x...)
end

function _pass(::Type{<:HasBranchingCtx}, reflection::Cassette.Reflection)
    ir = reflection.code_info

    if any(x -> isa(x, GotoIfNot), ir.code)
        printbranch && ccall(:jl_safe_printf, Cvoid, (Cstring,), "GotoIfNot detected in $(reflection.method)\nir = $ir\n")
        Cassette.insert_statements!(
            ir.code, ir.codelocs,
            (stmt, i) -> i == 1 ? 3 : nothing,
            (stmt, i) -> Any[
                Expr(:call, Expr(:nooverdub, GlobalRef(Base, :getfield)), Expr(:contextslot), QuoteNode(:metadata)),
                Expr(:call, Expr(:nooverdub, GlobalRef(Base, :setindex!)), SSAValue(1), true, QuoteNode(:has_branching)),
                stmt,
            ],
        )
        Cassette.insert_statements!(
            ir.code, ir.codelocs,
            (stmt, i) -> i > 2 && isa(stmt, Expr) ? 1 : nothing,
            (stmt, i) -> begin
                callstmt = Meta.isexpr(stmt, :(=)) ? stmt.args[2] : stmt
                Meta.isexpr(stmt, :call) || Meta.isexpr(stmt, :invoke) || return Any[stmt]
                callstmt = Expr(callstmt.head, Expr(:nooverdub, callstmt.args[1]), callstmt.args[2:end]...)
                return Any[
                    Meta.isexpr(stmt, :(=)) ? Expr(:(=), stmt.args[1], callstmt) : callstmt,
                ]
            end,
        )
    end
    return ir
end

const pass = Cassette.@pass _pass

function hasbranching(f, x...)
    metadata = Dict(:has_branching => false)
    Cassette.overdub(Cassette.disablehooks(HasBranchingCtx(; pass, metadata)), f, x...)
    return metadata[:has_branching]
end

Cassette.overdub(::HasBranchingCtx, ::typeof(+), x...) = +(x...)
Cassette.overdub(::HasBranchingCtx, ::typeof(*), x...) = *(x...)
Cassette.overdub(::HasBranchingCtx, ::typeof(Base.materialize), x...) = Base.materialize(x...)
Cassette.overdub(::HasBranchingCtx, ::typeof(Base.literal_pow), x...) = Base.literal_pow(x...)
Cassette.overdub(::HasBranchingCtx, ::typeof(Base.getindex), x...) = Base.getindex(x...)
Cassette.overdub(::HasBranchingCtx, ::typeof(Core.Typeof), x...) = Core.Typeof(x...)


text:
@testset "Day4.jl" begin
    # Write your tests here.
    lines =["ecl:gry pid:860033327 eyr:2020 hcl:#fffffd",
    "byr:1937 iyr:2017 cid:147 hgt:183cm"]

    pp = Dict(
        "ecl"=>"gry",
        "pid"=>"860033327",
        "eyr"=>"2020",
        "hcl"=>"#fffffd",
        "byr"=>"1937",
        "iyr"=>"2017",
        "cid"=>"147",
        "hgt"=>"183cm"
        )

    @test AdventOfCode2020.topassports(lines) == [pp]

    @test solve(Day4_1, "$(@__DIR__)/assets/day4_test1.txt") == 2

    lines_invalid = ["eyr:1972 cid:100",
    "hcl:#18171d ecl:amb hgt:170 pid:186cm iyr:2018 byr:1926",
    "",
    "iyr:2019",
    "hcl:#602927 eyr:1967 hgt:170cm",
    "ecl:grn pid:012533040 byr:1946",
    "",
    "hcl:dab227 iyr:2012",
    "ecl:brn hgt:182cm pid:021572410 eyr:2020 byr:1992 cid:277",
    "",
    "hgt:59cm ecl:zzz",
    "eyr:2038 hcl:74454a iyr:2023",
    "pid:3556412378 byr:2007"]

    @test solve(Day4_2, lines_invalid[1:2]) == 0
    @test solve(Day4_2, lines_invalid[4:6]) == 0
    @test solve(Day4_2, lines_invalid[8:9]) == 0
    @test solve(Day4_2, lines_invalid[11:13]) == 0
    @test solve(Day4_2, lines_invalid) == 0

    lines_valid = ["pid:087499704 hgt:74in ecl:grn iyr:2012 eyr:2030 byr:1980",
    "hcl:#623a2f",
    "",
    "eyr:2029 ecl:blu cid:129 byr:1989",
    "iyr:2014 pid:896056539 hcl:#a97842 hgt:165cm",
    "",
    "hcl:#888785",
    "hgt:164cm byr:2001 iyr:2015 cid:88",
    "pid:545766238 ecl:hzl",
    "eyr:2022",
    "",
    "iyr:2010 hgt:158cm hcl:#b6652a ecl:blu byr:1944 eyr:2021 pid:093154719"]

    @test solve(Day4_2, lines_valid[1:2]) == 1
    @test solve(Day4_2, lines_valid[4:5]) == 1
    @test solve(Day4_2, lines_valid[7:10]) == 1
    @test solve(Day4_2, lines_valid[12:12]) == 1
    @test solve(Day4_2, lines_valid) == 4
end