text:
using Serialization
using .OSRM
using ProgressBars

struct TransitNetwork
    stops::Vector{Stop}
    stopidx_for_id::Dict{String,Int64}

    routes::Vector{Route}
    routeidx_for_id::Dict{String, Int64}

    services::Vector{Service}
    serviceidx_for_id::Dict{String, Int64}

    trips::Vector{Trip}
    tripidx_for_id::Dict{String, Int64}

    patterns::Vector{TripPattern}

    # todo could this have less indirection?
    transfers::Vector{Vector{Transfer}}

    # NB could also put this inside stop and pattern structs, but would prevent stop packing
    patterns_for_stop::Vector{Vector{Int64}}
    trips_for_pattern::Vector{Vector{Int64}}
end

TransitNetwork() = TransitNetwork(
    Vector{Stop}(),
    Dict{String,Int64}(),
    Vector{Route}(),
    Dict{String,Int64}(),
    Vector{Service}(),
    Dict{String,Int64}(),
    Vector{Trip}(),
    Dict{String,Int64}(),
    Vector{TripPattern}(),
    Vector{Vector{Transfer}}(),
    Vector{Vector{Int64}}(),
    Vector{Vector{Int64}}()
    )

# ensure that stop times are ordered by stop stop_sequence in each trip
function sort_stoptimes!(net::TransitNetwork)
    for trip in net.trips
        sort!(trip.stop_times, by=st -> st.stop_sequence)
    end
end

function hash_tp(stop_times::Vector{StopTime}, service::Int64)
    hash::Int64 = 0

    for (i, st) in enumerate(stop_times)
        hash += st.stop * primes[i % length(primes)]
    end

    hash += service * primes[length(primes)]

    return hash
end

# create trip patterns
function find_trip_patterns!(net::TransitNetwork)
    # maps from trip pattern hash to trip patterns with that hash
    tp_hashes = Dict{Int64, Vector{Tuple{Int64, TripPattern}}}()

    trips_with_patterns = Vector{Trip}()
    sizehint!(trips_with_patterns, length(net.trips))

    for trip in net.trips
        hash = hash_tp(trip.stop_times, trip.service)
        if !haskey(tp_hashes, hash)
            tp_hashes[hash] = Vector{Tuple{Int64, TripPattern}}()
        end

        stops::Vector{Int64} = map(st -> st.stop, trip.stop_times)

        # check for a matching trip pattern
        found_pattern = false
        for (tpidx, tp) in tp_hashes[hash]
            # short circuit and means vector equality won't be tested unless vectors are same length
            if ((tp.service == trip.service) && length(tp.stops) == length(stops) && all(tp.stops .== stops))
                found_pattern = true
                new_trip = Trip(
                    trip.stop_times,
                    trip.route,
                    trip.service,
                    tpidx
                )
                push!(trips_with_patterns, new_trip)
            end
        end

        if !found_pattern
            # create a new trip pattern
            tp = TripPattern(stops, trip.service)
            push!(net.patterns, tp)
            tpidx = length(net.patterns)
            push!(tp_hashes[hash], (tpidx, tp))
        end
    end

    # clear the entire trips array
    empty!(net.trips)
    append!(net.trips, trips_with_patterns)

    @info "created $(length(net.patterns)) trip patterns"
end

# find transfers based on crow-flies distance
function find_transfers_distance!(net::TransitNetwork, max_distance_meters::Real)
    empty!(net.transfers)
    total_transfers = 0
    sizehint!(net.transfers, length(net.stops))
    for stop in net.stops
        # find nearby stops
        # could use a spatial index if this is slow
        lat_diff = meters_to_degrees_lat(max_distance_meters)
        lon_diff = meters_to_degrees_lon(max_distance_meters, stop.stop_lat)

        # bbox query for nearby stops
        candidate_stops = filter(t -> (
            (t[2].stop_lat > stop.stop_lat - lat_diff) &&
            (t[2].stop_lat < stop.stop_lat + lat_diff) &&
            (t[2].stop_lon > stop.stop_lon - lon_diff) &&
            (t[2].stop_lon < stop.stop_lon + lon_diff) &&
            (t[2] !== stop)
            ), collect(enumerate(net.stops)))

        candidate_xfers = map(t -> Transfer(
            t[1],
            distance_meters(stop.stop_lat, stop.stop_lon, t[2].stop_lat, t[2].stop_lon)
        ), candidate_stops)

        # some overselection possible in corners of bbox
        xfers = filter(xfer -> xfer.distance_meters <= max_distance_meters, candidate_xfers)
        total_transfers += length(xfers)
        push!(net.transfers, xfers)
    end

    @assert length(net.transfers) == length(net.stops)

    @info "Created $total_transfers transfers from $(length(net.stops)) stops"
end

# Find transfers through the street network using OSRM
function find_transfers_osrm!(net::TransitNetwork, osrm::OSRMInstance, max_distance_meters::Real)
    empty!(net.transfers)
    total_transfers = 0
    sizehint!(net.transfers, length(net.stops))
    for stop in ProgressBar(net.stops)
        # find nearby stops
        # could use a spatial index if this is slow
        lat_diff = meters_to_degrees_lat(max_distance_meters)
        lon_diff = meters_to_degrees_lon(max_distance_meters, stop.stop_lat)

        # bbox query for nearby stops
        candidate_stops = collect(filter(t -> (
            (t[2].stop_lat > stop.stop_lat - lat_diff) &&
            (t[2].stop_lat < stop.stop_lat + lat_diff) &&
            (t[2].stop_lon > stop.stop_lon - lon_diff) &&
            (t[2].stop_lon < stop.stop_lon + lon_diff) &&
            (t[2] !== stop)
            ), collect(enumerate(net.stops))))

        # destinations has same order as candidate_stops
        destinations = collect(map(s -> Coordinate(s[2].stop_lat, s[2].stop_lon), candidate_stops))

        dmat = distance_matrix(osrm, [Coordinate(stop.stop_lat, stop.stop_lon)], destinations)

        xfers = Vector{Transfer}()

        for dest_idx in 1:length(destinations)
            # convert the index in destinations back to the index in net.stops
            stop_idx = candidate_stops[dest_idx][1]
            netdist = dmat.distances[1, dest_idx]
            if netdist < max_distance_meters
                push!(xfers, Transfer(stop_idx, netdist))
            end
        end

        total_transfers += length(xfers)
        push!(net.transfers, xfers)
    end

    @assert length(net.transfers) == length(net.stops)
    @info "Created $total_transfers transfers from $(length(net.stops)) stops"
end

function index_network!(net::TransitNetwork)
    # first, patterns for stops
    sizehint!(net.patterns_for_stop, length(net.stops))
    for stop in net.stops
        push!(net.patterns_for_stop, Vector{Int64}())
    end

    for (i, pattern) in enumerate(net.patterns)
        for stop in pattern.stops
            push!(net.patterns_for_stop[stop], i)
        end
    end

    # now, trips for pattern
    sizehint!(net.trips_for_pattern, length(net.patterns))
    for pattern in net.patterns
        push!(net.trips_for_pattern, Vector{Int64}())
    end

    for (i, trip) in enumerate(net.trips)
        push!(net.trips_for_pattern[trip.pattern], i)
    end
end

function _interpolate_segment_times(net::TransitNetwork, stop_times::Vector{StopTime})::Vector{StopTime}
    first_time = stop_times[1].departure_time
    last_time = stop_times[length(stop_times)].arrival_time

    distances = Vector{Float64}()
    sizehint!(distances, length(stop_times) - 1)

    for i in 1:(length(stop_times) - 1)
        fr = net.stops[stop_times[i].stop]
        to = net.stops[stop_times[i + 1].stop]
        push!(distances, distance_meters(
            fr.stop_lat,
            fr.stop_lon,
            to.stop_lat,
            to.stop_lon
        ))
    end

    distfrac = cumsum(distances) ./ sum(distances)

    timedelta = last_time - first_time

    new_stop_times = Vector{StopTime}()
    sizehint!(new_stop_times, length(stop_times))
    push!(new_stop_times, stop_times[1])

    for i in 2:(length(stop_times) - 1)
        interp_time::Int32 = first_time + round(distfrac[i - 1] * timedelta)
        push!(new_stop_times, StopTime(
            stop_times[i].stop,
            stop_times[i].stop_sequence,
            interp_time,
            interp_time
        ))
    end

    push!(new_stop_times, stop_times[length(stop_times)])

    @assert length(stop_times) == length(new_stop_times)

    return new_stop_times
end

function interpolate_stoptimes!(net::TransitNetwork)
    for trip in net.trips
        last_bona_fide_stidx = 0
        # collect should avoid concurrent modification issues
        for (i, st) in collect(enumerate(trip.stop_times))
            is_bona_fide_stidx = false
            if (st.arrival_time == INT_MISSING && st.departure_time != INT_MISSING)
                st.arrival_time = st.departure_time
                is_bona_fide_stidx = true
            elseif (st.departure_time == INT_MISSING && st.arrival_time != INT_MISSING)
                st.departure_time = st.arrival_time
                is_bona_fide_stidx = true
            elseif (st.departure_time != INT_MISSING && st.arrival_time != INT_MISSING)
                # do nothing, continue to accumulate
                is_bona_fide_stidx = true
            end

            if is_bona_fide_stidx
                if last_bona_fide_stidx != i - 1
                    trip.stop_times[last_bona_fide_stidx:i] = _interpolate_segment_times(net, trip.stop_times[last_bona_fide_stidx:i])
                end
                last_bona_fide_stidx = i
            end
        end

        # check our work
        @assert all(diff(map(st -> st.arrival_time, trip.stop_times)) .>= 0)
        @assert all(diff(map(st -> st.departure_time, trip.stop_times)) .>= 0)
    end
end

function save_network(network::TransitNetwork, filename::String)
    serialize(filename, network)
end

function load_network(filename::String)::TransitNetwork
    return deserialize(filename)
end


text:
#! /usr/bin/julia

# Rosetta Code, Bitmap

type Color
    r::Uint8
    g::Uint8
    b::Uint8
end

type Image
    pic::Array{Color,2}
end

Image(w::Int, h::Int) = Image(Array(Color, w, h))

Base.fill!(a::Image, c::Color) = Base.fill!(a.pic, c)

function splat!(a::Image, x::Int, y::Int, c::Color)
    a.pic[x, y] = c
    nothing
end

function color(a::Image, x::Int, y::Int)
    a.pic[x, y]
end

function showpixel(a::Image, x::Int, y::Int)
    c = color(a, x, y)
    hex(c.r, 2)*hex(c.g, 2)*hex(c.b, 2)
end

w = 5
h = 7
a = Image(w, h)

purple = Color(0xff, 0, 0xff)
green = Color(0, 0xff, 0)

fill!(a, purple)
splat!(a, 2, 3, green)

for i in 1:h
    for j in 1:w
        print(showpixel(a, j, i), " ")
    end
    println()
end

ext: jl

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
######################
# Logit and Logistic #
######################
using StatsFuns: logit, logistic

struct Logit{T<:Real} <: Bijector{0}
    a::T
    b::T
end

(b::Logit)(x) = @. logit((x - b.a) / (b.b - b.a))
(ib::Inversed{<:Logit{<:Real}})(y) = @. (ib.orig.b - ib.orig.a) * logistic(y) + ib.orig.a

logabsdetjac(b::Logit{<:Real}, x) = @. - log((x - b.a) * (b.b - x) / (b.b - b.a))


text:
"""
	hasnode(network::AbstractNetwork, idx::Integer)

Check whether the node exists in the network
"""
hasnode(network::AbstractNetwork, idx::Integer) = (0 < idx <= network.N)

"""
	hasconnection(network::AbstractNetwork, idx_origin::Integer, idx_dest::Integer)

Check whether there exists a connections from node `idx_origin` to `idx_dest`
"""
hasconnection(network::AbstractNetwork, idx_origin::Integer, idx_dest::Integer) =
	Bool(adjMat(network, sparse=true)[idx_dest, idx_origin])

"""
	allEdges(network::AbstractNetwork; first_index::Integer=1)
	allEdges(adjacencyMatrix::AbstractMatrix, directed::Bool; first_index::Integer=1)

Return an array with all existing connections as a tuple (origin,destination).

# Arguments
- idx_origin: Which of 0 or 1 should be used as the first index.
- both_directions: (undirected only) if true will return both directions of an undirected
connections. i.e.: i -> j and j -> i.
"""
function allEdges(mat::AbstractMatrix, directed::Bool; first_index::Integer=1, both_directions::Bool=false)
	(first_index in [1,0]) || throw(ArgumentError("first_index must be 0 or 1"))
	@assert size(mat, 1) == size(mat, 2)

	N = size(mat, 1)

	if directed || both_directions
		edges = [(i,j) for i in 1:N for j in 1:N if mat[j,i] != 0]
	else
		edges = [(i,j) for i in 1:N for j in i+1:N if mat[j,i] != 0]
	end

	(first_index == 0) && (edges = [(e[1]-1, e[2]-1) for e in edges])
	return edges
end

allEdges(network::AbstractNetwork; first_index::Integer=1, both_directions::Bool=false) =
	allEdges(adjMat(network, sparse=true), isdirected(network); first_index=first_index, both_directions=both_directions)

"""
	neighbors(network::AbstractNetwork, idx_node::Integer; directed_behaviour::Symbol=:any)

Return a vector with the indices of the neighbors of the passed node.

If the network is directed the `directed_behaviour` can be used to specify
the considered neighbors.
`:origin` for nodes whose connections arive at the selected node,
`:destination` for nodes in which connections from the selected node arive,
and `:any` for the union of both previous groups.
"""
function neighbors(network::AbstractNetwork, idx_node::Integer; directed_behaviour::Symbol=:any)
	hasnodeOrError(network, idx_node)
	mat = adjMat(network, sparse=true)

	if isdirected(network)
		(directed_behaviour in [:origin, :destination, :any]) ||
			throw(ArgumentError("directed_behaviour must be one of: :origin, :destination, :any"))

		(directed_behaviour == :origin) && (return [i for i in 1:network.N if mat[idx_node, i] != 0])
		(directed_behaviour == :destination) && (return [i for i in 1:network.N if mat[i, idx_node] != 0])
		(directed_behaviour == :any) &&
			(return [i for i in 1:network.N if mat[i, idx_node] != 0 || mat[idx_node, i] != 0])
	else
		return [i for i in 1:network.N if mat[i, idx_node] != 0]
	end
end

"""
	adjVetToMat(vet::Vector, N::Integer)

Calculate the adjacency matrix based on the adjacency vector.

See also: `adjMatToVet`
"""
function adjVetToMat(vet::Vector{<:Integer}, N::Integer; sparse::Bool=false)
	mat = sparse ? SparseArrays.spzeros(Int, N, N) : BitArray(0 for i in 1:N, j in 1:N)

	let i = 0
		for idx in vet
			j = idx - i*N

			while j >= N
				i += 1
				j -= N
			end

			mat[i+1, j+1] = 1
		end
	end

	return mat
end

"""
	adjMatToVet(mat::AbstractMatrix)

Calculate the adjacency vector based on the adjacency matrix.

See also: `adjVetToMat`
"""
function adjMatToVet(mat::SparseArrays.AbstractSparseMatrix)
	vet = Array{Int}(undef, length(mat.rowval))
	idx_vet = 1
	colptr = mat.colptr
	rowval = mat.rowval

	for idx_col in eachindex(mat.colptr[1:end-1])
		idxs_row = @view rowval[colptr[idx_col] : colptr[idx_col+1]-1]

		for idx_row in idxs_row
			vet[idx_vet] = (idx_row-1) * size(mat,1) + (idx_col-1)
			idx_vet += 1
		end
	end

	return sort(vet)
end

function adjMatToVet(mat::BitArray{2})
	@assert size(mat,1) == size(mat,2)

	N = size(mat)[1]
	vetPos = 1
	vet = Vector{Int}(undef, sum(mat))

	for i in 1:N, j in 1:N
		if mat[i,j] != 0
			vet[vetPos] = (i-1)*N + (j-1)
			vetPos += 1
		end
	end

	return vet
end

adjMatToVet(mat::Array{Bool,2}) = adjMatToVet(BitArray(mat))
function adjMatToVet(mat::AbstractArray{<:Real,2})
	try
		BitArray(mat)
	catch e
		throw(ArgumentError("Matrix must be composed only of ones and zeros!"))
	end

	return adjMatToVet(BitArray(mat))
end

"""
	equivalentRandomNetwork(network::AbstractNetwork; seed, numReshuffle, maxTentatives)

Return a random network with the same degree distribution of the original.

# Keyworkd Arguments
- seed (Integer, default=-1): Seed for random creation
- numReshuffle (Integer, default=10): Number of reshuffling cycles.
Each cycle makes `numconnections(network)` rewires.
- maxTentatives (Integer, default=numconnections(network) * numReshuffle * 300) :
Maximum number of rewiring tentatives before ending the algorithm.
"""
function equivalentRandomNetwork(
	network::AbstractNetwork;
	seed::Integer = -1,
	numReshuffle::Integer = 10,
	maxTentatives::Integer = (numconnections(network) * numReshuffle * 300)
)
	(seed < 0) && (seed = rand(1:99999999))
	rng = Random.MersenneTwister(seed)
	numShuffle = numReshuffle * numconnections(network)

	edges = allEdges(network, both_directions=true)
	mat = adjMat(network, sparse=true)

	numRewires = 0
	numTentatives = 0

	while numRewires < numShuffle && numTentatives < maxTentatives
		numTentatives += 1

		# Choose random pair of edges
		edge1 = rand(rng, edges)
		edge2 = rand(rng, edges)
		while edge1 == edge2 || edge1[1] == edge2[2] || edge2[1] == edge1[2]
			edge2 = rand(rng, edges)
		end

		# Check if the rewiring is valid
		((edge1[1], edge2[2]) in edges) && continue
		((edge2[1], edge1[2]) in edges) && continue

		# Rewire
		mat[edge1[2], edge1[1]] = 0
		mat[edge2[2], edge2[1]] = 0
		mat[edge1[2], edge2[1]] = 1
		mat[edge2[2], edge1[1]] = 1

		edges = filter(x -> x != edge1 && x != edge2, edges)
		edges = vcat(edges, [(edge1[1], edge2[2]), (edge2[1], edge1[2])])

		numRewires += 1
	end

	return isa(mat, SparseArrays.AbstractSparseMatrix) ?
		CustomNetwork(SparseArrays.dropzeros(mat), directed=true) :
		CustomNetwork(mat, directed=true)
end

"""
	equivalentLatticeNetwork(network::AbstractNetwork; seed, numReshuffle, maxTentatives)

Return a network with the same degree distribution of the original but with higher
clustering coefficient and connections closer to the diagonal.

# Keyworkd Arguments
- seed (Integer, default=-1): Seed for random creation
- numReshuffle (Integer, default=10): Number of reshuffling cycles.
Each cycle makes `numconnections(network)` rewires.
- maxTentatives (Integer, default=numconnections(network) * numReshuffle * 300) :
Maximum number of rewiring tentatives before ending the algorithm.
"""
function equivalentLatticeNetwork(
	network::AbstractNetwork;
	seed::Integer = -1,
	numReshuffle::Integer = 10,
	maxTentatives::Integer = (numconnections(network) * numReshuffle * 300)
)
	isregular(network) && (return network)

	function distToDiag(i::Int, j::Int)
		dist = abs(i - j)
		return min(dist, network.N - dist)
	end

	(seed < 0) && (seed = rand(1:99999999))
	rng = Random.MersenneTwister(seed)
	numShuffle = numReshuffle * numconnections(network)

	edges = allEdges(network, both_directions=true)
	mat = adjMat(network, sparse=true)
	newnet = CustomNetwork(mat)

	numTriedNets = 0

	while clusteringcoefficient(newnet) <= clusteringcoefficient(network) && numTriedNets < 10
		numRewires = 0
		numTentatives = 0
		numTriedNets += 1

		while numRewires < numShuffle && numTentatives < maxTentatives
			numTentatives += 1

			# Choose random pair of edges
			edge1 = rand(rng, edges)
			edge2 = rand(rng, edges)
			while edge1 == edge2 || edge1[1] == edge2[2] || edge2[1] == edge1[2]
				edge2 = rand(rng, edges)
			end

			# Check if the rewiring is valid
			((edge1[1], edge2[2]) in edges) && continue
			((edge2[1], edge1[2]) in edges) && continue

			# Check if the new matrix will be closer to diagonal
			dist_old = distToDiag(edge1[1], edge1[2]) + distToDiag(edge2[1], edge2[2])
			dist_new = distToDiag(edge1[1], edge2[2]) + distToDiag(edge2[1], edge1[2])
			(dist_new >= dist_old) && continue

			# Rewire
			mat[edge1[2], edge1[1]] = 0
			mat[edge2[2], edge2[1]] = 0
			mat[edge1[2], edge2[1]] = 1
			mat[edge2[2], edge1[1]] = 1

			edges = filter(x -> x != edge1 && x != edge2, edges)
			edges = vcat(edges, [(edge1[1], edge2[2]), (edge2[1], edge1[2])])

			numRewires += 1
		end

		newnet = isa(mat, SparseArrays.AbstractSparseMatrix) ?
			CustomNetwork(SparseArrays.dropzeros(mat), directed=true) :
			CustomNetwork(mat, directed=true)
	end

	return newnet
end

"""
	saveAdjMat(network::AbstractNetwork, filepath::String)

Write the adjacency matrix of a network to the specified files.

See also: `saveAdjVet`
"""
saveAdjMat(network::AbstractNetwork, filepath::String) =
	DelimitedFiles.writedlm(filepath, Int.(adjMat(network)))

"""
	saveAdjVet(network::AbstractNetwork, filepath::String)

Write the adjacency vector of a network to the specified files.

See also: `saveAdjMat`
"""
saveAdjVet(network::AbstractNetwork, filepath::String) =
	DelimitedFiles.writedlm(filepath, adjVet(network))


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

using Markdown
using InteractiveUtils

# ╔═╡ f1b53f4d-6d7c-42b8-bd4b-b77e3b3a1055
using PlutoTest

# ╔═╡ e9086930-afeb-11eb-35ad-dba62b2459a9
md"""
# [Day 5: Doesn't He Have Intern-Elves For This?](https://adventofcode.com/2015/day/5)
Santa needs help figuring out which strings in his text file are naughty or nice.

A nice string is one with all of the following properties:

-   It contains at least three vowels (`aeiou` only), like `aei`, `xazegov`, or `aeiouaeiouaeiou`.
-   It contains at least one letter that appears twice in a row, like `xx`, `abcdde` (`dd`), or `aabbccdd` (`aa`, `bb`, `cc`, or `dd`).
-   It does not contain the strings `ab`, `cd`, `pq`, or `xy`, even if they are part of one of the other requirements.

For example:

-   `ugknbfddgicrmopn` is nice because it has at least three vowels (`u...i...o...`), a double letter (`...dd...`), and none of the disallowed substrings.
-   `aaa` is nice because it has at least three vowels and a double letter, even though the letters used by different rules overlap.
-   `jchzalrnumimnmhp` is naughty because it has no double letter.
-   `haegwjzuvuyypxyu` is naughty because it contains the string `xy`.
-   `dvszwmarrgswjxmb` is naughty because it contains only one vowel.

**How many strings are nice?**
"""

# ╔═╡ c7430382-47c7-4f2e-a260-fa4de8a8648c
puzzle_input = "zgsnvdmlfuplrubt
vlhagaovgqjmgvwq
ffumlmqwfcsyqpss
zztdcqzqddaazdjp
eavfzjajkjesnlsb
urrvucyrzzzooxhx
xdwduffwgcptfwad
orbryxwrmvkrsxsr
jzfeybjlgqikjcow
mayoqiswqqryvqdi
iiyrkoujhgpgkcvx
egcgupjkqwfiwsjl
zbgtglaqqolttgng
eytquncjituzzhsx
dtfkgggvqadhqbwb
zettygjpcoedwyio
rwgwbwzebsnjmtln
esbplxhvzzgawctn
vnvshqgmbotvoine
wflxwmvbhflkqxvo
twdjikcgtpvlctte
minfkyocskvgubvm
sfxhhdhaopajbzof
sofkjdtalvhgwpql
uqfpeauqzumccnrc
tdflsbtiiepijanf
dhfespzrhecigzqb
xobfthcuuzhvhzpn
olgjglxaotocvrhw
jhkzpfcskutwlwge
zurkakkkpchzxjhq
hekxiofhalvmmkdl
azvxuwwfmjdpjskj
arsvmfznblsqngvb
ldhkzhejofreaucc
adrphwlkehqkrdmo
wmveqrezfkaivvaw
iyphmphgntinfezg
blomkvgslfnvspem
cgpaqjvzhbumckwo
ydhqjcuotkeyurpx
sbtzboxypnmdaefr
vxrkhvglynljgqrg
ttgrkjjrxnxherxd
hinyfrjdiwytetkw
sufltffwqbugmozk
tohmqlzxxqzinwxr
jbqkhxfokaljgrlg
fvjeprbxyjemyvuq
gmlondgqmlselwah
ubpwixgxdloqnvjp
lxjfhihcsajxtomj
qouairhvrgpjorgh
nloszcwcxgullvxb
myhsndsttanohnjn
zjvivcgtjwenyilz
qaqlyoyouotsmamm
tadsdceadifqthag
mafgrbmdhpnlbnks
aohjxahenxaermrq
ovvqestjhbuhrwlr
lnakerdnvequfnqb
agwpwsgjrtcjjikz
lhlysrshsmzryzes
xopwzoaqtlukwwdu
xsmfrfteyddrqufn
ohnxbykuvvlbbxpf
bbdlivmchvzfuhoc
vtacidimfcfyobhf
tinyzzddgcnmiabd
tcjzxftqcqrivqhn
vgnduqyfpokbmzim
revkvaxnsxospyow
ydpgwxxoxlywxcgi
wzuxupbzlpzmikel
nscghlafavnsycjh
xorwbquzmgmcapon
asmtiycegeobfxrn
eqjzvgkxgtlyuxok
mmjrskloposgjoqu
gceqosugbkvytfto
khivvoxkvhrgwzjl
qtmejuxbafroifjt
ttmukbmpoagthtfl
bxqkvuzdbehtduwv
gvblrpzjylanoggj
cltewhyjxdbmbtqj
fbkgedqvomdipklj
uxvuplhenqawfcjt
fkdjmayiawdkycva
gnloqfgbnibzyidh
kyzorvtopjiyyyqg
drckpekhpgrioblt
tvhrkmbnpmkkrtki
khaldwntissbijiz
aoojqakosnaxosom
xfptccznbgnpfyqw
moqdwobwhjxhtrow
chfwivedutskovri
gprkyalfnpljcrmi
pwyshpwjndasykst
xuejivogihttzimd
bugepxgpgahtsttl
zufmkmuujavcskpq
urybkdyvsrosrfro
isjxqmlxwtqmulbg
pxctldxgqjqhulgz
hclsekryiwhqqhir
hbuihpalwuidjpcq
ejyqcxmfczqfhbxa
xljdvbucuxnnaysv
irqceqtqwemostbb
anfziqtpqzqdttnz
cgfklbljeneeqfub
zudyqkuqqtdcpmuo
iuvhylvznmhbkbgg
mpgppmgfdzihulnd
argwmgcvqqkxkrdi
pdhrfvdldkfihlou
cbvqnjrvrsnqzfob
lkvovtsqanohzcmm
vxoxjdyoylqcnyzt
kurdpaqiaagiwjle
gwklwnazaxfkuekn
rbaamufphjsjhbdl
tzbrvaqvizhsisbd
pbcqlbfjvlideiub
hiwoetbfywaeddtx
fjirczxtuupfywyf
omeoegeyyospreem
ozbbpupqpsskvrjh
pzvcxkvjdiyeyhxa
odclumkenabcsfzr
npdyqezqdjqaszvm
yodkwzmrhtexfrqa
rjcmmggjtactfrxz
mioxfingsfoimual
aqskaxjjborspfaa
wientdsttkevjtkf
tdaswkzckmxnfnct
voucjhzvkkhuwoqk
boaaruhalgaamqmh
iufzxutxymorltvb
pfbyvbayvnrpijpo
obztirulgyfthgcg
ntrenvhwxypgtjwy
ephlkipjfnjfjrns
pkjhurzbmobhszpx
gqbnjvienzqfbzvj
wjelolsrbginwnno
votanpqpccxqricj
bxyuyiglnmbtvehi
qyophcjfknbcbjrb
anoqkkbcdropskhj
tcnyqaczcfffkrtl
rsvqimuqbuddozrf
meppxdrenexxksdt
tyfhfiynzwadcord
wayrnykevdmywycf
mhowloqnppswyzbu
tserychksuwrgkxz
xycjvvsuaxsbrqal
fkrdsgaoqdcqwlpn
vrabcmlhuktigecp
xgxtdsvpaymzhurx
ciabcqymnchhsxkc
eqxadalcxzocsgtr
tsligrgsjtrnzrex
qeqgmwipbspkbbfq
vzkzsjujltnqwliw
ldrohvodgbxokjxz
jkoricsxhipcibrq
qzquxawqmupeujrr
mizpuwqyzkdbahvk
suupfxbtoojqvdca
ywfmuogvicpywpwm
uevmznxmsxozhobl
vjbyhsemwfwdxfxk
iyouatgejvecmtin
tcchwpuouypllcxe
lgnacnphdiobdsef
uoxjfzmdrmpojgbf
lqbxsxbqqhpjhfxj
knpwpcnnimyjlsyz
fezotpoicsrshfnh
dkiwkgpmhudghyhk
yzptxekgldksridv
pckmzqzyiyzdbcts
oqshafncvftvwvsi
yynihvdywxupqmbt
iwmbeunfiuhjaaic
pkpkrqjvgocvaxjs
ieqspassuvquvlyz
xshhahjaxjoqsjtl
fxrrnaxlqezdcdvd
pksrohfwlaqzpkdd
ravytrdnbxvnnoyy
atkwaifeobgztbgo
inkcabgfdobyeeom
ywpfwectajohqizp
amcgorhxjcybbisv
mbbwmnznhafsofvr
wofcubucymnhuhrv
mrsamnwvftzqcgta
tlfyqoxmsiyzyvgv
ydceguvgotylwtea
btyvcjqhsygunvle
usquiquspcdppqeq
kifnymikhhehgote
ybvkayvtdpgxfpyn
oulxagvbavzmewnx
tvvpekhnbhjskzpj
azzxtstaevxurboa
nfmwtfgrggmqyhdf
ynyzypdmysfwyxgr
iaobtgubrcyqrgmk
uyxcauvpyzabbzgv
fbasfnwiguasoedc
mgmjoalkbvtljilq
szgkxiqkufdvtksb
xgfzborpavdmhiuj
hmuiwnsonvfgcrva
zolcffdtobfntifb
mvzgcsortkugvqjr
pbbpgraaldqvzwhs
zvsxegchksgnhpuv
kdpdboaxsuxfswhx
jdfggigejfupabth
tpeddioybqemyvqz
mxsntwuesonybjby
tzltdsiojfvocige
ubtdrneozoejiqrv
fusyucnhncoxqzql
nlifgomoftdvkpby
pyikzbxoapffbqjw
hzballplvzcsgjug
ymjyigsfehmdsvgz
vpqgyxknniunksko
ffkmaqsjxgzclsnq
jcuxthbedplxhslk
ymlevgofmharicfs
nyhbejkndhqcoisy
rjntxasfjhnlizgm
oqlnuxtzhyiwzeto
tntthdowhewszitu
rmxyoceuwhsvfcua
qpgsjzwenzbxyfgw
sumguxpdkocyagpu
ymfrbxwrawejkduu
hetgrtmojolbmsuf
qzqizpiyfasgttex
qnmoemcpuckzsshx
ddyqiihagcmnxccu
oirwxyfxxyktgheo
phpaoozbdogbushy
uctjdavsimsrnvjn
aurbbphvjtzipnuh
hpbtrubopljmltep
pyyvkthqfsxqhrxg
jdxaiqzkepxbfejk
ukgnwbnysrzvqzlw
lfkatkvcssnlpthd
ucsyecgshklhqmsc
rwdcbdchuahkvmga
rxkgqakawgpwokum
hbuyxeylddfgorgu
tbllspqozaqzglkz
rqfwizjlbwngdvvi
xuxduyzscovachew
kouiuxckkvmetvdy
ycyejrpwxyrweppd
trctlytzwiisjamx
vtvpjceydunjdbez
gmtlejdsrbfofgqy
jgfbgtkzavcjlffj
tyudxlpgraxzchdk
gyecxacqitgozzgd
rxaocylfabmmjcvt
tornfzkzhjyofzqa
kocjcrqcsvagmfqv
zfrswnskuupivzxb
cunkuvhbepztpdug
pmpfnmklqhcmrtmf
tfebzovjwxzumxap
xpsxgaswavnzkzye
lmwijdothmxclqbr
upqxhmctbltxkarl
axspehytmyicthmq
xdwrhwtuooikehbk
tpggalqsytvmwerj
jodysbwnymloeqjf
rxbazvwuvudqlydn
ibizqysweiezhlqa
uexgmotsqjfauhzp
ldymyvumyhyamopg
vbxvlvthgzgnkxnf
pyvbrwlnatxigbrp
azxynqididtrwokb
lwafybyhpfvoawto
ogqoivurfcgspytw
cinrzzradwymqcgu
sgruxdvrewgpmypu
snfnsbywuczrshtd
xfzbyqtyxuxdutpw
fmpvjwbulmncykbo
ljnwoslktrrnffwo
ceaouqquvvienszn
yjomrunrxjyljyge
xpmjsapbnsdnbkdi
uetoytptktkmewre
eixsvzegkadkfbua
afaefrwhcosurprw
bwzmmvkuaxiymzwc
gejyqhhzqgsrybni
gjriqsfrhyguoiiw
gtfyomppzsruhuac
ogemfvmsdqqkfymr
jgzbipsygirsnydh
zghvlhpjnvqmocgr
ngvssuwrbtoxtrka
ietahyupkbuisekn
gqxqwjizescbufvl
eiprekzrygkncxzl
igxfnxtwpyaamkxf
soqjdkxcupevbren
fspypobyzdwstxak
qstcgawvqwtyyidf
gsccjacboqvezxvd
bfsblokjvrqzphmc
srezeptvjmncqkec
opmopgyabjjjoygt
msvbufqexfrtecbf
uiaqweyjiulplelu
pbkwhjsibtwjvswi
xwwzstmozqarurrq
nytptwddwivtbgyq
ejxvsufbzwhzpabr
jouozvzuwlfqzdgh
gfgugjihbklbenrk
lwmnnhiuxqsfvthv
bzvwbknfmaeahzhi
cgyqswikclozyvnu
udmkpvrljsjiagzi
zzuhqokgmisguyna
ekwcdnjzuctsdoua
eueqkdrnzqcaecyd
lnibwxmokbxhlris
fdrbftgjljpzwhea
iabvuhhjsxmqfwld
qgogzkynrgejakta
mfcqftytemgnpupp
klvhlhuqhosvjuqk
gdokmxcgoqvzvaup
juududyojcazzgvr
fyszciheodgmnotg
yfpngnofceqfvtfs
cahndkfehjumwavc
dxsvscqukljxcqyi
cqukcjtucxwrusji
vevmmqlehvgebmid
ahswsogfrumzdofy
ftasbklvdquaxhxb
tsdeumygukferuif
ybfgbwxaaitpwryg
djyaoycbymezglio
trzrgxdjqnmlnzpn
rumwchfihhihpqui
ffrvnsgrnzemksif
oizlksxineqknwzd
cirqcprftpjzrxhk
zrhemeqegmzrpufd
kqgatudhxgzlgkey
syjugymeajlzffhq
nlildhmgnwlopohp
flcszztfbesqhnyz
ohzicmqsajyqptrw
ebyszucgozsjbelq
enxbgvvcuqeloxud
ubwnvecbsmhkxwuk
noifliyxvlkqphbo
hazlqpetgugxxsiz
ihdzoerqwqhgajzb
ivrdwdquxzhdrzar
synwycdvrupablib
mqkdjkntblnmtvxj
qmmvoylxymyovrnq
pjtuxskkowutltlq
gchrqtloggkrjciz
namzqovvsdipazae
yfokqhkmakyjzmys
iapxlbuoiwqfnozm
fbcmlcekgfdurqxe
ednzgtczbplwxjlq
gdvsltzpywffelsp
oaitrrmpqdvduqej
gseupzwowmuuibjo
dfzsffsqpaqoixhh
tclhzqpcvbshxmgx
cfqkptjrulxiabgo
iraiysmwcpmtklhf
znwjlzodhktjqwlm
lcietjndlbgxzjht
gdkcluwjhtaaprfo
vbksxrfznjzwvmmt
vpfftxjfkeltcojl
thrmzmeplpdespnh
yafopikiqswafsit
xxbqgeblfruklnhs
qiufjijzbcpfdgig
ikksmllfyvhyydmi
sknufchjdvccccta
wpdcrramajdoisxr
grnqkjfxofpwjmji
lkffhxonjskyccoh
npnzshnoaqayhpmb
fqpvaamqbrnatjia
oljkoldhfggkfnfc
ihpralzpqfrijynm
gvaxadkuyzgbjpod
onchdguuhrhhspen
uefjmufwlioenaus
thifdypigyihgnzo
ugqblsonqaxycvkg
yevmbiyrqdqrmlbw
bvpvwrhoyneorcmm
gbyjqzcsheaxnyib
knhsmdjssycvuoqf
nizjxiwdakpfttyh
nwrkbhorhfqqoliz
ynsqwvwuwzqpzzwp
yitscrgexjfclwwh
dhajwxqdbtrfltzz
bmrfylxhthiaozpv
frvatcvgknjhcndw
xlvtdmpvkpcnmhya
pxpemuzuqzjlmtoc
dijdacfteteypkoq
knrcdkrvywagglnf
fviuajtspnvnptia
xvlqzukmwbcjgwho
bazlsjdsjoeuvgoz
nslzmlhosrjarndj
menvuwiuymknunwm
uavfnvyrjeiwqmuu
yrfowuvasupngckz
taevqhlrcohlnwye
skcudnogbncusorn
omtnmkqnqedsajfv
yqmgsqdgsuysqcts
odsnbtyimikkbmdd
vuryaohxdvjllieb
dhaxldeywwsfamlo
opobvtchezqnxpak
pzfnegouvsrfgvro
rzkcgpxdslzrdktu
ksztdtqzxvhuryam
ctnqnhkcooqipgkh
pyqbbvrzdittqbgm
koennvmolejeftij
rvzlreqikqlgyczj
xrnujfoyhonzkdgd
mmsmhkxaiqupfjil
ypjwoemqizddvyfd
qgugcxnbhvgahykj
cviodlsrtimbkgmy
xbfbbechhmrjxhnw
psuipaoucfczfxkp
hdhwcpeuptgqqvim
gsxlruhjeaareilr
vgyqonnljuznyrhk
eewezahlumervpyu
iiolebrxfadtnigy
tdadlrodykrdfscn
ocvdtzjxrhtjurpo
gidljbuvuovkhhrf
qwfcpilbjwzboohd
xzohxonlezuiupbg
vslpbkkqgvgbcbix
pivzqrzfxosbstzn
fyqcfboevcqmbhhs
yqsrneacnlxswojx
heicqpxxyrwcbsjz
yzynmnnoumkmlbeh
bncadbjdvvmczylw
hlnjskgfzbgmigfn
fphpszymugpcykka
zbifcktanxpmufvy
saklpkhoyfeqbguy
nqtqfcfxmpivnjyo
locygrwerxlsvzqm
qqflecydqvlogjme
njklmixvgkzpgppf
ugzkpjwjflaswyma
lriousvkbeftslcy
nsvsauxzfbbotgmh
tblcpuhjyybrlica
hqwshxcilwtmxrsf
xojwroydfeoqupup
tikuzsrogpnohpib
layenyqgxdfggloc
nqsvjvbrpuxkqvmq
ivchgxkdlfjdzxmk
uoghiuosiiwiwdws
twsgsfzyszsfinlc
waixcmadmhtqvcmd
zkgitozgrqehtjkw
xbkmyxkzqyktmpfi
qlyapfmlybmatwxn
ntawlvcpuaebuypf
clhebxqdkcyndyof
nrcxuceywiklpemc
lmurgiminxpapzmq
obalwqlkykzflxou
huvcudpiryefbcye
zlxbddpnyuyapach
gqfwzfislmwzyegy
jhynkjtxedmemlob
hmrnvjodnsfiukex
pstmikjykzyavfef
wuwpnscrwzsyalyt
hksvadripgdgwynm
tvpfthzjleqfxwkh
xpmrxxepkrosnrco
qjkqecsnevlhqsly
jjnrfsxzzwkhnwdm
pehmzrzsjngccale
bsnansnfxduritrr
ejzxkefwmzmbxhlb
pceatehnizeujfrs
jtidrtgxopyeslzl
sytaoidnamfwtqcr
iabjnikomkgmyirr
eitavndozoezojsi
wtsbhaftgrbqfsmm
vvusvrivsmhtfild
qifbtzszfyzsjzyx
ifhhjpaqatpbxzau
etjqdimpyjxiuhty
fvllmbdbsjozxrip
tjtgkadqkdtdlkpi
xnydmjleowezrecn
vhcbhxqalroaryfn
scgvfqsangfbhtay
lbufpduxwvdkwhmb
tshipehzspkhmdoi
gtszsebsulyajcfl
dlrzswhxajcivlgg
kgjruggcikrfrkrw
xxupctxtmryersbn
hljjqfjrubzozxts
giaxjhcwazrenjzs
tyffxtpufpxylpye
jfugdxxyfwkzqmgv
kbgufbosjghahacw
xpbhhssgegmthwxb
npefofiharjypyzk
velxsseyxuhrpycy
sglslryxsiwwqzfw
susohnlpelojhklv
lfnpqfvptqhogdmk
vtcrzetlekguqyle
jlyggqdtamcjiuxn
olxxqfgizjmvigvl
cyypypveppxxxfuq
hewmxtlzfqoqznwd
jzgxxybfeqfyzsmp
xzvvndrhuejnzesx
esiripjpvtqqwjkv
xnhrwhjtactofwrd
knuzpuogbzplofqx
tihycsdwqggxntqk
xkfywvvugkdalehs
cztwdivxagtqjjel
dsaslcagopsbfioy
gmowqtkgrlqjimbl
ctcomvdbiatdvbsd
gujyrnpsssxmqjhz
nygeovliqjfauhjf
mmgmcvnuppkbnonz
bhipnkoxhzcotwel
wkwpgedgxvpltqid
mliajvpdocyzcbot
kqjhsipuibyjuref
zqdczykothbgxwsy
koirtljkuqzxioaz
audpjvhmqzvhzqas
cxyhxlhntyidldfx
iasgocejboxjgtkx
abehujmqotwcufxp
fmlrzqmazajxeedl
knswpkekbacuxfby
yvyalnvrxgstqhxm
sjnrljfrfuyqfwuw
ssaqruwarlvxrqzm
iaxbpeqqzlcwfqjz
uwyxshjutkanvvsc
uxwrlwbblcianvnb
nodtifgrxdojhneh
mloxjfusriktxrms
lkfzrwulbctupggc
gcrjljatfhitcgfj
tkdfxeanwskaivqs
ypyjxqtmitwubbgt
ssxbygzbjsltedjj
zdrsnoorwqfalnha
xlgmissaiqmowppd
azhbwhiopwpguiuo
fydlahgxtekbweet
qtaveuqpifprdoiy
kpubqyepxqleucem
wlqrgqmnupwiuory
rwyocktuqkuhdwxz
abzjfsdevoygctqv
zsofhaqqghncmzuw
lqbjwjqxqbfgdckc
bkhyxjkrqbbunido
yepxfjnnhldidsjb
builayfduxbppafc
wedllowzeuswkuez
gverfowxwtnvgrmo
tpxycfumxdqgntwf
lqzokaoglwnfcolw
yqsksyheyspmcdqt
vufvchcjjcltwddl
saeatqmuvnoacddt
dxjngeydvsjbobjs
ucrcxoakevhsgcep
cajgwjsfxkasbayt
hknzmteafsfemwuv
xxwhxwiinchqqudr
usfenmavvuevevgr
kxcobcwhsgyizjok
vhqnydeboeunnvyk
bgxbwbxypnxvaacw
bwjzdypacwgervgk
rrioqjluawwwnjcr
fiaeyggmgijnasot
xizotjsoqmkvhbzm
uzphtrpxwfnaiidz
kihppzgvgyoncptg
hfbkfrxwejdeuwbz
zgqthtuaqyrxicdy
zitqdjnnwhznftze
jnzlplsrwovxlqsn
bmwrobuhwnwivpca
uuwsvcdnoyovxuhn
nmfvoqgoppoyosaj
hxjkcppaisezygpe
icvnysgixapvtoos
vbvzajjgrmjygkhu
jinptbqkyqredaos
dpmknzhkhleawfvz
ouwwkfhcedsgqqxe
owroouiyptrijzgv
bewnckpmnbrmhfyu
evdqxevdacsbfbjb
catppmrovqavxstn
dqsbjibugjkhgazg
mkcldhjochtnvvne
sblkmhtifwtfnmsx
lynnaujghehmpfpt
vrseaozoheawffoq
ytysdzbpbazorqes
sezawbudymfvziff
vrlfhledogbgxbau
bipdlplesdezbldn
ermaenjunjtbekeo
eyaedubkthdecxjq
gbzurepoojlwucuy
rsiaqiiipjlouecx
beqjhvroixhiemtw
buzlowghhqbcbdwv
ldexambveeosaimo
fpyjzachgrhxcvnx
komgvqejojpnykol
fxebehjoxdujwmfu
jnfgvheocgtvmvkx
qmcclxxgnclkuspx
rsbelzrfdblatmzu
vexzwqjqrsenlrhm
tnfbkclwetommqmh
lzoskleonvmprdri
nnahplxqscvtgfwi
ubqdsflhnmiayzrp
xtiyqxhfyqonqzrn
omdtmjeqhmlfojfr
cnimgkdbxkkcnmkb
tapyijgmxzbmqnks
byacsxavjboovukk
awugnhcrygaoppjq
yxcnwrvhojpuxehg
btjdudofhxmgqbao
nzqlfygiysfuilou
nubwfjdxavunrliq
vqxmmhsbmhlewceh
ygavmcybepzfevrp
kgflmrqsvxprkqgq
iaqyqmcaedscmakk
cvbojnbfmrawxzkh
jjjrprbnlijzatuw
lcsudrrfnnggbrmk
qzgxbiavunawfibc
gnnalgfvefdfdwwg
nokmiitzrigxavsc
etzoxwzxqkkhvais
urxxfacgjccieufi
lqrioqhuvgcotuec
dydbaeyoypsbftra
hhrotenctylggzaf
evctqvzjnozpdxzu
tbpvithmorujxlcp
pllbtcbrtkfpvxcw
fzyxdqilyvqreowv
xdleeddxwvqjfmmt
fcldzthqqpbswoin
sgomzrpjfmvgwlzi
axjyskmtdjbxpwoz
hcvaevqxsmabvswh
lfdlsfcwkwicizfk
isjbwpzdognhoxvm
oqnexibqxlyxpluh
zqfbgodsfzwgcwuf
kvmnwruwsjllbldz
kghazimdyiyhmokj
uiktgpsxpoahofxn
zkdwawxargcmidct
ftbixlyiprshrjup
nofhmbxififwroeg
mcdaqrhplffxrcdt
fbjxnwojcvlawmlb
rizoftvwfdhiwyac
eduogrtyhxfwyars
zoikunqxgjwfqqwr
zxwbbpmvctzezaqh
nghujwyeabwdqnop
vcxamijpoyyksogn
jnckdbuteoqlsdae
jurfqqawafmsiqwv
inepmztrzehfafie
tznzkyvzodbrtscf
xewbavjeppflwscl
ucndzsorexjlnplo
jpxbctscngxgusvu
mfmygcllauzuoaok
oibkuxhjmhxhhzby
zjkslwagmeoisunw
avnnxmopdgvmukuu
jmaargejcwboqhkt
yacmpeosarsrfkrv
iqhgupookcaovwgh
ebjkdnxwtikqzufc
imdhbarytcscbsvb
ifyibukeffkbqvcr
aloighmyvwybtxhx
yszqwrutbkiwkxjg
xyholyzlltjhsuhp
gykhmrwucneoxcrf
badkdgqrpjzbabet
sunaucaucykwtkjj
pumqkglgfdhneero
usgtyuestahlydxq
xmfhflphzeudjsjm
knywgmclisgpootg
mtojnyrnvxtweuzb
uuxufbwfegysabww
vobhwwocqttlbsik
yuydfezeqgqxqmnd
wbqgqkwbibiilhzc
sfdmgxsbuzsawush
ilhbxcfgordyxwvp
ahqoavuysblnqaeg
plwgtvpgotskmsey
ewjcmzkcnautrrmp
tyekgzbznlikcyqj
bqzctiuaxpriuiga
bimvbfjkiupyqiys
mpqtbcxfhwymxncw
htemlptvqhharjgb
mqbsmsruwzzxgcxc
zjyedjwhnvteuaid
pzoelkoidwglpttc
efydnsvlfimvwxhx
gfyhgoeiyjcgfyze
deqtomhwopmzvjlt
casafubtkoopuaju
yylsfarntbucfulg
mgjwsormkjsrrxan
lkkenpupgmjpnqqd
tegweszyohsoluot
lihsfdwxmxvwdxna
rrefrjjxerphejwb
guuazonjoebhymtm
ysofqzmfmyneziki
lmjgaliatcpduoal
qzthcpjwtgahbebr
wvakvephyukmpemm
simxacxxzfoaeddw
aetgqmiqzxbvbviz
jxlmhdmqggevrxes
mmuglnjmuddzgaik
svopsqhtrslgycgc
xnvcsiiqrcjkvecn
kkvumxtvashxcops
bduflsdyeectvcgl
vfrxbwmmytjvqnsj
eeqtdneiyiaiofxw
crtbgknfacjtwkfl
uuutuoxdsxolpbhd
lcrztwzreaswovtn
htorkvnvujmjdqzj
wttzuzvrzlyhfzyf
oraewznfwgdsnhuk
rctlkqqvkwbgrcgk
cfehrsrqhzyiwtmz
kbvxwcumjkhvjpui
xxlocexbmniiakfo
gtknkkzvykmlqghl
kcjuxvkuimhwqrtk
vohekwkuyuoacuww
vorctgughscysyfo
zmjevqplngzswxyq
qhswdrhrijnatkyo
joakcwpfggtitizs
juzlwjijcmtswdtq
icbyaqohpkemhkip
rpdxgpzxncedmvzh
rozkmimbqhbhcddv
wkkypomlvyglpfpf
jcaqyaqvsefwtaya
ghvmtecoxlebdwnf
lqrcyiykkkpkxvqt
eqlarfazchmzotev
vqwndafvmpguggef
dbfxzrdkkrusmdke
cmjpjjgndozcmefj
hbrdcwjuyxapyhlo
mmforetykbosdwce
zynfntqwblbnfqik
sodwujfwlasznaiz
yyvrivjiqnxzqkfp
uldbskmmjbqllpnm
fyhhrmrsukeptynl
hpfjekktvdkgdkzl
bozhkoekcxzeorob
uvpptyfrzkvmtoky
hkhfprmjdpjvfkcb
igxzwktwsqhsivqu
qceomwysgkcylipb
cglateoynluyeqgc
xcsdfkpeguxgvpfh
owjhxlcncdgkqyia
rpbmrpcesiakqpna
lueszxiourxsmezb
zelvsowimzkxliwc
vzxbttoobtvdtkca
pfxvzphzwscqkzsi
edsjorainowytbzu
ipsegdaluoiphmnz
mkhueokfpemywvuw
urxdnumhylpafdlc
ggluurzavsxkvwkl
ctclphidqgteakox
tfobosynxsktajuk
jzrmemhxqmzhllif
eemwekimdfvqslsx
yjkwpzrbanoaajgq
rlxghzanuyeimfhx
hozbgdoorhthlqpv
obkbmflhyanxilnx
xojrippyxjmpzmsz
ukykmbfheixuviue
qivlmdexwucqkres
rmyxxipqkarpjmox
fgaftctbvcvnrror
raawxozucfqvasru
dinpjbdfjfizexdh
gybxubwnnbuyvjcr
qrqitdvyoneqyxcg
jqzcfggayzyoqteo
cikqpvxizpdbmppm
stfpldgyhfmucjjv
slzbcuihmimpduri
aufajwfrsorqqsnl
iylmzraibygmgmqj
lcdyfpcqlktudfmu
pmomzzsdpvgkkliw
zpplirgtscfhbrkj
mvhyerxfiljlotjl
ofkvrorwwhusyxjx
xngzmvcgkqfltjpe
yxfxaqipmysahqqq
sdqafdzgfdjuabup
qcqajmerahcdgxfv
xqimrqtupbapawro
qfvkqwidzzrehsbl
himixxvueksiqfdf
vgtfqpuzxxmhrvvd
adiioqeiejguaost
jnzxuycjxvxehbvm
xedbpxdhphamoodk
jsrioscmwlsfuxrg
mtsynnfxunuohbnf
enamqzfzjunnnkpe
uwcvfecunobyhces
ciygixtgbsccpftq
ewjgcronizkcsfjy
wztjkoipxsikoimv
jrgalyvfelwxforw
imylyalawbqwkrwb
yflwqfnuuvgjsgcj
wkysyzusldlojoue
zopllxnidcffcuau
bscgwxuprxaerskj
zvnvprxxjkhnkkpq
nejwxbhjxxdbenid
chryiccsebdbcnkc
guoeefaeafhlgvxh
nzapxrfrrqhsingx
mkzvquzvqvwsejqs
kozmlmbchydtxeeo
keylygnoqhmfzrfp
srwzoxccndoxylxe
uqjzalppoorosxxo
potmkinyuqxsfdfw
qkkwrhpbhypxhiun
wgfvnogarjmdbxyh
gkidtvepcvxopzuf
atwhvmmdvmewhzty
pybxizvuiwwngqej
zfumwnazxwwxtiry
keboraqttctosemx
vtlzxaqdetbhclib
wjiecykptzexuayl
ejatfnyjjdawepyk
mpcrobansyssvmju
gqukndzganeueabm
ukzscvomorucdnqd
wfydhtbzehgwfazx
mtwqdzlephqvxqmx
dltmlfxbjopefibh
atcfrowdflluqtbi
vowawlophlxaqonw
vblgdjzvwnocdipw
uzerzksmkvnlvlhm
ytjwhpaylohorvxd
siprvfxvnxcdgofz
cbhjupewcyjhvtgs
apqtozaofusmfqli
tmssrtlxfouowqnr
ntutrvwnzzgmokes
zrsgpwdzokztdpis
nrobvmsxtfmrqdhv
kadkaftffaziqdze
yrovbgcyqtlsnoux
modheiwuhntdecqs
gzhjypwddizemnys
gaputpwpcsvzxjho
bgmouxwoajgaozau
oxuapfrjcpyakiwt
kntwbvhuaahdixzj
epqjdjbnkxdnaccx
dspltdvznhypykri
tdrgqmbnagrxdwtt
njfqawzjggmemtbg
chpemsgwpzjpdnkk
fpsrobmbqbmigmwk
flxptsrqaazmprnl
nzdunrxlcbfklshm
miuwljvtkgzdlbnn
xbhjakklmbhsdmdt
xwxhsbnrwnegwcov
pwosflhodjaiexwq
fhgepuluczttfvqh
tldxcacbvxyamvkt
gffxatrjglkcehim
tzotkdrpxkucsdps
wxheftdepysvmzbe
qfooyczdzoewrmku
rvlwikuqdbpjuvoo
bcbrnbtfrdgijtzt
vaxqmvuogsxonlgq
ibsolflngegravgo
txntccjmqakcoorp
vrrbmqaxfbarmlmc
dzspqmttgsuhczto
pikcscjunxlwqtiw
lwzyogwxqitqfqlv
gsgjsuaqejtzglym
feyeqguxbgmcmgpp
gmttebyebdwvprkn
mzuuwbhzdjfdryxu
fganrbnplymqbzjx
cvsrbdcvhtxxdmro
scmgkjlkqukoamyp
fkgrqbyqpqcworqc
hjsrvkdibdjarxxb
sztzziuqroeidcus
pxdfvcpvwaddrzwv
phdqqxleqdjfgfbg
cqfikbgxvjmnfncy"

# ╔═╡ 0cc15cb8-4a0f-4a04-a858-43088cb68ff3
is_vowel(c::Char) = contains("aeiou", c)

# ╔═╡ 255d9b2e-4bc3-4b52-a835-9609eacf0c38
vowel_count(str::AbstractString) = sum(map(c -> is_vowel(c) ? 1 : 0, collect(str)))

# ╔═╡ 8cd32bc3-661a-4e0a-b02b-0a51b13c1e6a
# It contains at least three vowels (aeiou only), like aei, xazegov, or aeiouaeiouaeiou.
has_enough_vowels(str::AbstractString) = vowel_count(str) >= 3

# ╔═╡ 47de64ba-6e35-402b-b222-c214f2801e5f
@assert !has_enough_vowels("dvszwmarrgswjxmb") # is naughty because it contains only one vowel.

# ╔═╡ 9a2f0a15-57b5-4d28-838d-0bdd7485f500
@assert has_enough_vowels("ugknbfddgicrmopn") # is nice because it has at least three vowels (u...i...o...)

# ╔═╡ a1634cc1-6e18-46d8-97fd-fb2bc38131e5
@assert has_enough_vowels("aaa")

# ╔═╡ 273e3edb-417e-4d47-9450-ad07c41c6281
# It contains at least one letter that appears twice in a row, like xx, abcdde (dd), or aabbccdd (aa, bb, cc, or dd)
has_double_letter(str::AbstractString) = begin
	for (idx, c) in enumerate(str[1:end-1])
		if c == str[idx+1]
			return true
		end
	end
	false
end

# ╔═╡ 12123c64-a5d4-42e8-ac89-12d9f0b8cd6e
@assert !has_double_letter("jchzalrnumimnmhp") # is naughty because it has no double letter.

# ╔═╡ c81aed85-fd24-4f1c-a2a7-be7e8ab8ad40
@assert has_double_letter("ugknbfddgicrmopn") # is nice because it has ... a double letter (...dd...)

# ╔═╡ 8972c6bf-bcad-481a-b763-d9a5b22bd358
# It does not contain the strings ab, cd, pq, or xy, even if they are part of one of the other requirements.
no_forbidden(str::AbstractString) = begin
	for bad in ["ab", "cd", "pq", "xy"]
		if contains(str, bad)
			return false
		end
	end
	true
end

# ╔═╡ f768f13a-d075-4ffb-9d81-b2bb3c5b9bab
@assert !no_forbidden("haegwjzuvuyypxyu") # is naughty because it contains the string xy.

# ╔═╡ 8a320756-41fa-4f49-8c72-225065d21a59
@assert no_forbidden("ugknbfddgicrmopn") # is nice because it has ... none of the disallowed substrings.

# ╔═╡ 29f607dd-133e-4cc8-87d4-b101b4017897
check1(str::AbstractString) =
	has_enough_vowels(str) &&
	has_double_letter(str) &&
	no_forbidden(str)

# ╔═╡ 064293dc-63ff-44a6-9ecb-5a0fba08df66
@assert check1("ugknbfddgicrmopn")

# ╔═╡ 69765b8f-6917-4781-bfc4-3f00843005a2
@assert check1("aaa")

# ╔═╡ c4811f7d-6801-4fcb-8065-205bbd5e2e03
part1 = split(puzzle_input, "\n") |>
	Base.Fix1(map, str -> check1(str) ? 1 : 0) |>
	sum

# ╔═╡ 1ac94504-2c3f-4908-a1c8-9ab132cafb48
@test part1 == 238

# ╔═╡ 12eeff68-cf8d-48c0-a76a-1b6404b32247
md"Your puzzle answer was `238`."

# ╔═╡ 95825a3c-e112-4451-94c4-c68929e28cf3
md"""
# Part Two
Realizing the error of his ways, Santa has switched to a better model of determining whether a string is naughty or nice. None of the old rules apply, as they are all clearly ridiculous.

Now, a nice string is one with all of the following properties:

-   It contains a pair of any two letters that appears at least twice in the string without overlapping, like `xyxy` (`xy`) or `aabcdefgaa` (`aa`), but not like `aaa` (`aa`, but it overlaps).
-   It contains at least one letter which repeats with exactly one letter between them, like `xyx`, `abcdefeghi` (`efe`), or even `aaa`.

For example:

-   `qjhvhtzxzqqjkmpb` is nice because is has a pair that appears twice (`qj`) and a letter that repeats with exactly one letter between them (`zxz`).
-   `xxyxx` is nice because it has a pair that appears twice and a letter that repeats with one between, even though the letters used by each rule overlap.
-   `uurcxstgmygtbstg` is naughty because it has a pair (`tg`) but no repeat with a single letter between them.
-   `ieodomkazucvgmuy` is naughty because it has a repeating letter with one between (`odo`), but no pair that appears twice.

**How many strings are nice under these new rules?**
"""

# ╔═╡ 542ba6ce-f781-4f10-9414-329177b35733
# It contains a pair of any two letters that appears at least twice in the string without overlapping, like xyxy (xy) or aabcdefgaa (aa), but not like aaa (aa, but it overlaps).
has_pair_twice_without_overlapping(str::AbstractString) = begin
	for idx1 in 1:length(str)-3
		for idx2 in idx1+2:length(str)-1
			first = str[idx1:idx1+1]
			second = str[idx2:idx2+1]
			if first == second
				return true
			end
		end
	end
	false
end

# ╔═╡ 8f6b4850-737c-4194-a0ee-6f49c98a7f49
@assert has_pair_twice_without_overlapping("xyxy")

# ╔═╡ 065377a0-b116-43aa-9b10-da16c2651d73
@assert !has_pair_twice_without_overlapping("aaa")

# ╔═╡ c8ecfeed-df4a-4d3f-9ce2-4583362d61a2
@assert has_pair_twice_without_overlapping("qjhvhtzxzqqjkmpb") # is nice because is has a pair that appears twice (qj)

# ╔═╡ 192ba648-8047-4a19-b813-b94f64fdbe73
@assert !has_pair_twice_without_overlapping("ieodomkazucvgmuy") # is naughty because it has ... no pair that appears twice.

# ╔═╡ 3c5f7512-2cbb-4d27-8ff4-9109ebb113d0
# It contains at least one letter which repeats with exactly one letter between them, like xyx, abcdefeghi (efe), or even aaa.
has_repeat_with_gap(str::AbstractString) = begin
	for idx in 1:length(str)-2
		if str[idx] == str[idx+2]
			return true
		end
	end
	false
end

# ╔═╡ 799c6130-ff6d-4baa-9ab3-6701464af1d7
@assert has_repeat_with_gap("xyx")

# ╔═╡ db901ec7-d5b8-49d5-96af-83f11a51a357
@assert has_repeat_with_gap("abcdefeghi") # efe

# ╔═╡ ea996734-921d-48a3-bbd4-b62c7b40901b
@assert has_repeat_with_gap("aaa")

# ╔═╡ 2ab72150-e617-4d50-a744-3056b954624f
@assert !has_repeat_with_gap("uurcxstgmygtbstg") # is naughty because it has ... no repeat with a single letter between them.

# ╔═╡ 51e8ef23-7d4c-4a6b-a477-b7fd17a388dd
check2(str::AbstractString) =
	has_pair_twice_without_overlapping(str) &&
	has_repeat_with_gap(str)

# ╔═╡ 06344808-f1bc-44d9-818a-6f95a257ae91
part2 = split(puzzle_input, "\n") |>
	Base.Fix1(map, str -> check2(str) ? 1 : 0) |>
	sum

# ╔═╡ 8d322661-54b4-4672-bb21-b83be802983d
md"Your puzzle answer was `69`."

# ╔═╡ af6ef30c-c753-472f-8fde-312e277f1c89
@test part2 == 69

# ╔═╡ Cell order:
# ╠═f1b53f4d-6d7c-42b8-bd4b-b77e3b3a1055
# ╟─e9086930-afeb-11eb-35ad-dba62b2459a9
# ╟─c7430382-47c7-4f2e-a260-fa4de8a8648c
# ╠═0cc15cb8-4a0f-4a04-a858-43088cb68ff3
# ╠═255d9b2e-4bc3-4b52-a835-9609eacf0c38
# ╠═8cd32bc3-661a-4e0a-b02b-0a51b13c1e6a
# ╠═47de64ba-6e35-402b-b222-c214f2801e5f
# ╠═9a2f0a15-57b5-4d28-838d-0bdd7485f500
# ╠═a1634cc1-6e18-46d8-97fd-fb2bc38131e5
# ╠═273e3edb-417e-4d47-9450-ad07c41c6281
# ╠═12123c64-a5d4-42e8-ac89-12d9f0b8cd6e
# ╠═c81aed85-fd24-4f1c-a2a7-be7e8ab8ad40
# ╠═8972c6bf-bcad-481a-b763-d9a5b22bd358
# ╠═f768f13a-d075-4ffb-9d81-b2bb3c5b9bab
# ╠═8a320756-41fa-4f49-8c72-225065d21a59
# ╠═29f607dd-133e-4cc8-87d4-b101b4017897
# ╠═064293dc-63ff-44a6-9ecb-5a0fba08df66
# ╠═69765b8f-6917-4781-bfc4-3f00843005a2
# ╠═c4811f7d-6801-4fcb-8065-205bbd5e2e03
# ╠═1ac94504-2c3f-4908-a1c8-9ab132cafb48
# ╟─12eeff68-cf8d-48c0-a76a-1b6404b32247
# ╟─95825a3c-e112-4451-94c4-c68929e28cf3
# ╠═542ba6ce-f781-4f10-9414-329177b35733
# ╠═8f6b4850-737c-4194-a0ee-6f49c98a7f49
# ╠═065377a0-b116-43aa-9b10-da16c2651d73
# ╠═c8ecfeed-df4a-4d3f-9ce2-4583362d61a2
# ╠═192ba648-8047-4a19-b813-b94f64fdbe73
# ╠═3c5f7512-2cbb-4d27-8ff4-9109ebb113d0
# ╠═799c6130-ff6d-4baa-9ab3-6701464af1d7
# ╠═db901ec7-d5b8-49d5-96af-83f11a51a357
# ╠═ea996734-921d-48a3-bbd4-b62c7b40901b
# ╠═2ab72150-e617-4d50-a744-3056b954624f
# ╠═51e8ef23-7d4c-4a6b-a477-b7fd17a388dd
# ╠═06344808-f1bc-44d9-818a-6f95a257ae91
# ╟─8d322661-54b4-4672-bb21-b83be802983d
# ╠═af6ef30c-c753-472f-8fde-312e277f1c89


text:
tol = 2e-2
tolP = 1e-3
maxnit = 2000

@testset "String with ODE" begin

preconI = SaddleSearch.localPrecon(precon = [I], precon_prep! = (P, x) -> P)

heading1("TEST: String type methods with ode12r solver")

heading2("Muller potential")
V = MullerPotential()
x = ic_path(V, :near, 15)
E, dE = objective(V)
precon = x-> hessianprecond(V, x)

path = StaticString(0.0009, tol, maxnit, preconI, serial(), false,
                        verbose)
PATHx, PATHlog, _ = run!(path, E, dE, Path(x))
@test PATHlog[:maxres][end] <= path.tol

path = ODEString(reltol=0.1, tol = tol, maxnit = maxnit,
                        precon_scheme = preconI,
                        path_traverse = serial(),
                        verbose = verbose)
PATHx, PATHlog, _ = SaddleSearch.run!(path, E, dE, Path(x))
@test PATHlog[:maxres][end] <= path.tol

path = StaticNEB(0.0007, 0.0002, 3, tol, maxnit, preconI, serial(),
                        false, verbose)
PATHx, PATHlog, _ = run!(path, E, dE, Path(x))
@test PATHlog[:maxres][end] <= path.tol

path = ODENEB(reltol=1e-2, k=0.0002, interp=3, tol = tol,
                        maxnit = maxnit,
                        precon_scheme = preconI, path_traverse = serial(),
                        verbose = verbose)
PATHx, PATHlog, _ = run!(path, E, dE, Path(x))
@test PATHlog[:maxres][end] <= path.tol

heading2("Double well potential")
c = 16.0
V = DoubleWell(diagm([.5, c*c*(0.02^(c-1)), c]))
x = ic_path(V, :near, 11)
E, dE = objective(V)

P = copy(V.A); P[1] = 1.0
precon = x->[P]

path = StaticString(1./c, tol, maxnit, preconI, serial(), false,
                        verbose)
PATHx, PATHlog, _ = run!(path, E, dE, Path(x))
@test PATHlog[:maxres][end] <= path.tol

path = ODEString(reltol=1e-3, tol = tol, maxnit = maxnit,
                        precon_scheme = preconI, path_traverse = serial(),
                        verbose = verbose)
PATHx, PATHlog, _ = SaddleSearch.run!(path, E, dE, Path(x))
@test PATHlog[:maxres][end] <= path.tol

path = StaticNEB(0.07, 0.1, 3, tol, maxnit, preconI, serial(),
                        false, verbose)
PATHx, PATHlog, _ = run!(path, E, dE, Path(x))
@test PATHlog[:maxres][end] <= path.tol

path = ODENEB(reltol=1e-1, k=0.1, interp=3, tol = tol,
                        maxnit = maxnit,
                        precon_scheme = preconI, path_traverse = serial(),
                        verbose = verbose)
PATHx, PATHlog, _ = run!(path, E, dE, Path(x))
@test PATHlog[:maxres][end] <= path.tol

preconP = SaddleSearch.localPrecon(precon = [P], precon_prep! = (P, x) -> precon(x))

path = StaticString(0.25, tolP, maxnit, preconP, serial(), false,
                        verbose)
PATHx, PATHlog, _ = run!(path, E, dE, Path(x))
@test PATHlog[:maxres][end] <= path.tol

path = ODEString(reltol=0.1, tol = tolP, maxnit = maxnit,
                        precon_scheme = preconP, path_traverse = serial(),
                        verbose = verbose)
PATHx, PATHlog, _ = SaddleSearch.run!(path, E, dE, Path(x))
@test PATHlog[:maxres][end] <= path.tol

path = StaticNEB(0.6, 0.00001, 1, tolP, maxnit, preconP, serial(),
                        false, verbose)
PATHx, PATHlog, _ = run!(path, E, dE, Path(x))
@test PATHlog[:maxres][end] <= path.tol

path = ODENEB(reltol=0.1, k=0.00001, interp=1, tol = tolP,
                        maxnit = maxnit,
                        precon_scheme = preconP, path_traverse = serial(),
                        verbose = verbose)
PATHx, PATHlog, _ = run!(path, E, dE, Path(x))
@test PATHlog[:maxres][end] <= path.tol

heading2("Vacancy migration potential")
V = LJVacancy2D(R = 3.1)
x = ic_path(V, :min, 9)
E, dE = objective(V)

precon = x->[copy(precond(V, xn)) for xn in x]

path = StaticString(0.001, tol, maxnit, preconI, serial(), false,
                        verbose)
PATHx, PATHlog, _ = run!(path, E, dE, Path(x))
@test PATHlog[:maxres][end] <= path.tol

path = ODEString(reltol=1e-2, threshold=1e-5, tol = tol,
                        maxnit = maxnit,
                        precon_scheme = preconI,
                        path_traverse = serial(),
                        verbose = verbose)   # allow failure on 600 iterations
PATHx, PATHlog, _ = SaddleSearch.run!(path, E, dE, Path(x))
@show PATHlog[:maxres][end]
println("[allowed test failure: target is $(path.tol)]")
# @test PATHlog[:maxres][end] <= path.tol

path = StaticNEB(0.001, 0.01, 3, tol, maxnit, preconI, serial(),
                        false, verbose)
PATHx, PATHlog, _ = run!(path, E, dE, Path(x))
@test PATHlog[:maxres][end] <= path.tol

path = ODENEB(reltol=1e-2, k=0.00001, interp=3, tol = tol,
                        maxnit = maxnit,
                        precon_scheme = preconI, path_traverse = serial(),
                        verbose = verbose)
PATHx, PATHlog, _ = run!(path, E, dE, Path(x))
@test PATHlog[:maxres][end] <= path.tol

preconP = SaddleSearch.localPrecon(precon = precon(x),
            precon_prep! = (P, x) -> precon(x))

path = StaticString(1.55, tolP, maxnit, preconP, serial(), false,
                        verbose)
PATHx, PATHlog, _ = run!(path, E, dE, Path(x))
@test PATHlog[:maxres][end] <= path.tol

path = ODEString(reltol=1e-1, tol = tolP, maxnit = maxnit,
                        precon_scheme = preconP,
                        path_traverse = serial(),
                        verbose = verbose)
PATHx, PATHlog, _ = SaddleSearch.run!(path, E, dE, Path(x))
@test PATHlog[:maxres][end] <= path.tol

path = StaticNEB(1.0, 0.001, 3, tolP, maxnit, preconP, serial(),
                        false, verbose)
PATHx, PATHlog, _ = run!(path, E, dE, Path(x))
@test PATHlog[:maxres][end] <= path.tol

path = ODENEB(reltol=1e-2, k=0.001, interp=3, tol = tolP,
                        maxnit = maxnit,
                        precon_scheme = preconP, path_traverse = serial(),
                        verbose = verbose)
PATHx, PATHlog, _ = run!(path, E, dE, Path(x))
@test PATHlog[:maxres][end] <= path.tol

end


text:
#!/usr/bin/env julia

# Compute probabilities that feature-level splicing changes are occuring
# monotonically across time points.
#
# As input it it takes a YAML formatted file listing timepoints. Each timepoint
# is a just a list of condition names, which should match the condition names in
# the experiment file fed to isolator. E.g.
#
#     - [h7-day20, empty-vector]
#     - [adult, h7-1yr, let-7g-oe, fetal-ventricle, fetal-atrium, IMR90]
#
# Any condition that was present in the experiment but not present here will be
# ignored when considering monotonicity.
#
# This will print a table containing any feature with > 0.5 probability of
# changing across timepoints in a specific pattern.

using DataStructures, DataFrames, YAML

if length(ARGS) != 1
    println(STDERR,
    """
        Usage: isolator summarize condition-feature-splicing-samples isolator-output.h5 | ./monotonicity.jl timepoints.yml
    """)
    exit(1)
end

function read_permutations(timepoints_filename)
    equiv_classes = YAML.load_file(timepoints_filename)

    header = split(strip(readline(STDIN)), '\t')
    condition_names = [replace(name, r"_splice_rate_samples$", "") for name in header[7:end]]

    equiv = Dict{Int, Int}()
    for (i, equiv_class) in enumerate(equiv_classes)
        for cond in equiv_class
            equiv[findfirst(condition_names, cond)] = i
        end
    end

    k = maximum(values(equiv))
    permutations = Vector{Uint64}[]

    for i in 1:length(condition_names)
        if !haskey(equiv, i)
            equiv[i] = 0
        end
    end

    features = Any[]
    for line in eachline(STDIN)
        entry = split(strip(line), '\t')
        push!(features, entry[1:6])
        rows = [split(row, ',') for row in entry[7:end]]

        M = [parse(Float32, rows[i][j])
             for i in 1:length(rows), j in 1:length(rows[1])]

        permutations_row = Array(Uint64, length(rows[1]))

        for j in 1:size(M, 2)
            ps = sortperm(M[:,j])
            for i in 1:length(ps)
                ps[i] = equiv[ps[i]]
            end

            u = Uint64(0)
            for p in ps
                if p == 0
                    continue
                end
                u *= k
                u += p
            end

            permutations_row[j] = u
        end
        push!(permutations, permutations_row)
    end

    return permutations, features, k
end

permutations, features, k = read_permutations(ARGS[1])

function decode_permutation(p, k, n)
    decoded = Array(Uint64, n)
    for i in 1:n
        decoded[n-i+1] = p % k
        p = div(p, k)
    end
    return decoded
end

function interesting_transcripts(permutations, features)
    num_samples = length(permutations[1])
    interesting = Any[]
    cutoff = 0.5
    tally = DefaultDict(Uint64, Uint64, () -> 0)
    for i in 1:size(permutations, 1)
        for key in keys(tally)
            tally[key] = 0
        end

        for j in 1:num_samples
            tally[permutations[i][j]] += 1
        end

        for (k, v) in tally
            freq = v / num_samples
            if freq > cutoff
                push!(interesting, (features[i], k, freq))
            end
        end
    end

    return interesting
end

interesting = interesting_transcripts(permutations, features)

println("gene_names\tgene_ids\tincluding_transcript_ids\texcluding_transcript_ids\tlocus\tfeature_type\tprobability")
for feature in interesting
    println(join(feature[1], '\t'), '\t', feature[3])
end

max_stars_repo_name: dcjones/isolator

max_issues_repo_name: dcjones/isolator

max_forks_repo_name: dcjones/isolator


text:
@testset "asymptotic.jl" begin
    @testset "estimate, statistic, and kernel matrix" begin
        kernel = (ExponentialKernel() ∘ ScaleTransform(0.1)) ⊗ WhiteKernel()
        biasedskce = BiasedSKCE(kernel)
        unbiasedskce = UnbiasedSKCE(kernel)

        for nclasses in (2, 10, 100), nsamples in (10, 50, 100)
            # sample predictions and targets
            dist = Dirichlet(nclasses, 1)
            predictions = [rand(dist) for _ in 1:nsamples]
            targets_consistent = [
                rand(Categorical(prediction)) for prediction in predictions
            ]
            targets_onlyone = ones(Int, length(predictions))

            # compute calibration error estimate and test statistic
            for targets in (targets_consistent, targets_onlyone)
                estimate, statistic, kernelmatrix = CalibrationTests.estimate_statistic_kernelmatrix(
                    kernel, predictions, targets
                )

                @test estimate ≈ unbiasedskce(predictions, targets)
                @test statistic ≈
                      nsamples / (nsamples - 1) * unbiasedskce(predictions, targets) -
                      biasedskce(predictions, targets)
                @test kernelmatrix ≈
                      CalibrationErrors.unsafe_skce_eval.(
                    (kernel,),
                    predictions,
                    targets,
                    permutedims(predictions),
                    permutedims(targets),
                )
            end
        end
    end

    @testset "consistency" begin
        kernel1 = ExponentialKernel() ∘ ScaleTransform(0.1)
        kernel2 = WhiteKernel()
        kernel = kernel1 ⊗ kernel2
        αs = 0.05:0.1:0.95
        nsamples = 100

        pvalues_consistent = Vector{Float64}(undef, 100)
        pvalues_onlyone = similar(pvalues_consistent)

        for nclasses in (2, 10)
            dist = Dirichlet(nclasses, 1)
            predictions = [Vector{Float64}(undef, nclasses) for _ in 1:nsamples]
            targets_consistent = Vector{Int}(undef, nsamples)
            targets_onlyone = ones(Int, nsamples)

            for i in eachindex(pvalues_consistent)
                # sample predictions and targets
                for j in 1:nsamples
                    rand!(dist, predictions[j])
                    targets_consistent[j] = rand(Categorical(predictions[j]))
                end

                # define test
                test_consistent = AsymptoticSKCETest(
                    kernel, predictions, targets_consistent
                )
                test_onlyone = AsymptoticSKCETest(kernel, predictions, targets_onlyone)

                # estimate pvalues
                pvalues_consistent[i] = pvalue(test_consistent; bootstrap_iters=200)
                pvalues_onlyone[i] = pvalue(test_onlyone; bootstrap_iters=200)
            end

            # compute empirical test errors
            ecdf_consistent = ecdf(pvalues_consistent)
            @test all(ecdf_consistent(α) ≤ α + 0.1 for α in αs)
            @test all(p < 0.05 for p in pvalues_onlyone)
        end
    end
end


text:
# Run project tests

TestDir = @__DIR__
cd(TestDir)

include("Coeftab_plot/test_coeftab_plot_turing.jl")
#include("Pk_plot/test_pk_plot.jl")


text:
using Test

using BasisFunctions,
    CompactTranslatesDict,
    DomainIntegrals,
    DomainSets,
    FrameFun,
    StaticArrays

using IntegralEquations2D

include("test_qbf.jl")
include("test_helmholtz.jl")
include("test_laplace.jl")

# @testset "QBF" begin
#     test_qbf()
# end

@testset "Helmholtz BEM" begin
    test_helmholtz()
end
@testset "Laplace BEM" begin
    test_laplace()
end


text:
"""
### Logger
out::String
levels::Dict
log::Function \
A Logger logs information with different levels. Holds the function log(),
connected to the function _log(). Methods:
- log(::Int64, ::String)
- log(::String)
- log(::HTTP.Stream, ::String)
##### example
------------------
##### field info
- out::String - Logfile output directory.
- log::Function -
Writes to HTML console, and also logs at level 1 with logger.
- levels::Dict
------------------
##### constructors
Logger(levels::Dict{level_count::Int64 => crayon::Crayons.Crayon};
                    out::String = pwd() * "logs/log.txt")
Logger(; out::String = pwd() * "/logs/log.txt")
"""
mutable struct Logger <: ServerExtension
    type::Symbol
    out::String
    levels::Dict
    log::Function
    function Logger(levels::Dict; out::String = pwd() * "logs/log.txt")
        log(level::Int64, message::String) = _log(level, message, levels, out)
        log(message::String) = _log(1, message, levels, out)
        log(c::Connection, message::String) = _log(c.http, message)
        new(:connection, out::String, levels::Dict, log::Function)
    end
    function Logger(; out = pwd() * "/logs/log.txt")
        if contains(out, "src")
            out = pwd() * "../logs/log.txt"
        end
        levels::Dict{Integer, Crayon} = Dict(
        1 => Crayon(foreground = :light_cyan),
        2 => Crayon(foreground = :light_yellow),
        3 => Crayon(foreground = :yellow, bold = true),
        4 => Crayon(foreground = :red, bold = true))
        Logger(levels; out = out)
    end
end

"""
### _log(level::Int64, message::String, levels::Dict, out::String) -> _
------------------
Binded call for the field log() inside of Logger(). See ?(Logger) for more
    details on the field log. All arguments are fields of that type. Return is a
    printout into the REPL as well as an append to the log file, provided by the
    out URI.
------------------
### example (Closure from Logger)
log(level::Int64, message::String) = _log(level, message, levels, out)
log(message::String) = _log(1, message, levels, out)
"""
function _log(level::Int64, message::String, levels::Dict, out::String)
    time = now()
    if level > 1
        open(out, "a") do o
            try
                write(o, "[" * string(time) * "]: $message\n")
            catch
                try
                    touch(out)
                    write(o, "[" * string(time) * "]: $message\n")
                catch
                    throw(ArgumentError("Cannot access logs."))
                end
            end
        end
    end
    println(Crayon(foreground = :light_gray, bold = true), "[", levels[level],
     string(time), Crayon(foreground = :light_gray, bold = true), "]: ",
     message)
end

"""
### _log(http::HTTP.Stream, message::String) -> _
------------------
Binded call for the field log() inside of Logger(). This will log both to the
    JavaScript/HTML console
------------------
### example (Closure from Logger)
log(http::HTTP.Stream, message::String) = _log(http, message)
"""
function _log(http::HTTP.Stream, message::String)
    write(http, "<script>console.log('" * message * "');</script>")
end

"""
### route_from_dir(dir::String) -> ::Vector{String}
------------------
Recursively appends filenames for a directory AND all subsequent directories.
------------------
### example
x::Vector{String} = route_from_dir("mypath")
"""
function route_from_dir(dir::String)
    dirs::Vector{String} = readdir(dir)
    routes::Vector{String} = []
    for directory in dirs
        if isfile("$dir/" * directory)
            push!(routes, "$dir/$directory")
        else
            if ~(directory in routes)
                newread::String = dir * "/$directory"
                newrs::Vector{String} = route_from_dir(newread)
                [push!(routes, r) for r in newrs]
            end
        end
    end
    routes::Vector{String}
end

"""
### File <: Servable
dir::String
f::Function
------------------
##### field info
- dir::String - The directory of a file to serve.
- f::Function - Function whose output to be written to http().
------------------
##### constructors
File(dir::String)
"""
mutable struct File <: Servable
    dir::String
    f::Function
    function File(dir::String)
        f(c::Connection) = begin
        open(dir) do f
            write(c.http, f)
        end
        end
        new(dir, f)
    end
end
write!(c::Connection, f::File) = f.f(c)

"""
### Files
type::Symbol
directory::String
f::Function
------------------
- type::Symbol - The type of extension. There are three different selections
you can choose from. **:connection :routing :func**. A :connection extension
will be provided in Connection.extensions. A :routing function is passed a
Dict of routes as an argument. The last is a function argument, which is just a
specific function to run from the top-end to the server.
- directory::String - The directory to route the files from.
- f::Function - The function f() called with a Connection.
------------------
##### constructors
Files(dir::String)
"""
mutable struct Files <: ServerExtension
    type::Symbol
    directory::String
    f::Function
    function Files(directory::String = "public")
        f(r::Dict, e::Dict) = begin
            l = length(directory) + 1
            for path in route_from_dir(directory)
                push!(r, path[l:length(path)] => c::Connection -> write!(c, File(path)))
            end
        end
        new(:routing, directory, f)
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""JobCondition describes current state of a job.

    IoK8sApiBatchV1JobCondition(;
        lastProbeTime=nothing,
        lastTransitionTime=nothing,
        message=nothing,
        reason=nothing,
        status=nothing,
        type=nothing,
    )

    - lastProbeTime::IoK8sApimachineryPkgApisMetaV1Time : Last time the condition was checked.
    - lastTransitionTime::IoK8sApimachineryPkgApisMetaV1Time : Last time the condition transit from one status to another.
    - message::String : Human readable message indicating details about last transition.
    - reason::String : (brief) reason for the condition's last transition.
    - status::String : Status of the condition, one of True, False, Unknown.
    - type::String : Type of job condition, Complete or Failed.
"""
mutable struct IoK8sApiBatchV1JobCondition <: SwaggerModel
    lastProbeTime::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: lastProbeTime
    lastTransitionTime::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: lastTransitionTime
    message::Any # spec type: Union{ Nothing, String } # spec name: message
    reason::Any # spec type: Union{ Nothing, String } # spec name: reason
    status::Any # spec type: Union{ Nothing, String } # spec name: status
    type::Any # spec type: Union{ Nothing, String } # spec name: type

    function IoK8sApiBatchV1JobCondition(;lastProbeTime=nothing, lastTransitionTime=nothing, message=nothing, reason=nothing, status=nothing, type=nothing)
        o = new()
        validate_property(IoK8sApiBatchV1JobCondition, Symbol("lastProbeTime"), lastProbeTime)
        setfield!(o, Symbol("lastProbeTime"), lastProbeTime)
        validate_property(IoK8sApiBatchV1JobCondition, Symbol("lastTransitionTime"), lastTransitionTime)
        setfield!(o, Symbol("lastTransitionTime"), lastTransitionTime)
        validate_property(IoK8sApiBatchV1JobCondition, Symbol("message"), message)
        setfield!(o, Symbol("message"), message)
        validate_property(IoK8sApiBatchV1JobCondition, Symbol("reason"), reason)
        setfield!(o, Symbol("reason"), reason)
        validate_property(IoK8sApiBatchV1JobCondition, Symbol("status"), status)
        setfield!(o, Symbol("status"), status)
        validate_property(IoK8sApiBatchV1JobCondition, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        o
    end
end # type IoK8sApiBatchV1JobCondition

const _property_map_IoK8sApiBatchV1JobCondition = Dict{Symbol,Symbol}(Symbol("lastProbeTime")=>Symbol("lastProbeTime"), Symbol("lastTransitionTime")=>Symbol("lastTransitionTime"), Symbol("message")=>Symbol("message"), Symbol("reason")=>Symbol("reason"), Symbol("status")=>Symbol("status"), Symbol("type")=>Symbol("type"))
const _property_types_IoK8sApiBatchV1JobCondition = Dict{Symbol,String}(Symbol("lastProbeTime")=>"IoK8sApimachineryPkgApisMetaV1Time", Symbol("lastTransitionTime")=>"IoK8sApimachineryPkgApisMetaV1Time", Symbol("message")=>"String", Symbol("reason")=>"String", Symbol("status")=>"String", Symbol("type")=>"String")
Base.propertynames(::Type{ IoK8sApiBatchV1JobCondition }) = collect(keys(_property_map_IoK8sApiBatchV1JobCondition))
Swagger.property_type(::Type{ IoK8sApiBatchV1JobCondition }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiBatchV1JobCondition[name]))}
Swagger.field_name(::Type{ IoK8sApiBatchV1JobCondition }, property_name::Symbol) =  _property_map_IoK8sApiBatchV1JobCondition[property_name]

function check_required(o::IoK8sApiBatchV1JobCondition)
    (getproperty(o, Symbol("status")) === nothing) && (return false)
    (getproperty(o, Symbol("type")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiBatchV1JobCondition }, name::Symbol, val)
end


text:
"""
    ModelBaseEcon

This package is part of the StateSpaceEcon ecosystem.
It contains the basic elements needed for model definition.
StateSpaceEcon works with model objects defined with ModelBaseEcon.
"""
module ModelBaseEcon

using MacroTools
using SparseArrays
using DiffResults
using ForwardDiff
using Printf

# The Timer submodule
include("Timer.jl")
using .Timer
export @timer, inittimer, stoptimer, printtimer

# The Options submodule
include("Options.jl")

# The "misc" - various types and functions
include("misc.jl")

# NOTE: The order of inclusions matters.
include("abstract.jl")
include("evaluation.jl")
include("variables.jl")
include("equation.jl")
include("steadystate.jl")
include("parameters.jl")
include("metafuncs.jl")
include("model.jl")
include("export_model.jl")
include("linearize.jl")

"""
    @using_example name

Load models from the package examples/ folder.
The `@load_example` version is deprecated - stop using it now.
"""
macro using_example(name)
    examples_path = joinpath(dirname(pathof(@__MODULE__)), "..", "examples")
    return quote
        push!(LOAD_PATH, $(examples_path))
        using $(name)
        pop!(LOAD_PATH)
        $(name)
    end |> esc
end

" Deprecated. Use `@using_example` instead."
macro load_example(name)
    Base.depwarn("Use `@using_example` instead.", Symbol("@load_example"))
    return esc(:(@using_example $name))
end
export @using_example, @load_example

end # module


text:
module Testcb
    f1() = println("called f1()")
    f2() = exit(99)
    cf1 = @cfunction f1 Cvoid ()
    cf2 = @cfunction f2 Cvoid ()

    push!(Base.DL_LOAD_PATH, pwd())
    dl = Base.Libc.dlopen("libmtcb")
    dls = Base.Libc.Libdl.dlsym(dl, :cb_setup)
    ccall(dls, Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}), cf1, cf2)
end


text:
module CarterMethodBLCoords

using ..CarterBoyerLindquist
using DocStringExtensions

"""
    $(TYPEDSIGNATURES)

From Bardeen et al. (1972) eq. (2.10):

```math
T = E \\left( r^2 + a^2 \\right) - L * a.
```
"""
@inline T(E, L, r, a) = E * (r^2 + a^2) - L * a

"""
    $(TYPEDSIGNATURES)

From Bardeen et al. (1972) eq. (2.10), for the special case of a null-geodesic ``\\mu = 0``:

```math
V_r = T^2 - \\Delta \\left[ (L - a E)^2 + Q \\right]
```
"""
@inline Vr(E, L, M, Q, r, a) = T(E, L, r, a)^2 - Δ(M, r, a) * ((L - a * E)^2 + Q)

"""
    $(TYPEDSIGNATURES)

From Bardeen et al. (1972) eq. (2.10), for the special case of a null-geodesic ``\\mu = 0``:

```math
V_\\theta =
    Q + \\cos^2 (\\theta) \\left[ a^2 E^2 - \\frac{L^2}{\\sin^2 (\\theta) } \\right].
```
"""
@inline Vθ(E, L, Q, a, θ) = Q + cos(θ)^2 * ((a * E)^2 - (L / sin(θ))^2)

"""
    $(TYPEDSIGNATURES)

The ``t`` compontent of the equation of motion for a photon around a black hole, multiplied
by ``\\Sigma``.

From Bardeen et al. (1972) eq. (2.9d):

```math
\\Sigma \\frac{\\text{d}t}{\\text{d}\\lambda} =
- a \\left( a E \\sin^2 \\theta - L \\right)
+ \\frac{\\left( r^2 + a^2 \\right) T}{\\Delta}.
```
"""
@inline function Σδt_δλ(E, L, M, r, a, θ)
    -a * (a * E * sin(θ)^2 - L) + (r^2 + a^2) * T(E, L, r, a) / Δ(M, r, a)
end

"""
    $(TYPEDSIGNATURES)

The ``r`` compontent of the equation of motion for a photon around a black hole, multiplied
by ``\\Sigma``.

Modified from Bardeen et al. (1972) eq. (2.9a):

```math
\\Sigma \\frac{\\text{d}r}{\\text{d}\\lambda} =
\\pm \\sqrt{\\lvert V_r \\rvert},
```

where, for implementation reason, the sign is always positive. Instead, the sign is applied
in [`δ`](@ref).
"""
@inline function Σδr_δλ(E, L, M, Q, r, a)
    V = Vr(E, L, M, Q, r, a)
    √abs(V)
end

"""
    $(TYPEDSIGNATURES)

The ``\\theta`` compontent of the equation of motion for a photon around a black hole,
multiplied by ``\\Sigma``.

Modified from Bardeen et al. (1972) eq. (2.9b):

```math
\\Sigma \\frac{\\text{d}\\theta}{\\text{d}\\lambda} =
\\pm \\sqrt{\\lvert V_\\theta \\rvert},
```

where, for implementation reason, the sign is always positive. Instead, the sign is applied
in [`δ`](@ref).
"""
@inline function Σδθ_δλ(E, L, Q, a, θ)
    V = Vθ(E, L, Q, a, θ)
    √abs(V)
end

"""
    $(TYPEDSIGNATURES)

The ``\\phi`` compontent of the equation of motion for a photon around a black hole,
multiplied by ``\\Sigma``.

From Bardeen et al. (1972) eq. (2.9c):

```math
\\Sigma \\frac{\\text{d}\\phi}{\\text{d}\\lambda} =
- \\frac{L}{\\sin^2 \\theta} - aE
+ \\frac{aT}{\\Delta}.
```
"""
@inline function Σδϕ_δλ(E, L, M, r, a, θ)
    (L / sin(θ)^2) - (a * E) + a * T(E, L, r, a) / Δ(M, r, a)
end

"""
$(TYPEDSIGNATURES)

From Fanton et al. (1997), eq. (76):

```math
S = 1 + \\alpha \\omega \\sin \\theta.
```
"""
S(θ, α, ω) = 1 + α * ω * sin(θ)

"""
$(TYPEDSIGNATURES)

Calculates and returns the observer's angles ``\\sin \\Theta`` and ``\\sin \\Phi``, where
the parameters in the function signature correspond to the Bardeen et al. (1972) paper.

From Fanton et al. (1997), eq. (74) and eq. (75):

```math
\\begin{align*}
\\sin \\Theta &=
\\frac{\\alpha \\Sigma}{\\sqrt{A}}
\\left\\{
        \\beta^2 + \\left( \\alpha + a \\sin \\theta \\right)^2
        + \\frac{
            A S^2 - \\left( r^2 + a^2 + a \\alpha \\sin \\theta \\right)
        }{\\Delta}
\\right\\}, \\\\
\\sin \\Phi &=
-\\frac{\\alpha \\Sigma \\sqrt{\\Delta}}{A S \\sin\\Theta}.
\\end{align*}
```
"""
function sinΦsinΨ(Σ₀, sinθ, A₀, Δ₀, S₀, r, a, α, β)
    # calc 1
    sinΦ =
        ((α * Σ₀) / √A₀) /
        sqrt(β^2 + (α + a * sinθ)^2 + (A₀ * S₀^2 - (r^2 + a^2 + a * α * sinθ)^2) / Δ₀)

    # calc 2
    sinΨ = -(α * Σ₀ * √Δ₀) / (S₀ * A₀ * sinΦ)

    # return
    (sinΦ, sinΨ)
end

function sinΦsinΨ(M, r, a, θ, α, β)
    # value cache
    Σ₀ = Σ(r, a, θ)
    sinθ = sin(θ)
    A₀ = A(M, r, a, θ)
    Δ₀ = Δ(M, r, a)
    ω = 2.0 * a * r / A₀
    S₀ = S(θ, α, ω)

    sinΦsinΨ(Σ₀, sinθ, A₀, Δ₀, S₀, r, a, α, β)
end

"""
$(TYPEDSIGNATURES)

Calculates conserved quantities
- angular momentum ``L``
- Carter parameter ``Q``

for a photon described with position described by `x` in a Kerr spacetime given by
`p`.

From Fanton et al. (1997), eq. (69):

```math
L = \\frac{\\Upsilon_1}{\\Upsilon_2},
```

where

```math
\\begin{align*}
\\Upsilon_1 &= \\sin \\theta \\sin \\Phi \\sin \\Theta, \\\\
\\Upsilon_2 &= \\frac{\\Sigma \\sqrt{\\Delta}}{A} + \\omega \\Upsilon_1,
\\omega = \\frac{2 a r}{A},
\\end{align*}
```

taken from eq. (72) and eq. (73).

From Fanton et al. (1997), eq. (70):

```math
Q =
\\frac{P^2}{\\Delta}
- \\left( \\lambda - a \\right)^2
- \\frac{\\Sigma^2}{A} \\left( \\frac{\\cos \\Phi}{\\Upsilon_2} \\right)^2,
```

and

```math
P = \\left( r^2 + a^2 - a L \\right),
```

taken from eq. (71).
"""
function LQ(M, r, a, θ, sinΦ, sinΨ)
    sinθ = sin(θ)
    Σ₀ = Σ(r, a, θ)
    A₀ = A(M, r, a, θ)
    Δ₀ = Δ(M, r, a)
    ω = 2.0 * a * r / A₀

    Υ₁ = sinθ * sinΦ * sinΨ
    Υ₂ = (Σ₀ * √Δ₀ / A₀) + ω * Υ₁

    Υ₁ = sinθ * sinΦ * sinΨ
    Υ₂ = (Σ₀ * √Δ₀ / A₀) + ω * Υ₁

    L = Υ₁ / Υ₂
    P = (r^2 + a^2 - (a * L))
    Q = (P^2 / Δ₀) - (L - a)^2 - ((Σ₀^2 / A₀) * (cos(asin(sinΨ)) / Υ₂)^2)
    L, Q
end

"""
$(TYPEDSIGNATURES)

From Bardeen et al. (1972) eq. (2.3):

```math
\\Sigma = r^2 + a^2 \\cos^2( \\theta ).
```
"""
Σ(r, a, θ) = r^2 + (a * cos(θ))^2

"""
$(TYPEDSIGNATURES)

From Bardeen et al. (1972) eq. (2.3):

```math
\\Delta = r^2 - 2 M r + a^2.
```
"""
Δ(M, r, a) = r^2 - 2 * M * r + a^2

"""
$(TYPEDSIGNATURES)

From Bardeen et al. (1972) eq. (2.3):

```math
A = (r^2 + a^2)^2 - a^2 \\Delta \\sin^2 ( \\theta ).
```
"""
A(M, r, a, θ) = (r^2 + a^2)^2 - a^2 * Δ(M, r, a) * sin(θ)^2

"""
$(TYPEDSIGNATURES)

From Bardeen et al. (1972) eq. (2.21):

```math
Z_1 =
1 + \\sqrt[3]{1 - \\frac{a^2}{M^2}}
\\left[
    \\sqrt[3]{1 + \\frac{a}{M}} + \\sqrt[3]{1 - \\frac{a}{M}}
\\right].
```
"""
Z₁(M, a) = 1 + ∛(1 - (a / M)^2) * (∛(1 + (a / M)) + ∛(1 - (a / M)))

"""
$(TYPEDSIGNATURES)

From Bardeen et al. (1972) eq. (2.21):

```math
Z_2 = \\sqrt{\\frac{3a^2}{M^2} + Z_1^2}.
```
"""
Z₂(M, a) = √(3(a / M)^2 + Z₁(M, a)^2)

end # module


text:
# # Routines for reading and writing VTK files
# function read_vtk_2d(filepath::String)::UnstructuredMesh_2D
#     @info "Reading $filepath"
#     file = open(filepath, "r")
#     name = "DefaultMeshName"
#     local F
#     local points
#     local cells
#     local cell_types
#     while !eof(file)
#         line_split = split(readline(file))
#         if length(line_split) > 0
#             if line_split[1] == "#"
#                 if line_split[2] == "vtk"
#                     name = readline(file)
#                 end
#             elseif line_split[1] == "DATASET"
#                 if line_split[2] != "UNSTRUCTURED_GRID"
#                     @error "DATASET type is $(line_split[2]). Only UNSTRUCTURED_GRID is supported."
#                 end
#             elseif line_split[1] == "POINTS"
#                 if line_split[3] == "float"
#                     F = Float32
#                 elseif line_split[3] == "double"
#                     F = Float64
#                 else
#                     @error "Unable to identify POINTS data type."
#                 end
#                 points = read_vtk_points_2d(file, line_split[2], F)
#             elseif line_split[1] == "CELLS"
#                 cells = read_vtk_cells(file, line_split[2])
#             elseif line_split[1] == "CELL_TYPES"
#                 cell_types::Vector{Int64} = read_vtk_cell_types(file, line_split[2])
#             end
#         end
#     end
#     close(file)
#
#     # Remove all cells that are not 2d
#     # UnstructuredMesh (UM) uses the same cell types as VTK.
#     delete_indices = findall(x->x ∉  UnstructuredMesh_2D_cell_types, cell_types)
#     deleteat!(cell_types, delete_indices)
#     deleteat!(cells, delete_indices)
#
#     faces = Vector{Vector{Int64}}(undef, length(cells))
#     for i in eachindex(cell_types)
#         faces[i] = vcat(cell_types[i], cells[i])
#     end
#     # We can save a lot on memory and boost performance by choosing the smallest usable
#     # integer for our mesh. The relationship between faces, edges, and vertices is
#     # Triangle
#     #   F ≈ 2V
#     #   E ≈ 3F/2
#     # Quadrilateral
#     #   F ≈ V
#     #   E ≈ 2F
#     # Triangle6
#     #   2F ≈ V
#     #   E ≈ 3F/2
#     # Quadrilateral8
#     #   4F ≈ V
#     #   E ≈ 2F
#     # We see that V ≤ F ≤ 2E, so we use 2.2F as the max number of faces, edges, or vertices
#     # plus a fudge factor for small meshes, where the relationships become less accurate.
#     # If 2.2*length(faces) < typemax(UInt), convert to UInt
#     if ceil(2.2*length(faces)) < typemax(UInt16)
#         faces_16 = convert(Vector{Vector{UInt16}}, faces)
#         return UnstructuredMesh_2D{F, UInt16}(name = name,
#                                               points = points,
#                                               faces = [ SVector{length(f), UInt16}(f) for f in faces_16]
#                                              )
#     elseif ceil(2.2*length(faces)) < typemax(UInt32)
#         faces_32 = convert(Vector{Vector{UInt32}}, faces)
#         return UnstructuredMesh_2D{F, UInt32}(name = name,
#                                               points = points,
#                                               faces = [ SVector{length(f), UInt32}(f) for f in faces_32]
#                                              )
#     else
#         return UnstructuredMesh_2D{F, UInt64}(name = name,
#                                               points = points,
#                                               faces = [ SVector{length(f), UInt64}(f) for f in faces]
#                                              )
#     end
# end
#
# function read_vtk_points_2d(
#         file::IOStream,
#         npoints_string::SubString{String},
#         F::Type{<:AbstractFloat}
#     )
#     npoints = parse(Int64, npoints_string)
#     points = Vector{Point_2D{F}}(undef, npoints)
#     for i in 1:npoints
#         xyz = parse.(F, split(readline(file)))
#         points[i] = Point_2D(xyz[1], xyz[2])
#     end
#     return points
# end
#
# function read_vtk_cells(
#         file::IOStream,
#         ncells_string::SubString{String},
#     )
#     ncells = parse(Int64, ncells_string)
#     cells = Vector{Vector{Int64}}(undef, ncells)
#     for i in 1:ncells
#         # Strip the number of points and account for base 1 indexing
#         pointIDs = parse.(Int64, split(readline(file))) .+ 1
#         cells[i] = pointIDs[2:length(pointIDs)]
#     end
#     return cells
# end
#
# function read_vtk_cell_types(
#         file::IOStream,
#         ncells_string::SubString{String},
#     )
#     ncells = parse(Int64, ncells_string)
#     cell_types = Vector{Int64}(undef, ncells)
#     for i in 1:ncells
#         cellID = parse(Int64, readline(file))
#         cell_types[i] = cellID
#     end
#     return cell_types
# end
#
# function write_vtk_2d(filename::String, mesh::UnstructuredMesh_2D{F, U}) where {F <: AbstractFloat,
#                                                                                U <: Unsigned}
#     @info "Writing $filename"
#     @warn "  - Face sets not currently supported for vtk files"
#     # Check valid filename
#     if !occursin(".vtk", filename)
#         @error "Invalid filename. '.vtk' does not occur in $filename"
#     end
#     file = open(filename, "w")
#     println(file, "# vtk DataFile Version 2.0")
#     println(file, mesh.name)
#     println(file, "ASCII")
#     println(file, "DATASET UNSTRUCTURED_GRID")
#
#     # Points
#     if F === Float64
#         point_type = "double"
#     elseif F === Float32
#         point_type = "float"
#     else
#         @error "Unrecognized point type."
#     end
#     npoints = length(mesh.points)
#     println(file, "POINTS $npoints $point_type")
#     for i in 1:npoints
#         x, y = mesh.points[i]
#         println(file, "$x $y 0.0")
#     end
#     println(file, "")
#
#     # Cells
#     ncells = 0
#     ncell_parts = 0
#     for cell in mesh.faces
#         ncells += 1
#         ncell_parts += length(cell)
#     end
#     println(file, "CELLS $ncells $ncell_parts")
#     for cell in mesh.faces
#         nverts = length(cell) - 1
#         write(file, "$nverts ")
#         for i in 1:nverts
#             vert = cell[i + 1] - 1 # 0 based index
#             if i < nverts
#                 write(file, "$vert ")
#             else
#                 println(file, "$vert")
#             end
#         end
#     end
#     println(file, "")
#
#     # Cell types
#     # UnstructuredMesh uses the same cell types as VTK
#     println(file, "CELL_TYPES $ncells")
#     for cell in mesh.faces
#         cell_type = cell[1]
#         println(file, "$cell_type")
#     end
#
#     close(file)
# end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
# linear constrained quadratica programming filter EXAMPLE

# using JuMP
# using PyPlot
# using Cbc

# include("LCQP_filter.jl")
using Gadfly

include("../src/CDtypes.jl")
include("../src/initializations.jl")
include("../src/CDmainAlg.jl")
include("../src/GetResult.jl")
include("../src/GramMatrix.jl")
include("../src/InnerProducts.jl")
include("../src/Normalization.jl")
include("../src/utils.jl")

op=1
srand(10)
if op == 0
y=rand(500,1)
n=size(1)
    isempty(n)
elseif op == 1
  n=200
A=triu(ones(n,n))
lines=sort!(collect(Set(round((n-1)*rand(20,1)+1, 0))))
    y=A[:,lines]*(40*rand(size(lines)[1],1)+1)+10*(randn(n,1))
draw(
    SVG("/Users/mariosouto/Dropbox/SAM/L1_Adaptive_Trend_Filter/WindPower_CaseStudy/teste.svg", 7inch, 4.5inch),
    plot(y=y,Geom.line)
     )
end

y=y[:,1]

@time BCD, y_path, β_best, y_best, λ_best, γ_best = l1_adaptive_trend_filter(
  y, [1, 3], numλ=100, numγ=4, verbose=true, lower_bounds=[-Inf, -Inf, -Inf], upper_bounds=[Inf, Inf, 0]
  )

draw(
  SVG("/Users/mariosouto/Dropbox/SAM/L1_Adaptive_Trend_Filter/WindPower_CaseStudy/teste2.svg", 7inch, 4.5inch),
  plot(
    layer(x=1:length(y_best), y=y_best, Geom.line, Theme(default_color=color("blue"))),
    layer(x=1:length(y), y=y, Geom.line, Theme(default_color=color("red"))),
    Guide.xlabel("Hours"), Guide.ylabel("MW"), Guide.title("Wind Power Fit"),
    Coord.Cartesian(xmin=0,xmax=340)
    ))

# components
step = β_best[1]
# spike = β_best[2]
slope = β_best[3]

draw(
  SVG("/Users/mariosouto/Dropbox/SAM/L1_Adaptive_Trend_Filter/WindPower_CaseStudy/step.svg", 14inch, 8inch),
  plot(
    x=1:length(step), y=step, Geom.point, Geom.line,
    Guide.xlabel("Hours"), Guide.ylabel("Step size"),
    Coord.Cartesian(xmin=0,xmax=length(y))
    ))
# draw(
#   SVG("/Users/mariosouto/Dropbox/SAM/L1_Adaptive_Trend_Filter/WindPower_CaseStudy/spike.svg", 14inch, 8inch),
#   plot(x=1:length(spike),y=spike, Geom.point, Geom.line, Guide.xlabel("Hours")))
draw(
  SVG("/Users/mariosouto/Dropbox/SAM/L1_Adaptive_Trend_Filter/WindPower_CaseStudy/slope.svg", 14inch, 8inch),
  plot(x=1:length(slope),y=slope, Geom.point, Geom.line)
  )

max_stars_repo_name: joaquimg/ComponentFilterCD

max_issues_repo_name: joaquimg/ComponentFilterCD

max_forks_repo_name: joaquimg/ComponentFilterCD


text:
using TransformVariables
using Bijectors

struct Transform{N,T} <: Bijector{N}
    t::T

    function Transform(t::T) where {T <:  TransformVariables.TransformTuple}
        new{t.dimension, T}(t)
    end
end

function (b::Transform{N,T})(x) where {N,T <: TransformVariables.AbstractTransform}
    t(x)
end

# function B.dimension(b::Transform{N,T}) where {N, T <:  TransformVariables.AbstractTransform}

function (b::Inversed{<: Transform{N,T}})(y) where {N,T <:  TransformVariables.AbstractTransform}
    inverse(b.t, y)
end

#################################

export logpdf_with_trans

function logpdfℝⁿ(m::JointDistribution{A0,A,B,M}, data::NamedTuple, x::AbstractVector) where {A0,A,B,M}
    _logpdfℝⁿ(M, m.model, m.args, data, x)
end

function logpdfℝⁿ(m::JointDistribution{A0,A,B,M},x::AbstractVector) where {A0,A,B,M}
    _logpdfℝⁿ(M, m.model, m.args, x)
end

@gg M function _logpdfℝⁿ(_::Type{M}, _m::Model, _args, _data, _x) where M <: TypeLevel{Module}
    Expr(:let,
        Expr(:(=), :M, from_type(M)),
        type2model(_m) |> sourceLogpdfℝⁿ(_data) |> loadvals(_args, _data))
end

function sourceLogpdfℝⁿ(_data=NamedTuple())
    function(_m::Model)
        _m = canonical(_m)
        _datakeys = getntkeys(_data)

        proc(_m, st :: Assign)     = :($(st.x) = $(st.rhs))
        proc(_m, st :: Return)     = nothing
        proc(_m, st :: LineNumber) = nothing
        function proc(_m, st :: Sample)
            x = st.x
            xname = QuoteNode(x)

            x ∈ _datakeys && return :(_ℓ += logpdf($(st.rhs), $x))

            q = @q begin
                _j += 1
                $x = _x[_j]

                println($x)
                _ℓ += logpdf(transformed($(st.rhs)), $x)
            end

            return q
        end

        wrap(kernel) = @q begin
            _ℓ = 0.0
            _j = 0
            $kernel
            return _ℓ
        end

        buildSource(_m, proc, wrap) |> flatten
    end
end


text:
module TestPkg

using SparseArrays

export test

function test()
    a = spzeros(Float64, 5, 5)
    a[1, 3] = 4
    a[2, 5] = 19
    return a
end

end


text:
# Confirm that foh alt gives same results as foh

### Solver options ###
dt = 0.1
opts = TrajectoryOptimization.SolverOptions()
opts.square_root = false
opts.verbose = true
opts.cache = false
# opts.c1 = 1e-4
opts.c2 = 5.0
opts.cost_tolerance_intermediate = 1e-5
opts.constraint_tolerance = 1e-5
opts.cost_tolerance = 1e-5
opts.iterations_outerloop = 50
opts.iterations = 500
# opts.iterations_linesearch = 50
opts.constraint_decrease_ratio = 0.25
opts.outer_loop_update_type = :individual
######################

### Set up model, objective, solver ###
# Model, objective (unconstrained)
model, obj_uncon = TrajectoryOptimization.Dynamics.pendulum!

# -Constraints
u_min = -2
u_max = 2
x_min = [-20;-20]
x_max = [20; 20]

# # -Constrained objective
# function cE(cdot,x,u)
#     cdot[1] = x[1] - pi
# end

obj_con = ConstrainedObjective(obj_uncon, u_min=u_min, u_max=u_max, x_min=x_min, x_max=x_max)#, cE=cE)

# Solver
opts.use_static = false
solver = Solver(model,obj_con,integration=:rk3_foh,dt=dt,opts=opts)

# -Initial state and control trajectories
U = 10.0*rand(solver.model.m,solver.N)

results = ConstrainedVectorResults(solver.model.n,solver.model.m,solver.obj.p,solver.N)

copyto!(results.U,U)
rollout!(results,solver)
update_constraints!(results,solver)
update_jacobians!(results,solver)

results

# zoh version
# # Build the Hessian of the Lagrangian, stacked: constraints, Jacobians, multipliers
# n = solver.model.n
# m = solver.model.m
# p = solver.obj.p
# pI = solver.obj.pI
# N = solver.N
# Q = solver.obj.Q
# R = solver.obj.R
# Qf = solver.obj.Qf
# dt = solver.dt
#
# # initialize
# L = zeros((n+m)*(N-1)+n,(n+m)*(N-1)+n)
# c = zeros(p*(N-1)+n)
# cz = zeros(p*(N-1)+n,(n+m)*(N-1)+n)
# λ = zeros(p*(N-1)+n)
# μ = zeros(p*(N-1)+n,p*(N-1)+n)
#
# # get inequality indices
# idx_inequality = []
# tmp = [i for i = 1:pI]
#
# for i = 1:N-1
#     idx_inequality = cat(idx_inequality,tmp .+ (i-1)*p,dims=(1,1))
# end
# idx_inequality
#
# # assemble constraints
# for i = 1:N-1
#     c[(i-1)*p+1:(i-1)*p+p] = results.C[i]
#     cz[(i-1)*p+1:(i-1)*p+p,(i-1)*(n+m)+1:(i-1)*(n+m)+(n+m)] = [results.Cx[i] results.Cu[i]]
# end
# c[(N-1)*p+1:(N-1)*p+n] = results.CN
# cz[(N-1)*p+1:(N-1)*p+n,(N-1)*(n+m)+1:(N-1)*(n+m)+n] = results.Cx_N
#
# # assemble lagrange multipliers
# for i = 1:N-1
#     λ[(i-1)*p+1:(i-1)*p+p] = results.LAMBDA[i]
# end
# λ[(N-1)*p+1:(N-1)*p+n] = results.λN
#
# # assemble penalty matrix
# for i = 1:N-1
#     μ[(i-1)*p+1:(i-1)*p+p,(i-1)*p+1:(i-1)*p+p] = results.Iμ[i]
# end
# μ[(N-1)*p+1:(N-1)*p+n,(N-1)*p+1:(N-1)*p+n] = results.IμN
#
# if !solver.state.second_order_dual_update
#     # first order multiplier update
#     λ .= λ + μ*c
#     λ[idx_inequality] = max.(0.0,λ[idx_inequality])
# end
# # second order multiplier update
# constraint_status = ones(Bool,p*(N-1)+n)
#
# for i in idx_inequality
#     if c[i] <= 0.0 #&& λ[i] == 0.0
#         constraint_status[i] = false
#     end
# end
#
# idx_active = findall(x->x==true,constraint_status)
#
# # build Hessian
# for i = 1:N-1
#     L[(i-1)*(n+m)+1:(i-1)*(n+m)+(n+m),(i-1)*(n+m)+1:(i-1)*(n+m)+(n+m)] = [dt*Q zeros(n,m); zeros(m,n) dt*R]
# end
# L[(N-1)*(n+m)+1:(N-1)*(n+m)+n,(N-1)*(n+m)+1:(N-1)*(n+m)+n] = Qf
#
# L .+= cz[idx_active,:]'*μ[idx_active,idx_active]*cz[idx_active,:]
#
# B = cz[idx_active,:]*inv(L)*cz[idx_active,:]'
#
# if solver.state.second_order_dual_update
#     λ[idx_active] = λ[idx_active] + inv(B)*c[idx_active]
#     λ[idx_inequality] = max.(0.0,λ[idx_inequality])
# end
#
# # store results
# for i = 1:N-1
#     results.LAMBDA[i] = λ[(i-1)*p+1:(i-1)*p+p]
# end
# results.λN .= λ[(N-1)*p+1:(N-1)*p+n]

# foh case
# Build the Hessian of the Lagrangian, stacked: constraints, Jacobians, multipliers
n = solver.model.n
m = solver.model.m
q = n+m
p = solver.obj.p
pI = solver.obj.pI
N = solver.N
Q = solver.obj.Q
R = solver.obj.R
Qf = solver.obj.Qf
dt = solver.dt
if solver.model.m != length(results.U[1])
    m += n
    p += n
end

# initialize
L = zeros((n+m)*N,(n+m)*N)
c = zeros(p*N+n)
cz = zeros(p*N+n,(n+m)*N)
λ = zeros(p*N+n)
μ = zeros(p*N+n,p*N+n)

# get inequality indices
idx_inequality = []
tmp = [i for i = 1:pI]

for k = 1:N
    idx_inequality = cat(idx_inequality,tmp .+ (k-1)*p,dims=(1,1))
end
idx_inequality

# assemble constraints
for k = 1:N
    c[(k-1)*p+1:(k-1)*p+p] = results.C[k]
    cz[(k-1)*p+1:(k-1)*p+p,(k-1)*(n+m)+1:(k-1)*(n+m)+(n+m)] = [results.Cx[k] results.Cu[k]]
end
c[N*p+1:N*p+n] = results.CN
cz[N*p+1:N*p+n,(N-1)*(n+m)+1:(N-1)*(n+m)+n] = results.Cx_N

# assemble lagrange multipliers
for k = 1:N
    λ[(k-1)*p+1:(k-1)*p+p] = results.LAMBDA[k]
end
λ[N*p+1:N*p+n] = results.λN

# assemble penalty matrix
for k = 1:N
    μ[(k-1)*p+1:(k-1)*p+p,(k-1)*p+1:(k-1)*p+p] = results.Iμ[k]
end
μ[N*p+1:N*p+n,N*p+1:N*p+n] = results.IμN

if !solver.state.second_order_dual_update
    # first order multiplier update
    λ .= λ + μ*c
    λ[idx_inequality] = max.(0.0,λ[idx_inequality])
end
# second order multiplier update
constraint_status = ones(Bool,N*p+n)

for i in idx_inequality
    if c[i] <= 0.0 #&& λ[i] == 0.0
        constraint_status[i] = false
    end
end

idx_active = findall(x->x==true,constraint_status)

# build Hessian
for k = 1:N-1
    # Unpack Jacobians, ̇x
    Ac1, Bc1 = results.Ac[k], results.Bc[k]
    Ac2, Bc2 = results.Ac[k+1], results.Bc[k+1]
    Ad, Bd, Cd = results.fx[k], results.fu[k], results.fv[k]

    xm = results.xmid[k]
    um = (U[k] + U[k+1])/2.0

    lxx = dt/6*Q + 4*dt/6*(I/2 + dt/8*Ac1)'*Q*(I/2 + dt/8*Ac1)
    luu = dt/6*R + 4*dt/6*((dt/8*Bc1)'*Q*(dt/8*Bc1) + 0.5*R*0.5)
    lyy = dt/6*Q + 4*dt/6*(I/2 - dt/8*Ac2)'*Q*(I/2 - dt/8*Ac2)
    lvv = dt/6*R + 4*dt/6*((-dt/8*Bc2)'*Q*(-dt/8*Bc2) + 0.5*R*0.5)

    lxu = 4*dt/6*(I/2 + dt/8*Ac1)'*Q*(dt/8*Bc1)
    lxy = 4*dt/6*(I/2 + dt/8*Ac1)'*Q*(I/2 - dt/8*Ac2)
    lxv = 4*dt/6*(I/2 + dt/8*Ac1)'*Q*(-dt/8*Bc2)
    luy = 4*dt/6*(dt/8*Bc1)'*Q*(I/2 - dt/8*Ac2)
    luv = 4*dt/6*((dt/8*Bc1)'*Q*(-dt/8*Bc2) + 0.5*R*0.5)
    lyv = 4*dt/6*(I/2 - dt/8*Ac2)'*Q*(-dt/8*Bc2)

    L[(k-1)*(n+m)+1:(k-1)*(n+m)+2*(n+m),(k-1)*(n+m)+1:(k-1)*(n+m)+2*(n+m)] = [lxx lxu lxy lxv; lxu' luu luy luv; lxy' luy' lyy lyv; lxv' luv' lyv' lvv]
end
L[(N-1)*(n+m)+1:(N-1)*(n+m)+n,(N-1)*(n+m)+1:(N-1)*(n+m)+n] = Qf

L .+= cz[idx_active,:]'*μ[idx_active,idx_active]*cz[idx_active,:]

B = cz[idx_active,:]*inv(L)*cz[idx_active,:]'

if solver.state.second_order_dual_update
    λ[idx_active] = λ[idx_active] + inv(B)*c[idx_active]
    λ[idx_inequality] = max.(0.0,λ[idx_inequality])
end

# store results
for k = 1:N
    results.LAMBDA[k] = λ[(k-1)*p+1:(k-1)*p+p]
end
results.λN .= λ[N*p+1:N*p+n]


text:
# The Computer Language Benchmarks Game
# https://salsa.debian.org/benchmarksgame-team/benchmarksgame/

# based on an OCaml program
# *reset*

using Printf

abstract type BTree end

mutable struct Empty <: BTree
end

mutable struct Node <: BTree
    left::BTree
    right::BTree
end

function make(d)
    if d == 0
        Node(Empty(), Empty())
    else
        Node(make(d-1), make(d-1))
    end
end

check(t::Empty) = 0
check(t::Node) = 1 + check(t.left) + check(t.right)

function loop_depths(d, min_depth, max_depth)
    for i = 0:div(max_depth - d, 2)
        niter = 1 << (max_depth - d + min_depth)
        c = 0
        for j = 1:niter
            c += check(make(d))
        end
        @printf("%i\t trees of depth %i\t check: %i\n", niter, d, c)
        d += 2
    end
end

function perf_binary_trees(N::Int=10)
    min_depth = 4
    max_depth = N
    stretch_depth = max_depth + 1

    # create and check stretch tree
    let c = check(make(stretch_depth))
        @printf("stretch tree of depth %i\t check: %i\n", stretch_depth, c)
    end

    long_lived_tree = make(max_depth)

    loop_depths(min_depth, min_depth, max_depth)
    @printf("long lived tree of depth %i\t check: %i\n", max_depth, check(long_lived_tree))

end

n = parse(Int,ARGS[1])
perf_binary_trees(n)


text:
export PrerecordedDriver

type PrerecordedDriver <: DriverModel{NextState, IntegratedContinuous}
    trajdata::Trajdata # log we pull from
    time_start::Float64   # starting time [s]
    t_current::Float64
    timestep::Float64
    nextstate::VehicleState
end
function PrerecordedDriver(trajdata::Trajdata, frame_start::Int, egoid::Int, timestep::Float64=get_mean_timestep(trajdata))
    time_start = get_time(trajdata, frame_start)
    t_current = NaN
    nextstate = get_vehiclestate(trajdata, egoid, frame_start)
    PrerecordedDriver(trajdata, time_start, t_current, timestep, nextstate)
end

get_name(::PrerecordedDriver) = "PrerecordedDriver"
action_context(model::PrerecordedDriver) = IntegratedContinuous(model.timestep, 1)
function reset_hidden_state!(model::PrerecordedDriver)
    model.t_current = model.time_start
    model
end
function observe!(model::PrerecordedDriver, scene::Scene, roadway::Roadway, egoid::Int)
    model.t_current += model.timestep

    i = searchsortedfirst(model.trajdata.frames, model.t_current, by=frame_or_time->isa(frame_or_time, TrajdataFrame) ? frame_or_time.t : frame_or_time) # index of first frame ≥ t_current

    if i == 1
        lo, hi = 1, 2
    elseif i ≤ nframes(model.trajdata)
        lo, hi = i-1, i
    else # t_current > all times in model.trajdata.frames
        lo, hi = i-2, i-1
    end

    t = model.t_current
    t_lo = get_time(model.trajdata, lo)
    t_hi = get_time(model.trajdata, hi)

    γ = (t - t_lo) / (t_hi - t_lo)

    if (γ == 0.0 && !iscarinframe(model.trajdata, egoid, lo)) ||
       (γ == 1.0 && !iscarinframe(model.trajdata, egoid, hi)) ||
       (!iscarinframe(model.trajdata, egoid, lo) || !iscarinframe(model.trajdata, egoid, hi))

        model.nextstate = VehicleState()
    else
        state_lo = get_vehiclestate(model.trajdata, egoid, lo)
        state_hi = get_vehiclestate(model.trajdata, egoid, hi)
        model.nextstate = lerp(state_lo, state_hi, γ, model.trajdata.roadway)
    end

    model
end

Base.rand(model::PrerecordedDriver) = NextState(model.nextstate)
Distributions.pdf(model::PrerecordedDriver, a::NextState) = a.s == model.nextstate ? Inf : 0.0
Distributions.logpdf(model::PrerecordedDriver, a::NextState) = a.s == model.nextstate ? Inf : -Inf


text:
module MomentConstraintTest

using Test
using ForneyLab
using ForneyLab: outboundType, isApplicable
using ForneyLab: SPMomentConstraintOutG
using ForneyLab: expectedValue, constrainedMoments

f(x::Float64) = x # Trivial
g(x::Float64) = 1.0*(x > 0.0) # Non-trivial

#--------
# Helpers
#--------

@testset "expectedValue" begin
    @test expectedValue(2.0, 0.0, 1.0, f) == 1.999999026302908
    @test expectedValue(2.0, 0.0, 1.0, g) == 0.8807970779778822
end

@testset "constrainedMoments" begin
    @test constrainedMoments(2.0, 0.0, 1.0, f) == (1.999999026302908, 0.9999990688765581)
    @test constrainedMoments(2.0, 0.0, 1.0, g) == (0.6155633430596962, 0.6210829316309023)
end

#-------------
# Update rules
#-------------

@testset "SPMomentConstraintOutG" begin
    @test SPMomentConstraintOutG <: SumProductRule{MomentConstraint}
    @test outboundType(SPMomentConstraintOutG) == Message{GaussianWeightedMeanPrecision}
    @test isApplicable(SPMomentConstraintOutG, [Nothing])

    @test ruleSPMomentConstraintOutG(Message(Univariate, GaussianMeanVariance, m=0.0, v=1.0), g, 0.95, 0.5) == Message(Univariate, GaussianWeightedMeanPrecision, xi=1.5449452063991789, w=1.123838185932649)
    @test ruleSPMomentConstraintOutG(Message(Univariate, GaussianMeanVariance, m=0.0, v=1.0), f, 2.0, 0.0) == Message(Univariate, GaussianWeightedMeanPrecision, xi=2.0000018622402607, w=9.31120130553964e-7)
end

@testset "averageEnergy" begin
    @test averageEnergy(MomentConstraint, ProbabilityDistribution(GaussianMeanVariance, m=0.0, v=1.0)) == 0.0
end

#------------
# Integration
#------------

@testset "MomentConstraint integration" begin
    fg = FactorGraph()
    @RV x ~ GaussianMeanVariance(0.0, 1.0)
    MomentConstraint(x; g=f, G=1.0, eta_init=0.0)
    pfz = PosteriorFactorization(fg)
    algo = messagePassingAlgorithm(x)
    algo_code = algorithmSourceCode(algo)

    @test occursin("ruleSPMomentConstraintOutG(messages[1], f, 1.0, 0.0)", algo_code)
end

end # module